張淑梅，張 彬，劉軍厚，劉洪波，蘇建忠，王 芳，張 巖

(哈爾濱醫(yī)科大學生物信息科學與技術(shù)學院，黑龍江哈爾濱150081)

目前，癌癥是嚴重威脅人類健康的三大殺手之一，對于這種嚴重危害人類健康的頑疾，現(xiàn)在的醫(yī)學界并不十分清楚它的發(fā)病機制。同時，人們對基因的本質(zhì)也漸漸有了更深入地認識。很長一段時間里，人們認為癌癥的形成只與基因突變有關(guān)［1－4］。但是，越來越多的證據(jù)表明，表觀遺傳修飾對癌癥的發(fā)生也起著十分重要的作用?；蛐蛄胁蛔?，而基因的表型發(fā)生了可遺傳的變化，稱為表觀遺傳［5］。這是由表觀遺傳修飾造成的。DNA甲基化是一種重要的表觀遺傳修飾，它制約著基因的表達。

據(jù)報道，人類的基因只是果蠅的2倍多。很難想象DNA的遺傳信息可以調(diào)控人類這樣復雜的生命體發(fā)育和生存的全過程［1－4］。維持細胞的功能，決定哪些基因表達、哪些基因不表達，是非常重要的，幾個基因的錯誤表達便會誘發(fā)正常細胞發(fā)生癌變［2］。

目前普遍認為，DNA甲基化與癌癥的發(fā)生有密切關(guān)系［6］。癌癥的甲基化異常表現(xiàn)為總體的甲基化水平降低與啟動子區(qū)域的甲基化水平升高［7］。例如，抑癌基因與修復基因的高甲基化會導致它們的失活，造成腫瘤抑制喪失與基因損傷增加。

由于涉及基因的“開”與“關(guān)”，DNA甲基化對腫瘤的產(chǎn)生起著重大的作用。同時研究表明，某些基因的異常甲基化與多種癌癥的產(chǎn)生有著顯著的關(guān)聯(lián)［8］。例如，基因P15的甲基化會使基因沉默，并使細胞過度激活與增殖，而這與白血病、淋巴瘤、鱗狀細胞癌、肺癌的發(fā)生都有重要的聯(lián)系［9］。是否存在一組甲基化異常的基因，與多種癌癥的發(fā)生有著重要的關(guān)聯(lián)以及這些基因在不同的癌癥中是否起著不同的作用，成為本文關(guān)心的問題。通過研究這些問題，會為癌癥的預測提供必要的方法，同時也增強了人們對癌癥與DNA甲基化關(guān)系更進一步的認識。

表觀遺傳標記可以在被割除的腫瘤和體液中探測到。例如，超甲基化的癌癥基因可以在尿斑中探測到，這在膀胱癌的檢測中很有意義［10］。DNA甲基化的生物標記物在疾病診斷和預后的領(lǐng)域正在興起，并且需要在臨床實踐中廣為應用和擴展。

本課題首先通過對不同癌癥DNA甲基化數(shù)據(jù)進行預處理，利用權(quán)重基因共表達網(wǎng)絡分析方法(WGCNA)篩選出甲基化基因模塊，并分析模塊向量基因，利用 DAVID(The Database for Annotation，Visualization and Integrated Discovery)進行功能注釋，然后對基因模塊進行功能分析，得到DNA甲基化與腫瘤間的關(guān)系。本課題有助于發(fā)現(xiàn)癌癥中DNA甲基化的生物標記物，為腫瘤的診斷及治療提供可能的靶點。

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

1.1 數(shù)據(jù)獲取

本課題所用的癌癥甲基化數(shù)據(jù)來自于GEO數(shù)據(jù)庫，包括乳腺導管癌甲基化數(shù)據(jù)(編號為GSE14865，平臺為 GPL4126，6 個樣本)［11］、胃癌甲基化數(shù)據(jù)(編號為GSE15291，平臺為GPL4126，7個樣本)、前列腺癌甲基化數(shù)據(jù)(編號為GSE15298，平臺為 GPL4126，20 個樣本)［12］、白血病甲基化數(shù)據(jù)(編號為GSE18400，平臺為 GPL4126，60個樣本，樣本為嬰兒期白血病數(shù)據(jù)和1個對照組)、食管鱗狀細胞瘤甲基化數(shù)據(jù)(編號為 GSE21238，平臺為GPL4126，6個樣本，其中包括有淋巴結(jié)轉(zhuǎn)移食道癌和無淋巴結(jié)轉(zhuǎn)移食道癌樣本以及轉(zhuǎn)移淋巴結(jié)細胞樣本)、肺鱗狀細胞瘤甲基化數(shù)據(jù)(編號為GSE9622，平臺為GPL4126，5個樣本)。

