嵇江淮 趙瀟瀟李乾鵬 安 奕 趙 磊 李冬果*

(1.首都醫(yī)科大學(xué)生物醫(yī)學(xué)工程學(xué)院生物醫(yī)學(xué)信息學(xué)系，北京 100069；2.首都醫(yī)科大學(xué)宣武醫(yī)院麻醉手術(shù)科，北京 100053；3. 國家老年疾病臨床研究中心，北京 100053)

膠質(zhì)母細胞瘤(glioblastoma, GBM)是人類最常見且致死率極高的一種中樞神經(jīng)系統(tǒng)腫瘤，發(fā)生率約占膠質(zhì)瘤的69%[1]。這種腦瘤具有高浸潤性，預(yù)后差，患者的中位生存期大約只有一年[2]。近年來，隨著大量的蛋白編碼基因(protein-coding genes, PCGs)被發(fā)現(xiàn)，編碼基因表達的失調(diào)被證明與疾病的進程有密切關(guān)系[3-4]。有研究[5]表明GBM中EGFR基因表達顯著升高，促進GBM患者細胞的凋亡。Meng等[6]發(fā)現(xiàn)TCTN1基因的過表達促進GBM細胞的增生，并且TCTN1的過表達可以作為預(yù)測GBM患者的獨立預(yù)后因素。這些研究表明PCGs在GBM進程中發(fā)揮非常重要的作用，但是這些研究主要關(guān)注PCGs的表達模式。目前，GBM中絕大部分PCGs的調(diào)控機制(特別是DNA甲基化調(diào)控機制)尚不清楚。

DNA甲基化是表觀遺傳調(diào)控的一種重要模式。在基因啟動子處的DNA甲基化對基因的表達具有重要的作用，并且參與到許多疾病的進程中[7]。一般來說，基因啟動子處的高甲基化往往抑制轉(zhuǎn)錄因子的組合并且下調(diào)甚至沉默癌癥抑制基因，基因啟動子的低甲基化往往激活致癌基因[8]。研究[9]表明GBM中MGMT基因啟動子甲基化和MGMT的表達具有很強的相關(guān)性。然而，這些研究僅僅分析一些特殊的PCGs的DNA甲基化模式，并沒有系統(tǒng)分析腫瘤中PCGs的DNA甲基化的模式以及全面評估腫瘤中DNA甲基化對基因表達的調(diào)控關(guān)系。

隨著高通量測序技術(shù)的發(fā)展，RNA-seq和Infnium 450k等眾多數(shù)據(jù)已經(jīng)應(yīng)用于癌癥分析。本研究系統(tǒng)整合多組學(xué)數(shù)據(jù)，通過重新注釋DNA甲基化陣列，系統(tǒng)分析GBM進程中PCGs的DNA甲基化調(diào)控模式。本研究將會幫助理解GBM中DNA甲基化的調(diào)控機制，為識別GBM生物標(biāo)志物和潛在的治療靶點提供新的見解。

1 數(shù)據(jù)來源及方法

1.1 數(shù)據(jù)的來源，整合和標(biāo)準(zhǔn)化

GBM的DNA甲基化數(shù)據(jù)來自Infinium HM450k平臺[10]，本研究從TCGA中下載GBM腫瘤樣本的HM450k數(shù)據(jù)，從GEO(GSE41826[11-12],其中包含58個正常膠質(zhì)細胞樣本)中下載正常樣本的HM450k數(shù)據(jù)。正常數(shù)據(jù)是通過同一平臺用同樣方法獲得的。GBM表達譜數(shù)據(jù)來自于TCGA數(shù)據(jù)庫，包括5個正常樣本數(shù)據(jù)和152個臨床樣本信息。為了保證數(shù)據(jù)的質(zhì)量，本研究挑選50個同時具有DNA甲基化數(shù)據(jù)和表達譜數(shù)據(jù)的樣本用作進一步的分析。人類全基因組的注釋數(shù)據(jù)(V19)收集于GENCODE數(shù)據(jù)庫[13]。

對下載的RNA-seq(counts)數(shù)據(jù)進行l(wèi)og2轉(zhuǎn)換并使用R包“edgeR”[14]對數(shù)據(jù)進行標(biāo)準(zhǔn)化，最終獲得PCGs的表達譜數(shù)據(jù)。由于HM450k數(shù)據(jù)存在一定的缺失值，所以，在構(gòu)建甲基化譜之前需要對含有缺失的數(shù)據(jù)進行處理。K最近鄰填充算法 (K-nearest neighbor, KNN) 是用數(shù)據(jù)集中缺失數(shù)據(jù)的k個最近鄰來估計缺失值[15]。本研究計算在所有腫瘤樣本中具有缺失值的探針數(shù)量，并用“DMwR”包中的knnImputation函數(shù)補全剩余的缺失值，以評估甲基化探針的質(zhì)量。最終有89 512個探針被移除，獲得392 867個甲基化位點探針。

