estimate或者CIBERSORT结果真的是很好的临床预后指标吗

生信技能树

发布于 2021-10-12 12:10:34

1.6K0

发布于 2021-10-12 12:10:34

文章被收录于专栏：生信技能树生信技能树

肿瘤免疫微环境我们讲了很多内容了，目录是：

都是依据肿瘤病人的转录组测序表达量矩阵进行的分析，也有几百篇类似的数据挖掘文章了，它们总是喜欢落脚到estimate或者CIBERSORT结果的预后意义。

那么，我们就来实际检验看看estimate或者CIBERSORT结果真的是很好的临床预后指标吗！

首先看看 estimate 的StromalSignature 和 ImmuneSignature

载入前面步骤的 estimate 算法得到的StromalSignature 和 ImmuneSignature，如果你不知道下面的rdata文件（ssgseaScore-for-ssGSEA_for_seurat.Rdata ）来源，建议再看看不同癌症内部按照estimate的两个打分值高低分组看蛋白编码基因表达量差异 ：

load('../expression/estimate_results/ssgseaScore-for-ssGSEA_for_seurat.Rdata')
head(ssgseaScore)

> head(ssgseaScore)
                 StromalSignature ImmuneSignature
TCGA-OR-A5JP-01A      -0.06619425     -0.13854919
TCGA-OR-A5JG-01A      -0.09059249     -0.07372133
TCGA-OR-A5K1-01A      -0.06014336     -0.02297215
TCGA-OR-A5JR-01A      -0.09859187     -0.01468964
TCGA-OR-A5KU-01A      -0.11097034     -0.16760435
TCGA-OR-A5L9-01A       0.01055877      0.10013901

ssgseaScore$pid=substring(rownames(ssgseaScore),1,12)
head(ssgseaScore)

接下来批量进行生存分析，区分各个癌症

仍然是32种癌症，内部的estimate 算法得到的StromalSignature 和 ImmuneSignature，各自根据中位值进行分组，然后进行生存分析看看estimate 算法得到的StromalSignature 和 ImmuneSignature指标是否有生存意义：

library(survminer) 
library(survival)
load(file = 'phe_stage.Rdata')
phe[1:4,1:4]

tp=unique(phe$type)
tp

sur_list <- lapply(tp, function(x){
  # x=tp[1]
  this_phe=phe[phe$type==x,]
  # 这里先看 os 
  survival_dat=as.data.frame(this_phe[,c('bcr_patient_barcode','OS','OS.time')])
  head(survival_dat)
  
  colnames(survival_dat)=c('pid','event','time')
  survival_dat=merge(survival_dat,ssgseaScore,by='pid') 
  survival_dat$time =  survival_dat$time/365
  
  survival_dat$group=ifelse(survival_dat$StromalSignature>median(survival_dat$StromalSignature),
                            'stromal_high','stromal_low')
  fit <- survfit(Surv(time, event) ~ group,
                 data = survival_dat)
  stromal_survp=ggsurvplot(fit,data = survival_dat, #这里很关键，不然会报错
                   legend.title = x,
                   pval = T, #在图上添加log rank检验的p值 
                   risk.table = F, #在图下方添加风险表
                   xlab = "Time in years", #x轴标题
                   xlim = c(0, 10), #展示x轴的范围
                   break.time.by = 1, #x轴间隔
                   size = 1.5#线条大小
  )
  survival_dat$group=ifelse(survival_dat$ImmuneSignature>median(survival_dat$ImmuneSignature),
                            'immune_high','immune_low')
  fit <- survfit(Surv(time, event) ~ group,
                 data = survival_dat)
  immune_survp=ggsurvplot(fit,data = survival_dat, #这里很关键，不然会报错
                           legend.title = x,
                           pval = T, #在图上添加log rank检验的p值 
                           risk.table = F, #在图下方添加风险表
                           xlab = "Time in years", #x轴标题
                           xlim = c(0, 10), #展示x轴的范围
                           break.time.by = 1, #x轴间隔
                           size = 1.5#线条大小
  )
  
  return(list(
    immune_survp=immune_survp,
    stromal_survp=stromal_survp
  ))
})

这里每个癌症都返回来了一个list，包含这stromal和immune的高低分组后的差异分析图！

对ImmuneSignature 出图：


x=6;y=6
all_plot <- arrange_ggsurvplots(lapply(sur_list, function(x) x[[1]]),
                                print = F,
                                ncol =x, nrow = y) 

x=30;y=30
ggsave(all_plot,filename = 'stromal_sur_plot.pdf',
       width = x,height = y)

可以看到 ImmuneSignature 中位数分组后在部分癌症是有统计学显著的生存意义哦！

ImmuneSignature 中位数分组生存分析

对 StromalSignature 出图：

x=6;y=6
all_plot <- arrange_ggsurvplots(lapply(sur_list, function(x) x[[2]]),
                                print = F,
                                ncol =x, nrow = y) 

x=30;y=30
ggsave(all_plot,filename = 'immune_sur_plot.pdf',
       width = x,height = y)

可以看到 StromalSignature 分组后在部分癌症是有统计学显著的生存意义哦！

StromalSignature 分组生存分析

这个时候其实我们无需看这么多图，仅仅是输出p值，还有HR值即可！

survdiff代替 survfit

代码如下所示：

# surv 构建对象
# survfit 拟合生存曲线
# survdiff 差异检验

