基于免疫细胞浸润评分实现膀胱癌分型及预后风险评估

殷桂草，郑生旗，张 伟，董 欣，祁乐中，李一帆

1. 扬州大学附属医院泌尿外科，江苏 扬州 225000

2. 扬州大学护理学院 公共卫生学院，江苏 扬州 225000

膀胱癌具有高度异质性，是世界上第五大常见恶性肿瘤［1］。2018年，全世界超过19.9 万人死于膀胱癌，新发膀胱癌病例超过54.9 万例［2］。约3/4 的膀胱癌在最初诊断时为非肌层浸润性膀胱癌［3］，可以通过经尿道膀胱肿瘤切除联合膀胱内灌注化疗进行治疗，预后一般较好；一旦病变累及膀胱肌层，或出现远处转移，患者的5年总存活率将从80%下降到15%以下［4-5］。因此，迫切需要开发新的预后评估策略和有效的治疗方法来改善膀胱癌患者的预后。

近年来，免疫疗法显示出积极的抗肿瘤作用［6-7］。研究证实，PD-1/PD-L1 检查点抑制剂对肌层侵袭性和转移性膀胱癌具有较好的疗效［8］。肿瘤微环境中的免疫成分在肿瘤细胞基因表达和临床预后中发挥重要作用［9］，不同免疫细胞浸润丰度的膀胱癌患者在基因组、转录组和生物学过程上存在显著异质性［10］。因此，量化肿瘤中关键免疫细胞成分可以为预测肿瘤患者的免疫治疗效果及预后评估提供新的方法。

随着高通量技术的成熟和肿瘤免疫研究算法的发展，转录组测序可以用来描述肿瘤免疫微环境［11］。然而，目前还没有研究系统地探究免疫浸润评分与膀胱癌患者治疗及预后的关系。本研究基于ssGSEA 算法获得16种免疫细胞浸润的评分，通过无监督聚类将患者分为高、低免疫浸润两个聚类，研究不同聚类之间免疫检查点相关基因的表达、化疗药物敏感性和免疫治疗效果的差异；通过WGCNA 发现与关键免疫细胞相关的核心基因，Cox 回归筛选预后相关基因，构建预后风险评分模型，并利用多因素分析发现该风险评分是膀胱癌患者独立的预后因素；最后，根据风险评分和临床病理参数建立列线图来预测患者1、3、5年的总存活率，为患者的预后评估及个性化治疗提供参考。

1 资料与方法

1.1 数据收集及处理

从TCGA 数据库（https：//portal.gdc.cancer.gov/repository）中获取19 份正常组织标本和411份膀胱癌组织标本的RNA 测序数据及相关临床数据。排除重复标本和患者总存活时间不足30 d 的标本。从GEO 数据库（https：//www.ncbi.nlm.nih.gov/gds/）下载GSE13507 中标本的mRNA表达数据和患者的临床数据。从TCIA 数据库（https：//tcia.at/home）下载膀胱癌患者对不同免疫治疗方案的免疫细胞阳性比例分数评分数据，免疫细胞阳性比例分数得分越高代表患者免疫治疗的效果可能越好。

1.2 采用ssGSEA分析免疫细胞浸润情况

利用“GSVA”包进行ssGSEA，计算正常组织和肿瘤组织16 种免疫细胞浸润的评分，并使用R中的“limma”、“ggplot2”和“ggpubr”包对正常组织和肿瘤组织免疫浸润评分的差异进行比较。

1.3 采用无监督聚类进行膀胱癌患者分型分析

为了进一步探讨免疫浸润评分对膀胱癌患者预后和治疗的价值，使用“ConsensusClusterPlus”包基于16 种免疫细胞浸润评分对膀胱癌患者进行聚类分型，设置maxK=9，reps=100，pItem=0.8，pFeature=1，distance=“manhattan”，clusterAlg=“pam”。同时，使用ggplot2 R 包生成无监督聚类和相应的代表性数据。之后，对不同聚类的患者进行预后分析，并进行主成分分析以验证分型的效果。使用“Heatmap”包绘制不同分型下免疫细胞浸润评分的热图。使用“limma”包、“ggplot2”包评估不同分型下免疫检查点相关基因的表达及患者对免疫治疗反应的差异。最后，使用“prophytic”包评估不同分型的膀胱癌患者对临床药物的敏感性。

