Prognostic stratification of colorectal cancer by immune profiling reveals SPP1 as a key indicator for tumor immune status

Introduction

Representing one of the major contributors to cancer-related mortality worldwide, colorectal cancer (CRC) accounts for over 1.9 million new cases and 0.9 million deaths annually [1]. Surgery is the only curative approach for CRC, and patients with advanced inoperable disease often receive a multidisciplinary, combined treatment model based on chemotherapy. A diverse array of therapeutic options can now be contemplated, inter alia, radiotherapy, chemotherapy, immunotherapy, and targeted therapy. Immunotherapy has significantly changed the therapeutic pattern for CRC by manipulating the patient’s intrinsic immune system to combat the tumor. Several inhibitors targeting different checkpoints, such as programmed cell death 1 (PD-1), programmed death-ligand 1 (PD-L1), and cytotoxic T lymphocyte antigen 4 (CTLA4), have been approved by the Food and Drug Administration (FDA) for immunotherapy of CRC [2]. In CRC, however, whether immune checkpoint blockades (ICBs) can be beneficial for CRC patients depends on the microsatellite instability (MSI) or DNA mismatch repair deficiency (dMMR) status of the tumor, while CRC with high MSI or dMMR accounts for only about 15 % of patients [3]. These tumors are capable of activating the immune system and responding to ICBs due to their high mutation burden and neoantigen expression [4].

Despite advancements in treatment for CRC, a high percentage of patients still do not experience tumor shrinkage or survival prolongation, with the average 5-year survivorship of advanced CRC far from satisfactory over the past several decades [5], 6]. Amidst this evolving landscape and expanded array of treatment options, a crucial, yet unanswered question pertains to the screening of patients who would benefit from immunotherapy. There are many novel strategies being explored to optimize immunotherapy for CRC, yet a comprehensive and in-depth understanding of the complex and dynamic immune regulatory mechanisms in the CRC tumor microenvironment (TME), as well as a precise characterization of the immune profile of CRC across different subtypes and stages, is still absent. There is a pressing need for effective and reliable prognostic and predictive biomarkers to guide individualized and precise immunotherapy selection. The American Joint Committee on Cancer (AJCC) TNM staging system has long been recognized as the primary clinical determinant of prognosis [7]. However, tumor heterogeneity due to differences in the molecular biology of CRC often leads to dramatically different outcomes in CRC patients with similar TNM characteristics [8], 9]. The current molecular classification systems for CRC are not sufficient to capture the complexity and variability of the immune microenvironment. Hence, the prediction of clinical outcomes and immunotherapeutic efficacy by a novel immune-associated molecular stratification method that can integrate multiple types of data and provide more accurate and comprehensive information remains an unmet need in CRC.

The use of large-scale gene sequencing technologies has improved the understanding of the genetic heterogeneity and genomic complexity of CRC, resulting in numerous molecular risk stratifications and transcriptomic subtypes characterizations of CRC [10]. Key genes in CRC have been demonstrated to influence the response to targeted drugs [11], 12]. As the understanding of the genomic and transcriptomic subtypes of CRC advances, the co-development of CRC biomarkers and their corresponding targeted agents is bound to have a transformative impact on the precision medicine of CRC. Since there are no generally accepted clinical risk stratification approaches, it may be of great significance to identify a unique molecular stratification for predicting the prognosis and response to immunotherapy in patients with CRC.

The aim of this study is to establish and validate an immune-associated molecular stratification model based on gene expression data from CRC samples. We will also explore the biological and functional characteristics of different immune-associated molecular subtypes and their associations with immune cell infiltration, immune checkpoint expression, tumor mutation burden, and neoantigen load. This study will provide new insights into the immune-related molecular mechanisms of CRC and facilitate the development of personalized immunotherapy strategies for CRC patients.

Materials and Methods

Results

NMF classifies CRC into two subclasses

