關于吉西他濱耐藥在浸潤性膀胱癌上的機制研究

江詩睿，吳張松，張璐，吳松，任力杰(通信作者*)

(1.安徽醫科大學深圳二院臨床學院，廣東 深圳 518035；2.深圳大學泌尿外科研究所，廣東 深圳 518116; 3.深圳大學第三附屬醫院，廣東 深圳 518001；4.深圳大學第一附屬醫院，廣東 深圳 518035)

0 引言

膀胱癌是全球最流行的癌癥之一。依據腫瘤惡性程度，膀胱癌可分為非肌層浸潤性膀胱癌和肌層浸潤性膀胱癌。大約有 25%的患者被診斷為肌層浸潤性膀胱癌，并伴有或以后發展為該疾病的轉移形式。全身化療常用于肌層浸潤性和轉移性膀胱癌中，以分別減少疾病的遠處傳播和作為姑息治療手段[1-2]。近來，各種基因共表達分析方法越來越廣泛地應用于生物信息學。基因共表達分析是指一組表達相關且相互調節的基因，使用大量基因表達數據構建基因間的相關性，從而挖掘基因功能的一類分析方法。在很多情況下，有著相似行為/變化的物質，會存在著一定的聯系。在不同條件下，可以根據相似的基因表達譜，數學特征和生物學意義將基因分為不同的模塊。同一模塊中的基因被認為具有相似的功能。在它們當中，WGCNA是一種系統的基因聚類方法，用于闡明疾病發生和發展的潛在分子機制。共表達網絡的構建有助于篩選生物標志物和治療靶標。這些方法已成功地用于各種研究中，例如頭頸鱗狀細胞癌的神經周圍浸潤的遺傳分析。WGCNA 揭示了每個模塊中高度連接的基因的分層組織，在模塊化網絡中可以發現參與其中的關鍵hub 基因[3]。

1 資料與方法

1.1 數據來源

