錢 晗, 曹夢菲, 呂書梅, 袁 偉

(江蘇大學附屬醫(yī)院 心血管內(nèi)科, 江蘇 鎮(zhèn)江, 212001)

在人口老齡化和代謝危險因素的雙重壓力下，中國居民的心血管疾病風險持續(xù)增加，心肌梗死患者的病死率亦逐漸上升[1]。心肌梗死是指心肌缺血性壞死，通常是在冠狀動脈病變的基礎上發(fā)生冠狀動脈血供減少或中斷，心肌細胞因嚴重且持續(xù)的缺血缺氧而發(fā)生壞死[2]。心肌細胞的死亡會引發(fā)一系列免疫炎癥反應，這些免疫炎癥反應介導的受損心肌修復在梗死后心室重構(gòu)、心功能和患者結(jié)局等方面起著關(guān)鍵作用[3]。急性心肌梗死的診斷主要依據(jù)臨床表現(xiàn)、血清生物學標志物和心電圖，其中臨床常用的生物學標志物包括肌紅蛋白、肌酸激酶-MB、心肌肌鈣蛋白I和心肌肌鈣蛋白T[4-5]。然而，這些標志物水平在心力衰竭、腎衰竭和甲狀腺疾病患者中也同樣升高，故結(jié)合多組學平臺探尋更多特異性的生物學標志物非常必要[6-7]。外周血單個核細胞(PBMCs)包括T細胞、B細胞、單核細胞、樹突狀細胞和自然殺傷細胞，在免疫應答、免疫監(jiān)測和免疫治療中發(fā)揮著重要作用[8]。本研究基于PBMCs分析心肌梗死和梗死恢復過程中潛在的關(guān)鍵基因和生物學改變，旨在為心肌梗死的診斷、評估、管理和預后分析提供新的見解。

1 材料和方法

1.1 數(shù)據(jù)集獲取

從GEO數(shù)據(jù)庫中下載關(guān)于心肌梗死PBMCs的轉(zhuǎn)錄組測序芯片數(shù)據(jù)集GSE59867[9-10], 芯片注釋平臺為GPL6244, 芯片信息為Affymetrix Human Gene 1.0 ST Array。該數(shù)據(jù)集收集了111例ST段抬高型心肌梗死(STEMI)患者入院、出院、心肌梗死后1個月、心肌梗死后6個月時的外周血樣本數(shù)據(jù)。對照組為46例無心肌梗死病史的穩(wěn)定型冠心病患者。基于芯片注釋平臺對數(shù)據(jù)集進行基因名轉(zhuǎn)換，用于后續(xù)共表達網(wǎng)絡的構(gòu)建。

1.2 加權(quán)基因共表達網(wǎng)絡分析(WGCNA)

通過WGCNA R包篩選與心肌梗死有關(guān)的易感基因模塊[11]。WGCNA是一種常用的生物信息學方法，可以識別具有相似表達模式的基因模塊，分析基因模塊與樣本表型之間的關(guān)系，繪制基因模塊中的調(diào)控網(wǎng)絡，識別關(guān)鍵模塊和調(diào)控基因[12-13]。首先，選取數(shù)據(jù)集中方差變異程度最高的前3 000個基因構(gòu)建共表達網(wǎng)絡; 其次，利用R函數(shù)pickSoftThreshold計算軟閾值1～20; 然后，對這些基因進行層次聚類和動態(tài)樹切割來劃分模塊，并合并相似的模塊(閾值為0.25); 最后，計算模塊與疾病之間的Spearman相關(guān)系數(shù)，進而確定易感基因模塊。為了探討所選模塊的潛在生物學功能，通過ClusterProfiler R包對模塊中的基因進行基因本體論(GO)分析和京都基因與基因組百科全書(KEGG)通路分析[14]。GO分析包括生物過程(BP)、分子功能(MF)和細胞成分(CC)。通過Ggplot2 R包進行結(jié)果展示。

1.3 基因差異分析和轉(zhuǎn)錄因子預測

使用Limma R包進行基因差異分析[15]。差異表達基因(DEGs)的篩選標準為矯正后P<0.05, |log2(倍數(shù)變化)|>0.5?；贑lusterProfiler R包對DEGs進行GO分析和KEGG通路分析,P<0.05表示顯著富集。對上述易感基因模塊和DEGs進行交集分析，通過ChEA3數(shù)據(jù)庫對交集基因進行轉(zhuǎn)錄因子預測。ChEA3是一種轉(zhuǎn)錄因子富集分析工具，整合了多個數(shù)據(jù)庫，如ENCODE、GTEx、ARCHS4等數(shù)據(jù)庫[16]。由于ChEA3數(shù)據(jù)庫并未對預測結(jié)果可靠度進行嚴格劃分，參考相關(guān)研究[16-17]結(jié)論，本研究對平均排名分≤30分的共同轉(zhuǎn)錄因子進行結(jié)果展示。基于STRING數(shù)據(jù)庫構(gòu)建交集基因的蛋白相互作用網(wǎng)絡，并將交互評分>0.4分的數(shù)據(jù)導入Cytoscape 3.8.2軟件進行可視化處理[18-19]。

