基于TCGA數(shù)據(jù)庫(kù)建立膽管癌自噬相關(guān)基因預(yù)后預(yù)測(cè)模型及其應(yīng)用

史華帝，左瑜芳，鐘富蘭，易小瓊，徐祖敏

廣東醫(yī)科大學(xué)附屬醫(yī)院，廣東湛江524000

多數(shù)膽管癌患者確診時(shí)已至晚期，失去手術(shù)機(jī)會(huì)，并且放化療效果不佳。既往25 年的數(shù)據(jù)表明，膽管癌的發(fā)病率和病死率在增加，然而膽管癌的預(yù)后和治療并沒有明顯改善，其5 年生存率僅為5%～10%，晚期膽管癌的中位生存期不超過(guò)12 個(gè)月［1］。因此，尋找影響膽管癌患者預(yù)后的相關(guān)因素具有重要的臨床意義。美國(guó)癌癥聯(lián)合委員會(huì)（AJCC）的分期手冊(cè)已經(jīng)成為癌癥患者分類、確定預(yù)后和確定最佳治療方法的基準(zhǔn)。然而隨著對(duì)癌癥研究的深入，人們發(fā)現(xiàn)TNM 分期系統(tǒng)并不足以評(píng)估預(yù)后，也不能反映癌癥的生物學(xué)異質(zhì)性。即使是同一分期的患者，其預(yù)后和治療反應(yīng)也有很強(qiáng)的異質(zhì)性，其他因素如年齡、表現(xiàn)狀況、腫瘤位置等也會(huì)影響患者的生存，因此它對(duì)臨床預(yù)后提供的信息有限［2］。因此，我們迫切需要開發(fā)可靠的預(yù)后生物標(biāo)志物，以提供更好的臨床預(yù)后價(jià)值。自噬是通過(guò)吞噬細(xì)胞質(zhì)蛋白或細(xì)胞器及其內(nèi)含物進(jìn)入囊泡，并與溶酶體融合，降解內(nèi)含物的內(nèi)容，實(shí)現(xiàn)細(xì)胞的代謝需要和某些細(xì)胞器的更新，從而形成自噬溶酶體的動(dòng)態(tài)穩(wěn)態(tài)過(guò)程［3］。這種分解過(guò)程涉及多種自噬相關(guān)基因（ARGs）［4］。在炎癥、神經(jīng)退行性變、衰老、腫瘤等病理?xiàng)l件下，自噬可以被抑制，以維持細(xì)胞內(nèi)穩(wěn)態(tài)［5］?；贏RGs構(gòu)建的預(yù)后模型已應(yīng)用于多種癌癥類型［6-7］。2020年8月—2020 年9 月，本研究通過(guò)對(duì)TCGA 數(shù)據(jù)庫(kù)中CHOL 隊(duì)列的轉(zhuǎn)錄組和臨床數(shù)據(jù)進(jìn)行生物信息的綜合分析，對(duì)膽管癌組織中差異表達(dá)的ARGs 進(jìn)行篩選，并利用Lasso 和Cox 回歸構(gòu)建基于ARGs 的預(yù)后預(yù)測(cè)模型，為膽管癌患者的預(yù)后風(fēng)險(xiǎn)分層、預(yù)后評(píng)估及治療策略的選擇提供新的參考指標(biāo)。

1 資料與方法

1.1 人類ARGs 集的獲取 從人類自噬數(shù)據(jù)庫(kù)（HADb，http：//autophagy. lu/clustering/index. html）中獲取232 個(gè)ARGs，從分子特征數(shù)據(jù)庫(kù)v6.2（MSig-DB，http//software. broadinstitute. org/gsea/msigdb）的GO_AUTOPHAGY 基因集中獲取394 個(gè)ARGs。將兩個(gè)數(shù)據(jù)庫(kù)的ARGs 合并，刪除重疊的ARGs，最后獲得531個(gè)ARGs。

