Cancer-immunity cycle-based molecular subtypes in breast cancer predict the response to immune checkpoint inhibitors

Patient cohort

Patients diagnosed with malignant breast cancer who were undergoing treatment in the Department of Breast Surgery at Fudan University Shanghai Cancer Center (FUSCC) (Shanghai, China) between 1 January 2013 and 31 December 2014 and willing to participate in this study were retrospectively selected. Detailed information on biospecimen collection, genomic data, transcriptomic profiles, and metabolomic detection were described in a previous study²⁵. In the current study samples with qualified transcriptome data were selected. A total of 752 patients were enrolled in this study. Missing values across different omics layers were processed following data-type-specific strategies. Missing metabolomic data values were filled by one-half of the minimum value. Mutation data were treated as wild-type when absent and copy number alteration (CNA) values were only retained for samples with complete genomic segmentation.

All procedures were performed in accordance with the Helsinki Agreement and with the approval of the independent Ethics Committee/Institutional Review Board of the FUSCC Ethical Committee. Ethical approval and informed consent were obtained [Approval No. 2011226-17 (2010-ZZK-31)].

Building of the CIC step enrichment matrix and the CIC score

An unbiased selection of genes relating to every step in the CIC was conducted based on the framework of published molecular signatures database²⁶. Several changes were made to better characterize the status of anticancer immunity. First, step 3 (priming and activation) was excluded in the analysis because step 3 is primarily defined as the activation of peripheral T cells and cannot be evaluated using in situ tumor indicators. Next, controversial genes used to evaluate a specific step were removed and manually collected genes involved in the CIC steps from recent published studies were included to update the evaluation gene set. Finally, 126 genes were identified and compiled to curate the CIC gene set collection (Table S1).

The activity levels of every step were calculated using single-sample Gene Set Enrichment Analysis (ssGSEA) based on gene expression profiles, as previous described²⁶. In brief, positive and negative scores were firstly calculated for each step. The final activation index was then obtained by subtracting the negative score from the positive score. A linear normalization for the score of each step was mapped to a 0–1 range to ensure that each step contributed equally to the final score. Ultimately, the CIC score was defined as the sum of the normalized scores from each step.

CIC clustering analysis

A correlation matrix was generated using Pearson coefficients based on the enrichment scores of every step for each patient to determine specific CIC step deficiencies. The pheatmap package (v1.0.12) in R was then used to create a hierarchically clustered heatmap to visually identify distinct patterns of enrichment and enable the unsupervised clustering of patients based on the similarity of the immune profiles. The gap statistic, elbow method, and Calinski–Harabasz Index were used to determine the optimal number of stable CIC subtypes. This analysis successfully identified and characterized three unique CIC subsets; each CIC subset was defined by a specific distribution of enrichment across the CIC steps.

Survival analysis

Survival analysis was performed using R packages [Survival (v3.4.0) and Survminer (v0.4.9)]. Kaplan–Meier survival curves were constructed and compared to the log-rank test to assess survival differences among the CIC subgroups. Cox proportional hazards models were used to identify independent prognostic variables. Initially, a univariate analysis was performed on variables (the CIC score, lymph node status, tumor size, histological grade, and Ki67 percentage). Subsequently, all variables with significant associations in the univariate analysis were included in a multivariate Cox regression model to determine the independent prognostic value.

Immune cell infiltration analysis

CIBERSORT was used to deconvolute immune cell fractions from bulk RNA-sequencing data²⁷. Specifically, the LM22 leukocyte signature matrix, which contains the expression of 547 genes across 22 immune cell subtypes, was utilized in the R environment with 1000 permutations to generate detailed leukocyte fractions for each patient. The ESTIMATE (v1.0.13)²⁸ and xCell (v1.1.0)²⁹ packages were used to assess the TME status by estimating levels of stromal cells, immune infiltration, and tumor purity from the transcriptome data.

Mutational signature and somatic CNA analysis

The deconstructSigs package (v1.9.0) was used to identify mutational signatures in samples with single nucleotide variants (SNVs)³⁰. Raw mutation counts were first normalized by callable genome size and sequencing depth to reduce technical confounding. Samples with insufficient tumor purity based on ABSOLUTE estimates were excluded prior to signature deconvolution. Patient mutational profiles were deconvoluted using a weighted combination of breast cancer-related signatures. In addition, differences in CNA frequency were evaluated with Pearson’s chi-squared test and Fisher’s exact test. P-values for these comparisons were corrected for multiple testing using the false discovery rate (FDR) method.

Differential abundance (DA) score

Metabolomic data were pre-processed following established practices from the extant literature and publicly available protocols, including the use of QC samples for signal correction, LOESS-based adjustment to remove within- and between-batch variation, and normalization procedures applied prior to pathway-level analyses. The DA score quantifies the propensity of a pathway to exhibit elevated levels of metabolites compared to other CIC groups. First, a non-parametric test for DA (the Benjamini–Hochberg corrected Mann–Whitney Wilcoxon test) was applied to all metabolites in each pathway to identify significant increases or decreases in abundance. The DA score was calculated with the following formula:

