基于生物信息學(xué)途徑篩選缺血性腦卒中關(guān)鍵基因及藥物預(yù)測(cè)

于 瑜，王鐘興

（中山大學(xué)附屬第一醫(yī)院麻醉科，廣東廣州 510080）

腦卒中作為全球第二大死亡原因，嚴(yán)重危害人類(lèi)健康。全球每年約有550 萬(wàn)人死于腦卒中，每4個(gè)成年人中就有1人發(fā)?。?］。發(fā)病后遺留的身體殘疾和神經(jīng)功能障礙給中風(fēng)幸存者的家庭及社會(huì)造成極大負(fù)擔(dān)。在腦卒中病例中，缺血性腦卒中（以下簡(jiǎn)稱腦缺血）約占70%，是腦卒中的主要類(lèi)型［1］。目前，臨床上腦缺血的主要治療方法包括腦血管再通和神經(jīng)保護(hù)。遺憾的是，這些措施收益有限。為了找到更有效的治療方式以及治療靶點(diǎn)，許多科學(xué)家從遺傳學(xué)角度出發(fā)，尋找與腦缺血密切相關(guān)的基因，希望為治療提供更多選擇。例如：有研究人員利用高通量組學(xué)技術(shù)發(fā)現(xiàn)CLDN20、GADD45G、RGS2、BAG5 和CTNND2 可以作為小鼠和人類(lèi)血液樣本中腦缺血的血液生物標(biāo)記物，協(xié)助腦缺血的臨床診斷［2］。MAP2K4（絲裂原活化蛋白激酶4）是MAPK 通路的重要樞紐，MAP2K4 基因多態(tài)性rs3826392 可能在腦缺血后的炎癥反應(yīng)中發(fā)揮調(diào)節(jié)作用［3］。這些從基因?qū)用娼庾x腦缺血病理生理過(guò)程的研究對(duì)治療方案的研發(fā)有一定啟迪作用。因此，本研究借助生物信息學(xué)手段，通過(guò)對(duì)公共數(shù)據(jù)庫(kù)的分析，尋找對(duì)缺血性腦卒中有重要作用的基因及相關(guān)通路，希望找到潛在治療藥物，能為缺血性腦卒中的進(jìn)一步研究和后續(xù)治療提供線索。

1 材料與方法

1.1 數(shù)據(jù)收集

大腦中動(dòng)脈栓塞（middle cerebral artery occlu?sion，MCAO）是研究缺血性腦卒中的常用動(dòng)物模型，可精確控制缺血時(shí)長(zhǎng)。本研究選擇GEO 數(shù)據(jù)庫(kù)（https：//www.ncbi.nlm.nih.gov/geo/）中 研 究MCAO 小鼠表達(dá)譜變化的芯片集GSE98319。該數(shù)據(jù)集包含MCAO 小鼠和假手術(shù)組小鼠大腦皮層的檢測(cè)結(jié)果。每組各3個(gè)樣本。

1.2 差異表達(dá)分析

用R 語(yǔ)言的ggplot 包繪制6 個(gè)樣本總體基因表達(dá)箱線圖，并對(duì)芯片集進(jìn)行多維度分析（multidi?mensional scaling，MDS）。Limma 包用于差異表達(dá)分析，根據(jù)測(cè)序平臺(tái)注釋信息及Ensemble 數(shù)據(jù)庫(kù)對(duì)基因類(lèi)型的注釋信息，將探針I(yè)D 轉(zhuǎn)換為相應(yīng)的基因名稱（gene symbol）并刪除LncRNA 的數(shù)據(jù)。選擇P<0.05 且|log2FC|>0.6 作為差異表達(dá)基因（dif?ferential expressed genes，DEG）篩選標(biāo)準(zhǔn)。差異分析結(jié)果以火山圖形式呈現(xiàn)。用R 語(yǔ)言繪制熱圖展示數(shù)據(jù)集的聚類(lèi)情況。

