# How do I resolve disagreements in the determination of significant genes?

I want to find some good predictors (genes). I have checked the Spearman correlation of the expression of each of 23 genes with dependent variables (responders Vs non-responders and I saw only the 5 genes identified initially showed significant correlation with the dependent variable.

This is my data, log transformed RNA-seq:

              TRG    CDK6 EGFR  KIF2C CDC20
Sample 1  TRG12  11.39 10.62  9.75 10.34
Sample 2  TRG12  10.16  8.63  8.68  9.08
Sample 3  TRG12   9.29 10.24  9.89 10.11
Sample 4  TRG45  11.53  9.22  9.35  9.13
Sample 5  TRG45   8.35 10.62 10.25 10.01
Sample 6  TRG45  11.71 10.43  8.87  9.44
...


23 differentially expressed genes between responder patients (TRG12) and non-responder patients (TRG45).

1- I tested each of 23 genes individually in this code and each of them gives p-value < 0.05 remained as a good predictor; For example for CDK6 I have done

> glm=glm(TRG ~ CDK6, data = df, family = binomial(link = 'logit'))
>
> summary(glm)

Call:
glm(formula = TRG ~ CDK6, family = binomial(link = "logit"),
data = df)

Deviance Residuals:
Min       1Q   Median       3Q      Max
-1.8286  -0.9317   0.6382   0.9500   1.5966

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)  26.8950     9.6627   2.783  0.00538 **
CDK6         -2.6932     0.9765  -2.758  0.00582 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 77.347  on 55  degrees of freedom
Residual deviance: 64.062  on 54  degrees of freedom
AIC: 68.062

Number of Fisher Scoring iterations: 5


Finally I obtained five genes and I put them in this model:

final <- glm(TRG ~ CDK6 + CXCL8 + IL6 + ISG15 + PTGS2 , data = df, family = binomial(link = 'logit'))


but

> summary(final)

Call:
glm(formula = TRG ~ CDK6 + CXCL8 + IL6 + ISG15 + PTGS2, family = binomial(link = "logit"),
data = df)

Deviance Residuals:
Min       1Q   Median       3Q      Max
-1.5696  -0.9090   0.2233   0.7125   2.0131

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)   3.2258    13.3877   0.241   0.8096
CDK6         -2.6101     1.0860  -2.403   0.0162 *
CXCL8         0.8702     0.9783   0.889   0.3737
IL6           0.9591     0.8883   1.080   0.2803
ISG15         0.9574     0.5101   1.877   0.0605 .
PTGS2        -0.4623     1.1005  -0.420   0.6744
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 77.347  on 55  degrees of freedom
Residual deviance: 53.917  on 50  degrees of freedom
AIC: 65.917

Number of Fisher Scoring iterations: 6


I am just seeing CDK6 as significant. I done cross validation

> cv.glm(df, final, K=nrow(df))$delta  0.2130050 0.2125748  Likely the delta values are not greatly differ (Model is right ???) Then I tried Sensitivity and Specificity > threshold=0.5 > predicted_values<-ifelse(predict(final,type="response")>threshold,1,0) > actual_values<-final$y
> conf_matrix<-table(predicted_values,actual_values)
> conf_matrix
actual_values
predicted_values  0  1
0 16  7
1 10 23
> sensitivity(conf_matrix)
 0.6153846
> specificity(conf_matrix)
 0.7666667


2- I went this way for step-wise method

full.model <- glm(TRG ~., data = df, family = binomial)

library(MASS)
step.model <- full.model %>% stepAIC(trace = FALSE)
coef(step.model)

> coef(step.model)
(Intercept)        CDK6        CHGA        CSF3       CXCL8       ERBB2         IL6       ISG15       KRT14       KRT17
-23770.6261   -469.3660    376.3528   -238.5697   1404.5905    590.2586    789.5295    842.2058    423.0402   -651.3403
MAGEA1         OSM       PTGS2
496.6933   -637.0583   -500.1040


In addition to 5 genes another genes remains as predictors but likely none of them are significant by Pr(>|z|) > 0.05

> summary(step.model)

Call:
glm(formula = TRG ~ CDK6 + CHGA + CSF3 + CXCL8 + ERBB2 + IL6 +
ISG15 + KRT14 + KRT17 + MAGEA1 + OSM + PTGS2, family = binomial,
data = df)

