坏死性凋亡识别新的分子表型及其对肺腺癌肿瘤免疫微环境的影响

Necroptosis Identifies Novel Molecular Phenotypes and Influences Tumor Immune Microenvironment of Lung Adenocarcinoma

坏死性凋亡识别新的分子表型及其对肺腺癌肿瘤免疫微环境的影响

发表期刊：Front Immunol

发表日期：2022 Jul 14

影响因子：8.786

DOI:  10.3389/fimmu.2022.934494

一、研究背景

        由于肺癌的高复发率和不恰当的诊断，肺癌是全世界癌症相关死亡的主要原因。其中大部分患者为非小细胞肺癌，可分为几种组织学类型，其中肺腺癌（LUAD）是常见的亚型，占了近60%的病例，而且这个比率还在增加。

        坏死性凋亡一种由RIP1激酶和RIP3介导的程序性细胞死亡，其特点是细胞内容物的大量渗漏和细胞器肿胀。越来越多的证据表明，坏死吞噬在癌症发病机制中起着不可或缺的作用。相反，坏死诱导与免疫检查点抑制剂（ICI）相结合可以增强抗肿瘤活性，甚至在ICI耐药的肿瘤中。一些研究路线已经验证了坏死吞噬可以招募炎症免疫细胞，致力于促进肿瘤的发生。

二、材料与方法

1、数据来源

1） 共有1064名患者被纳入分析，包括来自TCGA-LUAD（n=513）、GSE13213-GPL6480（n=117）、GSE26939-GPL9053（n=116）、GSE29016-GPL6947（n=72）和GSE31210-GPL570（n=246）

2）来自癌症药物敏感性基因组学GDSC 265个聚集的小分子

3） 泛癌症的多组学数据来自TCGA和UCSC Xena，并通过GSCALite网络服务器进行了一些相应的探索

2、分析流程

流程图

三、实验结果

01 - 坏死性凋亡相关的表观遗传和转录特征

        在TCGA-LUAD队列中共发现了30个差异表达的坏死相关基因。为了揭示LUAD中坏死调节因子的表观基因组学特征，作者分析了甲基化水平。结果表明，DNA启动子甲基化与TERT和GATA3的转录表达水平之间存在明显的负相关，而BACH2与其甲基化有微弱的正相关（图1A）。由于CNV可以在臂水平上影响整个基因组或在焦点水平上改变一个广泛的区域，引发对体细胞突变的更大影响，进一步描绘了这些基因的CNV景观。拷贝数值的探索显示了30个调节因子的异常改变，其中TERT和OTULIN有广泛的CNV扩增频率，而CDKN2A在LUAD中表现出显著的CNV缺失（图1B）。根据cBioPortal的泛癌突变分析，CDKN2A表现出最高的突变类型，有深度缺失，其次是TERT、TRIM和ZBP，扩增频率为11%。在561个LUAD样本中，111个出现了坏死相关基因的突变，频率为19.79%。然而，说明ALK的突变频率最高，其次是CDKN2A，而剪接位点只在CDKN2A和IDH1中观察到（图1C）。

图1 肺腺癌中坏死性凋亡因子的表观遗传学改变景观

        研究发现，大多数调节因子在LUAD的癌旁组织和肿瘤组织中都有明显的表达差异。有几个基因在肿瘤病例中的表达水平较低，如AXL、BAHC2、CFLAR和DNMT1，而其他基因则过度表达，如BNIP3、DIABLO、FADD和IDH1（补充图2A）。此外，BNIP3、DIABLO、DNMT1、FADD、LEF1、PANX1、PLK1、TNFRSF21和ZBP在I至IV的不同病理阶段的表达有明显变化（补充图2B）。对这30个基因进行了层次聚类，得到了四个聚类：聚类A（TERT、MYCN和ALK），聚类B（AXL、BNIP3、CFLAR、DIABLO、DNMT1、FADD、ID1、ITPK1、KLF9、MPG、PANX1、TLR3、TNFRSF1B、TRAF2、TRIM11），以及聚类D（BACH2、CDKN2A、GATA3、LEF1、PLK1和ZBP1）。为了充分了解LUAD中坏死相关基因的调控机制，收集了1,064个样本，并使用'ComBat'方法校正了批次效应。五个队列在校正前后的基线显示（图2A）。可以看出，成分1和2的总和从90.9%下降到14.9%。集合体分布也表明，批量效应被大致去除。用网络图描述了坏死凋亡分子与遗传相似性相互作用的综合景观及其对LUAD患者的预后价值( 补充图2C )。发现同一类别的坏死调节因子在转录表达水平上呈现出明显的正相关，并且在四种修饰模式中表现出明显的关联性。DIABLO、FADD、ID1、MPG、TRIM11、PLK1、TNFRSF21和TERT的过度表达与总生存率低有关，即风险因素，而ALK被认为是一个有利因素。单变量Cox模型显示，FADD、PANX1、PLK1、TRIM11、TRAF2、TERT和ID1可能作为风险因素，而ALK对预后表现出预防作用（补充图2D）。

