Identification and validation of tumor microenvironment-related prognostic biomarkers in breast cancer
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. Figure 1 presents the cells with a median score larger than the median of all the cell scores.
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 (Figure 2A). Survival analysis indicated that high risk patients had a HR of 1.75 (95% CI: 1.08−2.83; P=0.022) (Figure 2B). ROC curve analysis indicated that the AUC was 0.90, 0.63 and 0.58 for 3-, 5- and 10-year survival (Figure 2C).
The analysis of the test cohort corroborated the findings in training cohort (Figure 2D). The HR of high-risk patients was 1.68 (95% CI: 1.07−2.66, P=0.024; Figure 2E) compared to low-risk score patients; The AUC was 0.71, 0.66 and 0.65 for 3-, 5- and 10-year survival (Figure 2F).
LASSO bootstrapping showed that cDC was a promising biomarker surpassing more than 4,000 counterparts covering the memory B-cells, keratinocytes and NKT (Figure 2G). 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) (Figure 2H). 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 (Figure 3). Regarding TCGA, cDC score had a significant association with the molecular classification (P<0.001; Figure 3A) and stage (P=0.004; Figure 3B), but not the expression level of Ki67 (P=0.254; Figure 3C). With respect to METABRIC, cDC score was significantly associated with stage (P<0.001; Figure 3E) and the expression level of Ki67 (P<0.001; Figure 3F) but not the molecular classification (P=0.348; Figure 3D). From Figures 3A,3D, 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).
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 (Figure 4). GSE25066 indicated that high cDC score was significantly associated with pCR (P=0.008) (Figure 4C). GSE22358 (P=0.173) and GSE32646 (P=0.170) showed a borderline association (Figure 4B,4D). However, we failed to find a significant association between cDC score and pCR using GSE20194 (Figure 4A).
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. Figure 5A 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 Figure 5B). Spearman’s correlation analysis found that the expression level of IL-2 was positively correlated with cDC score (P<0.001) (Figure 5C). 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.
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.
Acknowledgments
Funding: This work was supported by
Footnote
Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://dx.doi.org/10.21037/tcr-21-1248
Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://dx.doi.org/10.21037/tcr-21-1248). The authors have no conflicts of interest to declare.
Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013). Institutional ethical approval and informed consent were waived.
Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.
References
- Siegel RL, Miller KD, Jemal A. Cancer statistics, 2019. CA Cancer J Clin 2019;69:7-34. [Crossref] [PubMed]
- Zeng H, Zheng R, Zhang S, et al. Female breast cancer statistics of 2010 in China: estimates based on data from 145 population-based cancer registries. J Thorac Dis 2014;6:466-70. [PubMed]
- Zhong S, Chen H, Yang S, et al. Identification and validation of prognostic signature for breast cancer based on genes potentially involved in autophagy. PeerJ 2020;8:e9621. [Crossref] [PubMed]
- Mao Y, Keller ET, Garfield DH, et al. Stromal cells in tumor microenvironment and breast cancer. Cancer Metastasis Rev 2013;32:303-15. [Crossref] [PubMed]
- Paget S. The distribution of secondary growths in cancer of the breast. 1889. Cancer Metastasis Rev 1989;8:98-101. [PubMed]
- Wang J, Zhang Q, Wang D, et al. Microenvironment-induced TIMP2 loss by cancer-secreted exosomal miR-4443 promotes liver metastasis of breast cancer. J Cell Physiol 2020;235:5722-35. [Crossref] [PubMed]
- Houthuijzen JM, Jonkers J. Cancer-associated fibroblasts as key regulators of the breast cancer tumor microenvironment. Cancer Metastasis Rev 2018;37:577-97. [Crossref] [PubMed]
- Jia D, Li S, Li D, et al. Mining TCGA database for genes of prognostic value in glioblastoma microenvironment. Aging (Albany NY) 2018;10:592-605. [Crossref] [PubMed]
- Chu DT, Phuong TNT, Tien NLB, et al. The Effects of Adipocytes on the Regulation of Breast Cancer in the Tumor Microenvironment: An Update. Cells 2019;8:857. [Crossref] [PubMed]
- Aran D, Hu Z, Butte AJ. xCell: digitally portraying the tissue cellular heterogeneity landscape. Genome Biol 2017;18:220. [Crossref] [PubMed]
- Tibshirani R. Regression Shrinkage and Selection via the Lasso. J R Stat Soc Series B Stat Methodol 1996;58:267-88.
- Tomczak K, Czerwińska P, Wiznerowicz M. The Cancer Genome Atlas (TCGA): an immeasurable source of knowledge. Contemp Oncol (Pozn) 2015;19:A68-77. [Crossref] [PubMed]
- Zhong S, Ma T, Zhang X, et al. MicroRNA expression profiling and bioinformatics analysis of dysregulated microRNAs in vinorelbine-resistant breast cancer cells. Gene 2015;556:113-8. [Crossref] [PubMed]
- Szklarczyk D, Gable AL, Lyon D, et al. STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res 2019;47:D607-13. [Crossref] [PubMed]
- Otasek D, Morris JH, Bouças J, et al. Cytoscape Automation: empowering workflow-based network analysis. Genome Biol 2019;20:185. [Crossref] [PubMed]
- Liu M, Zhou S, Wang J, et al. Identification of genes associated with survival of breast cancer patients. Breast Cancer 2019;26:317-25. [Crossref] [PubMed]
- Penault-Llorca F, Radosevic-Robin N. Ki67 assessment in breast cancer: an update. Pathology 2017;49:166-71. [Crossref] [PubMed]
- Nielsen TO, Leung SCY, Rimm DL, et al. Assessment of Ki67 in Breast Cancer: Updated Recommendations From the International Ki67 in Breast Cancer Working Group. J Natl Cancer Inst 2021;113:808-19. [Crossref] [PubMed]
- Redmond AM, Omarjee S, Chernukhin I, et al. Analysis of HER2 genomic binding in breast cancer cells identifies a global role in direct gene regulation. PLoS One 2019;14:e0225180. [Crossref] [PubMed]
- Trapani D, Rajasekar AKA, Mathew A. More options for adjuvant treatment of HER2-positive breast cancer: How to choose wisely? Int J Cancer 2019;145:2901-6. [Crossref] [PubMed]
- Wei T, Xing H, Wang H, et al. Bovine serum albumin encapsulation of near infrared fluorescent nano-probe with low nonspecificity and cytotoxicity for imaging of HER2-positive breast cancer cells. Talanta 2020;210:120625. [Crossref] [PubMed]
- Oshi M, Newman S, Tokumaru Y, et al. Plasmacytoid Dendritic Cell (pDC) Infiltration Correlate with Tumor Infiltrating Lymphocytes, Cancer Immunity, and Better Survival in Triple Negative Breast Cancer (TNBC) More Strongly than Conventional Dendritic Cell (cDC). Cancers (Basel) 2020;12:3342. [Crossref] [PubMed]
- Michea P, Noël F, Zakine E, et al. Adjustment of dendritic cells to the breast-cancer microenvironment is subset specific. Nat Immunol 2018;19:885-97. [Crossref] [PubMed]
- Kim CW, Kim KD, Lee HK. The role of dendritic cells in tumor microenvironments and their uses as therapeutic targets. BMB Rep 2021;54:31-43. [Crossref] [PubMed]
- Fu C, Peng P, Loschko J, et al. Plasmacytoid dendritic cells cross-prime naive CD8 T cells by transferring antigen to conventional dendritic cells through exosomes. Proc Natl Acad Sci U S A 2020;117:23730-41. [Crossref] [PubMed]
- Binnewies M, Mujal AM, Pollack JL, et al. Unleashing Type-2 Dendritic Cells to Drive Protective Antitumor CD4+ T Cell Immunity. Cell 2019;177:556-571.e16. [Crossref] [PubMed]
- Ferris ST, Durai V, Wu R, et al. cDC1 prime and are licensed by CD4+ T cells to induce anti-tumour immunity. Nature 2020;584:624-9. [Crossref] [PubMed]
- Waldmann TA. The biology of interleukin-2 and interleukin-15: implications for cancer therapy and vaccine design. Nat Rev Immunol 2006;6:595-601. [Crossref] [PubMed]
- Yang Y, Lundqvist A. Immunomodulatory Effects of IL-2 and IL-15; Implications for Cancer Immunotherapy. Cancers (Basel) 2020;12:3586. [Crossref] [PubMed]
- Liao W, Lin JX, Leonard WJ. Interleukin-2 at the crossroads of effector responses, tolerance, and immunotherapy. Immunity 2013;38:13-25. [Crossref] [PubMed]
- Paliard X, de Waal Malefijt R, Yssel H, et al. Simultaneous production of IL-2, IL-4, and IFN-gamma by activated human CD4+ and CD8+ T cell clones. J Immunol 1988;141:849-55. [PubMed]
- Granucci F, Vizzardelli C, Pavelka N, et al. Inducible IL-2 production by dendritic cells revealed by global gene expression analysis. Nat Immunol 2001;2:882-8. [Crossref] [PubMed]
- Sturm G, Finotello F, Petitprez F, et al. Comprehensive evaluation of transcriptome-based cell-type quantification methods for immuno-oncology. Bioinformatics 2019;35:i436-45. [Crossref] [PubMed]