1.2 原始資料的下載與整理 在癌癥基因組圖譜數(shù)據(jù)庫(kù)（TCGA，https：//portal. gdc. cancer. gov/）中選擇CHOL 隊(duì)列的轉(zhuǎn)錄組和臨床數(shù)據(jù)進(jìn)行下載。包含膽管癌組織36 例和正常膽管組織9 例。應(yīng)用Perl5.28.1 軟件將原始測(cè)序數(shù)據(jù)進(jìn)行合并，轉(zhuǎn)換基因的ID，并提取所有ARGs 的表達(dá)數(shù)據(jù)，若有多個(gè)探針檢測(cè)同一個(gè)基因的表達(dá)量，則取該基因表達(dá)量的均值作為該基因的表達(dá)值。本研究納入分析的臨床數(shù)據(jù)完整，未刪減患者的臨床信息。

1.3 膽管癌組織中差異ARGs 的篩選及其GO 功能富集和KEGG 信號(hào)通路分析 應(yīng)用R 3.6.1 軟件的“l(fā)imma”包對(duì)36 例膽管癌組織和9 例正常組織的ARGs 進(jìn)行表達(dá)差異分析，篩選標(biāo)準(zhǔn)為：FDR＜0.05，∣logFC∣≥2，其中FDR 為錯(cuò)誤發(fā)現(xiàn)率，F(xiàn)C 為差異倍數(shù)。利 用R 軟 件 的“clusterprofiler”、“org. HS. eg.db”、“enrichplot”、“ggplot2”包進(jìn)行基因本體（GO）功能富集，包括生物過(guò)程、細(xì)胞成分和分子功能。采用京都基因百科全書（KEGG）進(jìn)行信號(hào)通路分析，并進(jìn)行可視化。

1.4 膽管癌患者預(yù)后相關(guān)關(guān)鍵ARGs 的篩選 將膽管癌患者的生存信息與差異表達(dá)的ARGs 合并在一起，通過(guò)R 語(yǔ)言的“survival”包首先進(jìn)行單因素Cox 回歸分析，計(jì)算每個(gè)差異ARGs 與膽管癌患者生存的風(fēng)險(xiǎn)比（HR）和P 值，P＜0.05 的差異ARGs 被選出來(lái)進(jìn)行下一步分析。應(yīng)用Lasso 回歸分析以減少基因之間共線性的影響，防止后續(xù)構(gòu)建的預(yù)后風(fēng)險(xiǎn)模型變量過(guò)度擬合，通過(guò)構(gòu)建懲罰函數(shù)對(duì)自變量的回歸系數(shù)進(jìn)行壓縮，實(shí)現(xiàn)基因數(shù)據(jù)的降維，進(jìn)而獲得與膽管癌患者預(yù)后相關(guān)度更高的ARGs。最后，對(duì)Lasso 回歸篩選得到的ARGs 進(jìn)行多因素Cox 回歸分析，計(jì)算每個(gè)ARGs 的多因素回歸系數(shù)，得到與患者預(yù)后相關(guān)的關(guān)鍵ARGs。

1.5 風(fēng)險(xiǎn)預(yù)后模型的建立與分析 根據(jù)Cox 和Lasso回歸篩選得到的關(guān)鍵ARGs，構(gòu)建基于ARGs表達(dá)的風(fēng)險(xiǎn)評(píng)分方程［8-9］：風(fēng)險(xiǎn)評(píng)分Coefi×Xi，其中Coef是指基因在多因素Cox 回歸分析中的回歸系數(shù)，X 為基因的表達(dá)量，n 為預(yù)后相關(guān)ARGs 的數(shù)量。根據(jù)公式計(jì)算出每個(gè)患者的風(fēng)險(xiǎn)評(píng)分，取其中位數(shù)作為截?cái)嘀?，將膽管癌患者分為低風(fēng)險(xiǎn)評(píng)分組和高風(fēng)險(xiǎn)評(píng)分組。采用R 軟件“survival”包進(jìn)行Kaplan-Meier 生存曲線分析，應(yīng)用“timeROC”包繪制模型的ROC 曲線，計(jì)算曲線下面積（AUC），評(píng)價(jià)模型的敏感度和特異度。利用“rms”包繪制列線圖，為評(píng)估列線圖實(shí)際生存和預(yù)測(cè)生存的一致性，進(jìn)一步繪制校準(zhǔn)曲線以評(píng)價(jià)模型的可靠性。

