膠質母細胞瘤中異常甲基化調控的蛋白編碼基因識別研究

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

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

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

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

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

1 數據來源及方法

1.1 數據的來源，整合和標準化

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

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

1.2 PCGs甲基化譜的構建

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

1.3 不同甲基化模式PCGs功能的預測

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

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

1.4 生存分析

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

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

2 結果

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

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

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

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

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

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

對不同甲基化模式下的PCGs，采用富集分析方法進行功能和通路的顯著性分析。結果顯示，高甲基化的PCGs主要富集在神經元系統發展、細胞-細胞信號傳導等生物過程;低甲基化的PCGs富集在許多與腫瘤進程相關的生物過程，比如細胞黏附、細胞遷移、免疫相關的細胞增生和血管生成(圖1A、1B)。因此，低甲基化的PCGs可能是影響GBM進程的一個重要模塊。對于KEGG通路富集分析，高甲基化的PCGs主要富集在GABAergic突觸、胰島素分泌等通路中，而低甲基化的PCGs主要富集在金黃色葡萄球菌感染、補體和凝血級聯等通路中(圖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中預后良好的PCGs的鑒別

為了評估這些表達受對應啟動子異常甲基化影響的PCGs是否可以作為GBM的預后因子，本研究將表達譜數據與臨床信息相結合，進行基因生存分析(圖2)，基于生存分析的結果，獲得了51個與GBM整體存活顯著相關的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 預后良好的PCGsTab.1 PCGs with good prognosis

3 討論

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

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

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

綜上，本研究通過重新注釋DNA甲基化陣列，系統識別GBM中潛在的表達受異常甲基化調控的PCGs，加深對GBM中PCGs甲基化調控模式的理解，并對識別GBM風險標志物和潛在的治療靶點提出了新的認識。