引用本文： 唐垚, 陳文棟, 王燕瓊, 楊偉. 基于機器學習的體外循環后心肌缺血-再灌注損傷中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數據庫相關研究數據較少，未來仍需生物學實驗進行驗證。
利益沖突：無。
作者貢獻：唐垚、楊偉提出研究思路，負責研究設計；王燕瓊負責數據分析；唐垚負責論文撰寫；陳文棟、楊偉負責論文審閱與修改。