补充图2 肺腺癌的坏死调节分子和临床特征

02 - 通过无监督聚类确定LUAD的新亚型

        基于K-means的无监督学习方法应用29个坏死性凋亡相关基因的对患者进行分类。根据不同聚类的分布和累积分布函数（CDF），最终确定了最佳聚类和三种不同的坏死凋亡，其中聚类A为210例，聚类B为578例，聚类C为276例。PCA证实，三种聚类亚型可以通过29种坏死相关分子的表达水平来区分（图2B），热图直观地显示了坏死凋亡作用调节因子的表达变化。坏死凋亡相关基因在C群中表现出最高的表达水平，其次是B群，然后是A群（图2C）。还同时检查了三个坏死凋亡亚群的临床病理表型的分布，包括年龄、性别、吸烟史、复发情况和病理阶段，发现亚群B和C有更多的女性患者，群组A的患者更有可能吸烟，早期阶段（I期或II期）的患者大多来自B群。然而，患者的年龄在各表型之间的分布没有统计学意义。三种主要坏死性凋亡表型的KM生存分析显示，簇B的预后优势相对突出（图2D）。单变量和多变量的Cox比例危险回归模型表明，B群可以预测LUAD的总生存期（图2E）与C群一起，此外，年龄和阶段也是生存的基本预测因素。

图2 基于平均数链接和K-均值的无监督聚类，确定三种坏死性凋亡表型

03 - 肿瘤免疫微环境在不同的坏死性凋亡表型中的变异

        由于癌症免疫在肿瘤进展和细胞增殖中起着至关重要的作用，作者推测不同亚群亚型的肿瘤免疫微环境可能是与众不同的。首先调查了坏死性凋亡表型与多种免疫调节剂和免疫检查点之间的相关性。坏死凋亡对肿瘤免疫微环境的影响结果表明，亚群B和C中存在明显上调的趋化因子的表达水平，如CCL11、CCR1、CXCL10和PPBP（图3A，上部分）。白细胞介素、干扰素、受体和其他一些细胞因子的表达也在簇B和簇C中出现，而簇A的免疫分子水平强烈下调（图3A，中间和底部部分）。与簇 A 和 B 相比，簇 C 在七种免疫分子中表现出最高的表达，这些免疫分子具有最高的表达，这些免疫分子具有免疫治疗的潜在靶标，包括 CD274 （PD-L1）、PDCD1 （PD-1）、CD247、PDCD1LG2 （PD-L2）、CTLA-4、TNFRSF9 和 TNFRSF4（图 3B）。

        进一步探究了20种免疫细胞类型的浸润水平和免疫周期中每个步骤的水平。对基因特征的探索表明，免疫细胞的浸润在亚群B中明显较高，其次是亚群C，而在亚群A中发挥的作用明显较低，如CD4+幼稚T细胞、滤泡辅助T细胞和M1巨噬细胞。值得注意的是，活化的记忆型CD4+T细胞在簇B中的丰度最低，而在簇A和簇C中的丰度较高。通过ssGSEA算法计算出在消灭肿瘤细胞的每个步骤中发生的抗肿瘤免疫反应事件所定义的七个免疫周期信号的分数。在聚类B和C中都显示出较高的分数，但在聚类A中显示出特别低的水平（图3B）。此外，亚群A也具有明显的MeTIL弱点，这表明肿瘤浸润免疫细胞的比例很低（图3B）。

        还采用了TIDE算法和ESTIMATE方法来比较TME的免疫亚型中免疫和基质成分的异质性。因此，可以量化坏死凋亡和潜在的免疫疗效之间的相关性。研究发现，簇B的TIDE和功能障碍得分最高，其次是簇C，然后是簇A，而在排除分析中发现了一个相反的结果。聚类B和C的免疫和基质得分明显高于聚类A，表明聚类A显示出最高的肿瘤纯度（图3B上部）。此外，调查了几个免疫途径的生物标志物。结果显示，功能的归一化富集分数明显较高，主要位于簇C，包括共同抑制APC、检查点、促进炎症、MHC I类、副炎症、共同抑制/刺激T细胞和I型IFN反应（图3C）。聚类B显示出相对较高的APC共刺激、CCR、细胞溶解活性、HLA和II型IFN反应的得分，而聚类A的丰度最低（图3C）。通过相关策略，几个坏死调节基因，包括ALK、AXL、KLF9、LEF1、TNFRSF1B和ZBP1，与大多数免疫细胞呈正相关，而BNIP3、DIABLO、FADD、IDH1、IDH2、PLK1和SLC39A7等分子呈负相关（图3D）。

