Multi-omics comprehensive analysis reveals the predictive value of N6-methyladenosine- related genes in prognosis and immune escape of bladder cancer

2.1Data sources

In this work, transcriptome data, single nucleotide polymorphism (SNP) data, copy number variation (CNV) data and related clinical information were downloaded from the TCGA database (https://tcga-data. nci.nih.gov/tcga/) for 430 cases (411 BC cases and 19 normal cases). Firstly, raw data of each case was extracted and generated a matrix by using Perl software. Then, case IDs were conformed into a consistent layout, as well as annotated gene names accurately in the matrix. Next, a standardized transcriptomic matrix was generated by applying formula log2 (x + 1) for further investigation. The clinical information (Supplementary Table 1) was first examined to obliterate any duplicates or missing data, then matched with transcriptomic data of each case for further analysis. As for the SNPs data, it could be directly used after matching with the patient ID of transcriptome data. Meanwhile, transcriptome data and clinical information (Supplementary Table 2) of 188 BC cases were obtained from the GSE13507 dataset (https://www.ncbi. nlm.nih.gov/geo/query/acc.cgi?acc=GSE13507). The transcriptome data was annotated by using platform file GPL6102, deleting missing data and standardized these data by applying formula log2 (x + 1) for further investigation. The clinical information was examined to delete duplicating and missing data for further analysis.

2.2Identification of differentially expressed genes

Before analysis, we standardized all transcriptomics data by applying log2 (x + 1). Then, we applied the edge R package with FDR < 0.05 and |log2FC|⩾ 1 to identify differentially expressed genes (DEGs) of transcriptomics in both the TCGA and GEO cohorts.

2.3Protein-protein interaction network

We generated a protein-protein interaction network (PPI) for MRGs by using STRING (https://string-db. org/) and Cytoscape software.

2.4Enrichment analysis

We applied Gene Set Variation Analysis (GSVA), Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) for functional or pathway enrichment analysis. Briefly, expression profiles of DEGs or MRGs were put into R software to perform the enrichment analysis, calculating enrichment scores, estimating significance levels and performing multiple hypothesis, which finally indicated the potential gene functions and pathways.

2.5M6A scores were evaluated for patients

Principal components analysis (PCA) is an effective technique for reducing data dimensionality. Calculation steps (assuming that the sample data is n-dimensional, with a total of m): 1. Organize these sample data into a matrix Xnm by column; 2. Row-wise averaging X, which means first calculating the average value of each row, and then subtracting this average value from each element in that row; 3. Calculate the covariance matrix; 4. Find the eigenvalues of the covariance matrix, which are the corresponding eigenvectors; 5. Arrange the feature vectors into a matrix by row from top to bottom according to the corresponding feature values, and take the first k rows to form a matrix P; 6. Then Y = PX is the data set after dimensionality reduction to k. These steps are ultimately accomplished with the help of R package. In this work, we first used PCA to analyze the expression of 11 survival-related genes (CLIC4, KPNA2, PAQR4, MTHFD2, CDC25B, WDR62, TPX2, CRTAC1, KLHL3, SMAD6, ZNF552) and obtained the corresponding PCA score. Next, we calculated the m6A score of each patient according to the PCA score: m6A score = pca (1) + pca (2) + …+ pca (n). Finally, patients were divided into a high score (HS) group and a low score (LS) group according to the m6A score in both the TCGA and GEO cohorts for subsequent analysis.

2.6Cell culture

In this study, human bladder cancer cell lines UMUC3 and RT112 were cultured in RPMI-1640 (Gibco, USA) supplemented with 10% fetal bovine serum and 1% penicillin/streptomycin, as well as cultured at 37^∘C and in 5% CO2.

2.7Transfection assay

We obtained the overexpression plasmid of IGFBP1 and paired empty plasmid from GeneCopoeia (Guangzhou, China). Firstly, we planted BC cells into 6-well plates and waited for cell density reached at 80%–90%. Then, we changed the culture medium into 1.5 ml Opti-MEM (Gibco) for each well at two hours before transfection. Meanwhile, we set two groups for transfection: negative control group (NC) and knockdown group (SH). For each well in SH group, 10 ul IGFBP1 knockdown plasmid was added into 240 ul Opti-MEM, as well as 5 ul lipo2000 (Thermo, USA) was supplanted into 245 ul Opti-MEM and set still at room temperature for 5 minutes. Next, we mixture plasmid solution with lipo2000 solution and set still at room temperature for 10 minutes, following adding total 500 ul transfection solution into 1.5 ml Opti-MEM culture medium to form 2 ml volume transfection system in each well. The NC group followed the same protocol but changed the knockdown plasmid into empty plasmid. After 8 hours transfection, we changed the transfection culture medium into 2 ml RPMI-1640 supplemented with 10% serum. The transfected efficacy was confirmed by examining the fluorescence percentage. We chose the wells that had a minimum of 80% fluorescence percentage for further investigation.

2.8RT-PCR

The reactions were performed on the QuantStudio 5 Real-Time PCR instrument (Thermo Fisher, USA) using 2x SYBR Green Pro Taq HS Premix II (AgBio, China). Normalization of mRNA levels was against GAPDH levels. For calculating the relative gene expression, the 2-Δ⁢ΔCT approach was employed. (First, for all test samples and calibration samples, the CT value of the internal reference gene is normalized to the CT value of the target gene: ΔCT test = CT target, test – CT ref, test; ΔCT calibrator = CT target, calibrator CT ref, calibrator. Second, normalize the ΔCT value of the test sample with the ΔCT value of the calibration sample: Δ⁢ΔCT =ΔCT test – ΔCT calibrator. Finally, calculate the expression level ratio: 2-Δ⁢ΔCT = the ratio of the amount of expression.)

2.9CCK8 assay

This assay was followed by the protocol of Beyotime CCK8 kit. Briefly, we adjusted the cell density of transfected and control cells to 5000/100 ul. Then, we added 100 ul cell suspension into each well of 96-well plate, supplemented with 10 ul CCK8 solution, and cultured at 37^∘C for one hour. Afterward, we examined the optical density at 450 nm wavelength to examine the cell proliferation. This assay was repeated at least three times.

2.10Wound healing assay

We seeded the transfected cells into six-well plates and waited for the cell density reaching more than 90%. Then we scratched several vertical lines by using a 200 ul pipette tip. After three times washing with PBS solution, we changed the culture medium to 2 ml serum-free culture medium. Photos of this assay at 0 hour, 24 hour were taken to examine the wound healing rates.

2.11Transwell assay

The transfected cells were incubated in serum-free medium for 12 h and then adjusted to a density of 1*10⁶ cells/ml. Next, 200 mL of cell suspension was added to the transwell chambers (Corning, NO.3422), which were fit into the wells of 24-well plates. The wells of the plates also each contained 500 ml of RPMI-1640 medium supplemented with 10% FBS. After 24 hours of incubation, cells on the upper membrane were removed with cotton wool, whereas cells adhering to the lower surface were fixed in methanol for 30 min and then stained with 0.1% crystal violet for 10 minutes. After natural air drying, the migrating cells on the lower surface of the membrane were counted under an optical microscope at 200 folds magnification.

2.12Statistical analysis

We employed one-way ANOVA and t tests for comparisons between groups, and the comparison of two or more constituent ratios was performed by the chi-square test. Correlation analysis, heatmap, waterfall curve and box plots were completed by R software (version 3.5.1). We processed all statistical analyses by applying SPSS 19.0 software (SPSS. Inc., Chicago, IL, USA) or R software. In this work, we set the significance level to P< 0.05.

Figure 1.

m6A MRGs significantly affect the survival of bladder cancer patients. (A, B) Overall survival (OS) and disease-specific survival (DSS) prognostic networks were constructed according to 17 MRGs. (C–I) Kaplan-Meier analysis of OS for m6A MRG genes.

Figure 2.

m6A clusters were identified based on MRGs in BC patients. (A) BC samples (n= 596) were divided into four m6A clusters based on the consensus clustering matrix (k= 4) in both TCGA and GEO cohorts. (B) Heatmap of the expression status of 17 m6A MRGs in m6A clusters. (C) The differences in 23 immune cell infiltration levels in m6A clusters. (D, E) KEGG and GO functional analysis between m6A clusters A and B. (F, G) The results of GO and KEGG enrichment analysis for the intersection of differentially expressed genes between m6A clusters.

3.1Multi-omics biological characteristics of m6A MRGs and its relationship with survival in BC

First, we investigated the CNV status of 23 m6A MRGs in TCGA patients with BC. The results showed that 56.52% (13/23) of m6A-related genes gained more CNV frequencies than they lost, and VIRMA had the highest frequency. In contrast, 43.48% (10/23) of genes lost more CNV frequencies than they gained, and YTHDF2 had the top percentage (Supplementary Fig. 1A). Moreover, the localization of these genes in the chromosome was described as well (Supplementary Fig. 1B). Then, we found a good protein-protein interaction network, which revealed a potential regulatory relationship among these genes (Supplementary Fig. 1C). Next, we compared the expression statuses of these genes from 411 BC cases with 19 normal cases and found that HNRNPA2B1 and HNRNPC had the most significant differences (Supplementary Fig. 1D). After that, we found that these 23 genes were mutated in 24.51% of 411 TCGA BC patients, of which METTL3 was the most frequently mutated (Supplementary Fig. 1E). Furthermore, we investigated the expression status of 22 other genes when the METTL3 gene was mutated or not. The results showed that 17 genes were differentially expressed in METTL3 mutant or wild-type cells, and the results of RBM15 and YHDF2 are shown as examples (Supplementary Fig. 1F and G), which indicated that METTL3 was one of the key m6A MRGs in BC.

We combined 411 BC cases in TCGA and 188 GEO BC cases and finally identified 17 MRGs with complete data for subsequent analysis. Next, we constructed overall survival (OS) and disease-specific survival (DSS) prognostic networks according to these 17 MRGs (Fig. 1A and B). More precisely, the networks exhibited the detailed function of these MRGs in m6A methylation (erasers, readers or writers), as well as the correlation between genes and the role of the gene as a protective or risk factor in OS and DSS. Then, we investigated every single gene by Kaplan-Meier analysis and finally found that ALKBH5, ZC3H13, LRPPRC, IGFBP1, and RBM15B were associated with poor OS and that FMR1, YTHDF3, YTHDC1, RBMX and YTHDC2 were associated with better OS (Fig. 1C–L). Additionally, we performed Kapan-Meier analysis for DSS and found that 9 of 17 MRGs were statistically significant (Supplementary Fig. 2A–I). All these findings highlight the important roles of m6A MRGs in the prognosis of BC.

3.2m6A clusters were identified based on MRGs in BC patients

Based on the gene expression of these 17 m6A MRGs, a consensus clustering analysis was performed for 596 BC cases (408 from TCGA and 188 from GEO), and the results showed that our patients were well divided into 4 m6A clusters when k= 4 (Fig. 2A). After that, we found that these MRGs showed apparent differences among the four groups of patients (Fig. 2B). Then, a boxplot for immune infiltration status was generated based on the m6A clusters (Fig. 2C). According to the results, we noticed significant differences in immune cell infiltration among the 4 m6A clusters of patients. More specifically, m6A cluster A had a high level of immune cell infiltration (such as CD4+ T cells, CD8+ T cells, NK cells and B cells), suggesting good immune function. In contrast, m6A cluster D had a poor level of immune infiltration, suggesting poor immune function. After that, ssGSEA was performed to explore the underlying functions or pathways for these clusters. KEGG analysis showed that metabolism-related pathways, such as linoleic acid metabolism, were significantly downregulated in cluster A compared with cluster B, while tumor-related pathways, such as DNA replication pathways and cell cycle pathways, were apparently upregulated in cluster A (Fig. 2D). In addition, similar results were observed in KEGG analysis of comparisons between cluster A and other m6A clusters (Supplementary Fig. 3A and B). Accordingly, these enriched pathways could partly explain the reason for the worst prognosis and lowest m6A score of patients in cluster A based on the results from the subsequent analysis of our work. In addition, we conducted pairwise comparisons of different KEGG pathways between m6A clusters (Supplementary Fig. 3C–E).

Consistently, we also performed GO functional analysis between m6A clusters. From the results of the comparison between clusters A and B, we noticed that several GO biological processes (GOBP), such as DNA replication initiation, and GO cellular components (GOCC), such as replisome, were at higher levels in cluster A, while some GOBPs and GO molecular functions (GOMF) were downregulated in cluster A. Interestingly, we noticed that insulin-like growth factor (IGF)-related GOBP levels were lower in patients from cluster A, whose prognosis and m6A score were the worst, while IGF levels were commonly upregulated in advanced malignancies (Fig. 2E) [17, 18]. Nevertheless, comparisons of GO function analysis between cluster A and other m6A clusters were also performed (Supplementary Fig. 4A and B). In addition, we also conducted pairwise comparisons of GO functions between other m6A clusters (Supplementary Fig. 4C–E). After that, we distinguished the differentially expressed genes by comparing every two m6A clusters (clusters A and B, A and C, A and D, B and C, B and D, and C and D) and further intersected these differentially expressed genes among these comparisons to process Unicox analysis, which finally identified 11 OS-related genes (OSRGs) (CRTAC1, SMAD6, CDC25B, ZNF552, KLHL3, TPX2, WDR62, MTHFD2, PAQR4, KPNA2 and CLIC4). Additionally, we employed GO and KEGG enrichment analysis for these 11 OSRGs. We found that these genes were mainly enriched in microtubule cytoskeleton organization involved in mitosis, regulation of cellular response to growth factor stimulus and epithelial tube morphogenesis (Fig. 2F). At the same time, we noticed that OSRGs were mainly enriched in microRNAs in cancer, as well as TGF-β signaling pathways and cell cycle pathways (Fig. 2G), which revealed that these OSRGs were involved in tumorigenesis-related pathways and could apparently alter prognosis. Therefore, these findings reveal the important value of m6A MRGs in identifying immune cell infiltration and prognosis in BC patients and highlight OSRGs as tumor-related genes.

3.3m6A score was established to predict the prognosis of BC patients

Figure 3.

The m6A score was established to predict the prognosis of BC patients. (A) BC samples (n= 596) were separated into two gene clusters based on the consensus clustering matrix (k=2) in both TCGA and GEO cohorts. (B) Heatmap of 11 survival-related gene expression statuses between gene clusters. (C) Boxplot of 11 survival-related gene expression levels between gene clusters. (D, E) Kaplan-Meier analysis for OS (D) and DSS (E) between gene clusters. (F, G) Kaplan-Meier analysis of OS (F) and DSS (G) between the HS group and LS group. (H–J) Living status (H), T grade (I) and N grade (J) in the m6A score groups.

To further investigate the roles of these 11 OSRGs in our patients, a new consensus clustering analysis was performed. When k= 2, cases could be well separated so that these 596 BC cases were divided into 2 gene clusters (Fig. 3A). As shown in Fig. 3B, these 11 OSRGs were apparently different in gene expression between the 2 gene clusters, and the details are shown in Fig. 3C. To be more precise, we found that the expression levels of CRTAC1, SMAD6, ZNF552 and KLHL3 were higher in gene cluster B than in gene cluster A, while the others showed the opposite results. Then, we followed Kapan-Meier analysis for OS and DSS separately based on the gene clusters (Fig. 3D and E). The results showed that patients in gene cluster B had a better prognosis than patients in gene cluster A, which demonstrated the ability of these genes to predict the prognosis of BC patients. Furthermore, we calculated the m6A score of each case via principal component analysis (PCA) and divided them into a high-score (HS) group and a low-score (LS) group. Next, we ran Kapan-Meier analysis for OS and DSS based on m6A score groups separately (Fig. 3F and G) and found that patients in the HS group had better survival possibilities than those in the LS group. Surprisingly, we observed the same survival trends among groups with consistent clinical parameters (T grades, N grades, age and sex, etc.) (Supplementary Fig. 5A–H). In addition, we further investigated the clinical characteristics of these BCs between m6A score groups. As the results showed, patients in the HS group had less mortality than those in the LS group, and patients in the HS group showed lower T stage and N stage than those in the LS group (Fig. 3H–J). These findings demonstrated the ability of the m6A score to predict the prognosis of BC patients.

3.4m6A score is closely related to TMB and gene mutations in BC patients

Figure 4.

The m6A score is closely related to TMB and gene mutations in BC patients. (A) m6A score levels in m6A cluster groups. (B) m6A score levels in gene clusters. (C) Corresponding relations between the m6A cluster, gene cluster, m6A score and Fustat in our patients. (D) TMB levels in the m6A score groups. (E) The correlations between TMB levels and m6A score. (F) Kaplan-Meier analysis of OS in patients with high and low TMB levels. (G) Kaplan-Meier analysis for OS in cases with four combinations that consisted of different TMB levels and different m6A scores. (H, I) The top 20 mutant genes in the HS (H) and LS groups (I).

First, we examined the relationship between the m6A score and m6A clusters and gene clusters (Fig. 4A and B). From the results, we found that patients in m6A cluster B, m6A cluster D and gene cluster B had higher m6A scores, suggesting a better prognosis, which is basically consistent with our previous conclusion. Next, a more intuitive Sankey diagram showed the corresponding relationship among the m6A cluster, gene cluster, m6A score and fustat status among 596 BC cases (Fig. 4C). Next, we found that the m6A score was negatively correlated with the TMB level, and the HS group was associated with a lower TMB level than the LS group (Fig. 4D and E). Then, we investigated the prognosis of BC with high and low TMB levels via Kaplan-Meier analysis for OS. Surprisingly, the results showed that patients with higher TMB levels gained better survival possibilities than those with lower TMB levels (Fig. 4F). Moreover, patients with high TBM and high m6A scores had the best prognosis, while patients with low TMB and low m6A scores had the worst prognosis (Fig. 4G). Nevertheless, we know that the effect of TMB on cancer patients is double-sided. In our results, high m6A score and high TMB both showed better prognosis, while the inverse association between TMB and m6A score seems puzzling, but it is not contradictory since the overall outcome is usually caused by multiple biological factors. A single factor of TMB is not sufficient to alter the final outcome in our patients.

Additionally, we assessed the gene mutation status of 406 cases in the m6A score groups via Maftools (Fig. 4H and I). From the results we noticed that the total mutant frequencies between HS group and LS group were similar (93.33% vs 94.02%), but FGFR3 was the most common mutant gene in cases from HS group, while TP53 in LS group. The differences in gene mutation status between m6A score groups could partly explain the different outcomes of cases in different groups.

3.5The LS group predicts a higher risk of immune escape in BC patients

Figure 5.

The LS group predicts a higher risk of immune escape in BC patients. (A, B) The immune score of Treg cells (A) and CD8+ T (B) cells in the m6A score groups and correlations between immune score and m6A score. (C–E) The expression of PD-L1, PD-1 and CTLA4 between the HS and LS groups. (F–I) The immune dysfunction score (F), immune exclusion score (G), PD-1 and PD-L1 immunotherapy effect score (H), and negative PD-1 and PD-L1 immunotherapy effect score (I) between the HS and LS groups. (J) The MSI status between the HS and LS groups.

Furthermore, we assessed the correlation of the m6A score with tumor immunity in BC. According to the results, we found that the m6A score was closely related to immune cell infiltration, and the LS group had higher infiltration of Treg cells and lower infiltration of CD8+ T cells than the HS group (Fig. 5A and B). Meanwhile, iDC, monocyte, Tcm and NK CD56bright cell infiltration was lower in the LS group than in the HS group and was positively correlated with the m6A score (Supplementary Fig. 6A–D). These findings suggested that patients in the LS group might have poorer immune function. Then, we examined the levels of the immune checkpoints PDL1, PD1 and CTLA4 in the m6A score groups. The results showed that these checkpoints were upregulated in the LS group compared to the HS group, which indicated that cases in the LS group had a higher risk of immune escape (Fig. 5C–E). In addition, a higher immune dysfunction score and immune exclusion score were identified in the LS group, which further demonstrated that the LS group was associated with a higher risk of immune escape in BC (Fig. 5F and G). Surprisingly, we found that when the patients were treated with PD-1 and PD-L1 immunotherapy, the LS group achieved a better therapeutic effect than the HS group. In contrast, when patients were treated with negative PD-1 and PD-L1 immunotherapy, the HS group achieved a better therapeutic effect than the LS group (Fig. 5H and I). In addition, we found that cases in the HS group had higher microsatellite instability (MSI), which might suggest a possible method for the gene detection of BC (Fig. 5J). These results suggested that the LS group had a higher risk of immune escape than the HS group and might gain more benefit from ICIs.

Figure 6.

Experimental verification of m6A related gene and prognosis of bladder cancer. A: The gene expression of ALKBH5, ZC3H13, LRPPRC, IGFBP1, YTHDC1 and YTHDC2 in bladder cancer tissues and adjacent normal tissues was measured by qPCR. B: CCK8 was used to determine the growth activity in the IGFBP1 knockdown group and empty vector group of bladder cancer cells. C and D: Gene expression levels of IGFBP1 in different groups. E and F: A wound healing assay was used to measure the migration ability of UMUC3 and RT112 cells in the IGFBP1 knockdown group and empty vector group. G and H Transwell assays were used to measure the migration ability of UMUC3 and RT112 cells in the IGFBP1 knockdown group and empty vector group.

3.6Experimental verification of m6A related gene and prognosis of bladder cancer

In this work, we identified a total of six m6A related genes that were both associated with overall survival and disease-specific survival, of which ALKBH5, ZC3H13, LRPPRC, IGFBP1 were associated with poor survival and YTHDC1 and YTHDC2 were associated with better survival. Furthermore, qPCR showed that the gene expression of ALKBH5, ZC3H13, LRPPRC, IGFBP1 in bladder cancer tissues was significantly higher than that in normal tissues, and YTHDC1 and YTHDC2 in bladder cancer tissues was significantly lower than that in normal tissues (Fig. 6A).

Among the six m6A related genes, of which ALKBH5, LRPPRC and YTHDC1 have been clearly reported to be associated with bladder cancer in the literature, but the biological roles of ZC3H13, IGFBP1 and YTHDC2 in bladder cancer are poorly understood [18, 19, 20, 21, 22, 23, 24, 25, 26]. Furthermore, we analyzed ZC3H13, IGFBP1 and YTHDC2 in TCGA dataset, and found that IGFBP1 had a higher expression level in non primary bladder cancer and high lymph node metastasis patients (Fig. 6B and C), while ZC3H13 and YTHDC2 did not observe this trend. These findings suggest that IGFBP1 may be more related to the poor prognosis of bladder cancer. Therefore, we selected IGFBP1 for further experimental verification. The results showed that the growth activity of bladder cancer cells in the IGFBP1 knockdown group was significantly inhibited in UMUC3 and RT112 cells compared to that in the empty vector transfection group (Fig. 6E and F). In addition, in UMUC3 and RT112 cells, the metastatic ability of the IGFBP1 knockdown group was significantly inhibited compared with that of the empty vector transfection group (Fig. 6G and H). These results suggested the important role of m6A-related genes, especially IGFBP1, in the prognosis of bladder cancer.