基于TCGA數據庫建立膽管癌自噬相關基因預后預測模型及其應用

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

廣東醫科大學附屬醫院，廣東湛江524000

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

1 資料與方法

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

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

1.3 膽管癌組織中差異ARGs 的篩選及其GO 功能富集和KEGG 信號通路分析 應用R 3.6.1 軟件的“limma”包對36 例膽管癌組織和9 例正常組織的ARGs 進行表達差異分析，篩選標準為：FDR＜0.05，∣logFC∣≥2，其中FDR 為錯誤發現率，FC 為差異倍數。利 用R 軟 件 的“clusterprofiler”、“org. HS. eg.db”、“enrichplot”、“ggplot2”包進行基因本體（GO）功能富集，包括生物過程、細胞成分和分子功能。采用京都基因百科全書（KEGG）進行信號通路分析，并進行可視化。

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

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

2 結果

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

2.2 膽管癌組織差異表達ARGs 功能富集和信號通路分析結果 GO功能富集分析顯示，差異表達的ARGs 在生物學過程中包括自噬、利用自噬機制的過程、自噬的調節、線粒體的拆卸等，在細胞組分中包括內體膜、泡膜、晚期內體、自噬小體、溶酶體膜等，在分子功能中包括蛋白絲氨酸/蘇氨酸激酶活性、泛素樣蛋白連接酶結合、細胞黏附分子結合、鈣黏著蛋白的結合、蛋白激酶調節活性、蛋白磷酸酶的結合等，見表1。KEGG 信號通路分析顯示，差異表達的ARGs 主要涉及自噬—動物、線粒體自噬—動物、自噬—其他、志賀菌病、凋亡等信號通路，見表2。

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

表2 KEGG信號通路分析（前10位）

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

2.4 膽管癌患者預后預測模型的構建結果 將上述篩選出來的5 個ARGs 納入多因素Cox 回歸分析，最終得到這5 個關鍵ARGs 組成的膽管癌患者預后預測模型，根據5 個ARGs 的mRNA 表達水平和風險系數計算每個患者的風險評分。風險評分=（-3.739×VPS11）+（1.691×EVA1A）+（1.734×BNIP3）+（5.776×GABARAP）+（-1.310×VPS4B）。取風險評分的中位數0.96為截斷值，將36例患者分為低風險組（風險評分＜0.96，18 例）和高風險組（風險評分＞0.96，18例）。

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

圖1 Lasso篩選ARGs和交叉驗證結果

圖2 生存曲線與ROC曲線

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

3 討論

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

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

圖3 列線圖和校準曲線

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

綜上所述，膽管癌中差異表達的ARGs 主要涉及自噬、凋亡、mTOR 信號通路、鉑耐藥、PI3K-Akt 信號通路、HIF-1 信號通路等癌癥相關通路，通過篩選與預后相關的ARGs，我們成功構建了基于VPS11、EVA1A、BNIP3、GABARAP、VPS4B 共5 個ARGs 表達的風險預測模型，該模型具有較高敏感度和特異度，可有效預測膽管癌患者的預后，對指導臨床治療、制定精準治療方案具有重要價值。