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

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

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

1 材料與方法

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

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

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

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

2 結(jié)果

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

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

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

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

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

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

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

2.3 構(gòu)建風險評估模型 使用單變量Cox 回歸分析來檢驗訓練集中鐵死亡相關(guān)lncRNA 表達水平與OS之間的關(guān)系，從而篩選出符合條件的候選鐵死亡相關(guān)lncRNA（P<0.01）。此后，對高維數(shù)據(jù)進行LASSO 回歸，利用R 軟件中的“GLIMT”軟件包選擇意義最顯著的預后因素。最終得出9個與預后有關(guān)的鐵死亡相關(guān)lncRNA，分別為：Z97989.1、AL445524.1、AL391807.1、AL161431.1、AC024075.1、AC246787.2、FLG-AS1、LINC01352、AC090559.1，依據(jù)模型公式計算各lncRNA風險評分，風險評分值=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表達水平。

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

圖3 風險評分模型預后評估Fig.3 Risk scoring model prognostic assessment

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

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

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

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

圖6 預后風險評分與免疫細胞浸潤Fig.6 Prognostic risk score and immune cell infiltration

3 討論

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

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

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