肝細(xì)胞癌免疫相關(guān)基因和lncRNA聯(lián)合預(yù)后模型的構(gòu)建及驗(yàn)證

王威，程倩倩，周雪麗，季文斌，呂振宇，楊燕

作者單位：蚌埠醫(yī)學(xué)院第一附屬醫(yī)院腫瘤內(nèi)科，安徽 蚌埠233004

肝細(xì)胞癌（hepatocellular carcinoma，HCC）是常見消化系統(tǒng)惡性腫瘤，全球每年約有91萬(wàn)新發(fā)病例和83萬(wàn)死亡病例，在所有癌癥中居第6位及第3位［1-2］。盡管近年來(lái)HCC新藥不斷涌現(xiàn)，尤其是免疫治療或免疫聯(lián)合治療較大地改善了晚期HCC的整體預(yù)后，但有效療效預(yù)測(cè)和預(yù)后判斷生物標(biāo)志物的缺乏在一定程度上限制了人群的獲益。傳統(tǒng)臨床分期雖是指導(dǎo)HCC治療及判斷預(yù)后的重要因素，但一些臨床特征相同的病人，由于分子層面異質(zhì)性可能預(yù)后不同［3］。近年來(lái)，大量證據(jù)顯示包括HCC在內(nèi)的腫瘤免疫微環(huán)境與病人的生存關(guān)系密切［4］，許多學(xué)者試圖通過(guò)對(duì)免疫相關(guān)基因（immune-related genes，IRGs）或免疫相關(guān)lncRNA（immune-related lncRNAs，IRlncRNAs）的生物信息學(xué)分析來(lái)評(píng)估病人風(fēng)險(xiǎn)及預(yù)后，如Hong等［5］構(gòu)建了12個(gè)IRlncRNAs的HCC預(yù)后模型，劉航睿等［6］構(gòu)建了7個(gè)IRGs的食管癌預(yù)后模型，這些由IRGs或轉(zhuǎn)錄產(chǎn)物構(gòu)成的模型在病人預(yù)后評(píng)價(jià)中具有良好的應(yīng)用價(jià)值，但這些模型多只涉及一種類型的RNA，多類型RNA的聯(lián)合預(yù)測(cè)模型研究較少。考慮到HCC免疫學(xué)研究的重要性，以及癌癥發(fā)生發(fā)展中不同分子間的密切協(xié)同作用［7-8］，本研究首次構(gòu)建了IRGs及IRlncRNAs聯(lián)合的HCC預(yù)后模型，以期為HCC免疫學(xué)機(jī)制研究及預(yù)后預(yù)測(cè)提供一定參考。

1 資料與方法

1.1 數(shù)據(jù)獲取及處理 2022年6—8月通過(guò)TCGA官網(wǎng)下載HCC轉(zhuǎn)錄組及與對(duì)應(yīng)的臨床數(shù)據(jù)；https：//www.immport.org/網(wǎng)站下載IRGs；根據(jù)Ensembl基因轉(zhuǎn)移格式文件將基因表達(dá)矩陣區(qū)分為mRNA及LncRNA表達(dá)陣列，Perl軟件提取mRNA陣列中IRGs表達(dá)陣列；剔除生存信息缺失及存活天數(shù)小于30的病人；R軟件“WGCNA”函數(shù)對(duì)生存相關(guān)的核心免疫基因聚類，篩選與預(yù)后顯著相關(guān)的基因模塊（模塊中基因數(shù)＞40）；上述得到的核心基因納入單因素Cox回歸篩選出預(yù)后相關(guān)的IRGs（P＜0.001為過(guò)濾條件），在預(yù)后相關(guān)的IRGs和lncRNAs之間進(jìn)行共表達(dá)分析，得到IRlncRNAs（IRGs與lncRNA之間Pearson系數(shù)絕對(duì)值大于0.4，P＜0.001認(rèn)為是IRlncRNAs）。

