Expression profile information and relevant clinical information for a total of 622 patients with colon cancer (COAD) and rectal cancer (READ) were downloaded from TCGA (https://tcga-data.nci.nih.gov/tcga/ ) database. Thirty-two data without complete survival information were removed. Mutation data and expression profiles for pan-cancer were also downloaded via the TCGA database. The GSE81582 dataset containing expression profiles and relevant clinical information for 19 CRLM and 23 primary CRC samples was downloaded from GEO (https://www.ncbi.nlm.nih.gov/geo/ ) database.

Differential expression analysis was performed using the ‘limma’ package to select differentially expressed genes (DEGs) (24). The absolute value of logarithmically converted fold change (|Log2FC|) > 0.6 and utilized a false discovery rate (FDR) filtering level of 0.05.

Univariate Cox regression analysis identified genes associated with prognosis. least absolute shrinkage and selection operator (LASSO) analysis was performed using the ‘glmnet’ package (25) in R language to screen the genes further. Moreover, the following formula was used to determine the risk score: r i s k   s c o r e = ∑ i = 1 7 C o e f i * X i . The C o e f i is the risk factor calculated by the LASSO model for each gene, and the X i is the expression of each gene. The TCGA sample data were randomly divided into training and validation groups at a 1:1 ratio. Using the median risk score value, patients inside this training cohort were split into low- and high-risk groups. The survivor and ‘survminer’ package were used to generate overall survival (OS) Kaplan-Meier curves for high- and low-risk groups. To compare the survival curves, the log-rank test was used. The ‘timeROC’ package (26) developed receiver operating characteristic (ROC) curves to evaluate the predicting efficacy of the signature. Decision Curve Analysis (DCA) was performed using the ‘dca.R’ package.

To create a nomogram for predicting patient outcomes, we utilized the ‘rms’ package. We fitted a multivariate Cox proportional hazards regression model using the significant prognostic factors identified in our study. Subsequently, we employed the ‘cph’ function to estimate the regression coefficients and hazard ratios for each variable. We then used the ‘nomogram’ function to construct a graphical representation of the Cox model, facilitating an intuitive visualization and interpretation of the results. The performance of the nomogram was assessed using calibration plots, which were generated with the ‘calibrate’ function in the ‘rms’ package.

The proportion of 22 immune cell infiltrates in CRC patients was calculated using the ‘CIBERSORT’ package. Stromal score, immune score, and tumor purity of malignant tumors were assessed using the ‘ESTIMATE’ package. The expression of several immune checkpoint levels was compared using the Wilcoxon test (27). The immunophenoscore (IPS) was obtained through The Cancer Immunome Atlas (TCIA) (https://tcia.at/home ) database. The tumor immune dysfunction and exclusion (TIDE) score was calculated using the online website TIDE (http://tide.dfci.harvard.edu/ ). Neoantigen data for tumors were obtained from a previous study (28).

The response of each patient to each drug treatment was predicted by estimating the half-maximal inhibitory concentration (IC50) using the ‘oncoPredict’ package, and associations with 198 drugs were calculated using Spearman rank correlation coefficients.

The expression of HAMP in several cell lines was verified using the cancer cell line encyclopedia (CCLE, broadinstitute.org) database. HAMP expression in pan-cancer was observed using TIMER2.0 (cistrome.org). The MSI score for each tumor was obtained from a previous study (29). The TMB of each tumor was calculated using the ‘tmb’ function of the ‘maftools’ package. The MSI, TMB, and gene expression data of the samples were integrated separately for correlation analysis.

GO and KEGG enrichment analyses were performed on the identified DEGs using the ‘clusterProfiler package (30). The adjusted p-value< 0.05 was used as the criterion to filter significantly enriched GO terms and KEGG pathways. The c2.cp.kegg.v7.0.symbols.gmt genome was obtained through the MSigDB database (http://software.broadinstitute.org/gsea/msigdb ). GSEA analysis with the above data. The p-value< 0.05 was used as filtering criteria to screen significantly activated KEGG pathways.

Eight pairs of fresh CRC and adjacent normal tissues were obtained from CRC patients treated with radical surgery at the First Affiliated Hospital of Xi’an Jiaotong University. All subjects signed an informed consent form. In compliance with the Helsinki Declaration, the study was approved by the Medical Ethics Committee of the First Affiliated Hospital of Xi’an Jiaotong University. As previously reported, the immunohistochemical staining procedure was performed using the streptavidin peroxidase-conjugated method (SP-IHC) (31). The hepcidin antibody (cat. no. sc-100277) was purchased from (Santa Cruz Biotechnology, Texas, USA).

RPMI 1640 media (Gibco BRL, Carlsbad, CA, USA) supplemented with 10% FBS (Gibco BRL, Carlsbad, CA, USA) was used to maintain human colon cancer cells (DLD-1) and intestinal epithelial cells (FHCs) (Shanghai Institute of Cell Biology, Chinese Academy of Sciences, China). Cells were cultured at 37°C with 5% CO₂.

The target gene (HAMP) silenced RNA (si-RNA) was constructed by the Life Technologies (Thermo Fisher Scientific), and we transfected the DLD-1 cells with 100 nM siRNA using lipofectamine 2000 (Invitrogen, CA, USA) according to the manufacturer’s protocol (Invitrogen, CA, USA).

The TRIzol reagent was used to isolate the total RNA of cells (Invitrogen, Carlsbad, CA, USA). Then, the obtained RNAs were used for cDNA synthesis with the PrimeScript RT Reagent kit (TaKaRa, Osaka, Japan). cDNA was next used for qPCR with SYBR Green fluorescence signal detection assays (TaKaRa, Osaka, Japan) on an IQ5 instrument (Bio-Rad, CA, USA). The forward primer for HAMP is 5’- CTGACCAGTGGCTCTGTTTTC -3’, and the reverse primer is 5’- GAAGTGGGTGTCTCGCCTC-3’. Using the 2^−ΔΔCT approach, the level of mRNA expression was measured.

The RPMI 1640 medium supplemented with 10% FBS medium was placed in the bottom compartments of a transwell chamber, and RPMI 1640 serum-free medium containing 1*10⁵ transfected DLD-1 cells was added to the upper compartments. After incubation for 24h at 37°C and 5% CO₂, the upper chambers were emptied and fixed using methanol. The migrating cells were then stained with 0.1% crystal violet. Then the staining was observed under an inverted phase contrast microscope. The migrated cells were counted in five randomly selected high-power visual fields, and their number was averaged. The experiments were repeated to verify the consistency and reliability of the experimental results. Each experiment was performed three times under the same conditions and using the same methods and materials.

The chi-square test was used to compare the composition ratios of categorical variables between the two groups. The log-rank test completed the comparison of each Kaplan-Meier curve. The correlation was evaluated using the Pearson correlation analysis and Spearman correlation analysis. The Wilcoxon signed-rank test completed the comparison of different groups. The statistical analysis was based on the R (4.0.2) programming language. P value< 0.05 was the filtering criterion.

There were 1835 up-regulated genes and 9516 down-regulated genes in CRC samples compared to normal tissue in the TCGA dataset ( Figure 1A ), and 102 up-regulated genes and 33 down-regulated genes in CRLM samples compared to primary CRC samples in the GSE81582 dataset ( Figure 1B ). Details of these two differential gene sets were recorded in ( Supplementary Tables 1 , 2 ). Ninety-six genes were identified after extracting the intersection of the two differential gene sets ( Figure 1C ).

The TCGA sample data were randomly divided into training and validation groups at a 1:1 ratio ( Table 1 ). The training and validation groups were used to study the prognostic importance of CRLM-related genes by constructing and validating CRLM-related prognostic features together with the overall group. Using univariate Cox analysis, we discovered CRLM genes that are related to prognosis. Then, using LASSO analysis, 7 important CRLM-related genes (HAMP, COLEC11, MMP3, UGT2B7, C8G, SERPINA1, IFITM10) were further identified, and signatures were constructed. ( Figure 1D, E ). The following formula was used to compute the risk score: r i s k   s c o r e = ∑ i = 1 7 C o e f i * X i . The C o e f i for each gene and the expression of each gene in all samples can be found in Supplementary Table 3 . Patients were categorized based on their median risk score. This sample distribution was reasonable in both three groups according to the distribution of risk scores and the distribution of survival status. ( Figures 2A–C ). According to Kaplan-Meier analysis, the high-risk group’s OS was significantly lower in the training, validation, and overall groups compared to the low-risk group ( Figures 2D–F ). We displayed time-dependent ROC curves to measure the prediction model’s effectiveness and calculated the area under the curve (AUC). The Decision Curve Analysis (DCA) plots effectively demonstrated the clinical utility and net benefit of our model across a range of threshold probabilities ( Supplementary Figure 1 ). The results showed that this prediction model was valuable in the training, validation, and overall groups in predicting OS in the short and long term. ( Figures 2G–I ).

Risk Score, Age, and AJCC Stage were independent indicators of CRC prognosis by multivariate Cox regression ( Figure 3A ). Based on the findings, a nomogram ( Figure 3B ) was developed to predict OS by combining risk score with age and AJCC stage. Finally, calibration curves were used to assess the nomogram’s prediction capabilities. According to the statistics, actual survival times were quite close to the predictions of this nomogram for 1, 3, and 5 years ( Figures 3C–E ). The information above indicates that this nomogram has a strong predictive ability.

First, comparing the percentage of immune infiltrating cells in the different risk groups, statistically, significant differences were found for Plasma cells, Some T cells ( Figure 4A ). Although there was no difference in the immune scores, in the high-risk group, the score of stromal and ESTIMATE was significantly higher (P< 0.05), and the tumor purity was significantly lower in the high-risk group ( Figures 4B–E ). The clinical use of ICIs was guided by comparing the differences in immune checkpoint gene expression levels between the high- and low-risk groups ( Figure 4F ), and many checkpoint genes differed significantly between the two groups, including the potent immunotherapeutic targets PDCD1 (PD-1), CD274 (PD-L1), and CTLA4. To assess the efficacy of immune checkpoint inhibitor treatment, the IPS was used ( Figure 4G ). There were no statistically significant differences in IPS, IPS-CTLA4 blockers, IPS-PD1/PD-L1/PD-L2 blockers, or IPS-CTLA4 and PD1/PD-L1/PD-L2 blockers between the two groups. Interestingly, the TIDE scores were significantly higher in the high-risk group than in the low-risk group ( Figure 4H ). Furthermore, in the high-risk group, the MSI and Neoantigen levels were significantly higher ( Figure 4I, J ).

In CRC patients, we estimated the IC50 of 198 chemotherapy medications or inhibitors and determined the efficacy of the risk score as a predictor of chemotherapy response. The estimated IC50 for Axitinib, Cediranib, Dasatinib, Ibrutinib, Lapatinib, Osimertinib, Rapamycin, and Ribociclib in the high-risk group was significantly lower ( Figures 5A–H ). Additionally, compared to the high-risk group, the IC50 of Camptothecin, Cisplatin, Dabrafenib, Dactinomycin, Docetaxel, Tamoxifen, Gemcitabine, and Oxaliplatin was significantly lower in the low-risk group ( Figures 5I–P ). Finally, the correlation between the risk score and the above drugs was calculated using the Spearman rank correlation coefficient ( Figure 5Q ). These findings show that the risk score is connected to the drug’s sensitivity.

Of the 7 genes in this prognostic model, HAMP was the only gene that was highly expressed in both data sets of CRC and CRLM simultaneously ( Figure 6A, B , Supplementary Table 4 ). To verify the correlation between HAMP expression and CRC, we further explored the expression of HAMP in a variety of CRC cell lines in the CCLE database. HAMP expression was higher in the CRC cell line than in the normal cell line (p = 0.029) ( Figure 6C ). However, between samples of different genders, there was no significant difference in HAMP expression ( Figure 6D ). The expression of HAMP in cancer samples ascended with the worsening of the stage ( Figure 6E ). Moreover, the expression of HAMP in stage I was significantly lower than that in stage II (p = 0.014), III (p = 0.00042), and IV (p = 0.002) ( Figure 6E ). There was no significant expression difference in different MSI ( Figure 6F ). When compared to the high MSI (MSI-H) group, HAMP expression was considerably lower in the microsatellite stability (MSS) group (p = 0.0045) ( Figure 6F ). Demonstrated that HAMP is related to CRC and the progress of cancers. In the TCGA dataset, we divided the samples into the High HAMP group and the Low HAMP group according to the median of HAMP expression. According to the survival study, the high HAMP group had significantly lower rates of OS (P = 0.0023), disease-free interval (DFI) (P = 0.00015), and disease-specific survival (DSS) (P = 0.0031) ( Figures 6G–I ). Indicated that HAMP could function as an independent biomarker for OS prediction. To further verify HAMP expression in the CRC samples, In eight pairs of clinical samples, we employed qRT-PCR and IHC to confirm our findings. The results revealed that HAMP expression was significantly higher in CRC samples than in control samples at both the mRNA ( Figure 6J ) and protein ( Figure 6K ) levels (p< 0.05).

HAMP expression was significantly increased in most cancer tissues compared to normal tissues, including BRCA, COAD, ESCA, GBM, HNSC, KICH, KIRC, KIRP, LUAD, LUSC, and STAD ( Figure 7A ). Tumors regularly have altered DNA methylation states, and the DNA methyltransferase family is primarily responsible for DNA methylation. In the majority of cancer types, HAMP was significantly correlated with the four main methyltransferases ( Figure 7B ). Furthermore, ( Figure 7C ) shows the correlation of HAMP with DNA mismatch repair genes in pan-cancer. TMB of various tumor samples was counted separately. In BLCA, CESC, COAD, LGG, SARC, and THYM, HAMP was associated with TMB positively, but in THCA, HAMP was associated with TMB negatively. ( Figure 7D ). Finally, we noted that HAMP was positively correlated with MSI in COAD, whereas it was negatively correlated with MSI in GBM, LIHC, LUSC, OV, SKCM, STAD, TGCT, and UCEC ( Figure 7E ).

GO annotation findings indicated that HAMP-related genes were primarily involved in tumorigenesis and metastasis (e.g., positive regulation of cell activation, adhesion) and immune response-related processes (e.g., cytokine binding, immune receptor activity) ( Figure 8A ). KEGG analysis showed an association with tumorigenesis and metastasis pathways (e.g., PI3K-Akt signaling pathway, Cell adhesion molecules) and immune activation pathways (e.g., Cytokine-cytokine receptor interaction, Chemokine signaling pathway) ( Figure 8B ). Finally, GSEA was performed, and the same results showed that the tumorigenesis, metastasis pathways, and immune activation pathways were mainly involved ( Figure 8C ).

First, the percentage of immune infiltrating cells in the different HAMP expression groups was compared. HAMP expression was positively connected with macrophages (M0, M1, M2), whereas it was inversely related to dendritic cells activated, T cells CD4 memory resting, mast cells activated, plasma cells, and monocytes ( Figure 9A, B , Supplementary Figure 2 ). The high expression group had substantially higher immune, stromal, and ESTIMATE scores (P< 0.05). In the meantime, the tumor purity in the high-expression group was significantly lower (P< 0.05) ( Figures9C–F ). Clinical use of ICIs was guided by comparing differences in immune checkpoint gene expression levels between the high- and low-expression groups ( Figure 9G ), and many checkpoint genes differed significantly between the two groups, including powerful targets PD-1/L1 and CTLA4. Response to immune checkpoint inhibitor therapy was assessed with the IPS Score ( Figure 9H ). The IPS, IPS-PD1/PD-L1/PD-L2 blockers scores were significantly lower in the high expression group. Additionally, the group with solid expressiveness had significantly higher TIDE scores ( Figure 9I ). The levels of MSI and Neoantigen were considerably more significant in the group with high expression ( Figure 9J, K ). These findings imply an essential link between HAMP expression and immunotherapy.

We estimated the IC50 of 198 drugs or inhibitors in CRC patients to determine the validity of HAMP expression profile as a predictor of chemotherapy response. The estimated IC50 for Axitinib, Cyclophosphamide, Dasatinib, Rapamycin, Ribociclib, Sorafenib, Zoledronate, and SB216763 was noticeably lower in the high expression group ( Figures 10A–H ). Moreover, the estimated IC50 for Trametinib, PD0325901, SB505124, and SCH772984 in the low expression group was significantly lower compared to the high expression group ( Figures 10I–L ). Consequently, the connection of HAMP expression and the IC50 values of the above drugs was calculated using Spearman rank correlation coefficients ( Figure 10M ). These results reflect that the expression of HAMP is correlated with drug sensitivity.

We transfected the CRC cell line DLD-1 with si-RNA of HAMP (si-HAMP), and the qRT-PCR results showed that the expression of HAMP in the si-HAMP group was significantly lower than that in the control group ( Figure 11A ) in mRNA levels (p< 0.05). Then the transwell migration assay of the si-HAMP and control groups revealed that the number of migrated cells in the si-HAMP group was significantly lower than that in the control group ( Figure 11B ), indicating the silencing of HAMP in the CRC cells suppressed the cell migration. The statistical analysis of the transwell migration assay is shown in Figure 11C (p< 0.05).