引用本文： 唐垚, 陈文栋, 王燕琼, 杨伟. 基于机器学习的体外循环后心肌缺血-再灌注损伤中m6A相关基因聚类分析和免疫浸润分析. 中国胸心血管外科临床杂志, 2024, 31(10): 1475-1485. doi: 10.7507/1007-4848.202312072 复制
体外循环是目前临床心脏大血管手术中常用的一种技术。由于其过程中心脏与循环脱离，不可避免地导致心肌缺血、缺氧，且会引发全身炎症和应激反应，从而增加术后并发症的风险，如出血、低心输出量、心律失常、电解质失衡、溶血和全身炎症等[1]。并且体外循环后心脏复灌流也会导致心肌缺血-再灌注损伤，但目前并没有防治体外循环后心肌缺血-再灌注损伤的特效方法。
N6-甲基腺苷（N6-methyladenosine，m6A）RNA 甲基化是真核生物中最丰富的 mRNA修饰形式，这种动态和可逆的m6A修饰影响几乎所有重要的生物学过程，m6A修饰失调涉及多种疾病且参与疾病的发生和发展，包括癌症及潜在的治疗靶点[2]。m6A相关基因可分为编码、消码、读码3种功能基因[3]，其中m6A修饰过程由甲基转移酶催化，发挥“编码器”（writers）功能；另外，去甲基化酶发挥“消码器”（erasers）功能，可从RNA中除去甲基，使 m6A修饰保持动态平衡[4]；此外，甲基化识别酶则充当“阅读器”（readers）功能，参与介导mRNA翻译、降解、稳定等特定功能[5]。已有研究[6-8]表明m6A甲基化与心肌缺血-再灌注损伤相关，且m6A甲基化过程对心血管系统的发育、分化、稳态维持和应激反应至关重要。心脏中存在许多m6A修饰，一项高通量基因测序研究[9]发现，m6A甲基化RNA广泛参与人和小鼠心脏的心肌发育、能量代谢、应激反应和心肌重塑。在机器学习中，聚类是一种用于根据相似性将不同对象分组到分离聚类中的技术，即相似对象将位于同一聚类中，与相似对象的其他聚类分开，基于数据自身信息对数据进行分类，将数据分为若干类别，使属于同一类别的数据间相似性尽可能大，不同类型相似性尽可能小，以达到分类的目的，将若干基因分成不同簇，从而降低研究对象的复杂性。免疫浸润分析即使用算法预测免疫细胞的细胞丰度，分析浸润免疫细胞的种类和比例。本研究利用GEO数据库中行体外循环的心脏病患者术前、术后的数据集，筛选该数据集中体外循环前后具有表达差异的m6A相关基因，基于表达谱中m6A相关的差异基因对该数据集样本进行一致性聚类、对聚类后的样本进行功能富集和免疫分析，并利用机器学习的方法筛选特征基因，针对特征基因进一步分析。心肌缺血-再灌注损伤是心脏直视下的体外循环手术尚未攻克的难题。在本研究中我们主要从基因层面寻找体外循环后心肌缺血-再灌注损伤中与m6A相关的特异性生物标志物，期望能为探索体外循环后心肌缺血-再灌注损伤的发病机制提供有价值的参考，并为预防及治疗体外循环后患者的心肌缺血-再灌注损伤提供思路，使麻醉或手术过程中有效预防和治疗心肌缺血-再灌注损伤成为可能。
1 材料与方法
1.1 数据获取
以“cardiopulmonary bypass”为关键词在GEO数据库（
通过在线网站RM2Target（
1.2 方法
1.2.1 m6A基因差异分析
应用R软件（R 4.3.0版）对标准化基因表达矩阵进行分析，提取m6A相关基因的表达量矩阵。然后利用limma包对实验组及对照组中m6A相关基因进行差异分析。
1.2.2 模型选择
应用R软件对机器学习中支持向量机（support vector machines，SVM）模型和随机森林树（random forest，RF）模型进行比较，基于残差的箱线图、残差的反向累积分布图来确定后续疾病特征基因筛选所采用的方法，利用受试者工作特征（receiver operating characteristic，ROC）曲线进行验证。
1.2.3 随机森林树
基于m6A相关的差异基因，利用R包“randomForest”，筛选重要性评分＞3的基因，利用R包“rms”构建这些特征基因的列线图，根据列线图中患者各基因的得分计算综合评分，以预测疾病的发病率。
1.2.4 m6A分型
基于m6A相关的差异基因，利用R包“ConsensusClusterPlus”对行体外循环的患者进行一致性聚类分析，根据累积分布函数曲线、一致性得分和一致性矩阵来确定K亚型的最佳数量。利用“limma、pheatmap、ggpubr”等R包，将特征基因在不同的m6A簇中的表达可视化，绘制箱线图和热图。运用主成分分析验证m6A簇的可靠性。
1.2.5 m6A分型的差异分析
基于m6A基因簇，运用R包筛选出不同分型之间的差异基因，P≤0.05为差异具有统计学意义，并进行可视化。
1.2.6 不同m6A簇中差异基因的GO和KEGG分析
为了探索m6A簇中差异基因的信号通路和潜在功能，我们使用R包“clusterProfiler、org.Hs.eg.db”等对这些基因进行基因本体论（Gene Ontology，GO）和京都基因与基因组百科全书（Kyoto Encyclopedia of Genes and Genomes，KEGG）富集分析。
1.2.7 不同m6A簇中免疫细胞与特征基因的相关性
通过R软件，运用Spearman相关性分析，计算m6A相关基因与免疫细胞的相关系数，探讨免疫细胞与特征基因的关系。
2 结果
2.1 m6A相关基因的获取
获取26个m6A相关基因，其中编码基因有METTL3、METTL14、METTL16、WTAP、VIRMA、ZC3H13、RBM15、RBM15B、CBLL1， 读码基因有YTHDC1、YTHDC2、YTHDF1、YTHDF2、YTHDF3、HNRNPC、FMR1、LRPPRC、HNRNPA2B1、IGFBP1、IGFBP2、IGFBP3、RBMX、ELAVL1、IGF2BP1，消码基因有FTO、ALKBH5（表1）。同时，将26个m6A相关基因在染色体上的位置可视化（图1a）。
图1
m6A相关基因在GSE132176中的表达谱
a：26个m6A相关基因在染色体上的相对位置；b：宏观显示在GSE132176具有差异表达的m6A相关基因，红色表示高表达，蓝色表示低表达；c：显示在GSE132176具有差异表达的m6A相关基因的变化，**：
表1
m6A相关基因
2.2 m6A基因差异分析
应用R软件（R 4.3.0版）对标准化基因表达矩阵进行分析，提取m6A相关基因的表达量矩阵。然后利用limma包对m6A相关基因进行差异分析，并将差异分析结果可视化。热图显示了在体外循环前后具有表达差异的m6A相关基因（图1b），从箱线图可以看出METTL3、METTL14、RBM15B、YTHDF1、HNRNPC在实验组及对照组之间具有表达差异，其中METTL3、METTL14及HNRNPC在体外循环术后下调，而RBM15B、YTHDF1上调（图1c）。
2.3 模型选择
我们建立了RF和SVM模型，以便从m6A相关基因中筛选特征基因以预测体外循环后心肌再灌注损伤的发生。应用R软件验证，在残差的箱线图中，RF模型残差低，表示模型准确性高（图2a）。在残差的反向累积分布图中，RF的残差反向累积分布线主要位于SVM的残差线内，表明预测值与实际值差异较小（图2b）。ROC曲线图显示RF模型更为准确（图2c）。因此我们选择RF模型进行后续分析。
图2
筛选特征基因的机器学习模型的确定
a：残差的箱线图，残差越低表示模型准确性越高；b：残差的反向累积分布图；c：受试者工作特征曲线，曲线下面积越大，模型准确性越高
2.4 RF模型
应用RF模型筛选特征基因，筛选重要性评分＞3的基因，结果表明METTL3、YTHDF1、RBM15B、METTL14为心肌缺血-再灌注损伤中的特征基因（图3a～b）。并利用R包“rms”构建了这些特征基因的列线图，结果显示METTL3可能对预测体外循环后心肌再灌注损伤风险贡献最大，其次是METTL14（图3c）。列线图校准曲线显示校正后预测风险模型符合理想状态（图3d）；临床决策曲线表明METTL3、RBM15B、YTHDF1、METTL14作为标志基因时，患者获益更高，这表明RF模型值得使用（图3e）。此外，临床影响曲线表明该模型具有良好的分类能力（图3f）。以上说明基于列线图模型可能有利于预测体外循环后心肌缺血-再灌注损伤，且METTL3、RBM15B、YTHDF1、METTL14或许可以作为体外循环后心肌缺血-再灌注损伤的生物标志物。
图3
建立机器学习模型并验证模型预测能力
a：对照组误差（红色），正常组误差（绿色），总样本误差（黑色）交叉验证的最小误差点及对应的最优森林树；b：体外循环后心肌缺血-再灌注损伤的m6A差异基因的重要性评分；c：METTL3、RBM15B、YTHDF1、METTL14的列线图；d：列线图校准曲线，实线表示校正后曲线，虚线表示理想现实状态；e：临床决策曲线，红色曲线离All曲线越远，说明列线图模型准确性越高；f：临床影响曲线，红色曲线表示每个阈值概念下阳性（高风险）例数，蓝色曲线表示每个阈值下真阳性数量
2.5 m6A分型
基于5个与m6A相关的差异基因，我们将最大的k值设置为9，发现累积分布函数曲线一致性指数波动最小。并且当k=2时，一致性得分相对较高，因此我们将体外循环后样本分为两个聚类簇（A和B）（图4a～d）。此外，主成分分析结果表明，这些m6A相关的差异基因可以基本区分这两个聚类簇（图4e）。
图4
基于m6A差异基因的样本簇分类
a：当
2.6 m6A簇的差异分析
利用R包筛选了m6A簇之间的差异基因，根据热图发现，多数m6A相关基因在B簇中上调，少数下调（图5a）。从m6A相关差异基因的箱线图中可以看出METTL14、RBM15B、YTHDF1在B簇中上调，而METTL3和HNRNPC在B簇中下调，且METTL3下调具有显著性（图5b）。
图5
m6A簇差异分析及富集分析
a：基于m6A相关差异基因进行分簇后的差异表达热图；b：m6A簇中m6A相关基因的箱线图；c～e：m6A簇差异基因的GO富集分析相关图；f～g：m6A簇差异基因的KEGG富集分析图
2.7 m6A簇中差异基因的GO和KEGG分析
利用R包对m6A簇分型中差异基因进行了GO和KEGG的富集分析。GO结果显示，这些基因在生物过程方面主要参与细胞对氧气水平的反应、单核细胞分化、上皮细胞分化的调控、肌肉器官的发育、对脂多糖的反应（GO：0032496）及对细菌来源分子的反应（GO：0002237）；在细胞组分方面主要参与RNA聚合酶Ⅱ转录调控复合体与粘附体接合；在分子功能方面主要参与DNA结合转录因子结合、DNA结合转录激活剂活性及RNA聚合酶Ⅱ转录调节因子复合物（GO：0090575）（图5c～e）。KEGG富集分析结果显示，这些基因主要参与脂质与动脉粥样硬化、内质网中的蛋白质加工、NOD样受体信号通路、肿瘤坏死因子（tumor necrosis factor，TNF）信号通路、流体剪切应力与动脉粥样硬化、白介素（interleukin，IL）-17信号通路（图5f～g）。
2.8 免疫细胞与特征基因的相关性
通过R软件，运用Spearman相关性分析，对5个差异表达的m6A相关基因与免疫细胞进行浸润分析，并对浸润水平进行相关性分析，并计算这些基因与免疫细胞的相关系数（图6a）。METTL3与肥大细胞、1型辅助性T淋巴细胞（Th1）、17型辅助性T淋巴细胞（Th17）免疫细胞显著相关。在METTL3低表达组，肥大细胞、Th1、Th17免疫细胞均高表达（图6b）。METTL14与巨噬细胞显著相关，在METTL14低表达组，巨噬细胞高表达（图6c），由此可见这些基因的表达与多种免疫细胞浸润水平有关，这表明这些特征基因可能参与心肌缺血-再灌注损伤发病机制中的免疫调节。
图6
不同m6A簇中免疫浸润
a：基因与免疫细胞的相关系数；b：METTL3与免疫细胞的相关性；c：METTL14与免疫细胞的相关性
2.9 构建特征基因调控网络
此外，我们利用NetworkAnalyst在线分析网站，针对METTL3、YTHDF1、RBM15B、METTL14这4个特征基因构建了基因miRNA网络（图7a～b）与基因转录因子网络（图7c），结果表明有许多miRNA和转录因子参与了这些特征基因的调节，这为针对这些基因的治疗方法提供了方向。
图7
m6A相关基因的调控网络
a～b：m6A相关基因的miRNA调控网络；c：m6A相关基因的转录因子调控网络
3 讨论
体外循环技术为术中心脏停跳患者的组织器官提供所需的血液及氧气，但术后心肌再灌注会造成心肌缺血-再灌注损伤。研究[12]表明，心肌缺血-再灌注损伤主要是由于过度的炎症反应、过量的自由基损伤、钙超载、自噬、凋亡等造成。但针对这些原因的治疗并未解决根本问题，为了进一步揭示心肌缺血-再灌注损伤机制，寻找有效的心肌缺血-再灌注损伤防治靶点，仍需不断地探索和研究。近几年研究[13]发现m6A甲基化在多种细胞内调控基因的表达，参与细胞更新、凋亡、自噬、免疫、增殖等多种生物学过程；也有研究[6]表明m6A甲基化在心肌缺血-再灌注损伤中发挥关键调控作用。
在本研究中，我们使用生物信息学和机器学习的方法，分析了体外循环再灌注后m6A相关基因表达谱，通过聚类分析、免疫浸润分析，筛选出了m6A相关的预测标志基因，并建立了心肌缺血-再灌注风险模型。我们通过对机器学习中SVM和RF模型进行了比较，确定使用RF模型对特征基因进行筛选，共获得4个特征基因：METTL3、YTHDF1、RBM15B、METTL14，为心肌缺血-再灌注损伤提供了比较准确的诊断价值。此外，我们还绘制了包含这4个特征基因的列线图，将4个特征基因结合起来，希望获得更好的诊断价值；通过对体外循环前后再灌注的样本进行差异分析，我们得到了5个m6A相关的差异表达基因，基于这些差异基因又进行聚类分析，将样本分为两类，之后对这两类m6A簇中差异表达的基因进行富集分析和免疫浸润分析。
肥大细胞起源于骨髓造血祖细胞，并迁移到外周组织中，最终在此定居并完成分化和成熟。肥大细胞的生物学功能包括先天免疫、免疫系统的免疫调节以及组织修复和血管生成[14]；肥大细胞还参与心肌缺血-再灌注晚期的组织纤维化和重塑[15]；肥大细胞分泌的TNF-α启动细胞因子级联反应促进心肌缺血-再灌注的组织损伤[16-17]；肥大细胞相关的类胰蛋白酶水平升高可加重心肌再灌注损伤，阻断胰蛋白酶释放可能有利于预防心肌梗死和改善经皮冠状动脉介入治疗后心肌再灌注[18]；另外，一些研究员[19]还发现肥大细胞脱颗粒成分对心血管系统具有保护作用，可以延长心肌细胞的存活时间并诱导新生血管形成。
Th1与Th17在宿主免疫反应中起关键作用，但其是一把双刃剑，反应过度或失衡都会引发炎症和过敏反应[20]。Th1细胞产生干扰素-γ（IFN-γ）、IL-2和TNF-β，引起细胞介导的免疫和吞噬依赖性炎症[21]。Th17细胞属于CD4+ Th细胞亚群，可产生IL-17A、IL-17F、IL-22和IL-21，并参与机体的炎症反应[22]；Th17产生的IL-17导致缺血性心力衰竭的纤维化和心室重塑，造成心肌损伤[23]；并且有研究[22]表明，在缺血-再灌注模型中Th17细胞显著增加。
巨噬细胞是一种很重要的免疫细胞，在免疫功能中发挥至关重要的作用，其对心肌缺血-再灌注损伤的调节也有极其重要的作用。根据表面标记物和功能，巨噬细胞分为两大亚型：与炎症反应有关的经典活化巨噬细胞（M1）、与再生和损伤修复有关的选择性活化巨噬细胞（M2）[24]。在心肌损伤时，M1巨噬细胞是最初对清除死细胞和细胞外基质碎片产生反应的主要亚型[14，25]；在心肌损伤的增殖期，M2巨噬细胞逐渐占主导地位，促进受损心脏组织的修复和再生[26]。本研究也发现巨噬细胞与心肌缺血-再灌注特征基因高度相关，所以针对特征基因，正确及时地调节巨噬细胞极化或许是一个很有前途治疗心肌缺血-再灌注损伤的靶点。
本研究中所得到的体外循环后心肌缺血-再灌注损伤相关的m6A特征基因有METTL3、METTL14、RBM15B与YTHDF1。已有众多研究表明这些基因与心肌损伤相关。其中，METTL3通过m6A修饰促进 DGCR8与 pri-miR-143-3p 的结合，从而增强 miR-143-3p 表达以抑制 PRKCE 转录并进一步加重心肌细胞焦亡和心肌损伤[27]。另外，METTL3还可以和lncRNA-SNHG8以及PTBP1和ALAS2之间相互作用，通过m6A修饰促进lncRNA-SNHG8，从而调节ALAS2诱导氧化应激并加重心肌损伤[28]。关于METTL14，有研究[29]表明将其敲低可显著抑制心肌缺血-再灌注诱导的心肌细胞凋亡和坏死，抑制氧化应激和炎症因子分泌，并在体外和体内激活 Akt/mTOR 途径，同时Akt/mTOR 通路激活可显著减轻METTL14敲低后心肌缺血-再灌注损伤引起的心肌细胞凋亡。也有研究[30]表明WTAP/YTHDF1/m6A/FOXO3a轴共同调节心肌缺血-再灌注损伤进展，WTAP通过m6A阅读器YTHDF1触发了FOXO3a mRNA上m6A修饰，敲低WTAP可显著减少细胞凋亡和炎性细胞因子的释放，改善心肌缺血-再灌注损伤的进展。对于RBM15B，它可能通过调节m6A甲基化来调节替代或非法剪接的功能，抑制前体 mRNA 剪接[31-32]，但其在心肌缺血-再灌注中的作用尚无研究阐明。
综上所述，现有研究在一定程度上表明我们筛选特征基因是可靠的。我们对此进行了进一步分析，旨在为体外循环后心肌再灌注损伤的治疗提供指导方向。我们的研究可能为阐述体外循环后心肌缺血-再灌注损伤提供有价值的参考，并为个性化治疗及免疫治疗提供指导，但由于GEO数据库相关研究数据较少，未来仍需生物学实验进行验证。
利益冲突：无。
作者贡献：唐垚、杨伟提出研究思路，负责研究设计；王燕琼负责数据分析；唐垚负责论文撰写；陈文栋、杨伟负责论文审阅与修改。
体外循环是目前临床心脏大血管手术中常用的一种技术。由于其过程中心脏与循环脱离，不可避免地导致心肌缺血、缺氧，且会引发全身炎症和应激反应，从而增加术后并发症的风险，如出血、低心输出量、心律失常、电解质失衡、溶血和全身炎症等[1]。并且体外循环后心脏复灌流也会导致心肌缺血-再灌注损伤，但目前并没有防治体外循环后心肌缺血-再灌注损伤的特效方法。
N6-甲基腺苷（N6-methyladenosine，m6A）RNA 甲基化是真核生物中最丰富的 mRNA修饰形式，这种动态和可逆的m6A修饰影响几乎所有重要的生物学过程，m6A修饰失调涉及多种疾病且参与疾病的发生和发展，包括癌症及潜在的治疗靶点[2]。m6A相关基因可分为编码、消码、读码3种功能基因[3]，其中m6A修饰过程由甲基转移酶催化，发挥“编码器”（writers）功能；另外，去甲基化酶发挥“消码器”（erasers）功能，可从RNA中除去甲基，使 m6A修饰保持动态平衡[4]；此外，甲基化识别酶则充当“阅读器”（readers）功能，参与介导mRNA翻译、降解、稳定等特定功能[5]。已有研究[6-8]表明m6A甲基化与心肌缺血-再灌注损伤相关，且m6A甲基化过程对心血管系统的发育、分化、稳态维持和应激反应至关重要。心脏中存在许多m6A修饰，一项高通量基因测序研究[9]发现，m6A甲基化RNA广泛参与人和小鼠心脏的心肌发育、能量代谢、应激反应和心肌重塑。在机器学习中，聚类是一种用于根据相似性将不同对象分组到分离聚类中的技术，即相似对象将位于同一聚类中，与相似对象的其他聚类分开，基于数据自身信息对数据进行分类，将数据分为若干类别，使属于同一类别的数据间相似性尽可能大，不同类型相似性尽可能小，以达到分类的目的，将若干基因分成不同簇，从而降低研究对象的复杂性。免疫浸润分析即使用算法预测免疫细胞的细胞丰度，分析浸润免疫细胞的种类和比例。本研究利用GEO数据库中行体外循环的心脏病患者术前、术后的数据集，筛选该数据集中体外循环前后具有表达差异的m6A相关基因，基于表达谱中m6A相关的差异基因对该数据集样本进行一致性聚类、对聚类后的样本进行功能富集和免疫分析，并利用机器学习的方法筛选特征基因，针对特征基因进一步分析。心肌缺血-再灌注损伤是心脏直视下的体外循环手术尚未攻克的难题。在本研究中我们主要从基因层面寻找体外循环后心肌缺血-再灌注损伤中与m6A相关的特异性生物标志物，期望能为探索体外循环后心肌缺血-再灌注损伤的发病机制提供有价值的参考，并为预防及治疗体外循环后患者的心肌缺血-再灌注损伤提供思路，使麻醉或手术过程中有效预防和治疗心肌缺血-再灌注损伤成为可能。
1 材料与方法
1.1 数据获取
以“cardiopulmonary bypass”为关键词在GEO数据库（
通过在线网站RM2Target（
1.2 方法
1.2.1 m6A基因差异分析
应用R软件（R 4.3.0版）对标准化基因表达矩阵进行分析，提取m6A相关基因的表达量矩阵。然后利用limma包对实验组及对照组中m6A相关基因进行差异分析。
1.2.2 模型选择
应用R软件对机器学习中支持向量机（support vector machines，SVM）模型和随机森林树（random forest，RF）模型进行比较，基于残差的箱线图、残差的反向累积分布图来确定后续疾病特征基因筛选所采用的方法，利用受试者工作特征（receiver operating characteristic，ROC）曲线进行验证。
1.2.3 随机森林树
基于m6A相关的差异基因，利用R包“randomForest”，筛选重要性评分＞3的基因，利用R包“rms”构建这些特征基因的列线图，根据列线图中患者各基因的得分计算综合评分，以预测疾病的发病率。
1.2.4 m6A分型
基于m6A相关的差异基因，利用R包“ConsensusClusterPlus”对行体外循环的患者进行一致性聚类分析，根据累积分布函数曲线、一致性得分和一致性矩阵来确定K亚型的最佳数量。利用“limma、pheatmap、ggpubr”等R包，将特征基因在不同的m6A簇中的表达可视化，绘制箱线图和热图。运用主成分分析验证m6A簇的可靠性。
1.2.5 m6A分型的差异分析
基于m6A基因簇，运用R包筛选出不同分型之间的差异基因，P≤0.05为差异具有统计学意义，并进行可视化。
1.2.6 不同m6A簇中差异基因的GO和KEGG分析
为了探索m6A簇中差异基因的信号通路和潜在功能，我们使用R包“clusterProfiler、org.Hs.eg.db”等对这些基因进行基因本体论（Gene Ontology，GO）和京都基因与基因组百科全书（Kyoto Encyclopedia of Genes and Genomes，KEGG）富集分析。
1.2.7 不同m6A簇中免疫细胞与特征基因的相关性
通过R软件，运用Spearman相关性分析，计算m6A相关基因与免疫细胞的相关系数，探讨免疫细胞与特征基因的关系。
2 结果
2.1 m6A相关基因的获取
获取26个m6A相关基因，其中编码基因有METTL3、METTL14、METTL16、WTAP、VIRMA、ZC3H13、RBM15、RBM15B、CBLL1， 读码基因有YTHDC1、YTHDC2、YTHDF1、YTHDF2、YTHDF3、HNRNPC、FMR1、LRPPRC、HNRNPA2B1、IGFBP1、IGFBP2、IGFBP3、RBMX、ELAVL1、IGF2BP1，消码基因有FTO、ALKBH5（表1）。同时，将26个m6A相关基因在染色体上的位置可视化（图1a）。
图1
m6A相关基因在GSE132176中的表达谱
a：26个m6A相关基因在染色体上的相对位置；b：宏观显示在GSE132176具有差异表达的m6A相关基因，红色表示高表达，蓝色表示低表达；c：显示在GSE132176具有差异表达的m6A相关基因的变化，**：
表1
m6A相关基因
2.2 m6A基因差异分析
应用R软件（R 4.3.0版）对标准化基因表达矩阵进行分析，提取m6A相关基因的表达量矩阵。然后利用limma包对m6A相关基因进行差异分析，并将差异分析结果可视化。热图显示了在体外循环前后具有表达差异的m6A相关基因（图1b），从箱线图可以看出METTL3、METTL14、RBM15B、YTHDF1、HNRNPC在实验组及对照组之间具有表达差异，其中METTL3、METTL14及HNRNPC在体外循环术后下调，而RBM15B、YTHDF1上调（图1c）。
2.3 模型选择
我们建立了RF和SVM模型，以便从m6A相关基因中筛选特征基因以预测体外循环后心肌再灌注损伤的发生。应用R软件验证，在残差的箱线图中，RF模型残差低，表示模型准确性高（图2a）。在残差的反向累积分布图中，RF的残差反向累积分布线主要位于SVM的残差线内，表明预测值与实际值差异较小（图2b）。ROC曲线图显示RF模型更为准确（图2c）。因此我们选择RF模型进行后续分析。
图2
筛选特征基因的机器学习模型的确定
a：残差的箱线图，残差越低表示模型准确性越高；b：残差的反向累积分布图；c：受试者工作特征曲线，曲线下面积越大，模型准确性越高
2.4 RF模型
应用RF模型筛选特征基因，筛选重要性评分＞3的基因，结果表明METTL3、YTHDF1、RBM15B、METTL14为心肌缺血-再灌注损伤中的特征基因（图3a～b）。并利用R包“rms”构建了这些特征基因的列线图，结果显示METTL3可能对预测体外循环后心肌再灌注损伤风险贡献最大，其次是METTL14（图3c）。列线图校准曲线显示校正后预测风险模型符合理想状态（图3d）；临床决策曲线表明METTL3、RBM15B、YTHDF1、METTL14作为标志基因时，患者获益更高，这表明RF模型值得使用（图3e）。此外，临床影响曲线表明该模型具有良好的分类能力（图3f）。以上说明基于列线图模型可能有利于预测体外循环后心肌缺血-再灌注损伤，且METTL3、RBM15B、YTHDF1、METTL14或许可以作为体外循环后心肌缺血-再灌注损伤的生物标志物。
图3
建立机器学习模型并验证模型预测能力
a：对照组误差（红色），正常组误差（绿色），总样本误差（黑色）交叉验证的最小误差点及对应的最优森林树；b：体外循环后心肌缺血-再灌注损伤的m6A差异基因的重要性评分；c：METTL3、RBM15B、YTHDF1、METTL14的列线图；d：列线图校准曲线，实线表示校正后曲线，虚线表示理想现实状态；e：临床决策曲线，红色曲线离All曲线越远，说明列线图模型准确性越高；f：临床影响曲线，红色曲线表示每个阈值概念下阳性（高风险）例数，蓝色曲线表示每个阈值下真阳性数量
2.5 m6A分型
基于5个与m6A相关的差异基因，我们将最大的k值设置为9，发现累积分布函数曲线一致性指数波动最小。并且当k=2时，一致性得分相对较高，因此我们将体外循环后样本分为两个聚类簇（A和B）（图4a～d）。此外，主成分分析结果表明，这些m6A相关的差异基因可以基本区分这两个聚类簇（图4e）。
图4
基于m6A差异基因的样本簇分类
a：当
2.6 m6A簇的差异分析
利用R包筛选了m6A簇之间的差异基因，根据热图发现，多数m6A相关基因在B簇中上调，少数下调（图5a）。从m6A相关差异基因的箱线图中可以看出METTL14、RBM15B、YTHDF1在B簇中上调，而METTL3和HNRNPC在B簇中下调，且METTL3下调具有显著性（图5b）。
图5
m6A簇差异分析及富集分析
a：基于m6A相关差异基因进行分簇后的差异表达热图；b：m6A簇中m6A相关基因的箱线图；c～e：m6A簇差异基因的GO富集分析相关图；f～g：m6A簇差异基因的KEGG富集分析图
2.7 m6A簇中差异基因的GO和KEGG分析
利用R包对m6A簇分型中差异基因进行了GO和KEGG的富集分析。GO结果显示，这些基因在生物过程方面主要参与细胞对氧气水平的反应、单核细胞分化、上皮细胞分化的调控、肌肉器官的发育、对脂多糖的反应（GO：0032496）及对细菌来源分子的反应（GO：0002237）；在细胞组分方面主要参与RNA聚合酶Ⅱ转录调控复合体与粘附体接合；在分子功能方面主要参与DNA结合转录因子结合、DNA结合转录激活剂活性及RNA聚合酶Ⅱ转录调节因子复合物（GO：0090575）（图5c～e）。KEGG富集分析结果显示，这些基因主要参与脂质与动脉粥样硬化、内质网中的蛋白质加工、NOD样受体信号通路、肿瘤坏死因子（tumor necrosis factor，TNF）信号通路、流体剪切应力与动脉粥样硬化、白介素（interleukin，IL）-17信号通路（图5f～g）。
2.8 免疫细胞与特征基因的相关性
通过R软件，运用Spearman相关性分析，对5个差异表达的m6A相关基因与免疫细胞进行浸润分析，并对浸润水平进行相关性分析，并计算这些基因与免疫细胞的相关系数（图6a）。METTL3与肥大细胞、1型辅助性T淋巴细胞（Th1）、17型辅助性T淋巴细胞（Th17）免疫细胞显著相关。在METTL3低表达组，肥大细胞、Th1、Th17免疫细胞均高表达（图6b）。METTL14与巨噬细胞显著相关，在METTL14低表达组，巨噬细胞高表达（图6c），由此可见这些基因的表达与多种免疫细胞浸润水平有关，这表明这些特征基因可能参与心肌缺血-再灌注损伤发病机制中的免疫调节。
图6
不同m6A簇中免疫浸润
a：基因与免疫细胞的相关系数；b：METTL3与免疫细胞的相关性；c：METTL14与免疫细胞的相关性
2.9 构建特征基因调控网络
此外，我们利用NetworkAnalyst在线分析网站，针对METTL3、YTHDF1、RBM15B、METTL14这4个特征基因构建了基因miRNA网络（图7a～b）与基因转录因子网络（图7c），结果表明有许多miRNA和转录因子参与了这些特征基因的调节，这为针对这些基因的治疗方法提供了方向。
图7
m6A相关基因的调控网络
a～b：m6A相关基因的miRNA调控网络；c：m6A相关基因的转录因子调控网络
3 讨论
体外循环技术为术中心脏停跳患者的组织器官提供所需的血液及氧气，但术后心肌再灌注会造成心肌缺血-再灌注损伤。研究[12]表明，心肌缺血-再灌注损伤主要是由于过度的炎症反应、过量的自由基损伤、钙超载、自噬、凋亡等造成。但针对这些原因的治疗并未解决根本问题，为了进一步揭示心肌缺血-再灌注损伤机制，寻找有效的心肌缺血-再灌注损伤防治靶点，仍需不断地探索和研究。近几年研究[13]发现m6A甲基化在多种细胞内调控基因的表达，参与细胞更新、凋亡、自噬、免疫、增殖等多种生物学过程；也有研究[6]表明m6A甲基化在心肌缺血-再灌注损伤中发挥关键调控作用。
在本研究中，我们使用生物信息学和机器学习的方法，分析了体外循环再灌注后m6A相关基因表达谱，通过聚类分析、免疫浸润分析，筛选出了m6A相关的预测标志基因，并建立了心肌缺血-再灌注风险模型。我们通过对机器学习中SVM和RF模型进行了比较，确定使用RF模型对特征基因进行筛选，共获得4个特征基因：METTL3、YTHDF1、RBM15B、METTL14，为心肌缺血-再灌注损伤提供了比较准确的诊断价值。此外，我们还绘制了包含这4个特征基因的列线图，将4个特征基因结合起来，希望获得更好的诊断价值；通过对体外循环前后再灌注的样本进行差异分析，我们得到了5个m6A相关的差异表达基因，基于这些差异基因又进行聚类分析，将样本分为两类，之后对这两类m6A簇中差异表达的基因进行富集分析和免疫浸润分析。
肥大细胞起源于骨髓造血祖细胞，并迁移到外周组织中，最终在此定居并完成分化和成熟。肥大细胞的生物学功能包括先天免疫、免疫系统的免疫调节以及组织修复和血管生成[14]；肥大细胞还参与心肌缺血-再灌注晚期的组织纤维化和重塑[15]；肥大细胞分泌的TNF-α启动细胞因子级联反应促进心肌缺血-再灌注的组织损伤[16-17]；肥大细胞相关的类胰蛋白酶水平升高可加重心肌再灌注损伤，阻断胰蛋白酶释放可能有利于预防心肌梗死和改善经皮冠状动脉介入治疗后心肌再灌注[18]；另外，一些研究员[19]还发现肥大细胞脱颗粒成分对心血管系统具有保护作用，可以延长心肌细胞的存活时间并诱导新生血管形成。
Th1与Th17在宿主免疫反应中起关键作用，但其是一把双刃剑，反应过度或失衡都会引发炎症和过敏反应[20]。Th1细胞产生干扰素-γ（IFN-γ）、IL-2和TNF-β，引起细胞介导的免疫和吞噬依赖性炎症[21]。Th17细胞属于CD4+ Th细胞亚群，可产生IL-17A、IL-17F、IL-22和IL-21，并参与机体的炎症反应[22]；Th17产生的IL-17导致缺血性心力衰竭的纤维化和心室重塑，造成心肌损伤[23]；并且有研究[22]表明，在缺血-再灌注模型中Th17细胞显著增加。
巨噬细胞是一种很重要的免疫细胞，在免疫功能中发挥至关重要的作用，其对心肌缺血-再灌注损伤的调节也有极其重要的作用。根据表面标记物和功能，巨噬细胞分为两大亚型：与炎症反应有关的经典活化巨噬细胞（M1）、与再生和损伤修复有关的选择性活化巨噬细胞（M2）[24]。在心肌损伤时，M1巨噬细胞是最初对清除死细胞和细胞外基质碎片产生反应的主要亚型[14，25]；在心肌损伤的增殖期，M2巨噬细胞逐渐占主导地位，促进受损心脏组织的修复和再生[26]。本研究也发现巨噬细胞与心肌缺血-再灌注特征基因高度相关，所以针对特征基因，正确及时地调节巨噬细胞极化或许是一个很有前途治疗心肌缺血-再灌注损伤的靶点。
本研究中所得到的体外循环后心肌缺血-再灌注损伤相关的m6A特征基因有METTL3、METTL14、RBM15B与YTHDF1。已有众多研究表明这些基因与心肌损伤相关。其中，METTL3通过m6A修饰促进 DGCR8与 pri-miR-143-3p 的结合，从而增强 miR-143-3p 表达以抑制 PRKCE 转录并进一步加重心肌细胞焦亡和心肌损伤[27]。另外，METTL3还可以和lncRNA-SNHG8以及PTBP1和ALAS2之间相互作用，通过m6A修饰促进lncRNA-SNHG8，从而调节ALAS2诱导氧化应激并加重心肌损伤[28]。关于METTL14，有研究[29]表明将其敲低可显著抑制心肌缺血-再灌注诱导的心肌细胞凋亡和坏死，抑制氧化应激和炎症因子分泌，并在体外和体内激活 Akt/mTOR 途径，同时Akt/mTOR 通路激活可显著减轻METTL14敲低后心肌缺血-再灌注损伤引起的心肌细胞凋亡。也有研究[30]表明WTAP/YTHDF1/m6A/FOXO3a轴共同调节心肌缺血-再灌注损伤进展，WTAP通过m6A阅读器YTHDF1触发了FOXO3a mRNA上m6A修饰，敲低WTAP可显著减少细胞凋亡和炎性细胞因子的释放，改善心肌缺血-再灌注损伤的进展。对于RBM15B，它可能通过调节m6A甲基化来调节替代或非法剪接的功能，抑制前体 mRNA 剪接[31-32]，但其在心肌缺血-再灌注中的作用尚无研究阐明。
综上所述，现有研究在一定程度上表明我们筛选特征基因是可靠的。我们对此进行了进一步分析，旨在为体外循环后心肌再灌注损伤的治疗提供指导方向。我们的研究可能为阐述体外循环后心肌缺血-再灌注损伤提供有价值的参考，并为个性化治疗及免疫治疗提供指导，但由于GEO数据库相关研究数据较少，未来仍需生物学实验进行验证。
利益冲突：无。
作者贡献：唐垚、杨伟提出研究思路，负责研究设计；王燕琼负责数据分析；唐垚负责论文撰写；陈文栋、杨伟负责论文审阅与修改。
