膀胱癌m6A调节因子预后模型建立与分析

白洋洋，郭依琳，陈瑞廷，潘世杰，孙继建

1. 河南省中医院（河南中医药大学第二附属医院）泌尿外科（郑州 450002） 2. 郑州大学第二附属医院妇科肿瘤临床医学研究中心（郑州 450014）

膀胱癌（bladder cancer, BC）是常见的男性泌尿系统恶性肿瘤之一。2020年全球BC 新发病例数超过55 万，死亡病例数超过20 万[1]。在中国，BC 发病率位居全部恶性肿瘤的第13 位，严重威胁男性的生命健康[2-3]。根据是否侵犯膀胱肌层，BC 可分为非肌层浸润性和浸润性。非肌层浸润性BC 患者首选经尿道膀胱肿瘤切除术联合膀胱内治疗，术后生活质量较好，但约30%患者存在局部浸润或远处转移可能[4]。对于浸润性BC，约25%患者伴有淋巴结转移，5%出现远处转移，预后较差[5]。近年来，随着免疫治疗的应用，BC患者预后有所改善，但由于肿瘤异质性，并非所有BC 患者均能从中获益[6-7]。因此，探究BC 发生发展的分子机制对提高BC 患者治疗反应率意义重大。

N6-甲基腺嘌呤（N6-methyladenosine, m6A）是在RNA 腺嘌呤残基的第6 位含氮碱基位置上添加一个甲基基团，是真核生物中最丰富的表观转录组学修饰[8]。研究发现，m6A 修饰主要由m6A调节因子，即编码器-甲基转移酶复合物、消码器-去甲基化酶和读码器-m6A 读取蛋白调控[9]。有研究指出，m6A 调节因子具有免疫调节功能，可影响肿瘤微环境（tumor microenvironment, TME）中免疫细胞浸润，并对免疫治疗疗效产生影响[10]。目前，m6A 调节因子对BC 预后影响的相关研究较少。本研究利用预后相关的m6A 调节因子构建BC 预后预测模型，评估TME 中免疫细胞浸润情况，进一步预测免疫治疗和靶向治疗的疗效，为BC 的临床治疗提供参考。

1 资料与方法

1.1 RNA测序转录组数据和临床数据获取

利用癌症基因组图谱（The Cancer Genome Atlas, TCGA）数据库（https://www.cancer.gov/ccg/research/genome-sequencing/tcga），获取397 例BC组织的高通量转录组测序数据和对应的临床病理特征数据。397 例BC 患者的临床病理特征见表1。

表1 397例BC患者的临床病理特征Table 1. Clinical and pathological characteristics of 397 patients of BC

1.2 预后相关的m6A调节因子筛选

检索中国知网、Pubmed、Embase 数据库，中文检索词包括m6A 调节因子、m6A 修饰，英文检索词包括m6A RNA methylation、m6A regulatory factors。共筛选到26 个m6A 调节因子，包括9 个编码器（WTAP、ZC3H13、RBM15、METTL3、RBM15B、KIAA1429、METTL5、METTL14、METTL16）、14 个读码器（LRPPRC、HNRNPA2B1、IGF2BP1、IGF2BP2、IGF2BP3、YTHDF1、YTHDF2、YTHDF3、HNRNPC、FMR1、ELAVL1、YTHDC1、YTHDC2、RBMX）和3 个消码器（ALKBH5、FTO、ALKBH3）。通过R 软件corrplot 包计算26 个m6A 调节因子之间表达的相关性，利用survival 包对26 个m6A调节因子进行单因素Cox 回归分析，结果通过forestplot 包输出，筛选预后相关的m6A 调节因子。

1.3 预后预测模型构建及验证

利用最小绝对值选择与收缩算子（least absolute selection and shrinkage operator, LASSO）Cox 回归，筛选出5 个预后相关m6A 调控因子（P＜0.05），构建预后预测模型。每例BC 患者的风险评分等于5 个m6A 调控因子中每个基因表达量乘以系数之和（Coefi为风险系数，Xi 为基因表达量）。计算公式如下：风险评分=分别计算397例BC患者的风险评分，根据中位风险评分分为高风险组和低风险组，比较两组间5年生存差异。绘制受试者工作特征（receiver operating characteristic curve, ROC）曲线并计算ROC 曲线下面积（area under ROC curve, AUC），评估所构建预后预测模型的预测价值。利用单因素Cox 回归和多因素Cox 回归评估预后预测模型和临床病理特征对BC 预后的影响。