The DA scores range from −1 to 1. A score of −1 indicated that all metabolites in a pathway were decreased, whereas a score of 1 indicated that all metabolites were increased.

Cell lines

Cell lines (MDA-MB-231 and HCC1806) were obtained from the American Type Culture Collection [ATCC] (Manassas, VA, USA). MDA-MB-231 was cultured in DMEM (Basal Media, Shanghai, China) supplemented with 10% fetal bovine serum [FBS] (Gibco,Waltham, MA, USA), while HCC1806 was cultured in RPMI-1640 (Basal Media, Shanghai, China) with 10% FBS. Regular testing using a mycoplasma detection kit (Vazyme, Nanjing, Jiangsu, China) verified a lack of mycoplasma contamination in all cell lines. Standard culture conditions, including incubation at 37°C in 5% CO₂, were maintained for all cells. Furthermore, the authenticity of the cell lines was verified via short tandem repeat (STR) profiling.

Small interfering (si)RNA construction and transfection

PSAT1 siRNAs were synthesized by GenePharma Biotech (Shanghai, China). The sequences are listed in Table S7. MDA-MB-231 and HCC1806 cells were grown to 60% confluency, then transfected with control siRNA and PSAT1 siRNA utilizing Lipofectamine 2000 (Invitrogen, Carlsbad, CA, USA). The evaluation of knockdown efficiency was conducted 24–48 h post-transfection.

Plasmids and viral transduction

