Identification of chondroitin polymerizing factor as a biomarker for predicting immunotherapy response in breast cancer: a bioinformatics analysis of tumor microenvironment

Introduction

Breast cancer (BRCA) is one of the most prevalent malignancies worldwide and is characterized by substantial molecular and clinical heterogeneity (1). Although immune checkpoint blockade (ICB) therapy has shown promising efficacy in multiple cancer types (2), its therapeutic benefit in BRCA remains limited and highly heterogeneous (3). Increasing evidence indicates that the tumor microenvironment (TME) plays a critical role in determining immunotherapy response, and distinct immune subtypes are associated with different clinical outcomes. However, the mechanisms underlying immune suppression in BRCA and their impact on immunotherapy response remain incompletely understood. In particular, reliable biomarkers reflecting immune suppression status are still lacking, highlighting the need to identify novel immune-related biomarkers in BRCA (4,5).

TME is a multifaceted ecosystem, encompassing cancer cells, their surrounding stromal cells, and the extracellular matrix (ECM) (6). This intricate network significantly influences BRCA’s response to immunotherapy. Key elements of the TME, such as cancer-associated fibroblast (CAF), tumor-associated macrophages (TAMs), and myeloid-derived suppressor cells (MDSCs), are crucial for the emergence of resistance to immunotherapy (7). CAFs, in particular, have been implicated in tumorigenesis as they inhibit apoptosis and accelerate cancer cell invasion. Furthermore, CAFs are known to produce ECM proteins, which induce an immunosuppressive state in tumor cells (8). The multifaceted nature of TME complicates the treatment of tumors. To address this, efforts have been made to classify the TME. A notable contribution is the development of a transcriptomic analysis platform by Bagaev et al., in which TME is classified into four subtypes based on transcriptomic data, and the classification effectively predicts the response to immunotherapy across various types of tumors, such as bladder cancer, melanoma, cervical cancer, and gastric cancer (4). Specifically, TME is categorized into four distinct subtypes: immune-enriched (IE)/fibrotic (IE/F), IE/non-fibrotic (IE), fibrotic (F), and immune-depleted (D). Patients with subtype IE exhibit longer overall survival (OS) and disease-free survival (DFS) than those with subtypes F and D. Subtype F, in particular, shows the poorest OS (4). However, the classification cannot predict how BRCA patients will respond to immunotherapy.

The Estimation of STromal and Immune cells in MAlignant Tumours using Expression data (ESTIMATE) algorithm is a precise method that calculates gene expression signatures to evaluate stromal and immune cell infiltration (9). Weighted gene co-expression network analysis (WGCNA) is a method employed to investigate biological networks by analyzing pairwise correlations between variables (10). Both methods are instrumental in understanding TME. For instance, ESTIMATE can be utilized to assess stromal and immune cell infiltration in TME. Subsequently, WGCNA is employed to investigate the correlations between various gene sets and their associations with clinical features (11). The integration of these two methods may help develop effective biomarkers to predict immunotherapy responses.

Therefore, this study aimed to identify TME-related biomarkers associated with immunotherapy response in BRCA through integrated bioinformatics analyses and to evaluate their prognostic value and potential clinical significance. We present this article in accordance with the REMARK reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2026-1-0221/rc).

Methods

Datasets

Figure 1 was a flowchart that illustrates each step in a workflow. This study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. Transcriptomic profiles and clinical information of 1,050 BRCA patients were obtained from The Cancer Genome Atlas (TCGA) database. This retrospective study included patients with histologically confirmed BRCA, available transcriptomic expression data, and complete survival information, while cases lacking key clinical annotations or survival data were excluded from downstream analyses. Two microarray datasets, GSE42568 (n=104) and GSE58812 (n=107), were downloaded from the Gene Expression Omnibus (GEO) to analyze the correlation between chondroitin polymerizing factor (CHPF) level and OS of BRCA patients. Two ICB cohorts, IMvigor210 (12,13), including patients treated with the anti-PD-L1 agent atezolizumab, and Gide2019_PD1 (14), comprising melanoma patients receiving anti-PD-1 therapy, were utilized to confirm the roles of CHPF in predicting ICB response. Treatment response was evaluated according to the original study criteria reported for each cohort. A BRCA single-cell RNA sequencing (scRNA-seq) dataset (EGAS00001004809) was applied to assess CHPF levels in specific cell types and to validate the correlation of CHPF with ICB response (15). Because all data were obtained from publicly accessible databases with de-identified patient information, additional ethical approval or informed consent was not required. Patients with incomplete survival information or missing key clinical variables were excluded from the corresponding analyses, and all available eligible cases were included without prior sample size calculation due to the retrospective nature of this bioinformatics study, and no imputation for missing data was performed.

Figure 1 A flowchart that illustrates each step in a workflow. BRCA, breast cancer; CHPF, chondroitin polymerizing factor; ICB, immune checkpoint blockade; IS, immune score; scRNA-seq, single-cell RNA sequencing; TCGA, The Cancer Genome Atlas; TIDE, tumor immune dysfunction and exclusion; TMB, tumor mutation burden; WGCNA, weighted gene co-expression network analysis.

ESTIMATE algorithm and TME classification

The expression matrix was obtained from the 1,050 TCGA-BRCA samples, and then the immune cell infiltration (expressed as an ImmuneScore) and the presence of stromal cells (expressed as an StromalScore) in these samples were calculated by the “ESTIMATE” R package (16). In brief, ImmuneScore and StromalScore were obtained by single-sample Gene Set Enrichment Analysis (ssGSEA) based on the levels of two distinct sets of 141 genes, which indicate stromal and immune cell infiltration in TME. According to ImmuneScore and StromalScore, the TME of BRCA was classified into the following four subtypes: high ImmuneScore and high StromalScore (IM/S), high ImmuneScore and low StromalScore (IM), low ImmuneScore and high StromalScore (S), and low ImmuneScore and low StromalScore (D).

Transcriptomic signatures

