Ten data sets were downloaded from the public database, including clinical data and transcriptome data of 2,078 patients with CRC. GSE39582, GSE14333, GSE17536, GSE17537, and GSE72968 were all microarray data of the GPL570 platform and integrated as a training set (the combined GEO cohort). The CEL format data of the microarray was downloaded using the “GEOquery” package (Davis and Meltzer, 2007). ReadAffy function in the “affy” package was used to read data in CEL format (Gautier et al., 2004), background correction and standardization were carried out with RMA, and then the “SVA” package was used to remove batch effect among the data sets (Chakraborty et al., 2012). Probes corresponding to multiple genes were deleted, and the average expression level was taken when multiple probes corresponded to one gene. Clinical data and FPKM value (fragments per kilobase million) of TCGA-COAD were obtained from the UCSC website as an external validation set (Navarro Gonzalez et al., 2021). Then, FPKM was converted to TPM (transcripts per kilobase million) for subsequent analysis. GSE39582, GSE17536, and GSE17537 were used as the internal validation sets. The paired samples of GSE44076, GSE32323, GSE89076, and GSE113513 datasets have been reserved for verifying the gene expression level in the TMEIG score system. Detailed information on the data sets is shown in Supplementary Table S1.

CIBERSORT (Newman et al., 2015) and ESTIMATE (Becht et al., 2016) algorithms can infer the composition of 22 types of immune cells, immune score, stroma score, and tumor purity in the TME based on transcriptome data. In this study, transcriptome data from the combined GEO dataset (1,175 samples) and TCGA COAD dataset (471 tumors) were used for CIBERSORT and ESTIMATE analyses. After filtering out low expression immune cells, 15 types of immune cells were retained. ssGSEA is another algorithm for estimating immune cell composition in solid tumor TME and can also be used to calculate adaptive and innate immune components of samples. In the present study, the adaptive and innate immune scores of each sample were obtained with the “GSVA” package. ssGSEA parameters were set as follows: method = “ssgsea,” KCDF = “Gaussian.”

“ConsensusClusterPlus” is a re-sampling unsupervised clustering method to verify the rationality of clustering (Wilkerson and Hayes, 2010). In the combined GEO cohort, TME subtypes and TMEIG subtypes were obtained using “ConsensusClusterPlus” package. The parameters were set as follows: MaxK = 9, REPS = 1,000, pItem = 0.8, pFeature = 1, clusterAlg = “PAM,” distance = “Euclidean,” and seed = 1.

WGCNA analysis was used to identify gene modules most associated with traits (Langfelder and Horvath, 2008). Stromal score, immune score, estimate score, tumor purity, adaptive immune, innate immune, TME cluster 1, and TME cluster 2 were inputted as traits. The key parameters of WGCNA were set as follows: soft threshold power β = 4, scale-free R2 = 0.89. The relationship between modules and traits was analyzed using the Pearson correlation method. Gene significance (GS) and module membership (MM) are two important indicators in WGCNA analysis. GS is the correlation between the gene and trait. MM is defined as the correlation of the module eigengene and the gene expression profile. Genes with GS > 0.2 and MM > 0.8 are usually considered hub genes.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa and Goto, 2000) analyses were employed to explore the biological functions of the modules in WGCNA using the R package “clusterprofiler” (Yu et al., 2012). An adjusted p-value of less than 0.05 was regarded as statistically significant. In addition, Gene Set Enrichment Analysis (GSEA) was conducted (Subramanian et al., 2005). The gene sets “c2.cp.kegg.v6.2.symbols.gmt,” “c5.all.v7.0.symbols.gmt,” and “h.all.v7.2.symbols.gmt” on MSigDB website were chosen as the reference gene sets (Liberzon et al., 2015). The log fold change (FC) of differentially expressed genes between two groups was used as the input list for GSEA analysis. When analyzing the biological functions related to one gene, the Pearson correlation coefficient was used as the input list.

First, univariate Cox proportional hazards regression was employed to identify the prognostic genes using the “survival” R package. Genes with a p-value less than 0.05 were regarded as the candidates, input to least absolute shrinkage and selection operator (LASSO) regression (Friedman et al., 2010). After ten cross-validations, five prognostic genes and the corresponding coefficient remained when lambda = 0.0713387182. Then, TMEIG score was established as follows: TMEIG score = ∑ i Coefficient of  ( i ) × Expression of gene  ( i )

