基于TCGA和ICGC構(gòu)建和驗(yàn)證肝細(xì)胞癌腫瘤微環(huán)境相關(guān)基因的預(yù)后風(fēng)險(xiǎn)模型

李紹智 梁維仁 郭金友 鄭家平

肝細(xì)胞癌（hepatocellular carcinoma，HCC）是全球最常見的惡性腫瘤之一，其發(fā)病率占惡性腫瘤的第六位，病死率占第四位［1］。近年來有較多治療方式，包括手術(shù)治療和藥物治療等已經(jīng)取得巨大進(jìn)展，但HCC患者的預(yù)后仍不理想，其中主要原因是腫瘤的復(fù)發(fā)和轉(zhuǎn)移［2］。越來越多的研究表明腫瘤微環(huán)境（tumor microenvironment，TME）在腫瘤的發(fā)生、發(fā)展等過程中扮演著重要的角色［3-4］。腫瘤細(xì)胞可以通過分泌各種細(xì)胞因子、趨化因子等，改變TME并導(dǎo)致周圍細(xì)胞重新編程，使這些細(xì)胞能夠促進(jìn)或抑制腫瘤細(xì)胞的生存和發(fā)展［3，5］。研究報(bào)道TME相關(guān)基因在預(yù)測不同癌種患者的生存、預(yù)后具有重要價(jià)值［6-8］。在肝癌上，盡管一些研究已經(jīng)構(gòu)建了TME相關(guān)基因的預(yù)后風(fēng)險(xiǎn)模型，但這些研究缺乏獨(dú)立數(shù)據(jù)集的驗(yàn)證，并且預(yù)測的性能仍需提高［9-11］。生物信息學(xué)方法在探索和鑒定腫瘤預(yù)后相關(guān)生物標(biāo)志物上具有重要的價(jià)值［12］。本研究基于TCGA和ICGC數(shù)據(jù)庫篩選TME相關(guān)的關(guān)鍵基因，構(gòu)建和驗(yàn)證預(yù)后風(fēng)險(xiǎn)評估模型，以期為HCC患者的預(yù)后評估提供參考。

1 材料與方法

1.1 數(shù)據(jù)下載和預(yù)處理 HCC轉(zhuǎn)錄組數(shù)據(jù)、臨床信息、隨訪資料和體細(xì)胞突變信息從癌癥基因圖譜（the cancer genome atlas，TCGA）數(shù)據(jù)庫（https：//portal.gdc.cancer.gov/）中獲取，其中對隨訪信息有缺失的患者予以刪除。其次，獨(dú)立數(shù)據(jù)集國際腫瘤基因組協(xié)作組（international cancer genome consortium，ICGC）中的HCC隊(duì)列（LIRIJP）患者的轉(zhuǎn)錄組數(shù)據(jù)也被下載用于驗(yàn)證。根據(jù)先前已經(jīng)發(fā)表的TME相關(guān)研究［13-19］，4,061個(gè)TME相關(guān)的基因被選擇用于后續(xù)分析。

1.2 非負(fù)矩陣分解算法驗(yàn)證分子亞型 通過單因素Cox分析對TCGA中這4,061個(gè)TME基因表達(dá)水平進(jìn)行分析，獲得與預(yù)后相關(guān)的TME基因（P＜0.01）。通過非負(fù)矩陣分解（non-negative matrix factorization，NMF）對樣本進(jìn)行聚類［20］，聚類的數(shù)量k根據(jù)cophenetic和殘差平方和確定。最終選擇合適的聚類數(shù)量用于后續(xù)分析。為了比較不同聚類間微環(huán)境細(xì)胞群的差異，該方法［19］被用于計(jì)算每個(gè)樣本中各種免疫細(xì)胞的評分，包括CD8+T細(xì)胞、細(xì)胞毒性淋巴細(xì)胞、B細(xì)胞、中性粒細(xì)胞、單核細(xì)胞、骨髓樹突狀細(xì)胞等。不同免疫細(xì)胞在不同聚類之間的差異也被比較。

