m6A相關lncRNA 預測肺鱗癌患者預后和免疫反應分析①

毛宇澤,楊成鵬,姜騰蛟,駱 鵬,朱曉峰

(佳木斯大學附屬第一醫院心胸外科,黑龍江佳木斯 154003)

肺癌是世界上最常見的惡性腫瘤,非小細胞肺癌(NSCLC)約占所有肺癌的85%,包括腺癌、鱗癌和大細胞肺癌等[1]。肺鱗狀細胞癌(LUSC)是一種起源于支氣管鱗狀上皮且與吸煙密切相關的惡性腫瘤,約占NSCLC病例的20%～30% 。與DNA和蛋白質相類似,RNA經歷轉錄后修飾,其中包括N1-甲基化(m1A)、5-甲基化(m5C)、N6-甲基化(m6A)和7-甲基化(m7G)。m6A修飾是真核生物中調控基因表達的最豐富的RNA修飾。這些m6A修飾包括甲基轉移酶、甲基結合蛋白和去甲基化酶,分別被命名為“編寫器”、“閱讀器”和“橡皮擦”[2～4]。研究表明,m6A調節因子通過調節腫瘤的生長和進展以及治療反應,在NSCLC中發揮潛在作用[5]。m6A閱讀器YTH n6-甲基化RNA結合蛋白1(YTHDF1)抑制NSCLC中的腫瘤生長和集落形成,而YTHDF1的低表達水平與臨床結果差相關,并通過上調抗氧化系統促進癌細胞對順鉑的耐藥性[6]。因此,m6A的修飾在NSCLC中起著至關重要的作用。

本研究在LUSC中篩選了m6A相關的lncRNA建立并驗證了預后特征。進一步區分了高、低危組,并分析了免疫細胞在LUSC中的分布及免疫療效。有助于預測LUSC的預后并區分LUSC患者是否受益于免疫治療。

1 資料與方法

1.1 獲取LUSC RNA-seq轉錄組及臨床數據