首先進行數(shù)據(jù)的預處理和標準化，標準化的原則是對同一基因的不同探針對應的數(shù)值取平均值，并且只選擇對應于啟動子的探針。最后獲得包含4 029個基因的甲基化數(shù)據(jù)。

1.2 權(quán)重基因共表達網(wǎng)絡分析(WGCNA)的簡介

網(wǎng)絡分析在生物信息學中得到越來越多的應用。WGCNA(Weight Gene Co-express Network Analysis)是一種描述各個樣本的基因芯片相關(guān)的系統(tǒng)生物學方法。這種方法可以找到高相關(guān)的基因模塊，可以使用模塊特征基因(eigengene)或hub節(jié)點間的基因彼此間和外部采樣特征來聚類［13］。相關(guān)網(wǎng)絡促進了基于基因篩選的方法的發(fā)展，可以用于識別候選生物標記物或治療靶點。

1.3 甲基化基因模塊篩選的原理

本文通過構(gòu)建權(quán)重基因共表達網(wǎng)絡來識別癌癥中甲基化基因模塊。

為了便于把顯著差異的甲基化基因分類成模塊，鄰接矩陣被轉(zhuǎn)換成拓撲重疊矩陣。拓撲重疊矩陣不僅可以捕捉到xi，xj的直接互作，也可以捕捉到間接互作。這樣，定義了一個相似性測度:

其中，ki=代表點的連通性。1－TOMij是層次聚類的距離矩陣。

1.4 基因模塊的功能分析

通過WGCNA篩選出甲基化基因模塊并量化模塊與表型的關(guān)系。分析與癌癥表型顯著相關(guān)的基因模塊。挖掘出基因模塊的向量基因，并利用DAVID生物信息學分析工具對基因模塊進行GO功能注釋與KEGG通路富集研究。

2 結(jié)果

2.1 選擇合適的閾值:網(wǎng)絡拓撲結(jié)構(gòu)分析

構(gòu)建一個權(quán)重基因網(wǎng)絡，選擇一個合適的鄰接矩陣的閾值β，得到的閾值滿足網(wǎng)絡接近無于尺度的標準。通過WGCNA，選擇一組候選的閾值，并返回被檢測的網(wǎng)絡參數(shù)(見圖1)。從圖中可看出閾值選擇為5最合適，它既保證了網(wǎng)絡接近于無尺度網(wǎng)絡(模型指數(shù)大于0.9，完美無尺度網(wǎng)絡的模型適應指數(shù)是1)，同時也是使曲線趨于平滑的最小閾值，并且它也使得網(wǎng)絡的平均鏈接程度不會太小，這有利于網(wǎng)絡包含足夠的信息(例如，挖掘模塊)。

圖1 閾值分析Fig.1 Threshold value analysis

2.2 使用TOM(拓撲重疊矩陣)聚類

為了減少噪聲和偽關(guān)聯(lián)的影響，將鄰接矩陣轉(zhuǎn)換為拓撲重疊矩陣(TOM)。通過TOM，利用層次聚類產(chǎn)生一個基因的層次聚類樹(見圖2)。

在層次聚類樹中，通過各個分支的識別(即從樹圖“剪枝”)得到模塊。使用Dynamic Tree Cut的方法［14］，期望獲得較大較少的模塊，所以設(shè)定參數(shù)最小模塊大小(minModuleSize)為50，這樣從樹圖中剪枝共得到10個模塊，標簽為1至10，模塊大小依次遞減，從806至65個基因。模塊0保存著所有模塊外的基因。

樹圖中不同的深淺區(qū)域代表了不同的模塊。找到匹配的模塊，并返回各基因模塊的寬度(見表1)。

注:樹圖中每個葉節(jié)點代表一個基因，其中密集連接的分支代表了甲基化數(shù)值接近的基因;圖中不同的顏色代表不同的模塊。Notes:In the tree diagram，every leaf node represents a gene，and branches densely connected represent the genes which have the similar methylation values.The different colors represent different modules.

表1 各基因模塊中的向量基因Table1 Vector genes in the modules

2.3 量化模塊與表型的關(guān)聯(lián)

