Identification and diagnostic model construction of cuproptosis-related biomarkers in bone metastatic prostate cancer

Introduction

According to GLOBOCAN 2022, prostate cancer is the most frequently diagnosed cancer among men in 112 countries worldwide (1), and the emergence of bone metastases marks a crucial shift toward systemic disease and irreversible clinical decline. Approximately 90% of the patients with advanced prostate cancer develop bone metastases [bone metastatic prostate cancer (BMPC)]. Patients with BMPC often experience severe bone pain, pathological fractures, and functional impairment (2). This significantly reduces patient quality of life, shortens survival time, and leads to substantial medical costs (3). However, diagnosing bone metastases in clinical practice remains challenging (4). Current approaches depend largely on imaging studies and the detection of specific biomarkers. Early identification continues to be difficult because of the scarcity of highly sensitive and specific molecular indicators that substantially affect clinical decision-making and patient outcomes. Consequently, identification of novel and stable biomarkers for early detection is critical.

In recent years, the field of cell death-related research has continuously made breakthrough progress. In 2022, a new copper ion-dependent programmed cell death mode—cuproptosis—was systematically revealed and defined (5). This death mode is different from classic programmed death types such as apoptosis, necroptosis, pyroptosis, and ferroptosis, and is directly triggered by excessive accumulation of intracellular copper ions. Its core molecular mechanism is: copper ions (Cu⁺) accumulated in mitochondria can specifically bind to acylated metabolic enzymes in the tricarboxylic acid cycle (such as DLAT), induce abnormal aggregation of key proteins and cause instability of iron-sulfur cluster proteins, ultimately inducing lethal proteotoxic stress and mitochondrial respiratory chain failure (6). This discovery has opened up new research perspectives and intervention directions for tumor treatment, especially for tumor types sensitive to abnormal copper metabolism.

There is a close and specific intrinsic association between prostate cancer and copper metabolism disorders. Multiple clinical studies have confirmed that the copper levels in the serum and tumor tissues of prostate cancer patients are significantly elevated, and they show a significant positive correlation with high Gleason scores, highly aggressive phenotypes, and poor overall survival rates (7). The high dependence of tumor cells on copper ions to maintain their biological processes such as rapid proliferation, angiogenesis, invasion and metastasis, and drug resistance formation is defined as “cuproplasia” (6). The clarification of the mechanism of cuproptosis transforms the “addictive” characteristic of tumor cells to copper into a metabolically targetable vulnerability, that is, through artificial induction of cuproptosis, specifically eliminating copper-dependent tumor cells (8). In the specific pathological microenvironment of prostate cancer bone metastasis, microenvironmental characteristics such as hypoxia, extracellular matrix remodeling, and osteoclast-osteoblast interactions can deeply regulate copper homeostasis and the sensitivity of cells to cuproptosis, endowing this process with a high degree of complexity and potential targetability (9).

