膀胱癌m6A相关长链非编码RNA预后模型的构建和免疫分析

Construction of a model for bladder cancer prognosis and immune analysis based on m6A-related long non-coding RNAs

摘要

目的 探究膀胱癌中N6-甲基腺苷（N6-methyladenosine，m6A）相关长非编码RNA（long non-coding RNA，lncRNA）及其与患者预后的关系。 方法 通过癌症基因组图谱（the cancer genome Atlas，TCGA）数据库获取转录、突变和临床数据，获得m6A相关的lncRNA表达数据，采用共表达网络分析、Cox回归分析和最小绝对收缩和选择算法（least absolute shrinkage and selection operator，LASSO）回归分析鉴定与m6A相关的lncRNA，并构建风险预后模型和对其进行验证。然后，本研究还构建了列线图来预测膀胱癌患者的预后。通过基因本体（gene ontology，GO）、京都基因与基因组百科全书（kyoto encyclopedia of genes and genomes，KEGG）研究生物学功能的差异。使用ESTIMATE算法计算基质、免疫评分。采用单样本基因集富集分析（single sample gene set enrichment analysis，ssGSEA）算法定量分析免疫细胞浸润和免疫功能。采用oncoPredict算法预测潜在治疗药物。 结果 本课题组确定了1 739个m6A相关lncRNA，其中 14个与预后相关。再通过LASSO回归分析筛选了5个m6A相关的lncRNA（GRK5-IT1、AC008883.2、AC145207.5、AC103746.1、AC104564.3）来构建预后预测模型。Kaplan-Meier和受试者工作特征曲线（receiver operating characteristic，ROC）曲线显示，该特征在训练集、验证集和全集中具有良好的预测能力。与其他临床特征相比，m6A相关lncRNA模型具有更高的诊断效率。免疫细胞浸润和ssGSEA分析进一步证实，特征lncRNA与膀胱癌患者的免疫状态显著相关。此外，高风险组对多种药物的治疗反应与低风险组存在显著差异。 结论 5个m6A相关的特征lncRNA可能有助于评估膀胱癌患者的预后和免疫特征，并指导个性化治疗方案的制定。

Abstract

Objective To explore N6-methyladenosine（m6A）-related long non-coding RNAs（lncRNAs） in bladder cancer（BCa） and their correlations with prognosis. Methods Transcription，mutation，and clinical data of m6A-related lncRNAs were downloaded from The Cancer Genome Atlas（TCGA） database. Then m6A-related signature lncRNAs were identified by co-expression network analysis，Cox regression analysis，and LASSO regression analysis. A risk prognostic model was constructed and validated. A nomogram was constructed to predict the prognosis of BCa patients. Functional enrichment analysis was performed using GO and KEGG. Stromal and immune scores were evaluated using the ESTIMATE algorithm. Single sample gene set enrichment analysis algorithm was used for quantitative analysis of immune cell infiltration and immune functions. The oncoPredict algorithm was used to predict potential therapeutic drugs. Results We identified 1，739 m6A-related lncRNAs，of which 14 were correlated with prognosis. Five m6A-related lncRNAs（GRK5-IT1，AC008883.2，AC145207.5，AC103746.1，and AC104564.3） were selected using LASSO regression to construct a prognostic risk model. The Kaplan–Meier and receiver operating characteristic curves confirmed that the prognostic risk model had good predictive performance in training，test，and complete sets. The prognostic risk model based on m6A-related lncRNAs had higher diagnostic efficiency compared to models based on other clinical features. Immune cell infiltration and single sample gene set enrichment analysis further confirmed that these signature lncRNAs were significantly associated with the immune status of BCa patients. Furthermore，the high-risk group and the low-risk group showed significant differences in the responses to multiple drugs. Conclusion The five m6A-related signature lncRNAs may help assess the prognosis and immune characteristics of BCa patients and facilitate the determination of personalized treatment options.

Graphical abstract

关键词

膀胱癌 / N6-甲基腺苷 / 长链非编码RNA / 癌症基因组图谱 / 预后模型

