A novel model associated with tumor microenvironment on predicting prognosis and immunotherapy in Triple Negative Breast Cancer

doi:10.21203/rs.3.rs-2730381/v1

Download PDF

Research Article

A novel model associated with tumor microenvironment on predicting prognosis and immunotherapy in Triple Negative Breast Cancer

https://doi.org/10.21203/rs.3.rs-2730381/v1

This work is licensed under a CC BY 4.0 License

Journal Publication

published 23 May, 2023

Read the published version in Clinical and Experimental Medicine →

You are reading this latest preprint version

Background

TNBC (Triple Negative Breast Cancer) is the most aggressive and malignant subtype in breast cancer. Immunotherapy is a currently promising and effective treatment for TNBC, while not all patients are responsive. Therefore, it is necessary to explore novel biomarkers to screen sensitive populations for immunotherapy.

Methods

All mRNA expression profiles of TNBC from The Cancer Genome Atlas (TCGA) database were clustered into two subgroups by analyzing tumor immune microenvironment (TIME) with single sample gene set enrichment analysis(ssGSEA). A risk score model was constructed based on differently expressed genes (DEGs) identified from two subgroups using Cox and Least Absolute Shrinkage and Selector Operation (LASSO) regression model. And it was validated by Kaplan–Meier analysis and Receiver Operating Characteristic (ROC) analysis in Gene Expression Omnibus (GEO) and the Molecular Taxonomy of Breast Cancer International Consortium (METABRIC) databases. Multiplex immunofluorescence (mIF) and Immunohistochemical (IHC) staining were performed on clinical TNBC tissue samples. The relationship between risk score and immune checkpoint blockades (ICB) related signatures was further investigated, as well as the biological processes were performed by gene set enrichment analysis (GSEA).

Results

We obtained three DEGs positively related to prognosis and infiltrating immune cells in TNBC. Our risk score model could be an independent prognostic factor and the low risk group exhibited a prolonged overall survival (OS). Patients in low risk group were more likely to present a higher immune infiltration and stronger response to immunotherapy. GSEA revealed the model was associated with immune-related pathways.

Conclusion

We constructed and validated a novel model based on three prognostic genes related to TIME in TNBC. The model contributed a robust signature that could predict the prognosis and response to immunotherapy in TNBC.

Triple-negative breast cancer (TNBC)

TCGA

Tumor immune microenvironment (TIME)

Prognosis

Immunotherapy

Breast cancer is a tumor with a high prevalence rate in women^[1], and its incidence and mortality are showing an upward trend year by year. The latest global cancer burden data in 2020 shows that breast cancer surpasses lung cancer to become the largest cancer in the world^[2]. TNBC is one of the subtypes of breast cancer, which occupies about 15–20%. The prognosis of TNBC is extremely worst with the character of highest histological grade and rates of metastasis^{[3, 4]}. Due to the lack of expression of estrogen receptor, progesterone receptor and human epidermal growth factor receptor 2 (HER-2), TNBC is not sensitive to endocrine therapy and HER-2 targeted therapy.

With the emergence of large-scale sequencing and other multi-omics technologies such as gene expression profiling, proteomics, immune omics and other detection methods, it has been proved that TNBC has stronger heterogeneity^[5]. TNBC is no longer a kind of disease, but a collection of diseases^[6–10]. In response to the heterogeneous feature of TNBC and to achieve precise treatment of tumors, different subtypes should be given different treatments according to the biological characteristics of each subtype. At present, immunotherapy is a hot spot in cancer which treats tumors by regulating the immune system response. Studies showed that immunotherapy may have a better effect on highly immunogenic tumors (inflammatory tumors or "hot" tumors)^[11]. TNBC obtains unique biomolecular characteristics and higher immunogenicity than other types of breast cancer^{[12, 13]}. A large number of clinical trials confirmed that TNBC patients who received immunotherapy and responded well could have a better prognosis and long-lasting benefits. The classic results of KEYNOTE-522 study^[14] and Impassion130 study^[15] suggested that TNBC could benefit from immunotherapy in both early and advanced stages. Immunotherapy is emerging as a promising and effective treatment option for TNBC. However, not all TNBC patients were sensitive to immunotherapy. In order to maximize the benefits of treatment and reduce the side effects of drugs, it is necessary to screen out the patients with a sensitive response to immunotherapy^[16].

The tumor microenvironment (TME) plays an important role in tumorigenesis and tumor progression of many cancers, which is related to the prognosis of patients^[17]. Cancers develop through continuous interaction between tumor cells and various immune cells in TME. Stimulating a sustained antitumor immune response is the key factor to achieve the clinical benefits of immunotherapy^[18]. Tumor-infiltrating lymphocytes (TILs) as an important part of TME were considered closely related to efficacy of the immunotherapy and favorable prognosis in breast cancer^{[19, 20]}. But the accuracy and manifestation of TILs were found to be limited^[21].

In this study, we systemically portrayed the TIME in TNBC by analyzing the mRNA sequencing profiles from TCGA database. After clustering all samples into two subgroups, we distinguished three DEGs related to immune infiltration and prognosis of patients. A risk score model was constructed based on the DEGs to predict the survival of patients and the efficiency of immunotherapy in TNBC. We further validated the risk score model by GEO and MATABRIC cohorts. Furthermore, we determined the risk score to be an independent prognostic signature that could accurately predict the 3- and 5-year OS of patients with TNBC by ROC analysis. We also investigated the correlation between the model and proportions of infiltrating immune cells in tumor tissue slides from our hospital. Our finding clarified the relationship between risk score and ICB related signatures, which indicated the predictive significance of the model. The potential mechanism of risk score was also explored, which was related to immune-related signal ways. We believed that this model would improve risk stratification in TNBC and serve as a potential immunotherapy biomarker for patients to provide more exact judgment for individualized clinical management.