A preliminary classification of CRC subclasses was conducted through the NMF consensus clustering, with two subclasses being stratified reasonably. The clinicopathological parameters between the two clusters are illustrated in Supplementary Table S1. To find the optimal k value, cophenetic correlation coefficients were calculated. The consensus matrix heatmap displayed a distinct and well-defined boundary, underscoring the stable and robust clustering of the samples (Figure 1A). Data showed that the cophenetic correlation coefficient fell sharply when k=2 (Figure 1B). Additionally, the PCA (Figure 1C and D) analyses confirmed that CRC patients in different risk clusters were distributed in two directions. As shown in Figure 1E and F, Cluster 1 (C1) was found to exhibit a worse OS and progression-free survival (PFS) than Cluster 2 (C2) (p<0.001).

Figure 1:

Construction of an NMF subtype based on the differentially expressed immune-related genes in the colorectal cancer (CRC) cohort. (A) CRC subclasses classification through NMF consensus clustering. (B) The relationship between cophenetic coefficients concerning the number of clusters (k=2 to 10). (C, D) Plane and three-dimensional principal components analysis (PCA) between the two clusters. (E, F) Kaplan-Meier survival curves showing the (E) OS and (F) PFS of the two molecular subgroups in the CRC cohort.

Differences in potential biological functions, pathways, and immune infiltration between the two clusters

We conducted a more in-depth exploration of potential biological mechanisms and established pathways between C1 and C2 using GSVA, which enabled us to identify the biological pathways perturbed and to establish the relationships between these disrupted processes. Our findings revealed that genes within C1 were enriched in the chemokine and immune function pathway. We proceeded to calculate the absolute abundance of 10 distinct cell populations between the two clusters, encompassing both the TME score and immune cell types. As illustrated in Figure 2A, the results showed that C1 had a significantly higher enrichment in the KRAS signaling, inflammatory response, IL6/JAK/STAT3 signaling, IL2/STAT5 signaling, TNFA/NFKB signaling, TGFβ signaling, and other immune-related signaling pathways, while C2 had a significantly higher enrichment in the MYC, oxidative phosphorylation, fatty acid metabolism, and peroxisome pathways. To further characterize the immune landscape, we compared the TME scores and the infiltration levels of 21 immune cell types across the two clusters. The two clusters exhibited significant differences in stromal, immune, and ESTIMATE scores, with all differences reaching a high level of statistical significance (p<0.001, Figure 2B). Furthermore, a comparative analysis of immune cell infiltration revealed significant variations in the expression levels of key cell populations, including specific CD4 T cell subsets (resting and activated memory), delta gamma T cells, naive and memory B cells, activated NK cells, monocytes, M0, M1, and M2 macrophages, as well as resting and activated dendritic cells and resting mast cells (Figure 2C). C2 exhibited distinct immunophenotypic features, including enrichment in anti-tumor immune pathways, elevated infiltration of activated CD8⁺ T cells and memory CD4⁺ T cells, and overall lower stromal scores. These characteristics were associated with improved OS and PFS compared to C1.

Figure 2:

Differences in potential biological functions, pathways, and immune infiltration between the two clusters. (A) Heatmap of gene set variation analysis between the molecular subgroups in the colorectal cancer (CRC) cohort. (B) The variations in immune cell infiltration between the molecular subgroups. (C) The abundance of each TME infiltrating cells among the two clusters. *p<0.05, **p<0.01, ***p<0.001.

Gene expression difference between the two clusters

A heat map and a volcano plot were plotted to compare DEGs between the two clusters (Figure 3A). The results showed that NMF classification resulted in 56 significant gene expression changes in the two clusters (Figure 3B, including 43 up-regulated genes and 13 down-regulated genes listed in Supplementary Table S2). Leveraging the entire dataset, we employed four well-established machine learning approaches (SVM-RFE, Lasso, XGBoost, and RF) to identify significant genes related to CRC subclasses, with 31, 30, 20, and 16 genes being yielded, respectively. As illustrated by Venn diagrams in Figure 3C, the intersection of screening characteristic genes via four types of machine learning algorithms and univariate Cox regression analysis among DEGs all pertained to two genes, SPP1 and GUCA2A, implicating their potentially essential role in CRC. Interestingly, it showed that the AUC of SPP1 was 0.902, significantly higher than that of GUCA2A (Figure 3D and E).