图3 三种坏死性凋亡表型中的免疫调节剂和几种定量的肿瘤免疫微环境的改变

04 - 坏死性凋亡表型中基因组改变的划分

        作者进一步探讨了已确定的LUAD分子表型的基因组差异，包括突变特征、SNV和转录因子调节网络的活性。推断有四个突变特征与LUAD有关，即SBS1（年龄相关）、SBS2（APOBEC活性相关）、SBS4（吸烟相关）和SBS5（ERCC2突变相关）。然而，仅在SBS1中检测到统计学意义，显示与群组A相比，群组B有更多的年龄相关特征的突变，与群组C相比，群组A在吸烟相关特征的突变频率更高，其次是群组B（图4A）。因此，这些突变特征可以作为每个分子坏死性凋亡表型的扰动的解释。在前20个经常突变的基因中，发现C群在TP53、TTN、MUC16、CSMD3和RYR2的突变中更为丰富，而B群与A群和C群相比，突变频率相对较低（图4A）。B组和C组的拷贝数扩增和缺失的负担下降，但A组在病灶和手臂水平上都明显增加，B组的CNV负担明显降低（图4C）。探讨了CNV和GISTIC分数在不同染色体上的分布频率，验证了结果，并显示增益和损失之间最极端的频率差异位于8号和14号染色体上。

        为了进一步分析亚群之间的转录组差异，作者共调查了8个LUAD特异性转录因子和15个与癌症染色质重塑相关的基础调节因子。证实了三个坏死凋亡亚群的生物学相关性，因为调节子的活动与亚群密切相关（图4D，上部）。亚群A和C共享相似的调节子活性模式，而亚群A则与FOXA2和ATOH8的高活性不同，而亚群C则与激活的HOXA4明显相关。Cluster B包含高浓度的调节子，如HOXA4、DACH1、EPAS1、ETV5、FOXA2、ATOH8和SMAD6。与癌症染色质重塑有关的调节子活动信息强调了三种坏死性凋亡表型之间可供选择的不同调节模式，表明表观遗传驱动的转录网络也许是这些亚群的重要区分因素（图4D，底部分）。发现亚群A和C也有类似的调节子模式，而亚群A表现出更高的SIRT1、SIRT4和SIRT5的活性。同时，SIRT1、CLOCK、EP300、HDAC5、SIRT4、SIRT2、HDAC10和HDAC8的调节子在亚群B中明显富集。

图4 三种坏死性凋亡表型之间的基因组特征的区别

05 - 信号通路的功能紊乱

        作者探讨了坏死的三种表型中基因表达的差异，并分析了几种与癌症有关的信号通路的特点，这有助于了解潜在的调节机制和可药用的通路。坏死调节分子与激活的途径有关，如细胞凋亡、细胞周期A、EMT-A和激素AR-A，但在泛癌中一直抑制着细胞周期I、激素AR/ER-I和RTK-I等功能（图5A）。在LUAD中，大多数坏死相关基因激活了细胞周期、EMT和细胞凋亡等途径，而它们抑制了激素ER和RTK信号（图5B）。使用语义相似性来聚合紧密相关的类别，以筛选在坏死凋亡中起作用的30个调节因子中的中心基因。结果表明，TNFRSF21的权重最高，其次是DIABLO和TRAF2（图5C）。为了分析这些不同的坏死性凋亡表型的生物学行为，进行了GSVA富集的探索。如图5D所示，亚群A在RNA聚合酶、氨基酰tRNA生物合成、同源重组和DNA复制中明显富集。群组B包括与病毒性心肌炎、细胞粘附分子（CAMs）、自然杀伤细胞介导的细胞毒性和造血细胞系相关的生物途径。此外，C组在原发性免疫缺陷、产生IGA的肠道免疫网络、细胞周期、DNA复制、错配修复和碱基切除修复方面有突出的富集。我们还评估了不同聚类亚型中的其他癌变途径。NOTCH、NRF2、RAS和TGFβ的得分在聚类B中明显偏高，而细胞周期和TP53信号通路在聚类C中被激活（图5E）。

