Identifying robust DEGs via the RRA method
Using the selection criteria, 4 independently eligible BLCA datasets were enrolled for subsequent RRA analysis. A series of clinical traits, including GEO accession ID, Platform ID, as well as the number of genes for each platform are displayed in Table S Ⅰ(5-8). Based on RRA analysis data, we identified 111 up-regulated and 232 down-regulated remarkable DEGs (Supplementary file 1). Besides, the top 50 up-regulated, as well as down-regulated DEGs are depicted in the heatmap (Figure 1).
Functional enrichment analysis of DEGs
The biologically functioning DEGs were revealed via the GO and KEGG functional enrichment analysis using DAVID. We considered the results significant only if p<0.05, we have highlighted the three categories of the GO results in Figure 2A and Figure 2B. Results on the upregulated and downregulated DEGs in top 15 findings derived from the GO enrichment analysis are depicted in Table S Ⅱ and Table S Ⅲ. Of note, the upregulated genes were highly enriched in protein binding (ontology: MF), nuclear division during mitosis (ontology: BP), and cytoplasm (ontology: CC). Besides, the downregulated genes were highly abundant in, extracellular exosome (ontology: CC), and binding of calcium ions (ontology: MF) and cell adhesion (ontology: BP). As to KEGG pathway analysis, ECM-receptor interaction, Focal adhesion, P13K-Akt signaling cascade, Proteoglycans in cancer, as well as Vascular smooth muscle contraction, were mostly associated with these genes (Figure 2C).
WCGNA analysis and modules significance calculation
To reveal the key modules highly related to the clinical characteristics of BC, we analyzed the WGCNA on the TCGA-BLCA cohort by integrating the DEGs retrieved from the RRA analysis (Figure 3). Clinical information of BC sample from TCGA, including stage, age, grade, and TNM classification were retrieved (Figure 3A). We set the soft-thresholding power at 6 (scale free R2=0.9) and cut height as 0.25. Consequently, 4 modules were identified (Figure 3B-3D). Based on the heatmap showing module-trait correlations, the blue module shows the highest correlation with clinical symptoms (Figure 3E), particularly the stage (correlation coefficient=0.24, p=1E-06,). The blue module had 67 genes (see Supplementary file 2). We set the module membership (MM)>0.8 and gene significance (GS)>0.2 then identified 19 hub genes from the blue module: EDNRA, SERPINF1, COLEC12, FBLN5, DDR2, SFRP2, OLFML3, AEBP1, DCN, CDH11, TIMP2, LUM, DPT, COL6A3, COL16A1, EMILINN1, SPON1, OLFML1 and CRISPLD2. Through GO and KEGG analyses, we uncovered the prospective biological roles of the genes in the blue module. The highest remarkable GO terms for biological process, molecular function, and cellular component, as well as KEGG pathways, are depicted in Figure 4A-4D. Following this evaluation, genes within the blue modules were primarily linked to signal transduction, cell adhesion, and extracellular matrix organization.
Survival analysis and significant gene identification
We assessed whether the 19 hub genes in BC were clinically relevant. To achieve this, correlation assessment of the hub genes with prognosis outcome of BC patients in TCGA-BLCA data sets was performed. By optimizing the cut-off values for hub gene analysis, CDH11, COL6A3, EDNRA and SERPINF1 were highly expressed and were associated with poor prognosis (Figure 5A and Figure S2). Furthermore, receiver operating characteristics (ROC) curves demonstrated that they had high diagnostic potential as BC biosignatures (Figure S3, CDH11 AUC: 0.699, COL6A3 AUC: 0.697, EDNRA AUC: 0.833, SERPINF1 AUC: 0.804), suggesting the potential use of the genes as indicators in monitoring prognosis.
Differential expression of CDH11, COL6A3, EDNRA and SERPINF1
We compared the mRNA expression of CDH11, COL6A3, EDNRA and SERPINF1 between bladder tumor and neighboring healthy tissues, respectively. This was based on data for RNA-sequence obtained from the Oncomine and TCGA databases. Notably, the transcriptional levels of CDH11, COL6A3, EDNRA and SERPINF1 expressions were lowly expressed in BC tissues in comparison to healthy tissues (Figure 5B). Besides, there was a significant correlation of CDH11 mRNA expression and BC samples with a mild clinical stage (Figure 5C), whereas the lowest CDH11 mRNA expression was reported stage Ⅰ +Ⅱ. Similarly, we evaluated the association of CDH11 mRNA expression with different pathological grade, whereby it was revealed that mRNA expression of CDH11 is significantly correlated with lower pathological grades (Figure 5D). Additionally, mRNA levels of COL6A3, EDNRA and SERPINF1 were lower in BC tissues (Figure 5B). COL6A3, EDNRA and SERPINF1 mRNA expression in BLCA sample were significantly correlated with mild clinical staging, whereas the lowest COL6A3, EDNRA and SERPINF1 mRNA expression were detected in stage Ⅰ +Ⅱ (Figure 5C). Moreover, mRNA expression levels of COL6A3, EDNRA and SERPINF1 were related to lower clinicopathological grading (Figure 5D). Collectively, we demonstrated that the expressions of CDH11, COL6A3, EDNRA and SERPINF1 were lower in BC tissues compared to healthy tissues. Thus, the hub gene CDH11, COL6A3, EDNRA and SERPINF1 could play a pivotal role in bladder cancer progression. Overall, low expression of CDH11, COL6A3, EDNRA and SERPINF1 mRNA is significantly associated with mild clinical-pathological parameters in BC patients and is significantly lowered in the early disease stages. This may be vital in the early BC diagnosis.
Association of hub genes’ expression with tumor-infiltrating immune cells
Referring to the critical roles of invading immune cells within the tumor microenvironment, we comprehensively analyzed immune signatures plus immune infiltrates. From the TIMER web resource, the association between CDH11, COL6A3, EDNRA and SERPINF1 immune signatures and tumor purity or numerous vital immune cells was revealed. CDH11, COL6A3, EDNRA and SERPINF1 were all negatively correlated with tumor purity. The correlations (Cor>0.5 and p<0.05) were considered to be the strongest correlated. Although it was observed no or weak correlations of these genes with infiltration of CD8+ T cells, dendritic cells, CD4+ T cells, B cells, and neutrophils, CDH11, COL6A3 and SERPINF1 were significantly associated with macrophages. (Figure 6)
To assess the potential roles of CDH11, COL6A3, EDNRA and SERPINF1 in BC, GSEA was conducted for hallmark analysis of the genes on the TCGA-BLCA RNA-seq data. Genes in low expression CDH11, COL6A3, EDNRA and SERPINF1 groups were enriched in “MYC-TARGETS-V2” “MYC-TARGETS-V1”, and “OXIDATIVE-PHOSPHORYLATION” pathways (Figure 7). Meanwhile, the “DNA-REPAIR” gene set was abundant in low-expression groups of CDH11, COL6A3 and EDNRA, and “PEROXISOME” was enriched in the COL6A3 and EDNRA low-expression groups.