Figure 3:

Screening for differentially expressed hub genes by machine learning. (A) Heatmap and (B) volcano plot showing the gene expression differences between the two molecular subgroups in the colorectal cancer (CRC) cohort. (C) Venn diagram depicting the intersection of prognosis-related genes across different approaches. (D, E) Receiver operating characteristic curves of the hub gene to identify CRC subclasses.

High expression of hub gene SPP1 was associated with poor CRC prognosis

Our subsequent investigation focused on SPP1. To further validate the expression of SPP1 in CRC patients, tumor specimens derived from patients were collected and used for IHC staining of SPP1. As shown in the representative IHC staining images, SPP1 showed varying degrees of expression in CRC tumors. The overexpression rate of SPP1 was quantified to be 90.62 % (58/64), while 35.94 % (23/64) of the patients showed strong positive expression of SPP1 (Figure 4A). As depicted in Supplementary Table S3, the expression of SPP1 did not demonstrate any significant correlation with clinicopathologic parameters, including age, gender, tumor location, lymph node metastasis, neoadjuvant chemotherapy, AJCC stages, and tumor differentiation (p>0.05). As shown in Supplementary Figure S1, high expression of SPP1 was associated with more advanced staging of CRC. CRC cohorts were divided into high and low expression groups with respect to different SPP1 expression levels. As shown in the Kaplan-Meier survival curves, patients with high SPP1 expression in the tumor tissues demonstrated a poorer prognosis for OS (Figure 4B) and PFS (Figure 4C). As a validation, we also classified CRC patients into different expression groups based on the SPP1 protein expression. As illustrated in Figure 4D, high SPP1 expression was associated with a poor prognosis in CRC (p<0.01). These findings implied that SPP1 may play a crucial role and serve as a prognostic predictor for CRC patients.

Figure 4:

High expression of SPP1 was associated with poor colorectal cancer (CRC) prognosis. (A) Representative immunohistochemical staining images of SPP1 in sections derived from CRC patients. Scale bar: below (50 μm); above (200 μm). The bar graph indicates the proportion of patients with different expression levels. (B, C) Kaplan-Meier survival curves for OS and PFS of CRC cohorts with respect to different SPP1 expression levels. (D) Kaplan-Meier survival curves for PFS of CRC patients with respect to different SPP1 protein expression.

SPP1 expression correlates with CRC immune infiltration

To further investigate the relationship between SPP1 and immune components, we first estimated the immune cell fractions in the CRC cohort. Patients with high SPP1 expression were correlated with high InfiltrationScore (Figure 5A). In addition, correlation analysis revealed that SPP1 was positively correlated with neutrophil infiltration, whereas it was negatively correlated with CD8⁺ T cell infiltration (Figure 5B and C). To further corroborate this finding, IHC was carried out to determine the infiltration of CD8⁺ T cells and neutrophils in CRC tissues. As indicated by representative immunohistochemical staining, patients with high SPP1 expression were associated with more neutrophil infiltration and less CD8⁺ T cell infiltration (Figure 5D–G). Taken together, these results suggested that there is an association between SPP1 expression and immune cell fractions.

Figure 5:

SPP1 expression correlates with colorectal cancer (CRC) immune infiltration. (A) The abundance of each TME infiltrating cells among the two clusters. (B, C) Spearman correlation between SPP1 expression and immune infiltration in CRC. (D, E) Representative images of immunohistochemical staining of (C) neutrophils and (D) CD8⁺ T cells in tumor specimens from CRC patients with different SPP1 expression. Scale bars: above, 200 μm; below, 20 μm. (F, G) Immunohistochemistry scores (H-scores) quantification of immunohistochemical staining of (F) neutrophils and (G) CD8⁺ T cells in tumor specimens from CRC patients with different SPP1 expression. *p<0.05, **p<0.01, ***p<0.001.

Discussion