PSAT1 overexpression and control plasmids were purchased from YiXueSheng Biosciences, Inc. (Shanghai, China). HEK293T cells were transfected with each lentivirus expression vector and packaging plasmid mix using Lipofectamine 2000 (Invitrogen, Carlsbad, CA, USA) to generate stable expressing cell lines. The supernatant containing virus was collected 48 h after transfection, filtered, and used for infecting target cells in the presence of 8 μg/mL of polybrene (#H9268; Sigma, St. Louis, MO, USA).

RNA isolation and qRT-PCR

Complementary DNA (cDNA) was synthesized from the extracted RNA using the HiScript III All-in-one RT SuperMix (Vazyme, Nanjing, Jiangsu, China) after cell isolation with TRIzol reagent (Invitrogen). Primers were designed via the PrimerBank website (https://pga.mgh.harvard.edu/primerbank). ChamQ Universal SYBR qPCR Master Mix (Vazyme, Nanjing, Jiangsu, China) was used for qRT-PCR analysis on a QuantStudio™ 6 Flex System (Foster City, CA, USA). The qRT-PCR results were standardized to the level of actin beta (ACTB) gene expression, which served as an internal control. Detailed information on the gene-specific primers utilized for amplification is provided in Table S7.

Western blotting

Protein concentrations were determined using a BCA Assay Kit (Solarbio, Beijing, China) following extraction with NP-40 lysis buffer. Equal amounts of protein were resolved by SDS-PAGE and transferred onto polyvinylidene difluoride (PVDF) membranes (Millipore, Billerica, MA, USA). The PVDF membranes were incubated with primary antibodies against PSAT1 (1:1,000, #A14124; ABclonal, Wuhan, Hubei, China), PD-L1 (1:1,000, #A19135; ABclonal), and vinculin (1:5,000, #A2752; ABclonal), followed by HRP-conjugated Goat Anti-Rabbit IgG(H+L) (Proteintech, Wuhan, Hubei, China).Protein signals were detected with LumiBest ECL Substrate Solution Kit (ShareBio, Shanghai, China) and imaged with the ChemiDoc XRS System (Bio-Rad, Hercules, CA, USA).

Statistical analysis

The Mann–Whitney Wilcoxon and Kruskal–Wallis tests were used for continuous variables. Pearson’s chi-squared test and Fisher’s exact test were used for categorical variables. Pearson’s or Spearman’s correlation was used to create matrices. P-values were adjusted using the FDR correction for multiple comparisons. All tests were two-sided unless otherwise noted. All statistical analyses were performed using R (v4.2.1).

CIC features correlate with clinical characteristics and indicate prognosis in breast cancer patients

A CIC signature that included 126 genes representing 6 major steps in the CIC (Materials and Methods; Table S1) was first established to systematically characterize the anti-tumor immunity phenotypes in breast cancer patients. We then estimated the activation score of each step for every patient based on the gene expression profile described in our previous study²⁵. As a result, the normalized scores from each step were added to form the CIC score (Figure 1A and Table S2). The CIC score framework was performed in three external breast cancer transcriptomic cohorts the cancer genome atlas (TCGA), molecular taxonomy of breast cancer international consortium (METABRIC), and GSE103091 to validate the findings.

• Download figure
• Open in new tab
• Download powerpoint

• Download figure
• Open in new tab
• Download powerpoint

Figure 1

Correlations of CIC score with clinicopathologic features and prognosis in the breast cancer patients. (A) Genes for each step of the cancer-immunity cycle were selected from published signatures and recent studies. Step 3 (priming and activation) and controversial genes were excluded. A total of 126 genes were curated to form the final CIC gene set. The panel on the right displays the names of each step and the number of genes included. (B) Associations of CIC scores with clinicopathologic and molecular features in breast cancer patients. The stacked bar chart shows CIC scores of samples, ordered from left to right by increasing with different colors representing the various CIC steps. The heatmap represents the activation status of CIC steps for each patient. Bottom rows represent clinicopathologic and molecular features of breast cancer patients. (C–E) Comparisons of CIC scores with subgroups of clinical subtypes (C), PAM50 subtypes, (D) and histological subtypes (E). (F) Kaplan–Meier survival curve according to CIC score levels for RFS in the breast cancer cohort. Patients were stratified into CIC-high and -low groups using the cohort median CIC score as the cut-off value. (G–I) Comparisons of CIC scores with subgroups of survival status in TCGA-BRCA, METABRIC, and GSE103091 cohorts. (J, K) Univariable (J) and Multivariable (K) Cox regression analyses of CIC scores with clinical and molecular characteristics. (L) Comparisons of CIC scores with pCR status subgroups in the I-SPY2 pembrolizumab treatment arm. P values are displayed with two decimal places for formatting consistency. The exact P value for the comparison is P = 0.0495 (wilcoxon test). Boxplots of (C)–(E) show first and third quartiles, whiskers extend to the lowest and highest value within the 1.5× interquartile range, and the line indicates the median. P values were from the Kruskal–Wallis test. IDC, invasive ductal carcinoma; RFS, recurrence-free survival; OS, overall survival; pCR, pathologic complete response; Pembro, pembrolizumab.

The relationship between the CIC score and clinicopathologic characteristics in our breast cancer cohort was determined first. The CIC score was significantly associated with clinical subtypes, the PAM50 subtype, and the histological type (Figure 1B–E). Specifically, patients with HR-positive breast cancer had lower CIC scores, while patients with the basal subtype and invasive ductal carcinoma (IDC) had higher CIC scores. In addition, both the histological grade and Ki67 index, recognized markers of tumor aggressiveness, had a positive correlation with the CIC score, whereas lymph node status had no association (Figure S1A–C). This observation underscores the intricate interplay between anti-tumor immune activity and the malignant biological behavior of tumors. Interestingly, while age had no association with the CIC score in the overall breast cancer cohort, a significant negative correlation was detected in the TNBC subgroup, suggesting the presence of subgroup-specific immunologic features (Figure S1D).

The prognostic predictive value of the CIC score was then examined. Kaplan–Meier analysis revealed significant associations between higher CIC scores and better clinical outcomes [log-rank, recurrence-free survival (RFS), P = 0.045; Figure 1F]. Additional three breast cancer transcriptomic datasets confirmed these findings (Figures 1G–I and S1E). The CIC score was confirmed to be a protective factor for RFS based on Univariable Cox analysis (Figure 1J). Multivariable Cox regression analysis validated that the CIC score is an independent prognostic factor (Figure 1K). The predictive value of the CIC score for therapeutic efficacy in the I-SPY2 neoadjuvant therapy cohort was investigated. The CIC score was elevated in patients who achieved a pCR following neoadjuvant immunotherapy or chemotherapy, suggesting that anti-tumor immunity may modulate the efficacy of multiple therapeutic approaches (Figures 1L and S1F). Overall, the CIC score was associated with various clinical parameters and showed predictive value for both prognosis and treatment response, supporting the potential utility in patient stratification.

Identification of three subtypes by CIC deficiency features

A patient correlation matrix derived from the enrichment patterns of CIC step signatures was constructed across samples to delineate the subtypes of CIC deficiency based on the heterogeneous distribution of CIC signatures. Hierarchical clustering of the correlation matrix, representing the similarity of each patient to others in the cohort, was performed using Euclidean distance. As a result, distinct subtypes of the CIC based on distribution of the CIC steps scores were uncovered and termed cluster 1 [C1], 2 [C2], and 3 [C3] (Figure 2A and Table S3).

• Download figure
• Open in new tab
• Download powerpoint

Figure 2

Identification of three CIC subtypes. (A) Heatmap of Pearson correlation of patients with breast cancer across six CIC steps (up). Distance-based clustering revealed three distinct clusters across patients. Heatmap below showed the enrichment score of every step. (B) Box plots depicting the enrichment level of six CIC steps in the three CIC clusters. (C) Radar charts depicting enrichment or deficiency in CIC clusters at each CIC step. (D) Levels and trends of CIC scores among three CIC subtypes. (E) Comparisons of clinical subtypes fraction with CIC clusters. (F) Comparisons of Ki67 index with CIC clusters. (G) Correlations of T stage with CIC subtypes. Boxplots of (B) and (D) show first and third quartiles, whiskers extend to the lowest and highest value within the 1.5× interquartile range, and the line indicates the median. P values were from Kruskal–Wallis test. For (E) and (G), P values were obtained from two-sided Fisher’s exact test and adjusted by the Benjamini–Hochberg procedure.

The distinct CIC features associated with each subtype were further identified based on the classification. Analysis of the functional CIC status revealed downregulation in subtype C1 and upregulation in subtype C3 across most steps (Figure 2A, B). The overall CIC scores also rose progressively from C1–C3, which indicated a shift in CIC characteristics from inactive to active (Figure 2C, D). The stepwise increase in CIC enrichment scores from C1–C3 was most evident in steps 4 and 5, related to immune cell migration and infiltration. In contrast, step 7, representing T-cell cytotoxic function, was not significantly downregulated in the immune-cold subtype C1 and even slightly increased. These findings suggested that the absence of immune cell infiltration is a key contributor to CIC deficiency in C1 patients. C2 displays intermediate features positioned between C1 and C3, while accompanied by a notable deficiency in step 2 (Figure 2A, B). The defect in step 2, which is responsible for presenting tumor-associated antigens to APCs, points to a potentially distinct immunosuppressive mechanism that may underlie impaired anti-tumor immunity in C2 patients. The stability of CIC classification in the METABRIC breast cancer cohort was validated (Figure S2A). The clustering results in the METABRIC cohort confirmed the conservation of CIC subtypes and enrichment analysis across CIC steps reproduced the CIC features of C1, C2, and C3 noted in the discovery cohort (Figure S2B–D).

The associations between CIC subtypes and clinical characteristics of breast cancer patients were further analyzed. The CIC clusters also exhibited heterogeneity in clinical subtypes, like the CIC scores. Specifically, C1 was enriched with HR-positive patients, C2 with TNBC patients, and C3 with HER2-positive patients (Figure 2E). The C2 subtype also exhibited distinct characteristics representing associations with prognostic risk features that were primarily reflected by a higher Ki67 index and advanced tumor T stage (Figure 2F, G).

Together, this analysis provides an interpretable framework to classify CIC status in breast cancer patients. Specifically, CIC-based clustering identified C1 and C3 subtypes with classical “cold” and “hot” tumor features and a distinct C2 subtype with antigen presentation defects.

Immune characteristics of CIC subtypes

The immunologic features in the subgroup of breast cancer samples were then examined across CIC subtypes. First, the major functional pathways of T cells were analyzed (Figure 3A). Pathways related to T cell activation, differentiation, migration, and cytotoxicity were all upregulated in C3 patients. The increased expression of co-stimulatory molecules and the concurrent upregulation of exhaustion markers in C3 further confirmed the activation status of T cells (Figure S3A). Enrichment analysis of antigen presentation pathways showed that C2 exhibited downregulation like the immune-cold C1 subtype, underscoring the distinct CIC functional deficiency (Figure 3A). Next, hallmark pathway distributions were examined across subtypes to investigate potential drivers of CIC states (Figure 3B). C3 patients exhibited upregulation of interferon response and inflammatory pathways, indicating active anti-tumor immunity, which was consistent with the T cell functional analysis. Estrogen response pathways were upregulated in C1, likely reflecting the higher prevalence of HR-positive tumors in this subtype. Activation of DNA damage repair and mitosis-related pathways was noted in C2 patients. Together with the previously reported increased Ki67 index, C2 demonstrated a highly active proliferative phenotype (Figures 2F and 3B).

• Download figure
• Open in new tab
• Download powerpoint

• Download figure
• Open in new tab
• Download powerpoint

Figure 3

Immune characteristics of CIC subtypes. (A) ssGSEA enrichment analysis of T cell functional GO pathways across different CIC subtypes. (B) ssGSEA enrichment analysis of hallmark pathways across different CIC subtypes. (C) Comparisons of TILs level among CIC clusters. (D, E) Comparisons of immune score calculated by xCell (D) and ESTIMATE (E) among CIC clusters. (F) Comparisons of tumor purity calculated by ESTIMATE among CIC clusters. (G, H) Boxplots showing CIBERSORT deconvolution results for macrophages (G) and T cells (H) across CIC clusters. (I) Proportions of main cell types based on breast cancer scRNA-seq dataset among CIC subtypes. (J) Violin plots depicting Tpex and Tex signature enrichment among CD8+ T cells from three CIC subtypes. (K) Dot plots depicting selected genes expression among DCs from three CIC subtypes. Boxplots of (C)–(H) show first and third quartiles, whiskers extend to the lowest and highest value within the 1.5× interquartile range, and the line indicates the median. P values were from wilcoxon test adjusted by the Benjamini–Hochberg procedure.

Tumor-infiltrating lymphocytes (TILs) were evaluated in pathologic sections to assess immune cell enrichment across CIC subtypes. The increased TILs levels from C1-C3 highlighted enhanced immune infiltration across subtypes, mirroring the trends observed in steps 4 and 5 (Figure 3C). Algorithms were used to further validate the trends based on transcriptomics to assess immune infiltration in breast cancer patients. Results from the xCell and ESTIMATE suggested the enrichment of immune cells in C3, characterized by increased immune infiltration scores and reduced tumor purity (Figure 3D–F). In addition, these algorithms assessed stromal cell infiltration across subgroups, revealing a significant reduction in stromal score in C2 (Figure S3B, C).

CIBERSORT was applied to deconvolve the breast cancer transcriptome samples to further investigate the differences in infiltrating immune cell types between the clusters (Table S4). The observed variations were primarily concentrated in macrophages and T cells (Figure 3G, H). Patients in C1 had an enrichment of M2 macrophages and a downregulation of M1 macrophages (Figure 3G). Macrophages in C1 patients exhibited a more pro-tumor phenotype based on previous research involving M1 and M2 macrophages. T cells displayed a pattern like overall immune infiltration with CD8⁺ and activated CD4⁺ memory T cells progressively increasing from C1-C3 (Figure 3H). Notably, regulatory T cells (Tregs), which have a key role in immunosuppression, were enriched in C2 patients (Figure 3F). xCell and ESTIMATE were also used to assess the proportion of stromal cells; endothelial cells and fibroblasts were downregulated in C2, suggesting a looser extracellular matrix in this subtype (Figure S3D, E).

Single-cell RNA sequencing (scRNA-seq) offers a more accurate method to characterize the infiltration levels of various cellular subsets in the TME. Therefore, published paired bulk transcriptome and scRNA-seq data of breast cancer patients were used to validate differences in immune cell infiltration among clusters³¹. The patients were divided into three groups according to the CIC levels and the changes in cellular proportions among different clusters were further analyzed at the single-cell level based on the bulk transcriptome data (Figure S3F, G). The analysis of the single-cell data validated the enrichment of T cells in patients within C3 (Figure 3I). Tpex markers, which are derived from the functional enrichment analysis of tumor-infiltrating CD8+ T cells at the single-cell level, were associated with the immunotherapy response and shown to be upregulated in patients from C3 (Figure 3J). In contrast, Tex markers, which indicate T cell exhaustion, had the highest expression in C2 (Figure 3J). In addition, an unusually high proportion of dendritic cells (DCs) was present in C2 (Figure 3I). Therefore, a functional analysis of DCs was performed and the DCs in C2 patients showed low expression of activation and functional molecules, suggesting functional impairment of DCs (Figure 3K).

These results collectively suggested that C3 exhibits high infiltration of immune cells and C1 could be considered a classic immune-cold tumor with immunosuppressed macrophages. The high infiltration of Tregs and dysfunction of DCs led to CIC deficiency despite the moderate level of immune cell infiltration and loose matrix in C2 patients.

The CIC subtypes serve as a pan-cancer prognostic biomarker for ICI therapy

The prognostic value of CIC subtypes was first evaluated in our breast cancer cohort to determine the clinical relevance of CIC clustering. C1 represented the poorest prognosis, C3 the best prognosis, and C2 an intermediate prognosis based on Kaplan–Meier analysis [log-rank, distant metastasis-free survival (DMFS) P = 0.022; Figure 4A, B]. These findings were confirmed by the breast cancer transcriptomic dataset from METABRIC (Figure 4C).

• Download figure
• Open in new tab
• Download powerpoint

Figure 4

CIC subtypes are prognostic of response to ICI. (A, B) Kaplan–Meier survival curve (A) and Cox regression analysis (B) according to CIC clusters for DMFS in the breast cancer cohort. (C) Cox regression analysis according to CIC clusters for DMFS in METABRIC cohort. (D–F) Heatmap of Pearson correlation of patients receiving ICI therapy with resectable breast cancer (D; n = 69), metastatic melanoma (E; n = 121) and metastatic NSCLC (F; n = 27) across six steps of CIC. CIC subtypes are preserved in these cohorts. All analyses were performed using pre-treatment RNA-seq samples. (G, H) Box plots depicting the enrichment level of CIC scores and step 2 and 4 signatures in the three CIC clusters in the resectable breast cancer (G) and metastatic melanoma (H) cohorts. (I) pCR rate for patients in different CIC clusters upon treatment with ICI in resectable breast cancer patients. (J, K) Response rate for patients in different CIC clusters upon treatment with ICI in metastatic melanoma patients (J). Kaplan–Meier survival curve according to CIC clusters for OS in this cohort (K). (L) Benefit rate for patients in different CIC clusters upon treatment with ICI in metastatic NSCLC patients. Boxplots of (G) and (H) show first and third quartiles, respectively, whiskers extend to the lowest and highest value within the 1.5× interquartile range, and the line indicates the median. P values were from Kruskal–Wallis test. For (I), (J), and (L), P values were obtained from two-sided Fisher’s exact test and adjusted by the Benjamini–Hochberg procedure. DMFS, distant metastasis-free survival; NSCLC, non-small cell lung cancer; DCB, durable clinical benefit; NCB, non-clinical benefit.

• Download figure
• Open in new tab
• Download powerpoint

Figure 5

Genomic alterations among CIC subtypes. (A) Comparisons of tumor mutation burden (TMB) among CIC clusters. (B) Comparisons of ASCAT ploidy among CIC clusters. (C) Heatmap showing the comparison of immune receptor repertoires and tumor antigen production across CIC subtypes. (D) Bar plots showing the distribution of HLA-I LOH status among CIC subtypes. (E) Comparison of the somatic CNAs between every two clusters among the three CIC subtypes in breast cancer patients. Plots illustrates the −log₁₀ P value of each gene when compared between every two subtypes among the three subtypes in the amplification/gain-centric (red) or deletion/loss-centric (blue) calculations. (F) Somatic mutation profile among CIC clusters in breast cancer patients. Top genes with an overall frequency > 3% are shown. (G) Heatmap showing the contribution of breast cancer–related mutational signatures among CIC clusters. Boxplots of (A) and (B) show first and third quartiles, respectively, whiskers extend to the lowest and highest value within the 1.5× interquartile range, and the line indicates the median. P values were from the wilcoxon test were adjusted by the Benjamini–Hochberg procedure. For (D), P values were obtained from two-sided Fisher’s exact test and adjusted by the Benjamini–Hochberg procedure.

A significant challenge in optimizing ICI therapy lies in predicting treatment efficacy and stratifying patients, which remains a key area of research given the variable outcomes. Previous analyses have suggested that CIC status not only affects patient prognosis but is also closely related to the efficacy of ICI therapy (Figure 1F–K). Pre-treatment RNA-seq data from pan-cancer cohorts of ICI-treated patients were analyzed to assess the predictive value of CIC subtypes for ICI therapy. Six cohorts of patients were assembled, consisting of five different tumor types challenged with ICIs (Table S5). Cohorts of patients with early-stage breast cancer and early-stage non-small cell lung cancer (NSCLC) were assembled for neoadjuvant anti-PD1 therapy. Patients with metastatic melanoma, NSCLC, glioblastoma, and gastric cancer were also included.

The CIC subtypes were consistently maintained across all cohorts treated with anti-PD1 therapy, including patients with early-stage breast cancer³² (Figure 4D), early-stage NSCLC³³ (Figure S4A), metastatic melanoma³⁴ (Figure 4E), metastatic NSCLC³⁵ (Figure 4F), metastatic glioblastoma³⁶ (Figure S4B), and metastatic gastric cancer³⁷ (Figure S4C). The distribution of enrichment scores was largely consistent with pan-cancer analysis, showing that CIC scores progressively increased from C1-C3 with an upregulation of step 4 activation and a deficiency in step 2, corresponding to C2 (Figures 4G, H and S4D–G). C3 displayed a marked increase in response rate to anti-PD1 therapy in resectable and metastatic tumors. The improved response rate observed in C3 tumors was characterized by a significantly increased pCR in early-stage breast cancer (Figure 4I) and NSCLC patients (Figure S4J) receiving adjuvant ICI therapy. Similarly, C3 subsets showed the best treatment response and received the most benefit from ICI therapy among patients with metastatic cancer, while C1 patients had the worst response and the least benefit (melanoma Figure 4J, K; NSCLC, Figure 4L; glioblastoma, Figure S4K; gastric cancer, Figure S4L). In addition, CIC clustering was shown to have prognostic value for the survival of patients undergoing ICI therapy in the metastatic melanoma (Figures 4K and S4H) and lung cancer cohorts (Figure S4I) with survival data. Patients in C1 had a significantly shorter progression-free survival (PFS) and overall survival (OS) in the metastatic melanoma cohort. Likewise, a trend was observed in the metastatic NSCLC cohort, although not statistically significant. The aforementioned results demonstrated that CIC subtypes can be prognostic indicators of response and survival in patients treated with ICIs.

Genomic alterations underlie defects in the CIC

The key oncogenic mutations and CNAs were analyzed to identify factors driving the distinctions among CIC clusters in breast cancer patients. Analysis of global genomic alterations revealed that C2 harbored the highest tumor mutational burden (TMB) and aneuploidy, whereas C1 exhibited relatively low levels (Figure 5A, B). The immunogenomic landscape of the CIC clusters was further characterized. The immune repertoire of tumor-infiltrating immune cells was predictive of response to ICI therapy based on prior research. C2 showed a higher level of tumor neoantigen production, which may be associated with the high level of genomic alterations (Figure 5C). However, patients in C2 failed to mount stronger tumor antigen responses despite the higher TMB. Therefore, the status of HLA across different CIC subgroups was investigated. HLA-I loss of heterozygosity (LOH) and allele imbalance (AI) samples were conspicuously enriched within C2 patients, suggesting the reason why these patients presented antigen presentation defects (Figures 5D and S5A). The clonal distribution of T cell receptors (TCRs) and B cell receptors (BCRs) among CIC clusters was also elucidated (Figure 5C and Table S6). Patients in C1 had lower diversity and higher clonality of TCRs, indicating the lack of response to ICI therapy. However, BCR diversity and clonality followed a distinct pattern. Specifically, diversity is an indicator of ICI responsiveness and was markedly reduced in C2 patients. The CNA changes were then analyzed among CIC clusters as an exploratory analysis. CNA analysis demonstrated that the amplification of DNA fragment 10p14 was more frequent in C1 and the loss of 23p11 was more frequent in C2 (Figure 5E). 23p11 encodes a variety of tumor-associated antigens, such as the cancer/testis antigen G antigen (GAGE) and cancer testis antigen (CTA). Deletion of this region may lead to a reduction in the release of tumor-associated antigens, resulting in a deficiency in CIC function.

The differential prevalence of non-synonymous mutations with an overall frequency > 3% in the breast cancer cohort was investigated to explore genomic determinants of CIC features (Figure 5F). Among top mutations, GATA3 and NF1 showed statistically significant differences among CIC clusters (Figure S5B, C). Specifically, GATA3 mutations occurred more frequently in patients in C1 (Figure S5B), while patients in C2 had a higher proportion of NF1 mutations (Figure S5C). Previous studies have reported that GATA3 mutations promote breast cancer development³⁸ and may be associated with reduced immune infiltration³⁹. Therefore, a further analysis was performed to correlate the corresponding high-frequency mutated genes with changes in the microenvironment (Figure S5D, E). GATA3 mutations were associated with a lower proportion of CD8+ T cells and an increase in the proportion of M2 macrophages, like the changes in C1 (Figure S5D). NF1 mutations were associated with an increase in the proportion of Treg cells in C2 (Figure S5E). NF1-mutant tumors have been previously associated with DNA damage repair⁴⁰, which may explain the elevated TMB and HLA LOH observed in C2 patients. Interestingly, breast cancer-related mutational signature 1 (clock-like signature), signature 3 [homologous recombination deficiency (HRD)], and signature 13 (activity of APOBEC family of cytidine deaminases) were separately dominant in C1, C2, and C3, respectively (Figure 5G). The elevated HRD score in C2 further suggested a link between CIC deficiency and homologous recombination defects (Figure S5F).

To summarize, the current study has provided a comprehensive perspective on the associations between genomic events and CIC clusters. Specifically, HRD-induced genomic instability and HLA LOH may underlie the observed CIC deficiency in C2.

Sphingolipid metabolic reprogramming influences CIC deficiency

Metabolic reprogramming is a fundamental hallmark of cancer with a pivotal role in tumorigenesis and progression through intricate links with oncogenic signaling. Principal component analysis revealed that breast cancer has distinctive metabolome features from normal samples, while metabolic heterogeneity in CIC subtypes showed subtle differences compared to normal samples (Figure 6A). The divergence between normal and tumor metabolomes was quantified next using two metrics (Euclidean distance and correlation distance) to assess global metabolic shifts across clusters. Both measures showed that the metabolome distance between normal tissues and C2 tumors was markedly greater than observed for the other subtypes (Figure 6B–D), indicating that C2 exhibits the most pronounced metabolic reprogramming. Together with the immunophenotype analyses, these extensive metabolic alterations may underlie the impaired antigen-presentation capacity characteristic of C2.

• Download figure
• Open in new tab
• Download powerpoint

Figure 6

Patterns of metabolism in CIC clusters. (A) Principal component analysis of tumor samples (n = 442) and paired normal samples (n = 76) in the breast cancer cohort for the metabolomics data. (B–D) Global differences in metabolite level between normal tissues and tumors from each CIC cluster in the breast cancer cohort measured by the Euclidean distance (B) and the correlation distance (C). The inset summarizes the average distances between pairs of tissues as a percentage of the average distance between C1 and normal tissues (D). (E) Pathway-based analysis of metabolomic changes among CIC clusters. DA score reflects the overall changes in pathway metabolite levels. A score of 1 indicated increased metabolite levels in one CIC cluster tumors compared to other groups, while a score of −1 indicated decreased levels. (F) Comparison of ssGSEA enrichment scores of metabolism signatures among CIC clusters. (G) GSEA plots of the unsaturated fatty acid metabolic process pathway upregulated in C1 patients. (H) Log2FCs in lipid abundance between C1 tumor tissues and other CIC groups. A Log2FC of 0 (dashed black line) indicated no difference in lipid abundance between the two groups. (I) Volcano plots of enriched metabolism-associated genes in C1 compared to other groups. Boxplots of (B), (C), and (H) show first and third quartiles, whiskers extend to the lowest and highest value within the 1.5× interquartile range, and the line indicates the median. P values were from wilcoxon test adjusted by the Benjamini–Hochberg procedure. FA, fatty acyls; GP, glycerophospholipids; SP, sphingolipids; GL, glycerolipids; ST, sterol lipids.

Next, transcriptome- and metabolome-based metabolic pathway analyses were integrated to evaluate the specific metabolic changes in patients from different CIC clusters (Figures 6E, F and S6A) A comparison of the transcriptomic and metabolomic data revealed distinct enrichment patterns. Analysis of DA scores indicated that C1 patients, characterized by a low proportion of immune cells, exhibited metabolic inertia (Figure 6E), highlighting the contribution of immune cells to metabolic regulation. While carbohydrate metabolism pathways were relatively enriched in both C3 and C2, C3 had higher glycolysis activity than C2, showing a metabolic pattern of highly activated immune cells. C2 also had an upregulation of nucleotide metabolism pathways associated with tumor aggressiveness. In contrast, amino acid metabolism pathways were specifically enriched in each subtype, underscoring the diverse functional roles. Patients in C1 exhibited enrichment of lipid metabolism pathways based on ssGSEA analysis (Figure 6F, G). Correlation analysis further suggested that most lipid metabolism pathways are linked to inhibition of CIC steps (Figure S6A). Upregulated lipid metabolites in C1 were analyzed to identify lipid metabolic factors contributing to the immune-cold phenotype. Sphingolipids were most highly upregulated in C1 among lipid metabolites (Figures 6H and S6B). C1 also showed upregulation of lipid metabolism-related genes, including ACADSB, ELOVL2, and the sphingolipid metabolism gene (DEGS2; Figure 6I). Previous studies have demonstrated that sphingolipid metabolism is associated with tumor progression and may influence anti-tumor immunity via macrophages^41,42.

The heterogeneity of metabolic reprogramming across subtypes was investigated based on metabolic analyses, which highlighted sphingolipid metabolism as a potential contributor to CIC functional impairment.

Identification of the metabolic enzyme, PSAT1, as a potential target of antigen presentation

C2 was characterized as an atypical immune subtype exhibiting moderate CIC activation but marked defects in antigen presentation (Figure 2). Preliminary analyses revealed that C2 displays distinct immune and genomic features (Figures 3 and 5). Metabolomic profiling further showed a pronounced upregulation of amino acid and peptide metabolites in C2 (Figure 7A), accompanied by elevated expression of classic metabolic checkpoints, such as IL4I1 and IDO1 (Figure 7B). Notably, phosphoserine aminotransferase 1 (PSAT1), a rate-limiting enzyme in serine biosynthesis, was significantly upregulated in C2, which was aligned with the strong enrichment of the serine metabolic pathway (Figure 6E, F). Moreover, frequent copy number amplification of PSAT1 in C2 tumors underscored genomic dysregulation, highlighting PSAT1 as a key metabolic driver uniquely elevated in this subtype (Figure S6C).

• Download figure
• Open in new tab
• Download powerpoint

Figure 7

PSAT1 modulates CIC deficiency in C2. (A) Volcano plots of 594 annotated polar metabolites in C2 compared to other groups with metabolites color coded by category. (B) Volcano plots of enriched metabolism-associated genes in C2 compared to other groups. (C) Comparison of genes involved in serine metabolism pathways across different CIC clusters. (D) Correlations of PSAT1 expression with immunoregulatory genes, including IL10, CD274, and TGFB1. (E) Relative mRNA levels of immunoregulatory genes following siRNA-mediated knockdown of PSAT1 in HCC1806, n = 3 independent experiments. Graph bars represent mean ± SD. P values were from Student’s t test. (F) Kaplan–Meier survival curve according to PSAT1 levels for RFS and OS in breast cancer patients. 3-PG, 3-phosphoglycerate; 3-PHP, 3-phosphohydroxypyruvate; α-KG, α-ketoglutarate; 3-PS, 3-phosphoserine.

PSAT1 catalyzes the conversion of 3-phosphohydroxypyruvate and glutamate to 3-phosphoserine and α-ketoglutarate (α-KG) during serine synthesis^43,44. Upregulation of PSAT1 drives cell proliferation, metastasis, and therapy resistance in various cancers^45–47. The distribution of key metabolic enzymes in the serine synthesis pathway across the CIC clusters was next analyzed to evaluate the potential immunoregulatory function of PSAT1 (Figure 7C). Key enzymes involved in serine metabolism were significantly upregulated in C2 patients and the levels of the downstream product (serine) also showed an increasing trend in C2 (Figure S6D). Given that the TME of C2 patients exhibited an increase in Tregs and dysfunctional DCs, PSAT1 expression was significantly positively correlated with the expression of IL10, CD274, and TGFB1, which contribute to the immunoregulatory phenotype (Figure 7D). The associations between other serine metabolism-related genes and immunoregulatory molecules were also determined. The results revealed that PSAT1 had the most significant correlation with functional molecules related to Treg production (Figure S6E). Knocking down PSAT1 in breast cancer cell lines downregulated the expression of IL10, CD274, and TGFB1 (Figures 7E and S6F, G). Overexpression of PSAT1 in PSAT1-knockdown cell lines restored the corresponding phenotypes, further confirming the potential role of PSAT1 in causing CIC dysfunction (Figure S6H). PSAT1 also had prognostic value in breast cancer with high expression indicating poor outcomes (Figure 7F).

Overall, dysregulation of serine metabolism in C2 suggested that targeting PSAT1 might have potential benefits in overcoming CIC defects.