Identification of a six-gene signature associated with tumor mutation burden for predicting prognosis in patients with invasive breast carcinoma.

BACKGROUND: Breast cancer (BC) is one of the most common cancers with high mortality worldwide. In the present study, through bioinformatics analysis, we aimed to identify new biomarkers to predict the survival rate of BC patients.
METHODS: Differentially expressed genes (DEGs) between low- and high-tumor mutation burden (TMB) groups were identified by using The Cancer Genome Atlas (TCGA) dataset and integrated analysis. Gene Ontology (GO), Kyoto encyclopedia of genes and genomes (KEGG) enrichment analysis, and the protein-protein interaction (PPI) network, were applied to predict the function of these above DEGs. Then, the Cox proportional hazard model was developed to screen DEGs. Based on the prognostic signature, survival analysis was used on The Cancer Genome Atlas Breast Invasive Carcinoma (TCGA-BRCA) dataset. Finally, the single-sample gene set enrichment (ssGSEA) analysis was employed to estimate immune cells related to this signature.
RESULTS: To create a prognostic signature, 6 DEGs were identified. The results revealed that the survival time of patients with high-risk scores based on the expression of the six-gene signature was dramatically shorter than that of patients with low-risk scores in BC. Furthermore, survival analysis and multivariate cox analysis indicated that the six-gene signature was an independent prognostic factor of BC. Then, we built a nomogram that integrated the clinicopathological factors with the six-gene signature to predict the survival probability of BC patients. We eventually predicted the 20 most vital small molecule drugs by CMap, and Nadolol was considered as the most promising small molecule to treat BC. Moreover, ssGSEA analysis showed that the 6 genes were closely associated with immune cells.
CONCLUSIONS: We constructed a six-gene signature associated with TMB that can improve the prognosis prediction and could be seen as a biomarker for BC patients. 2020 Annals of Translational Medicine. All rights reserved.

Introduction

Breast cancer (BC) is one of the major malignancies among females, and mortality remains the second main reason for cancer deaths worldwide (1). In recent years, the diagnosis of BC mainly depends on pathological tests, imaging examinations, and assessments of tumor markers (2). Because of a high recurrence rate of BC, the age of BC onset has gradually become younger. Therefore, the early detection and diagnosis of BC play an essential role in preventing the genesis and development of BC (3). It is thus necessary to identify possible molecular biomarkers for the diagnosis, prognosis, and survival risk prediction for BC.

Tumor mutation burden (TMB) is a new biomarker for immune therapy and is defined as the number of somatic, coding, base substitution, and indel mutations per megabase (4). It is widely acknowledged that the TMB plays a vital role in the occurrence and development of cancer (5). In one study, cancer patients with high levels of TMB showed stronger immunotherapeutic responses than those with low levels of TMB (6). Exploring gene expressions that are closely correlated with TMB levels can, therefore, be of great importance, and TMB has also been linked to cancer prognosis (7). Thus, there is an urgent need to identify potential TMB-related gene signature for BC prognosis.

In the current study, we firstly identified DEGs based on TMB by RNA sequencing data from The Cancer Genome Atlas (TCGA). Then, functional enrichment analyses were further employed to analyze the biological roles of these DEGs. Finally, we constructed a six-gene signature for BC by using TCGA datasets and bioinformatics analysis and verified that the prognostic value of this signature was associated with some clinical factors.

Methods

Data source

We obtained the RNA-seq datasets with the clinic information of BC patients comprising 474 tumor samples and 481 normal samples from the The Cancer Genome Atlas Breast Invasive Carcinoma (TCGA-BRCA) datasets and analyzed them using the limma package (8) with the false discovery rate (FDR) <0.05 and |log(fold change) >1|. Subsequently, we downloaded the masked somatic mutation data from TCGA database. Then, we classified BC patients into the low-TMB subtype and the high-TMB subtype by using the median of TMB data as the cutoff value and compared the expression levels of TMB-related genes between both subtypes.

Functional enrichment analysis and protein-protein interaction (PPI) network analysis of differentially expressed genes (DEGs)

To investigate potential biological processes (BP), cellular components (CC), molecular functions (MF), and signaling pathways related to the DEGs, we applied the Gene Ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) pathway analyses by clusterProfiler package (9) in R. The potential connections between the DEGs were then visualised using the STRING database (10).

Cox proportional hazards model construction and model validation

