A Novel Peroxisome-Related Gene Signature Predicts Breast Cancer Prognosis and Correlates with T Cell Suppression

Introduction

Breast cancer (BC) is the leading cause of cancer-related deaths among women worldwide, and its diversity poses significant challenges for diagnosis and treatment.¹ Despite advancements, clinical outcomes and therapeutic responses remain unpredictable owing to the heterogeneity of the disease.² Gene expression profiling has categorized BC into four molecular subtypes: luminal A, luminal B, HER2-positive, and triple-negative breast cancer (TNBC), thereby improving the outcomes for certain groups.^3–5 However, this classification addresses only part of the heterogeneity, necessitating the exploration of additional dimensions to better understand BC and therapeutic resistance.

Metabolic reprogramming plays a crucial role in BC progression and heterogeneity, involving enhanced glycolysis, TCA cycle, pentose phosphate pathway, glutamate metabolism, and fatty acid metabolism.⁶ Key signaling pathways such as PI3K/AKT and AMPK, along with transcription factors such as c-MYC, p53, and HIF, drive these metabolic alterations.^7,8 These shifts could reveal new therapeutic targets, highlighting the need for more comprehensive studies on the metabolic aspects of breast cancer.

Although mitochondria are well known for their role in cancer cell metabolism, the metabolic functions of other organelles, such as peroxisomes, are poorly understood.⁹ Peroxisomes, which are involved in > 50 enzymatic activities crucial for lipid metabolism, play a significant role in cancer.¹⁰ Dysregulation of peroxisomal lipid metabolism, which is critical for cell membrane production, energy storage, and signaling, is linked to tumor progression.¹¹ Fatty acid oxidation within peroxisomes provides ATP and reduces oxidative stress in cancer cells, making it vital for tumor growth.¹² Enzymes involved in peroxisomal lipid processing are upregulated in various cancers, including prostate, colorectal, breast, liver, ovarian, brain, pancreatic, and bladder cancers.^13–24 As such, peroxisome degradation, targeting these peroxisome-related prognostic biomarkers, or chemically inhibiting peroxisomal lipid processing can significantly inhibit tumor growth.^24–26 There have not been any direct investigations into the role of peroxisome in breast cancer. Understanding the metabolic functions of peroxisomes in breast cancer could lead to novel therapeutic strategies targeting metabolic flexibility and resistance in breast cancer.

Methods and Materials

Peroxisome-Related Genes

Three sources of peroxisome-related genes (PRGs), previously utilized in prior studies, were identified.^27–29 These include the PeroxisomeDB 2.0 database (100 genes), the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (83 genes), and human liver peroxisomes (60 genes).^30–32 After eliminating redundant entries, a total of 113 unique PRGs were identified (Table S1).

Datasets and Processing

