All microarray datasets were downloaded from the TCGA database. RNA data were downloaded from the TCGA database (https://portal.gdc.cancer.gov/), including 41 control tissues and 421 CRC tissues with clinical data. The GSE14297 dataset was downloaded from the GEO database (http://www.ncbi.nlm.nih.gov/geo). Sixty-eight surgical samples, including tumor tissue and pericarcinomatous tissues, were obtained from patients who were clinically diagnosed at our cancer hospital from October 2020 to June 2021. Samples were collected after approval by the ethics committee of Cancer Hospital Institute.

Based on the mRNA profiles in the TCGA database, we used the “limma” and “edgeR” R packages to identify DEGs (13) between tumor and normal tissues following the criteria |log (fold change)| >1 and P value <0.05. Then, we downloaded immune-related genes from the IMMPORT database (https://www.immport.org/) (14) and explored the IRDEGs that overlapped both the DEGs in the TCGA database and the immune-related genes.

Then, we downloaded TFs from the Cistrome database (http://cistrome.org/) (15) and identified the genes overlapping TF genes and prognostic immune-related genes screened by univariate analysis. After that, we explored the relationship between IRDEGs and TF-DEGs and constructed a TF-IRDEG regulatory network by using Cytoscape software (version 3.6.1). Based on the miRNA and lncRNA matrix files of CRC downloaded from the TCGA database, we identified differentially expressed miRNAs (DEmiRNAs) and lncRNAs (DElncRNAs) with the R package “edgeR” with P<0.01 and |logFc|≥2 as cutoffs. Based on the prognostic genes identified by multivariate Cox proportional hazards analysis, we identified candidate mRNAs as targets of the DEmiRNAs that were recognized by three databases (TargetScan, miRTarBase, and miRDB) and explored the interactions between the DElncRNAs and DEmiRNAs. Based on the interactive regulatory relationships of the DEmiRNAs and DElncRNAs and the DEmiRNAs and DEmRNAs, we utilized Cytoscape software 3.6.1 to construct an lncRNA−miRNA−mRNA−ceRNA network.

We performed univariate and multivariate Cox proportional hazards analyses to screen the IRDEGs to identify hub IRDEGs with prognostic value and built a prognostic risk score prediction model using the “survival” R package and SangerBox online tool. We validated the model through receiver operating characteristic (ROC) curve analysis. We used CIBERSORT to explore the relative immune cell percentages in the hub genes in all tumor tissues (16). Immune information was downloaded from TIMER (17) (https://cistrome.shinyapps.io/timer/) to explore the correlation between the risk score and immune infiltration score through the R packages “limma”, “reshape2”, “tidyverse” and “ggplot2”. In addition, we utilized the online tool SangerBox to explore the relationship between the hub IRDEGs.

We used the “ESTIMATE” and “limma” R packages to estimate the proportions of structural components (stromal and immune cells) in the TME. Then, we calculated the immune score, stromal score, and ESTIMATE score of tumor tissue, divided the three scores into high- and low-score groups, and plotted the survival curves of each group with the “survminer” and “survival” R packages. We used CIBERSORT to explore the relative percentage of immune cells in CRC tumor tissue and used immunohistochemistry to validate the expression on CD4⁺ T cells and M0 macrophages. The antibodies used to detect these immune cells were CD4 rabbit mAb (48274S) and CD68 XP rabbit mAb (76437S), which were purchased from CST. The immunohistochemistry process is shown in the Supplemental Materials and Methods .

We used the “limma” R package to explore the DEGs between the high- and low-score groups based on the immune score and stromal score with the following criteria: |log (fold change)| >1 and P value <0.05. Then, we used the “VennDiagram” R package to identify genes shared in terms of the immune score and stromal score. After obtaining the intersecting DEGs, we performed WGCNA with the “WGCNA” R package to identify clinical trait-related modules (18).

We conducted Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses using the “org.Hs.e.g.db”, “enrichplot”, “ggplot2”, and “GOplot” R packages, and an adjusted P value of <0.05 was considered statistically significant.

We used the “VennDiagram” R package to identify immune hub genes associated with the TME and explored their prognostic value through the Gene Expression Profiling Interactive Analysis (GEPIA) (http://gepia.cancer-pku.cn/) online tool (19). We used the GEO online tool LOGpc (http://bioinfo.henu.edu.cn/DatabaseList.jsp) (20) to explore the protein expression and prognostic value of HAMP. The LOGpc online prognostic tool encompasses 209 expression datasets and provides 13 different survival terms for 31310 patients with 27 distinct malignancies.

We used immunohistochemistry to validate HAMP expression in our clinical samples, and each sample was assessed three times. The antibody used to detect HAMP was the anti-hepcidin+hepcidin-2 antibody (EPR 18937, Abcam), and the details of the process are provided in the Supplemental Materials and Methods . Based on the results above, we calculated the positive cells% and H score to explore the relationship through immunohistochemistry. Positive cells%, which reflects the number of positive cells, was calculated as follows: number of positive cells/total number of cells. The H score is a histological scoring method for immunohistochemistry. The positive number and staining intensity in each section were converted into corresponding values to achieve semiquantitative staining of tissues. The H score was calculated as follows: H score = (PI × I) = (percentage of weak-intensity cells ×1)+(percentage of moderate-intensity cells ×2)+(percentage of strong-intensity cells ×3), where I represents the grade of positive cells. The scores were assigned as follows: negative without staining, 0 points; weak-positive light-yellow staining, 1 point; medium-positive brownish-yellow staining, 2 points; and strong-positive tan staining, 3 points. PI indicates the percentage of positive cells. We used three standard immunofluorescence assessments to validate the results.

We used GSE14297 (21) to explore the potential value of HAMP in the liver metastasis of CRC, and we divided the data into three groups: tumor group, normal group and liver metastasis group. We downloaded the immunological therapy response data of CRC from The Cancer Imaging Archive (TCIA) database (https://tcia.at/) to explore the relationship between the expression of HAMP and therapeutic effects. We performed GSEA with gene sets that were downloaded from https://www.gsea-msigdb.org/gsea/msigdb to explore the biological functions of the hub genes and the risk score in the prognostic model. Then, the hub genes and risk score were divided into two groups: high- and low-expression groups based on the TCGA datasets. We used “c2.cp.kegg.v6.2.symbols.gmt” for analysis and selected the top five pathways. Then, we used the “plyr”, “ggplot2”, “grid”, and “gridExtra” R packages to integrate different significant pathways into a single diagram.

Based on the dataset we downloaded from TCGA-COAD and TCGA-READ, we explored the DEGs (mRNA) between the control tissues and tumor tissues. We used immune genes and DEGs in the TCGA database to explore immune-related genes, including 472 genes.

We used univariate Cox regression analysis to explore survival-related genes ( Figure S1 ) ( Supplementary File 1 ). Based on the genes screened by univariate Cox analysis, we explored TF-related immune genes (PIRDEGs) and constructed a regulatory network between TFs and TF-related immune genes, as shown in Figure 1A . Based on the genes screened by uniCox analysis, we continued to explore hub prognostic genes through multivariate Cox analysis. Thirty-two genes were identified. We explored the regulatory relationships of these genes. There were 9 downregulated DElncRNAs, 27 upregulated DElncRNAs, 1 upregulated DEmiRNA, and 3 candidate target mRNAs (TPM2, AKT3, and FGFR1) in the network. Based on the DEmiRNA-DElncRNA and DEmiRNA-DEmRNA interactions, we built a DEmiRNA-DElncRNA-DEmRNA network ( Figure 1B ).

Based on the survival-related genes identified through multivariate Cox analysis, we calculated the prognostic risk score, and then, we divided the CRC patients into low- and high-risk groups to build a risk score model. The risk score distribution, ROC curve and Kaplan–Meier curves were analyzed and are shown in Figures 2A–C . The area under the curve (AUC) of the ROC curve was 0.85 for 1 year, 0.82 for 3 years and 0.83 for 5 years ( Figure 2B ). As shown in the Kaplan–Meier curves of the relationship of the risk score with overall survival (OS) in the low-risk and high-risk groups, the low-risk-score group had a stronger positive association with OS (P<0.0001; Figure 2C ). We also explored their relationships with each other, and these genes had close relationships ( Figure 2D ). We analyzed the relationship with the risk score of the prognostic model and 22 kinds of immune cells: regulatory T cells, gamma delta T cells, follicular helper T cells, CD8 T cells, naive CD4 T cells, resting memory CD4 T cells, activated memory CD4 T cells, resting plasma cells, resting NK cells, NK-activated neutrophils, monocytes, resting mast cells, activated mast cells, M2 macrophages, M1 macrophages, M0 macrophages, eosinophils, resting dendritic cells, activated dendritic cells, naive B cells and memory B cells. As shown in Figure 2E , the risk score had a close correlation with follicular helper T cells (P<0.01) and had a small correlation with immune cells.

Based on the datasets we downloaded, we estimated the immune and stromal component proportions in the TME through scoring. The ESTIMATE score is a combination of the immune score and stromal score. We analyzed the relationships of these scores with clinical characteristics, including age, sex, TNM classification and stage. For the ESTIMATE score, there was a significant difference between M0 and M1 macrophages. For the immune score, M0 and M1 macrophages also showed significant differences, and the fourth stage was different from the first to third stages. For the stomal score, there was a large difference between the N0 and N2 stages ( Figure S2 ). We divided the scores into two groups: high-score and low-score groups. A high score indicates that there was a large number of immune or stromal components in the TME, and a low score indicates that there was a small number of immune or stromal components in the TME. We used Kaplan–Meier survival analysis to explore the relationship between the survival rate and TME scores. As shown in Figure S3 , all the scores had little correlation with the survival time of the patients. We used CIBERSORT to explore the relative percentage of immune cells in CRC tumor tissue based on the TCGA database, and T cells and macrophages accounted for a large portion of immune cells in CRC ( Figure S4A ). Immunohistochemical results of two kinds of immune cells confirmed this conclusion ( Figure S4B-F ).

To identify the DEGs in the TME, we first explored the DEGs in the stromal and immune components. As shown in Figures 3A, B , there were 14 genes in the low-score group and 729 genes in the high-score group of the stromal and immune groups. To further identify the key modules that were most correlated with CRC clinical traits, we performed WGCNA on the genes between the two groups. Clinical information such as age, TNM stage, and survival time was retrieved from TCGA. By setting a soft-thresholding power of 10, we eventually identified 7 modules ( Figures 3C, D ). From the heatmap of module-trait correlations, we found that the turquoise module, which contained 202 genes, was the most highly correlated with prognosis (P=1.7E-12; Figure 3E ).

We explored the biological functions of the DEGs in the TME and the hub survival-related immune genes through GO and KEGG enrichment analyses following the criteria of P value <0.05 and adj P value<0.05. For the DEGs in the TME, the top five GO enrichment terms were positive regulation of cytokine production (P=2.13E-23), neutrophil activation involved in the immune response (P=5.16E-22), neutrophil degranulation (P=4.07E-12), positive regulation of cell activation (P=3.58E-20) and T-cell activation (P=2.68E-19) ( Figure 4A ). The top five KEGG enrichment pathways were Staphylococcus aureus infection (P=3.98E-21), osteoclast differentiation (P=6.46E-17), phagosome (P=2.83E-15), rheumatoid arthritis (P=1.45E-12) and tuberculosis (P=9.22E-12) ( Figure 4B ).

For hub survival-related immune genes, the top five terms in GO enrichment were regulation of chemotaxis (P=2.08E-06), positive regulation of defense response (P=3.02E-06), positive regulation of chemotaxis (P=4.12E-06), positive regulation of positive chemotaxis (P=9.94E-06) and regulation of positive chemotaxis (P=1.12E-05)) ( Figure 4C ); the top five enriched KEGG pathways were regulation of lipolysis in adipocytes (P=1.64E-05), Ras signaling pathway (P=3.46E-05), ErbB signaling pathway (P=8.57E-05), chemokine signaling pathway (P=0.0001683) and growth hormone synthesis, secretion and action (P=0.0003152) ( Figure 4D ).

Two genes, HAMP and CCR8, overlapped between the survival-related immune genes and the TME-related genes ( Figure 5A ), and we explored their prognostic value, as shown in Figures 5B and S5. As a result, only HAMP was associated with survival time. We validated the prognostic value in CRC in the GEO database ( Figures 5C, D ). We continued to explore the expression of HAMP between tumor tissue and paracarcinoma tissue in clinical specimens and public databases ( Figures 6A, B and S6 ). The results indicated that HAMP is highly expressed in CRC tissue (P<0.001 in the public database). We also explored the clinical features related to this gene ( Figure S7 ). There were no differences across the different stages and M grades. For the T stage, HAMP was significantly differentially expressed between stages T1 and T4 and stages T1 and T3; for the N stage, N0 and N2 showed large differences. Based on the CIBERSORT results, we assessed the relationship between the expression of HAMP and CD4⁺ T cells and M0 macrophages, the results by immunohistochemistry ( Figures 6C, D ) and three-standard immunofluorescence ( Figures 6E, F ) and showed that HAMP had a negative correlation with CD4⁺ T cells and M0 macrophages ( Supplementary Table ).

We first performed GSEA to explore the biological functions of HAMP and the risk score of the prediction model in CRC in the TCGA-COAD and TCGA-READ datasets ( Figures 7A, B ). For the risk score of the prediction model, “amino sugar and nucleotide sugar metabolism”, “apoptosis”, “chemokine signaling pathway”, “intestinal immune network for IGA production” and “natural killer cell-mediated cytotoxicity” were enriched in the low-expression group, as shown in Figure 7A . The top 5 pathways with P<0.05 of the high and low HAMP expression groups are shown in Figure 7B . The high HAMP expression group was enriched in “cell adhesion molecules (CAMs)”, “chemokine signaling pathway”, “cytokine receptor interaction”, “intestinal immune network for IGA production” and “leukocyte transendothelial migration”. The low HAMP expression group was enriched in “lysine degradation”, “nucleotide excision repair”, “peroxisome”, “RNA degradation” and “ubiquitin-mediated proteolysis”. HAMP was also related to the effectiveness of CTLA4 immunotherapy, as shown in Figures 7C, D .

We also explored the value of HAMP in the liver metastasis of CRC. The P value between the liver metastasis and normal groups was <0.0001, and the P value between the liver metastasis and tumor groups was <0.0001 ( Figures 7E, F ).