First, we searched PubMed for studies on single-cell RNA sequencing of colorectal cancer from January 2017 to January 2022. According to the study design, only scRNA sequencing of whole tumor tissue was included, while studies analyzing flow cytometry isolated T cells or myeloid cells were excluded. To minimum the discrepancies and batch effect across sequencing platforms, only data generated from the 10×Genomics sequencing studies were included. Additionally, studies that provided raw sequencing data (including fastq file or raw reads count matrix) were selected, while studies only providing converted datasets (such as TPM files) were excluded. Finally, 78 samples from 3 studies were included in our study (GSE188711, GSE144735, GSE132257) ( Supplementary Figure 1A ).

Generally, the scRNA-seq data were performed according to the standard protocols of Seurat (version 3.0) (38). For each sample, the gene and count features were identified, and cells with less than 200 or more than 6000 features were filtered. Then, cells with mitochondrial RNA percentage > 15 or with ribosome RNA percentage < 3 were further removed. Then, the Doublet Analysis was performed using the DoubletFinder (version 2.0.3) (39) with default settings for each sample, and all potential doublets were removed. Subsequently, the scRNA datasets were merged into a larger Seurat file. The integrated data were normalized, scaled and processed for PCA analysis. The data were further visualized using the t-SNE and UMAP methods, respectively. The JackStraw method was used to identify the final PCA components for further analysis. The number of cell clusters were evaluated using the shared nearest-neighbor modularity optimization-based clustering algorithm set a resolution of 0.01, 0.05, 0.1, 0.3, 0.4, 0.5, 0.6, 0.7, 0,8, 0.9, 1.0 gradients, and the threshold of 0.01 was chosen, and all cells were divided into 9 highly distinct subpopulations. Similarly, in the subgroup analysis of T cells and myeloid cells, the resolution value was set as 0.3 and 0.2, respectively. The FindAllMarkers function was used to identified the marker genes of each cluster. To assign cell identities and lineages, differentially expressed genes from each cluster were analyzed to sets of previously reported cell-type markers. These 9 subgroups can be well characterized as nine major cell types in tumor. The harmony software (version 0.1.0) (40) was used to evaluated and remove the batch effect in each subgroup. Subsequently, the cells in each group were reanalyzed according to the standard protocol of Seurat. Differentially expressed genes across different groups were identified using the wilcoxauc function in presto package. The clusterProfiler package (41) was used for KEGG and GO pathway analysis.

We used Monocle2 (version 2.24.0) (42) and Monocle3 (version 1.0.0) (43) to determine the trajectory of T cells, myeloid cells, and epithelial cells. The data processing and pseudo-temporal ordering were all performed using Monocle 2, and graph-based trajectory inference function from Monocle 3 was used to generate the trajectory tree. Cell clustering and annotation were performed using Seurat (version 4.1.1) according to the standard protocol. The presto package wilcoxauc function was used to calculate the differential expressed genes of each cluster. The Clusterprofiler package was used for further GO or KEGG enrichment analysis. The dataset was then imported into Monocle2 for pre-processing, including quality control and selection of highly variable genes. For different cell clusters, the top 500~1000 highly variable genes were used for further analysis. The DDRTree (version 0.1.5) method is used for dimension reduction. Cell trajectory was visualized according to cell subtypes and cell states. The Basic Differential Analysis algorithm is used to identify differently expressed gene according to pseudotime function. BEAM (Branched Expression Analysis Modeling) was used to identify genes that are regulated in a branching-dependent manner. The Seuratwrapper package for Monocle 3 was used to cluster cells and develop learned graphs. In addition, Monocle3 was also used to construct the single-cell trajectories for macrophages. The function learn_graph was run with default parameters.