1.3 蛋白質(zhì)互作網(wǎng)絡(luò)構(gòu)建及核心基因篩選

利用STRING（http：//www.string-db.org/）在線數(shù)據(jù)庫(kù)對(duì)DEG 進(jìn)行分析，預(yù)測(cè)可能在腦缺血發(fā)病過(guò)程中發(fā)揮重要作用的基因所編碼蛋白質(zhì)之間的互作關(guān)系。選擇中度置信度（作用關(guān)系綜合得分0.4-0.7）作為顯著標(biāo)準(zhǔn)。將STRING 分析結(jié)果導(dǎo)入Cytoscape 3.7.1.軟件進(jìn)行可視化，使用Cytohubba和MCODE 插件對(duì)PPIN 進(jìn)行分析。通過(guò)Cytohubba 計(jì)算出每個(gè)節(jié)點(diǎn)的度（Degree）［4］，Degree得分前10的基因構(gòu)成的網(wǎng)絡(luò)可視作核心子網(wǎng)絡(luò)。通過(guò)MCODE（molecular complex detection）從蛋白質(zhì)互作網(wǎng)絡(luò)（protein-protein interaction network，PPIN）中篩出得分最高的3個(gè)功能模塊亦視作核心子網(wǎng)絡(luò)［5］。為了進(jìn)一步縮小核心基因范圍，將MCODE 評(píng)分最高模塊的基因及Degree 得分前10 的基因進(jìn)行整合，他們與DEG 中|log2FC|排名前20基因的交集即為核心基因（hub genes）。

1.4 小鼠大腦中動(dòng)脈梗塞模型建立及核心基因表達(dá)水平測(cè)定

1.4.1 大腦中動(dòng)脈栓塞模型建立及實(shí)驗(yàn)分組 3只6～8 周SPF 級(jí)雄性C57B/6J 小鼠購(gòu)自湖南斯萊克景達(dá)實(shí)驗(yàn)動(dòng)物有限公司，小鼠體質(zhì)量為18～25 g，許可證號(hào)：SCXK（湘）2019-0004。本實(shí)驗(yàn)經(jīng)中山大學(xué)實(shí)驗(yàn)動(dòng)物管理與使用委員會(huì)同意。3%異氟烷誘導(dǎo)后將小鼠仰臥固定于操作臺(tái)并戴上麻醉面罩，術(shù)中以1.5%異氟烷維持麻醉。頸部備皮消毒后頸正中切口，小心暴露右側(cè)頸動(dòng)脈鞘并將右側(cè)頸總、頸外及頸內(nèi)動(dòng)脈分離出來(lái)，在三角分叉處遠(yuǎn)端結(jié)扎頸外動(dòng)脈，近端穿一根線打活結(jié)，用動(dòng)脈夾夾閉頸總動(dòng)脈及頸內(nèi)動(dòng)脈，在頸外動(dòng)脈中間剪一小口，將線栓插入頸外動(dòng)脈后松開(kāi)頸內(nèi)動(dòng)脈夾，將線栓緩慢沿著頸內(nèi)動(dòng)脈插入，線栓上的黑色標(biāo)記點(diǎn)到達(dá)三角分叉處（插入約1 cm 左右），有明顯阻力時(shí)停止。固定線栓，松開(kāi)頸總動(dòng)脈夾，縫合傷口。將小鼠放回籠內(nèi)，側(cè)臥，插栓側(cè)朝上，維持體溫。缺血30 min后經(jīng)原切口拔出線栓，恢復(fù)血供，再灌注12 h 后留取小鼠雙側(cè)大腦皮質(zhì)。小鼠缺血側(cè)皮層為MCAO 組，對(duì)側(cè)為對(duì)照組（CTRL組）。

