随机森林构建肺腺癌的TME表型模型
Identification and validation of tumor environment phenotypes in lung adenocarcinoma by integrative genome-scale analysis
肺腺癌肿瘤环境表型的整合基因组分析鉴定与验证
发表期刊：Cancer Immunol Immunother
发表日期：2020 Jul
影响因子：5.442
DOI: 10.1007/s00262-020-02546-3
一、研究背景
肺癌是全球最常见的癌症，也是男性和女性癌症死亡的主要原因。继肺鳞状细胞癌之后，肺腺癌（LUAD）是目前最常见的非小细胞肺癌（NSCLC）亚型，占肺癌病例的40%以上。
目前正在研究包括PD-L1、肿瘤突变负担和预先存在的细胞毒性在内的生物标志物，这些生物标志物可以预测对免疫检查点阻断的反应性，以进一步提高精准免疫治疗。
实体肿瘤和浸润其中的免疫细胞在动态平衡中相互作用，形成疾病进展。越来越多的研究强调了肿瘤环境（TME）的关键作用，它不仅包括先天性和适应性免疫反应的不同免疫细胞类型，还包括一系列与免疫相关的细胞因子，在癌症进展和肿瘤逃出免疫监视的过程中起着关键作用。
二、材料与方法
1 数据来源
1）发现队列（n=514）：从UCSC Xena下载LUAD患者的4级基因表达数据(FPKM归一化)和TCGA的相应临床信息。数据用R Bioconductor软件包进行分析。
2）验证队列（n=1063）：GSE30219、GSE31210、GSE3141、GSE37745、GSE50081和GSE68465。使用sva软件包中的COMBAT方法消除批次效应，并用limma软件包进行背景调整和数据归一化。
3）回顾性地选取了2016年至2017年在复旦大学附属中山医院的34例LUAD患者。使用Illumina Hiseq 2500和BGI-500RNAseq平台对所有肿瘤样本进行RNA测序。
2 分析流程
流程图
三、结果展示
01 - LUAD中的微环境表型
为了选择最佳和稳定的聚类数，作者根据24个细胞群的GSVA分数使用ConsensusClusterPlus包进行无监督共识聚类（图S1B）。k-means聚类的结果表明，514例患者分为两个异质性聚类（基于欧氏距离和Ward's linkage），分别命名为TMEcluster A（242例患者）和B（272例患者），反映了每个患者的适应性免疫系统和先天性免疫系统24个免疫细胞的不同浸润模式（图1A）。可以观察到两个TME群之间OS的显著差异（图1B）。
TMEcluster A，即 "冷区"，其特点是微环境细胞浸润相对较低，因此，我们为该簇选择了蓝光来描绘免疫细胞的稀疏分布。相比之下，TME集群B，"热区域 "的特点是先天性和适应性免疫反应的激活。因此，作者为该群选择了一个 "红灯 "来表示丰富的免疫浸润（图1A，C）。通过Wilcoxon检验证实了两个TME簇之间免疫细胞浸润模式的显著差异。24个免疫细胞亚群的丰度之间几乎是全方位的正相关，见图S1C。
图S1
如森林图（图1D）所示，除Th2、中性粒细胞和aDCs外，大部分肿瘤浸润的免疫细胞与较好的OS相关。这一发现与TMEcluster B的患者预后较好一致，其中免疫细胞的丰度明显高于TMEclusterA。在临床特征方面，调查了TME簇的年龄、性别、阶段和病理亚型的分布情况（图1E）。
图1 TCGA队列中的LUAD TME群组
02 - 随机森林分类器模型的构建与验证
为了开发一种临床适用的方法，将不同人群的LUAD患者划分为TMEcluster A或B，作者建立了一个基于Breiman随机森林算法的分类模型。
首先，将TCGA中的514名患者随机分配到训练队列（n = 360）或测试队列（n = 154）中，同时将上述24个免疫细胞群的相对丰度设置为输入变量，并将无监督聚类确定的TME表型（A或B）（图1A）设置为响应因子。在训练队列中，经过1万次迭代，选择最佳mtry和ntrees分别为5和5000。将其拟合到训练队列和内部测试队列中验证其性能，发现正确预测率分别为100%和95.58%。
接下来，将该随机森林模型应用于两个外部验证队。计算24个免疫细胞子集的GSVA评分，并将数据拟合到如上所述的分类器模型中后，将患者分层为TMEcluster A和B，如热图和相应的box图所示，两个外部验证队列的TMEcluster B中也观察到先天性和适应性免疫系统的免疫细胞富集，这与TCGA患者定义的免疫浸润模式一致（图2A-D）。同样，TMEcluster B也与两个队列中更好的OS相关（图2E，F）。
图2 随机森林分类器模型在GEO和我院数据中的应用验证
03 - DEGs的识别和功能注释
为了确定两个TME簇的潜在生物学特征，对DEGs的表达谱进行了检查。在TME簇B中共有1220个基因上调，80个基因下调（图3A）。通过R clusterProfiler软件包对1300个DEG进行GO和KEGG富集分析，发现免疫相关途径在TME簇B中富集，这支持了TME簇B中 "热而丰富 "的免疫细胞浸润模式（图3B，C）。在TMEcluster A中，排名靠前的富集途径是由几个经典的代谢途径驱动的。由GSEA算法计算的富集图也显示了类似的结果（图S3）。
图3 差异表达基因(DEGs)和TME簇之间的功能注释
图S3 GSEA算法计算的富集图
作者选取了已发表的参与癌症发展的生物通路标志的综合集合，用ssGSEA推断其在两个TME簇中的差异激活模式。如图S4所示，除了与免疫相关的通路外，Wnt、Notch、RAS、EMT TGF-β、能量代谢等多个通路在TME簇B中也被上调，说明免疫浸润与肿瘤发生和侵袭相关的经典信号通路之间存在复杂的相互作用。
图S4
04 - TME集群的免疫微环境特征
为了进一步描绘两个TME集群的遗传和分子环境特征，作者全面研究了514个LUAD样本中几个免疫相关基因、细胞因子和miRNAs的表达谱。首先，采用了POPLAR试验中设计的7个基因面板，作为量化效应T细胞浸润水平和IFN-γ相关细胞毒性的代用指标。其次，还考察了Rooney等人定义的细胞溶解活性评分，该评分以PRF1和GZMA的几何平均值计算，以反映抗肿瘤反应的大小。
如图4A所示，8个参数在TMEcluster B中的表达水平均显著升高，说明这些患者的细胞毒功能得到了有效提高。比较两种表型中可能参与启动先天性免疫的分子，表现出相似的趋势（图4B，左）。此外，与TMEcluster A相比，TMEcluster B表现出MHC-I/ II相关抗原呈递分子表达的富集（图4B，右）。
TMEcluster B不仅有丰富的活跃的先天性和适应性免疫细胞，而且还有免疫抑制细胞（图1A和图S1C）。因此，测定了CD8+T细胞/Treg比值，以探究免疫激活和抑制的相对幅度（图4C）。该比值在TMEcluster B中显著较低，表明TME中效应T细胞的招募所启动的潜在反馈机制。
作者还研究了几种免疫调节剂在TMEclusters中的表达模式，包括15种免疫检查点分子（图4D，左）和20种共刺激分子（图4D，右）。如热图所示，在TMEcluster B中观察到了相对较高的共刺激分子和共抑制分子的表达（图4E），表明这些患者可能从免疫检查点抑制剂中获益。当调查两个外部验证队列中免疫相关分子的表达谱时，也观察到了类似的趋势（图S6）。免疫检查点、免疫细胞和突变负荷的相关矩阵如图4F所示。
图4 TCGA患者TME簇的免疫微环境特征
图S6 GEO和本院患者TME群的免疫微环境特点
05 - 与TME集群相关的肿瘤基因组改变
接下来作者研究了两个TME群中体细胞突变和拷贝数变异的分布模式。在TME簇A和B中，平均每个样本分别发现288.12和224.15个体性突变（图5A），说明突变负荷不能完全解释免疫浸润的不同模式。通过对测序数据的挖掘，在TMEcluster A和B中发现了735个不同的突变基因。
在LUAD中，TP53、EGFR、KRAS和ALK是重要的肿瘤驱动基因，也是最主要的突变基因之一。然而，虽然在TMEcluster A中略高，但这4个基因在两个TME群中的突变频率均无显著差异。与TME簇A相比，另外6个突出的癌症相关基因包括ROS1、KEAP1和STK11在TMEcluster B中具有较低的体细胞突变率，而其他的体细胞在两个群体之间具有一致的突变率(图5B)。
图5 与TME集群相关的肿瘤基因组改变
06 - TME集群的预后价值
最后，进行了单变量和多变量分析，评估TME簇和几个临床因素的预后价值。结果显示TME簇B独立预测OS优于TME簇A。通过计算Harrell's C-index，作者发现，与仅实施年龄和分期的模型相比，在比例危险模型中加入TME表型可显著提高其预后疗效，说明基于免疫细胞丰度的TME簇是不同人群中一个稳健且独立的预后因素。
四、结论
总的来说，本研究通过应用机器学习方法和多组学剖析，描述了LUAD免疫浸润模式的全面景观，并整合了之前提出的几种与不同免疫表型相关的生物标志物，从而揭示了肿瘤如何与免疫微环境相互作用，并可能指导LUAD患者更精确和个性化的免疫治疗策略。
