基于生物信息学和机器学习方法探索缺血性脑卒中关键风险基因☆

周嫱 柏娜 刘生刚 刘伟 张宏伟 柳华

越来越多的研究证实遗传因素在缺血性脑卒中（ischemic stroke,IS）发病中发挥着重要作用[1]，散发IS有50%的遗传风险[2]。近期一个大型研究对不同种族GWAS数据进行meta分析，确认了32个卒中风险位点，其中20个为IS风险位点（12个为新发现位点），多数位点涉及血管病变，如血压、心脏病、静脉血栓形成等[3]。然而，IS潜在的致病基因和分子途径尚未完全了解。因此，探索新的可能的关键风险基因将有助于发现新的治疗靶点，更详细阐明IS的病理生理机制，提高IS早期的诊断和治疗水平。近年来，基于高通量测序技术的生物信息学分析技术被广泛应用于探索疾病相关基因，进而从分子层面揭示疾病的发生和发展机制。机器学习作为人工智能的核心技术，在生物医学研究、个性化医疗、计算机辅助诊断等医学领域有广阔的前景。机器学习已经应用在IS诊断、预后等多个方面[4]。但把机器学习与生物信息学相结合以挖掘IS潜在靶基因的研究较少[5-6]。因此，本研究采用上述方法，从基因表达综合数据库（Gene Expression Omnibus，GEO）获取人类IS的转录组数据集，进行差异表达基因（differentially expressed gene,DEG）和功能富集分析，通过最小绝对值收敛和选择算子（least absolute shrinkage and selection operator,LASSO）和支持向量机-递归特征消除（support vector machines-recursive feature elimination,SVMRFE）2种机器学习算法从DEGs中筛选潜在的IS关键基因，以期为IS有效诊治提供依据。

1 资料与方法

1.1 数据下载和整理 从美国国家生物技术信息中心（National Center for Biotechnology Information，NCBI）GEO数据库（https://www.ncbi.nlm.nih.gov/geo/）下载两个IS转录组数据集GSE122709、GSE140275。其中，GSE122709作为训练数据集，包括10个IS患者样本和5个健康对照样本。患者均在起病4.5～24 h内到达医院，年龄50～75岁。IS定义为头颅MRI弥散加权成像（diffusion-weighted imaging,DWI）见新发病灶，伴有急性神经功能缺损，或者磁共振血管成像（magnetic resonance angiography，MRA）证实大脑前或中动脉闭塞，美国国立卫生研究院卒中量表（National Institute of Health Stroke Scale，NIHSS）评分中位数为12（6～21）。对照组由5名年龄、性别和血管危险因素（包括体重指数、高血压和高脂血症）匹配的健康成人组成。GSE140275作为验证数据集，包含3个前循环大血管闭塞型IS患者样本和3个健康对照样本。病例组样本随机选取自2018年9月至2019年6月在广西医科大学附属第一医院就诊的IS患者，对照样本随机选取自同期在同一医院进行体检的健康志愿者。

由于GSE122709数据集同时包含mRNA和lncRNA转录本数据，本研究使用从GENCODE数据库下载的人类参考基因组[Release 32（GRCh38.p13）]注释文件（gencode.v32.annotation.gtf）[7]对mRNA进行注释和提取。

1.2 差异表达基因（DEG）筛选 从数据集GSE122709中提取并整理得到mRNA表达矩阵后，使用“DESeq2”软件包以FDR＜0.05 和|log2FC|＞2为阈值进行差异分析，筛选出DEGs。结果用“pheatmap”程序包绘制火山图和热图可视化呈现。

1.3 DEGs富集分析 使用“ClusterProfiler”程序包对DEGs进行基因本体论（gene ontplogy，GO）、京都基因与基因组百科全书（Kyoto Encyclopedia of Genes and Genomes，KEGG）通路富集分析，对基因及其产物进行注释，并明确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 验证关键基因表达及诊断效能 使用独立的验证数据集GSE140275来判断筛选出的标志基因对IS和正常对照的潜在诊断价值。箱线图分析标志基因在IS和正常对照中表达差异。建立受试者工作特征（receiver operating characteristic curve，ROC）曲线，计算ROC曲线下面积（area under curve，AUC），评估关键基因的诊断效能[9]。

2 结果

2.1 IS的DEGs鉴定 GSE122709数据集含有15个样本（包括10个IS患者和5个健康对照），对其进行差异分析，共筛选出378个DEGs，包括176个上调基因和202个下调基因（图1）。

图1 GSE122709数据集中DEGs的热图（A）和火山图（B）

2.2 DEGs富集分析 对筛选出的378个DEGs进行GO、KEGG、DO富集分析。

GO功能富集分析发现，DEGs主要在血小板α颗粒、MHC II类蛋白复合物中显著富集，具有多糖结合、CXCR趋化因子受体结合、G蛋白偶联嘌呤核苷酸受体活性、趋化因子活性、受体配体活性、信号受体激活因子活性、细胞因子活性、碳水化合物结合、嘌呤核苷酸受体活性及核苷酸受体活性等分子功能，主要参与中心粒细胞迁移、中心粒细胞趋化性、骨髓白细胞游走、粒细胞迁移、粒细胞趋化性及G蛋白偶联嘌呤能核苷酸受体等生物学过程（图2A）。

KEGG通路富集分析发现DEGs主要富集在造血细胞谱系、移植物抗宿主病、病毒蛋白与细胞因子和细胞因子受体相互作用、I型糖尿病、COVID-19、肠道免疫网络用于IgA生产、同种异体移植物排异、甘氨酸、丝氨酸和苏氨酸代谢、趋化因子信号通路、流体剪切应力与动脉粥样硬化等通路（图2B）。

