A novel prognostic signatures based on metastasis- and immune-related gene pairs for colorectal cancer

Background: Metastasis remains the leading cause of mortality in patients diagnosed with colorectal cancer (CRC). The pivotal contribution of the immune microenvironment in the initiation and progression of CRC metastasis has gained significant attention.

Methods: A total of 453 CRC patients from The Cancer Genome Atlas (TCGA) were included as the training set, and GSE39582, GSE17536, GSE29621, GSE71187 were included as the validation set. The single-sample gene set enrichment analysis (ssGSEA) was performed to assess the immune infiltration of patients. Least absolute shrinkage and selection operator (LASSO) regression analysis, Time-dependent receiver operating characteristic (ROC) and Kaplan-Meier analysis were used to construct and validate risk models based on R package. CTSW and FABP4-knockout CRC cells were constructed via CRISPR-Cas9 system. Western-blot and Transwell assay were utilized to explore the role of fatty acid binding protein 4 (FABP4) / cathepsin W (CTSW) in CRC metastasis and immunity.

Results: Based on the normal/tumor, high-/low-immune cell infiltration, and metastatic/non-metastatic group, we identified 161 differentially expressed genes. After random assignment and LASSO regression analysis, a prognostic model containing 3 metastasis- and immune-related gene pairs was constructed and represented good prognostic prediction efficiency in the training set and 4 independent CRC cohorts. According to this model, we clustered patients and found that the high-risk group was associated with stage, T and M stage. In addition, the high-risk group also shown higher immune infiltration and high sensitivity to PARP inhibitors. Further, FABP4 and CTSW derived from the constitutive model were identified to be involved in metastasis and immunity of CRC.

Conclusion: In conclusion, a validated prognosis predictive model for CRC was constructed. CTSW and FABP4 are potential targets for CRC treatment.

Introduction

Colorectal cancer (CRC) is the third most common cancer globally, with approximately 0.7 million mortality cases in 2020 worldwide (1). Currently, distant metastasis remains the leading cause of CRC-related death. Metastatic CRC (m-CRC) is defined as metastatic cancer that has spread beyond the original CRC mass, with the most common sites of metastasis being the lymph nodes and liver (2). Surgical resection effectively cures most localized lesions of primary CRC, whereas up to 20% CRC patients have metastases with initial diagnosis, and approximately 25% of patients with AJCC stage I–II will develop metastases in the following years (3, 4). Therefore, targeted therapy and prognostic assessment of mCRC is a great challenge for global public health.

Although surgery can completely remove localized liver or lung metastases of mCRC, <20% of patients achieve a long-term cure through resection (5). Infiltration of lymph nodes and latent micrometastases are more common in mCRC and depend on systemic therapy with chemotherapy and immunotherapy combinations (6). In recent years, immunotherapy has demonstrated promising clinical results in the treatment of tumors, including mCRC patients (7, 8). In June 2020, the Food and Drug Administration (FDA) approved the immune checkpoint inhibitor pembrolizumab as a first-line treatment for the MSI-H/MMR-D metastatic CRC (9). Further exploration of the relationship between mCRC and immune microenvironment is helpful to explore potential prognostic indicators and guide clinical medication.

Relatively much attention has been given to the signature of genes, which are based on the phenotypes of cancer, to better predict tumor prognosis (10). Liang et al. identified a six-gene signature on stem cell characteristic to construct a novel prognostic marker for patients with colon adenocarcinoma (COAD) (11). Xu et al. constructed and validated a predictive model for lung adenocarcinoma based on the individualized characteristics of immune-related gene pairs to differentiate the response of lung adenocarcinoma patients to immunotherapy (12). Therefore, identifying tumor-specific biomarkers of prognosis or response to immunotherapy in cancer tissues would be of tremendous clinical value.

In this study, a three-gene-pair signature was developed based on public databases to predict prognosis of CRC, and its efficiency were achieved in multiple validation sets. In addition, the signature was adopted, which could avoid the difference in model effect caused by gene expression difference to a certain extent. This model may provide guidance for prognosis and targeted drug use in mCRC patients. Importantly, two of the genes, fatty acid binding protein 4 (FABP4) and cathepsin W (CTSW), among the prognostic model component genes, are closely associated with metastasis and immune microenvironment and may serve as targets for treating CRC.

Materials and methods

Data evaluation, extraction, and differentially expressed genes analysis

As the training set, transcriptomic profiling of 453 CRC patients containing complete clinical information with corresponding clinical data obtained from The Cancer Genome Atlas (TCGA) database was conducted. Subsequently, to validate the prognostic efficacy of the model, RNA-seq data from four independent cohorts, namely, GSE39582 (n=579), GSE17536 (n=177), GSE29621 (n=65), and GSE71187 (n=52). The R package edgeR was used to analyze differentially expressed genes. |logFC| ≥ 1.3, and p < 0.0001 were set as the thresholds to screen for CRC-associated genes, while |logFC| ≥ 1.0 and p < 0.05 were set as the thresholds to screen for immune- and metastasis-related genes.

Immune infiltration evaluation for patients