options(digits = 2)
km_results <- do.call(rbind,
        lapply(tp, function(x){
          # x=tp[1]
          this_phe=phe[phe$type==x,]
          # 这里先看 os 
          survival_dat=as.data.frame(this_phe[,c('bcr_patient_barcode','OS','OS.time')])
          head(survival_dat)
          
          colnames(survival_dat)=c('pid','event','time')
          survival_dat=merge(survival_dat,ssgseaScore,by='pid') 
          survival_dat$time =  survival_dat$time/365
          
          survival_dat$group=ifelse(survival_dat$StromalSignature>median(survival_dat$StromalSignature),
                                    'stromal_high','stromal_low')
          fit <- survdiff(Surv(time, event) ~ group,
                          data = survival_dat)
          stromal_p.val = 1 - pchisq(fit$chisq, length(fit$n) - 1)
          stromal_HR = (fit$obs[2]/fit$exp[2])/(fit$obs[1]/fit$exp[1])
          
          survival_dat$group=ifelse(survival_dat$ImmuneSignature>median(survival_dat$ImmuneSignature),
                                    'immune_high','immune_low')
          fit <- survdiff(Surv(time, event) ~ group,
                          data = survival_dat)
          immune_p.val = 1 - pchisq(fit$chisq, length(fit$n) - 1)
          immune_HR = (fit$obs[2]/fit$exp[2])/(fit$obs[1]/fit$exp[1])
          
          return(c(stromal_p.val,stromal_HR,
                   immune_p.val,immune_HR))
        }))
rownames(km_results)=tp
km_results

结果如下所示，可以看到BRCA这个癌症里面的 estimate 算法得到的StromalSignature 和 ImmuneSignature都是可以区分生存，因为p值都是0.05附近，结合前面的图表，可以看到：

其中StromalSignature 高，死得快。
而ImmuneSignature高死的慢，是保护因子。

而且可以看到下面的HR值也有可能是反过来了的，需要自行甄别！

而且真的是每个癌症都不一样，两个打分不一定是生存分析统计学显著，不一定是保护因子或者风险因子，唉！

# stromal_p.val,stromal_HR, immune_p.val,immune_HR 
> round(km_results,2)
     [,1] [,2] [,3] [,4]
ACC  0.70 1.16 0.03 2.34
BLCA 0.03 0.73 0.85 1.03
BRCA 0.06 0.77 0.04 1.33
CESC 0.47 1.18 0.06 1.56
CHOL 0.15 1.86 0.04 2.55
COAD 0.41 0.86 0.81 0.96
DLBC 0.87 1.11 0.63 0.74
ESCA 0.66 0.90 0.68 0.91
GBM  0.20 0.80 0.21 0.81
HNSC 0.23 1.16 0.04 1.30
KICH 0.76 0.84 0.64 1.31
KIRC 0.78 1.04 0.13 0.81
KIRP 0.06 0.59 0.83 1.06
LAML 0.80 0.95 0.01 0.59
LGG  0.00 0.57 0.01 0.63
LIHC 0.65 0.93 0.75 0.95
LUAD 0.16 1.21 0.07 1.28
LUSC 0.03 0.76 0.18 0.84
MESO 0.01 0.55 0.45 1.19
OV   0.63 0.94 0.36 1.13
PAAD 0.54 0.88 0.80 1.05
PCPG 0.17 2.92 0.64 1.39
PRAD 0.22 2.21 0.59 1.41
READ 0.16 0.61 0.80 1.10
SARC 0.06 1.47 0.01 1.69
SKCM 0.04 1.32 0.00 1.75
STAD 0.10 0.77 0.33 0.85
TGCT 0.04 0.00 0.19 0.25
THCA 1.00 1.00 0.75 1.16
THYM 0.92 0.94 0.73 1.26
UCEC 0.48 1.16 0.01 1.65
UCS  0.68 0.87 0.73 1.12
UVM  0.01 0.35 0.00 0.23

然后看CIBERSORT的22种免疫细胞比例的生存意义

同样的，载入前面步骤得到的CIBERSORT的22种免疫细胞比例，也是可以自己回顾使用CIBERSORT算法推断全部tcga样品的免疫细胞比例教程：

load('../expression/CIBERSORT_results/cibersort_raw-for-seurat.Rdata')
colnames(cibersort_raw)
[1] "Patients"                    
 [2] "B.cells.naive"               
 [3] "B.cells.memory"              
 [4] "Plasma.cells"                
 [5] "T.cells.CD8"                 
 [6] "T.cells.CD4.naive"           
 [7] "T.cells.CD4.memory.resting"  
 [8] "T.cells.CD4.memory.activated"
 [9] "T.cells.follicular.helper"   
[10] "T.cells.regulatory..Tregs."  
[11] "T.cells.gamma.delta"         
[12] "NK.cells.resting"            
[13] "NK.cells.activated"          
[14] "Monocytes"                   
[15] "Macrophages.M0"              
[16] "Macrophages.M1"              
[17] "Macrophages.M2"              
[18] "Dendritic.cells.resting"     
[19] "Dendritic.cells.activated"   
[20] "Mast.cells.resting"          
[21] "Mast.cells.activated"        
[22] "Eosinophils"                 
[23] "Neutrophils"                 
[24] "group"                       
[25] "type"

就作为学徒作业吧，我就不贴代码也不给结果啦！

每个细胞类型在每个癌症的生存情况，都需要km检验，都有一个p值和一个hr值，计算起来并不难，但是可视化就有点麻烦了！后面我们再细说这个可视化！

本文参与腾讯云自媒体分享计划，分享自微信公众号。

原始发表：2021-10-01，如有侵权请联系 cloudcommunity@tencent.com 删除

编程算法

本文分享自生信技能树微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体分享计划，欢迎热爱写作的你一起参与！

编程算法

登录后参与评论

0 条评论

热度