基于m7G相關(guān)miRNA的肝癌風(fēng)險(xiǎn)預(yù)后模型構(gòu)建

邱康,朱謙,金晨,孫傳東

(青島大學(xué),山東 青島 266071 1 醫(yī)學(xué)部; 2 附屬醫(yī)院肝膽胰外科; 3 附屬醫(yī)院小兒外科)

肝癌是世界上是第七常見的癌癥類型,5年生存率僅有18%,是致死率第二高的腫瘤,僅次于胰腺癌[1]。隨著環(huán)境中致癌因素的增加,肝癌在我國(guó)的發(fā)病率逐年上升[2]。盡管目前的臨床管理已經(jīng)極大地改善了肝癌病人的生存率,但由于轉(zhuǎn)移率高,肝癌病人的一般預(yù)后仍然極差[3]。因此,了解肝癌的發(fā)病機(jī)制,找到一種評(píng)估病人預(yù)后以及指導(dǎo)臨床治療的方法十分重要。microRNA(miRNA)作為短的單鏈RNA分子可以通過翻譯抑制調(diào)節(jié)蛋白質(zhì)編碼基因[4],肝組織miRNA在肝癌中表達(dá)失調(diào),參與了肝癌的發(fā)生和發(fā)展[5]。肝癌細(xì)胞通過多種方式影響miRNA的功能,包括染色體改變、聚腺苷基化修飾及轉(zhuǎn)錄因子的異常表達(dá)等[6]。N7-甲基鳥苷(m7G)具有特定修飾依賴性化學(xué)性質(zhì),普遍存在于哺乳生物中并且修飾特定轉(zhuǎn)錄位置。m7G存在于真核生物mRNA的5’端帽中,且在轉(zhuǎn)運(yùn)RNA(tRNA)、核糖體RNA(rRNA)、小核RNA(snRNA)和小核仁RNA(snoRNA)等其他種類的RNA中也發(fā)現(xiàn)了類似情況[7]。已有研究表明,介導(dǎo)內(nèi)部m7G修飾的最具特異性酶甲基轉(zhuǎn)移酶樣1(METTL1)和其輔助因子WDR4,參與調(diào)節(jié)miRNA的m7G修飾,從而影響生物發(fā)生和細(xì)胞遷移[8]。本研究基于生物信息學(xué)方法,應(yīng)用腫瘤基因組圖譜(TCGA)數(shù)據(jù)庫(kù)構(gòu)建了m7G修飾相關(guān)的miRNA的風(fēng)險(xiǎn)預(yù)后模型,并利用Estimate數(shù)據(jù)庫(kù)對(duì)免疫細(xì)胞浸潤(rùn)和免疫功能進(jìn)行分析,旨在構(gòu)建一個(gè)模型來(lái)評(píng)價(jià)肝癌病人的預(yù)后,為肝癌病人的治療和預(yù)后評(píng)價(jià)提供參考。現(xiàn)將結(jié)果報(bào)告如下。

1 材料和方法

1.1 數(shù)據(jù)下載與處理

