HOXC4 as a machine learning-derived immune biomarker for predicting treatment and prognosis in prostate cancer patients

Introduction

Prostate cancer (PCa) represents a leading cause of cancer-related morbidity and mortality among males worldwide. The incidence rate is rising, particularly in regions with high Human Development Index (HDI), which brings a heavy burden to patients, the healthcare system and society (1). Early detection and intervention are crucial for improving prognosis, with better long-term outcomes for patients in the localized stage. However, there are still a large number of cases diagnosed in the late or high-risk stage, with extremely poor prognosis. Despite advancements in imaging techniques, biopsies, and serum biomarkers such as prostate-specific antigen (PSA), existing diagnostic and prognostic tools still have shortcomings in risk stratification, efficacy prediction, and personalized treatment guidance, especially in early diagnosis and individual clinical trajectory assessment (2).

Recently, the molecular diversity of PCa has gained recognition. Genomic and transcriptomic analyses have uncovered intricate patterns of genetic alterations and dysregulation of gene expression that contribute to the onset, progression, and resistance to therapy of the disease (3). Specifically, several genetic anomalies—including fusions involving ETS family genes, MYC amplification, and the loss or mutation of PTEN and TP53—have been linked to tumor formation and disease advancement (4). However, translating these molecular findings into clinically relevant biomarkers remains a daunting challenge. Consequently, the quest for novel biomarkers that can enhance risk stratification, inform prognostic assessments, and predict treatment outcomes has emerged as a pivotal focus in PCa research (3). Particularly, immunological elements and components of the tumor microenvironment (TME) have gained attention as essential modulators of disease behavior and as promising avenues for biomarker discovery (5). Nonetheless, the identification and validation of immune-related gene signatures that are strongly correlated with PCa outcomes and therapeutic responses remain inadequately explored.