By combining these prognostic gene expression markers with the regression coefficient (β) from the multivariate Cox proportional hazards regression analysis, a risk score model was established based on the results of multivariate Cox regression analysis. Then, the formula of risk score was calculated as follows: expression of gene1 × β1gene1 + expression of gene2 × β2gene2 + expression of genen × βngenen. According to the median risk score, all patients were assigned into the high-risk and low-risk groups. Proportional assumptions for Cox proportional hazard model were examined by Kaplan-Meier analysis. Receiver operating characteristic (ROC) and the area under the curve (AUC) were plotted by using R package survival ROC. In addition, the distribution of risk score and survival status, along with the gene expression levels of each patient, was an­alyzed by using the R software.

Correlation analysis of immune infiltration analysis

The CIBERSORT method (11) was used to estimate immunocyte population fractions, which could sensitively and accurately distinguish 22 phenotypes of immune cells. Then, in subsequent review, those with the CIBERSORT P<0.05 were considered eligible. Moreover, we employed 29 immune-related gene sets, which represent different immune cell types, functions, and pathways. We performed the single-sample gene set enrichment analysis (ssGSEA) to calculate the activity or enrichment levels of immune cells, functions, or pathways in tumor samples.

Identification of small molecular drugs

We compared the six-gene signature with those in the Connectivity Map (cMap) database (12) to identify candidate small molecular drugs associated with BC. We then divided DEGs into the upregulated group and downregulated group. Next, we used these genes to query the cMap database. A connectivity score from -1 to 1 indicated the stimulation of the compound, while a negative score indicated the inhibition of the compound.

Statistical analysis

All the statistical analyses were performed using the R package. Survival curves were established by the Kaplan-Meier method and the log-rank test. Univariate and multivariate Cox proportional-hazard models were used to evaluate the hazard ratios of prognostic factors. We applied the ‘RMS’ package in R to plot the nomogram and calibration analysis. A P value <0.05 was viewed as statistically significant.

Results

Identification of DEGs between two TMB groups

A total of 343 DEGs were screened between low/high-TMB groups from the TCGA-BRCA dataset according to the screening criteria, including 168 down-regulated and 175 up-regulated genes (). Then, GO and KEGG analyses were applied to investigate the biological roles of the DEGs. The functional enrichment analysis revealed that these DEGs were mostly enriched in epidermis development, extracellular matrix, and receptor-ligand activity among BP, CC, and MF, respectively (). KEGG analysis demonstrated that these DEGs primarily participated in protein digestion and absorption ().

Figure 1

Differentially expressed TMB-related genes and enrichment analysis results of DEGs. (A) The heatmap of DEGs between the high‐TMB and low‐TMB groups in BRCA; (B,C) functional enrichment pathway analysis of TMB-related DEGs in BRCA, including biological processes (BP), cellular components (CC), molecular function (MF), and KEGG pathway; (D) comparisons of 22 immune cell infiltrations between the high‐TMB and low‐TMB groups.

To further explore the relationship between high and low levels of TMB and infiltration of immune cells, The CIBERSORT algorithm was used to measure the relative proportions of immune cells in low/high TMB groups. It revealed that the abundance levels of naive B cells, resting CD4 memory T cells, resting dendritic cells, and resting mast cells were significantly lower in the high-TMB group than the low-TMB group. Conversely, T cells follicular helper, Tregs, and macrophage subsets (M0/M1) were significantly higher in the high-TMB group (, P<0.05), indicating that TMB levels can affect the cancer progression. Moreover, the PPI network of the DEGs was further constructed by the STRING database with 273 nodes and 448 edges ().

Figure 2

The protein-protein interaction network (PPI) analysis of the DEGs. The PPI network among the DEGs was established by using the STRING database.

Construction of the Cox regression model and survival analysis

Multivariable Cox proportional hazards stepwise regression with backward selection was used to build a prognostic model. Then, we calculated the risk score for each patient in the datasets. Hence, we established a six-gene signature as follows: risk score = (‒0.2224 × IGFALS) + (0.2232 × PAX7) + (0.1243 × SPDYC) + (‒0.1224 × IGHA2) + (‒0.1172 × SERPINA1) + (‒0.3558 × ADRB1). Then, the distribution of risk score, survival status, and the corresponding heatmap of the gene expression levels of each patient in the signature was analyzed and plotted (). In addition, the survival risk curve showed that patients in the high-risk group had a significantly shorter overall survival time compared with patients in the low-risk group (). ROC analysis revealed reliable performance in the survival prediction of the model, and the AUC was 0.705 (). Moreover, survival analysis using Kaplan-Meier showed that these 6 genes were closely associated with the overall survival of BC patients (, P<0.05).

Figure 3

