基于TCGA數據庫膽管癌自噬-臨床預后模型的構建和評估

鄒文強，符廣華，楊艷梅，張新宇

（1.哈爾濱醫科大學附屬第二醫院 普通外科，2.哈爾濱醫科大學 腫瘤防治研究所，黑龍江 哈爾濱 150086）

膽管癌（cholangiocarcinoma，CCA）是起源于膽管上皮細胞的惡性腫瘤，根據病變的解剖部位可以分為肝內CCA、肝周CCA和遠端CCA[1]。CCA的發病率約占消化道腫瘤的3%，我國及東南亞地區是最高發的地區，而且近幾十年發病率呈上升趨勢[2]。由于CCA患者早期發病癥狀隱匿且惡性程度較高，一旦發現，多數腫瘤已進入晚期。研究表明根治性手術是最有效的治療方式，但由于淋巴結轉移早，易侵犯血管及神經，術后遠期預后較差[3]。CCA患者的預后評估方法主要是TNM分期和組織學診斷，但無法實現準確的個體化生存預測[4-6]；且臨床常規的腫瘤標志物往往敏感度和特異性不高，因此，尋找可靠的預后標志物來優化預后評估系統是必不可少的。

近年隨著對腫瘤的深入研究，發現自噬在腫瘤的發生發展中發揮著重要作用[7-9]。自噬是一種通過溶酶體自我分解細胞內的老化及受損的蛋白質和細胞器，實現細胞內穩態的途徑[10]。細胞自噬在正常生理狀態下可以促進機體抵御疾病，但當自噬異常時又會導致腫瘤、自身免疫性及神經退行性等疾病[11]，在腫瘤中的作用取決于組織類型、腫瘤分期和致癌突變類型等因素[12]。自噬對包括CCA在內的多種腫瘤，具有雙重調節作用。自噬通過消除受損的線粒體和減少ROS介導的染色體損傷，在腫瘤初期發揮抑癌功能。自噬通過抑制p53 的誘導和維持線粒體的代謝功能，在腫瘤進展期促進腫瘤細胞存活[13-15]。Sasaki等[16]的研究表明，CCA組織樣本中自噬相關蛋白（LC3、beclin-1）表達量顯著增加。CCA進展及轉移的重要驅動因素之一是上皮細胞向間質細胞的轉變（EMT），Chen等[17]發現調節自噬的信號通路（PI3K/AKT、p53和beclin1）同時對EMT有顯著調節的作用，促使CCA細胞存活。這些研究證實自噬和CCA存在密切關系，并表明ATGs作為CCA預后標志物的巨大潛力。既往研究表明依據自噬基因構建的預后風險模型，可以有效評估膀胱癌、乳腺癌及胃癌等患者的預后情況[18-20]。

本研究擬通過癌癥基因組圖譜數據庫（The Cancer Genome Atlas，TCGA）分析與CCA預后密切相關的ATGs，同時構建自噬基因預后風險模型，觀察其對CCA患者的預后價值，從而為臨床診斷和指導治療CCA提供個性化的模型。

1 材料和方法

1.1 數據下載及整理

在TCGA中獲取CCA的轉錄組數據和臨床數據，包括35 例CCA組織樣本和9 例正常組織樣本。使用Strawberry Perl軟件將各個獨立樣本的基因表達量匯總在一個文件中，并將基因的“ID”轉換為“Symbol”。自噬相關基因（ATGs）從人類自噬數據庫（Human Autophagy Database，HADb）中下載，包括232個ATGs的信息。在TCGA數據庫已經下載的CCA樣本信息中，利用Perl腳本將232個ATGs的基因表達量提取出來。

1.2 差異表達的自噬相關基因（DEATGs）的識別

通過R語言（R 4.0.2）使用“Willcox.test”來篩選DEATGs，篩選條件為：FDR＜0.05，|logFC|＞1。通過R語言繪制DEATGs的火山圖和箱線圖。

1.3 DEATGs的富集分析

基因本體功能（gene ontology，GO）富集分析和京都基因與基因組百科全書通路（Kyoto Encyclopedia of Genes and Genomes pathway，KEGG pathway）分析，可以明確DEATGs的主要生物學屬性。調用R語言中“org.Hs.eg.db”“enrichplot”和“gplot2”包，進行GO富集分析和KEGG通路分析；使用R語言中“GOplot”包，繪制GO氣泡圖和KEGG圈圖。

1.4 預后相關ATGs的識別和基于ATGs的預后風險評分模型構建

