First, we downloaded transcriptome RNA sequencing data, including 410 BCa and 19 normal tissue samples and 412 clinical data from The Cancer Genome Atlas (TCGA) database (https://tcga-data.nci.nih.gov/tcga/). The independent validation cohort (referred to as GSE32894) with 308 BCa samples was obtained from Gene Expression Omnibus (GEO) database. We further obtained the list of DDRGs from the GeneCards database (https://pathcards.genecards.org/). BCa was standardized according to the eighth edition of the American Joint Committee on Cancer (AJCC). We subsequently collected individuals or genes based on the following criteria: individuals with survival status and survival time more than 30 days; individuals with a gene expression matrix, and the mean and median of RNA read counts greater than one. At last, there were 392 samples and 18769 RNA from TCGA and 221 patients from GEO for further analysis. The descriptive statistics of clinical features of BCa patients involved in the current study are summarized in Table 1.

The transcriptome data were processed by using R x64 3.6.1 software (https://www.r-project.org/). The R package limma and Wilcox test were applied to obtain the differentially expressed genes. False-discovery rate (FDR) < 0.01 and log₂|fold change (FC)| > 1 were set as the cutoff value. Differentially expressed DDRGs were extracted from all of the differentially expressed genes, and it was performed as described previously [12].

Gene ontology analysis (GO) was applied to annotate differentially expressed DDRGs. The results of GO analysis were presented in three parts, namely, biological processes (BP), molecular functions (MF), and cellular components (CC). In addition, Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis was used to perform pathway enrichment analysis. Both GO analysis and KEGG analysis were conducted using R x64 3.6.1 software [12].

As described previously [13], a GSEA version 4.0.1 was first downloaded from the genomic enrichment analysis website (http:/Sofare.Broadstitute.org/GSEA/index.jsp). Then, the differentially expressed genes dataset was imported into the software with reference to the molecular signature database gene set (msigdb.v7.2.symbols gmt). Each analysis was repeated 1,000 times according to the default weighted enrichment statistical method using a genomic error detection rate (FDR) of <0.25 and a pedigree error rate (family‐wise error rate) < 0.05. GSEA analysis included four main statistical data: enrichment score (ES), standardized ES (normalized ES), error detection rate, and p-value.

In this step, we constructed the PPI network of prognosis-associated differentially expressed DDRGs that had been identified in previous analysis. The Search Tool for the Retrieval of Interacting Genes (STRING) database (version 11.0; https://string-db.org/cgi/input.pl) was used to evaluate the PPI information. Cytoscape software (version 3.8.1) was used to visualize the PPI networks and select hub genes for further discussion [12]. The CBioportal database was used to analyze mutations and copy number variations (CNVs) of differentially expressed DDRGs.

To identify survival-related genes, we integrated the expression of DDRGs with the OS of BCa patients. The relationships of DDRGs with OS were then analyzed by univariate Cox regression analysis. All samples in TCGA datasets were randomly assigned to the training dataset (60%) and the internal validation dataset (40%). These survival-related DDRGs were integrated into a stepwise multivariate Cox regression to select the optimal model for prognosis [14]. Finally, the prognostic model of BCa was established based on the multivariate co-efficiency multiplied by expression data. The formula was as follows: Risk score = coefficient(a) × gene expression(a) + coefficient(b) × gene expression(b) + ··· + coefficient(n) × gene expression(n). According to the median of risk scores, the patients were divided into high-risk and low-risk groups. The survminer package of R software was used to apply the Kaplan–Meier (KM) curve to investigate the connection between DDRGs and prognosis. The multivariate Cox analysis and stratification analysis were used to explore independent prognostic factors of BCa patients [15].

The Cistrome Cancer database (http://cistrome.org/CistromeCancer/), a comprehensive database of expression profiles and public ChIP-seq profiles from TCGA, predicts target genes and enhancer profiles of transcription factors (TFs) in TCGA cancer types [16]. Validated TFs were downloaded (318 in total, p < 0.05) with statistical relevance to the tumor. These data were combined and the differentially expressed TFs were used to draw a volcano map. Correlation tests between the TFs and prognosis-related DDRGs were performed, and the regulatory network of TFs-DDRGs was presented with Cytoscape [17].

The data were collected, analyzed, and presented using R software (version 3.6.1) and the corresponding software packages. The performance of the prognostic index was assessed by conducting the receiver operating characteristic (ROC) curve analysis.

We studied 410 BCa tissues and 19 normal tissues from the TCGA database. We identified 4,893 differentially expressed genes (DEGs) using the Wilcoxon signed-rank test; 3,468 of them were upregulated and 1,425 were downregulated. There were 123 differentially expressed DDRGs extracted from this set of genes, including 95 upregulated and 28 downregulated genes (Figures 1A,B).

According to gene ontology (GO) analysis in the biological process category, these differentially expressed DDRGs were chiefly enriched in “DNA replication,” “cell division,” and “mitotic nuclear division” (Supplementary Figure S1). GO analysis related to cellular components (CC) showed that they were mainly enriched in “nucleoplasm,” “nucleus,” and “cytoplasm” (Supplementary Figure S1B). GO analysis related to molecular functions (MF) demonstrated that they were involved in “protein binding,” “3′-5′DNA helicase activity,” and “chromatin binding” (Supplementary Figure S1C). KEGG pathway analysis revealed that “cell cycle,” “homologous recombination,” and “p53 signaling pathway” were the most enriched pathways (Supplementary Figure S1D). GSEA of the DEGs revealed that these genes were enriched in “SHEDDEN_LUNG_CANCER_POOR_SURVIVAL_A6”, “GOBERT_OLIGODENDROCYTE_DIFFERENTIATION_UP”, and “MARSON_BOUND_BY_E2F4_UNSTIMULATED” (Supplementary Figure S2).

To explore the prognostic potential of the differentially expressed DDRGs in BCa patients, we directed our efforts toward uncovering molecular biomarkers. First, we downloaded the clinical information of patients with BCa. By univariate Cox regression analysis, the top 22 differentially expressed DDRGs that were considered relevant with OS (Table 2). Most of these genes had positive hazard ratios, which indicated a higher risk for BCa patients.

Furthermore, the results of the enrichment analysis of the differentially expressed DDRGs showed that OS-related DDRGs were significantly enriched in similar GO terms and KEGG pathways, such as “negative regulation of transcription from RNA polymerase II promoter,” “G1/S transition of mitotic cell cycle,” and “positive regulation of smooth muscle cell proliferation” (Figure 2A). Also, “protein binding,” “protein kinase activity,” and “protein kinase binding” (Figure 2B) were the most significant MF. Furthermore, “nucleoplasm,” “nucleus,” and “cytoplasm” (Figure 2C) were the most significant CC. In addition, “Human papillomavirus infection,” “Pyrimidine metabolism,” and “MicroRNAs in cancer” were the most significant pathways (Figure 2D).

PPI network analysis of OS-associated differentially expressed DDRGs was performed using STRING (version 11.0) and visualized using Cytoscape (version 3.8.1). It showed that PBK, CDCA2, FOXM1, RRM2, ORC1, and HDAC4 were the hub genes (Supplementary Figure S3). According to mutation and CNV analysis for BCa performed using the cBioportal database, the three most common types of OS-related DDRGs included amplification, missense mutations, and deep deletions (Supplementary Figure S4). We also performed the regression analysis of risk score on tumor mutational burden with the p-value of 0.96. The scatter plot was shown in Supplementary Figure S7.

To investigate the possible molecular mechanisms regulating the OS-related DDRGs, we analyzed the TF-DDRG interaction. First, we observed the expression profiles of 318 TFs; there were 77 TFs that were differentially expressed between BCa tissues and non-tumor tissues. Among them, 41 were upregulated and 33 were downregulated (Figure 3A). We then constructed a regulatory network based on these 77 TFs and 22 OS-related DDRGs. We set cutoff thresholds as correlation scores >0.4 and p-values < 0.001. The TF-based regulatory network is shown in Figure 3B. We identified the key regulated factors of OS-related DDRGs using the TRRUST database. Seven key transcription factors (NCAPG, DNMT1, LMNB1, BRCA1, E2H2, CENPA, and E2F7) were found to be associated with the regulation of these DDRGs.

We separated the TCGA-BLCA cohort into training dataset (60%) and internal validation dataset (40%) by random selection. The following optimal model was chosen by a stepwise multivariate Cox analysis: [Expression level of CAD * (0.042)] + [Expression level of HDAC10*(−0.146)] + [Expression level of JDP2 * (0.069)] + [Expression level of LDLR * (0.022) + [Expression level of PDGFRA * (0.070)] + [Expression level of POLA2 * (0.108)] + [Expression level of SREBF1 * (0.008)] + [Expression level of STAT1 *(−0.009)] (Table 3). We built a prognostic index-risk score to divide patients with BCa into two groups, the high-risk group and the low-risk group, and then we made a risk curve. The KM analysis showed that there was a statistically significant difference in OS between the two groups with p-value < 0.0001 (Figure 4A). A higher risk score was associated with a shorter survival time. The area under the ROC curve (AUC) was 0.771, p = 2.34E-07 (Figure 4A), which indicated a high forecast ability of the DDRG-based risk score in survival surveillance. The prognostic index was confirmed in TCGA-BLCA internal validation dataset, with a similar KM plot (p-value <0.0001) and ROC curve (AUC = 0.743, p = 3.77E-06) (Figure 4B). The survival status and risk scores of each BCa patient are listed in Supplementary Figure S5A for the TCGA-BLCA training dataset, Supplementary Figure S5B for the internal validation dataset, and Supplementary Figure S5C for the GSE32894 cohort.

Furthermore, we validated the performance of the prognostic model in an independent dataset, GSE32894, through the KM survival plot and ROC curve analysis. A significant p-value < 0.0001 in the KM survival analysis and an area under the ROC curve of 0.827 verified the efficiency of our prognostic model, p = 2.29E-07 (Figure 4C). To compare the prediction ability of risk score with clinical features, we conducted ROC analyses for each clinical feature (Supplementary Figure S6A) and combined clinical features and risk scores (Supplementary Figure S6B). The results showed that the AUC of the prognostic index was larger than those based on age, count of examined lymph nodes, positive lymph nodes, and tumor stage (Supplementary Figure S6A). When combining clinical features and risk scores, we found that our prognostic index could improve the clinical prognostication (Supplementary Figure S6B).

Univariate Cox analysis revealed that age, number of lymph nodes, number of positive lymph nodes, stage, and risk scores were related to OS. After adjusting for age, gender, race, tobacco smoking history, number of lymph nodes, number of positive lymph nodes, and clinical stage by multivariate Cox analysis, the risk score was still associated with OS, with a p-value <0.001 (Table 4). However, we noticed that age, number of lymph nodes, and number of positive lymph nodes were also related to OS (p < 0.05). In addition, chemotherapy treatment is an important factor to OS. To further confirm the independence of the risk scores, we performed a stratification analysis. All of the patients were separated into two subgroups as shown in Table 1. We found that the risk scores still correlated with OS (p < 0.05, Figure 5).

We additionally accessed the correlation between DDRG index and subtypes in TCGA and GEO database (p < 0.001, Supplementary Figure S8). Except SREBF1, other DDRGs were correlated with diagnosis subtypes between Non-Papillary and Papillary (Table 5). Except JDP2, the expression of each gene was related to the molecular subtypes (Table 6).