1.4 關(guān)鍵基因的確定和受試者工作特征(ROC)曲線分析

LASSO回歸分析是常用的篩選變量的壓縮估計方法，通過構(gòu)造懲罰函數(shù)獲得較為精煉的模型，使得其壓縮一些回歸系數(shù)，同時設定一些回歸系數(shù)為0，因此保留了子集收縮的優(yōu)點，是一種處理具有復共線性數(shù)據(jù)的有偏估計，常被用于協(xié)變量篩選[20]?；趃lmnet R包對交集基因中存在蛋白相互作用的基因進行LASSO回歸分析，進一步篩選關(guān)鍵基因[21]。采用10倍交叉驗證法選擇懲罰項(λ), 并選擇最小二項偏差在1個標準誤差內(nèi)的最簡單模型的λ值(lambda.1se)。通過pROC R包對關(guān)鍵基因進行ROC曲線分析，計算曲線下面積(AUC), 評估關(guān)鍵基因?qū)π募」Ｋ赖脑\斷價值[22]。

1.5 基因集富集分析(GSEA)和外部數(shù)據(jù)集驗證

根據(jù)每個關(guān)鍵基因的表達量將數(shù)據(jù)集劃分為高表達組和低表達組，通過GSEA探討高表達組和低表達組之間潛在的生物學功能變化，將特征基因集作為預先定義的基因集。作為一種常見的生物信息學方法, GSEA可評估一個預先定義的基因集的基因在與表型相關(guān)度排序的基因表中的分布趨勢，從而判斷其對表型的影響[23]。從GEO數(shù)據(jù)庫中選取另一個數(shù)據(jù)集GSE123342, 將其用于驗證關(guān)鍵基因的表達。GSE123342數(shù)據(jù)集包括急性心肌梗死(65例)、梗死后30 d(64例)、梗死后1年(37例)、穩(wěn)定型冠心病(22例)樣本數(shù)據(jù)以及4例技術(shù)重復樣本數(shù)據(jù)。

2 結(jié) 果

2.1 易感基因模塊的確定

