聯(lián)合免疫分型構(gòu)建肺腺癌鐵死亡相關(guān)長鏈非編碼RNA風(fēng)險評分模型①

于海洋 潘躍銀 （中國科學(xué)技術(shù)大學(xué)附屬第一醫(yī)院腫瘤化療科，合肥 230001）

肺癌是我國最常見同時也是致死率最高的惡 性腫瘤，五年生存率僅為4%～17%，其中又以肺腺癌（lung adenocarcinoma，LUAD）所占比例最高，因此尋找LUAD 的治療新靶點是當(dāng)務(wù)之急［1］。鐵死亡是一種新型的非凋亡細(xì)胞死亡模式，其主要機(jī)制為鐵依賴性的脂質(zhì)過氧化損傷，這種導(dǎo)致線粒體功能障礙和細(xì)胞內(nèi)的毒性脂質(zhì)過氧化的機(jī)制在抑制癌癥生長和發(fā)展過程中起關(guān)鍵作用［2-3］。lncRNA 是長度大于200 個核苷酸，不具有或僅有少許蛋白編碼能力的RNA 序列［4］。lncRNA 還參與各種生物調(diào)控過程，包括與腫瘤的發(fā)生、發(fā)展和轉(zhuǎn)移有關(guān)的過程，如LINC00618 在人類白血病細(xì)胞中低表達(dá)，其通過提高脂質(zhì)ROS 和鐵水平，同時抑制SLC7A11 表達(dá)來促進(jìn)鐵死亡發(fā)生［5-6］。然而，鐵死亡相關(guān)lncRNA 在LUAD 中的意義尚未明確，本研究依據(jù)LUAD 患者鐵死亡與lncRNA 之間的相關(guān)性，結(jié)合免疫分型構(gòu)建新型風(fēng)險評分模型以評估LUAD患者的預(yù)后。

1 材料與方法

1.1 數(shù)據(jù)收集與處理 TCGA 數(shù)據(jù)庫（https：//portal.gdc.cancer.gov/）下載LUAD 患者的轉(zhuǎn)錄組數(shù)據(jù)和臨床信息，包括59 例非腫瘤組織和535 例腫瘤樣本。Perl 軟件對數(shù)據(jù)進(jìn)行預(yù)處理以及l(fā)ncRNA 表達(dá)信息和相關(guān)臨床病理參數(shù)的提取。下載GEO 數(shù)據(jù)庫中GSE31210 數(shù)據(jù)集中的表達(dá)數(shù)據(jù)及臨床信息以供外部驗證。采用R軟件通過共表達(dá)分析鐵死亡相關(guān)基因與lncRNA 之間的關(guān)系，以篩選出與鐵死亡相關(guān)的lncRNA。結(jié)合臨床數(shù)據(jù)，利用單因素Cox 回歸分析法及Kaplan-Meier（K-M）生存分析法篩選出與LUAD預(yù)后有關(guān)的鐵死亡相關(guān)lncRNA。

1.2 一致性聚類與免疫分型 基于得到的預(yù)后相關(guān)lncRNA 的綜合表達(dá)信息，使用R 軟件中的“ConensusClusterPlus”包對腫瘤樣本的不同亞組進(jìn)行了一致性聚類鑒定，用K-M 生存分析法繪制亞組間的生存曲線進(jìn)行判定，并繪制與臨床特征相關(guān)性的熱圖。運用CIBERSORT 軟件（https：//cibersort.stanford.edu/）對所下載的LUAD 患者mRNA 數(shù)據(jù)處理后可得到腫瘤微環(huán)境內(nèi)浸潤的免疫細(xì)胞，通過轉(zhuǎn)錄組的表達(dá)數(shù)據(jù)可反映免疫細(xì)胞的浸潤情況［6-7］；結(jié)合患者風(fēng)險評估值與浸潤免疫細(xì)胞進(jìn)行Pearson 相關(guān)分析，找出與患者腫瘤浸潤免疫細(xì)胞表達(dá)情況的相關(guān)性［8］。運用R軟件中“estimate”包評估每個患者腫瘤樣本免疫細(xì)胞和基質(zhì)細(xì)胞的浸潤水平（即免疫評分和基質(zhì)評分）。