1.3 預(yù)測模型構(gòu)建 對臨床數(shù)據(jù)整理過濾后，370例HCC患者（TCGA數(shù)據(jù)集）以3∶2的比例被隨機(jī)分為訓(xùn)練組（Training cohort）和測試組（Testing cohort），并比較兩組間臨床信息以確定是否存在組內(nèi)差異，進(jìn)一步通過Lasso、Cox回歸尋找合適的預(yù)測相關(guān)TME基因［21］。Lasso回歸能夠消除過于相關(guān)的基因，防止過度擬合，并能找到誤差最小的基因。基于Lasso回歸的結(jié)果，多因素Cox回歸分析被用于進(jìn)一步篩選重要的與HCC患者預(yù)后密切相關(guān)的TME樞紐基因。根據(jù)TME相關(guān)預(yù)后基因的表達(dá)及其多因素分析產(chǎn)生的回歸系數(shù)生成最終的預(yù)后特征：TME風(fēng)險(xiǎn)評分（TMEScore）=（βmRNA1×mRNA1的表達(dá)水平）+（βmRNA2×mRNA2的表達(dá)水平）+…+（βmRNAn×mRNAn的表達(dá)水平）。TMEScore的cutoff值是根據(jù)訓(xùn)練集中的中位TMEScore決定。根據(jù)TMEScore的cutoff值，將患者分為高風(fēng)險(xiǎn)組或低風(fēng)險(xiǎn)組。

1.4 預(yù)測模型性能評估 基于建立的TMEScore，進(jìn)一步比較訓(xùn)練集、測試集、整個(gè)人群以及獨(dú)立數(shù)據(jù)集（ICGC-JP）中人群的預(yù)后差異以及預(yù)測模型對于1、3、5年預(yù)后評估的預(yù)測性能。單因素、多因素Cox回歸分析用于研究TMEScore在整個(gè)人群中的預(yù)后價(jià)值，其中在單因素分析模型中具有顯著差異的參數(shù)被納入多因素分析模型中進(jìn)一步驗(yàn)證。亞組分析用來明確TMEScore在不同的亞組人群中的預(yù)后價(jià)值。值得一提的是，該模型的預(yù)測性能也與已發(fā)表的預(yù)后模型相比較［9-11］。

1.5 列線圖構(gòu)建和評估 列線圖是一種通過標(biāo)尺的長度來衡量不同變量或者特征對結(jié)果的影響，是一種直觀的風(fēng)險(xiǎn)模型可視化方法。因此，應(yīng)用TMEScore整合臨床變量構(gòu)建預(yù)后預(yù)測的列線圖用于評估1年、3年和5年的總體生存率，主要臨床特征包括年齡、性別、分級、分期和TMEScore。校準(zhǔn)曲線用于檢驗(yàn)列線圖的擬合的準(zhǔn)確性。

1.6 富集通路分析 為了明確高、低風(fēng)險(xiǎn)評分組的富集通路差異，從MSigDB數(shù)據(jù)庫中下載并獲得的“c2.cp.kegg.v7.4.symbols”基因集用于GSEA分析，并比較兩個(gè)風(fēng)險(xiǎn)組之間的富集通路評分差異，并用熱圖可視化。使用R包 “l(fā)imma”“org.Hs.eg.db”“clusterProfiler”和“enrichplot”來確定高、低風(fēng)險(xiǎn)評分組之間的富集通路差異。

1.7 統(tǒng)計(jì)學(xué)方法 采用R（3.6.1版）和Strawberry Perl（5.32.0.1版）軟件。計(jì)量資料用t檢驗(yàn)，計(jì)數(shù)資料用卡方檢驗(yàn)或Fisher確切概率法。預(yù)后分析采用Kaplan-Meier曲線。P＜0.05為差異有統(tǒng)計(jì)學(xué)意義。

2 結(jié)果

2.1 鑒定HCC不同的分子亞型 在數(shù)據(jù)預(yù)處理后，370例和240例HCC樣本的轉(zhuǎn)錄組芯片分別從TCGA和ICGC數(shù)據(jù)庫中獲得。通過單因素Cox分析對TCGA隊(duì)列中的TME相關(guān)基因表達(dá)分析，共有528/4,061個(gè)基因被篩選為TME相關(guān)的預(yù)后基因。通過NMF分析將所有患者分成C1和C2兩個(gè)亞群（見圖1A），結(jié)果發(fā)現(xiàn)患者的OS和PFS在兩個(gè)亞群中有明顯的差異，C1亞群患者的OS和PFS高于C2亞群的患者（見圖1B和1C）。進(jìn)一步比較兩個(gè)亞群之間的免疫細(xì)胞浸潤評分，發(fā)現(xiàn)CD8+T細(xì)胞、上皮細(xì)胞、單核細(xì)胞、骨髓樹突狀細(xì)胞、NK細(xì)胞和T細(xì)胞在兩個(gè)亞群之間有明顯的差異（見圖1D）。聚類分析也發(fā)現(xiàn)C1和C2兩個(gè)亞群分別匯聚與不同的免疫亞型（見圖1E）。

