Article Text

Original research
Gene panel predicts neoadjuvant chemoimmunotherapy response and benefit from immunotherapy in HER2-negative breast cancer
  1. Xunxi Lu1,2,
  2. Zongchao Gou3,4,
  3. Hong Chen2,
  4. Li Li2,
  5. Fei Chen2,
  6. Chunjuan Bao2 and
  7. Hong Bu1,2
  1. 1Department of Pathology, West China Hospital, Sichuan University, Chengdu, Sichuan, China
  2. 2Institute of Clinical Pathology, West China Hospital, Sichuan University, Chengdu, Sichuan, China
  3. 3Department of General Surgery, West China Hospital, Sichuan University, Chengdu, Sichuan, China
  4. 4Breast Disease Center, West China Hospital, Sichuan University, Chengdu, Sichuan, China
  1. Correspondence to Dr Hong Bu; hongbu{at}scu.edu.cn

Abstract

Background It is encountering the dilemma of lacking precise biomarkers to predict the response to neoadjuvant chemoimmunotherapy (NACI) and determine whether patients should use immune checkpoint inhibitors (ICIs) in early breast cancer (BC). We aimed to develop a gene signature to predict NACI response for BC patients and identify individuals suitable for adding ICIs.

Patients and methods Two I-SPY2 cohorts and one West China Hospital cohort of patients treated with NACI were included. Machine learning algorithms were used to identify key genes. Principal component analysis was used to calculate the ImPredict (IP) score. The interaction effects between biomarkers and treatment regimens were examined based on the logistic regression analysis. The relationship between the IP score and immune microenvironment was investigated through immunohistochemistry (IHC) and multiplex IHC.

Results The area under the curves of the IP score were 0.935, 0.865, and 0.841 in the discovery cohort, validation cohort 1, and in-house cohort. Marker-treatment interaction tests indicated that the benefits from immunotherapy significantly varied between patients with high and low IP scores (p for interaction <0.001), and patients with high IP scores were more suitable for immunotherapy addition.

Conclusions Our IP model shows favorable performance in predicting NACI response and is an effective tool for identifying BC patients who will benefit from ICIs. It may help clinicians optimize treatment strategies and guide clinical decision-making.

  • Breast Cancer
  • Immunotherapy
  • Immune Checkpoint Inhibitor
  • Biomarker
  • Neoadjuvant

Data availability statement

Data are available in a public, open access repository. Data are available upon reasonable request. All data relevant to the study are included in the article or uploaded as supplementary information.

Data availability statement

The raw sequencing data of the WCH-TNBC cohort in this study has been deposited to NCBI (Project ID: PRJNA1081147, SRA ID: SRP492770). Data generated by the TCGA and Gene Expression Omnibus (GEO) are available publicly from the https://portal.gdc.cancer.gov/ and https://www.ncbi.nlm.nih.gov/geo/ websites. Other data supporting the findings of this study are available within the article, its supplemental files, and from the corresponding author upon reasonable request.

http://creativecommons.org/licenses/by-nc/4.0/

This is an open access article distributed in accordance with the Creative Commons Attribution Non Commercial (CC BY-NC 4.0) license, which permits others to distribute, remix, adapt, build upon this work non-commercially, and license their derivative works on different terms, provided the original work is properly cited, appropriate credit is given, any changes made indicated, and the use is non-commercial. See http://creativecommons.org/licenses/by-nc/4.0/.

Statistics from Altmetric.com

Request Permissions

If you wish to reuse any or all of this article please use the link below which will take you to the Copyright Clearance Center’s RightsLink service. You will be able to get a quick price and instant permission to reuse the content in many different ways.

WHAT IS ALREADY KNOWN ON THIS TOPIC

  • Immune checkpoint inhibitors (ICIs) have revolutionized the treatment of breast cancer (BC) patients in the past few years. However, the clinical benefit of ICIs varies greatly among individuals. Current biomarkers are insufficient to predict the immunotherapy response for BC patients, especially for early BC. For example, both the KEYNOTE-522 and IMpassion031 show that the benefit gained from ICIs is not correlated with programmed cell death ligand 1 (PD-L1) positivity, a commonly used biomarker in clinical practice. Hence, developing precise predictive biomarkers remains an urgent medical need.

WHAT THIS STUDY ADDS

  • The ImPredict (IP) score shows robust performance for predicting neoadjuvant chemoimmunotherapy response in two cohorts of the I-SPY2 clinical trial and one in-house cohort. Importantly, the IP score is a promising tool to identify patients who are recommended for ICI combination therapy. Moreover, the IP model characterizes the immune microenvironment and provides additional prognostic information for triple-negative BC.

HOW THIS STUDY MIGHT AFFECT RESEARCH, PRACTICE OR POLICY

  • In clinical practice, PD-L1 fails to determine whether to use ICI drugs during the neoadjuvant chemotherapy phase in early BC. Our work proposes a novel and efficient model for the optimization of therapeutic strategies and precision oncology for early BC patients. Patients with high IP scores are strong candidates for ICI therapy combined with chemotherapy.

Introduction

