Identification and characterization of chemoimmunomodulatory induction trajectories in breast cancer via unsupervised transcriptomic analyses of pre- and post-treatment specimens

Chemotherapy is key to curative treatment strategies for aggressive breast cancer (BC). Despite this, the molecular features governing the immunomodulatory effects of chemotherapy (chemoimmunomodulation; CIM), which enable efficacy and synergy with other treatment modalities, are understudied, thereby limiting optimization of chemotherapeutic regimens to improve patient outcomes. A key roadblock to these efforts has been the lack of a robust framework for classifying CIM induction trajectories. To address this challenge, we developed the CIM Induction Classifier (CIMIC), an iterative unsupervised clustering pipeline that uses delta gene expression (Δlog2(TPM+1)) of 3,100 genes across 19 CIM pathways from paired pre- and post-treatment tumor transcriptomes to classify CIM induction trajectories. We applied CIMIC to data from two cohorts (GSE191127 (N = 20) and GSE12845 (N = 16)) of BC patients treated with neoadjuvant chemotherapy to identify distinct CIM trajectories and characterize the tumor-intrinsic and tumor-extrinsic programs underlying CIM heterogeneity. In both cohorts, the classifier consistently assigned patients to one of two CIM trajectories, functional CIM (Fun-CIM; N = 17) or dysfunctional CIM (Dys-CIM; N = 19). Overrepresentation analyses of induced genes showed an enrichment for immunostimulatory programs, including immune-effector differentiation, function, and cell killing (FDR ≤ 0.01) and tumor-intrinsic stress adaptation programs, including proteostasis, unfolded protein response, and mitochondrial homeostasis (FDR ≤ 0.01) for the Fun-CIM and Dys-CIM, respectively. Interestingly, Dys-CIM was overrepresented in the basal-like molecular subtype in both cohorts (p < 0.05) and, where data were available, was associated with recurrence (p = 0.041; N = 20). To further evaluate the clinical relevance of these induction states, we derived trajectory-specific signatures and applied single-sample gene set analysis to baseline tumors from chemotherapy-treated patients in the METABRIC (N = 412) and SCAN-B (N = 2,774) cohorts. In multivariable Cox models adjusted for treatment, clinicopathological features, and molecular subtype, the high Dys-CIM signature was associated with inferior overall survival (METABRIC, HR = 1.83, p = 0.041; SCANB, HR = 1.73, p = 0.014), poorer disease-specific survival (METABRIC, HR = 1.94, p = 0.035), and lower recurrence-free survival (METABRIC, HR = 1.68, p = 0.019), suggesting that a predisposition toward Dys-CIM represents a clinically adverse signature detectable at baseline. Taken together, our findings herein provide rationale for dissecting clinically relevant CIM trajectories and responses, supporting future efforts to link induced CIM trajectories to clinical outcomes and personalized immunomodulation strategies.