基于心力衰竭HMOX2診斷標志物的機器學習模型

武淑琴，王遠涵，曹紅艷，張巖波,3

1.山西醫(yī)科大學基礎醫(yī)學院數(shù)學教研室，山西太原 030001；2.山西醫(yī)科大學公共衛(wèi)生學院衛(wèi)生統(tǒng)計教研室，山西太原 030001；3.重大疾病風險評估山西省重點實驗室，山西太原 030001

心力衰竭（heart failure，HF）是心臟病發(fā)展到嚴重階段導致心臟機械活動功能障礙的一種高發(fā)病率和高死亡率的臨床綜合征[1]。HF 通常是由其他心臟疾病引起，可以是原發(fā)性疾病，也可繼發(fā)于其他疾病。原發(fā)或繼發(fā)性的心臟疾病導致心臟功能減退，使心臟無法有效泵血，最終引發(fā)心力衰竭。近年來，免疫和炎癥被認為是影響HF 的常見病理生理學因素，基于生物信息學發(fā)現(xiàn)潛在的特征基因已成為一種更新穎、更可靠的疾病診斷方法。對疾病采用免疫治療已成為一種較新穎的臨床治療策略。有研究表明，SERCA2a 等基因已成為最有可能治療HF 的基因之一[2-4]。本研究選取基因表達綜合（gene expression omnibus，GEO）數(shù)據(jù)庫中的微陣列表達數(shù)據(jù)集，尋找HF 與正常樣本之間的差異表達基因（differentially expressed genes，DEGs），并通過生物富集分析探索潛在的生物學功能。同時采用最小絕對收縮和選擇算子（the least absolute shrinkage and selection operator，LASSO）、支持向量機遞歸特征消除（ support vector machine recursive feature elimination，SVM-RFE）和隨機森林（random forest，RF）三種高效機器學習算法對差異基因進行篩選，并進行特征基因外部數(shù)據(jù)集驗證[5]。使用R 軟件包“CIBERSORT”探討HF 組織免疫浸潤的結果，并綜合分析診斷基因與免疫細胞之間的相關性。

1 資料與方法

1.1 資料來源與處理

在GEO 數(shù)據(jù)庫中使用關鍵詞“HF、正常”或“HF、健康”，通過美國國家生物技術信息中心（National Center for Biotechnology Information，NCBI）數(shù)據(jù)庫平臺進行數(shù)據(jù)搜索。選擇的數(shù)據(jù)集類型是數(shù)組表達譜，生物類型是智人，數(shù)據(jù)集的樣本量大于20。使用R 包sva 消除來自不同平臺的數(shù)據(jù)批次效應。對該數(shù)據(jù)集的log2 轉換分位數(shù)歸一化信號強度進行校正，并輸出校正結果。篩選后獲取人類HF 樣本基因表達譜公共數(shù)據(jù)集GSE5406 和GSE84796 作為聯(lián)合芯片數(shù)據(jù)集，而GSE21610 作為獨立驗證數(shù)據(jù)集，下載各矩陣數(shù)據(jù)集文件及相應的平臺文件，見表1。

表1 訓練/驗證數(shù)據(jù)集的信息

1.2 方法

1.2.1 差異基因鑒定和基因功能富集分析 合并兩個數(shù)據(jù)集為訓練集，使用R 包“l(fā)imma”對標準化的基因表達譜進行差異分析，調整后P<0.05 和logFC>0.5被確定為DEGs 的顯著性標準。使用ggplot2 繪制DEGs 火山圖，用pheatmap 繪制DEGs 熱圖。使用clusterProfiler 包對差異基因進行基因本體論（gene ontology，GO）、京都基因和基因組數(shù)據(jù)庫（Kyoto Encyclopedia of Genes and Genomes，KEGG）、疾病本體論（disease ontology，DO）和基因集富集分析（gene set enrichment analysis，GSEA）來探索和分析差異基因的生物學功能[11]。GO、KEGG 和DO 的截止標準設置為校正后P<0.05。GSEA 分析中，ggplot2 包用于繪圖，GSEA 的截止值設置為|NES|>1.0、校正后P<0.05。

1.2.2 機器學習算法模型構建 LASSO 回歸的目的是得到預測誤差最小的變量結果及其對應的回歸系數(shù)。通過約束回歸系數(shù)，得到最優(yōu)結果[6]。用“glmnet”包獲取最佳值，具體計算方法：Cost()；其中，為懲罰系數(shù)，為向量。SVM-RFE 是基于SVM 的最大間隔原理的序列后向選擇其算法，通過模型訓練樣本，然后對每個特征得分進行排序，去掉最小得分的特征，用剩余特征繼續(xù)訓練模型，進行下次迭代，直到選出最佳的特征數(shù)，其算法用“e1071”包實現(xiàn)[7]。

1.2.3 特征基因的驗證 為進一步測試機器學習基因信號篩查效果，使用外部數(shù)據(jù)集（GSE21610）作為驗證數(shù)據(jù)集進行驗證。其受試者操作特征曲線（receiver operating characteristic curve，ROC 曲線）是使用R 包“pROC”生成，并使用曲線下面積（area under the curve，AUC）和95%CI來驗證特征基因的診斷效率。

