基于DNA甲基化的分子亞型可預測肝細胞癌患者的預后

楊阮阮，李錦忠，龔曉兵

暨南大學附屬第一醫院感染科，廣東 廣州 510632

原發性肝癌是世界上第六大常見的癌癥，也是導致癌癥死亡的第三大原因，最常見的組織學是肝細胞癌(hepatocellular carcinoma，HCC)[1]。改善 HCC 患者預后一個重要阻礙是對HCC 的異質性不完全了解。目前的觀點認為，每個原發性腫瘤的變異可能由遺傳和表觀遺傳上不同的細胞組成，這導致每種腫瘤類型具有特定的表型異質性[2]，在這些腫瘤細胞中有一些被稱為腫瘤干細胞(cancer stem cells，CSC)，其負責產生異質性腫瘤病變，也有可能助于治療腫瘤的復發和轉移。因此，通過對表觀遺傳的分子亞型進行臨床特征及預后的分析對了解腫瘤異質性具有重要意義。

表觀遺傳是基因表達的可遺傳狀態而不改變DNA序列，甲基化是表觀遺傳學中最重要和最常見的修飾，也是調節基因組功能的重要手段[3]。全基因組甲基化水平降低是癌癥早期的一種顯著標志，并可能與癌癥的嚴重程度及轉移有明顯關聯。而特定基因啟動子區域的CpG 島異常高甲基化是癌癥另一重要現象，導致染色體結構的改變并使腫瘤抑制基因和其他癌癥相關基因沉默，使癌細胞能適應微環境并促進腫瘤發生和發展[4-5]。SHIN 等[6]分析了 HCC 和非腫瘤肝組織的甲基化譜，發現了HCC 中新的高甲基化基因，如TNFRSF10C 和IRF5 (參與細胞凋亡)，HOXA9(胚胎發育過程中受調控)和NPY(參與細胞運動和細胞增殖)。此外，在腫瘤組織69 個顯著低甲基化的基因中，已經鑒定出IL16 (趨化因子介導的免疫相關基因)，BLK 和PGR 癌基因以及組蛋白去乙酰化酶1(HDAC1)。然而，尚未確定這些基因的啟動子區域中的特定甲基化序列。因此，本文基于高通量組學數據整合了HCC的多個DNA甲基化生物標志物，開發了一種預后預測模型，以改善臨床預后評估和個性化治療。

1 資料與方法