1.3 建立風(fēng)險評估模型 篩選出預(yù)后信息不完整的31 例樣本后，LUAD 腫瘤樣本（n=504）被隨機(jī)分為訓(xùn)練集（n=252）和驗證集（n=252），隨后使用單變量Cox 回歸分析來檢驗訓(xùn)練集中鐵死亡相關(guān)lncRNA 表達(dá)水平與總生存期（OS）之間的關(guān)系。篩選出符合條件的候選鐵死亡相關(guān)lncRNA（P<0.01）。此后，對數(shù)據(jù)進(jìn)行LASSO 回歸，利用R 中的“GLIMT”包篩選意義最顯著的預(yù)后因素，得出與預(yù)后有關(guān)的鐵死亡相關(guān)lncRNA，并計算其風(fēng)險評分。最終得到構(gòu)成風(fēng)險評分模型的評分公式如下所示：風(fēng)險評分∑i N= 1（Ei*Ci），其中N 表示構(gòu)建風(fēng)險評分模型的鐵死亡相關(guān)lncRNA 數(shù)目，Ci 表示鐵死亡相關(guān)lncRNA 系數(shù)，Ei 表示鐵死亡相關(guān)lncRNA 的表達(dá)水平。

1.4 模型評價與驗證 依據(jù)風(fēng)險評分模型計算訓(xùn)練集中每個患者的風(fēng)險得分。隨后依據(jù)訓(xùn)練集風(fēng)險評分值中位數(shù)，將訓(xùn)練集與驗證集患者各分為高風(fēng)險組（風(fēng)險值>中位數(shù)）和低風(fēng)險組（風(fēng)險值≤中位數(shù)）。使用R 軟件繪制高低風(fēng)險組K-M 生存曲線以對比兩組患者生存差異，繪制ROC 曲線并用曲線下面積（AUC）來評估模型預(yù)測患者預(yù)后的準(zhǔn)確性。然后分別對訓(xùn)練集和驗證集患者的風(fēng)險評分與其他臨床因素（年齡、性別、分期）進(jìn)行單因素及多因素Cox 回歸驗證，以考察其是否為獨立預(yù)后因子。隨后通過臨床分組模型進(jìn)行進(jìn)一步驗證。最后利用來自GEO數(shù)據(jù)庫中的樣本進(jìn)行外部驗證。

2 結(jié)果

2.1 預(yù)后有關(guān)鐵死亡相關(guān)LncRNA 篩選 篩選出預(yù)后信息不完整的31例樣本后，本研究匯總了來自TCGA 數(shù)據(jù)庫中504 例腫瘤樣本、59 例非腫瘤樣本以及GEO 數(shù)據(jù)中226個LUAD 患者的臨床基線資料（表1）。利用R 語言共表達(dá)分析FerrDb 數(shù)據(jù)庫（http：//www.zhounan.org/ferrdb/）259 個鐵死亡相關(guān)基因與lncRNA，以篩選出與鐵死亡相關(guān)的lncRNA。同時結(jié)合LUAD 患者的臨床信息，單因素Cox 回歸分析最終篩選11 個與LUAD 預(yù)后有關(guān)的鐵死亡相關(guān)lncRNA（表2）。

表1 患者的臨床基線數(shù)據(jù)Tab.1 Clinical baseline data of patients

表2 單變量Cox 回歸鑒定11個鐵死亡相關(guān)lncRNA 的P 值和風(fēng)險比Tab.2 P values and hazard ratios of 11 ferroptosis-related lncRNA identified by univariate Cox regression

2.2 一致性聚類與免疫分析 基于R 軟件中的“ConensusClusterPlus”包對504 個腫瘤樣本的不同亞組進(jìn)行一致性聚類分析（圖1A～C）。結(jié)果輸出為K=2～9 個子組，當(dāng)K=2 時為最優(yōu)的劃分方式。根據(jù)最穩(wěn)定的K 值，所有患者被成功地分為兩個亞組。其中，亞組1（Cluster1）代表低水平的基因表達(dá)，而亞組2（Cluster2）代表較高水平的基因表達(dá)。差異表達(dá)基因的總體生存分析表明，Cluster1 的生存期明顯延長（P<0.001，圖1D）。

