{'fraction_non_alphanumeric': 0.04053098556211938, 'fraction_numerical': 0.022854275482973335, 'mean_word_length': 5.198784108581931, 'pattern_counts': {'":': 0, '<': 32, '<?xml version=': 0, '>': 3, 'https://': 4, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Integrative transcriptional characterization of cell cycle checkpoint genes promotes clinical management and precision medicine in bladder carcinoma

Background The aberrant regulation of cell cycle is significantly correlated with cancer carcinogenesis and progression, in which cell cycle checkpoints control phase transitions, cell cycle entry, progression, and exit. However, the integrative role of cell cycle checkpoint-related genes (CRGs) in bladder carcinoma (BC) remains unknown. Methods The transcriptomic data and clinical features of BC patients were downloaded from The Cancer Genome Atlas (TCGA), used to identify CRGs correlated with overall survival (OS) by univariate Cox regression analysis. Then, the multivariate and least absolute shrinkage and selection operator (LASSO) Cox regression analyses further developed a prognostic CRG signature, which was validated in three external datasets retrieved from Gene Expression Omnibus (GEO). The receiver operating characteristic curve (ROC) analysis was conducted for evaluating the performance of the CRG signature in prognosis prediction. RNA sequencing (RNA-Seq) was performed to explore the expression difference in the identified CRGs between tumor and normal tissue samples from 11 BC patients in the local cohort. Ultimately, genomic profiles and tumor microenvironment (TME), and the Genomics of Drug Sensitivity in Cancer (GDSC) were investigated to guide precision treatment for BC patients with different CRG features. Results The novel constructed 23-CRG prognostic signature could stratify BC patients into high-risk and low-risk groups with significantly different outcomes (median OS: 13.64 vs. 104.65 months). Notably, 19 CRGs were the first to be identified as being associated with BC progression. In three additional validation datasets (GSE13507, GSE31684, and GSE32548), higher CRG scores all indicated inferior survival, demonstrating the robust ability of the CRG signature in prognosis prediction. Moreover, the CRG signature as an independent prognostic factor had a robust and stable risk stratification for BC patients with different histological or clinical features. Then, a CRG signature-based nomogram with a better performance in prognostic prediction [concordance index (C-index): 0.76] was established. Functional enrichment analysis revealed that collagen-containing extracellular matrix (ECM), and ECM-related and MAPK signaling pathways were significantly associated with the signature. Further analysis showed that low-risk patients were characterized by particularly distinctive prevalence of FGFR3 (17.03% vs. 6.67%, p < 0.01) and POLE alterations (7.97% vs. 2.50%, p < 0.05), and enrichment of immune infiltrated cells (including CD8+ T cells, CD4+ naïve T cells, follicular helper T cells, Tregs, and myeloid dendritic cells). RNA-seq data in our local cohort supported the findings in the differentially expressed genes (DEGs) between tumor and normal tissue samples, and the difference in TME between high-risk and low-risk groups. Additionally, CRG signature score plus FGFR3 status divided BC patients into four molecular subtypes, with distinct prognosis, TME, and transcriptomic profiling of immune checkpoint genes. Of note, CRG signature score plus FGFR3 status could successfully distinguish BC patients who have a higher possibility of response to immunotherapy or chemotherapy drugs. Conclusions The CRG signature is a potent prognostic model for BC patients, and in combination with FGFR3 alterations, it had more practical capacity in the prediction of chemotherapy and immunotherapy response, helping guide clinical decision-making.


Introduction
Cell cycle checkpoint-related genes (CRGs) play pivotal roles in cell cycle progression (CCP), ensuring and control regulating cell cycle events (1). Generally, in eukaryotic cells, the mitotic cell cycle is composed of two stages, the interphase (G1, S, and G2) and the mitotic (M) phase. The gap phases of G1-to-S (2), Sto-G2 (3), and G2-to-M (4), likely as decision windows, can determine cell cycle entry and progression. In cancer cells, some CRGs preventing DNA damage are usually compromised, contributing to genetic alterations and genomic instability (5), but those CRGs involved in DNA replication stress are scarcely altered to endure the replication stress (6). On the other hand, cancer cells could potentiate DNA replication stress through transcriptional regulation of CRGs (7). Overall, CRGs and CRGrelated signaling pathways play a key role in regulating the phase transitions, CCP, and cell cycle entry and exit in cancer cells (8).
In previous studies, it has been revealed that aberrantly expressed CRGs might be an essential prerequisite for cells to become cancerous, leading to tumor development and progression. For instance, CRG ataxia telangiectasia mutated (ATM) plays a core role in responding to DNA damage and stimulating DNA repair signaling pathways, and its absence is highly prone to giving rise to carcinogenesis (9). Of note, its downregulation or inactivation is associated with the highly accumulated genomic aberrations, which is one of the hallmarks of cancer (10). Ataxia telangiectasia and Rad3-related (ATR), another CRG encoding protein kinase that regulates DNA damage response (DDR), is correlated with the polarization of M2 tumor-associated macrophages, lymph node metastasis, and poor prognosis of patients with nasopharyngeal carcinoma (11). The abnormal expression of CRGs has been found to be associated with the development and progression of multiple kinds of cancers, including melanoma (12), lung cancer (13), colorectal cancer (14), and hepatocellular carcinoma (15). A growing body of evidence has found that the dysregulation of CRGs could render cells to be cancerous and promote cancer cell proliferation, but most studies only disclose the role of a single CRG in cell cycle, tumor carcinogenesis, and progression. Currently, the role of integrated CRGs representing checkpoint mechanisms in the regulation of cell cycle in tumor carcinogenesis and progression remains to be fully delineated.
Previous studies have already expounded that dysregulation of CRGs is correlated with increased genomic instability and malignant progression in bladder carcinoma (BC) patients (16- 18), indicating that there is a great potential of CRGs to become prognostic or targeted biomarkers for BC patients. BC generally presents as non-muscle-invasive bladder cancer (NMIBC), muscle-invasive bladder cancer (MIBC), or metastatic BC, of which the MIBC subtype has a relatively worse prognosis and poor treatment responses (19). Moreover, BC is a molecularly heterogeneous cancer with divergent clinical outcomes (20). The heterogeneity of tumor is always a huge challenge for cancer management and could reduce the efficacy of molecularly targeted therapies; therefore, the dissection of molecular signatures is urgently needed. Based on transcriptomic profiling of CRGs, a novel signature was constructed in the present study with good performance in prognosis prediction and that could function as a biomarker for treatment response. Overall, this study was of guiding significance in the clinical management of BC patients and promoting precision treatment.

Data acquisition
In the current study, RNA-seq data (TPM: transcripts per million) and related clinical features of BC patients involved in the TCGA-BC cohort were collected (https://www.cbioportal.org/) as the training cohort (BC samples with no available RNA-seq data or survival information were excluded in subsequent analyses). The transcriptional profiles together with clinical features of GSE13507 (including 62 MIBC and 103 NMIBC patients), GSE31684 (including 79 MIBC and 14 NMIBC patients), and GSE32548 (including 131 BC patients; the muscle-invasive status was uncertain) were downloaded, as three independent validation groups, from GEO (https://www.ncbi.nlm.nih.gov/geo/). The tumor and normal tissue samples from the TCGA, GSE133624, GSE188715, and GSE13507 cohorts were retrieved to conduct the differentially expressed gene (DEG) analysis. The RNA-seq data were normalized by Log 2 (x+0.001). Cell cycle checkpoint gene sets were obtained from the database of Gene Ontology (GO), Biological Process (http://geneontology.org/), and Reactome (https:// reactome.org/). The term "cell cycle checkpoint" was used for the acquisition of CRGs. After merging, a total of 464 CRGs were included in the union set, which were then selected for the further processes of establishing a prognostic signature.

Identification of the prognostic CRG signature
In order to investigate whether transcriptomic characterization was associated with prognosis of BC patients, the unsupervised hierarchical clustering analysis was conducted by using the package "Fastcluster", dividing BC patients into different clusters. The principal component analysis (PCA) was performed via the "ggbiplot" package, further revealing the distinction between clusters. The Kaplan-Meier curve analysis showed the overall survival (OS) of patients between clusters. Subsequently, the univariate Cox proportional hazard regression analysis, by using the packages "rms" and "survival", was conducted to classify the relationship of OS and the expression of each CRG in the TCGA-BC cohort. The least absolute shrinkage and selection operator (LASSO) and multivariate Cox proportional hazard regression analyses were used to construct a prognostic CRG signature via the "glmnet" package. Subsequently, a risk formula was established, and the CRG score was generated for each BC patient with the following formula: CRG score = expression value of gene 1 × C 1 + expression value of gene 2 × C 2 + … + expression value of gene x × C x , where C x is the coefficient of gene x. The optimal CRG score was adapted as the cutoff value by using the package "maxstat" to divide BC patients into high-risk and low-risk groups, while the same method was utilized in the validation groups. The value of area under the receiver operating characteristic (ROC) curve (AUC) revealed the prognostic performance of specificity and sensitivity. All analyses were also applied in the validation groups.

RNA sequencing in the local BC cohort
Tumor and matched normal tissues were collected from 11 MIBC patients in our local cohort to perform RNA-seq. This study was approved by the ethics committee of the Chinese PLA General Hospital (S2019-302-01), and all enrolled patients have signed the informed consent. The total RNA of each sample was collected using a FastPure ® Cell/Tissue Total RNA Isolation Kit V2 (Vazyme, Jiangsu, China), and its concentration and RNA Integrity Number (RIN) were determined by using Qubit (Thermo Fisher Scientific, MA, United States) and an Agilent 2100 bioanalyzer (Agilent Technologies, CA, United States), respectively. One sample of normal tissue failed quality control and then was discarded. Library construction was conducted using the NEBNext ® Ultra ™ RNA Library Prep Kit for Illumina ® Kit (NEB, MA, United States) and finally sequenced on the Illumina Novaseq-6000 system (Illumina, MA, United States).

Construction of the predictive nomogram
The univariate and multivariate Cox proportional hazard regression analyses, by using the packages "rms" and "survival", were conducted to explore whether the CRG signature was an independent prognostic factor. Furthermore, the OS-related clinical features were selected for the construction of a CRG signature-based nomogram. The value of AUC was used to evaluate the performance of a novel constructed nomogram in prognosis prediction, and the calibration plots were built to perform the consistency between actual OS and predicted OS by using the package "rms".

Functional enrichment analysis
The Pearson correlation analysis was conducted to screen out the CRG score-associated genes by the cutoff criteria of | Pearson Correlation Coefficient| > 0.4 and p < 0.05, and the heatmap showing the expression of the identified genes was drawn by using the package "pheatmap". The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were conducted based on the expression of CRG score-associated genes.

Tumor microenvironment analysis
The stromal score and immune score could predict the infiltration levels of stromal and immune cells, respectively. Moreover, the ESTIMATE (Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data) algorithm was used to determine the tumor purity in tumor tissues (21). The CIBERSORT (Cell type Identification By Estimating Relative Subsets Of known RNA Transcripts) algorithm was a widely used method to characterize the tumor-infiltrated lymphocytes (TILs) inside tumor tissues based on the gene expression profile (22); thus, the abundance of TILs between different groups was compared by the evaluation of CIBERSORT. The T-cell dysfunction score, exclusion score, and Tumor Immune Dysfunction and Exclusion (TIDE) score were analyzed to estimate the tumor immune escape (23). Furthermore, the expression levels of immune checkpoint genes, such as PD-1, PD-L1, CTLA-4, TIM-3, and LAG3 (24), were investigated to explore the potential immune therapies for different BC patients.

Evaluation of therapeutic treatment responses
The data were downloaded from the Genomics of Drug Sensitivity in Cancer (GDSC, https://www.cancerrxgene.org/), a publicly available pharmacogenomic database, to predict the treatment response of chemotherapy for BC patients. The measuring parameter of half-maximal inhibitory concentration (IC 50 ) was used to estimate the chemotherapeutic treatment response and chemo-drug sensitivity between different groups via the "pRRophetic" package. In addition, the clinical information and transcriptomic data of 348 urothelial cancer (UC) patients involved in the IMvigor210 cohort receiving anti-PD-L1 therapy (Atezolizumab) were downloaded from the following website (25): http://research-pub.gene.com/ IMvigor210CoreBiologies/. The criteria of treatment response were defined as previously described: CR: complete response, PR: partial response, SD: stable disease, PD: progressive disease. In addition, the GSE176307 dataset (26), including 86 UC patients with wild-type FGFR3 and 17 patients with altered FGFR3 receiving anti-PD-1 or anti-PD-L1 treatment, was further employed to explore the response prediction ability of the CRG signature for BC patients with different molecular subtypes by immunotherapeutic treatment.

Statistical analysis
The statistical data were analyzed by the Kruskal-Wallis (K-W) test, Mann-Whitney U test, Chi-square test, and Fisher's exact test in R studio. Cluster analysis was performed by the unsupervised hierarchical clustering. The univariate and multivariate Cox proportional hazards regression models were employed to assess the hazard ratio of the signature and clinical features. The Kaplan-Meier curve analysis along with log-rank test was conducted to evaluate the clinical outcomes of BC patients. The statistically significant difference was determined by "p < 0.05".

Results
The expression profiling of CRGs correlated with OS of BC patients The design of this study was exhibited in a work flowchart ( Figure 1). Initially, by way of the unsupervised hierarchical clustering analysis, BC patients from the TCGA cohort were divided into two clusters (Figure 2A), which could be separated in the Dim2 axis ( Figure 2B), indicating that the expression profiling of CRGs between two clusters was noticeably distinct. Of note, BC patients in cluster 1 had worse OS (median OS: 32.02 vs. 64.80 months, p = 0.088, Figure 2C). The above results suggested that the expression profiling of CRGs potentially affected the clinical outcomes in BC.

Identification and validation of the prognostic CRG signature
Subsequently, a total of 398 BC patients (Table 1) were selected into the univariate Cox proportional hazard regression analysis to investigate the relationship between OS and the expression of 464 CRGs (Supplementary Table 1), and eventually, expression levels of 52 CRGs were significantly correlated with the OS of BC patients (Supplementary Table 2). Subsequently, the LASSO and multivariate Cox proportional hazard regression analyses determined a novel prognostic CRG signature, consisting of 23 CRGs ( Figure 2D, Table 2). The scatterplot demonstrated the distribution of CRG score and the corresponding survival status of BC patients in the training group, in which BC patients were divided into high-risk and low-risk groups ( Figure 2E). It was identified that BC patients in the high-risk group had a significantly worse OS than those in the low-risk group (median OS: 13.64 vs. 104.65 months, p < 0.0001, Figure 2F). The AUC values of the CRG signature for predicting OS at 1, 2, 3, 4, and 5 years were 0.78, 0.75, 0.76, 0.76, and 0.80, respectively ( Figure 2G), indicating good specificity and sensitivity of the CRG signature in predicting the prognosis of BC patients.
Three additional independent cohorts, namely, GSE13507, GSE31684, and GSE32548, were further employed to validate the predictive ability of the prognostic CRG signature (Figures 3A-F). In each validation cohort, the CRG signature had good risk stratification, and a higher CRG score indicated a significantly shorter OS (

Clinical association of the prognostic CRG signature
Underlying the CRG signature, a significantly higher proportion of female patients as well as patients diagnosed at age over 68.5 years old were observed in the high-risk group (p < 0.05, Figure 4A). Moreover, there was also a significantly increased proportion of BC patients with advanced T (tumor stage), N (lymph node status), M (metastasis), and clinical stages (p < 0.05) in the high-risk group. Overall, it was demonstrated that a higher CRG score indicated advanced T, N, M, and/or clinical stages. In addition, stratification by age, gender, TNM, or clinical stages in the TCGA cohort further revealed that the CRG signature had a good predictive prognosis ability, and the higher CRG score was invariably correlated with worse survival ( Figure 4B).

Construction of a nomogram based on the prognostic CRG signature
Apparently, the CRG signature was the most robust risk factor, by comparison with classical clinical features ( Figure 5A), and it was the only independent prognostic factor for BC patients (p < 0.01, Figure 5B). Subsequently, the CRG signature-based prognostic nomogram was constructed in combination with several clinical parameters together, including diagnosis age, and T, N, and M stage ( Figure 5C). The C-index of the novel constructed prognostic nomogram was 0.76, with 95% confidence interval ranging from 0.70 to 0.81. In addition, the AUC values of the nomogram for predicting OS at

Exploration of the CRG signature-related biological functions
The highly CRG score-associated genes were defined with the criteria of |Pearson Correlation Coefficient| > 0.4 and p < 0.05, and a total of 252 and 159 genes were found to be positively or negatively correlated with CRG score, respectively ( Figure 6A). The extracellular matrix (ECM)-related biological functions, focal adhesion, regulation of actin cytoskeleton, and MAPK signaling pathway were mainly enriched ( Figures 6B, C).

Genomic characteristics and TME underlying the CRG signature
The genomic alteration profiles, respectively in high-risk and low-risk groups, showed that the most prevalently altered genes were distinct (Figures 7A, B). Through statistical analysis, among the prevalently altered genes between high-risk and low-risk groups (the altered genes with prevalence ≤ 5.00% in both high-risk and low-risk groups were excluded), there was a total of 46 altered genes with higher alteration frequencies enriched in the low-risk group, including RYR2 (frequency: 20.29% vs. 11.67%, p < 0.05), FAT4 (frequency: 17.75% vs. 9.17%, p < 0.05), and FGFR3 (frequency: 17.03% vs. 6.67%, p < 0.01), whereas the high-risk group had a significantly higher prevalence of only 5 altered genes, namely, RB1 (25.00% vs. 15.94%, p < 0.05), FBXW7 (12.50% vs. 5.80%, p < 0.05), NFE2L2 The transcriptomic profiling of cell cycle checkpoint-related genes (CRGs) was correlated with overall survival (OS) of BC patients from the TCGA-BC cohort. (10.83% vs. 3.99%, p < 0.05), ASAP2 (5.83% vs. 1.45%, p < 0.05), and PCSK6 (5.83% vs. 1.09%, p < 0.05) ( Figure 7C). Notably, only one DDR-related gene POLE was found with higher alteration frequency in the low-risk group (frequency: 7.97% vs. 2.50%, p < 0.05, Figure 7D). Of note, BC patients with altered POLE had a trend of better clinical outcomes (median OS: 33.14 months vs. not reached, p = 0.13, Figure 7E). Subsequently, the molecular characterization of TME was presented, demonstrating that the high-risk group had a significantly increased stromal score and ESTIMATE score ( Figure 8A). Moreover, the high-risk group also had a significantly higher T-cell exclusion score and TIDE score ( Figure 8B). Additionally, more plasma B cells, CD8+ T cells, CD4+ naïve T cells, follicular helper T cells, Tregs, and activated myeloid dendritic cells were significantly more enriched in the low-risk group, whereas the high-risk group only had significantly higher infiltration levels of macrophage M0 and M2 ( Figure 8C).

Exploration of CRG signature or related genes in the local BC cohort
It was found that the DEGs between tumor and normal tissue samples significantly varied in public datasets and our local cohort, and remarkably, only 7 of the 23 discovered CRGs exhibited significant differences between tumor and normal tissue samples across the TCGA, GSE133624, GSE188715, and GSE13507 cohorts ( Figure 9A). Among which, SLC25A15, RAD9A, PRF19, THOC1, and TIPIN were significantly overexpressed in tumor tissues in both the TCGA ( Figure 9B) and our local cohorts ( Figure 9C). PPP2CB and FBXO31 were upregulated in normal tissues in the TCGA cohort, but our local samples did not differ significantly. Limited by the sample size, there was no statistically significant difference in TILs between high-risk and low-risk groups in our local cohort, whereas, a consistent trend of more CD8+ T cells, follicular helper T cells, Tregs, and activated dendritic cells were presented in the lowrisk group and relatively more macrophage M0 and M2 were presented in the high-risk group in our local cohort ( Figure 9D). In our local cohort, the ESTIMATE score was found to be relatively higher but not statistically significant in the high-risk group ( Figure 9E), but notably, high-risk individuals in our local BC cohort had significantly higher TIDE scores ( Figure 9F).

Clinical significance of CRG signature in BC patients with altered FGFR3
Regarding the CRG score and FGFR3 alteration status, BC patients were further divided into four different risk groups. It could be apparently seen that high-risk patients with wild-type FGFR3 had the highest CRG scores (p < 0.0001, Figure 10A). Moreover, patients with wild-type FGFR3 had significantly higher CRG scores than those with FGFR3 alterations (FGFR3mt-high vs. FGFR3wt-high and FGFR3mt-low vs. FGFR3wt-low, p < 0.0001). Expectedly, high-risk patients with wild-type FGFR3 had the significantly shortest OS ( Figure 10B), further demonstrating the robust ability of CRG signature in prognosis prediction.
Immunity of BC patients with altered FGFR3 underlying the CRG signature Compared to other groups, low-risk patients with altered FGFR3 had the significantly lowest stromal score, immune score, ESTIMATE score, and TIDE score ( Figures 10C, D), as well as a relatively higher enrichment of plasma B cells, activated NK cells, CD8+ T cells, and follicular helper T cells ( Figure 10E). However, among the other three groups, there was no statistically significant difference in ESTIMATE and TIDE  The prognostic CRG signature was validated in three independent datasets: GSE13507, GSE31684, and GSE32548. Kaplan-Meier curve analysis to compare OS between high-risk and low-risk groups in GSE13507 (A), GSE31684 (C), and GSE32548 (E). The receiver operating characteristic curve (ROC) analysis for validating the ability of CRG signature in prognosis prediction in GSE13507 (B), GSE31684 (D), and GSE32548 (F). Comparison analysis of the CRG scores between MIBC and NMIBC patients in GSE13507 (G) and GSE31684 (H).
scores ( Figures 10C, D), except that low-risk patients with wildtype FGFR3 had a significantly lower TIDE score than high-risk patients with wild-type FGFR3 (p < 0.05, Figure 10D). Remarkably, the distribution of TILs in high-risk patients with altered FGFR3 was similar to that of BC patients with wild-type FGFR3, regardless of whether they had a high or low CRG score ( Figure 10E). Interestingly, low-risk patients with altered FGFR3 had a lower expression of immune checkpoint genes, including PD-L1, CTLA4, TIM3, LAG3, PD-1, CD27, CD47, and IDO1 (p < 0.05); on the contrary, high-risk patients with altered FGFR3 had a higher expression of TIM3 and CD47 (p < 0.05), whereas the expression level of PD-L1 in BC patients with altered FGFR3 was quite low ( Figure 10F).

Estimate of chemotherapeutic treatment response
As investigated for commonly used chemotherapeutic drugs for BC patients (including cisplatin, docetaxel, paclitaxel, methotrexate, and doxorubicin), it was found that IC 50 values for the response prediction of chemotherapeutic treatment by cisplatin, docetaxel, paclitaxel, methotrexate, or doxorubicin significantly differed between patients with distinctive molecular subtypes ( Figures 11A-E). Regardless of BC patients (whether FGFR3 was altered or not), the higher CRG score indicated the relatively higher sensitivity to cisplatin and docetaxel ( Figures 11A, B). Moreover, compared to low-risk BC patients with altered FGFR3, those with wild-type FGFR3 seemed to be more sensitive to paclitaxel ( Figure 11C). Meanwhile, low-risk patients with altered FGFR3 had the highest sensitivity to methotrexate ( Figure 11D). For patients with wild-type FGFR3, those in the high-risk group were more sensitive to doxorubicin ( Figure 11E).

Prediction of immunotherapeutic treatment response
The CRG signature further exhibited a robust prognosis prediction ability, and CRG score was negatively correlated with the survival of UC patients treated with PD-1/L1 blockades (median OS in the IMvigor210 cohort: 7.92 months vs. 10.58 months, p = 0.028, Figure 12A; GSE176307 dataset: 4.57 months vs. 13.00 The functional enrichment analysis underlying the CRG signature. (A) The heatmap exhibited the expression levels of CRG signature-associated genes positively and negatively correlated with CRG score (|Pearson Correlation Coefficient| > 0.4, p < 0.05). The GO enrichment analysis (B) and KEGG pathway enrichment analysis (C) of the identified CRG signature-associated genes. months, p = 0.041, Figure 12B). When exploring the role of the CRG signature predicting the response of immunotherapeutic treatment in the IMvigor210 cohort, no statistically significant difference in CRG score was observed between patients who responded to PD-1/L1 blockades or not (CR/PR vs. SD/PD patients, Figure 12C). Regarding patients with altered FGFR3 in the GSE176307 dataset, the CRG scores were also equivalent between CR/PR and SD/PD patients ( Figure 12D). For patients with wild-type FGFR3 in the GSE176307 dataset, it was notably found that patients who completely/partially responded to immunotherapy have slight lower CRG scores than those with stable/progressive diseases (p = 0.097, Figure 12E).

Discussion
Bladder cancer, according to the latest cancer statistics worldwide, remains to be one of the most prevalent cancers, with approximately 550,000 new cases annually (27). Recently, next-generation sequencing technology has led to the significant advances in the research field of bladder cancer (28), discovering a number of genomic, transcriptomic, and proteomic biomarkers for predicting diagnosis or prognosis of patients, as well as promoting clinically favorable targeted therapeutics and effective immune therapies. Furthermore, the newly developed methods have been widely applied to explore the tumor immune microenvironment in various cancers, which expanded the repertoire of precision medicine, especially immunotherapy (29). In the present study, it was the first time to investigate the integrative transcriptional characterization of cell cycle checkpoint genes in BC progression, and eventually a novel CRG signature and CRG-based nomogram were established, with remarkably robust and stable capacity in prognosis prediction. The clinical association analysis underlying the CRG signature further disclosed that the highrisk group had more patients who were female, diagnosed at older age, and with more aggressive diseases. In addition, the molecular characterization, including functional differentiation analysis, genomic alteration profiling, and TME, between different subgroups promoted the potential strategies of precision treatment for BC patients.
The early retrospective study proposed a risk model, namely, CCP score, based on the expression signature of 31 CCP genes, which was created to predict the aggressiveness of prostate cancer (30). Afterwards, this 31-CCP gene expression signature exhibited significant prognostic values in various cancers, such as ductal carcinoma in situ of the breast (31), clear cell renal cell carcinoma (32), lung cancer (33), and bladder cancer (34). According to the study of bladder cancer, an optimized 12-CCP signature was established, of which the AUC values in predicting patient progression (non-progressor vs. progressor) were 0.70 in the Lindgren cohort and 0.68 in the CNUH cohort. In contrast, it seemed that the CRG signature might outperform the CCP score in prognosis prediction, with The characteristics of genomic alterations between high-risk and low-risk patients from the TCGA-BC cohort. The oncoprint plots exhibited the genomic alteration profile of the high-risk group (A) and low-risk group (B). (C) The genomic alteration enrichment analysis demonstrated the prevalently altered genes between high-risk and low-risk groups. (D) The alteration enrichment analysis of DNA damage response (DDR) genes between high-risk and low-risk groups. (E) Kaplan-Meier curve analysis for BC patients with or without POLE alterations.
the AUC values at 1, 2, 3, 4, and 5 years for OS all beyond 0.75. Furthermore, recent studies have revealed that the regulation of CCP, such as the initiation of replication, the overall rate of replication, and the recovery and resumption of replication forks, was strictly controlled in cancer cells by checkpoints (35-38). In brief, in cancer cells, cell cycle checkpoints could control CCP and further regulate the cancer cell divisions, unless uncontrolled CCP can drive more mutations and genomic instabilities causing the existence of cell cycle and cell apoptosis.
In past decades, a large body of studies focused on the regulation of CCP in tumor progression (39), in which the important role of checkpoints regulating the entire cell cycle processing was often neglected. In the present study, a total of 52 CRGs in cell cycle were found to be significantly correlated with the clinical outcomes of BC patients; meanwhile, the novel CRG signature was constructed as a prognostic model consisting of further selected 23 CRGs. Among these 23 selected CRGs, the expression of PPP2CB has already been proven to be involved in promoting BC cell proliferation and migration (40). The expression of checkpoint CUL4A could mediate the degradation of BECN1 protein to alleviate cell autophagy and enhance the growth of BC (41). The overexpressed PLK2 was detected in BC, and urinary PLK2 protein level was highly correlated with bladder transitional cell carcinoma (42). The suppressed expression of THOC1 could mediate BC cell apoptosis (43). Except for the in-depth insight into the biological functions of the above-mentioned four CRGs, 19 other identified CRGs involved in the signature, namely, The evaluation of tumor environment between high-risk and low-risk groups. (A) Comparison of stromal, immune, and ESTIMATE scores between high-risk and low-risk groups. (B) Comparison of T-cell dysfunction, exclusion, and TIDE scores between high-risk and low-risk groups. (C) The evaluation of 22 tumor-infiltrated lymphocytes between high-risk and low-risk groups. *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001, ns represented "not significant".
ANAPC4, MMAB, B9D2, NDEL1, SLC25A15, PSMA7, PSMB10, PSMB5, PSMB8, RAD9A, REXO2, ARID3A, CHMP4C, DDX39B, FBXO31, FBXO6, PRPF19, RGCC, and TIPIN, were first identified to play a potential role in BC progression, and their detailed functions remained enigmatic, which merited further experimental verification. Underlying the CRG signature, correlation analysis revealed that signature-associated genes were mainly responsible for the collagen-containing ECM, and ECM-related and MAPK signaling pathways, indicating that the ECM-specific heterogeneity and related signaling pathways between different risk groups mainly caused the significant difference in OS. It has been revealed that COL1A1 protein was markedly upregulated in MIBC, which could activate the epithelial-mesenchymal transition and TGF-b signaling pathway contributing to the proliferation and invasion of bladder cancer cells (44), whereas the significantly downregulated expression of collagen type IV-a1 and a2 (COL6A1 and COL6A2) was found to promote tumor progression in both the NMIBC and MIBC tissue samples (45). Moreover, the tumor-related macrophages could secrete the type I collagen via activating the PI3K/AKT signaling pathway to stimulate the development of bladder cancer; meanwhile, the number of macrophages and the expression of M2 macrophage-associated genes (ARG-1, IL-10, and TGF-b) were remarkably elevated in malignant bladder tumor tissue samples (46). Furthermore, TGF-b could prompt the tumor immune escape and the resistance of immune therapy (47), which has been validated in the IMvigor210 cohort in which BC patients not responding to the therapeutic treatment of anti-PD-L1 agent (Atezolizumab) were highly correlated with the TGF-b signature in fibroblasts and commonly had fibroblastand collagen-rich peritumoral stroma (48). In addition, a singlecell proteomic analysis further revealed the chaotic tumorassociated collagens in the TME of MIBC (49). Thus, therapeutic treatment targeting collagen-specific ECM and related signaling pathways could be a potential regimen for BC patients.
Genomic alteration analysis revealed that the high-risk group had higher prevalence of tumor suppressor genes RB1, FBXW7, and NFE2L2, consistent with the previous findings that the altered RB1 (50), FBXW7 (51), and NFE2L2 (52) were correlated with tumor progression and worse outcomes of BC patients. However, the potential contribution of frequently altered ASAP2 and PCSK6 in the progression of bladder cancer remains unclear, needing to be further investigated. As known, the FGFR3 alterations were mainly enriched in the luminal bladder cancer correlated with better prognosis; specifically, the low-risk group in the present study was also observed to have more patients with FGFR3 alterations. More impressively, it was uncovered that there was a significant difference in alteration frequency of only one DDR signaling pathway-related gene-POLE-between high-risk and low-risk The estimate of chemotherapy treatment response via the GDSC database. The violin plots exhibited the IC 50 value representing the response of chemotherapeutic treatment by cisplatin (A), docetaxel (B), paclitaxel (C), methotrexate (D), and doxorubicin (E) between four different molecular subgroups. *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001, ns represented "not significant".

A B
D E C FIGURE 12 The predictive role of the CRG signature for the response of immunotherapy treatment in the IMvigor210 cohort and the GSE176307 dataset. Kaplan-Meier curve analysis to compare overall survival of patients between high-risk and low-risk groups from the IMvigor210 cohort (A) and the GSE176307 dataset (B). Comparison of the CRG score between CR/PR and SD/PD groups in the IMvigor210 cohort (C) and in the groups of patients with altered FGFR3 (D) and with wild-type FGFR3 (E) from the GSE176307 dataset. groups, and Kaplan-Meier analysis revealed that MIBC patients with altered POLE exhibited a trend of better OS. Owing to the limited number of MIBC patients presenting POLE alterations in the present study, this result should be highlighted and verified in the future.
Taken together, the difference in TME, including tumor purity, T-cell viability, and the proportion of TILs, between high-risk and low-risk groups was further investigated. Macrophage M0 and M2 were consistently enriched in high-risk groups, which also had a lower tumor purity expressed by the higher ESTIMATE score and a higher level of exhausted T cells indicated by the higher TIDE score. Similarly, it has been reported that tumor-associated macrophages could restrain the infiltration level of T cells in tumor tissues (53). In addition, altered FGFR3 is highly enriched in luminal or its papillary subtypes, which have been characterized by immune-suppressive states (50). Controversially, a retrospective analysis revealed that there was no difference in immune checkpoint blockade response between patients with altered FGFR3 and those with no altered FGFR3 (54), and a recent study found the equivalent T-cell receptor diversity between these patients (26). Whether the status of FGFR3 could influence the treatment response by the immune checkpoint blockade is discussed heatedly. In the present study, it was found that the discrepancy of TME between high-risk and low-risk groups among BC patients with wild-type FGFR3 was a little small, except that low-risk patients had more activated T cells but a lower proportion of tumor-associated macrophages (especially for macrophage M0 and M2), while compared to the BC patients with wild-type FGFR3, high-risk patients with altered FGFR3 had nearly similar predicted infiltrated immune cell proportions and immune responses, but only the expression levels of TIM3 and CD47 were equal, such that these BC patients seemed to be more sensitive to immune checkpoint inhibitors targeting TIM3 and CD47, which needed further analysis. Additionally, the low-risk BC patients with altered FGFR3 exhibited the highest tumor purity and the most activated T cells, which were associated with better treatment responses of immunotherapies (23). However, in the present study, the expression levels of well-known immune checkpoint genes were found to be extremely low. Thus, in clinic, it was highly recommended that immunohistochemistry staining of immune checkpoint(s) should be detected first during immunotherapeutic treatment.
In addition, the different BC molecular subtypes were identified to have distinctive sensitivity to the chemotherapeutic drugs, such as cisplatin, docetaxel, paclitaxel, methotrexate, and doxorubicin, which were commonly used for BC patients (28). The CRG signature further demonstrated its predictive ability in therapeutic treatment response. Currently, cisplatin has been regarded as a typical chemotherapeutic strategy for UC patients, and it mainly functioned to induce cell cycle arrest and cell apoptosis (55, 56). Furthermore, a previous study has revealed that docetaxel could make antiproliferative and apoptotic effects on bladder cancer cells (57). In the present study, it was found that patients with a higher CRG score, regardless of FGFR3 status, seem to be more sensitive to cisplatin or docetaxel. Although clinical benefits of paclitaxel were limited owing to patients' resistance, the failure of first-line combination treatment of cisplatin and gemcitabine for advanced and/or metastatic UC patients provided an opportunity for paclitaxel because of its apoptotic effects (58). Through the evaluation analysis of chemotherapy treatment response, it was recommended that low-risk BC patients with wild-type FGFR3 might benefit more from paclitaxel. Additionally, the CRG signature also exerted its influence on the patients' sensitivity to methotrexate, especially suitable for low-risk BC patients with wildtype FGFR3; however, its applicability for BC patients in clinic was rare, which needed more verification. Doxorubicin has been proven to disturb the regulation of cell cycle and induce cell death, usually as a vital constituent of combination treatment for MIBC patients (59). Indeed, the alteration frequency of FGFR3 in MIBC patients was relatively lower, compared with NMIBC patients (60,61). In the present study, we further found that high-risk BC patients with wild-type FGFR3 were likely to have higher sensitivity to doxorubicin. Collectively, the checkpoint mechanisms in the regulation of cell cycle greatly influenced the response of chemotherapeutic treatment.
The CRG signature also demonstrated its ability in predicting the response of immunotherapeutic treatment. Of note, the CRG signature exhibited better response prediction capacity for BC patients with wild-type FGFR3, of whom those with lower CRG score might be more sensitive to the immunotherapeutic treatment with anti-PD-1 or anti-PD-L1 drugs. Regarding patients with altered FGFR3, the CRG signature could not differentiate the potential CR/PR and SD/PD patients. As described before, the alteration frequency of FGFR3 in lowgrade and/or early-stage bladder tumors was evidently higher comparatively, such that the CRG signature could perform better in the response prediction of immunotherapeutic treatment for higher-grade BC patients, especially for those with wild-type FGFR3, while only 17 BC patients with altered FGFR3 included in the GSE176307 dataset received immune therapy treatment hence, the ability of the CRG signature to predict the immunotherapeutic treatment response for BC patients with FGFR3 alterations remained to be explored.

Conclusions
The present study comprehensively delineated the integrative transcriptional profiling of cell cycle checkpoint genes in the BC progression. The novel prognostic CRG signature and nomogram exhibited good performance in prognosis prediction for BC patients; furthermore, the molecular characterization underlying the CRG signature provided deep insights into key risk factors leading to the aggressive bladder cancer, promoting the development of precision medicine. The predictive role of the CRG signature in treatment response offered potential precision treatment strategies for different BC molecular subtypes.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors.

Ethics statement
The studies involving human participants were reviewed and approved by The ethics committee of the Chinese PLA General Hospital (S2019-302-01). The patients/participants provided their written informed consent to participate in this study.

Author contributions
W-WS, J-ZG, YH, and BY proposed and designed this study. W-WS and J-ZG collected raw data. Y-PL, QS, QX, and B-YQ pre-treated raw data and performed the statistical analyses. W-WS, J-ZG, Y-PL, Z-QM, and BY wrote the original manuscript. YH and BY revised the manuscript. All authors have read and approved this study to be published.

Funding
This study was supported by the Department of Medical Oncology, Senior Department of Oncology, The Fifth Medical Center of PLA General Hospital.