引用本文: 張迎春, 何誠成, 段茜萍. 特發性癲癇患者外周血差異表達基因篩選及相關生物信息學分析. 癲癇雜志, 2020, 6(6): 498-506. doi: 10.7507/2096-0247.20200082 復制
癲癇(Epilepsy)是一種常見的嚴重慢性腦部疾病,全球共約 7 000 萬人受到該病的影響[1],流行病學調查顯示其發病率為 0.5%~0.7%,我國約有 650~910 萬癲癇患者,該病可見于各個年齡組,青少年和老年是癲癇發病的兩個高峰年齡段[2]。而特發性癲癇是由基因突變和某些先天因素所致,有明顯遺傳傾向的具體病因尚不清楚的癲癇[2]。人類遺傳物質的變異類型主要有基因突變(gene mutation)、染色體結構或數目變異,這些遺傳變異均可導致癲癇的發生[3]。丙戊酸鈉(VPA)作為一種廣譜抗癲癇藥物(AEDs),可用于多種發作類型的癲癇[2]。在治療全面性癲癇或未分類癲癇時,VPA是最為有效的 AEDs[4]。目前癲癇尚無好的可用于其早期診斷與預后判斷的分子標記物。大多數患者需要長期使用 AEDs 治療。因為疾病的反復發作性,癲癇患者多缺乏自信、情緒消極,嚴重影響學習、工作和生活。近年來,利用高通量測序等方法獲得了大量生物數據,而分析不同情況下的差異表達基因(Differentially expressed genes,DEGs),能夠很好篩選新的有效的分子靶點。目前數據庫中有使用基因芯片分析癲癇患者外周血基因改變的研究結果,我們通過未用AEDs 治療的特發性癲癇外周血樣本與健康對照樣本比對,及VPA治療有效的特發性癲癇外周血樣本與未用藥患者比對,綜合相關結果,試圖尋找癲癇患者外周血關鍵基因及通路改變,以進一步闡明特發性癲癇患者的發病及VPA治療相關機制,為臨床上可用的癲癇早期診斷、預后判斷和治療的靶點提供線索。
1 材料與方法
1.1 基因芯片數據
下載于 2020 年 1 月共享在 GEO 數據庫的 GEO 數據系列 GSE143272[5, 6]。使用了 GSE143272 中的 6 個未用(AEDs 治療)特發性癲癇患者的外周血樣本(A 組)、8 個丙戊酸鈉治療有效的特發性癲癇患者外周血樣本(B 組)、以及 10 個健康對照樣本(C 組)的基因芯片數據。GSE143272 是在 GPL10558(Illumina HumanHT-12 V4.0 expression beadchip)平臺上進行表達譜檢測。本研究中所有基因表達數據均從公共數據庫下載。本研究中使用的樣本信息見表 1。

