嘌呤代謝相關(guān)基因預(yù)測肺腺癌預(yù)后模型的構(gòu)建及驗(yàn)證

李杰 李強(qiáng)

肺腺癌是一種常見的惡性腫瘤，發(fā)病率和死亡率都很高[1]，盡管治療肺腺癌技術(shù)已取得了巨大進(jìn)展，但肺腺癌患者的5 年總生存期（overall survival，OS）仍低于15%[2]。因此，探索影響肺腺癌患者預(yù)后的生物標(biāo)志物非常有必要。嘌呤不僅能提供必要的能量和輔助因子以促進(jìn)細(xì)胞存活和增殖，而且還參與免疫反應(yīng)和宿主-腫瘤相互作用[3]。研究發(fā)現(xiàn)，嘌呤代謝與腫瘤的發(fā)生、發(fā)展密不可分[4]，其不僅在前列腺癌細(xì)胞的增殖和侵襲中發(fā)揮著重要的作用[5]，而且還影響著膀胱癌和乳腺癌的發(fā)生、發(fā)展[6]。但目前有關(guān)嘌呤代謝與肺腺癌關(guān)系的報道很少，因此，本研究基于癌癥基因組圖譜（the cancer genome atlas，TCGA）數(shù)據(jù)庫中肺腺癌患者的數(shù)據(jù)，通過生物信息學(xué)技術(shù)建立預(yù)后模型，探討嘌呤代謝相關(guān)基因與肺腺癌的預(yù)后關(guān)系。

1 材料和方法

1.1 數(shù)據(jù)收集 從TCGA 數(shù)據(jù)庫官網(wǎng)（https://www.cancer.gov/tcga）中獲取肺腺癌臨床數(shù)據(jù)和mRNA 轉(zhuǎn)錄組數(shù)據(jù)，其中肺腺癌mRNA 轉(zhuǎn)錄組樣本數(shù)據(jù)包含535個肺腺癌組織和59 個癌旁組織，從基因表達(dá)數(shù)據(jù)庫（gene expression database，GEO）官網(wǎng)（http://www.ncbi.nlm.nih.gov/geo）中獲取GSE26939 的臨床數(shù)據(jù)和mRNA轉(zhuǎn)錄組數(shù)據(jù)，排除沒有生存狀態(tài)或隨訪時間＜1 d 的樣本。此外，從人類基因（Genecards）數(shù)據(jù)庫官網(wǎng)（https://www.genecards.org/）獲取嘌呤代謝相關(guān)基因3 968 個。

1.2 差異基因篩選 使用R 4.1.5 軟件中的“l(fā)imma”包，以|log2差異倍數(shù)（fold change,FC）|＞2 且P＜0.05 為標(biāo)準(zhǔn)，篩選表達(dá)差異的嘌呤代謝相關(guān)基因，再使用“ggplot2”和“pheatmap”包，繪制成火山圖和熱圖進(jìn)行可視化。

1.3 預(yù)后差異基因富集分析 將上述差異基因與肺腺癌臨床數(shù)據(jù)合并，篩選出與肺腺癌預(yù)后有關(guān)的嘌呤代謝相關(guān)基因，然后使用R 軟件中的“clusterProfiler”和“enrichplot”包對預(yù)后相關(guān)的嘌呤代謝相關(guān)基因進(jìn)行京都基因與基因組百科全書（Kyoto encyclopedia of genes and genomes，KEGG）和基因本體論（gene ontology，GO）富集分析，P＜0.05 為差異有統(tǒng)計學(xué)意義[7]。

1.4 預(yù)后模型的建立及分析 首先，使用Perl 軟件將上述差異基因與肺腺癌臨床數(shù)據(jù)進(jìn)行整合，剔除無生存狀態(tài)記錄的患者；其次，使用R 軟件中的“survival”和“forestplot”包進(jìn)行單因素及多因素獨(dú)立預(yù)后Cox 回歸分析，篩選出P＜0.05 的嘌呤代謝差異基因；然后使用“survival”和“glmnet”包對上述篩出的差異基因進(jìn)行套索（the least absolute shrinkage and selection operator，LASSO）回歸分析，建立風(fēng)險評分預(yù)后模型。風(fēng)險評分公式為：風(fēng)險評分與預(yù)后有關(guān)的嘌呤代謝相關(guān)基因的表達(dá)量（i）（Coef 為基因在多因素Cox 回歸分析中的回歸系數(shù)，n為與預(yù)后有關(guān)的嘌呤代謝相關(guān)基因的總數(shù)目），計算每例患者的風(fēng)險評分，根據(jù)風(fēng)險評分中位值，將患者分為高風(fēng)險組和低風(fēng)險組，再使用“survivalROC”包繪制ROC 曲線，并計算約登指數(shù)，從而預(yù)測模型的準(zhǔn)確性[6]。