Key words

bladder cancer / m6A / lncRNA / TCGA / prognostic model

膀胱癌（bladder cancer，BC）是第十位常见的恶性肿瘤和第二常见的泌尿生殖系统恶性肿瘤，每年约有573 000例新发病例和213 000例死亡病例^[1]。在所有新发病例中，有大约四分之一被诊断为肌层浸润性膀胱癌（muscle invasive bladder cancer，MIBC)^[2]。MIBC比非肌层浸润性膀胱癌（non-muscle invasive bladder cancer，NMIBC）具有更高的侵袭性和更差的预后，其5年总生存率（overall survival，OS）低于50%^[3-4]。即使接受一线治疗，大约10%~30%的NMIBC患者也会进展为MIBC^[5-7]。因此，寻找新的经济可行且有效的预测生物标志物，对于改善膀胱癌的预后，实现个体化治疗具有重要意义。

长非编码RNA（long non-coding RNA，lncRNA）是1组长度超过200个核苷酸的非蛋白编码转录本^[8]。LncRNA能够参与基因表达调控以及多种病理过程，并在调节代谢相关基因的转录和翻译、充当竞争性内源性RNA（competitive endogenous RNA，ceRNA）等方面发挥着重要的功能作用^[9-11]。在癌症患者中，lncRNA的异常表达是一种常见的生物学现象，与患者的预后密切相关^[12]。

在众多RNA修饰中，N6-甲基腺苷（N6-methyladenosine，m6A）是最丰富的mRNA 修饰，平均每1 000个核苷酸中就含有1~2个m6A残基^[13-14]。m6A几乎参与RNA代谢的所有步骤，包括mRNA翻译、降解、剪接、输出和折叠^[15-16]。许多研究表明，m6A修饰在各种癌症中起着重要作用，通常是通过调节促癌或抑癌基因 mRNA中的m6A修饰水平来改变其表达水平，最终影响肿瘤的多种生物学行为^[17]。最近研究表明，在肿瘤发生发展过程中，m6A修饰在lncRNA的表达失调中起着不可或缺的作用^[18-19]，例如m6A诱导的lncDBET可通过FABP5介导的脂质代谢促进膀胱癌的恶性进展^[20]。

在本研究中，本课题组探究了m6A相关lncRNA在膀胱癌中的潜在应用价值。本研究表明m6A相关lncRNA可以用作膀胱癌预后预测和免疫治疗反应性预测的生物标志物，能为制定患者的个性化治疗方案提供参考依据，从而改善膀胱癌患者的临床获益。

1 材料与方法

1.1 数据获取

突变数据、转录组数据以及相应的临床信息均来自癌症基因组图谱（the cancer genome atlas，TCGA)（https://portal.gdc.cancer.gov/）数据库，包含412例肿瘤样本（tumor）和19例正常组织样本（normal）。从发表的文章中获得m6A调节因子：13个“读码器”基因（LRPPRC、RBMX、HNRNPA2B1、HNRNPC、FMR1、YTHDF3、YTHDF1、YTHDF2、IGF2BP2、IGF2BP3、IGF2BP1、YTHDC1和YTHDC2）、8个“编码器”基因（ZC3H13、METTL3、RBM15、RBM15B、VIRMA、WTAP、METTL16和METTL14）、2个“消码器”基因（FTO和 ALKBH5）^[21]。

1.2 LncRNA注释

lncRNA注释文件从GENCODE网站（https://www.gencodegenes.org）获得。在TCGA数据集中分别鉴定出16 876个lncRNA。

1.3 m6A相关lncRNA的鉴定

通过根据文献确定的23个公认的m6A调节因子，对每个数据集进行Pearson相关性分析，以识别与任一m6A调控因子具有显著相关性（根据|Pearson R|>0.4和P<0.001）的lncRNA。并使用R包“ggalluvial”^[22]进行可视化。

1.4 m6A相关lncRNA预后模型的构建