The gene expression profiles (FPKM values) of 1182 breast cancer patients (normal, n=113; tumor samples, n=1069) from The Cancer Genome Atlas (TCGA) breast adenocarcinoma (BRCA) cohort, generated using the Illumina HiSeq RNA-Seq platform, were downloaded from the TCGA Data portal (https://portal.gdc.cancer.gov/). Clinical data including patient age, sex, molecular subtype, cancer stage, and survival were obtained from the UCSC Xena database (https://xenabrowser.net/). For validation, a cohort comprising 327 samples from the Gene Expression Omnibus (GEO) GSE86166 database was used. Normalized expression data and clinical information for the GSE86166 dataset were downloaded from the GEO repository. Both datasets underwent normalization and were compared for common gene expression using the “limma” and “sva” R packages. Batch effects were addressed using the combat function of the Sva R package. A total of 111 PRGs common to both datasets were identified and used for downstream analysis. External validation was performed by merging four GEO transcriptomic datasets (GSE16446, n=120; GSE20685, n=327; GSE42568, n=121; and GSE58812, n=107) of patients with breast cancer. All the four datasets were normalized with the “limma” and batch effects were removed with “sva” R package. Survival data were provided for 640 patients with breast cancer; hence, only these patients were included in the analysis.

Single-Cell Analyses

The Tumor Immune Single-cell Hub (TISCH, accessed on 12 February 2024 at http://tisch.comp-genomics.org/) is a dedicated database for evaluating the tumor microenvironment (TME), providing cell-type annotations at the single-cell level.³³ TISCH utilizes a standard pipeline in MAESTRO to preprocess each single-cell dataset.³⁴ In this study, we investigated gene expression levels of the PRG risk signature in BRCA_GSE176078 (n = 26, single cells = 89,471), a single-cell dataset of primary breast cancer. The data were preprocessed and quality control steps were skipped. The “Seurat” R package was used to normalize the data, identify variable genes, and perform principal component analysis. The “NormalizeData” function was applied with the “LogNormalize” method and a scale factor of 10,000. The “FindVariableFeatures” function was employed to identify the top 2000 variable features using the “vst” method. These features were then used for principal component analysis through the “ScaleData” and “RunPCA” functions. Clusters and cell annotations were retained from the metadata files. For dimensionality reduction, we used the “RunTSNE” function, and visualized the results using t-SNE plots. Additionally, individual gene expression analyses of PRG risk genes were conducted on five breast cancer single-cell datasets using the “Gene” module of the TISCH website. The four additional breast cancer datasets included BRCA_GSE161529 (n = 52, single cells = 332,168), BRCA_GSE148673 (n = 6, single cells = 10,359), and BRCA_EMTAB8107 (n = 14, single cells = 33,043), and one metastatic dataset, BRCA_GSE150660 (n = 3, single cells = 10,605).

Generation of Risk-Gene Signature

A risk signature was generated by performing univariate Cox proportional hazards regression analysis using the R package “survival” to identify PRGs with significant prognostic impact in the training cohort. Significance was considered at p <0.05 (Wald test). The prognostic genes were subjected to a multivariate Cox proportional hazards regression model (R package “survival”) to obtain an optimal risk signature.³⁵ The risk score for each patient in both cohorts was calculated by multiplying the sum of the Cox regression coefficients for each signature gene with its corresponding mRNA expression value. The median risk score was adopted as the cut-off to categorize patients as high- or low-risk. We conducted a Kaplan-Meier analysis to compare the overall survival among the different risk groups. The best expression cutoff value was identified using the “surv_cutpoint” function from the “survminer” R package. The predictive ability of this signature was assessed using time-dependent ROC (receiver operating characteristic (ROC) and AUC (area under the curve) analyses. The “timeROC” R package was used to conduct 1-year, 5-year, and 10-year ROC analyses. This process was repeated for internal validation of the GEO cohort.

Independent Prognostic Value Assessment and Nomogram Construction

The independent prognostic value of the PRG risk score was assessed using univariate and multivariate Cox regression analyses, along with other clinical features such as age, molecular subtype, and cancer stage. Significant variables identified in the Cox regression analysis were used to construct a nomogram with the R packages “rms” and “survival”. Calibration plots were generated to evaluate the concordance between actual and predicted survival outcomes.

Genomic Aberrations

To calculate arm- and focal-level somatic copy number alterations (SCNAs) in tumor samples, we used the GISTIC2.0 (Genomic Identification of Significant Targets in Cancer) tool available on GenePattern (https://cloud.genepattern.org).³⁶ SNP6 files downloaded from the Genomic Data Commons (GDC) data portal (https://portal.gdc.cancer.gov/) were used as input files. We obtained Simple Nucleotide Variation data for the TCGA BRCA cohort from TCGA Data Portal (https://portal.gdc.cancer.gov/). To analyze gene mutations and create oncoplots, we used the R package “maftools”.

Tumor Microenvironment (TME) Annotations

To infer tumor purity and the admixture of immune and stromal cells from cancer tissue gene expression data, we utilized the “ESTIMATE” R package.³⁷ This package integrates gene expression data to estimate tumor purity as well as the levels of immune and stromal components in cancer tissues. To characterize the primary cellular constituents of the TME, we used the MCP-counter algorithm. This robust algorithm calculates the absolute abundance of ten cell types, including eight distinct immune cell types and two stromal cell types, by analyzing transcriptomic data from heterogeneous tissues.³⁸ Finally, for a detailed quantitative analysis of the relative abundance of 22 immune cell types within TCGA cohort, we employed the CIBERSORT algorithm.³⁹ CIBERSORT uses gene expression profiles from bulk cancer tissues to estimate cellular fractions, covering various subtypes of major cell lineages and providing a comprehensive view of the immune landscape.

Functional Enrichment Analysis

To explore the enrichment of hallmark cancer pathways across different risk groups, we employed the “clusterProfiler” package for gene set enrichment analysis (GSEA).⁴⁰ Additionally, the “gsva” package in R was used to estimate the enrichment of the KEGG pathways and Gene Ontology (GO) terms between the identified risk subgroups.⁴¹

Public Database Analysis

The expression levels of individual PRG risk genes across various molecular subtypes of breast cancer (BRCA) tissues and compared with normal breast tissues from the GTEx database were analyzed using the GEPIA 2 platform (http://gepia2.cancer-pku.cn/).⁴² The expression levels of individual PRG risk genes across breast cancer stages were assessed and compared using a Gene Set Cancer Analysis platform (http://bioinfo.life.hust.edu.cn/GSCA/).⁴³

Protein Expression Analysis

The Human Protein Atlas (HPA) was used to explore the expression of NMNAT2, QPRT, and NT5E in normal and breast cancer tissues.⁴⁴

Biomarker Ability to Predict Immunotherapy Response

To assess the response to immunotherapy, we employed a Tumor Immune Dysfunction and Exclusion (TIDE) algorithm (http://tide.dfci.harvard.edu).⁴⁵ Additionally, we compared the biomarker capability of our PRG risk signature with that of other established immunotherapy response biomarkers, such as MSI score, TMB, CD274, CD8, IFNG, T Clonality, B Clonality, and Merck18, using the Biomarker Evaluation module available on the TIDE website. We validated the immunotherapy response by examining the activity of the PRG risk signature in the IMvigor210 cohort, which comprised 348 patients with urothelial carcinoma treated with anti-PD-L1 immune checkpoint inhibitors. The participants in the IMvigor210 cohort were stratified into high- or low-risk groups based on the PRG risk signature. This stratification was achieved by multiplying the mRNA expression levels of the PRG risk signature genes by their respective coefficients derived from the TCGA data. To determine the optimal cutoff point for expression values, we utilized the “surv_cutpoint” function from the “survminer” R package, ensuring the precise categorization of patients into risk groups.

Drug Sensitivity Analysis

We conducted drug sensitivity analysis of 198 small molecules to determine their potential effectiveness in patients with breast cancer. Using the “oncoPredict” R package,⁴⁶ we calculated the half-maximal inhibitory concentration (IC50) for each drug. A drug was considered effective when the IC50 value was significantly lower in the high-risk group than in the low-risk group.

Cell Lines and Cell Culture

The human breast cancer cell lines MCF-7, MDA-MB-231, T47D, and BT549, along with the normal breast cell line MCF-10A and HMEC, were obtained from the Committee of Type Culture Collection of the Chinese Academy of Sciences (Shanghai, China). Cells were cultured in DMEM supplemented with 10% fetal bovine serum (FBS), 100 U/mL penicillin, and 100 mg/mL streptomycin. The cells were maintained in a humidified incubator at 37°C and 5% CO2. We regularly authenticated all the cell lines used in this investigation by assessing their morphology and conducting tests to ensure the absence of mycoplasma contamination.

Quantitative Real-Time PCR

Total RNA was extracted and purified using Trizol Reagent (Takara, Otsu, Japan). RNA was reverse-transcribed to synthesize cDNA. Quantitative real-time PCR (qRT-PCR) was performed using the SYBR Green PCR Kit (Takara, Otsu, Japan). The mRNA expression levels were normalized to the internal control, GAPDH, and the relative mRNA levels in the treated group were compared to those in the control group. Primers for quantitative PCR were designed using Primer Premier 5.0 and Beacon Designer 7.8 software, and synthesized by General Biosystems Co., Ltd. The primer sequences used in this study are listed in Table S2.

Statistical Analysis

Gene expression and risk score comparisons between two or more groups were conducted using Wilcoxon or Kruskal–Wallis tests. Categorical variables were compared using the chi-square test. qPCR results were analyzed using the Student’s t-test. Survival analysis was performed using the Kaplan-Meier method with Log rank test. Univariate and multivariate factor analyses were performed using the Cox regression hazard models. Correlations were assessed using Spearman’s or Pearson’s correlation test. All statistical analyses were conducted using R software version 4.3.1 (http://www.r-project.org).

Results

Overview of the Analytical Pipeline Used in This Study

The data processing and methodology used in this study are outlined in Figure 1. The analysis is divided into five main components: 1) construction of the peroxisome-related genes (PRG) risk signature, 2) comprehensive evaluation of clinical, functional, immunological, and mutational features linked to the PRG risk signature, 3) individual analysis of PRG genes, 4) development of a monogram and external validation, and 5) assessment of therapeutic response.

Figure 1 Study workflow summary. A peroxisome-related gene (PRG) signature was developed using univariate and multivariate Cox regression analyses on transcriptomic data from 1069 breast cancer (BRCA) patients in TCGA (training set) and validated in 327 tumor samples from the GEO database. The resulting 7-gene PRG risk signature enabled stratification of BRCA patients into high- and low-risk groups. Clinical, functional, immunological characteristics, and mutational landscapes were analyzed. Expression of each PRG was examined via RT-qPCR, single-cell analysis, and ProteinAtlas protein-level data. A nomogram was constructed, and the PRG signature was further validated across four external GEO datasets. Finally, the signature’s predictive value for immunotherapy and chemotherapy responses was assessed.

Identification of Peroxisome-Related Prognostic Gene Signature in Breast Cancer

To identify the prognostic gene signature comprising PRGs, we performed univariate Cox regression analysis of 111 PRGs, which indicated that 12 PRGs had a significant impact on the prognosis of breast cancer patients in the TCGA training cohort (Figure 2A). Among the 12 PRGs, 4 played a risky role with a hazard ratio greater than one (HR>1), while overexpression of the remaining 8 PRGs was protective for the prognosis of breast cancer patients. These 12 PRGs were further subjected to multivariate Cox regression analysis, which revealed 7 PRGs to be considered for the development of the risk signature (Figure 2B). A risk score was obtained for each patient by summing the mRNA expression values of each of the seven PRGs multiplied by their corresponding multi-Cox regression coefficients. The risk plot depicting the risk score for each individual in the training and validation datasets, along with the expression patterns of these seven PRGs, is illustrated in Figure 2C and D. Patients with high-risk scores demonstrated poor survival in both cohorts (Figure 2E and F). The predictive performance of the risk signature was assessed using a receiver operating characteristic (ROC) curve analysis. In TCGA BRCA cohort, the analysis showed an AUC of 0.649 at 1 year and 0.681 at 10 years (Figure 2G). A similar trend was observed in the validation cohort, with AUC of 0.651 and 0.573 at 1 and 10 years, respectively (Figure 2H).

Figure 2 Construction of PRG risk signature. (A) Univariate Cox regression analysis of peroxisome-related genes (PRGs) in the TCGA BRCA dataset. (B) Bar plots coefficient of multivariate-Cox regression analysis. (C) Risk plot and heatmap of PRGs in TCGA BRCA samples and (D) GEO (GSE86166) validation dataset. (E) Kaplan-Meier curves of survival analysis of TCGA BRCA samples, and (F) GEO validation dataset. (G) Recursive operating curve (ROC) analysis of TCGA BRCA samples and (H) the GEO validation dataset.

The PRG Risk Signature Primarily Affects the Prognosis of Luminal Subtypes

To ascertain the clinical implications of the PRG risk signature, we evaluated differences in the clinical features of patients with high- and low-risk breast cancer. The heatmap illustrates the pattern of clinical features of the patients and indicates a significant difference between the molecular subtype (P < 0.001) and cancer stage (P < 0.05) of the high- and low-risk patients (Figure 3A). Compared to the low-risk subgroup, the high-risk subgroup comprised a higher number of basal and luminal B subtypes and a lower number of luminal A subtypes (Figure 3B). Basal and luminal B subtype patients had a higher risk score than patients with the luminal A subtype (Figure 3C). Similarly, the number of patients with stage II and IV cancer was slightly higher and the number of patients with stage I cancer was lower in the high-risk subgroup (Figure 3D). The risk score progressively increased with cancer stage (Figure 3E). The risk signature only affected the prognosis of the luminal breast cancer subtypes (Figure 3F-J). High-risk patients with luminal A and B subtypes showed a poor prognosis compared to low-risk patients (Figure 3H and I). This difference in prognosis between the subgroups was not observed in the basal, Her2, and normal subtypes (Figure 3F-J). Consequently, the ROC curves indicated higher AUC values (1-year: 0.816 vs 0.484; 3-year: 0.729 vs 0.542; 5-year: 0.710 vs 0.582; 10-year: 0.720 vs 0.634) for the luminal subgroup than for the other subtypes (Figure 3K and L). However, the prognostic impact of the risk signature was evident in both cancer stage categories (stages I–II and III–IV) (Figure 3M and N).

Figure 3 Clinical evaluation of PRG risk signature. (A) Heatmap of distribution of clinical features among PRG risk subgroups. Chi-square test: Significance levels represented as: *P < 0.05; ***P < 0.001. (B) Bar plot depicting the distribution of BRCA subtypes and (D) cancer stage in each risk subgroup. (C) Box plots comparing the PRG riskScore among BRCA subtypes and (E) cancer stages in each risk subgroup. (F) Kaplan-Meier curves comparing overall survival between risk subgroups in Basal, (G) Her2, (H) Luminal A, (I) Luminal B and (J) Normal subtypes. P value < 0.05 was considered significant. (K) Recursive operating curve (ROC) analysis in TCGA BRCA luminal and (L) other subtypes. (M) Kaplan-Meier curves comparing overall survival between risk subgroups in stage I–II and (N) stage III–IV. P value < 0.05 was considered significant.

Functional and TME Characteristics of PRG Risk Signature Indicate Reduced Immune Response

Next, we evaluated the functional and immunological implications of PRG risk signature in breast cancer. Enrichment of hallmark cancer pathways indicated that the high-risk subgroup had a high proportion of cancer cells, as the expression of E2F targets, G2M checkpoints, and mitotic spindles were elevated (Figure 4A). In contrast, the low-risk subgroup experienced heightened inflammatory and immune response activity with the activation of pathways such as the interferon alpha/gamma response and TNFA signaling (Figure 4B). A similar outlook was evident in the enrichment analysis of Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways (Figure 4C and D). The inflammatory and immune response pathways were significantly upregulated in the low-risk subgroup, whereas terms and pathways associated with the cell cycle and cell division were enriched in the high-risk subgroup (Figure 4C and D). The ESTIMATE algorithm was used to confirm the immune activity level between risk subgroups by investigating the immune and stromal TME scores. The results indicated a high immune score for the low-risk subgroup, whereas there was no difference in the stromal score between the risk subgroups (Figure 4E). The MCP-counter algorithm further established that T-cell subsets were deficient in the high-risk subgroup (Figure 4F). Moreover, the CIBERSORT results showed a higher infiltration of CD8+ and regulatory T cells in the low-risk subgroup and highlighted differences in the infiltration of monocytic lineage cells between the subgroups (Figure 4G). The high-risk subgroup exhibited increased macrophage polarization toward the M2 phenotype, which is anti-inflammatory. Overall, these results indicate that peroxisome dysfunction in breast cancer patients results in the suppression of the immune response, resulting in poor prognosis.

Figure 4 Functional implications of PRG risk signature. (A) GSEA (Gene set enrichment analysis) plots of hallmark cancer pathways in high- and (B) low-risk subgroups. (C) Heatmap of gene set variation analysis (GSVA) of KEGG and (D) GO pathways in risk subgroups. (E) Violin plot of TME scores. (F) Heatmap showing infiltration of immune and stromal cells (MCP-counter) in BRCA samples. (G) Violin plot of difference in immune cell infiltration (CIBERSORT) between the risk subgroups. Significance levels represented as: *P < 0.05; **P < 0.01; ***P < 0.001.

Genomic Aberrations are Predominant in PRG High-Risk Subgroup

Previous studies have demonstrated that a high immune score is inversely correlated with genomic aberrations.⁴⁷ We investigated the mutation landscape of PRG risk subgroups and found that the high-risk subgroup, which demonstrated lower T-cell infiltration, had higher rates of focal synonymous copy number alterations (SCNA) than the low-risk subgroup (Figure 5A). A high frequency of arm-level amplifications and deletions was observed in the high-risk subgroup (Figure 5B). The high-risk subgroup showed frequent deletions at several arm levels compared with the low-risk subgroup. Mutations in BRCA-specific genes were also frequent in the high-risk subgroup (83.16% vs 78.53%) (Figure 5C and D). TP53 had 39% mutations in the high-risk subgroup compared to 25% in the low-risk subgroup. Mutations in PIK3CA showed an opposite trend (high-risk subgroup: 30% vs low-risk subgroup: 37%). Moreover, the high-risk subgroup demonstrated a high tumor mutation burden and a positive correlation with the PRG risk score (Figure 5E and F).

Figure 5 Mutational landscape of PRG risk subgroups. (A) Depiction of focal-level somatic copy number alteration (SCNA) between the PRG risk subgroups. (B) Arm-level amplification and deletion frequencies between the risk subgroups. (C) Oncoplot depicting the mutation frequency of top 10 mutated genes in the high- and (D) low-risk groups. (E) Box plot showing the comparison of tumor mutation burden (TMB) between the risk subgroups. (F) Pearson’s correlation between riskScore and TMB.

Individual Clinical and Immunological Evaluation of PRG Risk Signature

Understanding the PRG risk signature and investigating the functional profile of each PRG risk gene are crucial. Notable variations in individual gene expression were observed across subtypes. Specifically, ACBD5, PEX10, and SLC27A2 exhibited higher expression in breast cancer tissues (TCGA BRCA) than in normal tissues (GTEx) and was particularly higher in luminal than basal and HER2 subtypes (Figure 6A). In contrast, PEX3 and ACSL5 consistently showed lower expression levels in cancer tissues across all the subtypes. NOS2 expression was slightly elevated in the basal subtypes compared with that in the other subtypes. Among these genes, only ACSL5 and SLC27A2 demonstrated significant variations in expression with cancer stage (Figure 6B). A persistent decrease in expression from stage I to IV was observed for ACSL5, PEX3, SLC27A2, and PEX10, although PEX10 expression trend reversed to an increased after stage III (Figure 6C). Conversely, NOS2 and ACBD5 displayed a continuous upward trend from stage I to IV. Immunologically, ACSL5 exhibited a significant positive correlation with M1 macrophages and CD8 T cells, suggesting its involvement in proinflammatory and immune response functions (Figure 6D). CD8 + T cell infiltration was negatively correlated with ACBD5 and NOS2, indicating their potential roles in immune response suppression. An intriguing opposing functional relationship was observed between ACSL5 and NOS2 in macrophages and CD8 + T cells, highlighting their contrasting roles in the immune system.

Figure 6 Evaluation of individual PRG risk gene in breast cancer. (A) Box plots depicting the comparison of mRNA expression of individual PRG risk genes in TCGA breast cancer tissues to GTEx normal brain tissues at whole and across BRCA subtypes. (B) Box plots of comparison of mRNA expression of individual PRG risk genes across BRCA stages. (C) Trend plot depicting the mRNA expression pattern of individual PRG risk genes from stage I to IV. (D) Spearman’s correlation between individual PRG risk genes and fraction of infiltrated immune cells (CIBERSORT). Spearman correlation test; Statistical significance is indicated as follows: *P<0.05; **P < 0.01; ***P < 0.001.

PRG Risk Signature at Single Cell Resolution

To confirm the expression patterns of PRG risk genes at the cellular level, we performed single-cell analysis of the BRCA dataset (GSE176078) obtained from the TISCH database. The PRG risk score was visualized using t-SNE and violin plots, revealing slightly higher expression in malignant and stromal cells than in immune cells (Figure 7A-E). Certain immune cell subsets, such as monocytes/macrophages, dendritic cells, and proliferative T cells (Tprolif), also exhibited higher expression levels than the other immune cells (Figure 7C). Consistent with the BRCA bulk RNA subtype analysis, TNBC and HER2 subtypes showed elevated risk scores compared to ER+ samples (Figure 3C; Figure 7F and G). Specifically, ACBD5, PEX3, and PEX10 demonstrated enhanced expression in malignant cells, whereas ACSL5 showed increased expression in immune cells, such as monocytes/macrophages, dendritic cells, and Tprolif (Figure 7H and I). Variations in expression were also evident among the BRCA subtypes (Figure 7J). These findings corroborate the observations from the bulk RNA-seq analysis. We further validated these expression profiles using other BRCA single-cell datasets, confirming a consistent pattern for ACSL5, PEX3, and PEX10 (Figure 7K).

Figure 7 Visualization of the PRG risk signature at single-cell resolution. (A) t-SNE plot showing the main cell types in the breast cancer single-cell dataset (GSE176078). (B) t-SNE plot displaying the expression pattern of the PRG risk signature in the same dataset. (C) Violin plot illustrating the expression of the PRG risk signature across each cell type. (D) t-SNE plot demonstrating cellular lineage distribution in the dataset. (E) Violin plot depicting the expression pattern of the PRG risk signature in each cellular lineage. (F) t-SNE plot showing the distribution of cancer subtypes in the dataset. (G) Violin plot illustrating the expression of the PRG risk signature across cancer subtypes. (H) Bubble plot displaying the expression pattern of individual PRG risk genes across cellular lineages. (I) Bubble plot showing the expression pattern of individual PRG risk genes across cell types. (J) Bubble plot demonstrating the expression pattern of individual PRG risk genes across cancer subtypes. (K) Expression analysis of ACSL5, PEX3, and PEX10 across cell lineages and cell types in five breast cancer single-cell datasets using the “Gene” module from the TISCH database. Abbreviations: CD4Tcov, Conventional CD4 T cells; Tprolif, Proliferating T cells; CD8Tex, Exhausted CD8 T cells; DC, Dendritic cells; Mono/Macro, Monocytes/Macrophages; SMC, Smooth muscle cells; ER+, Estrogen receptor positive; HER2+, Human epidermal growth receptor 2 positive; TNBC, Triple negative breast cancer; NA, Not available.

Validation of PRG Risk Signature

To validate the mRNA expression of the PRG risk signature, we compared the expression levels of the PRG risk genes in normal breast cells (MCF-10A and HMEC) and breast cancer cells (MCF-7, MDA-MB-231, T47D, and BT549) (Figure 8). Overall, all the PRGs showed elevated expression in breast cancer cells, in both TNBC (MDA-MB-231 and BT549) and ER+ cells (MCF-7 and T47D), in comparison to both normal breast cells (MCF-10A and HMEC) (Figure 8). NOS2 was excluded from qPCR analysis because it is inducible and requires stimulation with IFN-γ or other inflammatory stimuli.⁴⁸

Figure 8 Validation of PRG risk genes in breast cancer. Qualitative PCR results showing mRNA expression level of PRG risk genes in normal breast cell (MCF-10A and HMEC) and breast cancer cells (MCF-7, T47D, BT549 and MDA-MB-231). The data represent the mean ± SEM (standard error of mean) of n=3 independent experiments (independent biological replicas) for each condition. Two-tailed unpaired T-test; *P<0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001. Abbreviation: ns, not significant.

Additionally, the protein-level expression of PRG risk genes in breast cancer and normal breast tissues was visualized using immunohistochemical data from The Human Protein Atlas database (https://www.proteinatlas.org/) (Figure 9A). In normal breast tissues, PRG risk genes were mainly observed in the glandular and myoepithelial cells (Figure 9B). Elevated levels of ACBD5, PEX3, DAO, and ACSL5 were observed in breast cancer tissues (Figure 9C).

Figure 9 Protein-level expression of PRG risk genes in breast cancer. (A) Representative Immunohistochemistry (IHC) images of individual PRG risk genes in normal and cancerous breast cancer tissues from HPA database. (B) Expression level of PRG risk genes in normal breast tissues at cellular level. (C) Expression level of PRG risk genes in breast cancer tissues.

Construction of PRG Nomogram

To enhance the clinical utility of the PRG risk signature, a nomogram was developed that integrates the PRG risk signature with various clinical characteristics of breast cancer patients, such as age, molecular subtype, and cancer stage. Comprehensive clinicopathological data were available for only 938 patients; therefore, these patients were included exclusively in the Cox regression analysis. Age, molecular subtype, cancer stage, and PRG risk score all exhibited independent prognostic value in both univariate and multivariate Cox regression analyses (Figure 10A). These variables were selected for inclusion in a nomogram, where each was assigned specific points according to its impact on survival risk (Figure 10B). For instance, a patient with a luminal A subtype, stage I cancer, age over 60 years, and a PRG risk score just than 0.3 would accumulate a total score of 75.4 points. This score translated into survival probabilities of 99.3%, 96.4%, 92.8%, and 81.7% at 1, 3, 5, and 10 years, respectively. The nomogram score revealed significant differences in overall survival between the low- and high-risk subgroups (Figure 10C). Its predictive capability was demonstrated using a calibration plot, which showed close alignment between the predicted and observed OS values at 1, 3, 5, and 10 years (Figure 10D). Additionally, the newly constructed nomogram exhibited a high predictive accuracy for OS, with area under the curve (AUC) values of 0.846, 0.785, 0.739, and 0.739 for 1-, 3-, 5-, and 10-year survival, respectively (Figure 10E).

Figure 10 Independent prognostic value of PRG risk signature and PRG nomogram. (A) Univariate and multivariate Cox-regression analysis. (B) A nomogram integrating PRG risk signature and clinical features to predict survival probability of the TCGA BRCA patients. (C) Kaplan-Meier curve depicting the OS between nomogram-stratified risk subgroups. (D) Calibration curves demonstrating the nomogram’s accuracy in predicting survival probabilities compared to observed survival outcomes. (E) Recursive operating curve (ROC) analysis of nomogram in TCGA BRCA samples.

External Validation

Four external independent breast cancer cohorts (GSE16446, n=120; GSE20685, n=327; GSE42568, n=121; and GSE58812, n=107) with complete transcriptomic and survival data were obtained to assess the robustness, general applicability, and limitations of the PRG risk signature in breast cancer. The four datasets were merged, batch effects were removed, and 640 patients with breast cancer were included (Figure 11A and B). The risk score was calculated and the patients were categorized into low- and high-risk subgroups by employing individual gene expression data and multi-Cox regression coefficients of each PRG risk gene. Figure 11C illustrates the plot of the risk score and survival outcomes, and the expression pattern of PRGs in the merged dataset. The expression pattern indicated an outlook similar to that observed in the training and internal validation datasets. The three-dimensional PCA (3dPCA) plot illustrated well-separated distinct spatial dimensions for the risk subgroups (Figure 11D). There was a significant survival difference (P < 0.051) between the high-risk subgroup (n=241) and low-risk subgroup (n=399) (Figure 11E).

Figure 11 External validation of PRG risk signature. (A) PCA plots illustrating the Principal Component Analysis of 4 GEO datasets before and (B) batch correction. (C) Risk plot and heatmap of PRGs in the merged dataset. (D) 3-dimensional Principal Component Analysis of risk subgroups in the merged dataset. (E) Kaplan-Meier survival curves for the merged dataset.

Implications for Immuno- and Molecular-Targeted Therapy

Tumor Immune Dysfunction and Exclusion (TIDE) models immune evasion in cancers by assessing T-cell dysfunction (high CTL infiltration) and exclusion (low CTL infiltration) and combining them into a TIDE score. Higher scores indicate greater immune evasion and lower response to immune checkpoint inhibitors (ICIs). In our study, high-risk patients with breast cancer had lower TIDE scores, suggesting a better response to immunotherapy (Figure 12A). We validated this using the IMvigor210 cohort (urothelial carcinoma patients on pembrolizumab), categorizing them into high- and low-risk subgroups based on the PRG risk signature activity. The high-risk subgroup showed improved survival rates, indicating an increased sensitivity to immunotherapy (Figure 12B). UC patients with an excluded immunophenotype had higher risk scores than those with a desert immunophenotype, suggesting that the PRG risk signature reflects CD8 T cell accumulation but is unable to infiltrate efficiently and remain in the periphery (Figure 12C). Overall, the PRG risk signature demonstrated comparable biomarker efficacy to other well-known immunotherapy response markers, such as MSI score, TMB, CD274, CD8, IFNG, T Clonality, B Clonality, and Merck18, as assessed using the Biomarker Evaluation module on the TIDE website (Figure 12D).

Figure 12 Biomarker efficacy of PRG risk signature. (A) Association between the PRG riskScore and TIDE score in TCGA BRCA patients. (B) Kaplan-Meier survival curves showing differences between PRG risk subgroups. (C) Box plot comparing immunophenotypes (desert, excluded, inflamed) based on PRG riskScore. (D) Bar plot displaying AUC values for predicting immunotherapy response across various pan-cancer immunotherapy cohorts for the PRG risk signature and other known immunotherapy biomarkers, along with the percentage of PRG risk signature genes mapping in the corresponding cohorts. (E) Bar plot depicting drug sensitivity analysis results.

Drug sensitivity analysis indicated resistance to the majority of 198 small anticancer molecules, including several chemotherapeutic drugs (Figure 12E). High risk was identified as sensitivity to BCL2 family protein inhibitors, such as WEHI-539 and Navitoclax (selective inhibitor of Bcl-XL), ABT-737 (Bcl-2 and Bcl-xL inhibitor), and UMI-77 (Mcl-1 inhibitor). In addition, several protein kinase inhibitors, such as polo-like kinase 1 (PLK1: BI-2536), Bruton’s tyrosine kinase (BTK: Ibrutinib), cyclin-dependent kinase 1 (CDK1: RO-330), epidermal growth factor receptor (EGFR: Erlotinib), and p38 MAPK (mitogen-activated protein kinase: BIRB 796), were also identified.

Discussion

Over the years, targeting cancer metabolism has been the cornerstone of anticancer therapeutic strategies.⁴⁹ However, a significant challenge in this approach is the metabolic plasticity of the cancer cells. During tumor initiation, progression, and response to conventional treatments, cancer cells undergo extensive metabolic rewiring, allowing them to activate alternative pathways and resist targeted therapies.¹² Therefore, more effective anti-metabolic cancer therapy may involve targeting multiple metabolic pathways simultaneously and employing personalized treatment strategies. In this study, we aimed to deepen our understanding of the underexplored aspects of cancer and peroxisome metabolism in breast cancer. We identified a novel seven-gene peroxisome-related risk signature that can independently predict the prognosis of patients with breast cancer. This risk score showed a progressive increase from stage I to stage IV, effectively differentiating the prognoses at each cancer stage. Notably, the Basal, Her2 and Luminal B subtypes had higher risk scores than the Luminal A subtype. A significant difference in prognosis was observed among the luminal subtypes, with a trend toward worse outcomes. To enhance predictive accuracy, we developed a nomogram based on the PRG risk signature. This tool offers improved predictive capability for breast cancer prognosis and provides a valuable resource for personalized treatment planning. An overview of the main aspects of PRG risk difference is summarized in Figure 13.

Figure 13 PRG-Based Risk Stratification of Breast Cancer. Breast cancer patients were classified into high- and low-risk subgroups based on the 7-gene PRG signature (ACBD5, ACSL5, DAO, NOS2, PEX3, PEX10, and SLC27A2). The high-risk subgroup exhibited elevated expression of ACBD5, DAO, NOS2, and PEX3, primarily in tumor cells, and these genes were linked to peroxisome biogenesis (ACBD5 and PEX3) and redox balance (DAO and NOS2). In contrast, the low-risk subgroup showed increased expression of ACSL5, PEX10, and SLC27A2, involved in fatty acid metabolism and primarily expressed by CD8+ T cells. The differential regulation of peroxisome activity between tumor and CD8+ T cells was associated with breast cancer prognosis. The high-risk subgroup demonstrated a higher frequency of genomic aberrations, correlated with immune suppression and reduced CD8+ T cell infiltration. This subgroup was resistant to chemotherapy but showed potential for responding to immunotherapy.

Peroxisomes are essential organelles that play key roles in lipid metabolism, including fatty acid oxidation, ether phospholipids, and bile acid synthesis, as well as in maintaining reactive oxygen species (ROS) homeostasis.⁵⁰ The oxidation of fatty acids within peroxisomes generates ATP, which is crucial for cancer cells to meet their increased energy demands and produce H2O2 as a by-product, contributing to oxidative stress.⁵¹ Peroxisomes act as ROS sinks because they contain certain enzymes such as catalase to prevent oxidative stress, which can contribute to cancer if uncontrolled.⁵² Additionally, the generation of nicotinamide adenine dinucleotide phosphate (NADPH), with its reducing power during catalytic reactions and fatty acid oxidation, forms another critical antioxidant pathway that helps counteract oxidative stress in cancer cells and prevents apoptosis.⁵³ Thus, peroxisomal activity is essential for cancer growth, balancing energy production and ROS management.

The PRG high-risk subgroup exhibited elevated expression of PEX3, ACBD5, DAO, and NOS2. PEX3 and ACBD5 are key components of peroxisome biogenesis and are indicative of upregulated peroxisome formation.^54,55 De novo peroxisome biogenesis commences with PEX3 binding to either PEX16 in the ER or PEX14 in the mitochondrial membrane.⁵⁴ These complexes containing PEX3/PEX16 or PEX3/PEX14 bud off from the organelles in the vesicles. In lymphoma cells, depletion of PEX3 significantly impairs peroxisome biogenesis.⁵⁶ Moreover, reducing PEX3 expression renders lymphoma cells more susceptible to ROS-induced apoptosis. Therefore, targeting PEX3 could serve as a crucial therapeutic approach in breast cancer treatment to impair peroxisomal activity, rendering cancer cells more susceptible to ROS-induced apoptosis. Our drug sensitivity analysis highlighted apoptosis inducers, particularly BCL2 family protein inhibitors (WEHI-539 and Navitoclax [selective inhibitor of Bcl-XL], ABT-737 [Bcl-2 and Bcl-xL inhibitors], and UMI-77 [Mcl-1 inhibitor]), as a key class of therapeutic agents.⁵⁷ These inhibitors promote apoptosis by disrupting the balance between pro- and anti-apoptotic proteins within the cell, thereby facilitating the induction of death in cancer cells with compromised peroxisomal activity.

NOS2, which is primarily elevated in the basal subtype and stage IV breast cancer patients, contributes to nitrosative stress by producing nitric oxide, and is increasingly recognized as a key player in cancer progression.⁵⁸ D-amino acid oxidase (DAO), a flavin adenine dinucleotide (FAD)-dependent oxidase, shows elevated expression in high-risk subgroups and is associated with a poor prognosis. DAO catalyzes the oxidation of neutral and polar d-amino acids with strict stereospecificity, resulting in the production of α-keto acids, ammonia, and hydrogen peroxide.⁵⁹ Therefore, the risk genes in the PRG risk signature are critically involved in or associated with redox homeostasis in breast cancer, highlighting the crucial role of peroxisomal activity in ROS balance in breast cancer prognosis. Additionally, excess peroxisomal ROS/RNS can directly inactivate peroxisomal matrix proteins, which correlates with another characteristic of the PRG high-risk subgroup: lower expression of PEX10. PEX10 is an integral component of the peroxisomal protein import machinery that facilitates the translocation of peroxisomal matrix proteins from the cytosol to the peroxisome.⁶⁰ Its suppression may indicate a dampened import of matrix proteins due to elevated ROS levels. These observations provide a hypothetical basis for in-depth functional investigations of these biomarkers.

Elevated levels of ACBD5 (the acyl-coenzyme A-binding domain protein) also suggest increased peroxisomal activity, especially in the metabolism of very long-chain fatty acids. ACBD5 interacts with VAPB (vesicle-associated membrane protein-associated protein B) in the endoplasmic reticulum, forming membrane contacts that are crucial for lipid metabolism and promoting peroxisome biogenesis.⁵⁵ Despite its importance in metabolism, ACBD5’s role in tumorigenesis has not been extensively studied. However, upregulation and genomic alterations in ACBD proteins (ACBD3/4/5) have been linked to cancer development. ACBD4 expression and polymorphisms predict overall survival in hepatitis B virus-related hepatocellular carcinoma patients after hepatectomy.⁶¹ Similarly, the ACBD5-RET rearrangement has been linked to papillary thyroid cancer via the enhanced phosphorylation of ERK proteins in the MAPK pathway.⁶² ACBD3 has been associated with breast cancer prognosis, showing variations in expression based on the hormone receptor status.⁶³ Similarly, our analysis found that ACBD5 was upregulated in the luminal subtypes compared to basal and Her2 subtypes. This indicated that ACBD5 could be a candidate biomarker for ER signal reprogramming in precancerous breast tissues. Overall, the upregulation of ACBD5 indicates increased peroxisomal activity and potential involvement in cancer, warranting further investigation of its role in tumorigenesis and its therapeutic potential.

Peroxisomes specialize in the oxidation of very long-chain fatty acids (VLCFAs) and branched-chain fatty acids, which are beyond mitochondrial capacity.⁶⁴ Enzymes such as ACSL4, SLC27A2, and SLC27A4 in peroxisomes convert fatty acids into their activated forms, fatty acyl-CoAs, which then undergo oxidation to produce ATP.⁶⁵ ACSL5, typically associated with the endoplasmic reticulum and mitochondrial outer membrane, is also localized in peroxisomes and catalyzes the formation of fatty acyl-CoAs from long-chain fatty acids.⁶⁶ High ACSL5 expression has been linked to the prognosis of various cancers, such as colorectal adenocarcinoma, breast cancer, and pancreatic cancer, indicating its potential as a prognostic biomarker.^15,18,67 In a previous study, high ACSL5 expression was associated with better prognosis in breast cancer, which aligns with our findings.⁶⁷ Through single-cell analysis, we identified unique ACSL5 expression in immune cells, particularly in T-cell subsets. This suggests that T cells expressing ACSL5 play a protective role in breast cancer progression, possibly through enhanced fatty acid metabolism that supports anti-tumor immune responses.^68,69 Similarly, SLC27A2, which functions like ACSL5, also showed a protective role in breast cancer, with single-cell analysis indicating its upregulation in proliferative T cells. Unlike PEX3 and ACBD5, which are directly associated with peroxisome biogenesis, ACSL5 expression may also indicate increased mitochondrial and the ER.⁶⁴ Therefore, further functional experiments are required to determine whether upregulated ACSL5 expression increases peroxisomal activity in T cells. Mitochondria and peroxisomes exhibit functional interdependence, particularly in fatty acid oxidation and cellular redox balance. For example, fatty acyl-CoAs undergo oxidation in peroxisomes and provide substrates for the citric acid cycle in mitochondria, supporting additional ATP production.⁷⁰ Dysfunction in peroxisomes, such as catalase inhibition, increases mitochondrial ROS and affects the mitochondrial membrane potential.^71,72 Peroxisomal disorders such as X-ALD result in mitochondrial defects and increased ROS levels.⁷³ These observations suggest a functional interplay between peroxisomes and mitochondria, indicating that they work together to maintain cellular redox balance.

Our study revealed that peroxisomal dysfunction is associated with suppressed immune activity and reduced infiltration of immune cells, particularly CD8 + T-cells. Cells in the tumor microenvironment (TME), including immune and stromal cells, often depend on lipid metabolism to meet their energy requirements.^68,69 We observed that the low-risk subgroup exhibited elevated fatty acid metabolism and increased infiltration of CD8 T cells, suggesting a potential role of enhanced peroxisome activity in these CD8 T cells. In contrast, the high-risk subgroup showed elevated expression levels of PEX3, ACBD5, NOS2, and DAO, indicating increased peroxisome biogenesis and ROS management.

Peroxisomal ROS can oxidize redox-sensitive cysteine residues, potentially disrupting cellular signaling pathways and transcription factors.⁷⁴ For example, hydrogen peroxide generated by fatty acyl CoA oxidase activates NF-κB, a transcription factor involved in inflammation and tumorigenesis, whereas peroxisomal catalase counteracts NF-κB activity.⁷⁵ This balance implies that peroxisomal regulation of oxidative stress in cancer cells may lead to reduced interferon signaling, decreased M1 macrophage polarization, and diminished CD8 + T cell infiltration. The PRG risk score effectively predicted the excluded immunophenotype of the TME in the IMvigor210 cohort, indicating the presence of CD8 T cells that did not efficiently infiltrate the tumor. The PRG risk signature also demonstrated comparable efficacy as a biomarker for predicting immunotherapy response. Disruption of peroxisomal activity may enhance this response by converting the excluded TME into an inflamed TME with increased CD8 T cell infiltration. Autophagy-mediated degradation of peroxisomes, also known as pexophagy, can disrupt peroxisomal activity.⁷⁶ Elevated intracellular ROS levels induce pexophagy in liver and breast cancer cells through ATM-mediated activation of ULK1 and inhibition of mTORC1.⁷⁷ Inhibition of heat shock protein family A member 9 (HSPA9) has also been documented to trigger pexophagy in cancer cells.⁷⁸ However, the effects of pexophagy induction on cancer cell viability and growth have not been thoroughly investigated. Other therapeutic targets were identified using drug sensitivity analysis. In addition to BCL2 inhibitors, several kinase inhibitors, such as polo-like kinase 1 (PLK1: BI-2536), Bruton’s tyrosine kinase (BTK: Ibrutinib), cyclin-dependent kinase 1 (CDK1: RO-330), epidermal growth factor receptor (EGFR: Erlotinib), and p38 MAPK (mitogen-activated protein kinase: BIRB 796), have also been reported as therapeutic drugs for the treatment of high-risk subgroups. These findings provide a foundation for investigating peroxisomal dysfunction as a potential therapeutic target to enhance immune response in cancer treatment. Given the pivotal role of peroxisomes in fatty acid metabolism and redox balance, along with their strong association with CD8+ T cells and immune activation, targeting peroxisomes may offer advantages over conventional therapies. Nevertheless, further functional studies are needed to confirm and refine these observations for therapeutic application.

The role of peroxisomes in cancer, particularly through bioinformatic analysis, remains underexplored. A review of the existing literature shows that three previous studies have employed similar methodologies to investigate the role of peroxisome-related genes (PRGs) in specific cancers, including hepatocellular carcinoma (HCC), lower-grade glioma (LGG), and renal clear cell carcinoma (RCC).^27–29 However, the prognostic PRGs identified in each study varied significantly, suggesting differential regulation of peroxisomes across cancer types. Specifically, 14 PRGs were identified in LGG, 9 in RCC, and 10 in HCC. Of the genes identified in our study, only three have previously been associated with risk in other cancers: DAO in RCC, and ABCD5 and ACSL5 in LGG. While our findings overlap only minimally with previously identified prognostic PRGs in other cancers, this variability underscores the distinct roles that peroxisomes could play in tumor biology. Future research is needed to further elucidate these differences, potentially opening new avenues for targeted therapies based on peroxisomal function.

This study has several limitations. The use of retrospective data necessitates prospective experimental validation in patients with breast cancer patients and mouse models to confirm the diagnostic and prognostic potential of the signature. The clinical application of the gene signature is currently limited because it requires preset risk-score thresholds and data normalization in a large, pre-collected cohort. These steps are essential for future research and validation.

Conclusions

Targeting cancer metabolism, particularly peroxisome metabolism, has emerged as a crucial strategy for combating breast cancer because of the metabolic plasticity of cancer cells that enables resistance to therapies. Our study identified a novel seven-gene peroxisome-related risk signature (PRG) that predicts breast cancer prognosis. The PRG risk signature, validated across multiple cohorts, effectively differentiates patient outcomes, particularly within luminal subtypes, and correlates with immune suppression and reduced CD8+ T cell infiltration in high-risk groups. Our findings suggest that peroxisome dysfunction in breast cancer may lead to poor prognosis by modulating the immune response. In addition, the PRG risk signature shows promise as a predictive biomarker for immunotherapy responses, potentially guiding therapeutic strategies. These results underscore the need for further investigation into the role of peroxisomes in breast cancer and their potential as therapeutic targets.