1.5 預(yù)后模型的評價 使用R 軟件中的“survival”包繪制高、低風(fēng)險組評分分布曲線、生存狀態(tài)圖及建?；虮磉_(dá)熱圖。使用“survminer”包繪制Kaplan-Meier生存曲線，評估高、低風(fēng)險組患者的OS。使用“timeR‐OC”包繪制時間依賴性ROC 曲線，計算1、3、5 年肺腺癌患者OS 的AUC，以評價模型的準(zhǔn)確性，AUC 越高代表模型越準(zhǔn)確。

1.6 獨(dú)立預(yù)后因素分析 將性別、年齡、腫瘤分期、腫瘤類型、吸煙和風(fēng)險評分作為主要變量進(jìn)行單因素及多因素Cox 回歸分析，使用R 軟件中的“forestplot”包繪制森林圖來分析影響肺腺癌患者的預(yù)后因素。使用“rms”包建立列線圖和決策曲線圖來估計1、3、5 年肺腺癌患者的總復(fù)發(fā)率。

1.7 預(yù)后模型驗(yàn)證 從GEO 數(shù)據(jù)庫中選擇GSE26939為外部檢驗(yàn)集，根據(jù)風(fēng)險評分公式，計算驗(yàn)證集中每個樣本的風(fēng)險評分。根據(jù)風(fēng)險評分中位值，將樣本分為高風(fēng)險組和低風(fēng)險組，使用Kaplan-Meier 生存曲線、ROC 曲線的AUC、生存狀態(tài)圖及建?；虮磉_(dá)熱圖來檢驗(yàn)預(yù)后模型的準(zhǔn)確性。

2 結(jié)果

2.1 差異基因分析 從3 968 個嘌呤代謝相關(guān)基因中，共篩選出355 個有表達(dá)差異的嘌呤代謝相關(guān)基因，其中155 個上調(diào)，200 個下調(diào)，見圖1（插頁）。

圖1 嘌呤代謝相關(guān)基因差異分析（A：差異基因的熱圖；B：差異基因的火山圖）

2.2 功能富集分析 KEGG 結(jié)果顯示，嘌呤代謝相關(guān)基因主要富集在細(xì)胞周期、黃體酮介導(dǎo)的卵母細(xì)胞成熟、卵母細(xì)胞減數(shù)分裂、細(xì)胞衰老、p53 信號通路、神經(jīng)活性配體-受體相互作用、人T 細(xì)胞白血病病毒1 感染、膽汁分泌、唾液分泌和造血細(xì)胞譜系等通路上（均P＜0.05），見圖2A（插頁）。分子功能富集分析顯示，嘌呤代謝相關(guān)基因與蛋白質(zhì)絲氨酸/蘇氨酸激酶活性、作用于DNA 的催化酶活性、組蛋白激酶活性、蛋白酪氨酸激酶活性、蛋白質(zhì)絲氨酸/蘇氨酸/酪氨酸激酶活性、ATP 酶活性、單加氧酶活性、G 蛋白偶聯(lián)肽受體活性、肽受體活性、DNA 依賴性ATP 酶活性等均有關(guān)（均P＜0.05），見圖2B（插頁）。生物學(xué)過程富集分析顯示，嘌呤代謝相關(guān)基因與核分裂、細(xì)胞器分裂、有絲分裂核分裂、有絲分裂核分裂的調(diào)節(jié)、核分裂調(diào)節(jié)、有絲分裂中期/后期轉(zhuǎn)變的調(diào)節(jié)、染色體分離的調(diào)節(jié)、細(xì)胞周期中期/后期轉(zhuǎn)變的調(diào)節(jié)、有絲分裂細(xì)胞周期的中期/末期過渡、細(xì)胞周期的中期/后期過渡等均有關(guān)（均P＜0.05），見圖2C（插頁）。細(xì)胞定位富集分析顯示，嘌呤代謝相關(guān)基因與染色體區(qū)域、染色體著絲粒區(qū)、濃縮染色體、濃縮染色體著絲粒區(qū)、濃縮染色體著絲粒、濃縮染色體外動著絲粒、凝聚核染色體、凝聚核染色體著絲粒區(qū)、著絲粒、紡錘體等均有關(guān)（均P＜0.05），見圖2D（插頁）。

