单细胞转录组高级分析六：TCGA生存分析
本文是参考学习单细胞转录组高级分析六：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文件都要下载，
image.png
生存下载页面
最后得到三个文件：
图片
生存分析
读取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)。
