3
$\begingroup$

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
[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

$\endgroup$
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.