Despite the unprecedented achievements of cancer immunotherapy in recent years, the benefits are not available to every patient, and the mechanisms that lead to non-response have yet to be fully elucidated [21]. The molecular heterogeneity of CRC predicts not only clinical outcomes but also treatment sensitivity, especially for tumor characterization of the TME and immune profile [8], 9]. TNM staging, relying solely on macroscopic information, falls short in capturing the tumor heterogeneity arising from molecular biological variations in CRC. A “one-size-fits-all” approach to treating tumors of identical pathological and histological type and clinical stage is no longer aligned with the demands of precision medicine [22]. It was only with the widespread adoption of sequencing technologies that we gained insight into the variations in genomic, transcriptomic, and epigenetic characteristics among tumors of the same classification [23]. As an illustration, the CRC Subtyping Consortium has systematically categorized CRC into four consensus molecular subtypes (CMS), incorporating diverse pathological features and molecular subtypes [10], 24]. This nuanced classification acknowledges tumor heterogeneity, thereby opening avenues for personalized prognostic and therapeutic strategies. Recognizing molecular distinctions has, in turn, expedited the adoption of targeted therapeutic regimens, leading to improved outcomes across a wide spectrum of tumors. In our study, we calculated CRC immune-related gene index scores to classify CRC patients into two immune phenotypes with different prognoses. We assessed the prognostic significance of the two clusters and their associations with enriched pathways, tumor immune scores, and differential gene expression. Notably, these features demonstrated robust prognostic value and exhibited strong correlations with clinical metrics. Through effective stratification of CRC patients, key genes selected based on machine learning demonstrated their potential in predicting prognosis and presented a significant promise in advancing the development of targeted drugs and enhancing treatment approaches for CRC patients.

To identify the most crucial genes related to CRC, four supervised machine learning algorithms, SVM-RFE, Lasso, XGBoost, and RF, were utilized here to pinpoint the hub genes associated with CRC. As a proven machine learning algorithm for classification, SVM-RFE facilitates the optimal feature selection through recursive feature elimination [25]. Lasso, as a regression analysis method that concurrently conducts feature selection and regularization, aggregates the relative risks associated with significant variables and maximizes their predictive power. Its efficacy has been demonstrated through successful applications in numerous clinical studies [26], [27], [28]. XGBoost has demonstrated exceptional performance across a spectrum of tasks in recent years, owing to its user-friendly interface, ease of parallelization, and remarkable predictive accuracy [29]. Similarly, RF is a supervised non-parametric machine learning method well-suited for addressing both classification and regression issues, encompassing tasks such as gene screening and disease diagnosis [30]. In this study, using univariate Cox regression, SVM-RFE, Lasso, XGBoost, and RF algorithms, GUCA2A and SPP1 were identified as the potential biomarkers associated with the prognosis of CRC. Each of the four algorithms demonstrates distinct strengths across various computational scenarios. In the context of medical research, particularly oncology, it remains challenging to determine a universally superior method in terms of accuracy and precision. However, integrated machine learning models that combine multiple algorithms have shown promise in generating clinically meaningful predictions and improving prognostic accuracy in cancer studies.

Although formal cross-validation or external test set validation was not implemented in this study, several design strategies were adopted to reduce the risk of overfitting. First, feature selection was based on the intersection of four different machine learning models and univariate Cox analysis, ensuring consistency across multiple independent algorithms. The identification of SPP1 as the most robust gene was further corroborated by its high AUC and biological relevance. Second, SPP1’s prognostic and immunological significance was validated at both transcriptomic and protein levels, providing functional support for its selection. Finally, the integration of transcriptomic data from two independent cohorts (TCGA and GSE39582) across different platforms (RNA-seq and microarray) offers an implicit level of cross-dataset generalizability [31]. Although we integrated TCGA and GSE39582 datasets using the Rank-In method to improve robustness, the clustering was performed on the merged dataset rather than on completely separate cohorts. This is also one of the limitations of our study. Future studies should validate the reproducibility and clinical relevance of these subtypes in independent datasets to confirm their generalizability and potential for clinical translation.

