1. 引言
膀胱癌(Bladder Cancer, BLCA)是常见的泌尿系统恶性肿瘤之一，具有较高的发病率和复发率。其中，非肌层浸润性膀胱癌约占全部病例的75%以上。尽管手术、化疗及免疫治疗不断发展，患者预后仍存在显著差异，其分子机制尚未完全阐明[1]。
肿瘤代谢重编程被认为是肿瘤发生发展的重要标志之一。嘧啶代谢(Pyrimidine Metabolism)作为核苷酸合成和降解的关键通路，在肿瘤细胞增殖、DNA修复及免疫调控中发挥重要作用。已有研究表明，嘧啶代谢相关基因在乳腺癌中可成功区分免疫“冷热肿瘤”，但其在膀胱癌中的作用仍缺乏系统研究。因此，本研究旨在基于嘧啶代谢相关基因对膀胱癌进行分型，并深入探讨其免疫微环境及临床预后意义[2]。
2. 材料与方法
2.1. 数据来源与预处理
本研究所用膀胱癌转录组数据来源于The Cancer Genome Atlas (TCGA)数据库中的膀胱尿路上皮癌项目(TCGA-BLCA)。共纳入416例膀胱癌患者的RNA-Seq表达数据及对应的临床信息。基因表达数据采用标准化处理后的转录本表达矩阵，用于后续分析。
嘧啶代谢相关基因集来源于Molecular Signatures Database (MSigDB)，选取与嘧啶代谢通路相关的基因作为研究对象。对基因表达矩阵进行整理，保留在TCGA-BLCA数据集中具有完整表达信息的嘧啶代谢相关基因用于后续分型分析[3]。
2.2. 嘧啶代谢相关分型的构建
基于嘧啶代谢相关基因的表达特征，采用ConsensusClusterPlus R包对膀胱癌样本进行一致性聚类分析。聚类过程中采用欧氏距离和PAM聚类算法，通过重复抽样1000次以评估聚类结果的稳定性。
根据一致性聚类生成的累积分布函数(CDF)曲线及deltaarea曲线，当聚类数k = 2时，CDF曲线斜率变化最为显著，且组内一致性较高，提示该分组方式具有较好的稳定性和生物学合理性。因此，最终将膀胱癌样本分为Cluster1和Cluster2两个嘧啶代谢相关亚型[4]。
2.3. 差异表达分析与功能富集分析
基于R语言limma包对两种嘧啶代谢亚型进行差异表达分析。差异基因筛选标准设定为|log2Fold Change| > 1且校正后P值(padj) < 0.05。筛选得到的差异基因用于后续功能注释分析。
对差异基因进行Gene Ontology (GO)功能富集分析，以探究其在生物过程(Biological Process)、细胞组成(Cellular Component)及分子功能(Molecular Function)层面的潜在生物学意义。富集结果以校正后P值 < 0.05作为显著性阈值[5]。
2.4. 肿瘤免疫微环境分析
为系统评估不同嘧啶代谢亚型的免疫细胞浸润特征，采用CIBERSORT算法对TCGA-BLCA样本中22种免疫细胞的相对丰度进行估计。仅保留CIBERSORT分析结果中P值 < 0.05的样本用于后续比较分析。
同时，采用xCell算法进一步对多种免疫及基质细胞类型进行评分，以提高免疫浸润评估的准确性和可靠性。通过多算法交叉验证，综合分析不同亚型在免疫细胞组成及肿瘤微环境状态方面的差异[6]。
2.5. 免疫检查点相关分析
对不同嘧啶代谢亚型中常见免疫检查点分子(如PD-1、PD-L1等)的表达水平进行比较分析。通过比较两组间免疫检查点基因的表达差异，评估不同亚型患者接受免疫检查点抑制剂治疗的潜在获益可能[7]。
2.6. 加权基因共表达网络分析(WGCNA)
采用WGCNA R包构建基因共表达网络，以识别与膀胱癌免疫微环境密切相关的基因模块。首先对基因表达矩阵进行样本聚类，剔除异常样本；随后选择合适的软阈值以满足无标度网络特性。
构建加权共表达网络后，通过动态剪切算法识别不同基因模块，并计算模块特征基因(Module Eigengene)。将各模块与ESTIMATE免疫评分进行相关性分析，筛选与免疫微环境显著相关的关键模块。进一步对目标模块中的基因进行蛋白互作网络(PPI)分析，以筛选潜在核心调控基因[8]。
2.7. 突变特征分析
基于TCGA提供的突变注释文件(Mutation Annotation Format, MAF)，对不同嘧啶代谢亚型的体细胞突变特征进行系统分析。比较两组在关键突变基因、突变类型及总突变负荷(Tumor Mutation Burden, TMB)方面的差异，并通过可视化方式展示突变谱特征[9]。
2.8. 药物敏感性预测分析
采用pRRophetic R包基于转录组表达数据预测膀胱癌患者对常用化疗药物及靶向药物的半数抑制浓度(IC50)。比较不同嘧啶代谢亚型在多种药物敏感性方面的差异，以评估其潜在治疗反应特征，为个体化治疗策略提供依据。
2.9. 预后模型构建与评估
基于两种嘧啶代谢亚型的差异表达基因，采用Lasso-Cox回归分析筛选与总体生存(Overall Survival，OS)显著相关的关键基因，并构建预后风险评分模型。Lasso回归通过引入惩罚参数以降低多重共线性的影响，仅保留回归系数不为0的基因用于模型构建。风险评分(Risk score)的计算公式如下：Risk score = β0 + Σ (βi × Expi)。其中，β0表示模型截距，βi表示第i个基因对应的回归系数，Expᵢ表示该基因在样本中的标准化表达水平。纳入模型的基因及其对应回归系数详见附录表S1。
根据风险评分的中位数将患者划分为高风险组和低风险组。通过Kaplan-Meier生存分析比较不同风险组的生存差异，并采用时间依赖性ROC曲线评估模型在1、2和3年生存预测中的准确性。同时，将风险评分与临床病理特征(包括年龄、性别及TNM分期)共同纳入单因素及多因素Cox回归分析，以评估风险评分的独立预后价值。此外，通过决策曲线分析(Decision Curve Analysis, DCA)评估风险模型在临床决策中的潜在应用价值[10]。
2.10. 统计学分析
所有统计分析均在R软件环境(R version 4.x)中完成。连续变量采用Wilcoxon秩和检验或t检验进行比较，分类变量采用卡方检验。双侧P值 < 0.05认为差异具有统计学意义[10]。
2.11. Human Protein Atlas数据库验证
为进一步验证预后模型中关键基因的表达特征，选取风险评分模型中回归系数绝对值较高的代表性基因，利用Human Protein Atlas (HPA)数据库对其在正常膀胱组织与膀胱癌组织中的蛋白表达情况进行免疫组化分析。通过对比不同组织中免疫组化染色强度及分布特征，评估关键基因在蛋白水平上的表达差异。
3. 结果
3.1. 嘧啶代谢相关分子分型的构建
基于TCGA-BLCA队列中嘧啶代谢相关基因的表达谱，采用ConsensusClusterPlus算法对膀胱癌样本进行一致性聚类分析。通过比较不同聚类数下的累积分布函数(CDF)曲线及delta area曲线变化情况，当聚类数k = 2时，曲线斜率变化最为显著，且组内一致性较高，提示该分型方式具有较好的稳定性。因此，最终将膀胱癌患者分为2个嘧啶代谢相关分子亚型，分别命名为Cluster1 (197例)和Cluster2 (206例)。
3.2. 不同嘧啶代谢亚型的差异表达特征及功能富集分析
对Cluster1和Cluster2两个亚型进行差异表达分析，共筛选获得2651个差异表达基因，其中上调基因1793个，下调基因858个。差异基因的筛选标准为|log2Fold Change| > 1且校正后P值(padj) < 0.05。
GO功能富集分析结果显示，这些差异基因在生物过程层面主要富集于T细胞活化、白细胞黏附、免疫应答调控等免疫相关通路，提示不同嘧啶代谢亚型在免疫功能状态方面存在显著差异(图1)。
Figure 1. Functional enrichment analysis of differentially expressed genes among pyrimidine metabolism-related subtypes
图1. 嘧啶代谢相关亚型间差异表达基因的功能富集分析
3.3. 不同嘧啶代谢亚型的免疫细胞浸润特征
采用CIBERSORT算法评估两种嘧啶代谢亚型的免疫细胞浸润情况。结果显示，Cluster2中CD8+ T细胞、滤泡辅助性T细胞及活化树突状细胞的浸润水平显著高于Cluster1；而Cluster1中M2型巨噬细胞的浸润水平显著升高。
进一步采用xCell算法对免疫及基质细胞进行分析，结果显示Cluster2中CD8+ T细胞、NK细胞及M1型巨噬细胞评分显著升高；而Cluster1中癌相关成纤维细胞(CAF)及内皮细胞评分明显升高。综合CIBERSORT和xCell分析结果，Cluster2表现为免疫细胞高浸润特征，而Cluster1以基质细胞富集为主(图2)。
Figure 2. Immune cell infiltration characteristics of pyrimidine metabolism-related subtypes
图2. 嘧啶代谢相关亚型的免疫细胞浸润特征
3.4. 不同嘧啶代谢亚型的免疫检查点表达特征
对两种嘧啶代谢亚型中免疫检查点相关基因的表达水平进行比较分析。结果显示，Cluster2中PD-1和PD-L1的表达水平均显著高于Cluster1，提示Cluster2肿瘤微环境处于相对免疫激活状态，可能具有更高的免疫治疗潜在获益(图3)。
Figure 3. Expression differences of immune checkpoint molecules in distinct pyrimidine metabolism-related subtypes
图3. 不同嘧啶代谢相关亚型中免疫检查点分子的表达差异
3.5. WGCNA分析及关键免疫相关模块的筛选
采用加权基因共表达网络分析(WGCNA)构建膀胱癌基因共表达网络，并与ESTIMATE免疫评分进行相关性分析。结果显示，Turquoise模块与免疫评分的相关性最强，被认为是与膀胱癌免疫微环境最密切相关的模块。
进一步对Turquoise模块中的基因进行蛋白互作网络(PPI)分析，筛选出多个可能参与免疫调控的关键基因。其中，CXCL10在网络中处于核心位置，其高表达可能与免疫治疗敏感性相关(图4)。
Figure 4. Weighted gene co-expression network analysis (WGCNA) associated with pyrimidine metabolism-related subtypes
图4. 与嘧啶代谢相关亚型相关的加权基因共表达网络分析(WGCNA)
3.6. 不同嘧啶代谢亚型的突变特征分析
基于TCGA提供的突变数据，对不同嘧啶代谢亚型的体细胞突变特征进行分析。结果显示，Cluster2中TP53、TTN和KDM6A的突变频率显著高于Cluster1，而Cluster1中FGFR3和ARID1A的突变更为常见。
在突变类型方面，两组均以错义突变为主要突变类型，但Cluster1中剪接位点突变比例相对较高。进一步分析发现，Cluster2的总突变数(包括SNV和INDEL)显著高于Cluster1，提示其具有较高的肿瘤突变负荷(TMB) (图5)。
Figure 5. Somatic mutation characteristics of pyrimidine metabolism-related molecular subtypes
图5. 嘧啶代谢相关分子亚型的体细胞突变特征
3.7. 不同嘧啶代谢亚型的药物敏感性差异
基于pRRophetic算法预测不同嘧啶代谢亚型对常用化疗药物及靶向药物的敏感性。结果显示，Cluster2对顺铂(Cisplatin)和吉西他滨(Gemcitabine)的预测IC50值显著低于Cluster1，提示其对上述化疗药物更为敏感。
在靶向药物方面，Cluster1对FGFR抑制剂(如Erdafitinib)的预测敏感性显著高于Cluster2，而Cluster2可能对PI3K抑制剂表现出更高敏感性。上述结果提示不同嘧啶代谢亚型在治疗反应方面存在明显差异。
3.8. 预后风险模型的构建与评估
基于两种嘧啶代谢亚型的差异表达基因，采用Lasso回归分析筛选与总体生存相关的关键基因，并构建预后风险评分模型。根据风险评分中位数将患者分为高风险组和低风险组。
Kaplan-Meier生存分析结果显示，高风险组患者的总体生存率显著低于低风险组(P < 0.0001)。时间依赖性ROC曲线分析显示，该模型在1、2和3年生存预测中的AUC分别为0.792、0.818和0.774，提示模型具有较好的预测性能。
进一步的单因素和多因素Cox回归分析显示，风险评分在校正年龄、性别及肿瘤分期等临床因素后，仍为影响患者生存的独立预后因素。基于风险评分构建的列线图可用于预测患者1、2和3年生存概率，校准曲线显示预测结果与实际观察结果具有较好一致性(图6)。
Figure 6. Kaplan-Meier survival curves and time-dependent ROC curves
图6. 高Kaplan-Meier 生存曲线以及时间依赖性ROC曲线
3.9. 关键基因的蛋白表达验证
HPA数据库免疫组化结果显示，预后模型中关键基因在膀胱癌组织中普遍呈现中到强阳性染色，而在正常膀胱组织中染色较弱或呈阴性表达。上述结果在蛋白水平上支持了转录组分析中所观察到的表达差异，进一步验证了所构建预后模型的生物学合理性。
4. 讨论
近年来，代谢重编程被认为是肿瘤的核心特征之一，不仅参与肿瘤细胞的能量供给和生物合成过程，还深刻影响肿瘤免疫微环境及治疗响应。本研究基于嘧啶代谢相关基因对膀胱癌进行系统分型，发现不同嘧啶代谢模式与免疫微环境状态、基因突变特征、药物敏感性及患者预后密切相关，提示嘧啶代谢在膀胱癌发生发展及免疫调控中具有重要作用。
4.1. 嘧啶代谢与膀胱癌免疫微环境的潜在联系
嘧啶代谢是核苷酸代谢的重要组成部分，直接参与DNA和RNA的合成、修复及细胞周期调控。肿瘤细胞通常通过增强从头合成途径以满足其快速增殖的需求，而这一过程往往伴随着代谢通路的异常激活。既往研究表明，嘧啶代谢紊乱可通过影响DNA损伤修复能力，增加基因突变负荷，从而间接调控肿瘤免疫原性。
本研究发现，基于嘧啶代谢相关基因的表达特征，可将膀胱癌样本稳定分为两种分子亚型。Cluster2表现为显著的免疫细胞浸润增强、免疫激活通路富集以及免疫检查点高表达，呈现典型“热肿瘤”特征；而Cluster1则以基质细胞富集、免疫抑制状态为主，符合“冷肿瘤”表型。这一结果提示，嘧啶代谢状态可能通过影响肿瘤细胞增殖压力、突变积累及免疫原性，进而塑造不同的肿瘤免疫微环境[11]。
4.2. 不同嘧啶代谢亚型的免疫细胞浸润特征及生物学意义
通过CIBERSORT与xCell等多种算法交叉验证免疫细胞浸润情况，本研究发现Cluster2中CD8+ T细胞、NK细胞、活化树突状细胞及M1型巨噬细胞显著升高，这些细胞被认为是抗肿瘤免疫反应的重要效应细胞。其中，CD8+ T细胞和NK细胞可直接介导肿瘤细胞杀伤，而活化树突状细胞在抗原呈递及T细胞启动中发挥关键作用。
相较之下，Cluster1中M2型巨噬细胞、癌相关成纤维细胞(CAF)及内皮细胞显著富集。M2型巨噬细胞和CAF被广泛认为可通过分泌免疫抑制因子、促进血管生成和基质重塑，抑制抗肿瘤免疫反应，从而促进肿瘤进展。这种免疫抑制型微环境可能是Cluster1患者预后较差的重要原因之一[12]。
4.3. 嘧啶代谢分型与免疫治疗潜在响应的关系
免疫检查点抑制剂(ICB)已成为膀胱癌的重要治疗手段，但临床疗效存在显著个体差异。本研究发现，Cluster2中PD-1/PD-L1表达水平显著升高，提示该亚型存在活跃的T细胞–肿瘤细胞相互作用，肿瘤微环境处于免疫激活状态。这一特征与其高免疫浸润及高突变负荷相一致，表明Cluster2患者可能更易从PD-1/PD-L1抑制剂治疗中获益[13]。
值得注意的是，尽管Cluster2具有较高的TMB，但单纯依赖TMB指标仍可能不足以准确预测免疫治疗疗效。本研究提示，将嘧啶代谢分型、TMB及免疫检查点表达水平相结合，可能有助于更精准地筛选免疫治疗获益人群[14]。
4.4. 突变谱特征及其对治疗策略的启示
在突变特征方面，Cluster2中TP53、TTN和KDM6A突变频率显著升高，整体突变负荷明显高于Cluster1。TP53突变通常与基因组不稳定性增加及免疫原性增强相关，这可能进一步解释Cluster2中高免疫浸润和免疫激活状态。
而Cluster1中FGFR3和ARID1A突变更为常见。既往研究表明，FGFR3突变在非肌层浸润性膀胱癌中较为高发，且常与免疫“冷”微环境相关。本研究中Cluster1对FGFR抑制剂表现出更高敏感性，提示该亚型患者可能更适合靶向治疗或联合治疗策略，以改善免疫抑制状态[15]。
4.5. 基于嘧啶代谢分型的药物敏感性差异
通过药物敏感性预测分析，本研究发现不同嘧啶代谢亚型对多种化疗及靶向药物存在显著差异。Cluster2对顺铂、吉西他滨等常用化疗药物敏感性更高，可能与其较高的DNA损伤及修复缺陷有关；而Cluster1对FGFR抑制剂更为敏感，与其FGFR3突变特征高度一致[16]。
上述结果提示，嘧啶代谢分型不仅具有预后评估价值，还可为膀胱癌的个体化用药提供重要参考。
4.6. 预后模型的临床价值与局限性
本研究基于差异基因构建的Lasso风险评分模型在生存预测中表现出良好的判别能力，且在多因素分析中被证实为独立预后因素。通过列线图整合临床特征与风险评分，可直观评估患者不同时间点的生存概率，具有一定的临床应用潜力。
然而，本研究仍存在一定局限性。首先，分析主要基于TCGA数据集，缺乏外部独立队列验证；其次，尚未结合单细胞或空间转录组数据对肿瘤微环境进行更精细解析；尽管本研究主要基于公共数据库进行分析，尚未进一步开展体内外功能实验，但通过补充多变量Cox回归分析、决策曲线分析以及Human Protein Atlas数据库的蛋白表达验证，在一定程度上增强了预后模型的稳定性、可重复性及生物学可信度。未来仍需结合独立外部队列及功能实验对相关机制进行深入验证。此外，相关机制尚需进一步通过体内外实验加以验证[17]。
5. 结论
本研究基于嘧啶代谢相关基因对膀胱癌进行了分子分型，成功识别出具有显著免疫微环境差异的两种亚型。不同嘧啶代谢亚型在免疫细胞浸润特征、免疫检查点表达、基因突变谱及药物敏感性方面存在明显差异，并与患者预后密切相关。基于差异基因构建的预后风险模型具有较好的预测性能，风险评分可作为膀胱癌患者的独立预后因素。上述结果提示，嘧啶代谢相关分型有助于揭示膀胱癌的免疫异质性，并为免疫治疗及个体化治疗策略的制定提供潜在依据。
NOTES
*通讯作者。