Development and validation of a herb-related gene signature for prognosis prediction and therapeutic response assessment in breast cancer

Introduction

Breast cancer remains the greatest threat to women’s health. Breast tumors are still the number one most prevalent tumor in women (1) and the leading cause of cancer-related deaths in women worldwide (2). Fortunately, as early screening methods and treatments for breast cancer continue to advance, the mortality rate from breast cancer continues to decline (2). Because breast cancer is highly heterogeneous, personalized and precise treatment is especially critical for each patient. Breast cancer is divided into four pathological subtypes based on expression of hormone receptors and human epidermal growth factor receptor 2 (HER2): luminal A, luminal B, HER2, and basal like. Each subtype has its own therapeutic strategies in current clinical practice (3). Compared to regular surgery, chemotherapy and radiotherapy, immunotherapy has become very popular in recent years. This is because the success of immune checkpoint inhibitors in various solid tumors has rekindled interest in breast cancer immunotherapy-based therapies and prevention (4-6).

Traditional Chinese medicine (TCM) has been used to treat tumors for hundreds of years. However, cultural differences and different sources of evidence have led to TCM often being overlooked in international oncology research. However, the exploration of certain TCM-related drugs may offer valuable insights for tumor treatment and diagnosis. For example, the Nobel Prize in Physiology or Medicine in 2015 was awarded to Tu Youyou because of her discovery of artemisinin, which has contributed greatly to the treatment of malaria. Artemisinin is derived from a very common TCM herb, Artemisia annua (7). There are also a number of herbs in TCM that have been used in the treatment of tumors. Astragalus polysaccharide is believed to have an inhibitory effect on HepG2 hepatocellular carcinoma cells (8). Another drug that has also demonstrated anti-liver cancer effects is Poria pachyman, the active component in Poria, which is believed to have inhibitory effects on liver cancer cells (9). Hedyotis diffusa injection induces ferroptosis via the Bax/Bcl2/VDAC2/3 axis in lung adenocarcinoma (10). These herbs have a variety of antitumor mechanisms, and various studies have shown that the mechanisms may all be closely related to the immune system. For example, the active components of several of these herbs have immunomodulatory effects. Astragalus polysaccharide stimulates T-cell immunity via dendritic cells (11). Pachyman has been shown to alter the levels of multiple immune cells and interferons in a mouse model of Kawasaki disease (12). Ethanol extract of Hedyotis diffusa will affect immune responses in normal Balb/c mice in vivo (13).

Much breast cancer research involves TCM, and various Chinese studies have documented many herbal prescriptions for its treatment. Many of these findings found in the Chinese literature are based on theoretical sources of TCM but lack the validation of modern medical methods. Due to its limited specificity, TCM is often accompanied by many side effects when used to treat breast cancer, such as liver and renal function damage and ingredient allergy. Nevertheless, with the development of computational and modern medical technologies, the underlying mechanisms of these ancient prescriptions are presented in a scientific way. For example, the mechanism of inhibition of spontaneous metastasis of breast cancer by Salvia and Ginseng was then validated by modern systemic pharmacology (14).

TCM literature contains numerous prescriptions for breast cancer, encompassing hundreds of distinct herbs. Each herb contains unique active ingredients with the potential to modulate multiple biological targets. However, the targets corresponding to these herbal drugs have not been studied systematically in breast cancer patients. Our team believes that there is great scientific value in these potential herbal therapeutic targets. For clinicians, the treatment of breast cancer is becoming increasingly precise, and more refined models are needed to determine the treatment effect and prognosis. Based on our team’s strong bioinformatics foundation, we explored a meaningful predictive model to determine the prognosis of breast cancer patients.

Assays based on quantitative real-time reverse transcription polymerase chain reaction (PCR) are the most commonly used predictive model or gene signature in clinical practice (15). In addition, Finak et al. identified a 26-gene prediction model. Several of these 26 genes are closely associated with immune responses (16). Current prognostic models for breast cancer, such as the Nottingham Prognostic Index (NPI), PREDICT, Oncotype DX, and MammaPrint, have significantly advanced personalized oncology. However, they are not without limitations. These include a reliance on traditional histopathological parameters with limited molecular insight, a lack of specificity in predicting responses to modern therapies like immunotherapy, and constrained applicability across diverse patient subgroups and molecular subtypes. Furthermore, many existing gene signatures are often derived from a single biological process, potentially overlooking the complex, holistic nature of tumor biology and therapeutic response. Thus, developing a precision treatment for breast cancer requires an understanding of the immune gene network associated with the disease, though this type of research is rare.

To overcome these limitations, a new generation of prognostic models is required. An ideal framework should be built on a solid biological rationale, integrate multi-omics data for a comprehensive perspective, and be robustly validated to ensure clinical reliability. Crucially, it must possess high specificity in predicting responses to diverse therapies, particularly immunotherapy. In this study, we construct a novel prognostic signature by systematically analyzing genes targeted by TCM herbs. We hypothesize that this model, derived from the holistic principles of TCM, will provide a powerful and clinically applicable tool for improving prognosis and guiding personalized treatment in breast cancer.

In our study, we collected a large number of TCM prescriptions for the treatment of breast cancer that could be found over 32 years. The data were analyzed for 221 available prescriptions. A 20-gene signature that can predict the prognosis of breast cancer patients was established by analysis of networks of pharmacological systems and bioinformatics. The model we developed using these genes was shown to have great predictive power and can be validated using various other datasets. It is also useful for exploration of new chemotherapeutic agents and prediction of response to immunotherapy. This predictive model not only provides clinicians with a new reference for treatment but also pioneers a new model of TCM research. This study provides more evidence for subsequent researchers to study and explore new possibilities for breast cancer treatment. We present this article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1259/rc).

Methods

Data collection

