周嬙 柏娜 劉生剛 劉偉 張宏偉 柳華

越來越多的研究證實(shí)遺傳因素在缺血性腦卒中（ischemic stroke,IS）發(fā)病中發(fā)揮著重要作用[1]，散發(fā)IS有50%的遺傳風(fēng)險(xiǎn)[2]。近期一個(gè)大型研究對不同種族GWAS數(shù)據(jù)進(jìn)行meta分析，確認(rèn)了32個(gè)卒中風(fēng)險(xiǎn)位點(diǎn)，其中20個(gè)為IS風(fēng)險(xiǎn)位點(diǎn)（12個(gè)為新發(fā)現(xiàn)位點(diǎn)），多數(shù)位點(diǎn)涉及血管病變，如血壓、心臟病、靜脈血栓形成等[3]。然而，IS潛在的致病基因和分子途徑尚未完全了解。因此，探索新的可能的關(guān)鍵風(fēng)險(xiǎn)基因?qū)⒂兄诎l(fā)現(xiàn)新的治療靶點(diǎn)，更詳細(xì)闡明IS的病理生理機(jī)制，提高IS早期的診斷和治療水平。近年來，基于高通量測序技術(shù)的生物信息學(xué)分析技術(shù)被廣泛應(yīng)用于探索疾病相關(guān)基因，進(jìn)而從分子層面揭示疾病的發(fā)生和發(fā)展機(jī)制。機(jī)器學(xué)習(xí)作為人工智能的核心技術(shù)，在生物醫(yī)學(xué)研究、個(gè)性化醫(yī)療、計(jì)算機(jī)輔助診斷等醫(yī)學(xué)領(lǐng)域有廣闊的前景。機(jī)器學(xué)習(xí)已經(jīng)應(yīng)用在IS診斷、預(yù)后等多個(gè)方面[4]。但把機(jī)器學(xué)習(xí)與生物信息學(xué)相結(jié)合以挖掘IS潛在靶基因的研究較少[5-6]。因此，本研究采用上述方法，從基因表達(dá)綜合數(shù)據(jù)庫（Gene Expression Omnibus，GEO）獲取人類IS的轉(zhuǎn)錄組數(shù)據(jù)集，進(jìn)行差異表達(dá)基因（differentially expressed gene,DEG）和功能富集分析，通過最小絕對值收斂和選擇算子（least absolute shrinkage and selection operator,LASSO）和支持向量機(jī)-遞歸特征消除（support vector machines-recursive feature elimination,SVMRFE）2種機(jī)器學(xué)習(xí)算法從DEGs中篩選潛在的IS關(guān)鍵基因，以期為IS有效診治提供依據(jù)。

1 資料與方法

1.1 數(shù)據(jù)下載和整理 從美國國家生物技術(shù)信息中心（National Center for Biotechnology Information，NCBI）GEO數(shù)據(jù)庫（https://www.ncbi.nlm.nih.gov/geo/）下載兩個(gè)IS轉(zhuǎn)錄組數(shù)據(jù)集GSE122709、GSE140275。其中，GSE122709作為訓(xùn)練數(shù)據(jù)集，包括10個(gè)IS患者樣本和5個(gè)健康對照樣本。患者均在起病4.5～24 h內(nèi)到達(dá)醫(yī)院，年齡50～75歲。IS定義為頭顱MRI彌散加權(quán)成像（diffusion-weighted imaging,DWI）見新發(fā)病灶，伴有急性神經(jīng)功能缺損，或者磁共振血管成像（magnetic resonance angiography，MRA）證實(shí)大腦前或中動(dòng)脈閉塞，美國國立衛(wèi)生研究院卒中量表（National Institute of Health Stroke Scale，NIHSS）評分中位數(shù)為12（6～21）。對照組由5名年齡、性別和血管危險(xiǎn)因素（包括體重指數(shù)、高血壓和高脂血癥）匹配的健康成人組成。GSE140275作為驗(yàn)證數(shù)據(jù)集，包含3個(gè)前循環(huán)大血管閉塞型IS患者樣本和3個(gè)健康對照樣本。病例組樣本隨機(jī)選取自2018年9月至2019年6月在廣西醫(yī)科大學(xué)附屬第一醫(yī)院就診的IS患者，對照樣本隨機(jī)選取自同期在同一醫(yī)院進(jìn)行體檢的健康志愿者。

