摘要
基于生物信息学研究中医非药物疗法治疗坐骨神经痛的物质基础与潜在作用机制，构建相关miRNA-mRNA。通过高通量基因表达数据库（Gene expression omnibus，GEO）检索坐骨神经痛相关芯片，分析得到差异表达基因。通过支持向量机特征消除算法（Support vector machine based recursive feature elimination，SVM-RFE）和随机森林模型算法（Random forest，RF）结合受试者工作特征（Receiver operating characteristic curve，ROC）曲线及箱式图筛选出具有较高诊断效能的坐骨神经痛核心基因。利用Mfuzz时序分析识别中医非药物疗法治疗坐骨神经痛过程中具有相同变化特征的基因群。利用蛋白质互作网络分析数据平台（String）进行蛋白质交互作用（Protein-protein interaction，PPI）分析，将PPI分析结果导入到Cytoscape3.10.2软件中，并通过cytoHubba插件挖掘核心基因。对中医非药物疗法治疗坐骨神经痛的核心基因进行基因本体论（Gene ontology，GO）分析、京都基因与基因组百科全书（Kyoto encyclopedia of genes and genomes，KEGG）及疾病富集分析(Disease ontology enrichment analysis，DO)通路富集分析。利用8个miRNA在线数据库检索坐骨神经痛核心基因及中医非药物疗法治疗坐骨神经痛核心基因上游miRNA，构建miRNA-mRNA调控网络。结果发现癌胚抗原细胞粘附分子4（CEA cell adhesion molecule 4，CEACAM4）、kelch 样家族成员 14（kelch like family member 14 Gene，KLHL14）、盘状蛋白（Discoidin CUB And LCCL domain containing 1，DCBLD1）、双特异性磷酸酶2（Dual specificity phosphatase 2，DUSP2）、IKAROS家族锌指2（IKAROS family zinc finger 2，IKZF2）及神经元细胞粘附分子（Neuronal cell adhesion molecule gene，NRCAM）有可能成为坐骨神经痛核心基因。中医非药物疗法可能通过成纤维细胞生长因子6（Fibroblast growth factor 6，FGF6）、过氧化物酶体增殖物激活受体γ（Peroxisome proliferator-activated receptor gamma，PPARG）、胰岛素样生长因子 1（Insulin like growth factor 1，IGF1）、载脂蛋白E（Apolipoprotein E，APOE）介导丝氨酸/苏氨酸激酶（Phosphatidylinositide 3-kinases/protein kinase B，PI3K/AKT）信号通路与腺苷酸激活蛋白激酶（Adenosine 5'-monophosphate-activated protein kinase，AMPK）信号通路治疗坐骨神经痛。坐骨神经痛核心基因与中医非药物疗法治疗坐骨神经痛核心基因经筛选后共募集到hsa-miR-27a-3p、hsa-miR-27b-3p等21个miRNA。本研究通过生物信息学方法发现了坐骨神经痛核心基因以及中医非药物疗法治疗坐骨神经痛的核心基因与作用机制，为坐骨神经痛的治疗提供了新的思路。本研究还构建了相关miRNA-mRNA调控网络，该调控网络可能在坐骨神经痛的发病及中医非药物疗法治疗中扮演重要角色。
Abstract
Based on bioinformatics, this study aimed to explore potential molecular mechanism of non-pharmacological therapy of traditional Chinese medicine in treating sciatica, as well as the the miRNA-mRNA regulatory network related to sciatica and the material basis. This study retrieved sciatica-related microarray data from the Gene Expression Omnibus (GEO) database and analyzed differentially expressed genes. Support Vector Machine based Recursive Feature Elimination (SVM-RFE) and Random Forest (RF) algorithms were applied to screen out core genes with high diagnostic performance by combining Receiver Operating Characteristic (ROC) curve and boxplot. Mfuzz was utilized to identify gene clusters with similar change characteristics under non-pharmacotherapy of traditional Chinese medicine in treating sciatica. The Protein-Protein Interaction (PPI) analysis was conducted on the String platform, and the results of PPI data were imported into Cytoscape 3.10.2 software. Core genes were then identified using the cytoHubba plugin. Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and Disease Ontology enrichment analyses were performed to obtain the mechanism of non-pharmacological therapy of traditional Chinese medicine in treating sciatica. Using 8 miRNA online databases, this study retrieved the core genes related to sciatica and the core genes in non-pharmacological therapy of traditional Chinese medicine in treating sciatica, and then constructed a critical miRNA-mRNA regulatory network. CEA Cell Adhesion Molecule 4(CEACAM4), kelch like family member 14 Gene(KLHL14)，Discoidin CUB And LCCL Domain Containing 1(DCBLD1),Dual Specificity Phosphatase 2(DUSP2), IKAROS Family Zinc Finger 2(IKZF2), and Neuronal Cell Adhesion Molecule Gene (NRCAM) may be biomarkers for the diagnosis of sciatica. Non-pharmacological therapy of traditional Chinese medicine in treating sciatica may be based on fibroblast growth factor 6(FGF6), Peroxisome Proliferator-Activated Receptor Gamma (PPARG), Insulin like growth factor 1(IGF1), Apolipoprotein E(APOE) mediates the Phosphatidylinositide 3-kinases/protein kinase B (PI3K/AKT) signaling pathway and Adenosine 5'-monophosphate-activated protein kinase(AMPK) signaling pathway in treating sciatica. A total of 21 miRNAs, including hsa-mi-27a-3p and hsa-miR-27b-3p, were recruited from the sciatica core gene and the core gene of non-pharmacological therapy of traditional Chinese medicine in treating sciatica.In this study, bioinformatics were used to successfully screen a number of core genes that may be associated with sciatica. At the same time, we also discovered the core genes and mechanism of action of non-pharmacological therapy of traditional Chinese medicine in treating sciatica, which provides a new idea for the treatment of sciatica. Finally, we constructed a miRNA-mRNA regulatory network, which may play an important role in the pathogenesis of sciatica and the therapeutic mechanism of non-pharmacological therapy of traditional Chinese medicine in treating sciatica.
坐骨神经痛是从臀部沿着坐骨神经的走行向下辐射的疼痛，通常是由椎管或神经孔内的坐骨神经受压引起[1]。研究表明，坐骨神经痛的患病率从1.6%到43%不等[2]，是现代社会的常见病之一，并呈现出逐渐年轻化的趋势，严重影响患者的生活质量，成为了当下临床亟待解决的重要问题[3]。
目前药物疗法在改善坐骨神经痛方面存在服药周期长、依从性差、潜在肝肾损伤等局限性，在一定程度上增加了坐骨神经痛患者的治疗负担[4]，而手术疗法治疗坐骨神经痛的方式也存在诸多后遗症[5]。据报道，许多方法可以改善坐骨神经痛，而90%的急性坐骨神经痛可以通过非手术治疗得到有效缓解[6]。长期以来，中医非药物疗法如针灸、推拿、艾灸等已被证明可有效治疗坐骨神经痛[7]，拥有着预后良好、成本低廉等优点，但具体作用机制尚不明确，仍需进一步探索。
MicroRNA（miRNA）是一种非蛋白质编码RNA，它通过与靶向信使mRNA特异性结合，降解靶向mRNA或抑制其翻译过程，从而参与疾病的发生、发展过程并协调广泛的生物过程[8]。许多miRNA被证明可通过靶向关键转录因子或信号通路影响坐骨神经痛过程，如胡椒碱被证明可通过miR-520a-P65调控网络治疗坐骨神经痛[9]，miR-338-3p-TRPV1调控网络可以诱导神经性疼痛，为神经性疼痛的治疗提供潜在的治疗靶点[10]。
为此，本文拟通过生物信息学研究中医非药物疗法对坐骨神经痛的干预机制，构建相关miRNA-mRNA调控网络，为未来中医非药物疗法治疗坐骨神经痛的研究提供坚实的理论参考。
1 材料与方法
1.1 数据来源及筛选差异表达基因
以“sciatica”为关键词检索高通量基因表达数据库（Gene expression omnibus，GEO）[11]并确定芯片号为GSE150408。该数据集包含了17名无腰痛或坐骨神经痛临床证据的健康志愿者外周血样本，17名经磁共振成像证实的坐骨神经痛患者外周血样本，25名接受了中医非药物疗法的坐骨神经痛患者外周血样本。同时，选择GSE124272数据集，其包含8名无腰痛或坐骨神经痛临床证据的志愿者外周血样本，8名经磁共振成像证实的腰椎间盘脱垂患者的信息[12]。利用R语言中limma包分析GSE150408数据集，将17名无腰痛或坐骨神经痛临床证据的健康志愿者外周血样本作为对照组，17名经磁共振成像证实的坐骨神经痛患者外周血样本作为实验组，结合P值（P.value）和差异倍数（Fold change，FC）进行筛选，将差异表达基因（Differential expressed genes，DEGs）的筛选条件设置为P.value＜0.05、|logFC|＞0.5。
1.2 挖掘坐骨神经痛核心基因
支持向量机递归特征消除（Support vector machine based recursive feature elimination，SVM-RFE）和随机森林模型（Random forest，RF）是生物医学领域中应用最为广泛的两种机器学习算法[13]，对疾病的预测具有良好的准确性[14]。故本研究分别利用SVM-RFE算法和RF算法[15]从DEGs中筛选出重要基因，将以上2种机器学习算法结果的交集基因确定为关键基因。SVM-RFE算法通过十字交叉验证后，保存前40个特征基因结果；RF算法选择ntree=500，对特征基因变量进行逐一筛选。基于GSE124272数据集以及GSE150408数据集，利用R语言中pROC包绘制关键基因受试者工作特征（Receiver operating characteristic curve，ROC）曲线，计算ROC的曲线下面积（Area under curve，AUC）。一般认为，AUC为 0.5表示诊断效能无歧视（即诊断患有和不患有疾病的患者的能力），0.7~0.8表示诊断效能良好，0.8~0.9表示诊断效能优秀，超过0.9为诊断效能杰出[16]。故本研究将AUC值设置为大于0.7，结合箱式图以筛选出坐骨神经痛核心基因。
1.3 Mfuzz时序分析基因调控
基因调控通常不遵循“开关”动态，而是更加细致和渐进，从而实现精细的基因控制功能，Mfuzz是一款基于R语言开发的工具包，其主要功能是对时间序列微阵列数据执行模糊聚类分析[17]。本研究从GSE150408数据集中提取“无坐骨神经痛-坐骨神经痛-中医非药物疗法治疗”作为不同的临床状态进行Mfuzz时序分析，以识别中医非药物疗法治疗坐骨神经痛过程中具有相同变化特征的基因群。
1.4 功能富集分析
利用R语言ClusterProfiler包、org.Hs.eg.db包、enrichplot包等[18]对Mfuzz时序分析结果进行基因本体论（Gene ontology，GO）分析、京都基因与基因组百科全书（Kyoto encyclopedia of genes and genomes，KEGG）通路富集分析以及疾病富集分析（Disease ontology enrichment analysis，DO）[19]，其中GO分析包括三个主要部分：生物学过程（Biological process，BP）、细胞组成（Cellular component，CC）和分子功能（Molecular function，MF）。
1.5 挖掘中医非药物疗法治疗坐骨神经痛核心基因
将Mfuzz时序分析结果导入到蛋白质互作网络（Protein-protein interaction，PPI）分析数据平台（String）（https//cn.string-db.org），进行PPI网络分析[20]，选择蛋白种类为“Homo sapiens”，隐藏掉与其他节点无交互作用的节点后得到PPI网络分析结果，将该结果导入Cytoscape3.10.2软件中，并进行CytoHubba分析[21]，得到中医非药物疗法治疗坐骨神经痛核心基因。
1.6 miRNA-mRNA调控网络构建
利用R语言检索8个在线数据库，分别为ENCORI、miRDB、miRWalk、RNA22、RNAInter、miRTarBase、TargetMiner、TargetScan，对坐骨神经痛核心基因及中医非药物疗法治疗坐骨神经痛核心基因的上游miRNA进行预测，将两者预测到的miRNA取交集，构建miRNA-mRNA调控网络。
2 结果
2.1 坐骨神经痛差异表达基因
将17名无腰痛或坐骨神经痛临床证据的健康志愿者外周血样本作为对照组，17名经磁共振成像证实的坐骨神经痛患者外周血样本作为实验组，通过R语言中limma包对GSE150408数据集进行处理，筛选得到462个DEGs（图1）。
2.2 挖掘坐骨神经痛核心基因
使用RF算法从462个DEGs筛选出了298个基因，使用SVM-RFE算法筛选出了35个基因，二者取交集后筛选出26个基因（图2）。为进一步评价诊断效能，通过GSE124272数据集及GSE150408数据集对上述26个基因进行验证。箱式图结果显示，癌胚抗原细胞黏附分子4（CEA Cell adhesion molecule4，CEACAM4）及盘状蛋白（Discoidin CUB And LCCL Domain containing1，DCBLD1）在两个数据集中的疾病组表达量均比控制组高，kelch 样家族成员 14（Kelch like family member 14 Gene，KLHL14）、双特异性磷酸酶2（Dual specificity phosphatase2，DUSP2）、IKAROS家族锌指2（IKAROS Family zinc finger 2，IKZF2）及神经元细胞粘附分子（Neuronal cell adhesion molecule gene，NRCAM）在两个数据集中的疾病组表达量均比控制组低（图3（a）、图4（a））。受试者工作特征曲线表明，KLHL14、DCBLD1、NRCAM、DUSP2、IKZF2、CEACAM4在两个数据集ROC曲线分析中AUC值均大于0.7，证明该6个基因具有良好的诊断价值，可作为坐骨神经痛核心基因（图3（b）、图4（b））。
图1GSE150408数据集差异表达基因分析结果
Fig.1Analysis results of differentially expressed genes in GSE150408 dataset
注：（a）火山图显示了|logFC|＞0.5，P＜0.05 时实验组和对照组之间的 DEGs，红色代表高表达，蓝色代表低表达；（b）热图显示了所有DEGs，红色代表高表达，蓝色代表低表达，在group轴中红色代表实验组，蓝色代表对照组.
图2两种机器学习算法结果
Fig.2Results of two machine learning algorithms
注：（a）RF算法结果中重要性排名前20基因；（b）SVM-RFE算法筛选结果，当选择35个基因时，模型预测的准确率达到89.2%；（c）SVM-RFE算法与RF算法交集结果.
2.3 Mfuzz时序分析
对GSE150408数据集进行Mfuzz时序分析，结果显示，基因被分为8个簇（Cluster）（图5）。而在原实验中，患者仅接受了为期2周的短期治疗，因此目标基因群的表达量不应恢复到正常人的水平，如簇2、5所示，也不应超出正常人的水平，如簇4、6所示，更不应在三个时间点持续处于降低或者升高的状态，如簇1、7所示，只有簇3、8符合正常的生理变化，因此保留簇3、8基因集以进行后续的分析。
图3GSE150408数据集验证结果
Fig.3Validation results of GSE150408 dataset
注：（a）箱式图；（b）ROC曲线；其中#代表P＜0.2，*代表P＜0.05，**代表P＜0.01，***代表P＜0.001.
2.4 功能富集分析
对簇3及簇8基因集进行GO分析、KEGG分析、DO分析。簇3基因集GO分析共得到572个条目，KEGG分析共得到21条通路，DO分析共得到271个条目。簇8基因集GO分析共得到406个条目，KEGG分析共得到20条通路，DO分析共得到59个条目。将排名前7的生物过程、排名前20的信号通路及排名前20的疾病分析富集结果以气泡图展示。可见，中医非药物疗法治疗坐骨神经痛过程中可能参与的生物学过程包括蛋白质分泌（Protein secretion）、成骨细胞分化（Osteoblast differentiation）、氧化应激反应（Response to oxidative stress）、DNA代谢过程的正调控（Positive regulation of DNA metabolic process），见图6中BP。细胞组分主要有细胞基膜（Apical plasma membrane）、高尔基体（Golgi lumen）、特定颗粒腔（Specific granule lumen）、质膜微囊（Cytoplasmic vesicle lumen），见图6中CC。分子功能主要包括DNA结合转录激活剂活性，RNA聚合酶II-特异性（DNA-binding transcription activator activity，RNA polymerase II-specific）、端粒DNA结合（telomeric DNA binding），见图6中MF。KEGG富集分析显示，中医非药物疗法治疗坐骨神经痛的通路与丝氨酸/苏氨酸激酶（Phosphatidylinositide3-kinases/protein kinase B，PI3K/AKT）信号通路、腺苷酸激活蛋白激酶（Adenosine5'-monophosphate-activated protein kinase，AMPK）信号通路，见图7。DO分析表明中医非药物疗法治疗坐骨神经痛的基因也可能参与调控视网膜疾病（Retinal disease）、肉瘤（Sarcoma）、卵巢癌（Ovarian cancer）、痴呆（Dementia）等疾病，见图8。
图4GSE124272数据集验证结果
Fig.4Validation results of GSE124272 dataset
注：（a）箱式图；（b）ROC曲线； 其中ns代表非显著性差异， # 代表 P＜0.2，* 代表P＜0.05，** 代表 P＜0.01，***代表P＜0.001.
图5Mfuzz分析
Fig.5Mfuzz analysis
注：Normal代表无坐骨神经痛，Disease代表坐骨神经痛，Cure代表中医非药物疗法治疗.
图6GO分析
Fig.6GO analysis
注：（a）簇3基因集分析结果；（b）簇8基因集分析结果.
2.5 挖掘中医非药物疗法治疗坐骨神经痛核心基因
将筛选出的簇3及簇8基因集分别输入String平台，进行PPI网络分析，将该结果导入到Cytoscape3.10.2软件中，采用Cytohubba插件中的Degree、MCC、MNC、DMNC 算法计算评分排名前10的基因，再取其交集。簇3基因集确定了4个核心基因，分别为成纤维细胞生长因子6（Fibroblast growth factor 6，FGF6）、过氧化物酶体增殖物激活受体γ（Peroxisome proliferator-activated receptor gamma，PPARG）、胰岛素样生长因子 1（Insulin like growth factor 1，IGF1）、载脂蛋白E（Apolipoprotein E，APOE），簇8基因集未确定出核心基因，如表1所示。
图7KEGG分析
Fig.7KEGG analysis
注：（a）簇3基因集分析结果；（b）簇8基因集分析结果.
图8DO分析
Fig.8DO analysis
注：（a）簇3基因集分析结果；（b）簇8基因集分析结果.
表1Cytohubba分析
Table1Cytohubba analysis
2.6 miRNA-mRNA调控网络构建
利用R语言检索ENCORI、miRDB、miRWalk等8个在线数据库对坐骨神经痛核心基因及中医非药物疗法治疗坐骨神经痛核心基因上游miRNA进行预测并取交集，结果显示，一共募集到了21个miRNA，用Cytoscape3.10.2软件构建miRNA-mRNA调控网络，如图9所示。
图9miRNA-mRNA调控网络
Fig.9miRNA-mRNA regulatory network
3 讨论
坐骨神经痛大多由椎间盘髓核移位或椎间盘间隙外环纤维化压迫腰神经根引起，最终产生坐骨神经径路及分布区域疼痛 [22-23]，越来越多的证据表明，突出的髓核诱导的氧化应激和炎症反应可能在坐骨神经痛过程中起到更为关键的作用[24]。手术治疗被认为能够有效缓解坐骨神经痛[25]，然而残余背痛和复发性椎间盘突出症是最为常见的术后并发症[26-27]。传统的保守治疗如静脉或皮下注射抗炎药物以及硬膜外类固醇注射在长期缓解坐骨神经痛方面并不理想[28]。故本研究利用生物信息学研究中医非药物疗法治疗坐骨神经痛作用机制及构建miRNA-mRNA调控网络，为精准治疗坐骨神经痛提供更好的分子靶标。
本研究结果显示，hsa-miR-27a-3p、hsa-miR-27b-3p等21个miRNA 在坐骨神经痛的发展及转归中发挥着重要的作用，CEACAM4、DUSP2、IKZF2等6个基因对坐骨神经痛具有较高的诊断价值。CEACAM4是中性粒细胞表达许多表面粘附分子之一，可促使中性粒细胞离开循环系统，参与炎症反应[29]。DCBLD1缺乏酶活性，可作为质膜信号转导中心的支架，能够激活自噬而使代谢异常[30]。在衰老过程中，DUSP2表达量将会提高，从而介导并放大炎症反应，对人体产生不利的功能变化[31]。有研究证实，IKZF2可编码锌指蛋白，该蛋白主要在淋巴细胞中充当转录阻遏蛋白而可导致免疫功能下降并增加促炎细胞因子的产生[32]。NRCAM在神经系统发育中起着重要作用，重点参与周围神经系统中髓鞘神经结构的形成[33]。KLHL14 属于 Kelch 样基因家族，可调控椎间盘中关键成骨转录因子的生成与转化[34]。Mfuzz时序分析及功能富集分析表明，中医非药物疗法可能通过FGF6、PPARG、IGF1、APOE介导PI3K/AKT信号通路及AMPK信号通路从而治疗坐骨神经痛，并且对于治疗视网膜疾病、肉瘤、卵巢癌、痴呆等疾病表现出了较高的潜在价值。有研究证实，FGF6可激活PI3K/AKT信号通路并维持和修复骨骼肌干细胞[34]。PPARG在近期的一项研究中被发现可通过AMPK信号通路参与氧化应激机制以及神经保护作用[35]。Khloud认为IGF1能够通过激活PI3K/AKT和AMPK信号通路参与热量限制，最终延缓衰老[36]。APOE是血浆脂质水平的关键调节剂，可通过激活PI3K/AKT信号通路发挥抗炎作用[37]。软骨细胞中脂肪因子也会由于APOE基因的缺失而异常生长，从而加速了椎间盘退行性改变，最终导致坐骨神经痛[38]。视网膜疾病、肉瘤、卵巢癌和痴呆等疾病被发现均与炎症反应和细胞死亡等基本生命活动有关[39-42]，而PI3K/AKT和AMPK信号通路是调控这些过程的关键途径[43-44]，侧面证实了本研究的准确性。传统的中医理论认为，人体的健康状态与气血、阴阳等概念密切相关，而这些概念在现代医学中往往难以找到对应的实体[45]。然而，通过将中医理论与现代生物学技术相结合，有可能揭示出中医非药物疗法背后的生物学机制，从而为其提供科学的依据[46]。
4 结论
1）成功预测出与坐骨神经痛相关的核心基因及中医非药物疗法治疗坐骨神经痛的核心基因及相关通路，以上基因和通路与衰老及炎症反应密切相关，提示除了直接针对神经的治疗外，调节衰老和炎症过程也将是治疗坐骨神经痛的有效策略。
`2）成功构建了坐骨神经痛miRNA-mRNA调控网络，这一发现为坐骨神经痛的研究提供了新的视角，而目前对于该miRNA-mRNA调控网络与坐骨神经痛相关性的研究仍处于空白，未来的研究工作应该进一步深入探讨其具体作用机制，以及它们与其他生物过程的相互关系，以期找到更多有效的治疗手段。
