基于GEO 數據庫和生物信息學分析篩選小鼠心肌缺血再灌注損傷相關的潛在樞紐基因

王建茹，彭廣操，朱明軍

河南中醫藥大學第一附屬醫院心內科，鄭州 450099

急性心肌梗死（acute myocardial infarction，AMI）是冠狀動脈粥樣硬化性心臟病的一種急危重表現，發病率和死亡率均較高。近年來，隨著人們生活方式的巨大變化，AMI 病例數量急劇增加，AMI 已成為我國居民住院和死亡的主要原因［1-2］。早期快速恢復冠狀動脈血流可減少心肌梗死的面積，隨著溶栓、經皮冠狀動脈介入治療等再灌注策略的發展，AMI 的死亡率有所降低。但大量證據表明再灌注會引起心肌的繼發性損傷，其所造成的心肌損傷可占最終心肌梗死面積的50%，這種現象叫作心肌缺血再灌注損傷 （myocardium ischemia-reperfusion injury，MIRI）；MIRI 已成為AMI 經介入或溶栓治療后最常見的臨床問題［3-4］。MIRI主要表現形式多樣，可進一步加重心臟結構和功能的損傷，最終引起多種不良心血管結局［5-6］。MIRI 的發生發展過程涉及炎癥反應、氧化應激、鈣超載、冠狀微循環障礙、免疫反應、細胞死亡、長鏈非編碼RNAs 表達水平變化等，但其詳細的機制仍不清楚［7-9］。因此，積極探索MIRI 的分子機制，挖掘關鍵基因，對于MIRI的防治至關重要。

近年來，利用生物信息學技術對芯片基因表達數據進行挖掘已廣泛地用于分析與疾病相關的差異表達基因（differentially expressed genes，DEGs），尋找關鍵基因，篩選與疾病診斷、治療和預后相關的生物標志物等方面。DEGs有可能為探索MIRI病理過程的發展及防治提供有價值的線索。現階段，研究人員可以從基因表達數據庫（gene expression omnibus，GEO）等渠道獲取芯片數據。但值得注意的是，在篩選DEGs 的過程中由于單項獨立研究有限的樣本數、樣本的異質性、技術平臺的差異、數據異常值及分析處理方法不恰當等因素的存在，研究結果可能出現一定的偏倚。目前，基于GEO 中基因表達芯片鑒定DEGs存在一些不足之處；因此，尋找更加有效的計算方法，可避免偏倚結果的出現，得到更加穩健的結果。穩健排序整合（robust rank aggregation，RRA）方法是一種可使多個數據集之間的偏差和誤差達到最小的標準方法，用于各種排序基因列表的比較，可有效避免上述問題［10］。RRA 基于每個基因是隨機排序的假設，如果一個基因在所有實驗中排名高，用RRA 計算所得P值越小，則其成為DEG 的可能性越大。RRA 已廣泛應用于多種腫瘤和非腫瘤疾病的多數據集的綜合分析，并表現出了良好的效果。目前，尚沒有關于RRA整合多個數據集分析MIRI的報道。

綜上所述，本研究擬從GEO 數據庫中獲取小鼠GSE61592、GSE83472 和GSE160516 數據集，利用limma 包篩選各數據集中的DEGs，然后采用RRA整合分析3 個數據集的DEGs 篩選穩健DEGs；一方面利用clusterProfiler 包對穩健DEGs 進行基因本體論（gene ontology，GO）和京都基因與基因組百科全書（Kyoto Encyclopedia of Genes and Genomes，KEGG）通路富集分析，另一方面利用String 數據庫構建穩健DEGs 的蛋白質-蛋白質相互作用（proteinprotein interaction，PPI） 網絡；利用MCODE 和cytoHubba 插件來篩選PPI 網絡中的子模塊和hub 基因，再利用clusterProfiler 包對最重要的子模塊基因和hub 基因進行GO 和KEGG 分析；最后，采用反轉錄- 定量聚合酶鏈反應（reverse transcriptionquantitative polymerase chain reaction，RT-qPCR） 技術檢測小鼠MIRI 模型中hub 基因的表達情況。本研究旨在進一步探索MIRI 病理過程中潛在的干預靶點和分子機制，為其治療提供新的策略和思路。

1 材料與方法

1.1 芯片數據的來源

本研究流程如圖1 所示。從GEO 數據庫中下載與小鼠MIRI相關的基因芯片數據集GSE61592［平臺ID 為GPL6887，假手術（sham）組3 例，MIRI 組3例］、GSE83472（平臺ID 為GPL6885，sham 組4 例，MIRI 組4 例）和GSE160516（平臺ID 為GPL23038，sham 組4 例，MIRI 組12 例），其中GSE83472 的樣本量信息是根據本次研究的需求提取的樣本數，3 個數據集檢測樣本均為心肌組織。