1.4 采用WGCNA 筛选与关键免疫细胞显著相关的基因集

WGCNA 的输入数据文件为：①TCGA 数据集中411 个膀胱癌肿瘤组织标本的转录组数据；②前述鉴定的关键免疫细胞的免疫浸润评分数据。首先，通过中位数绝对偏差测量411 个标本中基因表达的可变性，并确定在正常组织和肿瘤组织中差异表达的基因。构建带表型热图的样本树突图。然后，通过软阈值函数计算出最佳软阈值，并通过软阈值构造加权邻接矩阵。邻接矩阵进一步转化为拓扑重叠矩阵。然后，将具有相似表达的基因聚类成一个模块，称为共表达模块。计算P值和相关系数以确定共表达模块与表型之间的关联。最后，鉴定出与关键免疫细胞浸润显著相关的关键模块，并提取模块内的关键基因用于下一步研究。

1.5 功能富集分析研究关键基因的潜在功能

使用“ClusterProfiler”包进行富集分析，以探索关键基因潜在的功能，并选择GO 和KEGG 作为背景数据库。

1.6 构建基于关键基因的预测模型并验证

首先，采用“glmnet”包对1.4 筛选出来的基因进行LASSO Cox 回归分析，进一步确定与膀胱癌患者预后相关的基因。然后，从HPA 数据库下载膀胱癌预后相关基因在膀胱癌组织和正常膀胱组织中的免疫组织化学结果，验证其在组织中的蛋白质表达水平。针对无免疫组织化学结果的HOXB5和HOXB6采用qRT-PCR在临床标本中验证。具体方法：收集扬州大学附属医院8 份膀胱癌组织标本和邻近正常组织标本，其中正常组织标本距离肿瘤组织2 cm。所有操作规程均通过扬州大学附属医院伦理委员会审查（2020-YKL03-G042）。所有标本均存储在－80 ℃冰箱保存。使用TRIzol 试剂（美国Invitrogen 公司）从样品中提取总RNA，然后用RNeasyMaxi 试剂盒（德国Qiagen 公司）纯化获得RNA。完成上述步骤后，按照miScriptⅡRT试剂盒（德国Qiagen 公司）的说明书，逆转录合成cDNA。使用miScriptSYBR®GreenPCR 试剂盒（德国Qiagen公司）进行实验，反应所需仪器为LC480荧光定量PCR 仪。反应条件：95 ℃ 30 s，95 ℃10 s，55 ℃ 30 s，72 ℃ 30 s，共45 个循环。相关基因的引物序列见表1。然后，基于筛选并经验证的膀胱癌预后相关基因构建风险评分模型：风险评分=（Coef1×mRNA1的表达量）+（Coef2×mRNA2的表达量）+…+（Coefn×mRNAn 的表达量）。其中，Coef 为相应mRNA 的LASSO Cox 回归模型系数。通过中位风险评分将膀胱癌患者分为高、低风险组。使用R 中的“Kaplan-Meier”包评价两组患者的预后，使用R 中的“survivalROC”包生成ROC 曲线，并使用“timeROC”包计算模型的AUC值。同时，从GEO 数据库中下载GSE13507 数据集（包含165 例具有完整临床随访数据的膀胱癌患者），对模型进行外部验证。

表1 qRT-PCR引物设计Table 1 qRT-PCR primers

1.7 构建预测膀胱癌患者预后的列线图

首先，结合风险评分和临床病理特征进行单因素和多因素Cox 回归分析，以探索风险评分模型是否为膀胱癌患者预后的独立相关因素。然后，基于风险评分和独立的临床病理特征构建列线图，并使用C 指数和校准曲线验证列线图的预测效能。

1.8 统计学方法

