楊阮阮，李錦忠，龔曉兵

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

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

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

1 資料與方法

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

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

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

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

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

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

2 結(jié)果

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 討論

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

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

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

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