2 結(jié)果

2.1 膽管癌組織差異表達(dá)ARGs 的篩選 與正常膽管組織比較，膽管癌組織中有324 個(gè)存在差異表達(dá)的ARGs，其中表達(dá)上調(diào)的基因311 個(gè)、表達(dá)下調(diào)的基因13個(gè)。

2.2 膽管癌組織差異表達(dá)ARGs 功能富集和信號(hào)通路分析結(jié)果 GO功能富集分析顯示，差異表達(dá)的ARGs 在生物學(xué)過(guò)程中包括自噬、利用自噬機(jī)制的過(guò)程、自噬的調(diào)節(jié)、線粒體的拆卸等，在細(xì)胞組分中包括內(nèi)體膜、泡膜、晚期內(nèi)體、自噬小體、溶酶體膜等，在分子功能中包括蛋白絲氨酸/蘇氨酸激酶活性、泛素樣蛋白連接酶結(jié)合、細(xì)胞黏附分子結(jié)合、鈣黏著蛋白的結(jié)合、蛋白激酶調(diào)節(jié)活性、蛋白磷酸酶的結(jié)合等，見表1。KEGG 信號(hào)通路分析顯示，差異表達(dá)的ARGs 主要涉及自噬—?jiǎng)游铩⒕€粒體自噬—?jiǎng)游?、自噬—其他、志賀菌病、凋亡等信號(hào)通路，見表2。

表1 GO功能富集分析（前10位）

表2 KEGG信號(hào)通路分析（前10位）

2.3 膽管癌患者預(yù)后相關(guān)關(guān)鍵ARGs 的篩選結(jié)果將324個(gè)ARGs納入單因素Cox回歸分析，以P＜0.05為篩選條件，共發(fā)現(xiàn)12個(gè)基因可能與膽管癌患者的生存相關(guān)，分別是VPS25、VPS11、EVA1A、BNIP3、FAM13B、PPP1R15A、GABARAP、GMIP、VPS4B、UBC、FXR2、ATG101。Lasso 回歸將上述單因素Cox篩選出來(lái)的ARGs 進(jìn)行再次篩選，結(jié)果見圖1A，并使用交叉驗(yàn)證建立模型，見圖1B。最后納入下一步 分 析 模 型 的ARGs 有VPS11、EVA1A、BNIP3、GABARAP、VPS4B。

2.4 膽管癌患者預(yù)后預(yù)測(cè)模型的構(gòu)建結(jié)果 將上述篩選出來(lái)的5 個(gè)ARGs 納入多因素Cox 回歸分析，最終得到這5 個(gè)關(guān)鍵ARGs 組成的膽管癌患者預(yù)后預(yù)測(cè)模型，根據(jù)5 個(gè)ARGs 的mRNA 表達(dá)水平和風(fēng)險(xiǎn)系數(shù)計(jì)算每個(gè)患者的風(fēng)險(xiǎn)評(píng)分。風(fēng)險(xiǎn)評(píng)分=（-3.739×VPS11）+（1.691×EVA1A）+（1.734×BNIP3）+（5.776×GABARAP）+（-1.310×VPS4B）。取風(fēng)險(xiǎn)評(píng)分的中位數(shù)0.96為截?cái)嘀?，?6例患者分為低風(fēng)險(xiǎn)組（風(fēng)險(xiǎn)評(píng)分＜0.96，18 例）和高風(fēng)險(xiǎn)組（風(fēng)險(xiǎn)評(píng)分＞0.96，18例）。

2.5 膽管癌患者預(yù)后預(yù)測(cè)模型的預(yù)測(cè)價(jià)值 Kaplan-Meier生存分析顯示，高風(fēng)險(xiǎn)組總的生存時(shí)間低于低風(fēng)險(xiǎn)組（P＜0.05），見圖2A。ROC 曲線分析顯示，預(yù)測(cè)1、2、3年生存率的ROC曲線下面積（AUC）分別為0.973、0.989、0.941，均大于0.9。見圖2B。

