多基因模型在肝細胞癌預后中的應用

魏之菡 法博濤 俞章盛

（上海交通大學生命科學技術學院生物信息學與生物統(tǒng)計學系 上海交通大學-耶魯大學聯(lián)合生物統(tǒng)計與數(shù)據(jù)科學中心，上海 200240）

肝癌是全球惡性腫瘤患者死亡的主要原因之一，在我國2017 年腫瘤死亡排名中高居第二位，而肝細胞癌（Hepatocellular carcinoma，HCC）在原發(fā)性肝癌中占比高達70%-90%［1-3］。近年來，HCC 的治療措施有了很大的改善，包括手術切除、射頻消融術、肝移植、栓塞化療術、靶向藥物治療等，但是HCC患者的5 年生存率仍然偏低［4-5］。另一方面，處于同一腫瘤分期的HCC 患者可能有不同的預后表現(xiàn)，故僅基于臨床指標來預測HCC 患者預后可能不盡人意［6］。因此，引入分子水平信息對患者進行更加精準的風險分型十分必要，將促進精準醫(yī)療的開展。

目前，利用分子標志物進行HCC 風險分型預測的模型的大部分研究利用基因差異表達分析篩選出的mRNA 或miRNA 建立風險打分模型進行HCC 患者風險預測，每次通過中位數(shù)或軟件設置單數(shù)據(jù)集打分分界點將HCC 患者分成高低風險組來驗證模型的準確性，但用這些模型進行分型難以保證不同數(shù)據(jù)集的打分分界點的統(tǒng)一性［7-12］。Liu 等［13］整合多組學數(shù)據(jù)將HCC 分成5 種類型，部分類型卻呈現(xiàn)類似的總體生存期表現(xiàn)。有些模型雖然穩(wěn)健，但需要較多的分子信息來進行模型構建，如Chaudhary 等［2］用深度學習的方法整合mRNA、miRNA 及甲基化等多組學數(shù)據(jù)，F(xiàn)a 等［14］基于多條通路或通路子集將HCC 分成兩類生存風險有差異的組別。因此，如何利用較少的分子信息來建立預后模型，既能保證較高的預測穩(wěn)健性與準確性，又具備操作簡單、診斷成本低廉的特點，是本研究的主要目標。

本研究首先針對HCC 預后相關的多條通路，通過主成分分析篩選出41 個核心基因。然后通過k 均值聚類、支持向量機等機器學習的方法在TCGA 數(shù)據(jù)集上建立HCC 患者分組預后模型，并在TCGA 數(shù)據(jù)集內(nèi)部及3 個外部數(shù)據(jù)集上進行穩(wěn)健性、有效性驗證。最后通過HCC 患者組間基因差異表達分析和富集分析，篩選出與HCC 疾病進展有關的基因與通路，以期為個性化醫(yī)療提供理論依據(jù)，同時進一步驗證模型分類的有效性。

1 材料與方法

1.1 材料

肝細胞癌患者基因表達和臨床數(shù)據(jù)來自4 個數(shù) 據(jù) 集TCGA、LIRI-JP、GSE14520 及GSE54236。其中，TCGA 數(shù)據(jù)集通過TCGA-Assembler 2.0.6 軟件從癌癥和腫瘤基因圖譜（The Cancer Genome Atlas，TCGA）數(shù) 據(jù)庫下載（https：//tcga-data.nci.nih.gov/tcga/）［15］；LIRI-JP 數(shù)據(jù)集來自國際腫瘤基因組協(xié)作組（International Caner Genome Consortium，ICGC） 數(shù) 據(jù) 庫（https：//dcc.icgc.org/projects/LIRIJP）［16］；GSE14520 和GSE54236 兩個 數(shù) 據(jù)集 來 自NCBI 的GEO 數(shù) 據(jù) 庫（https：//www.ncbi.nlm.nih.gov/geo/）［17-18］。去除生存時間不正確或缺失及非HCC 患者數(shù)據(jù)，最終得到不同樣本量的HCC 患者數(shù)據(jù)集如表1 所示。其中，數(shù)據(jù)集GSE54236 缺乏生存時間刪失（Censoring）信息。