圖1 兩個(gè)聚類（C1和C2）生存、預(yù)后及免疫評分和亞型的比較。A. 通過非負(fù)矩陣分解算法進(jìn)行聚類。B～C. 兩個(gè)聚類之間的總生存期（Overall survival，OS）和無進(jìn)展生存期（Progression-free survival，PFS）。D. 腫瘤微環(huán)境中的不同免疫細(xì)胞評分在兩個(gè)聚類間的比較。E. 免疫亞型在兩個(gè)聚類間的比例

2.2 基于TME相關(guān)基因構(gòu)建預(yù)后風(fēng)險(xiǎn)評分 為了構(gòu)建TME相關(guān)的預(yù)后風(fēng)險(xiǎn)評分模型，隨機(jī)將TCGA隊(duì)列分成訓(xùn)練集（224例）和測試集（146例）兩組，兩組間的臨床基線資料差異無統(tǒng)計(jì)學(xué)意義（見表1）。為了驗(yàn)證與預(yù)后相關(guān)的TME基因，對訓(xùn)練集的224例患者行單因素Cox分析，共得到488個(gè)預(yù)后相關(guān)的TME基因。為了進(jìn)一步縮小TME基因的范圍，Lasso回歸和Akaike信息準(zhǔn)則也被用于篩選TME相關(guān)基因。最終，2個(gè)關(guān)鍵的TME預(yù)后基因（CDCA8和CBX2）被確認(rèn)，TMEScore的風(fēng)險(xiǎn)評分模型為CDCA8的表達(dá)水平×0.3533+CBX2的表達(dá)水平×0.5156。不同TMEScore對患者的預(yù)后影響以及該評分模型對預(yù)后的診斷性能被進(jìn)一步評估。單因素和多因素分析均揭示TMEScore是預(yù)后相關(guān)的風(fēng)險(xiǎn)因素（見表2）。生存分析顯示基于TMEScore的高、低風(fēng)險(xiǎn)分組能夠很好區(qū)分訓(xùn)練集（見圖2A）、測試集（見圖2B）、整個(gè)數(shù)據(jù)集（見圖2C）以及ICGC（見圖2D）中的人群預(yù)后，且高TMEScore患者比低TMEScore患者的預(yù)后更差。并且，TMEScore預(yù)測1、3和5年患者生存的AUC在訓(xùn)練集中為0.818、0.721和0.699（見圖2E）；在測試集中分別為0.671、0.652和0.633（見圖2F）；在整個(gè)TCGA中分別為0.762、0.704和0.664（見圖2G）。在ICGC數(shù)據(jù)集中，TMEScore預(yù)測1、3和5年患者生存的AUC為0.770、0.783和0.390（見圖2H）。該預(yù)測評分在ICGC數(shù)據(jù)集中的5年AUC較低的主要是由于其中位隨訪時(shí)間相對較短，較少的事件數(shù)影響了模型預(yù)測的效果。為了進(jìn)一步分析TMEScore在不同臨床特征的人群中的預(yù)測價(jià)值，亞組分析在TCGA的患者中開展。分層分析結(jié)果顯示TMEScore在年齡＞65歲、≤65歲、男性和女性、G1～2期、G3～4期、M0、N0、I～I(xiàn)I期、III～I(xiàn)V期、T1～2期、T3～4期人群中均顯示出預(yù)后預(yù)測價(jià)值（見圖2I-2T）。基于TMEScore風(fēng)險(xiǎn)評分在T1期患者及其他T分期患者中差異有統(tǒng)計(jì)學(xué)意義，提示該預(yù)后評分對區(qū)分早期和晚期患者具有一定的價(jià)值（見圖2U）。

表1 訓(xùn)練集和測試集患者的臨床基線特征［n（%）］

表2 單因素和多因素分析預(yù)后相關(guān)風(fēng)險(xiǎn)特征

圖2 基于TMEScore的風(fēng)險(xiǎn)評分的預(yù)測效能評估