The regression coefficient of the gene was designated (i) in the LASSO regression model. The combined GEO cohort was used as a training set, whereas GSE39582, GSE17536, and GSE17537 were used as the internal validation sets. In addition, the TCGA COAD cohort served as the external validation set.

Only GSE39582, GSE17536, GSE17537, and GSE72968 had overall survival data in the combined GEO cohort (Supplementary Table S2). The survival time was converted to months, and samples with a survival time of less than 1 month were removed during survival analysis. Finally, 864 samples in the combined GEO cohort and 435 samples in the TCGA COAD cohort were used for survival analysis (Supplementary Table S2). According to the best cutoff value determined using the “survminer” package, the patients were divided into high and low expression groups. Log-rank test was employed to evaluate statistical significance. Kaplan–Meier (KM) plots were visualized using the “survminer” R package. The risk factors diagrams were visualized using the “ggrisk” R package.

The mutation data of TCGA COAD were downloaded from the TCGA website and analyzed using the “maftools” package (Mayakonda et al., 2018). The tumor mutation burden (TMB) was calculated using the following formula: (total mutation/total covered bases) × 106. The driver genes in somatic alterations were also identified using the “maftool” package.

Tumor Immune Dysfunction and Exclusion (TIDE) algorithm was employed to predict ICB response based on the gene expression related to T cell dysfunction (Dysfunction) and T cell exclusion (Exclusion). The lower the TIDE score is reportedly associated with a better immunotherapy response (Jiang et al., 2018). Furthermore, the scores of cancer-associated fibroblasts (CAF), Dysfunction, Exclusion, M2 macrophages (M2), myeloid-derived suppressor cells (MDSC), and TIDE were obtained from the TIDE website. The IMvigor210 cohort is a large cohort of patients with metastatic urothelial cancer under anti–PD-L1 therapy (atezolizumab), which can be downloaded from the Creative Commons 3.0 license. GSE78220 is an anti–PD-1 therapy cohort containing mRNA expression data from pre-treatment melanomas. The two cohorts were used to validate the predictive power of the TMEIG score for ICB response.

The human CRC cell lines SW620, RKO, HCT116, HT29, and NCM460 (ATCC) were cultured in RPMI-1640 medium (Gibco, United States ) supplemented with 10% fetal bovine serum (FBS, Biological Industries, United States ) at 37°C in a humidified 5% CO₂ atmosphere.

Total RNAs of cell lines were extracted by TRIzol reagent (Invitrogen, United States ) and then was reversely transcribed as cDNA via PrimeScript™ RT Master Mix (Takara, Japan). Quantitive real-time PCR was performed using PowerUp™ SYBR™ Green Master Mix (Applied Biosystems, United States) in the StepOne™ Real-Time PCR System (Applied Biosystems). Each reaction was tested in triplicate. ACTB was used as the internal reference, and the 2(^−ΔΔCT) method was used for calculating the relative mRNA expression. The following primer sets were used:

Human FABP4: Forward: 5ʹ-GGG​CCA​GGA​ATT​TGA​CGA​AG-3ʹ, Reverse: 5ʹ-TCG​TGG​AAG​TGA​CGC​CTT​TC-3ʹ; Human SCG2: Forward: 5ʹ-GTG​AAG​CGA​GTT​CCT​GGT​CA-3ʹ, Reverse: 5ʹ-ATG​CTC​TTT​GAT​GGC​CTG​CT-3ʹ; Human CALB2: Forward: 5ʹ-GAA​GGC​AAG​GAA​AGG​CTC​TGG-3ʹ, Reverse: 5ʹ-GCC​ATC​TCG​ATT​TTC​CCA​TCT​G-3ʹ; Human SERPINE1: Forward: 5ʹ-CCT​GGT​TCT​GCC​CAA​GTT​CT-3ʹ, Reverse: 5ʹ-CCA​TGC​GGG​CTG​AGA​CTA​TG-3ʹ; Human HOXC6: Forward: 5ʹ-CAC​TAA​CCC​TTC​CTT​ATC​CTG​CC-3ʹ, Reverse: 5ʹ-TCA​TAG​GCG​GTG​GAA​TTG​AGG-3ʹ; Human ACTB: Forward: 5ʹ-GAT​TCC​TAT​GTG​GGC​GAC​GA-3ʹ, Reverse: 5ʹ-AGG​TCT​CAA​ACA​TGA​TCT​GGG​T-3ʹ.