圖1 研究流程圖Fig 1 Flow chart of the study

1.2 數據預處理和DEGs的篩選

利用perl腳本對表達數據進行注釋，將探針矩陣轉化為基因矩陣；利用impute 包對缺失值進行填充；利用avereps 函數對基因對應多個探針取均值；利用normalizeBetweenArrays 函數對數據進行矯正；若基因表達數值較大，則需轉換為以2 為底的對數。利用limma 包，以|log2差異倍數（fold change，FC）|＞1 和P＜0.05 為閾值對3 個數據集的表達數據進行差異分析，篩選DEGs。

1.3 RRA法整合分析芯片數據

利用RRA 包對獲得的DEGs 整合分析，按照|log2FC|＞1 和P＜0.05 為條件篩選穩健DEGs；然后利用pheatmap 包分別將P值排序前20 的上調和下調基因進行可視化。

1.4 穩健DEGs的功能富集分析

為了探究穩健DEGs 可能的功能作用，利用clusterProfiler包對其進行GO和KEGG富集分析。

1.5 PPI網絡的構建

將穩健DEGs 上傳至STRING 數據庫，相互作用得分設置為最高置信度≥0.400，進行PPI 分析。利用Cytoscape3.7.2軟件對PPI網絡進行可視化。

1.6 模塊分析和hub基因的篩選

將PPI 信息導入Cytoscape3.7.2 軟件中，利用MCODE 插件并以默認參數（度截止=2，節點評分截止=0.2，最大深度=100，k 評分=2）為標準對PPI 網絡進行模塊分析，篩選最重要的子模塊。CytoHubba插件可以通過其包含的12種算法分別對PPI網絡中的每個節點（基因）賦予一定的分數，然后再依據基因得分將排序靠前的基因作為hub 基因。本研究的hub基因由12 種算法排序前100 的基因的重疊基因組成。然后，再利用clusterProfiler包對最重要的子模塊基因和hub基因進行GO和KEGG富集分析。

1.7 Hub基因的驗證

1.7.1 MIRI小鼠模型構建 18 只6～8 周齡的C57BL/6 J雄性小鼠購于鄭州市惠濟區華興實驗動物養殖場，許可證號：SCXK-（豫）-2019-0002。實驗動物飼養于獨立通風籠盒動物房，自由飲食和飲水。本實驗已經河南中醫藥大學第一附屬醫院實驗動物福利倫理審查委員會批準（批準編號：YFYDW2020004）。小鼠被隨機分為sham 組和MIRI 組，每組9 只。戊巴比妥鈉（40 mg/kg）腹腔注射麻醉并固定小鼠后，手術區備皮，消毒，于左側第4、5 肋間切口，打開胸腔，撕開心包，輕壓胸廓擠出心臟；探查確定左心耳，于肺動脈圓錐與左心耳交界距離左心耳根部約3 mm 處結扎冠狀動脈左前降支，缺血30 min 后去除結扎線，關閉切口，縫合皮膚，24 h 后處死小鼠（即再灌注24 h）。Sham 組小鼠操作同MIRI 組小鼠，但不結扎冠狀動脈左前降支。

1.7.2 氯化三苯基四氮唑染色 再灌注24 h 后，迅速摘取心臟，用磷酸鹽緩沖液洗滌，從結扎線附近垂直于心臟長軸進行切片，將心臟平均切成5 片，然后置入2% 氯化三苯基四氮唑（triphenyl tetrazolium chloride，TTC）染色液（北京索萊寶科技有限公司）中37 ℃水浴15 min，結束后固定并拍照。TTC 染色后正常心肌組織呈紅色，梗死組織呈蒼白色。

1.7.3 RT-qPCR 各組小鼠取適量左心室組織（其中MIRI 組為心肌缺血危險區組織），加入TRIZOL（Invitrogen 公司，美國）進行裂解，振蕩混勻，冰上放置5 min；加入氯仿，振蕩靜置離心；將離心管內的上清液轉移至另一干凈的EP 管中，加入異丙醇，振蕩靜置離心；去上清液，加入75%乙醇洗滌RNA，離心后吸掉多余液體，自然控干RNA；用無RNase水溶解RNA，輕輕吹打助溶，4 ℃靜置10 min；用核酸蛋白濃度微量測定儀（Thermo Fisher 公司，美國）測RNA 濃度。按照反轉錄試劑盒（TOYOBO 公司，日本）說明書合成cDNA；按照2×Universal SYBR Green Fast qPCR Mix試劑（ABclonal公司，中國）的說明書配制反應液，在TL-988 熒光定量PCR 儀（西安天隆科技有限公司）上進行RT-qPCR；應用2-ΔΔCT對基因進行相對定量分析。引物由蘇州金唯智生物科技有限公司合成（序列見表1）。PCR 反應條件：95 ℃3 min預變性，95 ℃5 s變性，60 ℃34 s退火延伸，共45 個循環。采用SPSS 21.0 軟件對數據進行處理分析，數據以±s表示，2組間比較采用t檢驗，P＜0.05為差異有統計學意義。

