Gene expression datasets (GSE87211, GSE14333, GSE117536, GSE17537, GSE17538, GSE38832, and GSE103479) were obtained from the Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/) as of October 2023. Single-cell RNA sequencing (scRNA-seq) data were retrieved from GEO accession GSE132465. Bulk RNA-seq expression profiles and associated clinical data for rectal adenocarcinoma (READ) and colon adenocarcinoma (COAD) patients were acquired from The Cancer Genome Atlas (TCGA; https://portal.gdc.cancer.gov/). The TCGA pan-cancer dataset, comprising over 10,000 samples across 33 cancer types, was also included for external validation. A total of 10,598 macrophage polarization-related genes (MPGs) were retrieved from GeneCards (https://www.genecards.org/) using the keyword “macrophage polarization.” Mutation status of candidate genes was assessed using cBioPortal (https://www.cbioportal.org/) (32). Protein-level expression data of MPGs in normal and tumor tissues were accessed from the Human Protein Atlas (https://www.proteinatlas.org/) (33). Gene expression matrices were normalized using the NormalizeBetweenArrays function from the limma R package, and batch effect adjustment was performed using the “ComBat” algorithm from the sva package, with default parameters, to correct for potential batch effects across datasets. GSE87211 was designated the training cohort, while the TCGA-READ dataset served as the test cohort.

Weighted gene co-expression network analysis (WGCNA) was performed on the GSE87211 dataset using the WGCNA R package (34). An optimal soft-thresholding power (β) was selected to ensure scale-free network topology (power = 6, minimum module size = 30, with a module merging threshold of 0.25). An adjacency matrix was constructed and transformed into a topological overlap matrix (TOM) to measure gene connectivity. Genes were hierarchically clustered based on TOM dissimilarity, and distinct gene modules were identified using average linkage clustering. Module–trait correlations were calculated to identify modules most associated with clinical traits in the GSE87211 cohort. These modules were prioritized for downstream analysis.

Differential expression analysis was conducted using the limma package (35) in R. Genes with |log2FoldChange| ≥ 0.5, and adjusted p-value< 0.05 were considered differentially expressed. Differentially expressed genes (DEGs) were intersected with the curated macrophage polarization genes (MPGs) to identify a subset of differentially expressed MPGs (DEMPGs) relevant to rectal cancer.

Univariate Cox regression was performed to identify DEGs significantly associated with overall survival in the training and validation cohorts. Candidate hub genes were identified by intersecting WGCNA module genes, DEGs, and MPGs with prognostic significance. To refine this gene set, four machine learning algorithms—LASSO Cox regression (glmnet package), support vector machine (SVM; conducted using the “e1071” R package; Kernel function: Recursive Feature Elimination (RFE kernel)) (36), random forest (RF; 500 trees with default settings from the “randomForest” R package), and extreme gradient boosting (XGBoost; xgboost package) (37)—were applied. Genes selected by all four methods were defined as core DEMPGs for further modeling. The cross-validation strategy for machine learning models: 10-fold cross-validation, repeated 3 times to ensure model stability, using the “caret” R package

Multivariate Cox regression was applied to the core DEMPGs to construct a macrophage polarization gene signature (MPGS). The risk score for each patient was calculated as:

where βi represents the regression coefficient, and χi is the normalized expression value (FPKM) of each signature gene. Patients were stratified into high- and low-risk groups based on the median risk score. Kaplan–Meier survival curves and multivariate Cox models were used to evaluate the prognostic significance of MPGS, adjusting for clinical covariates. Receiver operating characteristic (ROC) curves were generated using the timeROC package (38), and calibration curves were plotted to compare predicted and observed survival. Validation was performed in the independent TCGA-READ dataset.

DEGs between high- and low-risk groups (|log2FoldChange| ≥ 0.5, and adjusted p-value< 0.05) were identified using limma. Gene ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and gene set enrichment analysis (GSEA) were conducted using the clusterProfiler package (39) to explore biological processes and pathways enriched in the high-risk group, false discovery rate (FDR)< 0.5 and normalized enrichment score (NES) > 1 were set at the cut-off criteria. Immune cell composition and tumor microenvironment (TME) scores were assessed using the CIBERSORT, QuanTIseq, and single-sample GSEA (ssGSEA) algorithms.

Protein-protein interaction (PPI) networks for hub genes were constructed using GeneMANIA (http://www.genemania.org) (40, 41), an integrative platform that incorporates data on co-expression, physical interaction, co-localization, and functional annotations. Functional enrichment analysis was conducted to identify biological processes and pathways potentially regulated by the candidate genes.

Drug response analysis was performed using the prophetic and oncoPredict R packages (42, 43). Half-maximal inhibitory concentration (IC50) values for various chemotherapeutic agents were predicted and correlated with the MPGS risk scores. Differences in drug sensitivity between risk groups were visualized using scatter plots, highlighting drugs with significant IC50 variation. The TIDE scores were calculated utilizing the Tumor Immune Dysfunction and Exclusion (TIDE, http://tide.dfci.harvard.edu/login/) database (44, 45). Moreover, immunophenoscore (IPS) of GC patients were obtained in The Cancer Immunome Atlas (TCIA, https://tcia.at/home) database (46).

Single-cell RNA-seq data were processed using Seurat v4.3.0 for quality control, normalization, and dimensionality reduction. Cells with fewer than 400 genes or mitochondrial content exceeding 20% were excluded. Doublets were removed using DoubletFinder v2.0.3. Integration across samples was performed with Harmony v1.2.3. Principal component analysis (PCA) and Uniform Manifold Approximation and Projection (UMAP) were used for dimensionality reduction and clustering and the top 30 PCs were retained for downstream analysis. Downstream analyses were based on integrated expression matrices.

Cell type annotation was performed using the Seurat FindAllMarkers function to identify cluster-specific marker genes (adjusted p< 0.05, min.pct > 0.25, |log2FC| > 0.25). Initial annotations were derived using the SingleR package and cross-validated against the CellMarker database. Manual curation was performed to confirm annotations based on canonical gene expression profiles from the literature.

Fifty hallmark gene sets were downloaded from the Molecular Signatures Database (MSigDB v7.5.1; https://www.gsea-msigdb.org/gsea/msigdb) (47) Metabolic activity scores were calculated at the single-cell level based on the mean scaled expression of all genes within each signature, as previously described (48). Differential expression of pathway scores between tumor and normal tissues was assessed using the FindAllMarkers function in Seurat, with an adjusted p-value threshold of< 0.05.

Chromosomal copy number variations (CNVs) were inferred using the R package “inferCNV”. Epithelial cells in normal tissues served as reference populations. For each cell subcluster, CNV scores were calculated by aggregating the CNV levels of all constituent cells. The threshold parameter was set to 0.1, while other settings remained at default.

Cell–cell communication analysis was carried out using the R package “CellChat” (version 1.1.3). To ensure consistent sampling across cell subclusters, 500 cells were randomly selected from each subpopulation using the subset function. The analysis incorporated three major signaling categories from the CellChat database: Secreted Signaling, ECM–Receptor, and Cell–Cell Contact. A minimum threshold of 10 cells per cluster was applied to filter out low-abundance populations (49).

A total of 40 paired rectal adenocarcinoma (READ) and adjacent normal tissue samples were collected from patients undergoing surgical resection at Yangpu Hospital, Tongji University, between November 2018 and November 2019. All procedures were approved by the Ethics Committee of Yangpu Hospital (Approval No. LL-2023-LW-012). Fresh specimens were fixed in 4% paraformaldehyde for immunohistochemistry and snap-frozen in liquid nitrogen for RNA and protein extraction.

Total RNA was extracted from paired tumor and adjacent tissues using TRIzol reagent and reverse-transcribed into cDNA using a commercial kit (Takara, Dalian, China). Quantitative real-time PCR (qRT-PCR) was performed using gene-specific primers ( Supplementary Table 1 ).

For protein extraction, tissues were lysed in RIPA buffer (Solarbio, China) with protease inhibitors (1:100, Thermo Scientific). Western blotting used the following primary antibodies: PYGM (1:1,000; ProteinTech, 19716-1-AP), β-actin (1:4,000; ProteinTech, 66009-1-Ig), Arg1(1:1,000; ProteinTech, 16001-1-AP), CD301(1:1,000; ProteinTech, 13590-1-AP), CD206(1:1,000; ProteinTech, 32647-1-AP), IL-10(1:1,000; ProteinTech, 60269-1-Ig).

Formalin-fixed, paraffin-embedded tissue blocks were sectioned at 4 μm thickness. Sections were dewaxed, rehydrated, and subjected to antigen retrieval using a pressure cooker for 30 minutes. Endogenous peroxidase activity was blocked using 3% hydrogen peroxide for 20 minutes. Non-specific binding was minimized with 5% BSA for 40 minutes. Sections were incubated overnight at 4°C with anti-PYGM primary antibody (1:100; ProteinTech, 19716-1-AP). Visualization was achieved using a DAB detection kit and counterstaining with hematoxylin.

Three human colorectal cancer (CRC) cell lines—HCT116, LOVO, and SW620—and the normal colonic epithelial cell line NCM460 were obtained from the Shanghai Institute of Biochemistry and Cell Biology. All cell lines were maintained in DMEM supplemented with 10% fetal bovine serum (FBS; Gibco, Carlsbad, CA, USA) at 37°C in a humidified incubator with 5% CO₂. THP-1 cells were cultured in RPMI-1640 medium (Gibco) supplemented with 10% fetal bovine serum (FBS). Lipofectamine 3000 (Invitrogen, Carlsbad, CA, USA) was used to transfect cells with an siRNA specific for PYGM and a control construct purchased from GeneChem (Shanghai, China) ( Supplementary Table 2 ). Cells were utilized for downstream assays at 48h post-transfection. Analyses were conducted in triplicate. PYGM overexpression plasmid was customized from GenePharma (Shanghai, China).

Migration and invasion assays were performed using 24-well Transwell chambers (Nest, China). Cells were seeded in serum-free DMEM (250μL) into the upper chamber, and 600μL of DMEM with 10% FBS was added to the lower chamber. For invasion assays, inserts were pre-coated with Matrigel (2mg/mL). After 24hours, non-migrated/invaded cells were removed, and cells on the lower membrane surface were fixed with 4% paraformaldehyde and stained with crystal violet for 10 minutes. Cells were quantified in five non-overlapping fields under a microscope (Nikon, Japan).

Confluent cells were scratched using a 10 μL pipette tip for the wound healing assay and cultured in a serum-free medium. Images were taken at 0 and 24hours using phase-contrast microscopy to assess wound closure.

To assess the rate of DNA synthesis, CRC cell lines were subjected to treatment with 5-ethynyl-2’-deoxyuridine (EDU) at a concentration of 50 μM, which was subsequently added to the cell culture plates. Following a 30-minute incubation, DNA was stained using Hoechst 33342, allowing for the visualization of positively stained cells under a microscope. HCT116 and SW620 cells, characterized by either PYGM overexpression or knockdown, were dissociated into single-cell suspensions using 0.25% trypsin. These cells were then stained with Annexin V-APC and 7-Aminoactinomycin D (7-AAD) to evaluate apoptosis rates. THP-1 monocytes were first differentiated into macrophage-like cells by treatment with 200 ng/mL phorbol 12-myristate 13-acetate (PMA) for 48 hours. Following differentiation, the PMA-treated THP-1 cells were gently washed with PBS to remove residual PMA and were then seeded into the chamber of the transwell system for indirect co-culture.

After 48 hours of co-culture, THP-1-derived macrophages were collected and stained with anti-CD301-APC and anti-CD206-APC, and the number of CD301 or CD206-positive cells in macrophages was analyzed by flow cytometry. Meanwhile, total RNA and protein of THP-1-derived macrophages were extracted for the detection of M2 macrophage markers.

All statistical analyses and visualizations were performed in R (v4.2.1). Visualization packages included ggplot2, ggpubr, and enrichplot. For comparisons between groups, the Wilcoxon rank-sum test was applied. A two-sided p-value< 0.05 was considered statistically significant.

The workflow for model construction and downstream analyses is summarized in Figure 1 . Weighted Gene Co-expression Network Analysis (WGCNA) was applied to the GSE87211 cohort to explore gene modules associated with rectal cancer. A soft-thresholding power of β = 18 was selected to ensure scale-free network topology ( Figure 2A ). Gene clustering yielded multiple expression modules, visualized via a dendrogram, each represented by a distinct color ( Figures 2B, C ). Among them, the dark red, dark grey, and brown modules demonstrated the strongest correlations with clinical traits (Pearson’s r = 0.88, –0.88, and 0.76, respectively; Figure 2D ). The grey module, comprising unassigned genes, was excluded. Differential expression analysis identified 3,555 downregulated and 3,750 upregulated genes in the GSE87211 cohort; integration with TCGA-READ yielded 5,961 downregulated and 356 upregulated genes, visualized as volcano plots and heatmaps ( Figures 2E, F ).

Univariate Cox regression revealed 734 and 632 survival-associated genes in the GEO and TCGA datasets. By intersecting prognostic DEGs, WGCNA-derived module genes, and macrophage polarization genes (MPGs) from GeneCards, 29 candidate genes were identified ( Figure 2G ).

LASSO regression narrowed this set to 11 genes (BRCA1, FAR1, GPSM2, KL, MAOB, POLA1, PTPRU, PYGM, SYP, TIMP1, TMOD1; Figure 2H ). Three additional machine learning methods—XGBoost, SVM-RFE, and Random Survival Forest (RSF)—identified the top 15 genes by feature importance. The intersection with LASSO output resulted in six shared hub genes ( Figures 2I, J ).

Multivariate Cox regression was performed on the six hub genes to refine the candidate genes. Three genes—TIMP1, MAOB, and PYGM—remained significant (p< 0.05) in both the GSE87211 and TCGA-READ cohorts ( Supplementary Figures S1A–S1F ).

These genes formed the basis of a prognostic risk model (MPGS), and the following formula was derived:

Patients in both cohorts were stratified into high- and low-risk groups based on the median risk score. Risk score distribution and corresponding survival status are displayed in Figure 3A , and gene expression heatmaps showed upregulation of all three genes in the high-risk group. Kaplan–Meier analysis confirmed that low-risk patients had significantly better overall survival (OS) in both datasets ( Figure 3B ). Receiver Operating Characteristic (ROC) analysis demonstrated good predictive performance of the MPGS for OS in both training and validation cohorts ( Figure 3C ).

Univariate and multivariate Cox regression analyses were conducted to evaluate whether the MPGS was an independent predictor of OS. In the GSE87211 dataset, both MPGS (p<0.001) and M (P=0,006) stage were significantly associated with OS in univariate analysis, and MPGS remained independently prognostic in the multivariate model ( Figure 3D ). Similar findings were confirmed in the TCGA-READ cohort ( Figure 3E ).

A prognostic nomogram integrating MPGS and clinical variables was developed to predict 1-, 3-, and 5-year survival probabilities. Calibration curves demonstrated good agreement between predicted and observed outcomes in both cohorts ( Figures 3F, G ). Stratified survival analysis within M-stage subgroups revealed that high-risk patients in the M0 group exhibited significantly poorer OS than low-risk patients, while no significant difference was observed in the M1 subgroup ( Figures 3H–I , Supplementary Figure S2 ). Additionally, both datasets’ Kaplan–Meier and ROC analyses showed that PYGM and MAOB exhibited strong diagnostic and prognostic performance ( Supplementary Figure S3 ).

To investigate the biological implications of the MPGS, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were conducted. GO analysis highlighted immune- and metabolism-related processes, including positive regulation of the MAPK cascade, macrophage activation, and epithelial cell proliferation. KEGG pathway analysis further identified enrichment in the TGF-β signaling pathway, oxidative phosphorylation, and other metabolism-associated pathways in both cohorts ( Figures 4A, B ).

GSEA showed that the genes in the high-risk group from both cohorts were significantly enriched in several hallmark pathways, including Jak-Stat signaling pathway, MAPK signaling pathway (KEGG), glycosaminoglycan metabolism, interleukin 4 and interleukin 13 signaling (Reactome), epithelial-to-mesenchymal transition in colorectal cancer, and the PI3K-AKT signaling pathway (WikiPathways) ( Figures 4C, D ). Additionally, several other pathways related to cancer progression, macrophage polarization, and metabolism were enriched in the GSE87211 and TCGA cohorts, respectively ( Figures 4E, F ).

Immune infiltration in the tumor microenvironment (TME) was assessed using CIBERSORT and QuanTIseq algorithms, revealing that M2 macrophages were enriched in the high-risk group ( Figures 5A, B )

We analyzed the IPS of patients stratified by MPGS ( Figure 5C , Supplementary Figures S4A, B ) and the expression levels of three core genes using data from the TCIA database. The results showed that patients with low MPGS and low PYGM expression exhibited significantly higher IPS values compared to those with high expression levels, suggesting that lower MPGS and PYGM expression may be associated with improved responsiveness to immunotherapy.

Subsequently, multiple datasets were analyzed using the TIDE algorithm to evaluate the immunotherapy response between high- and low-MPGS expression groups. In the TIDE model, higher scores indicate a greater likelihood of immune evasion and a lower probability of benefiting from immune checkpoint inhibitor (ICI) therapy (44). The analysis revealed that the high-MPGS group exhibited significantly higher TIDE and Exclusion scores ( Figures 5D–I ), suggesting reduced sensitivity to immunotherapy. Consistently, PYGM showed a similar trend in both the training and validation cohorts ( Supplementary Figures S4C, D ).

Drug sensitivity prediction using the oncoPredict and prophetic packages indicated that IC50 values for Cytarabine and Docetaxel were significantly lower in high-risk patients, suggesting higher drug sensitivity while Lenalidomide, GW441756, Bosutinib, Afatinib, Gefitinib, Methotrexate, and GW441756 shown the opposite trend ( Figure 5J ). These findings may inform patient stratification and personalized chemotherapy, pending clinical validation.

To extend the macrophage polarization gene signature (MPGS) ‘s prognostic utility, we assessed its performance in five additional GEO datasets containing survival information. Patients with higher MPGS-derived risk scores exhibited significantly worse overall survival in all cohorts. ROC curve analyses confirmed consistent predictive performance ( Figures 6A–F ).

Subsequently, we applied the MPGS to the TCGA pan-cancer dataset (n >11,000; 33 cancer types). Higher MPGS scores were significantly associated with worse disease-free survival (DFS), disease-specific survival (DSS), progression-free interval (PFI), and overall survival (OS) ( Figures 6G–J ). Furthermore, MPGS risk scores varied significantly across clinical stages ( Figure 6K ).

When evaluating individual cancer types, a hazard ratio (HR) > 1 for MPGS was observed in eight malignancies—BLCA, COAD, GBM, KIRC, LGG, LUSC, SKCM, and STAD—indicating a potential risk association. In contrast, MPGS was inversely associated with mortality (HR< 1) in 12 cancer types, suggesting a protective trend ( Supplementary Figure S5 ).

To further investigate the immune associations of the three hub genes (PYGM, MAOB, TIMP1), we performed immune cell infiltration analyses across 33 TCGA cancer types using ssGSEA and CIBERSORT algorithms. ssGSEA showed a positive correlation between hub gene expression and macrophage infiltration in multiple cancer types, including BLCA, COAD, ESCA, HNSC, LGG, LUSC, PCPG, PRAD, READ, SKCM, and UVM ( Figure 7A ). CIBERSORT analysis further revealed that M2 macrophages exhibited a consistent positive correlation with hub gene expression in cancers such as BLCA, READ, and TGCT, while M1 macrophage correlation was limited, observed mainly in LGG ( Figure 7B ). PYGM expression was also positively associated with canonical M2 macrophage marker genes ( Supplementary Figure S6 ).

Genomic profiling of MPGS genes using cBioPortal showed varying degrees of mutation frequency across cancer types ( Figure 7C ). Protein-protein interaction (PPI) analysis via GeneMANIA revealed that the hub genes, particularly PYGM, are functionally linked to glucose catabolism ( Figure 7D ). Matrix metalloproteinase 9 (MMP9) exhibited the strongest interaction, consistent with its established role in tumor invasion and metastasis (50).

To examine the relationship between MPGs and tumor metabolism at single-cell resolution, scRNA-seq data from 23 rectal cancer patients and 10 healthy donors were analyzed. After integration and quality control, 63,689 cells were retained. Based on canonical markers, seven major cell types were identified: plasma cells (TNFRSF17), B cells (CD79B, MS4A1), T cells (CD3D, CD3E), epithelial cells (KRT18, EPCAM), myeloid cells (CD68, LYZ), fibroblasts (ACTA2, TAGLN), and endothelial cells (PLVAP) ( Figures 8A , 9C ). Cell proportion analysis showed an increased abundance of epithelial and myeloid cells in tumor tissues ( Figure 8D ). Metabolic pathway enrichment analysis using hallmark gene sets from MSigDB revealed significant metabolic activation in epithelial cells, fibroblasts, and myeloid populations ( Figure 8E ).

Myeloid cells were further sub-clustered into 16 subsets and annotated into seven cell types, including M1/M2 macrophages, monocytes, cDC1, cDC2, and pDCs, based on marker expression and Spearman correlation with established cell-type profiles ( Figures 8F–I ). Tumor samples exhibited elevated M2 macrophages and monocytes ( Figure 8J ). Pathway analysis indicated significant metabolic enrichment—including glycolysis, fatty acid metabolism, and oxidative phosphorylation—in M2 macrophages ( Figure 8K ).

Given the significant enrichment of metabolic pathways in the epithelial cell subcluster, we applied the “inferCNV” R package to analyze this subcluster ( Figure 8L ), distinguishing malignant from normal epithelial populations. The results revealed that PYGM expression was markedly elevated in malignant epithelial cells ( Figure 8M ). Subsequent analyses of cell-cell communication and intercellular interactions demonstrated a strong association between malignant cells and M2 macrophages ( Figures 8N, O ). These findings suggest that PYGM, as a metabolism-related gene predominantly expressed by malignant epithelial cells, may play a regulatory role in modulating M2 macrophage activity within the tumor immune microenvironment, thereby contributing to tumor progression.

Due to the prominent biological functions and immunotherapy-specific relevance of PYGM, it was selected for further validation. Immunohistochemistry (n = 21), Western blotting (n = 26) and qRT-PCR(n=40) further confirmed increased PYGM protein expression in READ tissues compared to adjacent normal tissues ( Figures 9A–D ). In vitro, assays also demonstrated higher PYGM expression in colorectal cancer cell lines (HT29, SW620, LOVO) compared to normal epithelial cells (NCM460) ( Figures 9F, G ). Besides, the expression patterns of MAOB and TIMP1 were determined using qRT-PCR analysis of clinical samples and HPA database ( Supplementary Figure S7A, B ).

Logistic regression analysis revealed that elevated PYGM expression was significantly associated with advanced T stage (p = 0.03), tumor anatomical subdivision (p = 0.003), and male sex (p = 0.009) ( Figure 9H ). TCGA and clinical validation cohorts confirmed higher PYGM levels in patients with advanced stage and males, which was inconsistent in patients in COAD ( Figures 9I–J , Supplementary Figure S7C ). Survival analysis further validated that higher PYGM expression was associated with poor OS ( Figure 9K ).

In order to explore the functions of PYGM in RC, it was knocked down by siRNA in SW620 and overexpressed in HCT116, and the efficiency was verified by western blotting ( Figure 10A ).

EDU results showed that SW620 had a decreased proliferation, while HCT116 had an increased viability ( Figure 10B ). Furthermore, we detected the effect of PYGM on CRC apoptosis, and our study indicated that overexpression of PYGM significantly reduced the apoptosis rate of CRC cells, while knockdown of PYGM significantly increased the apoptosis rate of CRC cells ( Figure 10C ). The wound healing assay showed a marked decrease in cell migration following PYGM knockdown ( Figure 10D ) and an increase after its overexpression. Consistent with the results, transwell assays verified that PYGM knockdown inhibited SW620 invasion and migration, and its overexpression in HCT116 had the opposite trend ( Figures 10E, F ). THP-1 monocytes were induced into M2 macrophages using PMA and IL-4. PYGM-overexpressing or -silenced HCT116 and SW620 cells were co-cultured with these macrophages ( Figure 10G ). Flow cytometry revealed that CD206 and CD301 expression were elevated following PYGM overexpression, but reduced upon PYGM knockdown ( Figure 10H ). Similarly, RT-PCR showed corresponding changes in Arg1, CD206, CD301, and IL-10 mRNA levels in THP-1-derived macrophages ( Figure 10I ), which was further validated by western blotting ( Figure 10J ).