1.2 風(fēng)險(xiǎn)評(píng)分模型的建立及驗(yàn)證 套索算法（LASSO）是一種變量選擇方法，通過(guò)構(gòu)造一個(gè)懲罰系數(shù)得到比較精練的模型。使用R軟件“glmnet”包使用R“caret”包按1∶1的比例將總樣本隨機(jī)分為訓(xùn)練集和驗(yàn)證集，在訓(xùn)練集中使用R“glmnet”包對(duì)上述得到的IRGs和lncRNAs進(jìn)行Lasso回歸分析，創(chuàng)建預(yù)后模型。風(fēng)險(xiǎn)評(píng)分公式=∑Ni=1（Ei×Ci）（Ei為每個(gè)元素的表達(dá)量，Ci為相對(duì)應(yīng)的回歸系數(shù)）。根據(jù)評(píng)分中位值將所有樣本分為高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組，分析高風(fēng)險(xiǎn)組與低風(fēng)險(xiǎn)組HCC病人臨床病理特征的差異（χ2檢驗(yàn)）；再使用R軟件繪制訓(xùn)練集組、測(cè)試集組、總樣本組病人Kaplan-Meier曲線、時(shí)間依賴性受試者操作特征（ROC）曲線、風(fēng)險(xiǎn)因子圖對(duì)模型進(jìn)行評(píng)價(jià)；最后將HCC病人臨床病理特征及模型風(fēng)險(xiǎn)評(píng)分納入多因素Cox回歸分析，判斷模型預(yù)測(cè)病人預(yù)后的價(jià)值（P＜0.05為差異有統(tǒng)計(jì)學(xué)意義）。

1.3 不同風(fēng)險(xiǎn)組生存差異的機(jī)制分析 為了解不同風(fēng)險(xiǎn)組病人預(yù)后差異的分子機(jī)制，使用R語(yǔ)言limma包篩選高低風(fēng)險(xiǎn)病人差異表達(dá)的基因（FDR＜0.05，│Log2FC│＞1為過(guò)濾條件），并對(duì)差異基因進(jìn)行GO（基因本體論）和KEGG（京都基因與基因組百科全書）分析，繪制富集結(jié)果的柱狀圖和氣泡圖（GO分析篩選條件為FDR＜0.05，基因計(jì)數(shù)＞5；KEGG通路分析篩選條件FDR＜0.05）。

2 結(jié)果

2.1 預(yù)后相關(guān)的IRGs及IRlncRNAs 剔除預(yù)后信息缺失及生存天數(shù)＜30 d的樣本后共得到342個(gè)有效樣本進(jìn)行WGCNA分析，結(jié)果顯示Power值為3時(shí)為最優(yōu)無(wú)尺度網(wǎng)絡(luò)，動(dòng)態(tài)剪切法得到6個(gè)模塊（圖1A），黃色模塊及藍(lán)色模塊中的基因?yàn)轭A(yù)后相關(guān)的核心基因（圖1B），模塊中共提取388個(gè)預(yù)后相關(guān)的核心基因，上述基因與生存數(shù)據(jù)合并后使用R“survival”包和單因素Cox回歸進(jìn)行IRGs篩選；共得到42個(gè)預(yù)后相關(guān)的IRGs，共表達(dá)分析得到42個(gè)預(yù)后相關(guān)IRlncRNAs（圖2）。42個(gè)預(yù)后相關(guān)的IRGs均為HCC病人預(yù)后危險(xiǎn)基因；在42個(gè)預(yù)后相關(guān)IRlncRNAs中，F(xiàn)11-AS1、AC015908.3、TMEM220-AS1、AP001065.3、AC115619.1是HCC病人預(yù)后的保護(hù)IRlncRNAs，其余為預(yù)后危險(xiǎn)IRlncRNAs。

圖1 通過(guò)WGCNA篩選肝細(xì)胞癌（HCC）預(yù)后相關(guān)的基因模塊：1A為WGCNA聚類樹；1B為與HCC病人生存相關(guān)的基因模塊

圖2 HCC預(yù)后相關(guān)的IRGs（A）和IRlncRNAs（B）森林圖

2.2 LASSO回歸構(gòu)建模型 在訓(xùn)練集中對(duì)上述預(yù)后相關(guān)的IRGs及IRlncRNAs進(jìn)行LASSO回歸分析，選取誤差最小的λ值確定最佳模型（圖3），篩選出與預(yù)后顯著相關(guān)的6個(gè)IRGs（NFYC、PSMD1、PSMD6、BCL10、TXLNA、S100A16）和7個(gè)IRlncRNAs（MSC-AS1、AC015908.3、AP001065.3、AL139384.1、AL031985.3、ZFPM2-AS1、AL117336.2），計(jì)算各個(gè)元素的回歸系數(shù)，構(gòu)建預(yù)后模型。

圖3 基于IRGs和IRlncRNAs構(gòu)建的HCC預(yù)后模型LASSO回歸構(gòu)建模型：3A為L(zhǎng)ASSO回歸中的調(diào)諧平行度；3B為特征相對(duì)于log2(1ambda)的LASSO系數(shù)剖面

