Exploring the expression levels of SA-beta-galIn the Pan-cancer cohort, we found the expression levels of GLB1 in rectum adenocarcinoma (READ) and colon adenocarcinoma (COAD) ranked 3rd and 5th among 33 cancers, respectively (Fig. 1a). Besides, compared with normal tissue, the expression levels of GLB1 were higher in READ and COAD (Fig. 1b). For further validation of GLB1 expression in CRC, we performed IHC staining for GLB1 in TMA from the RJH cohort. The representative IHC photographs of normal and tumor tissue were displayed in Fig. 1c, d. The quantitative result showed that tumor tissue had significantly higher GLB1 expression than that in normal tissue (Fig. 1e). These results indicated that cell senescence (GLB1 encoding SA-beta-gal, a biomarker of cell senescence) may have a pivotal impact on CRC.Figure 1Exploring the expression of GLB1 encoding senescence-associated beta-galactosidase in Pan-cancer and the RJH cohort. (a) The expression levels of GLB1 in Pan-cancer. (b) Differences in GLB1expression levels between normal and tumor samples in Pan-cancer. (c) Representative IHC photographs of GLB1 in normal patient from the RJH cohort. Left: magnification: ×40; Right: magnification: ×200. (d) Representative IHC photographs of GLB1 in CRC patient from the RJH cohort. Left: magnification: ×40; Right: magnification: ×200. (e) Comparison of immunohistochemical staining intensity of GLB1 between normal and CRC patients from the RJH cohort.Enrichment analysis of differentially expressed SAGs (DeSAGs)Based on 256 SAGs after data normalization, differential expression analysis was conducted to identify 204 DeSAGs and 184 DeSAGs between normal and tumor tissue samples in TCGA and GEO cohorts, respectively. Figure S1a, b showed top 20 DeSAGs with absolute maximal FC in the form of heatmap in the TCGA and GEO cohorts, respectively. 155 overlapping DeSAGs from two cohorts were then collected for subsequent enrichment analysis (Fig. S1c). GO enrichment analysis revealed that these DeSAGs were enriched in pathways related with aging, cell aging, cell senescence both in TCGA and GEO cohorts (Fig. S1d, e). Similarly, KEGG enrichment analysis showed these DeSAGs were enriched in senescence pathways, in addition to CRC, immunity, and infection (Fig. S1f, g). These results indicate that DeSAGs may have an impact on TME characterization of CRC.Consensus clustering of senescence patternsTo clarify the impact of DeSAGs on CRC, 1226 CRC samples in the pooled cohort were classified into 4 clusters with different senescence patterns based on DeSAGs. There were 451 CRC samples in cluster1, 451 CRC samples in cluster2, 193 CRC samples in cluster3, 131 CRC samples in cluster4, respectively (Fig. 2a). The cumulative distribution function (CDF) curves for k = 2–9 were plotted, in which k represented the number of clusters (Fig. S2a). And the relative change in area under CDF curves were demonstrated in Fig. S2b. The result of principal component analysis (PCA) showed there were boundaries between CRC samples from different clusters, suggesting CRC samples with different senescence patterns might have significantly different characterizations (Fig. 2b). Survival analysis revealed significant differences in OS between 4 clusters, with CRC samples in cluster3 having the worst OS (p < 0.001; Fig. 2c). Similar to OS results, K-M curves showed that CRC samples in cluster3 also had the worst PFS with p < 0.001 (Fig. S2c).Figure 2Consensus clustering for different senescence patterns. (a) Consensus clustering of CRC samples based on DeSAGs. (b) The result of PCA on CRC samples with different senescence patterns. (c) The overall survival (OS) comparison between CRC samples from different clusters. (d) The comparison of cancer-related hallmark between CRC samples from different clusters. (e) The comparison of TIDE score between CRC samples from different clusters. (f) Differences in immune cell infiltration between clusters predicted by the CIBERSORT algorithm. *: p < 0.05; **: p < 0.01; ***: p < 0.001. (g) Differences in biological processes between clusters predicted by GSVA. *: p < 0.05; **: p < 0.01; ***: p < 0.001.In this study, we further explored the differences in biological processes, immune cells infiltration and stromal characteristics between 4 clusters. The result of GSVA showed that CRC samples from different clusters had notable distinctions in cancer-related biological processes (Fig. 2d). For example, the processes of inflammatory response, EMT, angiogenesis and TGF-β signaling in cluster3 were more active compared to other clusters. In the comparison of Tumor Immune Dysfunction and Exclusion (TIDE) score, CRC samples from cluster3 had higher TIDE score than other CRC samples with p < 2.2e−16 (Fig. 2e), indicating CRC samples in cluster3 may have tumor immune dysfunction and rejection. CIBERSORT algorithm was used to predict the infiltration of immune cells in TME of CRC samples. As shown in Fig. 2f, the infiltration of immune effector cells, such as CD8+ T cell, plasma cell, as well as activated NK cell, was lower in the CRC samples in cluster3 than that in other clusters. Inversely, macrophages (including M0, M1, M2) were highly infiltrated in TME of CRC samples from cluster3. Finally, some known biological processes were quantified and compared between different clusters (Fig. 2g). Although antigen processing machinery, CD8+ T effector were more active in CRC samples from cluster3, the level of stromal activities such as Pan-F-TBRS, angiogenesis, EMT 2/3 were also higher in CRC samples from cluster3. The above findings demonstrated CRC samples with different senescence patterns had significant differences in prognosis, biological activity, and TME, providing a new avenue to understand the heterogeneity of CRC.Construction of senescence-related modelA total of 611 differently expressed genes (DEGs) were identified from 4 clusters with different senescence patterns (Fig. S3a). In the training set, univariate Cox regression analysis identified 117 prognostic DEGs (p < 0.05) for further model construction (Fig. S3b). The hazard ratio (HR) of each prognostic DEG was calculated, HR < 1 meant protective factor for prognosis, and HR > 1 meant risk factor for prognosis. The coefficients of prognostic DEGs were calculated by LASSO regression analysis (Fig. S3c). And 10-Fold cross-validation was performed to determine the 23 prognostic DEGs for construction of senescence-related model (Fig. S3d). The CNV information and position of model genes were displayed in Fig. S3e. In terms of mutation, DUSP27 was the gene with the highest mutation frequency, and missense mutation was the most common type of mutation (Fig. S3f).Clinical implications of senescence-related modelAccording to the coefficients and expression level of model genes, SRS of CRC samples were calculated. We could find that CRC samples in higher AJCC-TNM stage had higher SRS (Fig. 3a–c). Meanwhile, CRC samples in pathologically advanced stage also had higher SRS than that in early stage (Fig. 3d). These results suggested that there may be a strong association between SRS and survival of CRC. Therefore, we further investigated the difference in prognosis among CRC samples with different SRS. In the training set, CRC samples with high SRS had poorer OS than that with low SRS (p < 0.001; Fig. 3e). Similar outcomes were observed in the internal and external validation sets (Fig. 3f-g). In the training set, the area under receiver operating characteristic (ROC) curves of SRS in predicting 1-, 3-, 5-year OS was 0.710, 0.706 and 0.691, respectively (Fig. 3h). And Fig. 3i-j displayed the AUC of 1-, 3-, 5-year OS prediction in the internal and external validation sets. In addition to OS, we also investigated the role of SRS in predicting PFS of CRC samples. As shown in Fig. S4a–c, there were significant differences of PFS between low- and high- risk groups in the training, internal validation and external validation sets. Figure S4d, f showed the ROC curves of SRS in predicting 1-year, 3-year and 5-year PFS of CRC samples in the training, internal validation and external validation sets. Based on the above results, we can infer that SRS can be used as a reliable and effective auxiliary tool for the prognosis prediction of CRC samples. Moreover, mutation landscapes of low- and high-risk groups were displayed in Fig. 3k, l. In pan-cancer, SRS still had the ability to distinguish the prognosis of different samples (Fig. S5).Figure 3Clinical relevance of senescence-related model. (a–d) The comparison of SRS between CRC samples in AJCC-T stages, AJCC-N stages, AJCC-M stages, as well as pathological stages, respectively. And p < 0.05 is considered statistically significant. (e–g) Survival analysis between CRC samples from low- and high-risk groups in the training set, internal and external validation sets. And p < 0.05 is considered statistically significant. (h–j) ROC curves of SRS in predicting OS at 1-, 3-, 5-year in the training set, internal and external validation sets, respectively. (k) The mutation landscape of genes in CRC samples from the low-risk group. Top 20 genes with highest mutation frequency displayed in the waterfall plotting. (l) The mutation landscape of genes in CRC samples from the high-risk group. Top 20 genes with highest mutation frequency displayed in the waterfall plotting.Development of nomogram for survival predictionAfter univariate/multivariate Cox regression analysis, we found that age, pathological stage, T stage, M stage and SRS were independent prognostic indicators (Fig. 4a, b). Based on these independent prognostic indicators, a nomogram for quantifying the 1-year, 3-year, and 5-year survival probabilities of CRC samples was developed (Fig. 4c). Compared with single prognostic indicator, the AUC of nomogram was higher, indicating that nomogram had better specificity and sensitivity in predicting prognosis of CRC samples (Fig. 4d). The results of DCA demonstrated net benefit of age, pathological stage, AJCC-T/M stage, SRS and nomogram within a range of clinically reasonable risk thresholds (Fig. 4e). The net reduction curves showed the number of reduced patients after the interventions based on the decision derived from age, pathological stage, AJCC-T/M stage, SRS and nomogram (Fig. 4f). By comparison, we found that the nomogram had the greatest net benefit and net reduction, indicating nomogram is superior to other indicators in making better clinical decisions. Finally, calibration curves at 1-, 3- and 5-year were drawn to illustrate the difference between nomogram-predicted and actual survival. As shown in Fig. 4g–i, the nomogram-predicted survival is basically consistent with the actual survival. Therefore, nomogram may be a complementary tool to better assess the prognosis of CRC. In addition, this further indicated the importance of SRS in tumor prognosis, inspiring us to explore the mechanism of CRC from the perspective of cellular senescence.Figure 4A nomogram for predicting CRC prognosis. (a) The forest plot of univariate Cox regression analysis of SRS and clinical indicators. (b) The forest plot of multivariate Cox regression analysis of SRS and clinical indicators. (c) A nomogram containing independent prognostic indicators, including age, AJCC-T stage, AJCC-M stage, pathological stage and SRS. (d) ROC curves of independent prognostic indicators. Higher AUC is more valuable in predicting prognosis. (e) Decision curve analysis (DCA) curves show net benefits of age, pathological stage, AJCC-T stage, AJCC-M stage, SRS and nomogram at a range of clinically reasonable risk thresholds. (f) The net reduction analysis indicates how many patient deaths could be prevented after the intervention, as determined by age, pathological stage, AJCC-T stage, AJCC-M stage, SRS and nomogram. (g–i) Calibration curves of nomogram at 1-, 3- and 5-year.Predicting consensus molecular subgroups (CMSs) of CRC samplesTo explore the relationship between senescence patterns, SRS and tumor characterization, we predicted CMS of each CRC sample. As demonstrated in Fig. S6a, 171 CRC samples were classified into CMS1; 312 CRC samples were classified into CMS2; 172 CRC samples were classified into CMS3; 311 CRC samples were classified into CMS4. There were significant differences of biological processes between 4 CMSs. For example, microsatellite instability (MSI) of CRC samples in CMS1 were higher than that of other CMSs. Oppositely, CRC samples in CMS2 had feature of microsatellite stability (MSS). Compared with other CMSs, the activities of differentiation, glycolysis, fatty acids were more active in CMS3. And CRC samples in CMS4 were characterized by TGF-β signaling, as well as EMT. This study further investigated the proportion of CMSs in each senescence patterns. As expected, CMS4 had the most prominent proportion in CRC samples from cluster3, which was consistent with the result of GSVA (Fig. S6b). At the same time, we found CMS4 was the largest proportion in the high-risk group (Fig. S6c). The alluvial figure demonstrated the relationship between SRS, CMSs and senescence patterns (Fig. S6d). And the specific number of low- and high-risk CRC samples in each CMS was shown in Fig. S6e. These results suggested significant heterogeneity in CRC, such as immune activity, metabolism status and tumor pathway activity. However, senescence patterns can distinguish CRC with different characteristics. This provides a new reference for the precision treatment of CRC.Characterization of TME between CRC samples with different SRSTo investigate the effect of SRS on CRC TME, we explored the relevance of SRS to cancer-related biological processes and 7 critical steps in the cancer-immunity cycles. SRS was negatively related with CD8+ T effector, antigen processing machinery and immune checkpoint. Inversely, SRS was positively related with Pan-F-TBRS, as well as angiogenesis (Fig. 5a). In terms of cancer-immunity cycles, SRS was negatively related with all steps (Fig. 5b). These results suggested CRC samples with higher SRS may have an anti-tumor immune deficiency. Although there was no significant association between immune score and SRS (r = − 9.19e−03, p = 0.75; Fig. 5c), there was a positive correlation with stromal score (r = 0.37, p = 7.63e−41; Fig. 5d). Subsequently, the result of comparison of immune biomarker genes between low- and high-risk groups showed CRC samples with low SRS had higher expression levels of immune genes than that with high SRS (Fig. 5e). And GSEA further validated the similar results, CRC samples in low-risk group was enriched in immune-related pathways, including interferon alpha response and interferon gamma response (Fig. 5f). As for CRC samples in high-risk group, they were characterized by EMT, hypoxia, TGF-β signaling (Fig. 5g). The infiltration of 22 immune cells also illustrated immune effector cells, such as plasma cell, CD8+ T cell, activated NK cell were highly infiltrated in TME of CRC samples with low SRS. Inversely, the abundance of macrophage M0 and macrophage M2 in high-risk group were significantly higher than that in low-risk group (Fig. 5h). Hence, we inferred senescence-related model genes played an important role in immune cells infiltration in CRC TME, which in turn affected the response of CRC to immunotherapy.Figure 5The heterogeneity of CRC samples with different SRS. (a) The correlation of SRS with several known biological processes. Orange represents the positive correlation and blue represents the negative correlation. (b) The correlation of SRS with 7 critical steps in the cancer-immunity cycles. Orange represents the positive correlation and blue represents the negative correlation. (c) The relationship between SRS and immune score calculated by ESTIMATE algorithm. (d) The relationship between SRS and stromal score calculated by ESTIMATE algorithm. (e) Comparison of immune biomarker gene expression levels between low- and high-risk groups. (f) Enrichment of cancer-related hallmarks in CRC samples from the low-risk group. (g) Enrichment of cancer-related hallmarks in CRC samples from the high-risk group. (h) The differences in the infiltration of immune cells between the low- and high-risk groups.Investigation of senescence-related model at the single-cell levelIn the RJH cohort, we first performed quality control on scRNA-seq of 4 CRC samples. Cells with RNA number > 50 and mitochondrial RNA percentage < 10% were retained for subsequent analysis (Fig. S7a). 1500 genes with the largest variances were selected for PCA analysis (Fig. S7b). PCA analysis was conducted on these genes and top 20 PCs were selected for further dimensional reduction (Fig. S7c, d). Finally, cells were divided into 18 clusters (0–17 clusters) by t-SNE (Fig. S7e). all of these cells were annotated as T cell, B cell, epithelial cell, monocyte, neutrophils, NK cell, endothelial cell, and smooth muscle cell (Fig. 6a). AUC score of each cell was then calculated and projected onto the t-SNE diagram (Fig. 6b). AUCell scoring algorithm was applied to assign all cells into low- and high-AUC groups. We found two peaks in the AUC score of the cells. When the AUC score threshold was set to 0.062, 6098 cells showed relatively higher AUC score (Fig. 6c). Compared with high-AUC group, the proportion of B cells in the low-AUC group was higher, while the proportion of epithelial cells was opposite (Fig. 6d). To understand the role of model genes on B cells, we further performed analysis of scRNA-seq of B cells (Fig. S7f, g). The evolutionary trajectory of B cells with different clusters and states over pseudotime was displayed in Fig. S7h–j. As shown in Fig. 6e, B cells were annotated as 4 subsets, including plasma cell, memory B cell, naïve B cell, and geminal center B cell. The density of 4 B cell subsets over pseudotime was displayed in Fig. 6f. At the beginning of evolution, germinal center B cell appeared, followed by memory B cell and naïve B cell. Finally, plasma cells evolved. The expression level of each model gene over pseudotime was displayed in Fig. 6g. The expression levels of model genes in different B cell states and subsets were shown in Fig. S8. The evolutionary trajectory demonstrated germinal center B cell had two distinct evolutionary directions. One direction evolved into memory cells, and the other direction evolved into plasma cells (Fig. 6h). As expected, the expression levels of senescence-related model genes in different cells states were significantly different, indicating model genes may play an important role in evolution of B cells (Fig. 6i). Exploring cell–cell communication found that epithelial cell highly expressed CXCL2, CXCL3, while EGR1, FOS, IER3, DUSP1, JUNB and NFKBIA were predicted as potential target genes (Fig. 6j). The ligand-receptor pairs with the top potential interactions were also predicted (Fig. 6k).Figure 6Exploration of senescence-related model in scRNA-seq of the Ruijin hospital cohort (RJH cohort). (a) Cell annotation of scRNA-seq from the Ruijin cohort. (b) The t-SNE plots based on the AUC score of cells. Cell subsets with high AUC score are highlighted. (c) AUC scoring of senescence-related model genes. The threshold was determined as 0.062 and the 6098 cells exceeded the threshold value. (d) The comparison of cell subset proportions between low and high AUC score groups. (e) Cell annotation of B cell subsets. (f) The density changes of B cell subsets over pseudotime. (g) The heatmap shows the expression level of model genes in B cells over pseudotime. (h) The evolutionary trajectory along the B cell subsets, which is started from germinal center B cell. (i) The role of model genes in evolution of B cell subsets. (j) Top-correlated ligand activities inferred in epithelial cells (tumor cells) according to NicheNet. Heatmap showing regulatory potential of top-correlated ranked ligands in epithelial cells (tumor cells) and the downstream target genes in surrounding cells. (k) The heatmap showing potential interaction of ligand-receptor pairs between epithelial cells (tumor cells) and surrounding cells.In GSE132465, cells from 23 CRC samples were classified into epithelial cell, stromal cell, T cell, B cell, mast cell, as well as myeloid cell (Fig. 7a). When the AUC score threshold was set to 0.054, 39,823 cells showed relatively higher AUC score (Fig. 7b). As shown in t-SNE plot, high-AUC score cells were mainly consisted of epithelial cell, stromal cell and myeloid cell (Fig. 7c), which indicated that CRC with high expression of model genes may be associated with activation of stromal and immune-suppressive activities. To this end, we investigated the cell–cell communication mechanisms of epithelial cell and surrounding cells according to the expression and downstream targets of ligand-receptor pairs. We found that epithelial cell showed high CCL3, CCL4 and ALCAM ligand activity (Fig. 7d). And the expressions of ligands in epithelial cell were shown in Fig. 7e. In addition, CCL3 and CCL4 encoding protein bound to receptor encoded by CCR1, whereas ligand encoded by ALCAM interacted with receptor encoded by CD6 (Fig. 7f), resulting in the expression of target genes such as CDKN1A, MMP9, and TGFB1 (Fig. 7g). These results indicated senescence-related model genes may play a non-negligible role in evolution of B cells. Meanwhile, epithelial cell developed interaction networks with surrounding cells, which potentially induced cellular senescence and promoted the remodeling of the TME. This also partly explained the difference in sensitivity of CRC to drug therapy, including chemotherapy and immunotherapy.Figure 7Validation of senescence-related model from a single-cell perspective. (a) Cell annotation of scRNA-seq from GSE132465. (b) AUC scoring of senescence-related model genes. The threshold was determined as 0.054 and the 39,823 cells exceeded the threshold value. (c) The t-SNE plots based on the AUC score of cells. Cell subsets with high AUC score are highlighted. (d) Top-correlated ligand activities inferred in epithelial cells (tumor cells) according to NicheNet. (e) Dot plots showing the intensity (dot intensity) and the expression percentage (dot size) of top- correlated ligands in epithelial cells (tumor cells). (f) Ligand-receptor pairs displaying potential interactions between epithelial cells (tumor cells) and surrounding cells. (g) Heatmap showing regulatory potential of top-correlated ranked ligands in epithelial cells (tumor cells) and the downstream target genes in surrounding cells.Prediction of response to drug treatment in CRC samples with different SRSBased on the positive correlation between SRS and AJCC-TNM stage and the effect of model genes on immune cells, this study further conducted mIF on CRC tissues of patients with different pathological stages. As shown in Fig. 8a, the infiltration of T cell (CD3) and B cell (CD19) in early CRC was higher than that in advanced stage. In addition, in the GEO cohort, CRC samples in high-risk group were more sensitive to anti-CTLA4 ICIs compared with that in low-risk group. On the contrary, CRC samples in low-risk group were more sensitive to anti-PD-1 ICIs compared with that in high-risk group (Fig. 8b). These results were validated in the TCGA cohort (Fig. 8c). To further validate the reliability of SRS in predicting response to immunotherapy, we compared the difference of SRS between non-responder and responder in the GSE3564, which contained immunotherapy information of melanoma patients. As expected, SRS of non-responders was significantly higher than that of responders (p = 0.033; Fig. 8d). The proportion of non-responders in high-SRS group was higher than that in low-SRS group (Fig. 8e). These results indicated CRC patients with low SRS were more likely to benefit from immunotherapy.Figure 8Responses of CRC samples with different SRS to drug therapy. (a) Multiplex immunofluorescence staining of CRC tissues in early stage and advanced stage. Yellow: CD3, far-red: CD19, red: PanCK, dark blue: DAPI. (b) The comparison of responses to anti-PD-1/ CTLA4 immunotherapy between low- and high-risk groups in the GEO cohort. (c) The comparison of responses to anti-PD-1/CTLA4 immunotherapy between low- and high-risk groups in the TCGA cohort. (d) The comparison of SRS between non-responders and responders. (e) The proportion of non-responders and responders in low- and high-SRS groups. (f) The comparison of SRS between CRC samples from chemotherapy-resistant and chemotherapy-sensitive groups. (g) The comparison of estimated IC50 of 5-FU between CRC samples from the low- and high-risk groups. (h) Top 6 potential compounds/drugs predicted to be positively or negatively associated with SRS in the PRISM database. (i) Top 6 potential compounds/drugs predicted to be positively or negatively associated with SRS in the CTRP database.In the GEO cohort, CRC samples with chemotherapy were classified into chemotherapy-resistant and chemotherapy-sensitive groups. As shown in Fig. 8f, SRS of CRC samples in chemotherapy-resistant group were higher compared with that in chemotherapy-sensitive group (p < 0.001). In the prediction of response to chemotherapy, IC50 of 5-FU in high-risk group was also higher than that in low-risk group (p < 0.001; Fig. 8g). These results suggested that CRC samples with high SRS may be difficult to benefit from chemotherapy.In order to complement drug treatment strategies of CRC, several compounds/drugs were predicted as potential treatment strategies for CRC samples with different SRS. Figure 8h showed top 6 compounds/drugs positively and negatively associated with SRS in the PRISM database. Similarly, we also predicted top 6 compounds/drugs positively and negatively associated with SRS in the CTRP database (Fig. 8i). Interestingly, we found the AUC of floxuridine (which can be rapidly bio-catabolized to 5-fluorouracil) and fluorouracil were positively related with SRS, suggesting that CRC patients with low SRS were more sensitive to fluorouracil. It further validated the result of Fig. 8g. All of these results indicated the potential value of SRS in the precision treatment of CRC and provided new potential therapeutic drugs for CRC.Validation of hub genes in senescence-related modelThis study performed protein–protein interaction network (PPI) analysis on 23 genes from senescence-related model, and network of these genes was displayed in Fig. 9a. Through cytoHubba, 6 genes were identified as hub genes, including CXCL13, MMP12, CXCL1, CXCL9, IL7R, and MMP3 (Fig. 9b). The relationship between the expression levels of hug genes and infiltrations of 22 immune cells was demonstrated in Fig. 9c.Figure 9PPI network of genes in senescence-related model. (a) PPI network of genes in model. Red line—indicates the presence of fusion evidence; Green line: neighborhood evidence; Blue line: cooccurrence evidence; Purple line: experimental evidence; Yellow line: text-mining evidence; Light blue line: database evidence; Black line: co-expression evidence. (b) Hub genes identified from genes in model. (c) Relationship between hub gene expression levels and immune cell infiltration. (d–i) Validation of expression levels of hub genes between normal cell line (NCM-460) and CRC cell lines (HCT 116, HT-29) by RT-qPCR. (d) CXCL13. (e) MMP12. (f) CXCL1. (g) CXCL9. (h) IL7R. (i) MMP3.In order to preliminarily validate the clinical practicality of senescence-related model, we compared the relative expression levels of these hub genes between normal (NCM-460) and CRC cell lines (HCT 116, HT-29) by RT-qPCR. Figure 9d–i showed the relative expression levels of CXCL13, MMP12, CXCL1, CXCL9, IL7R, as well as MMP3, respectively. As we can see, all hub genes are underexpressed both in HCT 116 and HT-29 compared to NCM-460, which is consistent with HR of hub genes in CRC prognosis calculated by univariate Cox regression analysis.