使用Strawberry Perl軟件將DEATGs表達量和生存狀態及時間，匯總在一個文件中。調用R軟件的“survival”包，篩選標準為P＜0.05，對DEATGs進行單因素Cox回歸分析，篩選出與CCA預后相關的自噬基因（ATGs）。再進行多因素Cox回歸分析，刪除相關性較高的基因來優化模型，最后得到用于構建預后風險評分模型的風險基因。使用R軟件獲取風險基因的評分系數（Coef）。模型的計算公式為：患者風險評分=（基因1的表達量×coef1+基因2的表達量×coef2+…）。按照患者風險評分，以風險值中位數為界限，將患者分為高風險和低風險組。

1.5 預后模型的驗證

使用R 4.0.2軟件進行生存分析、獨立預后分析和臨床相關性分析，繪制風險曲線和多指標ROC曲線。

2 結果

2.1 自噬相關基因的差異分析

從TCGA數據庫中下載了35 個腫瘤組織樣本和9 個正常組織樣本的轉錄組數據和臨床數據。在HADb數據庫中提取了232 個自噬相關基因，并提取自噬相關基因表達量。按照FDR＜0.05和log2FC＞1 的篩選標準，過濾得到49 個自噬相關差異基因（DEATGs），其中上調的基因有48個，下調的基因有1個（圖1）。

圖1 49個差異自噬基因的箱線圖

2.2 差異自噬基因的富集分析

對49個差異表達的自噬相關基因（DEATGs）進行GO富集分析和KEGG通路分析。經過GO富集分析，在生物學過程方面，基因主要富集在自噬、利用自噬機制的過程、宏觀自噬、自噬調控、自噬信號通路的調控、細胞生長及內在凋亡信號通路等；在細胞成分方面，基因主要富集在晚期內涵體、細胞的焦點粘連、核被膜及液泡膜上；在分子功能方面上，基因主要富集在細胞黏附分子結合、泛素樣蛋白連接酶結合、鈣黏素結合及蛋白絲氨酸/蘇氨酸激酶活性等（圖2）。KEGG通路分析如圖3，這些基因主要與自噬-動物、志賀氏癥、人乳頭瘤病毒感染、細胞凋亡、PI3K-AKT信號通路、鉑類耐藥、肝細胞癌、慢性粒細胞白血病、胰腺癌等通路顯著相關。

圖2 GO富集分析的氣泡圖

圖3 KEGG通路分析的圈圖和熱圖

2.3 預后相關自噬基因的篩選

將差異自噬基因的表達量與生存時間和生存狀態合并，按照P＜0.05的篩選標準，通過單因素Cox回歸分析，過濾得到5個與預后相關的差異自噬基因（RHEB、PPP1R15A、ATG101、BNIP3、NRG1），見圖4。

圖4 單因素Cox分析森林圖

2.4 自噬基因的預后模型的構建

將這5個預后相關基因通過多因素Cox回歸分析進行優化，刪除相關性較高基因，結果這5個基因都與患者的預后顯著相關，可以用來構建預后模型（見表1）。自噬基因的預后風險模型公式為（NRG1表達量×1.1207+BNIP3表達量×1.4002+ATG101表達量×1.3457-PPP1R15A表達量×1.1613-RHEB表達量×1.3279）。

表1 多因素Cox分析CCA預后影響因素

2.5 自噬預后模型的驗證

根據模型計算每個患者風險值，18 例屬于高風險組，18 例為低風險組。為了驗證自噬預后模型在預測臨床預后方面的可靠性，進行以下分析。Kaplan-Meier生存分析中P=8.52e-05（P＜0.05），因此高、低風險組的患者在生存率上有差異的可能性較大（圖5），高風險組3 年生存率為18.7%，低風險組3年生存率為70.3%。圖6A為風險曲線圖，從左到右患者的風險值逐漸增大；圖6B為生存狀態圖，描述了風險值與生存時間/生存狀態的關系，隨著風險值增大死亡患者增多，隨著風險值減小死亡患者減少；圖6C為風險熱圖，描述了風險值與基因表達量的關系，隨著風險值的增加，4個基因的表達量增加（RHEB、ATG101、BNIP3、NRG1），表明這4個基因為高風險基因；隨著風險值的增加，1個基因的表達量減少（PPP1R15A），表明這1個基因為低風險基因。在單因素獨立預后分析中，風險值與生存時間及生存狀態顯著相關（HR1.382，95%CI1.160～1.645，P＜0.001）（圖7A），在多因素獨立預后分析中風險值與生存時間及生存狀態也顯著相關（HR1.427，95%CI1.161～1.756，P＜0.001）（圖7B），表明風險值可以作為獨立預后因子，進一步驗證了Cox模型在預測臨床預后方面的可靠性。在多指標ROC曲線中，風險值的AUC值最大（圖8），表明風險得分有更好的預測臨床預后的能力。