1.2 PCGs甲基化譜的構(gòu)建

本研究用映射到PCGs啟動子區(qū)域探針的甲基化水平來構(gòu)建PCGs的甲基化譜。采用Zhi等[16]的重注釋方法將392 867個探針映射到PCGs的啟動子區(qū)域(轉(zhuǎn)錄開始位點到上游10 kb區(qū)域內(nèi))，使用最接近每個轉(zhuǎn)錄開始位點的探針來確定PCGs啟動子的DNA甲基化水平[17]，從而構(gòu)建GBM相關(guān)的PCGs甲基化譜。

1.3 不同甲基化模式PCGs功能的預(yù)測

本研究使用基于線性模型設(shè)計的R包“l(fā)imma”[18]計算腫瘤和正常樣本之間的差異甲基化和差異表達。通過Benjamini-Hochberg方法校正P值。為了全面分析表達受對應(yīng)啟動子異常甲基化調(diào)控的PCGs，分別將高甲基化的PCGs和低表達的PCGs及低甲基化的PCGs和高表達的PCGs取交集，得到兩種具有不同生物學(xué)意義的情況：高甲基化且低表達的PCGs及低甲基化且高表達的PCGs。

為了預(yù)測不同甲基化模式下PCGs的功能，采用富集分析方法，對挑選出的PCGs進行功能和通路的顯著性分析，使用“clusterProfiler”[19]包來預(yù)測不同甲基化模式下PCGs的功能。通過Benjamini-Hochberg方法校正P值，如果矯正后的P值≤0.05，該GO項和富集通路就認(rèn)為是顯著的。

1.4 生存分析

為了鑒別挑選出的PCGs是否具有良好的預(yù)后效果，本研究基于152個GBM患者的表達譜信息，依據(jù)中值將患者分為兩組。Kaplan-Meier生存分析和log-rank檢測被用來評估兩組患者的生存差異。

所有的研究均使用R3.5.1完成。

2 結(jié)果

2.1 重注釋Infinium HM 450K陣列并構(gòu)建PCGs的DNA甲基化譜

為了描繪PCGs的DNA甲基化調(diào)控模式，本研究將甲基化數(shù)據(jù)重新注釋到人類PCGs相關(guān)的啟動子區(qū)域，共有125 442個探針落在14 684個PCGs啟動子區(qū)域，通常每一個PCG都有幾個探針落在其啟動子區(qū)域，在此僅保留最接近每個轉(zhuǎn)錄開始位點的探針來確定PCGs啟動子的DNA甲基化狀態(tài)。

2.2 整合差異甲基化和差異表達的PCGs

為了能夠有效地識別GBM相關(guān)的風(fēng)險標(biāo)志物，研究差異甲基化可能的生物學(xué)意義，基于鑒別出的3 561個差異甲基化的PCGs，其中高甲基化的PCGs有873個(24.5%)，低甲基化的PCGs有2 688個(75.5%)。顯然低甲基化的PCGs的數(shù)量遠遠多于高甲基化的PCGs的數(shù)量。這些結(jié)果顯示，在GBM的發(fā)生發(fā)展中PCGs呈現(xiàn)為更多的低甲基化模式，并且這種全局的低甲基化可能導(dǎo)致致癌基因的激活和影響基因組的穩(wěn)定性。進一步地，本研究鑒別出 6 586 個差異表達的PCGs，其中有2 788個PCGs表達上調(diào)，3 798個PCGs表達下調(diào)。

為了分析異常甲基化對PCGs表達的調(diào)控機制，針對兩種不同生物學(xué)意義的情況，識別出240個高甲基化且低表達的PCGs及390個低甲基化且高表達的PCGs。結(jié)果表明，PCGs在腫瘤樣本和正常樣本中表現(xiàn)出差異甲基化和差異表達兩種模式。

2.3 不同甲基化模式下PCGs的富集分析