根據(jù)無標度網(wǎng)絡原理，選取10作為共表達網(wǎng)絡的軟閾值，見圖1A。通過層次聚類和動態(tài)樹切割，將前3 000個基因劃分為15個基因模塊，見圖1B。剔除聚類失敗的灰色模塊，棕色模塊的總體表達與心肌梗死的相關(guān)性最高(r=0.4), 且與出院至心肌梗死后6個月的演變過程呈負相關(guān)(r=-0.33), 見圖1C。棕色模塊中的共表達基因伴隨著心肌梗死的發(fā)生呈現(xiàn)高表達趨勢，而在心肌梗死恢復過程中則呈現(xiàn)低表達趨勢。由此提示，棕色模塊可能是易感基因模塊，其中包含243個基因，即ABHD5、ACER3、ACP3、ACSL1、ADAM17、ADAM9、ADM、AHR、ALCAM、ALDH2、ALPK1、ANKRD50、ANO5、AQP9、ARHGAP24、ARHGAP29、ARRDC4、ASPH、ATP6V0A1、ATP6V1A、B3GNT5、BACH1、BLVRB、BST1、C3AR1、C5、C9orf72、CALCRL、CAPG、CAPZA2、CARD6、CASP1、CASP5、CCDC88A、CCR2、CD14、CD163、CD1D、CD33、CD36、CD63、CLEC1A、CLEC4A、CLEC4D、CLEC4E、CLEC7A、CLMN、CNTLN、CPED1、CPM、CPNE8、CPVL、CR1、CR1L、CREG1、CRISPLD2、CSTA、CTNNA1、CTSH、CYBRD1、CYP1B1、CYP1B1-AS1、DACH1、DDIAS、DOCK4、DOCK5、DRAM1、DSC2、DSE、DUSP6、EDNRB、ENTPD1、ERLIN1、EVI5、F5、FAM114A1、FAM13A、FAM151B、FAM20A、FAM198B、FAR2、FBN2、FBP1、FCGR2A、FGD4、FGD6、FLT3、FLVCR2、FMO5、FPR1、FPR2、FRRS1、FUCA1、FUCA2、GAPT、GAS2L3、GCA、GIMAP8、GLT1D1、GPAT3、GPR141、HAL、HGF、HMGB2、HNMT、HORMAD1、HP、HPSE、IDH1、IFNGR1、IL13RA1、IL15、IL18、IMPA2、IRAK3、JAK2、KCNJ15、KCNJ2、KIF13A、KLHL8、KYNU、LHFPL2、LILRA5、LIN7A、LIPN、LMNB1、LPCAT2、LRMDA、LRRK2、LTA4H、LY86、LYVE1、MAP3K20、MCEMP1、ME1、METTL7A、METTL7B、MFSD1、MGST1、MGST2、MILR1、MNDA、MOSPD2、MS4A4A、MSRB1、MTARC1、MYOF、NAAA、NAIP、NETO2、NLN、NLRC4、NPL、OLFML2B、P2RX7、P2RY13、PGD、PIP4P2、PLA2G4A、PLBD1、PLD1、PLIN2、PLSCR1、PPARG、PPT1、PRRG4、PSTPIP2、PYGL、QPCT、RAB39A、RALB、RBP7、RGL1、RNASE2、RNASE6、RNF141、RNF217、RRAGD、S100A12、S100A9、S100Z、SAT2、SEMA3C、SEPTIN10、SERPINB1、SERPINB2、SESTD1、SGMS2、SH3PXD2B、SHTN1、SIPA1L2、SIRPD、SLC15A2、SLC1A3、SLC22A15、SLC22A4、SLC26A8、SLC31A1、SLC36A4、SLC7A7、SLC8A1、SLITRK4、SMAD1、SMPDL3A、SNX10、SOD2、SORT1、SPATA6、ST3GAL6、ST6GALNAC3、STEAP4、SULT1B1、TASL、TBC1D12、TCN2、TDRD9、TFEC、TGFBI、TLR1、TLR2、TLR4、TLR5、TLR6、TLR7、TLR8、TM6SF1、TMEM144、TMEM167A、TMTC2、TNFAIP6、TNFSF13B、TPST1、TSPO、UBE2D1、UGGT2、VNN2、VNN3、WASF1、WDFY3、WLS、ZC3H12C、ZFYVE16、ZNF438。

棕色模塊GO分析最顯著的結(jié)果分別是骨髓白細胞激活(BP)、囊泡(CC)和模式識別受體活性(MF), 見圖1D; KEGG分析結(jié)果則主要為感染、免疫炎癥等相關(guān)生物學過程的改變，包括Toll樣受體信號通路等，見圖1E。

2.2 基因差異分析

進一步對數(shù)據(jù)集進行基因差異分析，根據(jù)篩選標準共得到142個DEGs(圖2A), 其中上調(diào)的DEGs有77個，下調(diào)的DEGs有65個，見表1。BP分析結(jié)果主要包括免疫反應、免疫系統(tǒng)過程和免疫效應過程, CC分析結(jié)果包括分泌顆粒、分泌小泡和細胞質(zhì)小泡部分等, MF分析結(jié)果主要包括信號受體活性、分子傳感器活性和碳水化合物結(jié)合，見圖2B。KEGG分析提示，這些DEGs主要涉及免疫、炎癥和感染的相關(guān)生物學改變，見圖2C。DEGs與上述棕色模塊存在35個交集基因，見圖2D?；贑hEA3數(shù)據(jù)庫對35個交集基因進行共同轉(zhuǎn)錄因子預測，其中平均排名分≤30分的共同轉(zhuǎn)錄因子分別是CREB5、MTF1、NFE4、SPI1、ZNF467、NFE2、TFEC、MXD1、NR1H3、BORCS8MEF2B、NFIL3和CEBPE, 見表2。構(gòu)建蛋白互作網(wǎng)絡，其中有18個基因存在蛋白質(zhì)相互作用，見圖2E。

表1 差異基因與交集基因

表2 共同轉(zhuǎn)錄因子預測結(jié)果(平均排名分≤30分)

2.3 關(guān)鍵基因的確定

對18個存在蛋白相互作用的基因進一步行LASSO回歸分析，這些基因不同懲罰參數(shù)值所對應的系數(shù)見圖3A。本研究選取最小二項偏差在1個標準誤差內(nèi)的λ值(lambda.1se, 6個基因)，該λ值提供更精簡的模型，見圖3B。ROC曲線分析顯示,CD163、RNASE2、HP、FAM20A、MCEMP1和FAM198B基因表達水平對心肌梗死的發(fā)生均具有良好的診斷價值,AUC分別為0.831、0.798、0.775、0763、0.866和0.829, 見圖3C。與對照組穩(wěn)定型冠心病患者相比，這些關(guān)鍵基因的表達水平在入院時心肌梗死患者中顯著上升，并在出院時、心肌梗死后1個月、心肌梗死后6個月逐步下降，見圖3D。