1.4 基因集富集分析和免疫浸润分析

在上述构建的BC 预后预测模型中，通过MSigDB 数据库（http://www.broadinstitute.org/gsea/msigdb/index.jsp）获取京都基因与基因组百科全书（Kyoto Encyclopedia of Genes and Genomes,KEGG）相关通路基因集，利用R 软件GSVA 包对高风险组和低风险组进行富集分析，明确两组KEGG 相关信号通路。利用聚类分析和单样本基因集富集分析（single sample gene set enrichment analysis, ssGSEA）[12]评估高风险组和低风险组中23 种免疫细胞类型的频率差异，并计算免疫细胞的相对浸润丰度。利用估计恶性肿瘤组织中基质和免疫细胞（estimation of stromal and immune cells in malignant tumors using expression data,ESTIMATE）方法评估高风险组和低风险组中每个样本中免疫、基质成分及两者之和在TME 中的比例，即Immune 分数、Stromal 分数和Estinate 分数。

1.5 基因相关性分析

利用R 软件limma 包，比较高风险组和低风险组在免疫检查点相关基因PD-L1和CTLA- 4[13]、靶向治疗相关基因HER-2[14]、BRAF[15]、KRAS[16]和EGFR[17]的表达差异，绘制箱式图，以P＜0.05为差异具有统计学意义。

2 结果

2.1 预后相关m6A调节因子筛选

本研究纳入397 例BC 组织样本。利用R 软件分析26 个m6A 调节因子在BC 中表达水平之间的相关性，结果显示除IGF2BP2 与YTHDC1、IGF2BP2 与METTL3，以及ALKBH3 与FMR1 外，大部分m6A 调节因子之间的表达存在显著正相关性，其中HNRNPC 与HNRNPA2B1 的正相关性最强（r=0.66），见图1-A、图1-B。为进一步筛选与BC 预后相关的m6A 调节因子，参考Feng 等[18]的研究，通过单变量Cox 回归发现YTHDC1、IGF2BP3、LRPPRC、FTO 和ALKBH3 是BC 独立的预后因素，其中IGF2BP3、LRPPRC、FTO 和ALKBH3 为危险调节因子，其高表达与较差的预后有关，而YTHDC1 为保护调节因子，其高表达与较好的预后有关，见图2。

图1 BC中m6A调节因子表达相关性Figure 1. Correlation of m6A regulator expression in BC注：A. 26个m6A调节因子表达的相关性热图；B. 26个m6A调节因子表达的相关性圈图。

图2 BC中m6A调节因子单因素Cox回归分析Figure 2. Univariate Cox regression analysis of m6A regulators in BC

2.2 构建BC预后预测模型

利用上述筛选的YTHDC1、IGF2BP3、LRPPRC、FTO 和ALKBH3，通过LASSO Cox 回归方法建立BC 预后预测模型，见图3-A。在模型中，YTHDC1、IGF2BP3、LRPPRC、FTO 和ALKBH3 的系数分别为-0.08、0.12、0.03、0.19和0.03。根据模型计算每例BC 患者的风险评分，将患者分为高风险组和低风险组。低风险组的OS显著高于高风险组（P＜0.001），见图3-B。利用ROC 曲线分析预后模型预测BC 患者5年生存率的准确性，结果显示AUC 值为0.665。

图3 基于m6A调节因子构建BC风险预后模型Figure 3. Prognostic risk model of BC based on m6A regulators注：A. LASSO Cox回归模型构建预后预测模型；B. 高风险组和低风险组生存分析。

2.3 预后预测模型的临床价值

