15+分的纯生信：胃癌m6A和肿瘤免疫微环境
虽然MC目前2018年的IF是10.679，但是我们之所以说是15+分呢，是因为MC在6月就要公布在2019年的影响因子是15.96分，说不定还能在这段时间努努力，突破到16分！
(图片来源：桑格助手)m6A regulator-mediated methylation modification patterns and tumor microenvironment infiltration characterization in gastric cancer胃癌中m6A调节因子介导的甲基化修饰模式与肿瘤微环境浸润特征的研究一.研究背景越来越多的证据表明，m6A修饰通过与多种m6A调节因子的相互作用，在炎症、先天免疫、抗肿瘤等方面发挥着不可或缺的作用。由于大多数研究集中在单个TME细胞类型或单个调控因子上，因此由多个m6A调控因子综合作用介导的TME整体浸润特征并未得到全面认识。明确不同m6A修饰模式在TME细胞浸润中的作用，有助于加深对TME抗肿瘤免疫反应的认识，指导更有效的免疫治疗策略。二.分析流程
1.m6Ascore模型的建立
2.m6A修饰模式联系TME浸润特性对应“结果解读”中图1-2（部分）
3.m6Asocre联系TME浸润特性对应“结果解读”中图2-4
4.TCGA-STAD中m6Ascore对于免疫疗法的价值对应“结果解读”中图5-6
三.结果解读1.胃癌m6A调节因子的遗传变异情况首先先介绍下作者对胃癌数据集和预处理：本研究从GEO、TCGA共收集了7个符合条件的GC组。符合条件的GC数据集的基本信息汇总在表S1中。对于来自Affymetrix®的微阵列数据，下载原始CEL文件，采用多阵列平均方法，使用affy和simpleaffy包作背景调整和分位数标准化。对于TCGA数据集基因表达的RNAsequencing 数据(FPKM值)从GDC下载，使用R包TCGAbiolinks 进行整合分析。FPKM值转换为转录本每千碱基百万(TPM)值。利用R sva包对非生物技术偏差的bacth effect（批量效应)进行了校正。
表S1.数据集的基本信息本研究作者共确定了21个m6A调节因子，包括8个甲基化酶、2个去甲基化酶和11个结合蛋白。作者总结展示了调控因子介导的m6A RNA甲基化的动态可逆过程及其潜在的RNA生物学功能（图1A）。然后总结了GC中21个m6A调节因子拷贝数变异（CNV）和体细胞突变的发生率：在TCGA-STAD队列433个样本中，有101个m6A调节因子发生突变，频率为23.33%。发现ZC3H13突变频率最高，KIAA1429次之，而去甲基化酶(FTO和ALKBH5)和METTL3在GC样品中均未发生突变(图1B)。进一步分析发现，ELAVL1与KIAA1429、YTHDF1和ZC3H13、RBM15和YTHDC1之间存在显著的突变共现关系(图S1B)
图S1B.21个m6A调节因子的共现分析对GSE62717队列中CNV改变频率的调查显示，21个调控因子中CNV的改变较为普遍，主要集中在拷贝数的扩增上，而ELAVL1、YTHDF2和FMRI中CNV的缺失较为普遍(图1C)。m6A调节子在染色体上CNV改变的位置如图1D所示。根据21m6A调节因子的表达，区分出GC样品和正常样品(图1E)。为了确定上述基因变异是否影响胃癌患者m6A调节因子的表达，作者调查了正常与GC样品间调节因子mRNA表达水平，发现CNV的改变可能是干扰m6A调节因子表达的重要因素。与正常胃组织相比,在GC组织m6A调节因子与扩增 CNV显著高表达(例如CBLL1和FTO),反之亦然(例如ELAVL1和YTHDF2)(图1 C、F)。
图1.胃癌中m6A调节因子的遗传和表达变化2.m6A甲基化修饰模式由21个调控因子介导作者将 5个GEO数据集纳入一个meta队列。作单变量Cox回归模型，显示21个m6A调节因子在胃癌患者中的预后价值（图S1C）。
图S1C.21个m6A调节因子在胃癌患者中的预后价值图2A作了m6A调节因子网络分析：描述了m6A调节因子相互作用、调节因子连接及其对胃癌患者预后的影响。图S2A-H：作者writer、eraser基因表达的关联显著。而且发现writer基因高表伴随的eraser基因低表达的现象实际上取决于不同的writer基因和eraser基因：writer基因(WATP和RBM15)高表达的肿瘤，其eraser基因FTO的表达较低，而WATP和RBM15的高表达则不影响另一个eraser基因ALKBH5的表达（图S2A-B）。图S2l:与野生型相比，ZC3H13突变型肿瘤中ALKBH5表达明显上调，而FTO表达明显下调。
图S2.writer/eraser基因表达的关系（图中ALKBH5、FTO为eraser基因）将GEO五个队列1051个患者样本作一致性聚类分析，分为了三个亚群：A（n=389）、B（n=348）和C（n=322）。然后对三种主要m6A修饰亚型的预后分析显示，m6A亚群B修饰模式具有更好的OS(图2B)。为了探索这些不同的m6A修饰模式之间的生物学行为，作者进行了GSVA富集分析（图2C)。m6A亚群A显著富集了ECM受体相互作用、TGF -信号通路、细胞粘附、MAPK信号通路等间质激活致癌通路。m6A亚群B表达了与免疫完全激活相关的富集通路，包括趋化因子信号通路的激活、细胞因子-细胞因子受体相互作用、T细胞受体信号通路和Toll样受体信号通路(图2C)。而m6A亚群C与免疫抑制生物学过程显著相关。(图2D)
图2.m6A的甲基化修饰模式及各模式的生物学特性3.TME细胞浸润特征和转录组特征图3A：接着对TME细胞浸润的分析显示，亚群A明显富含天然免疫细胞浸润，包括自然杀伤细胞、巨噬细胞、嗜酸性粒细胞、肥大细胞、MDSC、浆细胞样树突状细胞。图3B：后续分析显示,亚群A的基质显著增强如激活EMT、转化生长因子β(TGFb)和血管生成途径的激活,证实了作者基于GSVA分析结果:A基质的激活可能抑制了免疫细胞的抗肿瘤作用。作者依据三个亚群的GSVA分析结果和免疫浸润情况，将三亚群表型归类：Cluster免疫表型特征A免疫排斥表型天然免疫细胞浸润和基质活化B免疫排斥表型适应性免疫细胞浸润和免疫激活C免疫荒漠表型免疫抑制然后，使用spearman的相关分析来检验每个TME浸润细胞类型和每个m6A调节因子之间的具体相关性(图S3A)。重点研究了调节因子KIAA1429：图S3B：ESTIMATE算法量化KIAA1429，KIAA1429低表达时免疫评分较高。图S4C：23种TME浸润免疫细胞在KIAA1429高、低表达患者中的特异性差异分析我们发现KIAA1429低表达的肿瘤在23个TME免疫细胞中浸润明显增加(图S3C)。图S3D：KIAA1429的表达降低导致MHC分子、共刺激分子和粘附分子的表达全面升高。图S3E：KIAA1429的通路富集分析。图S3F：在抗PD-L1免疫治疗队列中，KIAA1429低表达患者具有更好的OS。综上作者推测KIAA1429介导的m6A甲基化修饰可能促进TME树突细胞的激活，从而增强瘤内抗肿瘤免疫应答。
图S3.不同m6A修饰模式的TME细胞浸润特征和转录组特征图3D：作者进一步用ACRG队列探讨分析m6A修饰表型在不同临床特征和生物学行为中的特征（n=300）。ACRG队列的聚类分析显示了三种完全不同的m6A修饰模式(图3c-d)。三种不同的m6A修饰模式在m6A转录谱上存在显著差异。图3E：m6Acluster-A甲基化修饰模式以EMT分子亚型为特征，m6Acluster-B甲基化修饰模式以MSI亚型为特征。作者接着使用limma包测定了718个m6A表型相关的DEGs。然后用clusterProfiler包对DEGs作GO富集分析。结果显示这些基因表现出丰富的与m6A修饰和免疫显著相关通路，证实m6A修饰在肿瘤微环境的免疫调节中发挥着不可忽视的作用(图3F)。
图3.不同m6A修饰模式的TME细胞浸润特征和转录组特征4.m6A基因标记和功能注释的产生根据获得的718个m6A表型相关基因进行无监督聚类分析，将患者分为不同的基因组亚型。与m6A修饰模式的聚类结果一致，聚类也发现了三个具有不同的特征基因不同的m6A基因簇（图4A)。接着对三个基因簇的患者作了KM生存分析，显示m6A基因簇A具有更好的OS，其次是基因簇B以及基因簇C(图4B)。在三个m6A基因簇中，m6A调控因子的表达存在显著差异(图4C)。作者考虑到m6A修饰的个体异质性和复杂性，基于这些表型相关基因，构建了一套评分系统来量化胃癌个体患者的m6A修饰模式，称之为m6Ascore。用冲积图显示了m6Aclusters、ACRG分子亚型、基因簇和m6Ascore的变化(图4D)。为了更好地说明m6A特征，作者利用Spearman分析m6Ascore与ACRG队列中已知基因特征的相关性。负相关用蓝色表示，正相关用橙色表示(图4E)。作Kruskal-Wallis检测显示m6A基因簇在m6Ascore上存在显著差异。基因簇A的中位得分最低，而基因簇C的中位得分最高，说明m6Ascore含量低可能与免疫激活相关特征密切相关，而m6Ascore含量高可能与基质激活相关特征密切相关(图4F)。而对3个m6A亚群作Kruskal-Wallis，显示m6A亚群A较B、C亚群m6Ascore明显增加(图4G)。对基质相关通路活性的分析表明，高m6Ascore可能对基质激活通路有影响(图4H)。
图4.m6A特征结构与m6Ascore5.m6Ascore与TME表型特征的联系m6Ascore与ACRG亚型：图5A:采用Kruskal-Wallis检验比较四种ACRG分子亚型间的统计学差异，与其他三种ACRG分子亚型相比，EMT亚型患者的m6Ascore最低。图5B:作者进一步探究m6Ascore在预测患者预后方面的价值。将患者分为低m6Ascore组和高m6Ascore组后，作KM生存分析，低m6Ascore患者具有更好的OS。图5C：然后检测了m6Ascore预测胃癌患者辅助化疗疗效的能力。发现m6Ascore低且同时接受辅助化疗的患者具有更好的OS。图S6A：多变量Cox回归模型，证实m6Ascore是独立预后因素。图S6B：m6Ascore与临床特征的关联。
图S6.m6Ascore的预后价值及临床特征作者在TCGA-STAD中作了验证:TCGA-STAD将胃癌分为基因组稳定型(GS)、微卫星不稳定型(MSI)、EBV感染和染色体不稳定型(CIN) 4个分子亚型。作者评估了m6Ascore在这些分子亚型之间的差异:显示高m6Ascore明显集中于GS亚型，患者OS差；低m6Ascore则集中于MSI、EBV感染亚型，OS较好(图5D-E)。接着探究不同微卫星亚型间m6Ascore差异。微卫星高不稳定性亚型（MSI-H）的m6Ascore值较低，而MSI-L和MSS的m6Ascore值较高(图5F)。然后，利用maftools软件包分析了TCGA-STAD队列中低、高m6Ascore体细胞突变的分布差异。如图5G-H所示，低m6Ascore组（图5H）较高m6Ascore组（图5G）有更广泛的肿瘤突变负荷。根据高TMB状态的患者对抗pd -1/PD-L1免疫治疗有持久的临床反应。证实m6Ascore在预测免疫治疗结果中的价值。
图5.TCGA分子亚型m6A修饰及肿瘤体细胞突变的特点 微卫星状态和EBV感染：图S7A：从三种m6A修饰模式在MSI和MSS亚型中的比例进一步发现：cluster-B主要以MSI亚型为特征，而cluster-C以MSS亚型为特征。图S7B：MSS/MSI亚型中m6A调节因子的表达情况。图S7C：EB病毒感染在三个cluster的比例。图S7D：EBV阳性/阴性m6A调节因子的表达情况。图S7E：TMB定量分析：低m6Ascore呈现高TMB。图S7F：相关性分析，m6Ascore与TMB也呈显著负相关。
图S7.微卫星状态和EB病毒感染对三种m6A修饰模式和m6A调节因子的影响6.m6A甲基化修饰模式与肿瘤免疫表型及抗pd -1/L1免疫治疗反应显著相关作者将ACRG队列中建立的m6Ascore：图S8A-S8C：应用于其他独立的胃癌队列中，验证其预后价值。图S8D：对所有GEO队列的组合进行了验证。图S8E-S8F：评估m6Ascore预测RFS的能力。图S8G：继续扩展到TCGA的所有消化系统肿瘤。这些数据表明m6A修饰模式与更好的临床疗效相关。图S8H-S8I：ROC曲线显示m6A预测优势在老年患者中尤为明显。
图S8.m6Ascore对胃癌组和消化系统癌症组的预后价值作者又基于两个免疫治疗队列:抗PD-L1队列(IMvigor210)和抗PD-1队列(GSE78220)，研究m6A修饰签名是否可以预测患者对免疫检查点阻断治疗的反应。对两个队列作KM生存分析，两队列（图6A:IMvigor210；图6D:GSE78220）低m6Ascore患者均具有较好的OS。通过对两队列中低、高m6Ascore组患者对阻断PD-L1(图6b)和PD-1（图6E）免疫治疗反应比例分析，显示低m6Ascore比高m6Ascore患者抗PD-1/L1免疫治疗具有显著临床疗效优势(图6B-C，E-F，G-H)。对低m6Ascore组和高m6Ascore组间基质激活通路和调节性T细胞丰度分析，显示高m6Ascore组的调节性T细胞丰度显著提升，以及TME基质通路显著激活(图6I)。又将m6Ascore合并肿瘤新生抗原负荷作KM分析，显示低m6Ascore、高肿瘤新生抗原负荷的患者具有更好的OS(图6J)。对m6Ascore在抗pd -1/L1免疫治疗患者中的预后价值分析：ROC曲线（AUC=0.768）显示m6Ascore是一种可靠的免疫治疗预后和临床反应评估的生物标志物(图6K)。最后利用IMvigor210队列中已知的肿瘤免疫表型，来对m6Ascore在不同表型之间差异进行分析：显示免疫排斥表型、免疫沙漠表型的m6Ascore较高(图6L)。
图6.m6A修饰模式在抗pd -1/L1免疫治疗中的作用最后我们总结一下：作者以胃癌为例，在1938个胃癌样本mRNA表达谱的基础上，充分评估了21个m6A调节因子的状态，并将胃癌患者分为三组m6A甲基化修饰模式，通过对三组修饰模式的差异分析，发现分别与肿瘤的三种免疫表型高度对应，由此建立一套m6Ascore的评分系统。然后证实m6Ascore与胃癌的m6A修饰特征、微环境炎症水平、基质活性、分子亚型、遗传变异、微卫星状态以及患者预后显著相关。m6A对于抗PD-1/L1的免疫疗法有十分重要的价值。