基于TCGA数据集，通过单因素Cox回归分析，本组确定了14个m6A相关预后lncRNA。进一步通过使用R包“glmnet”^[23]、“survival”和“survminer”进行LASSO Cox回归分析以及多因素Cox回归分析，最终得到5个m6A相关的特征lncRNA，并将其用于构建预后风险模型。风险评分计算如下：

风险评分=∑(Coef∗EXP),

Coef表示系数，EXP表示公式中每个预后相关lncRNA的表达水平。根据中位风险评分，将患者分为高风险组和低风险组，使用Kaplan-Meier生存曲线评估风险模型的预后预测能力，并使用R包“pROC”^[24]绘制受试者工作特征（receiver operating characteristic，ROC）曲线以评估曲线下面积（area under the curve，AUC）值。

1.5 列线图的构建和验证

为了确定膀胱癌的独立预后因素，使用单因素和多因素Cox回归分析评价了风险评分与多个临床特征（年龄、性别、WHO等级和分期）之间的预后关系。随后基于上述分析结果，构建具有综合预后特征的列线图，以预测膀胱癌患者的1年、3年和5年生存率，并绘制了校准曲线。

1.6 差异表达基因（differential expressed genes，DEGs）的功能富集分析

使用R包“limma”^[22]确定高风险组和低风险组之间的DEGs（|log2（fold change|>1，P<0.05）。使用R包“clusterProfiler”^[25]对DEGs进行基因本体（gene ontology，GO）和京都基因组百科全书（kyoto encyclopedia of genes and genomes，KEGG）通路富集分析。此外，还进行了基因组富集分析（gene set enrichment analysis，GSEA），以确定与m6A相关预后lncRNA富集的KEGG通路和肿瘤标志物。从MSigDB数据库^[26]中获取“c2.kegg.v7.4.symbols”和“c5.go.v7.4.symbols”背景基因集，并使用R包“clusterProfiler”^[25]绘制GSEA图。

1.7 肿瘤免疫分析

基于ESTIMATE^[27]算法，本课题组探索了两风险组之间免疫和基质细胞的丰度，并计算了每组的StromalScore、ImmuneScore和ESTIMATEScore（StromalScore+ImmuneScore）^[28]。此外，本组还研究了免疫检查点在两风险组中的差异表达。随后，通过R包 “GSVA”和“GSEABase”^[29]对膀胱癌中浸润免疫细胞和免疫相关功能进行单样本基因组富集分析（single-sample Gene Set Enrichment Analysis，ssGSEA）评分。

1.8 肿瘤突变负荷（tumor mutation burden，TMB）分析和潜在药物预测

从TCGA网站下载体细胞突变数据后，使用R包“maftools”^[30]整合并分析不同风险分组中TMB与生存率之间的关系。然后分析评估了不同风险分组之间免疫检查点抑制剂治疗的反应差异。最后，再使用R包“oncoPredict”^[31]来预测针对不同风险组可用于治疗的潜在药物。

1.9 统计学方法

使用R语言进行数据计算和统计分析（https://www.r-project.org/，4.0.2版本）。基于预后m6A相关lncRNA的表达情况，采用Kaplan-Meier生存曲线分析和log-rank检验来比较不同亚组间的OS。此外，进行单因素和多因素Cox回归分析来评估特征在预测OS方面的独立预后价值。采用Wilcoxon检验来确定两组间变量的差异。检验水准α=0.05。

2 结 果

2.1 鉴定膀胱癌中m6A相关lncRNAs

根据获得的23个m6A调节因子，本研究共鉴定了1 739个具有共表达关系的m6A相关的lncRNA（|Pearson R|>0.4和P<0.001）（图1）。使用单因素Cox分析（P<0.05）分析，共选择了14个差异预后相关lncRNA（GRK5-IT1、AC008883.2、AC111182.2、AC145207.5、AL121957.1、AC103746.1、AC016831.6、CLIP1-AS1、AC068282.1、AL133325.3、AC104564.3、AL731559.1、AC087286.4、AC058791.1）（图2A）。

