Identification and validation of tumor microenvironment-related prognostic biomarkers in breast cancer.

BACKGROUND: Stromal cells and immune cells in tumor microenvironment (TME) have been reported to have significant value in the diagnosis and prognosis of cancers. We aimed to identify key biomarkers predicting survival in the TME of breast cancer.
METHODS: Cell type enrichment analysis was performed to estimate cell scores using the xCell method with gene expression data from public database. Least absolute shrinkage and selection operator (LASSO) regression was used to identify key signature from the cell scores.
RESULTS: Totally, 50 cells in TME had different scores between 1,078 breast cancer tissues and 112 adjacent normal tissues. We identified a 4-cell signature predicting breast cancer survival, including myocytes, natural killer T cell (NKT), conventional dendritic cell (cDC) and sebocytes, which was validated in the test set. Further analysis showed that cDC score was a key signature predicting prognosis of breast cancer. cDC score was significantly associated with molecular classification and stage of breast cancer, as well as expression level of Ki67. Spearman's correlation analysis found that cDC score was inversely correlated with the expression level of HER2. High cDC score may predicate better pathological complete response rate. Mechanism analysis indicated high cDC score was associated with elevated immune activity; IL-2 was a key gene associated with high cDC score; and Breast cancer patients with high IL-2 expression had a longer survival time.
CONCLUSIONS: In conclusion, cDC score was a key signature predicting prognosis for breast cancer. cDCs may exert antitumor effects by upregulating IL-2. 2021 Translational Cancer Research. All rights reserved.

Introduction

Breast cancer is commonest malignancy in women worldwide. Only in the United States, more than 250 thousand new cases are diagnosed with breast cancer each year and 66 thousand cases die from this malignancy (1). In China, it is estimated 208 thousand new cases and 55.5 thousand deaths in 2010 (2). Although significant progress has been made over the past decade in understanding of cancer biology, identification of risk factors, and new treatments of breast cancer, it is still insufficient in terms of early diagnosis, therapy, prevention and prediction of prognosis (3). Therefore, it is needed to find new markers for breast cancer.

The “seed and soil” hypothesis postulates that the tumor microenvironment provides fertile soil for the growth of tumor cells (4,5). Emerging evidence indicates that the cross-talk between tumor cells and tumor microenvironment exerts important effects on initiation, progression and metastasis of tumor (6). For example, tumors, as key drivers, control the differentiation of precursors of cancer-associated fibroblasts by secreting factors. Once present in the developing tumor, cancer-associated fibroblasts shape the tumor microenvironment to support tumor cell survival, dissemination, immune suppression, angiogenesis, and therapy resistance (7). Tumor microenvironment is a complex ecosystem of stromal cells and immune cells (8). Stromal cells and immune cells have been reported to have significant value in the diagnosis and prognosis of various cancers including breast cancer (9).

