肝癌細(xì)胞周期調(diào)控相關(guān)基因預(yù)后模型構(gòu)建與評(píng)估

孫東旭,朱文靜,金志朋,劉華元,朱朋程,史光軍

(1 大連醫(yī)科大學(xué)研究生院,遼寧 大連 116000； 2 青島大學(xué)附屬青島市市立醫(yī)院肝膽外科； 3 青島大學(xué)青島醫(yī)學(xué)院)

肝細(xì)胞癌是世界范圍內(nèi)發(fā)病率較高的惡性腫瘤，約占肝癌病人的90%[1]。盡管肝細(xì)胞癌的治療取得了一些進(jìn)展，但肝細(xì)胞癌病人的預(yù)后仍然很差[2]。既往生物信息學(xué)綜合性研究所構(gòu)建的肝癌預(yù)后模型等研究結(jié)果十分廣泛，包括基于免疫相關(guān)編碼基因集合[3]、p53相關(guān)的microRNA集合[4]等。但由于預(yù)后腫瘤標(biāo)志物和治療靶點(diǎn)尚未得到充分研究和臨床應(yīng)用，肝細(xì)胞癌病人的預(yù)后判斷和個(gè)體化診療仍是一大挑戰(zhàn)。本研究的目的是構(gòu)建預(yù)后模型，為肝癌病人的預(yù)后預(yù)測(cè)和個(gè)體化治療提供分子標(biāo)志物和新的方向。

1 資料和方法

1.1 肝癌轉(zhuǎn)錄表達(dá)數(shù)據(jù)的獲取和差異表達(dá)基因的篩選