使用R 4.0.3 和GraphPad 9.0 进行统计分析。连续性变量用均数±标准差（±s）描述，采用Student’st检验或Wilcoxon秩和检验进行比较；分类变量用频率（%）描述，采用χ2检验或Fisher 精确检验。绘制Kaplan-Meier 生存曲线和时间依赖性ROC 曲线，评估风险评分模型的预测价值。采用单因素和多因素Cox 回归分析确定膀胱癌患者预后的独立影响因素。P<0.05 为具有统计学意义。

2 结 果

2.1 正常组织与肿瘤组织免疫细胞浸润评分差异及与患者预后的相关性

与正常组织标本比较，膀胱癌组织标本中B 细胞、肥大细胞、中性粒细胞、辅助性T 细胞和肿瘤浸润淋巴细胞的浸润评分明显升高（均P<0.05），提示这五种免疫细胞可能在膀胱癌免疫微环境的改变中发挥重要作用，是膀胱癌的关键免疫细胞。见图1。

图1 癌症基因组图谱膀胱癌队列中正常组织与肿瘤组织16种免疫细胞浸润丰度比较Figure 1 Comparison of 16 kinds of immune cell infiltration abundance between normal tissues and tumor tissues in bladder cancer cohort from the Cancer Genome Atlas

聚类分析结果显示，411 份膀胱癌患者标本可初步分为三类（图2）。生存曲线分析发现，聚类2 中患者的预后显著优于聚类1 和聚类3（均P<0.05），而聚类1和聚类3患者的预后无显著差异（P>0.05），见图3A。基于此，本研究将聚类1 和聚类3 合并为一类（聚类1´）重新进行生存分析。结果显示，聚类1´和聚类2 中患者的存活率差异具有统计学意义（图3B）。进一步行主成分分析发现，当TCGA 样本被划分为两个簇时，样本可以明显分离（附图1），且聚类2 患者的免疫细胞浸润评分明显高于聚类1´患者（附图2）。

2.2 不同聚类患者治疗情况比较

比较不同聚类之间免疫检查点相关基因的表达发现，在聚类2 中免疫检查点相关基因的表达均显著升高（附图3A）。TCIA 数据库分析发现，聚类2 中患者对CTLA4 抑制剂、PD-1 抑制剂及两者联合更敏感，更能够从免疫治疗中获益（附图3B）。分析两种分型患者对常用化疗药物的敏感性发现，与聚类1´比较，聚类2患者对恩贝酸、多西他赛、环巴胺和阿卡地新更为敏感（附图3C）。

2.3 关键免疫细胞相关基因组筛选结果

从附图4可以看出，当软域值为5时，数据更符合幂律分布，平均连通性趋于稳定，数据适合进一步研究。确定了软阈值后，根据基因间表达量的相关性构建基因聚类树。将最小模块数设置为10，深度灵敏度设置为 2，使用动态剪切法识别模块。最后，将相似度小于0.8的模块合并，共得到25 个非灰色模块（图4）。在正相关性模块中，棕色模块中除肥大细胞外，其余相关系数均大于0.7，可见其与免疫细胞的相关性系数相对较大；在负相关性模块中，仅有黑灰色模块，见表2。本研究后续选择棕色模块和黑灰色模块，汇总两个模块共得到35 个基因。相关基因的表达见附图5。

图4 基于最优软阈值构建基因聚类树及共表达网络Figure 4 Construction of gene clustering tree and coexpression network based on optimal soft threshold

图5 风险评分模型在TCGA数据集中的验证Figure 5 Validation of the risk scoring model in the TCGA dataset

表2 不同模块与关键免疫细胞的相关性分析Table 2 Correlation analysis between different modules and key immune cells

2.4 关键基因富集分析结果

对35 个关键基因进行GO 分析结果显示，在生物过程及分子功能分析中，35 个关键基因显著参与单核细胞和淋巴细胞等免疫细胞的分化以及正向调节细胞因子的产生（附图6）；在细胞成分分析中，相关基因主要富集于质膜的外侧。KEGG 富集分析结果显示，这些基因主要参与T细胞受体信号通路、细胞黏附分子及细胞因子-受体相互作用（附图7）。

