泛癌水平的批量生存分析

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

都是依据肿瘤病人的转录组测序表达量矩阵进行的分析，也有几百篇类似的数据挖掘文章了，它们总是喜欢落脚到estimate或者CIBERSORT结果的预后意义。但是实际上我们也代码演示了：estimate或者CIBERSORT结果真的是很好的临床预后指标吗，这样做风险很大，后面留了一个思考题，就是CIBERSORT的22种免疫细胞比例的生存意义的全部癌症的探索，呼应我们的主题《泛癌水平的批量生存分析》。

但是呢，CIBERSORT的22种免疫细胞比例毕竟是算法推断，而且是从bulk转录组表达量数据里面推断出来的。目前我们已经有了单细胞技术，其实很容易收集整理几十个癌症的数据集，去看真实的各个细胞亚群的比例 :

 [1] "B.cells.naive"               
 [2] "B.cells.memory"              
 [3] "Plasma.cells"                
 [4] "T.cells.CD8"                 
 [5] "T.cells.CD4.naive"           
 [6] "T.cells.CD4.memory.resting"  
 [7] "T.cells.CD4.memory.activated"
 [8] "T.cells.follicular.helper"   
 [9] "T.cells.regulatory..Tregs."  
[10] "T.cells.gamma.delta"         
[11] "NK.cells.resting"            
[12] "NK.cells.activated"          
[13] "Monocytes"                   
[14] "Macrophages.M0"              
[15] "Macrophages.M1"              
[16] "Macrophages.M2"              
[17] "Dendritic.cells.resting"     
[18] "Dendritic.cells.activated"   
[19] "Mast.cells.resting"          
[20] "Mast.cells.activated"        
[21] "Eosinophils"                 
[22] "Neutrophils"       

但是单细胞研究有一个问题， 样品数量太少了，短时间内不太可能达到TCGA计划的每个癌症成百上千的样品数量。所以通常很少看比例，而是看各个细胞亚群的标记基因组成的特征基因集是否有特殊的生物学意义，比如生存分析的统计学显著。其中，2020的文章：《Identification of EMT signaling cross-talk and gene regulatory networks by single-cell RNA sequencing》，就做了一个很好的展示，它本来是单细胞研究，前面也是很标准的降维聚类分群：

标准的降维聚类分群

这个时候，每个细胞亚群，都会有自己的特异性高表达量基因，组成为了基因集。然后研究者拿这些基因集去TCGA数据库里面检验它们是否在各个癌症里面可以统计学显著的区分生存，而且判定它们是保护因子还是风险因子。如下所示，展示的是这15个单细胞亚群的标记基因在TCGA计划的BRCA癌症里面的HR值，关于HR值：

  • HR = 1: No effect
  • HR < 1: Reduction in the hazard
  • HR > 1: Increase in Hazard

Note that in cancer studies:

  • A covariate with hazard ratio > 1 (i.e.: b > 0) is called bad prognostic factor
  • A covariate with hazard ratio < 1 (i.e.: b < 0) is called good prognostic factor
判定它们是保护因子还是风险因子

如果是在多个癌症，上面的条形图就不适合了，所以热图展现 ：

热图看泛癌水平的批量生存分析

可以看到， 标红的C0和C10以及C13这3个单细胞亚群的特异性高表达量基因集，在多个癌症，都是 风险因子。接下来就可以挑选C0和C10以及C13这3个单细胞亚群，来展现它们特异性高表达量基因集，根据log2FC排序的：

 

其实C10和C13其实主要是细胞周期基因比如TOP2A和CDK1等等，它必然是肿瘤风险因子啦，在多个癌症里面都是表达量越高癌症病人预后越差。

文章的Pan-Cancer Survival Analysis的具体步骤是：

  • The Kaplan–Meier method was used to determine survival probability. The P values were determined by a log-rank test.
  • The signatures for each cluster were used to retrieve signature enrichment scores. The NES then categorized as high or low based on mean.
  • Univariate Cox proportional hazards models were fitted to calculate the hazard ratios using the coxph function in Sur- vival (v 2.44). P values less than 0.05 were considered to be statistically significant.

也就是说作者使用的是gsea方法来判定各个基因集在各个癌症的 enrichment scores，然后cox分析就是依据 enrichment scores啦，不过呢，作者并没有把 enrichment scores当做是连续值来看，而是依据中位值进行高低分组后再进行cox分析。

那，我们使用TCGA数据来演示：

###　首先查看 pan-cancer官网自带一个免疫细胞比例

这个 TCGA.Kallisto.fullIDs.cibersort.relative.tsv 文件，在前面的教程里面有给出下载地址：

