单细胞转录组高级分析六：TCGA生存分析

单细胞天地 2022-06-07

The following article is from 生信会客厅 Author Kinesin

上期专题我们介绍了单细胞转录组数据的基础分析，然而那些分析只是揭开了组织异质性的面纱，还有更多的生命奥秘隐藏在数据中等待我们发掘。本专题将介绍一些单细胞转录组的高级分析内容：多样本批次校正、转录因子分析、细胞通讯分析、基因集变异分析和更全面的基因集富集分析。不足之处请大家批评指正，欢迎添加Kinesin微信交流探讨!

前言

通过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包BiocManager::install(c("survival","survminer"), update = F)##读取HNSC数据集的表达矩阵和临床信息library(survival)library(survminer)library(tidyverse)dir.create("Survival")TCGAexpr <- read.delim("Data/TCGA-HNSC.htseq_fpkm.tsv", check.names=F)id.annota <- read.delim("Data/gencode.v22.annotation.gene.probeMap", check.names=F)id.annota <- id.annota[, 1:2]TCGAexpr <- merge(id.annota, TCGAexpr, by.x="id", by.y="Ensembl_ID")TCGAexpr <- aggregate(TCGAexpr[,-c(1,2)], list(TCGAexpr$gene), FUN=sum)TCGAexpr <- column_to_rownames(TCGAexpr,"Group.1")#过滤表达量过低的基因，减少分析数据量TCGAexpr <- TCGAexpr[rowSums(TCGAexpr)>10, ] %>% as.matrix()TCGAclin <- read.delim("Data/TCGA-HNSC.survival.tsv", check.names = F)TCGAclin <- TCGAclin[TCGAclin$sample %in% names(TCGAexpr),]TCGAclin <- TCGAclin[,-3]

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'))names(TFHset) <- 'TFHset'ES <- gsva(TCGAexpr, TFHset)

生存曲线

ES <- data.frame(t(ES), check.names=F)TCGAclin <- merge(TCGAclin, ES, by.x="sample", by.y=0)TCGAclin$TFH.set = ifelse(TCGAclin$TFHset>median(TCGAclin$TFHset),'high','low')png('Survival/survival.png', 800, 600)ggsurvplot(survfit(Surv(OS.time, OS)~TFH.set, data=TCGAclin), conf.int=F, pval=TRUE)dev.off()write.csv(TCGAclin, 'Survival/TCGAclin.csv', row.names = F)

生存曲线大体趋势与原文相符，差异没有他们的图那么明显。原文用于富集分析的基因数量远多于我演示用的基因数目，并且他们还对TCGA的数据做了过滤，这可能是他们效果更好的主要原因。

Cox回归分析

上述Kaplan-Meier方法只能针对分类变量做单因素生存分析，因此我将GSVA分析的结果转换成了只有high和low的二分类变量。如果不想转换可以用Cox比例风险回归模型分析：

res.cox <- coxph(Surv(OS.time, OS)~TFHset, data=TCGAclin)res.cox# Call:# coxph(formula = Surv(OS.time, OS) ~ TFHset, data = TCGAclin)

# coef exp(coef) se(coef) z p# TFHset -0.3702 0.6906 0.1371 -2.7 0.00694

# Likelihood ratio test=7.33 on 1 df, p=0.006786# n= 545, number of events= 252

Cox结果解读：