丙戊酸單藥治療癲癇患者療效的生物信息學(xué)分析

陳傲寒綜述， 黎銀潮， 陳樹(shù)達(dá)， 劉賢岳， 許曉偉， 趙 軻， 王誠(chéng)蔗， 周列民審校,2

癲癇是最常見(jiàn)的慢性神經(jīng)元性疾病之一，以神經(jīng)元持續(xù)性的異常放電導(dǎo)致中樞神經(jīng)系統(tǒng)功能失常為特征。任何年齡均會(huì)發(fā)病，全世界約7000萬(wàn)患者，其反復(fù)發(fā)作影響患者的身體健康，對(duì)患者及其家人的生活、心理及經(jīng)濟(jì)等方面造成負(fù)擔(dān)[1]。使用抗癲癇藥物(antiepileptic drugs，AEDs)是控制癲癇發(fā)作最基礎(chǔ)、最有效的治療方式，但是臨床上仍有30%的患者在足量、足療程治療后復(fù)發(fā)[2]。丙戊酸是一種非選擇性鈉離子通道阻滯劑，是治療癲癇的一線藥物，對(duì)控制癲癇全身發(fā)作和局部發(fā)作具有較好的效果，但是這種效果受基因變異的影響較大。Feng[3]等發(fā)現(xiàn)CYP3A4、SCN1A和UGT2B7基因多態(tài)性對(duì)丙戊酸治療癲癇的療效具有重大影響。為了從基因角度更深入的分析丙戊酸治療癲癇存在差異性的可能，本研究使用來(lái)自基因芯片公共數(shù)據(jù)庫(kù)GEO的一組癲癇患者與正常人的芯片，分析探討基因表達(dá)的差異性對(duì)丙戊酸治療癲癇療效的影響。

1 材料和方法

1.1 數(shù)據(jù)來(lái)源 本研究選取數(shù)據(jù)來(lái)自GEO數(shù)據(jù)庫(kù)的芯片，系列號(hào)GSE143272，采用商業(yè)化的離子通道芯片平臺(tái)GPL10558[4]。數(shù)據(jù)集包含16個(gè)癲癇患者丙戊酸單藥治療反應(yīng)好、9個(gè)癲癇患者丙戊酸單藥治療無(wú)效和50個(gè)正常人的外周血標(biāo)本。從這些外周血標(biāo)本中分離得到總RNA，再利用Illumina芯片技術(shù)和HumanHT-12_v.4 芯片進(jìn)行基因表達(dá)譜分析。

1.2 預(yù)處理數(shù)據(jù)與篩選差異表達(dá)基因 基于GPL10558平臺(tái)構(gòu)建的信息，將探針識(shí)別號(hào)轉(zhuǎn)化為正式基因名稱(保留mRNA探針，放棄其他非mRNA探針)。利用limma軟件包分析基因表達(dá)矩陣，經(jīng)affy、affyPLM軟件包進(jìn)一步處理，得到DEGs。采用Benjamini & Hochberg錯(cuò)誤發(fā)現(xiàn)率的方法降低假陽(yáng)性率[5～7]，取P值≤0.05，|logFC|>0作為截?cái)嘀档臉?biāo)準(zhǔn)，logFC>0和logFC<0分別對(duì)應(yīng)上調(diào)和下調(diào)的差異表達(dá)基因。用火山圖表示篩選出的DEGs，并作Venny圖取非交集進(jìn)一步選取DEGs。

1.3 信號(hào)通路富集和差異表達(dá)基因本體論(Gene ontology)分析 目前，基因本體論分析(GO)作為一種對(duì)基因產(chǎn)物及其功能特征進(jìn)行注釋的常用方法，已經(jīng)在分子生物學(xué)等領(lǐng)域得到廣泛應(yīng)用，在有效地鑒定高通量相關(guān)的遺傳數(shù)據(jù)的生物學(xué)屬性方面，GO分析包括細(xì)胞組分(cellular components，CC)、生物學(xué)過(guò)程(biological processes，BP)和分子功能(molecular function，MF)[8]。KEGG數(shù)據(jù)庫(kù)是來(lái)自日本的一個(gè)開(kāi)放性信息源數(shù)據(jù)庫(kù)，其數(shù)據(jù)的來(lái)源是基于RNA-seq實(shí)驗(yàn)和微陣列，常用于注釋所涉及的信號(hào)通路網(wǎng)絡(luò)和基因列表，而且是處理細(xì)胞、基因組、信號(hào)轉(zhuǎn)導(dǎo)通路和疾病等的常用數(shù)據(jù)庫(kù)[9]。利用DAVID數(shù)據(jù)庫(kù)對(duì)DEGs進(jìn)行KEGG/生物途徑富集和GO分析，該網(wǎng)站是功能注釋生物信息學(xué)微陣列分析的相關(guān)網(wǎng)站，本研究通過(guò)使用ggplot2包將療效好與差的DEGs的功能分析(細(xì)胞組分、生物學(xué)過(guò)程、分子功能和以及信號(hào)通路)可視化，其中錯(cuò)誤發(fā)現(xiàn)率(false discovery rate，F(xiàn)DR)<0.05說(shuō)明統(tǒng)計(jì)學(xué)具有顯著的差異性[10]。