風(fēng)險(xiǎn)評(píng)分=0.062×NFYC+0.127×PSMD1+0.061×PSMD6+0.024×BCL10+0.1×TXLNA+0.009×S100A16+0.031×MSC-AS1-0.334×AC015908.3-0.007×AP001065.3+0.04×AL139384.1+0.053×AL031985.3+0.083×ZFPM2-AS1+0.444×AL117336.2。在模型中，NFYC、PSMD1、PSMD6、BCL10、TXLNA、S100A16、MSC-AS1、AL139384.1、AL031985.3、ZFPM2-AS1、AL117336.2為病人風(fēng)險(xiǎn)評(píng)分的危險(xiǎn)變量，AC015908.3、AP001065.3為病人風(fēng)險(xiǎn)評(píng)分的保護(hù)變量。

2.3 模型的評(píng)價(jià) 根據(jù)評(píng)分中位值將342位病人分為高、低風(fēng)險(xiǎn)組。臨床相關(guān)性分析顯示：不同風(fēng)險(xiǎn)組病人的組織學(xué)分級(jí)Grade（P=0.001）、臨床分期（P=0.005）、T分期（P=0.010）差異有統(tǒng)計(jì)學(xué)意義（表1）；模型在訓(xùn)練集（171例）、測(cè)試集（171例）及總樣本集中分別進(jìn)行驗(yàn)證，生存分析一致表明，高風(fēng)險(xiǎn)組較低風(fēng)險(xiǎn)組病人總生存期明顯降低（均P＜0.05）（圖4A～4C）；1年時(shí)間依賴性ROC曲線提示，訓(xùn)練集、測(cè)試集、總樣本集曲線下面積95%CI分別0.85（0.79，0.91）、0.81（0.73，0.89）、0.83（0.78，0.88），提示該模型在預(yù)測(cè)HCC預(yù)后表現(xiàn)良好；風(fēng)險(xiǎn)熱圖顯示低風(fēng)險(xiǎn)組病人生存狀態(tài)更好（圖4D）；多因素分析顯示模型評(píng)分、臨床分期是HCC病人獨(dú)立的預(yù)后因子（均P＜0.05），見表2～4。

表1 模型風(fēng)險(xiǎn)評(píng)分與HCC病人臨床病理特征的相關(guān)性/例

表2 訓(xùn)練集預(yù)后的單因素及多因素回歸分析

表3 測(cè)試集預(yù)后的單因素及多因素回歸分析

表4 總樣本預(yù)后的單因素及多因素回歸分析

圖4 肝細(xì)胞癌（HCC）模型的評(píng)價(jià)及驗(yàn)證：4A為訓(xùn)練集組生存曲線；4B為測(cè)試集組生存曲線；4C為總樣本組生存曲線；4D為風(fēng)險(xiǎn)得分；4E為生存狀態(tài)；4F為模型變量的表達(dá)熱圖

2.4 高低風(fēng)險(xiǎn)病人差異基因的GO、KEGG分析 基于模型的不同風(fēng)險(xiǎn)組共存在1 241個(gè)差異表達(dá)基因。GO分析顯示：差異基因主要富集的功能為有絲分裂、染色體分離、DNA復(fù)制起點(diǎn)的結(jié)合、微管蛋白結(jié)合、花生四烯酸單加氧酶活性、膠原結(jié)合等（圖5A）；KEGG分析顯示：差異基因主要富集的通路有PI3K-AKT、細(xì)胞周期、人乳頭瘤病毒感染、黏著斑、ECM-受體相互作用、細(xì)胞衰老、蛋白質(zhì)消化吸收等（圖5B）。

圖5 肝細(xì)胞癌不同風(fēng)險(xiǎn)組差異基因的基因本體論（GO）（5A）及京都基因與基因組百科全書（KEGG）（5B）分析棒狀圖及氣泡圖

3 討論

盡管HCC的診斷和治療方式取得較大進(jìn)步，但具有不同分子特征的病人其治療及預(yù)后存在明顯差異，因此區(qū)分不同風(fēng)險(xiǎn)病人尤為重要。構(gòu)建預(yù)后模型，能有效評(píng)估病人潛在的臨床狀態(tài)及結(jié)局，并篩選出與HCC發(fā)生發(fā)展及預(yù)后相關(guān)的分子標(biāo)記，在個(gè)性化治療和預(yù)后預(yù)測(cè)方面彰顯優(yōu)勢(shì)。隨著技術(shù)的進(jìn)步，免疫治療給人類抗癌帶來(lái)了希望，同時(shí)在個(gè)性化的治療時(shí)代，關(guān)鍵生物標(biāo)志物及基因組的研究可引領(lǐng)腫瘤治療的發(fā)展。鑒于此，部分研究者開始通過(guò)對(duì)IRGs或轉(zhuǎn)錄產(chǎn)物的研究來(lái)評(píng)估病人風(fēng)險(xiǎn)，但既往研究通常只關(guān)注單一類型的RNA的表達(dá)，如mRNA、lncRNA等，考慮到生物調(diào)控的復(fù)雜性以及不同RNA分子在癌癥發(fā)展和進(jìn)展中的密切聯(lián)系，我們首次構(gòu)建了HCC多類型RNA的預(yù)后模型，并分析了不同風(fēng)險(xiǎn)病人的臨床特征，驗(yàn)證結(jié)果顯示該模型預(yù)測(cè)HCC病人生存性能良好；多因素Cox回歸提示該模型預(yù)測(cè)性能相對(duì)于傳統(tǒng)臨床特征是獨(dú)立的；最后我們探討了高低風(fēng)險(xiǎn)病人預(yù)后差異的可能分子機(jī)制，以期對(duì)HCC的臨床決策及深入研究提供一定幫助。

