基于TCGA數據庫的肝細胞癌鐵死亡相關LncRNA風險預測模型的構建

王少波 郭躍 曾昊 龐志剛

2020年肝癌的發病率位于全世界惡性腫瘤第六位，病死率位于第三位[1]。肝癌最主要的危險因素是乙型肝炎病毒感染和黃曲霉毒素暴露，此外，酒精、非酒精性脂肪肝、糖尿病、吸煙等危險因素也與肝癌有密切關系[1，2]。肝細胞癌（Hepatocellular carcinoma，HCC）是肝癌最常見的類型，占肝癌病例的75%～85%。HCC有手術切除、消融、化療栓塞和肝移植等治療方法[3]，但由于病因復雜和腫瘤的高度異質性[4]，腫瘤往往進展至晚期才被診斷，有效的治療方法有限，往往導致預后不良[5]。因此有必要開發一種新的HCC預后預測模型為HCC的診療策略提供新思路。

腫瘤耐藥的主要原因之一是癌細胞在死亡的執行過程中存在缺陷[6]，鐵死亡是2012年由Dixon等人發現的一種由脂質過氧化和活性氧累積所驅動的調節性細胞死亡方式[7]。近年來有研究發現誘導腫瘤細胞鐵死亡可用于癌癥治療，特別是對傳統療法有抵抗力的侵襲性惡性腫瘤[8]。

長鏈非編碼RNA（LncRNA）是長度超過200個堿基的非編碼RNA，在基因調控方面發揮重要作用[9]。此外，有研究發現LncRNA在腫瘤的發生、發展、侵襲、轉移過程中發揮重要作用[10]。目前關于鐵死亡相關LncRNA的研究較少，最近的一項研究發現核LncRNA LINC00336作為競爭的內源性RNA，在肺癌中上調并發揮癌基因作用，通過與RNA結合蛋白ELAVL1相結合抑制腫瘤細胞鐵死亡[11]。在一項研究中發現，LncRNA LINC00618在人類白血病細胞中通過增加脂質ROS和鐵的水平，同時降低SLC7A11的表達來促進鐵死亡，此外還能促進長春新堿誘導的鐵死亡和細胞凋亡[12]。然而，目前關于肝癌患者鐵死亡相關LncRNA特征及其與總生存期的相關性研究較少，本研究中，我們基于The Cancer Genome Atlas（TCGA）數據庫構建鐵死亡相關LncRNA預測模型，并對其功能富集和免疫相關功能進行分析，探討其可能存在的機制。

1 材料與方法

1.1 數據預處理 通過在TCGA數據庫下載FPKM標準化的374例HCC組織和50例正常組織的RNA-seq轉錄組數據，同時下載對應臨床數據，主要臨床信息包括年齡、性別、生存時間、生存狀態、組織學分級、腫瘤TNM分期。鐵死亡相關基因通過FerrDb數據庫下載。FerrDb數據庫是一個由人工收集管理的鐵死亡相關標志物和調控因子以及鐵死亡關聯疾病的數據庫[13]。通過FerrDb數據庫共下載鐵死亡相關基因382個，其中有246個鐵死亡相關基因在TCGA轉錄組數據中有表達。使用Pearson相關性分析法將鐵死亡基因和LncRNA進行共表達分析確定鐵死亡相關LncRNA，篩選條件為r>0.5，P<0.001。

1.2 差異基因分析 使用R軟件4.0.3版本“limma”包對374例HCC組織和50例正常組織的鐵死亡相關基因進行差異表達分析，篩選條件為P<0.05，|log2FC|>1。利用R軟件“clusterProfiler”包對鐵死亡相關差異表達基因進行基因本體論（GO）富集分析，包括生物學過程（BP）、細胞組成（CC）和分子功能（MF），同時進行京都基因與基因組百科全書（KEGG）通路富集分析。