1.2.4 免疫滲透評估 本研究探索免疫細胞在HF樣本中的浸潤，并分析診斷基因與免疫細胞之間的相關性，以尋找可能的病理生理過程。反卷積工具CIBERSORT 用于計算組織基因表達中的定量免疫細胞成分。比較HF 和正常人22 種免疫細胞的表達水平及免疫細胞之間的相關性。通過Spearman 秩相關分析得到診斷基因與免疫細胞亞型之間的關系，圖形用“ggplot2”包繪制。

1.3 統(tǒng)計學方法

數(shù)據(jù)預處理、模型的構建、驗證及評價均使用R 語言（版本4.2.2）及相關的R 程序包。所有統(tǒng)計學檢驗均為雙側檢驗，P<0.05 為差異有統(tǒng)計學意義。

2 結果

2.1 差異基因的識別

使用Affymetrix Human Genome U133 Plus 2.0 Array 分別分析GSE5406、GSE21610 數(shù)據(jù)集基因表達；使用Agilent-028004 SurePrint G3 Human GE 8x60K 微陣列分析GSE84796 數(shù)據(jù)集基因表達（表1）。根據(jù)篩選標準logFC>0.5 和調整后的P<0.05，獲得在正常樣本和HF 中差異表達的177個基因。采用火山圖描述數(shù)據(jù)集中所有差異基因的表達狀態(tài)，見圖1。上調和下調基因之間的差異明顯。

圖1 訓練數(shù)據(jù)集中所有差異基因的火山圖

圖3 三種機器學習算法示意圖

圖4 特征基因的相關信息

圖5 免疫細胞浸潤與通路活性分析

2.2 差異基因的功能富集分析

GO 分析3 個層面：生物過程（biological processes，BP）分析，差異基因主要參與細胞外基質組織、細胞外結構組織、外部封裝結構組織、細胞對細胞外刺激的反應、多機體繁殖等過程；細胞成分（cell component，CC）分析，差異基因主要富集在含膠原蛋白的細胞外基質、膠原蛋白三聚體、纖維狀膠原蛋白三聚體、帶狀膠原蛋白纖維和膠原蛋白三聚體復合體中；分子功能（molecular function，MF）分析，差異基因主要在細胞外基質結構成分、糖胺聚糖結合、肝素結合、膠原蛋白結合和細胞外基質結構成分賦予抗拉強度等方面發(fā)揮作用，見圖2A。在DO富集分析中，發(fā)現(xiàn)差異基因主要集中在一些心血管系統(tǒng)和炎癥性疾病中，包括心肌梗死、動脈硬化癥、肝炎、生殖細胞癌、腎功能不全、動脈硬化性心血管疾病等，見圖2B。KEGG 術語表明，差異基因主要參與PI3K-Akt 信號傳導、細胞因子–細胞因子受體相互作用、蛋白質的消化和吸收、cGMP-PKG 信號傳導等多條通路，從而在HF 的發(fā)生和發(fā)展中發(fā)揮重要作用，見圖2C。GSEA 富集分析發(fā)現(xiàn)，Th1-Th2細胞分化、MAPK 信號通路和B 細胞受體信號通路等與炎癥和免疫相關，且在HF 疾病和正常樣本中表現(xiàn)出顯著差異，見圖2D、2E。

2.3 機器學習確定HF 的診斷遺傳信號

LASSO 采用10 折交叉驗證確定16 個診斷基因信號，見圖3A；SVM-RFE 算法在特征數(shù)為6 時誤差最小，因此篩選出6 個最佳診斷信號，見圖3B；在RF 算法中，最終選擇26 棵樹作為參數(shù)，見圖3C，并把重要性評分>2 的5 個基因作為最佳診斷基因信號，見圖3D。整合三種機器學習算法確定HMOX2為HF 的診斷基因信號，見圖3E，其AUC 值為94.3%，見圖4A。

2.4 診斷基因信號的驗證

將基因信號HMOX2在數(shù)據(jù)集GSE21610中進行驗證，得出其在正常對照組中的表達高于HF 組，見圖4B，其AUC 值為83.3%，見圖4C。表明HMOX2具有較高的診斷效能，可作為潛在的生物標志物。

2.5 免疫細胞浸潤與通路活性分析

組織免疫細胞浸潤分析揭示227 個樣本的基因集合中的22 種免疫細胞亞型，見圖5A，發(fā)現(xiàn)3 種免疫細胞的差異：HF 疾病組CD8T 細胞（P<0.001）和中性粒細胞（P=0.008）顯著升高，而幼稚B 細胞（P=0.002）則為一組在HF 組織中低表達的細胞亞型，見圖5B。就免疫細胞亞型的相關性而言，巨噬細胞M0 與巨噬細胞M1 的正相關最顯著（r=0.35），活化的肥大細胞與靜息的肥大細胞負相關最顯著（r=–0.58），見圖5C。HF 特征基因與免疫細胞亞型的相關性分析顯示，HMOX2 與幼稚B 細胞呈最顯著正相關（r=0.24，P<0.001），與γ-δ 型T 細胞呈最顯著負相關（r=–0.15，P=0.027）。此外，HMOX2 還與漿細胞和巨噬細胞M1 相關，見圖5D。

