結腸腺癌能量代謝相關mRNA-LncRNA預后預測模型的構建與驗證*

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

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

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

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

1 材料與方法

1.1 數據收集

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

1.2 能量代謝相關mRNA與LncRNA的鑒定

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

1.3 功能富集分析

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

1.4 預后模型的構建

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

1.5 預后模型的驗證

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

1.6 統計學處理

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

2 結 果

2.1 能量代謝相關mRNA與LncRNA的鑒別

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

表1 r排名前10的堿基對

2.2 功能富集分析

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

2.3 預后預測模型的構建

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

續表2 多因素Cox回歸結果

2.4 預后模型的驗證

2.4.1能量代謝相關DEGs與DELs的生存分析

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

2.4.2風險得分與不同臨床病理特征的分析

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

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

2.4.3風險得分與患者臨床病理特征的生存分析

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

3 討 論

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

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

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

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