图6 风险评分模型在GSE13507数据集中的验证Figure 6 Validation of the risk scoring model in GSE13507 dataset

图7 基于风险评分构建膀胱癌患者预后相关列线图Figure 7 Prognostic nomogram of bladder cancer patients based on risk score

2.5 膀胱癌预后风险评分模型的构建与验证

基于之前筛选出来的35个基因，进一步通过LASSO Cox 回归筛选出4 个与患者预后相关的基因GPR171、HOXB3、HOXB5和HOXB6。通过HPA数据库比较相关基因的蛋白表达水平发现，GPR171 在肿瘤组织中高表达，HOXB3 在肿瘤组织中低表达（附图8）；通过临床组织标本进行验证发现，HOXB5和HOXB6mRNA 在肿瘤组织中均高表达（分别是正常组织的3.75 和2.87 倍，均P<0.01）。基于GPR171、HOXB3、HOXB5和HOXB6构建膀胱癌患者的预后风险评分模型：风险评分=（－0.3654×GPR171 表达水平）+（－0.2885×HOXB3表达水平）+（－0.2885×HOXB5 表达水平）+（0.3087×HOXB6表达水平）。采用上述公式计算样本的风险评分范围为－2.41~－0.03，评分中位数为－0.83，根据中位风险评分将膀胱癌患者分为高风险组和低风险组，四个基因在高低风险组的表达见附图9。

图8 列线图模型用于预测患者1、3、5年总存活率的校准图Figure 8 Calibration plots of nomogram models for predicting 1， 3， or 5 year overall survival rates

Kaplan-Meier 生存曲线提示，高风险组患者预后较差，风险评分与预后呈负相关（图5A）。ROC 曲线分析结果显示，风险评分预测患者1、3、5年总生存时间的曲线下面积分别为0.735、0.765、0.799，均超过0.7（图5B）。在外部数据集GSE13507 中验证评分模型的预测价值，结果与前述一致，低风险组预后好（图6A）；随时间变化的ROC 曲线显示，1年的AUC 为0.667，3年的AUC 为0.613，5年时AUC 为0.629（图6B）。提示基于关键免疫细胞相关基因构建的预后风险评分模型对于膀胱癌患者预后预测的准确性较好。

2.6 膀胱癌预后评估列线图的构建和验证

单因素Cox 回归分析结果显示，年龄、T 分期、病理分期和风险评分与患者的预后显著相关；多因素分析显示，风险评分、T 分期和病理分期是膀胱癌患者预后的独立危险因素。见表3。结合风险评分和临床参数构建列线图，C 指数为0.72（图7）；校准曲线显示，模型的预测结果与实际结果基本重合（图8），提示列线图对膀胱癌患者1、3、5年总存活率具有较好的预测效能。

表3 膀胱癌患者预后相关危险因素分析结果Table 3 Analysis of prognostic risk factors in patients with bladder cancer

3 讨 论

膀胱癌是世界范围内发病率较高且最致命的恶性肿瘤之一，肿瘤的高复发风险使得膀胱癌患者的预后较差，给社会和患者造成了沉重的经济负担［2］。目前，手术结合膀胱内灌注治疗是非肌层浸润性膀胱癌的主要治疗手段［12］。随着对膀胱癌的免疫病理机制和肿瘤免疫微环境的进一步了解，基于肿瘤免疫微环境调控的免疫治疗已成为膀胱癌患者新的治疗选择［13］。肿瘤免疫微环境由许多具有抗肿瘤或促肿瘤活性的不同免疫细胞亚群组成。既往研究发现，巨噬细胞作为先天免疫效应细胞，在免疫治疗中起着至关重要的作用，与免疫检查点阻断应答有关［14-15］。Jiang等［16］通过巨噬细胞相关基因构建并验证了预测膀胱癌患者预后的模型。除了巨噬细胞之外，CD8+T细胞、NK细胞和B细胞等多种免疫细胞参与机体免疫。这些细胞通过不同途径促进或抑制机体免疫微环境的形成，各个细胞间存在协同或抑制效应，从而在肿瘤免疫中发挥作用。因此，系统研究免疫细胞在膀胱癌的作用具有重要意义。

