m6A评分可作为胃癌的预后标志物并预测免疫治疗的效果，IF=10.679，实时IF=15.83

胃癌中m6A调控因子介导的甲基化修饰模式和肿瘤微环境浸润特征

m6A regulator-mediated methylation modification patterns and tumor microenvironment infiltration characterization in gastric cancer

文献精读汇报人：陈海影、徐欣（排名不分先后）

背景知识介绍

1.N6甲基腺苷（m6A）是一种RNA修饰，是哺乳动物细胞中一种动态可逆的过程，受甲基转移酶、去甲基化酶和结合蛋白的调节。深入了解m6A调节因子有助于揭示m6A修饰在转录后调控中的作用和机制，免疫应答的表观遗传调节作用已被证实，越来越多的证据也表明，m6A调节因子的表达失调和基因改变与多种生物学过程的紊乱有关，包括细胞的死亡和增殖失调、发育缺陷、肿瘤恶性进展、自我更新能力受损和免疫调节异常。

2.肿瘤部分由一个复杂的肿瘤微环境（TME）组成，新的证据揭示了其在肿瘤发展、免疫逃逸和免疫治疗反应中的重要作用，根据TME细胞浸润特征预测ICB的反应是提高现有ICB成功率和开发新的免疫治疗策略的关键步骤。

3.研究发现TME浸润性免疫细胞与m6A修饰之间存在着特殊的相关性。本研究中，作者以21种m6A调控因子为基础，对GEO、TCGA中的胃癌样本的m6A修饰模式进行综合评价，并将这些修饰模式与TME细胞浸润特性系统地关联起来。利用主成分分析法构建m6A评分，量化每个胃癌患者的m6A修饰模式。最终确定了3种不同的m6A修饰模式，三种模式下的TME细胞浸润特征与包括免疫排斥、免疫炎症和免疫荒漠在内的三种肿瘤免疫表型高度一致，表明M6A修饰在塑造个体肿瘤微环境特征方面发挥了不可忽视的作用。在单个肿瘤中评估m6A修饰模式可以预测肿瘤炎症的分期、亚型、TME间质活性、遗传变异和患者预后，m6A评分较低的患者显示出显著的治疗优势和临床效益。

  文献结果解读

NO.1 胃癌中m6A调节因子的遗传和表达变化

  图（1)

首先作者总结了调节因子介导的m6A RNA甲基化的动态可逆过程及其对RNA的潜在生物学功能。从图中可以看出m6A RNA在RBM15,WTAP,METTL14,METTL3组成的甲基转移酶介导下发生甲基化，在FTO,ALKBH5,ALKBH3脱甲基酶介导下发生去甲基化。YTHDFs、YTHDCs组成的一组特异性RNA结合蛋白可以识别m6A基序，从而发生RNA衰变或RNA翻译。

      图（2)

然后作者对TCGA-STAD队列中433例胃癌患者的21个m6A调控因子在GC(胃癌)中的体细胞突变的发生率进行了总结。在433例样本中，有101例发生了m6A调节因子突变，突变频率为23.33%。由图可见，ZC3H13的突变频率最高，其次是KIAA1429，而去甲基化酶（FTO和ALKBH5）和METTL3在GC样品中均未显示任何突变。进一步分析显示，ELAVL1与KIAA1429、YTHDF1和ZC3H13以及RBM15和YTHDC1之间存在显著的突变共现关系。

图（3）

接着作者对GSE62717队列中m6A调节子的CNV变异频率进行了总结。图C显示21个调节器中普遍存在CNV改变，主要集中在拷贝数扩增，而ELAVL1、YTHDF2和FMR1则具有广泛的CNV缺失频率。柱子的高度表示突变频率，缺失频率为蓝点；扩增频率为红点。图D显示了m6A调节因子在23对染色体上的CNV改变位置（采用RCircos R包绘制），同时基于这21个m6A调节因子的表达，可以完全区分GC样品和正常样品（如图E所示）。

图（4）

为了确定上述基因变异是否影响m6A调节因子在GC患者中的表达，作者研究了正常人和GC患者m6A调节因子的mRNA表达水平，发现CNV的改变可能是m6A调节因子表达受干扰的重要因素。从图C,F可知，与正常胃组织相比，CNV扩增的m6A调节因子在GC组织（如CBLL1和FTO）中的表达明显更高，反之亦然（如ELAVL1和YTHDF2），提示m6A调节因子的表达失衡在胃癌的发生发展过程中起着至关重要的作用。

NO.2 m6A甲基化修饰模式及各模式的生物学特性

图（5）