2.2 m6A相关lcnRNA预测模型的构建

本研究将单因素Cox分析得到的14个预后相关lncRNA进行LASSO Cox回归分析，使用最佳lambda值确认最终的5个特征lncRNA，使用公式建立诊断模型（图2B~D）。公式如下：

风险评分=GRK5-IT1×（0.425798805732769）+AC008883.2×（1.06390457325097）+AC145207.5×（0.295538921677524）+AC103746.1×（0.539718400881476）+AC104564.3×（-0.57014487458029）

在全集、训练集和验证集中，根据风险评分中位数值，将患者分为低风险组和高风险组，2组间的风险分布如图3A~F，2 组间预后相关特征lncRNA表达如图3G~I。生存分析显示，无论是在全集还是在训练集或验证集中，与低风险组相比，高风险组的总生存率显著降低（图3J~L）。为了确认风险评分是否可以用作膀胱癌患者的独立预后因素，单因素及多因素Cox分析显示，年龄、分期和风险评分均与生存率相关（P<0.001）（图3M~N）。1年、3年和5年ROC的曲线下面积（AUC）分别为0.767、0.705和0.734（图3O）。在模型的5年ROC中，风险评分的AUC为0.767，明显高于其他临床病理特征（图3P）。风险模型中的10年C-index也高于其他临床特征（图3Q）。基于m6A相关lncRNA的风险评分可以作为膀胱癌的独立预后因素，表明m6A相关lncRNA在评估临床预后的过程中具有潜在价值。

2.3 m6A相关预后lncRNA的列线图构建和验证

此外，本课题组根据风险评分与多种临床特征（年龄、性别、WHO等级、分期、T分期、N分期和M分期）的关系，绘制了列线图，以预测膀胱癌患者的预后（图4A），在本研究的列线图中，校准曲线具有良好的预测1年、3年和5年生存率的能力（图4B）。本课题组建立的风险评分特征可以为预测这些预后指标的生存提供最有用和最准确的指导。接下来本研究分析了分期和年龄对生存率的影响，发现虽然高风险组的生存率均显著低于低风险组，但分期越高、年龄越大的患者生存率相对更低（图4C~F）。

2.4 m6A相关预后lncRNA功能富集分析

为了分析膀胱癌相关的差异表达m6A相关lncRNA与生物学过程（biological process，BP）、分子功能（molecular function，MF）、细胞组分（cellular component，CC）、生物学通路以及疾病之间的关系，本研究首先对膀胱癌相关的差异表达m6A相关lncRNA进行了功能富集分析。相关基因主要富集在endopeptidase activity、cytokine activity、glycosaminoglycan binding、serine hydrolase activity、serine-type endopeptidase activity等生物学过程中，同时富集在collagen-containing extracellular matrix、apical part of cell、cornified envelope、endoplasmic reticulum lumen、intermediate filament cytoskeleton等细胞组分中和epidermis development、skin development、external encapsulating structure organization、epidermal cell differentiation、extracellular structure organization等分子功能（图5A）。KEGG富集分析，m6A相关预后lncRNA主要富集在cytokine-cytokine receptor interaction、PI3K-Akt signaling pathway、IL-17 signaling pathway、estrogen signaling pathway、rheumatoid arthritis等细胞通路之中（图5B）。此外，GSEA分析表明，cytokine-cytokine receptor interaction、ecm receptor interaction、focal adhesion、pathways in cancer和regulation of actin cytoskeleton主要富集在高风险组中（图5C），而drug metabolism cytochrome、metabolism of xenobiotics、PPAR信号通路、retinol metabolism和starch and sucrose metabolism则主要富集在低风险组中（图5D）。

2.5 肿瘤免疫分析