圖1 Lasso篩選ARGs和交叉驗(yàn)證結(jié)果

圖2 生存曲線與ROC曲線

2.6 膽管癌患者預(yù)后預(yù)測(cè)模型的可靠性 利用R軟件構(gòu)建基于5 個(gè)ARGs 的列線圖，見圖3A。該列線圖模型經(jīng)評(píng)價(jià)得到C指數(shù)為0.822，95%可信區(qū)間為0.721～0.924，表明該模型具有較好的區(qū)分能力。為評(píng)估列線圖實(shí)際生存和預(yù)測(cè)生存的一致性，我們繪制預(yù)測(cè)1年、2年、3年生存率的校準(zhǔn)曲線幾乎都落在了45°的對(duì)角線上（見圖3B、3C、3D），提示該模型準(zhǔn)確性較高。

3 討論

目前對(duì)于膽管癌發(fā)生發(fā)展的分子機(jī)制尚不明確。近年研究發(fā)現(xiàn)，自噬參與了膽管癌的發(fā)生發(fā)展［10-11］。然而自噬在腫瘤發(fā)生中的作用目前仍存在爭(zhēng)議，對(duì)于自噬是否影響膽管癌患者的預(yù)后，目前尚無(wú)定論。近年來(lái)，基因簽名經(jīng)常用于預(yù)測(cè)各種腫瘤的預(yù)后，在一定程度上甚至優(yōu)于TNM 分期和組織病理學(xué)診斷［12-13］。基于ARGs 表達(dá)的預(yù)后模型已被報(bào)道用于多種癌癥，如肺癌、胃癌、腸癌、乳腺癌、膠質(zhì)瘤等［6，14-17］。本研究中，我們通過(guò)Cox和Lasso回歸模型分析了531 個(gè)ARGs，最終篩選出5 個(gè)預(yù)后相關(guān)的關(guān)鍵ARGs，并成功構(gòu)建了膽管癌預(yù)后預(yù)測(cè)模型，結(jié)果顯示高風(fēng)險(xiǎn)組總的生存時(shí)間明顯低于低風(fēng)險(xiǎn)組，預(yù)測(cè)膽管癌患者1 年、2 年、3 年生存率的ROC 曲線下面積AUC 均大于0.9，提示該模型是具有較好敏感性和特異性的理想模型。

本研究中，我們?cè)赥CGA 數(shù)據(jù)庫(kù)中對(duì)膽管癌的ARGs 進(jìn)行差異表達(dá)分析，通過(guò)GO 功能富集和KEGG 信號(hào)通路分析發(fā)現(xiàn)，差異表達(dá)的ARGs主要涉及的生物學(xué)過(guò)程和信號(hào)通路有自噬、凋亡、線粒體、志賀菌病等。我們篩選出的5 個(gè)預(yù)后相關(guān)的關(guān)鍵ARGs，即 VPS11、EVA1A、BNIP3、GABARAP、VPS4B，均被報(bào)道與腫瘤的發(fā)生發(fā)展相關(guān)，涉及多種腫瘤細(xì)胞的自噬、凋亡、鉑類藥物耐藥等生物學(xué)過(guò)程。PETERSON 等［18］報(bào)道，VPS11 基因與黑色素瘤細(xì)胞的順鉑耐藥有關(guān)。SHEN 等［19］報(bào)道，EVA1A 誘導(dǎo)的自噬和細(xì)胞凋亡在抑制GBM 的發(fā)展中起作用，在3 種腦膠質(zhì)瘤細(xì)胞系（U251、U87 和SHG44）中，EVA1A 過(guò)表達(dá)通過(guò)激活自噬和誘導(dǎo)細(xì)胞凋亡，EVA1A 過(guò)表達(dá)也成功抑制了NOD/SCID 小鼠的腫瘤生長(zhǎng)。EVA1A 也被證實(shí)與肝癌、多發(fā)性骨髓瘤、乳腺癌有關(guān)［20-22］。LI 等［23］報(bào)道，lncRNA DGCR5/miR-27a-3p/BNIP3 軸可通過(guò)調(diào)節(jié)p38/MAPK 通路促進(jìn)胰腺癌細(xì)胞凋亡。BORTNIK 等［24］報(bào)道，GABARAP 高表達(dá)與所有乳腺癌亞型的侵襲性疾病表型的不良預(yù)后和臨床病理特征相關(guān)。SZYMANSKA 等［25］發(fā)現(xiàn)，結(jié)直腸癌組織中VPS4B mRNA 和蛋白表達(dá)下調(diào)，敲除VPS4B 基因可誘導(dǎo)細(xì)胞死亡。以上研究表明，這些關(guān)鍵ARGs 在腫瘤發(fā)生發(fā)展中起重要作用。然而，我們對(duì)這些ARGs 在膽管癌中的了解仍然很有限，對(duì)于自噬如何影響膽管癌患者預(yù)后的分子機(jī)制及其對(duì)膽管癌患者臨床治療的意義有待進(jìn)一步研究。