The absolute levels of immune and stromal cells within TME were assessed by microenvironment cell populations (MCP)-counter based on transcriptomic data in TCGA-BRCA (17). PROGENy (Pathway RespOnsive GENe) algorithm was applied to assess pathway activity based on transcriptomic data in TCGA-BRCA (18,19). The immunosuppression genes (ISGs) were downloaded from the DisGeNET database (https://www.disgenet.org), and immune suppression scores (ISSs) were calculated using these ISGs (20). The TIMER and CIBERSORT algorithms were employed for analyzing the associations of CHPF level with tumor-infiltrating immune cells (TIICs) and stromal cells in TME (21,22). Moreover, the correlation of CHPF level with molecular markers of TIICs was explored through correlation modules.

Prediction of immunotherapy response

Tumor mutation burden (TMB) (23), neoantigen load (24), and tumor immune dysfunction and exclusion (TIDE) (25), three promising immunotherapy biomarkers, were applied to predict ICB response (26). Transcriptomic data from TCGA-BRCA were processed by normalizing to transcripts per million, applying a log₂ transformation, and excluding low-expression genes, as previously described (12). TMB was determined by counting the number of nonsynonymous mutations per megabase of analyzed genomic sequence. TIDE score was obtained from the TIDE website (http://tide.dfci.harvard.edu/). Neoantigen load was determined by identifying new amino acid 9mers and 10mers arising from mutations in genes with a median expression TPM value >10 in the TCGA-BRCA cohort, as previously described (27).

WGCNA

A total of 15,951 genes were applied to build the co-expression network with the WGCNA R package (28,29). To create a scale-free network, the soft threshold power value (β) was determined using the pickSoftThreshold function. Subsequently, the gene expression similarity matrix was established by calculating Pearson’s association coefficient between gene pairs and transformed into an adjacent matrix with the power adjacent function. The adjacency matrix was then transformed into a topological overlap matrix (TOM), and hierarchical clustering was performed to group similar genes into a module based on TOM-based dissimilarity, with a minimum module size of 29 for the gene dendrogram. To evaluate the association between co-expression modules and clinical traits, module eigengenes (MEs), representing the first principal component of each module, were calculated. The relationship between these MEs and the four TME subtypes (IM/S, IM, S, and D) was assessed using Pearson’s correlation analysis. To further explore their relevance to ICB response, we further calculated the correlation between MEs and TMB and TIDE scores using Pearson’s correlation analysis.

Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) analysis

GO analyses were carried out by the ClusterProfiler R package to identify biological processes of each module. KEGG analyses were performed on a hypergeometric distribution by the R ‘phyper’ function. GO categories and KEGG pathways with a P value of <0.05 were enriched.

Cox regression analysis

Prognosis-related genes were identified by univariate Cox regression analysis using the “Limma” and “survival” R packages. Hazard ratio (HR) with a 95% confidence interval (95% CI) was applied to identify genes as protective (HR <1) or risky (HR >1), and HRs with corresponding 95% CIs were regarded as the primary measures for reporting prognostic effects in subsequent analyses.

Survival analysis

Survival curves for the TCGA-BRCA, GSE58812, GSE42568, IMvigor210, and Gide2019_PD1 cohorts were plotted using the Kaplan-Meier (K-M) method. OS was defined as the time from initial diagnosis or treatment initiation to death from any cause, and OS was regarded as the primary study endpoint, with follow-up time derived from the survival information provided in each original dataset. In the TCGA-BRCA, GSE58812, and GSE42568 cohorts, patients were divided into the high and low CHPF expression groups according to the median CHPF levels. In the IMvigor210 and Gide2019_PD1 cohorts, patients were divided into the high and low CHPF expression groups according to the optimal cutoff values of the CHPF levels determined using maximally selected rank statistics implemented in the “survminer” package. The “survminer” R package was employed for survival analyses.

Differentially expressed genes (DEGs)

The DEGs in RNA-seq data were identified by edgeR and DEseq2 R packages. Genes were deemed as differentially expressed if they had an FDR below 0.05 and a fold change (|FC|) greater than 2.

Data analysis of scRNA-seq

A scRNA-seq dataset, EGAD00001006608, was obtained from the European Genome-phenome Archive (EGA). The raw data were first preconditioned by the “Seurat” R package, followed by dimensionality reduction through principal components analysis (PCA) algorithm. Furthermore, the K-Nearest Neighbor method was utilized for classification processing. The data analysis was conducted as previously detailed (30).

Statistics analysis

R packages were primarily used for statistical analysis, bioinformatics, and image visualization. Survival curves were plotted using the K-M method, and their significance was assessed with log-rank tests. Univariate and multivariate Cox proportional hazards regression analyses were performed to evaluate the independent prognostic value of CHPF. Variables with P<0.05 in univariate analysis were included in the multivariate model. All statistical tests were two-sided. One-way analysis of variance (ANOVA) was applied to compare the differences among the four groups, with a P value below 0.05 considered statistically significant.

Results

Four distinct TME subtypes in BRCA

Previous studies have demonstrated that TME in many types of tumors can be categorized into four distinct subtypes: IE/F, IE, F, and D (4,12). Although the TME classifications can be used to predict immunotherapy efficacy, OS, and DFS across 20 different cancers (4), they are not associated with OS (Figure S1A) and patient responses to immunotherapy (Figure S1B-S1D) in BRCA. In the study, the TME of BRCA was classified into the following four subtypes by ESTIMATE algorithm: IM/S, IM, S, and D. As expected, subtypes IM/S and IM exhibited significantly higher immune scores compared with subtypes S and D, and there was no immune cell infiltration in subtypes S and D (Figure S1E). Subtypes IM/S and S exhibited significantly higher stromal scores compared with subtypes IM and D, and there was no stromal cell infiltration in subtypes IM and D (Figure S1F). It is noteworthy that the TME classification can be used to predict the 25-year OS of BRCA patients (Figure 2A). Patients with subtype IM had a significantly longer 25-year OS than those with subtypes S and D, while subtype IM/S showed the poorest OS. Although the 10-year OS did not differ between patients with subtypes IM/S and IM (Figure 2B), the survival probability beyond 10 years was lower in the patients with subtype IM/S compared to those with subtype IM (Figure 2C). To uncover the underlying reasons behind this, immunosuppression genes (ISGs) were obtained from the DisGeNET database, and immune suppression scores (ISSs) were calculated using these ISGs (Table S1) (20,31). Figure 2D showed that the ISSs of subtype IM/S were significantly higher than that of subtype IM. A total of 1,050 patients from the TCGA-BRCA cohort and additional patients from the external validation cohorts were included in the survival analyses, and the corresponding numbers of outcome events were obtained from the original survival data of each dataset.

Figure 2 Four distinct TME subtypes in BRCA. (A) 25-year OS of BRCA patients stratified by TME subtype classification in the TCGA-BRCA cohort. (B) 10-year OS between BRCA patients with subtypes IM/S and IM. (C) The survival probability beyond 10 years between BRCA patients with subtypes IM/S and IM. (D) ISSs were calculated using ISGs (https://www.disgenet.org). Violin plots illustrating variations in neoantigen load (E), TMB (F), and TIDE (G) among the four distinct subtypes. (H) The MCP-counter algorithm was used to analyze the differences of immune cell infiltration among the four distinct subtypes. (I) Heatmap displaying signaling pathway activity scores as determined by the PROGENy algorithm. (J) The differences in 12 TEM-related gene signatures were analyzed by ssGSEA algorithm among the four distinct subtypes. ANOVA, analysis of variance; BRCA, breast cancer; CAF, cancer-associated fibroblast; D, low ImmuneScore and low StromalScore; EMT, epithelial-mesenchymal transition; IM, high ImmuneScore and low StromalScore; IM/S, high ImmuneScore and high StromalScore; ISGs, immunosuppression genes; ISSs, immune suppression scores; MCP, microenvironment cell populations; MDSC, myeloid-derived suppressor cell; MHCI, major histocompatibility complex class I; NK, natural killer; OS, overall survival; S, low ImmuneScore and high StromalScore; ssGSEA, single-sample Gene Set Enrichment Analysis; TCGA, The Cancer Genome Atlas; TEM, T effector memory; TIDE, tumor immune dysfunction and exclusion; TMB, tumor mutation burden; TME, tumor microenvironment.

Subsequently, the roles of the TME classification in predicting ICB response were investigated based on the following biomarkers: neoantigen load, TMB, and TIDE. As shown in Figure 2E-2G, the TME classification can predict patient responses to immunotherapy in the TCGA-BLCA cohort. Especially, subtype IM/S had significantly lower neoantigen load (Figure 2E), less TMB (Figure 2F), and higher TIDE scores (Figure 2G) compared with subtype IM, indicating that the patients with subtype IM/S had worse ICB response than those with subtype IM.

Then, the MCP-counter algorithm was used to examine the differences in immune cells between subtype IM/S and subtype IM. Figure 2H revealed that more fibroblasts and endothelial cells were enriched in subtype IM/S than in subtype IM. The results from PROGENy pathway scores revealed that immunosuppressive and stemness-related pathways such as MAPK, TGF-β, and Wnt (32) were significantly activated in subtype IM/S compared with subtype IM (Figure 2I). Furthermore, the differences in 12 TEM-related gene signatures between subtype IM/S and subtype IM were analyzed (Figure 2J). The pro-tumor (cytokines, macrophages, MDSCs, etc.) and angiogenesis-related (angiogenesis, CAF, EMT, etc.) gene signatures were markedly elevated in subtype IM/S compared to subtype IM. Conversely, anti-tumor gene signatures [cytokines, effector cells, MHC-1, natural killer (NK) cells, etc.] were markedly decreased in subtype IM/S compared with subtype IM. These data demonstrate that stromal cell infiltration in subtype IM/S BRCA causes an immunosuppressive TME, which is associated with poor outcome and immunotherapy efficacy, compared with subtype IM BRCA.

Construction of co-expression modules using WGCNA

To reveal the molecular mechanisms responsible for immunosuppression, transcriptomic data were downloaded from the TCGA-BRCA dataset, and 15,951 genes (see table available online: https://cdn.amegroups.cn/static/public/tcr-2026-1-0221-1.xlsx) were used to construct co-expression networks with the WGCNA R software package. Herein, the soft threshold power value was selected as β=18 (scale-free R²=0.85 and mean connectivity =11.31) to assure a scale-free network (Figure 3A). Consequently, 22 co-expression modules, with gene counts varying from 37 to 4,235 (see table available online: https://cdn.amegroups.cn/static/public/tcr-2026-1-0221-2.xlsx), were yielded through clustering the genes with similar expression patterns (Figure 3B,3C). An adjacency heatmap displayed the TOM for all 15,951 genes, indicating that each module was independently validated in relation to the others (Figure 3D). The association of MEs with clinical traits was further analyzed. As shown in Figure 3E, the MEs in the green_M5 (r=0.49, P=4e−65), darkgreen_M11 (r=0.35, P=3e−31), darkred_M10 (r=0.47, P=1e−57), royalblue_M18 (r=0.26, P=2e−17), black_M4 (r=0.61, P=4e−109), and tan_M17 (r=0.4, P=4e−42) modules showed a significantly higher correlation with subtype IM/S BRCA compared with other modules. To further explore their relevance to immunotherapy, we next examined the correlations between MEs and TMB and TIDE scores. As shown Figure 3F,3G, several of these modules were significantly associated with TMB and TIDE. Notably, green_M5 showed a negative correlation with TMB and a positive correlation with TIDE, suggesting its potential role in reduced tumor immunogenicity and enhanced immune evasion.

Figure 3 Construction of co-expression modules using WGCNA. (A) Network topology analysis at various soft-thresholding powers showed that when β=18, the coexpression network exhibited a scale-free topology, with R²=0.85 and mean connectivity =11.31. (B) The dendrogram illustrating co-expression network modules derived from topological overlap. (C) The cluster dendrogram derived from module eigengenes. (D) An adjacency heatmap illustrating the topological overlap matrix for 15,951 genes. (E) Heatmap of module-trait relationships, with rows representing module eigengenes and columns representing clinical traits. Each cell contained the Pearson correlation coefficient (r) and P value. (F,G) The correlations between module eigengenes and TMB and TIDE scores. D, low ImmuneScore and low StromalScore; IM, high ImmuneScore and low StromalScore; IM/S, high ImmuneScore and high StromalScore; S, low ImmuneScore and high StromalScore; TIDE, tumor immune dysfunction and exclusion; TMB, tumor mutation burden; WGCNA, weighted gene co-expression network analysis.

A key module causing immunosuppression in the IM/S subtype

Then, genes in these modules were used to carry out KEGG pathway enrichment analysis. The anti-tumor and stemness-related pathways, including PI3K/Akt, TGF-β, Hippo, MAPK, and focal adhesion signaling (32-35) were enriched in green_M5 module (Figure 4A). The MAPK and PI3K/Akt pathways were enriched in the darkred_M10 module (Figure 4B). The chemokine and cell adhesion signaling pathways were enriched in the black_M4 module (Figure 4C). The scatterplot analysis revealed a significant correlation between GS and MM values across three modules, indicating a high correlation between key genes within these modules and poor outcomes (Figure 4D-4F). Subsequently, the number and proportion of genes significantly associated with OS in these modules and survival of BRCA patients were analyzed. As shown in Figure 4G, 22, 4, 9, 0, 0, and 0 genes in green_M5, darkred_M10, black_M4, darkgreen_M11, royalblue_M18, and tan_M17 modules were associated with patient survival, respectively, with P<0.05 and HR >1 (see table available online: https://cdn.amegroups.cn/static/public/tcr-2026-1-0221-3.xlsx). The green_M5 module was more significantly associated with poor prognosis and involved 33.8% of the genes. The green_M5 module was positively correlated with critical pro-tumor signaling, including CAF, EMT, protumor cytokines, angiogenesis, and MDSC infiltration (Figure 4H). Therefore, it was chosen for a subsequent detailed investigation.

Figure 4 A key module causing immunosuppression in IM/S subtype. KEGG enrichment analysis of genes within the M5 (A), M10 (B), and M4 (C) modules. The scatterplot illustrates the correlation of GS with MM values in the M5 (D), M10 (E), and M4 (F) modules. (G) The number and proportion of genes significantly associated with OS in these modules are shown. HR with a 95% CI was utilized to identify genes as protective (HR <1) or risky (HR >1). (H) The correlation of M5 module with pro-tumor signaling (CAF, EMT, protumor cytokines, angiogenesis, MDSC, etc.) was analyzed. CAF, cancer-associated fibroblast; CI, confidence interval; ECM, extracellular matrix; EMT, epithelial-mesenchymal transition; GS, gene significance; HR, hazard ratio; IM/S, high ImmuneScore and high StromalScore; KEGG, Kyoto Encyclopedia of Genes and Genomes; MDSC, myeloid-derived suppressor cell; MHCI, major histocompatibility complex class I; MM, module membership; NK, natural killer; OS, overall survival.

Up-regulated CHPF in the green_M5 module was correlated with poor outcome

To identify key genes in the green_M5 module associated with a worse prognosis, the association between the expression of 22 genes and OS of BRCA patients was analyzed. As shown in Figure 5A, the relationship of syndecan 1 (SDC1) and CHPF with poor prognosis was the most significant (HR: 1.62; 95% CI: 1.2–2.2; P<0.01), and multivariate Cox regression analysis further confirmed that CHPF remained an independent prognostic factor after adjustment for available clinical variables. SDC1 is crucial in regulating lymphocyte infiltration and severs as a predictor of immunotherapy response in BRCA (36,37). Interestingly, CHPF is instrumental in sustaining SDC1 expression and promoting micropinocytosis (38), which is important in supporting tumor cell growth under nutrient-deprived conditions. Thus, the role of CHPF in predicting prognosis and ICB response in BRCA was next investigated. In the TCGA-BRCA cohort, CHPF levels were increased in tumor tissues compared to normal controls (Figure 5B) and showed a positive correlation with disease progression (Figure 5C). CHPF levels in subtype IM/S were elevated compared to those in subtype IM (Figure 5D), indicating that upregulated CHPF expression is a potential mechanism for the poor prognosis of subtype IM/S BRCA. Moreover, high CHPF levels were markedly associated with poor OS and DFS in the TCGA-BRCA cohort (Figure 5E,5F). The correlation between CHPF level and OS in BRCA was further analyzed across multiple cohorts of BRCA patients. Indeed, the BRCA patients with high CHPF levels exhibited worse OS in the GSE42568 (Figure 5G) and GSE58812 (Figure 5H) datasets.

Figure 5 Up-regulated CHPF was correlated with poor outcome. (A) Univariate Cox regression analysis was used to determine the correlation between the levels of these 22 genes and OS of BRCA patients. All genes have an HR value >1 and a P value <0.05. (B) CHPF expression was analyzed in tumor tissues and normal controls in the TCGA-BRCA cohort. (C) CHPF expression was analyzed in patients with stage I/II and stage III/IV in the TCGA-BRCA cohort. (D) CHPF expression was analyzed among the four distinct subtypes in the TCGA-BRCA cohort. 25-year OS (E) and DFS (F) between BRCA patients with high CHPF expression and those with low CHPF expression. OS of between BRCA patients with high CHPF expression and those with low CHPF expression in the GSE42568 (G) and GSE58812 (H) datasets. *, P<0.05; **, P<0.01. BRCA, breast cancer; CHPF, chondroitin polymerizing factor; CI, confidence interval; D, low ImmuneScore and low StromalScore; DFS, disease-free survival; F, fibrotic; HR, hazard ratio; IM, high ImmuneScore and low StromalScore; N, normal; OS, overall survival; T, tumor; TCGA, The Cancer Genome Atlas.

CHPF was correlated with immunosuppressive TME in BRCA

The correlations of CHPF expression with immune and stromal cell infiltration, immune biomarker expression, and signaling pathway activation were next investigated. Figure 6A showed that CHPF levels were significantly negatively associated with NK cell, cytotoxic T lymphocyte, and CD8 T cell infiltration, while being positively associated with CAF, endothelial cell, and monocytic lineage infiltration. We further analyzed the association between CHPF level and the expression of TAM markers, including CD68, CD163, V-set immunoglobulin domain-containing 4 (VSIG4), membrane spanning four domains, subfamily A, member 4A (MS4A4A), CCL2, and IL-10. Figure 6B revealed that high CHPF expression was associated with the expression of these TAM markers. We further analyzed the correlations between CHPF expression and immune checkpoint molecules. As shown in Figure 6C,6D, CHPF expression was significantly associated with multiple immune checkpoint genes, including both stimulatory and inhibitory checkpoints. Specifically, CHPF expression was positively correlated with key inhibitory checkpoint molecules, including CD274 (PD-L1), TIGIT, CD80, HAVCR2, LAIR1, LGALS3, and VSIR (Figure 6D). The DEGs between the high and low CHPF expression groups (see table available online: https://cdn.amegroups.cn/static/public/tcr-2026-1-0221-4.xlsx) were used to reveal the signaling pathways and biological processes linked to CHPF. KEGG analysis displayed that upregulated CHPF was associated with PI3K/Akt, ECM dysregulation, and focal adhesion signaling pathways (Figure 6E). GO analyses revealed that low CHPF levels were markedly associated with anti-tumor immune responses such as lymphocyte activation, phagocytosis, complement activation, and B cell activation (Figure 6F). High CHPF expression had significantly higher TIDE scores (Figure 6G) and lower TMB scores (Figure 6H). These data demonstrate that low CHPF expression might predict good immunotherapy efficacy in BRCA.

Figure 6 CHPF was correlated with immunosuppressive TME in BRCA. (A) The correlations of CHPF level with immune and stromal cell infiltration in BRCA were analyzed using the TIMER and CIBERSORT algorithms. (B) The correlations of CHPF level with immunological markers in BRCA were analyzed using the TIMER algorithm. (C) The correlations of CHPF level with stimulatory checkpoints in BRCA. (D) The correlations of CHPF level with stimulatory checkpoints in BRCA. The TCGA-BRCA patients were divided into the CHPF high expression group and the CHPF low expression group, and then the DEGs between the two groups were identified (FDR <005 and |FC| >2) for KEGG (E) and GO (F) enrichment analysis. (G,H) Violin plot illustrating variations in TIDE (G) and TMB (H) between the CHPF high expression group and the CHPF low expression group in the TCGA-BRCA. (I) Eight clusters were annotated by the cell type-specific markers in the scRNA-seq dataset, EGAD00001006608, and CHPF was identified to be mainly enriched in fibroblasts. (J) CHPF expression was markedly higher in the NR group than in the R group, regardless of pre- or on- anti-PD-1 treatment. (K) In the EGAD00001006608 dataset, fibroblasts were divided into two groups: CHPF-positive and CHPF-negative. The DEGs between the two groups were identified (FDR <005 and |FC| >2) for KEGG enrichment analysis. BRCA, breast cancer; CHPF, chondroitin polymerizing factor; DEGs, differentially expressed genes; ECM, extracellular matrix; FC, fold change; FDR, false discovery rate; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; NR, non-responder; R, responder; TCGA, The Cancer Genome Atlas; TIDE, tumor immune dysfunction and exclusion; TMB, tumor mutation burden; TME, tumor microenvironment.

Furthermore, a BRCA scRNA-seq dataset was used to validate the correlation of CHPF with ICB response. Eight clusters were annotated by cell type-specific markers (Figure S2A) in accordance with a previous study (15), obtaining eight cell types, namely B cell, cancer cell, endothelial cell, fibroblast, mast cell, myeloid cell, plasmacytoid dendritic cell (pDC), and T cell (Figure S2B). CHPF was enriched in fibroblasts (Figure 6I and Figure S2C), indicating that CAFs, the major stromal component in TEM, were the main source of CHPF. More importantly, CHPF levels were markedly higher in the NR group compared to those in the R group, regardless of pre- or on- anti-PD-1 treatment (Figure 6J and Figure S2C).

To reveal the potential roles of CHPF in regulating fibroblast function in this dataset, the cells were divided into CHPF-positive fibroblasts and CHPF-negative fibroblasts, and then DEGs were identified between the two groups (see table available online: https://cdn.amegroups.cn/static/public/tcr-2026-1-0221-5.xlsx). KEGG analysis showed that CHPF levels were associated with ECM remodeling (Figure 6K), which accelerated tumor growth and angiogenesis (39). Additionally, ECM dysregulation was associated with CAF-driven TGF-β signaling, resulting in an immunosuppressive TME by inhibiting immune cell activity and resulted in the ineffectiveness of PD-1 blockade (40).

CHPF was a predictor of immunotherapy response

Finally, two independent cohorts of ICB-treated tumor patients, IMvigor210 and Gide2019_PD1, were used to validate the CHPF performance in predicting immunotherapeutic response. In the IMvigor210 cohort, CHPF expression was significantly higher in the NR group compared to the R group (Figure 7A). The role of CHPF in predicting immunotherapeutic outcomes was next analyzed. Figure 7B showed that high CHPF levels were markedly correlated with poor OS after ICB therapy. In the Gide2019_PD1 cohort, the high CHPF expression group exhibited a worse OS (Figure 7C) and PFS (Figure 7D) compared with the low CHPF expression group. These data demonstrate that CHPF may act as a promising predictor of immunotherapeutic response in BRCA.

Figure 7 CHPF was a predictor of immunotherapy response. (A) The correlation between CHPF level and ICB response was analyzed in the IMvigor210 cohort. (B) The correlation between CHPF level and 25-year OS in the IMvigor210 cohort was analyzed after ICB therapy. The correlations between CHPF level and OS (C) or PFS (D) in the Gide2019_PD1 cohort were analyzed after ICB therapy. CHPF, chondroitin polymerizing factor; ICB, immune checkpoint blockade; NR, non-responder; OS, overall survival; PFS, progression-free survival; R, responder.

Discussion

ICBs, particularly PD-1 (nivolumab, cemiplimab, pembrolizumab, etc.) (41), PD-L1 (atezolizumab, durvalumab, avelumab, etc.) (42), and CTLA-4 (ipilimumab, tremelimumab, botensilimab, etc.) (43), have revolutionized the therapeutic strategy of various malignancies. Nivolumab, the pioneer PD-1 inhibitor, has been empirically shown to enhance OS and DFS in conjunction with chemotherapy, outperforming the outcomes of chemotherapy alone (44). Complementarily, ipilimumab synergizes with nivolumab to stimulate de novo anti-tumor T cell responses (45). However, the response to ICB treatment varies from patient to patient.

Recent research indicates improved patient outcomes with personalized treatment, highlighting the potential of precision medicine in advancing cancer therapies (46,47). While genomic analysis is increasingly recognized as a key component in clinical decision-making, it often involves the use of targeted panels that encompass a limited number of genes, thereby capturing only a subset of cancer-causing alterations. RNA-seq presents an additional avenue to unravel the intricacies and diversity of tumors, and to identify new biomarkers for the creation of innovative therapeutic approaches. The technique provides a comprehensive perspective on tumor attributes, potentially facilitating the further recognition and refinement of personalized cancer treatments (48). By analyzing the transcriptomic data, a recent study classifies TME into four subtypes: IE/F, IE, F, and D (4). The classification can effectively forecast the response to immunotherapy in multiple types of tumors (4). Patients with subtype IE exhibit significantly longer OS and DFS compared to those with other subtypes. However, this predictive approach is not effective in BRCA.

These findings suggest that the biological heterogeneity of BRCA may limit the direct applicability of previously established TME classification systems. Therefore, it remains necessary to develop BRCA-specific TME stratification strategies and to identify key molecular determinants associated with prognosis and immunotherapy response.

The vast amount of data generated by large-scale exome and transcriptome sequencing can be daunting and impractical for routine treatment decision-making support. To address this, the ESTIMATE algorithm can be used to analyze immune and stromal cell infiltration in TME, which classifies TME into distinct subtypes (49-52). In the present study, the TME of BRCA is classified into four subtypes by ImmuneScore and StromalScore: IM/S, IM, S, and D. Interestingly, the classification can predict the response to immunotherapy in BRCA. Specifically, the BRCA patients with subtype IM exhibit significantly longer 25-year OS compared with the patients with subtypes S and D, with subtype IM/S showing the worst OS. Furthermore, while there was no difference in 10-year OS between patients with subtypes IM/S and IM, the survival probability beyond 10 years was lower for those with subtype IM/S compared to subtype IM. Given that the ImmuneScore between subtypes IM/S and IM is the same, the differences between the two subtypes should exist in the details of the StromalScore. To reveal the potential mechanisms responsible for immunosuppression, WGCNA was applied to identify the critical modules and the key genes in these modules. As a result, SDC1 and CHPF were identified as two key genes in the process, and they were most associated with poor prognosis in BRCA patients.

CHPF plays a crucial role in chondroitin sulfate chain elongation (53). Recent studies have suggested that CHPF functions as an oncogene in various types of tumors. Upregulated CHPF accelerates tumor cell proliferation and invasion in many types of tumors (54-56). In BRCA, CHPF overexpression is correlated with poor prognosis and has been shown to regulate the formation of chondroitin sulfate chains. In this study, we demonstrated that CHPF contributes to an immunosuppressive TME by decreasing the infiltration of NK cells and cytotoxic T lymphocytes while accelerating the recruitment of CAFs and endothelial cells. Mechanistically, this phenomenon may be explained by the role of CHPF in sustaining SDC1 expression (38). Previous reports have indicated that SDC1 regulates the TGF-β1/Smad pathway, which is a known driver of EMT and a potent suppressor of lymphocyte infiltration in BRCA (57). Furthermore, our scRNA-seq analysis revealed that CHPF is primarily enriched in CAFs and is associated with ECM remodeling. As suggested by existing research, an altered ECM can act as a physical barrier that restricts the infiltration and effector functions of NK cells and cytotoxic lymphocytes (40). Additionally, the positive association between CHPF and the activation of PI3K/Akt and focal adhesion signaling pathways further suggests that CHPF orchestrates a pro-tumor environment. The activation of these pathways has been shown to promote tumor cell survival and coordinate the secretion of immunosuppressive cytokines, ultimately leading to the observed resistance to ICB treatment (37). Collectively, these data demonstrate that high CHPF expression is significantly correlated with poor OS after ICB treatment.

Despite these findings, the precise mechanisms by which CHPF causes immunosuppressive TME in BRCA remain unclear. Future research should focus on elucidating the signaling pathways regulated by CHPF in the process. Moreover, several limitations should be acknowledged in the present study. First, this study was primarily based on retrospective analyses of publicly available transcriptomic datasets, which may introduce potential selection bias. Second, the biological functions and molecular mechanisms of CHPF in shaping the immunosuppressive TME were inferred from bioinformatics analyses and still require further experimental validation. Third, prospective clinical studies with larger and more homogeneous cohorts are needed to confirm the clinical utility of CHPF as a predictive biomarker for immunotherapy response in BRCA.

Conclusions

These data demonstrate that CHPF could potentially act as a predictor for immunotherapy response in BRCA. Furthermore, CHPF may serve as a promising biomarker for stratifying patients and guiding individualized immunotherapeutic strategies. However, prospective clinical validation and mechanistic investigations are still required before its clinical application. Overall, this study provides novel transcriptomic evidence linking CHPF to immunosuppressive TME remodeling and immunotherapy response in BRCA.

Acknowledgments

None.

Footnote

Reporting Checklist: The authors have completed the REMARK reporting checklist. Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2026-1-0221/rc

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2026-1-0221/prf

Funding: None.

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2026-1-0221/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. This study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Bray F, Laversanne M, Sung H, et al. Global cancer statistics 2022: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin 2024;74:229-63. [Crossref] [PubMed]
2. Zheng X, Fang Z, Liu X, et al. Increased vessel perfusion predicts the efficacy of immune checkpoint blockade. J Clin Invest 2018;128:2104-15. [Crossref] [PubMed]
3. Morad G, Helmink BA, Sharma P, et al. Hallmarks of response, resistance, and toxicity to immune checkpoint blockade. Cell 2021;184:5309-37. [Crossref] [PubMed]
4. Bagaev A, Kotlov N, Nomie K, et al. Conserved pan-cancer microenvironment subtypes predict response to immunotherapy. Cancer Cell 2021;39:845-865.e7. [Crossref] [PubMed]
5. Mei J, Cai Y, Xu R, et al. Conserved immuno-collagenic subtypes predict response to immune checkpoint blockade. Cancer Commun (Lond) 2024;44:554-75. [Crossref] [PubMed]
6. Ghosh S, Li X, Rawat K, et al. Chemokine-defined macrophage niches establish spatial organization of tumor immunity. Nat Immunol 2026;27:715-24. [Crossref] [PubMed]
7. Kundu M, Butti R, Panda VK, et al. Modulation of the tumor microenvironment and mechanism of immunotherapy-based drug resistance in breast cancer. Mol Cancer 2024;23:92. [Crossref] [PubMed]
8. Pei L, Liu Y, Liu L, et al. Roles of cancer-associated fibroblasts (CAFs) in anti- PD-1/PD-L1 immunotherapy for solid cancers. Mol Cancer 2023;22:29. [Crossref] [PubMed]
9. He M, Wu Y, Liu S, et al. Pancancer Analysis and the Oncogenic Role of UBTF in Breast Invasive Carcinoma. Int J Mol Sci 2026;27:2909. [Crossref] [PubMed]
10. Wang J, Li X, Chen S, et al. Identification of the role of MCM6 in bladder cancer prognosis, immunotherapy response, and in vitro experimental investigation using multi-omics analysis. Life Sci 2023;335:122253. [Crossref] [PubMed]
11. Zheng H, Liu H, Li H, et al. Weighted Gene Co-expression Network Analysis Identifies a Cancer-Associated Fibroblast Signature for Predicting Prognosis and Therapeutic Responses in Gastric Cancer. Front Mol Biosci 2021;8:744677. [Crossref] [PubMed]
12. Xiao Y, Dong Y, Yu T, et al. Characterization of the immune related lncRNAs in bladder cancer to aid immunotherapy. Front Immunol 2022;13:941189. [Crossref] [PubMed]
13. Mariathasan S, Turley SJ, Nickles D, et al. TGFβ attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells. Nature 2018;554:544-8. [Crossref] [PubMed]
14. Gide TN, Quek C, Menzies AM, et al. Distinct Immune Cell Populations Define Response to Anti-PD-1 Monotherapy and Anti-PD-1/Anti-CTLA-4 Combined Therapy. Cancer Cell 2019;35:238-255.e6. [Crossref] [PubMed]
15. Bassez A, Vos H, Van Dyck L, et al. A single-cell map of intratumoral changes during anti-PD1 treatment of patients with breast cancer. Nat Med 2021;27:820-32. [Crossref] [PubMed]
16. Dai H, Liu B, Rui H, et al. INPP5D/SHIP1-mediated immunometabolic remodeling of renal monocytes in idiopathic membranous nephropathy. Front Immunol 2026;17:1755723. [Crossref] [PubMed]
17. de Vries-Brilland M, Rioux-Leclercq N, Meylan M, et al. Comprehensive analyses of immune tumor microenvironment in papillary renal cell carcinoma. J Immunother Cancer 2023;11:e006885. [Crossref] [PubMed]
18. Reinscheid M, Weisser J, Pascual Maier N, et al. Acute and chronic infections drive distinct trajectories in human memory CD4(+) T cell formation. Immunity 2026;S1074-7613(26)00119-6.
19. Schubert M, Klinger B, Klünemann M, et al. Perturbation-response genes reveal signaling footprints in cancer gene expression. Nat Commun 2018;9:20. [Crossref] [PubMed]
20. Piñero J, Ramírez-Anguita JM, Saüch-Pitarch J, et al. The DisGeNET knowledge platform for disease genomics: 2019 update. Nucleic Acids Res 2020;48:D845-55. [Crossref] [PubMed]
21. Gu Z, He C, Cheng S, et al. Mature tertiary lymphoid structures tumor microenvironment-based risk model to assess patients with pancreatic ductal adenocarcinoma. Transl Cancer Res 2026;15:128. [Crossref] [PubMed]
22. Newman AM, Liu CL, Green MR, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods 2015;12:453-7. [Crossref] [PubMed]
23. Zhao X, Pan X, Wang Y, Zhang Y. Targeting neoantigens for cancer immunotherapy. Biomark Res 2021;9:61. [Crossref] [PubMed]
24. Cristescu R, Mogg R, Ayers M, et al. Pan-tumor genomic biomarkers for PD-1 checkpoint blockade-based immunotherapy. Science 2018;362:eaar3593. [Crossref] [PubMed]
25. Zheng K, Hai Y, Chen H, et al. Tumor immune dysfunction and exclusion subtypes in bladder cancer and pan-cancer: a novel molecular subtyping strategy and immunotherapeutic prediction model. J Transl Med 2024;22:365. [Crossref] [PubMed]
26. Carmona J. Biomarkers for checkpoint blockade. Nat Med 2018;24:1637. [Crossref] [PubMed]
27. Rooney MS, Shukla SA, Wu CJ, et al. Molecular and genetic properties of tumors associated with local immune cytolytic activity. Cell 2015;160:48-61. [Crossref] [PubMed]
28. Yin M, Zhang D, Wang T, et al. Identification and evaluation of glutamine-related gene characteristics based on multi-omics to predict the prognosis of patients with colorectal cancer. J Transl Med 2026;24:239. [Crossref] [PubMed]
29. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 2008;9:559. [Crossref] [PubMed]
30. Zhang X, Du Y, Xiong W, et al. Combined single-cell RNA-seq and bulk RNA-seq to analyze the expression and role of TREM2 in bladder cancer. Med Oncol 2022;40:23. [Crossref] [PubMed]
31. Li S, Mai Z, Gu W, et al. Molecular Subtypes of Oral Squamous Cell Carcinoma Based on Immunosuppression Genes Using a Deep Learning Approach. Front Cell Dev Biol 2021;9:687245. [Crossref] [PubMed]
32. Shi H, Yang Y. Identification of inhibitory immune checkpoints and relevant regulatory pathways in breast cancer stem cells. Cancer Med 2021;10:3794-807. [Crossref] [PubMed]
33. Gu M, Liu Y, Xin P, et al. Fundamental insights and molecular interactions in pancreatic cancer: Pathways to therapeutic approaches. Cancer Lett 2024;588:216738. [Crossref] [PubMed]
34. Li A, Chang Y, Song NJ, et al. Selective targeting of GARP-LTGFβ axis in the tumor microenvironment augments PD-1 blockade via enhancing CD8+ T cell antitumor immunity. J Immunother Cancer 2022;10:e005433. [Crossref] [PubMed]
35. Shen J, Cao B, Wang Y, et al. Hippo component YAP promotes focal adhesion and tumour aggressiveness via transcriptionally activating THBS1/FAK signalling in breast cancer. J Exp Clin Cancer Res 2018;37:175. [Crossref] [PubMed]
36. Zhong Y, Li F, Zhang S, et al. Syndecan-1 as an immunogene in Triple-negative breast cancer: regulation tumor-infiltrating lymphocyte in the tumor microenviroment and EMT by TGFb1/Smad pathway. Cancer Cell Int 2023;23:76. [Crossref] [PubMed]
37. Li C, Yang L, Zhang Y, et al. Integrating single-cell and bulk transcriptomic analyses to develop a cancer-associated fibroblast-derived biomarker for predicting prognosis and therapeutic response in breast cancer. Front Immunol 2023;14:1307588. [Crossref] [PubMed]
38. Yen HR, Liao WC, Chen CH, et al. Targeting chondroitin sulfate suppresses macropinocytosis of breast cancer cells by modulating syndecan-1 expression. Mol Oncol 2024;18:2569-85. [Crossref] [PubMed]
39. Wu F, Yang J, Liu J, et al. Signaling pathways in cancer-associated fibroblasts and targeted therapy for cancer. Signal Transduct Target Ther 2021;6:218. [Crossref] [PubMed]
40. Chakravarthy A, Khan L, Bensler NP, et al. TGF-β-associated extracellular matrix genes link cancer-associated fibroblasts to immune evasion and immunotherapy failure. Nat Commun 2018;9:4692. [Crossref] [PubMed]
41. Yi M, Zheng X, Niu M, et al. Combination strategies with PD-1/PD-L1 blockade: current advances and future directions. Mol Cancer 2022;21:28. [Crossref] [PubMed]
42. Wu M, Huang Q, Xie Y, et al. Improvement of the anticancer efficacy of PD-1/PD-L1 blockade via combination therapy and PD-L1 regulation. J Hematol Oncol 2022;15:24. [Crossref] [PubMed]
43. Farhangnia P, Ghomi SM, Akbarpour M, et al. Bispecific antibodies targeting CTLA-4: game-changer troopers in cancer immunotherapy. Front Immunol 2023;14:1155778. [Crossref] [PubMed]
44. Janjigian YY, Shitara K, Moehler M, et al. First-line nivolumab plus chemotherapy versus chemotherapy alone for advanced gastric, gastro-oesophageal junction, and oesophageal adenocarcinoma (CheckMate 649): a randomised, open-label, phase 3 trial. Lancet 2021;398:27-40. [Crossref] [PubMed]
45. Shitara K, Ajani JA, Moehler M, et al. Nivolumab plus chemotherapy or ipilimumab in gastro-oesophageal cancer. Nature 2022;603:942-8. [Crossref] [PubMed]
46. Qiu J, Li H, Tang L, et al. Development and validation of a herb-related gene signature for prognosis prediction and therapeutic response assessment in breast cancer. Transl Cancer Res 2026;15:48. [Crossref] [PubMed]
47. Raskin L, Fullen DR, Giordano TJ, et al. Transcriptome profiling identifies HMGA2 as a biomarker of melanoma progression and prognosis. J Invest Dermatol 2013;133:2585-92. [Crossref] [PubMed]
48. Koeppel F, Bobard A, Lefebvre C, et al. Added Value of Whole-Exome and Transcriptome Sequencing for Clinical Molecular Screenings of Advanced Cancer Patients With Solid Tumors. Cancer J 2018;24:153-62. [Crossref] [PubMed]
49. Yoshihara K, Shahmoradgoli M, Martínez E, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun 2013;4:2612. [Crossref] [PubMed]
50. Lou X, Zhao K, Xu J, et al. CCL8 as a promising prognostic factor in diffuse large B-cell lymphoma via M2 macrophage interactions: A bioinformatic analysis of the tumor microenvironment. Front Immunol 2022;13:950213. [Crossref] [PubMed]
51. Yue Y, Zhang Q, Sun Z. CX3CR1 Acts as a Protective Biomarker in the Tumor Microenvironment of Colorectal Cancer. Front Immunol 2021;12:758040. [Crossref] [PubMed]
52. Xu Z, Zheng T, Zheng Z, et al. TAGAP expression influences CD4+ T cell differentiation, immune infiltration, and cytotoxicity in LUAD through the STAT pathway: implications for immunotherapy. Front Immunol 2023;14:1224340. [Crossref] [PubMed]
53. Silagi ES, Shapiro IM, Risbud MV. Glycosaminoglycan synthesis in the nucleus pulposus: Dysregulation and the pathogenesis of disc degeneration. Matrix Biol 2018;71-72:368-79. [Crossref] [PubMed]
54. Hou XM, Zhang T, Da Z, et al. CHPF promotes lung adenocarcinoma proliferation and anti-apoptosis via the MAPK pathway. Pathol Res Pract 2019;215:988-94. [Crossref] [PubMed]
55. Sun W, Zhao F, Xu Y, et al. Chondroitin polymerizing factor (CHPF) promotes development of malignant melanoma through regulation of CDK1. Cell Death Dis 2020;11:496. [Crossref] [PubMed]
56. Duan X, Yang J, Jiang B, et al. Identification of chondroitin polymerizing factor (CHPF) as tumor promotor in cholangiocarcinoma through regulating cell proliferation, cell apoptosis and cell migration. Cell Cycle 2021;20:591-602. [Crossref] [PubMed]
57. Cui J, Jin S, Jin C, et al. Syndecan-1 regulates extracellular matrix expression in keloid fibroblasts via TGF-β1/Smad and MAPK signaling pathways. Life Sci 2020;254:117326. [Crossref] [PubMed]