Developing a prognostic risk model based on circulating tumor cell genes to predict prognosis and provide potential therapeutic strategies in colorectal cancer

Introduction

According to statistics from the Global Cancer Epidemiology Database (GLOBOCAN 2020), it was estimated that there were 1.9 million new cases of colorectal cancer (CRC) and 935,000 deaths globally in 2020, making them the third and second most common malignant tumors, respectively (1). Furthermore, some studies indicate that metastasis accounts for the primary cause of mortality among CRC patients, with the liver serving as the most primary and common site of metastatic CRC (mCRC) (2-4). Consequently, the prognosis of CRC patients remains a clinically significant issue worth exploring.

In recent years, with the advancement of bioinformatics technology and the improvement of public databases, an ever-increasing number of prognostic markers for CRC have been discovered, and the prognostic models constructed on the basis of these markers can effectively predict the overall survival of CRC patients (5,6). However, the precise molecular mechanisms underlying CRC metastasis remain poorly understood. Circulating tumor cells (CTCs) are the direct cause of cancer metastasis, and these cells drive metastasis by disseminating from primary tumors to seed metastases in distant organs (7). Some studies suggest that these rare cells can serve as cancer biomarkers to predict prognosis or help identify potential therapeutic targets (8-12). Moreover, the differences in gene expression between primary CRC cells and circulating tumor cells are closely related to tumor metastasis (13). Therefore, further exploration of the differentially expressed genes (DEGs) in CTCs may provide a more accurate assessment of patient prognosis and facilitate precision therapy for CRC patients.

In this study, we downloaded gene datasets for CRC and CTCs from the public databases The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO). Through comprehensive analysis, we identified the DEGs between CRC and normal, as well as between CTCs and primary CRC, and extracted the intersecting genes for inclusion in the study. Subsequently, we used univariate Cox regression and least absolute shrinkage and selection operator (LASSO) regression to obtain eight CTCs’ DEGs. Based on these eight genes, we developed a prognostic model for CRC, which demonstrated favorable results in both the receiver operating characteristic (ROC) curves and the Kaplan-Meier curve analyses. We present this article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2024-2268/rc).

Methods

Data collection

