最新m6A的7+分新思路！你也能学会！
Comprehensive Analysis of the PD-L1 and Immune Infifiltrates of m6 A RNA Methylation Regulators in Head and Neck Squamous Cell CarcinomaPD-L1和头颈部鳞状细胞癌m6A甲基化调节因子的免疫浸润的综合分析
一、研究背景
头颈部鳞状细胞癌(HNSCC)是世界上第六大最常见的恶性肿瘤。免疫系统的失衡在以免疫抑制为特征的HNSCC中起重要作用。
腺嘌呤第六N原子处的甲基化修饰（N6-甲基腺苷(m6A)）是mRNA上最常见的转录后修饰。m6A水平的异常与干细胞分化、免疫应答、胚胎发育、miRNA编辑有关，m6A的甲基化通过调节癌基因或抑癌基因mRNA表达，影响肿瘤的进展。m6A的甲基化水平主要取决于m6A甲基化调节因子的水平，m6A甲基化调节因子包括细胞内甲基转移酶、去甲基化酶、与甲基化位点结合的蛋白质，它们的异常可导致肿瘤的生长、迁移等。m6A甲基化调节因子可作为HNSCC患者的预后预测物。然而，m6A甲基化在肿瘤免疫微环境（TIME）中的作用尚不清楚，许多m6A甲基化调节因子（如METTL16和YTHDF3）在HNSCC中的预后作用尚不清楚，针对它们的治疗效果以及它们与PD-L1的相关性也不清楚。
常见的m6A甲基化调节因子：
甲基转移酶（“Writers"): METTL3, METTL14, WTAP,METTL16, KIAA1429,RBM15
去甲基化酶（”Erasers"):FTO、ALKBH5
与甲基化位点结合的蛋白质（“Reader")：YTHDC1, YTHDC2, YTHDF1, YTHDF2,YTHDF3, HNRNPC
二、分析流程
三、结果解读
1、HNSCC中m6A甲基化调节因子表达上调
图1 HNSCC中m6A甲基化调节因子表达水平上调（TCGA队列， *p <0.05, **p < 0.01, ***p < 0.001）
为了评估m6A甲基化调节因子在HNSCC启动和进展中的作用，作者在TCGA数据库中将样本分为HNSCC组织和癌旁正常组织，分析两组中15种m6A甲基化调节因子的表达水平。作者从TCGA中下载RNA-seq数据，获得mRNA表达谱，并分为43对肿瘤-正常组织的配对组，比较两种组织中m6A甲基化调节因子的表达水平的差异。图1A和图1B表明，m6A甲基化调节基因在肿瘤组织中表达上调。其中，"writers"（KIAA1429, METTL3, RBM15,WTAP, ZC3H13, METTL16)、“readers” ( HNRNPC, YTHDC1,YTHDF1, YTHDF2, YTHDF3), 以及“erasers” (ALKBH5 and FTO)在HNSCC组织中表达显著上调（p<0.01）。而对于”readers"YTHDC2，则在HNSCC组织中表达显著下调（p<0.05)。“writer” METTL14在正常和癌组织中表达无显著差异（p>0.05)。
以上结果表明m6A甲基化调节因子在HNSCC的进展中起重要作用。
2、基于m6A甲基化调节因子的聚类与HNSCC患者临床特征和生存的相关性
图2 TCGA队列中HNSCC cluster1/2亚型的临床病理特征差异及生存率差异
根据m6A甲基化调节因子的表达水平以及PAC（proportion of ambiguous clustering）方法，通过具有最低PAC的k值来推断最佳簇数，最终选择簇数k=2。图2A为m6A甲基化调节因子表达水平，将499名HNSCC患者分为cluster1(n=203)亚型和cluster2(n=296)。两组的m6A甲基化调节因子表达情况和临床特征的比较如图2B所示，个体m6A甲基化调节因子在cluster2中的表达低于cluster1，特别是METTL3和HNRNPC的表达水平。cluster2有更多的女性患者。cluster2的肿瘤WHO分级更低，cluster2中免疫评分显著升高。根据K-M曲线比较两组患者的生存率，如图2C和2D所示，cluster2患者的OS和DFS比cluster1患者的长。以上结果表明， 由m6A甲基化调节因子的表达水平定义的聚类亚型与HNSCC患者的异质性密切相关。
为了进一步研究15种m6A甲基化调节因子的相互作用，作者分析了它们之间的关系， 附图结果表明，YTHDC2、YTHDC1、YTHDF2、METTL14和RBM15中，两两的表达呈正相关， KIAA1429、ZC3H13、YTHDF3和FTO中，两两的表达也呈正相关。HNRNPC的表达水平与METTL3和WTAP的表达水平呈正相关。通过PCA，作者发现两个亚型的基因表达谱具有清晰的边界和显著的表达模式。
3、m6A RNA甲基化与PD-L1的关系以及HNSCC免疫图景
图3 PD-L1与m6A RNA甲基化的关系及HNSCC免疫细胞浸润的图景
为探究m6A甲基化与PD-L1的相关性，作者分析了两个亚型中m6A甲基化水平的差异以及两组中PD-L1和m6A甲基化的关系的差异。图3A表示，与正常组织相比，PD-L1在HNSCC中表达显著上调（p<0.05)。图3B表明，cluster2 中PD-L1的表达水平与cluster1相比显著升高（p<0.01)。如图3C所示，PD-L1与METTL3, RBM15, WTAP, HNRNPC, YTHDC2和YTHDF3呈显著正相关，PD-L1与KIAA1429, YTHDF1和FTO呈负相关，且以GSE65858作为外部验证集验证时，得到的结果一致。图3D展示了在TCGA的cluster1和2两组中，22种免疫细胞类型的浸润水平。
4、两个亚型中免疫细胞浸润的差异
图4 TCGA队列中两个亚型免疫细胞浸润的差异
图4A表明两个亚型在免疫评分上存在较大差异，其中cluster2免疫评分更高(p<0.05），且cluster2具有较好的预后。随后分析了两个亚型中22种免疫细胞的比例。通过图3D以及图4B-F，可看出cluster1具有较多的活化CD4初始T细胞，记忆T细胞、嗜酸性粒细胞；而cluster2具有较多的记忆B细胞、Tregs、NK细胞、中性粒细胞、M2巨噬细胞、活化的DC、肥大细胞。
作者进行GSEA富集分析以阐述两种亚型产生不同TIME的潜在机制。通过图4G发现，与肿瘤恶性特征有关的途径（G2M检查点、mTORC1信号通路、PI3K/AKT/mTOR1通路）都与cluster1有关。这表明，PI3K/AKT/mTOR通路可以解释两种亚型具有不同的TIME。
图5 TCGA队列中以去甲基化酶进行聚类生成的1/2亚型间生存和免疫浸润的逆转现象
为了验证作者的假设，对“Eraser”FTO和ALKBH5进行了一致性聚类。图5A可看出cluster1中的OS明显长于cluster2（p=0.0012）。然而，图5B和5C表明，cluster1的免疫评分和PD-L1表达水平低于cluster2。此外，图5D表明，22种免疫细胞类型的组成在两个集群之间没有显著差异。因此，去甲基化酶的聚集对患者生存和免疫浸润有相反的影响。
5、m6A甲基化调节因子预后信号的构建以及验证
图6 TCGA队列中，基于m6A甲基化调节因子的预后预测模型的构建与验证
作者为探究m6A甲基化调节因子对HNSCC患者的预后作用，将从TCGA数据集收集的样本分为训练集（n=351)和验证集(n=148)。附表提示，在TCGA训练和验证集中，年龄、性别、TNM分期、分级等基线特征没有统计学差异。
作者基于TCGA训练集中15种m6A甲基化调节因子的表达水平构建LASSO回归模型以筛选最具有预后价值的基因，共确定了7个调节因子，即METTL3、YTHDC2、METTL14、RBM15、YTHDC1、ZC3H13和HNRNPC（附图）。通过LASSO算法获得的系数，得到风险评分公式（风险评分=系数之和*m6A调节因子表达水平），并计算TCGA训练集和验证集的风险评分，根据评分的中位数将患者分为高风险和低风险两组。图6A和6B展示了两组患者中风险评分、OS、7个基于m6A调节因子表达谱的signature的差异。热图结果表明，高风险m6A甲基化调节因子，包括METTL3、YTHDC2、RBM15和YTHDC1，在高风险组中高表达，而保护性m6A甲基化调节因子的表达水平，包括METTL14、ZC3H13和HNRNPC在低风险组中高表达。通过绘制图6C&D的生存曲线，可知在TCGA训练和验证集中，低风险组的OS比高风险组长(p<0.0001）。为了评估7个风险signature的预后准确性，作者进行了3年和5年ROC曲线分析，发现TCGA训练集中3年和5年的AUC值分别为0.798和0.825，TCGA验证集中3年和5年的AUC值分别为0.775和0.814（附图）。AUC值表明，7个m6A调节因子signature对HNSCC患者的预后具有良好的鉴别性能，意味着基于7个风险signature的风险评分可准确预测HNSCC患者的预后。
6、 预后风险评分与HNSCC的分级、免疫评分和放疗有关
图7 TCGA队列中，预后评分与临床特征、免疫评分有关
接下来，作者讨论风险评分与临床特征的相关性。图7A显示了7个m6A在TCGA训练队列中的高风险和低风险组的表达水平。高风险组METTL3和YTHDC2的表达水平通常低于低风险组，高风险组HNRNPC和ZC3H13的表达水平高于低风险组。聚类亚型(p<0.001)、分级(p<0.05)和免疫评分(p<0.001)方面，高风险和低风险组之间差异显著。
作者还进一步研究了风险评分与亚型、分级和免疫评分之间的关系。图7B表明，cluster1的风险评分明显高于cluster2(p<0.001)，风险评分随着组织学分级的增加而增加，高免疫评分组的风险评分较低。这些发现表明，HNSCC患者的风险评分与亚型、分级和免疫评分显著相关。
图8 基于m6A甲基化调节因子风险评分模型的独立验证
为了验证m6A甲基化调节因子在免疫调节和患者生存方面的作用，作者以GSE65858作为外部验证集，根据风险评分公式计算GEO样本风险评分，并分为高低风险两组。图8A显示，与TCGA分析一致，低风险评分患者的OS比高风险评分患者更长。图8B显示了3年和5年ROC分析的结果，3年和5年的AUC值分别为0.854和0.846。图8C结果与上文一致，低风险组患者具有显著较高的免疫评分。图8D-F展示了在三种队列中，高风险评分组患者PD-L1的表达均显著下调。
作者还在TCGA训练集和验证集中，通过单因素和多因素Cox回归分析，探究风险评分是否为预测HNSCC患者预后的独立因素。附表结果表明，基于7个m6A甲基化调节因子signature建立的风险评分为预测预后的独立因素。
图9 在接受放射治疗和高风险评分的HNSCC患者中，生存获益显著
作者评估了接受放疗的HNSCC患者中风险评分的预后价值，以探究放疗是否会影响风险评分的预测性能。从图9A&C可知，接受放疗患者的OS显著延长。图9B&D表示低风险放疗生存效益相对不显著。通过比较图9A&B以及图9C&D可知，高风险患者接受放疗对OS的改善更为明显，即高风险患者接受放疗获益更多。
7、m6A调节因子signature的基因突变对免疫浸润细胞的影响
图10 风险评分与6种免疫细胞浸润丰度的相关性
作者探究风险评分和6种免疫浸润细胞的关系，以评估基于7个m6A调节因子signature对HNSCC肿瘤免疫微环境的影响。由图10A-C得知，风险评分与B细胞、CD4+T细胞、以及CD8+T细胞呈负相关；由图10D-F得知，风险评分与中性粒细胞、巨噬细胞、以及DC呈正相关。
图11 m6A调节因子基因突变对免疫细胞浸润的影响
最后，作者分析了基于m6A调节因子的体细胞拷贝数突变(CNAs)对免疫细胞浸润的影响，以阐明风险评分与免疫细胞浸润之间潜在的机制。由图11获悉，已鉴定的m6A调节因子特征的CNAs（包括arm-level缺失和arm-level增加），可显著影响B细胞、CD4T细胞、CD8T细胞、中性粒细胞、巨噬细胞和树突状细胞的浸润水平。
小结
基于TCGA训练和验证集以及GEO外部验证集，本研究系统地评估了m6A RNA调节因子在HNSCC中的预后价值、PD-L1的相关性、在TIME上的作用以及潜在的调控机制。PD-L1上调与m6A甲基化有关。通过对m6A调节因子的一致聚类，确定了两种HNSCC亚型（cluster1/2），这两种亚型对HNSCC患者的预后分层，同时呈现出明显不同的TIME。通过LASSO回归模型，得到7个具有预后特征的m6A调节因子的风险signature，并构建风险评分公式，且风险评分为预测患者预后的独立因素。高风险评分的患者可能从放疗中获益更多。m6A调节因子的风险signature与免疫细胞浸润有关。通过GSEA得知，m6A RNA甲基化可能与PI3K/AKT/mTOR信号通路协同参与HNSCC免疫微环境的调控。