关于吉西他滨耐药在浸润性膀胱癌上的机制研究

江诗睿，吴张松，张璐，吴松，任力杰(通信作者*)

(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 浸润性膀胱癌中吉西他滨耐药性涉及的三基因风险特征和潜在的信号通路