圖1 一致性聚類分析與鑒定Fig.1 Consistent cluster analysis and identification

腫瘤微環(huán)境差異分析顯示（圖2A～C），在Cluster1 中免疫細(xì)胞和基質(zhì)細(xì)胞含量明顯高于Cluster2，Cluster1 腫瘤純度顯著低于Cluster2，這一結(jié)果與上述生存分析結(jié)果也相互佐證。通過CiberSort 方法評估了一致性聚類中每個患者22 種不同免疫細(xì)胞的相對比例，并對這兩個風(fēng)險組產(chǎn)生的CiberSort 輸出進(jìn)行對比匯總（圖2D），包括靜息樹突狀細(xì)胞、靜息肥大細(xì)胞、靜息CD4+記憶T 細(xì)胞、幼稚B 細(xì)胞、中性粒細(xì)胞、CD8+T 細(xì)胞、輔助濾泡T 細(xì)胞在內(nèi)的免疫細(xì)胞在不同的風(fēng)險組中富集。結(jié)果表明靜息樹突狀細(xì)胞、靜息肥大細(xì)胞、靜息CD4+記憶T 細(xì)胞在Cluster1 中顯著高表達(dá)，幼稚B 細(xì)胞、中性粒細(xì)胞、CD8+T 細(xì)胞和輔助濾泡T 細(xì)胞在Cluster2 中明顯高表達(dá)（圖2E）。

圖2 腫瘤微環(huán)境與免疫細(xì)胞浸潤Fig.2 Tumor microenvironment and immune cell infiltration

2.3 構(gòu)建風(fēng)險評估模型 使用單變量Cox 回歸分析來檢驗訓(xùn)練集中鐵死亡相關(guān)lncRNA 表達(dá)水平與OS之間的關(guān)系，從而篩選出符合條件的候選鐵死亡相關(guān)lncRNA（P<0.01）。此后，對高維數(shù)據(jù)進(jìn)行LASSO 回歸，利用R 軟件中的“GLIMT”軟件包選擇意義最顯著的預(yù)后因素。最終得出9個與預(yù)后有關(guān)的鐵死亡相關(guān)lncRNA，分別為：Z97989.1、AL445524.1、AL391807.1、AL161431.1、AC024075.1、AC246787.2、FLG-AS1、LINC01352、AC090559.1，依據(jù)模型公式計算各lncRNA風(fēng)險評分，風(fēng)險評分值=EZ97989.1 ×（-0.940 467 387 287 318）+EAL445 524.1×0.000 879 959 066 730 395+EAL391807.1×（-0.657 561 513 468 988）+ EAL161431.1×0.00 416 410 759 392 408+EAC024 075.1×（-0.0 111 194 380 084 745）+EAC246787.2×（-0.137 501 955 859 741）+EFLGAS1×1.94901939831152+ELINC01352×（-1.83 848 358 454 151）+EAC090559.1×（-0.159 643 403 496 863），其中E表示鐵死亡相關(guān)lncRNA表達(dá)水平。

2.4 風(fēng)險評估模型的評價與驗證 利用“survival”包對高低風(fēng)險組患者的生存進(jìn)行差異分析，并繪制風(fēng)險評分曲線及生存狀態(tài)關(guān)系圖（圖3A）。隨著風(fēng)險評分的升高，相較于低風(fēng)險組而言，高風(fēng)險組患者的病死人數(shù)明顯增多。生存分析表明，與高風(fēng)險組相比，低風(fēng)險組患者OS 顯著延長（圖3C）。為了進(jìn)一步評估此模型的特異度和靈敏度，通過“survivalROC”程序包來繪制模型的ROC 曲線。ROC 曲線顯示：訓(xùn)練集中1 年和5 年AUC 值分別為0.708、0.758（圖3E），表明該模型具有較好的預(yù)測效能。通過單因素以及多因素Cox回歸分析探討了風(fēng)險評分模型是否是LUAD 患者預(yù)后的獨立預(yù)測因子。在除其他臨床病理特征（如腫瘤分期、年齡以及性別）的干擾后，此模型具有較好的預(yù)測效能。