1.4.2 標(biāo)本收集以及qRT-PCR 驗(yàn)證基因表達(dá)變化 按照奕杉生物組織RNA 提取試劑盒說(shuō)明書(shū)提取小鼠雙側(cè)皮層RNA，用翊圣生物Hifair? Ⅲ1st Strand cDNA Synthesis SuperMix for qPCR 試劑盒將RNA 逆轉(zhuǎn)錄成cDNA 后按照Hieff? qPCR SYBR?Green Master Mix 配置出實(shí)時(shí)熒光定量PCR（quanti?tative real-time polymerase chain reactionq，qRTPCR）體系進(jìn)行檢測(cè)。以Actin 為內(nèi)參，用2-ΔΔCT相對(duì)定量法分析待測(cè)基因的轉(zhuǎn)錄水平。所用引物序列見(jiàn)表1。

1.4.3 統(tǒng)計(jì)方法 采用SPSS 17.0 統(tǒng)計(jì)軟件進(jìn)行統(tǒng)計(jì)分析。RNA 相對(duì)表達(dá)量為連續(xù)性資料，以均數(shù)±標(biāo)準(zhǔn)差（±s）表示，組間比較采用t檢驗(yàn)（雙側(cè)），P<0.05認(rèn)為差異有統(tǒng)計(jì)學(xué)意義。

1.5 GSEA算法進(jìn)行功能富集分析

本研究選擇GSEA C2 數(shù)據(jù)集中的KEGG 基因集（c2.cp.kegg.v7.2.）和C5 數(shù)據(jù)集中的GO 基因集（c5.go.mf.v7.2.entrez、c5.go.bp.v7.2.entrez、c5.go.cc.v7.2.entrez）為預(yù)設(shè)基因集，利用R 語(yǔ)言對(duì)差異分析結(jié)果進(jìn)行GSEA 分析。FDR（false discovery rate）<0.05 作為富集通路的篩選標(biāo)準(zhǔn)。富集結(jié)果圖由R 語(yǔ) 言的clusterprofiler 包［6］和GSEA 函數(shù) 進(jìn)行繪制。

1.6 缺血性腦卒中藥物預(yù)測(cè)

Connectivity Map（CMap）數(shù)據(jù)庫(kù)包含1 309 種小分子化合物作用于多種細(xì)胞系引起的基因表達(dá)譜變化數(shù)據(jù)，可用于比較藥物誘導(dǎo)的基因圖譜與基因表達(dá)之間的相似性，并得到一個(gè)從-100 到100 的連通性評(píng)分：該評(píng)分大于0，表示化合物引起和上傳基因相似的變化；評(píng)分小于0，表示該化合物引起和上傳基因相反的變化，即該化合物可能對(duì)疾病有治療作用［7］。由于CMap 將上傳基因數(shù)目限制在150 個(gè)，本文將|log2FC|前150 上調(diào)DEG和全部下調(diào)DEG（100 個(gè)）分別上傳至Cmap 在線數(shù)據(jù)庫(kù)（https：//clue.io/）進(jìn)行藥物預(yù)測(cè)。連通性評(píng)分<-80 的小分子化合物作為本次預(yù)測(cè)的結(jié)果。

2 結(jié)果

2.1 缺血性腦卒中差異表達(dá)基因的篩選

GSE98319 數(shù)據(jù)集的箱線圖提示各樣本基因表達(dá)水平均衡（圖1A），具有可比性。MDS（圖1B）分析結(jié)果提示MCAO 組和假手術(shù)組能進(jìn)行有效區(qū)分。圖1C 顯示GSE98319 芯片集差異表達(dá)分析的結(jié)果；不同變化方向的基因表達(dá)基因用不同顏色進(jìn)行展示。將這些差異基因中l(wèi)ncRNA 刪去，根據(jù)篩選條件，篩出521 個(gè)差異mRNA，其中421 個(gè)上調(diào)、100 個(gè)下調(diào)。圖1D 為該數(shù)據(jù)集的熱圖，該數(shù)據(jù)集樣本聚類(lèi)效果較好，可信度較高。

