周丁杰 戈 偉 曹德東

武漢大學(xué)人民醫(yī)院腫瘤中心，湖北武漢 430060

肝細(xì)胞癌（hepatocellular carcinoma，HCC）是導(dǎo)致全球癌癥死亡的第四大原因[1]，每年造成80 多萬人死亡。HCC 患者的5 年生存率僅為5%～30%[2]。盡管目前靶向治療和免疫療法已進(jìn)入臨床，但肝癌患者的預(yù)后仍然較差。尋找新的治療靶點及檢測預(yù)后生物標(biāo)志物仍是當(dāng)下研究的重點。衰老相關(guān)分泌表型（senescence-associated secretory phenotype，SASP）是指由衰老細(xì)胞分泌的細(xì)胞因子、生長因子和酶等分子，這些分子可引發(fā)慢性炎癥反應(yīng)，促進(jìn)腫瘤的形成和發(fā)展[3]。研究表明，SASP 在肝癌中的作用與其成分有關(guān)。例如，SASP 中的白細(xì)胞介素（interleukin，IL）-6、IL-8 等因子被認(rèn)為是肝癌的重要促進(jìn)因子，能夠促進(jìn)肝癌細(xì)胞的增殖和侵襲，抑制免疫系統(tǒng)的反應(yīng)[4]；而轉(zhuǎn)化生長因子β 等因子則能夠促進(jìn)肝癌干細(xì)胞的增殖和轉(zhuǎn)移[5]。在肝癌中，SASP 參與了腫瘤微環(huán)境的形成和轉(zhuǎn)化過程[6]，并且lncRNA 作為SASP 中的重要組成部分，也在肝癌中發(fā)揮了重要作用。研究顯示，GAS5 參與了肝癌中的SASP，可以通過抑制IL-6 的表達(dá)來抑制腫瘤細(xì)胞的增殖和侵襲，并提高肝癌患者的生存率[7]。另外，GAS5 還可以通過調(diào)節(jié)核因子-κB 信號通路來抑制肝癌細(xì)胞的轉(zhuǎn)移和侵襲[8]。lncRNA-H19 可以通過介導(dǎo)細(xì)胞內(nèi)活性氧水平的升高，促進(jìn)IL-6 和IL-8 等SASP 相關(guān)因子的分泌，從而促進(jìn)肝癌細(xì)胞的增殖和侵襲[9]。上述結(jié)果顯示，SASP 相關(guān)的lncRNA 有望成為肝癌診斷和治療的新靶點?；赟ASP 相關(guān)的lncRNA，本研究建立了一個包含6個lncRNA 的風(fēng)險模型，旨在優(yōu)化肝癌臨床治療，改善肝癌患者預(yù)后[9-16]。

1 材料與方法

1.1 實驗材料

在癌癥基因組學(xué)圖譜（the cancer genome stlas，TCGA）數(shù)據(jù)庫[10]中檢索并下載了424 例HCC 患者樣本的RNA 測序數(shù)據(jù)，臨床資料包括年齡、性別、生存時間、生存狀態(tài)、臨床分期及TNM 分期，排除有模糊值和臨床資料不完整以及隨訪時間＜30 d 的樣本后，最終納入241 例HCC 樣本。并按1∶1 的比率將納入樣本按計算機隨機分為實驗組（120 例）和驗證組（121例）。通過Ensembl 網(wǎng)站[11]（http://asia.ensem-bl.org/）對轉(zhuǎn)錄組數(shù)據(jù)注釋以區(qū)分mRNA 和lncRNA 數(shù)據(jù)，提取HCC 的lncRNA 表達(dá)矩陣。

1.2 實驗方法