本研究于2019年10月從TCGA數據庫中下載膀胱癌表達譜數據，該數據包括412浸潤性膀胱癌樣本[4]。下載以RPKM呈現形式原始表達譜數據，用于篩選差異表達的基因及后續WGCNA分型分析。另外，對應臨床相關資料，包括無病生存(disease free progression，DFS)數據可從網頁下載(http：//www.cbioportal.org)。浸潤性膀胱癌細胞系的 IC50 值來自腫瘤藥物敏感性基因組學數據庫(Genomics of Drug Sensitivity in Cancer，DSC)。

1.2 納入和排除標準

根據下載的412浸潤性膀胱癌的臨床及用藥數據，有83例接受了吉西他濱治療，而239例患者未使用該藥治療。根據患者中有無進展生存數據記錄的，排除使用吉西他濱之前就發生腫瘤進展的病例(N=15)。去除吉西他濱治療期間接受放射治療的病例(N=3)。排除診斷浸潤性膀胱癌后一年內未接受吉西他濱治療的病例(N=5)。最終，有60例符合入組條件。最后按照DFS的中位數，分別把接受吉西他濱標準化治療組(N=60)及未接受吉西他濱治療組(N=239)分成高低DFS組進行差異分析。

1.3 GO功能富集和KEGG信號通路分析

利用R軟件的“clusterprofiler”程序包GO及KEGG功能富集，包括基因執行的分子功能(Molecular Function，MF)，基因所處的細胞組分(Cellular Component，CC)，基因以及參與的生物學過程(Biological Process，BP)這三部分[5]。

1.4 加權基因共表達網絡分析(WGCNA)

利用 R 軟件的“WGCNA” 程序包構建加權基因共表達網絡[6]。首先，通過方差分析確定的 TOP 20%的變異基因(4101個基因)作為候選基因，用作后續 WGCNA 分析。繪制患者的層次聚類樹以查找和刪除異常值。然后，根據標準無標度網絡選擇軟閾值β=7。最后，根據每個基因之間的相關性，進行平均連鎖層次聚類，將變異基因分為不同的模塊。

1.5 對吉西他濱有特殊反應的基因功能模塊的選擇

根據每個模塊中基因表達水平的高低，將患者分為高表達組和低表達組。根據無進展生存時間繪制生存曲線。選擇具有顯著意義的基因功能模塊(P＜0.05)作為吉西他濱特殊反應模塊。

1.6 共表達網絡的可視化和 hub 基因的定義

使用Cytoscape 3.7.2軟件對2.6部分篩選目標功能模塊行基因共表達網絡可視化[7]。hub基因是由Cytoscape的插件CytoHubba計算的，它提供了四種類型的基于本地的方法和七種類型的基于全局的方法[8]。每個模塊中的高度連接的基因是通過四種基于本地的方法來計算的，包括Degree，Maximum Neighborhood Component (MNC)，Density of Maximum Neighborhood Component (DMNC)，and Maximal Clique Centrality (MCC)。用以上四種方法的交集獲得的基因被鑒定為hub基因或關鍵控制基因，用韋恩圖可視化。

1.7 基于 mRNA 表達的預后標記的構建和準確性驗證

對所選的 24 個 hub 基因進行Lasso回歸分析和單變量Cox 回歸分析。由此，我們確定了三個與DFS顯著相關的基因(P＜0.05)。使用以下公式計算簽名的風險評分：其中 Coefi是系數，xi 是每個選定基因的 z 分數轉換的相對表達值。利用 R 軟件中的“survival”和“survminer”程序包探索最佳風險分值并繪制 Kaplan-Meier 生存曲線。重要的是使用“survminer”R軟件包的“surv_cutpoint”功能來將患者分為高風險和低風險組的最佳臨界值，并使用R包“FactoMineR”和“factoextra”行主成分分析來驗證分組結果。此外，使用Cox回歸分析評估風險評分與DFS之間的關聯，其中年齡，階段和性別為協變量。

1.8 基因功能驗證和富集分析

從GDSC數據庫獲得基因表達數據和藥物相關信息，并探索吉西他濱的IC50與浸潤性膀胱癌細胞系中風險基因表達的相關性[9]。進行GSEA分析來研究與浸潤性膀胱癌不同亞組相關的功能[10]。

1.9 數據分析

所有統計分析使用GraphPadPrism8.0和使用R軟件(version 3.5，www.r-project.org)進行。所有數據均表示為(±SD)。本研究中使用了 Student’st-test 和 Pearson’s χ2test來分析數據，當P＜0.05具有顯著性差異。

2 結果

2.1 分組情況

在生存分析中，死亡或腫瘤轉移以及任何原因引起的復發均被視為終末事件。在無進展生存數據記錄的患者中，經過吉西他濱治療的患者共83例，未使用該藥的患者為239例。排除使用吉西他濱之前發生腫瘤進展的病例(N=15)，去除吉西他濱治療期間接受放射治療的病例(N=3)，排除診斷后一年內未接受吉西他濱治療的病例(N=5例)。最終，篩選出60例符合條件患者。根據DFS的中位數，將60例吉西他濱治療的患者分為兩組，在較長生存期組中有30例患者，在較短生存期組中有30例患者。其中，無病生存期相對較長的組中，有8人的生存期超過5年，有3人的生存期超過10年(圖1A)。根據中位無病生存期，將未經吉西他濱治療的239例患者分為兩組，較長生存期組120例，較短生存期組119例。其中，無病生存期相對較長的組中，有28人的生存期超過5年，有2人的生存期超過10年(圖1B)。

圖1 吉西他濱給藥功效的差異具有其獨特的差異表達基因

圖2 加權基因共表達網絡分析(WGCNA)

2.2 吉西他濱共表達網絡的構建

通過R軟件的差異分析，我們發現了數量有限的差異表達基因。富集的KEGG途徑與免疫細胞，造血細胞和細胞因子有關。這些途徑對此問題都有一定的解釋，但是這些結果仍然相對分散。為了更好地解釋該問題，我們使用了WGCNA來構建對吉西他濱有特殊反應的基因共表達網絡。先前選擇的接受吉西他濱治療的患者的分層聚類樹(N=60)顯示，所有樣本均被聚集成一個大的聚類，而沒有需要刪除的離群值。樣本之間最變異的基因的前20%通過平均連鎖層次聚類被分組為基因功能模塊。如圖2C所示，總共鑒定出15個基因功能模塊。未分配的基因被分類為灰色模塊。在所有模塊中，最大的模塊是綠色模塊，其中包含338個基因。最小的模塊是深藍色模塊，僅包含42個基因。不同模塊基因表達的KM生存分析表明，黃色模塊(P=0.0084)和綠色模塊(P=0.032)與DFS的相關性最高。我們通過Cytoscape插件CytoHubba識別了24個hub基因(分別在綠色模塊中選擇了10個，在黃色模塊中選擇了14個)，其特征是對網絡中的節點進行排名。對于黃色模塊，GO和KEGG分析顯示主要途徑顯著富集了粘著斑，ECM受體相互作用，跨膜受體蛋白，細胞外結構組織，膠原生物合成過程等。hub 基因是 COL3A1，HTRA1，THY1，ADAM12，COL5A2，ITGA11，OLFML2B，SP ARC，FBN1，ANGPTL2，CDH11，COL1A1，CTSK，PXDN。該模塊包含差異表達的基因 PCOLCE 和 EFEMP2，它們在圖3D 中以黑色標記。對于綠色模塊，GO 和 KEGG 分析顯示主要途徑顯著富集于血管平滑肌收縮，肌動蛋白細胞骨架的調節，粘著斑，組織遷移，平滑肌收縮。hub 基因是 CNN1，ALDH1B1，ACTG2，KANK2，ACTA2，MYLK，TAGLN，CSRP1，TPM1，FLNA。該模塊包含差異表達的基因 FLNC，其在圖 4F 中以黑色標記。

圖3 對吉西他濱的特殊應答基因共表達網絡的構建

2.3 吉西他濱在浸潤性膀胱癌治療的風險評估的預測價值

為了更好地了解尿路上皮癌中與吉西他濱耐藥相關的潛在調節因子，我們鑒定了 24個 hub 基因。接下來，我們對24 個 hub 基因進行了單變量 Cox 回歸分析。結果表明，浸潤性膀胱癌患者中 KANK2，ITCA11，FBN1，PXDN，COL5A2和 THY1 的高表達具有較差的 DFS(圖 4A)。此外，我們還將 LASSO Cox 回歸算法應用于 24 個 hub 基因，以便準確識別涉及臨床結果的調節基因。根據最低標準選擇了五個基因 (ALDH1B1，COL5A2，FLNA，ITGA11 和 KANK2)(圖 4B，C)。最終，我們選擇 KANK2，ITCA11 和 COL5A2 進行多元Cox回歸分析，并使用從多元 Cox 回歸算法獲得的系數來計算風險評分，以建立浸潤性膀胱癌患者吉西他濱反應的預測模型(圖4D)，并指出 KANK2 可以作為不良 DFS 的獨立危險因素。

2.4 預后風險評分模型與浸潤性膀胱癌臨床病理表現出顯著相關性

為了研究三基因風險因子的預后模型的有效性，我們根據風險評分中位數將浸潤性膀胱癌患者分為高風險(吉西他濱耐藥)和低風險(吉西他濱敏感)組，三基因預測模型風險評分分布如圖5A 所示。此外，計算了風險評分與浸潤性膀胱癌狀況之間的關系(圖5B)。此外，結果表明，高風險組浸潤性膀胱癌患者的生存率較差(圖 5C)。為了判斷我們的分類是否正確，我們將通過 PCA 分析這兩個子類，結果顯示高風險亞組可以聚集在一起，低風險亞組也可以聚集在一起(圖5D)，并且在 3D 模型中得到驗證。接下來，風險評分模型與其他臨床數據結合，我們進行了單因素和多因素 Cox 回歸分析，以確定風險特征是否為獨立的預后指標。根據證據，危險評分與 DFS 顯著相關 (圖5E，F)。

2.5 浸潤性膀胱癌中吉西他濱耐藥的三基因風險分子特征和潛在的信號通路

我們首先顯示了高風險和低風險組(圖6A)中三種基因風險分子(KANK2，ITCA11和COL5A2)的表達水平。與DFS較短的組相比，DFS較短的組患者通常更容易表達的KANK2(圖6B)。為了研究 KANK2 對膀胱癌細胞化療反應的影響，我們利用高等級膀胱癌細胞系癌癥項目數據集(GDSC)中的藥物敏感性基因組學，分析 KANK2 表達水平與對化療敏感性分析。有趣的是，如圖6C 所示，KANK2 水平顯示出與吉西他濱的 IC50呈正相關的趨勢，但是P值不顯著。此外，為揭示吉西他濱在浸潤性膀胱癌中低應答的潛在機制，GSEA 揭示了幾種眾所周知的信號通路可能與吉西他濱耐藥相關，包括上皮間充質轉化 (EMT)(NES=1.95，FDRq＜0.001)，干細胞分化 (NES=1.47，FDRq=0.021)，hedgehog(Hh)信號傳導 (NES=1.57，FDRq=0.031)和Wnt信號通路(NES=1.67，FDRq=0.006)與高危亞組顯著相關(圖6D)。所有，這些發現為吉西他濱在浸潤性膀胱癌中的耐藥性提供了潛在的靶點，需要將來進一步研究。

3 討論

圖5 預后風險評分顯示浸潤性膀胱癌與臨床病理特征密切相關

我們根據來自TCGA數據庫的基因表達譜數據，通過WGCNA鑒定了與吉西他濱耐藥相關的兩個不同模塊，即綠色和黃色模塊。通過利用Lasso回歸分析和Cox回歸分析，我們還從綠色和黃色模塊中選擇了三個hub基因(例如，COL5A2，ITGA11和KANK2)得出了預后風險特征，從而將浸潤性膀胱癌患者分為高風險和低風險亞組。此外，高風險和低風險亞組不僅與患者疾病預后相關，而且與惡性膀胱癌中吉西他濱耐藥相關的關鍵信號通路密切相關。我們當前工作的可以開發專家決策支持系統，該系統可以特別預測患者對吉西他濱的敏感性，并指導臨床醫生選擇吉西他濱進行化學治療的最佳人選。此外，我們推斷靶向EMT，Hh和Wnt信號通路可以用來改善患者的吉西他濱治療結局，這顯然是未來研究的方法。這項研究是從臨床到分子水平全面分析吉西他濱在浸潤性膀胱癌患者中療效差異的因素。這是從組學角度首次驗證影響吉西他濱化療藥物治療浸潤性膀胱癌療效的各種因素。但是，考慮到所涉及患者的數量有限以及缺乏相應實驗數據來支持我們的結論，未來需要進一步研究以證明其可行性。更精確地來講，這些治療方案可以通過首先建立患者來源的 3D 類器官和/或在異種移植模型中進行驗證，顯示對個體化給藥的持續監測及療效驗證。

圖6 浸潤性膀胱癌中吉西他濱耐藥性涉及的三基因風險特征和潛在的信號通路