從TCGA數(shù)據(jù)庫(kù)(https://portal.gdc.cancer.gov)中檢索并下載肝癌的miRNA表達(dá)數(shù)據(jù)和臨床特征數(shù)據(jù),其中正常組織50例,腫瘤組織375例,臨床數(shù)據(jù)371例。利用軟件Strawberry Perl 5.32.1.1將基因表達(dá)數(shù)據(jù)進(jìn)行轉(zhuǎn)換整理為表達(dá)矩陣,用于后續(xù)分析。GSEA數(shù)據(jù)庫(kù)(http://www.gsea-msigdb.org/gsea/index.jsp)下載得到m7G相關(guān)基因并進(jìn)行單因素Cox回歸篩選,得到差異基因用于后續(xù)分析。在Estimate數(shù)據(jù)庫(kù)中下載肝癌相關(guān)的373個(gè)樣本的免疫評(píng)分,用于后續(xù)免疫細(xì)胞浸潤(rùn)和免疫功能的分析。

1.2 m7G差異基因分析

使用R語(yǔ)言中的“edgeR”和“Limma”包,分析腫瘤組織和正常組織間差異基因的表達(dá)情況。差異基因的篩選標(biāo)準(zhǔn):|log2FC|≥1,P<0.05。

1.3 構(gòu)建m7G風(fēng)險(xiǎn)評(píng)分模型

利用“survival”包和“glmnet”包對(duì)差異表達(dá)miRNA進(jìn)行單因素和多因素回歸分析,得到預(yù)后相關(guān)miRNA;同時(shí)結(jié)合臨床數(shù)據(jù)計(jì)算模型的風(fēng)險(xiǎn)評(píng)分,構(gòu)建m7G風(fēng)險(xiǎn)評(píng)分模型。公式為:風(fēng)險(xiǎn)評(píng)分(Riskscore)=∑回歸系數(shù)×基因表達(dá)量。然后利用“Rtsne”包執(zhí)行PCA分析。依據(jù)風(fēng)險(xiǎn)評(píng)分,繪制Kaplan-Meier生存曲線以及ROC曲線,進(jìn)一步對(duì)該模型準(zhǔn)確性進(jìn)行評(píng)估。

1.4 風(fēng)險(xiǎn)評(píng)分與臨床因素分析

利用“survival”包構(gòu)建具有風(fēng)險(xiǎn)評(píng)估和臨床因素的Cox模型。在臨床肝癌病人中,通過單因素和多因素Cox回歸分析,評(píng)估影響其預(yù)后的獨(dú)立危險(xiǎn)因素。

1.5 功能富集分析

使用R語(yǔ)言中的“clusterprofile”“enrichplot”“ggplot2”包進(jìn)行基因本體(GO)和KEGG通路分析,以確定與差異表達(dá)的m7G相關(guān)miRNA的作用位點(diǎn)和信號(hào)通路。

1.6 風(fēng)險(xiǎn)評(píng)分模型與免疫浸潤(rùn)的相關(guān)性分析

利用“Limma”包和“edgeR”包對(duì)肝癌病人高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組進(jìn)行基因表達(dá)差異分析(篩選標(biāo)準(zhǔn):|log2FC|>1,FDR<0.05)。采用“GSVA”包中的基因集富集分析算法(ssGSEA)計(jì)算肝癌病人免疫細(xì)胞和免疫功能表達(dá)的基因豐度并繪制熱圖,進(jìn)行免疫細(xì)胞和功能的相關(guān)性分析,比較高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組的免疫細(xì)胞表達(dá)及免疫功能。

1.7 統(tǒng)計(jì)學(xué)處理

利用R語(yǔ)言軟件(4.1.3)進(jìn)行統(tǒng)計(jì)學(xué)分析。應(yīng)用單因素和多因素Cox回歸分析風(fēng)險(xiǎn)模型的獨(dú)立預(yù)后價(jià)值,Kaplan-Meier法進(jìn)行生存分析,Mann-Whitney檢驗(yàn)比較兩組免疫浸潤(rùn)和免疫通路激活情況。所有生存結(jié)局均以P值和風(fēng)險(xiǎn)比(HR)表示,置信區(qū)間(CI)為95%。P<0.05為差異有顯著性。

2 結(jié) 果

2.1 肝癌病人與m7G相關(guān)的miRNA差異基因

從GSEA數(shù)據(jù)庫(kù)下載得到34個(gè)可能與m7G相關(guān)的基因,通過單因素Cox分析肝癌病人預(yù)后數(shù)據(jù),結(jié)合相關(guān)文獻(xiàn)結(jié)果[1]篩選出12個(gè)與m7G相關(guān)的基因。利用TargetScan數(shù)據(jù)庫(kù)預(yù)測(cè)m7G調(diào)控的相關(guān)miRNA,并得到相應(yīng)表達(dá)矩陣。結(jié)合TCGA數(shù)據(jù)庫(kù)中肝癌病人的miRNA數(shù)據(jù)與得到的m7G相關(guān)miRNA進(jìn)行差異分析,根據(jù)相關(guān)篩選標(biāo)準(zhǔn)(|log2FC|>1、FDR<0.05)獲得363個(gè)miRNA差異基因并繪制火山圖(圖1A)。其中上調(diào)基因320個(gè),下調(diào)基因43個(gè),差異表達(dá)最顯著的前20個(gè)上游miRNA熱圖見圖1B。

A:火山圖;B:熱圖。

2.2 Lasso回歸篩選m7G相關(guān)miRNA

對(duì)363個(gè)差異表達(dá)miRNA與生存時(shí)間取交集進(jìn)行單因素 Cox回歸分析,結(jié)果顯示共有32個(gè)miRNA與肝癌病人的預(yù)后顯著相關(guān)。將候選的32個(gè)miRNA納入Lasso回歸構(gòu)建模型,選取誤差最低的數(shù)值,進(jìn)一步篩選出11個(gè)有意義的特征基因miRNA,分別為miR-4661-5p、miR-301a-3p、miR-760、miR-9-3p、miR-7156-5p、miR-561-5p、miR-3911、miR-513c-5p、miR-4652-3p、miR-346、miR-548aq-5p(圖2A、B)。根據(jù)臨床生存時(shí)間進(jìn)行多因素Cox回歸分析顯示,結(jié)果顯示has-miR-5p、has-miR-9-3p、has-miR-3911、hsa-miR-513c-5p、has-miR-4652-3p、has-miR-538aq-5p差異有統(tǒng)計(jì)學(xué)意義(HR=1.004～1.553,P<0.05)。見圖2C。

A:Lasso回歸系數(shù)分析;B:Lasso回歸交叉驗(yàn)證圖;C:m7G相關(guān)miRNA多因素獨(dú)立預(yù)后分析。

2.3 m7G相關(guān)miRNA的風(fēng)險(xiǎn)模型構(gòu)建

應(yīng)用上述11個(gè)與m7G相關(guān)的miRNA構(gòu)建風(fēng)險(xiǎn)模型,并且基于風(fēng)險(xiǎn)評(píng)分將肝癌病人分為高、低風(fēng)險(xiǎn)組。PCA聚類分析結(jié)果表明,二者區(qū)分明顯(圖3A)。根據(jù)風(fēng)險(xiǎn)評(píng)分構(gòu)建Kaplan-Meier曲線,經(jīng)Cox回歸分析顯示,低風(fēng)險(xiǎn)組的整體生存水平明顯高于高風(fēng)險(xiǎn)組(HR=2.23,P<0.001) (圖3B)。ROC曲線分析結(jié)果顯示,肝癌病人的模型風(fēng)險(xiǎn)評(píng)分對(duì)1、3、5年生存評(píng)估的AUC分別為0.674、0.725和0.697(圖3C)。

2.4 風(fēng)險(xiǎn)評(píng)分與臨床因素的模型構(gòu)建

為了評(píng)估風(fēng)險(xiǎn)水平模型在臨床中的預(yù)測(cè)效果,將風(fēng)險(xiǎn)水平與臨床因素中的年齡、性別、腫瘤分期(G0=T1+T2,G1=T3+T4)和風(fēng)險(xiǎn)評(píng)分納入分析。根據(jù)風(fēng)險(xiǎn)評(píng)分以及臨床因素繪制列線圖,對(duì)肝癌病人的預(yù)后進(jìn)行預(yù)測(cè),結(jié)果顯示,該預(yù)測(cè)模型的一致性指數(shù)(C-index)為0.611(0.563～0.658),肝癌病人第1、3、5年的生存率分別為78.1%、51.8%和37.1%(圖4A)。為了進(jìn)一步評(píng)估預(yù)測(cè)模型是否符合實(shí)際情況,采用Bootstrap自抽樣法進(jìn)行檢測(cè),并繪制了預(yù)測(cè)模型的校準(zhǔn)曲線。結(jié)果表明,3年和5年擬合線重合度較高,表明該預(yù)測(cè)模型具有較好的生存預(yù)測(cè)功效,其3年和5年生存期預(yù)測(cè)值與實(shí)際值相似(圖4B)。為了鑒別預(yù)測(cè)模型的區(qū)分度,繪制了預(yù)測(cè)模型的ROC曲線,結(jié)果顯示1、3、5年的AUC分別為0.613、0.701和0.673(圖4C)。

A: 風(fēng)險(xiǎn)評(píng)分和臨床因素預(yù)后相關(guān)列線圖;B:預(yù)測(cè)模型標(biāo)定曲線;C:風(fēng)險(xiǎn)評(píng)分和臨床因素模型對(duì)預(yù)后區(qū)分度的ROC曲線。

2.5 風(fēng)險(xiǎn)評(píng)分與免疫評(píng)分相關(guān)基因的功能分析

對(duì)風(fēng)險(xiǎn)模型與免疫功能的相關(guān)性分析顯示,在風(fēng)險(xiǎn)模型中的高風(fēng)險(xiǎn)和低風(fēng)險(xiǎn)組中發(fā)現(xiàn)了539個(gè)差異表達(dá)的mRNA。將來(lái)自Estimate數(shù)據(jù)庫(kù)的經(jīng)免疫評(píng)分分組的960個(gè)基因與風(fēng)險(xiǎn)模型的差異表達(dá)基因相交,產(chǎn)生183個(gè)與風(fēng)險(xiǎn)模型和免疫浸潤(rùn)相關(guān)的mRNA(圖5A)。GO分析顯示,這些基因主要參與了消化、內(nèi)分泌過程、攝食行為等生物學(xué)過程(BP);細(xì)胞成分(CC)主要富集在含膠原蛋白的細(xì)胞外基質(zhì)、分泌顆粒腔、細(xì)胞質(zhì)囊泡腔中;分子功能(MF)富集的主要是受體配體活性、信號(hào)受體激活劑活性、生長(zhǎng)因子活性等(圖5B)。KEGG分析顯示,神經(jīng)活性配體-受體相互作用、胰腺分泌、蛋白質(zhì)的消化吸收是主要參與的富集通路(圖5B)。在眾多基因中篩選出了與免疫功能和免疫細(xì)胞相關(guān)的mRNA,其中S100A8、KRT4、FCGBP與免疫浸潤(rùn)呈正相關(guān),而SOHLH1、CRISP2、CTSV則與免疫浸潤(rùn)呈負(fù)相關(guān)(圖5C)。

A:風(fēng)險(xiǎn)評(píng)分與免疫評(píng)分差異基因韋恩圖;B:GO和KEGG分析;C:差異基因與免疫浸潤(rùn)相關(guān)性熱圖。

2.6 風(fēng)險(xiǎn)模型與免疫微環(huán)境之間的關(guān)系

應(yīng)用ssGSEA分析肝癌病人風(fēng)險(xiǎn)評(píng)分與免疫細(xì)胞表達(dá)豐度以及免疫功能相關(guān)性結(jié)果顯示,高風(fēng)險(xiǎn)和低風(fēng)險(xiǎn)組中輔助性T細(xì)胞、人類白細(xì)胞抗原(HLA)和主要組織相容性復(fù)合物Ⅰ類(MHCⅠ)在腫瘤免疫微環(huán)境中顯著表達(dá)(圖6A)。腫瘤浸潤(rùn)淋巴細(xì)胞(TIL)和CD8+T細(xì)胞在肝癌病人腫瘤免疫微環(huán)境中相關(guān)性最高(r=0.870),濾泡輔助性T細(xì)胞(Tfh)和肥大細(xì)胞的相關(guān)性最低(r=-0.004),在免疫功能中T細(xì)胞共抑制途徑和免疫檢查點(diǎn)顯示最高的正相關(guān)性(r=0.930)(圖6B)。高風(fēng)險(xiǎn)組活化的樹突狀細(xì)胞(aDCs)、巨噬細(xì)胞相較于低風(fēng)險(xiǎn)組富集更明顯(P<0.05),而低風(fēng)險(xiǎn)組中B細(xì)胞、肥大細(xì)胞、中性粒細(xì)胞、NK細(xì)胞、輔助性T細(xì)胞、細(xì)胞溶解活性、亞炎癥以及Ⅰ、Ⅱ型干擾素反應(yīng)均顯著升高(P<0.05)(圖6C)。因此,肝癌高風(fēng)險(xiǎn)組病人可能更傾向于應(yīng)用抗巨噬細(xì)胞或抗aDCs細(xì)胞進(jìn)行免疫治療。

A:免疫細(xì)胞與免疫功能表達(dá)豐度熱圖;B:免疫細(xì)胞和免疫功能之間的相關(guān)性熱圖;C:不同風(fēng)險(xiǎn)組中免疫浸潤(rùn)相關(guān)表達(dá)豐度的箱式圖。

3 討 論

m7G是最普遍的RNA修飾之一,有研究表明m7G修飾顯著參與腫瘤的發(fā)生發(fā)展[9]。此外,m7G修飾不僅發(fā)生在mRNA、rRNA、tRNA上,還通過修飾miRNA介導(dǎo)RNA的代謝和功能[8]。但是針對(duì)m7G修飾miRNA是否同樣影響了肝癌的發(fā)生發(fā)展,尚未見報(bào)道。既往研究表明,miRNA在轉(zhuǎn)錄后水平上調(diào)節(jié)基因的表達(dá),在人的肝癌和正常肝臟中的表達(dá)存在顯著差異,miRNA的失調(diào)提示肝臟可能發(fā)生早期癌變[10]。已有研究結(jié)果顯示,m7G修飾相關(guān)基因METTL1和WRD4參與了對(duì)miRNA的調(diào)控,推測(cè)可能有更多的m7G修飾基因參與了miRNA的調(diào)控。本研究提取肝癌中具有預(yù)后價(jià)值的m7G相關(guān)基因,篩選下游相關(guān)miRNA構(gòu)建風(fēng)險(xiǎn)模型,評(píng)估模型的預(yù)測(cè)價(jià)值,為臨床治療提供參考。

本研究通過肝癌的預(yù)后數(shù)據(jù)篩選得到了12個(gè)m7G相關(guān)基因,并通過Cox和Lasso回歸分析構(gòu)建模型,篩選出11個(gè)受m7G修飾影響的miRNA。其中血清外泌體中的miR-4661-5p可以作為早期肝癌的潛在診斷標(biāo)志物[11],miR-301a-3p可以通過靶向VGLL4抑制肝癌的惡性進(jìn)展,miR-760可以調(diào)控下游的NACC-1促進(jìn)肝癌的發(fā)生發(fā)展[12],miR-561-5p通過影響CX3CR1/NK細(xì)胞的浸潤(rùn)和功能促進(jìn)轉(zhuǎn)移性肝癌的發(fā)生發(fā)展[13],海綿化的miR-346在circMDK促進(jìn)肝癌進(jìn)展的過程中發(fā)揮作用[14],而miR-9-3p則通過下調(diào)TAZ在肝癌細(xì)胞中發(fā)揮腫瘤抑制作用[15]。目前肝癌研究中尚未見miR-7156-5p、miR-3911、miR-513c-5p、miR-4652-3p、miR-548aq-5p報(bào)道,其功能和機(jī)制有待進(jìn)一步研究。

本文通過12個(gè)相關(guān)miRNA與生存時(shí)間構(gòu)建風(fēng)險(xiǎn)模型,計(jì)算風(fēng)險(xiǎn)評(píng)分并進(jìn)行PCA聚類分析,結(jié)果顯示風(fēng)險(xiǎn)評(píng)分能夠較好地將病人區(qū)分為高低風(fēng)險(xiǎn)組;Kaplan-Meier分析顯示,高風(fēng)險(xiǎn)和低風(fēng)險(xiǎn)組病人生存率存在差異,低風(fēng)險(xiǎn)組存活率顯著高于高風(fēng)險(xiǎn)組;ROC曲線分析顯示,該模型的1、3和5年生存率的AUC均大于0.6,表明預(yù)后模型具有預(yù)測(cè)價(jià)值,尤其是對(duì)肝癌病人3年生存率的預(yù)測(cè)。本文根據(jù)臨床因素建立預(yù)后模型列線圖分析顯示,風(fēng)險(xiǎn)評(píng)分對(duì)于病人預(yù)后的評(píng)價(jià)具有顯著意義,根據(jù)風(fēng)險(xiǎn)評(píng)分構(gòu)建的預(yù)后模型可以更好地預(yù)測(cè)臨床病人的預(yù)后情況,為臨床治療提供指導(dǎo)意見。

RNA甲基化修飾除了直接影響腫瘤發(fā)展外,還通過調(diào)節(jié)腫瘤免疫力來(lái)影響腫瘤治療的有效性[16]。

作為m7G的關(guān)鍵基因,METTL1在肝癌放療中可能通過m7G修飾參與了免疫微環(huán)境的調(diào)節(jié)[17]。本文利用已經(jīng)構(gòu)建的風(fēng)險(xiǎn)模型與免疫相關(guān)基因分析得到m7G修飾相關(guān)miRNA的下游基因,并分別進(jìn)行了功能富集分析和相關(guān)免疫功能和免疫細(xì)胞的分析,結(jié)果顯示,受體配體活性、信號(hào)受體激活劑活性與免疫浸潤(rùn)的抗腫瘤免疫力相關(guān)[18];ssGSEA分析結(jié)果顯示,在m7G相關(guān)的免疫環(huán)境中,輔助性T細(xì)胞、HLA和MHCⅠ的表達(dá)水平更高。不同風(fēng)險(xiǎn)組免疫浸潤(rùn)相關(guān)表達(dá)豐度分析顯示,高風(fēng)險(xiǎn)組肝癌病人可能更傾向于抗巨噬細(xì)胞或抗aDCs細(xì)胞的免疫治療,提示m7G相關(guān)的miRNA可能與免疫浸潤(rùn)有關(guān),這可為免疫治療提供可能的指導(dǎo)方案。

綜上所述,本文研究構(gòu)建了與m7G修飾相關(guān)miRNA的風(fēng)險(xiǎn)預(yù)后模型,將肝癌病人分為高低風(fēng)險(xiǎn)組,基于風(fēng)險(xiǎn)評(píng)分構(gòu)建預(yù)測(cè)模型對(duì)于肝癌病人的預(yù)后評(píng)價(jià)良好。風(fēng)險(xiǎn)評(píng)分可以作為預(yù)測(cè)病人預(yù)后的獨(dú)立預(yù)后因素,對(duì)于病人的預(yù)后預(yù)測(cè)有著顯著意義。免疫浸潤(rùn)相關(guān)分析提示高風(fēng)險(xiǎn)組病人的預(yù)后與巨噬細(xì)胞和aDCs細(xì)胞顯著相關(guān),可能作為免疫治療的潛在方向。該模型的建立有助于為肝癌病人的生存和免疫治療提供預(yù)測(cè)和指導(dǎo),但是其實(shí)際價(jià)值還有待于更大樣本量和相關(guān)實(shí)驗(yàn)的驗(yàn)證。