基于DNA甲基化的分子亞型可預(yù)測(cè)肝細(xì)胞癌患者的預(yù)后

楊阮阮，李錦忠，龔曉兵

暨南大學(xué)附屬第一醫(yī)院感染科，廣東 廣州 510632

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

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

1 資料與方法

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

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

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

1.4 生存和臨床特征分析 Kaplan-Meier 圖用于說明DNA 甲基化譜圖定義的LIHC 亞組的總體存活率。使用對(duì)數(shù)秩檢驗(yàn)來評(píng)估集群之間差異的顯著性。使用R中的生存軟件包進(jìn)行生存分析。使用R的ggplot2 包進(jìn)行腫瘤分型與臨床特征的相關(guān)性分析。在所有測(cè)試中，P＜0.05被認(rèn)為差異具有統(tǒng)計(jì)學(xué)意義。

1.5 功能富集分析和基因組注釋 將R 中的clusterProfiler 包[10]與 KEGG 結(jié)合使用，對(duì)預(yù)后位點(diǎn)所在的基因組中的基因本體論，生物途徑和調(diào)控基序進(jìn)行基因富集分析。

1.6 預(yù)后預(yù)測(cè)模型的構(gòu)建和評(píng)估 R 中生存包的coxph函數(shù)用于基于CpG位點(diǎn)的甲基化圖譜和預(yù)后信息的組合來構(gòu)建Cox比例風(fēng)險(xiǎn)模型。模型的風(fēng)險(xiǎn)評(píng)分=CpG1*風(fēng)險(xiǎn)系數(shù)CpG1+CpG2*風(fēng)險(xiǎn)系數(shù)CpG2+CpGN*風(fēng)險(xiǎn)系數(shù)CpGN(“*”代表乘法)。Suivival、ROC 曲線和Calibration plot (校準(zhǔn)圖)判斷風(fēng)險(xiǎn)模型的臨床應(yīng)用價(jià)值。

2 結(jié)果

2.1 篩選與OS相關(guān)的潛在預(yù)后甲基化位點(diǎn) 經(jīng)對(duì)患者數(shù)據(jù)進(jìn)行預(yù)處理后，確定了16 381個(gè)甲基化位點(diǎn)。然后使用單變量Cox 回歸分析，在16 381 個(gè)甲基化位點(diǎn)中，有7 917 個(gè)CpG 位點(diǎn)被確定為肝細(xì)胞癌(LIHC)患者OS 相關(guān)的潛在DNA 甲基化生物標(biāo)記物(P＜0.05)。多因素獨(dú)立預(yù)后分析篩選出具有獨(dú)立預(yù)后能力的2 248個(gè)CpG位點(diǎn)。

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

圖1 DNA甲基化分類的選擇標(biāo)準(zhǔn)和一致性矩陣

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

圖3 每個(gè)DNA甲基化亞型的生存曲線

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

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

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

圖5 對(duì)2 248個(gè)CpG位點(diǎn)的注釋基因進(jìn)行KEGG通路富集分析

表1 異常CpG位點(diǎn)的注釋基因的GO功能富集

圖6 每個(gè)DNA甲基化簇的特定的高/低甲基化CpG位點(diǎn)

表2 多變量COX回歸得到的8個(gè)與OS顯著相關(guān)的CpG位點(diǎn)

圖7 高低風(fēng)險(xiǎn)組的熱圖分析

圖8 模型對(duì)臨床生存及預(yù)后的預(yù)測(cè)價(jià)值

3 討論

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

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

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

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