丙戊酸钠抗癫痫疗效相关枢纽基因鉴定

杨 炯，张明华，耿 淼，艾伦娜，田 磊，范 皎

(1.中国人民解放军总医院医疗保障中心，2.中国人民解放军总医院第二医学中心老年医学研究所，国家老年疾病临床医学研究中心，北京 100853)

癫痫是一种多因素参与的神经系统疾病，全世界约有6 000多万人受到影响。当前主要依靠药物对症治疗，通过控制患者癫痫发作频率改善患者的生活质量。患者对于抗癫痫药物(antiepileptic drug，AED)的反应存在个体差异，约有40%～50%的患者对第一次AED单药治疗无效[1]，其中30%的患者有抗药性[2]。有报道认为，早期发病、治疗前高发作频率、隐原性癫痫、脑神经解剖学异常等临床因素是药物反应性差的原因。也有研究发现，AED药物的代谢酶、转运体和靶点等分子的遗传多态性可用于预测药物治疗的效果[3]。神经环路由众多未知细胞和物质类型构成，癫痫的多因素特性涉及基因组和环境的复杂相互作用，任何遗传或环境因素的改变最终都会影响到基因的表达和功能，需要超越单一通路的范畴去探究疾病的信号传导机制，进一步为开发有效的药物靶标提供必要依据[4]。

目前，用于治疗癫痫的药物有二十多种，包括苯妥英钠、卡马西平和丙戊酸钠等。其中丙戊酸钠无论单用还是多药联用对多种类型的癫痫都有较好的疗效，目前是抗癫痫的一线用药。丙戊酸钠缓释片属于广谱抗癫痫药物，抗癫痫作用可能是通过竞争性抑制γ-氨基丁酸转移酶，提高脑内γ-氨基丁酸的含量，抑制神经元异常放电实现的。该药的安全性比较高，能够在维持药物浓度的情况下，防止患者出现药物不良反应。但是丙戊酸钠的抗癫痫与其他抗癫痫药物一样，其药效和药代动力学都存在个体差异，原因尚不清楚。

加权基因共表达网络分析(weighted gene co-expression network analysis, WGCNA)[5-6]依据基因的表达模式相似性构建权重网络，进一步通过分析基因模块与性状数据之间的相关性得到与表型高度相关的基因集合，最后通过模块内部基因关联分析得到模块内部的关键基因。WGCNA通过选择合适的加权系数对基因之间的相关系数加权，使之满足基因网络近似服从无尺度网络分布，拥有更好的统计功效，避免了大量的多重校正导致的假阴性结果。本研究通过对未接受过药物治疗的癫痫患者及接受丙戊酸钠治疗并出现不同疗效患者的外周血样本表达谱数据进行WGCNA分析，为探寻造成丙戊酸钠疗效差异相关的基因提供理论依据。

1 材料与方法

1.1 数据获取与预处理从GEO数据库(Gene Expression Omnibus)中下载表达谱芯片数据，数据编号为GSE143272，测序平台为Illumina Human HT-12 V4.0 expression beadchip。从该样本集中选取34例此前未接受过药物治疗的癫痫患者(drug-free组)以及25例服用丙戊酸钠并随访一年以上的患者(Valproate组)的外周血RNA样本。其中又根据药物治疗随访过程中癫痫是否发作将25例患者分为丙戊酸钠治疗有效组患者(VA responders组)16例，无效组患者(VA non-responders组)9例。有效组指治疗过程中无癫痫发作，无效组指用药期间至少经历3次以上癫痫发作。数据集共包括13 165个基因探针表达数据，经过合并和转换筛选生成 9 924 个基因和59个样本的表达矩阵，从9 924个基因的表达数据中筛选出SD值前5 000的基因进行后续分析。

