Role of ADME genes in breast cancer prognosis: an analysis of risk scoring models based on multi-omics data

Background: Breast cancer (BC) is a significant malignancy characterized by a high global incidence and a propensity for recurrence. Absorption, distribution, metabolism, and excretion (ADME) genes comprise a collection of genes that participate in the drug ADME. Understanding the role and prognostic value of ARGs (ADME related genes) in BC advancement is critical for personalized therapy. Therefore, an ARPS (ADME related prognostic signature) was created in this study to examine the clinical implications of ARGs in patients with BC.

Methods: A multi-omics investigation of ADME-related genes in BC was conducted using bulk RNA sequencing, single-cell RNA sequencing, and spatial transcriptome data. According to the expression profiles of ADME-related differentially expressed genes (DEGs), the ARPS was calculated, and all patients were stratified based on their risk scores. A prediction model was then created using Cox regression and stepAIC analyses. This model divided all patients into HR (High risk) and LR (Low risk) groups following the median risk score. Bioinformatics analyses were conducted to estimate the risk signature’s predictive capacity.

Results: This study identified five ARGs (SLC7A5, HSD11B1, ADHFE1, GSTM2, and TAP1) correlated with BC prognosis. The risk signature in the TCGA-BRCA, METABRIC, and GSE58812 cohorts revealed robust predictive accuracy for 1-, 3-, and 5-year survival. Compared to the gene signature alone, the nomogram integrating the ARPS and clinical parameters demonstrated improved prognostic performance. Immune infiltration analysis revealed a high level of immune checkpoint related gene expression and immune score in patients with ARPS LR, suggesting potential implications for immunotherapy responses.

Conclusion: The findings highlight the prognostic significance of ARPS in BC and its potential utility in guiding personalized treatment strategies. Combining ARPS with clinical parameters enhances prognostic accuracy and may help patients with BC make clinical decisions.

Introduction

In 2020, breast cancer (BC) replaced lung cancer as the leading cause of global cancer incidence while remaining the fifth most prevalent cause of cancer mortality worldwide (1). There are different types of BC, and even though many patients may have better outcomes than those with other solid tumors after different treatments, such as radical surgery, chemotherapy, radiation therapy, and targeted therapy, some patients with BC still have poor outcomes (2). Consequently, there is a continued need for novel biomarkers to identify patients with BC.