圖2 功能富集分析（A：通路富集分析；B：分子功能富集分析；C：生物學(xué)過程富集分析；D：細(xì)胞定位富集分析）

2.3 肺腺癌預(yù)后模型預(yù)測患者預(yù)后 對上述355 個有表達(dá)差異的嘌呤代謝相關(guān)基因進(jìn)行單因素Cox 回歸分析，得到24 個與肺腺癌患者預(yù)后相關(guān)的基因，見表1。使用LASSO 回歸對這24 個基因進(jìn)行交叉驗(yàn)證分析，得到7 個與預(yù)后密切相關(guān)的基因，見圖3（插頁）。將這7個基因進(jìn)行多因素Cox 回歸分析，最終得到5 個與肺腺癌患者預(yù)后顯著相關(guān)的基因（CD19、CYP17A1、KH‐DRBS2、INHA、PLK1），其中INHA 和PLK1 基因HR值＞1，為高風(fēng)險基因，表示INHA 和PLK1 高表達(dá)對肺腺癌的預(yù)后較差；CD19、CYP17A1 和KHDRBS2 基因HR值＜1，為低風(fēng)險基因，表示CD19、CYP17A1 和KH‐DRBS2 低表達(dá)對肺腺癌的預(yù)后較差，見表2。根據(jù)5個嘌呤代謝相關(guān)基因的β值和基因的表達(dá)量計算每例患者的風(fēng)險評分，風(fēng)險評分=（-0.542×CD19 表達(dá)量）+（-0.542×CYP17A1 表達(dá)量）+（-0.463×KH‐DRBS2 表達(dá)量）+（0.134×INHA表達(dá)量）+（0.129×PLK1表達(dá)量）。

表1 單因素Cox 回歸分析結(jié)果

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

圖3 預(yù)后風(fēng)險評分模型評估（A、B：套索回歸分析）

2.4 預(yù)后模型的預(yù)測價值 5 個基因的生存分析見圖4（插頁）。根據(jù)風(fēng)險評分的中位值（1.77），將患者分為低風(fēng)險組和高風(fēng)險組。Kaplan-Meier 生存曲線（紅色為低風(fēng)險組，藍(lán)色為高風(fēng)險組）顯示，低風(fēng)險組患者OS高于高風(fēng)險組（HR=3.85，95%CI：2.79～5.31，P＜0.01），見圖4A（封三）。生存狀態(tài)圖及建?；虮磉_(dá)熱圖（紅色為高風(fēng)險組，藍(lán)色為低風(fēng)險組）顯示，高風(fēng)險組患者預(yù)后較差，見圖4B（封三）。ROC 曲線分析顯示，1、3、5年肺腺癌患者OS 的AUC 分別為0.76、0.74、0.77，見圖4C（封三）。綜上表明，此預(yù)后模型對肺腺癌患者的預(yù)后預(yù)測準(zhǔn)確性很高。

圖4 預(yù)后模型的預(yù)測價值（A：Kaplan-Meier 生存風(fēng)險曲線；B：生存狀態(tài)圖和建?；虮磉_(dá)熱圖；C：時間依賴性ROC 曲線）

2.5 獨(dú)立預(yù)后因素分析 單因素Cox 回歸分析顯示腫瘤分期和風(fēng)險評分均是肺腺癌患者預(yù)后的獨(dú)立危險因素（均P＜0.05），見圖5A。 將上述兩個危險因素納入多因素Cox 回歸分析，結(jié)果顯示腫瘤分期和風(fēng)險評分均是肺腺癌患者預(yù)后的獨(dú)立危險因素（均P＜0.05），見圖5B。列線圖得分可推測患者未來1、3、5 年的生存率，見圖5C。決策曲線顯示其預(yù)測效能較好，見圖5D。

圖5 獨(dú)立預(yù)后生存分析（A：單因素Cox 回歸分析；B：多因素Cox 回歸分析；C：列線圖；D：決策曲線）

