基于DNA甲基化模式篩選評估結腸腺癌預后的標記物

劉 陽，王麗茹，張 巖

(哈爾濱工業大學 生命科學與技術學院計算生物學研究中心,哈爾濱 150001)

結直腸癌(Colorectal cancer,CRC)是一種發生率極高的惡性腫瘤[1]。結腸腺癌(Colon adenocarcinoma,COAD)是結直腸癌最常見的病理類型之一(約占95%)[2]。然而目前現有的結腸腺癌的治療手段(包括手術和化療等綜合治療)仍然不能達到令人滿意的效果，結腸腺癌的五年生存率依然不容樂觀[3]。大部分患者發現和初診為結腸腺癌時已處于癌癥晚期，因此，現階段的醫療診治過程中晚就診、缺乏可靠的生物標志物和不精準的治療靶點已成為治療結腸腺癌的主要障礙[4]，而早期診斷和治療對于改善患者的預后和生活質量至關重要。最新的研究表明，除了遺傳變異以外，擾亂的、不平衡的表觀遺傳基因組也是癌癥發生發展的重要原因[5]。DNA甲基化是一種重要的調控基因表達的表觀遺傳修飾，與癌癥的發生和發展緊密相關[6]。異常的DNA甲基化通過激活癌基因和/或使腫瘤抑制基因失活誘導癌癥的發生發展[7]。通過因此本課題旨在通過結腸腺癌DNA甲基化譜，分析患者的DNA甲基化模式并識別出與預后相關的結腸腺癌亞型。

1 材料與方法

1.1 DNA甲基化數據的獲取和預處理