在本模型納入的6個(gè)IRGs和7個(gè)IRlncRNAs中，有部分已被報(bào)道與HCC的發(fā)展和預(yù)后有關(guān)。如Tan等［9］發(fā)現(xiàn)PSMD1的高表達(dá)與HCC的不良預(yù)后顯著相關(guān)；Cheng等［10］發(fā)現(xiàn)BCL10在HCC組織中存在高頻點(diǎn)突變，但具體致病機(jī)制仍不清楚；S100A16則可以參與HCC病人的體液免疫，并影響病人預(yù)后［11］。至于IRlncRNA方面，MSC-AS1被發(fā)現(xiàn)可預(yù)測(cè)HCC病人的預(yù)后、免疫細(xì)胞浸潤(rùn)水平及免疫治療反應(yīng)［12］；Kou等［13］發(fā)現(xiàn)MSC-AS1的下調(diào)可抑制HCC細(xì)胞增殖、遷移、侵襲和促進(jìn)HCC細(xì)胞凋亡；ZFPM2-AS1可通過(guò)靶向調(diào)節(jié)miR-653/GOLM1軸［14］、miR-576-3p/HIF-1α軸［15］等途徑促進(jìn)HCC細(xì)胞增殖、侵襲，并與HCC不良預(yù)后相關(guān)。另外值得注意的是，AC015908.3的腫瘤干細(xì)胞相關(guān)功能新近在HCC中被鑒定出來(lái)［16］；AL031985.3則在HCC免疫、自噬等功能中發(fā)揮一定作用［17］，提示這些lncRNAs可能與HCC的發(fā)生發(fā)展密切相關(guān)。而關(guān)于NFYC、PSMD6、TXLNA、AP001065.3、AL139384.1、AL117336.2在HCC中的作用研究則較少，具體功能有待研究進(jìn)一步揭示。最后GO分析顯示高低風(fēng)險(xiǎn)病人風(fēng)險(xiǎn)差異基因主要參與細(xì)胞遺傳物質(zhì)的復(fù)制等事件，KEGG分析提示這些基因與某些腫瘤相關(guān)通路密切聯(lián)系，如“PI3K-AKT”“細(xì)胞周期”等，表明HCC病人預(yù)后差異涉及多途徑、多靶點(diǎn)，為HCC的多靶點(diǎn)聯(lián)合治療提供了理論支持。

HCC已進(jìn)入免疫治療時(shí)代，篩選免疫相關(guān)的生物標(biāo)志物用于預(yù)測(cè)HCC病人預(yù)后非常重要和必要［18-19］。本研究是首個(gè)IRGs及IRlncRNAs聯(lián)合的HCC預(yù)后模型的研究，該模型風(fēng)險(xiǎn)評(píng)分能夠作為HCC病人的獨(dú)立預(yù)后因素，將有助于臨床醫(yī)生對(duì)HCC病人的風(fēng)險(xiǎn)管理和決策制定。然而，我們的研究仍有一些不足。首先，本研究是基于國(guó)外數(shù)據(jù)庫(kù)的分析，基于東西方HCC病人特征的不同，該模型在中國(guó)人群的應(yīng)用是否同樣具有價(jià)值尚未可知；其次，本模型是在公共數(shù)據(jù)的回顧性分析基礎(chǔ)上構(gòu)建，可能存在選擇性偏倚；最后HCC具有高度異質(zhì)性及病人個(gè)體差異，我們關(guān)注的隊(duì)列是否具有代表性尚需在前瞻的獨(dú)立臨床隊(duì)列中進(jìn)行驗(yàn)證?？傊狙芯勘砻?，由IRGs及IRlncRNAs聯(lián)合構(gòu)建的預(yù)后模型在預(yù)測(cè)HCC病人的生存及臨床特征上有一定價(jià)值，可能為HCC的免疫學(xué)研究及預(yù)后預(yù)測(cè)提供有用工具。