作者研究了m6A调节因子的相互作用，并绘制了一幅图（如图所示）展示调节因子相互作用的关联。从图中可以看出，作者研究了21种m6A调节因子，并将其分为4个调节器簇，研究了m6A调节器相互作用、相互连接及其对GC患者预后意义的综合情况。

         图（6)

然后利用5个具有OS信息的数据集GSE15459、GSE34942、GSE57303、GSE62254/ACRG和GSE84437，作者根据21个m6A调节因子的表达情况，使用R包ConsensusClusterPlus对具有不同性质的m6A修饰模式的患者进行分类，最终通过无监督聚类识别出三种不同的修饰模式，包括m6Acluster-A 中的389例、m6Acluster-B 中的348例和m6Acluster-C 中的322例。基于三种修饰模式进行生存分析，由生存曲线可知，m6Acluster-B的患者显示出显著的生存优势。

   图（7)

接着作者进行了GSVA富集分析，由热图可见，m6Acluster-A在基质和致癌激活途径中显著富集，m6Acluster-B呈现与免疫充分激活相关的富集途径，m6Acluster-C与免疫抑制的生物学过程密切相关。

NO.3 不同m6A修饰模式下TME细胞浸润特征和转录组特征

   图（8）

作者研究了三种m6A修饰方式下TME浸润细胞的丰度，图中显示m6Acluster-A在NK细胞、巨噬细胞、嗜酸性粒细胞、肥大细胞、MDSC、浆细胞样树突状细胞等天然免疫细胞浸润方面非常丰富。

然而，m6Acluster-A修饰模式的患者并没有显示出与丰富免疫细胞浸润相匹配的生存优势。以往的研究表明，在免疫排斥表型的肿瘤也存在丰富的免疫细胞，而这些免疫细胞停留在肿瘤细胞巢周围的基质中，而不是穿透肿瘤细胞的实质。TME间质的激活被认为是T细胞抑制。GSVA分析的结果表明，m6Acluster-A修饰模式与基质激活显著相关。

因此，我们推测m6Acluster-A的基质激活抑制了免疫细胞的抗肿瘤作用。随后的分析表明（如图B所示），m6Acluster-A的间质活性显著增强，如上皮-间充质转化(EMT)、转化生长因子β(TGFb)和血管生成途径的激活，这证实了作者的推测。基于以上分析，发现三种M6A修饰模式具有明显不同的TME细胞渗透特性。m6Acluster-A为免疫排斥表型，以天然免疫细胞浸润和基质激活为特征；m6Acluster-B群为免疫炎症表型，特征为适应性免疫细胞浸润和免疫激活；m6Acluster-C群为免疫荒漠表型，特征为免疫抑制。

然后，作者使用CIBERSORT算法来确定肿瘤中的免疫细胞类型，以比较三种M6A修饰模式中免疫细胞组成的差异,发现三种M6A修饰模式在TME细胞类型组成上没有显著差异（文章补充文件 Figure S2K），提示M6A甲基化修饰不改变肿瘤的TME浸润细胞类型。

  图（9）

作者研究了三种不同m6A修饰模式的间质激活途径的差异：包括EMT、TGF-β和血管生成途径等(如图B所示)，用单因素方差分析检验三种修饰方式的差异，从图中可见，A簇中基质活性显著增强，而在CD8 T effector中相对减弱，证实了A簇基质激活抑制了免疫细胞的抗肿瘤作用。以上是针对5个GEO队列合并(用SVA包的 ComBat 算法去批次效应的)后的分析结果。

接着作者又将注意力集中在GSE62254(ACRG)队列，这次同样是基于21个m6A调节因子的表达情况对GSE62254(ACRG) 队列进行无监督聚类，与上面的5个GC数据集的聚类结果相类似，在GSE62254(ACRG) 队列也发现了三种完全不同的m6A修饰模式。主成分分析结果显示该队列这三种不同的m6A修饰模式之间，存在显著的m6A转录谱差异（如图D所示）。如图E所示，作者描绘了GSE62254(ACRG) 队列中三种不同m6A修饰模式间分子亚型的比例。

  图（10）

图C，为利用这21个m6A调节因子对GSE62254(ACRG)队列进行无监督聚类的热图，作者还将肿瘤分期、分子亚型、生存状态和年龄一并画了出来。其中黄色表示调控因子高表达，蓝色表示低表达。