圖3 風(fēng)險評分模型預(yù)后評估Fig.3 Risk scoring model prognostic assessment

依據(jù)同一風(fēng)險評分值，將驗證集患者分為高、低風(fēng)險兩組。與從訓(xùn)練集得出的結(jié)果一致，隨著風(fēng)險評分的升高，相較于低風(fēng)險組，高風(fēng)險組患者的預(yù)后明顯更差（圖3B）；且高風(fēng)險評分患者組的中位OS 與低風(fēng)險評分組比要更短（P=0.001，圖3D）。驗證集中1 年和5 年AUC 值為0.659、0.636（圖3E、F），這表明本預(yù)測模型在驗證集也具有良好的預(yù)測效能。對臨床因素和風(fēng)險評分值與預(yù)后關(guān)系的分析表明，與LUAD 患者的腫瘤分期、N 分期、免疫分組均與風(fēng)險評分相關(guān)（P<0.05），且風(fēng)險評分隨著患者腫瘤淋巴結(jié)分期增加，這說明此風(fēng)險評分預(yù)后模型與LUAD 淋巴結(jié)轉(zhuǎn)移存在一定關(guān)聯(lián)。R 軟件“survminer”包根據(jù)臨床分組從而對模型進(jìn)行進(jìn)一步驗證，如圖4 所示，以年齡、性別、分期及T、N、M 分期作為分組依據(jù)，本風(fēng)險預(yù)測模型同樣適用于臨床分組后的模型。依據(jù)同一風(fēng)險評分值，將來自GEO數(shù)據(jù)庫226 例LUAD 患者分為高、低風(fēng)險兩組，結(jié)果表明高風(fēng)險組患者的預(yù)后明顯差于低風(fēng)險組（P<0.001，圖5A），1 年、3 年及5 年ROC 曲線AUC 值分別為0.61、0.69、0.77（圖5B），表明本預(yù)測模型在外部驗證集同樣具有良好的預(yù)測效能。

圖4 風(fēng)險評分值與臨床因素Fig.4 Risk score values and clinical factors

圖5 GSE31210數(shù)據(jù)集外部驗證Fig.5 External validation of GSE31210 dataset

對免疫細(xì)胞表達(dá)情況和免疫預(yù)后RiskScore 風(fēng)險值進(jìn)行Pearson 相關(guān)性檢驗，由圖6 可知，巨噬細(xì)胞M0、靜息肥大細(xì)胞、活化肥大細(xì)胞、靜息樹突狀細(xì)胞和活化樹突狀細(xì)胞等腫瘤浸潤免疫細(xì)胞與患者免疫預(yù)后RiskScore風(fēng)險值密切相關(guān)（P<0.05），活化樹突狀細(xì)胞、巨噬細(xì)胞M0、活化肥大細(xì)胞、活化NK細(xì)胞和輔助濾泡T 細(xì)胞呈正相關(guān)（R>0）；記憶B 細(xì)胞、靜息樹突狀細(xì)胞、巨噬細(xì)胞M2、靜息肥大細(xì)胞、單核細(xì)胞和靜息CD4+記憶T細(xì)胞呈負(fù)相關(guān)（R<0）。

圖6 預(yù)后風(fēng)險評分與免疫細(xì)胞浸潤Fig.6 Prognostic risk score and immune cell infiltration

3 討論

