結(jié)腸腺癌能量代謝相關(guān)mRNA-LncRNA預(yù)后預(yù)測(cè)模型的構(gòu)建與驗(yàn)證*

黃浩東，劉 穎,劉小株,龔 軍,段敏捷,王惠來(lái),向天雨△

(1.重慶醫(yī)科大學(xué)醫(yī)學(xué)數(shù)據(jù)研究院,重慶 400016;2.重慶市急救醫(yī)療中心藥劑科,重慶 400014;3.重慶醫(yī)科大學(xué)附屬大學(xué)城醫(yī)院信息中心,重慶 401331)

結(jié)直腸癌(colorectal cancer,CRC)是常見(jiàn)的惡性腫瘤之一。其發(fā)病率居全球癌癥的第3位,病死率居全球癌癥的第2位[2]。在我國(guó),CRC的發(fā)生率和病死率也逐年上升[3]。結(jié)腸癌最常見(jiàn)的病理組織學(xué)分型為結(jié)腸腺癌(COAD)。信使RNA(messenger RNA,mRNA)是指導(dǎo)功能性蛋白質(zhì)合成的直接模板[5]。長(zhǎng)鏈非編碼RNA(long non-coding RNA,LncRNA) 是一種長(zhǎng)度大于200 bp且不具備明顯蛋白編碼能力的RNA轉(zhuǎn)錄物[6]。本研究旨在研究COAD中代謝相關(guān)RNA對(duì)預(yù)后的影響,尋找關(guān)鍵分子和調(diào)控通路,以更好地指導(dǎo)COAD的預(yù)防、診斷和治療。

本研究通過(guò)癌癥基因組圖譜(TCGA)獲取了COAD患者的轉(zhuǎn)錄本數(shù)據(jù)和臨床病理數(shù)據(jù)。采用相關(guān)性分析與Cox回歸相結(jié)合的方法構(gòu)建COAD患者能量代謝相關(guān)mRNA-LncRNA預(yù)后預(yù)測(cè)模型,并通過(guò)受試者工作特征(AUC)和Kaplan-Meier曲線等分析方法,在高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組中進(jìn)行預(yù)后預(yù)測(cè)模型的評(píng)估。以便為基于mRNA和LncRNA的COAD分子機(jī)制研究提供新視角。

1 材料與方法

1.1 數(shù)據(jù)收集

