基于壞死性凋亡相關長鏈非編碼RNA構建肝細胞癌預后模型及藥物治療反應分析

張玉俊，朱琳，趙璇，地力亞爾·吾斯曼江，王巖

作者單位：1新疆醫科大學公共衛生學院，新疆維吾爾自治區 烏魯木齊 830054；

2新疆醫科大學附屬腫瘤醫院，新疆維吾爾自治區 烏魯木齊 830011

肝細胞癌（hepatocellular carcinoma，HCC）是全球第三大最常見的惡性腫瘤，也是中國癌癥相關死亡的第二大主要原因［1］。目前，手術切除是早期肝細胞癌的主要治療方法。然而，由于缺乏特異性癥狀大多數肝細胞癌病例在診斷時已處于晚期階段手術切除率相對較低［2］。盡管近年來在治療方面取得了巨大進展，但由于其侵襲性的生物學特性，復發率高達70%，手術切除或肝移植成功后5 年生存率低于60%［3］。因此，確定新的診斷和預后生物標志物為HCC 的準確診斷、個體化治療和預后預測具有重要意義。

壞死性凋亡（necroptosis）是由受體相互作用蛋白激酶RIPK1，RIPK3 和混合譜系激酶結構域樣蛋白（mixed lineage kinase domain-like protein recombi?nant protein，MLKL）介導的程序性細胞死亡形式，其主要特征為質膜完整性的早期喪失，細胞內容物泄漏和細胞器腫脹［4］。據報道，壞死性凋亡在多種疾病的發生和發展中起著關鍵作用，例如神經退行性疾病，缺血性心血管疾病和癌癥轉移［5］。此外，壞死性凋亡通過在腫瘤微環境中釋放損傷相關分子、細胞因子和趨化因子，從而產生腫瘤促進或抗腫瘤作用［6］。長鏈非編碼RNA（long Noncoding RNA， ln?cRNA）是具有超過200 個核苷酸轉錄本的非編碼RNA，在體內發揮轉錄調節、mRNA 剪接和mRNA 轉錄后調節作用［7］，與腫瘤發展、轉移和免疫途徑密切相關，其可通過與miRNA 相互作用來調節壞死性凋亡相關基因產物，從而調節腫瘤的進程［8］。研究表明壞死性凋亡相關長鏈非編碼RNA（necroptosis-re?lated long noncoding RNAs，NRLs）與多種癌癥的發生發展密切相關，包括胃癌［9］、乳腺癌［10］和肺腺癌［11］，而NRLs 與HCC 間的關系尚未見廣泛報道。因此，進一步闡明NRLs 與 HCC 間的關系對于探索HCC新型治療靶點和改善病人預后至關重要。

在這項研究中，通過共表達網絡分析分析從TCGA 肝細胞癌（TCGA-LIHC）轉錄本數據中挖掘NRLs。然后，通過單變量Cox 和LASSO-Cox 回歸分析，鑒定出了用于預測HCC 病人總生存期（overall survival，OS）的4 個NRLs 特征。此外，基于4 個NRLs 預后特征的藥物敏感性分析為HCC 病人的臨床個體化治療提供了理論基礎。

1 資料與方法

1.1 一般資料本研究所有公共數據均符合以下納入標準：（1）納入病人經HCC 病理確診；（2）納入病人的臨床資料完整；（3）納入病人有OS 的隨訪信息，且OS 大于30 d。從TCGA 數據庫中獲得374 份HCC 組織和50 份正常組織樣本的RNA 測序數據和相應臨床信息。根據壞死性凋亡基因集M24778.gmt和文獻檢索［12］，共獲得67 個壞死性凋亡相關基因。從GENCODE（https：//www.gencodegenes.org/，截至2022 年4 月25 日）網站檢索人類LncRNA 的GTF 注釋文件。從以前的文獻中獲得47 個與免疫檢查點相關的基因［11］。由于TCGA 數據是公開的，本研究不需要當地倫理委員會的批準。然而，本研究嚴格遵守了TCGA數據訪問政策和出版指南。

1.2 方法

1.2.1 壞死性凋亡相關LncRNA 的鑒定使用從GENCODE 網站檢索到的人類LncRNA 的GTF 注釋文件，從TCGA-LIHC RNA-seq數據中鑒定HCC相關LncRNA。隨后，通過Pearson 相關分析HCC 樣本中壞死性凋亡相關基因與鑒定出的LncRNA 的共表達關系，過濾標準：（1）|Pearson 相關系數|>0.4；（2）P<0.001。為了證明NLRs 和相應的靶標mRNA 之間的相互調控連接，使用“igraph”R 包可視化lncRNAmRNA 網絡。最后，采用“limma”R 包鑒定差異表達NRLs，標準為：（1）|log2 差異倍數（fold change，FC）|>1；（2）錯誤發現率（FDR，P值由Benjamini& Hoch?berg 校正進行調整）<0.05，并采用火山圖進行可視化。

1.2.2 預后模型的建立與驗證首先，應用單變量Cox 回歸確定預后相關的NRLs（P<0.05），為防止模型的過度擬合使用 “glmnet”和“survminer”R 包進行LASSO-Cox 回歸（使用10 倍交叉驗證估計的懲罰參數）篩選出與HCC 預后相關的最佳NRLs，并構建預后模型。接下來，根據NRLs 的標準化表達水平和從LASSO-Cox 回歸分析中得出的回歸系數計算每個HCC 病人的生存風險評分。計算公式如下：風險評分=expNRL1×β1 + expNRL2×β2+ expNRL3×β3+… + expNRLn ×βn，“exp”代表已鑒定的NRL 的表達水平，“β”為通過LASSO-Cox 回歸確定的系數。根據風險評分中位值將病人劃分為高風險組和低風險組，使用“survival”R包分析高低風險組間OS是否存在差異（P<0.05），并采用“survminer”R 包分析臨床病理特征分層高低風險組間OS 的差異（P<0.05）以了解NRLs 預后特征的適用性。隨后，進行單變量-Cox 和多變量-Cox 回歸分析，以探索風險評分特征是否為HCC 病人的獨立預后指標（P<0.05）。最后，使用 “timeROC”“survival”R 包繪制了時間依賴的受試者操作特征（receiver operating characteristic curve，ROC）曲線以評估風險評分模型對1 年、3 年和5 年OS 的預測能力，并采用“ROCR” R 包繪制1年ROC曲線比較風險評分與傳統臨床特征指標的預測能力。

1.2.3 列線圖的構造和校準使用“rms” R 包，建立整合臨床病理學特征和壞死性凋亡相關長鏈非編碼RNA 風險評分（NRLs risk-score）的列線圖以預測HCC 病人的1 年、3 年和5 年OS。隨后，開發了校準曲線來說明實際結果與預測結果之間的一致性。

1.2.4 基因集富集分析使用基因集富集分析（gene set enrichment analyses ，GSEA）軟件4.1.2（http ：//www.gsea-msigdb.org /gsea/index.jsp）鑒定兩組中差異表達的KEGG 途徑，顯著富集的生物過程和途徑閾值設置為P<0.05，FDR<0.25。結果由“gri?dExtra”“grid”和“ggplot2” R 包進行可視化。

1.2.5 腫瘤微環境分析腫瘤微環境（tumor micro?environment，TME）由腫瘤、基質和免疫細胞組成。基于7 種平臺（XCELL、TIMER、QUANTISEQ、MCP?COUNTER、EPIC、CIBERSORT-ABS 和CIBERSORT）估計了TCGA-LIHC 數據集中免疫細胞的表達，并采用Spearman 相關性分析評估了免疫細胞亞群與風險評分的關系（P<0.05）。然后，使用“estimate”R 包計算TME 評分，包括基質評分（StromalScore）、免疫評分（ImmuneScore）以及總的估計評分（ESTIMATE?Score），并采用Wilcoxon 秩和檢驗比較腫瘤TIME 在高低風險組間的差異（P<0.05）。通過單樣本基因集富集分析（single-sample GSEA ，ssGSEA）計算16 個免疫細胞和13 種免疫相關通路活性在高低風險組間的浸潤評分差異（通過barplot 可視化）。最后，使用“ggpubr”R 包對高低風險組間的免疫檢查點激活進行比較。

1.2.6 基于NRLs 特征篩選疾病潛在藥物為了預測基于NRLs 特征兩種不同風險群體的HCC 病人對化療藥物的反應，使用“pRRophetic”R 包來評估第8屆美國癌癥聯合委員會（AJCC）指南強烈推薦的20種普通化療藥物的半最大抑制濃度（half-maximal inhibitory concentration，IC50）。并采用Wilcoxon 符號秩檢驗分析高低風險組間IC50值的差異（P<0.05）。

1.3 統計學方法本研究中使用的統計分析工具是R 軟件4.0.2。采用Mann-WhitneyU檢驗進行高低風險組間的差異分析，運用Kaplan-Meier 方法和對數秩檢驗比較高低風險組病人的OS。雙側檢驗，P<0.05認為差異有統計學意義。

2 結果

本研究納入來自TCGA 數據庫的343 例HCC 病人以1∶1 的比例隨機分配到訓練集（n=172）或測試集（n=171）。研究流程如圖1所示。

圖1 研究流程圖

2.1 鑒定HCC 病人中壞死性凋亡相關LncRNA經篩選從TCGA 數據庫中共獲得343 個HCC 組織和50個正常組織的樣本。根據人類LncRNA的GTF注釋文件，從TCGA-LIHC 基因表達文件中鑒定出16 774 個LncRNA。同時，根據壞死性凋亡基因集共鑒定出67 個壞死性凋亡基因。通過Pearson 相關分析最終獲得989 個NRLs（圖2A）。差異分析結果顯示，760 個NRLs（9 個下調，751 個上調）在HCC 和正常組織中差異表達（圖2B）。

圖2 鑒定HCC病人中壞死性凋亡相關LncRNA:A為壞死性凋亡基因與 lncRNA 之間的相關網絡；B為760個差異表達NRLs的火山圖

2.2 模型的構建與驗證在TCGA 訓練集中使用單變量Cox回歸分析，獲得了47個與OS顯著相關的NRLs（P<0.05），并制作了熱圖（圖3A、3B）。為了避免過度擬合并提高預后特征的準確性，對這些NRLs進行了LASSO 懲罰的Cox 分析（圖3C、3D），并在交叉驗證誤差最小點選擇了4 個NRLs，使用以下公式計算風險評分：NRLs risk-Score=（0.315 1×ZFPM2-ASI）+（0.863 7×MKLN1-AS）+（0.462 1×LINC0111 6）+（1.055 6×AP003390.1）。

圖3 識別HCC中NRLs預后特征：A為單變量Cox回歸分析森林圖；B為NRLs差異表達熱圖；C為NRLs的LASSO系數曲線；D為LASSO算法中變量選擇的10倍交叉驗證

在訓練集、測試集和完整的集合中按風險評分中位值將病人劃分為高風險組和低風險組，比較了HCC 病人的風險評分分布、生存狀態、生存時間以及NRLs 在兩組病人中的表達分布（圖4A-J）。這些結果都表明高風險組病人的預后較差。同時，按年齡（Age）、性別（Gender）、分級（Grade）、分期（Stage）進行分組，研究HCC 病人在臨床病理變量中的風險特征與預后的關系。結果顯示，除年齡>65 歲和女性外，其余的臨床病理學特征分層中低風險組的OS均顯著高于高風險組，驗證了風險評分模型的普遍適用性。

圖4 不同集合中風險模型的預后：A~C為訓練集、驗證集合完整集合的風險曲線；D~F為訓練集、驗證集和完整集合的生存狀態圖；G~I為訓練集、驗證集和完整集合中4個NRLs表達熱圖

在單變量-Cox 和多變量-Cox 回歸中，風險評分的風險比（HR）分別為1.06 和1.07（P<0.05）（圖5A、B）。NRLs 預后特征的1 年、3 年和5 年ROC 曲線下面積（AUC）分別為0.74、0.66 和0.67（圖5C）。在模型的1 年ROC 中，風險評分的AUC 為0.74，與年齡（Z=10.90，P<0.01）、性別（Z=15.25，P<0.001）、分級（Z=11.92，P<0.001）、分期（Z=14.45，P<0.001）間差異有統計學意義，表明其具有較強的預測能力（圖5D）。

圖5 驗證預后特征的預測值：A為臨床病理學因素和風險評分的單變量-Cox回歸；B為臨床病理學因素和風險評分的多變量-Cox回歸；C為基于風險評分特征的1年、3年和5年ROC曲線；D為風險評分模型的預測準確性與臨床病理特征的比較

2.3 列線圖的構造和校準根據風險評分和臨床病理學因素，開發了用于預測 1 年、3 年和5 年OS的列線圖（圖6A）。校準曲線顯示出列線圖的預測能力與實際觀察值具有良好的一致性（圖6B）。

圖6 模型的列線圖和校準曲線：A為用于預測HCC病人 1 年、3 年和 5 年OS 的列線圖；B為列線圖的校準曲線

2.4 GSEA 富集分析為了研究高低風險組間富集途徑的差異，采用GSEA 進行KEGG 途徑富集分析。結果顯示，NOD 樣受體信號通路、Notch 信號通路、細胞周期、RIG-I 樣受體信號通路、癌癥通路在高風險組中顯著富集；脂肪酸代謝、初級膽汁酸代謝、藥物代謝細胞色素P450、PPAR 信號通路、過氧化物酶體在低風險組中顯著富集（圖7）。

圖7 高低風險組間GSEA富集分析

2.5 腫瘤微環境分析首先，我們研究了風險評分與不同平臺腫瘤浸潤免疫細胞間的相關性。結果顯示，更多的免疫細胞與風險評分存在正相關關系（表1）。同時，免疫微環境分析結果顯示，高風險病人的ImmuneScore（Z=?3.53，P<0.001）和ESTIMATE?Score（Z=?2.52，P<0.05）明顯高于低風險病人（圖8A~8C）。為了進一步探索風險評分與免疫狀態之間的相關性，ssGSEA結果顯示（圖8D~8E），7種免疫細胞在高低風險組中存在差異包括：抗原呈遞過程中像成熟樹突狀細胞（aDCs）、樹突狀細胞（DCs）、未成熟樹突狀細胞（iDCs）、漿細胞樣樹突狀細胞（pDCs）、Th1 細胞（Th1 cells）、Th2 細胞（Th2 cells）、調節性T 細胞（Treg）且均在高風險組高表達。10 種免疫相關途徑在高低風險組間差異有統計學意義且均在高風險組中上調，包括抗原提呈細胞共抑制（APC co inhibition）、抗原提呈細胞共刺激（APC co stimulation）、細胞因子-細胞因子受體（CCR）、檢查點（check point）、人類白細胞抗原（HLA）、主要組織相容性復合體Ⅰ（MHC class Ⅰ）、副炎癥（Parainfla?mation）、T 細胞共抑制（T cell co inhibition）、T 細胞共刺激（T cell co stimulation）、Ⅱ型干擾素反應（typeⅡ IFN Reponse）。此外，通過比較不同風險組之間的免疫檢查點激活，我們發現幾乎所有免疫檢查點在高風險組中表達更多的活性（圖8F）。

表1 不同平臺免疫細胞與風險評分的相關性

圖8 腫瘤微環境分析：A~C為高低風險組間 StromalScore、ImmuneScore 和 ESTIMATEScore 的箱線圖；D~E為風險組中免疫細胞和免疫功能的ssGSEA評分；F為風險組間免疫檢查點表達的差異

2.6 基于NRLs 特征篩選疾病潛在藥物我們計算了20 種常規化療藥物在高低風險組中的IC50值，以評估HCC 病人對化療的反應。結果顯示，兩個風險組對16種化療藥物的反應差異有統計學意義（P<0.05）。其中AICAR、Axitinib、Cyclopamine、DMOG、Docetaxel、Gefitinib、Lapatinib、Nilotinib、Vinblastine藥物低風險組的IC50 值較低，表明這些藥物對低風險組病人的療效較好；而Cisplatin、Doxorubicin、Gemcitabine、Imatinib、Paclitaxel、Salubrinal、Tipifarn?ib 藥物在高風險組中IC50 值較低表明對高風險組病人療效較好。其余4個化療藥物的IC50值在高低風險組間不存在差異。

3 討論

HCC 有較高的發病率（28.33/10 萬）和死亡率（24.69/10 萬）對人類健康構成極大威脅，為HCC 病人確定特異和準確的預后特征對于改善預后至關重要。先前的研究表明，壞死性凋亡在許多類型癌癥的遷移和侵襲中發揮重要的作用被認為是消除癌細胞有希望的療法［13］。lncRNA 已被證明廣泛參與癌癥相關細胞通路，在預后和診斷方面具有良好的預測能力［14］。在本研究中構建的風險評分模型由4個NRLs組成，根據風險評分中位值將病人劃分為高低風險兩組。生存分析、Cox 回歸分析、ROC 曲線均表明風險評分可作為獨立預后因子，構建的用于預測HCC 病人OS 的列線圖顯示出較好的預測性能。腫瘤微環境分析顯示高低風險組間存在差異。最后，比較了高低風險組間的藥物敏感性，共篩選出16種化療藥物的IC50值在兩組間存在差異。

在這項研究中，獲得了760 個差異表達的NRLs來探索預后功能。通過單變量Cox、LASSO 和多變量Cox 回歸分析，最終確定了4 個NRLs（ZFPM2-AS1、MKLN1-AS、LINC01116、AP003390.1）與HCC病人的OS 顯著相關，以構建NRLs 預后特征。ZF?PM2-AS1 位于染色體8q23 上，與多種癌癥的發生發展密切相關。例如，肺腺癌、甲狀腺癌和胃癌等［15］。研究［16］表明，ZFPM2-AS1 在HCC 組織中過表達與較差的總生存期相關，其可充當 miRNA海綿通過多個軸促進HCC細胞增殖，凋亡，遷移和侵襲。Yan等［17］在基于TCGA 的HCC 病人預后特征研究中ZFPM2-AS1 被確定為預后lncRNA。與此結果一致，本研究發現與正常組織相比HCC 組織中ZFPM2-AS1 表達顯著升高。MKLN1-AS 在HCC 中可起到致癌調節劑的作用。Pan 等［18］體外細胞實驗表明沉默MKLN1-AS 可抑制HCC 的HuH7 和LM3 細胞增殖、血管生成、遷移和侵襲。其可通過作為miR-654-3p的分子海綿來增加HDGF 表達，在HCC 過程中產生促癌作用，并且MKLN1-AS 的下調可抑制HCC 細胞的侵襲性表型［19］。MKLN1-AS 在HCC 中被報道可作為診斷和預后生物標志物［20］，進一步證實了本研究的可靠性。LINC01116 異常表達與多種癌癥相關，包括肺癌，胃癌，結直腸癌，膠質瘤和骨肉瘤［21］。其在HCC中作為癌基因起作用，參與EMT和免疫調節，過表達可促進腫瘤細胞增殖，誘導EMT 相關因子的表達，與HCC 細胞的增殖、細胞周期進展和遷移密切相關［22］。Jiang 等［23］研究表明，LINC01116 過表達的HCC 病人通常具有較低的OS，可作為HCC的預后生物標志物。AP003390.1 與HCC 關系的研究報道較少，在本研究中AP003390.1在HCC病人中高表達與病人的預后密切相關。

TME 中的免疫細胞浸潤通常隨著腫瘤的發生和發展而變化。在這項研究中，基于TME 評分，高風險群體中的樣本具有較高的免疫評分。進一步分析表明，aDCs、DCs、iDCs、pDCs、Th1 cells、Th2 cells、Treg 均在高風險組高表達，表明這些異常浸潤的免疫細胞可能與HCC 病人的腫瘤啟動和發展有關。同時，APC co inhibition、APC co stimulation、CCR、Check point、HLA、MHC class Ⅰ、Parainflama?tion、T cell co inhibition、T cell co stimulation、Type ⅡIFN Reponse 在高風險組病人中具有較高的活性，表明HCC 的發生和發展可能與異常激活的免疫相關途徑存在關聯。隨后，分析了常見免疫檢查點表達與NRLs 預后特征之間的相關性。研究［24］指出，免疫檢查點基因的表達水平與免疫療法的療效高度相關。結果顯示，與低風險組相比，高風險HCC 病人中的大多數免疫檢查點表達都升高，表明風險較高的HCC 病人接受免疫治療時療效更好。最后，比較了20種常見化療藥物在高低風險組中的敏感性。結果顯示，9 種化療藥物對低風險組病人療效較好，7 種化療藥物對高風險組病人療效較好，為指導臨床醫生為HCC 病人選擇合適的抗癌藥物提供理論基礎。

這項研究仍然存在一些局限性。首先，單因素-Cox 結果顯示AC138356.1 是唯一對HCC 發病具有保護作用的NRLs，然而僅考慮了統計學意義將其排除在模型之外。其次，分析結果未經體內和體外驗證，生物學功能需要進一步闡明。因此，在接下來的工作中將收集和擴展臨床樣本，并嘗試通過外部實驗來驗證該模型的準確性。

（本文圖2見插圖8-7，圖3，4見插圖8-8，圖5，6見插圖8-9，圖7，8見插圖8-10）