1.2 加权基因共表达网络的构建使用R语言中的WGCNA程序包构建共表达网络。对样本进行聚类分析，根据样本聚类的距离鉴定是否存在离群样本。共表达网络符合无尺度网络，即出现连接度为k的节点的对数lgk与该节点出现的概率的对数lg [P(k)]呈负相关，且相关系数应>0.8。按照无尺度网络标准筛选合适的连接函数加权参数，即软阈值β。由软阈值计算网络拓扑重叠矩阵(topological overlap matrix, TOM)，通过动态剪切树法进行模块识别，设置每个基因网络模块最少的基因数目为30。

1.3 共表达模块与临床表型的相关性分析采用分层聚类方法识别基因模块，并用不同颜色来表示。计算每个模块的特征向量基因ME(module eigengene)，降维处理模块聚类，并合并相似模块，grey模块是无法聚集到其他模块的基因集合。通过计算模块和临床表型性状相关系数，给出模块和性状之间的相关系数热图。本研究中关注的临床表型有年龄(age)，性别(sex)，癫痫种类(epilepsy type)、是否接受丙戊酸钠治疗(treatment)、丙戊酸钠治疗是否有效(VA response)。通过模块的特征向量与性状的相关系数以及模块显著性P<0.05筛选出与临床表型显著相关的基因模块。计算模块内的基因表达与性状的相关性(gene significance，GS)和某个基因表达与模块内基因主成分表达的相关系数(module membership，MM)，通过设置GS、MM、q.weighted的取值范围对networkScreening函数计算得到的基因列表进行筛选，从而识别和鉴定出关键枢纽基因(hub genes)。

1.4 网络可视化及功能富集分析选取与临床表型显著相关的基因模块，利用Cytoscape软件绘制网络图。基于基因本体论(gene ontology, GO)数据库对与性状高度相关的模块分别进行基因功能注释分析，采用Fisher检验筛选出模块基因显著性富集的GO条目；基于Kyoto Encyclopedia of Genes and Genomes(KEGG)数据库对核心模块中基因进行信号通路(pathway)富集分析。

2 结果

2.1 加权基因共表达网络的构建首先对所有样本的基因表达值进行聚类分析并作图(Fig 1A)。通过WGCNA算法，根据无尺度网络拟合指数和平均连接度，计算并选取β=6作为本数据集的软阈值(Fig 1B)，计算基因间的邻接矩阵和TOM矩阵，并根据TOM矩阵构建基因间的分层聚类树，基于动态剪切树的方法把基因分成12个模块，分别用12种颜色矩形表示，纵坐标为基因占比(Fig 1C)。各个模块包含基因个数如下：black(177)，blue(675)，brown(540)，green(199)，greenyellow(59)，magenta(136)，pink(175)，purple(86)，red(181)，turquoise(1897)，yellow(322)，grey(553)，其中无法聚类到其他任何模块中的基因放入grey模块，在后续分析中将grey模块移除。

Fig 1 Construction of co-expression network

2.2 共表达模块与临床表型的相关性分析计算不同基因模块与临床表型之间的关系，绘制共表达模块与表型的相关性热图(Fig 2A)。与是否接受丙戊酸钠治疗(treatment)表型相关性最强的模块是yellow模块(r=0.57，P<0.000 1)，yellow模块中的基因总体上与治疗方式表型正相关，即yellow模块整体上在服用丙戊酸钠治疗组患者(valproate组)中高表达，在此前未接受过药物治疗的患者(drug-free 组)中低表达。与丙戊酸钠治疗有效(VA response)表型相关性最强的模块是blue模块(r=-0.53，P<0.000 1)，blue模块中的基因总体上与丙戊酸钠治疗有效表型负相关，即blue模块整体上在丙戊酸钠治疗有效组患者(VA responders组)中低表达，在无效组患者(VA non-responders组)中高表达。对各个模块进行层次聚类和用热图分析各模块之间的相关性(Fig 2B)。

Fig 2 Correlation between co-expression modules and clinical traits