1.1 數據獲取和預處理 從TCGA 數據(https://cancergenome.nih.gov/)下載 level 3 級的 RNA-seq 數據(FPKM)，包括374例肝細胞癌(Liver hepatocellular carcinoma，LIHC)患者樣本和50 例癌旁樣本。從UCSC(https://xenabrowser.net/)下載 430 例 LIHC 的 Illumina Infinium HumanMethylation450 陣列的甲基化數據。用 Perl(5.26)語言(http://www.perl.org/)的合并及 ID 轉換腳本提取DNA甲基化矩陣和轉錄組矩陣，本研究僅包括來自臨床隨訪時間超過30 d 的樣本數據。每個位點的甲基化水平由β值表示，其范圍為0(未甲基化)至1(完全甲基化)。使用sva R package[7]中的ComBat算法，通過整合所有DNA甲基化陣列數據，合并batch和patient信息來去除batch效應。

1.2 預后相關的甲基化位點篩選 首先，使用單變量Cox模型來計算每個異常甲基化基因的甲基化水平與患者總體存活(overall survival，OS)之間的關聯，當P＜0.05時，可以認為基因在單變量Cox分析中具有統計學意義[8]。然后，使用TMN 分期，病理分期(Stage)，組織學分級(Grade)，年齡和性別將從單變量COX比例風險回歸模型中獲得的顯著CpGs引入多變量COX比例風險回歸模型中，這些變量在單變量模型中也很重要。最后，選擇在單變量和多變量Cox 回歸分析中均顯著的CpG位點作為特征性CpG位點。

1.3 甲基化腫瘤分型 使用R 中的Consensus Cluster Plus軟件包[9]進行一致性聚類。根據最多可變的CpG位點識別LIHC子組。對應于一致性聚類的熱圖由pheatmap R 包生成。使用顏色漸變表示從0(白色)到1(深藍色)的共識值；排列矩陣，使得屬于同一簇的項目彼此相鄰。在這種布置中，對應于完美共識的矩陣將顯示顏色編碼的熱圖，其特征是沿著白色背景上對角線的藍色塊。

1.4 生存和臨床特征分析 Kaplan-Meier 圖用于說明DNA 甲基化譜圖定義的LIHC 亞組的總體存活率。使用對數秩檢驗來評估集群之間差異的顯著性。使用R中的生存軟件包進行生存分析。使用R的ggplot2 包進行腫瘤分型與臨床特征的相關性分析。在所有測試中，P＜0.05被認為差異具有統計學意義。

1.5 功能富集分析和基因組注釋 將R 中的clusterProfiler 包[10]與 KEGG 結合使用，對預后位點所在的基因組中的基因本體論，生物途徑和調控基序進行基因富集分析。

1.6 預后預測模型的構建和評估 R 中生存包的coxph函數用于基于CpG位點的甲基化圖譜和預后信息的組合來構建Cox比例風險模型。模型的風險評分=CpG1*風險系數CpG1+CpG2*風險系數CpG2+CpGN*風險系數CpGN(“*”代表乘法)。Suivival、ROC 曲線和Calibration plot (校準圖)判斷風險模型的臨床應用價值。

2 結果

2.1 篩選與OS相關的潛在預后甲基化位點 經對患者數據進行預處理后，確定了16 381個甲基化位點。然后使用單變量Cox 回歸分析，在16 381 個甲基化位點中，有7 917 個CpG 位點被確定為肝細胞癌(LIHC)患者OS 相關的潛在DNA 甲基化生物標記物(P＜0.05)。多因素獨立預后分析篩選出具有獨立預后能力的2 248個CpG位點。

2.2 不同的DNA甲基化預后亞組和簇間預后分析 使用2 248個潛在的預后甲基化位點的一致性聚類來識別LIHC 的不同DNA 甲基化分子亞組以進行預后分析。根據以下標準確定簇的數量：簇內的一致性較高，變異系數相對較低，CDF 曲線下的面積沒有明顯增加(圖1A、1B)。根據類別數計算了平均聚類一致性和聚類之間的變異系數。累積分布函數(CDF)曲線下的面積在7 個類別后開始趨于穩定(圖1B)，最后確定最佳簇數為7，并制作了相關樹狀圖表示定義明確的7 塊結構及其一致性(圖1C)，另外TMN 分期，Stage，Grade，年齡，性別和 DNA 甲基化亞組作為注釋，利用heatmap 函數生成與上述k=7 樹圖相對應的heatmap 圖譜(圖 2)。Kaplan-Meier 生存分析顯示 7 個組之間的預后差異有顯著統計學意義(P＜0.001)，聚類5和2的預后最佳，而聚類7和6的預后最差(圖3)。然后，分別根據TMN 分期，Stage，Grade 和年齡、性別分析了7 個群集的群集內比例(圖4)。相關臨床特征與不同組聚類間的關聯趨勢如下：聚類4、6 和7 具有較高的病理分期階段(Stage)；T分期較低的集群為1、2和5；N 等級較高的集群為4 和6；性別總體以男性為主。本研究發現，不同臨床特征在7 個群集中的所占比例各不相同，而且與集群的預后特征也相匹配(如聚類5的T分期較低而預后較好)。

圖1 DNA甲基化分類的選擇標準和一致性矩陣

圖2 以DNA 甲基化分級、TNM 分期、病理分期、組織學分級、性別和年齡為注釋的heatmap圖譜

圖3 每個DNA甲基化亞型的生存曲線

2.3 基因功能富集分析和基因組注釋 利用上述2 248個CpG位點的基因組注釋，共鑒定出2 438個相應的啟動子基因。然后對這2 438個基因進行功能富集分析，結果表明：在生物途徑(biology process，BP)，這些基因主要涉及ncRNA代謝過程、有絲分裂細胞周期相變的調控、細胞周期相變的調節、核糖核酸分解過程、信使核糖核酸分解過程。分子功能(molecular function，MF)主要參與催化活性，作用于RNA、核酸外切酶活性、單鏈DNA結合、核酸酶活性、催化活性，作用于tRNA。此外，細胞組成(cellular component，CC)主要涉及染色體區域、核斑點、染色體端粒的地區、濃縮的染色體和中心體(表1)。KEGG共發現24個相關的PATHWAYs(P＜0.01)，其中最集中且有意義的途徑是RNA 運輸、細胞周期、p53 信號通路和剪接體等(圖5)。

2.4 構建和評估LIHC 預后預測模型 聚類5包含大量的樣本和最多的特異性甲基化位點(圖6)，并且預后良好，因此選擇它作為種子聚類。多變量Cox 回歸用于構建腫瘤預后風險模型并發現8 個預后相關的甲基化位點：cg05489143、cg09600437、cg19165652、cg19569208、cg22732432、cg22958262、cg24153171、cg25455598 (表 2)，由它們構建的風險模型可用于預后評估。獲得的腫瘤預后風險模型是：風險評分=4.98*cg05489143-21.18*cg09600437+3.50*cg19165652+459*cg19569208+11.08*cg22732432+5.07*cg22958262-1.6.02*cg24153171 + 4.75*cg25455598。此外，根據上述模型計算的中位風險評分截止值，共有204 例患者被分為高風險組(n=101)和低風險組(n=102)，LIHC的高低風險與甲基化程度的關系如圖7，可以看到隨著風險評分的增高，高風險組的生存時間雖然沒有明顯下降，但是死亡率明顯升高。基于高、低風險分組采用Kaplan-Meier 方法繪制生存曲線(圖8A)，高低風險兩組的預后差異具有顯著統計學意義(P＜0.001)。同時，ROC曲線的AUC值為0.822(圖8B)，提示該模型可以很好地預測患者的存活率。

圖4 DNA甲基化簇之間的年齡、性別、TNM分期、病理學分期和組織學分級占比

圖5 對2 248個CpG位點的注釋基因進行KEGG通路富集分析

表1 異常CpG位點的注釋基因的GO功能富集

圖6 每個DNA甲基化簇的特定的高/低甲基化CpG位點

表2 多變量COX回歸得到的8個與OS顯著相關的CpG位點

圖7 高低風險組的熱圖分析

圖8 模型對臨床生存及預后的預測價值

3 討論

近年來HCC的治療取得了較大的突破，但肝癌的累積復發率分別在3年、5年時仍高達50%～60%、70%～100%[11]。盡管索拉非尼是目前HCC 最有效的靶向藥物，也是唯一一種在晚期HCC使用的全身性靶向腫瘤治療藥物，但其療效非常有限[12]。在這個前提下，HCC 中出現的表觀遺傳改變可能會成為一種新的治療靶點[13]。DNA 甲基化作為表觀遺傳學中最重要和最常見的修飾，它是一種共價化學修飾，表現為胞嘧啶環的5位碳上加入甲基(CH3)，這常見于5'CG3'這個基因序列中[14]。

本文研究了LIHC患者癌癥樣本和癌旁樣本之間的異常甲基化位點，首先通過單因素和多因素COX比例風險回歸分析最終篩選出具有獨立預后能力的CpG 位點，接著利用一致性聚類方法進行腫瘤分型，最后根據分型差異分析構建甲基化位點的預后模型。根據該模型計算可以看到隨著風險評分的增高，高風險組的生存時間雖然沒有明顯下降，但是死亡率明顯升高。另外，高低風險兩組的預后具有顯著差異，同時，ROC 曲線的AUC 值為 0.822，表明模型可以很好地預測患者的存活率。有研究表明HCC 腫瘤會表現出與危險因素、腫瘤分期、分化程度和癌癥治療后存活相關的特異性DNA甲基化特征[15]，并且大多數CpG 位點會傾向于從T1 期到T3 期逐漸高甲基化[16]，本研究也發現臨床分期等臨床參數與集群的預后特征相匹配。近年來，已經有很多關于基因甲基化對各種癌癥如肺癌[17]，乳腺頸癌[18]和卵巢癌[19]診斷，治療和預后評估的實驗及臨床研究。也有關于肝癌的相關研究，如CARM1介導的GAPDH甲基化是肝癌中葡萄糖代謝的關鍵調節機制[20]。另外，有實驗證明P14ARF mRNA 水平受原發性肝癌中的DNA 甲基化調節，P14ARF基因DNA甲基化可能與HCC的發生及TNM分期有關[21]。與基因突變不同，表觀遺傳的改變具有可逆性，特別是DNA甲基化和組蛋白修飾[22]。通過使用去甲基化劑如5-氮雜胞苷(阿扎胞苷)和5-氮雜-2'-脫氧胞苷，可以實現DNA甲基化基因在癌細胞系中的重新表達，臨床上在骨髓異常增生綜合征和急性髓性白血病中的應用較常見和成熟[23-24]。因此，對基因甲基化的研究對未來肝癌的診療意義重大。

近年來，使用綜合基因組工具的大規模基因組學和全基因組研究重塑了對癌癥進化和異質性的理解。例如，RAN TAO 等[25]通過 HBV 相關 HCC 的全基因組甲基化譜發現7 個新基因(WNK2、EMILIN2、TLX3、TM6SF1、TRIM58、HIST1H4F 和 GRASP)在HCC 中高甲基化，在成對的相鄰肝組織中低甲基化。SHEN等[26]也發現在HCC組織中有684個CpG位點顯著高甲基化，這些基因中的5 個(CDKL2、CDKN2A、HIST1H3G、STEAP4、ZNF154)在高達63%的患者血漿中具有可檢測的高甲基化DNA。鑒定的甲基化基因組可以是用于早期診斷的潛在生物標志物。但是以上這些研究只是從單個基因或多個基因的DNA 甲基化程度監測去評估在肝癌中的價值，肝癌的發病機制復雜，所以聯合基因表達和DNA甲基化位點進行綜合分析，并通過甲基化位點構建出了肝細胞癌的預后風險模型，把評估進行了量化分析，通過這個模型可以比較準確肝癌患者的生存情況，臨床實用性更強。

綜上，本研究在癌組織和癌旁組織中成功篩選出2 248個差異甲基化位點，DNA甲基化水平的差異與T分期、N 分期、M 分期、年齡、Stage、Grade 和預后的差異相關。接著利用一致性聚類方法得到7 個腫瘤亞組，亞組之間的預后有顯著差異，同時對差異甲基化位點相應的啟動子基因進行功能富集分析，它們主要涉及RNA運輸、細胞周期、p53信號通路和剪接體，為肝癌發生機制提供了理論基礎。最后根據分型差異分析構建甲基化位點的預后模型，該模型可以很好地預測患者的存活率。