在上面将GSE62254(ACRG) 队列患者分为三种不同的m6A修饰模式后，接着作者用limma R包来进行差异分析，以确定GSE62254(ACRG) 队列三种模式之间的差异基因（DEGs），adjusted P value < 0.001 的基因被认定为是DEGs。最终确定了718个m6A表型相关DEGs(见文献的补充图S4G和表S7)。并使用Y叔的clusterProfiler 包对这些DEGs进行GO富集分析，文献补充表S8总结了显著富集的生物过程(BP，biological processes )，图3F显示了这些基因富集于m6A修饰和免疫显著相关的生物学过程，再次证实了m6A修饰在肿瘤微环境的免疫调节中起着不可忽视的作用。

NO.4 m6A基因标签的构建

 图（11）

为了进一步验证m6A修饰在肿瘤微环境的免疫调控机制，作者这次根据这718个m6A表型相关基因（DEGs）对GSE62254(ACRG) 队列进行了无监督聚类，将患者分为不同的基因组亚型，这里使用了一致性聚类算法（consensus clustering algorithm ）来确定基因簇的分类个数及其稳定性。

如图A所示，与m6A修饰模式的聚类分组相似，基于718个差异基因的无监督聚类也将GSE62254(ACRG) 队列患者分为三种不同的m6A修饰基因组亚型，分别命名为m6A基因簇A-C。这也表明胃癌中确实存在三种不同的m6A甲基化修饰模式。

由图可见，三个不同的基因簇具有不同的特征基因。M6A基因簇C模式的胃癌患者表现出较差的分化，并且在弥漫性组织学亚型中富集，在m6A基因簇A和基因簇B中观察到相反的模式。Alive或MSI亚型的患者主要集中在M6A基因簇A，而临床IV期或EMT分子亚型的患者主要表现为M6A基因簇C模式。

图（12）

作者还对GSE62254(ACRG)队列3个基因簇的患者进行了生存分析，由生存曲线可见，基因簇C的患者预后最差，而基因簇A的患者预后最好，这个结果可能与二者不同的基因富集结果有关。

图（13）

作者采用单因素方差分析（one-way ANOVA test）检验GSE62254(ACRG) 队列三个基因簇中21个m6A调节因子表达的统计学差异。如图C所示，21个m6A调节因子在三个基因簇的表达存在显著差异，这与m6A甲基化修饰模式的预期结果一致。

（m6Ascore计算公式，i是m6A表型相关基因的表达量）

PS: 这里计算m6Ascore实在没想明白，查阅了一下资料，目前的猜想是：假设基因A表达量是10，其在主成分1的系数是1，在主成分2的系数是2。基因B的表达量是20，其在主成分1的系数是3，在主成分2的系数是4。那分数等于=（10×1+10×2）+（20×3+20×4）？？不确定这个猜想对不对，欢迎大家一起来研究、探讨，谢谢！！

M6A基因标签的产生：为了量化每个胃癌患者的M6A修饰模式，作者构建了一套评分系统来评估单个胃癌患者的M6A修饰模式-M6A基因标签，称之为m6Ascore。

建立m6A基因标签的步骤如下：

首先对GSE62254(ACRG) 队列所有样本进行归一化，然后两两差异分析，提取三种m6A修饰模式的718个交集的DEGs。接着对这718个DEGs进行单因素Cox回归分析(补充表S7提供了这718个DEGs的单因素cox分析结果)，进一步提取有预后意义（单因素cox分析P<0.05）的DEGs做主成分分析(PCA)来构建m6A相关基因标签，主成分1和主成分2都被选作构建签名分数。

    图（14)

图D的三冲积图显示了GSE62254(ACRG)队列不同的m6A修饰模式簇、分子亚型、基因簇和m6Ascore(高、低组)之间的关联。几乎所有EMT亚型(ACRG队列中的分子亚型)患者(41/45，91%)都被归入基因簇C，这与较差的生存结果有关。

    图（15）

为了研究m6A基因标签与一些相关的生物学途径之间的联系，作者用Spearman相关性分析，研究了GSE62254(ACRG)队列中m6Ascore与已知基因特征(补充表S9，来自Mariathasan等人的研究)之间的关系，蓝色为负相关，橙色为正相关。

图（16）

接着，作者用Kruskal-Wallis检验比较了GSE62254(ACRG)队列中三个基因簇之间m6Ascore的差异（P<0.001）。如F图所示，基因簇A的中位数得分最低，而基因簇C的中位数得分最高，这表明低m6Ascore可能与免疫激活（基因簇A）相关信号密切相关，而高m6Ascore可能与基质激活(基因簇C)相关信号密切相关。更重要的是，与其他m6A修饰模式簇相比（P< 0.001, Kruskal-Wallis test)），m6Acluster-A显示出显著的m6Ascore增加，而m6Acluster-B表现出最低的中位数得分（图4G）。