Prognostic analysis based on risk score model of the 6 genes. (A) Patients’ survival status distribution by the risk score; (B) patient survival status distribution of the low-risk group and the high-risk group; (C) heatmap of the 6 genes for low- and high-risk groups; (D) Kaplan-Meier curves for the low- and high-risk groups; (E) the receiver operating characteristic (ROC) curve validation of prognostic value by the risk score.

Figure 4

Evaluation of the prognostic value of the six-gene signature using Kaplan-Meier analysis. (A) ADRB1, (B) IGFALS, (C) IGHA2, (D) PAX7, (E) SERPINA1, (F) and SPDYC.

Association between the risk scores of the six-gene signature and clinicopathological features

For observing the association between the risk scores and each clinicopathological factors, we separately calculated the distribution of risk scores. We examined the expression of the 6 genes in high/low-TMB and high/low-risk patients in the TCGA-BRCA dataset (). We observed significant differences between the high- and low-risk groups concerning tumor grade (T, P<0.001) and lymph node (N, P<0.05). Furthermore, we found risk scores were significantly different between patients stratified by T grade (P<0.05), age (P<0.001), and status (P<0.001) (). Moreover, univariate and multivariate Cox analysis was applied to determine whether the six-gene signature can be employed to independently predict the survival of BC patients. The results showed that the six-gene risk score could be used to independently predict the overall survival of the patient (P<0.001, , ). Then, we built a nomogram that integrated the clinicopathological factors with the six-gene signature to predict the survival probability of BC patients (). Meanwhile, calibration plots for the probability of 3- and 5-year survival were consistent between predicted values by the nomogram and measured values ().

Figure 5

The relationship between the six-gene risk score and clinical information. (A) The heatmap showed the expression levels of the 6 TMB-related genes in the low- and high-TMB groups; (B) the heatmap indicated the expression levels of the 6 genes in the low- and high-risk groups. The distribution of clinicopathological features was compared between the low/high-TMB and low/high-risk groups. (C,D) forest plot of the univariate and multivariate Cox analysis for the independent six-gene signature; (E) a nomogram was used to predict the overall survival at 1 year, 3 years, and 5 years, with a six-gene risk score and clinical information. (F,G) the calibration plots for predicting overall survival at 3 years and 5 years. *, P<0.05; ***, P<0.001.

Table 1

Cox proportional hazards model analysis of prognostic factors

Variables	Univariate analysis				Multivariate analysis
	HR	95% CI	P value		HR	95% CI	P value
Six-gene risk score	1.453	1.314–1.605	2.40E–13*		1.424	1.276–1.589	2.61E−10*
Age (<60 vs. ≥60)	1.033	1.017–1.049	3.69E–05*		1.032	1.015–1.049	0.0001*
TNM stage (III + IV vs. I + II)	2.123	1.653–2.727	3.71E–09*		1.700	0.951–3.040	0.0733*
T (tumor)	1.555	1.235–1.956	0.0001*		1.030	0.735–1.441	0.862
N (lymph node)	1.707	1.396–2.088	1.84E–07*		1.174	0.847–1.628	0.334
M (metastasis)	5.821	3.104–10.914	3.96E–08*		1.116	0.451–2.760	0.810

*, P value <0.05 was considered statistically significant.

To further explore the relationship between the 6 TMB-related genes and the infiltration of immune cells, we conducted the ssGSEA algorithm to estimate the coefficient of the association of 6 genes in immune cell subsets. It was found that these 6 genes were closely associated with tumor-infiltrating immune cell subset, indicating that these genes mainly participate in immune response (, P<0.05). Furthermore, we applied the cMap database to identify potential therapeutic drugs for BC, and the top 20 most distinct small molecule drugs with their enrichment scores were identified (). Nadolol was related to highly significant negative scores (‒0.949). The potential small molecules could reverse the expression levels of the 6 TMB-related genes induced by BC, and represent new directions for research in BC treatment.

Figure 6

The 6 TMB-related prognostic genes were associated with immune-cell subset. Blue boxes indicate positive correlation; red boxes indicate negative correlation. *, P<0.05; **, P<0.01.

Table 2

List of the 20 most significant small molecule drugs that can reverse the tumoral status of BRCA