Breast cancer (BC) is the most prevalent malignant tumor among women globally and is classified into hormone receptor-positive/human epidermal growth factor receptor 2-negative (HR+/HER2−), HR+/HER2+, HR−/HER2+, and triple-negative breast cancer (TNBC).1 2 BC was previously considered “weakly immunogenic” among solid tumors. In recent years, omics researchers have uncovered that TNBC exhibits higher chromosomal instability, frequent copy number alterations, and complex structural rearrangements than the other BC subtypes.3 4 It also has more lymphocyte infiltration and a higher programmed cell death ligand 1 (PD-L1) positivity rate, thus immunotherapy is mainly tested and applied in TNBC at first.5 The Food and Drug Administration has approved pembrolizumab for the treatment of early high-risk TNBC and PD-L1 Combined Positive Score (CPS) ≥10 advanced or metastatic TNBC. Besides, some ongoing clinical trials, such as the CheckMate 7FL and KEYNOTE-756 studies, suggest that early high-risk HR+/HER2− BC can also benefit from immune checkpoint inhibitors (ICIs).6 7

Although immunotherapy has shown potential in treating BC, it is encountering the dilemma of lacking precise predictive biomarkers for ICI treatment efficacy. One significant issue is that only part of patients respond to ICI, and it needs urgent solutions in patient selection for ICI therapy, especially for early BC.8 PD-L1 immunohistochemistry (IHC) is the most extensively adopted biomarker in clinical practice. Nevertheless, it remains controversial due to the high heterogeneity, different detective standards, and various positive criteria.9 Moreover, KEYNOTE-522 and IMpassion031 studies have demonstrated that early TNBC patients could benefit from ICI, irrespective of the PD-L1 positivity.10 11 Some other biomarkers, such as tumor mutation burden (TMB) and microsatellite instability, are not appropriate biomarkers for BC. Evidence for these two markers mainly focuses on solid tumors such as lung cancer and colorectal cancer.12 13 It is rare for BC patients to exhibit TMB ≥10 mut/Mb or be detected as mismatch repair deficiency.14 15

Some studies have attempted to develop molecular models like IFNγ-6, GEP, and Pan_F_TBRs in melanoma and urothelial cancer.16 17 Given the complexity and heterogeneity of the tumor microenvironment (TME), these molecular markers were not necessarily suitable for BC. Considering the lack of ideal biomarkers for immunotherapy in early BC, we constructed a novel gene model for predicting ICI response, ImPredict (IP), using the gene array data from the early HER2− BC immunotherapy cohort, and validated it in two independent cohorts. IP model can predict the pathological complete response (pCR) outcomes of HER2− BC patients who receive neoadjuvant chemoimmunotherapy (NACI) and is expected to screen patients suitable for ICI addition. Furthermore, we investigated the relationship between the IP score and the TME of TNBC and conducted an exploratory analysis to figure out whether the IP score could provide additional prognostic information for TNBC.

Methods

Data collection

The discovery cohort (N=69; TNBC: N=29, HR+: N=40) was derived from the pembrolizumab arm in the I-SPY2 clinical trial,18 with gene array data and clinical information downloaded from GSE194040. Validation cohort 1 (N=71; TNBC: N=21, HR+: N=50) was the durvalumab/olaparib (DO) arm of the I-SPY2 clinical trial,19 and its gene chip data and clinical information were downloaded from GSE173819. Cohort A for the marker-treatment interaction test included data from 246 patients (neoadjuvant chemotherapy (NACT) alone: N=177, NACT plus pembrolizumab: N=69), and cohort B included data from 105 patients (NACT alone: N=34, NACT plus DO: N=71). The TNBC cohort of The Cancer Genome Atlas (TCGA) (N=147) was obtained from the GDC Data Portal. The immune phenotype cohort (N=43) was downloaded from GSE177043.20

Validation cohort 2 (N=55) and TNBC RNA-sequencing (RNA-seq) cohort (N=53) were collected from the West China Hospital (WCH). Validation cohort 2 included TNBC patients treated with NACI at WCH from January 2020 to June 2023. Pathological responses were reviewed, and pCR is defined as no invasive cancer in both the primary lesion (ductal carcinoma in situ may be present) and regional lymph nodes (ypT0/Tis ypN0). The RNA-seq cohort was collected from the treatment-naïve TNBC patients before surgery in WCH from January 2022 to October 2022.

Differentially expressed genes analysis

Differentially expressed genes (DEGs) were identified using the R package “limma” according to the cut-off of |Log2FC|>0.5 and adjusted p value <0.05. The adjusted p value was calculated by the Benjamini-Hochberg correction.

Weighted gene correlation network analysis

The gene co-expression network of discovery cohort was generated by the R package “WGCNA.”21 For each pair of genes, their similarity was calculated based on the expression patterns. The similarity matrix was converted into an adjacency matrix by raising the similarity to a positive power β (soft threshold) to satisfy the scale-free property. Then, the weighted adjacency matrix was transformed into a topological overlap matrix (TOM), and accordingly, the dissimilarity (1−TOM) was introduced to generate the systematic clustering tree between genes. Furthermore, module identification was conducted by a hybrid dynamic shear tree with the criterion of 30 as the minimum number of genes in each gene module. Relationships between modules and sample traits were analyzed by calculating the correlation between the traits and module eigengenes, which are the first principal components of the expression data of all genes in the modules and are considered representative of the modules. Finally, the gene modules with p value <0.05 were recognized as those that were significantly associated with pCR, and the members in these modules were selected for candidate genes.