1.3 鐵死亡預后相關LncRNA模型的構建與驗證 對上述鐵死亡相關LncRNA使用R軟件“survival”包進行單因素COX回歸分析，獲取與預后相關的LncRNA（P<0.01），再對預后相關的LncRNA進行多因素COX回歸分析，建立預后模型，計算每例患者的風險分數，根據中位風險值將病例分為低風險組和高風險組。使用KM法分析低風險組和高風險組的生存時間差異，從而確定風險評分對HCC患者預后的預測作用。對患者年齡、性別、組織學分級、臨床分期、風險評分分別進行單因素和多因素COX獨立預后分析，從而觀察風險評分是否可以獨立于其他臨床特征作為新的預后因子。通過受試者工作特征（ROC）曲線評價預后模型的預測性能。

1.4 通路富集分析和免疫相關功能研究 通過GSEA富集分析探索鐵死亡相關LncRNA預測模型下高、低風險組富集通路特征，通過ssGSEA算法評估高、低風險組的免疫細胞成分，并對高、低風險組腫瘤浸潤性免疫細胞亞群、免疫功能和免疫檢查點之間的差異進行評估。

1.5 統計學方法 使用R軟件4.0.3版本進行數據分析和圖形可視化處理，采用 COX 回歸模型篩選與HCC患者預后相關的LncRNA并構建風險評分模型。采用時間依賴的ROC曲線評估模型的預測性能，根據此模型對HCC患者進行Kaplan-Meier生存分析。P<0.05為差異有統計學意義。

2 結果

2.1 鐵死亡相關基因的富集分析 通過對TCGA轉錄組數據中的鐵死亡相關基因進行差異分析，共得到84個鐵死亡相關基因，其中上調基因13個，下調基因71個。GO富集分析結果顯示鐵死亡相關基因參與BP中細胞化學應激和氧化應激反應，在CC中主要參與煙酰胺腺嘌呤二核苷酸磷酸（NADPH）氧化酶復合體和次級溶酶體的構成，在MF方面與NADPH氧化酶活性有關。KEGG分析顯示高表達基因主要涉及鐵死亡、癌細胞碳代謝、5-羥色胺能突觸和VEGF信號通路等。見圖1。

圖1 鐵死亡相關差異表達基因的GO（A）和KEGG（B）分析

2.2 構建鐵死亡相關LncRNA預后模型 通過基因共表達分析，共得到319個鐵死亡相關LncRNA，單因素COX分析得到了58個與預后顯著相關的LncRNA，多因素COX分析最終確定了22個LncRNA為獨立預后因素，風險得分=-0.409×AC022007.1表達值+0.786×MKLN1-AS表達值+0.870×SNHG12表達值+0.635×Z95115.1表達值+1.077×PRRT3-AS1表達值+1.586×SREBF2-AS1表達值+1.094×MIR210HG表達值+1.248×AL050341.2表達值+0.899×LUCAT1表達值+1.342×NRAV表達值+1.396×LINC01224表達值+1.108×MIR4435-2HG表達值+1.220×LINC00205表達值+0.487×PTOV1-AS1表達值+0.320×AL671710.1表達值×1.069×ZFPM2-AS1表達值+0.526×AL357079.1表達值+1.639×AC145 207.5表達值+1.073×AC099850.3表達值+1.846×SNHG4表達值+2.689×KDM4A-AS1表達值+0.157×AC091057.1表達值。

2.3 驗證鐵死亡相關LncRNA預后模型 以風險評分的中位值作為截斷值將患者分為高風險組和低風險組，見圖2，Kaplan-Meier分析顯示高風險組患者的生存時間顯著小于低風險組（P<0.001）。見圖3。將患者年齡、性別、組織學分級、臨床分期、風險評分進行獨立預后分析，結果顯示臨床分期和風險評分在單因素和多因素COX分析中差異均具有統計學意義（P<0.001），ROC曲線對LncRNA預測模型的1年的生存預測性能評價顯示AUC=0.811，同時得到該模型3年和5年的生存預測性能AUC值分別為0.778和0.752。見圖4。

圖2 風險評分與相關LncRNA表達熱圖及生存時間點圖

圖3 高、低風險組的Kaplan-Meier分析

圖4 鐵死亡相關LncRNA單因素（A）和多因素COX分析（B）以及模型對1、3、5年生存率預測性能的AUC值（C）