DO富集分析显示DEGs的基因功能与卵巢癌、女性生殖器官癌、恶性卵巢表面上皮-间质瘤、卵巢上皮癌、生殖系统疾病等疾病相关（图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 验证关键基因表达及诊断效能 在独立的验证数据集GSE140275对B3GAT1、TVP23C进行外部验证，发现B3GAT1、TVP23C的表达量在IS组和对照组具有统计学差异（P＜0.05）。同时，ROC曲线显示B3GAT1、TVP23C在验证数据集的ROC曲线下面积（AUC）均接近1，即筛选出的2个关键基因在验证数据集内对区分IS患者与健康对照者具有较高的诊断效能（图4）。

图4 关键基因验证

3 讨论

本研究通过分析从GEO数据库获取的IS患者与健康对照组的mRNA转录本，鉴定出378个DEGs。GO富集分析发现DEGs参与的生物学过程主要集中在炎症反应方面，如中心粒细胞迁移、中心粒细胞趋化性、骨髓白细胞游走、粒细胞迁移和粒细胞趋化性等；KEGG通路富集分析发现DEGs除了与传统IS危险因素如糖尿病、动脉粥样硬化相关外，主要参与炎症和免疫反应通路，如造血细胞谱系、移植物抗宿主病、病毒蛋白与细胞因子和细胞因子受体相互作用、肠道免疫网络用于IgA生产、同种异体移植物排异、甘氨酸、丝氨酸和苏氨酸代谢、趋化因子信号通路等通路。应用机器学习的LASSO和SVM-RFE算法对378个DEGs进行进一步筛选，并对筛选的特征基因进行重叠，最终得到2个关键基因：B3GAT1、TVP23C。

B3GAT1（GlcAT-P）基因是葡萄糖醛酸转移酶基因家族的成员，该基因产物GlcAT-P在碳水化合物表位HNK-1的生物合成过程中作为葡萄糖醛酸转移反应中的关键酶发挥作用[10]。HNK-1表位在神经系统中高度表达[11]，GlcAT-P是大脑中主要的HNK-1合成酶，GlcAT-P基因敲除小鼠大脑中HNK-1抗原表位几乎完全消失，从而导致突触可塑性、记忆和学习障碍[12-13]。

有研究发现，AMPA受体亚基GluA2含有HNK-1表位[14]。AMPA型谷氨酸受体是由GluA 1-4亚基组成的离子型谷氨酸受体，它与NMDA谷氨酸受体以及红藻氨酸受体的激活介导了中枢神经系统中大部分的兴奋性突触传递[15]。大脑中动脉闭塞模型中，抑制AMPA受体激活能够抑制小胶质细胞激活、促炎细胞因子表达和氧化应激，显著减少梗死体积[16]。AMPA受体存在2种形式：含GluA2、Ca2+不可渗透的AMPAR，缺乏GluA2、Ca2+可渗透的AMPAR，大多是AMPA受体是含有GluA2亚基，对Ca2+不可渗透的[17]。在缺血/再灌注模型中，AMRA受体被堆积的谷氨酸激活，GluA2被显著内吞，突触后膜上的AMPA受体经历了从Ca2+不可渗透、含GluA2的AMPA受体到Ca2+可渗透、缺乏GluA2的AMPA受体的亚基组成转换，使Ca2+持续进入细胞内，导致细胞内钙超载，触发神经元兴奋性毒性诱导细胞凋亡[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病理生理过程中发挥重要作用。但B3GAT1在IS中的特异性表达和作用机制仍有待于进一步深入挖掘。

TVP23C基因产物TVP23C是GA膜蛋白TVP23同源物。在脑缺血-再灌注模型中可以观察到损伤细胞中GA碎裂和细胞凋亡[20]。众多研究表明GA参与离子稳态、细胞凋亡等氧化应激过程，被称为“GA应激”[21]，在IS中同样存在GA应激[22]。

GA通过维持细胞Ca2+稳态在氧化应激过程中起保护作用。此类保护作用的关键在于分泌途径中质膜相关Ca2+ATP酶-SPCA1。细胞受到刺激后SPCA1可以转运细胞质内Ca2+到GA腔内，从而缓解细胞质钙超载[21,23]。增加SPCA1的表达可以抑制GA应激减少脑缺血损伤[24-25]。SPCA1敲除小鼠中则表现出神经细胞凋亡增加[26-27]。同时，GA本身在细胞凋亡中也起着重要作用，许多凋亡调控成分已被鉴定并定位于GA，如Fas、Hippi蛋白、肿瘤坏死因子受体-1、Bcl-2家族成员，以及含泛素连接酶的杆状病毒IAP重复序列、GA抗凋亡蛋白等[21]。因此，本研究推测TVP23C是通过GA活动进而参与IS的病理生理机制。

本研究存在着一定的局限性。首先本研究数据来源于GEO数据库，样本方面有一定的局限性。其次，本研究通过绘制数据集的B3GAT1、TVP23C的ROC曲线，发现AUC值均接近于1，说明这2个基因在诊断IS方面具有良好的效能，但这个结果可能与验证集样本量太少、没有细致进行TOAST分型等因素有关，所以未来需要更大的样本进一步验证关键基因的诊断价值。

综上所述，本研究通过生物信息学和机器学习结合的方法筛选出TVP23C、B3GAT1可能为IS的关键风险基因，对于IS诊断方面有较好的效能。结合B3GAT1的表达分析，推测B3GAT1可能通过调控AMRA谷氨酸受体参与IS缺血性损伤的病理生理过程，这为诊断IS的潜在生物标志物和治疗靶点提供新的思路。