2.3 整合臨床參數(shù)的列線圖預(yù)測預(yù)后 通過整合臨床參數(shù)和TMEScore，一個(gè)預(yù)測患者1、3和5年生存評估列線圖被建立（見圖3A）。擬合曲線提示該列線圖具有良好的擬合性（見圖3B）。進(jìn)一步比較構(gòu)建的列線圖與其他臨床參數(shù)的預(yù)測性能，結(jié)果提示該列線圖預(yù)測患者1年生存的AUC比其他臨床參數(shù)高（AUC=0.777）（見圖3C）。決策曲線分析也提示構(gòu)建的列線圖的預(yù)測性能比其他臨床參數(shù)更好（見圖3D）。為進(jìn)一步驗(yàn)證構(gòu)建的TMEScore在預(yù)測患者預(yù)后上的價(jià)值，本研究搜集了已發(fā)表的與TME相關(guān)的風(fēng)險(xiǎn)模型，并比較這些模型的預(yù)測性能。相比于其他3個(gè)已發(fā)表的TME相關(guān)預(yù)測模型，本研究的模型在1、3和5年的AUC上均高于這三個(gè)模型（見圖3E）。生存分析提示Deng等、Tian等及本研究構(gòu)建的模型能夠良好區(qū)分高、低風(fēng)險(xiǎn)人群的預(yù)后（見圖3F）。以上結(jié)果也提示本研究構(gòu)建的TMEScore以及列線圖模型能夠有效評估患者生存預(yù)后風(fēng)險(xiǎn)，對具有高死亡風(fēng)險(xiǎn)的患者進(jìn)行有效的監(jiān)督隨訪和必要的干預(yù)對提升這些患者的生存具有重要意義。

2.4 富集通路分析和相關(guān)性分析 進(jìn)一步對高、低風(fēng)險(xiǎn)評分組內(nèi)的樣本行基因富集分析以明確高、低風(fēng)險(xiǎn)評分組之間涉及的信號通路差異，結(jié)果提示高風(fēng)險(xiǎn)評分組與KEGG_CELL_CYCLE，KEGG_DNA_REPLICATION，KEGG_HOMOLOGOUS_RECOMBINATION，KEGG_OOCYTE_MEIOSIS和KEGG_PRIMARY_IMMUNODEFICIENCY通路相關(guān)（見圖4A）。低風(fēng)險(xiǎn)評分組與KEGG_COMPLEMENT_AND_COAGULATION_CASCADES，KEGG_DRUG_METABOLISM_CYTOCHROME_P450，KEGG_FATTY_ACID_METABOLISM，KEGG_PPAR_SIGNALING_PATHWAY和KEGG_RETINOL_METABOLISM通路相關(guān)（見圖4B）。進(jìn)一步研究發(fā)現(xiàn)TMEScore與免疫檢查點(diǎn)（PDCD1、CD274、CTLA4）、DNA復(fù) 制（POLE2、FEN1、MCM6、POLD3）、錯(cuò)配修復(fù)（MSH6、MSH2）和上皮-間質(zhì)轉(zhuǎn)化（FAP、LOXL2）均呈正相關(guān)（見圖4C）。由于TMEScore與免疫檢查點(diǎn)成正相關(guān)，進(jìn)一步分析其與免疫細(xì)胞之間的相關(guān)性。結(jié)果顯示TMEScore與T細(xì)胞、CD8+T細(xì)胞、細(xì)胞毒性淋巴細(xì)胞、單核細(xì)胞、髓系樹突狀細(xì)胞、纖維母細(xì)胞呈正相關(guān)，與中性粒細(xì)胞呈負(fù)相關(guān)（見圖4D）。

圖3 列線圖構(gòu)建可視化風(fēng)險(xiǎn)評分模型。A. 構(gòu)建整合TMEScore的風(fēng)險(xiǎn)值及臨床參數(shù)的列線圖。B. 患者1年、3年、5年的擬合曲線。C. 不同臨床參數(shù)之間預(yù)測患者1年生存的AUC比較。D. 決策曲線評估列線圖的臨床獲益。E. 基于TMEScore的風(fēng)險(xiǎn)評分模型預(yù)測性能與先前發(fā)表的TME相關(guān)預(yù)測模型的比較。F. 基于TMEScore的風(fēng)險(xiǎn)評分模型與先前已發(fā)表的TME相關(guān)預(yù)測模型對患者的預(yù)后評估比較

圖4 基于TMEScore的風(fēng)險(xiǎn)評分值在樣本中的分子機(jī)制。A-B. 高、低風(fēng)險(xiǎn)評分組中的樣本參與的重要通路。C. HCC樣本風(fēng)險(xiǎn)評分值與腫瘤學(xué)重要通路的靶點(diǎn)的代表性基因表達(dá)之間的相關(guān)性。D. HCC樣本風(fēng)險(xiǎn)評分值與免疫細(xì)胞評分之間的相關(guān)性

3 討論