由于GSE122709數(shù)據(jù)集同時(shí)包含mRNA和lncRNA轉(zhuǎn)錄本數(shù)據(jù)，本研究使用從GENCODE數(shù)據(jù)庫下載的人類參考基因組[Release 32（GRCh38.p13）]注釋文件（gencode.v32.annotation.gtf）[7]對mRNA進(jìn)行注釋和提取。

1.2 差異表達(dá)基因（DEG）篩選 從數(shù)據(jù)集GSE122709中提取并整理得到mRNA表達(dá)矩陣后，使用“DESeq2”軟件包以FDR＜0.05 和|log2FC|＞2為閾值進(jìn)行差異分析，篩選出DEGs。結(jié)果用“pheatmap”程序包繪制火山圖和熱圖可視化呈現(xiàn)。

1.3 DEGs富集分析 使用“ClusterProfiler”程序包對DEGs進(jìn)行基因本體論（gene ontplogy，GO）、京都基因與基因組百科全書（Kyoto Encyclopedia of Genes and Genomes，KEGG）通路富集分析，對基因及其產(chǎn)物進(jìn)行注釋，并明確IS相關(guān)DEG的具體生物學(xué)途徑。其中GO富集分析包含細(xì)胞組分（cell components，CC），分 子功 能（molecular function,MF），生物學(xué)過程（biological process，BP）3個(gè)部分[8]。使用“DOSE”程序包進(jìn)行疾病本體（disease ontology，DO）富集分析。P＜0.05作為基因顯著性富集的標(biāo)準(zhǔn)。

1.4 機(jī)器學(xué)習(xí)篩選關(guān)鍵基因 采用“glmnet”程序包用于LASSO回歸分析，“e1071”和“caret”程序包用于SVM-RFE。最終確定兩種算法之間共同識別的基因?yàn)殛P(guān)鍵基因。

1.5 驗(yàn)證關(guān)鍵基因表達(dá)及診斷效能 使用獨(dú)立的驗(yàn)證數(shù)據(jù)集GSE140275來判斷篩選出的標(biāo)志基因?qū)S和正常對照的潛在診斷價(jià)值。箱線圖分析標(biāo)志基因在IS和正常對照中表達(dá)差異。建立受試者工作特征（receiver operating characteristic curve，ROC）曲線，計(jì)算ROC曲線下面積（area under curve，AUC），評估關(guān)鍵基因的診斷效能[9]。

2 結(jié)果

2.1 IS的DEGs鑒定 GSE122709數(shù)據(jù)集含有15個(gè)樣本（包括10個(gè)IS患者和5個(gè)健康對照），對其進(jìn)行差異分析，共篩選出378個(gè)DEGs，包括176個(gè)上調(diào)基因和202個(gè)下調(diào)基因（圖1）。

圖1 GSE122709數(shù)據(jù)集中DEGs的熱圖（A）和火山圖（B）

2.2 DEGs富集分析 對篩選出的378個(gè)DEGs進(jìn)行GO、KEGG、DO富集分析。

GO功能富集分析發(fā)現(xiàn)，DEGs主要在血小板α顆粒、MHC II類蛋白復(fù)合物中顯著富集，具有多糖結(jié)合、CXCR趨化因子受體結(jié)合、G蛋白偶聯(lián)嘌呤核苷酸受體活性、趨化因子活性、受體配體活性、信號受體激活因子活性、細(xì)胞因子活性、碳水化合物結(jié)合、嘌呤核苷酸受體活性及核苷酸受體活性等分子功能，主要參與中心粒細(xì)胞遷移、中心粒細(xì)胞趨化性、骨髓白細(xì)胞游走、粒細(xì)胞遷移、粒細(xì)胞趨化性及G蛋白偶聯(lián)嘌呤能核苷酸受體等生物學(xué)過程（圖2A）。

KEGG通路富集分析發(fā)現(xiàn)DEGs主要富集在造血細(xì)胞譜系、移植物抗宿主病、病毒蛋白與細(xì)胞因子和細(xì)胞因子受體相互作用、I型糖尿病、COVID-19、腸道免疫網(wǎng)絡(luò)用于IgA生產(chǎn)、同種異體移植物排異、甘氨酸、絲氨酸和蘇氨酸代謝、趨化因子信號通路、流體剪切應(yīng)力與動(dòng)脈粥樣硬化等通路（圖2B）。

DO富集分析顯示DEGs的基因功能與卵巢癌、女性生殖器官癌、惡性卵巢表面上皮-間質(zhì)瘤、卵巢上皮癌、生殖系統(tǒng)疾病等疾病相關(guān)（圖2C）。

