Machine Learning-Based Identification of Tamoxifen Resistance-Associated Genes and Their Application in Prognostic Modeling of Breast Cancer

Summary

We developed a machine learning–derived six-gene tamoxifen resistance signature that stratifies breast cancer patients by survival risk and may support personalized prognostic assessment and therapeutic decision-making in clinical practice.

Abstract

Tamoxifen is a key endocrine therapy for estrogen receptor-positive (ER⁺) breast cancer, but acquired resistance limits long-term efficacy. The molecular mechanisms remain complex, and predictive biomarkers are lacking. Gene expression data related to tamoxifen resistance were obtained from GEO (GSE67916), and differentially expressed genes (DEGs) were identified using the limma algorithm. Functional enrichment analyses (GO and KEGG) revealed involvement in immune processes, antiviral responses, endocytosis, lysosome pathways, and estrogen signaling. Three machine learning algorithms (LASSO, SVM-RFE, and RF) identified six hub genes (CAMK1D, CHAC1, KIAA0513, MED13, NDRG1, STXBP5). A prognostic risk model based on these genes was constructed using TCGA-BRCA data, effectively stratifying patients into high- and low-risk groups with significantly different overall survival. The model demonstrated good predictive accuracy (AUC = 0.70) and stable performance in time-dependent ROC analyses, validated in an independent cohort. This study provides a robust tamoxifen resistance–related gene signature and a multigene prognostic model, offering novel insights into resistance mechanisms and potential guidance for individualized prognosis and therapy in ER⁺ breast cancer.

Introduction

Breast cancer is the most commonly diagnosed cancer among women globally and remains a major contributor to cancer-associated deaths¹. Based on molecular profiling, approximately 60-70% of patients present with the estrogen receptor-positive (ER⁺) subtype, which generally responds well to initial endocrine therapy. Among various endocrine agents, tamoxifen—a selective estrogen receptor modulator (SERM)—has long been the standard of care for adjuvant and metastatic treatment of ER⁺ breast cancer, significantly improving survival outcomes in large-scale randomized trials².

Despite its established clinical efficacy, a substantial proportion of patients eventually develop acquired resistance, leading to disease recurrence and progression³. Consequently, tamoxifen resistance represents a critical bottleneck in the long-term management of breast cancer. Previous studies indicate that the mechanisms driving resistance are highly heterogeneous, involving aberrant ER signaling, activation of alternative growth factor pathways (e.g., PI3K/AKT, MAPK), dysregulation of apoptosis, and metabolic reprogramming⁴.

Recently, the tumor microenvironment (TME), particularly the immune microenvironment, has garnered increasing attention for its role in therapeutic response and resistance. Immune cell infiltration and inflammatory signaling are closely linked to breast cancer progression and sensitivity to treatment. Emerging evidence suggests that immune-related processes may modulate cellular stress responses and promote immune evasion, thereby contributing to endocrine resistance⁵. However, the specific immune landscape associated with tamoxifen resistance remains to be fully elucidated.

With the advent of high-throughput sequencing and bioinformatics, analyzing transcriptomic data offers a powerful approach to decoding the molecular basis of drug resistance. Compared to single-gene biomarkers, multi-gene signatures capture tumor heterogeneity more effectively, offering superior stability and accuracy in prognosis. Machine learning algorithms, such as LASSO regression, Support Vector Machines (SVM), and Random Forest (RF), have become essential tools for identifying robust biomarkers and constructing predictive models in precision medicine⁶. Although several gene signatures have been proposed^7,8,9, robust validation is often lacking.

While previous studies have explored tamoxifen resistance genes, few have systematically integrated multiple machine learning strategies to build a robust prognostic model validated in independent cohorts. In this study, we integrated transcriptomic data from GEO and TCGA to screen for tamoxifen resistance-associated genes. By applying an ensemble of machine learning algorithms, we identified core resistance genes and constructed a prognostic risk model. We further evaluated the model's association with the immune microenvironment and validated its predictive value in an independent cohort, aiming to provide new theoretical evidence for personalized treatment strategies.

Protocol

All data used in this study were obtained from publicly accessible databases (TCGA, GEO, and METABRIC). No human participants or animals were involved; therefore, institutional review board approval and informed consent were not required.

Data acquisition and preprocessing

Gene expression data related to tamoxifen resistance were retrieved from the Gene Expression Omnibus (GEO) database¹⁰. The GSE67916 dataset (Affymetrix Human Genome U133 Plus 2.0 Array) includes 18 breast cancer cell samples: 8 untreated control samples and 10 tamoxifen-resistant samples generated through long-term drug exposure. RNA-sequencing data and corresponding clinical follow-up information for the Breast Invasive Carcinoma cohort (TCGA-BRCA) were downloaded from The Cancer Genome Atlas (TCGA)¹¹.

