Identification and functional enrichment of differentially expressed genes (DEGs) associated with CRC-related PANoptosisRaw CRC mRNA expression data were downloaded and integrated from the TCGA database, resulting in a total of 701 samples, which included 51 samples from the normal group and 650 samples from the tumor group. Differential genes between the two groups were identified using the Limma package (P < 0.05 and |log2FC|> 0.585). This resulted in a total of 5048 differentially expressed genes (DEGs), including 2639 upregulated genes and 2409 downregulated genes (Fig. 1A,B). PRGs were obtained from previous literature13. The intersection of the 5048 DEGs and PRGs resulted in 151 intersecting genes (Fig. 1C). Functional annotation of these 151 genes revealed enrichment mainly in signaling pathways such as necroptosis, apoptosis, and the NOD-like receptor signaling pathway (Fig. 1D,E).Fig. 1Identification and functional enrichment of DEGs associated with CRC-related PANoptosis. (A) Volcano plot of differential gene expression. Blue and pink indicate downregulated and upregulated differential expression, respectively (screening conditions: P < 0.05 and |Log2FC|> 0.585). (B) Heatmap of differential gene expression. (C) Venn diagram of differentially expressed genes (DEGs) in CRC and PANoptosis-related genes (PRGs). (D,E) Functional enrichment of intersecting genes. GO-KEGG enrichment analysis of intersection genes from the Metascape database (D). A cluster network of enriched pathways in which nodes that share the same cluster are often located close to each other (E).Screening of key PRGs associated with CRC progression and survival analysisTo further determine the key genes that affect progression of CRC, random survival forest analysis was performed on the 151 intersecting genes, and genes with a relative importance greater than 0.2 were selected as feature genes, revealing the importance sequence of 12 genes (Fig. 2A,B). Finally, associations between survival and expression of these 12 key genes were analyzed, and BCL10, CDKN2A, DAPK1, PYGM and TIMP1 were found to be significantly associated with survival, indicating that they are key genes (Fig. 2C–N).Fig. 2Screening of key PRGs associated with CRC progression and survival analysis. (A) Random survival forest analysis plot of 151 differentially expressed genes. (B) Importance ranking of 12 feature genes that met the criteria among intersecting genes. (C–N) Survival analysis of 12 characteristic genes. Kaplan–Meier analysis of the overall survival difference between BCL10 (C), CDKN2A (D), CLU (E), CTNNB1 (F), DAPK1 (G), GPX3 (H), GSN (I), GSTM1 (J), LEF1 (K), PIK3CD (L), PYGM (M) and TIMP1 (N) high- and low-expression groups.Identification of disease-causative regions for key PRGs in CRCCRC GWAS data were analyzed to confirm the disease-causing regions for the 5 key genes. A Q-Q plot shows the significant disease-associated single-nucleotide polymorphisms (SNPs) identified in the GWAS data (Supplementary Fig. 1A–F). Precise mapping of the GWAS data revealed the important SNP loci in the enriched regions. The disease-causing regions corresponding to the key genes are shown, with BCL10 located on chromosome 1, CDKN2A on chromosome 9, DAPK1 on chromosome 9 and PYGM on chromosome 11. The significant SNP loci corresponding to the key genes are shown in the Supplementary material (Supplementary GWAS data).Relationships between key genes and the immune microenvironmentThe microenvironment, which significantly affects disease diagnosis, survival outcomes and clinical treatment sensitivity, is composed mainly of fibroblasts, immune cells, the extracellular matrix, various growth factors, inflammatory factors and specific physicochemical characteristics. Distributions of immune infiltration levels and immune cell correlations are displayed in different forms in Fig. 3A,B. Furthermore, differences in immune cells between the two groups were explored, and significant differences were found in the numbers of naive B cells, plasma T cells, CD8 + T cells, activated memory CD4 + T cells, gamma delta T cells, resting NK cells, activated NK cells, monocytes, M0 macrophages, M1 macrophages, M2 macrophages, resting dendritic cells, resting mast cells, activated mast cells, eosinophils and neutrophils (Fig. 3C). Relationships between key genes and immune cells were further investigated, revealing that multiple key genes correlated highly with immune cells (Fig. 3D–H).Fig. 3The landscape of immune infiltration between CRC and normal groups. (A) Relative percentages of 22 immune cells across all samples. (B) Heatmap showing the correlation of infiltration of 22 immune cell types. The colored squares represent the strength of the correlation; red represents a positive correlation, whereas purple represents a negative correlation. The deeper the color is, the stronger the correlation is. (C) Differences in immune cell content between normal patients (blue) and patients with CRC (yellow). P < 0.05 was considered to indicate statistical significance. (D–H) Spearman correlations of BCL10 (D), CDKN2A (E), DAPK1 (F), PYGM (G), and TIMP1 (H) gene expression with immune cell content.Analysis of drug sensitivity for key genesThe effectiveness of surgery combined with chemotherapy for early CRC treatment is well established. We used drug sensitivity data from the GDSC database and the R package “pRRophetic” to predict the chemotherapy sensitivity of each tumor sample to explore the relationship between key genes and common chemotherapy drugs. The results showed that BCL10 expression correlated significantly with sensitivity to AKT inhibitors (Supplementary Fig. 2A). CDKN2A expression correlated significantly with sensitivity to AKT inhibitor VIII, dasatinib and gefitinib (Supplementary Fig. 2B). DAPK1 expression correlated with sensitivity to cisplatin, dasatinib, erlotinib and gefitinib (Supplementary Fig. 2C), and that of PYGM correlated with sensitivity to AKT inhibitor VIII, cisplatin, dasatinib, erlotinib, gefitinib and gemcitabine (Supplementary Fig. 2D). Sensitivity to AKT inhibitors VIII, cisplatin, dasatinib and gefitinib correlated significantly with TIMP1 expression (Supplementary Fig. 2E).Association of key genes with GSVA and GSEA in CRCOur research evaluated the specific signaling pathways in which the 5 key genes are enriched to explore the potential molecular mechanisms by which these genes affect progression of CRC. GSVA results revealed enrichment of signaling pathways such as PI3K_AKT_MTOR_SIGNALING and IL2_STAT5_SIGNALING with high expression of BCL10 (Fig. 4A), P53_PATHWAY and WNT_BETA_CATENIN_SIGNALING with high expression of CDKN2A (Fig. 4B), IL2_STAT5_SIGNALING and HEME_METABOLISM with high expression of DAPK1 (Fig. 4C), HEDGEHOG_SIGNALING and BILE_ACID_METABOLISM with high expression of PYGM (Fig. 4D), and TGF_BETA_SIGNALING and DNA_REPAIR with high expression of TIMP1 (Fig. 4E). In addition, GSEA showed that pathways such as the cAMP signaling pathway, necroptosis pathway and T-cell receptor signaling pathway are related to BCL10 (Fig. 4F,K), the Hippo signaling pathway, IL-17 signaling pathway and Wnt signaling pathway to CDKN2A (Fig. 4G,L), the PI3K-Akt signaling pathway, Rap1 signaling pathway and TGF-beta signaling pathway to DAPK1 (Fig. 4H,M), the cGMP-PKG signaling pathway, Hedgehog signaling pathway and oxytocin signaling pathway to PYGM (Fig. 4I,N), and apoptosis, cholesterol metabolism and the relaxin signaling pathway to TIMP1 (Fig. 4J,O). These findings suggest that progression of CRC may be affected by the key genes through these pathways.Fig. 4Functional and pathway enrichment analysis of key genes. (A–E) Correlation analysis results of GSVA for BCL10 (A), CDKN2A (B), DAPK1 (C), PYGM (D) and TIMP1 (E) in CRC. (F–O) Correlation analysis results of GSEA for BCL10, CDKN2A, DAPK1, PYGM and TIMP1 in CRC and molecular interaction networks between various pathways.Construction of a clinical prognostic prediction model for key PRGs in CRCMultivariate regression analysis revealed that in all of our samples, the values of different clinical indicators of CRC and the distribution of key gene expression contributed to the entire scoring process to varying degrees (Fig. 5A). Moreover, predictive analysis of OS was performed for three-year and five-year periods (Fig. 5B), and the results showed that the key gene-related nomogram model had good predictive performance.Fig. 5Establishment and validation of the prognostic nomogram. (A) Nomogram based on the key gene signature and clinical information for predicting 3- and 5-year OS in patients with CRC. (B) Calibration curves were used to verify the consistency of the predicted and actual 3- and 5-year outcomes (x‐axis: predicted survival probabilities; y‐axis: actual observed survival probabilities). OS, overall survival.Enrichment analysis of transcription factors for key genes and the miRNA networkThe 5 key genes were used as the gene set for this analysis, and it was found that they are regulated by common mechanisms such as multiple transcription factors. Therefore, enrichment analysis of transcription factors was performed using a cumulative recovery curve. The motif with the highest standardized enrichment score (NES: 9.09) was cisbp__M4892, as shown by motif-TF annotation, and the results of selection analysis for important genes and all the motifs enriched in the key genes and their corresponding transcription factors are displayed in Supplementary Fig. 3A–E. In addition, 71 miRNAs were obtained through reverse prediction of the 5 key genes using the miRcode database, resulting in a total of 146 mRNA-miRNA relationships, which were visualized using Cytoscape (Supplementary Fig. 3F).Correlation analysis between key genes and disease-regulating genesGenes involved in development and progression of diseases are called disease-regulating genes. In this study, tumor progression genes related to CRC were obtained from the GeneCards database (https://www.genecards.org/). Expression levels of the top 20 genes based on the relevance score and intergroup differences in expression of tumor genes were analyzed, revealing significant differences in expression of MSH6, BRCA2, BRCA1, MSH2, POLE, PMS2, APC, TP53, PALB2, CHEK2, POLD1, BRIP1, CDH1, AXIN2, BARD1, MUTYH, PTEN and KRAS between the two groups of patients. Furthermore, significant correlations were found between expression levels of the 5 key genes and those of multiple tumor progression genes, with DAPK1 and AXIN2 showing a significant negative correlation (r = − 0.504) and BCL10 and KRAS a significant positive correlation (r = 0.644) (Fig. 6A,B). Analysis of the impact of tumor immunotherapy revealed differences in expression of key genes between high- and low-expression groups (Supplementary Fig. 4A–E). (For responder analysis, the Hexp group refers to samples with gene expression greater than the median value, and the Lexp group refers to samples with gene expression not greater than the median value).Fig. 6Correlation analysis of CRC-related genes. (A) Differential analysis of disease-regulating genes. Differences in expression of CRC-related genes between normal (yellow) and with CRC (blue) samples. (B) Correlation analysis of key genes and differentially expressed regulatory genes. The first plot indicates that DAPK1 correlated significantly negatively with AXIN2, the second plot indicates the Pearson correlation between regulatory genes and key genes, and the third plot indicates that BCL10 correlated significantly positively with KRAS. Purple and red indicate negative and positive correlations, respectively. The Pearson correlation coefficients and P values are shown at the top of the graphs (P < 0.01 indicates a significant correlation).Single-cell sample subtyping and annotationSingle-cell data from GSE217774 were downloaded, and single-cell analysis was performed using the Seurat package. Cells were clustered using the tSNE algorithm, resulting in 30 subtypes obtained through tSNE (Fig. 7A). Then, annotation of each cluster was performed using the R package SingleR. All clusters were annotated into 7 cell categories: Monocyte, Macrophage, Epithelial_cells, Endothelial_cells, T_cells, Tissue_stem_cells and B_cell (Fig. 7B). In addition, expression of BCL10, CDKN2A, DAPK1, PYGM and TIMP1 in Monocyte, Macrophage, Epithelial_cells, Endothelial_cells, T_cells, Tissue_stem_cells and B_cell was explored (Fig. 7C–E).Fig. 7Single-cell analysis. (A) Cell clustering yielded 30 subtypes. (B) Annotating each subtype by “SingleR” from the R package. (C–E) Scatter plot (C), violin plot (D) and bubble plot (E) shows expression of five key genes in Monocyte, Macrophage, Epithelial_cells, Endothelial_cells, T_cells, Tissue_stem_cells and B_cell.Validation of the differential expression of key genes in CRC cells and tissues, and clinical significance verificationTo verify the results of our data analysis, we extracted total RNA from CRC cell lines (HCT116, SW480 and RKO) and the normal epithelial colon cell line NCM460 and measured mRNA expression levels of BCL10, CDKN2A, DAPK1, PYGM and TIMP1. qRT-PCR assays showed that the levels of BCL10, CDKN2A and TIMP1 were significantly increase in CRC cells than in NCM460 cells but that DAPK1 was expressed at low levels in cancer cells (Fig. 8A–C,E). In addition to those in RKO cells, mRNA levels of PYGM in HCT116 and SW480 cells were significantly increased (Fig. 8D).Fig. 8Validation of the differential expression of BCL10, CDKN2A, DAPK1, PYGM and TIMP1 in CRC cells and tissues. (A–E) mRNA expression levels of key genes in different cell lines (NCM460, HCT116, SW480 and RKO) were measured by qRT-PCR. The results were normalized to the reference gene GAPDH. Data are shown as means ± SDs, and two-tailed unpaired t tests were used for statistical analysis of each marker; n = 4 independent experiments. (F–J) The expression levels of key genes in CRC tissues (n = 30) and adjacent normal tissues (n = 30). (K–T) Analysis of the overall survival (OS) and disease-specific survival (DSS) of the key genes. (*P < 0.05, **P < 0.01, ***P < 0.001).Then, we detected expression levels of BCL10, CDKN2A, DAPK1, PYGM and TIMP1 in CRC tissues and paired adjacent normal tissues (Fig. 8F–J). Compared with adjacent normal tissues, the expression levels of BCL10, CDKN2A and TIMP1 were significantly increased in CRC tissues, whereas DAPK1 and PYGM expression levels were decreased (Fig. 8F–J). As expected, the expression changes of these key genes in CRC tissues were consistent with the cell results. Combined with clinicopathological features of CRC patients, BCL10 levels in CRC tissues were significantly associated with differentiation (P = 0.022) and perineural invasion (P = 0.043) (Supplementary Table 2), CDKN2A were associated with age (P = 0.001) and differentiation (P = 0.048) (Supplementary Table 3), DAPK1 were associated with tumor location (P = 0.001), dukes’ stage (P = 0.049) and lymphatic metastasis (P = 0.049) (Supplementary Table 4), PYGM were associated with differentiation (P = 0.020) (Supplementary Table 5), TIMP1 were associated with age (P = 0.010) and perineural invasion (P = 0.028) (Supplementary Table 6).More interestingly, the expression of 5 key genes in CRC tissue are also associated with the overall survival (OS) and disease-specific survival (DSS) of patients. Kaplan–Meier survival plots were used to analyze patients’ survival data from our study, TCGA-COAD and TCGA-READ projects. Our results showed that CRC patients with low expression levels of CDKN2A, PYGM and TIMP1 had longer OS and DSS times than those in high group (Fig. 8L, N,O,Q,S,T). On the contrary, CRC patients with high BCL10 expression level had longer OS and DSS times than those in low group (Fig. 8K,P). Additionally, DAPK1 expression level was associated with OS in CRC patients (Fig. 8M), but not associated with DSS (Fig. 8R).