表1 RT-qPCR引物序列Tab 1 Primer sequences for RT-qPCR

Continued Tab

2 結果

2.1 GEO數據集預處理和DEGs的篩選

差異分析結果顯示，GSE61592 篩選出1 684 個DEGs，其中上調的832個，下調的852個；GSE83472篩選出55 個DEGs，其中上調的51 個，下調的4 個；GSE160516篩選出1 061個DEGs，其中上調的733個，下調的328個。圖2為3個數據集的火山圖。

圖2 各數據集的火山圖Fig 2 Volcano plots for each data set

2.2 RRA法篩選穩健DEGs

RRA 方法對3 個數據集進行了整合分析后，共鑒定出294 個穩健DEGs，其中上調的209 個，下調的85 個。圖3 為P值排序前20 個上調和下調穩健DEGs 的熱圖。

圖3 前20個上調和下調穩健DEGs的熱圖Fig 3 Heat map of the top 20 up-and down-regulated robust DEGs

2.3 穩健DEGs的功能富集分析

利用clusterProfiler 包、org.Mm.eg.db 包將篩選出的294 個穩健DEGs 進行功能富集分析，如圖4所示。GO 分析結果顯示，依據P＜0.05 確定了949個GO 條目；生物學過程（biological process，BP）條目為867 條，主要涉及細胞的遷移、趨化作用、對外界刺激反應的正性調節、調節腫瘤壞死因子超家族細胞因子的產生、細胞因子介導的信號通路等；分子功能（molecular function，MF）條目為62 條，主要涉及受體配體活性、細胞因子受體結合、細胞因子活性、趨化因子活性、細胞因子受體活性、Toll 樣受體結合等；細胞組分條目（cellular component，CC）為20 條，主要涉及細胞外基質、線粒體外膜、肌質網、溶酶體等。KEGG分析結果顯示，依據P＜0.05 共映射出5 條信號通路，包括細胞因子受體相互作用、病毒蛋白與細胞因子及細胞因子受體的相互作用、阿米巴病、白介素-17 和核因子κB （nuclear factor kappa B，NF-κB）信號通路。

圖4 穩健DEGs的功能富集分析Fig 4 Functional enrichment analysis of robust DEGs

2.4 穩健DEGs的PPI網絡構建和模塊分析

利用STRING 數據庫對篩選得到的穩健DEGs 構建PPI網絡（圖5A）。該網絡包括1 248 條邊和239 個節點，其中上調的174個，下調的65個。模塊分析共篩選出14 個子模塊，其中模塊1 最重要（圖5B）。模塊1 中的穩健DEGs 主要富集在急性炎癥反應、炎性細胞的遷移、細胞趨化、細胞黏附分子結合、趨化因子受體活性、病毒蛋白與細胞因子及細胞因子受體的相互作用、NF-κB 信號通路、趨化因子信號通路等（圖5C、D）。

圖5 穩健DEGs的PPI網絡構建和模塊分析Fig 5 PPI network construction and modular analysis of robust DEGs

2.5 Hub基因的篩選

共獲得17 個hub 基因（圖6A、B），相關信息詳見表2。功能富集分析的結果顯示，17 個hub 基因主要富集在白細胞和中心粒細胞的遷移、淋巴細胞趨化性、對外界刺激反應的正性調節、趨化因子的活性、趨化因子受體結合、細胞因子受體結合等功能；涉及病毒蛋白與細胞因子及細胞因子受體的相互作用、細胞因子受體相互作用、趨化因子信號通路、NF-κB信號通路等通路（圖6C、D）。

表2 Hub基因的相關信息Tab 2 Information about the hub genes

圖6 Hub基因的篩選及其功能富集分析Fig 6 Screening of hub genes and analysis of their functional enrichment

2.6 Hub基因的驗證