To screen genes associated with immune infiltration in CRC, we performed single-sample Gene Set Enrichment Analysis (ssGSEA). The R packages “GSVA” and “GSEABase” were used to evaluate the types of immune cells and the abundance of immune cell infiltration in the TCGA-CRC expression profile. Then, we divided the CRC samples into high-, medium-, and low-immune infiltration groups based on their immune infiltration levels. This was achieved through clustering and calculation of stromal cell and immune cell scores using the R packages “Sparcl” and “Estimate.” The R package “Pheatmap” was used for the presentation of immune infiltration group.

Gene pair construction and risk model establishment

Gene pairing methods are described in Hong et al. (13). Briefly, incorporated genes are paired as the form of A|B. If the expression level of A is higher than that of B, the pair is recorded as 1; otherwise, it is defined as 0. When the expression level of 0 or 1 is >20%, the gene pair is considered effective.

To construct prognostic risk model, univariate Cox regression combined with clinical data was first used to screen for prognostically relevant gene pairs (p < 0.05), followed by LASSO regression with Cox proportional risk regression analysis to fit the best predictive formulas. The R package “survivalROC” was used to calculate the time-dependent receiver operating characteristic (ROC) and corresponding area under the curve (AUC) of the model at 1, 3, and 6 years. The optimal cut-off value with ROC curve of 6 years was obtained as the cut-off point to distinguish high risk from low risk.

Identification and validation of risk model

To validate the prognostic efficacy of the risk model, Kaplan–Meier analysis based on the R package “survival” and “survminer” was used to assess the prognostic differences between high- and low-risk groups in the training and validation sets. In addition, correlation analysis and chi-square test were used to evaluate the correlation between risk scores and clinicopathological features, which were visualized by “ggpubr” package and “ComplexHeatmap” package.

Estimation of immune cells infiltrating

TCGA-CRC immune infiltration was assessed by a variety of methods including XCELL, TIMER, MCPcounter, EPIC, CIBERSORT-ABS, and CIBERSORT to explore the correlation between risk scoring and immune cells and presented by bubble plots. Wilcoxon signed-rank test was used to evaluate the difference in immune cell infiltration between the high- and low-risk groups, with a significance threshold of p-value < 0.05.

GSEA analysis

Gene set enrichment analysis (GSEA) was used to explore the association between FABP4/CTSW and molecular signatures associated with CRC malignant progression. TCGA-CRC data were sorted according to CTSW or FABP4 expression, and the top 25% and bottom 25% of patients’ transcriptome data were included in analysis. The GSEA Desktop Application and Molecular Signature Database (MSigDB) was acquired from https://www.gsea-msigdb.org/gsea/index.jsp.

Cell lines and culture

Human CRC cell lines HT-29 and HCT-116 were purchased from ATCC, and both were maintained in McCoy’s 5A medium containing 10% fetal bovine serum (FBS) and 1% streptomycin–penicillin at 37°C in a 5% CO₂ atmosphere.

CTSW and FABP4 knockout via CRISPR-Cas9 system

The sgRNA sequence specifically targeting CTSW or FABP4 (listed in Supplementary Table S1) was predicted by the GPP Web Portal and synthesized by Sangon Biotech (Shanghai) and subsequently ligated to the corresponding CRISPR-Cas9 vector lentiCRISPRv2 to form a complete plasmid. Subsequently, 293T cells, target plasmids, packaging plasmids, and membrane plasmids constituted the lentivirus packaging system and generated the respective knockout lentivirus. CRC cells were inoculated in six-well plates and transfected with the lentivirus and polybrene. After 7 days of puromycin screening, stable CTSW/FABP4 KO cells were constructed. Control cells were then transfected with a virus consisting of an empty plasmid lentiCRISPRv2. The knockdown efficiency was verified by Western blotting.

Western blotting

