丙戊酸鈉抗癲癇療效相關樞紐基因鑒定

楊 炯，張明華，耿 淼，艾倫娜，田 磊，范 皎

(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信號通路、免疫調控等生物學過程影響癲癇患者治療效果。本研究為揭示癲癇的耐藥機制、進一步提升丙戊酸鈉的臨床療效、減少患者的不良反應提供了藥效學理論依據。