摘要
本研究基于生物信息学和机器学习方法探究多发性硬化症（Multiple sclerosis, MS）的关键基因，从GEO数据库获取MS基因表达谱GSE21942和GSE32988，GSE32988作为验证数据集，采用主成分分析（Principal component analysis, PCA）评估样本聚类情况，筛选差异表达基因（Differential expression analysis, DEGs）并进行基因本体论分析（Gene ontology, GO）和京都基因与基因组百科全书（Kyoto encyclopedia of genes and genomes, KEGG）通路富集分析，使用加权基因共表达网络分析（Weighted gene co-expression network analysis, WGCNA）筛选与MS密切相关的基因模块，将基因模块与DEGs取交集获得候选基因，采用最小绝对收缩和选择算子算法（Least absolute shrinkage and selection operator, LASSO）及随机森林（Random forest, RF）对候选基因筛选得到潜在关键基因，采用第三方数据集GSE32988验证潜在关键基因的差异表达。通过进行受试者工作特征曲线（Receiver operating characteristic curve, ROC）验证，得到关键基因。采用MS动物模型验证关键基因表达水平。结果显示GSE21942具有良好重复性和相关性，共获得506个DEGs，富集分析结果表明DEGs主要富集于B细胞活化、谷氨酸（Glutamic acid, GLU）代谢、氧化应激（Oxidative stress, OS）等生物功能、内质网（Endoplasmic reticulum, ER）等细胞组分以及EB病毒感染（Epstein-barr virus, EBV）等。29个候选基因经机器学习算法筛选得到5个潜在关键基因，用GSE32988验证后得到GLUD1、VDAC1、DDX3X、LAMP1共4个关键基因。RT-qPCR鉴定DDX3X、LAMP1、GLUD1和VDAC1的表达水平与mRNA芯片的生物信息学分析结果一致；因此，DDX3X、LAMP1、GLUD1和VDAC1有可能成为MS治疗的新靶标。
Abstract
In order to explore the key genes of multiple sclerosis (MS) based on bioinformatics and machine learning methods; MS gene expression profiles GSE21942 and GSE32988 were obtained from the GEO database. GSE32988 was used as a validation dataset, Sample clustering was assessed using PCA to screen for differentially expressed genes (DEGs) and analyzed for GO and KEGG enrichment, the gene modules closely related to MS were identified using weighted gene co-expression network analysis (WGCNA), the intersection of the gene modules and DEGs was analyzed to obtain candidate genes. Candidate genes were screened to obtain potential key genes using machine learning algorithms, which include the Least Absolute Shrinkage Operator Algorithm (LASSO) and the Random Forest Algorithm (RF). The third-party dataset GSE32988 was used to validate the differential expression of potential key genes. key genes were obtained by performing subject operating characteristic curve (ROC) validation. An MS animal model was used to verify the expression levels of key genes; The results show that GSE21942 showed good repeatability and correlation, and a total of 506 DEGs were obtained. Enrichment analysis showed that DEGs were mainly enriched in biological functions such as B cell activation, glutamic acid (GLU) metabolism, oxidative stress (OS), as well as in EBV infection and the B cell receptor signaling pathway, etc. The 29 candidate genes were screened by a machine learning algorithm to obtain five potential key genes, and a total of four key genes, GLUD1, VDAC1, DDX3X, and LAMP1, were obtained after validation with GSE32988. RT-qPCR identified the expression levels of DDX3X, LAMP1, GLUD1, and VDAC1 in accordance with the results of bioinformatics analysis of mRNA microarrays; Consequently, DDX3X, LAMP1, GLUD1, and VDAC1 may become new targets for MS therapy.
Keywords
多发性硬化症（Multiple sclerosis，MS）是一种中枢神经系统（Central nervous system，CNS）自身免疫性脱髓鞘疾病[1]，截至2020年，全球约280万人（每100 000人中35.9人）患有MS，且自2013年以来，世界MS发病率和患病率均有所增加。MS发病可能与多种因素相关[2]，如遗传易感性、生活方式、病毒感染及免疫机制等，但致病机理尚不明确[3]，其病因及致病机理仍是未来研究的重点。目前研究认为，小胶质细胞与MS的发病密切相关[4]。MS的特征是白质和灰质的炎症性脱髓鞘病变[5]，其发病主要累及大脑、脑干、脊髓和视神经，可引起肢体无力瘫痪、感觉麻木疼痛、二便障碍以及情绪障碍[6]，而且有着极高的复发率和致残率，给社会造成了巨大的经济及医疗压力，其现有的治疗靶点主要包括树突状细胞、小胶质细胞、B细胞表面CD20等[7]。
机器学习作为高效处理生物信息的方法，侧重于通过非线性模型推断变量中的重要风险因素，已成为数据挖掘的重要手段[8]。本研究将机器学习及生物信息学技术结合，并利用MS动物模型进行实验验证，旨在探索MS关键基因，这将有助于深入了解疾病的发生发展，为其干预治疗提供新思路。
1 材料与方法
1.1 数据集获取与PCA分析
从NCBI GEO数据库（http: //www.ncbi.nlm.nih.gov/geo）获取MS基因表达谱。筛选标准：关键词为“Multiple sclerosis”，筛选条件为“Homo sapiens”和“Expression profiling by array”，获得GPL570平台由Kemppinen A等提供的芯片数据集GSE21942[9]和GPL6480平台由Sousa AP等提供的芯片数据集GSE32988[10]。
GSE21942芯片数据集为主数据集，包含54 675个基因和29个来自外周血单核细胞的样本，分为对照组与MS组，对照组含有15个健康样本，MS组含有14个MS疾病样本。GSE32988为验证数据集，包含41 000个基因和48个来自T细胞的样本，分为对照组与MS组，对照组含有8个健康样本，MS组含有40个疾病样本。对数据集GSE21942进行PCA，以评估样本聚类情况。
1.2 鉴定DEGs
利用GEO2R将基因数据集GSE21942按照MS与正常人进行分组，采用筛选标准adj.P. Value＜0.05、|log2 FC|≥1获得DEGs。
1.3 加权基因共表达网络分析
在R软件（http: //www.r-project.org）中，利用“WGCNA”包对GSE21942进行加权基因共表达网络分析（Weighted gene co-expression network analysis，WGCNA），构建基因共表达网络，基于无标度网络R2=0.85和平均连通度（Mean connectivity）＜100，获得合适的软阈值（β）。设置最小模块30，模块合并阈值0.25鉴定共表达模块，得到特征模块基因。使用“VennDiagram”包将特征模块基因与DEGs取交集，作为候选基因进行后续分析。
1.4 差异表达基因功能富集分析
使用“clusterProfiler”包对差异表达基因（Differential expression analysis，DEGs）进行基因本体论（Gene ontology，GO）富集分析，在DAVID数据库（https: //david.ncifcrf.gov/home.jsp）中对DEGs进行京都基因与基因组百科全书（Kyoto encyclopedia of genes and genomes，KEGG）通路富集分析。P＜0.05视为显著富集。
1.5 机器学习算法筛选潜在关键基因
最小绝对收缩和选择算子算法（Least absolute shrinkage and selection operator，LASSO）是一种采用了L1正则化的线性回归方法，通过参数缩减拟合广义线性模型的同时进行变量筛选，从而达到降维和获得最优模型的目的[11]。利用R软件中的“glmnet”包，对候选基因进行LASSO分析，通过十倍交叉验证选择合适的惩罚系数进一步筛选特征基因。
随机森林（Random forest，RF）算法对于分析高维数据潜在关系具有优势。使用R包“randomforest”对候选基因进行RF算法分析筛选特征基因。随后将两种特征基因取交集作为潜在关键基因。
1.6 外部数据集验证和ROC诊断分析关键基因
首先利用外部基因集GSE32988对潜在关键基因进行验证，以P ＜ 0.05为筛选标准得到关键基因。随后通过R软件中的“pROC”包进行受试者工作特征曲线（Receiver operating characteristic curve，ROC）诊断分析，在AUC＞0.5的情况下，越接近于1，说明该变量的预测结局诊断效果越好。AUC＞0.9时视为有较高的准确性，0.7≤AUC≤0.9时视为具有一定的准确性。
1.7 关键基因的定量RT-qPCR验证
SPF级C57BL/6小鼠6只，雌性，6~8周龄，体重16~18 g，购买自浙江维通利华实验动物技术有限公司，随机分为对照组和MS组，每组3只。对照组采用含有0.2% Cuprizone（Sigma公司，美国）的饲料喂养5周以诱导MS模型，对照组小鼠采用正常饲料进行饲养，处死后取小鼠大脑胼胝体组织，液氮冷冻后储存在-80℃中保存备用。采用TRIzol法提取组织总RNA。使用Fastking RT Kit试剂盒（天根生化科技有限公司，北京）进行逆转录，使用KAPA SYBR FAST试剂盒（KAPA公司，美国）检测mRNA水平。实验条件：预变性95℃、20 s，然后95℃、3 s和60℃、20 s进行40个循环。将GAPDH作为内参进行数据处理，利用2-△△Ct法计算mRNA的表达量。所有实验均经河南中医药大学实验动物伦理委员会批准（DWLL202103124）。本实验未使用多余的实验动物。引物序列如表1所示。
2 结果与讨论
2.1 主成分分析结果
主成分分析（Principal component analysis，PCA）结果显示主数据集GSE21942的两组样本之间区分度高，组内样本呈现出良好的重复性和相关性，第一主成分（PCA1）占比16.58%，第二主成分（PCA2）占比11.05%，如图1所示。
表1实时荧光定量PCR引物序列
Table1Real-time fluorescence quantitative PCR primer sequences
图1样本分布PCA散点图
Fig.1Sample distribution PCA scatter plot
2.2 筛选MS的DEGs
差异分析结果显示，共有506个DEGs，包括298个上调基因，208个下调基因，如图2所示。
2.3 DEGs的所有功能富集分析
对所有MS相关DEGs采用进行GO富集分析，GO富集结果按照FDR＜0.01进行筛选，KEGG通路富集结果按照P＜0.05进行筛选，选取最显著富集的条目进行可视化（图3（a））。KEGG通路富集分析结果表明（图3（b））这些基因主要参与EB病毒（Epstein-Barr virus，EBV）感染等通路，并影响氧化应激（Oxidative stress，OS）反应、谷氨酸（Glutamic acid，GLU）代谢等生物功能。
2.4 通过加权基因共表达网络分析筛选与MS病变相关的模块
基因表达量矩阵来自于GSE21942，根据拟合指数（图4（a））和平均连接程度（图4（b）），采用无标度拓扑指数R2 ＞ 0.85，软阈值β=6确保构建无标度网络。通过基因聚类分析，获得基因树状图（图4（c）），具有高表达相关性的基因往往被归入同一模块，共获得20个基因模块。
分析基因表达模块与临床特征的关联性，绘制基因网络中拓扑重叠的热图。结果显示蓝色（BLUE）模块（|R|=0.89，p=1.7×10-10）与MS相关性最高（图5），此模块基因被认为是进一步分析的关键模块，因此作为与MS相关的特征模块进行后续分析。共获得385个模块特征基因。分析模块特征基因与临床特征的相关性（图6，）相关性系数为0.9，表明BLUE模块特征基因与MS相关性密切。
图2GES21942差异分析
Fig.2Difference analysis in GSE21942
注：（a）与MS相关的DEGs可视化火山图;（b）与MS相关的DEGs可视化热图.
图3DEGs的GO功能富集分析和KEGG通路富集分析
Fig.3GO and KEGG analysis of DEGs
注：（a）GO功能富集分析；（b）KEGG通路富集析析.
图4基因共表达网络构建
Fig.4Gene co-expression network construction
注：（a）拟合指数；（b）平均连接程度；（c）基因树状图.
图5各模块与MS的关联性
Fig.5Heat map of Module-trait associations
图6BLUE模块基因与MS的关联性
Fig.6Association plot of BLUE module genes with MS
2.5 获取候选基因
将DEGs与BLUE模块特征基因取交集，获得29个基因（图7），将其作为候选基因进一步筛选关键基因。
2.6 通过机器学习算法筛选潜在关键基因
对29个候选基因进行LASSO分析，随着λ值增加（图8），一些基因的系数降为0。通过十倍交叉验证（图9），当λ=0.04时，建立了由8个基因组成的最佳模型：VCAD1、GLUD1、ST8SIA4、NONO、DDX3X、LAMP1、EIF5、ZBTB20，将此作为MS特征基因。
图7DEGs与BLUE模块基因的维恩图
Fig.7Venn diagram of DEGs-BLUE
图8交叉基因的LASSO系数曲线
Fig.8LASSO coefficient curves for crossover genes
图9LASSO回归系数分析
Fig.9LASSO screening for potential key candidate genes
对29个候选基因进行RF算法分析，选择Mean Decrease Accuracy ＞ 2且Mean Decrease Gini ＞ 0.05的基因作为区分对照组和MS组的重要基因（图10），一共获得9个特征基因（CYBB、HIPK3、TOMM2、ZFP36L1、EIF5、GLUD1、VDAC1、DDX3X、LAMP1）。
将RF结果与LASSO结果取交集（图11），得到5个潜在关键基因（DDX3X、EIF5、GLUD1、LAMP1、VDAC1）进行后续分析。
图10MS相关性基因的RF重要性排序
Fig.10Importance ranking of MS-associated genes
图11RF与LASSO的维恩图
Fig.11Venn diagram of LASSO-RF
2.7 与潜在关键基因表达一致的外部数据集验证和ROC诊断分析
采用外部数据集GSE32988验证5个潜在关键基因的表达（图12），其中4个基因（DDX3X、VDAC1、LAMP1、GLUD1）存在差异表达， 1个基因（EIF5）表达无显著差异。随后，将筛选出的与外部数据集表达一致的4个潜在关键基因进行ROC分析（图13），ROC曲线下方面积大小（Area Under Curve，AUC）分别为GLUD1（AUC=0.7995%， CI:1.00~0.58）、DDX3X（AUC=0.8195%，CI：0.93~0.69）、LAMP1（AUC=0.7695%，CI：0.95~0.97）、VDAC1（AUC=0.8195%，CI：0.99~0.64）。结果表明DDX3X、VDAC1、LAMP1对于诊断MS具有较高的准确性，GLUD1对于诊断MS具有一定的准确性。
2.8 在MS小鼠模型中验证关键基因
采用RT-qPCR检测CPZ小鼠胼胝体组织中关键基因DDX3X、VCAD1、LAMP1和GLUD1 mRNA表达水平。结果显示，与对照组比较，MS组小鼠VDAC1、LAMP1和GLUD1的mRNA表达显著增高，DDX3X的mRNA表达显著降低，差异均有统计学意义（均P＜0.05）（图14）。结果与外部数据集验证和ROC诊断分析（2.7小节）的结果一致。
图12潜在关键基因的外部基因集验证
Fig.12Key genes set validation of key genes
图13评估GSE32988中四个关键基因预测准确性的ROC曲线图
Fig.13ROC curve to evaluate prediction accuracy of the four hub genes in GSE32988
图14关键基因mRNA表达验证
Fig.14Validation of mRNA expression of key genes
3 讨论
本研究将机器学习算法与生物信息学相结合筛选并鉴定MS的关键基因，对于寻找MS的潜在治疗靶点具有重要意义。
GO富集结果表明，DEGs主要与OS反应、GLU代谢等生物过程和内质网（Endoplasmic reticulum，ER）等细胞组分有关。OS反应是机体内氧化产物和抗氧化防御机制之间平衡紊乱的应激状态[12]，可以导致少突胶质细胞的损伤[13]，甚至诱导细胞死亡导致髓鞘脱失，从而导致脱髓鞘症状[14]。ER是存在于所有真核细胞中的细胞器[15]，已有研究证明，ER的功能失常与神经退行性疾病的发生和进展联系密切[16]，当ER应激时，可导致线粒体功能障碍，从而进一步影响MS等神经退行性病变[17]。ER的小管分布受到溶酶体调节，这种调节作用在神经元发育中起着关键作用[18]，溶酶体对ER的降解调控途径是ER自噬[19]，一些神经系统疾病是由参与ER吞噬的基因突变引起的[15]，目前ER自噬参与控制神经炎症疾病已经被越来越多的证据证明[16]。GLU是CNS中主要的兴奋性神经递质，在认知功能中起关键作用[20]，当过量GLU出现时，会出现兴奋性毒性[21]，从而导致MS认知障碍等常见症状发生，这可能与其导致的神经元功能障碍和变性有关[22]。同时，GLU代谢可影响髓鞘的形成与再生[23]，且抑制神经元释放GLU导致髓鞘形成显著下降约39%[24]。
KEGG通路富集分析表明，DEGs主要在EBV感染和B细胞受体等信号通路中被激活。EBV感染是MS致病的重要危险因素[25]，机体在感染EBV后，患MS的风险会增加30倍以上[26]。B细胞作为EBV的主要宿主[27]，可通过多种途径影响MS的发病及转归，如进入CNS，影响异位淋巴组织生成等[28]。
电压依赖性阴离子通道（Voltage-dependent anion channel，VDAC），是线粒体外膜上存在的一种孔蛋白，是线粒体内代谢物交换的主要通道。哺乳动物线粒体表现出三种VDAC亚型（VDAC1、VDAC2和VDAC3），其中VDAC1表达最广泛[29]。VDAC1可以维持线粒体的正常功能并参与调控OS反应[30]。而线粒体已被证明与多种神经退行性疾病密切相关[31]。有研究表明，MS患者的免疫细胞中线粒体活动、代谢、线粒体DNA水平和线粒体形态存在异常[32]，VDAC1作为线粒体通透性转换孔的重要组成部分，它不仅调控着线粒体外膜的通透性，还紧密关联着细胞的凋亡调节过程[33]。OS反应是活性物质的产生和抗氧化防御机制之间的不平衡状态，被认为是多发性硬化症病理生理学的关键因素之一[34]，在MS的疾病进展过程中发挥着重要作用[35]，它与炎症反应关系密切[36]，其应激水平与MS病灶的敏感性之间存在正相关[37]。溶酶体相关膜蛋白（Lysosome-associated membrane protein，LAMP）家族是一类约由200个氨基酸组成的糖蛋白，其中LAMP1（也称为CD107a）是一种Ⅰ型跨膜糖蛋白，在维持溶酶体的完整性、pH和分解代谢上发挥重要作用[38]。作为溶酶体的标记物[39]，LAMP1缺失时会强烈影响溶酶体的功能，如胆固醇的运输[40]。胆固醇作为髓鞘的重要组成成分，对于大脑髓鞘膜生长有着重要作用[41]，它的代谢异常会损害髓鞘形成，从而导致多种神经退行性疾病的发生[42]，虽然目前MS活动与脂质代谢异常之间的因果关系仍未具体阐明，但该领域的进展有可能提高对MS的理解[43]。DDX3X是一种ATP依赖性的RNA解旋酶，DDX3X的突变或失调不仅与多种神经退行性和神经发育障碍有关[44]，而且还可诱发EBV感染[45]，EBV感染已被证实成为MS的重要诱导因素[46]。GLUD1是一种线粒体谷氨酸脱氢酶，在GLU的分解与合成中扮演着重要角色[47]。GLU作为神经递质，其传递异常与神经退行性疾病关系密切[48]。有实验表明，在转基因小鼠中，当GLUD1功能出现异常时，会导致神经元损伤加剧和GLU释放增强[49]，即GLU兴奋性毒性，它会加剧神经元的损伤，可能引发原发性或继发性髓鞘破坏[50]。
基于生物信息学分析结果，通过RT-qPCR进一步鉴定了关键基因的表达水平。RT-qPCR结果显示，DDX3X、LAMP1、GLUD1和VDAC1的表达水平与mRNA芯片的生物信息学分析结果一致，这些被证实的基因与前人的研究佐证，在MS的发生发展过程中起着重要的作用。
综上所述，MS的发病可能与DDX3X、VDAC1、LAMP1、GLUD1的表达，参与EBV感染等信号通路，并影响OS和GLU代谢等生物功能有关。不足之处在于数据样本量较少，且仅从正常人与MS患者之间的基因差异角度进行了研究，在MS患者的不同发病阶段情况可能不同。总之，本研究应用生物信息学技术探讨了MS发病的相关分子机制，为今后MS的治疗提供了新的靶点视角。