In the present study, we used xCell to perform cell type enrichment analysis from gene expression data (10). xCell is a gene signatures-based method learned from thousands of pure cell types from various sources and infers 64 immune and stromal cell types. xCell signatures have been validated extensively, and were shown to outperform other methods. We predicted cell scores using xCell with the data from The Cancer Genome Atlas (TCGA), Molecular Taxonomy of Breast Cancer International Consortium (METABRIC) and Gene Expression Omnibus (GEO). Then, we selected and validated prognostic biomarkers from cell scores for breast cancer using the least absolute shrinkage and selection operator (LASSO) (11). We present the following article in accordance with the TRIPOD reporting checklist (available at https://dx.doi.org/10.21037/tcr-21-1248).

Methods

Data acquisition

We downloaded RNA-Seq data and clinical data for 1,097 female breast cancer patients from the data portal for TCGA (accessed October 2020) (12). The data of 1,904 breast cancer patients were obtained from METABRIC database. Datasets of GSE96058, GSE20194, GSE22358, GSE25066 and GSE32646 were downloaded from GEO. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Cell type estimation

Cell type enrichment analysis was performed using the xCell (R package version 1.1) (10) with gene expression data. We compared the cell scores of tumors to those of adjacent normal tissues. Violin plot was drawn using vioplot package (version 0.3.5) in R software. The difference between two groups in cell scores was assessed using Wilcoxon test, and a P value less than 0.05 was considered significant.

Analysis of the cell score according to breast cancer pathological features

We investigated the relationships between cell scores and pathological features of breast cancer, including molecular classification [luminal (HER2−/ER+), HER2+(HR+), HER2+(HR−), and HER2−/HR−], stage (stage 1, stage 2, stage 3 and stage 4) and expression level of Ki67. R package of “beeswarm” (version 0.2.3) was used to illustrate the distribution of cell score in different groups.

Integration of protein-protein interaction (PPI) network

The PPI network was constructed as previously described (13). Briefly, we mapped a set of genes to Search Tool for the Retrieval of Interacting Genes (STRING) (14). Then, cytoscape software (version 3.8) was used to construct the network (15). The hub gene was the gene with highest connectivity.

Statistical analysis

We performed all the analyses using R software (version 4.0). For LASSO penalized regression analysis, we firstly calculated the cell scores using HTSeq-FPKM data from TCGA with the xCell. Only 1,063 breast cancer patients with a survival time longer than 0 day were included in the present analysis. Then the patients were randomly separated into two sets, training set and test set. We performed LASSO Cox regression with cell scores of the training-set patients. Depending on the regulation weight λ, all regression coefficients are shrunken to towards zero in LASSO, and the irrelevant features are set exactly to zero. Risk scores were calculated by as our previously study (3,16). We used “glmnet” package (R package version 4.0-2) to conduct the LASSO analysis and a P value <0.05 was considered statistically significance. Receiver operating characteristic (ROC) curve was drawn and the corresponding area under the ROC curve (AUC) was calculated to evaluate the prognostic value of the risk score by using ROCR package (R package version 1.0-11).

Regarding survival analysis, the patients were separated into two groups according to the expression level of a gene or the risk scores, and the median was used as cut-off. Then, the log-rank test was used to assess the overall survival (OS) with survival package (R package version 3.1-7). We calculated hazard ratios (HRs) and their 95% confidence intervals (CIs) using Cox proportional hazards.

We performed enrichment analyses using Gene Set Enrichment Analyses (GSEA) software (version 4.1). False discovery rate (FDR) <0.05, nominal P value <0.05 and normalized enrichment score >1 were used as criteria of statistical significance.

Results

Assessment of cell scores in tumor microenvironment

To identify signature predicting breast cancer survival in tumor microenvironment, we investigated the gene expression profile of 1,097 female breast cancer patients from TCGA. We firstly estimated the cell scores for each sample, and then compared the scores between tumors and adjacent normal tissues. Figure S1 presents the scores of all the cells. presents the cells with a median score larger than the median of all the cell scores.

Figure 1

Cell score difference between breast cancer tissues and adjacent normal tissues in the TCGA cohort. Only the cells with a median score larger than the median of all the cell score are presented. Red, breast cancer tissues; and Green, adjacent normal tissues. aDC, activated dendritic cell; cDC, conventional dendritic cell; CLP, common lymphoid progenitor; mv, microvascular; ly, lymphatic; HSC, hematopoietic stem cell; iDC, immature dendritic cell; MEP, megakaryocyte-erythroid progenitor; NKT, natural killer T cell; MSC, mesenchymal stem cell; TCGA, The Cancer Genome Atlas.

Conventional dendritic cell (cDC) predicts survival of breast cancer patients

The patients were randomly allocated to two sets, training set and test set. The characteristics of the two set patients were similar (Table S1). Key biomarkers associated breast cancer survival were screened from cell scores in training cohort using LASSO regression. We obtained four biomarkers including myocytes, natural killer T cell (NKT), cDC and sebocytes, and their coefficients were 0.098, −0.131, −0.021 and 0.012. The constructed risk score model divided patients into low- and high-risk groups according to their risk scores. Low-risk group had longer survival and less deaths (). Survival analysis indicated that high risk patients had a HR of 1.75 (95% CI: 1.08−2.83; P=0.022) (). ROC curve analysis indicated that the AUC was 0.90, 0.63 and 0.58 for 3-, 5- and 10-year survival ().

Figure 2

LASSO Cox regression of the cell scores with TCGA cohort. (A) The characteristics of the patients in training cohort order by the risk score calculated from the 4-cell signature. Dotted line: the median of the risk score. Top panel: the distribution of the risk score for the patients. Middle panel: survival time and status of the patients in the high- and low-risk groups, as defined by the risk score. Bottom panel: heatmap of the 4 cell scores. (B) The high-risk group has significantly shorter survival times than low-risk group in training cohort (HR =1.75, 95% CI: 1.08−2.83, P=0.022). (C) Receiver operating characteristic (ROC) curve of the 4-cell signature in test cohort. The area under the ROC curve (AUC) was 0.90, 0.63 and 0.58 for 3-, 5- and 10-year survival. (D) The characteristics of the patients in test cohort order by the risk score. (E) The high-risk group has significantly shorter survival times than low-risk group in test cohort (HR =1.68, 95% CI: 1.07−2.66, P=0.024). (F) ROC curve of the 4-cell signature in test cohort. The AUC was 0.71, 0.66 and 0.65 for 3-, 5- and 10-year survival. (G) Frequency of gene signatures selected by Lasso-bootstrapping highlighted the prominence of cDC contributing to predictive model. (H) Higher cDC score associated with longer overall survival in overall cohort (HR =0.61, 95% CI: 0.44−0.85; P=0.003). cDC, conventional dendritic cell; CI, confidence interval; CMP, common myeloid progenitor; HR, hazard ratio; HSC, hematopoietic stem cell; LASSO, least absolute shrinkage and selection operator; TCGA, The Cancer Genome Atlas.

The analysis of the test cohort corroborated the findings in training cohort (). The HR of high-risk patients was 1.68 (95% CI: 1.07−2.66, P=0.024; ) compared to low-risk score patients; The AUC was 0.71, 0.66 and 0.65 for 3-, 5- and 10-year survival ().

LASSO bootstrapping showed that cDC was a promising biomarker surpassing more than 4,000 counterparts covering the memory B-cells, keratinocytes and NKT (). Kaplan-Meier survival analysis indicated that higher cDC score was associated with longer OS in overall TCGA cohort (HR =0.61, 95% CI: 0.44−0.85; P=0.0033) (). We further validated it in the other two larger datasets, METABRIC (n=1,904) and GSE96058 (n=3,409). Breast cancer patients with high cDC score in METABRIC had a HR of 0.65 (95% CI: 0.58−0.73; P<0.001), and a similar result was found in GSE96058 (HR =0.70, 95% CI: 0.57−0.87; P<0.001) (Figure S2).

cDC score is associated with pathological features in breast cancer

We investigated the relationships between cDC score and pathological features of breast cancer with the data from TCGA and METABRIC (). Regarding TCGA, cDC score had a significant association with the molecular classification (P<0.001; ) and stage (P=0.004; ), but not the expression level of Ki67 (P=0.254; ). With respect to METABRIC, cDC score was significantly associated with stage (P<0.001; ) and the expression level of Ki67 (P<0.001; ) but not the molecular classification (P=0.348; ). From , we noted that HER2-positive breast cancer had lower cDC score than HER2-negative breast cancer. Therefore, we explored whether a difference in cDC score exists between HER2-positive and HER2-negative breast cancer. TCGA (P<0.001) but not METABRIC (P=0.712) showed a significant association (Figure S3A,S3B). Spearman’s correlation analysis found that cDC score was inversely correlated with the expression level of HER2 (also known as ERBB2) in both TCGA (P<0.001) and METABRIC (P<0.001) (Figure S3C,S3D). After adjusted for ERBB2 in the survival analysis, cDC score still had a significant association with breast cancer survival (TCGA: HR =0.61, 95% CI: 0.44−0.86, P=0.004; METABRIC: HR =0.67, 95% CI: 0.59−0.75, P<0.001).

Figure 3

cDC score is associated with pathological features in breast cancer. (A) cDC score in different molecular subtypes of breast cancer in TCGA. (B) cDC score in different stages of breast cancer in TCGA. (C) Distribution of cDC score according to the expression of Ki67 in TCGA. The median was used as cut-off. (D) cDC score in different molecular subtypes of breast cancer in METABRIC. (E) cDC score in different stages of breast cancer in METABRIC. (F) Distribution of cDC score according to the expression of Ki67 in METABRIC. The median was used as cut-off. cDC, conventional dendritic cell; TCGA, The Cancer Genome Atlas.

High cDC score may predicate better pathological complete response (pCR) rate

We also explored the association between cDC score and pCR rate using datasets from GEO including GSE20194, GSE22358, GSE25066 and GSE32646 (). GSE25066 indicated that high cDC score was significantly associated with pCR (P=0.008) (). GSE22358 (P=0.173) and GSE32646 (P=0.170) showed a borderline association (). However, we failed to find a significant association between cDC score and pCR using GSE20194 ().

Figure 4

High cDC score may predicate better pathological complete response (pCR) rate. (A) The pCR rate in GSE20194. (B) The pCR rate in GSE22358. (C) The pCR rate in GSE25066. (D) The pCR rate in GSE32646.

Tumors with high cDC score have elevated immune activity

In order to uncover the mechanisms of cDC score on outcome of breast cancer patients, GSEA was performed to compare the transcriptomes of high- and low-score tumors at the level of Hallmarks gene sets. We noted that all the 19 significant Hallmark gene sets enriched in high-score tumors (Table S2). The top six gene sets were IL-2 STAT5 signaling, complement, inflammatory response, KRAS signaling up, allograft rejection, and IL-6 JAK STAT3 signaling (Figure S4).

IL-2 was a key gene associated with high cDC score

To identify the key genes associated with high cDC score from Table S2, we firstly analyzed the associations between the genes and survival of breast cancer patients, and obtained 135 genes with a significant association (Table S3). Then, we constructed a network with the 135 genes using STRING. shows that IL-2 was the hub gene of the network. Breast cancer patients with high IL-2 expression had a HR of 0.56 (95% CI: 0.40−0.78; P<0.001) (Table S2 and ). Spearman’s correlation analysis found that the expression level of IL-2 was positively correlated with cDC score (P<0.001) (). We further analyzed the relationship between IL-2 and other cells. All the cells with a P value less than 0.05 are listed in Table S4.

Figure 5

IL-2 is a key gene associated with high cDC score. (A) IL-2 is the hub gene of the protein-protein interaction (PPI) network. (B) High IL-2 expression is associated with a better overall survival (HR =0.56; 95% CI: 0.40–0.78; P<0.001). (C) The expression level of IL-2 was positively correlated with cDC score. cDC, conventional dendritic cell; CI, confidence interval; HR, hazard ratio.

Discussion

In present study, we investigated the association between cell enrichment in microenvironment and survival of breast cancer. We estimated cell scores of each sample, and found that 50 cells had different scores between tumors and adjacent normal tissues, suggesting cells in tumor microenvironment may have important roles in breast cancer. Then, we used LASSO regression analysis to identify signatures predicting breast cancer survival from cell scores, and obtained a 4-cell signature. In addition, cDC score was found to be a key signature associated with longer survival time. Further analyses showed that cDC score was significantly associated with stage of breast cancer and expression level of Ki67. Since Ki67 is a nuclear marker expressed in all phases of the cell cycle other than the G0 phase, Ki67 is usually used to assess cell proliferation (17,18). Our results indicated that low cDC score may be associated with high cell proliferation. We also found that cDC score was inversely correlated with the expression level of HER2. HER2 is a member of the epidermal growth factor family of receptor tyrosine kinases (19). Presence of HER2 protein overexpression or ERBB2 gene amplification has a prognostic and predictive role in breast cancer (20). Breast cancer with HER2 overexpression is associated with malignant tumor growth and poor clinical outcome (21). Finally, we explored the association between cDC score and pCR, and high cDC score may predicate better pCR rate. Taking together, high cDC score was associated with better prognosis of breast cancer.

Dendritic cells (DCs) are major antigen-presenting cells in tumor immune microenvironment and exert anti-cancer activity (22,23). DCs are divided into two principal cell populations, including cDCs and plasmacytoid DCs (24,25). cDCs are subdivided into two subtypes, type 1 cDC (cDC1) and type 2 cDC (cDC2) (23,24). Both cDC1 and cDC2 had antitumor activities. cDC1s are involved in cross-presenting tumor-associated antigens to CD8+ T cells which are important effector cells involving tumor cell elimination (24). cDC2s at least exert antitumor effects by inducing the proliferation of antitumor CD4+ T cells (26). Ferris et al. reported that both CD4+ and CD8+ T cells are needed for tumor rejection, selective deletion of major histocompatibility complex (MHC) class II in cDC1s, and inhibition of early CD4+ T-cell priming, indicating that cDC1s and cDC2s are necessary for CD4+ T-cell priming (27). Therefore, cDCs play immunoregulatory roles in the tumor microenvironment. However, Oshi et al. did not find an association between cDC and disease-specific survival (DSS) with data from TCGA, but METABRIC showed a similar result as ours (22). Using different outcome may contribute to this difference. We used OS but Oshi et al. used DSS for TCGA. Combining the findings in our study, cDC score may be a key signature predicting OS of breast cancer.

We explored the potential mechanisms underlining the cDCs. GSEA analysis showed that tumors with high cDC score have elevated immune activity. We further found that IL-2 was a key gene associated with high cDC score. IL-2, the first of a series of lymphocytotrophic hormones to be recognized, is released from activated T lymphocytes and rapidly cleared from human circulation with a half-life of 3–22 min (28). IL-2 is pivotal for the generation and regulation of the immune response (28,29). IL-2 is predominately secreted by CD4+ and CD8+ T cells (30,31), and DCs also secret a lesser amount of IL-2 (32). As the first cytokine approved by the FDA to be used in cancer treatment, IL-2 exerts antitumor effects through immunomodulating T and NK cells in the tumor microenvironment.

When interpreting the results, the limitations of the present study need to be considered. We did not validate our results using experimental method. The mechanisms underlying cDCs are needed to further confirmed. We raised cDCs as a key signature, however, myocytes, NKT and sebocytes do indeed have their effects on breast cancer patients. Although our results suggested that cDC score was a key signature predicting prognosis for breast cancer, it is difficult to specifically profile cDCs from transcriptomic data (33).

Conclusions

In conclusion, cDC score was a key signature predicting prognosis for breast cancer and high cDC score was associated with elevated immune activity and better prognosis of breast cancer. cDCs may exert antitumor effects by upregulating IL-2.