1.2.1 差異表達(dá)基因及l(fā)ncRNA 的篩選 在GeneCards網(wǎng)站[12]（https://www.genecards.org/）中獲取相關(guān)系數(shù)＞7 的201 個SASP 基因，隨后利用R 軟件中的“l(fā)imma”包[13]對HCC 樣本和正常樣本中的SASP 相關(guān)基因進(jìn)行差異分析，篩選標(biāo)準(zhǔn)為|logFC|≥1 且FDR＜0.05，得到差異基因。隨后設(shè)定相關(guān)系數(shù)為|r|＞0.4 且P＜0.001，對HCC 樣本中的lncRNAs 和差異基因進(jìn)行Spearman 相關(guān)分析，獲取SASP 相關(guān)lncRNAs。

1.2.2 預(yù)后風(fēng)險模型的建立 利用Perl 語言腳本將臨床數(shù)據(jù)和對應(yīng)的lncRNA 矩陣樣本對應(yīng)后合并。使用R 軟件中“survival”包（https://github.com/therneau/survival）對SASP 相關(guān)lncRNAs 進(jìn)行單因素Cox 回歸分析，篩選出HCC 預(yù)后相關(guān)的lncRNAs。采用“glmnet”軟件包[14]進(jìn)行的lasso 回歸分析對篩選出的lncRNA 進(jìn)一步限制，采用十倍交叉驗證來獲得lambda 懲罰參數(shù)的理想值。經(jīng)過lasso 回歸解決共線性問題，再經(jīng)過多因素Cox 回歸分析后，將得到的lncRNAs 納入風(fēng)險模型中。lncRNA 模型的計算公式為：風(fēng)險得分=（coefi×xi）。系數(shù)值用“Coef”表示，所選lncRNAs 的表達(dá)值用“x”表示。

1.2.3 風(fēng)險模型的評估 應(yīng)用風(fēng)險模型計算出每個HCC 患者的風(fēng)險分?jǐn)?shù)。根據(jù)風(fēng)險評分的中位數(shù)，將HCC 患者分為低風(fēng)險組和高風(fēng)險組。對高風(fēng)險組和低風(fēng)險組用Kaplan-Meier 曲線評估生存差異，繪制受試者操作特征（receiver operating characteristic，ROC）曲線評估風(fēng)險模型靈敏度與特異度。單因素和多因素Cox回歸分析比較風(fēng)險模型和其他臨床病理變量（年齡、性別、臨床分期、TNM 分期）的關(guān)系，運用R 軟件“ggpubr”包（https://rpkgs.datanovia.com/ggpubr/）繪制了不同臨床病理參數(shù)下的高、低風(fēng)險組之間的生存差異。

1.2.4 風(fēng)險模型的免疫圖譜分析 利用Perl 語言將矩陣與含風(fēng)險值的樣本文件合并，得到高、低風(fēng)險組中不同樣本基因的表達(dá)量新矩陣。使用“GSVA”工具[15]探索樣本的風(fēng)險評分與樣本中免疫細(xì)胞浸潤和經(jīng)典的免疫檢查點相關(guān)標(biāo)志物之間的關(guān)系。

1.3 統(tǒng)計學(xué)方法

所有分析皆由R 軟件（4.1.3）和Perl 完成。散點圖展示樣本的風(fēng)險評分與生存狀態(tài)的關(guān)系，生存曲線采用Kaplan-Meier（K-M）法生成，使用單因素和多因素Cox 比例風(fēng)險回歸模型評估SASP 相關(guān)lncRNAs預(yù)后模型，采用ROC 曲線評價預(yù)測模型的可靠性和準(zhǔn)確性。計量資料采用t 檢驗，計數(shù)資料采用χ2檢驗或Fisher確切概率法。以P＜0.05 為差異有統(tǒng)計學(xué)意義。

2 結(jié)果

2.1 HCC 中SASP 相關(guān)差異表達(dá)基因（differential gene expression，DEGs）及l(fā)ncRNA 的篩選