從TCGA數據庫(https://cancergenome.nih.gov/)和TCGA-LUSC(n=551)下載LUSC的標準化RNA-seq數據和臨床信息(年齡、性別、TNM分期、生存時間和狀態),其中包括501個腫瘤樣本和49個正常樣本。我們在分析中排除了缺失生存信息的樣本。

1.2 m6A相關lncRNA的鑒定

從TCGA-LUSC中提取30個m6A調控因子的RNA表達矩陣,包括“編寫器”(METTL3,METTL14,METTL16,METTL5,WTAP,VIRMA,RBM15,RBM15B,ZC3H13,CBLL1,和ZCCHC4),“閱讀器”(YTHDF1,YTHDF2,YTHDF3,YTHDC1,YTHDC2,HNRNPA2B1,HNRNPC,FMF31,IGF2BP1,IGF2BP2,IGF2BP3EPLAVL1,G3BP1G3BP2,PRRC2A,RBMX),和“橡皮擦”(FTO和ALKBH5)。此外,采用Pearson相關分析對m6A相關的lncRNA進行篩選(|R|>0.4,P<0.001)。

1.3 預后m6A相關lncRNA特征的建立和驗證

采用單因素Cox回歸分析篩選預后后與m6A相關的lncRNA(P<0.05)。采用R包“glmnet”進行10倍交叉驗證,以進一步縮小范圍。最后,利用λmin(λmin = -4.5)鑒定了7個lncRNA并建立了與m6A相關的lncRNA預后特征,使用以下公式計算了風險評分:Risk score = coef (lncRNA1) × expr (lncRNA1) + coef (lncRNA2) ×expr (lncRNA2) +……+ coef (lncRNA) × expr (lncRNAn), 其中coef為系數,expr(lncRNAn)為lncRNAs的表達量。

將數據集以7:3的比例被劃分為訓練集和測試集。根據每個患者的風險評分,使用R軟件進行生存分析、風險曲線和繪制熱圖。使用ROC評估了模型對預測預后的敏感性。AUC使用R包“cesivivalROC”計算。采用單因素和多因素Cox回歸分析來評估該模型作為一個獨立的預后因素的潛力。采用分層分析來評估該模型對不同臨床病理特征的預測能力。

1.4 功能分析

使用GSEA軟件揭示高危組和低危組之間差異表達的功能通路。從Molecular Signatures Database數據庫下載并分析兩個廣泛應用的基因集(h.all.v7.2.symbols.gmt [cancer hallmarks]和c7.all. v7.4.symbols.gmt [immunologic signatures])。為了獲得標準化富集分數,進行1000次基因集排列。其中具有|NES|>1, nominalP<0.05,and q<0.25的通路被認為是顯著富集的。

1.5 腫瘤突變負荷(TMB)和腫瘤浸潤性免疫細胞(TICs)分析

采用CIBERSORT算法(http://cibersort.stanford.edu)在LUSC樣本中獲得了22個不同的TICs。采用TIDE(http://tide.dfci.harvard.edu/)來預測各亞組中免疫治療的療效。此外,從TCGA數據庫中下載了LUSC患者資料,并使用“maftools”軟件包進行分析和可視化。最后,計算每個患者的TMB。

1.6 統計學方法

采用R版本4.1.0和Perl編程語言(版本5.34.0)進行。基于預后m6A相關lncRNA的表達情況,采用K-M生存曲線分析和log-rank檢驗來比較不同亞組間的OS。此外,進行單因素和多因素Cox回歸分析來評估特征在預測OS方面的獨立預后價值。采用Wilcoxon檢驗來確定兩組間變量的差異。P<0.05為差異有統計學意義。

2 結果

2.1 LUSC中m6A相關lncRNA的鑒定

為識別LUSC中m6A相關的lncRNA,從癌癥基因組圖譜(TCGA)-LUSC數據集中提取30個m6A調控因子和14086個lncRNA表達矩陣進行皮爾遜相關分析,篩選出373個m6A相關的lncRNA(|R|>0.4,P<0.001)。與m6A相關的lncRNA共表達網絡見圖1a。為了篩選LUSC中重要的lncRNA,將臨床預后和lncRNA的表達整合到模型中。使用單變量Cox回歸分析篩選了9個與預后相關的lncRNA,見圖1b。其中 AP001189.3和HORMAD2-AS1(HR>1、P<0.05)與預后不良相關,而AL122125.1、AP001469.3、AC138035.1、AC002、AC002467.1、AC020658.4、AC243919.2、PRC1-AS1與預后良好相關(HR <1、P<0.05)。

圖1 m6A相關LncRNAs和m6A相關預后LncRNAs的鑒定a.共表達網絡; b.森林圖。

2.2 m6A相關lncRNA特征的建立和驗證

基于上述與m6A相關的lncRNA,完成LUSC的lncRNA特征的建立,對9個預測的m6A相關的lncRNA進行了Lasso回歸分析。確認了AL122125.1、AP001469.1、AC138035.1、AC002467.1、AC243919.1、AC243919.2、AP001189.3、HORMAD1-As1構建風險模型。為了探究lncRNA簽名的作用,將整個集合(n=495)分為訓練集(n=347)和測試集(N = 148),并對這些集合進行了Kaplan-Meier生存分析。結果顯示,高危組的預后比低危組差(訓練組和測試組的P=0.003和P=0.002;圖2a,b)。此外,兩組患者的風險曲線顯示,隨著風險評分的增加,LUSC的預后進一步惡化(圖2c,d)。

構建一個受試者工作特征曲線(ROC)來評估我們的模型的敏感性。訓練組1年、3年和5年總生存率(OS)的曲線下面積(AUC)分別為0.581、0.581和0.608。在測試集中,1-、3-和5年OS的AUC分別為0.655、0.655和0.555。這些結果表明,已建立的LUSC的預后m6A相關特征具有一般的敏感性(圖2e,f)。K-M生存曲線顯示,AP001189.3和HORMAD2-AS1表達水平較高與預后不良相關,AC002467.1、AC243919.2、AC138035.1、AL122125.1和AP001469.3表達水平較低,與LUSC良好預后相關(圖3a～g)。

圖2 m6A相關的lncRNA的預后特征a-b.訓練集和測試集中的Kaplan-Meier生存曲線;c-d.訓練集和測試集中,LUSC患者的風險評分和生存狀態的分布; e-f .訓練集和測試集中,與m6A相關的預后LncRNAs特征的ROC曲線。

圖3 Kaplan-Meier生存曲線a.AP001189.3; b.HORMAD2-AS1; c.AC002467.1; d.AC243919.2; e.AC138035.1; f.AL122125.1; g.AP001469.3(P<0.05)。

2.3 m6A相關LncRNA特征的分層和獨立預后分析

為了更好地評估已建立的m6A相關預后特征是否可以作為一個獨立的預后因素,基于臨床病理特征對訓練集和測試集進行了單變量和多變量Cox回歸分析。基于預后m6A相關lncRNA特征的風險評分可以預測訓練和測試集中LUSC的預后(單變量: HR=3.437、P<0.001和HR=6.960、P=0.005;多變量: HR=3.475、P<0.001和HR=10.386、P=0.003;圖4a～d)。K～M生存分析結果顯示,低危組在年齡(≤65歲、>65歲)、性別(男、女)、TNM期(I～II期、III～IV期)、T期(T1～2)、N期(N0、N1～3)、M期(M0)(P<0.05;圖4e-l)方面均有更好的預后。這些結果表明,高危評分的模型與預后不良相關,可作為LUSC的獨立預后因素。

圖4 m6A相關預后特征的獨立預后和分層分析a-b.訓練集中m6A相關預后特征的單變量和多變量回歸分析;c-d.測試集中m6A相關預后特征的單變量和多變量回歸分析;e-l .高危組和低危組的Kaplan-Meier生存曲線(P-value<0.05)。

圖5 高低風險組的功能富集分析a.在低風險組中具有豐富的腫瘤特征:細胞周期;b-d.高危組中富集的腫瘤特征: MAPK信號通路、細胞凋亡、JAK/STAT信號通路;e-h.免疫相關信號通路富集(e:低危組,f-h:高危組)。

2.4 m6A預后相關lncRNA特征的功能分析

為了檢驗高風險組和低風險組之間的功能差異,我們進行了基因集富集分析(GSEA),發現有幾種不同的腫瘤特征在兩個風險組之間有顯著差異富集。細胞周期通路在低危組中富集(圖6a),而MAPK、凋亡和JAK/STAT信號通路在高危組中富集(圖6b～d)。值得注意的是,兩個危險組之間許多差異表達的免疫相關信號通路與T細胞、B細胞和單核細胞相關(圖6e中的低風險組,圖6f～h中的高危組)。因此,研究結果表明,預后后與m6A相關的lncRNA與腫瘤發生和免疫通路相關。

2.5 腫瘤突變負荷(TMB)和腫瘤浸潤性免疫細胞(TICs)分析結果

我們分析了高低風險組中每個LUSC患者的突變譜。在總數中,顯著突變的前20個基因為TP53、TTN、TTN、CSMD3、MUC16、MUC16、RYR2、LRP1B、USH2A、SYNE1、ZFHX4、FAM135B、KMT2D、XIRP2、SPTA1、CDH10、NAV3、PCDH15、PAPPA2、RYR3、DNAH5、PKHD1(圖6a低危組,圖6b高危組)。低危組的TMB顯著高于高危組(圖6c),且風險評分與TMB呈負相關(圖6d;R=-0.13,P<0.05),提示低危組可能表現出更好的免疫應答。

圖6 TMBa.低風險組突變基因;b.高危組突變基因;c.低危組的TMB高于高危組(P<0.05);d.風險評分與TMB呈負相關(R=-0.13,P=0.0032)。

比較兩組之間的免疫細胞的差異,發現低風險組比高危組中幼稚B細胞,濾泡輔助T細胞,和激活NK細胞更豐富,而CD4記憶休息T細胞,M2巨噬細胞和中性粒細胞相對較少(圖7)。然后我們進行了腫瘤免疫功能障礙和排斥(TIDE)分析,探討兩個亞組之間的差異,發現高危組的TIDE評分高于低危組,提示高危組的腫瘤免疫逃逸和免疫治療反應較低(圖7b)。此外,我們還進行了免疫細胞與風險評分之間的相關性分析。我們發現M2巨噬細胞、中性粒細胞和CD4記憶靜息T細胞與風險評分呈正相關,而活化的NK細胞、濾泡輔助T細胞和幼稚B細胞與風險評分呈負相關(P<0.05)。這些結果表明,與m6A相關的預后lncRNA與腫瘤免疫相關,可能成為新的免疫治療的生物標志物。

圖7 高、低危組的免疫學特征a.TICs差異;b.箱線圖。

3 討論

肺癌是一種高度異質性的疾病,由于不同的組織學和分子特征,不同亞型的肺癌的治療效果也有所不同[7]。晚期肺鱗癌盡管其治療方法有限,但免疫治療已成為LUSC的標準治療方法。然而,只有一小部分表達高水平免疫檢查點如PD-L1的LUSC患者從免疫治療中獲益[8]。因此,探索其他生物標志物來識別可能受益于免疫治療的LUSC患者是至關重要的。

N6-甲基化(m6A)是一種在真核生物中調控基因表達RNA修飾[9]。有研究表明,m6A調控因子的異常表達通過在腫瘤生長、增殖、侵襲和轉移中發揮作用,促進腫瘤的發展和進展[10]。m6A相關基因在腫瘤中起著關鍵作用,但在LUSC中的功能與作用知之甚少。

在本研究中,建立并驗證了一個基于LUSC中與m6A相關的7個lncRNA的預后特征。在高低風險組中確定了特異性的癌癥和免疫相關通路,并研究了預后m6A相關lncRNA與腫瘤免疫之間的關系。M1巨噬細胞具有促炎和吞噬腫瘤細胞的作用,而M2巨噬細胞具有抗炎作用并促進腫瘤增殖[11]。結果顯示,高危組中M2巨噬細胞的數量高于低危組,提示高危組患者有更高的腫瘤進展風險。此外,TMB被廣泛稱為基于免疫檢查點抑制劑(ICI)治療應答的替代或補充生物標志物[12],高TMB患者可能比低TMB患者從ICI治療中獲益更多[13]。在我們的研究中,低危組的TMB評分高于高危組,TIDE評分較低,高危組相反,說明低危組對免疫治療的反應可能高于高危組。

綜上,本研究是對LUSC中預后m6A相關lncRNA的綜合分析,建立并驗證了一個與m6A相關的lncRNA預后特征,可用于區分受益于免疫治療的LUSC患者。本研究也加深了對m6A修飾在腫瘤中重要作用的理解,并希望對未來研究有指導作用。