Data collection and processing

In this study, we integrated bulk RNA-seq, single-cell RNA-seq (scRNA-seq), and spatial transcriptomics data. Bulk RNA-seq data from CRC patients were obtained from The Cancer Genome Atlas (TCGA; https://portal.gdc.cancer.gov/) as well as the GSE1433312 and GSE3883213 (Jorissen et al. 2009, Tripathi et al. 2014) datasets. These bulk transcriptomic datasets were subjected to quality control and normalization, and the corresponding clinical information was curated. Single-cell RNA-seq data were acquired from the GSE17831814 and GSE22585715 (Che et al. 2021, Wang et al. 2023) datasets, encompassing 10 primary CRC patients and 8 patients with colorectal liver metastases (Home - GEO - NCBI). Spatial transcriptomics data for all three cohorts were downloaded from the 10X Genomics official database (https://www.10xgenomics.com/).

Single-cell RNA sequencing data processing and analysis

Single-cell RNA sequencing (scRNA-seq) data from primary CRC and CRLM samples were processed using Seurat (v5.0). Cells expressing fewer than 200 genes, cells with extremely high UMI counts suggestive of doublets, and cells with a high proportion of mitochondrial transcripts (> 20%) were excluded. Potential doublets were further identified and removed using computational doublet-detection algorithms.

To minimize technical variability across patients and datasets, individual samples were first normalized independently and subsequently integrated using the Harmony algorithm. Batch correction was performed based on patient identity and sequencing dataset, generating a harmonized low-dimensional representation while preserving biological variation. Principal component analysis (PCA) was conducted on highly variable genes, followed by Uniform Manifold Approximation and Projection (UMAP) for visualization. A shared nearest-neighbor (SNN) graph was constructed, and graph-based clustering was performed to identify transcriptionally distinct cell populations.

Cell-type annotation was conducted using canonical marker genes reported in previous colorectal cancer single-cell studies. Differentially expressed genes were identified for each cluster and used to validate cluster assignments. Major cell populations, including malignant cells, T cells, B cells, NK cells, macrophages, MDSCs, endothelial cells, fibroblasts, and hepatocytes, were annotated according to established lineage-specific markers.

To infer large-scale chromosomal copy number variations (CNVs) and distinguish malignant from non-malignant cells, inferCNV was applied to the integrated and batch-corrected single-cell dataset. To minimize bias introduced by reference selection, multiple non-malignant cell populations, including immune and stromal cells, were used collectively as reference cells rather than relying on a single reference population. These reference cells were selected based on canonical lineage markers and the absence of malignant transcriptional features. Gene expression values were ordered according to chromosomal position, and inferCNV was performed using sliding-window smoothing and denoising procedures to reduce technical noise and enhance signal robustness. The resulting CNV profiles were visualized as heatmaps to evaluate chromosomal instability patterns across cell populations.

Importantly, inferCNV analysis was used primarily to support malignant cell identification and characterize genomic instability patterns rather than to infer functional mechanisms. Downstream biological conclusions were derived from integrated transcriptomic, trajectory, cell–cell communication, spatial transcriptomic, and experimental validation analyses.

To investigate intercellular communication within the tumor microenvironment, CellPhoneDB was applied to annotated cell populations. Significant ligand–receptor interactions were identified through permutation-based statistical testing, and interaction strengths were aggregated to construct cell–cell communication networks. Pathway-specific analyses focused on MIF, SPP1, TGFβ, and GALECTIN signaling pathways associated with immune suppression and metastatic progression.

Cell differentiation trajectories were reconstructed using Monocle 3. Highly variable genes were used for trajectory inference and pseudotime ordering. Dynamic gene expression changes along differentiation trajectories were analyzed to identify transcriptional programs associated with metastatic progression and PIGT expression.

Functional enrichment analysis

To investigate gene functions and pathway characteristics, functional enrichment analysis was performed on differentially expressed genes or cell type–specific genes. Gene Ontology (GO) and KEGG databases were used to identify significantly enriched biological processes and signaling pathways, highlighting key functional modules. Furthermore, Gene Set Variation Analysis (GSVA) was applied to quantify pathway activity at the single-cell level or across samples, enabling the detection of dynamic changes in critical signaling pathways.

Feature selection and machine learning modeling

In this study, stratified sampling was used to construct training and testing cohorts. Eight supervised learning algorithms (Random Forest, XGBoost, GLMNET, PLS, LDA, Decision Tree, CatBoost, and LightGBM) were trained and evaluated using repeated cross-validation. Feature importance was quantified uniformly across models using permutation-based importance implemented in the DALEX framework. The final top-ranked genes were identified by integrating intra-model mean importance and inter-model consensus, while early stopping and stability resampling were applied to mitigate overfitting and ensure both biological interpretability and statistical robustness.

Tumor immune infiltration analysis

To characterize the composition and functional states of immune cells in the tumor microenvironment, integrated single-cell RNA-seq data were analyzed. Immune cell populations, including T cells, B cells, NK cells, macrophages, and dendritic cells, were first identified and further subdivided into subpopulations based on canonical marker genes, with their proportions and transcriptional profiles quantified. Gene Set Variation Analysis (GSVA) was then applied to score key immune-related pathways, such as cytokine signaling, T cell receptor signaling, and inflammatory pathways, revealing active immune states and potential immunosuppressive mechanisms. Additionally, CIBERSORTx was used to infer immune infiltration from bulk RNA-seq data, validating the robustness of the single-cell findings. The results were visualized using heatmaps and network diagrams to illustrate differences in immune composition and pathway activity across cell populations and samples, providing molecular insights into the tumor immune microenvironment and informing potential immunotherapeutic strategies.

Spatial transcriptomics data processing

Spatial transcriptomics data were obtained from human CRC and matched adjacent mucosal tissues generated using the 10x Genomics Visium HD platform on formalin-fixed paraffin-embedded (FFPE) sections. Raw sequencing data were processed with the Space Ranger pipeline (10x Genomics), which performs alignment to the human reference transcriptome, assigns reads to spatial barcodes, and generates gene expression matrices with corresponding tissue coordinates and histology images.

The count matrices and spatial metadata were imported into Seurat. Low-quality spatial spots were filtered based on sequencing depth, gene detection counts, and mitochondrial gene content. Normalization and variance stabilization were performed using SCTransform to correct for technical variation while preserving biological heterogeneity. Principal component analysis (PCA) was applied for dimensionality reduction, and a shared nearest neighbor (SNN) graph was constructed to identify spatial clusters. Cluster labels were mapped back to the tissue coordinates for visualization of spatial organization.

Spatial domains and putative cell populations were annotated by integrating the spatial transcriptomic data with publicly available single-cell RNA-seq reference datasets from CRC. Reference-guided label transfer and canonical marker genes were used to assign cell types, including epithelial tumor cells, stromal cells, endothelial cells, T cells, B cells, neutrophils, and macrophage subpopulations.

To quantitatively evaluate spatial colocalization between tumor regions and immune cell populations, we employed a neighborhood enrichment analysis. Spatial neighborhoods were defined based on physical adjacency of spots, and the observed frequency of co-occurrence between two cell types or spatial domains was compared against a null distribution generated by randomly permuting cell-type labels across the tissue while preserving spatial coordinates. Enrichment z-scores were calculated to determine whether a given cell population was preferentially localized near tumor-enriched regions, tumor invasive margins, or other spatial domains. Additionally, radial distance-based co-occurrence analysis was performed to characterize the spatial proximity of macrophage subpopulations relative to tumor regions.

Cell migration, invasion, wound-healing, and cell viability assays

Human colorectal cancer cell lines (HCT116) were maintained in DMEM or RPMI-1640 medium supplemented with 10% fetal bovine serum (FBS) and antibiotics under standard conditions (37 °C, 5% CO₂). All cell lines were routinely authenticated and tested negative for mycoplasma contamination.

For the wound-healing assay, cells were seeded into 6-well plates and grown to full confluence, followed by serum starvation for 4–6 h. A uniform linear scratch was generated using a sterile 200-µL pipette tip, and detached cells were removed with PBS. The medium was replaced with 1% FBS to suppress proliferation. Images were captured at 0, 6, 12, 24, and 48 h using phase-contrast microscopy at predefined positions. The wound area was quantified using ImageJ software, and the closure rate was calculated as the percentage of the initial wound area.

Transwell migration assays were performed using 8-µm pore inserts (Corning, 24-well format). A total of 1 × 10⁵ serum-starved cells in 200 µL serum-free medium were seeded into the upper chamber, while 600 µL of medium containing 10% FBS was added to the lower chamber as a chemoattractant. After incubation for 12–24 h, cells that migrated to the lower surface of the membrane were fixed with methanol, stained with 0.1% crystal violet, and counted in at least five randomly selected fields under a light microscope.

For the invasion assay, inserts were pre-coated with Matrigel (BD Biosciences) to mimic the extracellular matrix barrier. Serum-starved cells (1–2 × 10⁵ per insert) were seeded in the upper chamber with serum-free medium, while 10% FBS medium was added to the lower chamber. After 24–48 h, non-invading cells on the upper surface were gently removed with a cotton swab. Invaded cells on the underside of the membrane were fixed, stained with 0.1% crystal violet, and quantified as described for the migration assay.

Cell viability was determined using the CCK-8 assay. Briefly, 3–5 × 10³ cells were seeded per well in 96-well plates and cultured for the indicated durations. Then, 10 µL of CCK-8 reagent (Dojindo) was added to each well and incubated for 1–2 h at 37 °C. Absorbance was measured at 450 nm using a microplate reader.

All experiments were independently repeated at least three times, and results were statistically analyzed to ensure reproducibility and minimize confounding effects from cell proliferation.

Multiplex immunofluorescence staining

We collected paired primary CRC tissues and CRLM specimens from our institution to systematically evaluate the differential expression of PIGT between primary and metastatic sites and to further characterize metastasis-associated alterations in the tumor immune microenvironment. To assess immune infiltration, CD8 (cytotoxic T-cell/immune effector marker), FOXP3 (regulatory T-cell/immunosuppressive marker), and CD163 (M2-like tumor-associated macrophage/immunosuppressive marker) were selected as representative indicators. The levels and relative proportions of immune effector and immunosuppressive cell infiltration were compared between primary tumors and liver metastases.This study was approved by the Institutional Ethics Committee of our hospital, and written informed consent was obtained from all participating patients.

Statistical analysis

All quantitative data are presented as mean ± standard deviation (SD) or median with interquartile range (IQR), depending on the data distribution. Differences between two groups were assessed using Student’s t-test or the Mann–Whitney U test, while comparisons among multiple groups were performed using one-way analysis of variance (ANOVA) or the Kruskal–Wallis test, with multiple comparison corrections applied as appropriate. Correlation analyses were conducted using Pearson or Spearman correlation coefficients, selected according to variable type and distribution. Survival analysis was performed using the Kaplan–Meier method with log-rank testing. All statistical analyses were conducted in R (version 4.3.3), and two-sided P values < 0.05 were considered statistically significant.