在肿瘤微环境（tumor micro-environment，TME）评分方面，高风险组患者的免疫评分、ESTIMATE评分和基质评分均显著高于低风险组患者（图6A~C）。此外，本研究分析了高风险组和低风险组之间免疫检查点相关基因的表达，值得注意的是，大多数免疫检查点在高风险组患者中的表达较高，这可能解释了高风险组较差的生存率（图6D）。在免疫细胞气泡图中，本研究发现，高风险组的样本与Th2细胞、CD8⁺ T细胞、肿瘤相关成纤维细胞（cancer associated fibroblasts，CAF）、巨噬细胞、单核细胞等的浸润显著正相关，与中央记忆型T细胞、CD4⁺ T细胞等的浸润呈负相关（图6E)。接下来，本研究进一步进行ssGSEA分析，免疫细胞浸润方面，高风险组样本中的活化的树突状细胞（active dendritic cells，aDCs）、树突状细胞（dendritic cells，DCs）、巨噬细胞、中性粒细胞和浆细胞样树突状细胞（plasmacytoid dendritic cells，pDCs）等多种免疫细胞浸润均比低风险组样本更多（图6F）；同时免疫功能方面，高风险组中的抗原呈递细胞共抑制（APC_co_inhibition）、抗原呈递细胞共刺激（APC_co_stimulation）、CC趋化因子受体（CCR）、检查点（Check-point）和细胞溶解活性（Cytolytic_activity）等多种免疫功能均显著低于低风险组（图6G）。

2.6 TMB分析和潜在药物预测

从TGCA数据库下载了体细胞突变数据，并分析了高风险组和低风险组组体细胞突变的变化。突变程度前十的基因分别是：TP53、TTN、KMT2D、MUC16、ARID1A、KDM6A、PIK3CA、SYNE1、RYR2和KMT2C（图7A、B）。此外，更低的TMB和更高的风险评分的患者的预后相比于其他组更差（图7C、D）。

本研究使用R包“oncoPredict”评估了常见药物的半抑制浓度（half maximal inhibitory concentration，IC₅₀），IC₅₀值越低，药敏度越高，治疗效果越好。阿立塞替（Alisertib）、顺铂（Cisplatin）、达沙替尼（Dasatinib）、厄洛替尼（Erlotinib）、沙普替尼（Sapitinib）、星形孢菌素（Staurosporine）、赛沃替尼（Savolitinib）、他拉唑帕尼（Talazoparib）和曲美替尼（Trametinib）在高风险组中的IC₅₀明显低于低风险组，提示这些药物对于高风险组有更好的疗效（图8）。

3 讨 论

膀胱癌作为泌尿系统最常见的恶性肿瘤之一，其中约75%为非肌层浸润性膀胱癌，其余则大多为肌层浸润性膀胱癌^[2]。非肌层浸润性膀胱癌不仅具有较高的复发率，并且约10%～30%的患者最终会进展为肌层浸润性膀胱癌^[5-7]。而肌层浸润性膀胱癌患者5年生存率较低^[32]。膀胱癌的诊断和监测仍然依赖于膀胱镜活检，因此，探索新的治疗靶点和潜在的生物标志物至关重要。

LncRNA参与包括膀胱癌在内的多种肿瘤的生物学过程，例如生长、转移、耐药性和肿瘤免疫微环境^[33]。此外，研究证据表明lncRNA的m6A修饰与癌症预后和治疗反应有关^[34-35]，同时m6A相关lncRNA也被证明具有预后价值，可以预测患者的总生存期^[36-37]。然而，目前报道的m6A 修饰以lncRNA依赖性方式在膀胱癌中发挥作用的研究尚少。

本研究首先从TCGA数据库中鉴定出1 739个m6A相关的lncRNA，再通过单因素Cox分析得到了14个差异预后相关lncRNA，最终通过LASSO Cox回归分析确定了5个特征lncRNA，并构建了预后预测模型，该模型在训练集和测试集中均表现出良好的预测能力，多因素Cox回归分析结果显示，该模型可作为膀胱癌的独立预后因素。此外，通过列线图整合的独立危险因素，这些因素也能显示出足够的预后评估能力。以上结果表明，m6A相关lncRNA可用作潜在的预后生物标志物，人们可以更早发现高危患者，并通过更积极的监测和管理来实施有针对性的预防，也有助于医生针对不同患者制订个体化治疗方案。

