Development and validation of a hypoxia-related prognostic signature for breast cancer.

Hypoxia, an important component of the tumor microenvironment, plays a crucial role in the occurrence and progression of cancer. However, to the best of our knowledge, a systematic analysis of a hypoxia-related prognostic signature for breast cancer is lacking and is urgently required. Therefore, in the present study, RNA-seq data and clinical information were downloaded from The Cancer Genome Atlas (TCGA) and served as a discovery cohort. Cox proportional hazards regression analysis was performed to construct a 14-gene prognostic signature (PFKL, P4HA2, GRHPR, SDC3, PPP1R15A, SIAH2, NDRG1, BTG1, TPD52, MAFF, ISG20, LALBA, ERRFI1 and VHL). The hypoxia-related signature successfully predicted survival outcomes of the discovery cohort (P<0.001 for the TCGA dataset). Three independent Gene Expression Omnibus databases (GSE10886, GSE20685 and GSE96058) were used as validation cohorts to verify the value of the predictive signature (P=0.007 for GSE10886, P=0.021 for GSE20685, P<0.001 for GSE96058). In the present study, a robust predictive signature was developed for patients with breast cancer, and the findings revealed that the 14-gene hypoxia-related signature could serve as a potential prognostic biomarker for breast cancer. Copyright: © Wang et al.

Introduction

Breast cancer remains one of the most prevalent malignancies and leading causes of cancer-associated mortality for women globally, and it is estimated that there were approximately 2,000,000 new cases and 612,000 breast cancer-associated deaths in 2017 worldwide (1). Presently, the main treatments for breast cancer include systemic (chemotherapy, endocrine therapy, targeted therapy and immunotherapy) and locoregional (surgery and radiation therapy) therapeutic approaches (2). Therapeutic modalities and clinical prognosis mainly rely on classical pathological risk and prognostic factors, such as age, stage, tumor grade and tumor type (3,4). However, advancements in genetic profiling provide a precise molecular complement of each patient's tumor in order to better understand the heterogeneity of breast cancer (5).

Hypoxia, a significant microenvironmental factor, predicts a poor patient outcome due to the increased expression of related genes (6). Rapidly growing and highly invasive tumors exhibit hypoxic phenotypes as a result of insufficient and/or irregular blood perfusion (7). A previous study has demonstrated that almost 50% of locally advanced breast cancers exhibit hypoxia, resulting in resistance to irradiation and chemotherapeutic failure (8). Although these studies have confirmed the role of hypoxia in breast cancer, systematic studies of hypoxia-related gene expression and prediction of overall survival are rare (9).

In the present study, the differential expression of hypoxia-related genes in breast cancer samples was identified to investigate the biological functions of the enriched pathways. A hypoxia-related risk signature, including 14 hypoxia-related genes that can be used to accurately predict the clinical outcome of patients with breast cancer, was investigated. The data suggest that autophagy-related genes play a crucial role in breast cancer and therefore may be promising prognostic biomarkers and targets for breast cancer.

Materials and methods

Hypoxia-related gene set

The hypoxia-related genes were extracted from the hallmark hypoxic gene set in the Molecular Signatures Database v7.0 (10), which includes a list of 200 genes involved in hypoxic regulation pathways directly or indirectly based on the literature available in PubMed and biological public databases.

Gene expression data

The RNA-seq breast cancer data were downloaded from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov) on September 18, 2019. The RNA-seq data of 1,109 breast cancer samples and 113 adjacent normal samples were combined into a gene expression matrix using Strawberry perl (version 5.30.0.1).

Furthermore, clinical data from 1,097 cases were obtained from TCGA database. A total of 113 adjacent normal samples and 1,039 breast cancer patients who were followed-up for ≥1 month were screened for further analysis.

Three independent datasets were extracted from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) (accession nos. GSE10886, GSE20685 and GSE96058) to validate the prognostic signature of the genetic risk score (11–13). The series_matrix.txt format of the three datasets was obtained and all probes were transformed into gene symbols. The follow-up period was described in months, and samples of <1 month were excluded from the analysis.

Extraction of hypoxia-related gene matrix and identification of differentially expressed genes

TCGA gene expression matrix was selected to extract hypoxia-related gene expression patterns. The ‘limma’ (version 3.40.6) R package (14) was used to identify differentially expressed mRNAs between normal breast samples and breast cancer samples, with the threshold of |log2 fold-change (log2 FC)| >1 and false discovery rate (FDR) <0.05.

Functional enrichment analysis

