Multi-omics Data Analyses Construct TME and Identify the Immune-Related Prognosis Signatures in Human LUAD
多组学数据分析构建TME并识别人类LUAD中与免疫相关的预后特征
发表期刊：Mol Ther Nucleic Acids
发表日期：2020 Sep 4
影响因子：7.032
DOI: 10.1016/j.omtn.2020.07.024
一、研究背景
肺癌是所有癌症类型中发病率和死亡率最高的癌症。许多流行病学调查和实验研究认为，LUAD的发生和发展主要与环境因素和基因改变有关。到目前为止，与基因相关的治疗策略主要有两类，即靶向治疗和免疫治疗。
TME通常被定义为肿瘤周围的环境，其中包括细胞外基质、血管以及免疫细胞和神经元等细胞角色，所有这些都与肿瘤的进展和治疗结果有很大的关系。越来越多的研究通过实验阐明了TME浸润在各种癌症类型的免疫治疗反应和耐药中的贡献作用，并探讨了它们对患者预后的影响。
二、材料与方法
1 数据来源
TCGA：535个LUAD样本和59个正常对照样本的RNA-seq图谱数据，561个LUAD样本的WES数据，以及504个Illumina 450k DNA甲基化阵列的图谱数据
2 分析流程
1）TME构建：采用ESTIMATE算法构建TME免疫评分的中位数分为高免疫组（n = 268）和低免疫组（n = 267）；使用CIBERSORT计算每个样本的浸润免疫含量；对于WES和甲基化谱的LUAD样本，通过映射RNA-seq谱的样本ID构建了一个高免疫队列和一个低免疫队列
2）多组学数据分析：分析高免疫队列和低免疫队列之间的基因表达、体细胞突变和DNA甲基化的差异；WES数据用于检测SNVs、SNPs和INDELs；CoMEt算法识别共存和互斥突变，包maftools用于创建体细胞突变的可视化；R包ChAMP用于处理甲基化阵列数据，筛选出缺失值超过20%的样本，使用451个样本，进一步分为高免疫队列219个样本和低免疫队列232个样本；limma软件包和Bumphunter算法分别识别差异甲基化探针和区域；利用Pearson相关性研究探针信号与基因表达水平之间的相关性，并从所有探针中随机选取与真DMP集相同数量的探针，构建100个随机集
3）功能丰富性分析：R软件包clusterProfiler
4）临床相关性：Kaplan-Meier生存分析
5）预后预测模型的建立和评价：单变量Cox、lasso回归模型、多变量Cox；R软件包timeROC绘制ROC曲线评价模型
三、结果展示
01 - 在LUAD建造TME并剖析浸润性免疫内容
为了评估浸润的基质细胞和免疫细胞的肿瘤相关效应，通过使用ESTIMATE算法，根据TCGA表达谱初步建立TME。ESTIMATE产生一个被称为 "estimate score"的指数，全面推断肿瘤的纯度。如图1A所示，与正常样本相比，LUAD样本的estimatescore分布明显偏低，基质得分和免疫得分也是如此。
接下来，研究肿瘤纯度与临床因素之间的关系。如结果所示（图S1A），estimate score在肿瘤大小、远处转移和肿瘤阶段上有显著差异。对比显示，I期样本的estimate score明显高于III期和IV期样本。此外，与高分（高于中位数）队列中的患者相比，低estimate score（低于中位数）队列中的患者预后较差（图S1A）。
从TNM分期系统方面来看，不同肿瘤大小的免疫评分有显著差异，但淋巴结和远处转移没有差异(图1D-1F)。综合TNM分期分类方面，如图1B所示，早期和晚期的免疫评分有显著差异，其中I期的免疫评分显著高于III期和IV期。相比之下，基质评分仅与远处转移显著相关，而非肿瘤大小、淋巴结和分期（图S1B）。较高的免疫评分与较长的总生存时间显著相关（图1C），而基质评分与患者预后无显著相关性（图S1B）。
由于免疫浸润水平和细胞组成与肿瘤进展和患者结局密切相关，将LUAD样本利用其免疫评分中位数分为高免疫队列和低免疫队列，并利用CIBERSORT进一步表征细胞组成，探讨免疫细胞亚群与临床特征之间的关系。剔除CIBERSORT产生的P值大于0.05的样本。共保留468个样本，其中203个属于低免疫力队列，265个属于高免疫力队列。
将每个样本的免疫含量剖析为22种免疫细胞成员。高免疫队列中记忆B细胞、CD8 T细胞、活化的记忆CD4T细胞、M1巨噬细胞、静止的树突状细胞、活化的肥大细胞和γdelta T细胞的比例明显较大，而浆细胞、M0巨噬细胞和活化的树突状细胞的比例较小 (图1G)。在考虑免疫评分与上述检测到的临床因素之间的关系的同时，推测不同的免疫细胞成员可能对患者的结果有不同的贡献，检测到记忆B细胞、CD8 T细胞、M0巨噬细胞、M2巨噬细胞和活化树突状细胞在TNM阶段、肿瘤大小或淋巴结上有显著差异，而其他细胞亚群在所有临床因素中似乎没有统计学意义（图1H）。此外，尽管组合显示出辨别力(图1C)，单一类型的细胞成员对LUAD的患者总体生存率贡献不大(图1H)。
02 - 免疫浸润依赖性差异化表达基因的鉴定
利用TCGA的LUAD样本表达谱来鉴别高免疫组和低免疫组之间的表达变化，其中高免疫队列中分别有611个和164个基因上调和下调(图2A)。发现CXCR4和CCL8等29个趋化因子显著上调（图2B），这些因子能够调节多种免疫细胞对肿瘤的招募。
使用clusterProfiler进行功能富集分析，发现上调的基因富集在免疫相关的生物过程中（图2C），表明它们在增强肿瘤相关免疫力方面具有积极作用。此外，一些上调的基因如LILRB4、RUNX3和CXCR3已在先前的研究中被实验验证为调节T细胞活化和支持肿瘤浸润。另外，下调的基因主要富集在代谢过程中（图2D），推测一些下调基因利用代谢开关调控免疫细胞和肿瘤细胞的活动。
03 - 不同免疫渗透水平下的体细胞突变比较
在检测到上述部分的转录改变后，进一步研究了是否有证据表明高免疫力和低免疫力队列的基因组层存在差异。基于TCGA门户网站的WES数据，如图3A和图S3A所示，在高免疫力和低免疫力人群中，大多数基因变异都是错义突变（约60%）。从全局的角度来看，低免疫队列的样本持有的变异数量明显多于高免疫队列的样本。
SNV方面，高免疫和低免疫队列中所有样本共检测到64,344和88,708个SNV，其中C>A是高免疫队列和低免疫队列中最常见的类型。无论SNV的类型如何，低免疫队列中的突变数都显著高于高免疫队列中的突变数（图3B）。所有SNVs的转座（Tv）和转座（Ti）之间的比例约为2:1，并且在两个队列中保持稳定（图S3D）。此外，高免疫队列中的SNPs、INSs和DELs也被低免疫队列中的SNP超过（图3C）。相比高免疫力队列中的样本，低免疫力队列中的样本具有明显较高的变异等位基因分数（VAFs）水平（图3D）。尽管两个免疫队列中4种类型的体细胞突变的变异数量存在显著差异，但所有变异中所占据的每个突变类型的内成分比几乎保持不变(图S3A-3C)，这表明所观察到的突变数量的差异不是由类型转移引起的。
在低免疫力队列中，129个基因在10%以上的样本中发生了突变，而在高免疫力队列中，只有62个基因符合这一标准，其中有56个基因发生了重叠。相应队列中最常突变的15个基因见图3E，TP53、TTN和MUC16在两个队列中都占据了前三位的位置，而且它们之间存在相互作用。
接下来，利用CoMEt算法研究了前25个最常突变基因的共现和排他性突变情况，与普遍的共现情况相比，有两个队列（KRAS-TP53、KRAS-TNR和STK11-TP53）中有3个独特的病例在表现出排他性突变（图3J），这表明它们可能在同一途径中产生冗余效应，并且它们之间具有选择优势，可以保留一个以上的突变拷贝。
有些基因在两组之间有差异性突变频率。检测到268个差异突变基因，按p值升序排序，前10名如图3F所示。不同的变异可能会对患者的其他基因改变甚至临床结果产生不同的影响。没有SNPs的KRAS的表达水平在高免疫和低免疫队列之间有显著差异，但当SNP rs121913530（C>A）存在时，则相反（图3I）。此外，STK11是另一个典型的例子，以证明两个队列之间不同的突变点（图3G）和预后影响差异的合理连锁反应（图3H）。
04 - TIME描绘LUAD的DNA甲基化模式
未能维持正常的DNA甲基化，其中包括CpG岛的低甲基化和CpG贫乏区域的低甲基化，增加了触发肿瘤形成和恶化的敏感性。因此，接下来使用来自TCGA的甲基化数据检测和比较不同免疫队列中DNA甲基化模式的影响。
451个样本中共鉴定了5,764个免疫相关的差异甲基化探针DMPs（图4A）。与低免疫力人群相比，高免疫力人群共检测到高甲基化位点5647个，涉及2386个基因，其中2221个位点位于1687个CpG岛上。相比之下，低甲基化位点的数量大大超过了68个基因相关的117个位点，位于56个CpG岛上。因此，高免疫群整体上倾向于有低甲基化的位置，但低甲基化只发生在少数基因上。此外，还发现许多DMP相关基因在两个队列之间有差异表达。从高免疫队列中的2386个低甲基化基因中，有63个上调和32个下调的DEG（图4B）。然而，从低甲基化基因组中，只检测到7个上调的DEGs。
在基因本体分析的基础上，对DMP相关基因的功能进行了研究。前15位富集的生物过程中，FDR最低的GO项表明它们在细胞分化和发育中的潜在作用（图4G）。DMP相关基因的基因集富集分析（GSEA）表明，具有高度正β差异的高甲基化基因对肿瘤相关神经生物学过程有更重要的贡献（图4D），说明异常甲基化诱导的肿瘤免疫攻击行为是通过对神经通路的识别和参与实现的。
考虑到先前关于DNA甲基化和基因表达水平之间相关性的发现，推测在本研究中是否存在类似的现象，以及在不同的免疫水平上这种趋势是否稳定。结果表明，在2441个DMP相关基因中，高免疫组有329个正相关基因和926个负相关基因，低免疫组有346个正相关基因和939个负相关基因。与使用随机选择的探针构建的100个随机集的相关系数的相对平衡分布相比，与DMP相关基因相关的探针信号容易与表达水平负相关（图4C）。
免疫水平不影响甲基化水平和表达水平之间的相关性，这由图4E所示的两个队列间相关系数的高度一致性所支持，两个队列之间负相关（图5B中的Venn图）或正相关（图5A中的Venn图）基因集的大量重叠。在低免疫组和高免疫组中，这些一致的正相关和负相关基因分别在免疫系统和细胞增殖中富集，而趋势不一致的基因则有其独特的功能。例如，低免疫力队列中独特的正相关基因可以参与突触相关功能（图5A）。此外，正相关基因的探针更常位于基因体和3′UTR区域，而负相关基因的探针更倾向于与启动子相邻的区域（图4F），说明DNA甲基化对表达的影响存在区域差异。
05 - 多组学特征提供准确的预后预测方法
为了从众多的基因改变中找出免疫相关的预后信号，采用了基于lasso回归和Cox比例危险回归的策略。还分别研究了三种改变的联合效应和单独效应，以确定哪种模型的性能最好。
首先，对于联合效应，将所有的基因改变进行合并，其中由DEGs、5个突变和217个DMPs组成的337个变量，采用单变量Cox比例危险模型确定对患者的总生存时间的显著独立影响。采用lasso回归模型删除贡献较小的变量，在最优参数 (图6A)下，保留52个变量，建立多变量Cox比例危害回归模型。随机将TCGA样本分为训练集和独立测试集（n=142）。在结果中（图S6），训练模型的平均协整指数（C-index）等于0.839。其次，根据建立的模型计算每个样本的风险评分，训练集上1年、3年、5年预后预测的AUC均值分别达到0.871、0.875、0.928。在测试集的预测方面, 1年、3年和5年生存期的AUC均值分别为0.796、0.786和0.777。此外，按风险评分中位数将样本分为高风险和低风险队列（图6F）。生存分析（图6G）显示，与低风险队列相比，高风险队列的总体生存率较差。
考虑到上述构建的模型具有很好的鲁棒性和有效性，再结合所有TCGA样本，生成一个由27个变量组成的整体预测模型（图6D）。此外，与上述发现一致，根据生存分析的结果，高危队列的预后比低危队列差（图5F和5G）。不管是1年、3年还是5年的生存率，风险评分都有很高的辨别力，AUC值分别等于0.861、0.850和0.916（图6B）。
采用与上述相同的策略，判断每类基因改变的单独效应是否具有同等甚至更优的效益，以及高免疫和低免疫队列之间的差异性突变是否能够取代频繁突变基因在预测模型中的作用。从图S9所示的结果来看，无论表达变化、体细胞突变、差异性DNA甲基化，与联合模型相比，没有任何单一特征能够提供足够强大的预后预测。
此外，当将频繁的体细胞突变替换为差异性突变时，得到了一个基于23个变量的预测模型，但性能没有改善（图S10）。另外，考虑到这些差异突变的频率相对较低，最终没有保留这个模型。
除了基因改变外，一些临床因素也可能对总生存时间有预测价值。虽然阶段与总体生存时间显著相关，但与上述仅基于多组学改变构建的总体模型相比，纳入临床因素的新模型的辨别能力并没有提高，如图6E所示。此外，这个新模型在其1年、3年和5年生存预测上并没有取得更好的表现（图6C）。只有由上述27种改变组成的多组学特征已经可以产生准确的预后预测。
四、结论
在本研究中，旨在根据TCGA的表达谱，估计LUAD的TME浸润模式，尤其是肿瘤相关免疫系统，然后通过分析多组学数据(RNA-seq、全外显子组测序和DNA甲基化阵列)，将免疫状态与遗传或表观遗传特征相关联，最终从显著的改变中建立预后预测模型。期待本研究结果能够为人类LUAD提供更全面的免疫基因组学图谱，并有可能找到更好的预后预测因子。