本研究还通过对肿瘤免疫进行分析评分，高风险组的免疫评分、ESTIMATE评分和基质评分均高于低风险组，且多种免疫检查点基因在高风险组中显著高表达，此外如CD8⁺ T细胞、肿瘤相关成纤维细胞、巨噬细胞等多种免疫细胞均在高风险组中存在较高的浸润；本课题组又通过ssGSEA分析发现，高风险组中16种免疫细胞的浸润水平显著高于低风险组。

近年来，针对程序性细胞死亡的免疫检查点抑制剂（ICI），包括PD-1、程序性死亡配体1（PD-L1）和细胞毒性T淋巴细胞相关蛋白4（CTLA4），已成为癌症治疗的另一大支柱^[38]。然而ICI治疗并不是一刀切的治疗方法，其疗效高度依赖于肿瘤的肿瘤免疫微环境特征^[39]。TMB已在多项临床试验中进行了研究，在PURE-01试验中，TMB与帕博利珠单抗治疗的反应率呈正相关；在ABACUS试验（阿替利珠单抗）、NABUCCO试验（伊匹木单抗联合纳武单抗）中，与治疗无反应患者相比，反应患者的TMB数值更高^[40-43]。此外在PURE-01试验中，14名患者在接受帕博利珠单抗治疗后TMB较治疗前基线数值有明显的下降^[42]。TME中的免疫细胞浸润可能会影响肿瘤患者的存活、转移和治疗耐药性^[44-46]。有研究表明，肿瘤的中性粒细胞浸润与不良预后相关^[47-48]，而Treg浸润水平高的患者预后更好^[49]。

然而本研究仍然存在局限性，虽然本课题组使用了1个测试集来评估风险评分模型，但该测试集中的样本较少。其次，在未来的研究中，还需要体内和体外实验来进一步确定膀胱癌中m6A 相关lncRNA发挥作用的机制和具体信号通路。

参考文献

[1] Sung H， Ferlay J， Siegel RL，et al. Global cancer statistics 2020：GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries[J]． CA Cancer J Clin，2021，71（3）：209-249．

[2] Siegel RL， Miller KD， Jemal A. Cancer statistics，2016[J]. CA Cancer J Clin，2016，66（1）：7-30．

[3] Patel VG， Oh WK， Galsky MD. Treatment of muscle-invasive and advanced bladder cancer in 2020[J]. CA A Cancer J Clinicians，2020，70（5）：404-423．

[4] Ghandour R， Singla N， Lotan Y. Treatment options and outcomes in nonmetastatic muscle invasive bladder cancer[J]. Trends Cancer，2019，5（7）：426-439．

[5] Kamat AM， Hahn NM， Efstathiou JA，et al. Bladder cancer[J]. Lancet，2016，388（10061）：2796-2810．

[6] Sanli O， Dobruch J， Knowles MA，et al. Bladder cancer[J]. Nat Rev Dis Primers，2017，3：17022．

[7] Mari A， Campi R， Tellini R，et al. Patterns and predictors of recurrence after open radical cystectomy for bladder cancer：a comprehensive review of the literature[J]. World J Urol，2018，36（2）：157-170．

[8] Romero-Barrios N， Legascue MF， Benhamed M，et al. Splicing regulation by long noncoding RNAs[J]. Nucleic Acids Res，2018，46（5）：2169-2184．

[9] Dhamija S， Diederichs S. From junk to master regulators of invasion：lncRNA functions in migration，EMT and metastasis[J]. Int J Cancer，2016，139（2）：269-280．

[10] Peng WX， Koirala P， Mo YY. LncRNA-mediated regulation of cell signaling in cancer[J]. Oncogene，2017，36（41）：5661-5667．

[11] Quinn JJ， Chang HY. Unique features of long non-coding RNA biogenesis and function[J]. Nat Rev Genet，2016，17：47-62．