绘制YTHDC1、IGF2BP3、LRPPRC、FTO 和ALKBH3 在高风险组和低风险组中的表达量及与临床特征的热图，结果显示其表达量在临床分期中的差异有统计学意义，见图4。单因素Cox 回归和多因素Cox 回归分析结果显示，年龄、临床分期和风险评分与OS 相关（P＜0.001），均为BC 独立的预后影响因素，见表2。

图4 基于5个m6A调节因子构建的预后模型与临床特征的关系热图Figure 4. Heatmap of the relationship between prognostic model based on five m6A regulators and clinical characteristics注：***P＜0.001，**P＜0.01。

表2 BC预后预测模型与临床特征的关系Table 2. Relationship between BC prognostic prediction model and clinical characteristics

2.4 基因集富集分析

基因集富集分析结果显示，高风险组在趋化因子、NOD 样受体、嘌呤代谢、丙酮酸代谢等信号通路富集，而低风险组在药物代谢细胞色素P450、亚油酸代谢、α-亚麻酸代谢等信号通路富集，见图5。

2.5 免疫浸润分析

采用ssGSEA 评估高风险组和低风险组中23种免疫细胞类型的频率差异，结果发现高风险组免疫细胞浸润程度显著高于低风险组，特别是激活的CD4+T 细胞、激活的CD8+T 细胞、激活的B细胞、嗜酸性粒细胞、激活的树突细胞、Th1 细胞、Th17 细胞、中性粒细胞等，提示高风险组和低风险组具有不同的免疫细胞浸润特征，见图6-A。ESTIMATE 计算高风险组和低风险组中每个样本中免疫、基质成分及两者之和在TME 中的比例，发现高风险组Immune 分数、Stromal 分数和Estinate 分数均高于低风险组，见图6-B。上述结果表明高风险组的TME 为免疫炎症型，具有丰富的免疫细胞浸润。

图6 BC预后预测模型的免疫浸润分析Figure 6. Immune infiltration analysis of prognostic risk model in BC注：A. ssGSEA分析；B. ESTIMATE分析；***P＜0.001。

2.6 基因相关性分析

进一步分析免疫检查点相关基因PD-L1和CTLA-4、靶向治疗相关基因HER-2、BRAF、KRAS和EGFR在高风险组和低风险组之间的表达差异，结果发现PD-L1、CTLA-4、EGFR和KRAS在高风险组中表达更高，HER-2和BRAF在低风险组表达更高，差异均有统计学意义（P＜0.05），见图7。以上结果提示高风险组可能对PD-L1 抑制剂、CTLA-4 抑制剂、EGFR 抑制剂和KRAS 抑制剂更敏感，而低风险组可能对HER2 抗体和BRAF 抑制剂更敏感。

图7 高低风险组间免疫检查点相关基因和靶向治疗相关基因表达差异Figure 7. Differences in the expression of immune checkpoint related genes and targeted therapy related genes between high and low risk group注：A. PD-L1基因；B. CTLA-4基因；C. KRAS基因；D. EGFR基因；E. HER-2基因；F. BRAF基因。

3 讨论

BC 是一种来源于膀胱黏膜的恶性肿瘤。手术和放化疗是其主要的治疗方式，然而即使接受标准治疗，仍有50%～70%患者出现复发或转移，晚期和转移性BC 的治疗效果不理想[5]。目前BC仍缺乏特异性预后标志物，因此，探究新型预后预测标志物并建立预后预测模型对治疗BC 具有重要的临床价值。

m6A 调节因子的异常表达可能影响肿瘤的发生发展[9]。随着m6A 修饰研究领域的不断突破，越来越多的新m6A 调节因子被发现。本研究分析了26 个m6A 调控因子表达的相关性，发现BC中m6A 调节因子之间的表达存在显著相关性。进一步筛选与BC 预后相关的m6A 调节因子，单变量Cox 回归结果显示YTHDC1、IGF2BP3、LRPPRC、FTO 和ALKBH3 是BC 独立的预后因素。IGF2BP3 和IGF2BP1 同属IGF2BP 家族，是m6A阅读蛋白的一种。有研究证实IGF2BP3 在BC 组织内的高表达，可通过调节HMGB1 信号通路促进BC 的发展，进而影响BC 患者预后[19]。周子健等通过临床标本验证了BC 组织内IGF2BP1 的高表达可提示BC 患者预后不良，转染IGF2BP1 siRNA 后BC 细胞的增殖和生长被显著抑制[20]。脂肪量和肥胖相关蛋白FTO 在多种癌症组织中表达，刺激肿瘤细胞代谢，诱导肿瘤细胞浸润和进展。苏甲林通过检测24 例BC 组织及癌旁正常组织FTO 的表达，发现与对应癌旁组织相比，BC组织中FTO 相对低表达，且与病理分期分级呈负相关[21]。