1.4 蛋白質(zhì)-蛋白質(zhì)相互作用(Protein-Protein interaction，PPI)網(wǎng)絡(luò)分析 通過(guò)數(shù)據(jù)庫(kù)STRING(The Search Tool for the Retrieval of Interacting Genes，STRING，version 11.0)可以評(píng)估蛋白質(zhì)之間的相互作用，該數(shù)據(jù)庫(kù)覆蓋了約2460萬(wàn)個(gè)蛋白質(zhì)，其顯示各種蛋白質(zhì)之間相互作用，可用于發(fā)現(xiàn)蛋白質(zhì)之間的直接(物理)和間接(功能)關(guān)聯(lián)[11]。將界值標(biāo)準(zhǔn)設(shè)定為最大相互作用數(shù)=0且置信度得分≥0.4，進(jìn)而評(píng)估DEGs的潛在蛋白質(zhì)-蛋白質(zhì)相互作用。隨后，用Cytoscape的插件Cytohubba分析PPI網(wǎng)絡(luò),使用MCC算法，從PPI網(wǎng)絡(luò)中挑選出與周圍基因具有高度連通性前5個(gè)基因，作為核心基因[12]。

1.5 差異表達(dá)基因功能模塊分析 使用Cytoscape中的分子復(fù)合物檢測(cè)(MCODE)插件，其根據(jù)拓?fù)潢P(guān)系對(duì)給定網(wǎng)絡(luò)進(jìn)行聚類，找到密集連接的區(qū)域。本研究分析基于DEGs的PPI網(wǎng)絡(luò)以篩選出基因功能模塊，node score cutoff=0.2，degree cutoff=2，max depth=100以及k-core=2。另外將篩選出的功能模塊映射到STRING數(shù)據(jù)庫(kù)中，進(jìn)行KEGG和GO分析以注釋信號(hào)傳導(dǎo)途徑和功能模塊組分的基因。

1.6 統(tǒng)計(jì)學(xué)方法 采用t檢驗(yàn)進(jìn)行DEGs的篩選以及采用Fisher精確檢驗(yàn)分析GO和KEGG注釋富集。所有統(tǒng)計(jì)分析在R版本4.0.3軟件中執(zhí)行。

2 結(jié) 果

2.1 差異表達(dá)基因 通過(guò)對(duì)GEO數(shù)據(jù)集GSE157797進(jìn)行差異基因篩選，以調(diào)整P值≤0.05，|logFC|>0為標(biāo)準(zhǔn)得出DEGs，其中在丙戊酸單藥治療有效組與正常人組對(duì)比的1836個(gè)差異表達(dá)基因中，分別包括972下調(diào)基因和864個(gè)上調(diào)基因(見(jiàn)圖1A)。在丙戊酸單藥治療無(wú)效組與正常人組對(duì)比的751個(gè)差異表達(dá)基因中，分別包括486下調(diào)基因和265個(gè)上調(diào)基因(見(jiàn)圖1B)。用火山圖表示如圖1A、圖1B,紅色表示上調(diào)的DEGs，藍(lán)色表示下調(diào)的DEGs。然后做Venny圖取非交集進(jìn)一步篩選出療效差異的相關(guān)基因。

圖1 GSE157797數(shù)據(jù)集基因表達(dá)水平

圖2 差異表達(dá)基因Venny圖

2.2 差異表達(dá)基因的功能分析 為了進(jìn)一步了解丙戊酸單藥治療有效和無(wú)效的DEGs的生物學(xué)功能，分別將對(duì)VPA治療反應(yīng)好的1508個(gè)基因和對(duì)VPA治療反應(yīng)差的423個(gè)基因進(jìn)行GO和KEGG/生物途徑富集分析。