2.4 GSEA結(jié)果和驗證

GSEA結(jié)果提示，這些關(guān)鍵基因涉及的生物學改變主要與糖脂代謝、活性氧、免疫炎癥等有關(guān)，見圖4。這些基因在外部數(shù)據(jù)集中大多也存在差異表達。基于外部驗證數(shù)據(jù)集GSE123342, 本研究同樣發(fā)現(xiàn)，與穩(wěn)定型冠心病相比,HP、FAM198B、CD163、FAM20A、MCEMP1表達水平在急性心肌梗死發(fā)生時顯著上升，而RNASE2僅表現(xiàn)出上升趨勢，見圖5。

3 討 論

本研究基于PBMCs轉(zhuǎn)錄組學的變化探討心肌梗死及梗死恢復過程中可能存在的生物學改變和潛在標志物, WGCNA和DEGs分析結(jié)果顯示，心肌梗死過程伴隨著免疫炎癥紊亂?；贚ASSO分析，本研究從上述2種分析方法得到的共同基因中鑒定出6個關(guān)鍵基因，即CD163、RNASE2、HP、FAM20A、MCEMP1和FAM198B。ROC曲線分析結(jié)果提示，這些關(guān)鍵基因?qū)π募」Ｋ赖陌l(fā)生均具有較高的診斷價值。本研究還發(fā)現(xiàn)，在梗死恢復過程中，這些基因的表達呈現(xiàn)明顯下降趨勢。GSEA分析結(jié)果表明，這些關(guān)鍵基因涉及的生物學改變主要與糖脂代謝、活性氧、免疫炎癥等有關(guān)。

CD163分子是Ⅰ型膜蛋白，蛋白表達限于單核細胞/巨噬細胞系, CD163抗原特異性釋放機制可能在炎癥調(diào)節(jié)過程中起重要作用。在動脈粥樣硬化過程中, CD163+巨噬細胞能促進血管生成和增加血管通透性，并伴隨炎癥反應；此外，破裂的冠狀動脈斑塊中CD163的表達增加了心肌梗死和冠心病發(fā)生風險[24]。RNASE2編碼的蛋白質(zhì)是非分泌型核糖核酸酶，屬于胰核糖核酸酶家族[25]。FAM20A是激酶編碼基因家族的成員，本身不具有激酶活性，通過與家族成員FAM20C形成復合物，增強FAM20C的激酶活性，從而使分泌通路內(nèi)的蛋白磷酸化[26]。MCEMP1又稱C19ORF59, 可編碼表達于肥大細胞、巨噬細胞等的跨膜蛋白[27]。研究[27]表明,MCEMP1基因可能是卒中診斷和預后評估的新生物標志物。生物信息學研究[28-29]提示,RNASE2、FAM20A和MCEMP1可能是心肌梗死的關(guān)鍵基因，然而這些基因在心肌梗死及梗死恢復過程中的具體作用機制有待進一步研究。HP基因編碼一種前蛋白，經(jīng)處理后產(chǎn)生α鏈和β鏈，結(jié)合為四聚體產(chǎn)生觸珠蛋白，其基因型與急性心肌梗死的發(fā)病風險密切相關(guān)，并且能決定心肌梗死面積[30-31]。FAM198B(GASK1B)是目前未知功能的新基因，可能編碼定位在高爾基體上的膜結(jié)合糖蛋白，被認為參與癌癥的轉(zhuǎn)移與進展[32]。本研究分析結(jié)果顯示,FAM198B可能是心肌梗死和反映心肌梗死恢復的關(guān)鍵基因。

本研究尚存在一定局限性: 首先，本研究基于GEO數(shù)據(jù)庫中下載的心肌梗死患者PBMCs的測序數(shù)據(jù)集進行分析，還需要進一步開展分子生物學實驗來驗證這些基因在心肌梗死中的表達及其可能的生物學機制; 其次，本研究僅分析了現(xiàn)有的樣本數(shù)據(jù)，未來還需基于更大的樣本規(guī)模和更詳細的樣本類型進一步深入研究; 最后，由于缺乏臨床資料，本研究未分析這些關(guān)鍵基因的表達與心肌梗死患者年齡、性別、射血分數(shù)等臨床指標的相關(guān)性。

綜上所述，本研究基于PBMCs確定了6個可能與心肌梗死密切相關(guān)的關(guān)鍵基因，即CD163、RNASE2、HP、FAM20A、MCEMP1和FAM198B, 為心肌梗死的診斷、評估、管理和預后分析提供了新的思路。