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

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

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

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

1 材料與方法

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

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

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

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

1.3 功能富集分析

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

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

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

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

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

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

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

2 結(jié) 果

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

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

表1 r排名前10的堿基對

2.2 功能富集分析

將共表達(dá)分析得到mRNA進(jìn)行KEGG與GO通路分析。KEGG通路結(jié)果顯示共表達(dá)分析得到的LncRNA可能通過各種代謝途徑參與了COAD的形成過程(圖2A)。GO生物學(xué)過程富集結(jié)果顯示與氧化還原過程、脂質(zhì)代謝過程和病毒轉(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ù)測模型的構(gòu)建

將上述的125個能量代謝相關(guān)DEGs和451個能量代謝相關(guān)DELs與臨床預(yù)后數(shù)據(jù)相結(jié)合,單因素Cox回歸分析得到 15個mRNA與22個LncRNA與預(yù)后相關(guān)；將其納入多因素Cox回歸分析,結(jié)果顯示LRP2、CTC-428G20.6、LINC02257、PRR7-AS1、RP11-29G8.3、RP11-677M14.3這6個堿基與COAD患者預(yù)后相關(guān)(表2),根據(jù)此結(jié)果,得到風(fēng)險得分=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ù)后模型的驗證

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

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

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

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

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

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

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

3 討 論

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

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

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

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