Raw microarray CEL files were processed in R (version 4.4.2) using the affy package. Background correction and normalization were performed with the Robust Multi-array Average (RMA) algorithm, including log2 transformation and quantile normalization. Probe IDs were mapped to gene symbols using platform annotation files; for genes with multiple probes, the average expression value was used. For TCGA RNA-seq data, transcripts per million (TPM) values were log2-transformed [log2(TPM + 1)]. Samples with incomplete survival information or missing clinical variables were excluded.

Identification of differentially expressed genes

Differential expression between tamoxifen-resistant and control samples was assessed using the limma R package¹² with empirical Bayes moderation. Genes with |log2 fold change| > 1 and an adjusted P value < 0.05 (Benjamini–Hochberg FDR) were defined as differentially expressed genes (DEGs).

Functional enrichment analysis

Gene Ontology (GO)¹³and Kyoto Encyclopedia of Genes and Genomes (KEGG)¹⁴ analyses were performed using the clusterProfiler¹⁵ R package. GO categories included biological process (BP), cellular component (CC), and molecular function (MF). Adjusted P values < 0.05 were considered significant.

Machine learning–based feature selection

Three machine learning algorithms were applied to identify hub genes: (1) LASSO regression¹⁶ (glmnet package) with 10-fold cross-validation to select the optimal penalty parameter (lambda.min); (2) Support vector machine–recursive feature elimination (SVM-RFE)¹⁷ (e1071 package) with fivefold cross-validation to identify the minimal gene subset with lowest classification error; (3) Random forest (RF)¹⁸ (randomForest package) with 500 trees (ntree = 500); genes were ranked by MeanDecreaseGini. Genes identified by all three methods were defined as hub genes.

Construction of the prognostic risk model

A multigene prognostic risk model was constructed using TCGA-BRCA gene expression and survival data. Survival-associated genes were screened using univariate Cox regression, followed by multivariate Cox regression to develop the final signature. The risk score formula was calculated as: Risk score = (0.01297 × CAMK1D) + (0.03021 × CHAC1) + (0.02018 × KIAA0513) + (0.00647 × MED13) + (0.00108 × NDRG1) + (0.04551 × STXBP5). Patients were stratified into high- and low-risk groups based on the median risk score.

Evaluation and validation of the prognostic model

Overall survival differences between groups were assessed with Kaplan–Meier analysis and the log-rank test. Predictive performance was evaluated using ROC curves (pROC package) and time-dependent ROC analysis (timeROC package). A nomogram integrating risk scores and clinical variables was constructed using the rms package. Calibration curves assessed agreement between predicted and observed survival probabilities. External validation was conducted in the independent Molecular Taxonomy of Breast Cancer International Consortium (METABRIC) cohort¹⁹ using the same formula and cutoff.

Immune infiltration analysis

Immune cell infiltration was estimated using CIBERSORT²⁰ with 1,000 permutations based on TCGA data. Samples with P < 0.05 were included. Differences in immune cell composition between high- and low-risk groups were assessed with the Wilcoxon rank-sum test, and correlations between hub gene expression and immune cell abundance were evaluated using Spearman's rank correlation.

Statistical analysis

All analyses were performed in R. Continuous variables were compared using the Wilcoxon rank-sum test, and categorical variables using the chi-square test. Two-sided P < 0.05 was considered statistically significant.

Results

Identification of differentially expressed genes associated with tamoxifen resistance

A total of 854 differentially expressed genes (DEGs) were identified between tamoxifen-resistant and control samples, including 556 upregulated and 298 downregulated genes. The distribution of DEGs showed a predominance of upregulated genes in resistant samples, suggesting extensive transcriptional activation associated with tamoxifen resistance (Figure 1).

<...

Discussion

Tamoxifen resistance remains a pivotal obstacle in the treatment of ER⁺ breast cancer. While mechanisms such as ER mutations are well known, the systemic molecular adaptations to long-term therapy remain less well understood. In this study, we integrated transcriptomics and machine learning to identify a robust six-gene signature (CAMK1D, CHAC1, KIAA0513, MED13, NDRG1, STXBP5) that predicts both tamoxifen resistance and patient prognosis.

Our functional analysis revealed that resist...

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
Affy package	Bioconductor	-	Microarray data processing
clusterProfiler package	Bioconductor	v4.4.2	GO and KEGG enrichment analysis
CIBERSORT algorithm	Newman et al.	-	Immune cell infiltration estimation
e1071 package (SVM-RFE)	CRAN	-	Support vector machine recursive feature elimination
Gene Expression Omnibus (GEO)	NCBI	GSE67916	Tamoxifen resistance dataset
glmnet package	CRAN / Bioconductor	-	LASSO regression
limma package	Bioconductor	-	Differential expression analysis
METABRIC dataset	cBioPortal / Curtis et al.	-	Validation cohort
pROC package	CRAN	-	ROC curve analysis
R software	R Core Team	v4.4.2	Statistical computing environment
Random forest package	CRAN	-	Random Forest algorithm
rms package	CRAN	-	Nomogram construction
TCGA-BRCA cohort	The Cancer Genome Atlas (TCGA)	-	Training cohort RNA-seq data
timeROC package	CRAN	-	Time-dependent ROC analysis