RNA sequencing profiles of 408 BLCA samples were obtained from The Cancer Genome Atlas (TCGA) database (https://www.cancer.gov/tcga), and 209 BLCA samples with complete demographic and prognostic characteristics, especially bone metastasis information integrity, were eventually selected. Among the 209 BLCA samples, 11 had bone metastasis. Additionally, information on cancer-related transcription factors (TFs) was accessed through the Cistrome database (http://cistrome.org) [13]. In addition, 2498 immune-related genes and their gene signatures were retrieved from the ImmPort database (https://immport.org) [14] and the Molecular Signatures Database (MSigDB) v7.0 (https://www.gsea-msigdb.org) [15], respectively. What is more, the single-cell RNA sequencing profiles of the tumor sample, as well as CD4 and CD8 samples of BLCA patients were downloaded from the Gene Expression Omnibus database (https://www.ncbi.nlm.nih.gov/geo/) (accession no. GSE164046, no. GSE149652).

Differentially expressed genes (DEGs) were identified between bone-metastatic and non-metastatic samples with the cutoff value of false discovery rate (FDR) 0.05 and $|$log2(fold change)$|$ $>$1 through the edgeR package. And DEIGs were subsequently recognized by matching DEGs with the 2498 immune-related genes downloaded from the ImmPort database. Gene expression values were pre-processed by z-score before fold change and p value were calculated. In addition, Gene Oncology (GO) and KEGG enrichment analysis was conducted. Subsequently, univariate Cox proportional hazards model was constructed to identify DEIGs significantly correlated with overall survival. Disease-free survival analysis of these identified DEIGs was also obtained from the Gene Expression Profiling Interactive Analysis (GEPIA) database (http://gepia.cancer-pku.cn/) [16].

Using univariate and multivariate Cox regression analyses, DEIGs with statistical significance were included in the proportional hazards model. Risk score for prognosis was calculated according to the following formula:

$${\mathit{\text{Risk score}}}_i={\sum }_{a=1}^n\beta a\times (\mathit{\text{expression level of gene }}a)$$

In the formula, “i” stands for the number of the patient, “n” stands for the number of significant DEIGs included in the model, and “$\beta$a” stands for the regression coefficient of the DEIGs. The median cut point of the risk scores was used as the cutoff value for high and low risk grouping. Model predictions were evaluated by plotting the Kaplan-Meier curve. Predictive performance was quantified using area under the curve (AUC). In order to further estimate the independent prognostic value of risk score, the univariate and mulitivariate Cox regression analyses were utilized. Moreover, 196 BLCA patients cohort from the Therapeutically Applicable Research To Generate Effective Treatments (TARGET) database (https://ocg.cancer.gov/programs/target) were extracted for further validation. Differential enrichment of tumor infiltrating immune cells and immune pathways in BLCA was evaluated using the CIBERSORT and ssGSEA methods.

Differentially expressed TFs were identified by the edgeR package. FDR was set to be less than 0.05 while $|$log2(fold change)$|$ greater than 1. Then, correlation analysis was performed to co-analyze DEIGs and TFs, illustrating their regulatory relationships with correlation coefficient $>$0.300 and p 0.001, in which case if the variables conformed to normal distribution and homogeneity of variance, Pearson correlation analysis was conducted, or else Spearman correlation analysis was used.

Differential enrichment of tumor infiltrating immune cells/pathways in BLCA was evaluated using the CIBERSORT and ssGSEA tools. The correlation between key DEIGs and immune cells/pathways was evaluated by Pearson/Spearman correlation coefficient. In addition, GSVA was applied to score KEGG pathways in each sample. To explore the potential KEGG pathways correlated with prognosis and metastasis in BLCA patients, Pearson/Spearman correlation analysis and univariate Cox regression were conducted. Moreover, the Venn diagram was plotted using the VennDiagram package [17].

Focusing on the immune genes, TFs and downstream immune cells/pathways and KEGG pathways were identified which were involved in the bone metastasis in BLCA. Therefore, we constructed a regulatory network to further elucidate the mechanism of bone metastasis in BLCA by exploiting the Cytoscape package [18].

To further validate the transcriptional regulation pattern between EBF1 and RBP7, data related to EBF1 acquired from chromatin immunoprecipitation sequencing (ChIP-seq) of in vitro human cell lines were downloaded from Gene Expression Omnibus (GEO) (GEO assession: samples GSM1958038 and GSM1958039 from GSE75503, GSM935375 from GSE31477, and GSM803386 from GSE32465) [19, 20]. ENCODE Transcription Factor Targets [21] and JASPAR Predicted Transcription Factor Targets [22] algorithms were then applied to re-predict the transcriptional regulation pattern between EBF1 and RBP7. Additionally, ATAC-seq data of BLCA samples were obtained from the TGCA project of chromatin accessibility landscape of primary human cancers (https://gdc.cancer.gov/about-data/publications/ATACseq-AWG), which were then utilized to explore the chromatin accessibility in specific locations of EBF1 and RBP7. What is more, the scRNA-seq profiles of BLCA samples from GEO were obtained to successively go through quality control, normalization, cell cycle assignment, feature selection, principal component analysis (PCA), Uniform Manifold Approximation and Projection (UMAP)/t-distributed Stochastic Neighbor Embedding (t-SNE) for dimensionality reduction and visualization, unsupervised clustering, differential expression analysis, and cell communication analysis.

Multiple online databases were used to access the external validation of our scientific hypothesis. CellMarkers (http://biocc.hrbmu.edu.cn/CellMarker/) [23] and PathCards (https://pathcards.genecards.org/) [24] were applied to find biomarkers related to immune genes/pathways and KEGG pathways. What is more, UALCAN (http://ualcan.path.uab.edu/index.html) [25], GEPIA (http://gepia.cancer-pku.cn/) [16], LinkedOmics (http://linkedomics.org/) [26], cBioportal (https://www.cbioportal.org/) [27], TISIDB (http://cis.hku.hk/TISIDB/) [28], UCSC xena (https://xena.ucsc.edu/) [29], OncoLnc (http://www.oncolnc.org/) [30], and PROGgeneV2 (http://www.progtools.net/gene/) [31] were employed for the external validation for our scientific hypothesis.

Survival profiles were analyzed by Kaplan–Meier survival analysis. Additionally, correlation analysis were performed using Pearson/Spearman correlation analysis. In addition, Biomarkers with prognostic value were distinguished by Cox univariate and multivariate analyses. All bioinformatic and statistical analyses were implemented with R software (Version 3.6.1; Institute for Statistics and Mathematics, Vienna, Austria; https://www.r-project.org). A significance level of a two-sided p value less than 0.05 was adopted.

The analysis procedure of this study was demonstrated in Fig. 1. Clinical information (Table 1) and RNA-seq data of 209 BLCA samples were accessed through TCGA database. Between bone-metastatic and non-metastatic samples, 1090 DEGs were identified (Fig. 2A,B). Additionally, GO and KEGG enrichment analyses, visualized as bubble plots, were performed to uncover that these DEGs actively participated in pathways such as axonogenesis, metal ion transmembrane transporter activity, and neuroactive ligand-receptor interaction which might play a role in the bone metastasis of BLCA patients (Supplementary Fig. 1A,B). In addition, using differential analysis, 86 of 1090 DEGs were found in the list of 2498 immune genes, which were defined as DEIGs and displayed in a heatmap plot (Fig. 2C) and a volcano plot (Fig. 2D). Moreover, in univariate Cox regression analysis, 31 DEIGs were exhibited to have a significant association with overall survival, in which 29 risk factors and 2 protective factors were identified in a forest plot (Fig. 2E). What is more, the differences of disease free survival between low and high expression of the 31 DEIGs were displayed in survival plots (Supplementary Fig. 2).

Fig. 1.

The analysis flowchart.

Fig. 2.

Differential gene expression analysis and univariate Cox proportional hazards model. (A) Heatmap demonstrated DEGs between bone-metastatic and non-metastatic BLCA samples. (B) Volcano plot showed over-expressed and under-expressed DEGs with p 0.05. (C) Heatmap demonstrated DEIGs between bone-metastatic and non-metastatic BLCA samples. (D) Volcano plot showed over-expressed and under-expressed DEIGs with p 0.05. (E) Forest plot illustrated 31 DEIGs correlated with overall survival by univariate Cox proportional hazards model. Abbreviation: DEGs, differentially expressed genes; DEIGs, differentially expressed immune genes.

Based on the 31 prognosis-related DEIGs, we constructed a multivariate prognostic prediction model. The receiver operating characteristic (ROC) curve of the model was exhibited in Fig. 3A, in which the AUC was 0.793, indicating a good predictive value. Sub-group survival analysis was conducted according to the median risk score, displaying a significant difference (p = 1.04 $\times$ 10${}^{-7}$) between the low risk and high risk groups (Fig. 3B). The clinical features of risk score were represented by a risk curve (Fig. 3C) and a scatter plot (Fig. 3D). The expression levels of the 31 DEIGs between low-risk and high-risk groups were demonstrated by a heatmap (Fig. 3E). Together with four clinical features (age, gender, M and stage), the risk score in our model proved to be an independent prognostic predictor by univariate (HR, 1.011; 95% CI, 1.005–1.017; p 0.001) and multivariate (HR, 1.009; 95% CI, 1.004–1.015; p 0.001) Cox regression analyses (Fig. 4A,B). Moreover, significant differences of censor (p 0.001), grade (p 0.05), T (p 0.05), N (p 0.001), and stage (p 0.001) between low risk and high risk groups in a 196 BLCA patients cohort from the TARGET database (Fig. 4C). Additionally, the distribution of different immune cell components in low risk and high risk groups was demonstrated in Supplementary Fig. 3A. What is more, the comparison of different immune cell fractions (Supplementary Fig. 3B) and ssGSEA scores (Supplementary Fig. 3C) between low risk and high risk group were displayed, which revealed that the risk score was correlated with expression of CD8 T cells (p 0.001), T cells follicular helper (p 0.05), and macrophages M0 (p 0.01), along with immune functions of NK cells (p 0.001) and Th2 cells (p 0.001).

Fig. 3.

Construction of a multivariate prognostic model. (A) ROC curve revealed a good prediction reliability (AUC = 0.793). (B) Kaplan-Meier survival curve showed a significant difference (p = 1.04 $\times$ 10${}^{-7}$) between low risk and high risk groups. (C) Risk curve displayed the distribution of patients with increasing risk scores in the low risk and high risk groups. (D) Scatter plot represented the distribution of survival time and survival status of low risk and high risk groups. (E) Heatmap demonstrated differential expression of the 31 DEIGs between low risk and high risk groups. Abbreviations: DEIGs, differentially expressed immune genes.

Fig. 4.

Independent prognostic analysis. (A) Univariate Cox regression analysis with age, gender, M, stage, and risk score. (B) Multivariate Cox regression analysis with age, gender, M, stage, and risk score. (C) Heatmap visualized differences of censor, age, gender, grade, TNM stage, clinical stage, and bone metastasis status between low risk and high risk groups in 196 TARGET patients (*: p 0.05; ***: p 0.001).

8 differentially expressed TFs ($|$log2FC$|$ $>$1 and FDR 0.05) were illustrated by a heatmap (Fig. 5A) and a volcano plot (Fig. 5B). Additionally, a venn diagram was utilized to identify 5 DEIGs which were KLRK1, IL9R, CNTFR, RBP7, and TNFSF4 related to metastasis, stage, and TF co-expression simultaneously (Fig. 5C). Aiming to identify the regulatory relationships between the identified TFs and key DEIGs, Pearson correlation analysis was utilized and 30 regulatory relationships were identified with $|$correlation coefficient$|$ $>$0.300 and p value 0.001. Eventually, EBF1 positively regulating RBP7 was selected as the most potential regulatory relationship (R = 0.677, p value 0.001). In addition, the significant differences of RBP7 expression between metastatic and non-metastatic BLCA samples (p = 0.018), as well as among stage i, ii, iii, and iv of BLCA samples (p = 0.018) were demonstrated (Fig. 5D,E).

Fig. 5.

Identification of differentially expressed TFs and co-expression analysis. (A) Heatmap revealed differentially expressed TFs between metastatic and non-metastatic BLCA samples. (B) Volcano plot showed over-expressed and under-expressed TFs with p 0.05 and $|$log2(fold change)$|$ $>$1. (C) Venn diagram presented the overlapped number of DEIGs related to bone metastasis, stage and co-expression. (D) Bee-swarm plot illustrated the differential RBP7 expression (p = 0.018) between metastatic and non-metastatic BLCA samples. (E) Bee-swarm plot displayed the differential RBP7 expression (p = 0.018) among stages i, ii, iii, and iv of BLCA samples. Abbreviations: TFs, transcription factors; DEIGs, differential expressed immune genes.

For the purpose of exploring the potential signaling pathways related to bone metastasis and BLCA prognosis, CIBERSORT, ssGSEA, and GSVA algorithms were utilized. The correlation between expression of immune cells/pathways and RBP7 expression was analyzed by Pearson’s correlation (Supplementary Fig. 4A). Moreover, 5 immune cells/pathways were identified to be statistically significant (Supplementary Fig. 4B), including macrophages M0 (R = 0.212, p value = 0.007), Type 2 helper T (Th2) cells (R = 0.23, p value 0.001), parainflammation (R = –0.201, p value = 0.003), neutrophils (R = –0.179, p value = 0.009) and MHC class I (R = –0.148, p value = 0.032). What is more, 9 prognosis-related KEGG pathways were characterized to be significantly correlated with RBP7 expression (Supplementary Fig. 5A,B).

Using the GSEA, 23 metastasis-related KEGG pathways were identified. Thus, as illustrated in the Venn plot (Fig. 6A), the oocyte meiosis pathway was identified to be the key KEGG pathway both significant in GSEA and GSVA co-expression (Fig. 6B). What is more, the correlation relationship between the oocyte meiosis pathway and RBP7 (R = 0.14, p = 0.042) was shown in Fig. 6C. Besides, GSEA result of the oocyte meiosis was demonstrated by the enrichment plot in Fig. 6D.

Fig. 6.

Co-expression and correlation analysis of the KEGG oocyte meiosis pathway. (A) Venn plot for the KEGG pathway selected using GSVA and GSEA algorithms. (B) GSEA analysis for the oocyte meiosis pathway. (C) RBP7 expression was positively correlated with the oocyte meiosis pathway. (D) GSEA analysis showed a positive correlation between the oocyte meiosis pathway and RBP7 expression. Abbreviation: KEGG, Kyoto Encyclopedia of Genes and Genomes; GSVA, Gene-Set Variation Analysis; GSEA, Gene-Set Enrichment Analysis.

The regulatory network based on TFs, RBP7, immune genes/pathways and KEGG pathways was constructed in Fig. 7A. Taken together, we proposed that the EBF1-RBP7 axis (representing TF-immune gene relationship) played an important role in bone metastasis of BLCA through Th2 cells and the oocyte meiosis pathway.

Fig. 7.

Construction of the regulatory network. (A) Regulatory network based on TFs transcription factors, RBP7, the 5 filtered immune cells/pathways, and 9 filtered KEGG pathways. (B) Correlation analysis between EBF1 and RBP7 (R = 0.71, p 0.001). (C,D) In ATAC-seq data of BLCA samples, multiple binding peaks were identified in EBF1 and RBP7 sequences. Abbreviations: TFs, transcription factors; KEGG, Kyoto Encyclopedia of Genes and Genomes.

ATAC-seq analysis was performed to explore the correlation between EBF1 and RBP7. Open chromatin loci on different chromosomes were displayed in Supplementary Fig. 6A. Additionally, the venn pie and upsetplot showed the distribution and intersection of different pick types such as genic, intergenic, intron, and exon in Supplementary Fig. 6B. Moreover, the distribution of binding loci relative to TSS was visualized in Supplementary Fig. 6C. Additionally, the correlation between EBF1 and RBP7 was analyzed, which turned out that the expression of EBF1 was positively correlated with RBP7 (R = 0.71, p 0.001) (Fig. 7B). What is more, multiple binding peaks in BLCA samples at promoters of EBF1 and RBP7, as well as various regulatory elements binding areas in the introns and in introns of neighboring genes were identified, which indicated that these regions might function as potential regulatory elements on upstream of EBF1 and RBP7 sequences (Fig. 7C,D).

To further illustrate the co-expression patterns between EBF1 and RBP7, ENCODE and JASPAR algorithms were utilized. In EBF1 ChIP-seq data of activated GM12878, LCL and Mutul cell line (homo sapiens), binding peaks were all found in RBP7 sequences (Fig. 8A), and 2 motifs of EBF1 were identified (Fig. 8B,C). Consequently, the transcriptional regulation pattern between EBF1 and RBP7 could be determined.

Fig. 8.

The results of ChIP-seq validation. In EBF1 ChIP-seq data of activated GM12878, LCL and Mutul cell lines (homo sapiens), binding peaks were all found in RBP7 sequences (A). Two EBF1-binding motifs (B,C).

Furthermore, in the scRNA-seq analysis of BLCA tumor cell lines, 12 seurat clusters and bulk labels were shown in t-SNE dimension reduction analysis (Fig. 9A). Additionally, the top 5 marker genes of each cluster was demonstrated in a heatmap (Fig. 9B). In addition, up- and down-regulated marker genes of each cluster were revealed in a point plot (Fig. 9C). More importantly, the distribution of 5 marker genes expression including EBF1, RBP7, CD24, MKI67, and CD44 was illustrated in Fig. 9D, which indicated that EBF1 and PBP7 had a co-expression relationship in clusters 0, 1, 2, 3 and 4, and they might also be co-expressed with tumor stemness and proliferation markers CD24, MKI67, and CD44. What’r more, the distribution of different stages in a cell cycle and relevant scores were visualized in Fig. 9E, which suggested that BLCA cells with high EBF1 and RBP7 expression had a relative high fraction of stages G2M and S cells, including clusters 2, 3 and 4. Last but not least, the cell communication among different clusters presented that clusters 1, 2 and 4 had the strongest interaction with other clusters (Fig. 9F). Moreover, in the scRNA-seq analysis of BLCA CD4 and CD8 T cell lines (Supplementary Fig. 7), a co-expression relationship between EBF1 and RBP7 was also identified.

Fig. 9.

Single-cell RNA sequencing of BLCA tumor cell lines. (A) t-SNE dimension reduction analysis showed 12 seurat clusters and bulk labels in BLCA tumor cell lines. (B) Top 5 marker genes of each cluster in a heatmap. (C) Marker genes of each cluster in a point plot. (D) The distribution of 5 marker genes including EBF1 and RBP7. (E) UMAP dimension reduction analysis showing the distribution of different stages in a cell cycle as G1, G2M and S, along with G2M.score and S.score in violin plot. (F) Ligand-receptors plot among 12 clusters.

In order to access the external validation of our scientific hypothesis, multiple online databases were employed. Three biomarkers of Th2 cells (IL1RL1, GATA3, and IL33) and the top five genes related to the oocyte meiosis pathway (ADCY3, ADCY5, ADCY9, MAPK3, and IGF1) were presented using the CellMarker and PathCards databases, respectively. All the external validation results are summarized in Table 2.

In UALCAN database validation (Supplementary Figs. 8,9), high ADCY3 (Supplementary Fig. 8F) and MAPK3 (Supplementary Fig. 8I) expression levels were characterized as being significantly correlated with worse BLCA patients’ prognosis. Additionally, expression levels of EBF1 (Supplementary Fig. 9A), RBP7 (Supplementary Fig. 9B), IL33 (Supplementary Fig. 9E), ADCY5 (Supplementary Fig. 9G), ADCY9 (Supplementary Fig. 9H), and IGF1 (Supplementary Fig. 9J) were identified to be elevated in normal tissues compared to BLCA primary tumor tissues, while GATA3 expression was escalated in BLCA primary tumor tissues (Supplementary Fig. 9D).

In GEPIA database validation (Supplementary Figs. 10,11), it was discovered that elevated EBF1, RBP7, IL1RL1, IL33, ADCY3, ADCY5, ADCY9, MAPK3, and IGF1 (Supplementary Fig. 10A–H,J) expression levels were associated with unfavorable prognoses in BLCA patients, while over-expression of GATA3 (Supplementary Fig. 10I) indicated better overall survival. In addition, ADCY3 (Supplementary Fig. 11A), IGF1 (Supplementary Fig. 11B), and EBF1 (Supplementary Fig. 11G) expression were recognized to be positively correlated with RBP7 expression, while GATA3 (Supplementary Fig. 11I) expression was negatively correlated with RBP7 expression.

In LinkedOmics database validation (Supplementary Figs. 12,13), it was shown that escalated expression of ADCY3 (Supplementary Fig. 12A), IGF1 (Supplementary Fig. 12B), ADCY9 (Supplementary Fig. 12E), MAPK3 (Supplementary Fig. 12H), and RBP7 (Supplementary Fig. 12J) were significantly correlated with poor BLCA patients’ prognosis, while over-expressed GATA3 (Supplementary Fig. 12I) suggested good prognosis. Moreover, ADCY3, IGF1, ADCY5, IL1RL1, ADCY9, IL33, EBF1, MAPK3, and GATA3 expression (Supplementary Fig. 13A–H) were illustrated to be positively correlated with RBP7 expression, while GATA3 (Supplementary Fig. 13I) expression was negatively correlated with RBP7 expression.

In TISIDB database validation (Supplementary Fig. 14), up-regulation of ADCY3 (Supplementary Fig. 14A), IGF1 (Supplementary Fig. 14B), ADCY9 (Supplementary Fig. 14E), and MAPK3 (Supplementary Fig. 14H) were significantly linked to worse BLCA patients’ prognosis, while up-regulated GATA3 (Supplementary Fig. 14I) expression was significantly related to better prognosis of BLCA patients.

In cBioportal database validation (Supplementary Fig. 15), ADCY3, IGF1, ADCY5, IL1RL1, ADCY9, IL33, EBF1, and MAPK3 (Supplementary Fig. 15A–H) expression were positively correlated with RBP7 expression, while expression of GATA3 (Supplementary Fig. 15I) was negatively associated with expression of RBP7.

In UCSC Xena database validation (Supplementary Fig. 16), down-regulation of ADCY3, IGF1, ADCY5, IL1RL1, ADCY9, IL33, and EBF1 (Supplementary Fig. 16A–G) were significantly correlated with better BLCA patients’ prognosis, while down-regulated GATA3 (Supplementary Fig. 16I) expression was significantly correlated with better BLCA patients’ prognosis.

In OncoLnc database validation (Supplementary Fig. 17), lower expression of ADCY3, IGF1, ADCY5, IL1RL1, ADCY9, IL33, EBF1, and RBP7 (Supplementary Fig. 17A–G,J) were revealed to be significantly correlated with favorable prognosis of BLCA patients, while down-regulated GATA3 (Supplementary Fig. 17I) expression was significantly correlated with favorable prognosis of BLCA patients.

In PROGgeneV2 database validation (Supplementary Figs. 18,19), low-expressed IGF1 (Supplementary Fig. 18B), ADCY9 (Supplementary Fig. 18E) and RBP7 (Supplementary Fig. 18J) were demonstrated to be significantly associated with better BLCA patients’ prognosis. What is more, down-regulated combined gene expression of EBF1, RBP7, IL1RL1, GATA3, IL33, ADCY3, ADCY5, ADCY9, MAPK3, and IGF1 was significantly correlated with favorable prognosis of BLCA patients (Supplementary Fig. 19).