2.2 缺血性腦卒中蛋白互作網(wǎng)絡(luò)的構(gòu)建以及核心基因的篩選

通過(guò)STRING預(yù)測(cè)521個(gè)DEG所編碼蛋白質(zhì)之間的相互作用關(guān)系。將預(yù)測(cè)結(jié)果導(dǎo)入Cytoscape 軟件構(gòu)建出蛋白質(zhì)互作網(wǎng)絡(luò)（圖2）。此外，用Cyto?scape 的插 件CytoHubba 和MCODE 篩選出PPIN 中的核心子網(wǎng)絡(luò)：度（Degree）排名前10的基因構(gòu)成的子網(wǎng)絡(luò)（圖3A）以及MCODE 評(píng)分最高的子網(wǎng)絡(luò)（圖3B）。為了進(jìn)一步縮小核心基因的范圍，Degree 前10 的基因及構(gòu)成MCODE 得分最高子網(wǎng)絡(luò)的基因與DEG中上調(diào)倍數(shù)最高的20個(gè)基因求交集（圖4），篩選出Drd4和Hcar22 個(gè)核心基因，表2 為他們的差異分析結(jié)果。

表1 qRT-PCR使用引物序列Table 1 Primers used for qRT-PCR

圖1 GSE98319芯片數(shù)據(jù)集箱線圖、MDS圖、差異表達(dá)基因火山圖以及熱圖Fig.1 Plots of GSE98319 gene chips

2.3 實(shí)時(shí)熒光定量PCR 驗(yàn)證核心基因

MCAO 組與CTRL 組各3 個(gè)樣本，經(jīng)統(tǒng)計(jì)分析發(fā)現(xiàn)腦缺血后Hcar2基因的相對(duì)表達(dá)量為（6.01±0.54）與對(duì)照組（1.00±0.09）相比有所上升，差異有統(tǒng)計(jì)學(xué)意義（t=15.79，P=0.000 094）。Drd4在腦缺血再灌注損傷組的相對(duì)表達(dá)量為（4.96±0.29）與對(duì)照組（1.00±0.02）相比，差異有統(tǒng)計(jì)學(xué)意義（t=23.74，P=0.000 019；圖4）。qRT-PCR 結(jié)果與前述芯片差異分析結(jié)果一致。

表2 芯片數(shù)據(jù)核心基因差異分析結(jié)果Table 2 Results of analysis of microarray data for differential expression of hub genes

圖2 差異表達(dá)基因的蛋白質(zhì)互作網(wǎng)絡(luò)Fig.2 PPIN of DEGs

圖3 核心子網(wǎng)絡(luò)Fig.3 Sub-networks of PPIN

圖4 Hcar2 and Drd4 qRT-pcr驗(yàn)證結(jié)果Fig.4 Hcar2 and Drd4 expression level by qRT-PCR validation

2.4 GSEA功能富集分析結(jié)果

根據(jù)篩選條件，本次分析共富集到GO 分析中的12 個(gè)細(xì)胞組分（cellular component，CC）、16 種分子功能（molecular function，MF）和158 個(gè)生物過(guò)程（biological process，BP）以及11條KEGG 通路。圖5和圖6 分別顯示了每一類(lèi)富集結(jié)果中P值最小和富集分?jǐn)?shù)最高的10條富集通路。

