大家好~我是汇然，今天给大家带来一篇肿瘤免疫相关的5分+干湿结合生信文章。2020年3月发表在“Journal of Pathology”上，题名为“Molecular subtyping reveals immune alterations in IDH wild-type lower-grade diffuse glioma”。全文包含6个Figure和1个Table，接下来我们一起康康免疫相关的生信文章套路是如何应用的吧！
本篇范文研究的主要问题是弥漫性低度胶质瘤(LGG)。与胶质母细胞瘤(GBM, IV级)相比，其病程更为缓慢，且具有高侵袭性。异柠檬酸脱氢酶(Isocitrate dehydrogenase,IDH)野生型LGGs通常与预后不良有关。在WHO的新分类中，弥漫性LGGs根据IDH突变和1p/19q 共同缺失状态分为三个诊断和预后亚组。IDH野生型肿瘤以高级别(Ⅲ级)为主，预后最差，其次为IDH突变、无1p/19q共同缺失的肿瘤，其次为IDH突变、1p/19q共缺失的肿瘤(少突胶质细胞瘤)。但临床上，IDH野生型LGGs的结果与IDH野生型GBMs的结果难以区分，且比IDH突变型GBMs的预后更差。与此同时，还有研究证明了并非所有IDH野生型LGGs的存活率都很低。目前对IDH野生型LGGs的定义和分型仍有争议。因此，进行进一步分层分析，找到能够划分IDH野生型LGGs预后的标志物是非常有必要的。有很多研究表明在肿瘤微环境中，胶质瘤不同分子亚型间存在差异，其中固有基因表达亚型的进化与微环境的免疫变化相关。且有研究发现具有间充质基因表达特征的GBMs可能对免疫治疗更敏感。目前，神经胶质瘤的免疫治疗受到广泛关注，更好地了解肿瘤免疫微环境对于提高免疫治疗的疗效至关重要。作者首先在基因表达谱比较的基础上，利用R包“consensus susclusterplus”对从CGGA数据库获得的49个IDH野生型弥漫性LGGs队列中高度可变表达(MAD)> 1的2501个基因进行无监督聚类的亚型鉴定，并以热图的形式展示两个主要亚型(Fig.1A)。接下来，使用R包“princomp”进行主成分分析(PCA)验证这两种亚型的表达谱存在显著差异((Fig.1B)。并分析了这两种亚型聚类的总生存期(OS)和无进展生存期(PFS)，具有显著差异(Fig.1C)。结果显示Sub1患者的临床结果明显较差，Sub2的总体生存期更长。且Sub2中包含更多的II级肿瘤。接下来作者以从TCGA数据库中收集的73个IDH野生型弥漫性LGG表达谱作为验证集，以类似的基因排序再现了CGGA训练集中确定的亚组聚类、PCA以及生存分析验证(Fig.D-F)。作者对IDH野生型弥漫性LGG的两种亚型进行GO及GSEA富集分析以确定差异有统计学意义的基因集(Fig.2A-C)，结果显示，Sub1中高表达基因主要富集于免疫、炎症反应、抗原加工提呈等信号通路，并且该组患者免疫和炎症反应显著增强。为了进一步解释这两种亚型之间的免疫异质性，作者使用“estimate”R包计算基质细胞和免疫细胞的比例，推导每例患者的肿瘤纯度，从而检测各组间基质和免疫内容的分布。以箱线图的形式展示了基质、免疫、纯度、细胞溶解和炎症评分的比较，以及人类白细胞抗原(HLA)基因的表达(Fig.2D-G)。结果表明，Sub1含有更多的免疫细胞。由于天然抗肿瘤免疫需要溶细胞性免疫应答，于是作者通过量化颗粒酶和PRF1的平均表达来计算免疫细胞介导的溶细胞活性，以评估细胞毒性免疫细胞活性。与免疫评分一致，Sub1肿瘤的溶细胞评分明显高于Sub2。进一步，我们计算了肿瘤炎症信号算法，发现Sub1呈阳性。随后为探究微环境中免疫细胞的比例，采用CIBERSORT算法研究了各组间浸润性免疫细胞的组成。发现Sub1中M2巨噬细胞的富集程度较高，而Sub2肿瘤细胞中淋巴细胞、幼稚B细胞和浆细胞的富集程度较高(Fig.2H)。随后在TCGA队列中进行验证，发现具有一致的富集结果及免疫差异（Fig.3）。由于已知T细胞和NK细胞衰竭已被确定为癌细胞逃避宿主免疫的重要机制，作者随后检测了几种明星分子的表达(LAG3、CTLA4、PD1、PD-L1和HAVCR2)。结果发现大部分明星分子在训练集和验证集的Sub1中都有高表达(Fig.S5A,B)，这表明这些肿瘤的免疫衰竭水平升高。并进一步在12例来自CGGA队列的石蜡包埋组织中，针对一些经典标记物进行IHC验证(M2巨噬细胞的CD163；免疫溶细胞活性的GZMA；免疫衰竭的PD1和HAVCR2)。同样，这些基因在sub1中蛋白表达水平较高。总的来说，这些结果在IDH野生型弥漫性LGGs中发现了广泛的免疫异质性。联&靠——联合免疫细胞标记，分析与IDH野生型弥漫性LGGs预后相关的因素
作者接下来探讨临床意义，分别以溶细胞信号和炎症信号评分对训练集和验证集进行分层预后分析，结果表明炎症评分高与不良预后相关，而细胞溶解活性可能在此肿瘤中不起作用（Fig.4A-B）。另一方面，尽管具有较高的免疫和溶细胞活性，但Sub1中免疫衰竭基因的高表达意味着免疫抑制状态，这可能导致较差的预后。因此，作者接下来以相关性热图展示了炎症和细胞溶解特征与M1、M2巨噬细胞、CD8+ T细胞、中性粒细胞和记忆激活CD4+ T细胞的相关性，以及炎症评分、溶细胞评分和免疫抑制基因(PD-1和CTLA4)之间可能存在相关性(Fig.4C-D)。这些发现提示炎症微环境和浸润性溶细胞性细胞免疫应答可能影响IDH野生型弥漫性低级别瘤变患者的预后。第一步：TCGA队列中Sub1和Sub2之间基因组改变的差异分析作者随后分析了TCGA数据库中的体细胞拷贝数改变（SCNVs）数据，以探讨这两种亚型之间基因组改变的差异。突变特征提供了肿瘤发展的机制，结果发现两组之间在突变频率上没有明显的富集(Fig.5A)。但通过采用GISTIC 2.0分析两组间拷贝数变化则具有显著差异，Sub1显示了更多缺失位点，如CDKN2A/CDKN2B、DMRTA1、C9orf53和MTAP(Fig.5B)。通过上述分层聚类分析，表明以上新的分类方法可能影响疾病预后。于是，作者利用Cox比例风险模型识别预后相关的信号，筛选基于TCGA队列中Sub1和Sub2之间差异基因集。①采用SAM分析鉴定了1453个差异表达基因(FDR<0.05)；②对差异基因进行单变量Cox回归分析，筛选出与生存相关的326个基因（Fig.6A）；③采用Cox风险比例回归模型分析确定具有预后意义的基因，使用“glmnet”R包构建具有最佳预后价值的基因集(Fig.6B)。最终获得5个基因签名，采用基因表达加权回归系数(Coeffs)的线性组合来计算患者的风险得分(Fig.6C)。Ps：这里简单介绍一下“SAM”，是基因芯片著性分析方法的简称，一个由Standford大学开发的免费软件，用来筛选差异基因。其原理是将基因表达数据与临床特征数据关联，可以找到与该特征相关性强的基因，再结合其他方法或数据进一步分析，以揭示共内在的生物学意义。并且SAM软件在筛选得到较多差异基因的同时，FDR值还能够保持在较低的水平，因而是一种非常理想的差异基因筛选工具，有兴趣的小伙伴可以多多尝试哦~根据风险评分的中位数将肿瘤分为低风险组和高风险组。Kaplan-Meier分析显示，高危组患者预后明显差于低危组(P<0.001，Fig.6D)。并在CGGA队列中验证这些基因签名，得到一致的结果(Fig.6F-G)。并且对这些基因签名和病理特征进行联系，以风险评分将患者排序，发现在TCGA集合中按亚型、分级和EM/PM分级分层的病例的分别分析风险评分的分布情况（Fig.6E）。接下来，作者通过单因素和多因素Cox回归分析进一步评估了炎症信号的预后独立性。如预期，获得的基因签名与患者总体生存显著相关，且不受其他因素影响，表明上述分类在IDH野生型弥漫性LGGs中具有预后意义(Table 1)。随后，作者通过计算风险评分、年龄和分组的曲线下面积(AUC)来评估预测准确性。风险评分AUC明显高于其他因素(Fig.6I)。这些结果证明了该特征在预测预后方面的优越性能。本篇范文中，作者通过分层聚类的方法提出了一种新分类，用于实现基因型和表型的结合，从而把握不同分子特征患者临床转归的同质性，以及胶质瘤诊断的异质性。通过分子病理学诊断弥补传统组织病理学的不足，能够更加精准地反映肿瘤特点。作者采用从CGGA和TCGA获得的122例IDH野生型弥漫性LGGs的转录组数据，将其进行聚类分析，定义了两种具有不同临床和生物学特征的亚型。并将该分类在另一个独立的验证集中重复验证，观察到一致的表型，以证明这个分层具有临床应用价值。进一步的GO和GSEA分析发现，Sub2显著富集于化学突触传递、神经递质和谷氨酸等部分。因有研究表明神经回路中的突触和电整合促进胶质瘤进展，由此提示作者提出的分层在肿瘤进展中具有潜在应用价值。接下来进入重头戏——免疫预后分析。已知侵袭性胶质瘤的特征是基质、免疫评分和细胞溶解活性较高，但纯度较低，通过作者研究发现预后较差的Sub1病例在功能注释上表现出相似的免疫表型，且免疫浸润细胞分析富集于M2巨噬细胞（在胶质瘤中具有原核细胞活性）；而Sub2中表现出幼稚B细胞和浆细胞的高水平浸润（与胶质瘤患者良好预后相关）。同时发现无监督聚类的两个IDH野生型LGGs的稳健亚型，具有不同的预后和免疫特征。随后，作者根据免疫相关评分确定生存预后相关的分子标志物。并进行如Fig.6所示的预后分层分析，证明自己提出的这种分类在不同的治疗策略及预后中具有临床意义，并可能影响IDH野生型弥漫性LGG患者的靶向治疗效果。下面我给大家进行了免疫分子分型类生信文章套路总结：学习过后，我们同样可以将以上生信套路用于分析多种类型的数据，包括临床，SCNAs，基因突变，DNA甲基化和mRNA，miRNA，长非编码RNA和蛋白质的表达等，有待小伙伴们进一步探索发现哦~计算肿瘤纯度和免疫评分
接下来给大家介绍两款计算肿瘤纯度的R包，一款是利用表达数据计算基质与免疫评分，预测肿瘤纯度的必杀器ESTIMATE包；另一个是能够定量肿瘤细胞纯度、倍性、绝对拷贝数的ABSOLUTE包。下面我们就分别来康康这两个计算肿瘤纯度的R包有什么区别，又是如何运用的。
我们都知道鉴于肿瘤的异质性，单细胞转录组是揭示细胞异质性的的有力武器。在肿瘤样本的单细胞转录组不只是无监督聚类再分析那么简单，目前对肿瘤样本的基因表达谱数据库已经非常丰富了，我们可以从更多角度计算并分析肿瘤样本的异质性。ESTIMATE就是一种预测肿瘤纯度的工具，使用表达数据估计恶性肿瘤组织中的基质细胞和免疫细胞并且可以利用基因表达数据预测肿瘤组织中浸润的基质/免疫细胞的存在。该算法基于单个样本基因集合的富集分析，产生三个得分：下面我们来演示一下如何实战应用这个包。ESTIMATE其实是基于单样本基因集富集分析ssGSEA算法，因此这种算法与目前单细胞转录组中有些基于富集的细胞类型定义还是比较相似的。通过对基质细胞和免疫细胞的两个基因集在各肿瘤样本表达矩阵中打分，进而推断出细胞的类型。该功能可测量原始的ESTIMATE得分，得分与基于DNA拷贝数的肿瘤纯度呈正相关。另外，该算法使用的基因list是HUGO Gene Symbol或Entrez ID，我们再输入数据中的行名称必须是其中的一种。安装包
首先我们输入以下三行命令来安装ESTIMATE包，这一过程种如果大家提示有安装超时，可能是由于网络导致的。##安装“estimate”包library(utils)rforge <- 'http://r-forge.r-project.org'install.packages('estimate', repos=rforge, dependencies=TRUE)
示例数据的整理与读取
ESTIMATE包的作者使用的是Affymetrix U133Plus2.0平台的10个卵巢癌样本信息，该平台包含17,256个基因。随后将每个微阵列平台的不同基因数量统一为10,412个共同基因。我们通过下图所示的代码导入R包后读取数据。##运行“estimate”包>library(estimate)##以Affymetrix U133Plus2.0平台的10个卵巢癌样本信息为示例> OvarianCancerExpr <- system.file('extdata', 'sample_input.txt',+ package='estimate')> read.table(OvarianCancerExpr)[1:4,1:4]s516 s518 s519 s520C9orf152 4.881540 4.575656 3.739469 3.695996ELMO2 7.298054 7.555440 7.533202 7.382355CREB3L1 5.569164 5.700406 5.959730 5.770007RPS11 13.389937 13.848820 13.642862 13.654622> #最后，计算基质和免疫评分，分别代表肿瘤组织中基质和免疫细胞的存在。> filterCommonGenes(input.f=OvarianCancerExpr,+ output.f='OV_10412genes.gct',+ id='GeneSymbol')[1] 'Merged dataset includes 10412 genes (0 mismatched).'接下来我们需要先对文本文档中由10412个共同基因（行）和10个样本（列）组成的基因表达矩阵读入，并转为gct格式。完成数据格式转换后，使用estimateScore对gct格式的input表达矩阵进行计算，得到基质评分、免疫评分，以及估计分数，并保存到本地文件中。> estimateScore(input.ds = 'OV_10412genes.gct',+ output.ds='OV_estimate_score.gct',+ platform='affymetrix')[1] '1 gene set: StromalSignature overlap= 141'[1] '2 gene set: ImmuneSignature overlap= 141'
最后一个 plotPurity函数，即Plot tumor purity，用于绘制肿瘤纯度。该函数需要estimateScore方法输出的包含每个样品的基质，免疫，估计值和肿瘤纯度GCT文件。这一步可以将每个样本生成散点图，绘制肿瘤纯度与ESTIMATE得分的相关性图，并导出PNG文件。然而目前这个函数仅支持Affymetrix平台。> plotPurity(scores='OV_estimate_score.gct', samples='s516',+ platform='affymetrix')> scores=read.table('OV_estimate_score.gct',skip = 2,header = T)> rownames(scores)=scores[,1]> scores=t(scores[,3:ncol(scores)])> 可以看到很简单的代码，首先把txt文档里面的表达矩阵读入R里面转为gct格式，然后对gct格式的input表达矩阵使用estimateScore得到计算好的3个score值并且保存到本地文件。值如下：> scoresStromalScore ImmuneScore ESTIMATEScore TumorPuritys516 -281.81487 171.5411 -110.2737 0.8316075s518 -426.14692 105.3890 -320.7580 0.8483668s519 -57.14977 -365.2374 -422.3871 0.8561698s520 1938.82379 2339.0707 4277.8944 0.3314725s521 -671.64710 147.6183 -524.0288 0.8637832s522 1458.13837 1176.8159 2634.9543 0.5472110s523 -268.89216 -928.4953 -1197.3875 0.9092887s525 973.42289 1320.0869 2293.5098 0.5884565s526 552.64161 2162.4612 2715.1029 0.5373262s527 -709.33568 1312.8416 603.5059 0.7689656> ptm <- proc.time()最后重要的事情说三遍，引用！引用！记得引用以下三个参考文献。[1] Subramanian A, et al.'Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles.'Proc Natl Acad Sci U S A 2005, 102:15545-15550.
[2] Barbie DA, et al.'Systematic RNA interference reveals that oncogenic KRAS-driven cancers require TBK1.'Nature 2009, 462:108-112.
[3] Verhaak RG, et al.'Integrated genomic analysis identifies clinically relevant subtypes of glioblastoma characterized by abnormalities in PDGFRA, IDH1, EGFR, and NF1.'Cancer Cell 2010, 17:98-110.
可能大家都听说过ABSOLUTE软件，这是Broad研究所开发的一款癌症CNV分析软件，运用的就是R包“ABSOLUTE”，能够定量肿瘤细胞纯度与倍性，绝对拷贝数，如果提供突变数据，还能探究克隆与亚克隆等。我们都知道当从癌细胞和正常细胞的混合种群中提取DNA时，混合中会丢失每个癌细胞的绝对拷贝数信息。比如，实验中很常见的步骤，从混合样品中提取 DNA 进行测序后，得到的也是一个混合样品的结果。肿瘤因其本身异质性高，存在部分突变部分不突变的情况，不一定是单纯的二倍体，直接分析拷贝数变异，得到的结果并不准确，评估肿瘤倍性也更加必要。这时候就需要ABSOLUTE算法来进行肿瘤细胞纯度与倍性定量评估。ABSOLUTE的目的是从混合的DNA群体中重新提取这些数据。从生成分段拷贝数开始，我们此时将其与复发性癌症核型的预先计算模型，以及体细胞点突变的等位基因分数值（可以选择不输入）一起输入，随后ABSOLUTE会重新提取局部DNA片段的绝对细胞拷贝数以及针对点突变的突变等位基因数的信息，并进行输出和整合。那么是什么因素影响着绝对拷贝数和排序呢？这里涉及到样本异质性、核型模型以及等位基因分数三个方面。② 来自内置于ABSOLUTE算法参考面板的核型模型③ 来自突变数据的等位基因部分（可以选择不提供此项数据）一是要明确肿瘤纯度表示肿瘤细胞相对于正常细胞（二倍体）和污染样本的血细胞的比例二是由于正在进行的亚克隆进化，每个肿瘤群体按倍性分组，倍性以正常单倍体基因组的基因组片段定义。
为了更好地对这些解决方案进行排名，ABSOLUTE以核型模型的形式引用了外部数据。核型模型不会影响单个解决方案的计算，只会影响其排名。来自SCNA，SSNVS和全癌核型模型的可能性被合并以产生排名。基因组倍增的频率因肿瘤类型而异，反映了疾病组织的特异性生物学。ABSOLUTE用拷贝数数据和突变数据推断样品的纯度。或者通过估计细胞多样性，即每个癌细胞的平均等位基因拷贝数，以揭示潜在的亚克隆种群。考虑到癌症中不同类型基因组事件的可能存在很多差异，仅使用SCNA在推断肿瘤异质性方面条件有限。序列突变信息能够为ABSOLUTE提供在肿瘤进展中更多的增量信息，从而可以对肿瘤异质性进行更全面的建模。安装包这一步要较ESTIMATE的安装稍微复杂一点，我们首先进入这个网站，注册并登录，才能获得安装包。注册后跳转到如下界面。我们根据提示步骤下载“ABSOLUTE_1.0.6.tar.gz ”包，以及“HAPSEG_1.1.1.tar.gz”准备好主要的安装包后我们需要先安装两个依赖包'numDeriv'&'DNAcopy'，使用下列命令语句安装，如果提示错误或者安装不成功请使用R 4.0以下版本，或者直接从Bioconductor官网中下载安装包，进行本地安装即可。> install.packages('numDeriv')> if (!requireNamespace('BiocManager', quietly = TRUE))> install.packages('BiocManager')
然后我们将之前下载的两个包进行本地安装，或者直接使用以下命令语句：
> install.packages('HAPSEG_1.1.1.tar.gz', repos = NULL)> install.packages('ABSOLUTE_1.0.6.tar.gz',repos = NULL)如果癌症组织混合物中癌症细胞的比例是x，那么正常细胞的比例是1-x。我们将基因组中的每个基因座位点记为α，那么癌症细胞中基因座的整数拷贝数：q(α)，癌细胞的平均倍性：τ
那么，在混合组织中，位点α的绝对拷贝数是x q(α) +2(1 − x)；因此，位点α的拷贝数是：R(α) = (a q(α) + 2(1 − x))/D = (a /D) q(α) + (2(1 − x)/D)；体细胞的点突变计算式是：F(α) = (a sq(α))/Ds = (x /Ds )sq (α)。这个算法的步骤呢是先生成copy number的片段化信息（segment）文件，再利用癌细胞核型模型，输出结果中提供了DNA局部片段的细胞绝对拷贝数信息。这里如果我们提供了体细胞点突变的数据，结果会输出突变的等位基因频率，这些都是可选择性的。测试数据
我们使用官方（https://www.genepattern.org/modules/docs/ABSOLUTE/2）提供的测试数据，包含拷贝数变异的segments文件和体细胞突变的maf文件（该包不能批量处理，因此两份数据记录的是同一个样本），然后用 ABSOLUTE 的 RunAbsolute 函数进行分析。下图展示的是官网（http://software.broadinstitute.org/cancer/cga/absolute_run）提供的示例代码，这是使用来自先前HAPSEG运行的输入数据和捆绑数据对混合实验数据进行ABSOLUTE调用的示例。该文件应在bundle目录中运行。通过取消注释registerDoMC，调用并指定要使用的内核数，可以在多核系统上运行它。这里还需要使用foreach和可选的doMC。DoAbsolute <- function(scan, sif) {registerDoSEQ()library(ABSOLUTE)plate.name <- 'DRAWS'genome <- 'hg18'platform <- 'SNP_250K_STY'primary.disease <- sif[scan, 'PRIMARY_DISEASE']sample.name <- sif[scan, 'SAMPLE_NAME']sigma.p <- 0max.sigma.h <- 0.02min.ploidy <- 0.95max.ploidy <- 10max.as.seg.count <- 1500max.non.clonal <- 0max.neg.genome <- 0copy_num_type <- 'allelic'seg.dat.fn <- file.path('output', scan, 'hapseg',paste(plate.name, '_', scan, '_segdat.RData', sep=''))results.dir <- file.path('.', 'output', scan, 'absolute')print(paste('Starting scan', scan, 'at', results.dir))log.dir <- file.path('.', 'output', 'abs_logs')if (!file.exists(log.dir)) {dir.create(log.dir, recursive=TRUE)}if (!file.exists(results.dir)) {dir.create(results.dir, recursive=TRUE)}sink(file=file.path(log.dir, paste(scan, '.abs.out.txt', sep='')))RunAbsolute(seg.dat.fn, sigma.p, max.sigma.h, min.ploidy, max.ploidy,primary.disease,platform, sample.name, results.dir, max.as.seg.count,max.non.clonal,max.neg.genome, copy_num_type, verbose=TRUE)sink()}arrays.txt <- './paper_example/mix250K_arrays.txt'sif.txt <- './paper_example/mix_250K_SIF.txt'## read in array namesscans <- readLines(arrays.txt)[-1]sif <- read.delim(sif.txt, as.is=TRUE)library(foreach)## library(doMC)## registerDoMC(20)foreach (scan=scans, .combine=c) %dopar% {DoAbsolute(scan, sif)}obj.name <- 'DRAWS_summary'results.dir <- file.path('.', 'output', 'abs_summary')absolute.files <- file.path('.', 'output',scans, 'absolute',paste(scans, '.ABSOLUTE.RData', sep=''))library(ABSOLUTE)CreateReviewObject(obj.name, absolute.files, results.dir, 'allelic',verbose=TRUE)## At this point you'd perform your manual review and mark up the file## output/abs_summary/DRAWS_summary.PP-calls_tab.txt by prepending a column with## your desired solution calls. After that (or w/o doing that if you choose toaccept## the defaults, which is what running this code will do) run the followingcommand:calls.path = file.path('output', 'abs_summary', 'DRAWS_summary.PPcalls_tab.txt')modes.path = file.path('output', 'abs_summary', 'DRAWS_summary.PPmodes.data.RData')output.path = file.path('output', 'abs_extract')ExtractReviewedResults(calls.path, 'test', modes.path, output.path, 'absolute','allelic')
- 第三个图显示了前面3个解对应的拷贝轮廓，如纯度，倍性，异质性片段的比例。
- 第四个图是等位基因比例图，输入突变数据才会有。可以用于解释每个肿瘤细胞平均的等位基因拷贝数，以显示多样性，还可以揭示亚克隆突变。
[1] https://www.genepattern.org/modules/docs/ABSOLUTE/2
[2] http://software.broadinstitute.org/cancer/cga/absolute_run
[3] http://software.broadinstitute.org/cancer/software/genepattern/analyzing-absolute-data
[4] Carter SL, Cibulskis K, Helman E, McKenna A, Shen H, Zack T, Laird PW, Onofrio RC, Winckler W, Weir BA, Beroukhim R, Pellman D, Levine DA, Lander ES, Meyerson M, Getz G. Absolute quantification of somatic DNA alterations in human cancer. Nat Biotechnol. 2012;30(5):413-21.
这两款计算肿瘤纯度的R包，ESTIMATE通过使用表达数据估计恶性肿瘤组织中的基质细胞和免疫细胞，并且可以利用基因表达数据预测肿瘤组织中浸润的基质/免疫细胞的情况。ABSOLUTE可以直接从体细胞突变和拷贝数变异的结果中推断肿瘤纯度和倍性，如果提供突变数据，还可以探究克隆与亚克隆等。两者的原理类似，前者侧重于免疫细胞的分析，后者侧重于肿瘤细胞的倍性分析。我比较推荐ESTIMATE包，因为这个包不太需要很多代码经验也能够轻松上手，而且分析得到的结果比较适用于初步接触R语言和肿瘤方向生物信息学研究的，最重要的是又双叒叕靠上了“免疫”这个大热门，好发文章呀哈哈~小伙伴们快快尝试起来吧！我们之前讲过计算免疫微环境的两种主流方法为CIBERSORT与ssGSEA，那么今天我们来学习一下ssGSEA为何方神圣吧。single sample Gene Set Enrichment Analysis(ssGSEA)?虽然不认识他，但是他的兄弟GSEA总还是熟悉的吧。与他兄弟不同的是，它主要针对单样本无法做GSEA而提出的一种实现方法。那么他又是如何估计每个样本中的免疫细胞含量的呢？此处也用到了一种算法——随机漫步（此处略过，毕竟咱们也不是研究算法的嘛，这些工作就交给科研工作者就好）。这种算法精简版就是可以计算出每个基因集在这个样本中的富集分数。
大家之前应该看过CIBERSORT的代码操作了，那么可能各位小伙伴也会去选择CIBERSORT。CIBERSORT固然有他的优势：R包简便，也有网页在线操作优势显而易见。但是各位看官莫要忘记，“LM22”文件中只有22种免疫细胞，要是再想多一些那就只能sorry了。ssGSEA其中最明显的特征就是可以高度定制化。根据输入的gmt文件即可对样本数据进行随机游走，之后对每个样本计算出其免疫细胞的富集得分。话不多说，马上进入实操环节。##输入文件准备 首先我们需要准备输入文件。其中文件主要为3个： | |
| |
| |
| 需要运算的样本文件（因文件过大，可在TCGA官网自行下载） |
“symbol.txt”（仅展示部分）:数据为TCGA-LUSC的转录本FPKM文件“immune cell.txt”（仅展示部分）:免疫细胞基因集大家可以点开Rstudio,按照如下方式导入资料脚本” ssGSEA.R”接下来大家将光标放在第一行后，点击“Run”即可运行每行代码，此处不建议刷住多行代码直接点击Run,对于初学者应当感悟每行代码所代表的意义，而且万一出错会导致二次检查，费时费力rm(list=ls())setwd('E:\rxys\03') #设置工作目录#引用分析所用包。此处需要提前下载欧library(GSVA)library(limma)library(GSEABase)#读取输入文件，并对输入文件处理symbol=read.table('symbol.txt',sep=' ',header=T,check.names=F)#读取输入文件，此处为TCGA-LUSC的FPKM标准化数据symbol=as.matrix(symbol)rownames(symbol)=symbol[,1]exp_data=symbol[,2:ncol(symbol)]dimnames=list(rownames(exp_data),colnames(exp_data))exp_data1=matrix(as.numeric(as.matrix(exp_data)),nrow=nrow(exp_data),dimnames=dimnames)exp_data1=avereps(exp_data1)exp_data1=exp_data1[rowMeans(exp_data1)>0,]geneSet=getGmt('immune cell.gmt', geneIdType=SymbolIdentifier())#读取我们所提供的免疫细胞基因文件ssGSEA_Score=gsva(exp_data1, geneSet, method='ssgsea', kcdf='Gaussian', abs.ranking=TRUE)#ssGSEA计算normalize=function(x){return((x-min(x))/(max(x)-min(x)))}#定义ssGSEA_Score矫正函数norm_ssGSEA_Score=normalize(ssGSEA_Score)#对ssGSEA_Score进行矫正norm_ssGSEA_Score=rbind(id=colnames(norm_ssGSEA_Score),norm_ssGSEA_Score)write.table(norm_ssGSEA_Score,file='norm_ssGSEA_Score.txt',sep=' ',quote=F,col.names=F)#此处的输出文件即为ssGSEA富集得分文件
此时会在文件夹中出现结果文件“norm_ssGSEA_Score.txt”，我们将其打开看一下内容。此时在文件中会对应展示每个样品中免疫细胞的富集分数，为了更好的展示，我们使用excel打开该文件（如下图：仅显示部分截图） 我们在论文中往往会对数据进行可视化操作，接下来就让我们对数据进行一下绘图操作。在此之前要对于结果文件进行处理，因为我们在ssGSEA中大多数情况下只对于肿瘤样本进行计算，因此我们要先筛掉那些正常样本screen=sapply(strsplit(colnames(norm_ssGSEA_Score),'\-'),'[',4)screen=sapply(strsplit(screen,''),'[',1)screen=gsub('2','1',screen)data=norm_ssGSEA_Score[,screen==0]#对肿瘤样品分组，TCGA的命名规则可以复习一波h_cluster = hclust(dist(t(norm_ssGSEA_Score)))#计算样品距离并聚类cluster_result=cutree(h_cluster,3)#将聚类结果分成低，中，高免疫三组library(pheatmap)norm_ssGSEA_Score1=read.table('norm_ssGSEA_Score.txt',sep=' ',header=T,row.names=1,check.names=F)#读取结果文件V1<-names(cluster_result)V2<-as.integer(cluster_result)category<-data.frame(V1,V2)category$V1<-as.character(category$V1)category[,2]=paste0('Cluster',category[,2])category=category[order(category[,2]),]norm_ssGSEA_Score1=norm_ssGSEA_Score1[,as.vector(category[,1])]cluster=as.data.frame(category[,2])row.names(cluster)=category[,1]colnames(cluster)='Cluster'pdf('heatmap.pdf',height=5,width=9)color.key<-c('#3300CC','#3399FF','white','#FF3333','#CC0000')pheatmap(norm_ssGSEA_Score1, annotation=cluster, color =colorRampPalette(color.key)(50), cluster_cols =F, fontsize=8, fontsize_row=8, scale='row', show_colnames=T, fontsize_col=3, main = 'Example heatmap')dev.off()
从图中可以看到低、中、高三组有较为明显的区分，接下来我们就可以利用这些分组完成一些自己课题的后续研究了。想必到了这里，大家应该对ssGSEA有了一个更为深刻的了解，趁热打铁，我们直接来到应用部分。我们今天介绍的是2019年发表在Oncology reports上的一篇文章“Tumor‑infiltrating M2 macrophages driven by specific genomic alterations are associated with prognosis in bladder cancer”这篇文章采用了ssGSEA与CIBERSORT两种方法进行免疫浸润的相关研究。最终通过膀胱癌的组织学分级、病理分期与M2巨噬细胞浸润水平的关联性，得出膀胱癌特异性基因组改变可能是M2巨噬细胞浸润的重要驱动力 在这里我们放两张论文中的图片：是不是与我们的可视化结果较为相似，不过这个是利用CIBERSORT的结果做出来的。不过我们可以看一下下面这幅图在这幅图中作者采用ssGSEA的方法对样品中的富集分数的计算（也就是我们的结果文件啦）再配合临床分期进行可视化操作。从而得出与低度膀胱癌组织相比，B细胞，巨噬细胞，嗜中性粒细胞等细胞在高度膀胱癌中富集更为显著。文献的全文我们已经打包在资料中，感兴趣的小伙伴可以复现一下这篇文献。无论是在线网站分析或者R语言代码分析，我们都不要忘记引用算法文献呀！！！[1] GSVA: The Gene Set Variation Analysis package for microarray and RNA-seq data
[2] Xue Y, Tong L, LiuAnwei Liu F, et al. Tumor‑infiltrating M2 macrophages driven by specific genomic alterations are associated with prognosis in bladder cancer. Oncol Rep. 2019;42(2):581-594. doi:10.3892/or.2019.7196
[3] Jia Q, Wu W, Wang Y, et al. Local mutational diversity drives intratumoral immune heterogeneity in non-small cell lung cancer. Nat Commun. 2018;9(1):5361. Published 2018 Dec 18. doi:10.1038/s41467-018-07767-w
最后给小伙伴们留一个小作业：大家看到的这幅图可以依据我们的结果文件并利用pheatmap包完成其中热图的部分绘制，大家可以动手尝试一下了，看一看谁做的比较像吧:
