銅死亡相關長鏈非編碼RNA預測喉癌的預后和調控機制

李淑楨，昝雨欣

湖北醫藥學院生物醫藥研究院（湖北十堰 442000）

喉癌屬于頭頸鱗癌的一種，多表現為進食嗆咳和呼吸窘迫，對人類健康構成嚴重威脅[1]。喉癌的治療手段包括手術、放療和化療[2]。晚期喉癌生存率下降，癌細胞轉移會導致嚴重的健康損害。近年來，腫瘤風險預測評分模型作為一種非侵入性方法，已被證實能夠有效預測患者的存活率，并在臨床實踐中得到逐步推廣[3-4]。因此，創建喉癌預測模型有助于在臨床治療前評估患者的生存率，從而進行精準治療。

當機體受到刺激后，細胞內銅的積累會導致Fe-S 簇蛋白的不穩定，這是一種獨特細胞死亡類型[5]。既往臨床試驗表明，與正常組織相比，各種惡性腫瘤中的銅含量更高，且銅積累與腫瘤的增殖和生長有關[6]。因此，研究銅死亡在癌癥中的作用具有較大的臨床應用潛力。長鏈非編碼RNA（long non-coding RNA，lncRNA）在多種腫瘤的發生與進展中扮演著重要角色，并在調控腫瘤細胞的增殖、遷移和侵襲等方面發揮作用[7]。有研究顯示，外泌體lncRNALINC02191通過靶向miR-204-5p/RAB22A軸及調節磷脂酰肌醇3-激酶 /蛋白激酶B/哺乳動物雷帕霉素靶蛋白信號通路（PAM 信號通路）促進喉鱗癌進展[8]。本研究應用生物信息學分析喉癌患者的銅死亡相關lncRNA 及其生物學和預后，探討其在預測喉癌患者預后中的生物學功能及作用。

1 資料與方法

1.1 銅死亡相關lncRNA的數據處理與鑒定

從癌癥基因組圖譜（The Cancer Genome Atlas，TCGA） 數據庫（https://portal.gdc.cancer.gov/） 下載基因表達的RNA-Seq 數據和喉癌的相應臨床和突變數據，其中包含12 例正常個體和111 例喉癌患者。根據TCGA 的基因注釋獲得RNA-Seq 數據，以區分lncRNA。為了鑒定與銅死亡相關的lncRNA，使用limma、dplyr、ggalluvial 和ggplot2軟件包進行Pearson 相關性分析，確定銅死亡相關基因表達水平與lncRNA 表達之間的關系。使用的閾值為|r|＞0.4 和P＜0.001。

1.2 預后銅死亡相關lncRNA模型構建

銅死亡相關基因列表來自銅死亡相關研究[9-11]。基于1 000 倍交叉驗證對顯著表達的lncRNA 進行Lasso Cox 回歸，以識別銅死亡相關的lncRNA。基于多變量Cox 回歸分析確定最佳預后lncRNA，然后計算風險評分[12]。將每個lncRNA 的表達水平（Exp 值）及其對應的回歸系數（β）相乘，然后將所有乘積項相加，得到綜合風險評分，即風險評分=Exp lncRNA1×β lncRNA1+Exp lncRNA2×β lncRNA2+Exp lncRNA3×β lncRNA3 +…+Exp lncRNAn×β lncRNAn（n為樣本數）。

1.3 風險評分分組能力評估

根據風險評分的中值將患者分為高、低風險組，使用pheatmap 包繪制基于風險評分的基因表達熱圖；采用Kaplan-Meier 生存分析，使用Survival 包繪制各組喉癌患者的總生存期（overall survival，OS）、無進展生存期（progression-free survival，PFS）和患者腫瘤分級下生存曲線；使用Survival 和ROC 包計算風險曲線的時間依賴性受試者工作特征（receiver operating characteristic，ROC）曲線下的1 年、3 年和5 年面積（area under the curve，AUC）；使用R 的rms 和pec 包生成一致性指數（concordance index，C-index）曲線。

1.4 主成分分析

使用R 語言中limma 和scatterplot3d 軟件包進行主成分分析（principal component analysis，PCA），由此產生的輸出以矩陣形式展示了高風險和低風險群體的空間分布，用于評估模型的準確性。

1.5 腫瘤突變負荷分析