In this regard, the HOXC4 gene, which belongs to the homeoprotein transcription factor family, has garnered increasing interest due to its established functions in morphogenesis and cellular differentiation, alongside its aberrant expression in various cancers (6,7). For instance, research has indicated that HOXC4 is upregulated and functionally pertinent in hepatocellular carcinoma and breast cancer, where it plays a role in processes like epithelial-mesenchymal transition and tumor progression (6,8). Despite these findings, the specific biological functions and clinical relevance of HOXC4 in PCa have not been comprehensively defined. Evidence indicates that HOXC4, along with other related HOX family members, is frequently overexpressed in aggressive forms of PCa and may interact with known oncogenic transcription factors, potentially modifying the PCa transcriptome and impacting disease progression (9). However, the mechanisms through which HOXC4 influences the immune microenvironment, its association with immune cell infiltration, and its prognostic or predictive significance in PCa remain largely unexamined. Due to the limitations associated with traditional biomarker discovery methods and the increasing accessibility of high-dimensional transcriptomic datasets, the incorporation of sophisticated bioinformatics and machine learning techniques has emerged as a highly effective approach in cancer research (10,11). Machine learning algorithms, especially those tailored for feature selection and predictive modeling, possess the ability to explore intricate gene expression patterns, reveal subtle correlations linked to clinical characteristics, and identify molecular attributes with significant predictive capacity (12). Notably, the utilization of machine learning in the analysis of immune-associated gene signatures in cancer has promoted the advancement of personalized prognostic tools, while also showing considerable potential in enhancing patient stratification and optimizing clinical therapeutic decision-making (5). In the present study, our goal is to identify core molecular biomarkers for PCa and explore their potential as feasible therapeutic targets. We integrated transcriptomic data from four publicly available Gene Expression Omnibus (GEO) datasets, which included a total of 258 clinical samples from PCa patients, and employed a series of machine learning algorithms to build stable and reliable predictive models. The predictive performance of the constructed models was further validated in independent patient cohorts (including an external MSKCC dataset) to assess their predictive accuracy, clinical subgroup applicability, and efficacy. To enhance the interpretability of our optimal-performing model, we integrated SHapley Additive exPlanations (SHAP) analysis, combined with CIBERSORT-based immune infiltration profiling and gene set enrichment analysis (GSEA) for functional enrichment assessment, to systematically characterize the association between the screened key genes and the immune microenvironment of PCa. This research is intended to identify hub genes related to the prognostic outcomes of PCa and clarify the clinical relevance of HOXC4 in the context of tumor immunity. Its core objective is to bridge the gap between molecular characteristic profiling and individualized patient management, and ultimately establish a theoretical basis for the rational development of HOXC4-focused strategies for cancer diagnosis, prognostic evaluation, and clinical treatment. We present this article in accordance with the TRIPOD reporting checklist (available at https://tau.amegroups.com/article/view/10.21037/tau-2026-0231/rc).

Methods

Acquisition, batch adjustment and integration of datasets

In this study, we systematically screened and retrieved 4 PCa-related gene expression datasets from the GEO public functional genomics database (https://www.ncbi.nlm.nih.gov/geo/), namely GSE32448, GSE32571, GSE46602, and GSE69223. The integrated cohort contained a total of 258 tissue samples, including 108 normal prostate tissue samples and 150 PCa tissue samples. The ComBat function was used to eliminate the batch effect between the training and testing cohorts. Principal component analysis (PCA) was performed to evaluate the consistency and comparability of gene expression profiles after batch correction.

To ensure the generalizability and clinical applicability of our findings, we also incorporated independent validation cohorts. Transcriptome data and corresponding clinical information for the The Cancer Genome Atlas-Prostate Adenocarcinoma (TCGA-PRAD) cohort were downloaded from the UCSC Xena platform. Furthermore, an independent clinical validation cohort from the MSKCC study (13) was retrieved via the cBioPortal database for prognostic validation. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Differential expression analysis to identify differentially expressed genes (DEGs)

In the training cohort, the limma R package (version 3.62.2) was used for differential expression analysis between PCa tissues and normal prostate tissues. The empirical Bayes method was adopted to screen statistically significant genes, and genes meeting the criteria of |log₂ fold change (FC)| >1 and false discovery rate (FDR) <0.05 were identified as DEGs. This screening criterion can effectively balance the sensitivity and specificity of DEG identification. Heatmaps were drawn using the ggplot2 package to visualize the results of differential expression analysis.

Analysis of Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG)

To investigate the biological functions and PCa-related mechanisms linked to the detected DEGs, GO enrichment analysis was carried out using the clusterProfiler package (version 4.14.6). This analysis covered biological processes (BP), cellular components (CC), and molecular functions (MF) (6). The corresponding bar charts and bubble charts were constructed with the RichGraph (version 1.26.6) and ggplot2 (version 3.5.2) packages (7,8). KEGG pathway enrichment analysis was implemented following the same procedure.

Machine learning-based feature selection for intersected DEGs

To identify stable diagnostic signature genes from the previously identified DEGs, we utilized three complementary machine learning-based feature selection methods. First, we performed least absolute shrinkage and selection operator (LASSO) logistic regression analysis using the glmnet R package, retaining genes with non-zero regression coefficients through 10-fold cross-validation. Next, we applied the random forest (RF) algorithm via the randomForest package, ranked genes according to their MeanDecreaseGini values, and retained the top 32 genes for further analysis. Finally, we employed the support vector machine recursive feature elimination (SVM-RFE) algorithm using the e1071 R package, selecting the gene subset that resulted in the lowest cross-validation error as the optimal gene panel. To ensure the stability of the final feature set, we selected genes that were consistently identified by all three methods for subsequent model development and functional analysis.

Model construction and feature interpretability

In this study, genes consistently identified across three machine learning algorithms were selected as the final signature gene set for constructing a diagnostic model. Using the transcriptomic profiles of these genes, we developed diagnostic prediction models based on the GEO training cohort and evaluated their performance in an independent external testing cohort. Ten widely used machine learning algorithms in tumor biomarker research were trained and systematically compared, including Partial Least Squares (PLS), Support Vector Machine with Radial Kernel (SVM-Radial), Logistic Regression, K-Nearest Neighbors (KNN), Extreme Gradient Boosting (XGBoost), Gradient Boosting Machine (GBM), Neural Networks (NeuralNet), RF, Decision Trees, and Generalized Boosted Regression Models (GlmBoost). Model performance in the testing cohort was assessed using receiver operating characteristic (ROC) curves and area under the curve (AUC) metrics. The algorithm achieving the highest AUC was selected as the optimal diagnostic classifier.

To improve the interpretability of the optimal diagnostic model, we performed SHAP analysis using three R packages: DALEX, kernelshap, and shapviz. SHAP analysis can effectively balance computational efficiency and variance reduction when processing high-dimensional gene expression data. The contribution of each feature gene to the local and global predictions of the model was visualized through three complementary methods: summary bar chart, beeswarm plot, and individual-level waterfall chart.

Although SHAP analysis provides intuitive interpretations of our model’s predictive results, we acknowledge several intrinsic limitations of this method. First, SHAP is based on the assumption of feature independence, which may not be strictly maintained for gene expression data due to the widespread occurrence of gene co-expression patterns. Second, in cases with small sample sizes, SHAP values may exhibit sensitivity to slight disturbances in input data. Third, SHAP only captures statistical correlations instead of definitive causal associations. For this reason, SHAP results were mainly utilized to gain biological insights, and all interpretations were carried out in combination with well-established biological domain knowledge.

Expression validation, prognostic assessment, and clinical subgroup analysis

The expression levels of candidate genes were compared between normal prostate and PCa tissues using the TCGA-PRAD dataset. Prognostic significance was evaluated using Kaplan-Meier survival analysis for overall survival (OS) and progression-free survival (PFS), with patients stratified according to optimal cutoff values. Genes significantly associated with survival outcomes were selected for further analysis. To assess the robustness of HOXC4 as a prognostic biomarker, external validation was performed using the MSKCC cohort (13). In addition, clinical and genomic data from cBioPortal and TCGA-PRAD were integrated to perform stratified survival analyses across age and clinical risk subgroups.

Functional enrichment analysis of signature genes

To clarify the potential biological mechanism of the screened core genes, we performed GSEA based on the transcriptomic data of the GEO training cohort. Samples were divided into high-expression and low-expression groups according to the expression level of the core gene. GSEA was performed using the clusterProfiler R package (14), and the c2.cp.kegg.v7.5.1.symbols.gmt gene set from the Molecular Signatures Database (MSigDB) was used as the reference annotation set. The significance threshold was set at P<0.05 to screen statistically significantly enriched pathways.

Immune cell infiltration and immune correlation analysis of key genes

The CIBERSORT algorithm was used to deconvolute gene expression profiles and estimate the relative proportions of 22 immune cell types based on the LM22 signature matrix with 1,000 permutations. Immune infiltration landscapes between normal prostate and PCa tissues in the GEO training cohort were compared and visualized using the ggplot2 and ggpubr packages. To explore the immunological relevance of the identified signature genes, Spearman correlation analysis was performed to assess their associations with the infiltration fractions of the 22 immune cell subsets. The resulting gene-immune cell associations were further visualized as a co-association network using the LinkET R package.

To further evaluate the immunological relevance of HOXC4, we extended the immune infiltration analysis to the TCGA-PRAD cohort. Patients were divided into HOXC4-high and HOXC4-low expression groups according to the median HOXC4 expression level, and differences in immune cell infiltration between the two groups were assessed using the Wilcoxon rank-sum test. In addition, Spearman correlation analysis was performed to examine the associations between HOXC4 expression and selected immunosuppressive cell populations, including regulatory T cells (Tregs) and M2 macrophages. We also analyzed the correlations between HOXC4 expression and representative immune checkpoint molecules, including B7-H3, CTLA-4, TIM-3, LAG-3, ICOS, CD86, TIGIT, PD-1, PD-L1, SIGLEC-15, and CD40. This analysis was designed to comprehensively assess the potential relationship between HOXC4 expression and the immunosuppressive tumor microenvironment in PCa.

Pan-cancer expression analysis of HOXC4

The TIMER3 database (http://compbio.cn/timer3/), a comprehensive public resource for systematic analysis of tumor immune cell infiltration and gene expression profiles across various tumor types, was used to explore the expression pattern of HOXC4 in different human malignant tumors. We compared the differential expression of HOXC4 between tumor tissues and paired adjacent normal tissues in the pan-cancer dataset to verify its potential as a broad-spectrum tumor biomarker.

Statistical analysis

All statistical analyses were performed using R software (version 4.5.0). Differential expression analysis was conducted using the limma package with thresholds of |log₂FC| >1 and FDR <0.05. Feature selection was performed using LASSO logistic regression, RF, and SVM-RFE. Ten machine learning models were trained and evaluated using repeated five-fold cross-validation, and model performance was assessed by ROC curves and the corresponding AUC. SHAP values were calculated using the DALEX and kernelshap packages. Group differences were assessed using Student’s t-test or the Wilcoxon rank-sum test, as appropriate. Categorical variables were compared using Fisher’s exact test. Kaplan-Meier curves were compared using the log-rank test, and Cox proportional hazards regression was used to estimate hazard ratios (HRs) and 95% confidence intervals (CIs). All statistical tests were two-sided, and P<0.05 was considered statistically significant.

Results

Dataset integration and cohort baseline characteristics

A total of 258 PCa tissue samples from 4 independent GEO datasets were included in this study. After batch effect elimination using the ComBat algorithm, PCA results showed that samples from different datasets were clustered together after correction, indicating that the batch effect was effectively removed, and the gene expression profiles of the training and testing cohorts had good consistency and comparability (Figure 1A,1B).

Figure 1 Data integration and batch effect correction. (A) PCA plot of the uncorrected GEO training cohort, showing that samples are clearly clustered according to the original dataset. (B) PCA plot of the training cohort after ComBat correction, showing overlapping clustering patterns of samples, which confirms the effective elimination of batch effects. GEO, Gene Expression Omnibus; PCA, principal component analysis.

Screening of DEGs

In the GEO training cohort, differential expression analysis between normal prostate tissues and PCa tissues identified a total of 123 DEGs, including 81 downregulated genes and 42 upregulated genes (Figure 2A,2B). Compared with normal prostate tissues, the top three most significantly downregulated genes in the PCa group were NEFH, SLC14A1, and KRT15, while the top three most significantly upregulated genes were OR51E2, AMACR, and HOXC6.

Figure 2 Profiling results of DEGs. (A) Heatmap showing the top 100 DEGs in the GEO training cohort. (B) Volcano plot of DEGs in PCa samples from the GEO training cohort (42 up-regulated DEGs and 81 down-regulated DEGs). DEG, differentially expressed gene; FC, fold change; GEO, Gene Expression Omnibus; PCa, prostate cancer.

Enrichment analysis results

GO enrichment analysis showed that the screened DEGs were mainly enriched in BP including intermediate filament organization and intermediate filament-based movement. In terms of CC, DEGs were significantly enriched in the basement membrane and collagen-containing extracellular matrix, while in terms of MF, DEGs were mainly concentrated in fatty acid binding and peroxidase activity (Figure 3A,3B). KEGG pathway enrichment analysis showed that DEGs were significantly enriched in Staphylococcus aureus infection, focal adhesion, ECM-receptor interaction, muscle contraction, and drug metabolism-cytochrome P450 pathways (Figure 3C,3D). These results indicate that DEGs play a crucial role in maintaining cell structural integrity, regulating cell adhesion, and mediating metabolic processes related to tumor progression.

Figure 3 Enrichment analysis. (A,B) GO enrichment analysis for DEGs. (C,D) KEGG pathway analysis for DEGs. BP, biological process; CC, cellular component; DEG, differentially expressed gene; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, molecular function.

Feature gene selection based on machine learning

Three complementary machine learning algorithms were used for feature selection from the 123 identified DEGs. LASSO regression analysis identified 26 genes with non-zero regression coefficients (Figure 4A,4B). The SVM-RFE algorithm screened out a 21-gene subset with the lowest cross-validation error (Figure 4C,4D), while the RF algorithm retained the 32 top-ranked genes based on MeanDecreaseGini (Figure 4E,4F; the top 30 are shown in Figure 4F).

Figure 4 Feature selection using machine learning techniques. (A) A 10-fold cross-validation approach to identify the optimal λ value in LASSO regression, where the model with the minimum binomial deviance was selected as the best performer. (B) LASSO coefficient trajectory curves for 38 candidate genes, demonstrating that most coefficients shrink to zero as the L1 norm decreases. (C) Accuracy curve from 10-fold cross-validation of the SVM-RFE algorithm, achieving the highest accuracy (0.926) when 21 features were included. (D) Error curve from 10-fold cross-validation of the SVM-RFE algorithm, with the lowest error (0.0735) obtained using 21 features, confirming it as the optimal feature set. (E) Out-of-bag error rate as a function of the number of trees in the random forest model. (F) The top 30 genes identified by the RF algorithm, ranked by their variable importance score. CV, cross-validation; LASSO, least absolute shrinkage and selection operator; RF, random forest; SVM-RFE, support vector machine recursive feature elimination.

Model construction, performance evaluation and interpretability analysis

The overlapping genes among the results of LASSO, RF, and SVM-RFE algorithms (RAB17, FAM107A, COL9A1, HOXC4, and TRPM4) were determined as the final feature set for subsequent model construction (Figure 5A). Based on the expression profiles of these five genes, we constructed 10 machine learning models using the training cohort and verified their predictive efficacy in the independent testing cohort. Among the 10 evaluated algorithms, the PLS model showed the best diagnostic performance in the testing cohort, with an AUC value of 0.977 (Figure 5B).

Figure 5 Predictive performance of the model and SHAP-based model interpretability. (A) Venn diagram illustrating the overlap of genes identified by the LASSO, RF, and SVM-RFE algorithms, with five final signature genes (RAB17, FAM107A, COL9A1, HOXC4, TRPM4) highlighted. (B) ROC curves for 10 classifiers in the testing cohort, demonstrating that the PLS model achieved the highest AUC value of 0.977. (C) SHAP summary bar plot for the PLS model, with genes ranked by their mean absolute SHAP values. (D) SHAP beeswarm plot displaying the direction and magnitude of each gene’s effect across individual patients. (E) SHAP waterfall plot for a representative patient, outlining the gene-specific contributions to the prediction threshold of the PLS model. AUC, area under the curve; DTS, decision trees; GBM, gradient boosting machine; glmBoost, generalized boosted regression model; KNN, K-nearest neighbors; LASSO, least absolute shrinkage and selection operator; NeuralNet, neural network; PLS, partial least squares; RF, random forest; ROC, receiver operating characteristic; SHAP, SHapley Additive exPlanations; SVM, support vector machine; SVM-RFE, support vector machine recursive feature elimination; XGBoost, extreme gradient boosting.

Subsequently, we performed SHAP analysis to clarify the interpretability of the optimal PLS model. The SHAP summary bar plot showed that RAB17 made the most significant contribution to the predictive performance of the model, while TRPM4 had the smallest impact on the model output (Figure 5C). The SHAP beeswarm plot (Figure 5D) further illustrated the strength and direction of the effect of each gene on the prediction at the individual patient level. High expression of RAB17, HOXC4, and TRPM4 was associated with positive SHAP values, indicating that the upregulation of these genes would shift the model prediction to the direction of positive PCa diagnosis, acting as positive predictive markers. In contrast, high expression of FAM107A and COL9A1 was negatively correlated with SHAP values, suggesting that these genes play a protective role in PCa and can be used as negative predictive markers.

The SHAP waterfall plot further clarified the contribution of each gene to the prediction results in representative patient samples (Figure 5E). The baseline expected prediction value of the PLS model was E[f(x)] = 0.58. Taking a representative sample as an example, RAB17 (expression = 5.01) reduced the predicted probability of PCa by 0.182; TRPM4 (expression = 5.47) further reduced it by 0.0756; FAM107A (expression = 6.97) caused an additional reduction of 0.0427; COL9A1 (expression = 3.77) increased the predicted probability by 0.0737; HOXC4 (expression = 5.66) increased it by 0.0515. The combined effect of these genes resulted in a final prediction value of f(x) = 0.405 for this sample.

Expression profiles, prognostic evaluation and clinical subgroup analysis of signature genes

We first evaluated the expression and prognostic significance of the identified signature genes in the TCGA-PRAD cohort. Violin plots showed that RAB17, HOXC4, and TRPM4 were significantly upregulated, whereas FAM107A and COL9A1 were downregulated in PCa tissues compared with normal controls (P<0.0001; Figure 6A). Survival analysis further showed that FAM107A was significantly associated with both OS (P=0.04) and PFS (P<0.001), whereas high HOXC4 expression was associated with poor PFS (P=0.02; Figure 6B-6D). Considering its prognostic relevance and the limited evidence regarding HOXC4 in PCa, HOXC4 was selected for subsequent comprehensive analysis.

Figure 6 Expression profiles and prognostic value of the identified feature genes. (A) Violin plot comparing gene expression levels between normal prostate and prostate tumor tissues from the TCGA cohort. ****, P<0.0001. (B) Kaplan-Meier overall survival curves for FAM107A (log-rank P=0.04). (C,D) The progression-free survival curves of FAM107A (log-rank P<0.001) and HOXC4 (log-rank P=0.02). (E) Kaplan-Meier DFS curves stratified by HOXC4 expression levels in the independent external validation cohort from MSKCC. (F) Forest plot summarizing the subgroup survival analysis of HOXC4 expression across different clinical strata based on age and clinical risk classification. CI, confidence interval; DFS, disease-free survival; MSKCC, Memorial Sloan Kettering Cancer Center; TCGA, The Cancer Genome Atlas.

External validation in the independent MSKCC cohort demonstrated that high HOXC4 expression was associated with shorter disease-free survival (DFS) (log-rank P=0.041; Figure 6E). Stratified survival analysis further suggested that elevated HOXC4 expression remained associated with unfavorable prognosis across selected clinical strata, particularly in patients aged ≤65 years (HR: 2.57, P<0.001) and in the low-risk subgroup (HR: 2.69, P=0.001; Figure 6F). These findings support HOXC4 as a potential prognostic biomarker with relevance across different clinical contexts.

Functional enrichment analysis of HOXC4

GSEA results showed that the HOXC4 high-expression group (Figure 7A) was significantly enriched in oxidative phosphorylation, proteasome, protein export, and ubiquitin-mediated protein degradation pathways. These pathways are mainly involved in cellular energy metabolism, protein processing and degradation, suggesting that the upregulation of HOXC4 expression may drive metabolic reprogramming and protein turnover during PCa progression. In contrast, the HOXC4 low-expression group (Figure 7B) showed an enrichment trend in cytokine-cytokine receptor interaction, drug metabolism-cytochrome P450, and xenobiotic metabolism pathways, which are involved in immune regulation and metabolic response BP.

Figure 7 GSEA for HOXC4. (A) Top 5 KEGG pathways significantly enriched in the HOXC4 high expression group. (B) Top 5 pathways enriched in the HOXC4 low expression group. GSEA, gene set enrichment analysis; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Immune cell infiltration profile in normal and PCa groups and its association with HOXC4

The CIBERSORT algorithm was used to evaluate the immune infiltration landscape of the samples, which were divided into the normal group and the PCa group (Figure 8A). Comparative analysis showed that there were significant differences in the immune infiltration profiles between the two groups: compared with the normal group, the proportion of Tregs and M2 macrophages in the PCa group was significantly increased, while the proportion of follicular helper T (Tfh) cells, monocytes, and resting mast cells was significantly decreased (Figure 8B). These results indicate that the PCa TME presents a typical immunosuppressive phenotype.

Figure 8 Immune infiltration landscape and its association with HOXC4 expression in normal and PCa groups. (A) Stacked bar plot showing the relative fractions of immune cell subsets inferred by CIBERSORT for normal and PCa samples. (B) Box plots showing the differential immune cell infiltration between normal and PCa groups. (C) Spearman correlation heatmap among the 22 immune cell subsets. (D,E) Scatter plots illustrating the direct Spearman correlation between HOXC4 expression levels and the infiltration fractions of (D) Tregs and (E) M2 macrophages in the TCGA-PRAD cohort. (F) Box plot comparing the infiltration fractions of 22 immune cell types between HOXC4-high and HOXC4-low expression groups. (G) Lollipop plot displaying the Spearman correlation coefficients between HOXC4 expression and various classical immune checkpoint molecules. *, P<0.05; **, P<0.01; ***, P<0.001; ns, not significant. PCa, prostate cancer; TCGA-PRAD, The Cancer Genome Atlas-Prostate Adenocarcinoma; Tregs, regulatory T cells; TPM, transcripts per million.

Spearman correlation heatmap showed a strong negative correlation (r=−0.58) between M1 and M2 macrophages, a strong positive correlation (r=0.69) between CD8⁺ T cells and M1 macrophages, and a negative correlation (r=−0.36) between regulatory T cells and M1 macrophages (Figure 8C).

To further evaluate the relationship between HOXC4 and the immunosuppressive microenvironment, we analyzed the correlations between HOXC4 expression and the infiltration fractions of regulatory T cells (Tregs) and M2 macrophages. Spearman correlation analysis showed that HOXC4 expression was positively correlated with Treg infiltration (r=0.16, P<0.001; Figure 8D) and M2 macrophage infiltration (r=0.17, P<0.001; Figure 8E) in the TCGA-PRAD cohort. These findings suggest that HOXC4 expression may be associated with an immune-suppressive microenvironment in PCa. Furthermore, patients from the TCGA-PRAD cohort were stratified into HOXC4-high and HOXC4-low expression groups to characterize the divergent immune landscapes. The CIBERSORT boxplot revealed distinct infiltration patterns of the 22 immune cell types between the two strata, highlighting significant variations in multiple immune cell subsets (Figure 8F). Given the prominent association of HOXC4 with immune suppression, we extended our analysis to explore its relationship with key immune checkpoint molecules. Spearman correlation analysis demonstrated that HOXC4 expression was associated with several immune checkpoints. Significant negative correlations were observed with CD40 and SIGLEC-15 (both P<0.001), PD-L1 (P<0.01), PD-1 (P<0.05), and TIGIT (P<0.05), whereas B7-H3 showed a significant positive correlation (P<0.05). No significant associations were observed with CTLA-4, TIM-3, LAG-3, ICOS, or CD86 (Figure 8G). These results underscore the potential clinical utility of HOXC4 as a candidate bio-indicator for predicting response or resistance to immune checkpoint blockade (ICB) therapies.

Analysis of HOXC4 expression in pan-cancers

TIMER analysis revealed that HOXC4 exhibits highly heterogeneous expression across diverse cancers. It functions as an oncogene (upregulated in tumors) in cancers, such as BLCA (bladder urothelial carcinoma), BRCA (breast cancer), CHOL (cholangiocarcinoma), HNSC (head and neck squamous cell carcinoma), LIHC (hepatocellular carcinoma), LUAD (lung adenocarcinoma), LUSC (lung squamous cell carcinoma), PCPG (pheochromocytoma and paraganglioma), and PRAD (prostate adenocarcinoma). Conversely, it functions as a tumor suppressor gene (downregulated in tumors) in other cancers, including KIRC (kidney renal clear cell carcinoma), READ (rectal adenocarcinoma), THCA (thyroid carcinoma), and UCEC (endometrial adenocarcinoma). Notably, no significant expression differences were observed in some cancers. These variations provide a preliminary landscape of HOXC4 function in pan-cancer initiation and progression (Figure 9).

Figure 9 Pan-cancer analysis of HOXC4. The expression level of HOXC4 in different cancers. **, P<0.01; ***, P<0.001; ns, not significant. TPM, transcripts per million.

Discussion

In this study, we systematically explored the core value of HOXC4 as an immune-related biomarker for the treatment and prognostic evaluation of PCa by integrating a variety of interpretable machine learning algorithms. We first integrated multiple PCa transcriptomic datasets, eliminated the batch effect for standardized processing, and screened out DEGs between tumor and normal tissues. Then, we used three complementary machine learning algorithms, LASSO, RF, and SVM-RFE, to screen and identify five core signature genes including HOXC4. Among the diagnostic models constructed based on these signature genes, the PLS model showed excellent diagnostic performance, and SHAP analysis further confirmed that HOXC4 was a key risk factor for PCa. We found that HOXC4 was significantly highly expressed in PCa tissues and correlated with shorter PFS of patients, suggesting that it has important potential in the prognostic evaluation of PCa. In addition, GSEA analysis showed that the high HOXC4 expression group was significantly enriched in energy metabolism and protein processing-related pathways, while the low expression group was mainly enriched in immune regulation-related pathways. CIBERSORT analysis confirmed that PCa has a typical immunosuppressive microenvironment, and the expression of HOXC4 is significantly correlated with the infiltration level of key immune cells, suggesting that HOXC4 may participate in the progression of PCa by regulating the tumor immune microenvironment. Pan-cancer analysis also confirmed the heterogeneous expression pattern of HOXC4 in different tumor types, and it showed an oncogene-like high expression characteristic in PCa.

The discovery of DEGs between PCa and normal tissues reveals a complex transcriptional regulatory network involving upregulated and downregulated genes in carcinogenesis. Mechanistically, upregulated genes such as HOXC6 and AMACR exert carcinogenic effects by regulating cell proliferation and metabolic adaptation (9). Importantly, abnormal expression of homeobox genes such as HOXC4 can alter the prostate transcriptome by competing with transcription regulatory factors such as HOXB13, FOXA1, and androgen receptors, affecting key pathways in tumor growth (9). Downregulation of genes may indicate inhibition of tumor suppressor pathways, which is consistent with the dual role of gene expression regulation in previous transcriptome studies (15,16). These DEGs promote malignant phenotypes by enriching BP and signaling networks, highlighting the necessity of adopting integrated methods to analyze the interaction between oncogenes and tumor suppressor genes.

Functional enrichment analysis further clarified the role of DEGs in tumor biology. GO and KEGG analysis showed that DEGs are mainly involved in cytoskeleton assembly, extracellular matrix composition, and metabolic regulation, which are key links in cancer cell invasion and metastasis (16,17). The enrichment of DEGs in the intermediate filament assembly and basement membrane interaction pathways provides a mechanistic basis for the high migration and invasion ability of PCa cells. In addition, the identification of MF such as fatty acid binding and peroxidase activity is consistent with the research conclusion that metabolic reprogramming and redox homeostasis are characteristics of malignant tumors (18,19). Other cancer studies have also confirmed that the enrichment of DEGs can reveal metabolic and immune escape vulnerabilities, validating the association of this discovery (20,21).

LASSO, RF and machine learning techniques such as SVM-RFE helped screen out five feature genes containing HOXC4. Although the mathematical principles of each algorithm are different, they all identify overlapping gene sets with biological significance (20,22). The penalty mechanism of LASSO is applicable to high-dimensional transcriptome data, while RF and SVM-RFE are adept at handling complex nonlinear relationships between variables (21,23). Multi-algorithm consistent recognition of HOXC4 highlights its core position. This integration method overcomes the problem of feature omission caused by bias or overfitting in a single algorithm (21), and improves the accuracy of biomarker screening.

Survival analysis shows that there are differences in the prognostic value of characteristic genes: high expression of HOXC4 is associated with shortened PFS, while high expression of FAM107A is associated with prolonged OS. This duality is consistent with the conclusion that dysregulation of homeobox genes in other cancers suggests poor prognosis, as it can promote cell cycle progression, resist apoptosis, and promote metastasis (7,24). Specifically, HOXC4 has been proven to enhance glycolysis and cell proliferation of pancreatic cancer and liver cancer, confirming its carcinogenic metabolic reprogramming effect (6,24). On the contrary, as a tumor suppressor gene, FAM107A—whose deletion is associated with increased tumor invasiveness (14). The complexity of this genetic prognostic association highlights the importance of interpreting it in conjunction with cellular and molecular backgrounds. Importantly, the prognostic robustness of HOXC4 was rigorously confirmed in an independent MSKCC cohort, mitigating the risk of cohort-specific bias. Furthermore, our stratified survival analysis revealed that HOXC4 acts as a highly prominent risk factor particularly in younger (≤65 years) and low-risk patients. This implies that HOXC4 possesses unique clinical utility in identifying patients with early-stage or seemingly low-risk profiles who are covertly predisposed to aggressive tumor progression, thereby aiding in more precise and personalized clinical surveillance.

GSEA further provided mechanistic clues for the biological role of HOXC4 in PCa. HOXC4-high samples were significantly enriched in oxidative phosphorylation, proteasome, protein export, and ubiquitin-mediated proteolysis pathways, indicating that HOXC4 may be involved in the regulation of energy metabolism and protein homeostasis. Oxidative phosphorylation is a major source of mitochondrial energy production and may support the elevated bioenergetic requirements of tumor cells. Meanwhile, proteasome- and ubiquitin-mediated proteolysis pathways are essential for protein turnover and cellular stress adaptation, which may facilitate malignant cell survival under unfavorable microenvironmental conditions. These findings are consistent with previous studies showing that HOXC4 promotes tumor progression and metabolic reprogramming in pancreatic cancer and hepatocellular carcinoma (6,24), and with pan-cancer evidence suggesting that HOX family dysregulation is broadly associated with tumor progression and immune microenvironment remodeling across solid tumors (25). Importantly, the enrichment of these metabolic and protein-processing pathways was observed together with increased infiltration of Tregs and M2 macrophages in HOXC4-high PCa, suggesting that HOXC4 may participate in metabolic-immune crosstalk within the TME. Given that nutrient competition and metabolic stress are important contributors to T cell dysfunction and immune escape (26,27), HOXC4 may link metabolic remodeling with immune suppression during PCa progression. However, this hypothesis remains to be validated by further experimental studies.

CIBERSORT immune cell infiltration deconvolution analysis confirmed that PCa has a typical immunosuppressive microenvironment and is significantly correlated with high expression of HOXC4. The proportion of Tregs and M2 macrophage infiltration increased in the HOXC4 high expression group, while Tfh, monocytes, and resting mast cells infiltration decreased, providing a basis for immune escape. As core immunosuppressive cells, Tregs inhibit T cell activation and proliferation by secreting IL-10, TGF-β, and expressing CTLA-4. The association between HOXC4 and Tregs infiltration suggests that it may upregulate the Tregs recruitment/proliferation pathway. This is consistent with the research pattern of “upregulation of HOXC4 promotes immune escape through the PD-L1 axis” in colorectal cancer (28), suggesting that HOXC4 in PCa may synergistically regulate Tregs infiltration and PD-L1 expression to form a dual inhibitory barrier. M2 macrophages promote tumor growth and inhibit T cells by secreting VEGF and expressing Arg-1. High expression of HOXC4 is associated with increased M2 infiltration, negative correlation between M1 and M2, and positive correlation between CD8⁺ T cells and M1, suggesting that it may promote macrophage polarization towards M2 and indirectly weaken CD8⁺ T cell function (29). Beyond altering the cellular composition of the TME, HOXC4 expression was associated with multiple immune checkpoint molecules. HOXC4 showed negative correlations with CD40, SIGLEC-15, PD-L1, PD-1, and TIGIT, but a positive correlation with B7-H3; no significant correlation was observed with CTLA-4. SIGLEC-15 has emerged as an important immune suppressor comparable to PD-L1 and is often expressed in tumors lacking PD-L1. Together, these findings suggest that HOXC4 is associated with a distinct immune-regulatory state in PCa rather than a generalized upregulation of immune checkpoints. Further mechanistic and clinical studies are needed to determine whether these correlations have predictive relevance for ICB. In addition, similar immune subtype classifications also exist in urothelial carcinoma and lung cancer (29,30), suggesting that HOXC4 may be a cross-cancer immune microenvironment molecular regulatory factor, providing mechanistic support for its position as an immune biomarker and immunotherapy sensitization target.

Pan-cancer investigations showed that the expression and function of HOXC4 were tissue-specific: it was overexpressed as an oncogene in bladder cancer, breast cancer, cholangiocarcinoma, liver cancer and lung cancer, and was associated with poor prognosis and invasion (7); low expression in renal clear cell carcinoma and rectal adenocarcinoma suggests tissue-specific regulation or compensatory pathways (7). Similar differences also exist among other members of the HOX family, whose oncogene/tumor suppressor functions are influenced by the microenvironment and co-regulatory factors (25). In terms of mechanism, this heterogeneity is related to lineage-specific transcription factors, chromatin modifying factors and non-coding RNAs. For example, lncRNA/miRNA in breast cancer and uveal melanoma regulates HOXC4 (8,31), suggesting that it is necessary to integrate multi-omics and tissue-specific analysis to analyze the cancer biological function of HOXC4.

Several limitations should be acknowledged. First, this study was primarily based on retrospective analyses of publicly available transcriptomic datasets. Although internal model evaluation and external validation using the MSKCC cohort were performed, the expression and clinical relevance of HOXC4 still require confirmation in independent, large-scale clinical PCa tissue cohorts. Second, the current framework provides predictive and correlative evidence rather than direct mechanistic validation. The biological mechanisms by which HOXC4 may regulate metabolic reprogramming and immune microenvironment remodeling remain to be clarified. Future in vitro and in vivo experiments, together with prospective multicenter clinical validation, are needed to determine the functional role and translational value of HOXC4 in PCa.

This study confirms that HOXC4 is a potential immune biomarker for PCa treatment response and prognosis, associated with mediating immune suppression and metabolic reprogramming. The interaction between the intrinsic characteristics of tumors and the host immune microenvironment is crucial for the progression and treatment sensitivity of PCa. Data show that tumors with high expression of HOXC4 tend to form an immunosuppressive microenvironment and enrich energy metabolism and protein processing pathways, while tumors with low expression exhibit immune regulatory responses. Targeted regulation of HOXC4 holds the potential to improve the prognosis of PCa patients, and thus deserves further clinical translational evaluation. Incorporation of HOXC4 expression status into machine learning-based predictive models is expected to guide individualized treatment decisions for PCa patients. The findings of this study highlight the application value of combining machine learning-driven feature screening with systematic biological validation in the precise diagnosis and treatment of PCa.

Conclusions

Identified through interpretable machine learning and externally validated clinical cohorts, HOXC4 may serve as a promising diagnostic and prognostic biomarker for PCa. High HOXC4 expression was associated with unfavorable clinical outcomes, metabolic reprogramming, and an immunosuppressive TME characterized by increased regulatory T cell and M2 macrophage infiltration. These findings suggest that HOXC4 may be involved in metabolic-immune crosstalk during PCa progression and may provide a potential basis for improved risk stratification, prognostic assessment, and individualized therapeutic decision-making. Further experimental and prospective clinical studies are warranted to validate the mechanistic role and translational applicability of HOXC4 in PCa.

Acknowledgments

The authors would like to thank the TCGA and GEO databases for providing the public data, as well as the developers of all R packages used in this study.

Footnote

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

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

Funding: This study was supported by the Students’ Platform for Innovation and Entrepreneurship Training Program of Jining Medical University (Grant No. S202510443022), the Breeding Program for High-level Scientific Research Projects at Jining Medical University (Grant No. JYGC2023KJ017), the Medical and Health Technology Project of Shandong Province (Grant No. 202501061008), and the Research Fund for Lin He’s Academician Workstation of New Medicine and Clinical Translation in Jining Medical University (Grant No. YHL.2025ZD03).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tau.amegroups.com/article/view/10.21037/tau-2026-0231/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. Prostate cancer. Nat Rev Dis Primers 2021;7:8.
2. Moon SK, You MW, Lim JW. Granulomatous Prostatitis Mimicking Prostate Cancer. Urology. 2020;137:e3-e5. [Crossref]
3. Krušlin B, Škara L, Vodopić T, et al. Genetics of Prostate Carcinoma. Acta Med Acad 2021;50:71-87. [Crossref] [PubMed]
4. Rebello RJ, Oing C, Knudsen KE, et al. Prostate cancer. Nat Rev Dis Primers 2021;7:9. [Crossref] [PubMed]
5. Peng Y, Song Y, Ding J, et al. Identification of immune-related biomarkers in adrenocortical carcinoma: Immune-related biomarkers for ACC. Int Immunopharmacol 2020;88:106930. [Crossref] [PubMed]
6. Yang T, Zhang XB, Li XN, et al. Homeobox C4 promotes hepatocellular carcinoma progression by the transactivation of Snail. Neoplasma 2021;68:23-30. [Crossref] [PubMed]
7. Xiao J, Li Y, Liu Y, et al. The involvement of homeobox-C 4 in predicting prognosis and unraveling immune landscape across multiple cancers via integrated analysis. Front Genet 2022;13:1021473. [Crossref] [PubMed]
8. Zhao S, Song C, Chen F, et al. LncRNA XIST/miR-455-3p/HOXC4 axis promotes breast cancer development by activating TGF-β/SMAD signaling pathway. Funct Integr Genomics 2024;24:159. [Crossref] [PubMed]
9. Luo Z, Farnham PJ. Genome-wide analysis of HOXC4 and HOXC6 regulated genes and binding sites in prostate cancer cells. PLoS One 2020;15:e0228590. [Crossref] [PubMed]
10. Kumar P, Paul RK, Roy HS, et al. Big Data Analysis in Computational Biology and Bioinformatics. Methods Mol Biol 2024;2719:181-97. [Crossref] [PubMed]
11. Nalina V, Prabhu D, Sahayarayan JJ, et al. Advancements in AI for Computational Biology and Bioinformatics: A Comprehensive Review. Methods Mol Biol 2025;2952:87-105. [Crossref] [PubMed]
12. Hajjo R, Sabbah DA, Bardaweel SK, et al. Identification of Tumor-Specific MRI Biomarkers Using Machine Learning (ML). Diagnostics (Basel) 2021;11:742. [Crossref] [PubMed]
13. Taylor BS, Schultz N, Hieronymus H, et al. Integrative genomic profiling of human prostate cancer. Cancer Cell 2010;18:11-22. [Crossref] [PubMed]
14. Qiu Z, Du X, Chen K, et al. Gene signatures with predictive and prognostic survival values in human osteosarcoma. PeerJ 2021;9:e10633. [Crossref] [PubMed]
15. Wang S, Zhang Y, Hu C, et al. Shiny-DEG: A Web Application to Analyze and Visualize Differentially Expressed Genes in RNA-seq. Interdiscip Sci 2020;12:349-54. [Crossref] [PubMed]
16. Yin H, Duo H, Li S, et al. Unlocking biological insights from differentially expressed genes: Concepts, methods, and future perspectives. J Adv Res 2025;76:135-57. [Crossref] [PubMed]
17. Yang P, Li Y, Li T, et al. Screening differentially expressed genes and the pathogenesis in atopic dermatitis using bioinformatics. Cell Mol Biol (Noisy-le-grand) 2023;69:73-8. [Crossref] [PubMed]
18. Nair AR, Kaniyala H, Vardhan MH, et al. Differentially Expressed Genes (DEGs) in Umbelliferone (UMB) Producing Endophytic Fusarium oxysporum (ZzEF8) Following Epigenetic Modification. J Basic Microbiol 2025;65:e2400582. [Crossref] [PubMed]
19. Zhou H, Jiang J, Chen X, et al. Differentially expressed genes and miRNAs in female osteoporosis patients. Medicine (Baltimore) 2022;101:e29856. [Crossref] [PubMed]
20. Zhou Q, Lan L, Wang W, et al. Identifying effective immune biomarkers in alopecia areata diagnosis based on machine learning methods. BMC Med Inform Decis Mak 2025;25:23. [Crossref] [PubMed]
21. Ma K, Nakajima H, Basak N, et al. Integrating explainable machine learning and transcriptomics data reveals cell-type specific immune signatures underlying macular degeneration. NPJ Genom Med 2025;10:48. [Crossref] [PubMed]
22. Fang D, Lin J, Wang J, et al. CEACAM6 as a machine learning derived immune biomarker for predicting neoadjuvant chemotherapy response in HR+/HER2- breast cancer. Front Immunol 2025;16:1662004. [Crossref] [PubMed]
23. Xu Z, Hao Q, Yan B, et al. Immuno-transcriptomic analysis based on machine learning identifies immunity signature genes of chronic rhinosinusitis with nasal polyps. Sci Rep 2025;15:19393. [Crossref] [PubMed]
24. Zhang H, Han B, Tian S, et al. HOXC4 promotes proliferation of pancreatic cancer cells by increasing LDHA-mediated glycolysis. Aging (Albany NY) 2024;16:11103-16. [Crossref] [PubMed]
25. Wang Y, Gao J, Ren Z, et al. A pan-cancer analysis of homeobox family: expression characteristics and latent significance in prognosis and immune microenvironment. Front Oncol 2025;15:1521652. [Crossref] [PubMed]
26. Zhao S, Peralta RM, Avina-Ochoa N, et al. Metabolic regulation of T cells in the tumor microenvironment by nutrient availability and diet. Semin Immunol 2021;52:101485. [Crossref] [PubMed]
27. Sharma P, Guo A, Poudel S, et al. Early methionine availability attenuates T cell exhaustion. Nat Immunol 2025;26:1384-96. [Crossref] [PubMed]
28. Liu L, Yu T, Jin Y, et al. MicroRNA-15a Carried by Mesenchymal Stem Cell-Derived Extracellular Vesicles Inhibits the Immune Evasion of Colorectal Cancer Cells by Regulating the KDM4B/HOXC4/PD-L1 Axis. Front Cell Dev Biol 2021;9:629893. [Crossref] [PubMed]
29. Peng M. Immune landscape of distinct subtypes in urothelial carcinoma based on immune gene profile. Front Immunol. 2022;13:970885. [Crossref]
30. Song Y, Yan S, Fan W, et al. Identification and Validation of the Immune Subtypes of Lung Adenocarcinoma: Implications for Immunotherapy. Front Cell Dev Biol 2020;8:550. [Crossref] [PubMed]
31. Wu S, Chen H, Zuo L, et al. Suppression of long noncoding RNA MALAT1 inhibits the development of uveal melanoma via microRNA-608-mediated inhibition of HOXC4. Am J Physiol Cell Physiol 2020;318:C903-12. [Crossref] [PubMed]