We employed the inferCNV (version 1.12.0) software (44) to determine the copy number variation (CNV) in each cell. The inferCNV algorithm is an analysis tool for single-cell copy number variant analysis. The core idea of the algorithm is to compare the genome-scale gene expression intensity of each cell from the experimental group to the average gene expression intensity of reference ‘normal’ cells, and the relative expression intensity of genes on each chromosome is obtained and presented in a heatmap. Massively over-abundant or less-abundant genes can be observed in malignant cells as compared to that of normal cells. The cells of interest and reference cells can be further integrated and analyzed using unsupervised clustering analysis, and the potential malignant in the cells of interest can be identified by the results of unsupervised clustering. First, we extracted the epithelial cell subsets, and used epithelial cells from normal colorectal tissues as reference to determine single-cell CNV of epithelial cells from tumor and adjacent tissues. Next, based on the CNV data, we used unsupervised clustering to classify all epithelial cells into 7 subgroups. We distinguished normal and malignant epithelial cells in tumor and adjacent tissues based on the unsupervised clustering data (groups of cells clustered with normal tissue epithelium without significant copy number amplification or deletion). To verify the results, we next calculated a CNV_score of each cell. We found that cell subgroups 4 and 7 had the lowest CNV_scores, indicating that they were most likely normal epithelial cells. Then, we generated heatmaps according to the CNV data. A hidden Markov model (HMM) was used to identify potential CNV sites, and the output was denoised using a Bayesian latent mixture model. The phylogenetic tree of tumor evolution was draft using Uphyloplot2 software (45) based on the CNV data, and the phylogenetic tree of each sample was classified according to tumor location and mutation status.

The CellphoneDB (46) algorithm and R ‘iTALK’ (version 0.1.0) (47) package were used to evaluate the crosstalk across the immunological and non-immunological cellular components in normal mucosa, tumor adjacent and tumor tissue. The quantification of cell-cell communications by calculating the expression of immune checkpoint genes using ‘iTALK’ has been reported in previous studies (48, 49). The expression level of immune checkpoints is a key factor in the remodeling of tumor immune microenvironment and immunotherapy response, and the quantile or custom-defined values of PD-L1 expression were associated with patient outcomes after immunotherapy in lung cancers and colorectal cancers (50, 51). Thus, in our analysis, the ligand or receptor expression density below the first quartile value was thought to be a false positive connection.

The Htseq-counts matrix and clinical data of the TCGA-COAD cohort (n=478) were acquired using the TCGAbiolinks package (version 2.20.1). Subsequently, a custom-drafted script was used to transform the Htseq-counts matrix to TPM (Transcripts Per Kilobase of Exon Model per Million mapped reads) data. The canonical gene expression marks of 24 different types of tumor-infiltrating immune cells as well as cell marker identified using scRNA-seq data were set according to a prior publication with minor modifications. Then, the R GSVA package (version 1.40.1) was used for single-sample gene set enrichment analysis. The R ConsensusClusterPlus Package (version 1.58.0) and factoextra package (version 1.0.7) were used for unsupervised clustering analysis.

