3
$\begingroup$

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

$\endgroup$
  • 1
    $\begingroup$ 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" $\endgroup$ – llrs Jan 15 '19 at 11:40
  • 2
    $\begingroup$ 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). $\endgroup$ – Devon Ryan Jan 15 '19 at 12:49
  • $\begingroup$ 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. $\endgroup$ – Angel Jan 15 '19 at 12:53
  • 2
    $\begingroup$ We can't substitute a proper advisor or professor. We can help you with specific problems or give advice to your research but we cannot take those decisions (and responsibility) from you. Personally (and this might come as harsh) I would advise to find someone on your department or center to help you in your research, because you need lot of help to perform analysis, or perhaps simply you don't trust your decisions. Either way I am sorry I cannot help more than I do. $\endgroup$ – llrs Jan 15 '19 at 13:54
  • 2
    $\begingroup$ @FereshTeh Then speak to who have hired you, that you need some courses or to go to meet someone or speak about your troubles. The sooner you solve it, the better for you. $\endgroup$ – llrs Jan 15 '19 at 14:10

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.