# Several models identify different which genes are significant

I want to find some good predictors (genes). 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 [1] 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)
[1] 0.6153846
> specificity(conf_matrix)
[1] 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

• There isn't a single significant value for a gene. It depends on the models and data used. Here with several models you get different values. So you need to decide which is the right model (and why) and act accordingly. Please clarify what do you mean by "horrible confusion" – llrs Jan 15 '19 at 11:40
• This is really going into pure statistics and at this point you should sit down with a local statistician who can explain how coefficients and probabilities in the various models are related to each other (e.g., they genes themselves have some correlation and increasing numbers of coefficients will suck up degrees of freedom). – Devon Ryan Jan 15 '19 at 12:49
• Thanks a lot, I have checked the Spearman correlation of the expression of each of 23 genes with dependent variable (responders Vs non-responders and I saw only the 5 genes identified initially showed significant correlation with dependent variable. However I am not sure as no one here to help me face to face at the time. – Angel Jan 15 '19 at 12:53