利用maftools 軟件包，通過讀取高風險和低風險樣本的基因突變文件，對突變頻率前15 個基因可視化，檢測兩組間的遺傳變異。此外，使用ggpubr、survival 和survminer 包比較高風險組和低風險組樣本之間的腫瘤突變負荷（tumor mutational burden，TMB）和存活率。

1.6 免疫相關功能分析

采用R 語言中limma、GSVA 包分析喉癌患者免疫相關功能的差異，pheatmap 包用于可視化結果。從腫瘤免疫功能障礙和排除（tumor immune dysfunction and exclusion，TIDE）分析網站（http://tide.dfci.harvard.edu/）探索風險曲線和免疫功能障礙和排斥的關系[13]。

1.7 喉癌干性指數

通過UCSC Xena 數據庫（https://xenabrowser.net/） 下載TCGA 泛癌兩類干性評分數據（StemnessScores_DNAmeth_20170210.ts 和StemnessScores_RNAexp_20170127.2.tsv），使用limma、ggExtra、ggplot2 軟件包進行可視化分析。

1.8 功能富集

下載MSigDB 數據庫中c2.cp.kegg.v7.5.1.symbols.gmt 和c5.go.v7.5.1.symbols.gmt 文件作為功能基因集，進行基因本體論（gene ontology，GO）和京都基因與基因組百科全書（Kyoto Encyclopedia of Genes and Genomes，KEGG）富集，隨機組合次數設為1 000 次、錯誤發現率＜0.05、標準化富集得分＞1、P＜0.05。

1.9 統計分析

采用R 4.0.3 軟件對數據進行分析。通過ROC 曲線、C 指數和PCA 曲線評估模型的預測性能，Spearman 相關檢驗分析基因表達間的相關性，生存曲線比較采用Kaplan-Meier 法和Log rank 檢驗。以P＜0.05 為差異有統計學意義。

2 結果

2.1 銅死亡相關lncRNA的鑒定與預后模型構建

以|r|＞0.4 和P＜0.001 為標準，從16 773 個lncRNA 和19 個銅死亡相關基因中鑒定出570 個銅死亡相關lncRNA，兩者之間的共表達關系見圖1。單變量Cox回歸分析篩選出17個lncRNA（圖2-A），Lasso Cox 回歸分析確定與銅死亡相關的lncRNA（圖2-B）。多變量Cox 分析確定5 個lncRNA作為獨立預后因素，并建立風險評估模型，其風險評分=（AL158166.2×4.594 224 454 551 14）+（HOXB-AS3×2.83871554342233）-（AC091729.2×5.560 000 728 450 9）+（AC121764.1×4.888 326 291 011 31）+（MAP3K2-DT×2.092 706 823 954 72）。熱圖顯示了5個lncRNA在不同風險亞組中的相對表達情況：AL158166.2、HOXB-AS3、AC121764.1、MAP3K2-DT為高風險lncRNA，而AC091729.2為低風險lncRNA（圖2-C）。同時相關熱圖還顯示了銅死亡相關基因PDHB、LIPT2、GCSH與5個lncRNA較為密切（圖2-D）。

圖1 銅死亡相關基因和喉癌相關lncRNA桑基圖Figure 1.Sankey map of genes associated with cuproptosis and lncRNAs associated with laryngeal cancer

圖2 銅死亡相關lncRNA預后模型Figure 2.Prognosis model of lncRNA associated with cuproptosis

2.2 不同風險組的臨床病理特征分層預后

風險曲線反映了喉癌患者的風險評分與生存狀態之間的關系，生存分析結果表明，低風險組患者中位OS 和PFS 約為6 個月，而高風險組患者中位OS 和PFS 約為3 個月（P＜0.05），高風險組患者的OS 和PFS 明顯短于低風險組患者（圖3-A、圖3-B）。進一步分析兩組在臨床分期方面的差異，結果顯示，低風險組患者在不同臨床分期下的中位生存期約為6 個月，高風險組患者的生存可能性僅2～3 個月（P＜0.05），不同分級（1～2 和3～4 期）高風險組患者生存期均低于低風險組（圖3-C、圖3-D），表明風險評分較高的患者預后較差。

圖3 喉癌患者的Kaplan-Meier生存分析Figure 3.Kaplan-Meier survival analysis of patients with laryngeal cancer

2.3 風險模型的預后能力驗證