圖3 列線圖和校準(zhǔn)曲線

準(zhǔn)確判斷癌癥患者的預(yù)后有助于指導(dǎo)臨床醫(yī)師進(jìn)行治療決策。近年來(lái)，列線圖作為腫瘤預(yù)后評(píng)估的工具在醫(yī)學(xué)研究和實(shí)踐中得到了廣泛的應(yīng)用［26-28］。列線圖的原理是通過(guò)多因素回歸模型（例如Cox、Logistic 回歸等），把復(fù)雜的回歸方程轉(zhuǎn)變?yōu)榭梢暬膱D形后，臨床醫(yī)師可以根據(jù)模型中每個(gè)風(fēng)險(xiǎn)變量的評(píng)分，然后相加后得到總評(píng)分，很方便的計(jì)算出疾病發(fā)生的概率及判斷患者預(yù)后的好壞。模型的區(qū)分能力常常用C 指數(shù)表示，其小于0.65 表示區(qū)分能力較差，0.65～0.75 表示模型的區(qū)分能力一般，大于0.75 則表示模型的區(qū)分能力較好［29］。我們模型的C指數(shù)為0.822，提示該模型對(duì)預(yù)后有較好的區(qū)分能力。除了區(qū)分能力外，校準(zhǔn)曲線也是評(píng)價(jià)模型準(zhǔn)確性的重要指標(biāo)，它能夠反映模型的預(yù)測(cè)風(fēng)險(xiǎn)與實(shí)際發(fā)病風(fēng)險(xiǎn)的一致程度，理想模型的校準(zhǔn)曲線剛好落在45°度的對(duì)角線上，校準(zhǔn)曲線越匹配對(duì)角線，預(yù)測(cè)準(zhǔn)確性越高［30］。從圖3B可以看出，我們構(gòu)建模型的校準(zhǔn)曲線幾乎都落在了對(duì)角線上，提示該模型的準(zhǔn)確性較高。然而，本研究也存在一定的局限性，我們的模型缺少獨(dú)立的外部數(shù)據(jù)集進(jìn)行驗(yàn)證。另外，模型真實(shí)的評(píng)估性能需要在實(shí)踐中進(jìn)行驗(yàn)證。

綜上所述，膽管癌中差異表達(dá)的ARGs 主要涉及自噬、凋亡、mTOR 信號(hào)通路、鉑耐藥、PI3K-Akt 信號(hào)通路、HIF-1 信號(hào)通路等癌癥相關(guān)通路，通過(guò)篩選與預(yù)后相關(guān)的ARGs，我們成功構(gòu)建了基于VPS11、EVA1A、BNIP3、GABARAP、VPS4B 共5 個(gè)ARGs 表達(dá)的風(fēng)險(xiǎn)預(yù)測(cè)模型，該模型具有較高敏感度和特異度，可有效預(yù)測(cè)膽管癌患者的預(yù)后，對(duì)指導(dǎo)臨床治療、制定精準(zhǔn)治療方案具有重要價(jià)值。