表1 HCC 數(shù)據(jù)集介紹

1.2 方法

1.2.1 基于信號通路的基因篩選及HCC 患者預后分組模型的建立 根據(jù)文獻中提到的與肝細胞癌生存風險顯著相關的13 個信號通路或通路子集［14］，我們分別對TCGA、LIRI-JP 及GSE14520 三個HCC 數(shù)據(jù)集相應信號通路中的基因表達數(shù)據(jù)進行主成分分析（Principle component analysis，PCA），分別截取主成分載荷（Loading）大于0.3 及總方差的累積解釋比例（Cumulative proportion）大于0.8 的主成分后取交集，得到通路中的關鍵基因。PCA 分析使用R語言內(nèi)置的prcomp 函數(shù)。

基于上一步得到的與生存相關的關鍵基因，首先對整個TCGA數(shù)據(jù)集做無監(jiān)督K均值聚類（K-means clustering，k=2），將HCC 患者分為兩組。用logrank 檢驗（R 包 survival［19］）比較兩組間生存風險差異，將高風險組記為S1（Subgroup 1），低風險組記為S2（Subgroup 2）。然后使用支持向量機（Support vector machine，SVM）進行有監(jiān)督建模。模型構建與驗證采用R 語言內(nèi)置的kmeans 函數(shù)和R 包e1071中svm 函數(shù)（https：//cran.r-project.org/web/packages/e1071/index.html），其中SVM 模型的最優(yōu)參數(shù)使用交叉驗證的方法選擇。

為了驗證該模型的穩(wěn)定性，我們用整個TCGA數(shù)據(jù)集進行交叉驗證，即將數(shù)據(jù)集隨機分成5 個大小一致的組，任意選擇3 組做訓練集，其余2 組做測試集，共得到10 個數(shù)據(jù)組合進行驗證。為了評價模型的預測能力，我們將整個TCGA 數(shù)據(jù)集作為訓練集，驗證模型在LIRI-JP、GSE14520 及GSE54236數(shù)據(jù)集中的預測能力。

對組別間生存風險差異評價使用以下3 種衡量指標：一致性指數(shù)（Concordance index，C-index，R包survcomp［20］）、log-rank 檢 驗 的P 值（Log-rank P value of Cox-PH model）及Brier 分數(shù)（Brier score，R包survcomp）。一致性指數(shù)由Harrell 提出［21］，該值范圍在0 到1 之間，越大表示模型預測結(jié)果與實際生存風險排序越一致?？偵嫫冢∣verall survival，OS）的KM（Kaplan-Meier）曲 線 繪 制 使 用R 包survminer 中的ggsurvplot 函數(shù)（https：//cran.r-project.org/web/packages/survminer/index.html），并用log-rank檢驗比較組間生存差異。Brier 分數(shù)衡量預測結(jié)果與真實結(jié)果間的差異，該值越小表示模型預測能力 越好。

在建模之前，對訓練集及測試集進行了兩步歸一化處理。首先，對所有數(shù)據(jù)集分別進行中位數(shù)標準化，即利用該數(shù)據(jù)集所有腫瘤樣本對應基因表達值的中位數(shù)（Median）和絕對中位差（Median absolute deviation，MAD）對該基因表達值進行標準化。其次，利用TCGA 訓練集的各基因歸一化后的均值及方差對測試集數(shù)據(jù)進一步標準化。需要指出的，無論是測序數(shù)據(jù)還是芯片數(shù)據(jù)，我們均用各數(shù)據(jù)內(nèi)最小值（近似儀器檢測下限）來填補缺失值，未log2 處理的數(shù)據(jù)進行對數(shù)化。針對芯片數(shù)據(jù)多個探針對應一個基因的情況，我們?nèi)∑淦骄怠?/p>