2.3 機(jī)器學(xué)習(xí)篩選關(guān)鍵基因 用LASSO和SVMRFE 2種算法對378個(gè)IS的DEGs進(jìn)行進(jìn)一步篩選。LASSO回歸分析篩選出7種特征基因（TBC1D3L、COLGALT2、TVP23C、MUC20、TAS2R4、CTRC、B3GAT1），SVM-RFE篩選得到14種特征基因（CDHR5、ZNF841、CCDC126、SPX、TAS2R3、KCNT2、DNAH14、B3GAT1、GMFB、TVP23C、SLC27A4、COPS2、GPR34、RAPGEFL1）。兩種算法共同識別到2個(gè)關(guān)鍵基因：B3GAT1、TVP23C（圖3）。

圖3 機(jī)器學(xué)習(xí)篩選關(guān)鍵基因

2.4 驗(yàn)證關(guān)鍵基因表達(dá)及診斷效能 在獨(dú)立的驗(yàn)證數(shù)據(jù)集GSE140275對B3GAT1、TVP23C進(jìn)行外部驗(yàn)證，發(fā)現(xiàn)B3GAT1、TVP23C的表達(dá)量在IS組和對照組具有統(tǒng)計(jì)學(xué)差異（P＜0.05）。同時(shí)，ROC曲線顯示B3GAT1、TVP23C在驗(yàn)證數(shù)據(jù)集的ROC曲線下面積（AUC）均接近1，即篩選出的2個(gè)關(guān)鍵基因在驗(yàn)證數(shù)據(jù)集內(nèi)對區(qū)分IS患者與健康對照者具有較高的診斷效能（圖4）。

圖4 關(guān)鍵基因驗(yàn)證

3 討論

本研究通過分析從GEO數(shù)據(jù)庫獲取的IS患者與健康對照組的mRNA轉(zhuǎn)錄本，鑒定出378個(gè)DEGs。GO富集分析發(fā)現(xiàn)DEGs參與的生物學(xué)過程主要集中在炎癥反應(yīng)方面，如中心粒細(xì)胞遷移、中心粒細(xì)胞趨化性、骨髓白細(xì)胞游走、粒細(xì)胞遷移和粒細(xì)胞趨化性等；KEGG通路富集分析發(fā)現(xiàn)DEGs除了與傳統(tǒng)IS危險(xiǎn)因素如糖尿病、動(dòng)脈粥樣硬化相關(guān)外，主要參與炎癥和免疫反應(yīng)通路，如造血細(xì)胞譜系、移植物抗宿主病、病毒蛋白與細(xì)胞因子和細(xì)胞因子受體相互作用、腸道免疫網(wǎng)絡(luò)用于IgA生產(chǎn)、同種異體移植物排異、甘氨酸、絲氨酸和蘇氨酸代謝、趨化因子信號通路等通路。應(yīng)用機(jī)器學(xué)習(xí)的LASSO和SVM-RFE算法對378個(gè)DEGs進(jìn)行進(jìn)一步篩選，并對篩選的特征基因進(jìn)行重疊，最終得到2個(gè)關(guān)鍵基因：B3GAT1、TVP23C。

B3GAT1（GlcAT-P）基因是葡萄糖醛酸轉(zhuǎn)移酶基因家族的成員，該基因產(chǎn)物GlcAT-P在碳水化合物表位HNK-1的生物合成過程中作為葡萄糖醛酸轉(zhuǎn)移反應(yīng)中的關(guān)鍵酶發(fā)揮作用[10]。HNK-1表位在神經(jīng)系統(tǒng)中高度表達(dá)[11]，GlcAT-P是大腦中主要的HNK-1合成酶，GlcAT-P基因敲除小鼠大腦中HNK-1抗原表位幾乎完全消失，從而導(dǎo)致突觸可塑性、記憶和學(xué)習(xí)障礙[12-13]。