The treated CRC cells were collected and washed twice with cold PBS, lysed on ice by adding appropriate doses of lysis buffer (containing protease inhibitor with phosphatase inhibitor), and then centrifuged at 14,000 rpm for 15 min at 4°C to remove debris. The primary antibodies anti-FABP4 (Abcam, ab92501), anti-Cathepsin W (Abcam, ab191083), anti-E-cadherin (Proteintech, 20874-1-AP), anti-N-cadherin (Proteintech, 22018-1-AP), and anti-MMP9 (Proteintech, 10375-2-AP) were used at a 1:1,000 dilution. A human-reactive STING Pathway Antibody Sampler Kit (CST #38866) used for primary antibodies associated with the STING pathway.

Transwell assay

For transwell assay, 5×10⁴ CRC cells were resuspended in 200 μl of serum-free MEM and inoculated in a 24-well plate, and 500 μl of complete medium containing 10% FBS fetal bovine serum was added to the lower layer. Cells were incubated at 37°C with 5% CO₂ for 24 h and counted under microscopic staining with crystal violet.

Human tissues and immunohistochemistry

Thirty-five patients with CRC were included in the study (Table 1). Tissue microarray sections were dewaxed in xylene, hydrated in graded alcohol, and finally in a closed solution for immunohistochemical staining. Sections were exposed to anti-FABP4 (Abcam, ab92501) and anti-CD8 (Proteintech, 66868-1-Ig) primary antibodies overnight.

TABLE 1

Table 1 Clinical characteristics of the patients.

Results

Identification of differentially expressed mRNAs

The workflow of this study is represented in Figure 1. First, we performed differential expression analysis on 612 cases of transcriptome profiling in TCGA-CRC public database (N=44, T=568) and identified 7,980 differential expressed genes Figure 2A. Subsequently, we analyzed metastasis-related genes of CRC samples (non-m=385, m=73) and identified 469 differentially expressed genes (Figure 2B). To obtain the degree of immune infiltration, CIBERSORT and ESTIMATE algorithms were utilized to quantify the activity or enrichment level of immune cells in CRC tissues. CRC patients was classified into high-, medium-, and low-immune group according to the immune score and stromal score (Figure 2C). We further analyzed the differentially expressed genes between the high- and the low-immune group (H-immu=48, L-immu=136), and 3,177 differentially expressed genes were obtained (Figure 2D). The intersection of the above differentially expressed genes was taken, and 161 differentially expressed genes were finally identified, which were associated with CRC metastasis and immunity (Figure 2E).

FIGURE 1

Figure 1 Flowchart of the study. First, we screened the TCGA-CRC database for differentially expressed genes associated with immunity and metastasis; second, we paired the differentially expressed genes and performed machine learning to fit a prognostic signature containing three gene pairs for prognostic validation in the training set and four independent validation sets; third, we explored the correlation between the prognostic signature and clinicopathological characteristics, immune infiltration, and drug sensitivity; finally, we focused on the role of two genes, FABP4 and CTSW.

FIGURE 2

Figure 2 Prognostic signature construction. (A, B) Differential expression analysis for N vs. T, non-metastatic vs. metastatic of TCGA-CRC data; (C) ssGSEA analysis showing immunity cluster of CRC patients; (D) differential expression analysis for high- vs. low-immunity groups; (E) differential expression genes intersection of the three groups; (F, G) differentially expressed gene pairing and LASSO regression for prognostic signature construction; (H) forest plot for the three gene pairs comprising the risk signature.

Construction and validation of risk models

The 161 differential genes were converted into gene pairs by iterative loop and expression-based 0 or 1 assignment, and then, 35 prognostic gene pairs were screened by univariate cox regression. The subsequent LASSO regression with multifactor cox regression analysis had a total of three gene pairs included in the Cox proportional risk model (Figures 2F–H). The formula of the model-based risk score was as follows: risk score = 0.397×C6orf15|PCSK1- 0.749×CTSW|FABP4 + 0.679×SPRR1B|PCSK1 (Table 2).

TABLE 2

Table 2 Information of the three-mRNA pair.

To evaluate the efficacy of the model, we calculated the areas under curve (AUCs) of the receiver operating characteristic (ROC) curves at 1, 2, and 6 years and took the cut-off value corresponding to the most optimal AUC (AUC = 0.711, value_cut-off = 0.703) as the point to assess the risk of patients (Figure 3A). Based on the cut-off values, patients in the TCGA-CRC dataset were divided into high- and low-risk groups and subjected to Kaplan–Meier (K-M) analysis to investigate the prognostic performance of the model. The survival time of patients in the high-risk group was significantly shorter than that in the low-risk group (p < 0.001), suggesting that the three-gene-pair model had good prognostic prediction efficiency in the training dataset (Figure 3B). In addition, the prognostic performance of the model was verified by K-M analysis in the independent validation sets GSE17536 (n = 177, p = 0.009), GSE29621 (n = 65, p = 0.005), GSE39582 (n = 579, p = 0.004), and GSE71187 (n = 52, p = 0.045) (Figures 3C–F). These results suggested that the three-gene-pair model had good prognostic prediction efficiency in the training set and multiple independent validation sets. Simultaneously, ssGSEA was applied to immunologically assess the four validation sets and to explore the relationship between their risk scores and immunity. We found that GSE39582, GSE17536, and GSE29621 had higher infiltration of Antigen-presenting cells (APC) co-stimulation, macrophages, and other related immune factors in the high-risk group (Supplementary Figure S1).

FIGURE 3

Figure 3 Prognostic validation of risk signature. (A) Time-dependent ROC curve analysis for the TCGA-CRC dataset; (B–F) Kaplan–Meier survival analysis of the risk signature in the TCGA-CRC training set and four independent validation sets (GSE17536, GSE29621, GSE39582, and GSE71187).

Clinical evaluation of risk model

Next, to explore the relationship between the model and clinical characteristics in the validation set, Wilcoxon signed-rank test was used to analyze the correlation between risk groups and patient age, gender, tumor stage, and TNM stage (Figure 4A). Chi-square test results reveal association between risk score and metastasis (p = 0.002) (Figure 4B). We found that model-based risk scores were not associated with patient age (p = 0.64) (Figure 4C) and gender (p = 0.087) (Figure 4D), while tumor grade (p < 0.05) (Figure 4E), T stage (p < 0.05) (Figure 4F), N stage (p < 0.05) (Figure 4G), and M stage (p < 0.05) (Figures 4B, H) were significantly associated with risk. Considering that genes identified from datasets were involved in immune response, we investigated whether the risk model was associated with the tumor microenvironment by Spearman correlation analysis (Figure 5A). The high-risk group had a significantly higher degree of immune infiltration (Figure 5B), with more CD8⁺T cells, tumor-associated fibroblasts, macrophages, monocyte infiltration, and less CD4⁺T cell infiltration, as compared with the low-risk group (Figures 5C–G). Furthermore, in the analysis of the correlation between risk score and cancer-targeted drugs, we found that the IC50 of PARP inhibitors Rucaparib, Olaparib, and FH535 in the high-risk group was significantly lower than that in the low-risk group, suggesting that the high-risk group was more sensitive to PARP inhibitors (Figure 5H). The high-risk score was associated with lower IC50 of Cisplatin, Imatinib, and Erlotinib, which are commonly used in the clinic (Figure 5I).

FIGURE 4

Figure 4 Correlation of risk signature with clinicopathological features. (A) Heat map reveals correlation between high and low risk score groupings and clinicopathological characteristics; (B) Chi-square test reveals significant correlation between risk score and metastasis; scatter plot indicates that (C) age and (D) gender are not correlated with risk score; scatter plot indicates significant correlation between (E) stage, (F) T, (G) N, and (H) M staging and risk score. ^**p < 0.01; ^***p < 0.001.

FIGURE 5

Figure 5 Analysis of immune cell infiltration and drug sensitivity between high- and low-risk subgroups. (A) Six commonly used methods such as XCELL, TIMER, and EPIC to assess the relevance of immune-related cells to risk scores. (B–G) Differences in CD4⁺ T cells, CD8⁺ T cells, macrophages, and others in high- and low-risk subgroups. (H, I) Differences in IC₅₀ between high- and low-risk subgroups for PARP inhibitors and commonly used first-line CRC drugs.

CTSW and FABP4 as potential therapeutic targets for CRC

Among the five genes comprising the risk model, we focused on two genes, CTSW and FABP4, as potential targets for CRC immunity and metastasis based on preliminary experiments and available reports. CTSW is a cysteine protease reported to associate with the membrane in the endoplasmic reticulum of natural killer cells and cytotoxic T cells. Wilcoxon signed-rank test was used to analyze the correlation between the CTSW expression (group by median) and patient age, gender, tumor stage, and TNM stage. The results showed that CTSW expression was negatively correlated with stage and M and N stage in TCGA (Supplementary Figure S2A). TCGA-CRC data were ranked according to CTSW expression, and the top 25% high expression data and the bottom 25% low expression data were included in GSEA enrichment analysis. Intriguingly, in addition to being enriched in immune response and epithelial–mesenchymal transition (EMT)-related processes, differentially expressed genes were also found to be associated with DNA damage and double-strand breaks (Figure 6A). Furthermore, we predicted the correlation of CTSW with immune cells through the online website TIMER, and the results showed that CTSW was significantly associated with CD4⁺T, CD8⁺T, macrophages, and dendritic cells (Figure 6B). To experimentally validate our bioinformatics analysis results, we generated CTSW-KO HT-29 and HCT-116 cells by CRISPR-Cas9, and the knockout efficiency was verified by Western blotting assay (Figure 6C). Next, the relationship between EMT-related proteins and CTSW was explored. Compared with negative control (NC), the expression of N-cadherin and MMP9 was increased in CTSW-KO cells, while E-cadherin was opposite, which was also consistent with our previous screening results (Figure 6D). Transwell assay also confirmed the enhanced migratory ability of CTSW-KO CRC cells compared with control (Figure 6F; Supplementary Figure S2C). Given that our GSEA analysis results suggested an association between CTSW expression levels and DNA damage pathways within CRC tumors, we next explored potential relationships between CTSW and the cyclic GMP-AMP synthase (cGAS) stimulator of interferon genes (STING) pathway, which is an important innate immune response pathway that can be activated by tumor DNA spill. Western blotting results showed that CTSW knockout significantly decreased the protein expression of cGAS, p-STING, p-TBK1, and p-IRF3, suggesting the important role of CTSW in activating the immune response (Figure 6E). The above results suggest that CTSW is associated with CRC metastasis except for immune infiltration, suggesting its important role in the development of CRC.

FIGURE 6

Figure 6 The biological function of CTSW in CRC immunity and metastasis. (A) The TCGA-CRC dataset was divided into CTSW-high and CTSW-low expression group based on CTSW expression. GSEA analysis demonstrated that the differentially expressed genes were significantly associated with immunity, EMT, and DNA damage repair pathways; (B) TIMER analysis showed a significant positive correlation between CTSW and CD8⁺ T cells, CD4⁺ T cells, and neutrophils in CRC; (C) Western blot assay validated the knockout efficiency of CTSW in HT-29 and HCT-116 cells; (D) Western blot assay probed the expression of EMT-related proteins in CTSW-KO cells; (E) Western blot assay showed the expression of cGAS-STING-related genes in CTSW-KO cells; (F) The migratory capability of CRC cells after CTSW knockout was tested by transwell assay. ^*p < 0.05; ^**p < 0.01; ^***p < 0.001.

FABP4 is a fatty acid binding protein related to fatty acid uptake, transport, and metabolism. Wilcoxon signed-rank test was used to analyze the correlation between the FABP4 expression (group by median) and patient age, gender, tumor stage, and TNM stage in TCGA. The results showed that FABP4 expression was correlated with T stage (Supplementary Figure 2B). We ranked the TCGA-CRC samples according to FABP4 expression levels and took each high and low 25% of transcriptome data for GSEA analysis. Similar to the results of CTSW, FABP4-related genes were not only associated with metastasis and immune response but also with DNA double-strand breaks and mismatch repair (Figure 7A). The results of online website TIMER also showed that FABP4 was strongly associated with CD4⁺T, macrophages, neutrophils, and dendritic cells (Figure 7B). Subsequently, FABP4-knockout CRC cells were constructed by CRISPR-cas9 system and verified by Western blotting assay (Figure 7C). It was found that E-cadherin expression was upregulated, while N-cadherin and MMP9 expression was downregulated in FABP4-knockdown CRC cells (Figure 7D). FABP4-knockout attenuated the migratory and invasion ability of CRC cells compared with NC, suggesting that FABP4 may promote CRC metastasis (Figure 7F; Supplementary Figure S2C). Considering that FABP4 is relevant to DNA damage repair, the association between FABP4 and cGAS-STING pathway was explored. The expression levels of cGAS, p-STING, p-TBK1, and p-IRF3 was significantly decreased in FABP4-KO cells (Figure 7E). Strikingly, FABP4 and CD8 expression in tissue microarray samples from 35 CRC patients were examined by immunohistochemistry, indicating that CD8 infiltration was increased in patients with raised FABP4 expression (Figure 7G). These findings suggest that FABP4 is closely related to immune response and metastasis and could be a potential therapeutic target for CRC.

FIGURE 7

Figure 7 The biological function of FABP4 in CRC immunity and metastasis. (A) The TCGA-CRC dataset was divided into FABP4 high- and low-expression group based on FABP4 expression. GSEA analysis showed that the differentially expressed genes were significantly associated with immunity, EMT, and DNA damage repair pathways; (B) TIMER analysis showed a significant positive correlation between FABP4 and CD8⁺ T cells, CD4⁺ T cells, and neutrophils in CRC; (C) Western blot validated the knockout efficiency of FABP4 in HT-29 and HCT-116 cells; (D) Western blot probed the expression of EMT-related proteins in FABP4-KO cells; (E) Western blot showed the expression of DNA damage-repair-related pathway (cGAS-STING)-associated genes in FABP4-KO cells; (F) transwell assay to explore the migratory ability of CRC cells after FABP4 knockout; (G) Immunohistochemical staining based on tissue microarrays from CRC patients to investigate the expression of FABP4 and CD8. **p < 0.01; ***p < 0.001.

Discussion

Previous studies have shown that the vast majority of genes associated with CRC metastasis and prognosis are expressed by cells in the tumor microenvironment, among which immune cells occupy an important position (14–16). In this study, we first used the ESTIMATE method to predict the degree of immune cell infiltration based on the TCGA-CRC data and combined with differential expression analysis to screen out genes associated with both immune and metastasis. Next, the differential genes were assigned in pairs for the prognostic prediction model construction followed by Cox regression and LASSO regression, and the prognostic prediction ability was validated in the training cohort and four independent validation cohorts. Subsequently, we evaluated the relationship of the risk score with clinicopathological features, immune infiltration, and immune cells. Finally, CTSW and FABP4 were explored as potential therapeutic targets for CRC.

Immunity is an essential ingredient of the tumor microenvironment. Tumor cells metastasize by evading monitoring by the immune system, which is the main cause of tumor death (17). EMT is an indispensable phenotype associated with metastasis, and tumor cells with EMT characteristics are present at the frontline of invasion (18). During EMT, specific cell surface proteins and cytoskeletal proteins are altered resulting in loss of epithelial cell polarity and increased invasiveness (19). Therefore, in order to obtain immune-related genes, we used the ssGSEA algorithm to evaluate the immune and stromal scores of the included samples based on 29 immune cell types and grouped samples by unsupervised clustering. The differentially expressed genes of high- versus low-immune group, metastatic versus non-metastatic group, and normal versus CRC tissue group were intersected. Intersecting genes were paired and constructed to obtain a prognostic prediction model prompting immune cell infiltration. It can be observed that patients in the model-based high-risk group had a worse prognosis and a higher degree of immune infiltration than those in the low-risk group.

In recent decades, advances in genomic research and the development of precision targeted therapies have significantly improved the prognosis of CRC patients including those with advanced disease (20). CRC drugs targeting VEGF, EGFR, BRAF V600E, PDL1, and others are currently being used in first-line treatment of CRC or included in clinical studies (21–23). Thus, we evaluated the IC50 of commonly used and targeted drugs for CRC in the high- and low-risk groups and found that the high-risk group was more sensitive to the first-line drugs cisplatin and oxaliplatin, DNA-repair-related PARP inhibitors (Rucaparib, FH535), and EGFR inhibitors (Erlotinib). Interestingly, our data revealed that CTSW and FABP4 in the model are both related to DNA damage repair pathways and involved in regulating the cGAS-STING pathway, which may partially explain the different sensitivity of high- and low-risk group to DNA damage-related inhibitors, but the specific mechanism needs to be further explored. Taken together, these results suggest that patients in the high-risk group may have a better response to immunotherapy and target therapy than low-risk patients.

Based on the expansion of cancer databases and the advancement of sequencing technology, increasing numbers of prognostic models are being established. Okuno et al. have reported a model based on eight-miRNAs that can robustly predict the risk of early recurrence of gastric cancer (24). Based on cuproptosis-associated regulators, a five-gene signature was constructed to predict outcome and responsiveness to immunotherapy in CRC patients (25). Sun et al. constructed and validated a programmed necrosis-related signature by analyzing the expression profile of necrotizing apoptosis-related genes in CRC (26). Nevertheless, these models require the exact expression values of genes, which are susceptible to different detection methods and individual variation. Here, we refer to the method used by Hong et al. to construct an immune-related lncRNA model for hepatocellular carcinoma prognosis by randomly pairing two genes and assigning an overall value of 0 or 1 (13). This approach alternatively requires the relative expression of the two genes that constitute the gene pair.

In this study, all five genes included in the model play important roles in cancer. C6orf15 was reported to be associated with liver cancer prognosis, lymphoma susceptibility, and systemic lupus erythematosus (27–29). PCSK1 is a member of the chymotrypsin-like preprotein convertase family and associated with obesity and diabetes, while the family is involved in the regulation of immune cells in the tumor immune response (30–32). SPRR1B, a cytosolic protein of keratinocytes, is a marker of highly differentiated epithelial cells and has been reported as a potential target for predicting immunotherapeutic response in pan-cancer (33, 34). CTSW is a cysteine protease that is closely associated with NK cells and cytotoxic T cells (35). Currently, CTSW is reported to be a characteristic gene pertinent to the immune microenvironment of breast and endometrial cancers (36, 37), but neither its specific role nor its mechanism in cancer has been revealed. FABP4 is a novel adipokine that regulates inflammation and angiogenesis and plays a central role in controlling lipolysis and the development of diabetes (38). In CRC, FABP4 has been reported to be involved in metabolic reprogramming and has been associated with TNM staging, differentiation, and metastatic tropism to the liver or lung (39, 40). Nonetheless, the definite mechanism of FABP4 and its relevance to immunity remains to be explored in CRC. Through a review of the previous studies, we finally spotlighted the role of CTSW and FABP4 in CRC. Our results uncovered that both CTSW and FABP4 were positively associated with DNA damage repair and immune response (including the cGAS-STING pathway). CTSW was negatively correlated with the expression of EMT-associated genes (N-cadherin and MMP9) and cell migration ability, while the opposite result was observed for FABP4 in CRC. These results suggest the essential role of CTSW and FABP4 in CRC metastasis and immunity.

In conclusion, a prognostic model consisting of three metastasis- and immune-related gene pairs was identified and validated on four external datasets in CRC. FABP4 and CTSW may act as critical regulators during CRC progression.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Ethics statement

Written informed consent was obtained from the individual(s) for the publication of any potentially identifiable images or data included in this article.

Author contributions

The research project was designed by BP, MX and SW. BP, YY, WD, LS and MX conducted the experiments. Statistical analysis was performed by BP, YY, WD and LS. The first draft of the manuscript was written by BP and MX. All authors were involved in revising the manuscript critically. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by grants from the National Natural Science Foundation of China (Grant Nos. 82272629 and 82203489), Jiangsu Provincial Key Research and Development Plan (Grant No. BE2019614), Jiangsu Provincial Medical Key Discipline Cultivation Unit (JSDW202239), Key Project of Science and Technology Development of Nanjing Medicine (ZKX21042), Elderly Health Research Project of Jiangsu Province (Grant No. LR2021017), Specialized Cohort Research Project of Nanjing Medical University (NMUC2021013A), Nanjing Medical University Science and Technique Development Foundation Project to MX (No. NMUB20210199), and Project of Science and National Nature Science Foundation of China (No. 82002230).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2023.1161382/full#supplementary-material

Supplementary Figure 1 | The relationship between immune cells and risk scores in validation cohorts. After immune assessment of the validation set by ssGSEA, multiple immune cells including antigen-presenting cells were highly infiltrated in the high-risk score group in the (A) GSE39582, (B) GSE17536, and (C) GSE29621.

Supplementary Figure 2 | Relationship between CTSW/FABP4 and clinicopathological features and cell invasion ability. TCGA-CRC samples were divided into high and low groups according to the median expression of CTSW/FABP4, and the relationship between age, gender, stage, TNM stage and grouping were analyzed (A, B). Transwell assay to explore the invasion ability of CRC cells after CTSW/FABP4 knockout (C). ^*p < 0.05; ^**p < 0.01; ^***p < 0.001.

References

1. Sung H, Ferlay J, Siegel R, Laversanne M, Soerjomataram I, Jemal A, 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(3):209–249. doi: 10.3322/caac.21660

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Bien J, Lin A. A review of the diagnosis and treatment of metastatic colorectal cancer. Jama (2021) 325(23):2404–5. doi: 10.1001/jama.2021.6021

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Amin MB, Greene FL, Edge SB, Compton CC, Gershenwald JE, Brookland RK, et al. The eighth edition ajcc cancer staging manual: continuing to build a bridge from a population-based to a more "Personalized" approach to cancer staging. CA Cancer J Clin (2017) 67(2):93–9. doi: 10.3322/caac.21388

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Ciardiello F, Ciardiello D, Martini G, Napolitano S, Tabernero J, Cervantes A. Clinical management of metastatic colorectal cancer in the era of precision medicine. CA Cancer J Clin (2022) 72(4):372–401. doi: 10.3322/caac.21728

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Benson AB, Venook AP, Al-Hawary MM, Azad N, Chen YJ, Ciombor KK, et al. Rectal cancer, version 2.2022, nccn clinical practice guidelines in oncology. J Natl Compr Cancer Net. JNCCN (2022) 20(10):1139–67. doi: 10.6004/jnccn.2022.0051

CrossRef Full Text | Google Scholar

6. Cervantes A, Adam R, Roselló S, Arnold D, Normanno N, Taïeb J, et al. Metastatic colorectal cancer: esmo clinical practice guideline for diagnosis, treatment and follow-up. Ann Oncol Off J Eur Soc Med Oncol (2022) 34(1):10–32. doi: 10.1016/j.annonc.2022.10.003

CrossRef Full Text | Google Scholar

7. Casak SJ, Marcus L, Fashoyin-Aje L, Mushti SL, Cheng J, Shen YL, et al. Fda approval summary: pembrolizumab for the first-line treatment of patients with msi-H/Dmmr advanced unresectable or metastatic colorectal carcinoma. Clin Cancer Res an Off J Am Assoc Cancer Res (2021) 27(17):4680–4. doi: 10.1158/1078-0432.Ccr-21-0557

CrossRef Full Text | Google Scholar

8. Overman MJ, McDermott R, Leach JL, Lonardi S, Lenz HJ, Morse MA, et al. Nivolumab in patients with metastatic DNA mismatch repair-deficient or microsatellite instability-high colorectal cancer (Checkmate 142): an open-label, multicentre, phase 2 study. Lancet Oncol (2017) 18(9):1182–91. doi: 10.1016/s1470-2045(17)30422-9

PubMed Abstract | CrossRef Full Text | Google Scholar

9. André T, Shiu KK, Kim TW, Jensen BV, Jensen LH, Punt C, et al. Pembrolizumab in microsatellite-Instability-High advanced colorectal cancer. New Engl J Med (2020) 383(23):2207–18. doi: 10.1056/NEJMoa2017699

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Du S, Zeng F, Sun H, Liu Y, Han P, Zhang B, et al. Prognostic and therapeutic significance of a novel ferroptosis related signature in colorectal cancer patients. Bioengineered (2022) 13(2):2498–512. doi: 10.1080/21655979.2021.2017627

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Liang Y, Su Q, Wu X. Identification and validation of a novel six-gene prognostic signature of stem cell characteristic in colon cancer. Front Oncol (2020) 10:571655. doi: 10.3389/fonc.2020.571655

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Xu F, Zhan X, Zheng X, Xu H, Li Y, Huang X, et al. A signature of immune-related gene pairs predicts oncologic outcomes and response to immunotherapy in lung adenocarcinoma. Genomics (2020) 112(6):4675–83. doi: 10.1016/j.ygeno.2020.08.014

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Hong W, Liang L, Gu Y, Qi Z, Qiu H, Yang X, et al. Immune-related lncrna to construct novel signature and predict the immune landscape of human hepatocellular carcinoma. Mol Ther Nucleic Acids (2020) 22:937–47. doi: 10.1016/j.omtn.2020.10.002

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Calon A, Espinet E, Palomo-Ponce S, Tauriello DV, Iglesias M, Céspedes MV, et al. Dependency of colorectal cancer on a tgf-B-Driven program in stromal cells for metastasis initiation. Cancer Cell (2012) 22(5):571–84. doi: 10.1016/j.ccr.2012.08.013

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Calon A, Lonardo E, Berenguer-Llergo A, Espinet E, Hernando-Momblona X, Iglesias M, et al. Stromal gene expression defines poor-prognosis subtypes in colorectal cancer. Nat Genet (2015) 47(4):320–9. doi: 10.1038/ng.3225

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Isella C, Terrasi A, Bellomo SE, Petti C, Galatola G, Muratore A, et al. Stromal contribution to the colorectal cancer transcriptome. Nat Genet (2015) 47(4):312–9. doi: 10.1038/ng.3224

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Garner H, de Visser KE. Immune crosstalk in cancer progression and metastatic spread: a complex conversation. Nat Rev Immunol (2020) 20(8):483–97. doi: 10.1038/s41577-019-0271-z

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Ganesh K, Massagué J. Targeting metastatic cancer. Nat Med (2021) 27(1):34–44. doi: 10.1038/s41591-020-01195-4

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Massagué J, Ganesh K. Metastasis-initiating cells and ecosystems. Cancer Discovery (2021) 11(4):971–94. doi: 10.1158/2159-8290.Cd-21-0010

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Bando H, Ohtsu A, Yoshino T. Therapeutic landscape and future direction of metastatic colorectal cancer. Nat Rev Gastroenterol Hepatol (2023). doi: 10.1038/s41575-022-00736-1

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Heinemann V, von Weikersthal LF, Decker T, Kiani A, Vehling-Kaiser U, Al-Batran SE, et al. Folfiri plus cetuximab versus folfiri plus bevacizumab as first-line treatment for patients with metastatic colorectal cancer (Fire-3): a randomised, open-label, phase 3 trial. Lancet Oncol (2014) 15(10):1065–75. doi: 10.1016/s1470-2045(14)70330-4

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Venook AP, Niedzwiecki D, Lenz HJ, Innocenti F, Fruth B, Meyerhardt JA, et al. Effect of first-line chemotherapy combined with cetuximab or bevacizumab on overall survival in patients with kras wild-type advanced or metastatic colorectal cancer: a randomized clinical trial. Jama (2017) 317(23):2392–401. doi: 10.1001/jama.2017.7105

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Benson AB, Venook AP, Al-Hawary MM, Arain MA, Chen YJ, Ciombor KK, et al. Colon cancer, version 2.2021, nccn clinical practice guidelines in oncology. J Natl Compr Cancer Net. JNCCN (2021) 19(3):329–59. doi: 10.6004/jnccn.2021.0012

CrossRef Full Text | Google Scholar

24. Okuno K, Watanabe S, Roy S, Kanda M, Tokunaga M, Kodera Y, et al. A liquid biopsy signature for predicting early recurrence in patients with gastric cancer. Br J Cancer (2023) 128(6):1105–16. doi: 10.1038/s41416-022-02138-1

PubMed Abstract | CrossRef Full Text | Google Scholar

25. He R, Zhang H, Zhao H, Yin X, Lu J, Gu C, et al. Multiomics analysis reveals cuproptosis-related signature for evaluating prognosis and immunotherapy efficacy in colorectal cancer. Cancers (2023) 15(2):387. doi: 10.3390/cancers15020387

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Sun M, Ji X, Xie M, Chen X, Zhang B, Luo X, et al. Identification of necroptosis-related subtypes, development of a novel signature, and characterization of immune infiltration in colorectal cancer. Front Immunol (2022) 13:999084. doi: 10.3389/fimmu.2022.999084

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Oh S, Kwon HJ, Jung J. Estrogen exposure causes the progressive growth of sk-Hep1-Derived tumor in ovariectomized mice. Toxicol. Res (2022) 38(1):1–7. doi: 10.1007/s43188-021-00100-6

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Chung SA, Brown EE, Williams AH, Ramos PS, Berthier CC, Bhangale T, et al. Lupus nephritis susceptibility loci in women with systemic lupus erythematosus. J Am Soc Nephrol JASN (2014) 25(12):2859–70. doi: 10.1681/asn.2013050446

CrossRef Full Text | Google Scholar

29. Skibola CF, Bracci PM, Halperin E, Conde L, Craig DW, Agana L, et al. Genetic variants at 6p21.33 are associated with susceptibility to follicular lymphoma. Nat Genet (2009) 41(8):873–5. doi: 10.1038/ng.419

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Mahmoud R, Kimonis V, Butler MG. Genetics of obesity in humans: a clinical review. Int J Mol Sci (2022) 23(19):11005. doi: 10.3390/ijms231911005

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Rose M, Duhamel M, Rodet F, Salzet M. The role of proprotein convertases in the regulation of the function of immune cells in the oncoimmune response. Front Immunol (2021) 12:667850. doi: 10.3389/fimmu.2021.667850

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Zhou Y, Liu S, Liu C, Yang J, Lin Q, Zheng S, et al. Single-cell rna sequencing reveals spatiotemporal heterogeneity and malignant progression in pancreatic neuroendocrine tumor. Int J Biol Sci (2021) 17(14):3760–75. doi: 10.7150/ijbs.61717

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Zhang Z, Wang ZX, Chen YX, Wu HX, Yin L, Zhao Q, et al. Integrated analysis of single-cell and bulk rna sequencing data reveals a pan-cancer stemness signature predicting immunotherapy response. Genome Med (2022) 14(1):45. doi: 10.1186/s13073-022-01050-w

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Parikh AS, Puram SV, Faquin WC, Richmon JD, Emerick KS, Deschler DG, et al. Immunohistochemical quantification of partial-emt in oral cavity squamous cell carcinoma primary tumors is associated with nodal metastasis. Oral Oncol (2019) 99:104458. doi: 10.1016/j.oraloncology.2019.104458

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Brinkworth RI, Tort JF, Brindley PJ, Dalton JP. Phylogenetic relationships and theoretical model of human cathepsin W (Lymphopain), a cysteine proteinase from cytotoxic T lymphocytes. Int J Biochem Cell Biol (2000) 32(3):373–84. doi: 10.1016/s1357-2725(99)00129-6

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Wang Z, Liu H, Gong Y, Cheng Y. Establishment and validation of an aging-related risk signature associated with prognosis and tumor immune microenvironment in breast cancer. Eur J Med Res (2022) 27(1):317. doi: 10.1186/s40001-022-00924-4

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Chen P, Yang Y, Zhang Y, Jiang S, Li X, Wan J. Identification of prognostic immune-related genes in the tumor microenvironment of endometrial cancer. Aging (2020) 12(4):3371–87. doi: 10.18632/aging.102817

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Prentice KJ, Saksi J, Robertson LT, Lee GY, Inouye KE, Eguchi K, et al. A hormone complex of Fabp4 and nucleoside kinases regulates islet function. Nature (2021) 600(7890):720–6. doi: 10.1038/s41586-021-04137-3

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Montero-Calle A, Gómez de Cedrón M, Quijada-Freire A, Solís-Fernández G, López-Alonso V, Espinosa-Salinas I, et al. Metabolic reprogramming helps to define different metastatic tropisms in colorectal cancer. Front Oncol (2022) 12:903033. doi: 10.3389/fonc.2022.903033

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Zhang Y, Zhang W, Xia M, Xie Z, An F, Zhan Q, et al. High expression of Fabp4 in colorectal cancer and its clinical significance. J Zhejiang Univ Sci B (2021) 22(2):136–45. doi: 10.1631/jzus.B2000366

PubMed Abstract | CrossRef Full Text | Google Scholar