Identification of key genes associated with pCR

To determine the key genes associated with pCR, we performed gene selection using four machine learning algorithms. Boruta is an algorithm employed for feature selection based on the Random Forest classification. We obtained the gene importance using the “Boruta” function embedded in the “Boruta” R package with maximum iterations up to 300 times, and then confirmed variables were selected as key genes for model construction. Least absolute shrinkage and selection operator (LASSO) regression incorporates the L1 regularization term, realizing sparsity by limiting the sum of the absolute values of the parameter vectors. This enables LASSO to carry out feature selection by estimating numerous parameters to be zero. We choose the variables with non-zero parameters as the model genes. In this study, we conducted LASSO regression with the “glmnet” R package. The optimal value of λ was determined using 10-fold cross-validation. In the end, λ was selected based on a one-SE criterion. Support Vector Machine-Recursive Feature Elimination (SVM-RFE) is a sequential backward selection algorithm based on the maximum margin principle of SVM. It trains the model with samples, then scores and ranks each feature, removes the feature with the lowest score, trains the model again with the remaining features, iterates the process, and finally selects the required number of features. SVM-RFE was conducted using the “e1071” R package and codes from https://github.com/johncolby/SVM-RFE. eXtreme Gradient Boosting (XGBoost) is a type of gradient-boosting algorithm, which builds a model by progressively constructing numerous decision trees. XGBoost assigns an importance score for every feature in the input data, aiding in identifying the most significant features in the model. The “xgboost” R package was used to achieve key gene selection.

IP signature construction

First, we analyzed the DEGs between pCR and residual disease (RD) samples and identified gene modules associated with pCR through weighted gene correlation network analysis (WGCNA). Next, the overlapping genes between DEGs and WGCNA genes were retained for subsequent analysis. Four machine learning algorithms were conducted to select the key genes correlated with pCR, respectively. Genes highly expressed in pCR samples were defined as gene set P (positive genes), and genes highly expressed in RD samples were defined as gene set N (negative genes). The P and N gene sets were extracted and their expression was then z-score transformed. Finally, principal component analysis (PCA) was performed and principal component 1 (PC1) was calculated to serve as the gene signature score. This approach has the advantage of focusing the score for the set on the largest block of well-correlated (or anti-correlated) genes in the gene set while down-weighting contributions from genes that do not track with other set members.17 IP score was defined using a method similar to the Genomic Grade Index22:

Embedded Image

where PC1P and PC1N represent the first principal component of gene set P and gene set N, respectively.

Transcriptomic-based predictive biomarkers

PD-L1 is the gene expression of CD274. IFNγ-6 gene signature was an average gene expression-based biomarker. T cell-inflamed GEP score was calculated as the weighted sum of signature gene expressions.16 Antigen-presenting machinery (APM) signature was calculated using the single-sample gene set enrichment analysis (ssGSEA),23 and pan-fibroblast TGF-β response signature (Pan-F-TBRs) was calculated based on PCA using the “calculate_sig_score” function implemented in the “IOBR” R package.24 The detailed information of these signatures was summarized in online supplemental table S3.

Supplemental material

RNA extraction