Deviance Residuals:
Min          1Q      Median          3Q         Max
-1.317e-04  -2.100e-08   2.100e-08   2.100e-08   1.396e-04

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)  -23770.6  2922188.9  -0.008    0.994
CDK6           -469.4    64920.6  -0.007    0.994
CHGA            376.4    46191.2   0.008    0.993
CSF3           -238.6    30857.1  -0.008    0.994
CXCL8          1404.6   172558.7   0.008    0.994
ERBB2           590.3    71989.5   0.008    0.993
IL6             789.5   100130.7   0.008    0.994
ISG15           842.2   102934.0   0.008    0.993
KRT14           423.0    52565.5   0.008    0.994
KRT17          -651.3    79489.6  -0.008    0.993
MAGEA1          496.7    60319.6   0.008    0.993
OSM            -637.1    82121.8  -0.008    0.994
PTGS2          -500.1    73541.7  -0.007    0.995

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 7.7347e+01  on 55  degrees of freedom
Residual deviance: 1.4741e-07  on 43  degrees of freedom
AIC: 26

Number of Fisher Scoring iterations: 25


3- I went to lasso regression

> library(glmnet)
> lassoModel <- glmnet(
+   x=data.matrix(df[,-1]),
+   y=df$$TRG, + standardize=TRUE, + alpha=1.0, + family="multinomial") > #Derive coefficients for each gene to each subtype > co <- coef(lassoModel, s=idealLambda, exact=TRUE) > co$$TRG12
24 x 1 sparse Matrix of class "dgCMatrix"
1
(Intercept)  4.441322326
ALB         -0.010867331
AQP9        -0.051908908
CALML5       0.126359364
CCL4        -0.059869284
CDK6         0.283781437
CHGA        -0.088283426
CREB3L3      0.109988758
CSF3        -0.012311176
CXCL5       -0.030532485
CXCL6       -0.063610112
CXCL8       -0.084385732
ERBB2       -0.102303095
FGFRL1       0.148759702
IL1B        -0.072225606
IL6         -0.122671026
ISG15       -0.297386587
KLK5         0.039891341
KRT14       -0.004536929
KRT17        0.041693665
MAGEA1      -0.154075458
OSM          0.028987648
PTGS2       -0.034960405
S100A7A     -0.008190105

$TRG45 24 x 1 sparse Matrix of class "dgCMatrix" 1 (Intercept) -4.441322326 ALB 0.010867331 AQP9 0.051908908 CALML5 -0.126359364 CCL4 0.059869284 CDK6 -0.283781437 CHGA 0.088283426 CREB3L3 -0.109988758 CSF3 0.012311176 CXCL5 0.030532485 CXCL6 0.063610112 CXCL8 0.084385732 ERBB2 0.102303095 FGFRL1 -0.148759702 IL1B 0.072225606 IL6 0.122671026 ISG15 0.297386587 KLK5 -0.039891341 KRT14 0.004536929 KRT17 -0.041693665 MAGEA1 0.154075458 OSM -0.028987648 PTGS2 0.034960405 S100A7A 0.008190105 >  Likely none of genes is significant so I just tried 5 genes from first glm > finalLasso <- glm(df$TRG ~ CDK6 + IL6 + ISG15 + CXCL8 + CALML5 , data=df, family=binomial(link="logit"))
> summary(finalLasso)

Call:
glm(formula = df\$TRG ~ CDK6 + IL6 + ISG15 + CXCL8 + CALML5, family = binomial(link = "logit"),
data = df)

Deviance Residuals:
Min       1Q   Median       3Q      Max
-1.6780  -0.6733   0.1643   0.6925   1.9182

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)  -0.4178    14.7456  -0.028   0.9774
CDK6         -2.2672     1.1537  -1.965   0.0494 *
IL6           1.2524     0.9197   1.362   0.1733
ISG15         1.1763     0.5479   2.147   0.0318 *
CXCL8         0.5413     0.6443   0.840   0.4008
CALML5       -0.7890     0.4219  -1.870   0.0615 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 77.347  on 55  degrees of freedom
Residual deviance: 49.839  on 50  degrees of freedom
AIC: 61.839

Number of Fisher Scoring iterations: 6

>


That 2 genes are significant; really confusing :(

Please somebody save me from this horrible confusion, finally which genes are significant remaining in model and which not

First confusion: Why when I am putting each of 23 genes individually in glm model, 5 of them returns p-value < 0.05 but when summary of final model containing these 5 genes says only p-value of one of them is < 0.05

Second confusion: When I am using stepwise regression, summary of model says none of my genes is significant while some of them had been remained as predictors by coef(step.model)

Third confusion: CDK6 as my significant genes in two of model shows negative coefficient while I know that show be correlated with dependent variable positively