近年来，有研究利用m6A 调节因子建立头颈鳞癌[22]、肝癌[23]、结直肠癌[24]的预后预测模型，并验证发现其具有较好的预测价值，可以指导后续临床治疗。Wu 等分析了23 个m6A 调节因子与临床病理特征和预后的相关性，提示m6A 调节因子可以作为BC 预后评估的潜在分子靶点[25]。Chen 等分析了13 个m6A 调节因子在BC 中的表达情况，并选择FTO、YTHDC1 和WTAP 构建预后预测模型[26]。Zheng 等利用FTO 和YTHDF2评估BC 预后预测模型的应用价值[27]。本研究纳入了更多m6A 调控因子，更全面地探讨了m6A修饰对BC 预后预测的价值。通过Cox 回归模型选择5 个m6A 调节因子（YTHDC1、IGF2BP3、FTO、ALKBH3和LRPPRC）建立BC预后预测模型，ROC 曲线检测风险评分预测BC 患者5年生存率的准确性，结果显示AUC 值为0.665。本研究与Zheng 等[27]建立的模型（5年AUC=0.614）相比，5年AUC 值更高，进一步分析发现Zheng 等[27]利用单因素Cox 筛选预后相关m6A 调节因子时设置筛选标准为P＜0.15，最后纳入预后模型构建的YTHDF2（P=0.112）对BC 患者预后的影响较小，本研究将筛选标准设置为P＜0.05，以确保参与模型构建的每个m6A 调节因子都与BC 预后相关。另外，按照风险评分将BC 患者分为高风险组和低风险组，经过单因素Cox 回归及多因素Cox 回归发现年龄、临床分期和风险评分与OS 相关，是BC 独立的预后影响因素。

m6A 调节因子具有免疫调节功能，可影响TME 中免疫细胞浸润，并对免疫治疗疗效产生影响[10]。利用ssGSEA 和ESTIMATE 法比较高风险组和低风险组中免疫细胞浸润水平的差异，结果发现，两组具有不同的免疫细胞浸润特征，高风险组具有丰富的免疫细胞浸润，其TME 为免疫炎症型，即热肿瘤对免疫治疗疗效更好[28]。

靶向治疗和免疫治疗的快速发展极大程度提高了BC 患者的OS，并改善了其预后。但由于肿瘤异质性，并非所有的BC 患者均能从靶向治疗和免疫治疗中获益。为了探讨建立的预后预测模型对靶向治疗和免疫治疗的指导作用，本研究选择了免疫检查点相关基因PD-L1和CTLA- 4、靶向治疗相关基因HER-2、BRAF、KRAS和EGFR，比较高风险组和低风险组中这些基因表达的差异，结果发现PD-L1、CTLA-4、EGFR和KRAS在高风险组中表达更高，HER-2和BRAF在低风险组中表达更高，提示高风险组可能对PD-L1 抑制剂、CTLA-4 抑制剂、EGFR 抑制剂和KRAS 抑制剂更敏感，而低风险组可能对HER2抗体和BRAF 抑制剂更敏感，间接预测了高低风险组患者对免疫治疗和靶向治疗的疗效差异。本研究存在一定局限性，仅分析了TCGA 数据库，未进行临床试验验证。

综上所述，本研究利用5 个m6A 调节因子在BC 中构建预后预测模型并验证其临床应用价值，评估预后预测模型中免疫细胞浸润程度，该模型可能成为评估BC 患者是否适合免疫治疗和靶向治疗的参考指标，为临床医师的治疗选择提供依据。