分析模塊與模型的顯著關(guān)聯(lián)。由于已有每個模塊的特征基因(eigengene)，使特征基因(eigengene)與表型相關(guān)聯(lián)，并找到最大相關(guān)性。由于已有模塊與表型，可以可視化這種關(guān)聯(lián)，用顏色標注相關(guān)性。

圖3中可以清晰地看到模塊與癌癥表型的相關(guān)性。模塊0的基因是樹圖中剔除的基因，從圖3中也可看出它與各癌癥表型的相關(guān)性較差，因此不予考慮。

可以看到，胃癌(gastric cancer)與模塊1、4、10有較強的相關(guān)性(p≤0.05);前列腺癌(prostate cancer)與模塊2有較強的相關(guān)性(p≤0.05)。說明這兩種癌癥在上述模塊中甲基化程度較高。

同時，可以看到各種嬰兒期白血病(ALL)，如MLL－AF4白血病、MLL－ENL白血病、未擴散白血病(Untranslocated infant ALL)在模塊7、9都有較強的相關(guān)性(p≤0.05)，而正常人體細胞在這兩個模塊中p值都大于0.05，沒有顯著的相關(guān)性。由于選用的基因來自于啟動子區(qū)域，可以得出結(jié)論:在上述基因模塊中，白血病對應基因的甲基化程度要比正常細胞的甲基化程度高。

圖3 模塊與癌癥表型關(guān)系圖Fig.3 Module-trait relationship

食道鱗狀細胞瘤(Esophageal Squamous Cell Carcinomas)在模塊1、4、10都有顯著的相關(guān)性。同時還發(fā)現(xiàn)，有淋巴結(jié)轉(zhuǎn)移食道鱗狀細胞瘤(ESCC with metastasis)和轉(zhuǎn)移淋巴結(jié)(Metastatic lymph node)比無淋巴結(jié)轉(zhuǎn)移食道鱗狀細胞瘤在上面的模塊中具有更高的相關(guān)性?？梢缘贸鼋Y(jié)論:在上述模塊中，有淋巴結(jié)轉(zhuǎn)移食道鱗狀細胞瘤比無淋巴結(jié)轉(zhuǎn)移食道鱗狀細胞瘤對應基因的甲基化程度高。

同時，也可看到胃癌與食管鱗狀細胞瘤在模塊1、4、10都有顯著的相關(guān)性，說明在這幾個基因模塊中，兩者甲基化程度較高。

2.4 基因模塊的功能分析

為了進一步了解上面的基因模塊與癌癥發(fā)生與發(fā)展的關(guān)系，挖掘上面得到的模塊的向量基因。并對這些向量基因進行基因本體功能分類及生物學通路分析。

首先，對與胃癌與食管鱗狀細胞瘤顯著相關(guān)的模塊1、4、10的向量基因進行功能注釋。這些模塊中共得到1 148個基因。

通過DAVID分析，1 148個基因有647個注釋到了189類生物學過程，其余為未知功能基因。設(shè)定閾值為p≤0.05，則基因注釋到96類生物學過程。這些生物學過程主要包括:基因沉默，蛋白質(zhì)降解過程，己糖降解，Wnt受體信號通路，蛋白激酶活性負調(diào)節(jié)等(見表2)。

表2 與胃癌和食管鱗狀細胞瘤顯著相關(guān)的模塊向量基因的功能富集聚簇Table 2 Functional annotation for module vector genes significantly associated with gastric cancer and ESCC

同理，對與前列腺癌顯著相關(guān)的模塊2的向量基因進行功能注釋。注釋的372個基因有204個注釋到了127類生物學過程。設(shè)定閾值為p≤0.05，則基因注釋到79類生物學過程。這些生物學過程主要包括:蛋白激酶活性負調(diào)節(jié)，細胞增殖調(diào)節(jié)，調(diào)控細胞死亡，參與細胞形態(tài)分化等(見表3)。

表3 與前列腺癌顯著相關(guān)的模塊向量基因的功能富集聚簇Table 3 Functional annotation for module vector genes significantly associated with prostate cancer

同理，對與白血病顯著相關(guān)的模塊7、9的向量基因進行功能注釋。注釋的148個基因有141個注釋到了27類生物學過程。設(shè)定閾值為p≤0.05，則基因注釋到14類生物學過程。這些生物學過程主要包括:磷代謝過程;mRNA代謝過程;轉(zhuǎn)錄調(diào)控;磷酸化蛋白質(zhì)氨基酸等(見表4)。

