李淑楨，昝雨欣

湖北醫(yī)藥學(xué)院生物醫(yī)藥研究院（湖北十堰 442000）

喉癌屬于頭頸鱗癌的一種，多表現(xiàn)為進(jìn)食嗆咳和呼吸窘迫，對(duì)人類(lèi)健康構(gòu)成嚴(yán)重威脅[1]。喉癌的治療手段包括手術(shù)、放療和化療[2]。晚期喉癌生存率下降，癌細(xì)胞轉(zhuǎn)移會(huì)導(dǎo)致嚴(yán)重的健康損害。近年來(lái)，腫瘤風(fēng)險(xiǎn)預(yù)測(cè)評(píng)分模型作為一種非侵入性方法，已被證實(shí)能夠有效預(yù)測(cè)患者的存活率，并在臨床實(shí)踐中得到逐步推廣[3-4]。因此，創(chuàng)建喉癌預(yù)測(cè)模型有助于在臨床治療前評(píng)估患者的生存率，從而進(jìn)行精準(zhǔn)治療。

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

1 資料與方法

1.1 銅死亡相關(guān)lncRNA的數(shù)據(jù)處理與鑒定

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

1.2 預(yù)后銅死亡相關(guān)lncRNA模型構(gòu)建

銅死亡相關(guān)基因列表來(lái)自銅死亡相關(guān)研究[9-11]?；? 000 倍交叉驗(yàn)證對(duì)顯著表達(dá)的lncRNA 進(jìn)行Lasso Cox 回歸，以識(shí)別銅死亡相關(guān)的lncRNA?；诙嘧兞緾ox 回歸分析確定最佳預(yù)后lncRNA，然后計(jì)算風(fēng)險(xiǎn)評(píng)分[12]。將每個(gè)lncRNA 的表達(dá)水平（Exp 值）及其對(duì)應(yīng)的回歸系數(shù)（β）相乘，然后將所有乘積項(xiàng)相加，得到綜合風(fēng)險(xiǎn)評(píng)分，即風(fēng)險(xiǎn)評(píng)分=Exp lncRNA1×β lncRNA1+Exp lncRNA2×β lncRNA2+Exp lncRNA3×β lncRNA3 +…+Exp lncRNAn×β lncRNAn（n為樣本數(shù)）。

1.3 風(fēng)險(xiǎn)評(píng)分分組能力評(píng)估

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

1.4 主成分分析

使用R 語(yǔ)言中l(wèi)imma 和scatterplot3d 軟件包進(jìn)行主成分分析（principal component analysis，PCA），由此產(chǎn)生的輸出以矩陣形式展示了高風(fēng)險(xiǎn)和低風(fēng)險(xiǎn)群體的空間分布，用于評(píng)估模型的準(zhǔn)確性。

1.5 腫瘤突變負(fù)荷分析

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

1.6 免疫相關(guān)功能分析

采用R 語(yǔ)言中l(wèi)imma、GSVA 包分析喉癌患者免疫相關(guān)功能的差異，pheatmap 包用于可視化結(jié)果。從腫瘤免疫功能障礙和排除（tumor immune dysfunction and exclusion，TIDE）分析網(wǎng)站（http://tide.dfci.harvard.edu/）探索風(fēng)險(xiǎn)曲線(xiàn)和免疫功能障礙和排斥的關(guān)系[13]。

1.7 喉癌干性指數(shù)

通過(guò)UCSC Xena 數(shù)據(jù)庫(kù)（https://xenabrowser.net/） 下載TCGA 泛癌兩類(lèi)干性評(píng)分?jǐn)?shù)據(jù)（StemnessScores_DNAmeth_20170210.ts 和StemnessScores_RNAexp_20170127.2.tsv），使用limma、ggExtra、ggplot2 軟件包進(jìn)行可視化分析。

1.8 功能富集