大量文獻(xiàn)表明TME在HCC的病理生成及復(fù)發(fā)、轉(zhuǎn)移中起關(guān)鍵作用，這表明HCC的生物學(xué)行為不僅與腫瘤細(xì)胞有關(guān)，也與TME密切相關(guān)［22-23］。因此，為了明確HCC中與TME相關(guān)的預(yù)后生物標(biāo)志物，本研究根據(jù)先前發(fā)表的研究整理和過濾了4,061個(gè)TME相關(guān)的基因。通過Lasso Cox回歸結(jié)果確定了2個(gè)關(guān)鍵的TME預(yù)后基因（CDCA8和CBX2）。其中，CDCA8是染色體乘客復(fù)合體（Chromosomal passenger complex，CPC）的一個(gè)重要組成部分，與有絲分裂過程中的細(xì)胞動態(tài)定位調(diào)節(jié)有關(guān)。JOEN等［24］報(bào)道，沉默CDCA8能夠恢復(fù)ATF3腫瘤抑制因子的功能和以及失活致癌信號通路AKT/β-catenin，從而抑制HCC的生長和干性，提示靶向CDCA8可能是治療原發(fā)性HCC復(fù)發(fā)、轉(zhuǎn)移的潛在靶點(diǎn)。CUI等［25］發(fā)現(xiàn)高表達(dá)的CDCA8與HCC患者更差的總生存率（P=0.0054）和無進(jìn)展生存（P=0.0009）相關(guān)。體外實(shí)驗(yàn)顯示抑制CDCA8能減緩細(xì)胞增殖，阻斷G0/G1階段的細(xì)胞周期。體內(nèi)實(shí)驗(yàn)表明抑制CDCA8能抑制腫瘤的生長。這些發(fā)現(xiàn)也提示CDCA8可能是HCC的一個(gè)有效治療靶點(diǎn)。另一個(gè)關(guān)鍵基因CBX2，是核心蛋白復(fù)合物1（polycomb repressive complex 1，PRC1）的一個(gè)重要組成部分。MAO等［26］發(fā)現(xiàn)高表達(dá)的CBX2與HCC患者的不良預(yù)后有關(guān)，且敲除CBX2能抑制HCC細(xì)胞的增值和促進(jìn)HCC細(xì)胞的凋亡，提示CBX2是一個(gè)可能的治療靶標(biāo)。

基于這兩個(gè)關(guān)鍵的預(yù)測子，本研究構(gòu)建TMEScore并發(fā)現(xiàn)高TMEScore與患者更差的預(yù)后有關(guān)，且在獨(dú)立數(shù)據(jù)集和亞組人群中也驗(yàn)證了其預(yù)后價(jià)值，提示該風(fēng)險(xiǎn)評分特征能作為預(yù)測患者預(yù)后的獨(dú)立因素。基于TMEScore構(gòu)建的列線圖能很好評估患者的預(yù)后風(fēng)險(xiǎn)。盡管先前研究也報(bào)道TME相關(guān)基因在預(yù)后預(yù)測上的價(jià)值［9-11］，但通過比較后發(fā)現(xiàn)本研究的模型在預(yù)測性能上高于已發(fā)表的模型。這些結(jié)果也驗(yàn)證了本研究的風(fēng)險(xiǎn)模型在患者預(yù)后預(yù)測上的價(jià)值。

TMEScore與免疫檢查點(diǎn)（PDCD1、CD274、CTLA4）、DNA復(fù)制（POLE2、FEN1、MCM6、POLD3）、錯(cuò)配修復(fù)（MSH6、MSH2）和上皮-間質(zhì)轉(zhuǎn)化（FAP、LOXL2）均呈正相關(guān)。這與GSEA富集分析的結(jié)果相一致。這些重要的機(jī)制也已經(jīng)被報(bào)道與HCC的發(fā)生、發(fā)展及療效有關(guān)［27］。免疫相關(guān)分析也揭示TMEScore與多種免疫細(xì)胞相關(guān)，提示該風(fēng)險(xiǎn)評分可能與腫瘤浸潤免疫細(xì)胞（tumor-infiltrating lymphocytes，TILs）水平相關(guān)。TILs水平也被認(rèn)為是重要的免疫治療生物標(biāo)志物，能夠反映免疫治療療效［28］。因此，本研究構(gòu)建的風(fēng)險(xiǎn)評分模型可能對預(yù)測免疫治療療效也有一定的價(jià)值，但仍需要后續(xù)的研究。

本研究通過生物信息學(xué)的方法構(gòu)建和驗(yàn)證了一個(gè)TME相關(guān)的雙基因預(yù)后風(fēng)險(xiǎn)評估模型（CDCA8和CBX2）。該模型能夠有效預(yù)測HCC患者預(yù)后風(fēng)險(xiǎn)水平，對于具有高死亡風(fēng)險(xiǎn)的人群，臨床醫(yī)師能夠采取必要的監(jiān)督隨訪和合理的干預(yù)措施提升患者的生存。