2.4 基因富集分析 GSEA結果提示鐵死亡相關LncRNA預測模型的高風險組中富集多條免疫和腫瘤相關通路，如錯配修復、Wnt信號通路、mTOR信號通路、MAPK信號通路、TGF-β信號通路、p53信號通路、VEGF信號通路、RIG-Ⅰ樣受體信號通路、T細胞受體信號通路、FcγR介導的吞噬作用、Toll樣受體信號通路、自然殺傷細胞介導的細胞毒作用等，而在低風險組中則富集了多條與合成代謝相關通路，如脂肪酸代謝、膽汁酸合成、異生物質代謝和部分氨基酸代謝等。見圖5。

圖5 鐵死亡相關LncRNA的GSEA

2.5免疫相關分析 基于ssGSEA算法對高、低風險組腫瘤免疫浸潤程度進行評估，結果顯示免疫細胞和免疫功能中B細胞、NK細胞、中性粒細胞、肥大細胞、輔助性T細胞、細胞毒活性、MHCⅠ類分子、Ⅰ型和Ⅱ型干擾素反應在高、低風險組間的差異有統計學意義（P<0.05）。由于免疫檢查點在免疫治療方面發揮重要作用，我們進一步探討了高、低風險組間免疫檢查點的差異，結果顯示CTLA4、CD44、HHLA2、TMIGD2、TNFSF15等在兩組患者中的表達差異有統計學意義（P<0.05）。見圖6。

圖6 免疫細胞（A）與免疫相關功能（B）在高、低風險組中的差異及免疫檢查點在高、低風險組中的表達（C）

3 討論

腫瘤靶向治療過程中產生的耐藥機制已經成為晚期腫瘤治療失敗的重要原因[14]，有研究表明，以FSP1/CoQ10軸為靶點的鐵死亡療法對于由細胞可塑性引起的耐藥機制有巨大的應用潛力[15]，可用于預防各種來源的惡性腫瘤的治療抵抗或轉移。由于鐵死亡在腫瘤發生發展中的重要作用，其相關的LncRNA也受到了特別關注[16]。本研究基于TCGA數據建立了一種新的鐵死亡相關LncRNA預后預測模型，該模型能較好地區分高低風險人群，與傳統的臨床預后因素相比，該模型對于HCC預后預測的效能明顯提高。本模型共納入22個LncRNA，根據LncRNA在高、低風險組間差異表達的熱圖可以看出，模型中所有LncRNA都傾向于在高風險組中高表達，其中部分LncRNA被報道與各種腫瘤的發生發展密切相關，如有研究表明，MIR210HG在非小細胞肺癌的組織和細胞中過表達，沉默MIR210HG可抑制腫瘤細胞增殖、遷移和侵襲，同時增加腫瘤細胞的放射敏感性和化學抗性[17]。MIR210HG高表達還與HCC的臨床分期晚、腫瘤體積大、有血管侵犯及組織分化不良有關，是影響患者總體生存的獨立不良預后因素[18]。另有研究表明，LINC01224可通過LINC01224/miR-330-5p/CHEK1軸參與HCC的發生發展[19]。在胃癌中LINC01224還可通過介導miR-193a-5p調控CDK8影響胃癌的生物學行為[20]。

在HCC腫瘤微環境中，由不同免疫抑制細胞亞群介導的促腫瘤免疫反應在驅動免疫逃避中起關鍵作用。通過對高、低風險組免疫評分分析發現，高風險組人群中的NK細胞比例和Ⅰ型和Ⅱ型干擾素反應較低，提示高風險組中的腫瘤免疫抑制程度較高并導致預后較差。近年來基于免疫檢查點抑制的免疫療法已被證實能提升晚期腫瘤患者的生存率，本研究中高、低風險組間免疫檢查點的顯著差異提示兩組人群可能對免疫療法的敏感程度不同。

本研究通過對TCGA數據庫HCC基因表達譜數據的挖掘，構建了鐵死亡相關的LncRNA預測模型，可以較好地預測患者預后，相比于單個LncRNA位點具有更好的預測效果。然而本研究仍存在一定的不足之處，如該預測模型是基于公共數據庫TCGA的臨床信息和基因表達信息所建立的，還需要更多的臨床數據對其作用進行驗證，另一方面，未對該模型納入的LncRNA進行腫瘤機制實驗研究，對LncRNA如何影響腫瘤增殖轉移及預后尚不清楚，有待于進一步研究。