图5 对坏死和癌症相关信号通路的调查

        此外，研究了这三种表型与cGAS-STING信号通路之间的联系。它表明，cGAS、IRF3和TBK1在C簇中表达水平最高，但在A簇中下调最多（补充图7A-D），这意味着坏死细胞在塑造肿瘤免疫微环境和免疫反应中的意义，并进一步提示坏死细胞在化疗和免疫疗法中治疗LUAD的潜在可行性。

补充图7 通过Kruskal-Wallis和Welch单因素方差分析试验检查的三种坏死性凋亡表型中cGAS-STING相关基因的水平

06 - NecroScore的临床价值

        作者根据与坏死性凋亡相关的分子开发了一个名为NecroScore的评分系统，以衡量LUAD患者的坏死性凋亡水平。在TCGA和GEO-meta队列中，患者被独立分为低和高NecroScore组，其中TCGA-LUAD的阈值为-3.80，GEO的阈值为-0.22。通过比较三种聚类亚型的NecroScore水平，显示NecroScore在聚类C中最高，其次是聚类A，然后是聚类B，通过Welch单因素方差检验（图6A）。根据TCGA（图6B）和GEO队列（图6C）的KM生存分析，NecroScore高的患者比NecroScore低的患者预后更差。不仅如此，不同临床分层的患者，包括年轻人和老年人、女性和男性、早期阶段和晚期阶段也有生存意义。在高NecroScore组中发现了更高的不吸烟率和更多的晚期患者，而年龄和性别分布没有统计学差异（图6D）。单变量Cox模型显示年龄、病理分期和NecroScore与生存结果显著相关（图6E，上部），多变量Cox回归验证了年龄，阶段和NecroScore可以作为独立的预后因素（图6E，底部）。通过t检验，低NecroScore患者与高NecroScore患者相比有明显的治疗优势（图6F）。在三个独立的预后因素中，Stage表现出最高的净效益，其次是年龄，然后是NecroScore（图6G）。此外，根据决策分析，三个变量的组合在3年、5年和10年内具有最高水平的净效益（图6G）。作者建立了一个nomogram模型来预测LUAD患者的1年、3年和5年的生存率（补充图9A）。校准曲线验证了nomogram模型对总体生存结果的预测是乐观的（补充图9B）。还发现，低TMB和高NecroScore组合的患者预后不佳，而高TMB和低NecroScore的患者则表现出明显的优势（图6H）。

图6 NecroScore的临床意义 补充图9   综合诊断模型

        进一步研究了免疫细胞浸润水平之间的相关性。基于TIMER、CIBERSORT、CIBERSORT-ABS、ESTIMATE、QUANTISEQ、MCP-counter和XCELL的免疫反应热图，使用ssGSEA方法，如补充图10A所示。它显示，几种细胞类型的浸润丰度与NecroScore呈负相关，如骨髓树突细胞、癌症相关的成纤维细胞、造血干细胞和单核细胞。从另一项研究中收集的免疫细胞的不同基因标志物被用来验证这一结果。大多数标志物呈现负相关，特别是在内皮细胞、静止的树突状细胞和静止的肥大细胞中，在激活的记忆CD4+T细胞和成纤维细胞中可以发现微弱的关联（补充图10B）。就免疫表型而言，高新陈代谢分数组中C1、C2和C6的比例较高，而低新陈代谢分数组中C3和C4的样本较多(补充图10C)。来自TCIA的证据显示，低NecroScore患者对PD-1无反应的抗CTLA-4阳性或阴性状态的IPS明显较低(补充图10D )，表明抗CTLA-4抑制剂也可能对低分的PD-1阴性状态患者有影响。此外，发现几个免疫检查点，如CD274、PDCD1、PDCD1LG2、LAG-3和TMIGD2，在低NecroScore组中被下调。有趣的是，在低NecroScore组中，TIDE评分、MSI和功能障碍评分都明显增加。同时，低NecroScore患者的排除分数较低，表明免疫检查点负表达或弱表达的个体也可能从免疫检查点阻断（ICB）治疗中表现出有希望的结果。