采用單因素和多因素生存回歸分析，探討喉癌風險模型的獨立預后模型。單變量Cox 回歸結果顯示，患者性別風險比（hazard ratio，HR）為0.292，與OS 獨立相關（P＜0.05）。多變量Cox回歸結果顯示，性別[HR=0.236（0.097，0.572），P＜0.05]與OS 獨立相關，表明風險模型是喉癌患者的獨立預后因素（圖4-A、圖4-B）。采用ROC 曲線評估風險評分的預測準確性，在風險模型中，1 年、3 年和5 年OS 分別為0.809、0.699和0.745（圖4-C），風險評分的AUC 為0.809，優于年齡（0.516）、性別（0.377）、分級（0.498）和分期（0.484）（圖4-D）。該風險模型的C 指數值高于其他臨床特征（圖4-E），表明該模型具有良好的預測能力。

圖4 風險模型的預后能力驗證Figure 4.Validation of prognostic proficiency of risk models

2.4 主成分分析

PCA 分析結果表明，在喉癌患者中，無論是所有基因（圖5-A）、銅死亡相關基因（圖5-B），還是銅死亡相關lncRNA（圖5-C），均未能顯著區分不同高低風險患者。相對而言，基于風險模型的分析揭示了兩種風險聚類間重疊較小，表明模型能夠有效區分低風險和高風險患者群體（圖5-D），從而驗證了該風險模型在構建過程中的可靠性。

圖5 主成分分析結果Figure 5.Results of PCA

2.5 高低風險人群的腫瘤突變負荷及相關預后

使用maftools 算法觀察高風險組和低風險組的突變，發現對于大多數基因，高風險組的突變頻率高于低風險組（TP53：低風險，75%；高風險：84%；TTN：低風險，52%；高風險，54%；CSMD3：低風險，33%；高風險，37%），見圖6-A、圖6-B。此外，高低風險組之間的TMB差異并不顯著（P=0.93），見圖6-C。高突變負荷組和低突變負荷組的存活率無顯著差異（P=0.21），見圖6-D。綜合生存分析結合風險評分顯示，四組患者的生存時間差異有統計學意義（P＜0.05），高突變和低風險患者的生存時間最長，低突變和高風險患者的生存時間較短（圖6-E）。

圖6 高低風險人群的腫瘤突變負荷及相關預后Figure 6.TMB and associated outcomes in high - and low-risk populations

2.6 腫瘤干細胞和免疫分析

風險模型僅對基質細胞評分具有正相關作用（r=0.24，P=0.013）（圖7-B），對免疫細胞評分基因則無（r=0.05，P=0.6）（圖7-A、圖7-C）。TIDE 算法是近年開發的用于確定腫瘤免疫檢查點治療效果的工具[13]。使用TIDE 算法進一步研究高風險組和低風險組患者對免疫治療敏感性的差異，發現除了MSI Expr Sig 低風險組的TIDE評分與高風險組相比有差異外，其他免疫指標均無統計學意義（圖7-D）。干性指數是評估腫瘤細胞與干細胞相似性的指標，反映了腫瘤細胞中的干細胞特征，如自我更新能力和多向分化潛能。通過UCSC Xena 數據庫探索高低風險組的與干細胞相關性，結果顯示，風險曲線與RNAss呈負相關（r=-0.21，P=0.025）（圖7-E），與DNAss 無明顯相關性（P=0.069）（圖7-F）。

圖7 免疫相關功能分析Figure 7.Analysis of immune-related functions

2.7 功能富集分析

基因集富集分析（gene set enrichment analysis，GESA）結果顯示，在高風險組中，GO 富集分析的前五位顯著富集項包括適應性免疫應答、骨骼系統發育、外部包膜結構、高爾基體內腔及免疫球蛋白復合體（圖8-A）。同時，KEGG 通路富集分析指出，前五位顯著富集的途徑包括細胞外基質（extracellular matrix，ECM）受體相互作用、刺猬（Hedgehog）信號通路、WNT 信號通路、基底細胞瘤、黑素瘤（圖8-B）。

圖8 風險曲線在喉癌中功能富集分析Figure 8.Functional enrichment analysis of risk curves in laryngeal carcinoma

3 討論