圖5 生存曲線分析

圖6 風險值與生存時間/生存狀態、基因表達量的關系

圖7 單因素（A）和多因素（B）分析患者預后影響因素

3 討論

CCA作為臨床常見的惡性腫瘤，預后極差，可能與肝血管的侵犯、遠處轉移、較早淋巴結的轉移及術后早期復發有關，因此，尋找可靠的早期診斷標志物來完善預后評估系統，對于CCA患者預后有極大幫助。本研究在TCGA數據庫中篩選出49個差異表達的自噬相關基因DEATGs，且通過富集分析證實自噬參與CCA發生發展及耐藥的過程，本研究應用單因素及多因素Cox回歸分析最終篩選出5個與CCA預后顯著相關的自噬基因（RHEB、PPP1R15A、ATG101、BNIP3、NRG1）。

RHEB基因屬于RAS家族的一員，編碼的蛋白屬于GTP結合蛋白，在與GTP結合的情況下Rheb被激活，當GTP被水解時Rheb受到抑制。在TSC1-TSC2-TBC1D7/RHEB/mTORC1經典的信號通路調節中，扮演重要角色。在參與調節自噬的發生發展、調控細胞分化及抑制腫瘤發面發揮重要作用[21]。在結腸癌組織中，p27KIP1受到RHEB的調控，RHEB/p27 激活mTORC1 信號通路來增強自噬，促進腫瘤細胞的存活，加速結腸癌的進展[22]。He等[23]的研究表明，FADD基因在乳腺癌中顯著上調，FADD基因敲減降低了乳腺癌細胞系中Rheb的表達，抑制了mTORC1的活性，抑制腫瘤細胞的存活。Chen等[24]在宮頸癌組織樣本中的研究表明，circMYLK與miR-1301-3p結合可以促宮頸癌細胞生長，其機制是RHEB/mTORC1信號通路的激活。

PPP1R15A基因又稱為GADD34，在應激生長停止條件和DNA損傷劑治療后其轉錄水平增加。PPP1R15A編碼的蛋白質可以調節電離輻射后的凋亡，在促進細胞死亡和未折疊蛋白反應中起著重要的作用。Holczer等[25]在肝癌細胞中的研究發現，GADD34的表達被抑制后，導致mTOR的快速激活和自噬水平下調，隨后凋亡細胞死亡。GADD34表達的減少顯著抑制腫瘤，并導致MDSCs和T細胞的積累減少，抑制GADD34減少了MDSCs分泌的血管上皮生長因子α和轉化生長因子β[26]。

ATG101基因位于染色體12q13.13，編碼的蛋白是ULK1 自噬復合物的重要組成部分，與ATG13結合維持ULK1復合物的穩定性，進而調控自噬[27]。研究表明，PTCH1是已知的腫瘤抑制因子，抑制腫瘤細胞增殖和誘導細胞凋亡，減少PTCH1的表達量會增加ATG101介導的自噬通量來促進癌細胞增殖[28]。

BNIP3基因的全稱為BCL2相互作用蛋白3基因，該基因編碼一種線粒體蛋白，含有BH3 結構域，可以充當促凋亡因子。細胞內存在很多蛋白可以誘導線粒體自噬和凋亡，但BNIP3能將線粒體自噬和凋亡串聯在一起，功能更加特殊[29]。m6A去甲基化酶（FTO）在乳腺癌中顯著高表達，BNIP3是其下游靶點，促進腫瘤細胞增殖和轉移[30]。Li等[31]在胰腺癌的研究中發現，BNIP3甲基化抑制了線粒體介導的胰腺癌細胞凋亡的誘導。

NRG1基因編碼一種膜糖蛋白，介導多種細胞信號通路，該基因有多種異構體，該基因的失調與癌癥等疾病密切相關。關于前列腺癌中一項研究發現，腫瘤微環境下產生的NRG1，可以促進腫瘤細胞的耐藥性。研究表明，阻斷NRG1會減少其介導的HER3的激活，抑制NRG1信號通路的誘導，起到抑制胰腺腫瘤細胞（PC）和腫瘤相關成纖維細胞（CAFS）相關的促腫瘤細胞增殖和轉移的作用[32]。Trombetta 等[33]研究發現，肺腺癌中NRG1的高表達會通過PI3K-AKT和MAPK信號通路異常激活ErbB2/ErbB3，使肺腺癌細胞侵襲性增強。

本研究基于上述5 個自噬基因構建預后風險模型，發現風險模型的AUC值（0.906）要明顯優于其他指標，因此可以作為CCA的獨立預后指標。當然本研究還存在一定的局限性，譬如本研究屬于純數據庫分析，缺少相關的功能實驗進行驗證。