随着免疫检查点抑制剂在临床中的运用，肌层浸润性膀胱癌患者能够获得更好的预后。PD-1（PDCD1）与肿瘤组织产生的PD-L1（CD274）结合，导致有限的宿主免疫应答。PD-L1 抑制剂通过增加浸润的CD8+T 细胞水平，显示出有效的抗肿瘤免疫应答［17］。T 淋巴细胞在癌症的“免疫监视”中发挥着重要作用，不仅通过表达多态性抗原受体以识别特异性抗原，还拥有效应功能以及记忆特征［18-19］。因此，本研究以16 种免疫细胞为基础，通过ssGSEA 算法量化了TGCA 膀胱癌队列中免疫细胞的丰度，系统地探讨了免疫细胞评分与膀胱癌患者治疗及预后的关联。通过无监督聚类将膀胱癌基于免疫细胞评分分为两个亚型，其中聚类2 患者的预后明显优于聚类1´患者。有报道称，PD-L1 的高表达水平与膀胱癌肿瘤组织的恶性程度及不良预后显著相关，且此类患者术后复发率较高［15］。本文资料显示，PDCD1、CD28和CTLA4等重要免疫检查点相关基因在聚类2 显著上调，而免疫检查点相关基因的表达与肿瘤免疫微环境显著相关，低表达的免疫检查点相关基因可能促进肿瘤免疫逃逸［11］。阿特珠单抗是目前唯一批准用于膀胱癌治疗的PD-L1 抑制剂，但只有少数患者从免疫治疗中受益［20］。本研究开发的免疫评分聚类分型方法根据膀胱癌患者对PD-1 或CTLA4 抑制剂的敏感性进行分层发现，聚类2 对PD-1 抑制剂和CTLA4 抑制剂联合治疗反应更强，表明聚类2 是免疫有利亚群，聚类2 聚类中的患者更可能从阿特珠单抗治疗中获益。此外，本研究还发现聚类2 对于部分化疗药物的半抑制浓度小于聚类1´，提示聚类2 的膀胱癌患者更容易从这些化疗药物中获得临床疗效，避免更高剂量的化疗药物所致各种不良反应。