The gene expression level data and clinical information for patients were obtained from the TCGA database (http://gdc.cancer.gov). The inclusion criteria for patient samples were as follows: (I) the dataset provides complete gene expression profiles, and (II) comprehensive clinical information, including age, sex, TNM staging, and overall survival. A total of 600 CRC samples were included. 49 normal samples were obtained by TCGA database. Additionally, the GSE72970 dataset from the GEO (https://www.ncbi.nlm.nih.gov/) database was utilized as a validation set, incorporating 124 CRC samples. From the GSE82198 dataset, three CTC samples and three primary CRC samples were included. Single-cell data were sourced from the GSE146771 dataset. Immunotherapy data were derived from the GSE91061 dataset, which includes complete clinical information for 101 patients with malignant melanoma, encompassing overall survival data, follow-up information, and immunotherapy efficacy data. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Screening for significantly DEGs

Based on the TCGA database and the GSE82198 datasets, differential analysis was conducted between CRC vs. normal and CTC vs. primary CRC groups using the “limma” package of R. Furthermore, the Benjamini & Hochberg method was employed for multiple testing corrections. The two groups of DEGs were screened with a threshold of adjusted P value <0.05 & |log2FC| ≥1. Then, the “clusterProfiler” package of R was used for Gene Ontology (GO) functional and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of the two groups of DEGs, with multiple testing corrections using the Benjamini & Hochberg method. DEGs were screened in both groups using thresholds of adjusted P value <0.05 and count >10. Finally, significantly DEGs were created by intersecting the two sets of screened DEGs.

Construction and validation of prognostic risk model

Significantly DEGs were subjected to univariate Cox regression analysis using the “survival” package of R. A threshold of P value <0.05 was chosen to select genes that are significantly associated with prognosis at the expression level. LASSO regression algorithm was performed on the above prognostically significant genes using the “glmnet” package of R. The penalty parameter was adjusted through 10-fold cross-validation to select key genes. Subsequently, the “survminer” package of R was used to conduct stepwise Cox regression analysis on the key genes to select the model genes. We multiplied the expression levels of the model genes by the corresponding regression coefficients obtained from the stepwise Cox regression analysis and performed an exponential operation to derive a comprehensive risk score. The specific formula is as follows:

Riskscore=h(t,X)=h0(t)∗exp(β1X1+β2X2+...+βnXn)[1]

In this formula, β represents the regression coefficients, h(t, X) represents the hazard rate associated with the covariate [gene expression level (X)] at a specific time (t), while h0(t) denotes the baseline hazard rate, reflecting the underlying risk in the absence of covariate effects. Using the risk score calculation formula, we calculated the risk score for each sample in the training set, which is the TCGA dataset. Subsequently, based on the median risk score, we divided the samples into a high-risk group (risk score greater than the median) and a low-risk group (risk score less than or equal to the median). To evaluate the prognostic predictive ability of different risk groups, we plotted Kaplan-Meier survival curves and time-dependent ROC (tROC) curves, and calculated their concordance index (C-index) values and area under the curve (AUC) values. Furthermore, we utilized samples from the GEO dataset as a validation set. Within the validation samples, we employed the stepwise Cox regression analysis to compute the regression coefficients corresponding to the model genes derived from the training set. Thereafter, the regression coefficients corresponding to the model genes and their expression levels were substituted into the risk score calculation formula to calculate the risk score for each sample. Based on the median risk score, we once again divided the samples into a high-risk group and a low-risk group, and plotted Kaplan-Meier survival curves and tROC curves to further validate the effectiveness of the model.

Immunological analysis

This study used the cell-type identification by estimating relative subsets of RNA transcripts (CIBERSORT) algorithm (14), microenvironment cell populations counter (MCPcounter) algorithm (15), and single-sample gene set enrichment analysis (ssGSEA) algorithm (16) to respectively assess the proportions of immune cell types in the analyzed tumor samples. Subsequently, based on the estimation of stromal and immune cells in malignant tumors using expression data (ESTIMATE) algorithm, we utilized the “ESTIMATE” package in R to calculate the ESTIMATE score, immune score, and stromal score. The differences in infiltration among different risk groups were assessed using the Wilcoxon test. Additionally, immune checkpoint gene expression data were extracted based on TCGA expression data. The expression differences of immune checkpoint genes among different risk groups were compared using the Wilcoxon test. Furthermore, the ssGSEA algorithm calculated hallmark gene sets enrichment scores. The correlation between the risk score and hallmark gene set enrichment scores was then analyzed.

Construction of nomogram

To further investigate the prognostic independence between clinical prognostic factors and risk group, we included clinical factors and risk group from the training set in univariate and multivariate Cox regression analyses. We used a P value threshold of less than 0.05 to screen for independent prognostic factors and to generate forest plots. Based on those above-screened independent prognostic factors, we constructed a nomogram using the “rms” package of R. Based on the nomogram model, calibration curves for one, three, and five years were plotted.

Drug sensitivity analysis

We estimated how sensitive each patient is to chemotherapy drugs using the Genomics of Drug Sensitivity in Cancer (GDSC; https://www.cancerrxgene.org/) database. We quantified the half-maximal inhibitory concentration (IC50) using the “pRRophetic” package of R. Then, we compared the IC50 differences of 138 chemotherapy drugs between different risk groups using the Wilcoxon test.

Prediction of immunotherapy response

We utilized the Tumor Immune Dysfunction and Exclusion (TIDE) database(http://tide.dfci.harvard.edu/) to predict each patient’s response to immune checkpoint therapy and represented it using TIDE scores. Based on the response of patients to immunotherapy in the GSE91061 dataset, patients were classified into progressive disease (PD), stable disease (SD) disease group, partial response (PR) group, and complete response (CR) group. We computed the risk score for patients in the GES91061 dataset using normalized raw counts. Subsequently, we stratified the samples into high and low-risk groups using the median risk score as the threshold. Finally, we analyzed the relationship between the efficacy of programmed death-ligand 1 (PD-L1) inhibitors and the risk grouping.

Correlation between prognostic models and epithelial-mesenchymal transition (EMT)

To explore the association between prognostic models and EMT, we used the Msigdb database (http://www.gsea-msigdb.org/gsea/msigdb/index.jsp) to identify relevant gene sets related to EMT using “HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION” as the search term. Using the ssGSEA algorithm, we calculated the EMT score, and compared the differences in EMT scores between high and low-risk groups. We then divided the groups into high and low EMT score groups using the optimal cutoff value. Finally, we utilized the Kaplan-Meier analysis of R to assess the association between the grouping of high and low EMT scores and actual prognosis information.

Single-cell analysis

We leveraged the Tumor Immune Single Cell Hub 2 (TISCH2), a single-cell RNA sequencing database focused on the tumor microenvironment (TME), to analyze the expression of model genes obtained from TCGA across different cell types. This further revealed variations in the TME among patients with CRC, thereby explaining the heterogeneity of CRC to some extent.

Statistical analysis

We identified the DEGs of CTCs through differential analysis. Subsequently, we obtained eight model genes using univariate and multivariate Cox regressions as well as LASSO regression. We used Kaplan-Meier curves and ROC curves at one-year, three-year, and five-year time points to evaluate the prognostic predictive ability of the model in the training cohort, and we further validated the feasibility of the model in the validation cohort using the same methods. Additionally, we constructed a nomogram that incorporated clinical risk factors and further explored the application of this prognostic risk model in the context of the immune microenvironment, immunotherapy, chemotherapy drug sensitivity, EMT levels, and single-cell analysis. Graph plotting and statistical analysis were conducted using R Studio 4.0.4. The distinction between the two groups was determined through either the paired two-tailed Student’s t-test or the Mann-Whitney-Wilcoxon test. Statistical significance was established at P<0.05.

Results

Clinical characteristics of patients.

As presented in Table 1, this study encompassed data from 600 CRC patients obtained from the TCGA database. The mean age of the cohort was 66.1 years, comprising 289 males and 311 females. In terms of the pathologic T classification, the distribution of patients was as follows: 20 patients were classified as T1, 118 as T2, 401 as T3, and 61 as T4. Regarding the Pathologic N classification, the patient distribution was as follows: 339 patients were categorized as N0, 153 as N1, and 108 as N2. According to the pathologic M classification, the study included 515 patients with M0 and 85 patients with M1. For the stage classification, the distribution was detailed as follows: 113 patients were in Stage I, 221 in Stage II, 181 in Stage III, and 85 in Stage IV.

Identification of DEGs and enrichment analysis

We identified 1,727 DEGs between CTCs and primary CRC samples from the TCGA database (Figure 1A), as well as 3,564 DEGs between CRC and normal samples from the GSE82198 dataset (Figure 1B). Subsequently, GO and KEGG pathway enrichment analyses were conducted for the DEGs in each group to explore the functional terms associated with the key genes. The enrichment analysis results displayed the top seven genes ranked by adjusted P value (Figure 1C,1D). Additionally, 409 intersecting genes were obtained between the two sets of DEGs (Figure 1E). Finally, a heatmap of the intersecting genes was generated, as shown in Figure 1F,1G.

Figure 1 Comparative analysis and functional enrichment of DEGs in CRC and CTCs. (A) Volcano plots comparing CTCs versus primary CRC group; (B) volcano plots comparing CRC versus normal group, examining the effect size (log₂FC)—log₁₀(P value), with blue and red points representing downregulated and upregulated factors, respectively; (C) enrichment bubble plots for DEGs from CTCs versus primary CRC group; (D) enrichment bubble plots for DEGs from CRC versus normal group, with the vertical axis indicating pathway names and the color intensity reflecting the significance of the P values; (E) Venn diagram depicting the intersection of DEGs between the two groups; (F) heatmap of intersected genes based on CRC versus normal group; (G) heatmap of intersected genes based on CTCs versus primary CRC group. BP, biological processes; CC, cellular components; CRC, colorectal cancer; CTCs, circulating tumor cells; DEGs, differentially expressed genes; FC, fold change; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, molecular functions.

Construction and validation of the prognostic model

We conducted a univariate Cox regression analysis on the intersecting genes, thereby identifying 49 genes significantly associated with the prognosis of CRC (Figure 2A). Using the LASSO regression algorithm, we filtered down to 24 key genes (Figure 2B,2C). Moreover, we utilized a stepwise Cox regression algorithm to derive the optimal gene combination, ultimately obtaining eight model genes (Figure 2D). We constructed a risk model based on eight model genes. In the training set, we calculated the risk score for each sample and divided the patients into high-risk and low-risk groups using the median risk score value of 0.95 as a threshold (available online: https://cdn.amegroups.cn/static/public/tcr-2024-2268-1.pdf). The distribution of patients’ risk scores and survival times is shown in Figure 3A. We plotted Kaplan-Meier curves to assess the association between the classification of high-risk and low-risk groups and the actual prognosis of CRC patients (Figure 3B). Based on the risk model, we plotted ROC curves for one year, three years, and five years, with AUC values of 0.686 (95% CI: 0.611–0.754), 0.712 (95% CI: 0.646–0.772), and 0.739 (95% CI: 0.654–0.822), respectively (Figure 3C). We also calculated the C-index for the risk model, which is 0.688. In the validation set, we also calculated the risk score for each sample and divided the samples into high-risk and low-risk groups using the median risk score value of 1.01 as a threshold (Table S1). The distribution of patients’ risk scores and survival times is shown in Figure 3D. We also plotted Kaplan-Meier curves to further validate the model (Figure 3E). Additionally, we generated ROC curves for one year, three years, and five years, with AUC values of 0.666 (95% CI: 0.564–0.762), 0.718 (95% CI: 0.615–0.819), and 0.737 (95% CI: 0.578–0.862), respectively (Figure 3F). The C-index value is 0.669 for the risk model.

Figure 2 Selection of model genes using LASSO and Cox regression. (A) Forest plot of genes significantly associated with the prognosis of colorectal cancer; (B) distribution of LASSO coefficients; (C) likelihood deviance of LASSO coefficients distribution, with two vertical dashed lines representing lambda.min (left line) and lambda.1se (right line); (D) forest plot of the multivariate Cox regression for the eight model genes. *, P<0.05; **, P<0.01; ***, P<0.001. AIC, Akaike information criterion; CI, confidence interval; LASSO, least absolute shrinkage and selection operator.

Figure 3 Prognostic analysis of colorectal cancer using risk scores in TCGA and GEO datasets. (A) Distribution of risk scores (upper plot) and survival status (lower plot) in the TCGA dataset; (B) Kaplan-Meier curves related to the prognosis of colorectal cancer based on the risk score prediction model in the TCGA dataset; (C) ROC curves for the gene prognostic features at one year, three years, and five years in the TCGA dataset; (D) distribution of risk scores (upper plot) and survival status (lower plot) in the GEO dataset; (E) Kaplan-Meier curves related to the prognosis of colorectal cancer based on the risk score prediction model in the GEO dataset; (F) ROC curves for the gene prognostic features at one year, three years, and five years in the GEO dataset. AUC, area under the curve; GEO, Gene Expression Omnibus; ROC, receiver operating characteristic; TCGA, The Cancer Genome Atlas.

Analysis of the immune microenvironment

In the training set, we employed the CIBERSORT algorithm to calculate the proportions of 22 types of immune cells in each sample. We then compared the differences in the proportions of various immune cells among different risk groups, identifying ten types of differentially infiltrating immune cells (DICs, Figure 4A). The ssGSEA algorithm was used to calculate the proportions of 28 types of immune cells in each sample. Similarly, we compared the differences in immune cell proportions among different risk groups, identifying 11 types of DICs (Figure 4B). Using the Estimate algorithm, we computed immune and stromal scores and analyzed the differences in these scores regarding infiltration levels among different risk groups (Figure 4C). The MCPcounter algorithm revealed significant differences in four types of immune cells (Figure 4D). Additionally, we compared the differences in the expression of immune checkpoint genes among different risk groups, identifying seven types of significantly different immune checkpoint genes (Figure 4E). As shown in Figure 4F, we computed the correlation between risk scores and enrichment scores of hallmark gene sets.

Figure 4 Immune analysis between different risk groups. (A) Comparison of immune cell types with significant differences between different risk groups using CIBERSORT algorithms; (B) comparison of immune cell types with significant differences between different risk groups using ssGSEA algorithms; (C) comparison of immune cell types with significant differences between different risk groups using ESTIMATE algorithms; (D) heatmap of immune cells with significant differences between different risk groups based on the MCPcounter algorithm; (E) comparison of immune checkpoint genes with significant differences between different risk groups; (F) heatmap of the correlation between risk scores and gene sets. -, P>0.05; ., P≈0.05; *, P<0.05; **, P<0.01; ***, P<0.001; ****, P<0.0001. CIBERSORT, cell-type identification by estimating relative subsets of RNA transcripts; ESTIMATE, estimation of stromal and immune cells in malignant tumors using expression data; MCPcounter, microenvironment cell populations counter; ssGSEA, single-sample gene set enrichment analysis.

Nomogram model based on independent prognostic factors

We conducted univariate and multivariate Cox regression analyses on clinical factors and Risk Group of CRC samples in the training set to identify significant independent prognostic factors, as shown in Figure 5A,5B. To further analyze the correlation between these independent prognostic factors (age, pathologic T, stage, and risk group) and survival prognosis, we incorporated them into the construction of the nomogram model, as depicted in Figure 5C. The calibration curves for one, three, and five years were plotted to validate the nomogram (Figure 5D).

Figure 5 Screening independent risk factors in the TCGA dataset and constructing and validating a nomogram model. (A) The forest plot of univariate Cox regression analysis. (B) The forest plot of multivariate Cox regression analysis. (C) Nomogram for predicting survival rates based on independent prognostic factors. (D) Calibration plot comparing the one-year, three-year, and five-year survival rates predicted by the nomogram with the actual survival rates. CI, confidence interval; TCGA, The Cancer Genome Atlas.

Drug sensitivity analysis

This study compared the sensitivity of patients in the training set to five common chemotherapy drugs, including cisplatin, docetaxel, gemcitabine, paclitaxel and vinblastine (Figure 6). The results showed that the IC50 of the high-risk group was significantly lower than that of the low-risk group, indicating that the high-risk group had a stronger sensitivity to chemotherapy drugs.

Figure 6 Comparison of IC50 levels for five chemotherapeutic agents—cisplatin, docetaxel, gemcitabine, paclitaxel, and vinblastine—among different risk groups. IC50, half-maximal inhibitory concentration.

Prediction of immunotherapy response

We compared the differences in TIDE scores between different risk groups, which revealed that patients in the high-risk group exhibited higher TIDE scores than those in the low-risk group (Figure 7A). Based on the GSE91061 dataset, we analyzed the relationship between the risk model and the response to immunotherapy. The results showed a significant association between high-risk patients and poor prognosis (Figure 7B). Additionally, the risk scores of the SD/PD (no response) group were significantly higher than those of the CR/PR (response) group (Figure 7C). As shown in Figure 7D, in the low-risk group, 30% of patients belonged to the CR/PR group, while 70% belonged to the SD/PD group. In the high-risk group, 10% of patients were in the CR/PR group, and 90% were in the SD/PD group.

Figure 7 Assessment of immune therapy responses in the GSE91061 cohort based on the risk model. (A) Comparison of TIDE scores among different risk groups; (B) Kaplan-Meier curves of immune therapy responses in the GSE91061 cohort; (C) box plots of immune therapy responses in the GSE91061 cohort; (D) bar charts of immune therapy responses in the GSE91061 cohort. CR, complete response; PD, progressive disease; SD, stable disease; PR, partial response; TIDE, Tumor Immune Dysfunction and Exclusion.

Estimation and clinical analysis of EMT scores

Our analysis compared the EMT scores between different risk groups (Figure 8A), demonstrating a significant elevation in EMT scores within the high-risk group compared to the low-risk group. Furthermore, Kaplan-Meier curves constructed based on EMT scores revealed a distinct trend: Patients with higher EMT scores exhibited notably poorer prognoses (Figure 8B).

Figure 8 Analysis of EMT score differences among different risk groups and survival analysis based on EMT score. (A) Comparison of EMT scores among different risk groups; (B) Kaplan-Meier curves for high-risk and low-risk groups based on EMT score classification; EMT, epithelial-mesenchymal transition.

Single-cell analysis

Utilizing the single-cell dataset GSE146771 from the TISCH database, we analyzed the expression of eight characteristic genes within the immune microenvironment. The GSE146771 dataset comprises 20 cell clusters and 13 types of immune cells (Figure 9A,9B), with Figure 9C depicting the distribution and quantity of various cell types. Analysis of the immune microenvironment indicated that eight model genes were predominantly expressed in various immune cells, with the expression of the sushi repeat containing protein x-linked (SRPX) gene being particularly prominent (Figure 9D,9E).

Figure 9 Single-cell analysis of 8 model genes based on the single-cell dataset GSE146771. (A) Annotation of all cell clusters; (B) annotation of all immune cell types; (C) percentage of all immune cell types; (D) the expression of 8 model genes; (E) expression patterns of the eight model genes. Exp. Level, expression level.

Discussion

Given the significance of CTCs in tumor metastasis, we innovatively linked CTC-related genes to the prognosis, immunotherapy, and chemotherapy of CRC patients. In this study, we identified two distinct sets of DEGs by analyzing TCGA (CRC and normal samples) and GSE82198 (CTC and primary CRC samples). Subsequently, we performed GO functional and KEGG signaling pathway enrichment analysis on the two groups of DEGs, ultimately obtaining 409 intersecting genes. We used univariate Cox regression, LASSO regression, and stepwise Cox regression analyses to identify eight genes. Subsequently, we established a risk score model based on these genes. Based on the median risk score, patients were divided into high-risk and low-risk groups. In the training set, Kaplan-Meier analysis revealed that the prognosis of the high-risk group was significantly worse than that of the low-risk group. The AUC values of the ROC curves at one year, three years, and five years were 0.686, 0.712, and 0.739, respectively, and the C-index of the model was 0.688. These results indicate that the model can effectively predict the prognosis of CRC patients and has potential value for personalized treatment. In the validation set, Kaplan-Meier analysis confirmed that the high-risk group had a worse prognosis. The AUC values of the ROC curves at one year, three years, and five years were 0.666, 0.718, and 0.737, respectively, and the C-index of the model was 0.669, which further validates the model’s feasibility. We performed univariate Cox regression analysis on clinical factors and risk groups in the training set, identifying age, pathological T stage, cancer stage, and risk group for the construction of the nomogram model. The results of the calibration curve indicate that the Nomogram model, which incorporates clinical risk factors, has a good ability to predict the prognosis of CRC patients, further demonstrating the broad applicability of risk models in clinical practice.

The significant achievements of immunotherapy in clinical trials have positioned it as a new cornerstone in cancer treatment. Immunotherapy has made significant strides in CRC patients with the deficient mismatch repair-microsatellite instability-high (dMMR-MSI-H) phenotype, but patients with proficient mismatch repair-microsatellite instability-low (pMMR-MSI-L) or microsatellite stable (MSS) phenotypes, which constitute the majority of patients with CRC, have shown little or no benefit (17). This condition can be attributed to different resistance mechanisms, such as the low tumor mutation burden (TMB) of pMMR/MSS tumors and the shift towards an immunosuppressive TME (18). In CRC, the heterogeneous environment consists of numerous cellular components that can facilitate an immunosuppressive microenvironment, including immune, stromal, and endothelial cells (19). TME is a dynamic and complex system that continuously and comprehensively influences different stages of tumor development (20). Moreover, various immune cell types of the TME are reported to affect cancer-cell progression and metastasis (21). Hence, in the training set, we analyzed the immune status of patients in different risk groups. We found that patients in the high-risk group had high stromal scores, and a high stromal level in tumors often indicates a more invasive tumor, leading to a poorer prognosis (22). The high-risk group had significantly more CD8⁺ T cells, Tregs, M0, and CD56⁺ NK cells. Some studies have shown that Tregs can promote the formation of an immunosuppressive microenvironment in CRC and inhibit effector T cells through a cyclooxygenase-2-prostaglandin e2-dependent mechanism, thereby promoting tumor growth (23,24). In contrast, in the low-risk group, there were significantly more CD4⁺ memory T cells, activated dendritic cells, activated CD4⁺ T cells, eosinophils, type 2 and type 17 T helper cells, memory B cells, and neutrophils. Memory cells can enhance surveillance against tumor cells, making it difficult for tumor cells to escape recognition and clearance by the immune system. We also found that the low-risk group had high expression levels of immune checkpoint genes. These findings indicate that the low-risk group has better immune function, aligning with the risk model’s predictions.

We calculated the TIDE scores for all CRC patients in the training set and found that the low-risk group exhibited lower TIDE scores. This indicates a relatively favorable tumor immune microenvironment in the low-risk group, where immune cells are more adept at recognizing and clearing tumor cells, potentially leading to enhanced efficacy of immune checkpoint inhibitor (ICI) therapy. In the GSE91061 dataset, we observed that patients in the low-risk group showed greater sensitivity to PD-L1. This result demonstrates the robust validation of our risk model derived from CRC in melanoma patients, suggesting its broad applicability and providing potential guidance for the immunotherapy of other cancers. Based on the GDSC database analysis, we evaluated patients’ sensitivity to five commonly used chemotherapy drugs. It was observed that high-risk patients exhibit a higher sensitivity to chemotherapy during treatment. These five chemotherapy drugs disrupt the nucleic acid metabolism or cell division processes of cancer cells, thereby inhibiting their proliferation and dissemination, ultimately leading to apoptosis. Moreover, tumors with a higher degree of malignancy typically exhibit a higher proliferation rate, rendering them more sensitive to chemotherapy drugs that impact cell proliferation. Consequently, this result further validates the importance and practicality of risk models in guiding clinical chemotherapy strategies.

EMT is an evolutionarily conserved developmental program that endows cancer cells with metastatic properties by enhancing their mobility, invasiveness, and resistance to apoptotic stimuli (25). Furthermore, the ease of survival and metastasis of CTCs is attributed to the acquisition of more mesenchymal traits through EMT (26). Therefore, we quantitatively analyzed the progression of EMT. The results showed that in the low-risk group, the EMT score was significantly lower than in the high-risk group, indicating that patients in the high-risk group had a significantly higher level of EMT progression than those in the low-risk group, further confirming the effectiveness of the risk model. Additionally, the Kaplan-Meier curve for the high EMT score group exhibited a significantly worse prognosis, consistent with the results of multiple studies (27,28). Finally, at the single-cell level, we discovered that these eight model genes are expressed in various immune cells, indicating their crucial role in immune response and maintenance of immune function. Additionally, the expression of the SRPX gene in immune cells is particularly prominent. The SRPX gene is a potential influencing factor in the liver metastasis of CRC, and inhibiting the expression of the SRPX gene can reduce the migratory ability of CRC cells (29). Therefore, the SRPX gene could be a potential target site. Modulating the expression or activity of the SRPX gene in immune cells may affect the function of immune cells and their immune response to tumors, thus potentially serving as a critical strategy for cancer immunotherapy.

Despite our comprehensive validation of the risk model, there are inherent limitations. Primarily, our findings are based solely on bioinformatics analysis and lack experimental validation. Secondly, we employed various algorithms to assess the sensitivity of chemotherapy drugs and immunotherapy to validate the prognostic performance of the model, yet further prospective research and clinical data are imperative for validation. Furthermore, our study only revealed a correlation between model gene expression and prognosis through statistical analysis without elucidating its underlying mechanisms.

Conclusions

In summary, this study established a risk prognosis model for CRC based on CTCs’ DEGs, enabling effective prognosis prediction for CRC patients. Furthermore, this model demonstrates accurate predictive capabilities in determining sensitivity to chemotherapy and immunotherapy, and it has also shown promising predictive abilities in the context of melanoma. The outcomes of this study may offer improved personalized treatment approaches for CRC and provide valuable insights for other tumor research endeavors.

Acknowledgments

None.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2024-2268/rc

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2024-2268/prf

Funding: This research was supported by Medical and Health Research Project of Zhejiang Province (grant number: 2022KY1093).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2024-2268/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Sung H, Ferlay J, Siegel RL, et al. Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J Clin 2021;71:209-49. [Crossref] [PubMed]
2. Okholm C, Mollerup TK, Schultz NA, et al. Synchronous and metachronous liver metastases in patients with colorectal cancer. Dan Med J 2018;65:A5524. [PubMed]
3. Tang J, Chen H, Wong CC, et al. DEAD-box helicase 27 promotes colorectal cancer growth and metastasis and predicts poor survival in CRC patients. Oncogene 2018;37:3006-21. [Crossref] [PubMed]
4. Engstrand J, Nilsson H, Strömberg C, et al. Colorectal cancer liver metastases - a population-based study on incidence, management and survival. BMC Cancer 2018;18:78. [Crossref] [PubMed]
5. Ye L, Zhang T, Kang Z, et al. Tumor-Infiltrating Immune Cells Act as a Marker for Prognosis in Colorectal Cancer. Front Immunol 2019;10:2368. [Crossref] [PubMed]
6. Chen L, Lu D, Sun K, et al. Identification of biomarkers associated with diagnosis and prognosis of colorectal cancer patients based on integrated bioinformatics analysis. Gene 2019;692:119-25. [Crossref] [PubMed]
7. Au SH, Storey BD, Moore JC, et al. Clusters of circulating tumor cells traverse capillary-sized vessels. Proc Natl Acad Sci U S A 2016;113:4947-52. [Crossref] [PubMed]
8. Yu M, Stott S, Toner M, et al. Circulating tumor cells: approaches to isolation and characterization. J Cell Biol 2011;192:373-82. [Crossref] [PubMed]
9. Hong B, Zu Y. Detecting circulating tumor cells: current challenges and new trends. Theranostics 2013;3:377-94. [Crossref] [PubMed]
10. Krebs MG, Metcalf RL, Carter L, et al. Molecular analysis of circulating tumour cells-biology and biomarkers. Nat Rev Clin Oncol 2014;11:129-44. [Crossref] [PubMed]
11. Paterlini-Brechot P, Benali NL. Circulating tumor cells (CTC) detection: clinical impact and future directions. Cancer Lett 2007;253:180-204. [Crossref] [PubMed]
12. Hou JM, Krebs MG, Lancashire L, et al. Clinical significance and molecular characteristics of circulating tumor cells and circulating tumor microemboli in patients with small-cell lung cancer. J Clin Oncol 2012;30:525-32. [Crossref] [PubMed]
13. Alix-Panabières C, Cayrefourcq L, Mazard T, et al. Molecular Portrait of Metastasis-Competent Circulating Tumor Cells in Colon Cancer Reveals the Crucial Role of Genes Regulating Energy Metabolism and DNA Repair. Clin Chem 2017;63:700-13. [Crossref] [PubMed]
14. Chen B, Khodadoust MS, Liu CL, et al. Profiling Tumor Infiltrating Immune Cells with CIBERSORT. Methods Mol Biol 2018;1711:243-59. [Crossref] [PubMed]
15. Becht E, Giraldo NA, Lacroix L, et al. Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. Genome Biol 2016;17:218. [Crossref] [PubMed]
16. Xiao B, Liu L, Li A, et al. Identification and Verification of Immune-Related Gene Prognostic Signature Based on ssGSEA for Osteosarcoma. Front Oncol 2020;10:607622. [Crossref] [PubMed]
17. Al Zein M, Boukhdoud M, Shammaa H, et al. Immunotherapy and immunoevasion of colorectal cancer. Drug Discov Today 2023;28:103669. [Crossref] [PubMed]
18. Ghiringhelli F, Fumet JD. Is There a Place for Immunotherapy for Metastatic Microsatellite Stable Colorectal Cancer? Front Immunol 2019;10:1816. [Crossref] [PubMed]
19. Egan H, Treacy O, Lynch K, et al. Targeting stromal cell sialylation reverses T cell-mediated immunosuppression in the tumor microenvironment. Cell Rep 2023;42:112475. [Crossref] [PubMed]
20. Chen F, Zhuang X, Lin L, et al. New horizons in tumor microenvironment biology: challenges and opportunities. BMC Med 2015;13:45. [Crossref] [PubMed]
21. Liu Y, Zhang Q, Xing B, et al. Immune phenotypic linkage between colorectal cancer and liver metastasis. Cancer Cell 2022;40:424-437.e5. [Crossref] [PubMed]
22. Huijbers A, Tollenaar RA, v Pelt GW, et al. The proportion of tumor-stroma as a strong prognosticator for stage II and III colon cancer patients: validation in the VICTOR trial. Ann Oncol 2013;24:179-85. [Crossref] [PubMed]
23. Yaqub S, Henjum K, Mahic M, et al. Regulatory T cells in colorectal cancer patients suppress anti-tumor immune activity in a COX-2 dependent manner. Cancer Immunol Immunother 2008;57:813-21. [Crossref] [PubMed]
24. Mougiakakos D, Johansson CC, Trocme E, et al. Intratumoral forkhead box P3-positive regulatory T cells predict poor survival in cyclooxygenase-2-positive uveal melanoma. Cancer 2010;116:2224-33. [Crossref] [PubMed]
25. Mittal V. Epithelial Mesenchymal Transition in Tumor Metastasis. Annu Rev Pathol 2018;13:395-412. [Crossref] [PubMed]
26. Wei C, Yang C, Wang S, et al. Crosstalk between cancer cells and tumor associated macrophages is required for mesenchymal circulating tumor cell-mediated colorectal cancer metastasis. Mol Cancer 2019;18:64. [Crossref] [PubMed]
27. Li S, Zhang J, Qian S, et al. S100A8 promotes epithelial-mesenchymal transition and metastasis under TGF-β/USF2 axis in colorectal cancer. Cancer Commun (Lond) 2021;41:154-70. [Crossref] [PubMed]
28. Li Q, Lai Q, He C, et al. RUNX1 promotes tumour metastasis by activating the Wnt/β-catenin signalling pathway and EMT in colorectal cancer. J Exp Clin Cancer Res 2019;38:334. [Crossref] [PubMed]
29. Wei B, Xiao S, Lou W. In silico whole-transcriptome analysis reveals a potential hsa_circ_0000375-miR-424-5p-TPM2/SRPX/SRGAP1 regulatory network related to liver metastasis of colorectal cancer. Heliyon 2023;9:e21688. [Crossref] [PubMed]