Study Design

The flow chart of the research was shown in Fig. 1. We analyzed the mRNA profiles of TNBC samples from TCGA databases and clustered them into two different subgroups according to the immune characteristics of TME. Then we performed difference analysis of prognosis and immune between the two groups. By Cox and LASSO analysis, the DEGs were distinguished. Then we evaluated the immune-related and prognostic significance of DEGs from TIMER and K-M plotter website. At the same time, we constructed a risk score model and verified it in GEO and METABRIC database. We also validated the relationship between the model and infiltrating immune cells of TNBC clinical tissue samples using mIF and IHC staining. We acclaimed the predictive value of the model for immunotherapy and prognosis in TNBC. The correlations between the risk score and ICB related signatures were further explored in the TCGA datasets. GSEA was conducted for biological functional exploration of risk score model.

Database Collection

TNBC datasets were downloaded from TCGA database (https://gdc.cancer.gov/), including clinical data and mRNA-seq profiles measured by Illumina HiSeq2000 RNA Sequencing platform. The GSE58812 datasets and the METABRIC datasets were used for validation. The clinical characterizes and mRNA-seq data were downloaded from the GEO portal (https://www.ncbi.nlm.nih.gov/geo/) and cBioPortal (https://www.cbioportal.org/datasets), respectively. In the end, 112 patients with TNBC were enrolled this research from TCGA database, as well as 107 patients were obtained from GSE58812 datasets and 209 patients were from METABRIC datasets.

Landscape Of Immune Infiltration In Tumor Microenvironment

GSVA package was conducted to realize ssGSEA^[22] and calculate the enrichment scores of 58524 genes representing 29 immune associated gene sets of TCGA cohorts. The TNBC samples were stratified into two subgroups, which were high and low immunogenicity groups, based on the immune characteristic clustering using R package, sparcl^[23]. We defined them as Immunity-H and Immunity-L group, respectively. Estimation of Stromal and Immune Cells in Malignant Tumors using the Expression Data (ESTIMATE) algorithm was used to calculate the tumor purity, ESTIMATE score, immune score and stromal score of TNBC samples by using R package, estimate. Levels of infiltrating immune cells and stromal cells were predicted, which formed the basis for inferring tumor purity^[24]. The CIBERSORT algorithm was used to identify the relative fractions of 22 tumor infiltrating immune cells in each tumor samples, using R package, CIBERSORT. The bar charts of each sample were drawn using R package, barplot^[25]. We used multivariate Cox regression analysis to explore the factors associated with TNBC with immunity subgroup and migration status as covariates. Survival comparisons between Immunity-H and Immunity-L group were plotted with R-packages, survival and survminer.

Degs Between Subgroups And Construction Of Risk Score Model

Based on the previous clustering subgroups, DEGs were determined by using R package, limma. The threshold of DEGs were |log2 fold change (log2FC)| > 2 and false discovery rate (FDR) < 0.05. All results were visualized by heatmap with the R-package, pheatmap. Volcano plot was draw by R package, ggplot. To quantify the effects of immune and prognostic status and minimize overfitting, Cox and LASSO regression model using glmnet package in R software were adopted to screen out the optimal DEGs signatures and construct a risk score model^[26]. The risk score was calculated according to the expression level of each gene and its corresponding regression coefficient. The following formula was performed ^[27]. Using cox regression model with R package, survival to confirm the prognostic significance of the risk score model, compared with clinical information of TNBC patients.

The Immune-related And Prognostic Significance Of Degs

TIMER (https://cistrome.shinyapps.io/timer/) was used to confirm the correlation between DEGs and six tumor infiltrating immune cells. Kaplan-Meier Plotter (http://kmplot.com/analysis/) was adopted to analyze the prognosis of DEGs based on Affymetrix microarrays. Hazard ratios (HR) with 95% confidence intervals and log rank P-value were determined.

Validation Of Risk Score Model

TCGA datasets were set as training group, while GSE58812 and METABRIC datasets were set as validation group to prove the efficiency of risk score model using R package, ggrisk. According to the median risk score, we divided patients into two groups. Then, survival comparisons between high and low risk score groups were visualized by Kaplan–Meier survival curves and examined by log rank test using the R-packages, survival and survminer. ROC curves were developed to estimate the sensitivity and predictive ability of risk score model by R package, timeROC. After the approval of the ethics committee and with the informed consent of the patients, 32 Formalin-fixed paraffin-embedded (FFPE) tissue sections diagnosed as TNBC were collected from the First Affiliated Hospital of Xi'an Jiaotong University during 2015 to 2020 to analyse the correlation between risk score and immune cells in TME.

Detection Of The Fractions Of Immune Cells And Expression Of Degs

Various immune cells were tested by multiplex immunofluorescence (mIF) in PerkinElmer Mantra Quantitative Pathology Workstation/Quantitative Pathology Analysis Platform. The detected marker and its corresponding cell type was presented in Table 1. Immunohistochemistry (IHC) was performed on FFPE sections using CXCL13 antibody(rabbit, 1: 200, proteintech), CCL5 antibody༈rabbit, 1༚200, Affinity༉and GZMB antibody༈rabbit, 1༚200, proteintech༉ according to the manufacturer’s instructions to detect DEGs. The staining was confirmed by two independent investigators. The following method was used to quantify the expression of genes. Firstly, three fields of view were randomly selected and 20 cells were counted in each field of view. The staining results were judged according to the staining intensity and the percentage of positive cells. Staining intensity (SI): negative (no staining) 0 score, weak positive (light yellow) 1 score, medium positive (brown) 2 score, strong positive (brown) 3 score; and the positive cell percentage (PP): no positive cells 0 score, 1–50% 1 score, 51–100% 2 score. Secondly, multiplying the SI and PP is the final staining score of the field, and the sum of the three fields was taken as the average value.

Table 1

Detected marker and its corresponding cell type
Detected Marker	PANCK+	CD8+	CD56+	CD68+
Detected Marker	PANCK+	CD8+	CD56+	HLA-DR+	HLA-DR-
Cell type	Tumor cells	CD8 + T cells	Natural Killer (NK) cells	Macrophages M1	Macrophages M2

Correlation Between Immune Checkpoint Genes And Ips

Six key genes of immune checkpoint blockade therapy were enrolled to analyze the correlation with risk score in TNBC, which were CD274 (PD-L1, programmed death ligand 1), programmed death 1 (PD-1, also known as PDCD1), cytotoxic T-lymphocyte antigen 4 (CTLA-4), T-cell immunoglobulin domain and mucin domain-containing molecule-3 (TIM-3, also known as HAVCR2), indoleamine 2,3-dioxygenase 1 (IDO1) and lymphocyte activating gene 3 (LAG3). The Immunophenoscore (IPS) were calculated on a 0–10 scale based on four main parts (effector cells, immunosuppressive cells, MHC (major histocompatibility complex) molecules, and immunomodulators) related gene expression z-scores using machine learning to determine the immunogenicity. The higher IPS was associated with increased immunogenicity and it has been verified that IPS could predict the response of patients to immunotherapy. The IPS of breast cancer patients were downloaded from The Cancer Immunome Atlas (TCIA) (https://tcia.at/home)^[28].

Gene Set Enrichment Analysis

The GSEA (version 4.1.0) software (https://www.gsea-msigdb.org/gsea/index.jsp) was performed to explore the potential molecular mechanisms between the RiskScore_H and RiskScore_L groups in TNBC. A NOM p value < 0.05 and FDR q value < 0.25 were considered statistically significant.

Statistical analysis

Statistical analyses were performed with R software (version 4.1.2), IBM SPSS Statistic 25.0 and GraphPad Prism 8.0. Differences between two groups were compared using two independent sample t-tests. Pearson analysis was performed to evaluate the correlation between risk score with ICB related genes. The correlation of risk score with OS was performed by Kaplan–Meier analysis and log rank test.

Landscape of immune cell infiltration of TME in TNBC

Using ssGSEA algorithm, we calculated 29 immune-related gene sets of TNBC datasets from TCGA database and described the TIME status of TNBC patients. By unsupervised clustering, the samples were divided into two different immune-related subgroups, including Immunity_H group (52 samples) with a high expression level of immune-related gene sets and Immunity_L group (60 samples) with a low expression level of immune-related gene sets. We also calculated the immune scores (ESTIMATEScore, ImmuneScore, StromalScore) and tumor purity between two subgroups using the ESTIMATE algorithm. The involvement of the subgroups with the immune scores and tumor purity was explored and depicted in the comprehensive heatmap (Fig. 2a). It suggested that the immune scores of Immunity_H group were generally higher than those of Immunity_L group, and the tumor purity was lower than that of Immunity_L group (P < 0.05; Fig. 2b).

The Immune Status And Prognostic Difference Between Two Subgroups

In order to further clarify and verify the accuracy of the immune differences between Immunity_H and Immunity_L groups, firstly, we compared some immune-related gene expression levels between Immunity_H and Immunity_L groups. The expression levels of MHC gene family, CD4, CD8, CD274 (PD-L1) and CTLA4 in Immunity_H group were higher than Immunity_L group (P < 0.05; Fig. 3a-b). Meanwhile, the relative subpopulations of infiltrating immune cells were estimated with CIBERSORT approach and the distribution of the proportion of immune cells in TNBC were characterized by bar plot. By comparison, the proportion of immune cells differed between Immunity_H and Immunity_L groups. Patients from Immunity_L group exhibited a remarkably higher infiltration of Macrophages M0, Macrophages M2 and activated Mast cells. Immunity_H group was marked by a significant increase in the subpopulation of the antitumor lymphocyte cell subsets, such as CD8 T cells, CD4 memory activated T cells, gamma delta T cells and Macrophages M1 (P < 0.05; Fig. 3c). Based on above results, it could be regarded that Immunity_H group was characterized an elevated inflammation response than that of Immunity_ L group. Interestingly, Immunity_H group was characterized by a higher immune score, which indicated an immunologically “hot” phenotype. Previous reports indicated that the cell composition of TME had an important impact on the prognosis of patients, including the degree of immune cells and stromal cells infiltration. In addition, the median overall survival time of Immunity_H group was longer than Immunity_L group (P = 0.035, HR = 0.34[0.13–0.93]); Fig. 3d).

Screening Of Degs And Construction Of Risk Score Model

In order to further explore the reasons of the difference in immunogenicity between Immunity_H and Immunity_L group, we analyzed the mRNA-seq data of TNBC patients from TCGA database. And the volcano plot and heat map were drawn to distinguish the DEGs between the two groups (Fig. 3e-f). Combined with LASSO regression (Fig. 3g) to reduce the risk of overfitting, three genes with the most significant different expression level were obtained, namely CXCL13, CCL5 and GZMB. And the risk score model was constructed with the following formula: risk score= -0.0060×CXCL13 expression-0.0046×CCL5 expression-0.0606×GZMB expression. Then we calculated the risk score of each sample and combined with its clinical characteristics to analyze the relationship between these various factors and OS by Cox risk regression analysis. Based on the results of the univariate Cox analysis (Fig. 3h), it showed that tumor stage (P < 0.001, HR = 9.615[3.368–27.450]), tumor size (P = 0.010, HR = 4.147[1.397–12.313]), lymph node status (P < 0.001, HR = 2.873[1.821–4.533]) and risk score (P = 0.017, HR = 4.674[1.312–16.645]) had a statistically significant impact on patient survival. The results of multivariate Cox analysis (Fig. 3i) revealed that lymph node status (P < 0.001, HR = 3.457[1.701–7.026]) and risk score (P = 0.037, HR = 5.347[1.107–25.829]) were favorable prognostic factors for survival of TNBC patients. The results suggested that risk score could be an independent prognostic factor in both univariable and multivariable Cox regression analyses.

The Correlation Of Degs With Immune Infiltration And Prognosis

We used the Kaplan-Meier Plotter to evaluate prognostic value of DEGs in TNBC. The results showed that the expression of DEGs was related to OS of TNBC patients, and their expression level was positively correlated with better OS among all TNBC patients (P < 0.05; CXCL13, HR = 0.31[0.2–0.48]; CCL5, HR = 0.36[0.24–0.54]; GZMB, HR = 0.31[0.2–0.48]; Fig. 4a). The TIMER database was performed to investigate whether DEGs were related to immune infiltration in breast cancer. The results revealed that the expression level of DEGs was positively correlated with tumor immune infiltration in breast cancer. Except for macrophages, other immune cells had significant statistical differences (P < 0.05; Fig. 4b). These results suggested that DEGs may serve as an immune-related tumor marker in breast cancer which could predict the prognosis in TNBC.

Validation Of Risk Score Model Based On Degs In Tnbc

In order to further verify the accuracy of the prognostic prediction of risk score model in TNBC, we confirmed the reliability in three TNBC cohorts, including 112 cases from TCGA database as training set (Fig. 5a-c), 107 cases from GEO database (Fig. 5d-f) and 209 cases from METABRIC database as validation sets (Fig. 5g-i). We calculated the risk score of each sample using the risk score model formula (risk score=-0.0060×CXCL13 expression-0.0046×CCL5 expression- 0.0606×GZMB expression). According to the median value of risk score, the TNBC patients were divided into RiskScore_L group with low risk score and RiskScore_H group with high risk score. Kaplan-Meier survival curve was performed. The results showed that the prognosis of RiskScore_L group was significantly better than that of RiskScore_H group in all three databases (P < 0.05; TCGA, HR = 0.24[0.09–0.63]; GEO, HR = 0.28[0.13–0.59]; METABRIC, HR = 0.67[0.47–0.97]; Fig. 5a, d, g). The ranking was based on the risk score values of the DEGs from low to high, the risk score distribution (Fig. 5b), patient survival status (Fig. 5e) and risk gene expression (Fig. 5h) in sets were shown, respectively. The time-dependent ROC curve analysis of the risk score in three databases revealed the robustness of the prognostic capability of the risk score for OS (TCGA, 3 year survival, AUC = 0.703; 5 year survival, AUC = 0.785; GEO, 3 year survival, AUC = 0.757; 5 year survival, AUC = 0.720; METABRIC, 3 year survival, AUC = 0.608; 5 year survival, AUC = 0.64; Fig. 5c, f,i). The above results indicated that the risk score model based on three immune-related DEGs was validated to be a strong predictive signature on the survival of TNBC patients.

We collected the clinical samples from our hospital to identify our study. We detected the immune cells in tumor microenvironment of 32 TNBC pathological tissue specimens using mIF, which including marker PanCK, CD8, HLA-DR, CD68, CD56 (Fig. 6a). The expression levels of three DEGs (CXCL13, CCL5 and GZMB) were performed by IHC staining on TNBC tissue samples (Fig. 6b). The risk score and immune cells fractions of 32 TNBC tumor tissue specimens were calculated (Supplementary Table 1). And the samples were divided into two groups based on the median of risk score. We analyzed the correlation between risk score and the proportions of immune cells in TME. The results displayed that the fractions of immune cells were different in two groups, especially about NK cells and macrophages. Low risk score group was associated with higher proportions of NK cells and macrophages M1/M2 ratio (Supplementary Fig. 1).

Predictive Significance For Immunotherapy And Immune-related Pathways

We checked the correlation between the risk score and six key ICB related genes, including CD274 (PD-L1), PDCD1, CTLA4, HAVCR2, IDO1 and LAG3. We observed that risk score had a remarkably negative relationship to the expression of these ICB related genes (P < 0.05; Fig. 7a), which showed that risk score might act as a significant role in the prediction of responsiveness to ICB treatment in TNBC. Furthermore, we confirmed the IPS to predict efficacy of ICB with PD-1 and/or CTLA-4 inhibitors, including the IPS, IPS-PD1 blocker, IPS-CTLA4 blocker and IPS- PD1 + CTLA4 blocker. In our study, the predictive value of risk score for ICB treatment were checked. RiskScore_H group was marked by a significantly lower IPS than RiskScore_L group (P < 0.05; Fig. 7b) which suggested that RiskScore_L group might be more suitable for immunotherapy. These results suggested the risk score model might be a potential predictive signature for response to immunotherapy in TNBC. At the same time, we performed GSEA to explore the biological processes related to risk score in TNBC. The top twenty significant pathways with comprehensive details were recorded in Supplementary Table 2. KEGG significant pathways associated with low risk score were mainly involved in immune-related pathways, including antigen processing and presentation, B cell receptor signaling pathway, cytokine-cytokine receptor interaction, chemokine signaling pathway, intestinal immune network for IgA production, leukocyte transendothelial migration, natural killer cell mediated cytotoxicity, T cell receptor signaling pathway (Fig. 7c). GO biological processes significant pathways associated with low risk score were also mainly enriched in immune-related signal pathways, such as IL-1 signaling pathway, IL-6 signaling pathway, leukocyte proliferation, lymphocyte activation involved in immune response, mononuclear cell differentiation, regulation of immune response, T cell activation, T cell differentiation and so on (Fig. 7d). Through GSEA results of risk score, the immune-related pathways were observed obvious enrichment in low risk score group which suggested that risk score was closely linked to immune cell infiltration, inflammatory reactions and TME modification of TNBC.

Immunotherapy has started being considered one of the emerging and effective treatments in breast cancer, especially for the patients with TNBC^[29]. TNBC has high expression levels of TILs, PD-L1, tumor mutational burden, microsatellite instability, and mismatch repair genes and other characteristics. These indicators cause the immunogenicity of TNBC which are different from other types of breast cancer. And these might be the reasons why TNBC patients are more sensitive and responsive to immunotherapy. We might can use these biological indicators as efficient predictors of beneficial population in breast cancer immunotherapy^{[5, 30]}. The interaction between tumor cells and immune cells is related to cancer development and progression^[31]. In recent years, TILs as an important part of TME have been confirmed as an prominent biomarkers in the prediction of immunotherapy ^[29]. Therefore, it is necessary to explore more signatures from TIME in TNBC. Recently, the studies of Jérôme GALON group discussed the meaning of the individual immune parameters and combination of multiple immune factors, such as immune cell populations, immune genes, cytokines and so on. They introduced a novel immune-based assay named the ‘Immunoscore’ based on the main immune parameter in cancer with classification, prognostic and predictive value^[32].

In this study, we depicted the immune infiltration status of TME with heatmap based on 29 immune-related gene sets of TNBC datasets from TCGA database which were evaluated by ssGSEA algorithm. It is suggested that TNBC had internal immunological heterogeneity which has already been reported ^[33]. Based on the immune characteristic clustering, all TNBC samples were clustered into two subgroups which were Immunity_L and Immunity_H group. Through ESTIMATE algorithm, we calculated tumor purity, stromal, and immune scores of each sample according to immune or stromal components in TME. The patients with TNBC in Immunity_H group were characterized with the abundant of immune scores which presented higher immune infiltration levels and associated with stronger immune activity^[34]. Immune cell compositions were analyzed with CIBERSORT algorithm^[25]. The infiltration of Macrophages M0, Macrophages M2 and activated Mast cells were enriched in Immunity_L group. On contrary, Immunity_H group obtained more infiltration of CD8 T cells, CD4 memory activated T cells, gamma delta T cells and Macrophages M1. Above all, it indicated that Immunity_H group is an immunologically “hot” phenotype with an elevated inflammation response. We also compared the expression levels of some confirmed immune-related genes (MHC gene family, CD4, CD8, CD274 and CTLA4)^[35] between two subgroups. The results showed the significant up-regulated expression of immune-related genes in Immunity_H group than Immunity_L group. As reported before, TME was related to tumorigenesis, which affected the outcome of patients with TNBC ^[36]. We explored the prognostic value between Immunity-H and Immunity-L group, and we found patients in Immunity_H group were accompanied with a prolonged median OS. Recent studies also indicated that “immune-inflamed” cluster in TNBC was abundant with adaptive and innate immune cells infiltration and high expression of immune checkpoint molecules. And it had significantly better progression free survival (RFS) and OS than other clusters^[11] which was similar to our research. In order to discover the cause of the difference between two subgroups, we compared the mRNA expression level of TNBC and the DEGs were distinguished using limma R package. After the application of Cox and LASSO regression analysis, three optimal DEGs were determined which were CXCL13, CCL5 and GZMB.

CXCL13 was a chemokine ligand playing an important role in recruiting lymphocytes (such as naive and effector CD4 + T cells, memory CD4 + T cells, CD8 + T cells, B cells and natural killer cells, etc.) to the site of inflammation^[37]. In breast cancer tissues, CXCL13 gene was overexpressed at mRNA level and protein level. And CXCL13 immunohistochemistry marker was related to the degree of tumor immune infiltration^[38]. Studies have proved that the expression of CXCL13 gene can predict the increase of TILs and the improvement of disease free survival and OS in early TNBC^{[39, 40]}. CCL5 was a chemokine involved in immune regulation and inflammation^[41], and its mRNA expression was significantly related to immune activation and pathologic complete response increase of neoadjuvant chemotherapy in breast cancer (P < 0.001, OR = 1.41[1.23–1.62])^[42]. In a retrospective cohort of 72 TNBC patients and a publicly available dataset analysis, there was a direct correlation between TILs and CCL5 expression in TNBC. In addition, CCL5 gene expression was significantly correlated with better RFS (P = 0.0012; HR = 0.39[0.22–0.71])^[43]. The GZMB gene encoded a member of protein granzyme subfamily and was a part of S1 family of serine protease peptidases. It was secreted by natural killer cell cells and cytotoxic T lymphocytes which was related to the apoptosis of target cells, chronic inflammation and wound healing. B cells could kill tumor cells through GZMB^{[20, 44]}. GZMB + TILs were markers of immune escape from primary tumors and were associated with metastatic tumors^[45]. In summary, it could be seen that the three genes CXCL13, CCL5 and GZMB were all immune-related prognostic genes in breast cancer. Most of the conclusions tended to be pro-inflammatory. The higher the gene expression level, the better the patient’s prognosis and response to treatment^[46]. Our study further proved the expression level of three DEGs were higher in breast tumor tissue than normal breast tissue. Meanwhile, we also confirmed the favorable effect of CXCL13, CCL5 and GZMB on prognosis and the positive correlation with immune activity in TNBC.

Based on above results, a risk score model was constructed. The formula was as follows: risk score= -0.0060×CXCL13 expression-0.0046×CCL5 expression-0.0606×GZMB expression. We confirmed the prognostic significance of the risk score model in univariable and multivariable Cox regression model. It suggested that the risk score could be regard as an independent favorable prognostic factor in TNBC. The predictive capability of the risk score model was also successfully validated in TCGA, GEO and METABRIC databases. The risk score distribution, patient survival status, risk gene expression of three DEGs and ROC curve analysis using risk score were also performed showed a better prognostic value. In all cohorts, patients in low risk score group showed prolonged survival time than high risk score group which verified the model could provide a valuable prognostic signature for patients with TNBC. We suggested the risk score was related to the fractions of different immune cells in TME. Low risk score group had higher proportions of NK cells and macrophages M1/M2 ratio which might be related to a better prognosis. Because of the small clinical enrolled sample size, we did not get the statistical difference. In subsequent research, we planned to enroll external patients with TNBC for further exploration at a larger cohort.

Our risk score model derived from TIME and it was related to the survival outcomes in TNBC, so we further explored whether risk score correlated with immunotherapy of TNBC. Refer to existing studies, the expression level of ICB related genes might be related to the efficiency of ICB treatment^[47]. Therefore, we analyzed the relationship between risk score and several important ICB related genes, such as CD274, PDCD1, CTLA4, HAVCR2, IDO1 and LAG3. The result showed that risk score was negatively correlated with the expression of ICB related genes. It suggested that the risk score may could predict the response to ICB treatment. Low risk group patients were more likely to respond to immunotherapy. At the same time, we downloaded the IPS which could predict the response of the cancer patients to immunotherapy with anti-PD-1, anti-PD-L1 and/or anti-CTLA-4 treatment from TCIA^[48]. We found that IPS of high risk score group was significantly lower than that of low risk score group. In addition, GESA enrichment analysis exhibited that a low risk score showed a mainly significant enrichment in immune-related pathways from KEGG and GO biological processes. Above all, the risk score model might be used as a predictor of immunotherapy in TNBC. We could say that the patients with low risk score showed sensitive responses to ICB treatment.

In conclusion, TNBC has internal heterogeneity. We distinguished three immune-related prognostic genes from clustering subgroups correlated with TIME. And we built and verified a novel prognostic model based on the three genes. The model could be used in molecular subgrouping and accurately predict the prognosis and the efficacy of immunotherapy in TNBC.

CTLA-4	cytotoxic T-lymphocyte antigen 4
DEGs	differentially expressed genes
ER	estrogen receptor
ESTIMATE	Estimation of Stromal and Immune Cells in Malignant Tumors using the Expression Data
FDR	false discovery rate
FFPE	Formalin-fixed paraffin-embedded
GEO	Gene Expression Omnibus
GSEA	gene set enrichment analysis
HER-2	human epidermal grow the factor receptor 2
HR	hazard ratios
ICB	immune checkpoint blockades
IDO1	indoleamine 2,3-dioxygenase 1
IHC	Immunohistochemistry
IPS	Immunophenoscore
LAG3	lymphocyte activating gene 3
LASSO	Least Absolute Shrinkage and Selector Operation
METABRIC	Molecular Taxonomy of Breast Cancer International Consortium
MHC	major histocompatibility complex
mIF	multiplex immunofluorescenc
OR	odd ratio
OS	overall survival
PD-1	programmed cell death-1
PD-L1	programmed cell death-ligand 1
PP	positive cell percentage
RFS	recurrence free survival
ROC	Receiver Operating Characteristic
SI	staining intensity
TCGA	The Cancer Genome Atlas
TCIA	The Cancer Immunome Atlas
TILs	tumor-infiltrating lymphocytes
TIM-3	T-cell immunoglobulin domain and mucin domain-containing molecule-3
TIME	tumor immune microenvironment
TME	tumor microenvironment
TNBC	triple negative breast cancer

Ethics approval and consent to participate

The studies involving human participants were reviewed and approved by the Ethics Committee of the First Affiliated Hospital of Xi’an Jiaotong University. The patients/participants provided their written informed consent to participate in this study.

Consent for publication

Not applicable.

Availability of data and materials

The datasets generated and/or analysed during the current study are available in the TCGA repository, [https://gdc.cancer.gov/], GEO portal, [https://www.ncbi.nlm.nih.gov/geo/] and cBioPortal, [https://www.cbioportal.org/datasets]. In addition, The original contributions presented in the study are included in the article/supplementary material. Further inquiries can be directed to the corresponding authors.

Competing interests

The authors declare that they have no competing interests

Funding

This study was supported by International Cooperation Foundation Project of Shaanxi Province (No. 2019KW-033) and National Natural Science Foundation of China (No. 82002803).

Authors' contributions

All authors contributed to the study conception and design. JZ and QT contributed to the study design. MZ and QT contributed to data collection. JZ and MZ performed statistical analysis and interpretation. JZ and JY drafted the manuscript. All authors read and approved the final manuscript.

Acknowledgments

We would like to thank all the databases and software for providing and analyzing transcriptomics and clinicopathological data in our research.

FERLAY J, FORMAN D, MATHERS C D, et al. Breast and cervical cancer in 187 countries between 1980 and 2010 [J]. The Lancet, 2012, 379(9824): 1390–1.
World Health Organization. (2020). [J]. https://www.iarc.fr/faq/latest-global-cancer-data-2020-qa/.
GARRIDO-CASTRO A C, LIN N U, POLYAK K. Insights into Molecular Classifications of Triple-Negative Breast Cancer: Improving Patient Selection for Treatment [J]. Cancer Discov, 2019, 9(2): 176–98.
SUN X, LUO H, HAN C, et al. Identification of a Hypoxia-Related Molecular Classification and Hypoxic Tumor Microenvironment Signature for Predicting the Prognosis of Patients with Triple-Negative Breast Cancer [J]. Frontiers in Oncology, 2021, 11(3253):
BIANCHINI G, BALKO J M, MAYER I A, et al. Triple-negative breast cancer: challenges and opportunities of a heterogeneous disease [J]. Nat Rev Clin Oncol, 2016, 13(11): 674–90.
LEHMANN B D, BAUER J A, CHEN X, et al. Identification of human triple-negative breast cancer subtypes and preclinical models for selection of targeted therapies [J]. J Clin Invest, 2011, 121(7): 2750–67.
WU S Y, WANG H, SHAO Z M, et al. Triple-negative breast cancer: new treatment strategies in the era of precision medicine [J]. Sci China Life Sci, 2021, 64(3): 372–88.
BURSTEIN M D, TSIMELZON A, POAGE G M, et al. Comprehensive genomic analysis identifies novel subtypes and targets of triple-negative breast cancer [J]. Clin Cancer Res, 2015, 21(7): 1688–98.
HE Y, JIANG Z, CHEN C, et al. Classification of triple-negative breast cancers based on Immunogenomic profiling [J]. Journal of Experimental & Clinical Cancer Research, 2018, 37(1): 327.
JIANG Y Z, MA D, SUO C, et al. Genomic and Transcriptomic Landscape of Triple-Negative Breast Cancers: Subtypes and Treatment Strategies [J]. Cancer Cell, 2019, 35(3): 428 – 40.e5.
XIAO Y, MA D, ZHAO S, et al. Multi-Omics Profiling Reveals Distinct Microenvironment Characterization and Suggests Immune Escape Mechanisms of Triple-Negative Breast Cancer [J]. Clinical Cancer Research, 2019, 25(16): 5002.
JIA H, TRUICA C I, WANG B, et al. Immunotherapy for triple-negative breast cancer: Existing challenges and exciting prospects [J]. Drug Resistance Updates, 2017, 32(1–15.
MCARTHUR H L, PAGE D B. Immunotherapy for the treatment of breast cancer: checkpoint blockade, cancer vaccines, and future directions in combination immunotherapy [J]. Clin Adv Hematol Oncol, 2016, 14(11): 922–33.
SCHMID P, CORTES J, PUSZTAI L, et al. Pembrolizumab for Early Triple-Negative Breast Cancer [J]. N Engl J Med, 2020, 382(9): 810–21.
SCHMID P, ADAMS S, RUGO H S, et al. Atezolizumab and Nab-Paclitaxel in Advanced Triple-Negative Breast Cancer [J]. N Engl J Med, 2018, 379(22): 2108–21.
MARRA A, VIALE G, CURIGLIANO G. Recent advances in triple negative breast cancer: the immunotherapy era [J]. BMC Med, 2019, 17(1): 90.
QIU P, GUO Q, YAO Q, et al. Characterization of Exosome-Related Gene Risk Model to Evaluate the Tumor Immune Microenvironment and Predict Prognosis in Triple-Negative Breast Cancer [J]. Front Immunol, 2021, 12(736030.
KENNEDY L B, SALAMA A K S. A review of cancer immunotherapy toxicity [J]. CA Cancer J Clin, 2020, 70(2): 86–104.
HANAHAN D, COUSSENS L M. Accessories to the crime: functions of cells recruited to the tumor microenvironment [J]. Cancer Cell, 2012, 21(3): 309–22.
BYRNE A, SAVAS P, SANT S, et al. Tissue-resident memory T cells in breast cancer control and immunotherapy responses [J]. Nat Rev Clin Oncol, 2020, 17(6): 341–8.
WANG S, XIONG Y, ZHANG Q, et al. Clinical significance and immunogenomic landscape analyses of the immune cell signature based prognostic model for patients with breast cancer [J]. Brief Bioinform, 2021, 22(4):
HäNZELMANN S, CASTELO R, GUINNEY J. GSVA: gene set variation analysis for microarray and RNA-seq data [J]. BMC Bioinformatics, 2013, 14(7.
WITTEN D M, TIBSHIRANI R. A framework for feature selection in clustering [J]. J Am Stat Assoc, 2010, 105(490): 713–26.
YOSHIHARA K, SHAHMORADGOLI M, MARTíNEZ E, et al. Inferring tumour purity and stromal and immune cell admixture from expression data [J]. Nat Commun, 2013, 4(2612.
NEWMAN A M, LIU C L, GREEN M R, et al. Robust enumeration of cell subsets from tissue expression profiles [J]. Nat Methods, 2015, 12(5): 453–7.
TIBSHIRANI R. The lasso method for variable selection in the Cox model [J]. Stat Med, 1997, 16(4): 385–95.
ZHENG S, ZOU Y, LIANG J Y, et al. Identification and validation of a combined hypoxia and immune index for triple-negative breast cancer [J]. Mol Oncol, 2020, 14(11): 2814–33.
CHAROENTONG P, FINOTELLO F, ANGELOVA M, et al. Pan-cancer Immunogenomic Analyses Reveal Genotype-Immunophenotype Relationships and Predictors of Response to Checkpoint Blockade [J]. Cell Rep, 2017, 18(1): 248–62.
CHIC N, BRASó-MARISTANY F, PRAT A. Biomarkers of immunotherapy response in breast cancer beyond PD-L1 [J]. Breast Cancer Research and Treatment, 2022, 191(1): 39–49.
EMENS L A. Breast Cancer Immunotherapy: Facts and Hopes [J]. Clinical Cancer Research, 2018, 24(3): 511.
BAXEVANIS C N, FORTIS S P, PEREZ S A. The balance between breast cancer and the immune system: Challenges for prognosis and clinical benefit from immunotherapies [J]. Seminars in Cancer Biology, 2021, 72(76–89.
BRUNI D, ANGELL H K, GALON J. The immune contexture and Immunoscore in cancer prognosis and therapeutic efficacy [J]. Nat Rev Cancer, 2020, 20(11): 662–80.
METZGER-FILHO O, TUTT A, DE AZAMBUJA E, et al. Dissecting the heterogeneity of triple-negative breast cancer [J]. J Clin Oncol, 2012, 30(15): 1879–87.
XU M, LI Y, LI W, et al. Immune and Stroma Related Genes in Breast Cancer: A Comprehensive Analysis of Tumor Microenvironment Based on the Cancer Genome Atlas (TCGA) Database [J]. Front Med (Lausanne), 2020, 7(64.
LI N, WANG J, ZHAN X. Identification of Immune-Related Gene Signatures in Lung Adenocarcinoma and Lung Squamous Cell Carcinoma [J]. Front Immunol, 2021, 12(752643.
WANG X, SU W, TANG D, et al. An Immune-Related Gene Prognostic Index for Triple-Negative Breast Cancer Integrates Multiple Aspects of Tumor-Immune Microenvironment [J]. Cancers (Basel), 2021, 13(21):
WEINSTEIN A M, STORKUS W J. Therapeutic Lymphoid Organogenesis in the Tumor Microenvironment [J]. Adv Cancer Res, 2015, 128(197–233.
GU-TRANTIEN C, LOI S, GARAUD S, et al. CD4â࿽º follicular helper T cell infiltration predicts breast cancer survival [J]. J Clin Invest, 2013, 123(7): 2873–92.
CRISCITIELLO C, BAYAR M A, CURIGLIANO G, et al. A gene signature to predict high tumor-infiltrating lymphocytes after neoadjuvant chemotherapy and outcome in patients with triple-negative breast cancer [J]. Ann Oncol, 2018, 29(1): 162–9.
RAZIS E, KALOGERAS K T, KOTSANTIS I, et al. The Role of CXCL13 and CXCL9 in Early Breast Cancer [J]. Clin Breast Cancer, 2020, 20(1): e36-e53.
PINTO J A, ARAUJO J, CARDENAS N K, et al. A prognostic signature based on three-genes expression in triple-negative breast tumours with residual disease [J]. npj Genomic Medicine, 2016, 1(1): 15015.
DENKERT C, VON MINCKWITZ G, BRASE J C, et al. Tumor-infiltrating lymphocytes and response to neoadjuvant chemotherapy with or without carboplatin in human epidermal growth factor receptor 2-positive and triple-negative primary breast cancers [J]. J Clin Oncol, 2015, 33(9): 983–91.
ARAUJO J M, GOMEZ A C, AGUILAR A, et al. Effect of CCL5 expression in the recruitment of immune cells in triple negative breast cancer [J]. Sci Rep, 2018, 8(1): 4899.
ARABPOUR M, RASOLMALI R, TALEI A R, et al. Granzyme B production by activated B cells derived from breast cancer-draining lymph nodes [J]. Mol Immunol, 2019, 114(172–8.
SHIBUTANI M, MAEDA K, NAGAHARA H, et al. A comparison of the local immune status between the primary and metastatic tumor in colorectal cancer: a retrospective study [J]. BMC Cancer, 2018, 18(1): 371.
BEDOGNETTI D, HENDRICKX W, MARINCOLA F M, et al. Prognostic and predictive immune gene signatures in breast cancer [J]. Curr Opin Oncol, 2015, 27(6): 433–44.
SHI Z, SHEN J, QIU J, et al. CXCL10 potentiates immune checkpoint blockade therapy in homologous recombination-deficient tumors [J]. Theranostics, 2021, 11(15): 7175–87.
GUI C P, WEI J H, CHEN Y H, et al. A new thinking: extended application of genomic selection to screen multiomics data for development of novel hypoxia-immune biomarkers and target therapy of clear cell renal cell carcinoma [J]. Brief Bioinform, 2021, 22(6):

No competing interests reported.

Download PDF

Journal Publication

published 23 May, 2023

Read the published version in Clinical and Experimental Medicine →

Editorial decision: Major revision
17 Apr, 2023
Reviews received at journal
14 Apr, 2023
Reviewers agreed at journal
03 Apr, 2023
Reviewers invited by journal
01 Apr, 2023
Submission checks completed at journal
31 Mar, 2023
Editor assigned by journal
31 Mar, 2023
First submitted to journal
24 Mar, 2023

You are reading this latest preprint version

A novel model associated with tumor microenvironment on predicting prognosis and immunotherapy in Triple Negative Breast Cancer

Status:

Journal Publication

Version 1

Abstract

Background

Methods

Results

Conclusion

Figures

Background

Materials And Methods

Study Design

Database Collection

Landscape Of Immune Infiltration In Tumor Microenvironment

Degs Between Subgroups And Construction Of Risk Score Model

The Immune-related And Prognostic Significance Of Degs

Validation Of Risk Score Model

Detection Of The Fractions Of Immune Cells And Expression Of Degs

Correlation Between Immune Checkpoint Genes And Ips

Gene Set Enrichment Analysis

Statistical analysis

Results

Landscape of immune cell infiltration of TME in TNBC

The Immune Status And Prognostic Difference Between Two Subgroups

Screening Of Degs And Construction Of Risk Score Model

The Correlation Of Degs With Immune Infiltration And Prognosis

Validation Of Risk Score Model Based On Degs In Tnbc

Predictive Significance For Immunotherapy And Immune-related Pathways

Discussion

Conclusion

Abbreviations

Declarations

References

Additional Declarations

Supplementary Files

Status:

Journal Publication

Version 1