m6A regulator-based methylation modification patterns characterized by distinct tumor microenvironment immune profiles in colon cancer
基于m6A调节器的甲基化修饰模式，在结肠癌中表现出不同的肿瘤微环境免疫特征
发表期刊：Theranostics
发表日期：2021 Jan 1
影响因子：8.579
DOI: 10.7150/thno.52717
一、研究背景
N6腺苷（m6A）的甲基化，在mRNA、lncRNA和miRNA中广泛观察到，是最常见的RNA修饰类型，在多种生理过程和疾病进展中起着至关重要的作用。m6A修饰也是一种动态的、可逆的过程，它由不同类型的调节蛋白控制：甲基转移酶（"写入者"）、去甲基化酶（"擦除者"）和结合蛋白（"读取者"）。
结肠癌（CC）是最常见的恶性肿瘤之一，仍然是全球癌症死亡的主要原因，30％ - 50％的患者在治疗后5年内出现复发，转移，甚至死亡。随着对肿瘤微环境背景的多样性和复杂性的认识增加，通过全面解析肿瘤微环境的组成部分确定的肿瘤免疫表型将有助于指导和预测免疫治疗反应性。
二、材料与方法
1 数据来源
1）GEO数据库：GSE39582/CIT队列（N = 557）、GSE37892（N = 130）、GSE14333（N = 226）
2）TCGA：TCGA-COAD（N = 394）数据集，包括RNA测序数据、基因组突变数据
3）基于ICI的队列的基因组和临床信息：三个免疫治疗队列：用nivolumab（抗PD-1 mcAb）或ipilimumab（抗CTLA-4 mcAb）治疗的转移性黑色素瘤，以及用atezolizumab（抗PD-L1 mcAb）治疗的转移性尿道癌（mUC）
2 分析流程
1）聚类：用NMF进行共识聚类，根据23个m6A调控因子的表达确定不同的m6A修饰模式（R包'NMF'）。
2）功能分析：R软件包'GSVA'的GSVA分析，R包'clusterProfiler'中对m6A表型相关基因进行GO注释（FDR<0.01）。
3）ssGSEA量化肿瘤微环境中28种免疫细胞类型的相对丰度，CIBERSORT估计22个不同白细胞亚群的丰度与结肠癌的基因表达谱。
4）免疫表征、TIDE和ESTIMATE
5）显著突变的基因和肿瘤突变特征筛选：利用MutSigCV算法来识别显著突变的基因（SMGs），具有统计学意义（q <0.1），并在人类癌症细胞系百科全书（CCLE）认证的基因被确认为SMG。
6）鉴定不同m6A修饰表型之间的DEGs：R包'limma'（P<0.001）
7）m6Sig分数的构建：单变量Cox回归模型，使用随机森林的递归特征消除（RFE）和'caret'软件包中的10倍交叉验证，进行PCA分析提取主成分。m6Sig score= ∑(PC1i+PC2i)
三、结果展示
01 - 结肠癌中m6A调节因子遗传变异的景观
作者调查了23个m6A RNA甲基化调控基因在CC中的作用。图1A显示，这些m6A调节剂的动态可逆过程。进行GO富集和Metascape分析，显著富集的生物过程总结在图1B和图S1A中。测定了23个m6A调节因子在CC中的体性突变流行率（图1C），还研究了所有m6A调节剂的突变共存情况（图S1B）。
进一步分析23个m6A调节因子，发现CNV突变普遍存在（图1D）。23个m6A调节子的CNV改变在染色体上的位置如图1E所示。此外，基于配对的肿瘤-正常标本进行PCA，发现23个m6A调节子完全区分了CC样本和正常样本(图1F)。进一步分析表明，ALKBH5、METTL14、YTHDC2和YTHDF3在肿瘤样本中显著下调，而CBLL1、ELAVL1、HNRNPA2B1、HNRNPC、IGF2BP1、KIAA1429、LRPPRC、METTL3、RBM15和YTHDF1在肿瘤样本中显著上调（图1G）。
与正常对照样本相比，CC样本中具有CNV扩增的m6A调节因子的表达明显增加（图1D,1G）。此外，作者评估了这些m6A调节剂之间的相互调节（图S1C）。采用Cox回归分析来确定这些m6A调节剂与CC患者预后的关系。森林图显示，YTHDC1可视为保护性因素，与无复发生存期延长显著相关，而IGF2BP1被公认为危险因素（图S1D-E）。
图1 结肠癌中m6A调控因子基因改变
图S1 23种m6A调节剂的相关性和预后分析
02 - 鉴定23个调节剂介导的m6A甲基化修饰模式
在m6A调节器网络中全面展示了23个m6A调节器的相互作用、调节器连接及其在CC患者中的预后意义（图2A）。作者利用NMF算法的共识聚类分析，根据23个m6A调节器的表达情况，对具有质地不同的m6A修饰模式的样本进行分层。确定了三个不同的修饰模式集群，包括221例模式集群1，530例集群2和162例集群3（图2B，图S2A-C）。将这些集群称为m6A-C1、m6A-C2和m6A-C3，其中m6A-C1表现出突出的生存优势，而m6A-C3在meta-GEO队列中预后最差。
图2 m6A甲基化修饰模式及相关的生物途径
图S2 23个m6A调节器的无监督的聚类在meta-GEO和TCGA队列中的应用
此外，在TCGA-COAD队列中进行了相同的分析，也得到了类似的结果（图2C,图S2D）。多变量Cox比例危险性回归分析进一步证明，调整临床病理因素后，该修饰模型与患者生存结果相关（图S3A，B）。不同的m6A修饰模式之间的m6A调节剂的表达有显著差异（图S2C-D）。
图S3 m6A簇的预后价值以及三种m6A修饰模式中不同标志的富集
03 - 以不同免疫景观为特征的m6A修饰模式
为了探索三种不同的m6A修饰模式背后的生物分子变化，针对Hallmarker基因集进行了GSVA富集分析（图2D）。用ssGSEA比较28个免疫浸润细胞亚群在不同m6A修饰模式中的相对丰度（图3A）。还用CIBERSORT进一步表征了免疫浸润情况，在这个m6A甲基化修饰分层中也观察到了一致的结果（图3B）。
在m6A-C2亚型中，基质活化明显增强（图S3C-D）。使用ESTIMATE算法对三种修饰模式的免疫细胞整体浸润和肿瘤细胞纯度进行量化。m6A-C1表现出最高的免疫得分，其次是m6A-C2和m6A-C3（图3C，上图）。相反，m6A-C3的肿瘤纯度高于m6A-C2和m6A-C1，说明m6A-C2和m6A-C1亚型肿瘤周围有更多的非肿瘤成分（图3C，下图）。
Marisa等(CIT cohort/GSE39582)将CC患者分层为4种优势分子亚型(CIN、CSC、dMMR和KRASm)，CIN亚型的患者主要集中在m6A-C2和m6A-C3，而dMMR亚型主要集中在m6A-C1肿瘤内（图3D）。作者还比较了不同m6A修饰群中PD-L1的表达水平，并观察到m6A-C1亚型的PD-L1显著上调（图3E）。
图3 不同m6A修饰模式下的TME特征
基于这些发现，作者证实三种m6A修饰模式具有不同的免疫浸润特征。m6A-C1被识别为以免疫激活和丰富的免疫细胞浸润为特征的免疫炎症型；m6A-C2被识别为以基质激活和减弱的免疫细胞浸润为特征的免疫排斥型；m6A-C3被识别为以免疫抑制为特征的免疫荒漠型。
通过使用Spearman相关分析（图S4A）进一步研究了每个m6A调节因子与免疫细胞浸润之间的具体关联。IGF2BP1与YTHDF1之间存在较强的正相关关系（图2A）。IGF2BP1表现出与活化DC，未成熟DC和CD8 + T细胞的浸润水平显著负相关（ 图S4A-B）。
将IGF2BP1分层为高表达与低表达亚组，GSEA分析表明，IGF2BP1低表达的样本显示出参与干扰素γ/α反应、TNFα通过NF-κB和炎症反应信号通路的基因富集（图S4C）。ESTIMATE算法在IGF2BP1低表达亚组中表现出更高的免疫评分，证实了上述发现（图S4D）。四类免疫原性分析显示，IGF2BP1低表达亚组的特点是交叉呈递相关的MHC分子和免疫效应细胞上调，但免疫抑制细胞下调（图S4E）。此外，具有IGF2BP1突变的肿瘤比没有IGF2BP1突变的肿瘤藏有明显更多的TML（图S4F）。
图S4 TME浸润和m6A调节剂与IGF2BP1在结肠癌中的作用的关系
04 - 结肠癌中m6A表型相关的DEGs
共524个代表三种m6A修饰模式的关键区分指标的DEGs被认为是m6A相关的签名（图4A）。对这些签名基因的GO富集分析显示，与RNA修饰、转录和免疫调节相关的生物过程相关（图4B）。进行了无监督共识聚类分析，得到了三种稳定的转录组表型（图S5A-B）。这些分层将患者分为三个不同的m6A基因特征亚组，这些亚组具有不同的临床病理特征，被定义为m6A基因-S1、m6A基因-S2和m6A基因-S3（图4C）。
临床分期较晚的患者以m6A基因-S3亚组为代表，CIN亚型和PD-L1表达下调的患者主要集中在m6A基因-S2和S3亚组（图4C，图S5C）。生存分析表明，CC样本中的三种m6A基因标志有明显的预后差异。证明m6A基因-S1签名与较好的预后相关，而m6A基因-S3与较差的生存结果相关（图4D）。在考虑年龄、性别、MMR和阶段后，m6A基因标志与生存期的关联仍有统计学意义（图4E）。还比较了三个基因特征亚组中23个m6A调节因子的表达水平，并在图S5D中显示。
图4 构建m6A基因的差异表达标志和功能注释
图S5 结肠癌队列中524个m6A表型相关基因的无监督聚类
05 - m6Sig评分的构建及其临床意义
作者开发了一种被称为m6Sig评分的评分方案，该方案基于已识别的m6A相关的标志性基因，来量化CC患者个体的m6A修饰模式，用桑基图（图5A）说明了m6Sig评分构建的工作流程。具有CIN亚型的m6A基因-S2与较高的m6Sig得分有关，而m6A基因-S1表现出较低的m6Sig得分。m6A-C3表现出最高的m6Sig得分，其次是m6A-C2和m6A-C1（图S6A）。m6Sig得分与免疫激活过程和DNA修复标志明显负相关，但与EMT和基质相关标志正相关（图5B）。m6Sig评分与免疫评分之间也存在显著的反向相关性（图S6B）。与高m6Sig评分的患者亚组相比，低m6Sig评分亚组的MHC分子和免疫效应细胞的比例较高，但抑制细胞和检查点分子的比例较低（图S6C）。
图5 m6Sig评分的构建及临床特征的相关性探讨
通过将患者分为低分或高分亚组，来确定m6Sig评分在预测生存结果方面的预后能力。在CIT队列中，低m6Sig评分的患者与较好的预后显著相关（图5C），ROC曲线分析结果验证了建立的风险模型的预测优势（图S6D）。多变量Cox回归模型分析证实，m6Sig评分是评价患者结局的稳健和独立的预后生物标志物（图S6F）。探讨了m6Sig评分和分子亚型之间的关系，发现dMMR亚型比其他CC亚型的m6Sig评分更低（图5D）。调查了PD-L1的表达水平，低m6Sig评分组有明显的升高（图5E）。
通过整合TCGA-COAD数据库的临床特征和基因组信息，对构建的m6A评分系统进行了验证，发现m6Sig评分在TCGA队列中显示出潜在的预后预测价值（图S6E）。m6Sig评分在TCGA队列中显示出潜在的预后预测价值（图S6E），低m6Sig评分的患者表示有突出的生存获益（图5F）。TCGA队列的多变量分析也证实，m6Sig评分可以作为CC的独立预后生物标志物（图S6G）。TCGA分析揭示了CC的全面分子特征，并建议根据MSI状态将肿瘤细分为四个亚型。发现MSI-H亚型样本的m6Sig评分明显低于其他3个亚型样本（图5G）。
图S6 m6Sig评分与TIM免疫调节和生存结果相关
作者研究了不同m6Sig评分组的肿瘤突变负荷分布模式，发现低m6Sig评分组具有较高的突变频率（图5H）。还注意到低m6Sig评分亚组的体细胞拷贝数改变（SCNA）频率更高，这与之前SCNA与免疫逃避和肿瘤细胞增殖正相关的发现一致（图5I）。进一步对低m6Sig评分亚组与高评分亚组的CC样本进行了显著突变基因（SMG）分析。为了进一步了解CC中操作的突变过程，从TCGA-COAD突变谱中提取了三个突变信号。低m6Sig得分亚型表现出显著较高比例的突变签名6（图5J）。
06 - m6Sig评分在预测免疫治疗效果中的作用
分析发现，TIDE在低m6Sig评分组中显著降低，IPS在低m6Sig评分组中显著升高(图6A-B，图6C-D)。接下来，研究m6A修饰特征是否可以预测患者对ICI治疗的反应。在抗PD-1队列和抗CTLA-4队列中，低m6Sig评分组患者表现出显著的临床优势，生存期明显延长(图6E;图6G)。与高m6Sig评分的患者相比，低m6Sig评分的患者对ICI治疗的显著治疗效果和免疫反应得到了证实（图6H）。
图6 m6Sig评分可以预测免疫治疗的效果
四、结论
本研究基于23个m6A调控因子，全面评估了1307个结肠癌样本中的m6A修饰模式，并将这些修饰模式与TME细胞浸润特征进行了系统关联。这一综合分析表明，RNA甲基化的失调为理解肿瘤免疫的调控奠定了重要基础。更广泛地讲，评估个体肿瘤的m6A修饰模式将有助于提升对TME浸润特征的认知，并为免疫治疗疗效提供重要的启示。
网友评论