對不同甲基化模式下的PCGs，采用富集分析方法進行功能和通路的顯著性分析。結(jié)果顯示，高甲基化的PCGs主要富集在神經(jīng)元系統(tǒng)發(fā)展、細胞-細胞信號傳導(dǎo)等生物過程;低甲基化的PCGs富集在許多與腫瘤進程相關(guān)的生物過程，比如細胞黏附、細胞遷移、免疫相關(guān)的細胞增生和血管生成(圖1A、1B)。因此，低甲基化的PCGs可能是影響GBM進程的一個重要模塊。對于KEGG通路富集分析，高甲基化的PCGs主要富集在GABAergic突觸、胰島素分泌等通路中，而低甲基化的PCGs主要富集在金黃色葡萄球菌感染、補體和凝血級聯(lián)等通路中(圖1C、1D)。

圖1 不同甲基化模式下PCGs的富集分析Fig.1 Enrichment analysis of PCGs under different methylation patternsA and B are the GO function analysis of hypermethylated-low expressed PCGs and hypomethylated-highly expressed PCGs,respectively. C and D are KEGG pathway analysis of hypermethylated-low expressed PCGs and hypomethylated-highly expressed PCGs,respectively.The depth of color represents the number of PCGs enriched in the GO item or pathway; PCGs: protein-coding gens.

2.4 GBM中預(yù)后良好的PCGs的鑒別

為了評估這些表達受對應(yīng)啟動子異常甲基化影響的PCGs是否可以作為GBM的預(yù)后因子，本研究將表達譜數(shù)據(jù)與臨床信息相結(jié)合，進行基因生存分析(圖2)，基于生存分析的結(jié)果，獲得了51個與GBM整體存活顯著相關(guān)的PCGs(P<0.05)，如表1所示。

圖2 SOCS1和AEBP1基因生存分析Fig.2 Survival analysis of SOCS1 and AEBP1 genesA and B represent Kaplan-Meier curves of SOCS1 and AEBP1 genes in GBM, respectively.The red and blue lines at the bottom represent the number of patients who survived.

表1 預(yù)后良好的PCGsTab.1 PCGs with good prognosis

3 討論

DNA甲基化的變化可能破壞癌癥中特定啟動子的調(diào)控[20]。近年來，隨著高通量測序技術(shù)的發(fā)展，表觀遺傳調(diào)控已經(jīng)成為生物醫(yī)學(xué)研究的一大熱點。一項研究[21]表明GBM中B3GNT5、FABP7等13個基因啟動子的甲基化和表達水平之間呈現(xiàn)負(fù)相關(guān)。

本研究對重注釋得到的PCGs甲基化譜進行差異甲基化分析，結(jié)果顯示,GBM中多呈現(xiàn)一種低甲基化的模式。為了驗證這個結(jié)果，本研究重新注釋了GBM中l(wèi)ncRNA的甲基化譜，結(jié)果發(fā)現(xiàn)在鑒別出的5 567個異常甲基化的lncRNAs中高甲基化的有1 214個(21.8%)，低甲基化的有4 353個(78.2%)。為了評估是否不同甲基化模式的PCGs對應(yīng)著不同的生物學(xué)功能，對這些基因進行功能富集分析。結(jié)果顯示低甲基化的PCGs參與了許多與腫瘤發(fā)生及進程相關(guān)的生物過程。由此推斷這些低甲基化的PCGs更有可能成為GBM診斷和治療的生物標(biāo)志物。有研究[22-23]表明AEBP1和SOCS1在GBM中過表達,沉默他們的表達可以抑制GBM細胞的增生。在本研究中，AEBP1和SOCS1均屬于低甲基化且高表達的基因，進一步地，AEBP1和SOCS1基因表達水平較低的GBM患者具有更加良好的預(yù)后。結(jié)果表明AEBP1和SOCS1在GBM中表達的上調(diào)可能是受它們基因啟動子的低甲基化調(diào)控，而抑制它們的表達可以提高GBM患者的生存概率。

當(dāng)然，由于GBM不同于一般的腫瘤，樣本的獲取是一個重大的難題。本研究需要進一步搜集大規(guī)模樣本進行重復(fù)性計算或者動物實驗來驗證挖掘出的表達受異常甲基化調(diào)控的PCGs以及潛在的治療靶點。

綜上，本研究通過重新注釋DNA甲基化陣列，系統(tǒng)識別GBM中潛在的表達受異常甲基化調(diào)控的PCGs，加深對GBM中PCGs甲基化調(diào)控模式的理解，并對識別GBM風(fēng)險標(biāo)志物和潛在的治療靶點提出了新的認(rèn)識。