從TCGA數(shù)據(jù)庫(https://portal.gdc.cancer.gov/)下載TCGA-LIHC肝癌數(shù)據(jù)集。TCGA數(shù)據(jù)庫肝癌數(shù)據(jù)集包含374例肝細(xì)胞癌腫瘤組織樣本和50例癌旁正常肝組織樣本的表達(dá)數(shù)據(jù)以及臨床數(shù)據(jù)。使用統(tǒng)計(jì)學(xué)軟件R軟件(3.6.1版)[5]和Bioconductor ‘edge’軟件包分析肝細(xì)胞癌樣本與正常組織間差異表達(dá)基因的表達(dá)差異[6-7]。|Log2FC|>2和校正后P值<0.05的基因被定義為差異表達(dá)基因。

從GEO數(shù)據(jù)庫(https://www.ncbi.nlm.nih.gov/geo/)[8]GPL10558平臺(tái)(Illumina HumanHT-12 V4.0 expression beadchip)下載肝癌數(shù)據(jù)集GSE36376。GSE36376數(shù)據(jù)集包含240例肝細(xì)胞癌組織樣本和193例癌旁組織樣本的表達(dá)數(shù)據(jù)和臨床數(shù)據(jù)。|Log2FC|>1和校正后P值<0.05的基因被鑒定為差異表達(dá)基因。使用維恩圖在線工具(http://bioinformatics.psb.ugent.be/webtools/Venn/) 繪制韋恩圖鑒定共同上調(diào)和下調(diào)基因。

從ICGC數(shù)據(jù)庫(https://dcc.icgc.org/projects/LIRI-JP/)LIRI-JP肝癌數(shù)據(jù)集下載231例肝癌樣本的表達(dá)數(shù)據(jù)和臨床數(shù)據(jù)。這些樣本主要來自日本乙型肝炎病毒(HBV)或丙型肝炎病毒(HCV)感染人群[9]。樣本數(shù)據(jù)使用了標(biāo)準(zhǔn)化的計(jì)數(shù)值。

1.2 肝癌差異表達(dá)基因的通路和功能富集分析

利用Metascape網(wǎng)站[10]對(duì)差異表達(dá)基因進(jìn)行通路和功能富集分析?；騁O功能注釋及基因參與通路來源于以下數(shù)據(jù)庫的并集:Kyoto Encyclopaedia of Genes and Genomes (KEGG) Pathway, Gene Ontology (GO) Biological Processes, Reactome Gene Sets, Canonical Pathways, CORUM, TRRUST, DisGeNET, PaGenBase, Transcription Factor Targets, COVID。將基因組中的所有基因作為富集背景。P值的計(jì)算基于累積超幾何分布，q值的計(jì)算采用Benjamin-Hochberg (BH)進(jìn)行多重檢驗(yàn)[11]。最后使用Cytoscape可視化網(wǎng)絡(luò)[12]。

1.3 肝癌細(xì)胞周期調(diào)控相關(guān)基因預(yù)后模型的構(gòu)建和驗(yàn)證

采用單因素Cox回歸分析細(xì)胞周期調(diào)控相關(guān)差異表達(dá)基因的預(yù)后價(jià)值。根據(jù)表達(dá)量的中位值將病人分為高表達(dá)組和低表達(dá)組，通過在線Kaplan-Meier plotter (http://kmplot.com/analysis/)進(jìn)行Kaplan-Meier生存曲線驗(yàn)證[13]。使用Lasso Cox回歸分析方法建立預(yù)后模型[14-15]。采用‘glmnet R’包使用LASSO算法進(jìn)行選擇和收縮自變量。根據(jù)中位風(fēng)險(xiǎn)評(píng)分將病人分為高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組?；谀Ｐ椭械幕虮磉_(dá)，采用‘stats’R包的‘prcomp’程序進(jìn)行主成分分析(PCA)；同樣基于模型中的基因表達(dá)，采用‘Rtsne’R包中的t-分布隨機(jī)相鄰嵌入分析(t-SNE)方法，分析不同風(fēng)險(xiǎn)組的分布，確定各風(fēng)險(xiǎn)組的區(qū)分顯著性。采用‘survminer’ R包的‘sur_cutpoint’程序來確定最佳截?cái)啾磉_(dá)值，進(jìn)行Kaplan-Meier生存分析確定高低風(fēng)險(xiǎn)組的病人生存情況差異。使用單因素和多因素Cox回歸分析確定模型風(fēng)險(xiǎn)評(píng)分是否為總生存期(OS)的獨(dú)立預(yù)后因素。應(yīng)用‘survival ROC’R包進(jìn)行時(shí)間依賴性受試者工作特征(ROC)曲線分析，以評(píng)估模型基因集的預(yù)測(cè)能力。生成用于模型可視化和臨床應(yīng)用的列線圖(Nomogram)，應(yīng)用校準(zhǔn)曲線(Calibration curve)評(píng)價(jià)列線圖的校準(zhǔn)度，應(yīng)用決策曲線分析(DCA)評(píng)價(jià)臨床適用度。

1.4 樣品采集和標(biāo)準(zhǔn)化處理

收集青島大學(xué)附屬青島市市立醫(yī)院肝膽外科3例確診為肝細(xì)胞癌病人的肝癌組織和癌旁組織，樣本采集和存儲(chǔ)采用標(biāo)準(zhǔn)化的方法。對(duì)組織樣本進(jìn)行基因轉(zhuǎn)錄水平二代測(cè)序(NGS)，對(duì)數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理，統(tǒng)計(jì)方法采用Mann-WhitneyU檢驗(yàn)。

1.5 統(tǒng)計(jì)學(xué)分析

所有統(tǒng)計(jì)分析均使用R軟件。除特殊標(biāo)注外，計(jì)量資料比較采用t檢驗(yàn)，計(jì)數(shù)資料比較采用χ2檢驗(yàn)。應(yīng)用Cox回歸估計(jì)危險(xiǎn)比(HR)和95%置信區(qū)間(CI)。生存分析采用Kaplan-Meier法，采用logrank檢驗(yàn)確定差異是否有統(tǒng)計(jì)學(xué)意義。使用BH法校正P值。采用雙側(cè)檢驗(yàn)，P<0.05為差異有統(tǒng)計(jì)學(xué)意義。

2 結(jié) 果

2.1 肝癌腫瘤組織和正常肝臟組織差異表達(dá)基因的篩選

TCGA數(shù)據(jù)庫TCGA-LIHC肝癌數(shù)據(jù)集共篩選出3 619個(gè)差異表達(dá)基因 (|log2FC|>2, FDR<0.05)，差異表達(dá)基因的熱圖和火山圖見圖1 A、B。 GEO數(shù)據(jù)庫GSE36376肝癌數(shù)據(jù)集共篩選出687個(gè)差異表達(dá)基因 (|log2FC|>1, FDR<0.05)。應(yīng)用韋恩圖共同鑒定了141個(gè)差異表達(dá)基因，其中70個(gè)基因表達(dá)顯著上調(diào)，71個(gè)基因表達(dá)顯著下調(diào)。見圖1 C、D和表1。

2.2 肝癌細(xì)胞周期調(diào)控相關(guān)預(yù)后基因的確定

通路及功能富集分析顯示，肝癌差異表達(dá)基因共參與了409個(gè)重要功能及通路(圖1)，其中有95個(gè)通路和功能與肝癌細(xì)胞周期調(diào)控密切相關(guān)，通過統(tǒng)計(jì)歸納，最后確定了28個(gè)與肝癌細(xì)胞周期調(diào)控相關(guān)基因。見表2。單因素Cox回歸分析顯示，與肝癌預(yù)后相關(guān)的細(xì)胞周期調(diào)控基因有24個(gè)，其中包括CDC20、AURKA、NUSAP1、HMMR、TP2A和MDK等(HR>1,FDR<0.05)(圖2A)；基因表達(dá)熱圖顯示了這些基因的表達(dá)水平(FDR<0.05)(圖2B)。應(yīng)用在線Kaplan-Meier Plotter分析驗(yàn)證肝癌病人細(xì)胞周期調(diào)控相關(guān)基因的預(yù)后價(jià)值，最終確定這24個(gè)細(xì)胞周期調(diào)控相關(guān)基因均與肝癌病人的預(yù)后顯著相關(guān)(圖2C)。

A、B：TCGA-LIHC肝癌數(shù)據(jù)集374例肝癌組織與50例癌旁正常組織差異表達(dá)基因熱圖和火山圖；其中行代表基因，列代表樣本；藍(lán)色代表低表達(dá)，紅色代表高表達(dá)。C、D：TCGA-LIHC肝癌數(shù)據(jù)集和GSE36376肝癌數(shù)據(jù)集中肝癌組織和癌旁組織差異表達(dá)基因的交集，其中70個(gè)基因表達(dá)上調(diào)，71個(gè)基因表達(dá)下調(diào)。

表1 TCGA-LIHC數(shù)據(jù)集與GSE36376數(shù)據(jù)集交集差異表達(dá)基因

A:單因素Cox回歸分析確定影響肝癌病人預(yù)后的細(xì)胞周期調(diào)控相關(guān)基因;B:影響肝癌病人預(yù)后的細(xì)胞周期調(diào)控相關(guān)基因在肝癌腫瘤組織和癌旁正常組織中表達(dá)熱圖；C:Kaplan-Meier Plotter在線分析驗(yàn)證細(xì)胞周期調(diào)控相關(guān)基因的預(yù)后價(jià)值。

2.3 肝癌細(xì)胞周期調(diào)控基因預(yù)后模型的構(gòu)建

基于TCGA數(shù)據(jù)庫TCGA-LIHC肝癌病人隊(duì)列，用Lasso Cox回歸分析建立預(yù)后模型?；趹土P參數(shù)的最優(yōu)值λ，確定了一個(gè)8個(gè)基因的基因集(圖3)。風(fēng)險(xiǎn)評(píng)分計(jì)算方法如下:風(fēng)險(xiǎn)評(píng)分=e(0.319×CDC20表達(dá)量-0.393×NUSAP1表達(dá)量+0.438×HMMR表達(dá)量+0.066×ARID3A表達(dá)量+0.068×RACGAP1表達(dá)量+0.123×NCAPG表達(dá)量-0.141×SPC24表達(dá)量+0.004×MELK表達(dá)量)。根據(jù)其中位截?cái)嘀?，將病人分為高風(fēng)險(xiǎn)組(n=182)和低風(fēng)險(xiǎn)組(n=183)(圖3A)。PCA和t-SNE分析顯示，高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組病人離散方向不同(圖3B、C)，高風(fēng)險(xiǎn)病人早期死亡的可能性高于低風(fēng)險(xiǎn)病人(圖3D)。Kaplan-Meier曲線分析顯示，高風(fēng)險(xiǎn)組的OS明顯低于低風(fēng)險(xiǎn)組(圖3E,P<0.001)，低風(fēng)險(xiǎn)評(píng)分的肝癌病人較高風(fēng)險(xiǎn)評(píng)分者有更好的預(yù)后。應(yīng)用ROC曲線評(píng)估模型的預(yù)測(cè)能力，生存時(shí)間1年的ROC曲線下面積(AUC)為0.800(95%CI=0.737～0.863)，2年為0.750(95%CI=0.687～0.813)，3年AUC為0.731(95%CI=0.659～0.804)，表明本文建立的預(yù)后模型具有良好的預(yù)后預(yù)測(cè)準(zhǔn)確度和特異度(圖3F)。利用TCGA隊(duì)列中多因素Cox回歸模型生成的系數(shù)，將風(fēng)險(xiǎn)評(píng)分與分期、分級(jí)、年齡和性別等重要的臨床變量整合在一起，以進(jìn)一步提高預(yù)后預(yù)測(cè)的準(zhǔn)確性，建立了模型可視化和臨床應(yīng)用的列線圖(圖4A)。校準(zhǔn)曲線檢測(cè)出列線圖預(yù)測(cè)與實(shí)際觀測(cè)之間的最佳預(yù)測(cè)閾值(圖4B)。最后，通過1、2和3年的DCA比較風(fēng)險(xiǎn)評(píng)分與其他臨床指標(biāo)的臨床凈效益(圖4C～E)，結(jié)果顯示，在上述閾值概率的大部分范圍內(nèi)，風(fēng)險(xiǎn)評(píng)分顯示出更大的凈收益，表明風(fēng)險(xiǎn)評(píng)分在預(yù)測(cè)肝癌病人預(yù)后方面具有較好的臨床應(yīng)用價(jià)值。

2.4 肝癌細(xì)胞周期調(diào)控基因預(yù)后模型的驗(yàn)證

為了檢驗(yàn)肝癌病人隊(duì)列模型的穩(wěn)健性，按照與TCGA數(shù)據(jù)庫TCGA-LIHC肝癌病人隊(duì)列構(gòu)建模型的相同公式，將ICGC數(shù)據(jù)庫LIRI-JP肝癌病人隊(duì)列分為高風(fēng)險(xiǎn)組(n=182)和低風(fēng)險(xiǎn)組(n=78)(圖5A)。PCA分析和t-SNE分析確定了病人在兩個(gè)亞組中離散方向的分布，見圖5B、C。與低風(fēng)險(xiǎn)組相比，高風(fēng)險(xiǎn)組病人早期死亡可能性更高(圖5D)，生存時(shí)間更短(圖5E,P<0.001)。ROC曲線分析顯示，生存時(shí)間1年的AUC為0.722(95%CI=0.584～0.861)，2年為0.739(95%CI=0.633～0.845)，3年為0.733(95%CI=0.627～0.839)，預(yù)后模型具有良好的預(yù)測(cè)準(zhǔn)確度和特異度(圖5F)。

2.5 肝癌細(xì)胞周期調(diào)控基因預(yù)后模型風(fēng)險(xiǎn)評(píng)分的獨(dú)立預(yù)后價(jià)值

單因素Cox回歸分析顯示，TCGA-LIHC肝癌病人隊(duì)列(構(gòu)建隊(duì)列)和LIRI-JP肝癌病人隊(duì)列(驗(yàn)證隊(duì)列)的風(fēng)險(xiǎn)評(píng)分與OS之間存在顯著相關(guān)性(構(gòu)建隊(duì)列：HR=3.767,95%CI=2.661～5.333,P<0.001;驗(yàn)證隊(duì)列：HR=3.752,95%CI=2.240～6.266,P<0.001)。多因素Cox回歸分析顯示，風(fēng)險(xiǎn)評(píng)分是OS的獨(dú)立預(yù)測(cè)因子(TCGA數(shù)據(jù)庫肝癌病人隊(duì)列:HR=3.436,95%CI=2.402～4.916,P<0.001;ICGC數(shù)據(jù)庫肝癌病人隊(duì)列:HR=3.264,95%CI=1.920～5.549,P<0.001)。見圖6。

2.6 肝癌細(xì)胞周期調(diào)控相關(guān)預(yù)后基因的轉(zhuǎn)錄表達(dá)水平鑒定

本文NGS結(jié)果顯示，包括CDC20、AURKA和NUSAP1等在內(nèi)的16個(gè)細(xì)胞周期調(diào)控相關(guān)預(yù)后基因在肝癌中表達(dá)顯著上調(diào)(圖7)。

3 討 論

肝癌等惡性腫瘤細(xì)胞的特點(diǎn)是無限增殖，這與細(xì)胞周期調(diào)控密切相關(guān)。盡管細(xì)胞周期調(diào)控的機(jī)制已經(jīng)成為腫瘤研究的核心領(lǐng)域，但其具體機(jī)制仍不明確，細(xì)胞周期調(diào)控的機(jī)制以及相關(guān)基因?qū)Ω伟┎∪祟A(yù)后的預(yù)測(cè)價(jià)值也尚不清楚。既往的研究結(jié)果表明，基于p53相關(guān)的microRNA集合[5]、免疫相關(guān)編碼基因集合[4]、CpG島甲基化表型(CIMP)相關(guān)基因[16]、控制胚胎發(fā)育的claudin基因家族[17]等構(gòu)建的肝癌預(yù)后模型顯示了優(yōu)秀的預(yù)測(cè)能力。與這些研究相比，本研究1、2、3年的ROC曲線及DCA曲線等結(jié)果均顯示本文構(gòu)建的預(yù)后模型具有良好的準(zhǔn)確性、特異性及臨床適用性，能夠準(zhǔn)確預(yù)測(cè)肝癌病人的預(yù)后。

A:病人風(fēng)險(xiǎn)評(píng)分分布及中位值;B:PCA確定高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組病人離散方向的分布;C:t-SNE確定高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組病人離散方向的分布;D:高風(fēng)險(xiǎn)組與低風(fēng)險(xiǎn)組病人生存時(shí)間、生存狀態(tài)散點(diǎn)圖;E:高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組病人OS的Kaplan-Meier生存曲線分析;F:ROC曲線評(píng)估模型的預(yù)后預(yù)測(cè)能力。

A:TCGA-LIHC肝癌病人隊(duì)列風(fēng)險(xiǎn)評(píng)分及臨床變量指標(biāo)1、2和3年OS的復(fù)合列線圖預(yù)測(cè);B:1、2、3年觀測(cè)OS概率和復(fù)合列線圖預(yù)測(cè)OS概率的校準(zhǔn)曲線；C～E:TCGA-LIHC肝癌病人隊(duì)列風(fēng)險(xiǎn)評(píng)分及臨床變量指標(biāo)1、2、3年的DCA。

A:病人風(fēng)險(xiǎn)評(píng)分分布及中位值;B:PCA確定高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組離散方向的分布;C:t-SNE確定高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組兩個(gè)亞組中離散方向的分布;D:高風(fēng)險(xiǎn)組與低風(fēng)險(xiǎn)組病人生存時(shí)間、生存狀態(tài)散點(diǎn)圖;E:高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組病人OS的Kaplan-Meier生存曲線分析;F:ROC曲線評(píng)估模型的預(yù)后預(yù)測(cè)能力。

A:構(gòu)建隊(duì)列風(fēng)險(xiǎn)評(píng)分單因素Cox回歸分析；B:驗(yàn)證隊(duì)列風(fēng)險(xiǎn)評(píng)分單因素Cox回歸分析;C:構(gòu)建隊(duì)列風(fēng)險(xiǎn)評(píng)分多因素Cox回歸分析;D:驗(yàn)證隊(duì)列風(fēng)險(xiǎn)評(píng)分多因素Cox回歸分析。

表2 28個(gè)肝癌細(xì)胞周期調(diào)控相關(guān)基因縮寫及英文全稱

ns：P>0.05;*P<0.05;**P<0.01;***P<0.001;****P<0.000 1。

本文構(gòu)建的預(yù)后模型中，參與模型的共有8個(gè)細(xì)胞周期調(diào)控相關(guān)基因，分別為RACGAP1、CDC20、NUSAP1、HMMR、ARID3A、NCAPG、SPC24和MELK。迄今為止的研究顯示，其中6個(gè)致癌基因CDC20[18]、NUSAP1[19]、RACGAP1[20-21]、NCAPG[22]、MELK[23]和SPC24[24]已經(jīng)在肝癌中被確定具有重要作用，但HMMR和ARID3A在肝癌中的作用尚不清楚。有生物信息學(xué)研究結(jié)果表明，HMMR可能是肝癌中較高表達(dá)的致癌基因[25]。本文研究表明，HMMR可能通過調(diào)控肝癌細(xì)胞周期影響病人的預(yù)后。此外，ARID3A基因在腫瘤中作用研究甚少，本文研究顯示ARID3A可能通過調(diào)控細(xì)胞周期影響肝癌病人的預(yù)后。為了驗(yàn)證本文篩選出的預(yù)后基因的表達(dá)水平，我們使用NGS技術(shù)檢測(cè)3例肝癌組織與癌旁組織基因表達(dá)，結(jié)果顯示16個(gè)細(xì)胞周期調(diào)控相關(guān)預(yù)后基因在肝癌中表達(dá)顯著上調(diào)，在轉(zhuǎn)錄水平上證明了細(xì)胞周期調(diào)控相關(guān)預(yù)后基因的作用。

綜上所述，本研究成功構(gòu)建了細(xì)胞周期調(diào)控相關(guān)基因的預(yù)后模型，為肝癌病人的預(yù)后及治療提供新的方向。本文測(cè)序分析為后續(xù)的模型驗(yàn)證提供了轉(zhuǎn)錄水平表達(dá)數(shù)據(jù)基礎(chǔ)，但仍需檢測(cè)更多的組織樣本進(jìn)行驗(yàn)證，并進(jìn)行更加深入的基礎(chǔ)實(shí)驗(yàn)研究。