m6A相關(guān)lncRNA 預(yù)測肺鱗癌患者預(yù)后和免疫反應(yīng)分析①

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

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

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

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

1 資料與方法

1.1 獲取LUSC RNA-seq轉(zhuǎn)錄組及臨床數(shù)據(jù)

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

1.2 m6A相關(guān)lncRNA的鑒定

從TCGA-LUSC中提取30個m6A調(diào)控因子的RNA表達(dá)矩陣,包括“編寫器”(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相關(guān)分析對m6A相關(guān)的lncRNA進(jìn)行篩選(|R|>0.4,P<0.001)。

1.3 預(yù)后m6A相關(guān)lncRNA特征的建立和驗證

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

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

1.4 功能分析

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

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

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

1.6 統(tǒng)計學(xué)方法

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

2 結(jié)果

2.1 LUSC中m6A相關(guān)lncRNA的鑒定

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

圖1 m6A相關(guān)LncRNAs和m6A相關(guān)預(yù)后LncRNAs的鑒定a.共表達(dá)網(wǎng)絡(luò); b.森林圖。

2.2 m6A相關(guān)lncRNA特征的建立和驗證

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

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

圖2 m6A相關(guān)的lncRNA的預(yù)后特征a-b.訓(xùn)練集和測試集中的Kaplan-Meier生存曲線;c-d.訓(xùn)練集和測試集中,LUSC患者的風(fēng)險評分和生存狀態(tài)的分布; e-f .訓(xùn)練集和測試集中,與m6A相關(guān)的預(yù)后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相關(guān)LncRNA特征的分層和獨立預(yù)后分析

為了更好地評估已建立的m6A相關(guān)預(yù)后特征是否可以作為一個獨立的預(yù)后因素,基于臨床病理特征對訓(xùn)練集和測試集進(jìn)行了單變量和多變量Cox回歸分析?；陬A(yù)后m6A相關(guān)lncRNA特征的風(fēng)險評分可以預(yù)測訓(xùn)練和測試集中LUSC的預(yù)后(單變量: 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生存分析結(jié)果顯示,低危組在年齡(≤65歲、>65歲)、性別(男、女)、TNM期(I～I(xiàn)I期、III～I(xiàn)V期)、T期(T1～2)、N期(N0、N1～3)、M期(M0)(P<0.05;圖4e-l)方面均有更好的預(yù)后。這些結(jié)果表明,高危評分的模型與預(yù)后不良相關(guān),可作為LUSC的獨立預(yù)后因素。

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

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

2.4 m6A預(yù)后相關(guān)lncRNA特征的功能分析

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

2.5 腫瘤突變負(fù)荷(TMB)和腫瘤浸潤性免疫細(xì)胞(TICs)分析結(jié)果

我們分析了高低風(fēng)險組中每個LUSC患者的突變譜。在總數(shù)中,顯著突變的前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),且風(fēng)險評分與TMB呈負(fù)相關(guān)(圖6d;R=-0.13,P<0.05),提示低危組可能表現(xiàn)出更好的免疫應(yīng)答。

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

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

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

3 討論

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

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

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

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