喉癌作為呼吸道常見惡性腫瘤，過去幾年在治療策略上雖有所進步，但由人類乳頭瘤病毒感染率檢測方法差異引起的研究間結果不一致性，仍對臨床治療構成挑戰[14]。鑒于此，開發一個精確的喉癌風險評估模型對于患者的預后判斷和生存期預測具有重要的臨床價值。lncRNA 是一類非蛋白質編碼RNA，通過轉錄調控等影響腫瘤發生發展[15]。研究發現，RASSF8-AS1通過靶向miR-664b-3p/TLE1軸在喉鱗狀細胞癌中顯示出抑制作用[16]。lncRNA-NEAT1通過海綿狀miR-577/miR-1224-5p阻斷細胞凋亡途徑[17]。這些發現表明lncRNA 有潛力成為新的喉癌生物標志物和治療靶點。銅死亡是與線粒體三羧酸循環密切相關的一種新型細胞死亡方式，盡管已有針對銅離子的小分子抗癌藥物依來昔洛爾的臨床試驗，但其治療效果并不理想[10,18]。因此，銅死亡與lncRNA 在喉癌中的相互作用及其調節機制亟待進一步研究。

本研究初篩17 個預后銅死亡相關lncRNA。通過多變量Cox 回歸分析，進一步確認了AL158166.2、HOXB-AS3、AC121764.1、MAP3K2-DT和AC091729.2作為構建預后模型的關鍵lncRNA。ROC、生存曲線、C-index 分析結果表明，5 種銅死亡相關lncRNA 的預后模型可準確區分高風險和低風險人群，能可靠地預測喉癌患者的預后。有研究表明，MAP3K2-DT表達與乳腺癌患者中的TGFβ1表達呈正相關[19]。而HOXB-AS3在多種癌癥中發揮重要作用，包括縮短頭頸鱗癌患者的生存期，并通過miR-378a-3p在卵巢癌中上調細胞外酸化率，促進腫瘤細胞的增殖、侵襲和遷移，以及通過與DNMT1的結合影響肝癌細胞的增殖和凋亡[20-22]。

腫瘤干細胞被廣泛認為是導致腫瘤復發的關鍵因素，而腫瘤干性指數的開發旨在精確識別并量化腫瘤中具有此類干細胞特性的細胞群，為癌癥干細胞的檢測與監控提供了一種有效的量化工具。Zhang 等研究發現，與親本細胞相比，喉癌干細胞樣細胞在小鼠體內具有更強的腫瘤形成能力和化療抵抗性[23]。GSEA 通路分析表明，銅死亡相關的lncRNA 可能與WNT 信號傳導途徑的發展有關。Tang 等報道了DGCR5通過激活人喉癌細胞中的WNT 信號通路，促使miR-506形成海綿狀結構，誘導腫瘤干細胞樣特性[24]，與本研究所得干性指數相符。但目前仍缺乏銅死亡與喉癌干細胞相關的基礎研究，未來或可以此為研究方向。

有研究表明，喉癌中腫瘤-宿主邊界處的基質腫瘤浸潤淋巴細胞與免疫系統激活相關基因的表達和預后顯著相關[25]。然而，在本研究風險模型中，高低風險曲線與免疫細胞及其功能之間未發現顯著關聯。TMB 作為免疫檢查點阻斷治療的預測生物標志物，在臨床上已被證實在耐輻射頭頸鱗癌組中較對照組更高，提示免疫治療可能是更合適的治療選擇[26-27]。在高風險患者中，TP53和TTN的表達水平上升，其中TP53是多種癌癥中常見的突變基因，影響多種癌癥的發展及患者生存率，而TTN的參與則與多種癌癥的發展相關，TTN-AS1的過表達與多種癌癥的不良預后相關[28-29]。這表明TP53和TTN可能成為喉癌免疫治療的新靶點。lncRNA 與免疫相關功能之間的確切關系尚未完全闡明，未來應評估喉癌患者中免疫細胞浸潤的預后價值，并進一步探討免疫細胞在銅死亡介導的喉癌患者靶向治療中的作用。本研究存在一些局限性：首先，研究結果基于公共生物數據庫，缺乏臨床樣本的直接驗證；其次，對于預后模型中的免疫作用及其機制仍需深入研究，這將是未來研究的重點。

綜上，本研究中為喉癌構建了銅死亡相關的lncRNA 模型，并分析了風險評分與TMB、免疫治療和腫瘤干性指數之間的關聯，從而為喉癌患者的預后預測和臨床治療提供參考。