GO分析結(jié)果顯示，癲癇患者VPA單藥治療效果反應(yīng)好的1508個(gè)DEGs的細(xì)胞組分在細(xì)胞質(zhì)、細(xì)胞核和細(xì)胞膜中均有富集；生物過(guò)程主要是參與蛋白翻譯，T細(xì)胞受體信號(hào)通路和核轉(zhuǎn)錄的mRNA分解過(guò)程；分子功能富集在蛋白質(zhì)以及RNA結(jié)合過(guò)程中。KEGG/生物途徑分析顯示，主要參與RNA轉(zhuǎn)運(yùn)和核糖體相關(guān)信號(hào)通路。423下調(diào)基因DEGs的細(xì)胞組分在細(xì)胞核中富集；生物過(guò)程主要是高爾基囊泡介導(dǎo)的轉(zhuǎn)運(yùn)和以DNA為模板的轉(zhuǎn)錄調(diào)控過(guò)程中；分子功能富集在DNA、蛋白質(zhì)的結(jié)合以及泛素蛋白連接酶活性。KEGG/生物途徑分析顯示，主要參與mRNA監(jiān)測(cè)途徑和溶酶體相關(guān)信號(hào)通路(見(jiàn)圖3、圖4)。

圖4 DEGs的GO分析

2.3 核心基因的功能模塊及表達(dá)水平分析

2.3.1 PPI 網(wǎng)絡(luò)的構(gòu)建及10個(gè)核心基因的篩選 為了篩選出對(duì)丙戊酸反應(yīng)存在差異的基因，首先構(gòu)建PPI網(wǎng)絡(luò)，再使用Cytoscape計(jì)算PPI網(wǎng)絡(luò)，最后分別參照MCC方法各選取前5個(gè)基因，并把這些基因命名為核心基因(見(jiàn)圖5)。

圖5 前5個(gè)核心基因

2.3.2 功能模塊分析 使用Cytoscape中的MCODE插件和DVAID網(wǎng)站分析DEGs功能模塊，其按照基因計(jì)數(shù)對(duì)結(jié)果進(jìn)行降序排序，P值<0.05(見(jiàn)表1)。顯示核心功能模塊中對(duì)VPA治療反應(yīng)好的102個(gè)DEGs(見(jiàn)圖6)主要富集在UBC13-MMS2 復(fù)合物和U6小核核糖核蛋白中，其涉及的生物學(xué)過(guò)程集中在多聚體結(jié)合的調(diào)節(jié)和蛋白質(zhì)去甲基化的負(fù)調(diào)控,生物分子功能分析顯示集中在泛素連接酶抑制劑活性，組蛋白前mRNA DCP結(jié)合和蛋白翻譯激活劑活性，KEGG信號(hào)通路分析核心模塊的DEGs主要參與核糖體和剪接體的相關(guān)信號(hào)通路。對(duì)VPA治療反應(yīng)差102個(gè)DEGs主要富集在UBC13-MMS2復(fù)合物和U6小核核糖核蛋白中，其涉及的生物學(xué)過(guò)程集中在多聚體結(jié)合的調(diào)節(jié)和蛋白質(zhì)去甲基化的負(fù)調(diào)控,生物分子功能分析顯示集中在泛素連接酶抑制劑活性，組蛋白前mRNA DCP結(jié)合和蛋白翻譯激活劑活性，KEGG信號(hào)通路分析核心模塊的DEGs主要參與核糖體和剪接體的相關(guān)信號(hào)通路，而在與平行對(duì)照的腦組織相比，這些基因的表達(dá)改變，顯示了探究差異表達(dá)基因的相關(guān)生物學(xué)功能及信號(hào)通路對(duì)研究VPA療效具有重要價(jià)值。

表1 DEGs的功能模塊相關(guān)分析

圖6 DEGs的功能模塊

3 總結(jié)與討論