本研究對TCGA 數(shù)據(jù)庫中的HCC 隊列的RNA測序數(shù)據(jù)清洗后共納入241 例樣本。通過“l(fā)imma”包對HCC 樣本和正常樣本進(jìn)行差異分析，得到58 個顯著差異的SASP 基因。其中有36 個上調(diào)基因及22 個下調(diào)基因（圖1A）。熱圖展示了所有的DEGs 在樣本中的表達(dá)情況（圖1B）。并利用Spearman 分析得到DEGs 的1 048 個lncRNAs。

圖1 腫瘤組織和正常組織間的DEGs 表達(dá)

2.2 構(gòu)建lncRNA 風(fēng)險模型

利用Perl 腳本從測序數(shù)據(jù)中分離出SASP 相關(guān)lncRNA 的表達(dá)矩陣，隨后聯(lián)合生存資料，利用單因素Cox 回歸分析確定了64 個預(yù)后相關(guān)且差異表達(dá)的lncRNAs（P＜0.001）。隨后對這些lncRNAs 進(jìn)行了Lasso 回歸，得到14 個SASP 相關(guān)的lncRNAs（圖2A），通過10 輪交叉驗證確定了最佳懲罰參數(shù)值（圖2B）。最后通過多因素Cox 回歸分析，從上述14 個lncRNAs中篩選出6 個lncRNAs，包括AL031985.3、AC124067.4、AC009283.1、AC107021.2、LINC00324、HPN-AS1，構(gòu)成風(fēng)險模型。公式：風(fēng)險評分=AL031985.3 expression×（0.768 307）+AC124067.4 expression×（0.301 32+AC009 283.1 expression×（-0.367 63）+AC107021.2 expression×（0.495 064）+LINC00324 expression×（-0.511 08）+HPNAS1 expression×（-0.751 87）。多因素Cox 回歸分析結(jié)果見表1。圖2C 展示了模型中6 個lncRNAs 在HCC樣本中的表達(dá)情況。本研究使用Cytoscape 來進(jìn)一步可視化lncRNAs 和對應(yīng)的mRNA 的關(guān)系（圖2D，|R2|＞0.4 和P＜0.001）。Sankey 圖顯示，AC009283.1、HPN-AS1 和LINC00324 是保護(hù)因素，而AC107021.2、AC124067.4 和AL031985.3 是風(fēng)險因素（圖2E）。

表1 多因素Cox 回歸分析篩選用于構(gòu)建預(yù)后風(fēng)險模型的SASP 相關(guān)lncRNA

2.3 風(fēng)險模型的評估

2.3.1 高、低風(fēng)險組間的生存分析 生存分析結(jié)果顯示，HCC 患者的預(yù)后隨著風(fēng)險值升高而惡化（圖3A）。K-M 曲線結(jié)果顯示，與高風(fēng)險組比較，低風(fēng)險組患者的預(yù)后更好（P＜0.001）（圖3B）。風(fēng)險值ROC 曲線顯示，1、3、5 年生存率的曲線下面積（AUC）分別為0.812、0.754、0.739（圖3C）。

2.3.2 單因素和多因素Cox 回歸分析 為了進(jìn)一步評估此風(fēng)險模型在HCC 中的應(yīng)用，本研究結(jié)合了其他病理參數(shù)，包括年齡、性別、臨床分期、TNM 分期，在實驗組和驗證組中分別進(jìn)行單因素和多因素Cox 回歸分析。結(jié)果顯示，在實驗組中，風(fēng)險評分的HR 在單變量Cox回歸分析中為1.437（95%CI：1.280，1.614）（圖4A，P＜0.001）。在多因素Cox 回歸分析中，HR 值為1.409（95%CI：1.242，1.599）（圖4B，P＜0.001）。在實驗隊列中，單變量、多變量Cox 的風(fēng)險評分HR 值分別為1.372（95%CI：1.209，1.557）（圖4C，P＜0.001）和1.366（95%CI：1.178，1.585）（圖4D，P＜0.001），提示此風(fēng)險模型可作為HCC 患者的獨立預(yù)后因素。