1.2.2 基于單因素Cox 回歸的基因篩選及HCC 患者預后分組模型的建立 為了與上述基于信號通路中關鍵基因的模型做對比，利用單因素Cox 回歸進行基因篩選。取TCGA、LIRI-JP 及GSE14520 3 個HCC 數(shù)據(jù)集的共有基因，然后在TCGA 數(shù)據(jù)集上，利用單因素Cox 回歸篩選上述基因中與生存風險最顯著相關的子集（Log-rank P 值最?。?。為保證可比性，篩選的基因數(shù)目、HCC 患者預后分組模型的建立方法及穩(wěn)定性、預測能力的檢驗與1.2.1 相同。

1.2.3 HCC 疾病進展相關的差異基因篩選 應用

1.2.1 中構建的風險預測模型可分別將3 個主要數(shù)據(jù)集TCGA、LIRI-JP 及GSE14520 中的HCC 患者分成兩個生存風險存在顯著差異的組別。我們對兩組間的腫瘤樣本（HCC 高風險與低風險組相比）進行差異表達基因（Differentially expressed genes，DEGs）篩選，其中測序數(shù)據(jù)、芯片數(shù)據(jù)分別使用R 包DESeq2 和limma 進行分析［22-23］。差異檢驗使用BH（Benjamini-Hochberg）方法控制假陽性率［24］。對3 個數(shù)據(jù)集中滿足|log2FoldChange|＞1、校正后P 值＜0.05 的差異表達基因取交集進行后續(xù)分析。

1.2.4 差異表達基因的功能富集及注釋 使用R 包clusterProfiler 對上一步得到的與HCC 疾病進展有關的差異基因進行基因本體（Gene ontology，GO）富集分析和京都基因與基因組百科全書（Kyoto encyclopedia of genes and genomes，KEGG）信號通路分析，探索與HCC 疾病進展有關的生物學過程［25］。篩選條件均為BH 校正后的P 值＜0.01，q 值＜0.05。

1.2.5 差異表達基因的蛋白質(zhì)互相作用網(wǎng)絡及子網(wǎng)絡構建 使用STRING 11.0 數(shù)據(jù)庫（https：//stringdb.org）分析與HCC 疾病進展有關的差異表達基因的蛋白質(zhì)相互作用網(wǎng)絡（Protein-protein interaction network，PPI），并用Cytoscape 軟件進行可視化處理［26-27］。其中蛋白質(zhì)節(jié)點的度（Degree）用插件CentiScaPe 計算，圖示節(jié)點大小與度成正比，并用插件MCODE（Molecular complex detection）分析子網(wǎng)絡［28-29］。對度排名前20 的蛋白質(zhì)，與KEGG 通路過程（P 值＜0.05）取交集，得到hub 基因（Hub genes）。使 用 插 件ClueGo+CluePedia 分 析hub 基因參與的重要KEGG 通路過程［30-31］。

1.2.6hub基因的預后相關性分析 將TCGA 數(shù)據(jù)集中355 例HCC 患者的總體生存時間與hub基因表達特征進行統(tǒng)計檢驗，探索這些hub基因?qū)Ω渭毎┑念A后價值。根據(jù)每個hub基因表達水平的中位數(shù)將HCC 患者分為高表達組（表達值＞中位數(shù)）和低表達組（表達值≤中位數(shù)），通過KM 曲線圖、logrank 檢驗來研究兩組之間的總生存期差異。

2 結(jié)果

2.1 基于信號通路的基因篩選結(jié)果及HCC患者預后分組模型的表現(xiàn)