由圖7A 可知，與sham 組相比，MIRI 組小鼠經過缺血再灌注處理后，心肌明顯出現了梗死，說明MIRI 模型構建成功。如圖7B 所示，與sham 組相比，MIRI 組小鼠心肌趨化因子（C-C 基序） 配體4［chemokine（C-C motif） ligand 4，Ccl4］、Ccl6、Ccl7、趨化因子（C-X-C 基序）受體4［chemokine（C-X-C motif） receptor 4，Cxcr4］、趨化因子（C-C基序） 受體2 ［chemokine（C-C motif）receptor 2，Ccr2］、信號調節蛋白β1（signal-regulatory protein β 1，Sirpb1）、低親和力免疫球蛋白γ Fc 區受體Ⅱb（low affinity immunoglobulin gamma Fc region receptor Ⅱb，Fcgr2b）、白細胞表面抗原Cd53（leukocyte surface antigen CD53，Cd53）、花生四烯酸5-脂氧合酶激活蛋白（arachidonate 5-lipoxygenase activating protein，Alox5ap）、髓樣分化初級反應基因88（myeloid differentiation primary response gene 88，Myd88）、巨噬細胞清道夫受體1 （macrophage scavenger receptor 1，Msr1）、 基質金屬肽酶14（matrix metallopeptidase 14，Mmp14）、髓樣細胞上表達的觸發受體2 （triggering receptor expressed on myeloid cells 2，Trem2）、樁蛋白（leupaxin，Lpxn）的mRNA 表達量均上調，而低親和力免疫球蛋白γ Fc 區受體Ⅲ（low affinity immunoglobulin gamma Fc region receptor Ⅲ，Fcgr3）、補體C1q亞組分亞單位B（complement C1q subcomponent subunit B，C1qb）、去整合素和含金屬蛋白酶結構域蛋白（a disintegrin and metalloproteinase domain-containing protein 8，Adam8）的mRNA表達未見明顯的差異。回顧文獻發現，目前僅Fcgr3、C1qb、Trem2、Lpxn、Cd53、Alox5ap、Sirpb1、Fcgr2b在MIRI 與正常組織中的差異表達未被報道。

圖7 Hub基因的驗證Fig 7 Validation of the hub genes

3 結論

本研究結合微陣列基因表達譜和生物信息學方法，探索了MIRI 病理過程中可能的新的核心基因及其潛在的分子機制。研究首先利用limma 包對小鼠GSE61592、GSE83472 和GSE160516 數據集中的DEGs 進行了篩選。其次，利用RRA 法對DEGs 進行合并分析，共鑒定出209 個上調和85 個下調穩健DEGs。然后，對穩健DEGs的PPI網絡進行功能模塊分析，共發現14 個子模塊，其中模塊1 最重要；同時，鑒定出了17個hub基因，而且這些基因大部分都在模塊1 中。最后，本文對穩健DEGs、模塊1 中的基因和17個hub基因分別進行了功能富集分析，結果發現這些基因大都是通過調控炎性細胞的遷移、趨化因子和細胞因子及其受體活性、細胞因子受體相互作用、趨化因子信號通路、Toll 樣受體信號通路、NFκB 信號通路等生物學功能和通路，來介導MIRI的病理過程。眾所周知，MIRI 是一個復雜的事件，涉及炎癥細胞（如中性粒細胞、單核/巨噬細胞、淋巴細胞） 和關鍵激活信號（如細胞因子、趨化因子）等［11］。研究報道，炎癥信號通路Toll 樣受體和NFκB 信號通路以及兩者之間的交互作用可加劇MIRI的嚴重程度［12-14］。此外，MIRI可強化趨化因子的表達，這些趨化因子可誘導炎性細胞的活化和浸潤，進而進一步產生強烈的炎癥反應，而抑制趨化因子可促進損傷心肌的愈合［11，15］。以上相關文獻的報道從側面證實了本研究富集分析結果的科學性。

為了進一步確證生物信息學分析預測的結果，本研究構建了小鼠MIRI 模型并用RT-qPCR 方法檢測了17 個hub 基因mRNA 的表達情況。結果顯示，在MIRI 過程中Ccl7、Sirpb1、Fcgr2b、Cxcr4、Cd53、Alox5ap、Ccl6、Myd88、Msr1、Ccl4、Mmp14、Ccr2、Trem2、Lpxn的表達上調，而Fcgr3、C1qb、Adam8未見明顯差異。既往其他研究也證實了在MIRI 的過程中Ccl7、Cxcr4、Ccl6、Myd88、Ccl4、Msr1、Ccr2、Mmp14的表達上調［16-22］，這與本研究的結果吻合，從側面佐證了本研究結果的可靠性和科學性。Adam8參與了MIRI 的病理過程，但在再灌注0.5 h 和24 h 時其表達并沒有增加，48 h 后才明顯增加［23］。目前，Fcgr3、C1qb、Trem2、Lpxn、Cd53、Alox5ap、Sirpb1、Fcgr2b在MIRI中的差異表達情況，均未被其他研究報道；但Fcgr3和C1qb在本研究中的生物信息學分析預測結果與動物實驗的驗證結果不一致。綜上所述，Trem2、Lpxn、Cd53、Alox5ap、Sirpb1、Fcgr2b這6 個hub 基因的高表達在MIRI 中發揮著重要作用，將有助于闡釋MIRI可能的分子機制。