下載MSigDB 數(shù)據(jù)庫(kù)中c2.cp.kegg.v7.5.1.symbols.gmt 和c5.go.v7.5.1.symbols.gmt 文件作為功能基因集，進(jìn)行基因本體論（gene ontology，GO）和京都基因與基因組百科全書(shū)（Kyoto Encyclopedia of Genes and Genomes，KEGG）富集，隨機(jī)組合次數(shù)設(shè)為1 000 次、錯(cuò)誤發(fā)現(xiàn)率＜0.05、標(biāo)準(zhǔn)化富集得分＞1、P＜0.05。

1.9 統(tǒng)計(jì)分析

采用R 4.0.3 軟件對(duì)數(shù)據(jù)進(jìn)行分析。通過(guò)ROC 曲線(xiàn)、C 指數(shù)和PCA 曲線(xiàn)評(píng)估模型的預(yù)測(cè)性能，Spearman 相關(guān)檢驗(yàn)分析基因表達(dá)間的相關(guān)性，生存曲線(xiàn)比較采用Kaplan-Meier 法和Log rank 檢驗(yàn)。以P＜0.05 為差異有統(tǒng)計(jì)學(xué)意義。

2 結(jié)果

2.1 銅死亡相關(guān)lncRNA的鑒定與預(yù)后模型構(gòu)建

以|r|＞0.4 和P＜0.001 為標(biāo)準(zhǔn)，從16 773 個(gè)lncRNA 和19 個(gè)銅死亡相關(guān)基因中鑒定出570 個(gè)銅死亡相關(guān)lncRNA，兩者之間的共表達(dá)關(guān)系見(jiàn)圖1。單變量Cox回歸分析篩選出17個(gè)lncRNA（圖2-A），Lasso Cox 回歸分析確定與銅死亡相關(guān)的lncRNA（圖2-B）。多變量Cox 分析確定5 個(gè)lncRNA作為獨(dú)立預(yù)后因素，并建立風(fēng)險(xiǎn)評(píng)估模型，其風(fēng)險(xiǎn)評(píng)分=（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個(gè)lncRNA在不同風(fēng)險(xiǎn)亞組中的相對(duì)表達(dá)情況：AL158166.2、HOXB-AS3、AC121764.1、MAP3K2-DT為高風(fēng)險(xiǎn)lncRNA，而AC091729.2為低風(fēng)險(xiǎn)lncRNA（圖2-C）。同時(shí)相關(guān)熱圖還顯示了銅死亡相關(guān)基因PDHB、LIPT2、GCSH與5個(gè)lncRNA較為密切（圖2-D）。

圖1 銅死亡相關(guān)基因和喉癌相關(guān)lncRNA?；鶊DFigure 1.Sankey map of genes associated with cuproptosis and lncRNAs associated with laryngeal cancer

圖2 銅死亡相關(guān)lncRNA預(yù)后模型Figure 2.Prognosis model of lncRNA associated with cuproptosis

2.2 不同風(fēng)險(xiǎn)組的臨床病理特征分層預(yù)后

風(fēng)險(xiǎn)曲線(xiàn)反映了喉癌患者的風(fēng)險(xiǎn)評(píng)分與生存狀態(tài)之間的關(guān)系，生存分析結(jié)果表明，低風(fēng)險(xiǎn)組患者中位OS 和PFS 約為6 個(gè)月，而高風(fēng)險(xiǎn)組患者中位OS 和PFS 約為3 個(gè)月（P＜0.05），高風(fēng)險(xiǎn)組患者的OS 和PFS 明顯短于低風(fēng)險(xiǎn)組患者（圖3-A、圖3-B）。進(jìn)一步分析兩組在臨床分期方面的差異，結(jié)果顯示，低風(fēng)險(xiǎn)組患者在不同臨床分期下的中位生存期約為6 個(gè)月，高風(fēng)險(xiǎn)組患者的生存可能性?xún)H2～3 個(gè)月（P＜0.05），不同分級(jí)（1～2 和3～4 期）高風(fēng)險(xiǎn)組患者生存期均低于低風(fēng)險(xiǎn)組（圖3-C、圖3-D），表明風(fēng)險(xiǎn)評(píng)分較高的患者預(yù)后較差。

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

2.3 風(fēng)險(xiǎn)模型的預(yù)后能力驗(yàn)證

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