将yellow和blue模块作为关键模块进行GS和MM分析，yellow模块的GS与MM的相关系数r=0.42，P<0.000 1(Fig 3A)，blue模块的GS与MM的相关系数r=0.29，P<0.000 1(Fig 3B)。进一步分析yellow模块和blue模块的eigengene表达热图(Fig 3C、3D)及yellow模块和blue模块基因与各表型之间的热图和聚类图(Fig 3E、3F)。

Fig 3 Correlation between yellow and blue modules with clinical traits

在两个模块中通过基因显著性(GS)和模块成员关系(MM)筛选yellow和blue模块中关键枢纽基因。设置3个筛选标准：|GS|>0.4，|MM|>0.2，q.weighted<0.05。yellow模块筛选出14个枢纽基因(Tab 1)，blue模块筛选出10个枢纽基因(Tab 2)。

Tab 1 Hub genes screened in yellow module

Tab 2 Hub genes screened in blue module

2.3 网络可视化及功能富集分析分别根据两个模块中筛选出的枢纽基因权重做出共表达网络基因间的相互作用关系图(Fig 4)，其中节点形状的大小代表与该节点连接的边的数量，边的宽度代表两节点间连接的权重大小。

Fig 4 Hub gene mapA: Hub gene map of yellow module; B： Hub gene map of blue module.

分别对两个模块中的枢纽基因进行GO富集分析(Fig 5A、5B)，发现yellow模块的基因功能主要富集于：免疫反应调控(GO:0050776)、细胞表面受体信号通路(GO:0007166)、T细胞受体信号通路(GO:0050852)、MAPK信号转导(GO:0000188)；blue模块的基因功能主要富集于：细胞增殖调控(GO:0042127)、脂质糖反应(GO:0032496)、GTPase活性激活(GO:0090630)、吞噬作用(GO:0006909)、神经元凋亡(GO:0051402)等生物学过程。此外，KEGG通路分析显示，yellow模块基因主要富集于：自然杀伤细胞介导的细胞毒性(hsa04650)、Ras信号通路(hsa04014)、细胞黏附分子(hsa04514)、初级免疫缺陷(hsa05340)(Fig 5C)；blue模块基因主要富集于：NF-κB信号通路(hsa04064)、细胞周期(hsa04110)、抗生素生物合成(hsa01130)等通路(Fig 5D)。

Fig 5 GO and KEGG enrichment analysis of yellow and blue modules

3 讨论

癫痫的发作是一种由多因素共同参与的复杂过程，机制较为复杂，遗传因素、脑部疾病等均可诱发该疾病。从该病临床特点上分析，发作一般没有征兆，有可能导致患者溺水、摔倒和烫伤等意外情况，威胁患者生命安全，影响生活质量。目前尽管已有20多种药物应用于抗癫痫临床治疗，但仍有约1/3的患者药物治疗无效。分析单个基因研究癫痫发病机制难以实现突破，基于网络的分析方法可能有助于发现疾病相关的基因网络。在传统的基因水平研究更多关注于强效应基因，本研究采用的WGCNA算法补充了对弱效应基因的分析。WGCNA算法构建的基因网络关系服从近似无尺度网络分布，因此相较于常规的聚类方法，更具有生物学数据特性，很好还原基因在生物学过程中的作用，有助于鉴定出与特定临床表型相关的重要基因模块和关键枢纽基因。