CMap name	Enrichment	P value
Nadolol	−0.949	0.00046
Finasteride	−0.743	0.00056
Gabexate	−0.859	0.0007
Primidone	−0.856	0.00074
Metronidazole	−0.751	0.00176
Metrifonate	0.757	0.002
Omeprazole	0.806	0.00265
Acenocoumarol	−0.735	0.0027
Medrysone	0.668	0.00371
Biperiden	−0.719	0.00376
Phensuximide	−0.791	0.00384
Quinidine	0.869	0.00409
Bupivacaine	0.783	0.00414
Prestwick-860	0.775	0.00483
Rilmenidine	−0.771	0.00559
Gentamicin	−0.767	0.00589
Betahistine	0.758	0.00647
Diphemanil metilsulfate	−0.684	0.00721
Doxylamine	0.683	0.00813
Verteporfin	0.834	0.00929

Discussion

BC is a heterogeneous disease with different biological properties caused by a variety of factors involving the accumulation of gene changes (13). Despite the improved diagnosis and therapy of BC, the prognosis is still poor (14). Finding accurate biomarkers for early diagnosis and a more accurate prognosis of BC will help to improve the efficiency of therapies for BC and provide molecular markers for targeted treatment (15). Thus, it is urgent to develop more sensitive and specific biomarkers for the prognosis of BC patients.

Currently, many studies are aimed at finding predictive biomarkers of immune responses. (16). TMB has been used to reveal the process of mutation accumulation in tumors and has proven to be effective in predicting the immune response in various cancers, including BC (17). Previous studies have reported that TMB is closely associated with the prognosis of cancer patients (18). In this study, a total of 343 TMB-related DEGs also showed a considerable difference between the 2 groups with dif­ferent TMB levels, which may be of great significance. To demonstrate the molecular function of these DEGs, we applied the GO and KEGG analyses. The results indicated that DEGs were mainly enriched in epidermis development, extracellular matrix, and receptor-ligand activity, in addition to protein digestion and absorption. Moreover, when we compared the differential abundance of immune cells between the low and high TMB groups via CIBERSORT algorithm, the results revealed that the abundance levels of naive B cells, resting CD4 memory T cells, resting dendritic cells, and resting mast cells were significantly lower in the high-TMB group compared to those in the low-TMB group. Similarly, a previous study showed that increased overall survival was related to relatively larger fractions of resting CD4 memory T cells (19).

Moreover, using univariate and stepwise multivariate Cox regression analyses, we established a prognostic signature for BC based on the expression of 6 genes, including IGFALS, PAX7, SPDYC, IGHA2, SERPINA1, and ADRB1. Among these, IGHA2 and SERPINA1 are considered risky prognostic genes in BC. Kang et al. reported that the differential expression of IGHA2 in the early and advanced stages of BC indicated that IGHA2 might be a novel marker for the progression of BC (20). Chan et al. showed that SERPINA1 expression is correlated with the outcome of ER+ and ER+/HER2+ BC (21). Moreover, Kang et al. reported that knockdown of ADRB can suppress tumor progression by inhibiting tumor growth, metastasis, and angiogenesis in BC, and that ADRB1 was abundantly expressed in several types of BC cells (22). As a member of the Spy1/RINGO family, SPDYC might also play a vital role in regulating both normal and abnormal growth processes in the breast (23).

Furthermore, some previous findings have shown that insulin-like growth factor binding protein (IGFALS) can act as a possible tumor-suppressor gene silenced by methylation in hepatocellular carcinoma (24). Also, He et al. showed that deregulated levels of PAX7 can contribute to muscle wasting in cancer cachexia (25). Previous studies also have suggested that SERPINA1 could inhibit caspase-3-mediated apoptosis and has been proven to predict the development of cutaneous squamous cell carcinoma (26). Finally, SERPINA1 has been shown to contribute to the development of metastatic features in epithelial ovarian cancer cells (27).

Through bioinformatics analyses, based on the six-gene signature, we found that BC patients with high-risk scores had a significantly decreased survival time compared with low-risk scores. ROC analysis demonstrated the relatively high sensitivity and specificity of the six-gene signature model. Similarly, in this study, the survival analysis of the six genes based on a prognostic signature all demonstrated a significant difference between the high and low expression groups. Then, we integrated the six-gene signature with other independent clinical variables to construct a comprehensive model for monitoring progression in BC, which suggested that the six-gene risk score could be used independently to estimate the overall survival of the patients. Based on the role of TMB in immune responses (28), ssGSEA analysis revealed that the 6 genes are closely associated with tumor-infiltrating immune cells. Nadolol was additionally identified as an important small molecular agent in BC’s development.

Conclusions

In this study, we integrated TMB with gene expression and built a prognostic prediction model that can improve survival prediction for BC patients and reveal the TMB-based prognostic signature which may be used as a prognostic biomarker for BC patients. Further experimental studies are still needed to confirm the established signature.

The article’s supplementary files as