Breast cancer patient transcriptome data and clinical information were downloaded from The Cancer Genome Atlas (TCGA: n=1,097, E =458) (https://portal.gdc.cancer.gov/repository). As a first validation set, gene expression information was obtained from the Molecular Taxonomy of Breast Cancer International Consortium (METABRIC: n=1,466, E =401) database (https://www.cbioportal.org/study?id=breastcancer_metabric). Additionally, GSE20685 (n=327, E =73) and GSE10886 (n=226, E =48) datasets from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) were used as validation across multiple independent cohorts (n= the total number of patients, E = the number of outcome events). Normalization of the gene expression profiles was performed using the scale method provided by the R package “limma”. All datasets included only patients with complete data and an overall survival (OS) time longer than 30 days. Additionally, data on gene-level proteomics were obtained from Clinical Proteomic Tumor Analysis Consortium (CPTAC) (https://proteomics.cancer.gov/programs/cptac). Infiltration scores for 24 immune cells were obtained from ImmuCellAI (https://guolab.wchscu.cn/ImmuCellAI#!/) (17). ImmuCellAI estimates the abundance of 24 immune cells based on gene expression datasets such as RNA sequencing (RNA-seq) and microarray data. Therefore, immune cell infiltration can be quantified between groups, and immune checkpoint blockade responses can be predicted. This method involves obtaining a reference expression profile for each cell type from the GEO database and curating a gene signature from a publication. Next, a total expression deviation was calculated based on the gene signatures in the input expression profiles compared with the reference expression profiles of 24 immune cell types. Based on the enrichment fraction of the genetic profiles of the immune cell types, deviations were assigned based on the single-sample gene set enrichment analysis (ssGSEA) algorithm (17). The Profiling Relative Inhibition Simultaneously in Mixtures (PRISM) Repurposing dataset (19Q4, released December 2019) and Cancer Therapeutics Response Portal (CTRP v.2.0, released October 2015) provide information on the drug sensitivity of cancer cell lines (CCLs). The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Study design

First, we collected all literature on TCM prescriptions for breast cancer treatment published in the CNKI (China National Knowledge Infrastructure; www.cnki.net), Wanfang (name of a Chinese academic database; www.wanfangdata.com.cn) and VIP (China Science and Technology Journal Database; www.cqvip.com) databases from 1990 to 2022. We strictly screened out papers that met the standards according to our inclusion criteria and exclusion criteria. Through the HERB database (a high-throughput experiment- and reference-guided database of TCM), we screened 73 potential targets with frequencies greater than or equal to 4. Excluding 13 targets that do not encode proteins, we finally included 60 protein-encoding candidate genes.

A regression algorithm based on the least absolute shrinkage and selection operator (LASSO) with 10-fold cross-validation was used for feature selection with the R package “glmnet” (version: 4.0.2; https://cran.r-project.org/web/packages/glmnet/index.html). The LASSO Cox regression model was fitted using 10-fold cross-validation to determine the optimal value of the penalty parameter (λ). The optimal λ was selected using the ‘lambda.1se’ criterion, which is the value of λ that results in the most regularized model whose cross-validated error is within one standard error of the minimum (the ‘lambda.min’). This criterion promotes model parsimony and enhances the generalizability of the derived gene signature by selecting a simpler model with fewer predictors. The model identified 20 herb-related genes with nonzero coefficients, and the samples were divided into high- and low-risk groups based on the median score. An herb-related risk score (HRS) was developed as follows: HRS = sum (each gene’s expression * corresponding coefficient). The high-low grouping was applied to TCGA as well as three validation sets (METABRIC, GSE20685, and GSE10886). To evaluate the predictive power of the gene signature, time-dependent receiver operating characteristic (tROC) curve analyses were conducted using the survivalROC R package. Independent risk factors for OS of patients were identified based on multivariate Cox regression analysis. The performance of the nomogram was evaluated by calibration plots. It was determined whether the model’s projected probability was consistent with the actual outcome by using the concordance index (C-index). The R package rms was used to plot the nomograms and calibration plots (Version: 4.0.2). To focus on cancer-specific mortality, we excluded patients who died within 30 days of diagnosis, as these early events are often related to non-oncological causes (e.g., post-surgical mortality) rather than tumor progression.

Estimation of immunotherapy and prediction of chemotherapy drugs

Six immune infiltration cell scores were downloaded from Tumor Immune Estimation Resource (TIMER) database (available at http://cistrome.org/TIMER) (18). ESTIMATE (Estimation of Stromal and Immune cells in Malignant Tumor tissues using Expression data), immune, and stromal scores of breast cancer were analyzed based on expression data using the “estimate” R package. Researchers can use this package to determine the purity of tumors, the presence of stromal cells, and the presence of immune cells in tumor tissues (19). ImmuCellAI, a gene set signature-based method for predicting immunotherapy responses with high accuracy by assessing gene expression, provides data regarding immunotherapy response (anti-Programmed cell death protein 1 or anti-Cytotoxic T-lymphocyte–associated protein 4). By using this tool, we were able to divide breast cancer patients into two groups based on their response to immunotherapy; the immunotherapy response group and HRS group were compared using the chi-square test.

We employed the CTRP dataset, which contains sensitivity data for 481 chemicals in 835 CCLs, and the PRISM Repurposing dataset (19Q4, published December 2019), which contains sensitivity data for 1,448 chemicals in 482 CCLs, to gather drug sensitivity data for CCLs. The AUC value is provided in both datasets as a measure of drug sensitivity, with lower values indicating greater sensitivity. The R package pRRophetic (Version: 4.15-1) is equipped with a built-in ridge regression model for predicting chemotherapy responses based on the expression profiles of TCGA samples (20,21). We identified compounds with negative correlation coefficients for CTRP and PRISM by using Spearman correlation (Spearman r=0.25 for CTRP and 0.30 for PRISM). Additionally, drug response was compared between high-HRS (highest decile) and low-HRS (lowest decile) groups.

Statistical analyses

All statistical analyses were conducted using R version 4.0.3. The Mann-Whitney U test and Pearson chi-square test were applied to compare continuous and categorical variables between the training cohort and validation cohort. To validate TCGA set findings using the “c5.go.v7.2.entrez.gmt” gene set, GSEA and METABRIC analysis were conducted using the clusterProfiler R package (version 3.18.0; https://bioconductor.org/packages/clusterProfiler/). To conduct Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses, the R package “clusterProfiler” was utilized based on differentially expressed genes (DEGs) [|log2 fold change (FC)| ≥1, FDR <0.05] between high-risk and low-risk groups. Kaplan-Meier survival curves were constructed using the Kaplan-Meier method, followed by a log-rank test, to determine significant differences. The “prcomp” function in the statistics R package was used to perform principal component analysis (PCA) based on expression of genes in the signature. Moreover, survival analysis was performed on each gene using the “surv_cutpoint” function of the R package “survminer”. Using the R package “timeROC” (Version 0.4; https://cran.r-project.org/web/packages/timeROC/index.html), a tROC curve was calculated to assess the predictive power of the gene signature. The R package “meta” (Version: 4.15-1; https://cran.r-project.org/web/packages/meta/index.html) was utilized to conduct a meta-analysis of the results of multivariate analysis. We used the “RMS” R package to establish the nomogram and calibration curve. All reported P values were derived from two-sided statistical tests. The significance level for all statistical analyses was set at P<0.05. Key estimates, including hazard ratios (HRs) from Cox regression models and area under the curve (AUC) values from tROC analysis, are presented alongside their 95% confidence intervals (CIs) to indicate the precision of the estimates.

Results

Construction of the flow chart

Our study design is illustrated in a schematic diagram in Figure 1. First, 60 genes related to herbs were screened from the published literature. Figure 2A shows the detailed criteria for the literature selection and gene screening procedures. The top 15 herbs that appeared most frequently in all TCM prescriptions are listed in Figure 2B. Detailed information about the screened herbs in the TCM prescription, the potential targets of each herb, and the frequency statistics of targets can be found in table available at https://cdn.amegroups.cn/static/public/tcr-2025-1259-1.xlsx. To predict survival, we used the LASSO algorithm to identify promising genes and to create a robust signature consisting of 20 herb-related genes (Figure 1B). After analyzing the cohort from TCGA for prognostic value, we validated the results with three independent datasets. Moreover, we validated the results of multivariate analysis by performing a meta-analysis (Figure 1C). To determine whether this signature is clinically applicable, we assessed differences in clinicopathological characteristics and immunotherapy responses between risk groups. Finally, a nomogram based on the HRS and other independent factors was developed to predict breast cancer risk assessment and survival probability (Figure 1D).

Figure 1 Schematic diagram of the study design. (A) Sixty potential candidate genes were confirmed by a literature search and HERB database screening. (B) LASSO regression was used to establish a novel and robust herb-related gene signature for breast cancer prognosis. (C) Characterization of herb-related genes and validation of the prognostic value of the gene signature in different cohorts. (D) Immune response prediction according to risk score and potential drug prediction. ****, P<0.0001. AUC, area under the curve; CDA, clinical decision analysis; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; LASSO, least absolute shrinkage and selection operator; METABRIC, Molecular Taxonomy of Breast Cancer International Consortium; PCA, principal component analysis; TCGA, The Cancer Genome Atlas.

Figure 2 Mining and analysis of TCM prescription data. (A) Flow chart of the screening of 60 potential candidate genes. (B) The top 15 herbs were found in TCM prescriptions after screening. TCM, traditional Chinese medicine.

Establishment of an herb-related gene signature for prognosis

A total of 1,097 patients from the TCGA-BRCA cohort were included in the analysis. LASSO Cox regression analysis was performed to construct a prognostic model based on the expression profiles of 60 herb-related genes. According to the optimal value of λ, a 20-gene signature was identified (Figure 3A,3B). The HRS was calculated as follows: risk score = −0.0182 × BCL2 − 0.0333 × PTGS2 − 0.1058 × CASP9 − 0.0042 × OXA1L + 0.0108 × CHRM1 − 0.0123 × CASP8 + 0.1561 × SLC6A3 − 0.0466 × ICAM1 − 0.1319 × ADRB1 − 0.059 × AKR1B1 + 0.0026 × DIO1 + 0.1922 × IL12B + 0.0024 × MMP9 + 0.2364 × NOS2 + 0.0117 × NOS3 + 0.2332 × PCYT1A + 0.0762 × PRKCB + 0.0181 × PRSS1 − 0.0354 × PTGER3 + 0.2284 × IL10.

Figure 3 Establishment of an herb-related gene signature with LASSO. (A,B) Gene markers were identified using the LASSO Cox regression model. (C) Patients were divided into low- and high-risk groups according to the median risk score. (D) Difference in survival analysis between high- and low-risk groups. (E) Scatter plot distribution of survival time and survival status according to the risk score of different samples. (F) ROC curve and AUC of the gene. The higher the AUC is, the greater the predictive power is. (G) Heatmap of gene expression of the signature. (H) Distribution of subgroups and clinical stages in different high- and low-risk groups. The X-axis in (C,E,G) is the risk score. AUC, area under the curve; CI, confidence interval; HR, hazard ratio; LASSO, least absolute shrinkage and selection operator; ROC, receiver operating characteristic.

Using the median risk score as the cutoff, patients in the TCGA cohort were stratified into high-risk (n=548) and low-risk (n=548) groups (Figure 3C). Kaplan-Meier survival analysis revealed that patients in the high-risk group had a significantly poorer OS than those in the low-risk group (P<0.001; Figure 3D). Differential expression patterns of the 20 herb-related genes between the two risk groups are presented in Figure 3E. The prognostic performance of the HRS was further assessed by tROC analysis (Figure 3F), yielding areas under the curve of 0.676 (95% CI: 0.593–0.760) at 1 year, 0.713 (95% CI: 0.655–0.772) at 3 years, and 0.709 (95% CI: 0.657–0.761) at 5 years. The expression distribution of individual genes within the signature is further illustrated in Figure 3G.

In addition, older age, postmenopausal status, advanced tumor node metastasis (TNM) stage, a higher mortality rate, and positive expression of estrogen receptor and HER2 were significantly associated with a higher risk score in the TCGA cohort (Figure 3H and Table 1, left).

Prognostic and correlation analysis of included herb-related genes

First, correlation analysis of 20 included herb-related genes was performed and revealed a strong correlation between them (Figure 4A). In univariate Cox regression analysis, most of the genes (18/20, 90%) were found to be differentially expressed between tumor tissues and adjacent nontumorous tissues (Figure 4B), and nine of them correlated directly with patient survival (Figure 4C). Based on survival analysis of these nine genes in the signature, NOS3, PCYT1A and SLC6A3 were associated with poor OS in the cohort from TCGA, and ADRB1, BCL2, CASP9, PTGER3, IL12B and PRSS1 were associated with better OS (all adjusted P<0.05, Figure 4D).

Figure 4 Characterization of individual genes. (A) The red color represents a positive correlation between 20 genes; the blue color represents a negative correlation. Correlations are stronger when the color is darker. (B) Differential expression of 20 genes between breast cancer samples and adjacent normal samples. (C) Univariate Cox regression analysis showed the relationship between expression of 20 genes and OS. (D) Nine genes were statistically significant by single-gene survival analysis. **, P<0.01; ***, P<0.001. CI, confidence interval; HR, hazard ratio; OS, overall survival.

Functional analysis in TCGA and METABRIC cohorts

By performing GO enrichment and KEGG pathway analyses on these herb-associated genes, we determined which biological functions and pathways are associated with the risk level. As shown in Figure 5A,5B, genes were mainly enriched in the cell cycle followed by metabolism-related pathways in GO enrichment analysis. Furthermore, the biological signaling pathways of many malignant tumors were found in KEGG analysis, including breast cancer and prostate cancer. As expected, the results indicated that several tumorigenesis-related molecular functions were enriched, such as the mTOR signaling pathway, Hippo signaling pathway and Ras signaling pathway (Figure 5C,5D). Interestingly, the TCGA and METABRIC cohort genes were obviously enriched in the endocrine resistance pathway.

Figure 5 Functional analysis of 20 herb-related genes. (A,B) GO enrichment terms with the highest significance in TCGA and METABRIC cohorts. (C,D) The most significant or shared KEGG pathways were identified in these cohorts. BP, biological process; CC, cellular component; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; METABRIC, Molecular Taxonomy of Breast Cancer International Consortium; MF, molecular function; TCGA, The Cancer Genome Atlas.

Validation of the herb-related signature

To monitor the robustness of the signature constructed from the cohort from TCGA, the GEO database (GSE20685, GSE10886) and the METABRIC cohort (GSE20685) were also categorized into high- or low-risk groups based on median values calculated according to the same formula as for TCGA. METABRIC baseline available characteristics based on risk level are displayed on the right in Table 1. Similarly, 20 herb-related genes differed in expression between the two groups (Figure 6A-6C). In addition, PCA confirmed distribution in distinct directions for the patients in the two subgroups (Figure 6D-6F), and patients in the high-risk group had a lower chance of surviving and died earlier (Figure 6G-6I). Furthermore, a meta-analysis was conducted to summarize the results of multivariate analysis of gene signatures across all cohorts. According to Figure 6J, patients with a higher HRS had a worse prognosis than those with a lower HRS (overall HR =2.51; 95% CI: 2.12–2.98; P<0.001).

Figure 6 Validation of the gene signature in external validation sets. An expression heatmap of the 20 herb-related genes in the METABRIC study between the high- and low-IRS groups (A), GSE20685 (B) and GSE10886 (C) cohorts. PCA of patients in the high-risk group and low-risk group in the METABRIC (D), GSE20685 (E) and GSE10886 (F) cohorts. Analysis of the Kaplan-Meier curves for the OS of patients in the METABRIC high-risk and low-risk groups (G) and GSE20685 (H) and GSE10886 (I) cohorts. (J) The training set from TCGA and three external validation sets were analyzed in a meta-analysis. CI, confidence interval; HR, hazard ratio; IRS, immunoreactive score; METABRIC, Molecular Taxonomy of Breast Cancer International Consortium; PCA, principal component analysis; TCGA, The Cancer Genome Atlas.

The gene signature serves as a valuable marker for response to immunotherapy and prediction of chemotherapy drugs

We identified immunotherapy targets and assessed response to immunotherapy in patients with high and low HRS. There was a negative correlation between expression of most immune cell markers and the risk score. Correlation analysis showed that B-cell expression (r=−0.26; P<0.001), T-cell CD4⁺ expression (r=−0.30; P<0.001), T-cell CD8⁺ expression (r=−0.23; P<0.001), neutrophil expression (r=−0.15; P<0.001), macrophage expression (r=0.11; P<0.001) and myeloid dendritic cell expression (r=−0.16; P<0.001) correlated with HRS (Figure 7A). Tumor purity, immune score, ESTIMATE score, and stromal score were also analyzed in the high- and low-risk groups. Immune, ESTIMATE, and stromal scores were higher in the low-risk group, which suggests that immunotherapy may be more effective in these patients. A higher tumor purity score further confirmed the poor prognosis of the high-risk group (Figure 7B-7E). Detailed data for Figure 7A-7E can be found in table available at https://cdn.amegroups.cn/static/public/tcr-2025-1259-2.xlsx.

Figure 7 Prediction of chemotherapy drugs and response to immunotherapy of high- and low-HRS patients. (A) Analysis of the correlation between expression of six immune cell types and risk scores. (B) Tumor purity, (C) immune score, (D) ESTIMATE score, and (E) stromal score are shown in boxplots. (F) Proportion of immune-responsive and nonimmune-responsive populations in the high- and low-HRS groups. (G,H) Graph showing differences in survival between the HRS-high and HRS-low groups in nonimmune-responsive (G) and immune-responsive (H) patients. (I,J) A bubble plot illustrates the negative correlation between AUC values and HRS for three CTRP-derived compounds (A) and two PRISM-derived compounds (B). Graphs with longer lines indicate a stronger negative correlation, indicating greater sensitivity to drugs. Boxplots show the results of differential drug response analysis for the three CTRP-derived compounds (I) and the three PRISM-derived compounds (J). It is important to note that lower values on the y-axis of boxplots indicate greater drug sensitivity. **, P<0.01; ***, P<0.001; ****, P<0.0001. AUC, area under the curve; CTRP, Cancer Therapeutics Response Portal; ESTIMATE, Estimation of Stromal and Immune cells in Malignant Tumor tissues using Expression data; HRS, herb-related risk score; NR, non-responder; PRISM, Profiling Relative Inhibition Simultaneously in Mixtures; R, responder.

In Figure 7F, the predictability of immunotherapy response is shown for patients in the HRS-high and HRS-low groups. It was predicted that 65% of patients in the high-risk group would respond to immunotherapy but that 82% of patients in the low-risk group would. According to these data, we speculate that immunotherapy may have a greater effect on patients in the low-risk group. Survival analysis was also performed between the high- and low-risk groups of patients who responded to immunotherapy versus those who did not. There was a relatively good prognosis for patients with low HRS, both among responders and nonresponders (Figure 7G,7H).

As we defined the specificity of two groups of patients with different HRS. We would also like to further understand whether new discoveries can identify drugs that have different effects on patients with different risks. Drug candidates with different drug sensitivity were screened in two different drug response databases (CTRP and PRISM). First, a Spearman correlation analysis was conducted using AUC values and HRS in order to identify compounds and drugs with negative correlation coefficients (Spearman r<−0.25 for CTRP or −0.30 for PRISM). In order to identify compounds with lower estimated AUC values in the high-HRS group (log2-FC >0.05), differences in drug response were compared between the high-HRS (highest decile) and low-HRS (lowest decile) groups. The lower the AUC, the more sensitive the drug is. Three CTRP-derived compounds (sirolimus, PI-103, and panobinostat) and two PRISM-derived compounds (napabucasin and temsirolimus) were found to be potentially sensitive to the high-risk group (Figure 7I for CTRP and Figure 7J for PRISM).

Establish and validate a nomogram for breast cancer patients

We conducted univariate and multivariate Cox regression analyses to determine whether HRS was an independent prognostic predictor of OS. According to univariate analysis, HRS was significantly associated with OS in both the TCGA and the METABRIC cohorts (HR =2.03, 95% CI: 1.25–3.30, P=0.004; HR =1.82, 95% CI: 1.59–2.10, P<0.001, respectively) (Figure 8A,8B, top). In the multivariate analysis, HRS remained an independent predictor of OS even after accounting for other confounding factors (TCGA: HR =1.77, 95% CI: 1.08–2.90, P=0.02; METABRIC: HR =1.55, 95% CI: 1.34–1.79, P<0.001) (Figure 8A,8B, bottom).

Figure 8 Clinicopathological features in univariate and multivariate Cox regression analyses of OS in TCGA training (A) and METABRIC validation (B) cohorts. CI, confidence interval; ER, estrogen receptor; HER2, human epidermal growth factor receptor 2; HRS, herb-related risk score; IDC, invasive ductal carcinoma; METABRIC, Molecular Taxonomy of Breast Cancer International Consortium; OS, overall survival; PR, progesterone receptor; TCGA, The Cancer Genome Atlas.

Since nomograms can directly predict a patient’s prognosis, we developed a nomogram based on a number of independent predictors, including age, ER status, PR status, HER2 status, tumor stage and HRS (Figure 9A). In addition, the calibration plots for the 3- and 5-year OS probabilities indicated that the predictions and actual observations were in good agreement (Figure 9B,9C). Moreover, to understand the clinical benefits of the nomogram, we conducted a decision curve analysis and compared it with previous reported models related to ferroptosis, pyroptosis and cuproptosis (22-24). Compared with previous models, the herb-related model had a higher clinical net benefit rate in predicting 3- and 5-year OS (Figure 9D,9E).

Figure 9 Combination of the HRS signature and clinicopathological features improves survival prediction. (A) A nomogram was constructed to predict the 3- and 5-year survival rates of breast cancer patients. (B,C) Calibration analysis indicated a high degree of accuracy for predicting survival at 3 and 5 years. (D,E) Herb-related gene signatures and ferroptosis-related, pyroptosis-related and cuproptosis-related gene signatures were compared at 3 and 5 years with decision curve analysis. ER, estrogen receptor; HER2, human epidermal growth factor receptor 2; HRS, herb-related risk score; IHC, immunohistochemistry; OS, overall survival; PR, progesterone receptor.

In summary, these above results demonstrated that this nomogram based on HRS had excellent clinical risk predictive ability in breast cancer patients.

Discussion

Complementary and alternative medicine (CAM) is recognized in the United States and in many parts of the world and is widely used for the treatment of disease and health management (25). According to reports, one in two Americans use alternative medicine, and 629 million people go to alternative medicine doctors for health problems each year (26). In Asian countries such as China, Mongolia, Korea, and Thailand, traditional medicine has a long history of making great contributions to the health of patients. Despite being widely used in these countries, CAM approaches, including TCM, have not been subjected to rigorous clinical trials. They are considered more empirical than scientific. An English-language literature of CAM found that out of more than 1,000 citations, only 17 had statistically significant randomized trials (27). TCM has a long history of treating breast tumors. Early Chinese medical literature described different types of breast tumors and discussed their clinical symptoms, pathophysiological changes, and prognosis (28-31). The most commonly used term for breast cancer in ancient Chinese medical records is “breast rock” in the Yellow Emperor’s Classic of Internal Medicine (written around 250 B.C.), which provides the first clinical description of breast cancer. Approximately 10 years after diagnosis, patients with breast cancer can expect progression, metastasis, and death. Although there are ancient medical records and some records of treatment success that are not verifiable as evidence, the lack of rigorous evidence-based medicine makes the results of TCM for breast cancer hard to believe. We cannot ignore the plants or herbs mentioned in these ancient prescriptions. Approximately 80% of the world’s population still relies on herbal medicine as their primary source of therapy, according to a study conducted by the World Health Organization (32). However, advances in phytopharmacology and related disciplines such as analytical chemistry have allowed us to explore the mechanisms behind this treatment with more evidence. Several databases have been developed to give us more access to the potential targets behind these TCM herbs with greater precision. For example, HERB is one such high-throughput database (http://herb.ac.cn/) (33). Therefore, we believe that re-analyzing the herbs in these ancient prescriptions using modern medical and bioinformatics approaches could provide meaningful insights. Our team has extensive experience in breast cancer treatment, and we believe that research in the diagnosis and treatment of breast cancer is of the greatest benefit to patients. Therefore, we wanted to discover a model similar to the 21-gene test to help doctors treat patients with precision.

In spite of this, previous studies have suffered from numerous unavoidable deficiencies. As a first point, many gene signatures do not have sufficient validation groups to confirm their predictive ability. Second, some gene signatures are based on relatively independent cellular mechanisms, and this lacks a holistic perspective, for example, ferroptosis (34), pyroptosis (35) and the epithelial–mesenchymal transition (36). Third, previous reports stop at the analysis of patient survival prognosis once the prognostic model is well established and rarely elaborate on chemotherapy and immunotherapy options.

To address the abovementioned shortcomings, we implemented the following improvements and explorations. First, we used three validation sets across multiple independent cohorts to confirm that the gene signature we built has very high predictive power. Second, TCM’s treatment philosophy emphasizes holistic and comprehensive treatments for maintaining the health of the entire body. Therefore, the candidate genes selected may differ considerably from previous independent mechanism gene signatures. We collected a large number of TCM prescriptions with rigorous criteria and performed high-throughput data analysis of the herbs documented in the prescriptions. Third, our model is capable of predicting not only the prognosis of breast cancer patients but also the response of these patients to immunotherapy based on their risk profile. In addition, the choice of new chemotherapeutic agents offered provides some suggestions for clinicians.

In our study, we screened effective TCM prescriptions available from 1990 to 2022 and identified 15 herbs with the highest frequency. Potential targets of these herbs were selected and analyzed through HERB, a high-throughput database, and 60 candidate genes related to herbs were obtained. Candidate genes were screened, and a gene signature was established by the LASSO algorithm. Next, the gene signature was tested using a cohort from TCGA to confirm whether it can predict the prognosis of breast cancer patients. Subsequently, the prognostic value of this gene signature was evaluated in three validation sets across multiple independent cohorts. The risk score derived from the herb-related gene signature is called the HRS in our study. We also performed a detailed analysis of these screened genes, including their expression and prognosis in breast cancer patients. The search for genetic relationships has been proven through many methods. As a result, we analyzed the clinicopathological characteristics of patients in different risk groups to predict their response to immunotherapy. The gene signature also played an important role in the prediction of new potential chemotherapy agents based on tumor heterogeneity in patients of different risk groups. Based on our signature, we found different outcomes in high- and low-risk patients with use of sirolimus, PI-103, panobinostat, napabucasin and temsirolimus. Sirolimus, also known as rapamycin, is a macrocyclic lactone antibiotic produced by Streptomyces hygroscopicus, isolated from the soil of the Vai Atari region in Rapa Nui (37). Studies have shown that sirolimus increases expression of A2M mRNA (38). In murine T-cell lymphoma 27658642, PI-103 inhibits PI3K-AKT signaling and induces apoptosis. Panobinostat is an oral deacetylase (DAC) inhibitor approved by the Food and Drug Administration (FDA) for treatment of multiple myeloma. Panobinostat inhibits cell growth by increasing reactive oxygen metabolism (39). Panobinostat also enhances Tax transcription in freshly isolated patient T cells (40). Napabucasin has been shown to inhibit cell growth in a variety of tumors (41,42) and temsirolimus to cause a decrease in cell proliferation in breast cancer cells by interfering with the cell cycle (43). These drugs are not common chemotherapy agents for breast cancer, and even the effects of these drugs are not well understood. Nonetheless, we found that they have different effects in patients with different risks, which provides more evidence for exploration of new chemotherapy regiments. Finally, according to our signature, a nomogram was built based on the risk score and other clinicopathological variables to quantify risk assessment and the survival probability of breast cancer patients. In conclusion, our data suggest that low expression levels of protective genes and high expression levels of risk genes render patients with high HRS and patients in the high-risk group less sensitive to immunotherapy.

Indeed, key molecules and their pathways, clinical feature assessments, and neoadjuvant therapy in breast cancer all play significant roles in influencing treatment outcomes (44-46), we focus on the relationship between breast cancer and chemical components. Our analysis of these herb-related genes and the risk score revealed that they are associated with apoptosis and reactive oxygen species. Since these genes are screened from various plant agents, we paid special attention to the relationship between them and chemical components in our analysis and discussion.

In breast cancer, the BCL2 protein inhibits mitomycin-induced apoptosis (47). PTGS2 affects celecoxib-induced apoptosis in breast cancer (48). The mutant form of CASP9 inhibits celecoxib and results in increased apoptosis (49). Pilocarpine can increase apoptosis in fibroblasts, and CHRM1 can affect this process (50). Methylselenic acid results in increased cleavage of the CASP8 protein, leading to increased apoptosis in breast cancer cells (51). Clenbuterol results in increased activity of the ADRB1 protein, which causes increased apoptosis (52). NOS2 has an effect on nicotine-inhibited apoptosis (53). PRKCB can promote docetaxel-induced apoptosis (54). The following genes are important in production of reactive oxygen species. Cotreatment of manganese with the SLC6A3 protein increases the response to oxidative stress in HEK293 cells (55). Myricetin and cyanidin attenuate the cellular response to reactive oxygen species, leading to decreased expression of ICAM1 and MMP9 mRNA, respectively (56). OXA1L, a mitochondrial inner membrane protein, may be involved in energy metabolism on mitochondrial membranes (57). We are currently unable to directly connect these genes with more or fewer pathways. In summary, the biological functions associated with tumor apoptosis and reactive oxygen species production of the novel gene signature in breast cancer still need further investigation.

Undeniably, our study has some limitations. First, as this study lacks validation of laboratory and clinical data, assessment of prognostic value and potential clinical application of herb-related gene markers needs to be validated in larger prospective trials. Second, detailed and standardized information on adjuvant therapies (e.g., specific chemotherapeutic regimens, duration of endocrine therapy, or use of targeted agents) was not uniformly available for the TCGA and validation cohorts. The absence of this critical confounding factor means that our prognostic model captures a composite signal reflecting both the tumor’s intrinsic aggressiveness and its response to the heterogeneous background of real-world treatments. Meanwhile, while the primary validation in the large METABRIC cohort was well-powered, the smaller sample sizes in the GEO cohorts (GSE20685 and GSE10886) result in greater uncertainty in the performance estimates, despite the consistent direction of the observed effects. Consequently, while the HRS demonstrates robust prognostic value, its independent contribution beyond established treatments cannot be precisely quantified in this analysis. This fundamental limitation underscores the necessity and defines the primary objective of our future work: a prospective validation study in a well-characterized cohort with meticulously curated treatment data, which will allow for direct adjustment of therapeutic modalities and solidify the clinical utility of the signature. Finally, while we did not conduct an internal validation (e.g., bootstrapping) on the TCGA development set to quantify the optimism of the model, the consistent and strong performance of the HRS across the METABRIC and GEO datasets, which are distinct in their patient populations and technical platforms, serves as a rigorous and ultimate test of model stability and effectively demonstrates its generalizability beyond the development cohort. Our preliminary data found that these herbs are related to the immune system, and it is our next task to continue to explore this relationship.

Conclusions

A novel herb-related gene signature from 221 TCM prescriptions was established. The prognostic value was validated using multiple external datasets. Using this model, we scored the prognostic risk of breast cancer patients and divided them into high- and low-risk groups. Comparing these groups, we found that they have different sensitivities to the immune response: a low HRS is associated with features of an immunoreactive tumor microenvironment and potential sensitivity of high-HRS tumors to agents like sirolimus. The establishment of this model contributes to precision treatment of breast cancer.

Acknowledgments

We thank Ms. Xin Yue for her help in the previous literature data statistics.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1259/rc

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1259/prf

Funding: This study was financially supported by the postdoctoral research start-up funding (Post Doc No. 333291 to H.Y.) from China and Postdoctoral Science Foundation of Chongqing (grants No. CSTB2023NSCQ-BHX0023).

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

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

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

References

1. Wyld L, Reed MWR, Morgan J, et al. Bridging the age gap in breast cancer. Impacts of omission of breast cancer surgery in older women with oestrogen receptor positive early breast cancer. A risk stratified analysis of survival outcomes and quality of life. Eur J Cancer 2021;142:48-62. [Crossref] [PubMed]
2. Siegel RL, Miller KD, Jemal A. Cancer statistics, 2020. CA Cancer J Clin 2020;70:7-30. [Crossref] [PubMed]
3. Sorlie T, Tibshirani R, Parker J, et al. Repeated observation of breast tumor subtypes in independent gene expression data sets. Proc Natl Acad Sci U S A 2003;100:8418-23. [Crossref] [PubMed]
4. Lipson EJ, Forde PM, Hammers HJ, et al. Antagonists of PD-1 and PD-L1 in Cancer Treatment. Semin Oncol 2015;42:587-600. [Crossref] [PubMed]
5. Emens LA, Ascierto PA, Darcy PK, et al. Cancer immunotherapy: Opportunities and challenges in the rapidly evolving clinical landscape. Eur J Cancer 2017;81:116-29. [Crossref] [PubMed]
6. Emens LA. Breast Cancer Immunotherapy: Facts and Hopes. Clin Cancer Res 2018;24:511-20. [Crossref] [PubMed]
7. Czechowski T, Branigan C, Rae A, et al. Artemisia annua L. plants lacking Bornyl diPhosphate Synthase reallocate carbon from monoterpenes to sesquiterpenes except artemisinin. Front Plant Sci 2022;13:1000819. [Crossref] [PubMed]
8. Jiao J, Yu J, Ji H, et al. Synthesis of macromolecular Astragalus polysaccharide-nano selenium complex and the inhibitory effects on HepG2 cells. Int J Biol Macromol 2022;211:481-9. [Crossref] [PubMed]
9. Qin L, Huang D, Huang J, et al. Integrated Analysis and Finding Reveal Anti-Liver Cancer Targets and Mechanisms of Pachyman (Poria cocos Polysaccharides). Front Pharmacol 2021;12:742349. [Crossref] [PubMed]
10. Huang F, Pang J, Xu L, et al. Hedyotis diffusa injection induces ferroptosis via the Bax/Bcl2/VDAC2/3 axis in lung adenocarcinoma. Phytomedicine 2022;104:154319. [Crossref] [PubMed]
11. An EK, Zhang W, Kwak M, et al. Polysaccharides from Astragalus membranaceus elicit T cell immunity by activation of human peripheral blood dendritic cells. Int J Biol Macromol 2022;223:370-7. [Crossref] [PubMed]
12. Chu MP, Wang D, Zhang YY, et al. Pachyman treatment improves CD4+CD25+ Treg counts and serum interleukin 4 and interferon γ levels in a mouse model of Kawasaki disease. Mol Med Rep 2012;5:1237-40. [Crossref] [PubMed]
13. Kuo YJ, Lin JP, Hsiao YT, et al. Ethanol Extract of Hedyotis diffusa Willd Affects Immune Responses in Normal Balb/c Mice In Vivo. In Vivo 2015;29:453-60.
14. Han H, Qian C, Zong G, et al. Systemic pharmacological verification of Salvia miltiorrhiza-Ginseng Chinese herb pair in inhibiting spontaneous breast cancer metastasis. Biomed Pharmacother 2022;156:113897. [Crossref] [PubMed]
15. Paik S, Shak S, Tang G, et al. A multigene assay to predict recurrence of tamoxifen-treated, node-negative breast cancer. N Engl J Med 2004;351:2817-26. [Crossref] [PubMed]
16. Finak G, Bertos N, Pepin F, et al. Stromal gene expression predicts clinical outcome in breast cancer. Nat Med 2008;14:518-27. [Crossref] [PubMed]
17. Miao YR, Zhang Q, Lei Q, et al. ImmuCellAI: A Unique Method for Comprehensive T-Cell Subsets Abundance Prediction and its Application in Cancer Immunotherapy. Adv Sci (Weinh) 2020;7:1902880. [Crossref] [PubMed]
18. Li T, Fan J, Wang B, et al. TIMER: A Web Server for Comprehensive Analysis of Tumor-Infiltrating Immune Cells. Cancer Res 2017;77:e108-10. [Crossref] [PubMed]
19. Yoshihara K, Shahmoradgoli M, Martínez E, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun 2013;4:2612. [Crossref] [PubMed]
20. Peng Y, Yu H, Jin Y, et al. Construction and Validation of an Immune Infiltration-Related Gene Signature for the Prediction of Prognosis and Therapeutic Response in Breast Cancer. Front Immunol 2021;12:666137. [Crossref] [PubMed]
21. Yang C, Huang X, Li Y, et al. Prognosis and personalized treatment prediction in TP53-mutant hepatocellular carcinoma: an in silico strategy towards precision oncology. Brief Bioinform 2021;22:bbaa164. [Crossref] [PubMed]
22. Zhang H, Yu X, Yang J, et al. Comprehensive analysis of pyroptotic gene prognostic signatures associated with tumor immune microenvironment and genomic mutation in breast cancer. Front Immunol 2022;13:933779. [Crossref] [PubMed]
23. Zhang Y, Liang Y, Wang Y, et al. A novel ferroptosis related gene signature for overall survival prediction and immune infiltration in patients with breast cancer. Int J Oncol 2022;61:148. [Crossref] [PubMed]
24. Song S, Zhang M, Xie P, et al. Comprehensive analysis of cuproptosis-related genes and tumor microenvironment infiltration characterization in breast cancer. Front Immunol 2022;13:978909. [Crossref] [PubMed]
25. Su D, Li L. Trends in the use of complementary and alternative medicine in the United States: 2002-2007. J Health Care Poor Underserved 2011;22:296-310. [Crossref] [PubMed]
26. Eisenberg DM, Kessler RC, Foster C, et al. Unconventional medicine in the United States. Prevalence, costs, and patterns of use. N Engl J Med 1993;328:246-52. [Crossref] [PubMed]
27. Jacobson JS, Workman SB, Kronenberg F. Research on complementary/alternative medicine for patients with breast cancer: a review of the biomedical literature. J Clin Oncol 2000;18:668-83. [Crossref] [PubMed]
28. Zhang D. Treatment of Cancer by Integrated Chinese-Western Medicine. 1st edition. Boulder: Blue Poppy Press; 1989.
29. Wu JN. Ling Shu, Or, the Spiritual Pivot. Washington: Univ of Hawaii Pr; 1993.
30. Fu QZ. Fu Qing-Zhu's Gynecology. Boulder: Blue Poppy Press; 1992.
31. Zhu DX. The Heart & Essence of Dan-Xi's Methods of Treatment: A Translation of Zhu Dan-Xi's Zhi Fa Xin Yao. 1st edition. Boulder: Blue Poppy Press; 1993.
32. Cohen I, Tagliaferri M, Tripathy D. Traditional Chinese medicine in the treatment of breast cancer. Semin Oncol 2002;29:563-74. [Crossref] [PubMed]
33. Fang S, Dong L, Liu L, et al. HERB: a high-throughput experiment- and reference-guided database of traditional Chinese medicine. Nucleic Acids Res 2021;49:D1197-206. [Crossref] [PubMed]
34. Ren X, Du H, Cheng W, et al. Construction of a ferroptosis-related eight gene signature for predicting the prognosis and immune infiltration of thyroid cancer. Front Endocrinol (Lausanne) 2022;13:997873. [Crossref] [PubMed]
35. Ling Y, Wang Y, Cao C, et al. Molecular subtypes identified by pyroptosis-related genes are associated with tumor microenvironment cell infiltration in colon cancer. Aging (Albany NY) 2022;14:9020-36. [Crossref] [PubMed]
36. Venugopal A, Michalczyk A, Khasraw M, et al. EMT Molecular Signatures of Pancreatic Neuroendocrine Neoplasms. Int J Mol Sci 2022;23:13645. [Crossref] [PubMed]
37. Yakupoglu YK, Kahan BD. Sirolimus: a current perspective. Exp Clin Transplant 2003;1:8-18.
38. Cui Y, Huang Q, Auman JT, et al. Genomic-derived markers for early detection of calcineurin inhibitor immunosuppressant-mediated nephrotoxicity. Toxicol Sci 2011;124:23-34. [Crossref] [PubMed]
39. Yaseen A, Chen S, Hock S, et al. Resveratrol sensitizes acute myelogenous leukemia cells to histone deacetylase inhibitors through reactive oxygen species-mediated activation of the extrinsic apoptotic pathway. Mol Pharmacol 2012;82:1030-41. [Crossref] [PubMed]
40. Schnell AP, Kohrt S, Aristodemou A, et al. HDAC inhibitors Panobinostat and Romidepsin enhance tax transcription in HTLV-1-infected cell lines and freshly isolated patients' T-cells. Front Immunol 2022;13:978800. [Crossref] [PubMed]
41. Li Y, Rogoff HA, Keates S, et al. Suppression of cancer relapse and metastasis by inhibiting cancer stemness. Proc Natl Acad Sci U S A 2015;112:1839-44. [Crossref] [PubMed]
42. Zhang Y, Jin Z, Zhou H, et al. Suppression of prostate cancer progression by cancer cell stemness inhibitor napabucasin. Cancer Med 2016;5:1251-8. [Crossref] [PubMed]
43. Fung AS, Wu L, Tannock IF. Concurrent and sequential administration of chemotherapy and the Mammalian target of rapamycin inhibitor temsirolimus in human cancer cells and xenografts. Clin Cancer Res 2009;15:5389-95. [Crossref] [PubMed]
44. Ren Z, Li Y, Yang X, et al. LSM4 as a potential prognostic indicator and therapeutic target in triple-negative breast cancer progression. Gland Surg 2025;14:1990-2004. [Crossref] [PubMed]
45. Lin GL, Zhang M, Du XJ, et al. Influencing factors of axillary lymph node metastasis and prognosis in patients with T2 breast cancer. Gland Surg 2025;14:1949-61. [Crossref] [PubMed]
46. Fennelly S, Shah B, Issac M, et al. Patient outcomes and clinician perspectives following one year of ad hoc implementation of neoadjuvant endocrine therapy in early breast cancer. Transl Breast Cancer Res 2025;6:34. [Crossref] [PubMed]
47. Pirnia F, Schneider E, Betticher DC, et al. Mitomycin C induces apoptosis and caspase-8 and -9 processing through a caspase-3 and Fas-independent pathway. Cell Death Differ 2002;9:905-14. [Crossref] [PubMed]
48. Sugimoto T, Bartholomeusz C, Tari AM, et al. Adenovirus type 5 E1A-induced apoptosis in COX-2-overexpressing breast cancer cells. Breast Cancer Res 2007;9:R41. [Crossref] [PubMed]
49. Jendrossek V, Handrick R, Belka C. Celecoxib activates a novel mitochondrial apoptosis signaling pathway. FASEB J 2003;17:1547-9. [Crossref] [PubMed]
50. Reina S, Sterin-Borda L, Passafaro D, et al. Muscarinic cholinoceptor activation by pilocarpine triggers apoptosis in human skin fibroblast cells. J Cell Physiol 2010;222:640-7. [Crossref] [PubMed]
51. Li Z, Carrier L, Rowan BG. Methylseleninic acid synergizes with tamoxifen to induce caspase-mediated apoptosis in breast cancer cells. Mol Cancer Ther 2008;7:3056-63. [Crossref] [PubMed]
52. Burniston JG, Tan LB, Goldspink DF. beta2-Adrenergic receptor stimulation in vivo induces apoptosis in the rat heart and soleus muscle. J Appl Physiol (1985) 2005;98:1379-86. [Crossref] [PubMed]
53. Argentin G, Cicchetti R. Evidence for the role of nitric oxide in antiapoptotic and genotoxic effect of nicotine on human gingival fibroblasts. Apoptosis 2006;11:1887-97. [Crossref] [PubMed]
54. Hung CH, Chan SH, Chu PM, et al. Docetaxel Facilitates Endothelial Dysfunction through Oxidative Stress via Modulation of Protein Kinase C Beta: The Protective Effects of Sotrastaurin. Toxicol Sci 2015;145:59-67. [Crossref] [PubMed]
55. Roth JA, Eichhorn M. Down-regulation of LRRK2 in control and DAT transfected HEK cells increases manganese-induced oxidative stress and cell toxicity. Neurotoxicology 2013;37:100-7. [Crossref] [PubMed]
56. Yi L, Chen CY, Jin X, et al. Differential suppression of intracellular reactive oxygen species-mediated signaling pathway in vascular endothelial cells by several subclasses of flavonoids. Biochimie 2012;94:2035-44. [Crossref] [PubMed]
57. Itoh Y, Andréll J, Choi A, et al. Mechanism of membrane-tethered mitochondrial protein synthesis. Science 2021;371:846-9. [Crossref] [PubMed]