根據(jù)文獻中提到的與HCC 生存風險顯著相關的13 條信號通路或通路子集，我們分別對TCGA、LIRI-JP 及GSE14520 三個HCC 數(shù)據(jù)集中相應信號通路中的基因表達數(shù)據(jù)進行主成分分析，取交集后得到41 個核心基因：CDK1、DUT、EGF、ENO1、FOXM1、G6PD、GMDS、GNAS、GNG4、GPI、GPX7、HEXA、IRAK1、ITPR1、PPP2R5A、MAPK1、MAPK9、MAPK13、MAP2K2、RAC1、RAD51、RAF1、RENBP、RRM2、TGFA、TK1、TKT、TYMS、UGDH、UCK2、IKBKG、SQSTM1、WASF1、GNPDA1、GNB5、RALBP1、PLCB1、CYB5R1、HKDC1、NPL和UAP1L1。

根據(jù)信號通路篩選的41 個基因及無監(jiān)督聚類得到的標簽對TCGA 基因表達數(shù)據(jù)建立HCC 風險預測模型，其中穩(wěn)健性及準確性分析結(jié)果如表2 所示。TCGA 數(shù)據(jù)集交叉驗證結(jié)果表明，模型的預測準確率在0.71 附近，HCC 患者組間總體生存期的logrank 檢驗P 值平均為0.03，且Brier 分數(shù)較低。在其他3 個數(shù)據(jù)集上，該模型區(qū)分的HCC 患者組間總體生存期均存在顯著差異（log-rank 檢驗P 值＜0.05），預測準確率均高于0.7（一致性指數(shù)＞0.7）。4 個數(shù)據(jù)集的KM 曲線（圖1）可進一步直觀表明，該模型可將HCC 患者分成生存風險存在顯著差異的兩組（S1 高風險組和S2 低風險組）。其中，TCGA 數(shù)據(jù)集log-rank 檢驗P 值為0.000 18，LIRI-JP 數(shù)據(jù)集log-rank 檢驗P 值為0.000 38，GSE14520 數(shù)據(jù)集logrank 檢驗P 值為0.002 1，GSE54236 數(shù)據(jù)集log-rank 檢驗P 值為0.012。

2.2 基于單因素Cox回歸的基因篩選結(jié)果及HCC患者預后分組模型的表現(xiàn)

針對TCGA、LIRI-JP 及GSE14520 三個HCC 數(shù)據(jù)集的共有基因，利用單因素Cox 回歸得到TCGA數(shù)據(jù)集中與生存風險最顯著相關的41 個基因：CDCA8、SFPQ、KIF20A、G6PD、PSRC1、UTP11、TTC26、VRK2、KIF2C、CENPA、NEIL3、CCNJL、UCK2、TTK、ZNF643、NCAPG、YBX1、NDC80、PDE6A、TXNRD1、TPX2、SEPHS1、HDAC2、DYNC1LI1、DBF4、HJURP、DLGAP5、CHORDC1、C19orf26、CEP85、MCM10、CPSF6、EZH2、CBX2、GLMN、NCAPD2、GNL2、ZNF239、NUP205、KIF18A及TRIP13。

根據(jù)單因素Cox 回歸篩選的41 個基因及無監(jiān)督聚類得到的標簽對TCGA 基因表達數(shù)據(jù)再次建立HCC 風險預測模型，其中穩(wěn)健性及準確性分析結(jié)果如表2 所示。TCGA 數(shù)據(jù)集交叉驗證結(jié)果表明，該模型的預測準確率也在0.71 附近，HCC 患者組間總體生存期的log-rank 檢驗P 值平均為0.02，且Brier分數(shù)較低。在其他3 個數(shù)據(jù)集上，該模型區(qū)分的HCC 患者組間總體生存期均存在顯著差異（log-rank檢驗P 值＜0.05），預測準確率均高于0.65（一致性指數(shù)＞0.65）。

表2 兩個模型在TCGA 數(shù)據(jù)集做交叉驗證及外部數(shù)據(jù)上的預測表現(xiàn)

2.3 HCC疾病進展相關的基因差異表達分析與富集分析的結(jié)果