有研究發(fā)現(xiàn)，AMPA受體亞基GluA2含有HNK-1表位[14]。AMPA型谷氨酸受體是由GluA 1-4亞基組成的離子型谷氨酸受體，它與NMDA谷氨酸受體以及紅藻氨酸受體的激活介導(dǎo)了中樞神經(jīng)系統(tǒng)中大部分的興奮性突觸傳遞[15]。大腦中動(dòng)脈閉塞模型中，抑制AMPA受體激活能夠抑制小膠質(zhì)細(xì)胞激活、促炎細(xì)胞因子表達(dá)和氧化應(yīng)激，顯著減少梗死體積[16]。AMPA受體存在2種形式：含GluA2、Ca2+不可滲透的AMPAR，缺乏GluA2、Ca2+可滲透的AMPAR，大多是AMPA受體是含有GluA2亞基，對Ca2+不可滲透的[17]。在缺血/再灌注模型中，AMRA受體被堆積的谷氨酸激活，GluA2被顯著內(nèi)吞，突觸后膜上的AMPA受體經(jīng)歷了從Ca2+不可滲透、含GluA2的AMPA受體到Ca2+可滲透、缺乏GluA2的AMPA受體的亞基組成轉(zhuǎn)換，使Ca2+持續(xù)進(jìn)入細(xì)胞內(nèi)，導(dǎo)致細(xì)胞內(nèi)鈣超載，觸發(fā)神經(jīng)元興奮性毒性誘導(dǎo)細(xì)胞凋亡[18]。目前，缺血后GluR2表達(dá)調(diào)控的分子機(jī)制尚不清楚，GluR2上的HNK-1由GlcAT-P和HNK-1磺基轉(zhuǎn)移酶在高爾基體（golgi apparatus,GA）合成，與攜帶HNK-1的GluA2相比，不攜帶HNK-1的GluA2被顯著內(nèi)吞并且在細(xì)胞表面表達(dá)較少[14，19]。

基于以上研究，我們推測B3GAT1可能通過調(diào)控AMRA谷氨酸受體GluR2亞基上的HNK-1表位影響Ca2+內(nèi)流，導(dǎo)致細(xì)胞內(nèi)Ca2+超載，從而在IS病理生理過程中發(fā)揮重要作用。但B3GAT1在IS中的特異性表達(dá)和作用機(jī)制仍有待于進(jìn)一步深入挖掘。

TVP23C基因產(chǎn)物TVP23C是GA膜蛋白TVP23同源物。在腦缺血-再灌注模型中可以觀察到損傷細(xì)胞中GA碎裂和細(xì)胞凋亡[20]。眾多研究表明GA參與離子穩(wěn)態(tài)、細(xì)胞凋亡等氧化應(yīng)激過程，被稱為“GA應(yīng)激”[21]，在IS中同樣存在GA應(yīng)激[22]。

GA通過維持細(xì)胞Ca2+穩(wěn)態(tài)在氧化應(yīng)激過程中起保護(hù)作用。此類保護(hù)作用的關(guān)鍵在于分泌途徑中質(zhì)膜相關(guān)Ca2+ATP酶-SPCA1。細(xì)胞受到刺激后SPCA1可以轉(zhuǎn)運(yùn)細(xì)胞質(zhì)內(nèi)Ca2+到GA腔內(nèi)，從而緩解細(xì)胞質(zhì)鈣超載[21,23]。增加SPCA1的表達(dá)可以抑制GA應(yīng)激減少腦缺血損傷[24-25]。SPCA1敲除小鼠中則表現(xiàn)出神經(jīng)細(xì)胞凋亡增加[26-27]。同時(shí)，GA本身在細(xì)胞凋亡中也起著重要作用，許多凋亡調(diào)控成分已被鑒定并定位于GA，如Fas、Hippi蛋白、腫瘤壞死因子受體-1、Bcl-2家族成員，以及含泛素連接酶的桿狀病毒IAP重復(fù)序列、GA抗凋亡蛋白等[21]。因此，本研究推測TVP23C是通過GA活動(dòng)進(jìn)而參與IS的病理生理機(jī)制。

本研究存在著一定的局限性。首先本研究數(shù)據(jù)來源于GEO數(shù)據(jù)庫，樣本方面有一定的局限性。其次，本研究通過繪制數(shù)據(jù)集的B3GAT1、TVP23C的ROC曲線，發(fā)現(xiàn)AUC值均接近于1，說明這2個(gè)基因在診斷IS方面具有良好的效能，但這個(gè)結(jié)果可能與驗(yàn)證集樣本量太少、沒有細(xì)致進(jìn)行TOAST分型等因素有關(guān)，所以未來需要更大的樣本進(jìn)一步驗(yàn)證關(guān)鍵基因的診斷價(jià)值。

綜上所述，本研究通過生物信息學(xué)和機(jī)器學(xué)習(xí)結(jié)合的方法篩選出TVP23C、B3GAT1可能為IS的關(guān)鍵風(fēng)險(xiǎn)基因，對于IS診斷方面有較好的效能。結(jié)合B3GAT1的表達(dá)分析，推測B3GAT1可能通過調(diào)控AMRA谷氨酸受體參與IS缺血性損傷的病理生理過程，這為診斷IS的潛在生物標(biāo)志物和治療靶點(diǎn)提供新的思路。