Alox5ap位于染色體13q12-13，與心肌梗死、動脈粥樣硬化等心血管疾病的易感性密切相關，其編碼的5-脂氧合酶激活蛋白，可將花生四烯酸催化為白三烯，進而參與炎癥反應介導多種心血管疾病［24］。Fcgr2b作為唯一的抑制性Fc 受體，在調控炎癥和免疫方面發揮著重要作用，與自身免疫性疾病的易感性密切相關［25］。Fcgr2b廣泛表達于巨噬細胞、樹突狀細胞等免疫細胞，可抑制激活FcγRs誘導的吞噬前信號，進而抑制吞噬作用［25］。研究［26］顯示，高效的吞噬作用可以及時清除MIRI過程中凋亡的心肌細胞，否則，將引起炎癥反應、細胞進一步凋亡等后續損傷。因此，Fcgr2b可能通過調控吞噬作用來參與MIRI的病理過程。Trem2為免疫球蛋白超家族受體的成員，在小膠質細胞、巨噬細胞、樹突狀細胞、髓系細胞等細胞中表達，可與相對分子質量為12 000的跨膜免疫接頭DNAX 激活蛋白（DNAX-activating protein of 12 kDa，DAP12）的酷氨酸活化基序結合，通過激活下游信號通路，在免疫炎癥等反應中發揮作用［27-30］。Trem2已被報道通過調控細胞焦亡及凋亡、炎癥反應等，參與了腦、肝、腎、肺缺血再灌注損傷的病理過程［27-30］。由于各器官缺血再灌注存在相似的病理生理過程，因此Trem2可能也介導了MIRI。研究［31-32］顯示，B 細胞被募集到了MIRI 的心肌組織，其加劇MIRI 的機制可能與其分泌IgM 進而激活補體系統和炎性細胞有關；同時，被募集到MIRI 心肌組織的還有T 細胞尤其是CD4+T 細胞，其主要表現為Th1 促炎癥表型，通過釋放γ 干擾素等，加劇MIRI。Cd53為四次穿膜蛋白（tetraspanin）超家族的成員，高表達于T 細胞、B 細胞等免疫細胞表面，在免疫系統中發揮著重要作用，主要表現為調控免疫細胞的黏附、遷移、增殖、活化以及信號轉導等［33］。Lpxn是樁蛋白（paxillin）家族的成員，表達于白細胞細胞系、B 細胞等細胞，可調控整合素的活性，在B 細胞免疫反應中發揮著不可或缺的作用［34］。Sirpb1是免疫球蛋白超家族的信號調節蛋白成員，胞內段缺乏信號轉導基序，但可通過與DAP12結合，轉導活化信號，參與細胞增殖、分化、免疫功能等［35］。值得注意的是，Sirpb1與B 細胞受體通路相關，參與了B 細胞的增殖［36］。也就是說，Cd53、Lpxn、Sirpb1可能通過調控免疫機制來介導MIRI 的病理過程。綜上所述，我們尚未發現有關這6 個hub 基因在MIRI 中作用的相關報道，但回顧文獻的結果提示這些hub 基因表現出了介導MIRI 病理過程的可能性，為探索它們在MIRI 中的作用機制提供了一定的依據。

本研究運用生物信息學的方法發現了在小鼠MIRI 過程中的一些hub 基因，但仍有一些局限性。雖運用RRA 法分析了多個數據集，但由于基因芯片數據的可獲得性有限以及每個數據集中滿足研究條件的樣本較少，最終導致本研究納入的樣本量有限。生物信息學預測結果雖進行了實驗驗證，但由于種屬的差異，結果仍需要進一步在臨床實踐過程中進行佐證。總之，本研究利用基因表達譜和生物信息學相結合的方法共挖掘出17個hub基因，通過回顧既往文獻和實驗驗證后共鑒定出6個潛在的hub基因（Fcgr2b、Alox5ap、Trem2、Sirpb1、Cd53、Lpxn），研究結果可為進一步探討MIRI 的分子機制和治療靶點提供新的思路和切入點。