KEGG 結(jié)果顯示細(xì)胞因子-細(xì)胞因子受體的相互作用、Toll 樣受體信號(hào)通路、抗原加工呈遞、胞質(zhì)DNA 傳感通路、RigⅠ樣受體信號(hào)、MAPK 信號(hào)通路、自噬、自然殺傷細(xì)胞活性、NOD 樣受體信號(hào)通路等與炎癥反應(yīng)密切相關(guān)的通路被富集，提示腦缺血再灌注損傷過(guò)程中炎癥反應(yīng)發(fā)揮了重要作用。GO富集分析得到了類(lèi)似結(jié)果：對(duì)病毒的防御反應(yīng)、病毒生命周期的變化、骨髓白細(xì)胞介導(dǎo)的免疫反應(yīng)、對(duì)γ 干擾素和Ⅰ型干擾素的反應(yīng)、對(duì)外部刺激反應(yīng)的正向調(diào)節(jié)和白細(xì)胞的增殖以及細(xì)胞因子產(chǎn)生的正向調(diào)控相關(guān)基因在腦缺血后被顯著富集，且全部處于激活狀態(tài)，可能在腦缺血再灌注損傷中發(fā)揮重要作用；圖5C、6C 顯示，整合子復(fù)合體、高爾基體順面膜囊、血液微粒、細(xì)胞外基質(zhì)、分泌顆粒膜、質(zhì)膜的外側(cè)面、富含M-纖維膠凝蛋白的顆粒體、構(gòu)成細(xì)胞外基質(zhì)的膠原蛋白、嗜天青顆粒和雙鏈斷裂部位可能是腦缺血再灌注損傷中重要的細(xì)胞組分。圖5D、6D 提示，腦缺血后腦組織在Ⅰ型干擾素受體結(jié)合、細(xì)胞因子活性、氣味結(jié)合、趨化因子的活性、G 蛋白偶聯(lián)受體結(jié)合、生長(zhǎng)因子結(jié)合、清道夫受體的活性、生長(zhǎng)因子的活性、絲氨酸水解酶活性、細(xì)胞因子受體結(jié)合等生物學(xué)功能方面發(fā)生了顯著變化。在上述所有富集到的通路中，除高爾基體順面膜囊的聚集為抑制外，其余通路均為上調(diào)或激活。

圖5 KEGG、GO功能富集分析結(jié)果圖Fig.5 Functional enrichment analysis results by GSEA algorithm

圖6 GSEA功能富集結(jié)果Fig.6 Functional enrichment analysis results by GSEA algorithm

2.5 CMap藥物預(yù)測(cè)結(jié)果

從CMap 下載預(yù)測(cè)藥物結(jié)果，基于對(duì)藥物的連通性評(píng)分進(jìn)行排序和篩選。連通性評(píng)分

表3 CMap預(yù)測(cè)結(jié)果Table 3 Prediction Results from CMap

圖7 4種預(yù)測(cè)藥物的簡(jiǎn)化結(jié)構(gòu)式Fig.7 Structure of 4 potential drugs

3 討論