從癌癥基因組圖譜(The Cancer Genome Atlas，TCGA，https://www.cancer.gov/about-nci/organization/ccg/research/structural-genomics/tcga)數據庫中下載了結腸腺癌樣本DNA甲基化數據和這些樣本對應病人的臨床信息，甲基化數據全部由Illumina Infinium HumanMethylation450 BeadChip芯片平臺產生。選取324個記錄了患者生存狀態的結腸腺腫瘤樣本數據進行接下來的分析。研究表明，發生在DNA啟動子區的異常甲基化對癌癥的發生發展起著重要調控作用。在本研究中啟動子區域定義為轉錄起始位點上游2 kb至下游0.5 kb。首先選取癌癥樣本數據中屬于啟動子區的位點，去除其中的SNP位點、對應于性染色體上的位點和一個位點對應多個基因的位點。接下來去除70%以上的樣本中有缺失的位點，最后使用k近鄰算法(knn)補全其余的缺失值，該算法使用“impute”R包中的impute.knn()函數實現，其中參數k取10、maxp取3 000。

1.2 提取分類特征

為尋找在癌癥和癌旁樣本中甲基化水平有差異且差異具有顯著性的CpG位點，計算每個位點在癌癥和癌旁樣本中的甲基化水平平均值，求出均值差并用t檢驗方法計算該差異的顯著性。本研究將均值差絕對差異>0.1且經Benjamini-Hochberg多重檢驗校正方法校正后p<0.05的CpG位點視作差異甲基化位點。將324個樣本按7∶3的比例劃分為訓練集與測試集，分組原則為隨機分組且死亡率相似。最后分為226個樣本的訓練集和98個樣本的測試集。

為獲得與預后相關的結腸腺癌分子亞型，將使用顯著影響生存的CpG位點作為分類特征。首先，使用每個樣本所對應患者的生存狀態和生存時間為每個差異甲基化CpG位點和患者的性別、年齡、臨床分期構建單變量COX比例風險回歸模型，分析每個CpG位點的甲基化水平、年齡、階段與預后的關系。再將單變量中顯著影響預后的臨床因素作為協變量引入多變量模型篩選得到獨立影響預后的CpG位點。對于每個CpG位點i，單變量和多變量COX比例風險回歸模型的公式定義如下：

h(t,x)i=h0(t)exp(βmethymethyi)

(1)

h(t,x)i=h0(t)exp(βmethymethyi+βstagestage)

(2)

式中h0(t)是基準風險方程，可以是任意一個針對時間t的非負方程，methyi是代表CpGi甲基化水平的矢量，stage是代表患者臨床分期的矢量，βmethy和βstage代表回歸系數。該模型使用Benjamini-Hochberg多重檢驗校正方法對P值進行校正。

1.3 識別預后相關的DNA甲基化亞組

在指定聚類算法和度量距離時，使用ConcensusClusterPlus包內嵌的多種聚類算法和距離進行嘗試，以探尋最優的聚類結果。根據輸出的結果文件，使用歐幾里得距離為度量的K均值聚類算法被納入考慮范圍。

分組數k的確定是本研究要解決的一個重要問題，對從k=2到k=10的不同分組數，應用上一步的聚類方法，除了固有的輸出結果外，對不同分組數下類別的平均一致性和一致性變異系數進行了計算。最終用于確定聚類方法和聚類數的標準是：在某一聚類數k下，組間的一致性數值相對較高、變異系數相對較低，并且CDF曲線下面積在相鄰的兩個類別數量間的改變趨勢較緩。變異系數根據下式計算：

(3)

式(3)中SD是同一分類數下不同分組間一致性的標準差，MN是同一分類數下不同分組間的平均一致性。根據上述標準，最終選擇使用K均值算法，以歐幾里得距離作為相似度指標的聚類方法，確定聚類數目k為7。

1.4 不同甲基化亞組間的預后分析

通過與預后相關的CpG位點區分出了七個亞組，自然應該考察這七個DNA甲基化亞組間的預后情況。對得到的甲基化亞組進行Kaplan-Meier生存分析。這一分析通過“survival”R包中的函數survfit()和survdiff()完成。并使用對數秩檢驗(log-rank test)分析生存差異的統計學顯著性。

1.5 各甲基化亞組間的特征描述

目前結腸腺癌還沒有通過臨床病理特征定義的細分亞型，在分析并確定了不同DNA甲基化亞型間的預后差異之后，進一步分析各個甲基化亞型與臨床特征的相關性，所有檢驗方法若無特殊說明均采取雙側檢驗。

根據臨床數據的情況，使用費希爾精確檢驗(Fisher’s exact test)分析各臨床特征在不同甲基化亞型間的分布是否有顯著差異。

在每個甲基化亞型內，對每個臨床特征和總體之間進行超幾何檢驗，以分析各臨床特征在各甲基化亞型內的富集是否具有顯著性。

1.6 篩選甲基化亞型的特異性標記

使用一種基于香農熵模型可從全基因組甲基化譜中鑒定差異甲基化區域(DMR)的軟件QDMR[12](Quantitative Differentially Methylated Regions)，從用于識別七個甲基化亞組的137個CpG位點中篩選每個亞組的特異性甲基化標記。由于在前述亞組確定的過程中采用了一致性聚類的方法，這意味著各個亞組內的樣本有著相似的甲基化模式。因此，本研究計算用于一致性聚類的137個特征位點在七個亞組中的甲基化均值，以該均值來表征各位點在不同亞組中的甲基化狀態，將得到285×6維矩陣作為QDMR的輸入。在衡量特異性的過程中，參數SD被設置為0.07。

1.7 構建和評估甲基化亞型預測模型

本研究的DNA甲基化亞型由前述無監督的一致性聚類方法得到，為了驗證DNA甲基化分型結果并建立更為便捷精確的結腸腺癌甲基化分型方法，使用有著7個甲基化亞型標簽的訓練集構建有監督的SMO分類模型。在98個樣本的測試集中使用1.2中得到的預測模型，根據樣本的甲基化數據將患者劃分到7個DNA甲基化亞型中，也就是使用SMO分類模型將訓練集得到的7個甲基化亞型標簽分配給測試集中的樣本。然后對測試集中7個甲基化亞型進行生存分析，驗證模型的穩定性和準確性。

2 結果和分析

2.1 結腸腺癌DNA甲基化特征的篩選

本研究將均值差絕對差異>0.1且p<0.05的CpG位點視作樣本對之間顯著差異的CpG位點。最終得到26 158個差異甲基化位點，這一結果可視化使用R包ggplot(見圖1)。

圖1 結腸腺癌差異甲基化位點火山圖Fig.1 Volcano plot of differential methylation sites of colon adenocarcinoma注：橫坐標是癌癥癌旁樣本間不同CpG位點的均值差，縱坐標是經過多重檢驗校正后的p值負對數。在均值差閾值和P值分別設置為0.1和0.05的情況下，紅色區域為差異且顯著的低甲基化位點，藍色區域是差異且顯著的高甲基化位點。

為獲得與預后相關的結腸腺癌分子亞型，使用顯著影響存活的CpG位點作為分類特征。首先將所有結腸腺癌樣本按照7∶3的比例劃分為訓練集和測試集，分別包含226個樣本和98個樣本。接著，為訓練集中的樣本基于生存時間和存活情況構建COX比例風險回歸模型。單變量COX模型獲得的2 838個與預后顯著相關位點被用于多變量COX模型構建，3個顯著的臨床因素“分期”(gender:p=0.017,age:0.004,stage:0.016)被作為協變量也引入其中，最終多變量的回歸模型獲得137個依舊顯著的位點，部分位點所對應的基因(見表1)。

表1 特異性位點基因對應表Table 1 Characteristic sites and corresponding gene symbols

P值使用Benjamini-Hochberg方法進行了多重檢驗校正，Hazard.Ratio代表著COX模型的風險比，風險比是相對而言的，風險比大于1的位點被認為與不良預后相關，小于1的位點被認為與良好預后相關。137個CpG位點將被用作分組特征。

2.2 基于一致性聚類的結腸腺癌亞組識別

進一步對所獲得的137 個CpG位點進行一致性聚類以獲得與預后相關的結腸腺癌DNA甲基化亞型。本研究計算每個類別的平均聚類一致性和聚類之間的變異系數，結合ConcensusClusterPlus包的輸出結果(見圖2)，根據1.3中描述的方法來確定選擇最優的聚類算法和類別數k。

圖2 基于K均值算法的結腸腺癌亞型分類標準Fig.2 Classification standards for colon adenocarcinoma subtypes based on k-means algorithm注：(a)表示每個類別號k與k -1 相比的累積分布函數(CDF)曲線下面積的相對變化。橫軸表示類別編號k，縱軸表示CDF曲線下面積的相對變化。(b)中紅線代表平均類一致性，藍線代表類間的變異系數。

對于K均值算法，CDF曲線下面積改變趨勢從聚類數目k=6開始明顯趨緩且分類數為7時，其平均聚類一致性曲線和變異系數曲線有一個明顯的拐點，擁有相對高的一致性系數和相對低的變異系數。因此，對于K均值算法，認為它在聚類數為7時具有較高的聚類穩定性，k=7是最適聚類數。最終確定使用聚類數k=7時K均值算法的聚類結果用于下一步的分析。

2.3 七個DNA甲基化亞組的預后分析

因為使用的是與預后相關的CpG位點區分的亞組，進一步本研究對7個亞組進行Kaplan-Meier生存分析，并使用對數秩檢驗(log-rank test)研究生存差異的統計學顯著性(見圖3)。與其他癌癥類型相比較而言，結腸腺癌屬于惡性程度高、預后狀況差的一類腫瘤，目前還沒有應用于臨床的可以區分腫瘤預后的分型方法。在本文的DNA甲基化亞組中，生存分析顯示這7個DNA甲基化亞組主要分為兩大組，之間的預后差異存在顯著性。三年(1 095 d)存活率分析結果表明，cluster1、cluster3、cluster4、cluster6生存時間明顯優于其他亞型。

圖3 訓練集中甲基化亞組的生存曲線Fig.3 Survival of methylation subgroups in training set注：橫坐標是存活時間(d)，縱坐標是存活概率.

2.4 不同甲基化亞型的臨床特征描述

對各個亞型中樣本的臨床特征進行分析，臨床特征年齡、腫瘤分期、T分期、N分期、M分期、組織學類型等較為經典的6個指標，同時，歷年來有多篇NCS文章表明微衛星不穩定性與結腸腺癌的發生發展密切相關，故也將其選出作為結腸腺癌的重要臨床特征。其中T分期是對原發腫瘤的評估，隨著腫瘤體積的增加和鄰近 組織受累范圍的增加，依次用T1～T4來表示；N分期是對區域淋巴結的評估，本研究中患者該項指標被納入分析的分期有N0和N1，分別表示著沒有區域淋巴結轉移和1～3枚區域淋巴結轉移；M分期是對腫瘤是否有遠端轉移的評估；腫瘤解剖學細分評估了腫瘤位置，不同位置的結腸腺癌腫瘤侵犯范圍、術式選擇等方面存在差異。

本研究在7個DNA甲基化亞型之間，使用Fisher’s exact test分析各臨床特征在不同亞型之間的分布差異是否具有顯著性(見表2)。結果顯示：年齡(p=0.017)、組織學類型((p=0.039)、微衛星不穩定性(p=0.003)、N分期((p=0.038)在不同亞型樣本中的分布差異是顯著的。這說明基于DNA甲基化譜分出的結腸腺癌亞型間的預后差異，一定程度上可由患者的年齡、發生部位、淋巴結受累情況、微衛星不穩定性去解釋。而臨床分期、腫瘤的大小和腫瘤轉移情況，對于甲基化亞型的區分沒有顯著的指導意義。

表2 甲基化亞型間Fisher’s exact test結果Table 2 Fisher’s exact test results among methylation subtypes

在確定了不同亞型間的臨床特征有著差異性后，進一步分析臨床特征在各亞型內的富集情況是否具有顯著性。對每個臨床特征在各亞型內的分布情況和在總體中的分布進行超幾何檢驗，檢驗結果(見表3)。

表3 超幾何檢驗各特征在亞型內富集的顯著性Table 3 Hypergeometric test for significance of feature enrichment in subtypes

根據上一步對亞型間進行Fisher’s exact test的結果，將7個甲基化亞型和基于這兩個臨床特征下的樣本分組進行比較(見圖4)。

圖4 DNA甲基化亞型和N分期類型、組織學類型、年齡、微衛星不穩定性相互富集Fig.4 Enrichment of DNA methylation clusters,N stage,histology type, age,and microsatellite instability注：A.N分期類型對應于DNA甲基化亞型。B.DNA甲基化亞型對應N分期。C.組織學類型對應于DNA甲基化亞型。D.DNA甲基化亞型對應組織學類型。E.年齡對應于DNA甲基化亞型。F.DNA甲基化亞型對應年齡。G.DNA甲基化亞型對應微衛星不穩定性。H.微衛星不穩定性對應于DNA甲基化亞型。

2.5 不同甲基化亞型特異性標記的確定

提取訓練集樣本中通過QDMR得到的79個特征CpG位點的甲基化數據，部分位點對應的基因(見表4)。生成特征位點甲基化矩陣。使用R包pheatmap為該矩陣繪制熱圖，并為它從0-1的甲基化值關聯上藍色到紅色。每個類別中特性低甲基化位點在熱圖中顯示出藍色，特異高甲基化位點顯示出紅色(見圖5)。這79個特征位點可以作為結腸腺癌中不同DNA甲基化亞型的特異性DNA甲基化標記，代表著每個亞型獨特的DNA甲基化模式(見圖5)。class1具有最大數量的特異性低甲基化位點，class2具有最大數量的特異性高甲基化CpG位點。

圖5 七個甲基化亞型中特征CpG位點甲基化水平Fig.5 CpG site methylation level among seven subtypes

表4 特異性位點基因對應表Table 4 Characteristic sites and corresponding gene symbols

2.6 構建DNA甲基化預測模型

為了驗證本研究獲得的甲基化亞型及更為準確高效地對樣本進行分類，使用訓練集樣本構建基于決策樹學習有監督的SMO分類模型(見表5)，所得到的模型有著87.61%的分類準確性，每一行代表模型預測的樣本類別，每列代表樣本真實的類別。C1-C7對應著class1-class7，表中對角線上的數字表示每個class中預測類別與實際類別相符的樣本數量。

表5 SMO分類器的混淆矩陣Table 5 Confusion matrix of SMO classifier

本研究的分類數量為7，而ROC曲線通常針對二分類的研究，故將7個亞型的真實值和模型的預測值轉化為7個二分類對象以適應ROC曲線特征(見圖6)，使用訓練集樣本訓練出的該模型對6個亞型的分類都具有較高的AUC值(ROC曲線下面積),而模型對class6的特異性和敏感性則相對較低。

圖6 預測模型的敏感性和特異性Fig.6 Sensitivity and specifity of prediction model

將訓練集中得到的模型應用到預留的檢驗集中，對檢驗集中未知類別標簽的患者進行DNA甲基化亞型的分類。檢驗集中98個樣本根據特征CpG位點的甲基化模式被分配到了7個不同的組內。對不同分組的3年(109 5 d)存活率進行分析，并使用對數秩檢驗(Log-rank test)分析發現生存差異具有較高的顯著性(p=0.052)。生存曲線(見圖7)。

圖7 檢驗集生存分析結果Fig.7 Survival analysis results of test set

3 討 論

誘發癌癥的因素復雜而多樣，眾所周知，目前的主流觀點認為癌癥由于基因突變引起[10]，是一個累積漸進的過程，遺傳改變的積累導致了惡性腫瘤的發生和發展。但是最新的研究表明，除了遺傳變異以外，異常的DNA甲基化也在癌癥的發生和發展中起到重要作用[11]。目前，異常的DNA甲基化在胃癌中是被廣泛研究的一種失調的表觀遺傳機制。例如，有研究證明胃癌中一些腫瘤抑制基因或者腫瘤相關基因(如p16、RUNX3、MLH1和CDH1等)會被啟動子的甲基化沉默[12]。這些研究都驗證了DNA甲基化作為腫瘤標志物的重要價值。在確定了DNA甲基化與癌癥的發生和進展相關聯后，甲基化模式的動態性可能有助于癌癥的早期診斷和評估：異常甲基化位點增加或減少的時間趨勢可以幫助預測惡性腫瘤轉化的速率和概率。因此，異常DNA甲基化被認為是可用來評估癌癥前期進展的潛在早期診斷生物標志物[13]，是臨床實踐中癌癥診斷的理想靶標。此外，識別癌癥特異或亞型特異的生物標志物也具有預測價值。

本研究中包含了大量的結腸腺癌樣本的 Infinium HumanMethylation450 BeadChip 陣列數據集從 TCGA 數據庫下載，這些 數據可以用于本文的結腸腺癌異質性分析。龐大的樣本量使得本研究能夠更全面地探索結腸腺癌的分子亞型和分子異質性。

許多研究表明，表觀遺傳修飾(DNA 甲基化)在結腸腺 癌中的早期檢測，以及改善分子分類，預后和輔助治療中發揮關鍵作用，體現了 基因組范圍的分子水平分析在精準醫學時代具有重要的生物學和臨床意義[14]。利用基因啟動子區域內篩選的預后相關的 CpG 位點進行無監督的聚類分析，通過一致性聚類獲得七個預后不同的分子亞型，深入分析發現不同亞型之間有著分子或臨床特征上的差異，這證實了結腸腺腫瘤的異質性及對結腸腺癌詳細分類的必要性。

一致性聚類與其他無監督聚類方法(如層次聚類)相比，提供了聚類的類別 數量的選取參考依據[15]。本研究首次提出建議根據 DNA 甲基化水平將結腸腺腫分為七個預后分子亞型。這種詳細程度帶來了較高的類內一致性，可更好地指導個性化醫療。在本研究中，SMO模型的甲基化亞型分類能力有著較高的穩定性。但是，第六亞型的AUC曲線太直，分析可能是由于樣本量過少，過擬合導致。同時，在模型運用到檢驗集時，在訓練集和檢驗集中，class6，class5，class3，class2的生存曲線差異較大，其可能是由于分類器對class6的靈敏性和特異性相對較低，導致預測類別有誤所致。觀察混淆矩陣可注意到實際類別為class6的不少樣本被錯誤劃分為class3，class5和class2，從而可以解釋檢驗集中class3，class5和class2生存率提升的現象，這說明在一定程度上本研究建立的分類模型對于結腸腺癌DNA甲基化亞型的預測與預后情況相關聯。值得注意的是，本研究建立的模型為各個甲基化亞型分配的樣本數量占總樣本數的相對數目多少在訓練集和檢驗集中是相似的，例如：分配到class3的樣本數目在訓練集和檢驗集樣本數均是最多。這也從一定程度上說明本研究建立的模型對甲基化亞型分類能力的穩定性。

在使用 QDMR的分析中，發現了79個在亞型中特異性高/低甲基化的 CpG 位點，對應于76個基因，它們定義了結腸腺癌的特定 DNA 甲基化亞組。這些位點可以被視為診斷結腸腺癌的精準醫療的靶標和生物標志物。在這些特異性的CpG位點對應的基因中，需多先前已報道與結腸腺癌有關。其中，Sfrp5基因的甲基化被認為與結腸腺癌中致癌物誘導過程WNT基因的響應有關[16]，Thbs4基因是大腸癌中一個與年齡相關的甲基化腫瘤抑制基因[17]，GDF10基因與結腸癌發生密切相關[18]。

總之，使用TCGA 數據庫中的結腸腺腫瘤數據識別出七種不同的 DNA 甲基化預后分子亞型，這些亞型在分子水平或臨床特征上也存在著顯著差異。這可以對單一腫瘤內部的異質性更詳細地解釋。同時，該研究方法也可以推廣到其他腫瘤的分型中。

4 結 論

使用TCGA中結腸腺癌的DNA甲基化數據，篩選到137個與預后相關的CpG位點。基于一致性聚類方法識別出7個與預后相關的亞型，這些亞型的預后差異顯著且亞型特征可由年齡、N分期、微衛星不穩定性、解剖學發生部位反映。在不同的結腸腺癌亞型中篩選出79個特異性高/低甲基化位點，這些位點代表著每個亞型獨特的甲基化模式。利用結腸腺癌亞型特異性DNA甲基化特征構建了基于序列最小最優化的分類模型，該模型對結腸腺癌亞型的分類準確性達到87.61%，提供了用于結腸腺癌分型的精確的DNA甲基化標記。