补充图10 免疫特征和NecroScore之间的相关性

07 - 低和高NecroScore组的突变状态

        为了分析LUAD的NecroScore相关机制，对TCGA-LUAD队列的体细胞突变进行了调查。比较两个评分组的突变频率，在高NecroScore组发现更多的突变（补充图12A），包括非同义和同义突变（补充图12B，C）。同时，在高NecroScore组中观察到TP53、TTN、PCDH15、LRP1B、ZFHX4、NALCN、CSMD3、MUC16、XIRP2、APOB、PAPPA2、PRDM9、ANK2和LRRC7的突变频率明显较高，而KEAP1在低NecroScore组中有更多的突变，由森林图显示（补充图12D）。此外，这些基因之间存在着明显的共现性（补充图12E）。还发现了高突变基因变异对NecroScore的影响，发现NecroScore与坏死相关基因CDKN2A突变频率呈正相关（补充图12F）。

补充图12 NecroScore与肿瘤突变状态相关联

08 - 药物发现和化疗反应分析

        作者分析了坏死性凋亡相关分子与LUAD治疗的临床效率之间的关联。根据GDSC的药物反应数据显示，一些基因，如BAHC2、LEF1、DNMT1、TERT、TNFRSF1B、TNFSF10、MPG和ALK与药物有协同作用，而PANX1、AXL、SLC39A7、TNFRSF21和ID1与药物有拮抗作用（图7A）。ALK表现出与CH5424802（艾乐替尼、ALK和RET抑制剂，图7A）的明显协同作用。同样明显的是，ID1与vorinostat（HDAC抑制剂）、PHA-793887（CDK2抑制剂）、PIK-93（PI4K抑制剂）和SNX-2112（HSP90抑制剂，图7A）等药物的负相关关系最强。

        鉴于坏死性凋亡调节分子在化疗中的重要性，进一步分析了NecroScore是否能准确预测LUAD患者的化疗敏感性。利用脊回归来预测低和高NecroScore组的IC50。结果显示，高NecroScore的几个LUAD特异性药物如顺铂、多西他赛和紫杉醇的IC50较低。相反，低NecroScore的患者可能对厄洛替尼（酪氨酸激酶抑制剂/TKI）、吉非替尼（EGFR-TK抑制剂）和拉帕替尼（TKI）反应更敏感（图7B），表明化疗是高NecroScore组的一个有希望的选择。相比之下，低NecroScore组的靶向治疗可能导致有利的结果。比较了基于坏死性凋亡表型的分类器得到的细胞系群内药物反应的AUC。只有在90%以上的LUAD细胞系上测试的274种化合物被用于分析（补充表5-7）。SB-590885（一种B-RAF抑制剂）和UNC0642（一种G9a/GLP抑制剂）的AUCs在第3组LUAD细胞系中明显降低，但681640和KN001-27的信息是不存在的（图7C）。

        此外，确定了可以针对干性特征和NecroScore差异的潜在化合物。对24个化合物的分析显示，19个MoA被这些化合物所共享（图7D）。四个化合物(氟苯尼考、氟苯尼考、美索里达嗪和匹莫齐特)共享多巴胺受体拮抗剂的MoA；甲氧苄啶和奥芬那定共享乙酰胆碱受体拮抗剂的MoA；金刚烷胺和美金刚共享谷氨酸受体拮抗剂的MoA。此外，还可以看到帕莫西汀（pyrvinium-pamoate）是AKT抑制剂的MoA。这些结果有助于根据NecroScore并参考突变特征，为每个分层的LUAD组采取更有针对性的化疗策略（补充图12）。

图7 以坏死为基础的亚型的药物敏感性和化疗潜力与NecroScore之间的关联性

四、结论

        这项研究描绘了LUAD肿瘤微环境中坏死分子的分子和免疫景观。三种坏死性凋亡表型可以区分LUAD的临床特征、免疫浸润和体细胞突变特征。在一定程度上，开发的NecroScore系统可以预测预后和化疗效率。作者用多组学方法对个体肿瘤坏死形态进行综合评估，将加深对肿瘤免疫微环境的理解，并制定更有效的免疫治疗策略。