丙戊酸是一種廣譜抗癲癇藥物，廣泛應(yīng)用于全身性和部分性癲癇的治療。其療效有顯著的個(gè)體差異性，臨床和分子因素是丙戊酸療效個(gè)體差異的基礎(chǔ)。臨床因素包括性別、年齡、發(fā)病機(jī)制、肝腎功能、藥物相互作用等。除去這些臨床因素，對(duì)丙戊酸的反應(yīng)在患者間也存在顯著的差異。潛在的分子因素，例如丙戊酸影響受體、轉(zhuǎn)運(yùn)體和藥物代謝酶的基因表達(dá)，單獨(dú)或結(jié)合編碼這些基因的多態(tài)性。既往的研究報(bào)道關(guān)于丙戊酸反應(yīng)差異的分子因素主要集中在藥物轉(zhuǎn)運(yùn)與代謝的差異基因上，例如UDP糖基轉(zhuǎn)移酶(UDP glycosyltransferase,UGT)和細(xì)胞色素(Cytochrome P450,CYP)等[13]。在本項(xiàng)研究中，發(fā)現(xiàn)了與以往關(guān)注點(diǎn)不同的影響基因，如KLHL21、FBOX30、LRR1、FBXL15、KCTD6、SKP1、MEX3C、UBR4、UBE2C、FBXW7，這些基因編碼蛋白參與人體內(nèi)泛素化過(guò)程。蛋白質(zhì)泛素化是一種由泛素激活酶(E1s)、泛素結(jié)合酶(E2s)和泛素連接酶(E3s)催化的高度有序多步酶級(jí)聯(lián)反應(yīng)。由Clu(Cullin)蛋白作為支架蛋白與SKP1和F-box蛋白構(gòu)成的SCF(Skp1-Cullin1-F-box)復(fù)合體也稱為CRL復(fù)合體，占人類E3s的絕大部分。F-box可特異性識(shí)別底物，F(xiàn)BOX30通過(guò)對(duì)RARγ(RA exerts its functions through RA receptors,RARs)的干預(yù)拮抗BMP(bone morphogenic protein,BMP)通路[14,15],與此同時(shí)FBXL15亦可負(fù)向調(diào)控Smurf1穩(wěn)定性來(lái)正向調(diào)控BMP[16],二者共同影響了胚胎神經(jīng)系統(tǒng)發(fā)育；FBXW7是一種腫瘤抑制因子，參與多種信號(hào)通路轉(zhuǎn)導(dǎo)，影響腫瘤發(fā)生[17]。SKP1(S-phase kinase-associated protein1,SKP1)連接Clu和F-box，其表達(dá)降低破壞SCF結(jié)構(gòu)穩(wěn)定性。在神經(jīng)退行性疾病如帕金森患者的神經(jīng)元中發(fā)現(xiàn)SKP1下降，推測(cè)損害神經(jīng)遞質(zhì)傳遞相關(guān)蛋白功能，加劇神經(jīng)元內(nèi)聚集物形成[18,19]。MEX3C作為一種RNA結(jié)合型E3s與某些mRNA后阻止其翻譯，并用泛素化的方式誘導(dǎo)其降解[20,21]。UBE2C(Ubiquitin-conjugating enzyme E2C)編碼另一類E3s蛋白成員，引導(dǎo)多聚泛素化到底物中的賴氨酸，在調(diào)控細(xì)胞周期中發(fā)揮重要作用[22]。UBR4(Ubiquitin protein ligase E3 component N-recognin 4)在人腦與脊髓組織中高度表達(dá)，尤其在海馬錐體神經(jīng)元富集，表明該蛋白在建立影響學(xué)習(xí)和記憶過(guò)程的神經(jīng)元網(wǎng)絡(luò)中起著關(guān)鍵作用[23]。KLHL21為Kelch-like基因家族的一員，特異性結(jié)合IKK復(fù)合物中的IKKβ亞基，促炎刺激時(shí)在巨噬細(xì)胞中表達(dá)下調(diào)，減少促炎因子表達(dá)[24]。此外，KLHL21和Clu3相互作用[25]，介導(dǎo)體外aurora B激酶的泛素化。KCTD6是鉀離子通道蛋白四聚體家族成員，其BTB(Bric-a-brack,Tram-track,Broad complex)可以高親和力的結(jié)合Clu3[26]。LRR1蛋白則與CRL2構(gòu)成CRL2LRR1復(fù)合物[27]參與泛素化的下一個(gè)過(guò)程。蛋白泛素化影響中樞神經(jīng)系統(tǒng)離子通道蛋白如鈉通道(Navs)，已證實(shí)E3s作用于Nav1相似的PY模序來(lái)修飾蛋白，而丙戊酸正是通過(guò)阻滯Navs來(lái)發(fā)揮抗癲癇作用。本研究篩選出的10個(gè)基因影響了大腦內(nèi)尤其是Navs蛋白泛素化修飾過(guò)程，故我們推測(cè)人群中蛋白泛素化相關(guān)基因差異是丙戊酸療效差異性的分子因素之一。我們分析癲癇患者與正常人外周血基因表達(dá)情況，利用生物信息學(xué)分析的方法，提示丙戊酸療效差異性大的另一種可能，為指導(dǎo)用藥及精準(zhǔn)治療提供新思路。