基于信號通路中41 個核心基因的模型可將HCC 患者分成高風險組S1 和低風險組S2，S1 組患者的總體生存時間顯著低于S2 組。通過對TCGA、LIRI-JP 及GSE14520 三個HCC 數(shù)據(jù)集各自S1 與S2組之間進行基因差異表達分析，最終篩選出61 個上調(diào)基因及122 個下調(diào)基因。3 個數(shù)據(jù)集所得的差異表達基因數(shù)目具體如圖2-A 所示。

對上述得到的HCC 疾病進展相關的差異表達基因進行GO 富集分析（前10 個顯著結(jié)果見圖2-B），發(fā)現(xiàn)上調(diào)基因主要富集在染色體分離、細胞核分裂、細胞器分裂、細胞核有絲分裂等生物學過程，下調(diào)基因主要富集在小分子分解代謝過程、類固醇代謝過程、羧酸生物合成過程等生物學代謝過程。對差異表達基因進一步進行KEGG 通路富集分析（圖2-C），發(fā)現(xiàn)上調(diào)基因主要涉及細胞周期、卵母細胞減數(shù)分裂、p53 信號通路、DNA 復制等信號通路，下調(diào)基因主要涉及視黃醇代謝、細胞色素P450 蛋白負責的藥物代謝等代謝通路。

圖1 四個數(shù)據(jù)集的兩組間總體生存期差異

2.4 差異表達基因的蛋白質(zhì)互作網(wǎng)絡及hub基因 篩選

利用STRING 11.0 數(shù)據(jù)庫與Cytoscape 3.7.2 軟件，構建如圖3 所示的162 個節(jié)點、1 172 條邊的差異表達基因的蛋白質(zhì)互作網(wǎng)絡。通過插件MCODE分析后發(fā)現(xiàn)兩個主要的子網(wǎng)絡模塊1（35 個節(jié)點，557 條邊，Score=34.353）和模塊2（28 個節(jié)點，174條邊，Score=12.889），圖示環(huán)狀。將編碼度排名前20 的結(jié)點蛋白質(zhì)的差異基因與KEGG 通路富集過程取交集，得到12 個hub基因：CCNB1、CDK1、

RRM2、BUB1B、CCNB2、TTK、CDC20、MCM4、RFC4、PTTG1、MCM2和MAD2L1。這12 個基因均屬于模塊1 的相互作用網(wǎng)絡，包含66 條邊的互作關系（圖4-A），經(jīng)過KEGG 富集分析發(fā)現(xiàn)主要集中在細胞周期、卵母細胞減數(shù)分裂等信號通路（圖4-B）。

2.5 hub基因的預后分析結(jié)果

根據(jù)TCGA 數(shù)據(jù)集中的355 例HCC 患者的生存數(shù)據(jù)，對12 個hub基因分別進行預后差異分析。如圖5 所示，除了基因MCM4，每個基因高表達組的生存風險率均顯著高于低表達組（Log-rank P 值＜0.05）。

3 討論

肝細胞癌是全球常見的惡性腫瘤之一，呈現(xiàn)出發(fā)病機制復雜、風險因素多樣的特點［4］。因此，融入分子水平信息對患者進行準確分類十分關鍵，能夠指導不同分型患者的個性化治療方案設計。

圖2 差異表達基因（S1 和S2 組相比）的篩選及功能注釋

圖3 差異表達基因的蛋白質(zhì)相互作用網(wǎng)絡分析（組S1 和S2 相比）

圖4 hub 基因的分析研究