精準(zhǔn)的療效預(yù)測和預(yù)后判斷方法是臨床防治LUAD 的重中之重。目前以臨床特征或血清分子標(biāo)志物（CEA、CA125 等）作為預(yù)后判斷的方式較為局限，高通量測序技術(shù)的廣泛應(yīng)用則為癌癥診斷、治療和預(yù)后研究提供了越來越多的測序數(shù)據(jù)［9-10］。lncRNA 指可通過影響鄰近的基因或其他染色體上的遠(yuǎn)距離基因來達(dá)到調(diào)控編碼基因表達(dá)目的并且長度超過200個核苷酸的非編碼RNA 序列［1］。既往研究提示在LUAD 中有很多異常表達(dá)的lncRNA（如lncRNA MIR31HG、lncRNA H19、lncRNA HOXCAS3 等），作為LUAD 患者特異性的生物學(xué)標(biāo)志物在診斷、預(yù)測預(yù)后等方面起一定作用［11-13］。CAI等［14］研究表明，在裸鼠肺癌腫瘤模型中，過表達(dá)轉(zhuǎn)鐵蛋白受體1（transferrin receptor 1，TFR1）可促進(jìn)肺癌細(xì)胞對鐵吸收，從而促進(jìn)腫瘤細(xì)胞的增殖。已有研究表明通過構(gòu)建鐵死亡相關(guān)基因風(fēng)險評分模型可以預(yù)測肺腺癌患者的預(yù)后［15］；此外，LUO 等［16］通過構(gòu)建鐵死亡相關(guān)lncRNA 預(yù)后模型可在一定程度上預(yù)測頭頸部鱗狀細(xì)胞癌患者的預(yù)后。然而，目前為止尚未有預(yù)測肺腺癌患者的鐵死亡lncRNA 風(fēng)險評分模型。因此，構(gòu)建相關(guān)風(fēng)險評分模型有助于LUAD 患者的預(yù)后評估。

本研究經(jīng)過對TCGA 數(shù)據(jù)庫的LUAD 樣本的匯總和篩選，最終得到504例符合標(biāo)準(zhǔn)的樣本數(shù)據(jù)，并將其隨機(jī)分為訓(xùn)練集與驗證集，最終得出由9 個鐵死亡相關(guān)的lncRNA 構(gòu)成的風(fēng)險評分模型，結(jié)果表明，AL445524.1、AL161431.1 和FLG-AS1 的高表達(dá)均與患者的 OS 呈正相關(guān)；Z97989.1、AL391807.1、AC024075.1、AC246787.2、LINC01352和AC090559.1 高表達(dá)與較差的預(yù)后相關(guān)。先前已有研究發(fā)現(xiàn)AL161431.1 及LINC01352 分別與肝癌和子宮內(nèi)膜癌的發(fā)生發(fā)展聯(lián)系緊密，且通過構(gòu)建lncRNA 預(yù)后模型可在一定程度上預(yù)測患者的預(yù)后［17-19］。無論是單因素分析還是多因素Cox 回歸分析均表明該風(fēng)險預(yù)后模型是影響LUAD 患者預(yù)后的獨立危險因素（P<0.001）。該風(fēng)險模型對訓(xùn)練集和驗證集中1年和5年AUC 值分別為0.708、0.758；驗證集中1 年和5 年AUC 值分別為0.659、0.636。通過一致性聚類使用多種不同的聚類方法，從而找到一種比每種單獨的方法更合適的聚類方法對腫瘤進(jìn)行分組，以便有助于為每個患者制定個性化的治療方法。根據(jù)基因表達(dá)水平，使用R 軟件的“ConensusClusterPlus”包 將 數(shù) 據(jù) 分 為Cluster1 與Cluster2。生存分析表明，Cluster1 的生存時間明顯延長，提示生存時間與鐵死亡相關(guān)基因的綜合表達(dá)水平有關(guān)。

然而，本研究仍存在一定程度的不足，本研究是一項基于TCGA 數(shù)據(jù)庫的回顧性研究，有可能產(chǎn)生一定范圍內(nèi)的偏差。此外，由于缺少合適的數(shù)據(jù)集，本研究尚未在其他數(shù)據(jù)庫進(jìn)行深入的驗證；且本研究是一項基于高通量測序結(jié)果的回顧性研究，并未進(jìn)行進(jìn)一步的功能研究和體內(nèi)外實驗。因此需要對本風(fēng)險預(yù)后模型進(jìn)行更加深入的探索和后續(xù)研究。