1.2 數據預處理和差異表達基因的篩選
在 R 軟件(4.0.2 版本)中對上述基因芯片數據進行了處理。在處理數據時使用了 Biobase 程序包、GEO query 程序包及 limma 程序包[7]。將 A 組樣本與 C 組樣本進行對比(組 1)、同時 B 組樣本與 A 組樣本(組 2)進行對比篩選 DEGs。通過將每個系列矩陣的探針名轉換成基因符號、建立比較模型、貝葉斯檢驗等方法在 R 軟件中對上述數據進行對比篩選出 DEGs。用 P 值<0.05 及│fold change│ > 2 的標準定義 DEGs。
1.3 層次聚類分析
層次聚類分析是一種聚類分析,其目的是建立一個聚類層次,并能在二元樹中對相似的元素進行分組。在 R 軟件中利用 gplots 程序包和 R ColorBrewer 程序包,提取每個數據集 DEGs 的表達值進行層次聚類分析,并在熱圖中進行可視化。
1.4 差異表達基因的功能注釋及通路富集分析
將組 1、組 2 的 DEGs 取并集進行分析。利用在線工具 DIVID[8, 9](6.8 版本,網址:https://david-d.ncifcrf.gov/)對這些 DEGs 的進行了基因本體論(GO)和通路富集分析,取 P 值<0.05、基因計數>2 的結果。
1.5 構建 PPI 網絡及關鍵基因分析
通過將所有 DEGs 導入 STRING 數據庫[10] (11.0 版本,網址:http://string-db.org)并在線計算,分析這些 DEGs 編碼的蛋白質之間的相互作用網絡。隨后,在 Cytoscape 加載 STRING 網絡,在 Cytoscape 中使用 CytoHubba 應用程序[11](版本 0.1)通過 MCC 方法計算出相應權重,繪出 PPI 網絡圖,以使這些關鍵基因(hub gene)更直觀可視化。
2 結 果
2.1 差異表達基因篩選
以 P < 0.05 和│fold change│>2 作為標準,從組 1 的對比得到了 46 個 DEGs(31 個上調和 15 個下調),從組 2 的對比得到了 32 個 DEGs(23 個上調和 9 個下調)。對兩組結果取交集得到 4 個 DEGs,分別為:TREML3P(Triggering receptor expressed on myeloid cells like 3,pseudogene;髓細胞表達的觸發受體樣 3,假基因)、KCNJ15(Potassium voltage-gated channel subfamily J member 15;鉀電壓門控通道亞族 J 成員 15)、ORM1(Orosomucoid 1;血清類黏蛋白 1)、RNA28S5(RNA,28S ribosomal 5;28S 核糖體 RNA 基因 5)。其中 TREML3P、KCNJ15、ORM1 在組 1 中表達上調在組 2 中下調,RNA28S5 在組 1 中表達下調,組 2 中上調,見表 2。

2.2 差異表達基因的層次聚類分析
在獲得 DEGs 的表達值之后,對 DEGSs 進行層次聚類分析。圖 1 為組 1 的熱圖,圖 2 為組 2 的熱圖,如圖所示,兩組的 DEGs 均可以清楚地區分出實驗樣本和對照樣本。


2.3 差異表達基因的功能注釋及通路富集分析
利用在線工具 DAVID 對兩組 DEGs 的并集進行分析。得到了這些 DEGs 的生物學過程、細胞成分、分子功能、KEGG(京都基因和基因組百科全書,Kyoto Encyclopedia of Genes and Genomes)通路富集結果。在表 3 中列出了結果的前 10 位。圖 3、4、5 分別顯示了生物學過程、分子功能、KEGG 通路結果。




2.4 PPI 網絡分析和關鍵基因分析
首先,通過將所有 DEGs 導入 STRING 數據庫并在 Cytoscape 軟件里加載 STRING 網絡,74 個 DEGs 中的 57 個被篩選進入到 DEGs PPI 網絡復合體中,包含 57 個節點和 260 條聯線,其中 17 個候選 DEGs 未進入該 DEGs PPI 網絡;隨后,在 Cytoscape 中使用 Cytohubba 應用程序通過 MCC 方法得到了 55 個 DEGs 的關鍵程度并以不同色階進行展示,如圖 6 所示。圖 7 則更進一步清晰地展示了前 26 位關鍵基因。


3 討論
在本研究中,通過將 A 組樣本與 C 組樣本對比(組 1)以及 B 組樣本與 A 組樣本對比(組 2)篩選 DEGs,發現 TREML3P、KCNJ15、ORM1 在組 1 中表達上調在組 2 中下調,而 RNA28S5 在組 1 中表達下調在組 2 中上調。提示 TREML3P、KCNJ15、ORM1、RNA28S5 可能在特發性癲癇發病機制及丙戊酸鈉治療機制上有一定作用。其中,鉀電壓門控通道類基因已有研究發現與癲癇有相關[12, 13]。這 4 個基因的具體功能尚待進一步研究。關鍵基因分析篩選出 26 位關鍵基因,按權重值由大到小分別為:中性粒細胞表達的彈性蛋白酶(Elastase,neutrophil expressed,ELANE)、抵抗素(Resistin,RETN)、精氨酸酶 1(Arginase 1,ARG1)、脂質運載蛋白 2(Lipocalin 2,LCN2)、分泌性白細胞肽酶抑制劑(Secretory leukocyte peptidase inhibitor,SLPI)、結合珠蛋白(Haptoglobin,HP)、肽聚糖識別蛋白 1(Peptidoglycan recognition protein 1,PGLYRP1)、殺菌/滲透性增加蛋白(Bactericidal/permeability-increasing protein,BPI)、防御素 α4(Defensin alpha 4,DEFA4)、運鈷胺素蛋白 1(Transcobalamin 1,TCN1)、ORM1、髓過氧化物酶(Myeloperoxidase,MPO)、基質金屬肽酶 9(Matrix metallopeptidase 9,MMP9)、組織蛋白酶 G(Cathepsin G;CTSG)、C-X-C 基序趨化因子配體 8(C-X-C motif chemokine ligand 8,CXCL8)、核糖核酸酶 A 家族成員 3(Ribonuclease A family member 3,RNASE3)、核糖核酸酶 A 家族成員 2(Ribonuclease A family member 2,RNASE2)、S100 鈣結合蛋白 A12(S100 calcium binding protein A12,S100A12)、防御素 α1B(Defensin alpha 1B,DEFA1B)、防御素 α1(Defensin alpha 1,DEFA1)、防御素 α3(Defensin alpha 3,DEFA3)、癌胚抗原相關細胞粘附分子 8(Carcinoembryonic antigen related cell adhesion molecule 8,CEACAM8)、四次跨膜蛋白結構域 A3(Membrane spanning 4-domains A3,MS4A3)、前列腺素-內過氧化物合成酶 2(Prostaglandin-endoperoxide synthase 2,PTGS2)、肽酶抑制劑 3(Peptidase inhibitor 3,PI3)、C-C 基序趨化因子配體 3(C-C motif chemokine ligand 3,CCL3)。
我們對上述關鍵基因與癲癇患者的關系進行的文獻檢索,對近十來年的相關文獻進行查閱,發現我們篩選的一些關鍵 DEGs 與癲癇發生發展及VPA的治療作用有相關,部分可作為生物標志物進行預后判斷,例如:有研究發現顳葉癲癇患者外周血中 RETN 水平相對于對照組更低[14];而另一項研究則發現 RETN 在VPA治療后的第 6、12 個月呈下降[15]。研究發現,結合珠蛋白(HP)與癲癇發作相關[16, 17]。還有研究發現,毛果蕓香堿誘發的癲癇小鼠的腦 MPO 表達顯著增加[18]。此外,通過與神經炎癥的相互作用,活性氧(Reactive oxygen species,ROS)包括 MPO 在創傷后癲癇的發生中有一定作用[19]。目前有研究認為 MMP9 和 S100 是血腦屏障功能障礙的生物標志物,發現它們的濃度均在癲癇發作后明顯增加[20],S100B 和 MMP9 濃度與癲癇發作在時間和解剖學上有相關,可作為神經元損傷的生物標志物并有助于顳葉癲癇的預后判斷[21]。另外,與健康對照相比,顳葉癲癇患者血清中組織蛋白酶 B(cathepsin B)有所增高[22]。在一些不同病因的癲癇中,CXCL8 呈現升高[23],而其他一些趨化因子 CCL3 等也參與了癲癇的發生發展[24]。相比VPA治療無效者,VPA治療有效者外周血中 PTGS2 表現出了下調[6]。現有研究數據支持我們的研究結果,對我們的這項生信分析結果做了部分驗證,證明了我們篩選關鍵DEGs的方法是有效的。同時,那些未被深入研究的關鍵分子在癲癇中的作用及其能否作為癲癇患者診治與預后判斷的生物標志物,均值得我們進一步研究。
而通過 GO 分析顯示:DEGs 生物學過程集中在免疫反應、防御反應、炎癥反應、趨化作用等,分子功能集中在過氧化物酶活性、趨化因子活性等,而 KEGG 通路富集分析 DEGs 主要參與細胞因子-細胞因子受體相互作用、Toll 樣受體信號通路、趨化因子信號通路等信號通路。這提示炎癥及免疫機制也與特發性癲癇密切相關。現有研究發現免疫反應[25]、趨化因子[24]、細胞因子及炎癥信號通路[23, 26]等均參與了癲癇的發生發展。而我們篩選的部分重要的關鍵 DEGs 可能通過多種信號通路及復雜的機制,從而參與特發性癲癇的發生發展。
綜上,現有研究報道支持我們該項生物信息學分析結果。同時,在本研究中仍有很多分子及其在癲癇發生發展中的作用尚待進一步的驗證和研究,以進一步闡明特發性癲癇的發病機制,為尋找臨床上可用的癲癇早期診斷、預后判斷和治療的靶點提供一定的線索。
癲癇(Epilepsy)是一種常見的嚴重慢性腦部疾病,全球共約 7 000 萬人受到該病的影響[1],流行病學調查顯示其發病率為 0.5%~0.7%,我國約有 650~910 萬癲癇患者,該病可見于各個年齡組,青少年和老年是癲癇發病的兩個高峰年齡段[2]。而特發性癲癇是由基因突變和某些先天因素所致,有明顯遺傳傾向的具體病因尚不清楚的癲癇[2]。人類遺傳物質的變異類型主要有基因突變(gene mutation)、染色體結構或數目變異,這些遺傳變異均可導致癲癇的發生[3]。丙戊酸鈉(VPA)作為一種廣譜抗癲癇藥物(AEDs),可用于多種發作類型的癲癇[2]。在治療全面性癲癇或未分類癲癇時,VPA是最為有效的 AEDs[4]。目前癲癇尚無好的可用于其早期診斷與預后判斷的分子標記物。大多數患者需要長期使用 AEDs 治療。因為疾病的反復發作性,癲癇患者多缺乏自信、情緒消極,嚴重影響學習、工作和生活。近年來,利用高通量測序等方法獲得了大量生物數據,而分析不同情況下的差異表達基因(Differentially expressed genes,DEGs),能夠很好篩選新的有效的分子靶點。目前數據庫中有使用基因芯片分析癲癇患者外周血基因改變的研究結果,我們通過未用AEDs 治療的特發性癲癇外周血樣本與健康對照樣本比對,及VPA治療有效的特發性癲癇外周血樣本與未用藥患者比對,綜合相關結果,試圖尋找癲癇患者外周血關鍵基因及通路改變,以進一步闡明特發性癲癇患者的發病及VPA治療相關機制,為臨床上可用的癲癇早期診斷、預后判斷和治療的靶點提供線索。
1 材料與方法
1.1 基因芯片數據
下載于 2020 年 1 月共享在 GEO 數據庫的 GEO 數據系列 GSE143272[5, 6]。使用了 GSE143272 中的 6 個未用(AEDs 治療)特發性癲癇患者的外周血樣本(A 組)、8 個丙戊酸鈉治療有效的特發性癲癇患者外周血樣本(B 組)、以及 10 個健康對照樣本(C 組)的基因芯片數據。GSE143272 是在 GPL10558(Illumina HumanHT-12 V4.0 expression beadchip)平臺上進行表達譜檢測。本研究中所有基因表達數據均從公共數據庫下載。本研究中使用的樣本信息見表 1。

1.2 數據預處理和差異表達基因的篩選
在 R 軟件(4.0.2 版本)中對上述基因芯片數據進行了處理。在處理數據時使用了 Biobase 程序包、GEO query 程序包及 limma 程序包[7]。將 A 組樣本與 C 組樣本進行對比(組 1)、同時 B 組樣本與 A 組樣本(組 2)進行對比篩選 DEGs。通過將每個系列矩陣的探針名轉換成基因符號、建立比較模型、貝葉斯檢驗等方法在 R 軟件中對上述數據進行對比篩選出 DEGs。用 P 值<0.05 及│fold change│ > 2 的標準定義 DEGs。
1.3 層次聚類分析
層次聚類分析是一種聚類分析,其目的是建立一個聚類層次,并能在二元樹中對相似的元素進行分組。在 R 軟件中利用 gplots 程序包和 R ColorBrewer 程序包,提取每個數據集 DEGs 的表達值進行層次聚類分析,并在熱圖中進行可視化。
1.4 差異表達基因的功能注釋及通路富集分析
將組 1、組 2 的 DEGs 取并集進行分析。利用在線工具 DIVID[8, 9](6.8 版本,網址:https://david-d.ncifcrf.gov/)對這些 DEGs 的進行了基因本體論(GO)和通路富集分析,取 P 值<0.05、基因計數>2 的結果。
1.5 構建 PPI 網絡及關鍵基因分析
通過將所有 DEGs 導入 STRING 數據庫[10] (11.0 版本,網址:http://string-db.org)并在線計算,分析這些 DEGs 編碼的蛋白質之間的相互作用網絡。隨后,在 Cytoscape 加載 STRING 網絡,在 Cytoscape 中使用 CytoHubba 應用程序[11](版本 0.1)通過 MCC 方法計算出相應權重,繪出 PPI 網絡圖,以使這些關鍵基因(hub gene)更直觀可視化。
2 結 果
2.1 差異表達基因篩選
以 P < 0.05 和│fold change│>2 作為標準,從組 1 的對比得到了 46 個 DEGs(31 個上調和 15 個下調),從組 2 的對比得到了 32 個 DEGs(23 個上調和 9 個下調)。對兩組結果取交集得到 4 個 DEGs,分別為:TREML3P(Triggering receptor expressed on myeloid cells like 3,pseudogene;髓細胞表達的觸發受體樣 3,假基因)、KCNJ15(Potassium voltage-gated channel subfamily J member 15;鉀電壓門控通道亞族 J 成員 15)、ORM1(Orosomucoid 1;血清類黏蛋白 1)、RNA28S5(RNA,28S ribosomal 5;28S 核糖體 RNA 基因 5)。其中 TREML3P、KCNJ15、ORM1 在組 1 中表達上調在組 2 中下調,RNA28S5 在組 1 中表達下調,組 2 中上調,見表 2。

2.2 差異表達基因的層次聚類分析
在獲得 DEGs 的表達值之后,對 DEGSs 進行層次聚類分析。圖 1 為組 1 的熱圖,圖 2 為組 2 的熱圖,如圖所示,兩組的 DEGs 均可以清楚地區分出實驗樣本和對照樣本。


2.3 差異表達基因的功能注釋及通路富集分析
利用在線工具 DAVID 對兩組 DEGs 的并集進行分析。得到了這些 DEGs 的生物學過程、細胞成分、分子功能、KEGG(京都基因和基因組百科全書,Kyoto Encyclopedia of Genes and Genomes)通路富集結果。在表 3 中列出了結果的前 10 位。圖 3、4、5 分別顯示了生物學過程、分子功能、KEGG 通路結果。




2.4 PPI 網絡分析和關鍵基因分析
首先,通過將所有 DEGs 導入 STRING 數據庫并在 Cytoscape 軟件里加載 STRING 網絡,74 個 DEGs 中的 57 個被篩選進入到 DEGs PPI 網絡復合體中,包含 57 個節點和 260 條聯線,其中 17 個候選 DEGs 未進入該 DEGs PPI 網絡;隨后,在 Cytoscape 中使用 Cytohubba 應用程序通過 MCC 方法得到了 55 個 DEGs 的關鍵程度并以不同色階進行展示,如圖 6 所示。圖 7 則更進一步清晰地展示了前 26 位關鍵基因。


3 討論
在本研究中,通過將 A 組樣本與 C 組樣本對比(組 1)以及 B 組樣本與 A 組樣本對比(組 2)篩選 DEGs,發現 TREML3P、KCNJ15、ORM1 在組 1 中表達上調在組 2 中下調,而 RNA28S5 在組 1 中表達下調在組 2 中上調。提示 TREML3P、KCNJ15、ORM1、RNA28S5 可能在特發性癲癇發病機制及丙戊酸鈉治療機制上有一定作用。其中,鉀電壓門控通道類基因已有研究發現與癲癇有相關[12, 13]。這 4 個基因的具體功能尚待進一步研究。關鍵基因分析篩選出 26 位關鍵基因,按權重值由大到小分別為:中性粒細胞表達的彈性蛋白酶(Elastase,neutrophil expressed,ELANE)、抵抗素(Resistin,RETN)、精氨酸酶 1(Arginase 1,ARG1)、脂質運載蛋白 2(Lipocalin 2,LCN2)、分泌性白細胞肽酶抑制劑(Secretory leukocyte peptidase inhibitor,SLPI)、結合珠蛋白(Haptoglobin,HP)、肽聚糖識別蛋白 1(Peptidoglycan recognition protein 1,PGLYRP1)、殺菌/滲透性增加蛋白(Bactericidal/permeability-increasing protein,BPI)、防御素 α4(Defensin alpha 4,DEFA4)、運鈷胺素蛋白 1(Transcobalamin 1,TCN1)、ORM1、髓過氧化物酶(Myeloperoxidase,MPO)、基質金屬肽酶 9(Matrix metallopeptidase 9,MMP9)、組織蛋白酶 G(Cathepsin G;CTSG)、C-X-C 基序趨化因子配體 8(C-X-C motif chemokine ligand 8,CXCL8)、核糖核酸酶 A 家族成員 3(Ribonuclease A family member 3,RNASE3)、核糖核酸酶 A 家族成員 2(Ribonuclease A family member 2,RNASE2)、S100 鈣結合蛋白 A12(S100 calcium binding protein A12,S100A12)、防御素 α1B(Defensin alpha 1B,DEFA1B)、防御素 α1(Defensin alpha 1,DEFA1)、防御素 α3(Defensin alpha 3,DEFA3)、癌胚抗原相關細胞粘附分子 8(Carcinoembryonic antigen related cell adhesion molecule 8,CEACAM8)、四次跨膜蛋白結構域 A3(Membrane spanning 4-domains A3,MS4A3)、前列腺素-內過氧化物合成酶 2(Prostaglandin-endoperoxide synthase 2,PTGS2)、肽酶抑制劑 3(Peptidase inhibitor 3,PI3)、C-C 基序趨化因子配體 3(C-C motif chemokine ligand 3,CCL3)。
我們對上述關鍵基因與癲癇患者的關系進行的文獻檢索,對近十來年的相關文獻進行查閱,發現我們篩選的一些關鍵 DEGs 與癲癇發生發展及VPA的治療作用有相關,部分可作為生物標志物進行預后判斷,例如:有研究發現顳葉癲癇患者外周血中 RETN 水平相對于對照組更低[14];而另一項研究則發現 RETN 在VPA治療后的第 6、12 個月呈下降[15]。研究發現,結合珠蛋白(HP)與癲癇發作相關[16, 17]。還有研究發現,毛果蕓香堿誘發的癲癇小鼠的腦 MPO 表達顯著增加[18]。此外,通過與神經炎癥的相互作用,活性氧(Reactive oxygen species,ROS)包括 MPO 在創傷后癲癇的發生中有一定作用[19]。目前有研究認為 MMP9 和 S100 是血腦屏障功能障礙的生物標志物,發現它們的濃度均在癲癇發作后明顯增加[20],S100B 和 MMP9 濃度與癲癇發作在時間和解剖學上有相關,可作為神經元損傷的生物標志物并有助于顳葉癲癇的預后判斷[21]。另外,與健康對照相比,顳葉癲癇患者血清中組織蛋白酶 B(cathepsin B)有所增高[22]。在一些不同病因的癲癇中,CXCL8 呈現升高[23],而其他一些趨化因子 CCL3 等也參與了癲癇的發生發展[24]。相比VPA治療無效者,VPA治療有效者外周血中 PTGS2 表現出了下調[6]。現有研究數據支持我們的研究結果,對我們的這項生信分析結果做了部分驗證,證明了我們篩選關鍵DEGs的方法是有效的。同時,那些未被深入研究的關鍵分子在癲癇中的作用及其能否作為癲癇患者診治與預后判斷的生物標志物,均值得我們進一步研究。
而通過 GO 分析顯示:DEGs 生物學過程集中在免疫反應、防御反應、炎癥反應、趨化作用等,分子功能集中在過氧化物酶活性、趨化因子活性等,而 KEGG 通路富集分析 DEGs 主要參與細胞因子-細胞因子受體相互作用、Toll 樣受體信號通路、趨化因子信號通路等信號通路。這提示炎癥及免疫機制也與特發性癲癇密切相關。現有研究發現免疫反應[25]、趨化因子[24]、細胞因子及炎癥信號通路[23, 26]等均參與了癲癇的發生發展。而我們篩選的部分重要的關鍵 DEGs 可能通過多種信號通路及復雜的機制,從而參與特發性癲癇的發生發展。
綜上,現有研究報道支持我們該項生物信息學分析結果。同時,在本研究中仍有很多分子及其在癲癇發生發展中的作用尚待進一步的驗證和研究,以進一步闡明特發性癲癇的發病機制,為尋找臨床上可用的癲癇早期診斷、預后判斷和治療的靶點提供一定的線索。