SPP1 (secreted phosphoprotein 1), also known as osteopontin, is implicated in a number of physiological and pathological events, including tissue remodeling, inflammation, carcinogenesis and autoimmunity [32]. A growing body of evidence has demonstrated that SPP1 is overexpressed in a variety of tumors and is associated with poor tumor prognosis [32], [33], [34], [35]. SPP1 has also been found to activate intracellular pathways leading to the modulation of immune cells in the tumor immune microenvironment. In our study, we provided a preliminary dissection of the potential mechanisms underlying the role of SPP1 in the CRC immune microenvironment. The study suggested that SPP1 expression significantly affects the abundance of immune cell infiltration. In particular, our analysis revealed that high SPP1 expression was associated with decreased CD8⁺ T cell infiltration and vice versa in neutrophils. One intriguing finding that was found to be significantly associated with SPP1 was the infiltration of macrophages (Figure 5). This study supported evidence from previous observations, which showed that SPP1 macrophages act as tumor-specific macrophages contributing to the poor prognosis of CRC patients through pro-angiogenesis and formation of an immune-excluded desmoplastic structure [36], 37].

Although direct mapping to CMS was not performed in this study, the immune and stromal characteristics observed in our clusters suggest potential correspondence: Cluster 1 resembles features of CMS1/CMS4 (immune-high and mesenchymal), while Cluster 2 aligns more closely with CMS2/CMS3 (immune-low, epithelial) [10]. Future studies incorporating CMS annotations may further validate and contextualize our immune-based classification for clinical translation. Collectively, the infiltration of immune cells in CRC serves as an indicator of the immune status of patients and could, to some extent, potentially explain the disparity in survival outcomes observed between the two groups. All these results revealed the relevance of SPP1 to the elaboration of immune cells in the CRC microenvironment and highlighted the potential value of identifying and developing therapeutic strategies against this target and other molecules involved in its interactions to overcome immunosuppression and increase the response of CRC to immunotherapy.

Although GUCA2A was also identified as a candidate gene, its relatively lower AUC and limited prognostic impact led us to focus on SPP1 for downstream validation. GUCA2A (Guanylate cyclase activator 2A), a peptide hormone secreted by intestinal epithelial cells, plays a key regulatory role in guanylate cyclase 2C (GUCY2C) signaling through both autocrine and paracrine pathways [38], [39], [40], [41]. Loss of GUCA2A expression at the mRNA and protein levels represents one of the most frequently observed molecular alterations in CRC [41]. Silencing of its receptor, GUCY2C, has been associated with tumorigenic processes such as cellular transformation, unchecked proliferation, and genomic instability [42]. Despite its apparent relevance, current research on GUCA2A remains limited, and the underlying mechanisms through which it contributes to CRC pathogenesis have yet to be fully elucidated. Further exploration of GUCA2A’s functional role in CRC immune modulation may be warranted.

Despite the encouraging results obtained, there are several critical issues that require attention in the present study. One of the limitations encountered in our study stems from the extended median survival of CRC patients, coupled with the majority of enrolled patients for immunohistochemical validation being hospitalized in 2019. Consequently, we were unable to assess the long-term survival of these patients, which constitutes a noteworthy limitation. In addition, the data analyzed in our study were sourced from a public database. However, it is important to emphasize that further validation through multi-center clinical trials is imperative for establishing the robustness and generalizability of our findings in future research. Although only one external dataset (GSE39582) was used, future studies incorporating multiple independent cohorts are warranted to further validate and refine the proposed immune-associated classification. While our study has shed light on the potential role of the key genes in CRC, further investigation through basic experimental exploration is warranted to gain a deeper understanding of their functional significance and mechanisms of action.

Conclusions

We employed the NMF algorithm to categorize CRC patients into two immunophenotypically distinct clusters with prognostic implications and validated the unique profile associated with each. The identification of SPP1, a pivotal prognostic gene identified through machine learning, emerged as a crucial indicator for assessing tumor immune status, providing novel insights to enhance prognostic assessments and therapeutic strategies tailored to individual CRC patients.