Genes implicated in drug absorption, distribution, metabolism, and excretion are designated as ADME genes (ARGs) (3, 4). The PharmaADME group has disclosed that ARGs comprise 266 extended genes and 32 core genes (http://www.pharmaadme.org), which include phase I and II drug-metabolizing enzymes, modifiers, and transporters that influence hepatic drug clearance and metabolism (5–7).

ARGs are extensively utilized in cancer research to understand their expression profiles in different cancer types and their impact on patient outcomes. Researchers have identified the differential expression of ARGs in tumors, with some genes related to favorable overall survival (OS) rates in certain cancers, while others are linked to unfavorable outcomes (8–10). Studies have revealed that ARGs may affect the survival of patients with cancer through various mechanisms related to drug metabolism and disposition. Tang et al. developed a novel ADME-related 14-gene prognosis model in HNSCC (11). The model assigned patients into two groups, LR or HR, and the results revealed that patients with LR have significantly improved OS and DFS and benefit more from immunotherapy and chemotherapy.

Moreover, ARGs are being explored as possible therapeutic targets and prognostic biomarkers for cancer treatment, highlighting their importance in personalized medication and patient management improvement. Wang et al. established a risk score signature based on ARGs that distinguishes HR from patients with LR sarcoma, demonstrating longer survival in the LR group and offering a direction for future targeted therapies (12). However, the biological roles and predictive value of ARGs in BC remain poorly understood.

In this study, DEARGs (Differentially Expressed ARGs) were explored in BRCA using data from The Cancer Genome Atlas (TCGA). Then, a 5-gene signature was established to predict survival outcomes in the TCGA training cohort using Cox regression and stepAIC analyses, and its prognostic usefulness was extensively validated using external cohorts. Furthermore, the underlying connotations between the signature and landscape of the tumor microenvironment (TME), namely the expression level of immune checkpoints, predictive enrichment of tumor-infiltrating immune cells, and the level of tumor mutation, were revealed, providing novel insights for personalized immunotherapy. Finally, the predictive value of the signature in patients with BC treated with diverse therapeutic modalities was verified. Figure 1 provides is a workflow diagram of our study.

Figure 1

Figure 1. Flowchart for comprehensive analysis of ADME related patterns in postoperative patients with breast cancer (BC).

Materials and methods

Preprocessing and dataset source

In this study, the training cohort consisted of currently accessible treatment and expression data of patients with BRCA obtained from TCGA (https://portal.gdc.cancer.gov/). The “TCGAbiolinks” (13) package in R was used to download RNA-seq data, the mutation profiles, and relevant clinical information. Additionally, information from the cBioPortal website database (https://www.cbioportal.org/) (14) and Gene Expression Omnibus database (http://www.ncbi.nlm.nih.gov/geo/) was used to obtain two independent validation cohorts: METABRIC and GSE58812. The human BRCA scRNA-seq dataset, GSE176078, was obtained from TISCH (http://tisch.comp-genomics.org/) (15). The BRCA spatial transcriptome dataset GSE203612-GSM6177603 was obtained from the Gene Expression Omnibus database (https://www.ncbi.nlm.nih.gov/geo/). A group of ARGs (n = 298) was derived from the PharmaADME Consortium (http://www.pharmaadme.org) (5).

Processing BC spatial transcriptome sequencing data

Seurat (16), an R program for analyzing spatial transcriptome data, was used. This entailed normalizing unique molecular identifier numbers, scaling the data, and identifying the most variable characteristics with “SCTransform.” Downscaling and unsupervised cluster analysis were then performed with “RunPCA.” For the cluster analysis, default parameters were used, focusing on the 30 most important principal components. The “SpatialFeaturePlot” function was also used for subgroup and gene visualization. The “AUCell” (17) R package is a useful tool for spatial transcriptome ADME-related gene analysis. Its primary goal is quantifying and exhibiting ADME-related activities at geographic transcriptome resolution.

scRNA-seq data analysis

The 10× scRNA-seq data GSE176078 was transformed to a Seurat object using the “Seurat” R package. Clusters with fewer than three cells, less than 50 genes, and more than 5% of the mitochondrial genes were deleted. Principal component analysis was applied to the top 1500 variable genes. The “FindNeighbors” and “FindClusters” routines were used to perform cell clustering analysis based on the top 15 principal components. The “FindAllMarkers” program was used to identify marker genes in various cell clusters, with FDR < 0.01 and |log2FC| > 1 as the criterion. Clusters were annotated using the “CellMarker (version 2.0)” (18) database to identify various cell types. The “ssGSEA” function from the Seurat package was used to quantify the activity of a specific gene set in each cell.

Identification of differential ADME related genes and functional analysis

For BRCA and normal cases in TCGA datasets, DEARGs were identified using the R package “limma” (19). Notably, the cutoff was defined as FDR < 0.05 and log2|FC| > 1. ARGs were characterized using the Kyoto encyclopedia of genes and genomes (KEGG) and gene ontology (GO) pathways, identified using the Metascape (20) website, to investigate their potential biological roles and signaling pathways. In the case of FDR < 0.1 and p < 0.05, the result was determined to be statistically significant.

Construction and validation of ARPS

The ADME-related prognostic genes were identified from DEARGs using univariate Cox proportional hazards analysis with a threshold of p < 0.05. The stepAIC analysis using the MASS package was then performed to identify the most predictive ARGs for BC prognosis in line with DEARGs. The most predictive ARGs were subsequently loaded for further analysis using the multivariate Cox proportional hazards regression model. The risk scores were computed by combining the expression of each DEARG and the relevant coefficient. All patients were divided into LR and HR groups based on median risk scores.

Additionally, time-dependent ROC and KM analyses were used to assess ARPS’s predictive performance. An external independent validation cohort used data from METABRIC and GSE58812 to determine the generalization degree of ARPS. Prognostic independence was assessed for clinical parameters, including ADME-related risk scores, in patients with BC using univariate and multivariate Cox regression analyses. Key risk variables were included to develop a nomogram for predicting survival. Calibration curves were plotted, and decision curve analysis (DCA) was performed to assess the nomogram’s precision.

Immunogenomic landscape assessment

The relative proportions of the 22 types of immune cells were estimated using the CIBERSORT (21) R package. The R package “estimate” was then selected to compute the score of ESTIMATE, stroma, and immunity to evaluate the tumor purity. The hypoxia score of BC was acquired from the cBioPortal (https://www.cbioportal.org/).

Analysis of gene mutations

A waterfall diagram displaying the distribution of genes with high somatic mutation frequency in patients with BC was created using the “maftools” (22) R package TCGA provided the copy number variation (CNV) data; the patients in several risk categories were examined using the GenePattern “gistic2” module (23). The output results were demonstrated using the ChromPlot feature in the R package “maftools.” Concurrent with this, the tumor mutation burden (TMB) of every sample was computed to investigate the correlation between the risk score and TMB.

Statistical analysis

The R software (version 4.3.1) was used for all statistical analyses. The Wilcoxon test was used for pairwise comparisons between two groups; the Kruskal–Wallis test ***p < 0.001; ****p < 0.0001 was used for multiple group comparisons. Survival analysis was performed using the KM approach and log-rank test. The outcome was considered statistically relevant at p = 0.05.

Results

ADME-related characteristics in spatial transcriptome and scRNA-seq

Following dimensionality reduction clustering, we used SCTransform’s method to adjust for spatial sequencing depth and performed a series of normalization procedures, thereby identifying nine different cell types in space (Supplementary Figure S1). We calculated ADME-related activity in every cell subgroup using the AUCell R package to assess the significance of ADME-associated genes in each cell subset (Figure 2A). ADME-related activity was more abundant in normal cells (Figure 2B). We then computed the association between cell content and ADME-related activity at all sites (Figure 2C) and between cell content and Spearman’s correlation analysis. Obtained from 26 patients with BC, scRNA-seq data consists of 89,471 cells. The cells were grouped into 11 main clusters using marker genes for distinct cell types: B cells, CD4Tconv, CD8Tex, dendritic cells, endothelial cells, fibroblast cells, malignant cells, monocytes/macrophages, plasma cells, SMC cells, and tprolif cells (Figure 2D). We computed the expression levels of ADME-related genes across all cells using the “ssGSEA” function in the Seurat package to estimate ADME activity in various cell types (Figure 2E). Of the 11 cell types, dendritic and fibroblast cells exhibited specifically more ADME ADME-related activity (Figure 2F).

Figure 2

Figure 2. ADME-related characteristics in the spatial transcriptome and scRNA-seq. (A) Spatial visualization of the ADME intensity. (B) Differential analysis of ADME-related activity in mixed malignant and normal regions. (C) Spearman’s correlation of ADME-related activity with microenvironmental components at spatial transcriptome resolution. (D) Types of single cells identified using marker genes. (E) ADME enrichment score (activity) in each cell. (F) Distribution of ADME scores in different cell types. ADME, Absorption, distribution, metabolism, and excretion.

Variant landscape of ARGs in patients with BC

First, ARS was calculated using the ssGSEA algorithm for 1050 patients with BRCA and 98 normal patients, derived from the TCGA database. Compared to the control group, a significantly lower ARS was observed in the BRCA group, which confirmed the model accuracy (Figure 3A). Then, the association between ARS and clinical parameters in patients with BRCA was further investigated. A higher ARS was associated with a higher stage (Figure 3B). In this study, we screened 36 DEGs for the BRCA and normal groups to investigate their distinct transcriptomic signatures (Figures 3C, D; Supplementary Table S1). The expression of each of the 36 DEGs in the TCGA-BRCA cohort is displayed in Figure 3E. Subsequently, the intricate relevance of DEG-associated proteins was clarified by constructing a network of interactions between proteins (Figure 3F).

Figure 3

Figure 3. Different landscapes of ARGs in patients with BC. (A) Box plot with Wilcoxon test displaying the variations in the SSGSEA ARS between BRCA and normal samples. (B) Box plot contrasting SSGSEA ARS across stages. FDR = 0.05; |log2FC| > 1. (C) Volcano plot of DEGs in BRCA (blue, downregulated DEGs; red, upregulated DEGs; and gray, unaltered genes). (D) Venn diagram of ARGs and BRCA-associated DEGs. (E) Differential analysis heatmap comparing BRCA with the normal group. The normal group is blue; the BRCA group is red; the blue square denotes low expression; the yellow square denotes strong expression. (F) PPI network, including DEGs connected to ADME. (G) Oncoplot of the ADME-related DEGs in the TCGA-BRCA cohort. (H) Frequencies of CNV gain, loss, and non-CNV among the DEGs associated with ADME. (I) ADME-related DEG GO keywords and KEGG pathway enrichment analysis. Several hues represent distinct terms or paths.

Additionally, the molecular alteration landscape of ADME-related DEGs in BRCA was examined, and the most prevalent variant was the missense mutation (Figure 3G). ABCC9, FMO2, ALDH1A3, AOX1, and SULF1 were the top five mutated genes. The top 20 mutated ADME-related DEGs exhibited significant CNV alterations based on CNV mutation frequency (Figure 3H). The Metascape website was employed for GO and KEGG enrichment analyses to examine the regulatory mechanisms of DEGS. According to the enrichment analysis, the most enriched functions included drug metabolism-cytochrome P450, detoxification, retinol metabolism, xenobiotic metabolic process, and active transmembrane transporter activity (Figure 3I).

Construction of ARPS in BC

A univariate Cox regression analysis was used to screen for ARGs with prognostic significance to build a prognostic gene model from the ADME-associated DEGs. Therefore, seven genes were found to have significant prognostic values (p < 0.05) (Figure 4A; Supplementary Table S2). StepAIC analyses were performed to minimize the gene count and streamline the model, resulting in a final set of five ARGs with coefficients for the prognosis model (Figures 4B, C; Supplementary Table S2). The 5-gene prognosis model was defined as follows: Risk score = (0.1279) * SLC7A5 + (-0.2116) * HSD11B1 + (-0.2388) * ADHFE1 + (-0.1828) * GSTM2 + (-0.2116) * TAP1. Patients with BRCA were stratified into HR (n = 525) and LR (n = 525) groups according to the median signature. Notably, in the TCGA cohort, the LR group outperformed the HR group regarding the OS rate (p < 0.0001, Figure 4D). The independent validation cohorts, METABRIC and GSE58812, were utilized to confirm the model’s robustness, and the results demonstrated that the LR group also outperformed the HR group regarding OS (METABRIC, p < 0.0001, Figure 4E; GSE58812, p = 0.029, Figure 4F). The survival status and risk scores distribution are illustrated in Figures 4G–I for TCGA-BRCA, METABRIC, and GSE58812 cohorts, respectively. These findings confirm the strong performance of the ADME-related prognosis model for predicting patients with BC outcomes across various datasets.

Figure 4

Figure 4. Development of an ARPS for patients with BC. (A) The forest plot illustrates the seven prognostic genes identified in the risk model derived from the univariate regression analysis. (B) The forest plot illustrates the final five prognostic genes identified in the risk model derived from stepAIC regression analysis. (C) Regression coefficients of five genes derived from stepAIC regression analysis. (D–F) OS in low- and high-risk patients in (D) TCGA-BRCA, (E) METABRIC, (F) GSE58812. (G–I) Distribution of risk scores based on survival status and time in (G) TCGA-BRCA, (H) METABRIC and (I) GSE58812. StepAIC refers to the stepwise application of the Akaike Information Criterion.

Development and evaluation of a nomogram

This study employed univariate and multivariate Cox regression analyses to evaluate the effect of ADME-related signatures on prognosis prediction. The univariate Cox regression analysis indicated that age, TNM.T, stage, TNM.M, grade, neoplasm status, and risk score were significantly associated with OS. The multivariate Cox regression analysis indicated that stage, age, risk score, and radiation were associated with OS in patients with BC (Figures 5A, B). In this study, we developed a nomogram for predicting 1-, 3-, and 5-year survival in patients with BC, utilizing the correlation between clinicopathological features and the ADME-related signature (Figure 5C). Using the same nomogram, the risk score was computed for all patients, and the patients were assigned in line with their risk scores. The nomogram model outperformed the gene signature model in terms of prognosis. The prognosis of LR and HR groups differed significantly (p < 0.001, Figure 5D).

Figure 5

Figure 5. Development and evaluation of the nomogram survival model. (A) Univariate study of the clinicopathological characteristics and risk score in TCGA-BRCA. (B) Multivariate analysis of clinicopathological characteristics and risk scores of TCGA-BRCA. (C) A nomogram for forecasting the prognosis of patients with BRCA. (D) Kaplan–Meier analysis of the two BRCA cohorts based on the nomogram score. (E) ROC curve study of the nomogram in TCGA-BRCA. (F) Calibration graphs illustrating the probabilities of 1-, 3-, and 5-year overall survival in TCGA-BRCA. (G) DCA of the nomogram forecasting 1-, 3-, and 5-year OS.

Besides, the combined model had AUCs for 1-, 3-, and 5-year survival rates of 0.822, 0.844, and 0.827, respectively (Figure 5E). Moreover, the nomogram’s predictive accuracy was revealed via the calibration curve (Figure 5F). Furthermore, a DCA (Figure 5G) was conducted to compare the nomogram’s clinical applicability regarding the 1-, 3-, and 5-year survival. The results demonstrated that the 3- and 5-year OS was better predicted by the nomogram, which provided more net clinical benefit than the 1-year OS. In general, when utilizing these essential clinical parameters to assess the prognosis of patients with BC, the nomogram demonstrated solid prediction power and clinical applicability.

Immune characteristics of ADME-related prognostic subgroups

The TCGA-BRCA cohort was used to investigate the relationship between the ARPS and the immune status of patients. Most normal cells in tumor tissue consist of infiltrating immune and stromal cells that modulate cancer biology besides disrupting tumor signaling. The ESTIMATE algorithm was employed to analyze the immune microenvironment to elucidate the relationship between the ADME-related risk score signature and its biological function in the immune response. This algorithm demonstrated that a high immune score was noted more frequently in patients with higher risk scores (Figure 6A). A notable reduction in the immune infiltration levels of CD8 T cells, gamma delta T cells, activated NK cells, and activated mast cells was observed in patients with HR.

Figure 6

Figure 6. Immunological attributes of ADME-associated prognostic groupings. (A) Violin plot illustrating reduced immune infiltration, stromal, and ESTIMATE scores, alongside increased tumor purity in patients with HR. (B) A box plot of the 22 invading immune cell types was generated using CIBERSORT. (C) The correlation between TME-infiltrating cells and genes in the ADME-related model. (D) Bubble plot illustrating the average and percentage expression of model genes across various cell subtypes. (E) Box plot illustrating the expression levels of genes related to immunological checkpoints. (F) Violin plot illustrating significant elevation in hypoxic scores among HR patients. *p < 0.05; **p < 0.01; ***p < 0.001; and ****p < 0.0001.

Conversely, an increase was noted in plasma cells, resting memory CD4 T cells, and M2 macrophages (Figure 6B). Moreover, five genes within the prognostic model were highly correlated with TIICs (Figure 6C). A violin plot revealed that SLC7A5 and TAP1 were highly expressed in dendritic and Tprolif cells (Figure 6D). ICIs are antitumor immunotherapies. They are frequently utilized in clinical practice. The differential expression of ICGs between HR and LR groups may lead to a distinct susceptibility to ICIs. As depicted in Figure 6E, most immune checkpoints, such as PDCD1, CD274, CTLA-4, IDO1, TIGIT, and PDCD1LG2, were significantly highly expressed in the LR group. These results indicate that ICI may be associated with ARPS. Additionally, patients with HR exhibited an elevated hypoxia score (Figure 6F) based on the analysis of hypoxia-responsive gene expression.

Genomic modifications in patients classified as low- or high-risk

A greater non-synonymous TMB was observed in the protein-coding regions of the genome in high-risk patients (Figure 7A). The top 20 genes exhibiting the highest mutation frequency were identified in both the risk groups (Figures 7B, C). Notably, the inverse frequency was recorded for PIK3CA (high/low-risk, 26%/39%) and TP53 (high/low-risk, 42%/26%) (Figures 7D, E). The fraction genome change (FGA) was significantly elevated in the high-risk group (Figure 7F). Mutations concentrated in the DNA-binding domain of the respective protein may significantly contribute to the decline in tumor suppression efficacy and reduction in patient survival rates.

Figure 7

Figure 7. Mutation profile among ADME-related prognostic groupings. (A) Analysis of TMB. (B, C) Waterfall plot depicting somatic mutation characteristics in the high- and low-risk score categories. (D, E) Comparison of several TP53 and PIK3CA mutation loci. (F) Analysis of differences in FGA among various risk score groups.

Discussion

BC represents a substantial global health threat due to its increasing prevalence and high mortality (24). Although therapeutic strategies for BC have progressed significantly, the prognosis of patients with BC remains unfavorable, mainly because of resistance to conventional therapies and late diagnosis (25). Several molecular markers for building predictive models for various cancers have been reported in recent studies (26–30). Integrating molecular biomarkers into prognosis models may offer a more comprehensive disease evaluation, improve the accuracy of predicting patient outcomes, and ultimately lead to better patient care and outcomes in BC management.

In this study, a novel ADME-related prognostic model was constructed based on the TCGA-BRCA cohort, and its robustness was validated using the external METABRIC and GSE58812 cohorts. The risk scores were calculated using Cox regression and stepAIC analysis to predict the patients’ prognosis in BC. The analysis indicated that patients in the high-risk group exhibited shorter survival times in the training and validation cohorts. We then identified whether the prognostic signature was correlated with immunomodulators, the tumor microenvironment, and gene mutations.

In this study, a novel prognostic model was created based on five ARGs: SLC7A5, HSD11B1, ADHFE1, GSTM2, and TAP1. SLC7A5, a key component of AA transporters, is essential for sodium-independent transport of large neutral amino acids across the membrane (31). This biological process is critical for the rapid growth and proliferation of tumor cells (32). Altered regulation of SLC7A5 has been observed in various cancer types, including ovarian and non-small-cell lung cancers (33–35). HSD11B1, identified among the five prognostic genes in this study, regulates glucocorticoid levels and is associated with unfavorable outcomes in patients diagnosed with clear cell renal cell carcinoma (36). Single nucleotide polymorphisms in HSD11B1 may correlate with BC risk in postmenopausal women (37). Alcohol dehydrogenase iron-containing 1 (ADHFE1) is an oncogene associated with BC that negatively affects patient survival rates (38).

Additionally, ADHFE1 facilitated metabolic reprogramming characterized by elevated levels of D-2HG and reactive oxygen species, a shift toward reductive glutamine metabolism, and alterations in the epigenetic landscape (38). The glutathione-S-transferase family (GST) represents a significant group of antioxidant enzymes in living organisms, with GSTM classified under the Mu subfamily of GSTs (39). Reports indicate that GSTM enzymes are involved in the metabolism of tumor chemotherapeutic drugs and in protecting organelles or cells from stress injuries (40). GSTM2 encodes a protein containing a GST structural domain integral to numerous biological functions (41). Guo E et al. demonstrated that GSTM2 low expression in colon cancer correlates with improved patient prognosis, as evidenced by bioinformatics analysis (42). Transporter associated with antigen processing 1 (TAP1) is a transporter protein that presents tumor antigens in the MHC I or HLA complex. Mutant TAP1 has been reported to influence MHC-I function in tumor surveillance (43). A prior study suggested that TAP may serve as a cancer treatment through immunotherapy, given its role in the peptide-MHC I complex and its ability to enhance the immune response (44).

A powerful tool for risk stratification of patients is the in-depth molecular characterization of tumor heterogeneity. This study revealed that compared to patients with LR, patients with HR have a lower immunological checkpoint level, hypoxia score, and a higher degree of tumor mutation. According to the findings, patients in the HR group responded poorly to immunotherapy. Moreover, in the HR group, tumor suppressor genes (for example, TP53) revealed a greater mutation frequency; conversely, tumor-promoting genes (PIK3CA) revealed reduced mutation frequency. The poor prognosis of high-risk patients might have their origin in these biological carcinogens.

The significant differences observed in the immunogenomic landscape between the high- and low-risk groups are particularly noteworthy for potential clinical applications. The finding that the low-risk group exhibited high immune scores and, critically, significantly higher expression levels of key immune checkpoint genes such as PDCD1, CD274, and CTLA-4, strongly suggests that patients in the low-risk group may be more likely to respond favorably to immune checkpoint inhibitor therapies. Conversely, the high-risk group’s association with elevated hypoxia, known to create an immunosuppressive environment, alongside reduced infiltration of cytotoxic T cells (like CD8 T cells and gamma delta T cells) and increased presence of M2 macrophages, indicates a potentially less ‘inflamed’ or more immunosuppressive tumor microenvironment, which is generally associated with poorer response to immunotherapy. These immunological distinctions identified by the ARPS highlight their potential as a predictive biomarker for stratifying patients who might benefit most from immune-based treatments.

Furthermore, the results demonstrated that the LR and HR groups differed significantly in ICI characteristics. Patients with HR exhibited a much lower immune infiltration level of T cells CD8, T cells gamma delta, increased plasma cell levels, and macrophages M2. The disease progression and patient prognosis were influenced by important factors such as the composition and dynamics of ICI in TME of BC, particularly the levels of CD8+ T cells. A higher CD8+ T cell infiltration level in BC generally correlates with a reduced possibility of metastatic disease and a better prognosis (45). Increased interest in using tumor-infiltrating gamma-delta T cells for cancer immunotherapy results from their ability to recognize stress antigens in an MHC-unconstrained manner and their correlation with favorable prognosis across many cancer types (46). A study linked poor prognosis in BC to the presence of invading plasma cells (47). Macrophages, the key elements of the TME, are mostly associated with a poor prognosis (48). By reacting to cancer cell-secreted factors, such as macrophage-CSF and GM-CSF, the recruited macrophages typically polarize toward the M2 phenotype during tumor growth (49–51). Crucially important for protumorigenic activities, M2 macrophages express arginase 1 and extensively generate cytokines, growth factors, and proteases. M2 macrophages also induce cancer cell migration, invasion, immunosuppression, and matrix remodeling (52, 53).

Limitations

While this study presents a novel and validated ADME-related prognostic signature for breast cancer, it is important to acknowledge its limitations, which also pave the way for essential future research. Our analysis, primarily based on bulk RNA-seq of overall cohorts, did not allow for a detailed investigation of how the ARPS performs or how the five genes function within distinct molecular subtypes of breast cancer (e.g., Luminal, HER2+, Triple-Negative), which have vastly different biologies and treatment responses. Similarly, a comprehensive assessment of the signature’s predictive value in patients receiving specific treatment modalities beyond general outcomes (such as different chemotherapy regimens, targeted therapies, or specific neoadjuvant/adjuvant settings) was limited by the scope and detailed treatment information available in the public databases, as noted. Furthermore, while we discussed potential roles based on known gene functions, delineating the precise underlying biological mechanisms and interactions of this specific 5-gene combination in breast cancer progression and response requires dedicated experimental validation. Investigating the signature’s performance and the function of these genes in diverse ethnic populations is also a critical area for future study to ensure broad clinical applicability. Therefore, while our ARPS shows promise for risk stratification and suggests implications for immunotherapy, future prospective studies, analyses within specific BC subtypes, and detailed mechanistic investigations are needed to fully understand its clinical utility and biological underpinnings.

Conclusions

This study generated a new ARPS for BC using a TCGA cohort and verified it in two outside cohorts in general. The risk model performed well in predicting patient survival. Besides, there was an association between ARS and ICI in BC. This study offers a promising ARPS to guide the personalized treatment of patients with BC. Furthermore, it provides new insights into possible immunotherapeutic and combined strategies for BC, as targeting ARGs may reverse ICI in BC.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Author contributions

JJ: Conceptualization, Data curation, Formal Analysis, Writing – original draft, Writing – review & editing. XZ: Writing – review & editing, Conceptualization, Formal Analysis. MD: Writing – review & editing, Formal Analysis, Data curation, Writing – original draft. LD: Writing – review & editing, Resources, Supervision.

Funding

The author(s) declare that no financial support was received for the research and/or publication of this article.

Acknowledgments

We thank all participants involved in this study. We also acknowledge the investigators who participated and provided data unselfishly in TCGA database (TCGA-BRCA), cBioPortal database (METABRIC) and GEO database (GSE58812, GSE176078 and GSE203612-GSM6177603).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2025.1582862/full#supplementary-material

Supplementary Table 1 | ADME related differentially expressed genes.

Supplementary Table 2 | Construction of an ADME-related prognostic signature for BC patients.

Supplementary Figure 1 | 9 distinct cell types in space of BC.

Abbreviations

BC, Breast cancer; ADME, Absorption, distribution, metabolism, and excretion; GO, Gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; LR, Low Risk; HR, High Risk; HGNC, HUGO Gene Nomenclature Committee; DEG, Differentially expressed genes.

References

1. Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataram I, Jemal A, et al. Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA: A Cancer J clinicians. (2021) 71:209–49. doi: 10.3322/caac.21660

PubMed Abstract | Crossref Full Text | Google Scholar

2. Hua X, Long ZQ, Zhang YL, Wen W, Guo L, Xia W, et al. Prognostic value of preoperative systemic immune-inflammation index in breast cancer: A propensity score-matching study. Front Oncol. (2020) 10:580. doi: 10.3389/fonc.2020.00580

PubMed Abstract | Crossref Full Text | Google Scholar

3. Tao L, Zhang P, Qin C, Chen SY, Zhang C, Chen Z, et al. Recent progresses in the exploration of machine learning methods as in-silico ADME prediction tools. Advanced Drug delivery Rev. (2015) 86:83–100. doi: 10.1016/j.addr.2015.03.014

PubMed Abstract | Crossref Full Text | Google Scholar

4. Gleeson MP, Hersey A, Montanari D, and Overington J. Probing the links between in vitro potency, ADMET and physicochemical parameters. Nat Rev Drug Discov. (2011) 10:197–208. doi: 10.1038/nrd3367

PubMed Abstract | Crossref Full Text | Google Scholar

5. Hu DG, Marri S, McKinnon RA, Mackenzie PI, and Meech R. Deregulation of the genes that are involved in drug absorption, distribution, metabolism, and excretion in hepatocellular carcinoma. J Pharmacol Exp Ther. (2019) 368:363–81. doi: 10.1124/jpet.118.255018

PubMed Abstract | Crossref Full Text | Google Scholar

6. Hovelson DH, Xue Z, Zawistowski M, Ehm MG, Harris EC, Stocker SL, et al. Characterization of ADME gene variation in 21 populations by exome sequencing. Pharmacogenetics Genomics. (2017) 27:89–100. doi: 10.1097/FPC.0000000000000260

PubMed Abstract | Crossref Full Text | Google Scholar

7. Jittikoon J, Mahasirimongkol S, Charoenyingwattana A, Chaikledkaew U, Tragulpiankit P, Mangmool S, et al. Comparison of genetic variation in drug ADME-related genes in Thais with Caucasian, African and Asian HapMap populations. J Hum Genet. (2016) 61:119–27. doi: 10.1038/jhg.2015.115

PubMed Abstract | Crossref Full Text | Google Scholar

8. De Iuliis F, Salerno G, Taglieri L, and Scarpa S. Are pharmacogenomic biomarkers an effective tool to predict taxane toxicity and outcome in breast cancer patients? Literature review. Cancer chemotherapy Pharmacol. (2015) 76:679–90. doi: 10.1007/s00280-015-2818-4

PubMed Abstract | Crossref Full Text | Google Scholar

9. Suthandiram S, Gan GG, Zain SM, Bee PC, Lian LH, Chang KM, et al. Effect of polymorphisms within methotrexate pathway genes on methotrexate toxicity and plasma levels in adults with hematological Malignancies. Pharmacogenomics. (2014) 15:1479–94. doi: 10.2217/pgs.14.97

PubMed Abstract | Crossref Full Text | Google Scholar

10. Zhang C, Ma Q, Shi Y, Li X, Wang M, Wang J, et al. A novel 5-fluorouracil-resistant human esophageal squamous cell carcinoma cell line Eca-109/5-FU with significant drug resistance-related characteristics. Oncol Rep. (2017) 37:2942–54. doi: 10.3892/or.2017.5539

PubMed Abstract | Crossref Full Text | Google Scholar

11. Tang X, Li R, Wu D, Wang Y, Zhao F, Lv R, et al. Development and validation of an ADME-related gene signature for survival, treatment outcome and immune cell infiltration in head and neck squamous cell carcinoma. Front Immunol. (2022) 13:905635. doi: 10.3389/fimmu.2022.905635

PubMed Abstract | Crossref Full Text | Google Scholar

12. Wang J, Wang G, Hu T, Wang H, and Zhou Y. Identification of an ADME-related gene for forecasting the prognosis and responding to immunotherapy in sarcomas. Eur J Med Res. (2024) 29:45. doi: 10.1186/s40001-023-01624-3

PubMed Abstract | Crossref Full Text | Google Scholar

13. Colaprico A, Silva TC, Olsen C, Garofano L, Cava C, Garolini D, et al. TCGAbiolinks: an R/Bioconductor package for integrative analysis of TCGA data. Nucleic Acids Res. (2016) 44:e71. doi: 10.1093/nar/gkv1507

PubMed Abstract | Crossref Full Text | Google Scholar

14. Gao J, Aksoy BA, Dogrusoz U, Dresdner G, Gross B, Sumer SO, et al. Integrative analysis of complex cancer genomics and clinical profiles using the cBioPortal. Sci Signaling. (2013) 6:l1. doi: 10.1126/scisignal.2004088

PubMed Abstract | Crossref Full Text | Google Scholar

15. Sun D, Wang J, Han Y, Dong X, Ge J, Zheng R, et al. TISCH: a comprehensive web resource enabling interactive single-cell transcriptome visualization of tumor microenvironment. Nucleic Acids Res. (2021) 49:D1420–30. doi: 10.1093/nar/gkaa1020

PubMed Abstract | Crossref Full Text | Google Scholar

16. Slovin S, Carissimo A, Panariello F, Grimaldi A, Bouché V, Gambardella G, et al. Single-cell RNA sequencing analysis: A step-by-step overview. Methods Mol Biol (Clifton N.J.). (2021) 2284:343–65. doi: 10.1007/978-1-0716-1307-8_19

PubMed Abstract | Crossref Full Text | Google Scholar

17. Lu Y, Li K, Hu Y, and Wang X. Expression of immune related genes and possible regulatory mechanisms in alzheimer’s disease. Front Immunol. (2021) 12:768966. doi: 10.3389/fimmu.2021.768966

PubMed Abstract | Crossref Full Text | Google Scholar

18. Hu C, Li T, Xu Y, Zhang X, Li F, Bai J, et al. CellMarker 2.0: an updated database of manually curated cell markers in human/mouse and web tools based on scRNA-seq data. Nucleic Acids Res. (2023) 51:D870–6. doi: 10.1093/nar/gkac947

PubMed Abstract | Crossref Full Text | Google Scholar

19. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. (2015) 43:e47. doi: 10.1093/nar/gkv007

PubMed Abstract | Crossref Full Text | Google Scholar

20. Zhou Y, Zhou B, Pache L, Chang M, Khodabakhshi AH, Tanaseichuk O, et al. Metascape provides a biologist-oriented resource for the analysis of systems-level datasets. Nat Commun. (2019) 10:1523. doi: 10.1038/s41467-019-09234-6

PubMed Abstract | Crossref Full Text | Google Scholar

21. Chen B, Khodadoust MS, Liu CL, Newman AM, and Alizadeh AA. Profiling tumor infiltrating immune cells with CIBERSORT. Methods Mol Biol (Clifton N.J.). (2018) 1711:243–59. doi: 10.1007/978-1-4939-7493-1_12

PubMed Abstract | Crossref Full Text | Google Scholar

22. Mayakonda A, Lin DC, Assenov Y, Plass C, and Koeffler HP. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res. (2018) 28:1747–56. doi: 10.1101/gr.239244.118

PubMed Abstract | Crossref Full Text | Google Scholar

23. Kuehn H, Liberzon A, Reich M, and Mesirov JP. Using GenePattern for gene expression analysis. Curr Protoc bioinformatics Chapter. (2008) 7:7.12.1–7.12.39. doi: 10.1002/0471250953.bi0712s22

PubMed Abstract | Crossref Full Text | Google Scholar

24. Katsura C, Ogunmwonyi I, Kankam HK, and Saha S. Breast cancer: presentation, investigation and management. Br J Hosp Med (London England: 2005). (2022) 83:1–7. doi: 10.12968/hmed.2021.0459

PubMed Abstract | Crossref Full Text | Google Scholar

25. Barzaman K, Karami J, Zarei Z, Hosseinzadeh A, Kazemi MH, Moradi-Kalbolandi S, et al. Breast cancer: Biology, biomarkers, and treatments. Int Immunopharmacol. (2020) 84:106535. doi: 10.1016/j.intimp.2020.106535

PubMed Abstract | Crossref Full Text | Google Scholar

26. Lv Y, Wu L, Jian H, Zhang C, Lou Y, Kang Y, et al. Identification and characterization of aging/senescence-induced genes in osteosarcoma and predicting clinical prognosis. Front Immunol. (2022) 13:997765. doi: 10.3389/fimmu.2022.997765

PubMed Abstract | Crossref Full Text | Google Scholar

27. Chen Y, He J, Jin T, Zhang Y, and Ou Y. Functional enrichment analysis of LYSET and identification of related hub gene signatures as novel biomarkers to predict prognosis and immune infiltration status of clear cell renal cell carcinoma. J Cancer Res Clin Oncol. (2023) 149:16905–29. doi: 10.1007/s00432-023-05280-2

PubMed Abstract | Crossref Full Text | Google Scholar

28. Hu G, Li J, Zeng Y, Liu L, Yu Z, Qi X, et al. The anoikis-related gene signature predicts survival accurately in colon adenocarcinoma. Sci Rep. (2023) 13:13919. doi: 10.1038/s41598-023-40907-x

PubMed Abstract | Crossref Full Text | Google Scholar

29. Chen Y, Lin QX, Xu YT, Qian FJ, Lin CJ, Zhao WY, et al. An anoikis-related gene signature predicts prognosis and reveals immune infiltration in hepatocellular carcinoma. Front Oncol. (2023) 13:1158605. doi: 10.3389/fonc.2023.1158605

PubMed Abstract | Crossref Full Text | Google Scholar

30. Wang Y, Xu J, Fang Y, Gu J, Zhao F, Tang Y, et al. Comprehensive analysis of a novel signature incorporating lipid metabolism and immune-related genes for assessing prognosis and immune landscape in lung adenocarcinoma. Front Immunol. (2022) 13:950001. doi: 10.3389/fimmu.2022.950001

PubMed Abstract | Crossref Full Text | Google Scholar

31. Kimball SR. Regulation of translation initiation by amino acids in eukaryotic cells. Prog Mol subcellular Biol. (2001) 26:155–84. doi: 10.1007/978-3-642-56688-2_6

PubMed Abstract | Crossref Full Text | Google Scholar

32. Yanagida O, Kanai Y, Chairoungdua A, Kim DK, Segawa H, Nii T, et al. Human L-type amino acid transporter 1 (LAT1): characterization of function and expression in tumor cell lines. Biochim Biophys Acta. (2001) 1514:291–302. doi: 10.1016/s0005-2736(01)00384-4

PubMed Abstract | Crossref Full Text | Google Scholar

33. Abd El-Rehim DM, Ball G, Pinder SE, Rakha E, Paish C, Robertson JF, et al. High-throughput protein expression analysis using tissue microarray technology of a large well-characterised series identifies biologically distinct classes of breast cancer confirming recent cDNA expression analyses. Int J Cancer. (2005) 116:340–50. doi: 10.1002/ijc.21004

PubMed Abstract | Crossref Full Text | Google Scholar

34. Curtis C, Shah SP, Chin SF, Turashvili G, Rueda OM, Dunning MJ, et al. The genomic and transcriptomic architecture of 2,000 breast tumours reveals novel subgroups. Nature. (2012) 486:346–52. doi: 10.1038/nature10983

PubMed Abstract | Crossref Full Text | Google Scholar

35. Kaira K, Oriuchi N, Imai H, Shimizu K, Yanagitani N, Sunaga N, et al. L-type amino acid transporter 1 (LAT1) is frequently expressed in thymic carcinomas but is absent in thymomas. J Surg Oncol. (2009) 99:433–8. doi: 10.1002/jso.21277

PubMed Abstract | Crossref Full Text | Google Scholar

36. Han D, Yu Z, Zhang H, Liu H, Wang B, and Qian D. Microenvironment-associated gene HSD11B1 may serve as a prognostic biomarker in clear cell renal cell carcinoma: a study based on TCGA, RT−qPCR, Western blotting, and immunohistochemistry. Bioengineered. (2021) 12:10891–904. doi: 10.1080/21655979.2021.1994908

PubMed Abstract | Crossref Full Text | Google Scholar

37. Feigelson HS, Teras LR, Diver WR, Tang W, Patel AV, Stevens VL, et al. Genetic variation in candidate obesity genes ADRB2, ADRB3, GHRL, HSD11B1, IRS1, IRS2, and SHC1 and risk for breast cancer in the Cancer Prevention Study II. Breast Cancer research: BCR. (2008) 10:R57. doi: 10.1186/bcr2114

PubMed Abstract | Crossref Full Text | Google Scholar

38. Mishra P, Tang W, Putluri V, Dorsey TH, Jin F, Wang F, et al. ADHFE1 is a breast cancer oncogene and induces metabolic reprogramming. J Clin Invest. (2018) 128:323–40. doi: 10.1172/JCI93815

PubMed Abstract | Crossref Full Text | Google Scholar

39. Pearson WR, Vorachek WR, Xu SJ, Berger R, Hart I, Vannais D, et al. Identification of class-mu glutathione transferase genes GSTM1-GSTM5 on human chromosome 1p13. Am J Hum Genet. (1993) 53:220–33.

PubMed Abstract | Google Scholar

40. Pljesa-Ercegovac M, Savic-Radojevic A, Matic M, Coric V, Djukic T, Radic T, et al. Glutathione transferases: potential targets to overcome chemoresistance in solid tumors. Int J Mol Sci. (2018) 19:3785. doi: 10.3390/ijms19123785

PubMed Abstract | Crossref Full Text | Google Scholar

41. Nebert DW and Vasiliou V. Analysis of the glutathione S-transferase (GST) gene family. Hum Genomics. (2004) 1:460–4. doi: 10.1186/1479-7364-1-6-460

PubMed Abstract | Crossref Full Text | Google Scholar

42. Guo E, Wei H, Liao X, Wu L, and Zeng X. Clinical significance and biological mechanisms of glutathione S-transferase mu gene family in colon adenocarcinoma. BMC Med Genet. (2020) 21:130. doi: 10.1186/s12881-020-01066-2

PubMed Abstract | Crossref Full Text | Google Scholar

43. Van Kaer L, Ashton-Rickardt PG, Ploegh HL, and Tonegawa S. TAP1 mutant mice are deficient in antigen presentation, surface class I molecules, and CD4-8+ T cells. Cell. (1992) 71:1205–14. doi: 10.1016/s0092-8674(05)80068-6

PubMed Abstract | Crossref Full Text | Google Scholar

44. Lou Y, Vitalis TZ, Basha G, Cai B, Chen SS, Choi KB, et al. Restoration of the expression of transporters associated with antigen processing in lung carcinoma increases tumor-specific immune responses and survival. Cancer Res. (2005) 65:7926–33. doi: 10.1158/0008-5472.CAN-04-3977

PubMed Abstract | Crossref Full Text | Google Scholar

45. Sun X, Zhai J, Sun B, Parra ER, Jiang M, Ma W, et al. Effector memory cytotoxic CD3+/CD8+/CD45RO+ T cells are predictive of good survival and a lower risk of recurrence in triple-negative breast cancer. Modern pathology: an Off J United States Can Acad Pathology Inc. (2022) 35:601–8. doi: 10.1038/s41379-021-00973-w

PubMed Abstract | Crossref Full Text | Google Scholar

46. Boufea K, González-Huici V, Lindberg M, Symeonides S, Oikonomidou O, and Batada NN. Single-cell RNA sequencing of human breast tumour-infiltrating immune cells reveals a γδ T-cell subtype associated with good clinical outcome. Life Sci alliance. (2020) 4:e202000680. doi: 10.26508/lsa.202000680

PubMed Abstract | Crossref Full Text | Google Scholar

47. Parkes H, Collis P, Baildam A, Ralphs D, Lyons B, Howell A, et al. In situ hybridisation and S1 mapping show that the presence of infiltrating plasma cells is associated with poor prognosis in breast cancer. Br J Cancer. (1988) 58:715–22. doi: 10.1038/bjc.1988.296

PubMed Abstract | Crossref Full Text | Google Scholar

48. Nasir I, McGuinness C, Poh AR, Ernst M, Darcy PK, and Britt KL. Tumor macrophage functional heterogeneity can inform the development of novel cancer therapies. Trends Immunol. (2023) 44:971–85. doi: 10.1016/j.it.2023.10.007

PubMed Abstract | Crossref Full Text | Google Scholar

49. Qian BZ and Pollard JW. Macrophage diversity enhances tumor progression and metastasis. Cell. (2010) 141:39–51. doi: 10.1016/j.cell.2010.03.014

PubMed Abstract | Crossref Full Text | Google Scholar

50. Su S, Liu Q, Chen J, Chen J, Chen F, He C, et al. A positive feedback loop between mesenchymal-like cancer cells and macrophages is essential to breast cancer metastasis. Cancer Cell. (2014) 25:605–20. doi: 10.1016/j.ccr.2014.03.021

PubMed Abstract | Crossref Full Text | Google Scholar

51. Sousa S, Brion R, Lintunen M, Kronqvist P, Sandholm J, Mönkkönen J, et al. Human breast cancer cells educate macrophages toward the M2 activation status. Breast Cancer research: BCR. (2015) 17:101. doi: 10.1186/s13058-015-0621-0

PubMed Abstract | Crossref Full Text | Google Scholar

52. Santoni M, Romagnoli E, Saladino T, Foghini L, Guarino S, Capponi M, et al. Triple negative breast cancer: Key role of Tumor-Associated Macrophages in regulating the activity of anti-PD-1/PD-L1 agents. Biochim Biophys Acta Rev Cancer. (2018) 1869:78–84. doi: 10.1016/j.bbcan.2017.10.007

PubMed Abstract | Crossref Full Text | Google Scholar

53. Gao J, Liang Y, and Wang L. Shaping polarization of tumor-associated macrophages in cancer immunotherapy. Front Immunol. (2022) 13:888713. doi: 10.3389/fimmu.2022.888713

PubMed Abstract | Crossref Full Text | Google Scholar