[12] Ji P， Diederichs S， Wang WB，et al. MALAT-1，a novel noncoding RNA，and thymosin beta4 predict metastasis and survival in early-stage non-small cell lung cancer[J]. Oncogene，2003，22（39）：8031-8041．

[13] He LE， Li HY， Wu AQ，et al. Functions of N6-methyladenosine and its role in cancer[J]. Mol Cancer，2019，18（1）：176．

[14] Krug RM， Morgan MA， Shatkin AJ. Influenza viral mRNA contains internal N6-methyladenosine and 5'-terminal 7-methylguanosine in cap structures[J]. J Virol，1976，20（1）：45-53．

[15] Liu N， Zhou KI， Parisien M，et al. N6-methyladenosine alters RNA structure to regulate binding of a low-complexity protein[J]. Nucleic Acids Res，2017，45（10）：6051-6063．

[16] Liu Q， Gregory RI. RNAmod：an integrated system for the annotation of mRNA modifications[J]. Nucleic Acids Res，2019，47（W1）：W548-W555．

[17] Wang TY， Kong S， Tao M，et al. The potential role of RNA N6-methyladenosine in Cancer progression[J]. Mol Cancer，2020，19（1）：88．

[18] Meng XR， Deng YY， He SH，et al. m⁶A-mediated upregulation of LINC00857 promotes pancreatic cancer tumorigenesis by regulating the miR-150-5p/E2F3 axis[J]. Front Oncol，2021，11：629947．

[19] Yang X， Zhang S， He CY，et al. METTL14 suppresses proliferation and metastasis of colorectal cancer by down-regulating oncogenic long non-coding RNA XIST[J]. Mol Cancer，2020，19（1）：46．

[20] Liu PH， Fan BY， Othmane B，et al. m⁶A-induced lncDBET promotes the malignant progression of bladder cancer through FABP5-mediated lipid metabolism[J]. Theranostics，2022，12（14）：6291-6307．

[21] Jin Y， Wang ZW， He D，et al. Analysis of m6A-related signatures in the tumor immune microenvironment and identification of clinical prognostic regulators in adrenocortical carcinoma[J]. Front Immunol，2021，12：637933．

[22] Ritchie ME， Phipson B， Wu D，et al. Limma powers differential expression analyses for RNA-sequencing and microarray studies[J]. Nucleic Acids Res，2015，43（7）：e47．

[23] Friedman J， Hastie T， Tibshirani R. Regularization paths for generalized linear models via coordinate descent[J]. J Stat Softw，2010，33（1）：1-22．

[24] Robin X， Turck N， Hainard A，et al. pROC：an open-source package for R and S+ to analyze and compare ROC curves[J]. BMC Bioinformatics，2011，12：77．

[25] Yu GC， Wang LG， Han YY，et al. clusterProfiler：an R package for comparing biological themes among gene clusters[J]. OMICS，2012，16（5）：284-287．

[26] Subramanian A， Tamayo P， Mootha VK，et al. Gene set enrichment analysis：a knowledge-based approach for interpreting genome-wide expression profiles[J]. Proc Natl Acad Sci U S A，2005，102（43）：15545-15550．

[27] Yoshihara K， Shahmoradgoli M， Martínez E，et al. Inferring tumour purity and stromal and immune cell admixture from expression data[J]. Nat Commun，2013，4：2612．

[28] Chen BB， Khodadoust MS， Liu CL，et al. Profiling tumor infiltrating immune cells with CIBERSORT[J]. Methods Mol Biol，2018，1711：243-259．

[29] Hänzelmann S， Castelo R， Guinney J. GSVA：gene set variation analysis for microarray and RNA-seq data[J]. BMC Bioinformatics，2013，14：7．

[30] Mayakonda A， Lin DC， Assenov Y，et al. Maftools：efficient and comprehensive analysis of somatic variants in cancer[J]. Genome Res，2018，28（11）：1747-1756．

[31] Maeser D， Gruener RF， Huang RS. oncoPredict：an R package for predicting in vivo or cancer patient drug response and biomarkers from cell line screening data[J]. Brief Bioinform，2021，22（6）：bbab260．