# pan-cancer官网自带一个免疫细胞比例
library(data.table)
cib = fread('TCGA.Kallisto.fullIDs.cibersort.relative.tsv',data.table = F)
cib[1:4,1:4]
codes=substring(cib$SampleID,14,16)
table(codes)
cib = cib[codes=='01A',]
cib$bcr_patient_barcode=gsub('[.]','-',substring(cib$SampleID,1,12))

可以看到这个 TCGA.Kallisto.fullIDs.cibersort.relative.tsv 文件 主要是记录了TCGA数据的全部样品的各自的22个免疫细胞打分情况  ：

> cib[1:4,1:4]
                      SampleID CancerType B.cells.naive B.cells.memory
1 TCGA.OR.A5JG.01A.11R.A29S.07        ACC  0.000000e+00     0.04852871
2 TCGA.OR.A5LG.01A.11R.A29S.07        ACC  7.168876e-03     0.01112538
3 TCGA.OR.A5JD.01A.11R.A29S.07        ACC  2.342245e-05     0.01460659
4 TCGA.OR.A5LH.01A.11R.A29S.07        ACC  4.729908e-02     0.03817952
> dim(cib)
[1] 9780   28

然后整理临床信息：

同样的 TCGA-CDR-SupplementalTableS1.xlsx 文件，在前面的教程给出来了下载地址：

library(readxl)
phe=read_excel('../MC3/TCGA-CDR-SupplementalTableS1.xlsx',sheet = 1)
colnames(phe)
phe[1:4,1:4]
phe=phe[,c(2,3,26:33)]
head(phe)
library(gplots)
balloonplot(table(phe$gender,phe$type))

colnames(cib)
colnames(phe)
dat=merge(as.data.frame(phe),cib,by='bcr_patient_barcode')
colnames(dat)

dat$OS.time =  dat$OS.time/365
boxplot(dat$OS.time)
dat=dat[dat$OS.time>0,]
table(dat$type)

接下来写一个函数

这个函数可以在前面的临床信息加上免疫细胞比例信息，里面提取容易免疫细胞再区分癌症进行批量生存分析，得到HR值和对应的p值：


library(survminer) 
library(survival)
sur_for_df <- function(df,choose ,by){
  # 挑选指定的癌症，做循环 
  tp=sort(unique(df[,choose]))
  tp
  sur_list <- lapply(tp, function(x){
    # x=tp[1]
    print(x)
    this_phe=df[ df[,choose] == x,]
    # 这里先看 os ，挑选指定的细胞比例来进行高低分组 
    survival_dat=this_phe[,c( by ,'OS','OS.time')]
    colnames(survival_dat)=c('group','event','time')
    survival_dat=na.omit(survival_dat)
    survival_dat$group = ifelse( survival_dat$group > median( survival_dat$group),
                                 'high','low')
    # 如果无法分组，就不用进行生存分析比较啦 
    if(length(unique(survival_dat$group))==2){
      fit <- survfit(Surv(time, event) ~ group,
                     data = survival_dat)
      # 这里无需出图，所以下面的代码注释掉即可,但是我不想删除它们，以后可以使用
     if(F){
       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(survp)
     }
      # 这个survdiff 函数，紧跟着前面的 survfit 函数， 来对两个分组进行统计检验
      fit <- survdiff(Surv(time, event) ~ group,
                      data = survival_dat)
      p.val = 1 - pchisq(fit$chisq, length(fit$n) - 1)
      HR = (fit$obs[1]/fit$exp[1])/(fit$obs[2]/fit$exp[2])
      return(c(p.val,HR))
    }

})
  return(sur_list)
  
}

所以我们可以对全部的免疫细胞比例独立运行上面的函数，并且保存信息：

# choose='type'
# by='Monocytes'
# df=dat
unlist(lapply(sur_for_df(dat,'type','Monocytes'), '[',1))
hr_df = do.call(rbind,
        lapply(colnames(dat)[13:34], function(celltype){
          unlist(lapply(sur_for_df(dat,'type', celltype ), '[',1))
        }))
dim(hr_df)
rownames(hr_df)=colnames(dat)[13:34]
colnames(hr_df)=sort(unique( dat$type))

pheatmap::pheatmap(hr_df)

希望大家可以看出来一些规律：

不同免疫细胞比例的生存情况

挺奇怪的，基本上没有显而易见的规律啊！

再次呼应了前面的结果：estimate或者CIBERSORT结果真的是很好的临床预后指标吗 ？

写在文末

欢迎加入我们pan-cancer数据挖掘的讨论群，因为已经满200人了，所以需要我们生信技能树的官方拉群小助手帮忙拉群哦！！！（名额有限，先到先得！！！）

这个时候请直接付款28元给小助手，就可以进群，或者你转发此推文到朋友圈然后截图给小助手，就可以仍然以18.8元进群！

一个简单的门槛，隔绝那些营销号！我们也会在群里共享生存分析相关的资料，仅此而已，考虑清楚哦！

赞 (0)

相关推荐