To investigate the functional correlation of these differentially expressed hypoxic genes, the present study used the ‘clusterProfiler’ (version 3.12.0) R package (15) to perform Gene Ontology (GO) functional annotations and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis. Those with P- and Q-values <0.05 were considered as significant categories.

Establishment of the prognostic signature

Univariate Cox proportional hazard regression analysis was performed to identify the differentially expressed hypoxia-related genes associated with overall survival time. P<0.05 was considered to indicate a statistically significant difference. By weighting the Cox regression coefficients (β), the hypoxia-related risk signature was constructed with multivariate Cox regression. The risk score of each patient was then calculated according to the following formula: Risk score=β gene(1) × expression gene(1) + β gene(2) × expression gene(2) + ···+ β gene(n) × expression gene(n). The patients were classified into the high- and low-risk groups using the median hypoxia-related signature risk score as the cut-off point. The Kaplan-Meier (KM) method was used to evaluate survival differences. Multivariate Cox regression analysis and stratified analysis were used to estimate the prognostic value of the genetic risk score model for discovery cohorts. A Receiver Operating Characteristic (ROC) curve was used to identify the accuracy of model prediction.

Validation of the 14-gene hypoxia-related prognostic signature

To validate the robustness of the prognostic signature, the same coefficients from the training set were used to validate the validation sets, including the GSE10886 (n=106), GSE20685 (n=327) and GSE96058 (n=3,273) datasets. KM survival analysis was performed to determine the risk score of each patient and analyze the overall survival (OS) time between risk groups by the log-rank test.

Association between prognostic risk score and clinical symptoms