表4 與白血病顯著相關(guān)的模塊向量基因的功能富集聚簇Table 4 Functional annotation for module vector genes significantly associated with leukemia

接著，對與胃癌與食管鱗狀細胞瘤顯著相關(guān)的模塊1、4、10的向量基因進行KEGG通路分析。這些模塊中共得到1 148個基因。

通過DAVID分析，對向量基因進行生物學通路富集分析。采用Fisher精確檢驗，p＜0.05表示一系列基因能代表與某些生物學通路相關(guān)的生物學功能發(fā)生了改變。本次分析中通路發(fā)生改變的主要有:產(chǎn)生癌癥(Pathways in cancer)，產(chǎn)生腎上皮細胞癌(Renal cell carcinoma)(見表5)。

表5 基因模塊1、4、10生物學通路中富集情況Table 5 Gene modules 1、4、10 biological pathway enrichment

通過以上的富集情況，發(fā)現(xiàn)基因模塊1、4、10的相關(guān)基因富集到了產(chǎn)生癌癥的通路。由于基因的啟動子區(qū)域甲基化程度較高，會產(chǎn)生抑制表達的作用。基因表達的缺失導致低氧誘導因子(缺氧誘導因子－α)的積累，從而產(chǎn)生多種生長因子，如血管內(nèi)皮生長因子和血小板衍生生長因子，使細胞運動、細胞轉(zhuǎn)化、防止細胞凋亡等生物學效應的調(diào)節(jié)功能缺失，造成了腫瘤的生成。

同時，模塊1、4、10的基因也富集到了產(chǎn)生腎上皮細胞癌的通路，這也說明了相關(guān)基因啟動子區(qū)域的甲基化程度較高，影響到多種癌癥的發(fā)生。

再對與前列腺癌顯著相關(guān)的模塊2的向量基因進行 KEGG通路分析。通過 DAVID分析，采用Fisher精確檢驗，本次分析中通路發(fā)生改變的是:細胞分裂周期(見表6)。

表6 基因模塊2生物學通路中富集情況Table 6 Gene module 2 biological pathway enrichment

通過上面的富集，基因模塊2的相關(guān)基因富集到影響細胞分裂周期的通路上。有絲分裂是一個重復序列的過程，細胞周期蛋白依賴性激酶(CDKs)是關(guān)鍵調(diào)節(jié)酶，它通過調(diào)控細胞基質(zhì)來控制細胞進程。細胞周期蛋白依賴性激酶抑制劑(CKIs)，如基因CDC7、MAD1L1、CCNB3 等，參與 CDKs的負調(diào)控，從而提供了一個通過該細胞周期負調(diào)控的通路。而它又反過來激活p53抑癌蛋白?；駽DC7、MAD1L1、CCNB3的高甲基化，抑制p53的表達，同時細胞不能進行正常分裂，從而造成腫瘤細胞的產(chǎn)生。

3 討論

目前的研究認為DNA甲基化與腫瘤密切相關(guān)。腫瘤的DNA甲基化改變表現(xiàn)為總體的甲基化水平降低與啟動子區(qū)域CpG島的甲基化水平升高。所篩選的基因模塊的向量基因的甲基化水平普遍較高，就是由于基因的啟動子區(qū)域CpG島的甲基化異常造成的。

通過對基因模塊進行GO功能注釋，發(fā)現(xiàn)了與各癌癥顯著相關(guān)的甲基化異常的基因模塊內(nèi)的相應基因注釋到了諸如基因沉默，Wnt受體信號通路;蛋白激酶活性負調(diào)節(jié);細胞增殖調(diào)節(jié);調(diào)控細胞死亡;參與細胞形態(tài)分化等生物過程，而這些生物學過程又與癌癥的發(fā)生有著顯著的關(guān)聯(lián)。說明這些甲基化異常的基因模塊對腫瘤的發(fā)生與發(fā)展起著重大的作用。

同時，對與胃癌與食管鱗狀細胞瘤顯著相關(guān)的模塊1、4、10的向量基因進行生物學通路富集分析，得到產(chǎn)生癌癥的通路。說明甲基化異常的基因模塊確實與腫瘤的生成有著重要的聯(lián)系。而對于與胃癌與食管鱗狀細胞瘤顯著相關(guān)的模塊1、4、10富集到產(chǎn)生腎上皮細胞癌的通路。也說明了甲基化異常的基因模塊同時與多種癌癥的發(fā)生有著千絲萬縷的聯(lián)系。