RNA was isolated from the pre-treatment FFPE specimen using the miRNeasy FFPE Kit (QIAGEN, Cat#217504) according to the manufacturer”s protocol. The quality and amount of RNA were assessed by NanoDrop OneC (Thermo Fisher Scientific). Specimens with adequate RNA which had an optical density of 1.8≤260/280 ratio ≤2.1 and a total yield ≥3 µg were eligible for reverse transcription.

Real-time quantitative PCR

We evaluated the IP score of the WCH immunotherapy cohort by quantifying the expression of 36 signature genes using the real-time quantitative PCR assay. Reverse transcription was performed using the HiScript III 1st Strand cDNA Synthesis Kit (R312, Vazyme). For expression analysis, all genes were detected in technical triplicates using the Taq Pro Universal SYBR qPCR Master Mix (Q712, Vazyme). QuantStudio 6 Flex Real-Time PCR System (Applied Biosystems) was applied for RNA quantification. Relative gene expression levels were quantified using the 2−ΔCt method with the mean level of GAPDH and ACTB as an internal reference. Primers for 36 model genes and 2 reference genes were listed in online supplemental table S7.

Statistics analysis

The receiver operating characteristic (ROC) curve was visualized and the area under the curve (AUC) was calculated to evaluate the predictive performance of biomarkers for therapy response using the “pROC” R package. The CI of AUC was obtained by the bootstrap resample embedded in the “roc” function of the “pROC” package. Marker-treatment interaction analysis was performed for each biomarker, which is a test for interaction based on the logistic regression models. A significant interaction test indicates a differential treatment effect according to the biomarker. The Kaplan-Meier method was applied to generate survival curves, and the log-rank test was used to determine the statistical significance of differences. Cox proportional hazard regression models were used to estimate the HRs with 95% CIs for variables. Time-dependent ROC analysis was performed using the R package “timeROC.” The nomogram model was generated via the R package “rms.” Statistical difference in the distribution of continuous variables between two groups was compared by the Wilcoxon rank-sum test, and that of three or more groups was examined by the Kruskal-Wallis test. Correlation coefficients were calculated by Spearman”s correlation analyses. The χ2 test or Fisher”s exact test was used to analyze categorical variables distribution. A two-sided p value <0.05 was considered statistically significant. All statistical analyses were performed using SPSS (V.25.0) and R software (V.4.1.1).

Results

Construction of the IP model

The overall design of this study was displayed in online supplemental figure S1. In the discovery cohort, 462 DEGs were identified between pCR (N=31) and RD ((N=38) patients with the cut-off of adj. p value <0.05 and |Log2FC|>0.5, including 265 downregulated and 197 upregulated genes in the pCR group (figure 1A). In the WGCNA process, the soft threshold β was set to 9 to construct a gene co-expression network (online supplemental figure S2A), and then 29 gene modules were determined, as represented by different colors. The heatmap revealed the eigengenes adjacency of these modules (online supplemental figure S2B). Subsequently, the correlations between gene modules and clinical traits, including MammaPrint status and pCR, were analyzed. Module-trait relationship indicated that five modules (dark red, green, gray, pink, salmon) were positively correlated with pCR, and five modules (yellow, midnight blue, light yellow, magenta, grey60) were negatively correlated with pCR according to the criterion of p value <0.05 (figure 1B). Next, 405 overlapping genes between 462 DEGs and 7114 WGCNA genes were extracted as candidates (figure 1C). GO analysis revealed that candidate genes highly expressed in the pCR group were enriched in immune response, lymphocyte activation, cytokine production, etc, whereas those highly expressed in the RD group were enriched in hormone regulation, cell migration, metabolic process, etc (online supplemental figure S3A, B). Furthermore, four different machine learning algorithms were used to reduce noise or redundant genes, and the retained genes were subjected to a PCA-based procedure to develop gene signatures. By comparing the AUCs of the four gene sets derived from four algorithms, the final prediction signature was determined (figure 1C).

Supplemental material

Figure 1

The construction process of IP signature. (A) Volcano plot of DEGs in the pCR group versus the RD group (|Log2FC|>0.5; p<0.05). (B) Relationships between the module eigengenes and clinical traits. The red and blue boxes indicate modules positively and negatively related to pCR, respectively (p<0.05). (C) IP signature is built using the key genes identified by four machine-learning algorithms. DEGs, differentially expressed genes; IP, ImPredict; LASSO, least absolute shrinkage and selection operator; pCR, pathological complete response; RD, residual disease; SVM-RFE, Support Vector Machine-Recursive Feature Elimination; WGCNA, weighted gene correlation network analysis; XGBoost, eXtreme Gradient Boosting.

After gene selection, 36, 29, 25, and 74 genes were identified by the Boruta, XGBoost, LASSO, and SVM-RFE (online supplemental figure S4A–D, online supplemental table S1), respectively. For all patients, TNBC subgroup, and HR+ subgroup in the discovery cohort, all four signatures exhibited significantly higher scores in the pCR group compared with the RD group (online supplemental figure S5A–D). In the validation cohort 1, only the signature score which was based on the 36 genes selected by the Boruta algorithm was significantly higher in the pCR group than in the RD group, regardless of the patient subgroups. Whereas, no significant difference was observed in the other three signature scores between the pCR and RD groups of TNBC patients (online supplemental figure S6A–D). Besides, only the Boruta signature score reached an AUC of above 0.8 in every subgroup of the discovery cohort (all: 0.935; TNBC: 0.989; HR+: 0.890) and validation cohort 1 (all: 0.865; TNBC: 0.833; HR+: 0.895), demonstrating that it is the most robust one among the four signatures (online supplemental figure S7A–D, online supplemental table S2). Therefore, we adopted the Boruta signature as the final IP model.

Predictive performances of the IP score and five gene expression-based signatures

We observed that TNBC exhibited a higher IP score than HR+ BC, which is consistent with the higher response rate of TNBC (online supplemental figure S8A,B). As illustrated in the bar plots, pCR patients showed a higher IP score than RD patients in the discovery cohort (figure 2A) and validation cohort 1 (figure 2B). In recent years, several predictive biomarkers for immunotherapy have been developed. To compare their prediction performances with the IP score, we calculated these biomarker scores based on the descriptions in the original publications (online supplemental table S3). ROC analysis displayed that the discrimination ability of the IP score was superior to PD-L1 mRNA, Pan_F_TBRs, APM, IFNγ-6, and GEP (figure 2C–F, online supplemental figure S9A, B, online supplemental table S4). Furthermore, we collected the data of STAT1, Mitotic, and ESR1_PGR_ave signatures in the validation cohort 1, which were considered to be associated with pCR in the biomarker analysis of the original citation. The IP score showed superior performance than the Mitotic signature and tended to perform better than the STAT1 and ESR1_PGR_ave signatures in all patients (online supplemental figure S10A–C, online supplemental table S5). Besides, DCA demonstrated that the IP score achieved higher overall net benefit compared with the other signatures (online supplemental figure S11A–D). Taken together, these results suggest that the IP score is a more excellent biomarker for immunotherapy response prediction.

Figure 2

Predictive performances of IP score and five biomarkers in the discovery set and validation set 1. (A–B) Distribution of the IP score in the discovery set (A) and validation set 1 (B). (C–D) ROC curves of the IP score and five biomarkers in the discovery set (C) and validation set 1 (D). (E–F) AUCs of the IP score and five biomarkers in the discovery set (E) and validation set 1 (F). APM, antigen-presenting machinery; AUC, area under the curve; HR, hormone receptor; IP, ImPredict; pCR, pathological complete response; PD-L1, programmed cell death ligand 1; RD, residual disease; ROC, receiver operating characteristic; TNBC, triple-negative breast cancer.

Validation in an in-house immunotherapy cohort

To further test the predictive power of the IP signature, we next collected a validation cohort 2 including 55 patients from WCH, all of whom were diagnosed with TNBC and treated with NACI, and their baseline information was summarized in online supplemental table S6. The distribution of IP scores in the validation cohort 2 is displayed in figure 3A. Consistent with the above results, pCR patients showed a significantly higher IP score than the RD patients (figure 3B), and there was a rising trend for the IP score as the Miller-Payne grading increased (figure 3C, Kruskal-Wallis test, p=6.3e-05). Moreover, ROC analysis demonstrated that the AUC of the IP score reached 0.841 (figure 3D). We assessed the predictive power of PD-L1 CPS in 24 patients. The AUC of PD-L1 CPS was 0.750, which was inferior to the IP score (figure 3E). Univariate logistic regression analysis indicated that the IP score, Ki-67, and tumor stage were statistically significant for the pCR outcome (figure 3F). Then the three variables were input into the multivariate logistic regression analysis, which suggested that the IP score is an independent predictor for pCR (figure 3G).

Figure 3

Validation of the IP score in an in-house immunotherapy cohort. (A) Distribution of the IP score in the WCH immunotherapy cohort. (B) pCR patients show a significantly higher IP score than RD patients (Wilcoxon rank-sum test, p<0.0001). (C) The IP score increases with the Miller-Payne upgrade (Kruskal-Wallis test, p=6.3e-05). (D–E) ROC curves of the IP score (D) and PD-L1 CPS (E) in the WCH immunotherapy cohort (IP: N=55, PD-L1: N=24). (F) Univariate logistic regression analysis indicates that the IP score, Ki-67, and tumor stage are statistically significant predictors for the pCR. (G) Multivariate logistic regression analysis including IP score, Ki-67, and tumor stage. AUC, area under the curve; CPS, Combined Positive Score; IP, ImPredict; pCR, pathological complete response; PD-L1, programmed cell death ligand 1; RD, residual disease; ROC, receiver operating characteristic; WCH, West China Hospital.

IP score is a promising biomarker to identify patients benefiting from ICI

To explore whether the IP score could predict the benefit of ICI addition, we collected the corresponding NACT alone arms of the discovery and validation set 1 in the I-SPY2 trial. The scores of IP and five molecular signatures were calculated, and for each signature, the “high” versus “low” groups were divided by the median of signature score or the optimal cut-off for predicting pCR in the immunotherapy arm. Further, we examined the interaction between treatment and biomarkers through the logistic regression analysis which incorporated the treatment-by-biomarker interaction term for pCR. In cohort A grouped by the median value of the IP score, the pCR rate of the NACT+pembro arm was significantly higher than that of the NACT arm in the IP-high group (overall: 81.8% vs 20.0%, p<0.001; TNBC: 79.2% vs 18.3%, p<0.001; HR+: 88.9% vs 26.3%, p<0.01), whereas there was no significant increase in the pCR rate of the NACT+pembro arm compared with the NACT arm in IP-low group (online supplemental figure S12A). In cohort A grouped by the median values of the other markers, the NACT+pembro arm exhibited an increasing trend in pCR rates than the NACT arm in both high and low groups (online supplemental figure S12B–E). Notably, the marker-treatment interaction test was significant for the IP score but not for PD-L1, APM, Pan_F_TBRs, IFNγ-6, and GEP (figure 4A). Besides, we also divided cohort A according to the optimal cut-off value of each signature in the NACT+pembro arm and performed analyses as above. Among these biomarkers, the difference in pCR benefit brought by pembrolizumab was most pronounced in the IP-high group (online supplemental figure 13A–E), and the p value of the IP score in the marker-treatment interaction test was the most significant (figure 4B, IP: p< 0.001, PD-L1: p=0.016, Pan_F_TBRs: p=0.719, APM: p=0.054, IFNγ-6: p=0.033, GEP: p=0.076). Moreover, these analyses were repeated in cohort B. It was discovered that regardless of the grouping way, the NACT+DO arm exhibited a significantly higher pCR rate than the NACT arm in the IP-high group, but not in the IP-low group (online supplemental figure S14A). However, in both high and low groups divided by the other biomarkers, the pCR rates were higher in the NACT+DO arm than in the NACT arm (online supplemental figure S14B–E). In the marker-treatment interaction test, the p value of the IP score remained the most significant, though it did not meet the 0.05 standard, which may be due to the limitation of the sample size (online supplemental figure S15A, B). To sum up, the addition of immunotherapy to NACT increases the pCR rates in IP-high patients but not the IP-low patients and the IP score is more effective than the other biomarkers in screening patients suitable for immunotherapy addition.

Figure 4

Marker-treatment interaction test in cohort A. (A) Marker-treatment interaction test when the patients are grouped by the medians of signature scores in cohort A. (B) Marker-treatment interaction test when the patients are grouped by the optimal cut-off of predicting pCR in the immunotherapy arm of cohort A. aORs are log10 transformed. bORs are HR-adjusted and log10 transformed. APM, antigen-presenting machinery; IP, ImPredict; NACT, neoadjuvant chemotherapy; PD-L1, programmed cell death ligand 1.

Associations of the IP score with immune microenvironment in TNBC

We examined the relationships between the IP score and extrinsic/intrinsic immune escape mechanisms of the TNBC tumors using the TCGA-TNBC (N=147) and WCH-TNBC (N=53) cohorts. The high and low IP groups were divided by the IP score median in these two cohorts. The IP-high group exhibited low tumor purity, prompting the diversity of TME in patients with high IP scores (online supplemental figure S16A). Both the CYT score and ImmuneScore were positively correlated with the IP score, demonstrating an activated immune response and more effector immune cell infiltration in IP-high patients (online supplemental figure S16B). Additionally, the IP-high group had a higher abundance of CD8+ T cells, B cells, NK cells, and mono/myeloid cells, while the IP-low group was infiltrated with more immunosuppression cells, including neutrophils, endothelial cells, and fibroblasts (online supplemental figure S16C).25 In line with these results, cytokine expression analysis revealed that numerous members of chemokines, interleukins, interferons, and receptors were highly expressed in the IP-high group, whereas an increased secretion of immunoinhibitory cytokines was observed in the IP-low group, such as platelet-derived growth factor A and D, transforming growth factor-β3, vascular endothelial growth factor A and C, etc (online supplemental figure S16D). TMB was reported to be associated with immunotherapy response in solid tumors, but there was no significant correlation between TMB and IP score in the TCGA-TNBC cohort (online supplemental figure S17A). Two innate immunity-sensing pathways, cGAS-STING and NLRP3 inflammasome, were more active in the IP-high group (online supplemental figure S17B). Almost all histocompatibility complex (MHC) molecules as well as co-stimulatory and co-inhibitory molecules showed higher expression in the IP-high group (online supplemental figure S17C, D). T-cell receptor and B-cell receptor (BCR) clonalities were significantly increased in the IP-high group, in which BCR evenness was significantly lower (online supplemental figure S17E). GSEA demonstrated that the IP-high group was enriched with pathways related to immune response, (figure 5A), which were corroborated by the GSEA based on KEGG and Reactome pathways (online supplemental figure S18A, B). We further validated the above discoveries in clinical specimens. H&E-stained pathological sections examination indicated that the IP-high TNBC was infiltrated with more stromal and intratumoral tumor-infiltrating lymphocytes (online supplemental figure S19A–D). We next performed IHC staining for representative markers in the WCH-TNBC cohort (N=52) and observed that the IP-high group showed weaker CDH2, CD31, CA9 expression and stronger MHC-I, CD8, CD4, CD20 levels (figure 5B–E).

Figure 5

Differences in the tumor microenvironment indicators between the high and low IP groups in the WCH-TNBC cohort. (A) Volcano plot displaying the differently enriched HALLMARK pathways identified by gene set enrichment analysis in the high and low IP groups. (B) Representative images of CD31, CA9, and CDH2 IHC. (C) Differences in the CD31, CA9, and CDH2 IHC levels between the high and low IP groups. (D) Representative images of MHC-I, CD8, CD4, and CD20 IHC. (E) Differences in the MHC-I, CD8, CD4, and CD20 IHC levels between the high and low IP groups. IHC, immunohistochemistry; IP, ImPredict; TNBC, triple-negative breast cancer; WCH, West China Hospital.

We further investigated the associations of the IP score with CD8+ T cell spatial localization and two CD8+ T cell subtypes. According to the CD8 IHC stain, three main spatial phenotypes were defined (figure 6A): inflamed (N=16, tumor with CD8+ T cell infiltration in tumor center and invasive margin), excluded (N=17, tumor with CD8+ T cell embedded surrounding tumor border) and desert (N=19, tumor devoid of CD8+ T cell infiltration). We found that the inflamed phenotype presented a higher IP score than the desert and excluded phenotypes, and the IP-high group included more inflamed phenotypes (figure 6B). Similar results were also found in the public GSE177043 cohort (figure 6C). In addition to immunophenotypes, we evaluated the abundances of cytotoxic CD8+ T cell (GZMB+CD8+) and stem/progenitor-like CD8+ T cell (TCF1+CD8+) using mIHC imaging of 50 TNBC samples (figure 6D). Quantitative analysis demonstrated that the IP-high group exhibited significantly higher abundances of cytotoxic CD8+ T cells and stem/progenitor-like CD8+ T cells (figure 6E).

Figure 6

Associations of the IP score with immunophenotypes and CD8 subgroups in TNBC. (A) Representative images of CD8+ T cell spatial phenotypes in the tumor center and border. (B–C) Differences in the IP levels of three immunophenotypes and the portions of three immunophenotypes in the high and low IP groups of the WCH-TNBC cohort (B) and GSE177043 (C). (D) Representative mIHC images of TCF1, CD8, and PanCK. (E) Comparison of the abundances of TCF1+CD8+ T cells and GZMB+CD8+ T cells in the high and low IP groups. GSEA, gene set enrichment analysis; IP, ImPredict; TNBC, triple-negative breast cancer; WCH, West China Hospital.

Association of the IP score with TNBC clinical characteristics

Finally, we attempted to investigate whether the IP score is associated with clinical characteristics and reflects the clinical prognosis of TNBC. Notably, the IP-high group included a higher proportion of patients with basal subtype (figure 7A), and comprised more patients with the (basal-like1) BL1 and immunomodulatory (IM) subtypes, while the IP-low group consisted of more luminal androgen receptor (LAR) and mesenchymal-like (MES) subtypes (figure 7B, C). The IP-high group showed more immune-enriched (IE) phenotypes, while the IP-low group included more immune-enriched with fibrotic (IE/F), fibrotic (F), and immune-depleted (D) phenotypes (figure 7D). Besides, the IP-low group exhibited significantly poorer disease-specific survival (DSS) and overall survival (OS) compared with the IP-high group in the TCGA-TNBC cohort (figure 7E, F). More importantly, the IP score was also a prognostic indicator for DSS and OS, independent of histology type, pathological tumor stage, and pathological lymph node stage (N stage) (figure 7G, H). Considering the equal effects of the N stage and IP score for evaluating prognosis, we performed survival analysis of high and low IP groups in different N stages. Kaplan-Meier analysis indicated that the IP-low patients showed poorer DSS and OS compared with the IP-high patients in the N0 and N1-3 groups (online supplemental figure S20A–D). IP score seemed to perform better in stratifying patients into high-risk and low-risk in the N1-3 group than in the N0 group. Next, we developed two nomograms to predict DSS and OS: the DSS nomogram integrating N stage and IP score (online supplemental figure S20E), and the OS nomogram integrating histology, N stage, and IP score (online supplemental figure S20F). The C-indices of combined DSS and OS nomogram models were greater than 0.8. ROC analysis showed that the combined models were superior to the IP score for predicting the 3-year and 5-year DSS/OS, but inferior to the IP score for predicting the 8-year DSS/OS (online supplemental figure S20G–N). Therefore, it appears that the IP score can provide additional information for the prognosis evaluation of TNBC.

Figure 7

Exploratory analysis for the associations of IP score with TNBC clinical characteristics in the TCGA-TNBC cohort. (A–D) Proportions of different molecular subtypes in the high and low IP groups of TNBC. (E–F) Kaplan-Meier survival curves of DSS (E) and OS (F) in the high and low IP groups. (G–H) Multivariate analysis for the DSS (G) and OS (H) in the TCGA-TNBC cohort. DSS, disease-specific survival; IP, ImPredict; OS, overall survival; TCGA, The Cancer Genome Atlas; TNBC, triple-negative breast cancer.

Discussion

TNBC patients are prone to earlier relapse and have a poorer prognosis than other molecular subtypes. Previously, chemotherapy was the only systemic treatment for TNBC.26 Recent studies have confirmed that the PD-1 inhibitor pembrolizumab is effective in TNBC, changing the pattern of TNBC treatment.10 11 27 Emerging evidence demonstrates that ICI combined with chemotherapy can generate a significant increase in pCR rate for HR+/HER2− patients.6 7 However, not all patients can benefit from immunotherapy. Due to the potential toxicity, uncertain efficacy, and economic burden of immunotherapy, it is necessary to discriminate potential responders through effective biomarkers. Applications of existing biomarkers are controversial or insufficient in predictive power, thus developing predictive markers for immunotherapy in BC remains an urgent priority.9 In this study, we used the public I-SPY2 immunotherapy cohort and machine learning algorithms to construct a predictive gene score. It is promising in determining the NACI response of HER2− BC patients and screening patients who can benefit from immunotherapy. For TNBC, a high IP score indicates positive immunotherapy response, abundant immune infiltration, and an inflammatory phenotype, while a low IP score is related to immunosuppression and poor prognosis.

Unlike the advanced TNBC, the benefit of immunotherapy in early TNBC patients is irrelevant with PD-L1 level. In the KEYNOTE-522 study concerning the early TNBC, event-free survival shows an increase after the addition of pembrolizumab to NACT, regardless of PD-L1 expression status.28 KEYNOTE-355 study indicates that the progression-free survival benefits from pembrolizumab are only observed in PD-L1 positive (CPS>10) advanced/metastatic TNBC.27 Likewise, when another ICI atezolizumab (PD-L1 inhibitor) is used for neoadjuvant treatment of early TNBC, the increase in pCR rate is also independent of PD-L1 positivity.11 While in the first-line treatment of advanced/metastatic TNBC, atezolizumab combination only brings benefit for the OS of PD-L1 positive patients.29 In our study, the p for interaction of IP score in cohort A is the most significant, regardless of whether it is based on the median or the optimal cut-off grouping. In cohort B, possibly due to the small sample size of the chemotherapy group, the p for interaction of IP score does not reach statistical significance, but it also remains the smallest among all signatures. These results indicate that the IP score is the most effective marker for identifying patients suitable for immunotherapy addition. Therefore, patients determined to have a high IP score are recommended to add ICI drugs, which might significantly increase the pCR rate.

As for our IP model, we observe that many positive genes are involved in lymphocyte activation (TSPAN33, DAPP1, GZMB), response to IFNγ (GBP1, IL12RB2), leukocyte migration (MCOLN2, CCL13), viral process (LAMP3), and antigen presentation (HLA-E), whereas the negative genes are enriched in G-protein signaling (GPRC5C, RADIL, RGS22), microtubule activity (CCDC74B, KIF3A, FAM179B), dehydrogenase activity (NQO1, DHRS2), extracellular matrix (UGDH), hypoxia (CA12), and angiogenesis (VEZF1), which are considered immunosuppressive.30–33 We construct the model using the PCA, which is a dimension reduction method. In brief, by employing PCA to extract PC1, we can obtain the main expression level of the gene set, thus eliminating the influence of some outlier genes on the whole.17 Unlike ssGSEA, PCA does not rely on the entire transcriptome data, making it applicable for small gene sets with good correlation. Considering that there are some challenges in terms of cost-effectiveness and technical aspects when obtaining the whole transcriptome, PCA is an ideal method for this IP model to be translated into a clinical-grade biomarker.

Our findings suggest that IP score is positively associated with immune infiltration and cytotoxic T-cell activity, but not with TMB. In line with our data, McGrail et al have demonstrated that in some type II cancers including BC and prostate cancer, CD8+ T cells do not increase with more new tumor antigens, and the former is usually the basis for effective immunotherapy response.34 These findings reflect the complexity of the immune system. Numerous mutations do not necessarily convert into immunogenic new antigens. The quality of mutations, the production, and identification of new antigens, along with the intricate immune microenvironment, are also pivotal factors in eliciting an immune response.35

GSEA reveals that immunotherapy response is positively related to the IFN signaling and immune activation, but is also negatively affected by some pathways such as steroid hormone biosynthesis, fatty acid metabolism, epithelial-mesenchymal transition, angiogenesis, etc. Interestingly, we also discover that the IM and IE subtypes, which have been confirmed to be more sensitive to ICIs,36 37 dominate in the IP-high group, while the LAR and MES types are more prevalent in the IP-low group. Previous studies have shown that the LAR TNBC exhibits enhanced hormone receptor signaling and fatty acid metabolism, and MES TNBC displays active EMT and angiogenesis activity,4 38 which are consistent with the aforementioned findings. In summary, our findings provide new insights into the mechanism link between the immunotherapy response and biological process, emphasizing the non-negligible role of some regulatory mechanisms in immunotherapy resistance. This may potentially facilitate the precision of immunotherapy and the development of drug combinations.

Some limitations in our study should be noted. First, the sample sizes of immunotherapy cohorts in our study are relatively small, and it needs future validation in a larger, prospective cohort to determine the cut-off of the IP score. Second, the heterogeneity of the baseline characteristics may affect the predictive performance. We are trying to collect more samples and develop an integrated model incorporating clinicopathological factors to improve accuracy.

In conclusion, we established a robust model for predicting the NACI response of HER2− BC and selecting patients suitable for adding immunotherapy to NACT. This model is a promising tool to facilitate the optimization of immunotherapeutic regimens for HER2− BC patients.

Supplemental material

Data availability statement

Data are available in a public, open access repository. Data are available upon reasonable request. All data relevant to the study are included in the article or uploaded as supplementary information.

Data availability statement

The raw sequencing data of the WCH-TNBC cohort in this study has been deposited to NCBI (Project ID: PRJNA1081147, SRA ID: SRP492770). Data generated by the TCGA and Gene Expression Omnibus (GEO) are available publicly from the https://portal.gdc.cancer.gov/ and https://www.ncbi.nlm.nih.gov/geo/ websites. Other data supporting the findings of this study are available within the article, its supplemental files, and from the corresponding author upon reasonable request.

Ethics statements

Patient consent for publication

Ethics approval

This study was approved by the Ethics Committee of West China Hospital, Sichuan University, China (No.20221410), and abided by the Declaration of Helsinki of using tissue samples for scientific research only. The written informed consent was exempted by the ethical committee for this retrospective study.

Acknowledgments

We would like to thank and express our heartfelt gratitude to Dr. Zhong Xiaorong (Breast Disease Center, West China Hospital) and Li Lingling (Data Manage Center, West China Hospital), who assisted with the data collection of the validation cohort in West China Hospital.

References

Supplementary materials

Footnotes

  • XL and ZG contributed equally.

  • Contributors XL and ZG contributed to the study design, data analysis, and paper writing. ZG and HC contributed to sample collection and data acquisition. LL, FC, and CB performed the experiments. HB contributed to project supervision and paper review. All authors read and approved the final version of the manuscript. HB takes full responsibility for the finished work and the conduct of the study, had access to the data, and controlled the decision to publish.

  • Funding This work was supported by the 1·3·5 Project for disciplines of excellence, West China Hospital, Sichuan University (ZYGD18012); Natural Science Foundation of Sichuan Province (2022NSFSC1385); Post-Doctor Research Project, West China Hospital, Sichuan University (2021HXBH080).

  • Competing interests None declared.

  • Provenance and peer review Not commissioned; externally peer reviewed.

  • Supplemental material This content has been supplied by the author(s). It has not been vetted by BMJ Publishing Group Limited (BMJ) and may not have been peer-reviewed. Any opinions or recommendations discussed are solely those of the author(s) and are not endorsed by BMJ. BMJ disclaims all liability and responsibility arising from any reliance placed on the content. Where the content includes any translated material, BMJ does not warrant the accuracy and reliability of the translations (including but not limited to local regulations, clinical guidelines, terminology, drug names and drug dosages), and is not responsible for any error and/or omissions arising from translation and adaptation or otherwise.