For the IHC experiment, we collected 16 pairs of CRC tissue (cancer and adjacent normal tissue) from patients who received surgery at the Department of Colorectal Cancer Surgery, the Second Affiliated Hospital of Harbin Medical University (Harbin, China) between January 2014 and December 2020. Ethics approval was also granted by the Ethics Committee of Harbin Medical University (No. KY 2022-063). The primary antibodies used in IHC were as follows: anti-FABP4 (Proteintech, #12802-1-AP, 1:200 dilution), anti-SCG2 (Proteintech, #20357-1-AP, 1:200 dilution), anti-CALB2 (Proteintech, #12278-1-AP, 1:200 dilution), anti-HOXC6 (Affinity, #DF3078, 1:150 dilution), and anti-PAI1 (SERPINE1) (Affinity, #AF5176, 1:200 dilution). Paraffin sections were incubated with primary antibodies at 4°C overnight, followed by treatment with HRP-conjugated secondary antibodies at 37°C temperature for 60 min following PBS rinse. Then, tissues were counter-stained with hematoxylin and further treated with DAB for 2 min. The IHC results were independently analyzed by two experienced pathologists. A staining scoring system was evaluated by both staining intensity (negative = 0, weak = 1, and strong = 2) and staining area (<5% = 0, 5%–25% = 1, 25%–50% = 2, 50%–75% = 3, and >75% = 4). The staining intensity score was computed, and the score of the staining area was the final staining score. A total score of ≤3 was considered a weak expression. A total score of >3 was considered a strong expression. The details of IHC performance and scoring system are described in Supplementary Table S3.

Heat maps were visualized with the “ComplexHeatmap” package (Gu et al., 2016). The “ggplot2” package was used to visualize boxplots, scatter plots, and Sankey plots. The log-rank test and Pearson method were used for KM survival and correlation analyses, respectively. The difference between the two groups was tested by the Wilcox test. It should be noted that * represented a p-value less than 0.05, ** represented a p-value less than 0.01, *** represented a p-value less than 0.001, and **** represented a p-value less than 0.0001. All analyses were performed in R 4.0.3.

The flow diagram describes the construction of TME subtypes and the TMEIG score in CRC (Figure 1). We integrated microarray data of 1,175 patients with CRC from the GPL570 platform and then used the combat function of the “SVA” package to remove batch effects. The principal component analysis (PCA) diagrams of five cohorts before and after batch effect removal are shown in Figures 2A,B. The results indicated that the batch effect was negated, and the combined cohort could be used for subsequent analysis. To fully dissect the heterogeneity of the TME in patients with CRC, the CIBERSORT algorithm was used to assess the proportion of immune cells in the TME. Macrophages and mast cells were the most abundant immune populations in the combined GEO cohort, followed by memory resting CD4 and CD8 T cells. Figures 2C,D shows the proportion of immune cells in each patient, which partly reflects the heterogeneity of immune cells in the TME. A total of 15 types of immune cells were retained after eliminating low expression cells (such as memory B cells, CD4 naive T cells, gamma delta T cells, activated NK cells, monocytes, resting mast cells, and eosinophils). The detailed results of the CIBERSORT analysis are shown in Supplementary Table S4. Then, the ESTIMATE algorithm was used to calculate patients’ immune scores and stromal scores. Collectively, CIBERSORT and ESTIMATE algorithms were used to comprehensively describe the correlations among the immune cells, immune score, and stromal score in the tumor microenvironment of patients with CRC (Figure 2E). Resting NK cells were inversely correlated with M0/M1/M2 macrophages (correlation values = −0.09, correlation values = −0.28, correlation values = −0.26; p-value < 0.05; Supplementary Table S4). Furthermore, CD8 T cells were negatively related to M0 macrophages (correlation values = −0.31, p-value < 0.05) and positively related to M1 macrophages (correlation values = 0.16, p-value < 0.05).

Based on these quantitative indicators describing the TME, we conducted unsupervised clustering in these 1,175 patients using the “ConsesusClusterPlus” package. As shown in Supplementary Figure S1, the clustering result was the most stable when K = 2. The PCA plot also demonstrated significant differences between the two clusters (Figure 3A). Then, survival analysis was employed to compare the prognosis between the two TME clusters (Supplementary Table S5). The OS in TME cluster 2 was significantly better than that in TME cluster 1 (Figure 3B log-rank test, p = 0.047). Furthermore, we explored 11 critical biological gene signatures between the two TME subtypes using a heat map (Mariathasan et al., 2018). The results indicated that cell cycle genes and DNA damage repair (DDR) genes were markedly decreased, and angiogenesis (Angio) markers, TGFβ receptor and ligand (TGFβ), antigen-processing machinery (APM), and F-TBRS genes were significantly increased in TME cluster 1 as compared to TME cluster 2 (Figure 3C). In addition, CD8 Teff cells and immune checkpoint signatures (ICI) were highly expressed in TME cluster 1 (Figures 3C–E). The low expression of cell cycle-associated genes may indicate that tumor cells in TME cluster 1 were in a dormant phase and were not easily cleared by the immune system. A comparison of the immune score and stromal score revealed that the immune score and stromal score of TME cluster 1 were higher than those of TME cluster 2 (Figure 3D). CIBERSORT analysis demonstrated that immunosuppressive cells (M0, M1, and M2) were significantly reduced, and immunoreactive cells (CD8 T cell, CD4 memory resting T cells, resting dendritic cells, and activated dendritic cells) were significantly increased in the TME cluster 2 (Figure 3F). Furthermore, the GSEA results indicated that immune-related functions (activation of immune response, positive regulation of cytokine production, cytokine–cytokine receptor interaction, and IL6-JAK-STAT3 signaling) significantly varied between TME clusters 1 and 2 (Figures 3G–I; Supplementary Table S6).

To identify key gene modules in TME clusters 1 and 2, WGCNA analysis was employed. The WGCNA analysis processes are shown in Supplementary Figures S2A–E. Adaptive immunity and innate immunity were derived from ssGSEA analysis. Stromal score, immune score, estimate score, tumor purity, adaptive immunity, innate immunity, and TME clusters 1 and 2 were used as traits. The heat map of module–trait relationships is shown in Figure 4A. Results indicated that blue (cor = 0.79, p < 1e−200), brown (cor = 0.53, p = 2.2e−63), and green (cor = 0.98, p < 1e−200) modules displayed the high correlations with adaptive immunity (Figures 4B–D).

Thus, the blue, brown, and green modules were identified as the key modules. We performed GO and KEGG analyses to explore the biological functions of genes within the key modules. As shown in Supplementary Figures S2F–H, GO and KEGG terms were closely related to the immune function, such as regulation of immune system process, cytokine production, TNF signaling pathway, TNF superfamily cytokine production, and TNF superfamily cytokine production. In the three modules, 223 genes with GS > 0.2 and MM > 0.8 were identified as candidate genes. We used the “Limma” package to obtain the differentially expressed genes (DEGs) between TME clusters 1 and 2, and the results are shown in the volcano map (Figure 4E). p-value < 0.05 and logFC >0.5 were set as parameters, and 719 DEGs were obtained (698 upregulated and 21 downregulated). Since there were very few downregulated genes, we mainly used the upregulated genes to compare with the candidate genes of WGCNA. A total of 202 TMEIGs were eventually identified after comparing candidate genes with upregulated genes (Figure 4F).

A total of 202 TMEIGs were again used for unsupervised clustering in the combined GEO cohort, and the clustering process is shown in Supplementary Figures S3A–E. The clustering result was most stable when k = 2. PCA plot also revealed significant differences between the TMEIG subtypes (Figure 4G; Supplementary Table S7). The KM plots revealed that patients in TMEIG cluster A exhibited better OS (Figure 4H, log-rank test, p = 0.047). Similarly, the heat map of tumor-related pathways showed that cell cycle and DDR signatures were significantly decreased in TMEIG cluster B. Angio, transforming growth factor-beta (TGFβ), antigen processing machinery (APM), TGF-beta response signatures (TBRS) of fibroblasts (F-TBRS), and immune checkpoint signatures were increased in TMEIG cluster B as compared to that in cluster A (Supplementary Figure S3J). ESTIMATE analysis showed that TMEIG cluster B had higher immune and stromal scores than TMEIG cluster A (Figure 4I). CIBERSORT analysis showed that immunosuppressive cells (M0/M1/M2 macrophages) increased significantly in TMEIG cluster B, whereas immunoreactive cells (CD8 T cell, CD4 memory resting T cells, resting dendritic cells, and activated dendritic cells) decreased significantly compared to TMEIG cluster A (Figure 4K). GSEA results indicated that immune-related functions (activation of immune response, cytokine-cytokine receptor interaction, and inflammatory response) significantly varied between TMEIG cluster A and TMEIG cluster B (Figures 4L–N). These results demonstrated that TMEIG subtype clustering accurately reflected the differences between TME subtypes.

Gene signature is a simple and effective model widely used in clinical practice (Paik et al., 2004; van 't Veer et al., 2002; Parker et al., 2009). To further facilitate the application of TME subtypes in CRC, we intended to construct a TMEIG score system. First, 662 genes were obtained by comparing 2,799 genes in blue, brown, and green modules with 698 DEGs (Supplementary Figure S4A). Univariate Cox regression analysis was performed in the combined GEO cohort and TCGA COAD cohort. With a p-value less than 0.05, 287 and 47 prognostic genes were obtained from the combined GEO and TCGA COAD cohorts, respectively. Then, Lasso regression was used to identify 27 common genes in the combined GEO cohort (Supplementary Figure S4B). Details of the Lasso regression are shown in Supplementary Figures S4C and D. After cross-validating the results ten times, five genes (SERPINE1, FABP4, SCG2, CALB2, and HOXC6) and their corresponding lambda coefficients were obtained when lambda = 0.0713387182. The TMEIG score was constructed based on the expression and coefficient of the five genes as described in the methods. According to the optimal cutoff value, patients were divided based on whether their TMEIG scores were high or low (Supplementary Table S8). OS analysis suggested that the high TMEIG scores were associated with poor prognosis in patients with CRC (Figure 5A, log-rank test, p < 0.0001).

Internal validation cohorts indicated that the OS of the patients with high TMEIG scores was poorer than those with low scores (Figures 5B–D, GSE39582, GSE17536, GSE17537, log-rank test, p < 0.0001). In addition, the PFS and DSS in the low TMEIG score group were superior to those of the high TMEIG score group (Figures 5E–G, GSE17536 DFS, GSE17536 DSS, GSE17537 DFS, log-rank test, p < 0.0001). Then, the TCGA COAD cohort also revealed that the overall survival of the high TMEIG score group was poorer (Figure 5H; Supplementary Table S8). When analyzing the relationship between clinicopathological parameters and the TMEIG score, we observed that the scores of patients exhibiting stage III & IV, T 3 & 4, lymphatic invasion, and venous invasion were significantly higher (Figures 5I–N; Supplementary Table S9), suggesting the high TMEIG score was associated with poor clinical prognosis. Furthermore, the risk factor diagrams of the combined GEO and TCGA COAD cohorts indicated that the high TMEIG score group had significantly higher mortality than the low TMEIG score group (Figure 5O; Supplementary Figure S4E). All the results demonstrated poor prognoses in patients with high scores. The distribution of patients in the TME clusters, TMEIG clusters, and TMEIG score groups are shown in Figure 5P, which indicates that most patients in TME cluster 1 belonged to TMEIG cluster B and the high TMEIG score group. Consistent with the results, the high TMEIG score group had a higher proportion of patients in TME cluster 1 and TMEIG cluster B (Figures 5Q,R). Moreover, patients in TME cluster 1 and TMEIG cluster B exhibited higher TMEIG scores than that in TME cluster 2 and TMEIG cluster A (Supplementary Figures S4F,G). This evidence demonstrated that the TMEIG score could effectively surrogate TME and TMEIG subtypes. Furthermore, TMEIG score, age, sex, and T, N, M, and TNM stages were included in stepwise Cox regression in the GSE39582 cohort, which possessed comprehensive clinical information. Results indicated that TMEIG score combined with age, sex, M stage, and TNM stage exhibited the best predictive power (Supplementary Figure 5S, AUC: 0.75, 0.69, 0.7, 0.72, and 0.73 at 1, 3, 5, 7, and 10 years, respectively), and was validated in the TCGA COAD cohort (Supplementary Figure S4H).

To evaluate whether the TMEIG score could predict the efficacy of ICB treatment, we analyzed the expression of crucial immune checkpoint molecules between high and low TMEIG score groups. The results showed that the expression of immune checkpoint molecules (PD-1 [PDCD1], PD-L1 [CD274], cytotoxic T-lymphocyte associated protein 4 [CTLA4], B- and T-lymphocyte attenuator [BTLA], T cell immunoglobulin and ITIM domain [TIGIT], hepatitis A virus cellular receptor 2 [HAVCR2], and lymphocyte-activation gene 3 [LAG3]) was significantly higher in the high score group (Figure 6A). Patients with CRC who have microsatellite instability-high (MSI-H) tumors are more likely to benefit from immune checkpoint inhibitors than patients with microsatellite stable (MSS)/MSI-low (MSI-L). To explore the relationship between TMEIG score and MSI status, the MSI status of TCGA COAD patients was downloaded from the supplements of previous studies focusing on MSI detection (Supplementary Table S10). There were 72 patients identified as MSI-H and 355 identified as MSI-L/MSS in TCGA COAD determined by MSI-PCR. We then analyzed whether the TMEIG score had prognostic value across MSI-H and MSI-L/MSS subgroups. KM plots demonstrated that patients with high TMEIG scores exhibited poor overall survival in MSI-H and MSI-L/MSS subgroups (Supplementary Figures S5A,B). Further analysis showed that patients with MSI-H possessed higher TMEIG scores (Figure 6B, p = 6.9e−06). Next, we explored the proportion of patients with MSI-H and MSI-L/MSS in high and low TMEIG score groups. We observed more MSI-H CRC patients in the high TMEIG score group (Supplementary Figure S5C, High: 26%, Low: 14%). Previous studies reported that TMB was a predictor of the efficacy of ICB therapy. When exploring the TMB of patients from the TCGA COAD cohort, results indicated no statistical difference between high and low TMEIG score groups (Figure 6C). In addition, the top six driver genes with the highest alteration frequency were further analyzed. The alteration frequency of APC, TP53, TTN, KRAS, PIK3CA, and MUC16 in high and low TMEIG score groups are displayed in Supplementary Figure S5D. The ESTIMATE algorithm revealed that immune, stromal, and estimate scores were significantly higher in the high TMEIG score group (Figure 6D). Furthermore, pathway heat maps demonstrated that EMT, angiogenesis (vimentin [VIM], Twist-related protein 1 [TWIST1], zinc finger E-box binding homeobox 1 and 2 [ZEB1 and ZEB2]), and T-cell factor-beta (TCFβ) signatures were significantly more activated in the high TMEIG score group. In contrast, cell cycle and DDR signatures were highly expressed in the low TMEIG score group (Figure 6E). These results indicated significant differences in the TME and biological pathways between high and low TMEIG score groups. To dissect the relationship between TMEIG score and ICB response, we used the TIDE algorithm to predict ICB response based on transcriptome signatures. The TIDE algorithm (Figure 6F) showed that the TMEIG score was positively correlated with CAF (R = 0.68, p < 2.2e−16), Dysfunction (R = 0.37, p < 1.7e−15), and Exclusion (R = 0.14, p < 0.0036), and negatively correlated with M2 macrophages (R = 0.51, p < 2.2e−16) and MDSC (R = −0.3, p < 2.7e−10). In addition, there was a strong positive correlation between the TMEIG and TIDE scores (Figure 6F, R = −0.15, p < 0.003).

Then, we explored the predictive power of the TMEIG score in two ICB therapy cohorts. In IMvigor210 and GSE78220 cohorts, the ICB response rates were significantly lower in the high TMEIG score group (Figures 6G,H). Notably, patients with high TMEIG scores exhibited worse overall survival (Figures 6I,J; IMvigor210: log-rank test, p = 0.018, GSE78220: log-rank test, p = 0.053). The insignificant result in GSE78220 (28 patients) can be attributed to the small sample size. Furthermore, the biomarker evaluation module on the TIDE website was used to assess the accuracy of the TMEIG score using multiple ICB cohorts as compared to other published biomarkers. The TMEIG score demonstrated an AUC of more than 0.5 in nine out of 16 ICB cohorts (Supplementary Figure S6), demonstrating its robustness as a predictive biomarker (Fu et al., 2020).

To further understand the functions of the biomarker genes consisting of the TMEIG score, we analyzed the expression levels of SERPINE1, FABP4, SCG2, CALB2, and HOXC6 in the TCGA-COAD cohort. The results demonstrated that the expression values of SERPINE1 and HOXC6 were significantly upregulated in tumors, whereas FABP4, SCG2, and CALB2 were highly expressed in normal tissues (Figure 7A). Furthermore, the same results were found in paired differential expression in multiple CRC cohorts (Figures 7B–E, GSE32323, GSE44076, GSE89076, and GSE113513). Our immunohistochemistry (IHC) results revealed that CALB2 exhibited relevant stronger staining in two cases of tumor (2/16), with the other 14 cases displaying low expression (14/16). For FABP4, seven and two cases in tumor and normal samples, respectively, exhibited stronger staining, whereas nine and 14 cases in tumor and normal samples, respectively, were with low staining intensity. For HOXC6, seven and nine samples out of 16 exhibited strong staining in normal and tumor samples, respectively. The protein level of SCG2 was high in 12 cases of CRC samples (12/16) and 11 cases of normal samples (11/16). SEP1NG1 was found strongly stained in 11 cases of tumor samples (11/16) and nine cases of normal samples (9/16). Representative immunohistochemical images are shown in Figure 7F, and the high-resolution images are shown in Supplementary Figure S7. The qPCR experiments (Supplementary Figure S8A) revealed that the expression of HOXC6 and SERPINE1 was significantly upregulated in RKO and HT29 cell lines, and FABP4 expression was downregulated in nearly all CRC cell lines analyzed. Although SCG2 and CALB2 were downregulated in patients with CRC in multiple cohorts, qPCR experiments showed that their expression was upregulated in several CRC cell lines (such as HCT116 and HT29). This discrepancy may be due to false positives in RNA sequencing or the heterogeneity between clinical tissues and tumor cells. Studies involving more clinical samples or cell lines may be needed to confirm the expression of the two genes at the RNA and protein levels in the future. KM plots of these genes are displayed in Supplementary Figure S9A. Results indicated that all genes were closely related to OS. GSEA analysis indicated that the five genes were involved in multiple cancer biological functions: cell motility, angiogenesis, cell migration, programmed cell death, MAPK signaling pathway, and PI3K-Akt signaling pathway. Notably, the immune-related pathway “cytokine-cytokine receptor interaction” was also significantly enriched in most of these genes (SERPINE1, FABP4, SCG2, and HOXC6) (Supplementary Figure S9B). Next, we summarized several immunological molecules from our previous studies, such as immune checkpoint genes and cytotoxic genes, and analyzed the correlation with five genes in TCGA COAD (Supplementary Figures S8B,C). Results demonstrated that most of the five genes were significantly correlated with immune checkpoint genes (BTLA, CD274, CTLA4, HAVCR2, LAG3, PDCD1, and TIGIT) and cytotoxic genes (granzyme A [GZMA], GZMB, GZMK, GZMM, interferon-gamma [IFNG], perforin 1 [PRF1], and tumor necrosis factor superfamily member 11 [TNFSF11]), which revealed that these genes might play an important role in tumor immunity. We then explored the correlation between the five genes and immune cells infiltrating the TME. As shown in the correlation heatmap, the five genes were positively related to macrophages (such as M0, M1, and M2), inversely correlated with resting NK cells and resting memory CD4⁺ T cells, CD8⁺ T cells, and activated memory CD4⁺ T cells (Figure 7G), which might explain the poor ICB response in the high TMEIG score group. Moreover, the gene set prioritization module on the TIDE website indicated that HOXC6 was the most appropriate target to treat TME resistance to ICB (Figure 7H). The expression of HOXC6 was positively associated with T cell dysfunction in GSE12417, METABRIC, and TCGA Endometria datasets (Figure 7H, left panel). In addition, high HOXC6 expression was also associated with poorer ICB outcomes in multiple cohorts treated with ICB (Figure 7H, second to left panel). Among the immune-suppressive cell types, HOXC6 was highly expressed on the MDSC and CAF (Figure 7H, right panel).