The transcriptome and clinical data of BCa were obtained from TCGA database (https://portal.gdc.cancer.gov), including 414 BCa cases and 19 normal cases [14]. The clinical features of BCa patients included age, gender, grade, T stage, N stage, M stage, and AJCC stage. A total of 412 BCa cases have clinical data; however, there was only a total of 396 BCa cases with unabridged transcriptome and clinical data simultaneously. Besides, a total of 165 BCa patients with complete mRNA expression profile and clinical data were downloaded from GSE13507 dataset in GEO database (https://www.ncbi.nlm.nih.gov). All these data were processed using Perl language and R language. The clinicopathologic feature of TCGA cohort and GEO cohort was showed in Table 1. We used the Ensemble database (http://asia.ensembl.org/signature.html) to convert Ensemble IDs into gene symbols.

The Molecular Signatures Database (MSigDB) (https://www.gsea-msigdb.org/gsea/msigdb) is a collection of annotated gene sets using GSEA software. We extracted 79 URGs from REACTOME PROTEIN UBIQUITINATION gene set (C2: curated gene sets; systematic name: M27742) from MSigDB database. All these URGs were related to protein ubiquitination. The detailed list of these 79 URGs were presented in Supplementary Table S1.

Firstly, the expression matrix of 79 URGs was extracted from BCa transcriptome and merged with overall survival (OS) time using Perl language. There was a total of 396 BCa cases with complete OS data and mRNA expression data in TCGA database. We then performed univariable Cox regression analysis to screen prognostic URGs associated with OS. The cut-off p value was set as 0.05. Next, consensus clustering analysis was performed for identifying ubiquitination-related molecular subtypes of BCa patients using R package “ConsensusClusterPlus”. The association between ubiquitination-related molecular subtypes and OS was explored using R package “survival” and “survminer”. The associations between ubiquitination-related subtypes and clinicopathologic characteristics (including gender, grade, T stage, N stage, stage, and age) were presented by utilizing R package “pheatmap”.

First of all, we divided randomly all BCa patients in TCGA database into training cohort and testing cohort (internal verification). All cases in GSE13507 dataset were used as validating cohort (external verification). Then, univariate and multivariate Cox regression analysis were conducted to establish a ubiquitination-based prognostic index for predicting OS of BCa using training cohort.

Next, all cases in training cohort, testing cohort and validating cohort were categorized into low-risk group and high-risk group based on the median risk score. We then performed survival analysis and the time-dependent receiver operating characteristic (ROC) curve to explore the performance of this ubiquitination-based prognostic index in training cohort, testing cohort and validating cohort. The expression heatmap, the distribution of risk score and survival time of training cohort, testing cohort and validating cohort were presented using “pheatmap” R package. Moreover, we performed univariate and multivariate independent prognostic analysis to demonstrate whether this ubiquitination-based prognostic index was an independent predictor of OS in BCa patients.

The ESTIMATE algorithm was utilized to evaluate tumor microenvironment (TME) scores while the CIBERSORT method was used to calculate the scores of 22 types of immune infiltrating cells [15,16]. Then, we investigated the associations of ubiquitination-related molecular subtypes with BCa tumor microenvironment and immune infiltrating cells. In addition, the association of the ubiquitination-related molecular subtypes with PD-L1 gene expression level [17] was also explored.

Single sample gene set enrichment analysis (ssGSEA) was conducted to calculate the infiltrating score of 16 types of immune infiltrating cell and the activity of 13 types of immune-related function in TCGA cohort using R package “gsva”. We then investigated the associations of the ubiquitination-based prognostic index with immune infiltrating cells and immune function activity. Besides, the associations of ubiquitination-based prognostic index with TME scores were explored.

UALCAN database (http://ualcan.path.uab.edu/) is a portal for facilitating tumor subgroup gene expression and survival analyses. The mRNA expression levels of risk URGs between normal and tumor tissues were demonstrated using UALCAN database. The prognostic value of risk URGs was validated in GEO cohort using univariate Cox regression analysis. Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment for high-risk group and low-risk group was performed using Gene Set Enrichment Analysis (GSEA) method in the whole TCGA cohort.

Statistical analysis was performed utilizing R programming language. Univariate and multivariate Cox regression analysis were performed to establish a ubiquitination-related prognostic index for predicting OS of BCa. Univariate and multivariate independent prognostic analysis were used to demonstrate whether this ubiquitination-related prognostic index was an independent predictor of OS. Survival analysis and the time-dependent receiver operating characteristic (ROC) curve were performed to explore the performance of ubiquitination-related prognostic index. Statistical significance was considered at the level of p value <0.05.

The flowchart of this study was presented in Figure 1. Firstly, the expression matrix of 79 URGs was extracted from TCGA database and merged with OS status and time. We then conducted univariable Cox regression analysis in the whole TCGA cohort to screen URGs associated with OS. The results showed that there was a total of six URGs associated with OS, including CDC73, PRKDC, RNF40, TMEM129, UBE2N and USP5. Next, the expression matrix of these six URGs was used to conduct consensus clustering analysis to identify ubiquitination-related molecular subtypes of BCa. Finally, a total of four ubiquitination-related molecular subtypes of BCa were identified, including 116 cases of subtype 1, 90 cases of subtype 2, 132 cases of subtype 3, 58 cases of subtype 4 (Figures 2A,B). The difference of OS among these ubiquitination-related molecular subtypes was statistically significant (p = 0.005, Figure 2C). As indicated by heatmap and histogram, these four molecular subtypes have significantly different grade (p < 0.05, Figures 2D,E ).

We explored the discrepant expression level of PD-L1 in TCGA database among these four molecular subtypes. The results showed that the expression level of PD-L1 in subtype 3 and subtype 4 was significantly increased compared with that in subtype 1 and subtype 2. However, the PD-L1 expression level was similar between subtype 1 and subtype 2, between subtype 3 and subtype 4, respectively. (Figure 2F).

We used the CIBERSORT method to calculate the scores of 22 types of immune infiltrating cells and the ESTIMATE algorithm to evaluate TME. According to the CIBERSORT method, the activated mast cell was the exclusively discrepant immune infiltrating cell among these four ubiquitination-related subtypes (Figure 3). The ESTIMATE scores, immune scores, stromal scores in subtype 2 and subtype 3 were significantly higher in comparison with that in subtype 1 and subtype 4; the tumor purity scores in subtype 1 and subtype 4 were significantly increased compared with that in subtype 2 and subtype 3. However, the ESTIMATE scores, immune scores, stromal scores and tumor purity scores were similar between subtype 1 and subtype 4, between subtype 2 and subtype 3, respectively (Figure 4).

There was a total of 200 cases in training cohort, 196 cases in testing cohort and 165 cases in validating cohort. Firstly, univariate and multivariate Cox regression analyses were conducted to develop a novel ubiquitination-based index utilizing six URGs (including HLA-A, TMEM129, UBE2D1, UBE2N, UBE2T and USP5) for prognosis prediction of BCa using training cohort. The calculation formula of risk score is shown as follows: Risk score= (−0.0010853) * HLA-A + (−0.03553995) * TMEM129 + (−0.08324461) * UBE2D1 + (−0.07059857) * UBE2N + (0.01191996) * UBE2T + (0.02752775) * USP5 (Table 2).

Then, we calculated the risk score of each patient using above calculation formula in training cohort, testing cohort, and validating cohort. All patients were divided into high-risk score group and low-risk score group according to the median risk score. The difference in OS between low-risk and high-risk group was statistically significant in training cohort (p < 0.001), testing cohort (p = 0.003), whole TCGA cohort (p < 0.001) and validating cohort (p = 0.035), respectively. Low-risk score was associated with significantly preferable OS in comparison with high-risk score in training cohort, testing cohort, whole TCGA cohort and validating cohort, respectively. The area under ROC curve (AUC) for OS prediction was 0.736, 0.723, 0.723 and 0.683 in training cohort, testing cohort, whole TCGA cohort and validating cohort, respectively, suggesting the promising value of this novel ubiquitination-based prognostic index for prognosis prediction of BCa (Figure 5). The expression heatmaps, the distributions of risk score and survival time of training cohort, testing cohort, whole TCGA cohort and validating cohort were presented in Figure 6.

Univariate and multivariate independent prognostic analysis were used to demonstrate whether this index was an independent predictor of OS in BCa patients. Univariate analysis showed that age (p = 9.97E-05), stage (p = 9.56E-08), T stage (p = 6.08E-05), N stage (p = 1.57E-07) and this prognostic index (p = 1.96E-10) were associated with OS of BCa. Multivariate analysis demonstrated that only age (p = 0.000189) and this prognostic index (p = 3.82E-06) were independent predictors for OS of BCa, indicating the great performance of this prognostic index (Table 3).

The difference of clinicopathologic features between high and low risk group were explored. The results demonstrated that ubiquitination-related molecular subtypes (p < 0.001), age (p < 0.001), grade (p = 0.002), T stage (p = 0.003), N stage (p = 0.002), AJCC stage (p < 0.001) and were significantly different between high-risk and low-risk patients. However, gender (p = 0.365) was similar between high-risk and low-risk patients (Figures 7A–G).

Subgroup survival analysis demonstrated that this prognostic index was especially suitable for older, male, high grade, AJCC stage III-IV, stage N0, stage T3-4 BCa patients. Low-risk score was associated with significantly preferable OS in comparison with high-risk score in patients of age >65 or male or high grade or AJCC stage III-IV or N0 stage or T3-4 stage. However, the difference in OS between low-risk and high-risk group was not statistically significant in patients of age ≤65 or female or AJCC stage I-II or N1-3 stage or T1-2 stage (Figures 7H–R). Besides, considering the extremely low case number of patients with low grade, we did not conduct subgroup survival analysis in this subgroup.

We calculated the risk score for each patient of each ubiquitination-related molecular subtype. The results revealed that this prognostic index was especially applicable in BCa patients of subtype 1 and subtype 3. Low-risk score was associated with significantly preferable OS in comparison with high-risk score in BCa patients of subtype 1 and subtype 3. The area under ROC curve (AUC) for OS prediction was 0.778, 0.713 in BCa patients of subtype 1 and subtype 3, respectively. However, the difference in OS between low-risk and high-risk group was not statistically significant in subtype 2 and subtype 4 patients (Figure 8).

The immune functions and immune cells infiltration were quantified using ssGSEA in TCGA cohort. The infiltrating proportions of CD8+ T cells, interdigitating dendritic (iDC) cells, natural killer (NK) cells, T helper 2 (Th2) cells in low-risk group were significantly increased compared with that in high-risk group; patients with low-risk score have higher scores of immune function or pathways than those with high-risk, including cytolytic activity, human leukocyte antigen (HLA), major histocompatibility complex (MHC) class I, para-inflammation, type I interferon (IFN) response (Figures 9A,B). Besides, the associations of ubiquitination-based prognostic index with TME scores were explored. The correlation analysis demonstrated that the risk score was negatively associated with immune score derived from ESTIMATE algorithm. However, the risk score was not associated with ESTIMATE score, stromal score, and tumor purity score (Figures 9C–F).

As indicated by UALCAN database, the mRNA expression levels of HLA-A, TMEM129, UBE2D1, UBE2N, UBE2T in BCa tissue were significantly increased compared with that in normal tissue; however, the mRNA expression level of USP5 was similar between normal and tumor tissue. We performed univariable Cox regression analysis to validate the prognostic value of risk URGs in validating cohort (GEO cohort). The results showed that TMEM129, UBE2D1, and UBE2T were successfully validated to be associated with OS in GEO cohort. However, the prognostic value of HLA-A, UBE2T, and USP5 was not found to be associated with OS in GEO cohort. KEGG functional enrichment for high-risk group and low-risk group was performed using GSEA method. Toll like receptor signaling pathway, cytosolic DNA sensing pathway, other glycan degradation, and lysosome were major four KEGG enriched pathways in low-risk patients. GAP junction, DNA replication, oocyte meiosis, and RNA polymerase were major four KEGG enriched pathways in high-risk patients (Figure 10).