2.6 外部數(shù)據(jù)集驗(yàn)證 從GEO數(shù)據(jù)庫中選擇GSE26939為外部驗(yàn)證集，GSE26939 數(shù)據(jù)集的風(fēng)險評分中位值為0.99，通過驗(yàn)證集的Kaplan-Meier 生存曲線（紅色為低風(fēng)險組，藍(lán)色為高風(fēng)險組）顯示，低風(fēng)險組患者OS 高于高風(fēng)險組（HR=8.94，95%CI：2.96～27.01，P＜0.01），見圖6A（封三）。生存狀態(tài)圖及建?；虮磉_(dá)熱圖（紅色為高風(fēng)險組，藍(lán)色為低風(fēng)險組）顯示，高風(fēng)險組患者的預(yù)后較差，見圖6B（封三）。驗(yàn)證集的ROC 曲線1、3、5 年肺腺癌患者OS 的AUC 分別為0.96、0.82、0.84，見圖6C（封三）。綜上表明該模型在外部驗(yàn)證集中有很好的預(yù)測性能。

圖6 預(yù)后模型的性能評估（A：Kaplan-Meier 生存分析曲線；B：生存狀態(tài)圖和建?；虮磉_(dá)熱圖；C：時間依賴性ROC 曲線）

3 討論

本研究以嘌呤代謝相關(guān)基因作為背景，最終建立了由5 個嘌呤代謝相關(guān)基因（CD19、CYP17A1、KH‐DRBS2、INHA、PLK1）組成的預(yù)后模型。結(jié)果顯示INHA 和PLK1 高表達(dá)對肺腺癌的預(yù)后較差，CD19、CYP17A1 和KHDRBS2 低表達(dá)對肺腺癌的預(yù)后較差。GO 和KEGG 富集分析表明這些基因主要富集在細(xì)胞周期和細(xì)胞分裂信號通路上，由此可大膽推測，嘌呤代謝相關(guān)基因可能通過促進(jìn)細(xì)胞增殖來影響肺腺癌患者的預(yù)后。

INHA 是重要的代謝相關(guān)基因之一，據(jù)報道，INHA的功能不僅與前列腺癌雄激素非依賴性轉(zhuǎn)移和卵巢腫瘤的血管生成有關(guān)[8]，還與肺腺癌的免疫浸潤有關(guān)[9]。PLK1 參與了多種細(xì)胞周期調(diào)節(jié)途徑，其可作為G2/M 檢查點(diǎn)并負(fù)責(zé)中心體、紡錘體組裝和染色體分離的調(diào)節(jié)[10]。既往研究發(fā)現(xiàn)，PLK1 的過表達(dá)與多種癌癥的發(fā)生、發(fā)展密切相關(guān)[10-11]。較高的PLK1 轉(zhuǎn)錄和蛋白水平對胃癌患者的預(yù)后有不良影響[12]。PLK1 高表達(dá)也與卵巢透明細(xì)胞癌密切相關(guān)[13]。還有研究證實(shí)，PLK1 的表達(dá)可以預(yù)測轉(zhuǎn)移性非小細(xì)胞肺癌患者的生存[14]，此外，研究發(fā)現(xiàn)PLK1 還與肺腺癌的免疫微環(huán)境有關(guān)[15]。CD19 是重要的B 細(xì)胞表面標(biāo)志物，目前被廣泛應(yīng)用于血液系統(tǒng)的腫瘤治療中[16]。據(jù)報道CD19 可能在肺腺癌的不同階段發(fā)揮雙重作用，并且還與免疫預(yù)后相關(guān)[17]。CYP17A1 可以將睪酮轉(zhuǎn)化為雌二醇，是非小細(xì)胞肺癌易感基因之一，但其多態(tài)性與亞洲人群中的非小細(xì)胞肺癌發(fā)展無關(guān)[18]。有研究表明CYP17A1較高的基因突變和拷貝數(shù)變異可能通過影響B(tài) 細(xì)胞功能來影響肺腺癌患者的易感性[19]。KHDRBS2 可以作為抑癌基因，既往研究報道，KHDRBS2 表達(dá)水平可以預(yù)測肺癌患者OS[20]，這與本研究結(jié)果一致。

綜上所述，嘌呤代謝相關(guān)基因可影響肺腺癌患者的預(yù)后，有潛在的臨床應(yīng)用價值。但本研究仍有不足之處，如本研究的全部數(shù)據(jù)均來源于公共數(shù)據(jù)庫，未進(jìn)行實(shí)驗(yàn)研究及臨床數(shù)據(jù)驗(yàn)證，因此，未來的研究需要多中心、大樣本來進(jìn)一步驗(yàn)證嘌呤代謝相關(guān)基因預(yù)測肺腺癌患者預(yù)后的準(zhǔn)確性，為今后的治療提供潛在的治療靶點(diǎn)。