從TCGA(https://portal.gdc.cancer.gov；level3)下載TCGA-COAD項(xiàng)目COAD組織與正常組織的RNA-seq數(shù)據(jù)集和相關(guān)臨床資料。 從GENCODE (https://www.gencodegenes.org/human/)數(shù)據(jù)庫(kù)下載基因的注釋信息(包括19 955個(gè)mRNA和16 888個(gè)LncRNA注釋信息)。將下載的轉(zhuǎn)錄本數(shù)據(jù)注釋分為mRNA與LncRNA兩個(gè)部分。

1.2 能量代謝相關(guān)mRNA與LncRNA的鑒定

利用edgeR包對(duì)差異表達(dá)基因進(jìn)行篩選COAD和正常對(duì)照組間的差異表達(dá)的mRNA(DEGs)與差異表達(dá)的(DELs),篩選條件為錯(cuò)誤發(fā)現(xiàn)率(FDR)<0.05,且基因表達(dá)差異倍數(shù)|logFC|>1,結(jié)果以火山圖的形式表示。從MSigDB (version 7.0;http://software.broadinstitute.org/gsea/msigdb/)數(shù)據(jù)庫(kù)收集能量代謝相關(guān)基因組,將DEGs與能量代謝相關(guān)基因組相交得到能量代謝相關(guān)DEGs。對(duì)DELs與能量代謝相關(guān)DEGs 的進(jìn)行共表達(dá)分析,選擇相關(guān)系數(shù)(r)>0.6且P<0.05的共表達(dá)對(duì)來(lái)得到能量代謝相關(guān)的LncRNA。

1.3 功能富集分析

使用 DAVID(https://david.ncifcrf.gov/)數(shù)據(jù)庫(kù)[9]將共表達(dá)分析得到的能量代謝相關(guān)LncRNA調(diào)控的mRNA進(jìn)行GO功能富集,以明確這些基因所涉及的生物學(xué)過(guò)程、分子功能及細(xì)胞組成。同時(shí)進(jìn)行KEGG通路分析,以明確涉及的特定能量代謝過(guò)程，結(jié)果以氣泡圖形式進(jìn)行展示。

1.4 預(yù)后模型的構(gòu)建

將共表達(dá)分析得到的能量代謝相關(guān)mRNA與LncRNA表達(dá)量與患者臨床信息合并,并進(jìn)行單因素Cox回歸和多因素Cox回歸,以P<0.05為差異具有統(tǒng)計(jì)學(xué)意義。最后風(fēng)險(xiǎn)得分(risk score)= (基因1的風(fēng)險(xiǎn)系數(shù)×基因1的表達(dá)量)+(基因2的風(fēng)險(xiǎn)系數(shù)×基因2的表達(dá)量)+(基因3的風(fēng)險(xiǎn)系數(shù)×基因3的表達(dá)量)+…+(基因n的風(fēng)險(xiǎn)系數(shù)×基因n的表達(dá)量)。

1.5 預(yù)后模型的驗(yàn)證

以中位風(fēng)險(xiǎn)評(píng)分作為分界點(diǎn),將COAD患者分為高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組。利用Kaplan-Meier 曲線和受試者工作特征 (ROC) 曲線評(píng)估能量代謝相關(guān)特征的預(yù)測(cè)能力,并與其他臨床特征的預(yù)測(cè)能力進(jìn)行對(duì)比。為了驗(yàn)證風(fēng)險(xiǎn)評(píng)分是否能獨(dú)立于其他臨床病理特征,分別采用：(1)非參數(shù)檢驗(yàn)檢測(cè)患者不同臨床病理特征分組中風(fēng)險(xiǎn)評(píng)分是否有差異;(2)單因素Cox回歸和多因素Cox回歸驗(yàn)證風(fēng)險(xiǎn)評(píng)分是否為患者不良預(yù)后的獨(dú)立危險(xiǎn)因素。

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

采用SPSS25.0和R語(yǔ)言(4.3.0版)進(jìn)行統(tǒng)計(jì)分析和畫(huà)圖,差異表達(dá)分析采用edgeR包實(shí)現(xiàn)；代謝相關(guān)基因與LncRNA 采用Pearson相關(guān)分析進(jìn)行共表達(dá)分析,通過(guò)cor.test函數(shù)實(shí)現(xiàn)；cox回歸分析采用survival包實(shí)現(xiàn)；通過(guò)Enhanced Volcano包繪制火山圖；通過(guò)survminer包繪制Kaplan-Meier 曲線,log-rank法檢驗(yàn)P值；通過(guò)survival ROC包繪制ROC曲線；通過(guò)ggplot2包繪制氣泡圖；通過(guò)forestplot包繪制Cox回歸森林圖；非參數(shù)檢驗(yàn)采用Mann-WhitneyU檢驗(yàn)。以P<0.05為差異具有統(tǒng)計(jì)學(xué)意義。

2 結(jié) 果

2.1 能量代謝相關(guān)mRNA與LncRNA的鑒別

共收集到TCGA-COAD項(xiàng)目中的471個(gè)COAD組織與41個(gè)正常組織,其中有完整臨床病理特征信息的患者387例,基于GENCODE數(shù)據(jù)庫(kù)共標(biāo)記19 551個(gè)mRNA 與14 036個(gè)LncRNA。edge包篩選出4 851個(gè)DEGs與1 468個(gè)DELs(圖1A)。通過(guò)搜索MsigDB數(shù)據(jù)庫(kù),共下載了1 384個(gè)能量代謝相關(guān)基因,其中碳水化合物代謝294個(gè),脂質(zhì)代謝742個(gè),氨基酸及衍生物代謝375個(gè)。將上述能量代謝相關(guān)基因與DEGs相交,結(jié)果顯示326個(gè)DEGs與能量代謝相關(guān)(圖1B)。將326個(gè)能量代謝相關(guān)DEGs與1 468個(gè)DELs進(jìn)行了326×1 468個(gè)Pearson相關(guān)分析,根據(jù)r>0.6且P<0.05,得到共表達(dá)堿基對(duì)2 079對(duì),其中125個(gè)mRNA,上調(diào)46個(gè),下調(diào)79個(gè)；451個(gè)LncRNA,上調(diào)261個(gè),下調(diào)190個(gè)。r排名前10的堿基對(duì)見(jiàn)表1。

表1 r排名前10的堿基對(duì)

2.2 功能富集分析

將共表達(dá)分析得到mRNA進(jìn)行KEGG與GO通路分析。KEGG通路結(jié)果顯示共表達(dá)分析得到的LncRNA可能通過(guò)各種代謝途徑參與了COAD的形成過(guò)程(圖2A)。GO生物學(xué)過(guò)程富集結(jié)果顯示與氧化還原過(guò)程、脂質(zhì)代謝過(guò)程和病毒轉(zhuǎn)錄等相關(guān)(圖2B)。GO細(xì)胞組成富集結(jié)果顯示與細(xì)胞質(zhì)、細(xì)胞外泌體和細(xì)胞膜等相關(guān)(圖2C)。GO分子功能富集結(jié)果顯示與蛋白質(zhì)同源二聚活性、鐵離子結(jié)合和磷酸吡哆醛結(jié)合等相關(guān)(圖2D)。

2.3 預(yù)后預(yù)測(cè)模型的構(gòu)建

將上述的125個(gè)能量代謝相關(guān)DEGs和451個(gè)能量代謝相關(guān)DELs與臨床預(yù)后數(shù)據(jù)相結(jié)合,單因素Cox回歸分析得到 15個(gè)mRNA與22個(gè)LncRNA與預(yù)后相關(guān)；將其納入多因素Cox回歸分析,結(jié)果顯示LRP2、CTC-428G20.6、LINC02257、PRR7-AS1、RP11-29G8.3、RP11-677M14.3這6個(gè)堿基與COAD患者預(yù)后相關(guān)(表2),根據(jù)此結(jié)果,得到風(fēng)險(xiǎn)得分=LRP2×4.037 98+CTC-428G20.6×(-1.937 56)+LINC02257×1.046 69+PRR7-AS1×0.869 12+RP11-29G8.3×(-1.459 44)+RP11-677M14.3×0.841 68。

表2 多因素Cox回歸結(jié)果

續(xù)表2 多因素Cox回歸結(jié)果

2.4 預(yù)后模型的驗(yàn)證

2.4.1能量代謝相關(guān)DEGs與DELs的生存分析

將上述得到的風(fēng)險(xiǎn)得分以中位數(shù)分為高風(fēng)險(xiǎn)組(風(fēng)險(xiǎn)得分>中位數(shù))與低風(fēng)險(xiǎn)組(風(fēng)險(xiǎn)得分 <中位數(shù)),Kaplan-Meier分析得到的生存曲線見(jiàn)圖3A,高、低風(fēng)險(xiǎn)組的中位數(shù)生存期分別為5.066年和8.334年。圖3B、圖3C和圖3D分別為1年、3年和5年生存時(shí)間對(duì)應(yīng)的ROC曲線?？梢钥闯鲲L(fēng)險(xiǎn)得分在預(yù)測(cè)1年生存期時(shí)對(duì)應(yīng)的AUC=0.767,與其他的臨床病理特征相比最高；在預(yù)測(cè)3年生存期時(shí),對(duì)應(yīng)的AUC=0.694,低于病理分期Stage AUC=0.745與M分期=0.704；在預(yù)測(cè)5年生存期時(shí),對(duì)應(yīng)的AUC=0.691,低于病理分期Stage AUC=0.706。綜合來(lái)說(shuō),風(fēng)險(xiǎn)得分在預(yù)測(cè)短期生存期時(shí)有一定的優(yōu)勢(shì),在預(yù)測(cè)長(zhǎng)期生存期時(shí),與病理分期預(yù)測(cè)性能相似。

2.4.2風(fēng)險(xiǎn)得分與不同臨床病理特征的分析

為比較風(fēng)險(xiǎn)得分是否與COAD患者的臨床病理特征具有相關(guān)性,以風(fēng)險(xiǎn)得分公式得出的值,以不同臨床病理特征為分組,進(jìn)行Mann-WhitneyU檢驗(yàn),分析不同組間的風(fēng)險(xiǎn)得分值的差異,結(jié)果如表3所示。從表中得出風(fēng)險(xiǎn)得分與COAD患者的不同臨床病理特征相關(guān)性不大。

表3 不同臨床病理特征與風(fēng)險(xiǎn)得分的Mann-Whitney U檢驗(yàn)

2.4.3風(fēng)險(xiǎn)得分與患者臨床病理特征的生存分析

如圖4所示,風(fēng)險(xiǎn)得分與患者的其他臨床病理特征進(jìn)行單因素Cox回歸與多因素Cox回歸檢驗(yàn),得出高風(fēng)險(xiǎn)得分是患者不良預(yù)后的獨(dú)立危險(xiǎn)因素。

3 討 論

近年來(lái),隨著生物信息學(xué)的發(fā)展,越來(lái)越多的研究利用mRNA或LncRNA的表達(dá)量預(yù)測(cè)肝癌、乳腺癌、胰腺癌和結(jié)直腸癌等患者的預(yù)后[10-13]。能量代謝相關(guān)mRNA和LncRNA的結(jié)直腸腺癌預(yù)測(cè)模型尚未構(gòu)建。能量代謝過(guò)程參與生命發(fā)生、發(fā)展的全過(guò)程,在COAD的進(jìn)展中起重要作用,糖、脂質(zhì)和氨基酸代謝過(guò)程產(chǎn)生三磷酸腺苷(adenosine triphosphate,ATP),而腫瘤細(xì)胞的惡性特征(快速增殖、侵襲和遷移)需依靠大量ATP維持[4,14]。并且一些研究表明,LncRNA可以通過(guò)調(diào)節(jié)能量代謝相關(guān)基因影響癌癥進(jìn)展。如：WANG等[7]研究發(fā)現(xiàn),LINRIS在總生存率較差的患者CRC組織中上調(diào),敲除LncRNA LINRIS減弱了CRC細(xì)胞中Myc介導(dǎo)的糖酵解途徑,從而抑制CRC細(xì)胞的生長(zhǎng)。TANG等[8]研究發(fā)現(xiàn),LncRNA GLCC1通過(guò)穩(wěn)定c-Myc的泛素化,從而重新編程葡萄糖代謝促進(jìn)CRC增殖。因此尋找出能量代謝相關(guān)生物標(biāo)志物并建立預(yù)后預(yù)測(cè)模型有利于COAD患者的個(gè)性化治療。

本研究通過(guò)對(duì)比COAD組織與正常組織,篩選出DEGs與DELs,然后通過(guò)DEGs與MsigDB數(shù)據(jù)庫(kù)提供的糖、脂質(zhì)、氨基酸三大類代謝基因組得到能量代謝相關(guān)DEGs,再將其與DELs進(jìn)行共表達(dá)分析進(jìn)一步篩選得到能量代謝相關(guān)DELs與DEGs。對(duì)這些能量代謝相關(guān)DEGs進(jìn)行富集分析發(fā)現(xiàn)它們除了在hsa01100:Metabolic pathways通路聚集較多外,Count數(shù)較高的通路還有hsa00564:Glycerophospholipid metabolism、hsa03320:PPAR signaling pathway和hsa04975:Fat digestion and absorption等。而生物學(xué)過(guò)程富集于氧化還原過(guò)程與脂質(zhì)代謝過(guò)程。然后,本研究進(jìn)一步采用單因素Cox回歸與多因素Cox回歸得到具有預(yù)后預(yù)測(cè)意義的RNA,其中包含1個(gè)mRNA和5個(gè)LncRNA,分別為L(zhǎng)RP2、CTC-428G20.6、LINC02257 、PRR7-AS1、RP11-29G8.3和RP11-677M14.3。其中LRP2是唯一1個(gè)mRNA,LRP2基因編碼的Megalin是一種配體結(jié)合的跨膜蛋白,通過(guò)和不同的配體結(jié)合發(fā)揮神經(jīng)及內(nèi)分泌調(diào)節(jié)、抗凋亡等作用[15-17]。Megalin/LRP2已在多篇文章中報(bào)道,如JAKOVAC等[18]推測(cè)其可能有助于鱗狀上皮中異常細(xì)胞的更好存活和腫瘤發(fā)生；ANDERSEN等[19]表明,黑色素瘤細(xì)胞中的Megalin對(duì)細(xì)胞的生存較重要,因?yàn)楹谏亓黾?xì)胞中的Megalin/LRP2表達(dá)的降低會(huì)導(dǎo)致其增殖和存活率降低。并且LRP2已被FEDIRKO等[20]證實(shí)參與CRC中維生素D的代謝過(guò)程,從而增加西歐人群CRC的患病風(fēng)險(xiǎn)。XIAO等[21]通過(guò)生物信息學(xué)方式發(fā)現(xiàn)LINC02257高表達(dá)會(huì)導(dǎo)致COAD預(yù)后不良,與本研究結(jié)果相符。其他4個(gè)LncRNA尚未見(jiàn)文獻(xiàn)報(bào)道。最后,本研究利用上述篩選出的6個(gè)能量代謝相關(guān)的RNA構(gòu)建了COAD預(yù)后預(yù)測(cè)模型。

為了驗(yàn)證這6個(gè)能量代謝相關(guān)RNA是否具有預(yù)后意義,本研究根據(jù)構(gòu)建的預(yù)后模型建立了一個(gè)風(fēng)險(xiǎn)得分公式,按照其中位數(shù)分為了低風(fēng)險(xiǎn)組與高風(fēng)險(xiǎn)組。在Kaplan-Meier生存分析中,低風(fēng)險(xiǎn)組與高風(fēng)險(xiǎn)組的中位數(shù)生存期分別是5.066年和8.334年(log-rank檢驗(yàn)顯示P<0.05),低風(fēng)險(xiǎn)組COAD患者生存期明顯高于高風(fēng)險(xiǎn)組生存期。而風(fēng)險(xiǎn)得分在1、3和5年的ROC曲線下面積分別是0.767、0.694和0.691。其中在1年的ROC曲線下面積相比于COAD患者其他臨床病理特征的ROC曲線下面積最高,而在3年與5年的ROC曲線下面積略低于COAD患者病理分期,因此,風(fēng)險(xiǎn)得分在預(yù)測(cè)短期生存期時(shí)有一定的優(yōu)勢(shì),在預(yù)測(cè)長(zhǎng)期生存期時(shí),與病理分期預(yù)測(cè)性能相似。此外本研究對(duì)風(fēng)險(xiǎn)得分與不同的臨床病理特征進(jìn)行非參數(shù)檢驗(yàn),結(jié)果顯示,不同的臨床病理特征與風(fēng)險(xiǎn)得分不相關(guān)。為了探討高風(fēng)險(xiǎn)得分是否為COAD患者不良預(yù)后的獨(dú)立危險(xiǎn)因素,本研究納入臨床病理特征與風(fēng)險(xiǎn)得分進(jìn)行Cox回歸分析,結(jié)果顯示高風(fēng)險(xiǎn)組生存時(shí)間明顯低于低風(fēng)險(xiǎn)組(HR=3.78,95%CI：2.29～6.22,P<0.05),高風(fēng)險(xiǎn)得分為COAD患者不良預(yù)后的獨(dú)立危險(xiǎn)因素。本研究也有一些不足之處,例如研究篩選出的6個(gè)RNA,其中4個(gè)LncRNA(CTC-428G20.6 、PRR7-AS1、RP11-29G8.3和RP11-677M14.3)還沒(méi)有相關(guān)報(bào)道,對(duì)其如何影響COAD的發(fā)展尚不清楚,因此需要進(jìn)一步的前瞻性實(shí)驗(yàn)研究。

本研究利用6個(gè)RNA(LRP2、CTC-428G20.6、LINC02257 、PRR7-AS1、RP11-29G8.3和RP11-677M14.3 )建立的COAD患者的能量代謝相關(guān)預(yù)后預(yù)測(cè)模型具有較好的性能。對(duì)COAD的個(gè)性化治療具有一定的積極意義,但仍然需要進(jìn)一步的研究。