Estrogen receptor—positive breast cancer survival prediction and analysis of resistance–related genes introduction

Data download

The raw counts of the RNA sequence data of the TCGA-BRCA dataset were downloaded from the TCGA website (https://portal.gdc.cancer.gov/). The clinical features and the disease-free survival (DFS) status of these patients were downloaded from the cBioportal database (http://www.cbioportal.org). We selected 473 ER+ and HER2- breast cancer patients in stage I-III or the subsequent analysis (Table 1). The R software (version 3. 6. 3) was used to standardize and process data.

GSE21653 (100 samples of ER+ and HER2- early invasive breast cancer), GSE17705 (298 samples of ER+ breast cancer) and GSE11121 (200 patients) datasets were acquired from GEO database (https://www.ncbi.nlm.nih.gov/geo/) as the validation cohorts. GSE21653 has DFS as the endpoint, GSE17705 has relapse-free survival (RFS) as the endpoint and GSE11121 has distant metastases-free survival (DMFS) as the endpoint.

Screening of differentially expressed genes

The resistance-related differentially expressed genes (DEGs) were detected with the “limma” package in R software from GSE67916 (DEGs between MCF7 tamoxifen resistant (MCF7/TamR) cell lines and MCF7 parental cell lines) and GSE24460 (DEGs between MCF7 doxorubicin resistant (MCF/ADR) cell lines and MCF7 parental cell lines). Defined genes with adjusted p values <0. 05 and —log2foldchange— values >1. 5 as DEGs. The genes both up-regulated (16 genes) and down-regulated (21 genes) in the two GEO datasets were finally screened (Fig. S1B). Volcano plots and heatmaps were performed using the R packages “ggpubr”, “ggthemes” and “pheatmap” to exhibit the differentially expressed genes (Figs. S1A, Fig. 1C).

Identification and validation of a risk signature associated with breast cancer

The data downloaded from TCGA datasets were used as the training cohort to analyze the correlation of genes’ expression with the DFS time based on univariate Cox regression analysis. Patients without survival for less than 0 days were removed to identify prognostic genes, only the genes with the p values < 0. 05 were selected as candidates (Fig. 2A). Then, the candidate genes were selected for stepwise multivariate Cox regression analysis (Fig. 2B). The Akaike information criterion (AIC) was used to avoid over-fitting, we chosed the best-fit predictive model with the lowest AIC (AIC =359. 44). Finally, a four-gene risk signature was established and calculated their risk coefficients use the formula: The riskScore = (−0. 175*expression of AMIGO2) + (−0. 210* expression of LGALS3BP) + (−0. 106* expression of SCUBE2) + (−0. 293* expression of WLS). Using the optimal riskScore cut-off value (−10.09) calculated by the R package “survminer”, we divided patients in both the training cohort and the validation cohort into high-risk and low-risk group. Then, Kaplan–Meier curve and log-rank test were used to analyze the different survival time between high-risk group and low-risk group. Subsequently, a time-dependent receiver operating characteristic (ROC) curve was calculated to compare the sensitivity and specificity of risk factors.

Development of the prognostic nomogram

To identify independent predictors of DFS, we tested clinical risk factors and the gene-based riskScore. Then we used the R package “rms” to establish the prognostic nomogram. The discrimination and calibration abilities of the nomogram were assessed with a concordance index (C-index) and calibration plot. For further evaluation of the predictive performance, we also used the time-dependent ROC curves.

Function and pathway enrichment analyses

The R package “clusterProfiler” was used to classify the 37 DEGs’ GO terms or KEGG pathways. Functional enrichment analyses were carried out for GO terms and KEGG pathways through the significance threshold of p value < 0. 05.

Immune infiltration analyses

The 22 kinds of immune cell infiltration levels of TCGA-BRCA dataset were quantified by CIBERSORT algorithm with 1,000 permutations (http://cibersort.stanford.edu/) (Chen et al., 2018a). CIBERSORT calculates a p value to confirm the accuracy of the results, only the samples with the p value < 0.05 were selected for analyzing the immune infiltration levels between the high- and low- risk groups. Moreover, to explore the correlation between the four genes’ expression and the abundance of immune cell infiltrates (B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages, and dendritic cells), the Tumor Immune Estimation Resource (TIMER, https://cistrome.shinyapps.io/timer/) was utilized.

Statistical analyses

All the statistical analyses were performed with R software (version 3. 6. 3). All the statistical analyses included univariate and multivariate Cox regression analysis, ROC curve analysis, Kaplan–Meier survival analysis and correlation analysis. The correlation analysis was based on the R package “corrplot”. Statistical differences were compared by Wilcoxon test between two groups and Kruskal-Wallis H test for multigroup comparison. P < 0.05 was considered statistically significant.

Screening DFS-related resistance-related DEGs

The flow chart was shown in Fig. 1. To screen the DFS-related resistance-related DEGs, two datasets was downloaded from GEO database. After “limma” filtering between MCF7 doxorubicin or tamoxifen resistant cell lines and MCF parental cell lines (Fig. S1A), the intersection of up-regulated and down-regulated genes was taken respectively (Fig. S1B). A total of 37 resistance-related DEGs were obtained. The expression profile of the 37 DEGs in GSE24460 and GSE67916 were shown in Fig. S1C.

Next, we investigated whether genes could be used as prognostic candidates in ER positive (ER+) and HER2 negative (HER2-) breast cancer patients in stage I-III. The univariate Cox regression analysis found that eight genes were associated with DFS time in ER+ and HER2- breast cancer, all of these genes were protective factors (Fig. 2A). The eight candidates were further analyzed by multivariate Cox regression analysis and examined four genes (AMIGO2, LGALS3BP, SCUBE2 and WLS) as the prognostic candidates (Table 2, Fig. 2B). In addition, we verified the prognosis of these four genes in the Kaplan–Meier Plotter website (https://kmplot.com/analysis/). Consistent with our results, WLS and SCUBE2 were protective factors, high expression of the two genes had longer RFS time (Fig. S2C). The results of other two genes were not statistically significant. A risk score formula (the formula shown in methods) derived from the expression levels of the four prognostic candidates was constructed by the multivariate Cox regression analysis.

Identification and validation of a four-gene based risk signature

All enrolled patients were divided into high- and low-risk groups with its optimal cut-off value among the training cohort by using the riskScore. The survival status plot and heatmap for the expression of the four genes were showed in Figs. 2C–2D. The prognostic capacity of the four-gene based riskScore was estimated by calculating the AUC of a time-dependent ROC curve (Fig. 3A). The AUCs of the four-gene risk signature were 0.686 and 0.714 for 3- and 5-year DFS respectively, suggesting that the risk signature had certain sensitivity and specificity. Kaplan–Meier curve showed that patients in the high-risk group had a lower DFS time than those in the low-risk group (Fig. 3A).

We also analyzed the correlation of the four genes in ER+ and HER2- breast cancer patients in TCGA-BRCA dataset. The correlation heatmap showed that the maximum R values was 0.14, indicating that these four genes were independent of each other (Fig. 2F).

To validate the prognostic prediction performance of riskScore, the same as training cohort, patients in the validation cohorts were classified into high- and low-risk groups using the same formula obtained from the training cohort. In all validation cohorts, patients with a high riskScore had a shorter survival time and a high riskScore was associated with worse prognosis (Figs. 3B–3D). Additionally, the time-dependent ROC analyses for the survival prediction of the risk signature showed AUCs of 0. 624 at 3 years, 0. 613 at 5 years in GSE21653 and 0. 618 at 3 years, 0. 616 at 5 years in GSE17705, demonstrating that this risk signature was capable of predicting survival in ER+ and HER2- or ER+ and HER2 equivocal patients (Figs. 3B–3C). Moreover, we also used the time-dependent ROC analysis to explore the prognostic prediction performance if the subtype of breast cancer was unknown. The AUCs were 0.730 at 3 years and 0.669 at 5 years in GSE11121 (Fig. 3D), suggesting that our risk signature also has good predictive ability when breast cancer subtype is unknown.

Stratified analyses with clinical factors

As shown in Fig. 4, the four-gene based risk signature could serve as a promising candidate for predicting the survival of breast cancer in different subtypes stratified by the clinical factors, including N0-1 type (p = 0.0014), N2-3 (p = 0.026), T1-2 (p = 0.0006) T3-4 (p = 0.04), stage I-II (p = 0.0022), stage III (p = 0.0078), PR positive (PR+, p < 0.0001) patients, respectively. But the PR negative (PR-) subtype was not statistically significant.

As shown in Fig. 5, the riskScore value was higher expression in PR- subtype than in PR+ subtype, SCUBE2 and WLS were higher expression in PR+ subtype than in PR-subtype, the expression of SCUBE2 was negatively correlated with T and Stage, these results were consistent with the clinical prognostic significance.

Functional enrichment analyses

To explain the function of the resistance-related DEGs, we performed enrichment analyses of the GO terms and KEGG pathways (Figs. S2A–S2B).

GO analyses showed that these resistance-related DEGs were enriched in several biological processes (BP), including antigen processing and presentation of peptide antigen via MHC class I, antigen processing and presentation of exogenous peptide antigen via MHC class I and epithelial structure maintenance (Fig. S2A).

KEGG analysis showed that the resistance-related DEGs were related to endocytosis, cellular senescence, antigen processing and presentation and adherens junction (Fig. S2B).

Revealing the correlation between the risk signature and tumor-infiltrating immune cells

Because TILs play a key role in chemotherapy response and our enrichment analyses results showed that the DEGs associated with antigen presentation, first, we evaluated the abundance of the 22 kinds of immune cells in each TCGA-BRCA patient (Fig. 6A). Next, we utilized the Wilcoxon test to compare the difference immune cell infiltration levels between the high- and low groups and found that B cells memory (p = 0.042) and T cells follicular helper (p = 0.016) were high infiltrated in the high-risk group, and B cells naive (p = 0.047) and mast cells resting (p = 0.042) were low infiltrated in the high-risk group (Fig. 6B). Unfortunately, due to the limitation of sample size, the relationship between the immune cell infiltration levels and survival is unknown. Whereafter, we assessed the correlation between the four resistance-related DEGs and six infiltrating immune cells by using the TIMER database. The results revealed that the genes had good correlation with immune cells (Fig. S3).

Building a prognostic nomogram combining riskScore with clinical risk factors

Nomograms are widely used in predicting cancer patients’ prognoses, principally because they have the ability to decrease the statistical prediction models into a single numerical assessment, which is tailored to the profile of an individual patient (Iasonos et al., 2008; Long et al., 2018).

In this article, we combined the four-gene based riskScore with clinical risk factors to construct a nomogram to assess the probability of 3- and 5- year DFS for patients with ER+ and HER2- early stage invasive breast cancer in TCGA-BRCA dataset. After univariate and multivariate Cox regression analyses filtering, the riskScore, T and Stage remained in the final Cox model for DFS (Figs. 7A–7B). A nomogram was established based on the final Cox model for DFS to evaluate of the 3- or 5-year DFS (Fig. 7C). The clinical stage and riskScore of the patient can be matched with this Nomogram to obtain the total points and estimate the survival rate of the patient. The nomogram achieved a C-index of 0. 740 (95% CI, 0. 653–0. 826), and the calibration plots demonstrated that the nomogram has good predictive power (Fig. 7D). Time-dependent ROC analyses also indicated that the nomogram was accurate in prediction of DFS compared to use riskScore or clinical risk factors alone (Fig. 7E). In conclusion, all the results suggest the clinical practical use of the nomogram for prediction of DFS.