[32] Miller KD， Nogueira L， Devasia T，et al. Cancer treatment and survivorship statistics，2022[J]. CA Cancer J Clin，2022，72（5）：409-436．

[33] Li Y， Li G， Guo X，et al. Non-coding RNA in bladder cancer[J]. Cancer Lett，2020，485：38-44．

[34] Wen SM， Wei YL， Zen C，et al. Long non-coding RNA NEAT1 promotes bone metastasis of prostate cancer through N6-methyladenosine[J]. Mol Cancer，2020，19（1）：171．

[35] Dong PX， Xiong Y， Yue JM，et al. Long non-coding RNA NEAT1：a novel target for diagnosis and therapy in human tumors[J]. Front Genet，2018，9：471．

[36] Luo YJ， Zhang ZH， Zheng B，et al. Comprehensive analyses of N⁶-methyladenosine-related long noncoding RNA profiles with prognosis，chemotherapy response，and immune landscape in small cell lung cancer[J]. Cancer Sci，2022，113（12）：4289-4299．

[37] Zhang M， Zhang J， Liu Y. Comprehensive analysis of molecular features，prognostic values，and immune landscape association of m6A-regulated immune-related lncRNAs in smoking-associated lung squamous cell carcinoma[J]. Front Genet，2022，13：887477．

[38] Peyrottes A， Ouzaid I， Califano G，et al. Neoadjuvant immunotherapy for muscle-invasive bladder cancer[J]. Medicina（Kaunas），2021，57（8）：769．

[39] van Dorp J， van der Heijden MS. The bladder cancer immune micro-environment in the context of response to immune checkpoint inhibition[J]. Front Immunol，2023，14：1235884．

[40] Powles T， Kockx M， Rodriguez-Vida A，et al. Clinical efficacy and biomarker analysis of neoadjuvant atezolizumab in operable urothelial carcinoma in the ABACUS trial[J]. Nat Med，2019，25（11）：1706-1714．

[41] Gao JJ， Navai N， Alhalabi O，et al. Neoadjuvant PD-L1 plus CTLA-4 blockade in patients with cisplatin-ineligible operable high-risk urothelial carcinoma[J]. Nat Med，2020，26（12）：1845-1851．

[42] Necchi A， Anichini A， Raggi D，et al. Pembrolizumab as neoadjuvant therapy before radical cystectomy in patients with muscle-invasive urothelial bladder carcinoma（PURE-01）：an open-label，single-arm，phase Ⅱ study[J]. J Clin Oncol，2018，36（34）：3353-3360．

[43] van Dijk N， Gil-Jimenez A， Silina K，et al. Preoperative ipilimumab plus nivolumab in locoregionally advanced urothelial cancer：the NABUCCO trial[J]. Nat Med，2020，26（12）：1839-1844．

[44] Joyce JA， Fearon DT. T cell exclusion，immune privilege，and the tumor microenvironment[J]. Science，2015，348（6230）：74-80．

[45] Hui LL， Chen Y. Tumor microenvironment：sanctuary of the devil[J]. Cancer Lett，2015，368（1）：7-13．

[46] Waniczek D， Lorenc Z， Śnietura M，et al. Tumor-associated macrophages and regulatory T cells infiltration and the clinical outcome in colorectal cancer[J]. Arch Immunol Ther Exp，2017，65（5）：445-454．

[47] Wu XY， Qu DF， Weygant N，et al. Cancer stem cell marker DCLK1 correlates with tumorigenic immune infiltrates in the colon and gastric adenocarcinoma microenvironments[J]. Cancers，2020，12（2）：274．

[48] Buonaguro L， Tagliamonte M. Selecting target antigens for cancer vaccine development[J]. Vaccines，2020，8（4）：615．

[49] Wu ZL， Zhu KJ， Liu QG，et al. Profiles of immune infiltration in bladder cancer and its clinical significance：an integrative genomic analysis[J]. Int J Med Sci，2020，17（6）：762-772．