在本研究中，对癫痫患者外周血样本的表达谱数据集进行了WGCNA分析，根据丙戊酸钠的治疗效果，将患者分为丙戊酸钠治疗有效组患者和无效组患者，同时和未经药物治疗的癫痫患者一起构建共表达网路分析基因表达差异。在共表达网络分析中，识别并聚类成12个颜色模块，并对各模块进行基因与临床表型的相关性分析，其中yellow模块中的基因与治疗方式表型的相关性最为明显，也即yellow模块中的基因在接受丙戊酸钠治疗后的患者中会发现明显变化；blue模块中的基因与丙戊酸钠疗效表型的相关性最为显著，也即blue模块的基因在丙戊酸钠有效患者和无效患者中表达存在明显差异。有研究报道，MAPK信号通路调控神经元的基因表达调控癫痫和认知障碍中的突触兴奋，与癫痫的发作有密切联系[7]。p38 MAPK主要由炎症因子和环境压力激活，抑制p38 MAPK可以减少癫痫诱导的锥体细胞缺失、细胞变性、神经元损伤和海马神经元退化[8]。本研究发现，yellow模块中的基因部分富集到MAPK信号失活通路，提示丙戊酸钠可能通过抑制MAPK信号通路来发挥抗癫痫的作用。NF-κB广泛存在于神经细胞中，参与免疫应激反应、炎症反应及细胞的增殖与凋亡。在癫痫患儿的外周血单个红细胞中NF-κB活化与癫痫发作的严重程度成正相关[9]。也有研究显示NF-κB的活化与癫痫易感性有关。本研究分析显示与丙戊酸钠疗效明显相关的blue模块基因富集在NF-κB通路，癫痫患者体内NF-κB通路的活化可能影响丙戊酸钠的药物敏感性。此外，yellow模块和blue模块基因明显富集在T细胞受体信号通路和T细胞激活机制，这可能与丙戊酸钠对组蛋白乙酰化的修饰能力有关；有研究在肿瘤细胞治疗过程中，利用丙戊酸钠的这种活性增强治疗效果[10]。综合来看，丙戊酸钠抗癫痫的药理作用与其对免疫反应的调控机制存在密切的关联。

根据共表达网络的显著性，本研究在yellow模块中筛选出包括S1PR5、GNLY、CD160在内的14个枢纽基因，在blue模块中筛选出包括MAGED1、FBXO31在内的10个下调枢纽基因。S1PR5是一种淋巴细胞表面的G-蛋白偶联受体，主要表达于NK细胞表面，S1PR5表达缺失影响稳态及炎症条件下小鼠体内自然杀伤细胞(natural killer, NK)的分布[11]。有研究发现，癫痫患者发作期间外周血中NK细胞活性低下，但发作后外周血中NK细胞水平增高[12]。早期研究报道，部分癫痫患者服用苯妥英钠后NK细胞活性降低，而服用卡马西平的癫痫患者NK细胞活性增加[13]。目前，癫痫及抗癫痫药物对NK细胞活性的研究结果尚不完全清楚。本研究提示，服用丙戊酸钠的癫痫患者相比于未经药物治疗的患者S1PR5基因水平出现明显变化，有可能会调控NK细胞的分布或功能。

MAGED1是黑色素瘤相关抗原家族成员之一，广泛表达于各组织中，参与细胞周期调控、细胞凋亡、分化、细胞黏附等多个生物学过程[14]。已有研究发现MAGED1调节脂肪前体细胞发育进而参与糖脂代谢调控途径[15]。FBXO31基因编码的蛋白属于F-box蛋白家族的一员，可以通过靶向介导细胞周期蛋白CyclinD1的泛素化和降解来抑制细胞增殖和细胞周期进程[16]。此外，FBXO31通过降解MDM2在DNA损伤修复过程中发挥重要作用[17]。

丙戊酸钠是抑制癫痫患者发作常用药物，但是患者之间存在较大的个体差异。临床上通常从药代动力学的角度，利用治疗药物检测技术，检测其稳态血药浓度调整给药剂量的方式，实现个体化的用药。如果仍然不能有效控制癫痫发作，一般联合使用拉莫三嗪、加巴喷丁、卡马西平、托吡酯、苯巴比妥等药物。目前较少从药效动力学角度考虑提高个体化治疗水平。本研究通过WGCNA分析，从药效学的角度发现了与丙戊酸钠疗效发生高度相关的模块，包含MAGED1、FBXO31等关键枢纽基因，可能通过调节NF-κB信号通路、免疫调控等生物学过程影响癫痫患者治疗效果。本研究为揭示癫痫的耐药机制、进一步提升丙戊酸钠的临床疗效、减少患者的不良反应提供了药效学理论依据。