本研究基于生物信息學(xué)途徑篩選了腦缺血關(guān)鍵基因并預(yù)測(cè)了潛在藥物。對(duì)結(jié)果進(jìn)行分析發(fā)現(xiàn)構(gòu)成核心子網(wǎng)絡(luò)的基因大多與炎癥反應(yīng)相關(guān)，他們包括Toll 樣受體（toll like receptor，TLR）家族（Tlr2、Tlr4）的成員、趨化因子CC-chemokine ligand（CCL）及相應(yīng)受體CC-chemokine receptor（CCR）家族（Ccl5、Ccl2、Cxcl16、Ccr7），白細(xì)胞分化抗原Cd86、Cd14和Cd68等。其中Tlr2和Tlr4均已被證實(shí)在大鼠腦缺血再灌注損傷后表達(dá)水平顯著升高［8］，同時(shí)甘草酸處理MCAO 大鼠可下調(diào)缺血腦內(nèi)TLR2、高遷移率族蛋白B1（high mobility group box 1，HMGB1）及基質(zhì)金屬蛋白酶9（matrix metallopro?teinase-9，MMP-9）從而發(fā)揮神經(jīng)保護(hù)作用［9］。腦缺血后TLR4/NF-κB 信號(hào)通路激活可能引起血小板聚集和星形膠質(zhì)細(xì)胞激活從而加重神經(jīng)損傷［10］。小鼠腦缺血72小時(shí)后調(diào)節(jié)T細(xì)胞表面CCR5表達(dá)增多，有利于減輕血腦屏障的損傷［11］。Georgakis等［12］基于全基因組關(guān)聯(lián)分析發(fā)現(xiàn)較高水平CCL2有更高的缺血性卒中遺傳傾向。這些結(jié)果均提示：炎癥反應(yīng)的調(diào)控對(duì)于腦缺血再灌注損傷的治療有著巨大的意義。CD68、CD86、CD14、TLR2和TLR4都與巨噬細(xì)胞的功能密切相關(guān)。腦卒中患者血液中CD14highCD16+單核細(xì)胞的數(shù)量上升，這類(lèi)單核細(xì)胞顯著高表達(dá)TLR2 蛋白，同時(shí)CD14highCD16+單核細(xì)胞的高水平與高死亡率和腦卒中早期臨床惡化相關(guān)［13］。結(jié)合本文分析結(jié)果，我們推測(cè)單核-巨噬細(xì)胞以及炎癥反應(yīng)的調(diào)控可能在腦缺血損傷中扮演了重要角色，可能是一個(gè)強(qiáng)有力的潛在治療靶點(diǎn)。此外，ISG15-/-的MCAO 小鼠死亡率高，梗死加劇，神經(jīng)系統(tǒng)恢復(fù)惡化，提示ISG15 可能作為內(nèi)源性神經(jīng)保護(hù)適應(yīng)因子能減輕腦卒中損傷［14］。

有文獻(xiàn)報(bào)道本文篩選出的核心基因Hcar2在小鼠腦缺血24 h后顯著升高，48 h開(kāi)始下降［15］。但其在腦缺血損傷中的作用以及機(jī)制仍不明確，另一核心基因Drd4編碼多巴胺受體D4（dopamine receptor 4，D4R），能被多巴胺、腎上腺素和去甲腎上腺素激活參與情緒控制、視網(wǎng)膜節(jié)細(xì)胞的晝夜節(jié)律的調(diào)節(jié)［16］。其在腦缺血的病理生理過(guò)程中的作用亦有待進(jìn)一步探索。艾司洛爾是一種選擇性β1腎上腺素能受體拮抗劑，對(duì)腦缺血有一定保護(hù)作用：缺血前給予艾司洛爾能有效減輕雙側(cè)頸動(dòng)脈閉塞引起的神經(jīng)功能障礙［17］。在沙鼠腦缺血模型中，甲巰咪唑可誘導(dǎo)甲狀腺功能減退以降低脂質(zhì)過(guò)氧化、增加超氧化物歧化酶（superoxide dismutase，SOD1），使海馬CA1 區(qū)神經(jīng)元的死亡得到緩解［18］。艾司洛爾和甲巰咪唑的腦保護(hù)作用目前已相對(duì)明確，但具體作用機(jī)制尚未闡明。吐根酚堿和水仙環(huán)素是否對(duì)腦缺血損傷有保護(hù)作用則未見(jiàn)報(bào)道。

本文進(jìn)行通路富集選用的GSEA 算法不需要事先篩選出差異基因，直接將全部基因與基因集進(jìn)行比較，能保留表達(dá)變化不大但功能重要的基因，分析結(jié)果更可靠。本研究雖然通過(guò)實(shí)時(shí)熒光定量pcr 對(duì)Hcar2及Drd4的表達(dá)水平進(jìn)行了驗(yàn)證，但未檢測(cè)篩選出的藥物的作用，仍須后續(xù)實(shí)驗(yàn)進(jìn)行驗(yàn)證。Hcar2及Drd4以及4種潛在藥物目前均未在腦缺血中得到充分的研究，可能是尚未發(fā)掘的新靶點(diǎn)。