15+分的纯生信：胃癌m6A和肿瘤免疫微环境
欢迎关注「科研菌」
今天跟大家分享的是2020年3月发表在Molecular Cancer杂志上的一篇文章m6Aregulator-mediated methylation modification patterns and tumor microenvironment infiltration characterization in gastric cancer.在文章中作者利用TCGA、GEO等多个数据库和IMvigor210队列研究胃癌中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模型的建立
m6Ascore模型的建立
2.m6A修饰模式联系TME浸润特性
对应“结果解读”中图1-2（部分）
m6A修饰模式联系TME浸润特性
3.m6Asocre联系TME浸润特性
对应“结果解读”中图2-4
m6Asocre联系TME浸润特性
4.TCGA-STAD中m6Ascore对于免疫疗法的价值
对应“结果解读”中图5-6
TCGA-STAD中m6Ascore对于免疫疗法的价值
三.结果解读
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分析结果和免疫浸润情况，将三亚群表型归类：
3种亚群免疫表型
然后，使用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特征结构与m6Ascore
5.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的免疫疗法有十分重要的价值。
还是和往常一样，后台回复「32b」，即可获取今天小编为大家解读的文献和『详细版』思路分析图。本期的分享就到这里啦，一起期待下一期『泛癌m6A与5mC表观遗传crosstalk』吧~
往期热门文章推荐
5. 生信杂志1|被“灌水”的Aging(IF=5.515)还能不能投？
6. 生信杂志2|推荐一本生信友好并且最快1月接受的5+分杂志
8. 生信杂志4|即将突破5分的Frontiers in Oncology是灌水杂志吗？
9. 生信杂志5|博士毕业神器3+分纯生信杂志：平均一审只要一个月，年刊量1000+
10. 生信杂志6|马上4分年刊近2000的杂志JCB：对纯生信友好还不要版面费！
本期的文献分享就到这里啦，一起关注「科研菌」，一起期待下期的精彩内容吧
扫上方二维码即可关注科研菌