However, the specific contribution of cuproptosis-related genes (CRGs) to BMPC progression and their potential as diagnostic markers remain largely unexplored. To address this gap, this study aimed to systematically explore the biomarkers related to cuproptosis in BMPC using multi-dimensional bioinformatics analyses. We screened the differentially expressed genes (DEGs) related to cuproptosis and evaluated their potential clinical applications. By identifying biomarkers associated with cuproptosis-related differentially expressed genes (CRDEGs), we aimed to construct a diagnostic model for BMPC, which was expected to effectively assist in the diagnosis of bone metastasis in patients before imaging evidence appears, and to preliminarily explore the specific regulatory mechanisms of CRDEGs in the bone metastatic microenvironment. The results of this study will provide a certain theoretical basis for the precise diagnosis of BMPC and the formulation of individualized treatment strategies, and promote the substantial transformation of cuproptosis, an emerging research field, from basic research to clinical application. We present this article in accordance with the TRIPOD reporting checklist (available at https://tau.amegroups.com/article/view/10.21037/tau-2026-1-0036/rc).

Methods

Data download

To identify DEGs, RNA-seq data were retrieved from the Gene Expression Omnibus (GEO) database utilizing the R package “GEOquery” (10) (Version 2.70.0). Specifically, the BMPC datasets GSE29650 (11) and GSE32269 (12) were downloaded from a public repository (https://www.ncbi.nlm.nih.gov/geo/) (13). Both cohorts comprised human prostate cancer bone metastasis tissues, with detailed demographic and clinical characteristics summarized in Table 1. Furthermore, we incorporated a comparative analysis of BMPC and primary localized prostate cancer (PLPC) groups.

To compile a comprehensive list of candidate genes, we integrated resources from the GeneCards database and peer-reviewed literature. Specifically, a total of 24 protein-coding CRGs were initially retrieved from the GeneCards database (14) (https://www.genecards.org/) using “Cuproptosis” as the keyword filter. Complementing this dataset, an additional 48 relevant CRGs were manually curated from PubMed-indexed publications (12). Following the removal of duplicate entries using Venn diagram analysis, a final non-redundant set comprising 64 distinct CRGs was established, and its detailed characteristics are summarized in Table S1.

To ensure data consistency across platforms, batch effects in the GSE29650 and GSE32269 datasets were corrected using the ComBat function implemented in the R package “sva” (v3.50.0) (15). This yielded a unified combined GEO dataset comprising 59 BMPC and 22 PLPC samples. Subsequently, quality control procedures were conducted via the R package “limma” (version 3.58.1) (16), encompassing probe annotation, standardization, and normalization steps. To assess the impact of batch correction, principal component analysis (17) was performed on gene expression matrices from both raw and corrected datasets.

CRDEGs associated with bone metastasis of prostate cancer

Based on the grouping of the combined GEO dataset, differential expression analysis was conducted between the BMPC and PLPC groups utilizing the R package “limma” (Version 3.58.1). DEGs were identified based on a threshold of |logFC| >0 and an adjusted P value <0.05. The Benjamini-Hochberg method was used to correct P values. Specifically, transcripts with logFC >0 were categorized as upregulated, whereas those with logFC <0 were classified as downregulated. The resulting data visualization, including volcano plots, was generated using the R package “ggplot2” (Version 3.4.4).

To identify CRDEGs specifically associated with BMPC, the intersection of DEGs derived from the combined GEO dataset and curated CRGs were determined using Venn diagram analysis. Subsequently, the top 20 CRDEGs were visualized using a heatmap generated with the R package “pheatmap” (version 1.0.12).

Functional and pathway enrichment analysis

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed to gain insight into the biological functions and pathways linked to CRDEGs. These analyses were implemented with the “clusterProfiler” package (version 4.10.0) (18) in R. Statistical significance for enriched terms was assessed using the Benjamini-Hochberg method, with an adjusted P value (q value) cutoff of <0.05 (19,20).

Gene set enrichment analysis (GSEA)

To further elucidate the biological pathways underlying BMPC progression, GSEA (21) was conducted utilizing the R package “clusterProfiler” (Version 4.10.0). Analyses were performed on ranked expression profiles derived from the combined GEO dataset. Specifically, the default parameters were adjusted as follows: a random seed of 2020 was set to ensure reproducibility, and gene set size constraints were defined with a minimum of 10 genes and a maximum of 500 genes per set. Significantly enriched gene sets were identified based on the Molecular Signatures Database collection “c2.all.v2023.2.Hs.symbols”, applying strict screening criteria of P<0.05 and false discovery rate value (q value) <0.25.

Gene set variation analysis (GSVA)

To evaluate whether distinct biological pathways exhibited differential activity between the BMPC and PLPC groups, this study retrieved the C2 curated pathway gene set “c2.cp.v2023.2.Hs.symbols.gmt” from the Molecular Signatures Database (22). Subsequently, the R package “GSVA” (Version 1.50.0) was employed to perform (23) on the entire combined GEO dataset. The resulting pathway activity scores were compared between the two groups, and statistical significance was defined as P<0.05.

Construction of diagnostic model for BMPC

A multi-step feature selection strategy was implemented to establish a robust diagnostic model for BMPC. Initially, logistic regression analysis was employed to identify the CRDEGs significantly associated with the BMPC versus PLPC groups (P<0.05). Subsequently, a support vector machine (SVM) algorithm (24) was utilized to refine the gene pool based on accuracy and error rate criteria. Following this, the least absolute shrinkage and selection operator (LASSO) regression analysis was conducted using the R package “glmnet” (25) (Version 4.1-10) with parameters set.seed [500] and family = ’binomial’. The resulting diagnostic model, visualized using trajectory plots, yielded a set of key genes critical for BMPC identification. A risk score was constructed by multiplying the expression level of each gene by its corresponding coefficient from the LASSO regression and summing the products. The risk score was calculated as follows (26).

RiskScore=∑iCoefficient(genei)∗mRNAExpression(genei)[1]

Validation of the diagnostic model for BMPC

A comprehensive evaluation framework was implemented to validate the diagnostic efficacy of the established model rigorously. Specifically, a nomogram was constructed utilizing the R package “rms” (Version 8.0-0) to visualize the prognostic impact of key genes. Calibration curves were generated to quantitatively assess the concordance between predicted probabilities and actual clinical outcomes. Furthermore, the R packages “ggDCA” and “pROC” were employed to perform decision curve analysis (DCA) and receiver operating characteristic (ROC) curve analysis, respectively. The diagnostic accuracy was quantified by calculating the area under the curve (AUC), where values closer to 1.0 indicate superior discriminatory power.

Protein-protein interaction (PPI) analysis

To elucidate the functional associations among the key genes, a PPI network was constructed utilizing the STRING database (27) (Version 11.5) with a confidence threshold of ≥0.150. The resulting network was visualized using the Cytoscape software. GeneMANIA (28) was used to predict functionally similar genes and expand the interaction network of key genes.

Construction of key gene regulatory network

To systematically investigate the regulatory mechanisms governing the key genes, a multilayered network was constructed by integrating transcriptional, post-transcriptional, and pharmacological interactions. First, transcription factors (TFs) were retrieved from the ChIPBase database (29) (http://rna.sysu.edu.cn/chipbase/) to analyze their regulatory effects on target genes, forming the mRNA-TF regulatory axis. Subsequently, miRNAs associated with these genes were identified using the TarBase database (30) (http://www.microrna.gr/tarbase) to construct an mRNA-miRNA interaction module. Furthermore, the target RNA-binding proteins (RBPs) were predicted using the StarBase v3.0 database (31) (https://starbase.sysu.edu.cn/) to complete the mRNA-RBP regulatory network. Finally, the Comparative Toxicogenomics Database (32) (https://ctdbase.org/) was used to predict direct and indirect drug targets, thereby establishing a comprehensive mRNA-drug regulatory framework. All resulting networks were visualized using the Cytoscape software (33).

Differential expression analysis, correlation analysis, and ROC curve analysis

To evaluate the diagnostic potential of key genes, a comparative grouping analysis was conducted between the BMPC and PLPC groups based on their expression profiles. The diagnostic efficacy was quantitatively assessed using the R package “pROC” (Version 1.18.5) to generate ROC curves and calculate the AUC. Generally, an AUC value closer to 1.0 indicates a superior discriminatory power. Spearman’s correlation analysis was performed on the expression levels of key genes within the combined GEO dataset. The resulting correlation matrix was visualized as a heatmap utilizing the R package “pheatmap” (Version 1.0.12). Correlation coefficients were categorized based on their absolute values: <0.3 denoted weak/non-correlation, 0.3–<0.5 indicated weak correlation, 0.5–<0.8 represented moderate correlation, and ≥0.8 indicated strong correlation.

Immune infiltration analysis

To systematically evaluate the tumor microenvironment, single-sample GSEA was used to quantify the relative abundance of immune cell infiltrates (34). The infiltration levels of distinct immune cell types, namely activated CD8⁺ T cells, activated dendritic cells, gamma-delta T cells, natural killer cells, and regulatory T cells, were estimated for each sample. This process generated an integrated immune infiltration matrix for the combined GEO datasets. Subsequently, differential expression patterns between low- and high-risk groups were visualized using the R package “ggplot2” (Version 3.4.4).

To further investigate the interplay between immune infiltration and biological processes, Spearman’s correlation analysis was conducted to quantify the associations among distinct immune cell types. The resulting correlation matrix was visualized as a heat map using the R package “pheatmap” (version 1.0.12). Furthermore, the relationships between key genes and immune cell infiltration were analyzed. A correlation bubble plot was generated via the R package “ggplot2” (Version 3.4.4) to illustrate the strength and direction of these gene-immune interactions.

Statistical analysis

All statistical analyses were conducted using R software (Version 4.2.2). For continuous variables, the independent Student’s t-test was used to compare normally distributed data between the two groups, whereas the Mann-Whitney U test (Wilcoxon rank-sum test) was used for non-normally distributed comparisons. The Kruskal-Wallis test was used to compare three or more groups. Spearman’s correlation analysis was used to determine associations between distinct molecular entities. Unless otherwise specified, all statistical P values were calculated as two-sided tests, with a threshold of P<0.05 considered statistically significant.

Ethical statement

The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. This study utilized publicly available datasets from the GEO database (GSE29650 and GSE32269). As the data were derived from previously published studies and are fully anonymized, no additional ethical approval or patient consent was required for this secondary analysis.

Results

The technology roadmap is presented in Figure 1.

Figure 1 Technology roadmap. CRDEGs, cuproptosis-related differentially expressed genes; CRGs, cuproptosis-related genes; DEGs, differentially expressed genes; GO, Gene Ontology; GSEA, gene set enrichment analysis; GSVA, gene set variation analysis; KEGG, Kyoto Encyclopedia of Genes and Genomes; LASSO, least absolute shrinkage and selection operator; PPI, protein-protein interaction; RBP, RNA-binding protein; ROC, receiver operating characteristic; SVM, support vector machine; TF, transcription factor.

Data collection and correction

Distribution boxplots (Figure 2A,2B) and principal component analysis plots (Figure 2C,2D) demonstrated that batch effects were effectively removed from the BMPC dataset.

Figure 2 Batch effects removal of GSE29650 and GSE32269. (A) Box plot of the GEO datasets distribution before batch effect removal. (B) Box plot of the combined GEO dataset distribution after batch effect removal. (C) PCA plot of the GEO datasets before batch effect removal. (D) PCA plot of the combined GEO dataset after batch effect removal. Red color represents BMPC dataset GSE29650, and blue color represents BMPC dataset GSE32269. BMPC, bone metastatic prostate cancer; GEO, Gene Expression Omnibus; PCA, principal component analysis.

CRDEGs associated with BMPC

Differential expression analysis of the combined GEO dataset revealed significant transcriptomic differences between the BMPC and PLPC groups. A total of 3,401 DEGs met the threshold criteria (|logFC| >0 and P<0.05), comprising 1,744 upregulated and 1,657 downregulated genes.

To visualize these findings, a volcano plot was generated to illustrate DEG distribution (Figure 3A). Furthermore, the intersection of these DEGs with CRGs was analyzed and 24 CRDEGs were identified, including MT1M, HSPD1, MT1E, MT1F, MT1H, GPX4, MT1G, GLS, MT1X, CYP1A1, AQP2, MT3, SLC25A3, CDKN2A, LIPT1, MTF1, ARID1A, MT2A, AQP1, PDHA1, SLC27A5, MT4, TAX1BP1, and DAXX (Figure 3B). Based on these intersection results, the CRDEG expression profiles were further characterized across different sample groups within the combined GEO dataset. The top 20 CRDEGs were visualized using hierarchical clustering via the R package “pheatmap” (Figure 3C).

Figure 3 Differential gene expression analysis between BMPC group and PLPC group in the combined GEO dataset. (A) Volcano plot of differential gene expression analysis between BMPC group and PLPC group in the combined GEO dataset. (B) Venn diagram of DEGs and CRGs in the combined GEO dataset. (C) Heatmap of CRDEGs in the combined GEO dataset. Red indicates BMPC group, and blue indicates PLPC group. In the heatmap, red represents high expression, while blue represents low expression. BMPC, bone metastatic prostate cancer; CRDEGs, cuproptosis-related differentially expressed genes; CRGs, cuproptosis-related genes; DEGs, differentially expressed genes; GEO, Gene Expression Omnibus; PLPC, primary localized prostate cancer.

GO and KEGG pathway enrichment analysis

GO and KEGG enrichment analyses were performed to elucidate the biological significance of the 24 CRDEGs. The detailed statistical results are summarized in Table 2. The analysis revealed that these genes were significantly enriched in specific biological processes and molecular functions associated with metal ion homeostasis and transport. Specifically, they are primarily involved in the cellular responses to Cu and Zn ions, including detoxification and stress responses. Functional annotations highlighted activities related to polyol transmembrane transporters, water channel activity, and carbohydrate transport across membranes. At the pathway level, CRDEGs were notably enriched in metabolic processes, such as mineral absorption and lipoic acid metabolism, as well as renal functions, such as proximal tubule bicarbonate reclamation. To visualize these findings, bar plots were generated to illustrate the distribution of enriched terms (Figure 4A). Additionally, network diagrams were constructed to depict the intricate relationships among the biological processes, molecular functions, and KEGG pathways (Figure 4B-4D).

Figure 4 GO and KEGG enrichment analysis for CRDEGs. (A) Bar plots showing the GO and KEGG enrichment results of CRDEGs: BP, MF and KEGG pathways, with the horizontal axis representing GO terms and KEGG entries. (B-D) Network diagrams of the GO and KEGG enrichment results for CRDEGs: BP (B), MF (C) and KEGG (D). The selection criteria for GO and KEGG enrichment analyses were P value <0.05 and FDR (q value) <0.05. Red nodes represent terms, blue nodes represent molecules, and edges indicate the relationships between terms and molecules. BP, biological process; CRDEGs, cuproptosis-related differentially expressed genes; FDR, false discovery rate; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, molecular function.

GSEA

To further investigate the biological pathways underlying the differential expression, GSEA was conducted based on logFC values derived from the combined GEO dataset. GSEA results are visualized in a mountain plot (Figure 5A), with detailed statistical data summarized in Table 3. The analysis revealed a significant enrichment of specific biological functions and signaling cascades. Notably, the dataset exhibited marked activation or suppression of pathways, such as Croonquist Il6 Deprivation Dn (Figure 5B), Tang senescence Tp53 target Dn (Figure 5C), and the Wp Pi3k Akt signaling pathway (Figure 5D).

Figure 5 GSEA for combined datasets. (A) Mountain plots illustrating three GSEA-derived biological functions of the combined GEO dataset. (B-D) GSEA results showing that the combined GEO dataset was significantly enriched in pathways/functions Croonquist Il6 Deprivation Dn (B), Tang senescence Tp53 Targets Dn (C) and Wp Pi3k Akt signaling pathway (D). The selection criteria for GSEA were P value <0.05 and an FDR (q value) <0.05. FDR, false discovery rate; GEO, Gene Expression Omnibus; GSEA, gene set enrichment analysis; NES, normalized enrichment score.

GSVA

To further characterize the functional differences at the pathway level, GSVA was performed on the entire gene set of the combined GEO dataset (Table 4). Pathways with statistically significant enrichment (P<0.05) were ranked according to the absolute values of logFC, and the top 20 pathways were selected for detailed visualization. The differential expression patterns of these 20 pathways between the BMPC and PLPC groups were characterized using a heat map (Figure 6A). Furthermore, a group comparison plot was generated to provide a comprehensive overview of the relative expression levels across samples (Figure 6B).

Figure 6 Gene set variation analysis for the combined GEO dataset. (A,B) Heatmap (A) and group‑comparison plot (B) of GSVA results between the BMPC and PLPC groups in the combined GEO dataset. ***, P<0.001, denoting a highly statistically significant difference. Red denotes the BMPC group, and blue denotes the PLPC group. The selection criterion for GSVA was P value <0.05. In the heatmap, blue represents low enrichment and red represents high enrichment. BMPC, bone metastatic prostate cancer; GEO, Gene Expression Omnibus; GSVA, gene set variation analysis; KEGG, Kyoto Encyclopedia of Genes and Genomes; PLPC, primary localized prostate cancer.

The GSVA analysis demonstrated that pathways P27 cell cycle G1 S, EBV EBNA3C to P27 cell cycle G1 S, RB pathway, gastric cancer network 1, ATR pathway, regulation of sister chromatid separation at the metaphase anaphase transition, phosphorylation of the APC/C, inhibition of the proteolytic activity of APC/C required for the onset of anaphase by mitotic spindle checkpoint components, cohesin dissociation in anaphase, PRC2 methylates histones and DNA, constitutive signaling by overexpressed ERBB2, DNA replication licensing, MCM pathway, unwinding of DNA, pre-IC formation, miRNA targets in extracellular matrix (ECM) and membrane receptors, miR-509-3p alteration of YAP1 ECM axis, omega-3 omega-6 fatty acid synthesis, bile acid biosynthesis, and response to metal ions showed biologically meaningful differences between the BMPC and PLPC groups (P<0.05).

Construction of diagnostic model for BMPC

Univariate logistic regression was performed based on the 24 CRDEGs, revealing 23 significant associations with bone metastasis (P<0.05; see Table S2). An SVM model was developed using these 23 CRDEGs and SVM algorithms to identify the gene set that yielded the lowest error rate (Figure 7A) and the highest classification accuracy (Figure 7B). The model achieved optimal accuracy when all 23 CRDEGs were included. These genes were: AQP2, SLC27A5, GLS, MT3, GPX4, CYP1A1, PDHA1, ARID1A, MT1H, LIPT1, MT1G, HSPD1, MT1E, CDKN2A, SLC25A3, MT2A, TAX1BP1, AQP1, MT1X, MT1M, MT1F, MTF1, and MT4. Subsequently, a LASSO regression model designated the BMPC diagnostic model was established based on the same 23 CRDEGs. Model selection was visualized using the LASSO variable trajectory diagram (Figure 7C) and the LASSO regression model diagram (Figure 7D). The final LASSO model retained the following 13 CRDEGs as key predictors: GPX4, HSPD1, MT1E, MT1F, ARID1A, LIPT1, GLS, MT1H, AQP1, SLC27A5, CDKN2A, MT3, and CYP1A1.

Figure 7 Construction of diagnostic model for bone metastatic prostate cancer. (A) Number of genes yielding the lowest error rate obtained with the SVM algorithm. (B) Number of genes achieving the highest accuracy obtained with the SVM algorithm. (C) LASSO variable trajectory diagram for CRDEGs in the combined GEO dataset. (D) LASSO regression model diagram. CRDEGs, cuproptosis-related differentially expressed genes; GEO, Gene Expression Omnibus; LASSO, least absolute shrinkage and selection operator; SVM, support vector machine.

Validation of the diagnostic model for BMPC

A nomogram was developed based on the 13 key genes included in the BMPC diagnostic model (Figure 8A). Among these, the expression level of the key gene SLC27A5 contributed substantially more to the diagnostic performance of the model than any other variable. The calibration curve for the BMPC model (Figure 8B) exhibited only a minor deviation from the ideal diagonal line, indicating good agreement between the predicted and observed probabilities. DCA (Figure 8C) shows that, across a clinically relevant threshold range, the model provided a higher net benefit compared to both the “treat all” and “treat none” strategies, reflecting its clinical utility. Furthermore, the ROC curve (Figure 8D) demonstrated that the logistic regression model derived from the combined GEO dataset achieved strong diagnostic accuracy.

Figure 8 Validation analysis of a diagnostic model for BMPC. (A) Nomogram of the combined GEO dataset showing the key genes in the BMPC diagnostic model. (B,C) Calibration curve (B) and DCA plot (C) of the BMPC diagnostic model based on the key genes from the combined GEO dataset. (D) ROC analysis of the linear predictors from the logistic regression model in the combined GEO dataset. In the DCA plot, the vertical axis represents net benefit and the horizontal axis represents the probability threshold (or threshold probability). AUC, area under the curve (when the AUC exceeds 0.9, it indicates high accuracy); BMPC, bone metastatic prostate cancer; CI, confidence interval; DCA, decision curve analysis; FPR, false positive rate; GEO, Gene Expression Omnibus; ROC, receiver operating characteristic; TPR, true positive rate.

GSEA

The combined GEO dataset was stratified into high- and low-risk groups using the median LASSO risk score from the BMPC model as the cut-off value. Differential expression analysis identified 560 DEGs that met the criteria |logFC| >0 and P<0.05. Among these, 276 genes were upregulated (logFC >0, P<0.05), and 284 were downregulated (logFC <0, P<0.05). GSEA was then performed across all genes in the combined GEO dataset using logFC values derived from the comparison between the high- and low-risk groups (Figure 9A). The analysis revealed significant enrichment in biologically relevant pathways and functional gene sets, including the Croonquist Il6 Deprivation Dn (Figure 9B), Foroutan Tgfb Emt Up (Figure 9C), Manalo Hypoxia Up (Figure 9D), and the Wp Pi3k Akt signaling pathway (Figure 9E).

Figure 9 GSEA for combined GEO datasets based on BMPC risk score. (A) Mountain plots showing the four biological functions identified by GSEA in the combined GEO dataset based on the prostate cancer bone metastasis risk score. (B-E) GSEA results indicating that the combined GEO dataset is significantly enriched in Croonquist Il6 Deprivation Dn (B), Foroutan Tgfb Emt Up (C), Manalo Hypoxia Up (D), and Wp Pi3k Akt signaling pathway (E). The selection criteria for GSEA were P value <0.05 and an FDR (q value) <0.05. BMPC, bone metastatic prostate cancer; FDR, false discovery rate; GEO, Gene Expression Omnibus; GSEA, gene set enrichment analysis; NES, normalized enrichment score.

GSVA

GSVA revealed that pathways mutation-caused aberrant Aβ to electron transfer in complex IV, electron transfer in complex IV, arsenic to electron transfer in complex IV, mitochondrial complex IV assembly, omega-3 omega-6 fatty acid synthesis, alpha-linolenic (omega-3) and linoleic (omega-6) acid metabolism, cholesterol biosynthesis pathway (WikiPathways), cholesterol biosynthesis (KEGG), white fat cell differentiation, ACE inhibitor pathway, scavenging by class A receptors, collagen chain trimerization, collagen biosynthesis and modifying enzymes, miRNA targets in ECM and membrane receptors, NCAM1 interactions, and anchoring fibril formation were significantly differentially enriched between the high- and low-risk BMPC groups (P<0.05). The top 20 pathways meeting the significance threshold (P<0.05) were selected and ranked in descending order based on the absolute value of their log-fold change (|logFC|). A heat map was used to visualize the differential enrichment of these 20 pathways between the two risk groups (Figure 10A). To validate these differences, the Mann-Whitney U test was applied, and the results are displayed in a group comparison plot (Figure 10B).

Figure 10 GSVA for combined GEO datasets based on BMPC risk score. (A,B) Heatmap (A) and group comparison plot (B) of the GSVA results between the high- and low-risk groups in the combined GEO dataset. ns, P≥0.05 (not statistically significant); *, P<0.05 (statistically significant); **, P<0.01 (highly statistically significant). The selection criterion for GSVA was P<0.05. In the heatmap, blue represents low enrichment and red represents high enrichment. BMPC, bone metastatic prostate cancer; GEO, Gene Expression Omnibus; GSVA, gene set variation analysis.

Differential expression analysis, correlation analysis, and genomic localization of key genes

The group comparison plot (Figure 11) shows the differences in the expression of the 13 key genes between the BMPC and PLPC groups in the combined GEO dataset. Among them, AQP1 showed significant differential expression (P<0.05); MT3, CDKN2A, LIPT1, and ARID1A exhibited highly significant differences (P<0.01); and CYP1A1, MT1H, GLS, MT1F, MT1E, HSPD1, and GPX4 were significantly differentially expressed (P<0.001).

Figure 11 Differential expression analysis of the 13 key genes. Group comparison map of the key genes in the BMPC group and PLPC group of the combined GEO datasets. ns, P≥0.05 (not statistically significant); *, P<0.05 (statistically significant); **, P<0.01 (highly statistically significant); ***, P<0.001 (extremely statistically significant). BMPC, bone metastatic prostate cancer; GEO, Gene Expression Omnibus; PLPC, primary localized prostate cancer.

A correlation heat map of the 13 key genes was generated using the combined GEO dataset (Figure 12A), which revealed strong positive correlations among several members of the MT1 family. Specifically, MT1H showed significant positive associations with MT1E, MT1F, and the other MT1 genes.

Figure 12 Correlation analysis and chromosomal localization analysis of the 13 key genes. (A) Correlation analysis among the key genes. (B) Chromosomal localization map of the key genes in the human genome. An absolute correlation coefficient (r value) below 0.3 indicates weak or no correlation; 0.3–<0.5 indicates weak correlation; 0.5–<0.8 indicates moderate correlation; and ≥0.8 indicates strong correlation.

Furthermore, we visualized the chromosomal positions of the 13 core genes using a circular ideogram (Figure 12B). Chromosome 1 harbors ARID1A; chromosome 2 contains LIPT1, GLS, and HSPD1; chromosome 7 encodes AQP1; chromosome 9 carries CDKN2A; chromosome 15 hosts CYP1A1; chromosome 16 accommodates MT1H, MT1F, MT1E, and MT3; and chromosome 19 bears GPX4 and SLC27A5.

Immune infiltration analysis

The abundances of 28 immune cell types in the BMPC samples were estimated using the single-sample GSEA algorithm. A group comparison plot (Figure 13A) revealed that four immune cell populations exhibited significant differences between the high- and low-risk BMPC groups (P<0.05): activated CD8⁺ T cells, CD56 dim natural killer cells (CD56 dim NK cells), memory B cells, and monocytes.

Figure 13 Immune infiltration analysis by ssGSEA algorithm. (A) Group comparison plot of immune cells in the low- and high-risk BMPC samples. (B,C) Correlation analysis results of immune cells infiltration abundance in the low-risk group (B) and the high-risk group (C) of BMPC samples. (D,E) Bubble plots showing the correlation between immune cells infiltration abundance and the key genes in the low-risk group (D) and the high-risk group (E) of BMPC samples. ns, P≥0.05 (not statistically significant); *, P<0.05 (statistically significant); **, P<0.01 (highly statistically significant). Interpretation of correlation coefficient (|r|) values: negligible (<0.3), weak (0.3–<0.5), moderate (0.5–<0.8), and strong (≥0.8). Blue bars/points represent the low-risk group; red bars/points represent the high-risk group. Red indicates positive correlation, and blue indicates negative correlation. The intensity of the color represents the strength of the correlation. BMPC, bone metastatic prostate cancer; MDSC, myeloid-derived suppressor cell; ssGSEA, single-sample gene set enrichment analysis.

The correlation patterns among the infiltration levels of these four immune cell types were visualized using correlation heat maps (Figure 13B,13C). In the low-risk group, most immune cell pairs showed strong negative correlations, with memory B cells and monocytes displaying the strongest significant inverse association (r=−0.384, P=0.01; Figure 13B). In contrast, the high-risk group exhibited predominantly positive correlations among immune cells, although memory B cells and monocytes still demonstrated the most pronounced negative correlation (r=−0.391; Figure 13C).

Furthermore, the associations between the 13 key genes and immune cell infiltration levels were illustrated using correlation bubble plots (Figures 13D,13E). In the low-risk group, most immune cells showed relatively strong correlations with key genes, among which ARID1A and activated CD8⁺ T cells had the most significant negative correlation (r=−0.417, P=0.007; Figure 13D). Similar strong correlations were observed in the high-risk group, with LIPT1 and CD56 dim NK cells exhibiting the strongest negative association (r=−0.628, P=0.004; Figure 13E).

Interaction network profiling of key genes

The PPI network of 13 key genes is shown in Figure 14A. Figure 14B illustrates the interaction networks and functional relationships between the 13 key genes and their homologous or co-expressed partners.

Figure 14 PPI network analysis and network construction of key genes. (A) PPI network diagram of the 13 key genes. (B) Interaction network between the 13 key genes and their predicted functionally similar genes. Circular nodes represent genes, with size reflecting the genes’ attributes and characteristics. Edges represent relationships between genes, such as interactions or functional connections; line thickness indicates the strength or importance of the association. PPI, protein-protein interaction.

Multilayer regulatory network of key genes in bone metastasis

The mRNA-TF regulatory network (Figure 15A) consisted of 12 key genes and 89 TFs; comprehensive details are listed in table online: https://cdn.amegroups.cn/static/public/tau-2026-1-0036-1.xlsx. The mRNA-miRNA regulatory network (Figure 15B) encompassed eight key genes and 188 miRNAs, with the corresponding annotations provided in table online: https://cdn.amegroups.cn/static/public/tau-2026-1-0036-2.xlsx. In the mRNA-RBP regulatory network (Figure 16A), five key genes interacted with 87 RBPs; the full information is available in table online: https://cdn.amegroups.cn/static/public/tau-2026-1-0036-3.xlsx. Finally, the mRNA-drug regulatory network (Figure 16B) linked nine key genes to 56 drugs or small-molecule compounds, as detailed in table online: https://cdn.amegroups.cn/static/public/tau-2026-1-0036-4.xlsx.

Figure 15 Regulatory network of key genes (mRNA-TF, mRNA-miRNA). (A) mRNA-TF regulatory network of key genes. (B) mRNA-miRNA regulatory network of key genes. Color yellow denotes mRNA, purple denotes TF, and blue denotes miRNA. miRNA, microRNA; TF, transcription factor.

Figure 16 Regulatory network of key genes (mRNA-RBP, mRNA-drug). (A) mRNA-RBP regulatory network of key genes. (B) mRNA-drug regulatory network of key genes. Color yellow represents mRNA, green represents RBP, and pink represents drug. RBP, RNA-binding protein.

Discussion

The occurrence of bone metastasis is a key marker of the progression of prostate cancer from local lesions to advanced systemic diseases. It is an important feature of prostate cancer progression, which is also closely related to the prognosis and treatment strategies of patients. Studies have shown that the specific crosstalk between prostate cancer cells and bone tissue is a key factor in the development of bone metastasis, which further complicates treatment and prognosis (35). Current diagnostic methods for prostate cancer lack specificity for early identification and diagnosis of BMPC, which limits the accuracy of clinical decision-making (36). Therefore, it is particularly important to study the mechanism of bone metastasis in prostate cancer and develop new biomarkers and diagnostic models to improve the early identification of BMPC.

This study aimed to investigate the role of cuproptosis-related biomarkers in BMPC and develop a diagnostic model using bioinformatic approaches. By integrating gene expression profiles from publicly available databases, we identified 24 CRDEGs in BMPC. Functional enrichment analyses revealed their involvement in key biological processes and signaling pathways; based on these findings, we constructed a 13-gene diagnostic model for BMPC. The identification of these candidate genes and the development of a diagnostic model offer novel insights into the molecular mechanisms underlying prostate cancer bone metastasis and hold promise for improving early detection strategies. Furthermore, these results may serve as a foundation to identify potential therapeutic targets and advance personalized treatment approaches for BMPC.

We identified 24 CRDEGs, including MT1M, HSPD1, and MT1E. These genes play pivotal roles in copper ion metabolism and modulate tumor cell survival and migratory processes, ultimately leading to pathological changes, such as bone metastasis, in prostate cancer. Metallothionein 1M (MT1M) belongs to the metallothionein family, a group of small, cysteine-rich cytoplasmic proteins that serve as primary intracellular buffers for transition metals such as copper and zinc. It can bind excess copper ions to maintain metal homeostasis and protect cells from heavy metal toxicity, DNA damage, and oxidative stress. In prostate epithelial cells, MT1M expression exhibits high sensitivity, being upregulated even at extremely low copper concentrations (37). In this study, MT1M showed lower expression changes in BMPC compared with PLPC, suggesting that this gene plays a tumor-suppressive role in the occurrence of prostate cancer bone colonization. Reviewing the literature, the potential mechanism by which low MT1M expression promotes prostate cancer bone metastasis may unfold through multidimensional synergistic actions. The core initiating factor in this process is hypermethylation of its gene promoter region. This epigenetic silencing mechanism has been confirmed in various cancers, including hepatocellular carcinoma and thyroid cancer (38,39), suggesting that loss of MT1M expression is an early and stable marker of prostate cancer progression, laying the foundation for subsequent mechanistic abnormalities. Within the specialized microenvironment of bone metastasis in prostate cancer cells, MT1M, acting as an intracellular copper ion buffer, exhibits a direct reduction in expression that leads to intracellular metal ion imbalance. Although uncontrolled copper ion concentrations may trigger oxidative stress, cancer cells can evade immune system-mediated copper-induced death surveillance through compensatory mechanisms such as upregulating copper transporters. This enables them to gain a survival advantage within the metal-stressed microenvironment of bone tissue (40). Concurrently, reduced MT1M expression releases inhibition of the NF-κB signaling pathway. As a potent pro-oncogenic TF, abnormally activated NF-κB initiates the expression of multiple pro-inflammatory and pro-oncogenic genes. This not only leads to uncontrolled tumor cell proliferation and suppressed apoptosis but also promotes the secretion of pro-inflammatory factors, constructing an inflammatory microenvironment conducive to metastasis. This mechanism has been validated in hepatocellular carcinoma research (41). Additionally, low MT1M levels promote epithelial-mesenchymal transition by activating the SOD1/PI3K-AKT pathway, enhancing cancer cells’ migration and invasion capabilities. This facilitates their detachment into the bloodstream or bone marrow cavity and successful skeletal metastasis (42). Simultaneously, reduced MT1M expression reshapes the tumor immune microenvironment (TIME) by inducing the accumulation of immunosuppressive cells such as tumor-associated macrophages. This suppresses antitumor immune responses, helping metastatic cells evade immune surveillance and further driving the bone metastasis process (43). Clinically, MT1M expression levels serve as a potential biomarker for assessing prostate cancer aggressiveness and metastatic risk. Restoring MT1M expression (e.g., via demethylating agents) or targeting downstream signaling pathways like NF-κB may represent novel therapeutic strategies for metastatic prostate cancer.

The Heat Shock Protein Family D Member 1 (HSPD1) (also known as Heat Shock Protein 60, HSP60) is another key participant in BMPC identified in this study. It belongs to the chaperone protein family and primarily participates in the proper folding of mitochondrial proteins. Under cellular stress (such as heat shock and oxidative stress), it is upregulated to maintain protein homeostasis. HSPD1 does not directly bind to copper ions or participate in the regulation of copper metabolism, but is responsible for the mitochondrial unfolded protein response. The main mechanism of cuproptosis involves the binding of copper to lipoylated proteins in mitochondria, leading to the aggregation of toxic proteins and the loss of iron-sulfur cluster proteins. This triggers mitochondrial stress, which in turn activates unfolded protein response through pathways such as Activating Transcription Factor 5, mobilizing key chaperone proteins like HSPD1 to repair damaged proteins and alleviate copper-induced proteotoxic stress. This mechanism has been supported in relevant literature (44). In this study, compared to PLPC, HSPD1 was upregulated in BMPC, which aligns with most of the previous literature, indicating it may be a marker of poor prognosis in prostate cancer patients, including the occurrence of bone metastasis.

Glutathione peroxidase 4 (GPX4), a central regulator of lipid peroxidation, emerged as a key gene identified in this study. A recent meta-analysis assessing the oncogenic relevance of GPX4 demonstrated that its elevated expression correlated with worse overall survival, increased lymph node involvement, metastatic potential, and more advanced clinical stages (45). Research has indicated that GPX4 expression is significantly higher in castration-resistant prostate cancer than in hormone-sensitive prostate cancer, and its mechanisms of action include supporting tumor metabolic reprogramming, participating in immune evasion, and ferroptosis (46). This highlights the delicate interplay between the pro-survival and cell death pathways in the context of copper metabolism. Dysregulation of GPX4 may lead to increased oxidative stress and subsequent tumor progression, making it a promising candidate for targeted therapies aimed at inducing ferroptosis in cancer cells. Our findings demonstrated that GPX4 expression was significantly lower in prostate cancer bone metastases than in primary prostate tumors, which contradicts the findings of other studies linking high GPX4 levels to poor prognosis. Therefore, we suggest that the role of GPX4 in the metastatic process may vary depending on the tumor type and environment. Future research could validate the specific function of GPX4 in the metastatic process through experiments involving gene knockdown or overexpression and by integrating multi-omics data (e.g., transcriptomics, proteomics, and metabolomics) to comprehensively analyze the involvement of GPX4 in the metastatic process.

GO and KEGG enrichment analyses revealed that CRDEGs were predominantly involved in the “cellular response to copper ion” and “detoxification of copper ion” pathways, highlighting the central role of copper homeostasis in BMPC. These pathways collectively reflect the cellular mechanisms that regulate copper uptake and distribution and are essential for maintaining mitochondrial activity, redox balance, and major oncogenic signaling programs (47,48). The cellular response to the copper ion pathway encompasses copper homeostasis, energy metabolism, and signal transduction pathways. Cu ions regulate key cancer processes by modulating multiple pathways (such as RTK, PI3K/AKT, and MAPK), mitochondrial function, and the tumor microenvironment (49,50). Disruption of copper homeostasis activates oncogenic signaling and promotes tumor progression, suggesting that targeting copper metabolism may offer therapeutic strategies for prostate cancer. The “detoxification of copper ions” pathway, primarily mediated by metallothioneins, copper chaperones, and ATP7A/ATP7B transporters, further underscores the importance of metal buffering in maintaining intracellular stability (51,52). Epidemiological and molecular evidence indicate that systemic copper overload and impaired detoxification correlate with prostate cancer risk and prognosis, although biological effects may vary depending on the efficiency of local buffering systems, such as metallothioneins (53).

Additionally, enrichment in pathways such as “mineral absorption”, “lipoic acid metabolism”, and “water channel activity” suggests broader metabolic and transport alterations accompanying copper dysregulation. Taken together, these results reinforce the idea that Cu-dependent metabolic imbalance is a key feature of BMPC and may provide opportunities for therapeutic intervention.

In this study, we constructed and validated a diagnostic model for BMPC based on a 13-gene signature. This model demonstrated excellent discriminative ability (AUC =0.97), good calibration, and, crucially, potential clinical net benefits confirmed by DCA in an independent cohort. Although direct comparison with standard clinical markers such as Prostate-Specific Antigen was not feasible due to the limitations of public datasets, the described above multidimensional internal validation provides robust support for its reliability. Future research directions include external validation in prospective cohorts with complete clinicopathological information, and integration of this molecular signature with existing clinical indicators to develop an integrated diagnostic model with greater clinical application value. In the SVM-based diagnostic model, the expression level of SLC27A5 exerted a disproportionately strong influence on the classification performance compared to other variables, suggesting its potential as a key discriminative biomarker for bone metastasis in prostate cancer. To the best of our knowledge, existing studies have not directly reported the expression or function of SLC27A5 in prostate cancer or its bone metastases, and existing studies have mainly focused on its role in hepatocellular carcinoma (54). Accordingly, this is the first study to propose the potential biological significance of SLC27A5 in BMPC and provide a new research direction for subsequent experimental validation. Fatty acid transport protein 5 (FATP5) encoded by SLC27A5 primarily mediates long-chain fatty acid uptake and participates in lipid metabolic regulation (55,56). Given that prostate cancer progression and metastasis are tightly linked to lipid reprogramming and enhanced fatty acid utilization, dysregulation of SLC27A5 may influence metabolic flux, oxidative stress, and ultimately, metastatic behavior. These findings suggest that SLC27A5 may represent a previously unrecognized metabolic node in BMPC and warrant further mechanistic and experimental validation.

In this study, we performed an integrated analysis of key genes identified in BMPC samples using the STRING and GeneMANIA platforms. The analysis revealed that these genes are closely associated with metal-ion homeostasis, antioxidant defense, cell-cycle regulation, metabolic reprogramming, and chromatin remodeling. They interact directly or indirectly through co-expression, shared protein domains, and predicted associations, suggesting a coordinated role in governing tumor initiation, bone metastasis, and remodeling of the immune microenvironment.

The bioinformatics analysis in this study revealed a significant association between risk stratification based on CRG signatures and the characteristics of the immune microenvironment in bone metastases. There were significant differences in the infiltration abundance of activated CD8⁺ T cells, CD56 dim NK cells, memory B cells, and monocytes between the high- and low-risk groups, and the correlation patterns among these immune cell subsets were completely different. These findings closely link the copper metabolism status of tumor cells with the local immune landscape, providing a new perspective for understanding the immune escape mechanism of bone metastases. Remodeling of TIME is tightly linked to the transcriptional regulation of key genes, which acts as a core molecular mechanism driving the progression of BMPC. In this study, correlation analysis between core immune cells and key genes in high- and low-risk groups revealed robust associations between most immune cells and key genes in both groups, with distinct specificities in the core regulatory targets: ARID1A was significantly negatively correlated with activated CD8⁺ T cells in the low-risk group (r=−0.417, P=0.007), while LIPT1 exhibited the most prominent negative correlation with CD56 dim NK cells in the high-risk group (r=−0.628, P=0.004). These findings uncover the specific gene-immune cell interaction patterns in BMPC with different risk grades and identify pivotal molecular targets for elucidating the remodeling mechanism of the immune microenvironment in bone metastatic foci.

The intersection of cuproptosis and TIME is particularly pronounced and clinically relevant in the context of bone metastasis, which remains a leading cause of morbidity and mortality for osteotropic malignancies (57). The bone marrow is not merely an inert structural reservoir; it serves as the primary site of hematopoiesis and functions as an immune-privileged secondary lymphoid organ harboring a highly specialized, inherently immunosuppressive microenvironment (57). When solid tumors—such as those originating in the breast, prostate, lung, or malignant melanoma—disseminate to the skeletal system, they systematically hijack local cellular machinery to establish a pre-metastatic niche (58). Since all the data used in this study were obtained from public transcriptome databases, with a limited sample size and a lack of more systematic immunomics information, we did not conduct further in-depth analysis of immune mechanisms. This design is consistent with the exploratory orientation of this study, but it also suggests that future research can combine multi-omics and experimental studies to further clarify the interaction between cuproptosis and immune regulation in the formation of bone metastasis.

There are some limitations in this study. This study conducted a preliminary exploration on the role of cuproptosis mechanism in the occurrence and development of metastatic prostate cancer and the possible targeted therapy in the later stage. The cuproptosis-related biomarkers identified in this study have not been experimentally validated, and their biological functions in prostate cancer bone metastasis remain to be confirmed using functional assays. The relatively small sample size, coupled with the absence of multicenter clinical data, may constrain the generalizability of the findings. Moreover, potential batch effects from integrating heterogeneous datasets may introduce technical variability, compromising the robustness and clinical translatability of the model. The PPI analysis and regulatory network analysis of key genes are only initial exploratory findings, which are speculative and require further experimental verification in the later stage. To address these gaps, future research should prioritize larger prospectively collected cohorts from multiple institutions and incorporate mechanistic studies to validate the biological relevance and clinical utility of the proposed biomarkers.

Conclusions

Employing an integrative bioinformatics strategy, this study identified a set of novel carcinogenesis-related biomarkers associated with bone metastasis in prostate cancer and developed a diagnostic model for prostate cancer bone metastasis with high discriminatory power (AUC =0.97). The analysis suggested that these genes are involved in regulating copper ion homeostasis and tumor metabolic reprogramming, and may have potential associations with the occurrence and development of BMPC. Although this work lacks experimental validation in vivo and in vitro, the findings provide novel hypotheses and research directions for further exploring the molecular mechanisms of cuproptosis in prostate cancer bone metastasis and its potential clinical applications.

Acknowledgments

None.

Footnote

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

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

Funding: None.

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tau.amegroups.com/article/view/10.21037/tau-2026-1-0036/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. This study utilized publicly available datasets from the GEO database (GSE29650 and GSE32269). As the data were derived from previously published studies and are fully anonymized, no additional ethical approval or patient consent was required for this secondary analysis.

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

References

1. Bray F, Laversanne M, Sung H, et al. Global cancer statistics 2022: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin 2024;74:229-63. [Crossref] [PubMed]
2. Li X, Wang Y, Fan Z, et al. High-level sustainable production of the characteristic protopanaxatriol-type saponins from Panax species in engineered Saccharomyces cerevisiae. Metab Eng 2021;66:87-97. [Crossref] [PubMed]
3. Wu YH, Gugala Z, Barry MM, et al. Optimization and Characterization of a Bone Culture Model to Study Prostate Cancer Bone Metastasis. Mol Cancer Ther 2022;21:1360-8. [Crossref] [PubMed]
4. Polavaram NS, Dutta S, Islam R, et al. Tumor- and osteoclast-derived NRP2 in prostate cancer bone metastases. Bone Res 2021;9:24. [Crossref] [PubMed]
5. Tralongo P, Ballato M, Fiorentino V, et al. Cuproptosis: A Review on Mechanisms, Role in Solid and Hematological Tumors, and Association with Viral Infections. Mediterr J Hematol Infect Dis 2025;17:e2025052. [Crossref] [PubMed]
6. Boaru DL, Leon-Oliva D, Castro-Martinez P, et al. Cuproptosis: Current insights into its multifaceted role in disease, cancer, and translational/therapeutic opportunities. Biomed Pharmacother 2025;190:118422. [Crossref] [PubMed]
7. Denoyer D, Clatworthy SA, Masaldan S, et al. Heterogeneous copper concentrations in cancerous human prostate tissues. Prostate 2015;75:1510-7. [Crossref] [PubMed]
8. Yao L, Jiang B, Xu D. Strategies to combat cancer drug resistance: focus on copper metabolism and cuproptosis. Cancer Drug Resist 2025;8:15. [Crossref] [PubMed]
9. Bakir M, Dabaliz A, Dawalibi A, et al. Microenvironmental and molecular pathways driving dormancy escape in bone metastases. International Journal of Molecular Sciences 2025;26:11893. [Crossref] [PubMed]
10. Davis S, Meltzer PS. GEOquery: a bridge between the Gene Expression Omnibus (GEO) and BioConductor. Bioinformatics 2007;23:1846-7. [Crossref] [PubMed]
11. Hörnberg E, Ylitalo EB, Crnalic S, et al. Expression of androgen receptor splice variants in prostate cancer bone metastases is associated with castration-resistance and short survival. PLoS One 2011;6:e19059. [Crossref] [PubMed]
12. Chen B, Zhou X, Yang L, et al. A Cuproptosis Activation Scoring model predicts neoplasm-immunity interactions and personalized treatments in glioma. Comput Biol Med 2022;148:105924. [Crossref] [PubMed]
13. Barrett T, Wilhite SE, Ledoux P, et al. NCBI GEO: archive for functional genomics data sets--update. Nucleic Acids Res 2013;41:D991-5. [Crossref] [PubMed]
14. Fishilevich S, Nudel R, Rappaport N, et al. GeneHancer: genome-wide integration of enhancers and target genes in GeneCards. Database (Oxford) 2017;2017:bax028. [Crossref] [PubMed]
15. Leek JT, Johnson WE, Parker HS, et al. The sva package for removing batch effects and other unwanted variation in high-throughput experiments. Bioinformatics 2012;28:882-3. [Crossref] [PubMed]
16. Ritchie ME, Phipson B, Wu D, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 2015;43:e47. [Crossref] [PubMed]
17. Ben Salem K, Ben Abdelaziz A. Principal Component Analysis (PCA). Tunis Med 2021;99:383-9. [PubMed]
18. Yu G, Wang LG, Han Y, et al. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 2012;16:284-7. [Crossref] [PubMed]
19. Mi H, Muruganujan A, Ebert D, et al. PANTHER version 14: more genomes, a new PANTHER GO-slim and improvements in enrichment analysis tools. Nucleic Acids Res 2019;47:D419-26. [Crossref] [PubMed]
20. Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res 2000;28:27-30. [Crossref] [PubMed]
21. Subramanian A, Tamayo P, Mootha VK, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A 2005;102:15545-50. [Crossref] [PubMed]
22. Liberzon A, Subramanian A, Pinchback R, et al. Molecular signatures database (MSigDB) 3.0. Bioinformatics 2011;27:1739-40. [Crossref] [PubMed]
23. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics 2013;14:7. [Crossref] [PubMed]
24. Sanz H, Valim C, Vegas E, et al. SVM-RFE: selection and visualization of the most relevant features through non-linear kernels. BMC Bioinformatics 2018;19:432. [Crossref] [PubMed]
25. Engebretsen S, Bohlin J. Statistical predictions with glmnet. Clin Epigenetics 2019;11:123. [Crossref] [PubMed]
26. Simon N, Friedman J, Hastie T, et al. Regularization Paths for Cox's Proportional Hazards Model via Coordinate Descent. J Stat Softw 2011;39:1-13. [Crossref] [PubMed]
27. Szklarczyk D, Gable AL, Lyon D, et al. STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res 2019;47:D607-13. [Crossref] [PubMed]
28. Franz M, Rodriguez H, Lopes C, et al. GeneMANIA update 2018. Nucleic Acids Res 2018;46:W60-4. [Crossref] [PubMed]
29. Zhou KR, Liu S, Sun WJ, et al. ChIPBase v2.0: decoding transcriptional regulatory networks of non-coding RNAs and protein-coding genes from ChIP-seq data. Nucleic Acids Res 2017;45:D43-50. [Crossref] [PubMed]
30. Vlachos IS, Paraskevopoulou MD, Karagkouni D, et al. DIANA-TarBase v7.0: indexing more than half a million experimentally supported miRNA:mRNA interactions. Nucleic Acids Res 2015;43:D153-9. [Crossref] [PubMed]
31. Li JH, Liu S, Zhou H, et al. starBase v2.0: decoding miRNA-ceRNA, miRNA-ncRNA and protein-RNA interaction networks from large-scale CLIP-Seq data. Nucleic Acids Res 2014;42:D92-7. [Crossref] [PubMed]
32. Grondin CJ, Davis AP, Wiegers JA, et al. Predicting molecular mechanisms, pathways, and health outcomes induced by Juul e-cigarette aerosol chemicals using the Comparative Toxicogenomics Database. Curr Res Toxicol 2021;2:272-81. [Crossref] [PubMed]
33. Shannon P, Markiel A, Ozier O, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res 2003;13:2498-504. [Crossref] [PubMed]
34. Xiao B, Liu L, Li A, et al. Identification and Verification of Immune-Related Gene Prognostic Signature Based on ssGSEA for Osteosarcoma. Front Oncol 2020;10:607622. [Crossref] [PubMed]
35. Liu X, Zhao C, Lin J, et al. Analysis of Factors Influencing Bone Metastasis in Prostate Cancer and Diagnostic Value of Serum PSA, CysC and D-D. Arch Esp Urol 2024;77:72-8. [Crossref] [PubMed]
36. Genitourinary Oncology Committee of Chinese Anti-cancer Association. Expert consensus on clinical diagnosis and treatment of bone metastases and bone-related diseases of prostate cancer (2021 edition). Zhonghua Zhong Liu Za Zhi 2021;43:1016-26. [PubMed]
37. Andrews GK. Regulation of metallothionein gene expression by oxidative stress and metal ions. Biochem Pharmacol 2000;59:95-104. [Crossref] [PubMed]
38. Ji XF, Fan YC, Gao S, et al. MT1M and MT1G promoter methylation as biomarkers for hepatocellular carcinoma. World J Gastroenterol 2014;20:4723-9. [Crossref] [PubMed]
39. Chen Y, Quan R, Bhandari A, et al. Low metallothionein 1M (MT1M) is associated with thyroid cancer cell lines progression. Am J Transl Res 2019;11:1760-70. [PubMed]
40. Krężel A, Maret W. The Functions of Metamorphic Metallothioneins in Zinc and Copper Metabolism. Int J Mol Sci 2017;18:1237. [Crossref] [PubMed]
41. Mao J, Yu H, Wang C, et al. Metallothionein MT1M is a tumor suppressor of human hepatocellular carcinomas. Carcinogenesis 2012;33:2568-77. [Crossref] [PubMed]
42. Li D, Peng W, Wu B, et al. Metallothionein MT1M Suppresses Carcinogenesis of Esophageal Carcinoma Cells through Inhibition of the Epithelial-Mesenchymal Transition and the SOD1/PI3K Axis. Mol Cells 2021;44:267-78. [Crossref] [PubMed]
43. Si M, Lang J. The roles of metallothioneins in carcinogenesis. J Hematol Oncol 2018;11:107. [Crossref] [PubMed]
44. Zhou G, Liu A, Bai J, et al. Decreased ATF5 level contributes to improved mitochondrial function in oocytes exposed to vitrification stress. Front Cell Dev Biol 2024;12:1431683. [Crossref] [PubMed]
45. Wu H, Liao X, Huang W, et al. Examining the prognostic and clinicopathological significance of GPX4 in human cancers: a meta-analysis. Free Radic Res 2025;59:239-49. [Crossref] [PubMed]
46. Fan Y, Wang Y, Dan W, et al. PRMT5-mediated arginine methylation stabilizes GPX4 to suppress ferroptosis in cancer. Nat Cell Biol 2025;27:641-53. [Crossref] [PubMed]
47. Puchkova LV, Broggini M, Polishchuk EV, et al. Silver Ions as a Tool for Understanding Different Aspects of Copper Metabolism. Nutrients 2019;11:1364. [Crossref] [PubMed]
48. Tang X, Yan Z, Miao Y, et al. Copper in cancer: from limiting nutrient to therapeutic target. Front Oncol 2023;13:1209156. [Crossref] [PubMed]
49. Feng Q, Huo C, Wang M, et al. Research progress on cuproptosis in cancer. Front Pharmacol 2024;15:1290592. [Crossref] [PubMed]
50. Ye T, Zheng JL, Song XY, et al. Research progress in cuproptosis-based cancer nanotherapy. Chinese Bulletin of Life Sciences 2025;37:147-57.
51. Bray F, Ferlay J, Soerjomataram I, et al. Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin 2018;68:394-424. [Crossref] [PubMed]
52. Sung H, Ferlay J, Siegel RL, et al. Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J Clin 2021;71:209-49. [Crossref] [PubMed]
53. Siegel RL, Miller KD, Jemal A. Cancer statistics, 2020. CA Cancer J Clin 2020;70:7-30. [Crossref] [PubMed]
54. Li X, Wang J, Guo Z, et al. Copper metabolism-related risk score identifies hepatocellular carcinoma subtypes and SLC27A5 as a potential regulator of cuproptosis. Aging (Albany NY) 2023;15:15084-113. [Crossref] [PubMed]
55. Lubiński J, Lener MR, Marciniak W, et al. Serum Essential Elements and Survival after Cancer Diagnosis. Nutrients 2023;15:2611. [Crossref] [PubMed]
56. Bobrowska-Korczak B, Skrajnowska D, Kiss AK, et al. Lipid peroxidation as a predictive biomarker of the early stage of cancer. J Biol Regul Homeost Agents 2019;33:799-810. [PubMed]
57. Hiraga T. Immune microenvironment of cancer bone metastasis. Bone 2025;191:117328. [Crossref] [PubMed]
58. Roato I. Interaction among cells of bone, immune system, and solid tumors leads to bone metastases. Clin Dev Immunol 2013;2013:315024. [Crossref] [PubMed]