在本課題中，首先下載了乳腺導管癌、胃癌、前列腺癌、白血病、肺鱗狀細胞瘤、食管鱗狀細胞瘤等6種癌癥及亞型的DNA甲基化數(shù)據(jù)，經(jīng)過預處理后利用WGCNA篩選出了甲基化基因模塊，通過量化模塊與癌癥表型的關(guān)系發(fā)現(xiàn)了與各癌癥顯著相關(guān)的6個基因模塊。然后，挖掘這些基因模塊的向量基因，對這些基因進行GO功能注釋和KEGG通路富集分析。通過GO功能注釋發(fā)現(xiàn)了基因模塊內(nèi)相應的基因與可能導致腫瘤產(chǎn)生的生物學過程有關(guān);利用KEGG數(shù)據(jù)庫對基因模塊的向量基因進行功能聚類，發(fā)現(xiàn)模塊內(nèi)的基因富集到產(chǎn)生癌癥的通路也說明甲基化異常的基因模塊與癌癥的發(fā)生有著顯著的內(nèi)在關(guān)聯(lián)。同時，也發(fā)現(xiàn)某些甲基化異常的基因模塊(模塊1、4、10)與多種癌癥的發(fā)生有著顯著的關(guān)聯(lián)?；诖耍疚挠兄诎l(fā)現(xiàn)癌癥中的DNA甲基化的生物標記物，為腫瘤的診斷及治療提供可能的靶點。

References)

［1］ WU C，MORRIS J R.Genes，genetics，and epigenetics:A correspondence［J］.Science，2001，293(5532):1103 －5.

［2］ WOLFFE A P.Chromatin remodeling:Why it is important in cancer［J］.Oncogene，2001，20(24):2988 －90.

［3］ PENNISI E.Behind the scenes of gene expression［J］.Science，2001，293(5532):1064－7.

［4］ VALLBOHMER D，BRABENDER J，YANG D，et al.DNA methyltransferases messenger RNA expression and aberrant methylation of CpG islands in non-small-cell lung cancer:association and prognostic value［J］.Clinical Lung Cancer，2006，8(1):39 －44.

［5］ ALVAREZ-VENEGAS R，AVRAMOVA Z.Methylation patterns of histone H3 Lys 4，Lys 9 and Lys 27 in transcriptionally active and inactive Arabidopsis genes and in atx1 mutants［J］.Nucleic Acids Res，2005，33(16):5199－207.

［6］ AHMAD I，RAO，D.N.Chemistry and biology of DNA methyltransferases［J］.Critical Reviews in Biochemistry and Molecular Biology，1996，31(5－6):361－380.

［7］ VERTINO P M，YEN R W，GAO J，et al.De novo methylation of CpG island sequences in human fibroblasts overexpressing DNA(cytosine-5-)-methyltransferase［J］.Mol Cell Biol，1996，16(8):4555 －65.

［8］ AHUJA N，LI Q，MOHAN A L，et al.Aging and DNA methylation in colorectal mucosa and cancer［J］.Cancer Res，1998，58(23):5489 －94.

［9］ WHEELER J M，BECK N E， KIM H C， et al.Mechanisms of inactivation of mismatch repair genes in human colorectal cancer cell lines:the predominant role of hMLH1［J］.Proc Natl Acad Sci U S A，1999，96(18):10296－301.

［10］ BAYLIN S B.Tying it all together:epigenetics，genetics，cell cycle，and cancer［J］.Science，1997，277(5334):1948－9.

［11］ MURAKAMI J，ASAUMI J，MAKI Y，et al.Influence of CpG island methylation status in O6-methylguanine-DNA methyltransferase expression of oral cancer cell lines［J］.Oncol Rep，2004，12(2):339－45.

［12］NEPHEW K P，HUANG T H.Epigenetic gene silencing in cancer initiation and progression ［J］.Cancer Lett，2003，190(2):125－33.

［13］ ZHANG B，HORVATH S.A general framework for weighted gene co-expression network analysis［J］.Stat Appl Genet Mol Biol，2005，4(1):1 －43.

［14］ LANGFELDER P，ZHANG B，HORVATH S.Defining clusters from a hierarchical cluster tree:the Dynamic Tree Cut package for R ［J］.Bioinformatics，2008，24(5):719－20.