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

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

1 資料與方法

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

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

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

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

1.4 機器學習篩選關鍵基因 采用“glmnet”程序包用于LASSO回歸分析，“e1071”和“caret”程序包用于SVM-RFE。最終確定兩種算法之間共同識別的基因為關鍵基因。

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

2 結果

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

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

2.2 DEGs富集分析 對篩選出的378個DEGs進行GO、KEGG、DO富集分析。

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

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

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

2.3 機器學習篩選關鍵基因 用LASSO和SVMRFE 2種算法對378個IS的DEGs進行進一步篩選。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個關鍵基因：B3GAT1、TVP23C（圖3）。

圖3 機器學習篩選關鍵基因

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

圖4 關鍵基因驗證

3 討論

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

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

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

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

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

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

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

綜上所述，本研究通過生物信息學和機器學習結合的方法篩選出TVP23C、B3GAT1可能為IS的關鍵風險基因，對于IS診斷方面有較好的效能。結合B3GAT1的表達分析，推測B3GAT1可能通過調控AMRA谷氨酸受體參與IS缺血性損傷的病理生理過程，這為診斷IS的潛在生物標志物和治療靶點提供新的思路。