圖4 風(fēng)險(xiǎn)模型的預(yù)后能力驗(yàn)證Figure 4.Validation of prognostic proficiency of risk models

2.4 主成分分析

PCA 分析結(jié)果表明，在喉癌患者中，無(wú)論是所有基因（圖5-A）、銅死亡相關(guān)基因（圖5-B），還是銅死亡相關(guān)lncRNA（圖5-C），均未能顯著區(qū)分不同高低風(fēng)險(xiǎn)患者。相對(duì)而言，基于風(fēng)險(xiǎn)模型的分析揭示了兩種風(fēng)險(xiǎn)聚類(lèi)間重疊較小，表明模型能夠有效區(qū)分低風(fēng)險(xiǎn)和高風(fēng)險(xiǎn)患者群體（圖5-D），從而驗(yàn)證了該風(fēng)險(xiǎn)模型在構(gòu)建過(guò)程中的可靠性。

圖5 主成分分析結(jié)果Figure 5.Results of PCA

2.5 高低風(fēng)險(xiǎn)人群的腫瘤突變負(fù)荷及相關(guān)預(yù)后

使用maftools 算法觀察高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組的突變，發(fā)現(xiàn)對(duì)于大多數(shù)基因，高風(fēng)險(xiǎn)組的突變頻率高于低風(fēng)險(xiǎn)組（TP53：低風(fēng)險(xiǎn)，75%；高風(fēng)險(xiǎn)：84%；TTN：低風(fēng)險(xiǎn)，52%；高風(fēng)險(xiǎn)，54%；CSMD3：低風(fēng)險(xiǎn)，33%；高風(fēng)險(xiǎn)，37%），見(jiàn)圖6-A、圖6-B。此外，高低風(fēng)險(xiǎn)組之間的TMB差異并不顯著（P=0.93），見(jiàn)圖6-C。高突變負(fù)荷組和低突變負(fù)荷組的存活率無(wú)顯著差異（P=0.21），見(jiàn)圖6-D。綜合生存分析結(jié)合風(fēng)險(xiǎn)評(píng)分顯示，四組患者的生存時(shí)間差異有統(tǒng)計(jì)學(xué)意義（P＜0.05），高突變和低風(fēng)險(xiǎn)患者的生存時(shí)間最長(zhǎng)，低突變和高風(fēng)險(xiǎn)患者的生存時(shí)間較短（圖6-E）。

圖6 高低風(fēng)險(xiǎn)人群的腫瘤突變負(fù)荷及相關(guān)預(yù)后Figure 6.TMB and associated outcomes in high - and low-risk populations

2.6 腫瘤干細(xì)胞和免疫分析

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

圖7 免疫相關(guān)功能分析Figure 7.Analysis of immune-related functions

2.7 功能富集分析

基因集富集分析（gene set enrichment analysis，GESA）結(jié)果顯示，在高風(fēng)險(xiǎn)組中，GO 富集分析的前五位顯著富集項(xiàng)包括適應(yīng)性免疫應(yīng)答、骨骼系統(tǒng)發(fā)育、外部包膜結(jié)構(gòu)、高爾基體內(nèi)腔及免疫球蛋白復(fù)合體（圖8-A）。同時(shí)，KEGG 通路富集分析指出，前五位顯著富集的途徑包括細(xì)胞外基質(zhì)（extracellular matrix，ECM）受體相互作用、刺猬（Hedgehog）信號(hào)通路、WNT 信號(hào)通路、基底細(xì)胞瘤、黑素瘤（圖8-B）。

圖8 風(fēng)險(xiǎn)曲線(xiàn)在喉癌中功能富集分析Figure 8.Functional enrichment analysis of risk curves in laryngeal carcinoma

3 討論

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

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

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

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

綜上，本研究中為喉癌構(gòu)建了銅死亡相關(guān)的lncRNA 模型，并分析了風(fēng)險(xiǎn)評(píng)分與TMB、免疫治療和腫瘤干性指數(shù)之間的關(guān)聯(lián)，從而為喉癌患者的預(yù)后預(yù)測(cè)和臨床治療提供參考。