本研究利用信號通路中41 個核心基因的表達數(shù)據(jù)，通過k 均值聚類、支持向量機等方法在TCGA數(shù)據(jù)集上建立了HCC 分類模型。該模型能夠?qū)CC患者分成生存風險存在顯著差異的兩組，交叉驗證結(jié)果表明模型具有較好的穩(wěn)健性，并且3 個外部數(shù)據(jù)集的測試結(jié)果證明該模型能夠?qū)CC 患者進行準確分類?；谛盘柾分?1 個核心基因的模型與基于單因素Cox 回歸篩選41 個基因的模型在TCGA 數(shù)據(jù)集上交叉驗證和LIRI-JP 數(shù)據(jù)集上準確性驗證的結(jié)果差別不大，甚至從某種程度上來說（log-rank 檢驗P 值），基于單因素Cox 回歸的模型表現(xiàn)更好。但是從GSE14520 和GSE54236 兩個數(shù)據(jù)集上的表現(xiàn)來看，基于信號通路核心基因的模型預測準確性（一致性指數(shù)＞0.7）遠好于基于單因素Cox 回歸的模型（一致性指數(shù)＞0.65）。綜上所述，這兩個模型均能對患者的生存風險進行有效區(qū)分，但是基于信號通路中核心基因的模型在多個數(shù)據(jù)集上的整體表現(xiàn)更好，更穩(wěn)定并準確預測HCC 患者的生存風險。本研究的工作也是圍繞信號通路中核心基因的模型而展開。與基于深度學習方法進行HCC 預后預測的文獻結(jié)果相比，模型在3 個相似的數(shù)據(jù)集上有不遜色的表現(xiàn)，而且需要的基因數(shù)更少，無需甲基化、miRNA 的數(shù)據(jù)，能夠進一步降低患者檢測成本［2］。

圖5 12 個hub 基因在肝細胞癌中的預后價值（總體生存期）

對分類后生存風險差異顯著的兩組HCC 患者的表達數(shù)據(jù)進行生物信息學分析，發(fā)現(xiàn)高風險組的上調(diào)基因主要富集在染色體分離、細胞核分裂、細胞器分裂、細胞核有絲分裂等生物學過程（GO 分析）及細胞周期、卵母細胞減數(shù)分裂、p53 信號通路、DNA 復制等信號通路（KEGG 分析）。其中染色體分離、細胞核分裂、細胞器分裂、細胞核有絲分裂、卵母細胞減數(shù)分裂及DNA 復制等都是發(fā)生在細胞周期中的生物學過程，p53 信號通路最終也誘導細胞周期的調(diào)控，而細胞周期調(diào)控和癌癥進展關系緊密，細胞周期失調(diào)在促進肝細胞癌發(fā)生中起核心 作用［32-33］。

本研究共發(fā)現(xiàn)12 個hub基因，除基因MCM4外，每個基因高表達組的生存風險率均顯著高于低表達組，而它們在S1 組中的相對高表達或許部分解釋了S1 組與S2 組HCC 患者間的生存風險差異。這些基 因CCNB1［4，34-37］、CCNB2［35-37］、CDK1［34-37］、R R M 2［35-37］、B U B 1 B［34，36-37］、T T K［35，37］、CDC20［4，34，36-37］、MCM2［36］、MCM4［36］、RFC4［36］、PTTG1［35-36］和MAD2L1［4，34，36-37］均 在HCC 癌 組織中表達顯著上調(diào)，該結(jié)果進一步表明分類模型的有效性，并說明這些基因可能在HCC 的發(fā)生、進展中起關鍵性作用，可以作為疾病預警的重要標志物或者潛在治療靶點。要說明的是，MCM2 和MCM4基因在HCC 預后中的作用與有些研究不一致［38-39］，可能與數(shù)據(jù)集的構成有關，需要進一步驗證。

4 結(jié)論

通過在TCGA 數(shù)據(jù)集上對信號通路中41 個核心基因建模，可以對HCC 患者的預后風險進行準確分類。交叉驗證表明該模型具備較高的穩(wěn)健性，且在另外3 個數(shù)據(jù)集上得到準確的預測結(jié)果。將生存風險差異顯著的兩組進行生物信息學分析，發(fā)現(xiàn)HCC高風險組上調(diào)基因富集在細胞周期信號通路中，并發(fā)現(xiàn)了11 個hub基因的上調(diào)與HCC 風險顯著相關，生物信息學分析也進一步驗證模型的有效性。