图（17）

进一步，对活化基质相关通路的分析表明，如图H中显示（GSE62254队列），高m6Ascore组和低m6Ascore组的间质激活途径存在差异，高m6Ascore与基质通路的激活显著相关。

NO.5 TCGA分子亚型与肿瘤体细胞突变的m6A修饰特征

                           图（18）

作者还研究了GSE62254(ACRG)队列不同分子亚型间m6Ascore的差异，用Kruskal-Wallis检验比较4种分子亚型的差异（p<0.0001）。接着，为进一步确定m6Ascore在预测患者预后方面的价值，根据Survminer R包确定的 cutoff 值 0.0291，将患者分为高、低m6Ascore组。采用Kaplan-Meier绘制生存曲线，对ACRG队列中低m6Ascore（157例）和高（143例）m6Ascore患者进行生存分析，如图B所示，低m6Ascore的患者显示出显著的生存优势(HR 3.0(2.12-4.21)，其5年生存率是高m6Ascore患者的两倍(69.4%比33.5%)。此外，包括年龄、性别、TNM状态、组织类型、MSI状态、TP53状态和分子亚型在内的多因素cox回归分析，证实了m6Ascore是评估患者预后的可靠和独立的预后生物标志物(HR2.54(1.71-3.8)，文献补充图S6A)。

作者还检测了m6Ascore预测胃癌患者接受辅助化疗疗效的能力，发现在同时接受辅助化疗的患者中，低m6Ascore组的患者显示出显著的治疗优势，5年生存率为77.5%对59.2%(图C)。而且辅助化疗不干扰m6Ascore的预测能力：因为无论是否接受化疗，低m6Ascore组总是显示出明显的生存优势(图C)。

此外，作者还发现年轻的患者、弥漫性组织学亚型和晚期患者与较高的m6Ascore显著相关，这意味着这些患者具有m6Acluster-A修饰模式和免疫排斥表型的特征，临床预后较差。这些结果表明，m6Ascore可以用来评估患者的某些临床特征，如MSI状态、分子亚型、组织学亚型以及临床分期等(文献补充图S6B)。

                                  图（19）

接着，作者又在TCGA-STAD队列中验证了m6A修饰模式的作用和特征。同样发现低m6Ascore组患者具有更显著的生存优势（如图D)。且不同分子亚型间m6Ascore的同样具有显著差异(如图E，Kruskal-Wallis试验，p<0.0001)。不同微卫星亚型间的m6Ascore的也具有差异(如图F）。较高的m6Ascore患者明显集中在GS亚型，患者的生存率较低，而较低的m6Ascore集中在MSI和EBV感染的亚型，这与较好的生存率有关(5年生存率，25.9%比43.3%；HR1.81(1.26-2.62)；图5D-E)。预后较好的MSI-H亚型与较低的m6Ascore相关，而MSI-L和MSS亚型的m6Ascore较高(图5F)。此外，对TCGA-STAD队列的多因素cox分析也证实了m6Ascore可以作为胃癌独立的预后生物标志物(文献补充图S6C)。

                              图（20）

作者使用maftools R包对TCGA-STAD队列高m6Ascore（图G）和低m6Ascore(图H)患者建立的肿瘤体细胞突变瀑布图中，每一列代表一个患者，上面的柱状图显示TMB，右边的数字表示对应的基因的突变频率，右边的条形图显示了每种突变类型的比例，通过图G、H可见低m6Ascore组患者比高m6Ascore组表现出更广泛的TMB。TMB定量分析证实m6Ascore低的肿瘤与高TMB显著相关(文献补充图S7E，Wilcoxon，p<2.2e−16)。m6Acore和TMB的相关性分析显示两者呈显著负相关(文献补充图S7F，Cor=−0.561，p=1.33e−29)。

越来越多的证据表明，高TMB状态的患者对抗PD-1/PD-L1免疫治疗有持久的临床反应。因此，上述结果间接证明肿瘤m6A修饰模式的差异可能是介导抗PD-1/PD-L1免疫治疗临床应答的重要因素。并间接证实了m6Ascore在预测免疫治疗结果中的价值。

NO.6  m6A修饰模式在抗PD-1/L1免疫治疗中的作用

                                   图（21）

为了进一步探索m6A修饰模式在抗PD-1/L1免疫治疗中的作用，作者又下载了数据集IMvigor210（晚期尿路上皮癌队列，抗PD-L1抗体干预）和GSE78220（转移性黑色素瘤队列，抗PD-1抗体（pembrolizumab）干预），分别分析证实m6Ascore对抗PD-L1、PD-1治疗效果的预测价值。抗PD-1/L1免疫治疗的患者，其生存曲线分别如图A和D显示，低m6Ascore的患者有着更好的预后和临床结局。图B和图E，在低m6Ascore中，获得治疗完全缓解或部分缓解的患者所占的比例比高m6Ascore组更大。

                                图（22）

由图6B,C,E,F,G可知，与高m6Ascore的患者相比，低m6Ascore的患者接受抗PD-1/L1免疫治疗，具有显著的治疗优势和临床反应。

                                     图（23）

此外，作者还研究了高、低m6Ascore组PD-L1的表达差异，如图H所示，低m6Ascore组的患者PD-L1的表达更高，这表明这类患者对抗PD-1/L1免疫治疗有潜在的反应。

                                  图（24）

进一步的研究表明，在抗PD-L1免疫治疗队列（IMvigor210）中，低m6Ascore组和高m6Ascore组间的基质激活途径和调节性T细胞丰度(被认为是免疫抑制)的存在显著差异。在高m6Ascore的患者中，调节性T细胞和TME基质被显著激活，这介导了肿瘤的免疫耐受(图6I)。

                                      图（25）

作者还评估了与免疫治疗效果密切相关的肿瘤新抗原负荷。如图J，采用Kaplan-Meier曲线对接受抗PD-L1（IMvigor210）免疫治疗的患者按m6Ascore和新抗原负荷进行分层生存分析，发现低m6Ascore和高新抗原负荷的患者显示出巨大的生存优势(P<0.0001，Log-rank test），H为高抗原负荷，L为低抗原负荷，NEO为新抗原负荷

                                   图（26）

图K，m6A修饰模式的定量（即m6Ascore）在接受抗PD-1/L1免疫治疗的患者中的预测价值，AUC = 0.768。表明m6Ascore是评估接受免疫治疗的患者预后和临床反应一个潜在的、强有力的生物标记物。

                                图（27）

此外，作者还研究了IMLIGN210队列中三种肿瘤免疫表型之间m6Ascore的差异（(p=0.015，Kruskal-Wallis检验），发现高m6Ascore者主要与免疫排斥和免疫荒漠表型显著相关，免疫检查点抑制剂在这些表型中很难发挥抗肿瘤作用(图6L)。

综上所述，m6A甲基化修饰模式与肿瘤免疫表型和对抗PD-1/L1免疫治疗的应答显著相关，建立m6A修饰标签将有助于预测对抗PD-1/L1免疫治疗的应答。

                    全文论证思路图

结论

本研究中揭示了m6A甲基化修饰对肿瘤微环境的广泛调节机制。m6A修饰模式的差异是导致个体肿瘤微环境异质性和复杂性的不可忽视的因素。对个体肿瘤m6A修饰模式的综合评价将有助于加深对肿瘤微环境细胞浸润特征的认识，指导更有效的免疫治疗策略。

特别鸣谢：在做文献汇报的过程中，学习和模仿了珠江肿瘤罗鹏师兄他们公众号作文献汇报的风格，非常感谢罗鹏师兄及其团队的伙伴们！！很开心可以邀请到王亮教授作为我们文献精读系列活动的指导老师，谢谢老师的支持和鼓励！！嘻嘻！！

参考文献：Bo Zhang1,2†, Qiong Wu1,2†, Ben Li3, Defeng Wang1, Lei Wang4and You Lang Zhou1. m6A regulator-mediated methylation modification patterns and tumor microenvironment infiltrationcharacterization in gastric cancer

                                                                                      编辑: 陈海影；徐欣

                                                                                                                                 初校审：梁晓杰

                                                                                                                                    校审专家：王亮  

校审专家介绍

王亮，医学博士，副主任医师，副教授，硕士研究生导师，现任北京同仁医院血液科副主任。从事血液系统肿瘤诊疗工作10余年，对各类淋巴瘤、骨髓瘤的诊断治疗具有丰富的临床经验。主要研究方向为淋巴瘤的靶向免疫治疗及耐药机制探索，主持多项国家自然科学基金，发表专业论著三十余篇，SCI收录二十余篇，发表在The Lancet Haematology, Journal of Hematology and Oncology, Clinical Cancer Research, Cancer, British Journal of Cancer等国际著名期刊。曾荣获“敬佑生命·2017荣耀医者-青年创新奖”、首批“广东省杰出青年医学人才”。学术兼职：中华医学会血液学分会青年委员；国家自然科学基金评审专家。担任《Aging Pathobiology and Therapeutics》副主编、《中国癌症防治杂志》编委、《协和医学杂志》青年编委。

相关阅读：