The standard protocol of the multi-color immunofluorescent staining has been reported in our previous studies. Generally, the colorectal samples were cut into 5mm×5mm pieces, and fixed overnight in 4% PFA. The fixed tissues were embedded in paraffin, and sliced into sequential 8-μm thick slides. Immunofluorescent-staining was performed according to a standard protocol on paraffin-embedded slides. Multicolor immunofluorescence staining was performed using the Tyramide SuperBoost Kits (Invitrogen, Cat#B40912-40926). The Panoramic MIDI was used to scan the IF-stained slides (3DHISTECH Digital Pathology Company, Budapest, Hungary).

All statistical analyses were performed using R software (version 4.1.0) or GraphPad Prism 8.0 (GraphPad software Inc.). The comparison of continuous variables between groups was performed using student t-test or welch t-test. The comparison of continuous variables across multiple groups were was performed using t-test with Bonferroni or Dunnett adjustment. The comparison of enumeration data was performed using chi-square test. The survival curve was drafted using Kaplan-Meier methods, and log-rank test was used to compare the survival difference across different groups. Two-sided p-value less than 0.05 was considered statistically significant.

We analyzed single-cell RNA (scRNA) sequencing data from 78 treatment-naive samples, including 42 colorectal cancer, 13 tumor-adjacent, and 23 normal mucosal samples ( Figures 1A, B ). The KRAS gene primarily included RAS activating mutations such as KRAS^G12D and KRAS^G12C . Whereas the TP53 gene had more diverse mutations sites. In the subsequent analysis, the tumors were categorized in to 4 groups, including KRAS (tumors only bear KRAS mutations), DUAL (tumors bear KRAS and TP53 dual mutations), TP53 (tumors only bear TP53 mutations) and WT (neither KRAS nor TP53 are mutated). The patient cohort information was provided to identify the origins of the patients. The detailed clinical information of the samples was presented ( Supplementary Table 1 ). After quality control, a total of 132,618 cells with 3,827 median gene count per cell were retained. Then, we categorized all cells into 9 distinct subgroups ( Figure 1C ) through a standard Seurat pipeline. As expected, diverging distribution of cells within each subgroup was observed when depicting the UMAPs based on cohort, tumor location and sample identity ( Figure 1C ). Using previously reported canonical cell markers, the cell subgroups can be identified as T lymphocytes (n= 43,740), Epithelial cells (n= 27,505), Myeloid cells (n=14,680), Fibroblasts cells (n=14,625), Plasma cells (n=13,416), B cells (n=11,971), Endothelial cells (n= 4,082), Master cells (n=1,406), and NK-like cells (n=1,193) ( Supplementary Figure 1B ). The distribution of the cells in normal mucosa, tumor adjacent and tumor tissues were visualized using UMAP plot ( Figure 1D ). Moreover, the cell proportion regarding the sample phenotype was also provided ( Supplementary 1C, D ). The violin plot ( Figure 1E ), feature plot ( Figure 1F ) and heatmap ( Figure 1G ) demonstrated that the expression of canonical marker genes was distinctly expressed across the 9 major cell types. The presence and spatial distribution of epithelial cells, endothelial cells, fibroblasts cells and T cells were also characterized in human colorectal cancers using immunofluorescent staining ( Figure 1H ). Finally, the proportions of cells across each sample were presented. The results showed that the proportions of Epithelial cells and Myeloid cells were increased in tumor tissue, while the proportions of fibroblasts, plasma cells and NK-like cells were decreased in tumor ( Figure 1I ).

Copy number variation (CNV) is a hallmark of malignant cells. Here, we performed CNV analysis to distinguish the malignant cells from the normal epithelial cells, and to evaluate the clonality of colorectal cancers. First, we constructed the CNV atlas of epithelial cells derived from tumor adjacent and tumor tissue with normal mucosa epithelial cells as a reference ( Figure 2A ; Supplementary Figures 2A, B ). The CNV scores of each epithelial cells (n= 27,505) were calculated, and according to unsupervised cluster analysis, the epithelial cells from tumor and tumor adjacent tissue with minimal CNV and clustered with normal epithelial cells were identified as normal mucosal epithelial cells (n=13,106) ( Figure 2B ; Supplementary Figure 2C ), while the remainders were identified as malignant epithelial cells (n=14,399). Then, the CNV atlas of all malignant epithelial cells were reconstructed with all normal mucosa epithelial cells as a reference ( Figure 2C ). The malignant cells showed significant copy number amplification on chromosomes 7, 8, 13, 20 and local copy number loss on regions on chromosomes 3, 5, 6, 14 ( Figure 2C ). Unsupervised cluster analysis categorized the malignant cells into 5 clusters ( Figure 2D ). On the t-SNE plot containing malignant cells, there was a substantial overlap between cell clusters and patient identity, indicating that the variation across patients was a main contribution factor for tumor heterogeneity ( Figure 2D ). We also observed tumor cells derived from KRAS/TP53 dual mutant cancers had the highest CNV score ( Figure 2E , left and middle), while tumor cells from left-sided colorectal cancers had a higher CNV score than right-sided tumors ( Figure 2E , right). Since metabolic reprogramming plays a key role in colorectal tumorigenesis, we further assessed the metabolic heterogeneity of tumor cells from a single-cell perspective. The metabolic activity of tumor cells was higher than that of normal mucosal epithelial cells, and the metabolic characteristics of the 5 tumor subclusters were significantly different ( Figure 2F , left). One of the mechanisms leading to this metabolic variability is the distinct expression of genes regulating metabolic pathways ( Figure 2F , right). The metabolic activity of tumor cells also showed high heterogeneity based on KRAS/TP53 mutation status ( Figure 2G ). Finally, based on the single-cell CNV data, we constructed the clonal trajectory of the colorectal cancers using the UPhyloplot2 algorithm ( Figure 2H ). The results showed that the most samples shared the similar trunk clones. This result is similar to the previous results of bulk sequencing that identified the CNV of colorectal cancer genome with a trunk clone.

T lymphocytes play an important role in cancer immunotherapy. Their variety in cell composition, transcriptomic patterns, and functional features has a profound impact on T cell-based immunotherapy. Despite the fact that the comprehensive T-cell landscape of colorectal cancer has been well characterized, there is a dearth of in-depth analysis demonstrating the impact of tumor phenotype, including tumor sideness and mutation status, on T cell activity. Here, we identified 43,740 lymphocytes from all 78 samples (normal mucosa = 8,459 cells; tumor adjacent = 4,507 cells, tumor = 30,774 cells). After assessment and correction of the batch-effect, we performed PCA analysis and dimension reduction, and categorized the lymphocytes into 16 subgroups, which were majorly recognized as functional CD4⁺ and CD8⁺ T cells based on previously reported T cell markers ( Figures 3A, B ; Supplementary Figures 3A, B ) (31). T cell subgroup characteristics based on tSNE clusters, cell origin, tumor location, and sample identity ( Figure 3C ), as well as the spearman correlation across subgroups ( Figure 3D ), were also presented. Interestingly, the gender group shows a bias phenomenon regarding CD8⁺ Runx3⁺, CD8⁺ CTSW⁺ cytotoxic, and CD8+ exhaustion. Consequently, we investigated the gender-specific gene expression feature of these T cells. We observed that the expression of NKG7, GZMH, GZMA and GZMB in tumor-infiltrating CD8+ CTSW+ cytotoxic cells were significantly increased in female patients, indicating that these T cells have a stronger antitumor cytotoxic effect ( Supplementary Figure 3C ). Notably, two groups of potentially exhausted CD8⁺ T cells (CD8⁺ Tex), LAG3⁺ CD8⁺ Tex and LAG3⁺ CTLA4⁺ CD8⁺ Tex, were identified. The expression of GZMA, GZMB, IFN-γ and CXCL13 were significantly reduced in the LAG3⁺ CTLA4⁺ CD8⁺ Tex ( Figure 3E ), suggesting the terminal state of T-cell exhaustion. Next, we assessed the ICL/ICRs (immune checkpoint ligand/receptors) landscape expressed on T cells. The expression level of clinically targetable immune checkpoint ligand or receptors were quantified in all the T cell subtypes ( Figure 3F ). Notably, tumor sideness had a minimal effect on the ICL/ICRs expression pattern ( Supplementary Figure 3D ), while the ICL/ICRs showed unique distribution patterns based on tumor KRAS/TP53 mutation status ( Figure 3G ). These results suggested that the tumor sideness is not an effective predictor for immune checkpoint blockade-based immunotherapies. For instance, the main contributor of CTLA-4 and TIGHT were Tregs across all groups, the CD8⁺ T subgroups showed decreased expression of PD-1 in KRAS/TP53 dual mutant tumors, and the main contributors of PD-1 in the mutant cancer appeared to be Follicular helper T cells (Tfh cells) ( Figure 3G ).

Next, gene set enrichment analysis (GSEA) of the DGEs was performed according to KRAS/TP53 mutation status, and significantly downregulated metabolic regulating pathways, including Glycolysis, Cholesterol homeostasis, and Myc targets, were observed in CD4⁺ and CD8⁺ T cells from KRAS/TP53 dual mutant cancers ( Supplementary Figures 4A–E ). In addition to the expression levels of cytokines and ICL/ICRs, the immune-metabolic condition of T cells is also tightly associated with antitumor immune function (52). Resting T cells usually have a higher level of oxidative phosphorylation, while T cell activation attenuated mitochondrial respiratory, oxidative phosphorylation and enhanced glucose and glutamine metabolism. We observed a global decrease in CD4⁺ T cell metabolism in KRAS/TP53 dual mutant colorectal cancers ( Figure 3H ). The oxidative phosphorylation and TCA cycle, the two most prominent energy sources for T lymphocytes, in CD4⁺ T cells were both impaired in KRAS/TP53 dual mutant colorectal cancers ( Figure 3I ), indicating the degraded function of these T cells. Interestingly, the CD4⁺ T cells in left-sided colorectal tumors had significantly higher levels of oxidative phosphorylation and the TCA cycle than those in right-sided tumors ( Figure 3J ), suggesting a higher T cell function in left-sided colorectal cancers. This may be one reason why left-sided colon cancer has a better prognosis than right-sided colon cancer in clinical practice. Similar results were also observed in CD8 ⁺ T cells ( Figures 3K–M ). Taken together, these results indicated that KRAS/TP53 mutation status had a significant impact on molecular and metabolic function of tumor-infiltrating T cells.

To better understand the effect of mutant KRAS/TP53 on T cell functions. We also analyzed the transcriptome trajectory using Monocle2, which reconstructs putative branching transcriptional trajectories to identify potential relationships across calculated states, to construct the putative branching trajectories based on pseudotime ( Figures 4A–J ). Pseudotime ordering of all CD4⁺ T cells yields a total of 8 cell states organized into 4 main branches ( Figure 4A ). We next analyzed the trajectories of CD4⁺ T cells using branched expression analysis modeling (BEAM) and hierarchical clustering to identify genes enriched across states ( Figures 4B, C ; Supplementary Table 2 ). The heatmap clearly showed the enrichment of 64 genes according to pseudotime-based cell states ( Figure 4C ). Additional analysis revealed that the T cell trajectories of KRAS/TP53 mutant tumors were markedly different from those of wild-type tumors, with more trajectory branches ( Figure 4D ). Notably, GO analysis revealed that the expression profiles of a subset of CD4⁺ cells in KRAS/TP53 mutant tumors were enriched for response to hypoxic, implying that they may be affected by the hypoxic microenvironment ( Figure 4E ). While functions of CD8⁺ T cells were impaired in KRAS/TP53 mutant cancers, their cell-state branching trajectory seems to be unaffected by tumor mutation status, suggesting a different phenotypic plasticity compared to CD4⁺ T cells ( Figures 4F-J ).

Myeloid cells are considered to be the immune cells with the most plasticity. A total of 14,680 myeloid cells were identified in our cohort. The myeloid-derived cells were categorized into 10 subgroups, including macrophages, typical tumor-associated macrophages (TAMs), dendritic cells (DCs) and neutrophils, according to canonical myeloid cell markers ( Figures 5A, B , Supplementary Figures 5A, C ) ( 32). The existence of SPP1⁺ TAMs and FLOR2⁺ TAMs in human colorectal cancers were also validated using immunofluorescent staining ( Figure 5C ). Exception of neutrophils, plasma DCs (pDCs), and SPP1⁺ TAMs, the transcriptome profile of the myeloid cells showed a similarly to some extent ( Figure 5D ). According to the clustering results, C1q⁺ TAMs, SPP1+ TAMs, and CCL20⁺IL1B⁺ macrophages were substantially more abundant in tumor tissue than in non-malignant colorectal tissues ( Figure 5E ). The metabolic features of these 10 subgroups of myeloid cells also differed significantly ( Figure 5F ). The high metabolic activity of Ki67⁺ C1q⁺ TAMs may be related to the characteristics of active proliferation, while the elevated glycolysis activity of SPP1⁺ TAMs is also consistent with the hypoxic tumor microenvironment ( Figure 5G ). The expression of major ICRs (including PD-1, LAG-3, TIGHT and TIM-3) and ICLs (including PD-L1, PD-L2, HAVCR2, LGALS9, CEACAM1, FGL1, NECTIN2, and PVR) were also quantified in the myeloid cell types ( Figure 5H ). The C1q⁺ TAMs were the predominant contributor of immunosuppressive ICLs, including PD-L1, PD-L2, HAVCR2, LGALS9, and CEACAM1. Similarly, to the results found in T cells, the ICL/ICRs expressed on myeloid cells also showed unique distribution patterns based on tumor KRAS/TP53 mutation status ( Figure 5I ). The GSEA analysis of DEGs showed a significant upregulated TNF-α pathway and downregulated inflammatory response pathway in myeloid cells in KRAS/TP53 dual mutant colorectal cancers ( Figure 5J ). Despite the fact that the energy metabolism function of myeloid cells was significantly hindered in KRAS/TP53 mutant cancers, the glycolysis function of TAMs did not appear to be substantially reduced ( Figure 5K ), suggesting adaptation of TAMs to the KRAS/TP53 mutant tumor microenvironment. Finally, trajectory analysis using monocle3 also revealed a more complex evolutionary trajectory of myeloid cells KRAS/TP53 in tumors ( Figure 5L ). In summary, our data indicated that the KRAS/TP53 status significantly affects the cellular states of tumor-infiltrating myeloid cells, particularly macrophages. Since macrophages are important providers of immunosuppressive ligands, patient stratification based on KRAS/TP53 status for individualized treatment may be a viable immunotherapeutic strategy for colorectal cancers.

Recent studies have underlined the significance of cell-to-cell communication in the development of diverse malignancies. The CellphoneDB algorithm and R ‘iTALK’ (version 0.1.0) package were used to evaluate the crosstalk across the immunological and non-immunological cellular components in normal mucosa, tumor adjacent and tumor tissue ( Figure 6A ). Compared with normal mucosal epithelial cells, the communication between tumor cells and immune cells was significantly enhanced ( Figure 6B , left). As expected, myeloid cells and T cells have also showed a strong ability to communicate with the tumor microenvironment ( Figure 6B , right). Further results revealed complex cell-cell communication networks across cancer and immune cells ( Figure 6C ). Next, we investigated the role of cell-cell communication in TME remodeling in colorectal tissue from the perspectives of immune checkpoints, cytokines and growth factors ( Figure 6D ). We found that the CD274-PDCD1 ligand-receptor pair’s connection was quite weak through the integration analysis. The ligand-receptor communication density identified by the iTALK algorithm is majorly determined by the expression levels of receptor and ligand genes in target cells. In contrast, some inhibitory ICR/ICL pairs such as CD24-SIGLEC10 and LGALS9-HAVCR2 were frequently observed between myeloid cells, master cells, B cells and tumor cells ( Figures 6C, D ). Of note, the immune and non-immune components showed distinct cell-cell communication patterns in normal mucosa, tumor adjacent and colorectal tumors. Tumor tissue presented more unique inhibitory ICR/ICL pairs such as CD80-CTLA4 and CD70-CD27 ( Figure 6D ). Endothelial cells are the primary involvers of cell crosstalk via growth factor in tumor adjacent region, while malignant epithelial cells are dominantly involved in the tumor tissues. In tumor tissue, the malignant epithelial cells not only produced growth factors such as HB-EGF, but also expressed CD44 and CD9 receptors to receive growth signaling, further promoting the tumor progression ( Figure 6D ).

Next, we detailed the cell-to-cell interactions across malignant epithelial cells, myeloid cells and T cells in colorectal tumors tissue ( Figures 6E, F ). The data demonstrated that Tregs were the primary source of CTLA4, while Tfh produced the CD40L signal, suggesting that CD4⁺ T cells play a complicated role in tumorigenesis ( Figure 6F ). For tumor cells, CD24 was the most important intrinsic molecular providing a “don’t eat me” signal, meanwhile, SDC4 was the most important receptor for cytokine signaling ( Figure 6F ). Finally, our data showed that tumor sideness had minimal effect on cell-to-cell communication, however KRAS/TP53 mutation status had a significant impact ( Figures 6G, H ). For instance, compared with KRAS/TP53 wildtype colorectal cancer, the CD24-SIGLEC10 immunosuppressive signaling and the HBEGF-CD9 tumor growth signaling across tumor and immune cells were dramatically increased in KRAS/TP53 dual mutant tumors ( Supplementary Table 3 ), which was highly consistent with the high malignancy and poor prognosis of KRAS cancers. Taken together, we demonstrated the cell interactions across immune cells and non-immune cells in colorectal tumors based on tumor sideness and mutation status, providing evidence for the combination therapy targeting immune checkpoints, growth factors and cytokine signaling, and providing additional support for the individualized treatment of colorectal cancer based on KRAS/TP53 mutation status.

Finally, we validated our findings using bulk sequencing data from TCGA-COAD cohort (n= 478) using ssGSEA analysis. We quantified the immune cell components in the bulk sequencing data using both canonical ssGSEA markers and scRNA-seq identified marker genes. The heatmap of the immune landscape and corresponding clinical information was presented ( Figure 7A ). The samples were categorized into two categories, the high immune infiltration and low immune infiltration group, using unsupervised clustering ( Figure 7A ). There was no correlation between T-cell and epithelial-cell infiltration score in the patient cohort, suggesting that the ssGSEA makers could efficiently distinguish the major cellular components in colorectal cancers ( Figure 7B , left). The infiltration score of the same cell type was evaluated using single-cell markers and canonical gene sets, and the results showed a high degree of congruence ( Figures 7B, C ). We further investigated the correlation between the immune infiltration score of these 59 immune cells and the prognosis of colorectal cancers ( Figure 7D ). The results showed that increased epithelial score and CD4⁺ unknown T cell score were associated with a poorer prognosis, while increased B cell score and CD8⁺ CTSW⁺ Cytotoxic T cell score were associated with a better prognosis in the TCGA-COAD cohort ( Figure 7D ). After stratifying patients based on tumor mutations, we observed that fibroblast enrichment predicted a poor prognosis, whereas CCL20⁺ IL1B⁺ macrophage enrichment indicated significantly prolonged survival in KRAS/TP53 wild-type tumors ( Figure 7E ; Supplementary Figures 6A–D ). Meanwhile, CD4⁺ unknown T cell and B cell infiltration score were optimal prognosis indicators in KRAS/TP53 dual mutant tumors ( Figure 7F ). Collectively, our results highlight the importance of diverse immune cells in the prognosis of colorectal cancer; we also discuss the variation in immune cell score in prognosis based on mutation status. The results provided evidence that stratification of the patients might improve the prognosis prediction and efficiency of immunotherapy.