The associations between the high- and low-risk groups of breast cancer and six clinical symptoms (age, sex, stage, T stage, M stage, N stage and molecular subtypes) were validated using the ‘Beeswarm’ (https://cran.r-project.org/web/packages/beeswarm/index.html;version 0.2.3) package. P<0.05 was considered to indicate a statistically significant difference.

Results

Hypoxia-related phenotype expression in normal and breast cancer tissues

According to the cut-off criteria of |log2 FC| >1 and FDR <0.05, 26 (38.81%) upregulated and 41 (61.19%) downregulated hypoxia-related genes were identified between the normal and breast cancer tissues. A boxplot of the 67 differentially expressed hypoxia-related genes is presented in Fig. 1. The expression heatmap and volcano plot are presented in Figs. S1 and S2. Green or red color represents significantly downregulated and upregulated genes, respectively.

Figure 1.

Boxplot of 67 differentially expressed hypoxia-related genes. Green represents the adjacent normal sample, and red represents the breast cancer sample. N, normal; T, tumor.

GO and KEGG terms enrichment analysis

GO and KEGG pathway analyses for differentially expressed hypoxia-related genes were performed. The top 10 most significantly enriched GO and KEGG terms were presented in Fig. 2A and B, respectively. In the GO analysis, 498 significantly enriched GO categories were detected, such as ‘carbohydrate biosynthetic process’, ‘gluconeogenesis’ and ‘cellular carbohydrate metabolic process’ (Fig. 2A). As presented in Fig. 2B, 21 KEGG pathways with P- and Q-values <0.05 were identified, such as the ‘HIF-1 signaling pathway (hsa04066)’ and ‘glycolysis/gluconeogenesis (hsa00010)’.

Figure 2.

Top ten categories of functional enrichment analysis. The colors inside the circle represent the z-score, and red dots represent upregulated and blue dots represents downregulated. (A) GO enrichment analysis. (B) KEGG pathway enrichment analysis. GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Identification of hypoxia-related genes for prognosis by univariate Cox regression

A univariate Cox regression was performed to investigate the hypoxia-related genes associated with OS in TCGA dataset. A total of 24 genes (PFKL, P4HA2, CDKN1C, GRHPR, B4GALNT2, SDC3, PPP1R15A, SIAH2, NDRG1, STC2, TNFAIP3, IRS2, BTG1, PGK1, JUN, TPD52, MAFF, TIPARP, BHLHE40, ISG20, LALBA, ERRFI1, BCL2 and VHL) were identified to have a significant prognostic value in patients with breast cancer (P<0.05; Fig. 3).

Figure 3.

A total of 24 hypoxia-related genes were significantly associated with overall survival according to univariate analysis. Red represents an upregulated gene, and green represents a downregulated gene.

Establishment of prognostic signature by multivariate Cox regression analysis

A total of 14 hub genes (PFKL, P4HA2, GRHPR, SDC3, PPP1R15A, SIAH2, NDRG1, BTG1, TPD52, MAFF, ISG20, LALBA, ERRFI1 and VHL) were identified from the multivariate Cox regression (Table I), and the following genetic risk score signature was established: Genetic risk score = 0.3276 × PFKL + 0.3437 × P4HA2-0.3935 × GRHPR-0.2794 × SDC3 + 0.2226 × PPP1R15A-0.2481 × SIAH2 + 0.1751 × NDRG1-0.1977 × BTG1+ 0.1796 × TPD52-0.2649 × MAFF-0.3585 × ISG20 + 0.1624 × LALBA-0.3556 × ERRFI1+ 0.3039 × VHL.

Table I.

A total of 14 genes significantly associated with overall survival according to multivariate Cox regression analysis.

Gene symbol	Coefficient	Hazard ratio	95% confidence interval	P-value
PFKL	0.3276	1.3876	1.0298–1.8699	0.0313
P4HA2	0.3437	1.4101	1.0327–1.9256	0.0305
GRHPR	−0.3935	0.6746	0.4445–1.0237	0.0643
SDC3	−0.2794	0.7561	0.5971–0.9576	0.0203
PPP1R15A	0.2226	1.2494	0.9345–1.6704	0.1328
SIAH2	−0.2481	0.7802	0.6403–0.9507	0.0138
NDRG1	0.1751	1.1914	1.0198–1.3918	0.0272
BTG1	−0.1977	0.8205	0.6317–1.0657	0.1381
TPD52	0.1796	1.1967	0.9777–1.4649	0.0816
MAFF	−0.2649	0.7672	0.5427–1.0844	0.1334
ISG20	−0.3585	0.6986	0.5411–0.9020	0.0059
LALBA	0.1624	1.1763	1.0296–1.3438	0.0168
ERRFI1	−0.3556	0.7007	0.5567–0.8820	0.0024
VHL	0.3039	1.3551	1.0266–1.7886	0.0318

With the 14 hypoxia-related genes, the genetic risk score signature and survival curve were constructed (Figs. 4 and 5). ROC curves for 1-, 3- and 5-year survival rates were created to evaluate the predictive performance of the 14-gene risk signature (Fig. 6A). The area under the curve (AUC) was 0.728 (1-year), 0.726 (3-year) and 0.736 (5-year). These results indicated that the 14-gene risk signature had an ability to predict the survival of patients with breast cancer.

Figure 4.

Establishment of prognostic risk score models. (A) Risk score plot. (B) Heatmap of prognostic genes among 14 hypoxia-related genes. (C) Survival time and status for each patients with breast cancer.

Figure 5.

KM survival curves for the hypoxia-related gene risk score for breast cancer in TCGA dataset. The KM survival curves demonstrated that the OS time is longer in the low-risk group compared with in the high-risk group in TCGA dataset. KM, Kaplan-Meier; OS, overall survival; TCGA, The Cancer Genome Atlas.

Figure 6.

ROC and Kaplan-Meier curves of the risk signature based on the 14 characteristic genes. (A) ROC analysis of overall survival for the 14-gene signature in The Cancer Genome Atlas cohort. Based on median risk score, the 14-gene risk signature divided patients into high-risk and low-risk groups with distinct prognosis in the (B) GSE10886, (C) GSE20685 and (D) GSE96058 cohorts. ROC, Receiver Operating Characteristic; AUC, area under the curve.

The patients in three datasets (GSE10886, GSE20685 and GSE96058) were then categorized at the median risk score (high risk vs. low risk) based on the genetic risk score. As was expected, the survival curves revealed a significant difference between the high- and low-risk group in the validation datasets (GSE10886, GSE20685 and GSE96058; P=0.00738, P=0.0212 and P<0.001, respectively; Fig. 6B-D).

Association between prognostic risk signature and clinical symptoms

By investigating the association between clinical symptoms and risk signature in TCGA database, it was found that the risk score was significantly associated with age, stage, N stage and molecular subtype (P=0.022, P=0.005, P=0.041 and P<0.0001, respectively; Fig. 7). However, the risk score was not associated with sex, T stage and M stage (P=0.375; P=0.217 and P=0.158, respectively; Fig. S3).

Figure 7.

Association between prognostic risk signature and clinical symptoms. The association between prognostic risk signature and (A) age, (B) stage, (C) N stage and (D) molecular subtypes. HER-2, human epidermal growth factor receptor-2.

Discussion

Breast cancer is a progressive and complex disease with a high morbidity and mortality (16). Previously, numerous studies have demonstrated that hypoxic phenotypes are associated with malignant progression and the development of resistance to the treatment of breast cancer (17,18). While the majority of studies on hypoxia phenotypes have focused on a single gene, high-throughput biological technologies were used in the present study to predict the progression and outcomes of patients with breast cancer by detecting the expression of hypoxia-related genes systematically and comprehensively.

Functional analysis revealed that the HIF-1 signaling pathway and glycolysis/gluconeogenesis signaling pathway were enriched by the differentially expressed hypoxia-related genes. HIF-1α is involved in angiogenesis and energy metabolism, and in the activation of cancer stem cells during tumor initiation and progression (19–21). The ‘Warburg effect’ refers to the aerobic glycolysis of tumors, and is a main characteristic of tumors (22). Cancer-associated fibroblasts are the main components of the microenvironment of breast tumors and play a role in tumor metabolism (23). However, in breast cancer metabolism, breast cancer cells induce the glycolysis of fibroblasts, and the products are transferred to cancer cells for oxidative phosphorylation, which significantly improves the metabolic efficiency of cancer cells (24). This may be the reason for the rise in breast cancer mortality.

Among the 14 hypoxia-related genes, P4HA2, SIAH2, NDRG1, BTG1 and TPD52 have been reported to be associated with breast cancer in previous studies (25–33). P4HA2 participates in the process of extracellular matrix collagen deposition as the downstream activated element of HIF-α (25). Furthermore, it plays a crucial role in the remodeling of the surrounding environment of breast cancer, which contributes to the invasiveness of breast cancer (25). The silencing of P4HA2 has been shown to inhibit the proliferation of breast cancer cells (26). In a previous study by Toss et al (27), P4HA2 was demonstrated to be an independent factor of a poor prognosis for patients with ductal carcinoma in situ. NDRG1, another downstream activator of the HIF pathway, plays a pivotal role in tumor tolerance to oxygen concentration fluctuations (28). SIAH2 mediates the ubiquitination and degradation of NRF1, which affects the stability of mitochondrial-related gene expression in tumor cells, and thus promotes the heterogeneity of tumor metabolism (29). TPD52 is expressed in a variety of tumor cells, which makes it a target for the treatment of various types of cancer, such as pancreatic cancer, prostate cancer and cervical cancer (30–32). TPD52 is overexpressed in breast cancer and is a significant target gene for regulating breast cancer proliferation and invasion (33).

BTG1 is well understood as a gene that regulates the cell cycle and apoptosis. The overexpression of BTG1 may inhibit cell proliferation and induce cell apoptosis via the downregulation of Bcl-12. By contrast, the expression level of BTG1 in breast cancer cells has been shown to be markedly lower compared with that in normal breast cells (34). VHL has been confirmed in a number of studies as a prognostic gene for clear cell carcinoma of the kidneys, due to its potential ability to degrade HIF-1α (35,36). As a target gene of miRNA-155 in triple-negative breast cancer, VHL is involved in tumor angiogenesis and is associated with a poor prognosis (37). ISG20, another potential marker of clear cell carcinoma of the kidneys, has also been reported to be associated with breast cancer (38). ISG20, with exonuclease activity, has also been found to promote tumor angiogenesis (38). LALBA, ERRFI1, PFKL and SDC3 are novel prognostic markers for breast cancer.

Further analyses demonstrated that the 14-gene hypoxia-related model gene signature could accurately predict the clinical prognosis of patients with breast cancer, which was further validated in three independent datasets (GSE10886, GSE20685 and GSE96058) in the GEO database. The expression of hypoxia-related genes was increased in the high-risk groups, suggesting that the level of hypoxia was increased in the high-risk groups with poor prognostic outcomes.

Studies have confirmed that the tolerance of MCF-7 cells to doxorubicin may be due to hypoxia-induced MDR1 expression under hypoxic conditions, which is the most widely studied mechanism for the cancer cells being resistant to chemotherapy (39,40). In addition, the inhibition of HIF-1 has been shown to downregulate the expression of MDR1 (41). These are also closely associated with a poor therapeutic effect; thus, the recognition of the hypoxia level is crucial for the prognosis of patients with breast cancer.

To investigate the association between hypoxia-related genes and the prognosis of patients, it was found that the hypoxia-related gene expression risk signature could independently predict the overall survival rate of patients with breast cancer. By performing multivariate analysis for survival, the 14-gene hypoxia-related signature in the present study was demonstrated to serve as a potential biomarker for the classification of the prognosis of patients with breast cancer, and may also guide the therapeutic strategies for these patients.

However, the current study was limited by the absence of experimental evidence. Further studies investigating cell lines and involving animal experiments will be needed to confirm the present findings and delineate the molecular mechanism. In conclusion, the current study demonstrated the prognostic predictive value of this 14-gene signature for breast cancer.