本文是参考学习单细胞转录组高级分析六:TCGA生存分析的学习笔记。可能根据学习情况有所改动。
前言
通过scRNA分析得到新的marker基因,gene signatures,以及解释研究结论的基因集和代谢通路,一般都会用公共数据库的bulkRNA数据验证一下。我们这个专题所用数据的来源文章就利用TCGA数据对TFH signatures做了生存分析。
We next utilized clinical data and bulk mRNA-seq data from The Cancer Genome Atlas (TCGA) (Liu et al., 2018) to determine whether CD4+ TFH signatures in HNSCC were related to survival (STAR Methods).
本篇文章将模拟原文的过程做一个生存分析,但是我在原文中找不到TFH signatures对应的200个基因,所以只能用其中的几个基因模拟一下。本次分析的结果与原文会有一定差异,大家关注生存分析如何操作即可。
UCSC Xena下载数据
UCSC Xena简介
UCSC是加利福尼亚大学圣克鲁兹分校(University of California, Santa Cruz)的简称,UCSC Xena是他们开发的一个癌症基因组学数据分析平台。该平台整合了多个癌症公共数据库的资源,比如来自TCGA, ICGC等大型癌症研究项目的数据,不仅可以下载各个数据集,还提供在线分析功能。
HNSC数据集下载示例打开https://xenabrowser.net/datapages/,进入数据集选择页面:
很容易就能找到我们需要的Head and Neck Cancer
大家应该会发现这个页面会有GDC TCGA和TCGA两个开头的数据,TCGA是原版的数据,GDC TCGA是GDC(美国癌症研究所开发维护的一个癌症数据库资源门户网站,Genomic Data Commons)整理的TCGA数据,两个数据大同小异。点击箭头所指链接,进入以下页面:
下载箭头所指的表达数据和生存数据
表达数据下载页面,
箭头所指表达矩阵和Gene Mapping文件都要下载,
生存下载页面
最后得到三个文件:
生存分析
读取TCGA数据
##安装生存分析R包
GSVA分析
原文富集分析用的200个基因找不到具体名称,我用了作者特异强调的8个上调基因代替,如下所示:
we first defined our TFH gene set based on the top 200 differentially upregulated genes between CD4+ Tconv cluster 1 and CD4+ Tconv cluster 7, which represent the two terminal-most differentiated states of CD4+ T cells. As a test statistic for enrichment in the bulk expression profile from each patient, we used the Kolmogorov-Smirnov (KS) test to compare genes in the gene set versus those not in the gene set.
TFHset = list(c('PDCD1', 'CXCL13', 'TIGIT', 'TOX', 'MAF', 'CXCR5', 'CD40LG', 'IL6ST'))
生存曲线
ES <- data.frame(t(ES), check.names=F)
生存曲线大体趋势与原文相符,差异没有他们的图那么明显。原文用于富集分析的基因数量远多于我演示用的基因数目,并且他们还对TCGA的数据做了过滤,这可能是他们效果更好的主要原因。Cox回归分析上述Kaplan-Meier方法只能针对分类变量做单因素生存分析,因此我将GSVA分析的结果转换成了只有high和low的二分类变量。如果不想转换可以用Cox比例风险回归模型分析:
res.cox <- coxph(Surv(OS.time, OS)~TFHset, data=TCGAclin)
Cox结果解读:
z值:Wald统计值,它对应于每个回归系数与其标准误差的比率(z = coef / se(coef))。
coef:回归系数,正数表示危险(死亡风险)较高,负数代表危险较低。
exp(coef):在此代表危险比(HR,hazard ratio)。