为了更好地将这些相关免疫浸润评分应用于临床，本研究通过WGCNA 筛选出与关键免疫细胞显著相关的基因，GO 富集分析显示这些基因主要参与免疫细胞的分化与细胞因子产生的正向调节。KEGG 富集分析表明这些基因主要参与T 细胞受体信号通路。通过LASSO Cox 回归进一步筛选出四个与膀胱癌患者预后相关的基因。其中，GPR171 是一种G 蛋白偶联受体，由神经内分泌途径蛋白衍生的BigLEN 是GPR171 的内源性配体。有研究显示，GPR171 可调节摄食和焦虑行为［21］，并在多种肿瘤中被证实是T 淋巴细胞的特征基因［22-23］。研究发现，GPR171在肿瘤浸润淋巴细胞中上调，且其在黑色素瘤患者中的表达可以预测患者接受免疫治疗的反应［24］。GPR171 在T淋巴细胞中的表达动力学与其他免疫检查点（如PD1、TIGIT、CD112R 和TIM36）非常相似，在T 淋巴细胞中GPR171 可上调并通过信号转导抑制T淋巴细胞介导的免疫反应。此外，GPR171 信号的破坏可促进T 淋巴细胞介导的抗肿瘤免疫［25］。HOXB3、HOXB5 和HOXB6 都属于HOXB 家族。HOXB基因家族是一个长度约为180 bp的同源片段，编码60 个氨基酸的同源结构域，与转录调控相关。HOXB 基因家族成员在肿瘤细胞增殖、迁移和凋亡以及肿瘤的形成和发展中均起着至关重要的作用［26］。已有研究表明，多个HOXB 基因家族成员在膀胱癌、胆管癌、子宫内膜癌和乳腺癌种差异表达［27］，如ISL LIM 同源盒蛋白1 和LIM 同源盒蛋白5 在膀胱肿瘤中表达上调，并与肿瘤的发展及转移密切相关［28］。HOXB3是一种蛋白编码基因，与急性髓系白血病的一个生物学亚群显著相关。乳腺癌和肺腺癌患者HOXB3 的表达与患者的预后相关。乳腺癌患者的HOXB3表达水平较低，往往预后较差［29］，而在肺腺癌中高表达HOXB3 的患者预后不良。研究发现，HOXB3 与恶性肿瘤的增殖、侵袭和转移显著相关［30］。在膀胱癌中，HOXB3与吉西他滨的耐药性相关，并可作为膀胱癌患者预后的生物标志物［29］。HOXB5 能够与特定的DNA 侧结合，从而影响器官发育和肿瘤发生中的关键信号通路［31］。据报道，HOXB5在肿瘤中表达失调，如HOXB5在视网膜母细胞瘤细胞表达上调，并通过增加基质金属蛋白酶的产生促进肿瘤细胞的迁移和侵袭［32］。此外，HOXB5还具有通过激活Wnt/β-连环蛋白途径促进非小细胞肺癌恶性进展的能力［31］。在乳腺癌中，Zhang 等［33］报道HOXB5 的过表达增强了乳腺癌细胞的侵袭能力。研究证实，HOXB5在膀胱癌组织和细胞系中过度表达，并促进膀胱癌细胞的增殖和迁移［34］。HOXB6 在结直肠肿瘤中差异表达并可作为患者预后的生物标志物［35］。研究表明，敲低HOXB6的表达能够抑制肝细胞癌细胞的迁移、侵袭和增殖［36］。此外，Huang等［37］报道了Hsa_circ_0007031 通过海绵miR-196a-5p 调节HOXB6 来促进骨肉瘤细胞的增殖和迁移。目前关于HOXB6 与膀胱癌的研究还甚少，需要进一步研究其机制。

本研究基于上述关键基因构建了风险评分模型，并通过外部验证证实了该风险评分模型的预测效能。结合临床特征构建的列线图也展示出对膀胱癌患者预后较好的预测价值。与其他基因集构建的膀胱癌预后风险评估模型比较，本研究构建的模型在预测膀胱癌患者5年存活率上有着更高的效能［38-39］。此外，相比Liang等［40］构建的15个基因的模型，本研究构建的模型纳入的基因数较少，检验难度小，便于临床应用。

综上所述，本文系统研究了膀胱癌患者免疫细胞浸润状态，并基于关键免疫细胞相关基因构建了预后风险评分模型，对膀胱癌患者的个体化治疗具有重要意义。但是，由于样本量有限，本研究中构建的模型还需通过前瞻性、多中心、更大样本的研究进一步验证，并通过开展相关功能体内外实验探索关键免疫细胞和核心基因在膀胱癌中发挥作用的相关机制。

本文附图见电子版。

志谢研究得到国家自然科学基金（82002675）、江苏省科技计划青年基金（BK2020938）、扬州市重点研发计划社会发展项目（YZ2020110）、扬州市软科学研究计划（YZ2022267）、江苏省博士后研究资助计划（2020Z268）、扬州大学高层次人才研究启动基金（2019LYF）支持

AcknowledgementsThis work was supported by the National Natural Science Foundation of China (82002675),Science and Technology Plan Youth Fund Project of Jiangsu Province (BK2020938), Key R&D Plan-Social Development Project of Yangzhou City (YZ2020110), Soft Science Research Plan of Yangzhou City (YZ2022267), Postdoctoral Research Funding Program of Jiangsu Province (2020Z268),and Yangzhou University High-level Talent Research Startup Fund (2019LYF)

利益冲突所有作者均声明不存在利益冲突

Conflict of InterestsThe authors declare that there is no conflict of interests

©The author(s) 2024. This is an open access article under the CC BY-NC-ND 4.0 License (https://creativecommons.org/licenses/by-nc-nd/4.0/)