3 討論

本研究表明慢性HF是一個全身性T細胞集落激活過程。在HF 進展中，T 細胞亞群（Th1-Th2）時空分布的變化對心室重構有較大的影響。Th2 細胞可能代表慢性HF 的治療靶點，有助于減少組織炎癥。此外，MAPK 信號通路通過血管緊張素Ⅱ與HF 發(fā)生關聯(lián)。通常認為抑制MAPK 信號通路可通過阻斷血管緊張素Ⅱ的分泌減輕HF 的發(fā)展。有研究發(fā)現(xiàn)B細胞通過依賴和獨立于抗體產生的機制在HF 的進展中發(fā)揮重要作用[9]。這也為B 細胞反應通路在HF中的作用提供新的見解。

HMOX2 又稱血紅素加氧酶-2，是血紅素加氧酶家族的組成性亞型成員，主要在大腦和睪丸中表達[10-12]。與長期以來一直是心血管研究重點的血紅素加氧酶-1（HMOX1）相比，對HMOX2 的研究仍處于起步階段。有報道稱，HMOX2 通過頸動脈中的BKCa2+通道在氧感知中發(fā)揮重要作用[13]。同時，HMOX2還通過調節(jié)細胞中的血紅素濃度及CO 和H2S 的水平影響多種生物過程。CO 可激活cGMP 信號通路、抑制強血管收縮劑內皮素-1 的產生，還可通過調節(jié)H2S 途徑調節(jié)腦血流量。這些生物過程大多是通過HMOX2 首先感知O2濃度而實現(xiàn)的[14]。據(jù)推測，HMOX2 能抑制缺氧疾病中的全身反應，但具體機制仍不清楚[15]。一般來說，氧化應激存在于大多數(shù)心血管疾病中[16]；其具體機制是大量的心肌細胞（心肌細胞、內皮細胞和中性粒細胞）可產生活性氧（reactive oxygen species，ROS）。在正常的生理條件下，心臟發(fā)揮防御性抗氧化功能，與ROS 的產生保持動態(tài)平衡。然而，在病理因素的刺激下，該平衡很快被改變，大量的ROS 被釋放，引起功能蛋白和脂質的過氧化及DNA 損傷，導致心肌收縮功能受損和細胞外基質重塑[17]。HMOX2 對清除細胞內的ROS 很重要，且在保護細胞免受ROS 引起的損傷方面也發(fā)揮重要作用[18]。本研究表明，HMOX2 的表達水平在HF 患者中降低，因此，根據(jù)HMOX2 活性表達情況，可對HF 診斷發(fā)揮一定的潛在作用。在HF治療方面，HMOX2 靶點的作用還在研究中，尚未得到充分的證據(jù)支持，需要更多的實驗和臨床研究確定HMOX2 靶點的潛在治療效果及其在HF管理中的具體作用機制。HMOX2 對HF 靶點基因的確定，可為HF 免疫治療提供新的思路。通過調控HMOX2 活性，可對HF 的分子治療具有潛在的作用。

本研究進一步探索HF 組織中免疫細胞的分布。其疾病組CD8T 細胞和中性粒細胞顯著升高，低表達的細胞亞型為幼稚B 細胞。研究發(fā)現(xiàn)HF 患者的CD8T 細胞、巨噬細胞、中性粒細胞的表達與正常人有顯著差異[19]。HF 后的組織損傷和纖維化伴隨復雜的免疫細胞反應，免疫細胞功能是HF 心肌損傷病理過程的中心環(huán)節(jié)，激活CD8T 細胞可抑制心室重塑。有研究發(fā)現(xiàn)，中性粒細胞在壓力超負荷引起的HF 實驗模型中發(fā)揮有害作用[20]。這一發(fā)現(xiàn)與中性粒細胞通過釋放脫顆粒和募集微泡參與多種心血管疾病發(fā)生、發(fā)展的基本作用機制相吻合。通過免疫細胞相關性棒棒糖圖探討HMOX2 與免疫細胞之間的相關性；免疫細胞相關性棒棒糖圖是一種熱圖的變種，其中每個細胞類型都用一列表示，每個基因用一行表示。通過觀察該圖，可判斷特征基因在免疫細胞中是否具有特異性表達，或某些免疫細胞類型中是否存在共同的表達模式。結果顯示：HMOX2 與幼稚B 細胞呈最顯著正相關，與γ-δ 型T 細胞呈最顯著負相關。此外，HMOX2 還與漿細胞和巨噬細胞M1 相關。這些結果表明，HMOX2 可通過調節(jié)免疫浸潤在HF 中發(fā)揮關鍵作用。

綜上所述，本研究通過LASSO 回歸、SVM-RFE和RF 三種機器學習算法篩選得到HMOX2診斷基因信號，可用作HF 診斷的特異生物標志物，為探索HF 的發(fā)病機制提供科學依據(jù)。