The workflow chart is shown in the Supplementary Data (Supplementary Figure S1). The public gene-expression data for transcriptome profiling and the corresponding clinical annotation were obtained from Gene Expression Omnibus (GEO) and The Cancer Genome Atlas (TCGA) database on May 1, 2021. There were four eligible CC cohorts of gene-expression data (GSE39582, GSE33113, and TCGA–Colon Adenocarcinoma [TCGA-COAD) (discovery data) and GSE17538 (independent validation data)]. We downloaded the raw microarray data form the Affymetrix Human Genome U133 Plus 2.0 Array of GEO database and the RNA sequencing data (fragments per kilobase of transcript million mapped reads (FPKM) value) of TCGA. We employed the “ComBat” algorithm in “SVA” package to adjust the batch effects from nonbiological technical biases among different CC RNA-seq data. And all of the RNA-seq data were adjusted for background adjustment and quantile normalization with robust multiarray averaging method in “affy” and “simpleaffy” packages. And the DNA sequencing of annotated somatic mutation of single-nucleotide polymorphisms (SNPs) and copy number variation (CNV) data for CC were also downloaded from TCGA. All CC samples were coded according to the third Edition of International Classification of Diseases for Oncology (ICD-O-3). And the exclusion criteria included patients with incomplete survival information and missing data on neoplasm histologic type.

From previous research, we identified a total of 27 PR genes presented in the Supplementary Data (Supplementary Table S1). All of PR genes were gathered from previous study and MSigDB database (Latz et al., 2013; Shi et al., 2015; Orning et al., 2018; Karki and Kanneganti, 2019; Li et al., 2021). For example, the previous study suggested that the caspase (CASP) family (CASP1, CASP3, CASP4, CASP5, CASP6, CASP8, and CASP9) was related to GSDMD, GSDMB, GSDMA, and GSDMC, which were significant for cancer cell pyroptosis (Shi et al., 2015; Li et al., 2021). And study showed that CASP3 and Granzyme B (GZMB) could help to convert cell apoptosis into pyroptosis (Orning et al., 2018). A protein–protein interaction (PPI) network for the differentially expressed genes (DEGs) was constructed with Search Tool for the Retrieval of Interacting Genes (STRING), version 11.0 (https://string-db.org/).

We built a novel PR molecular subtype based on the level of 27 PR genes identified from three CC cohorts. The unsupervised clustering analysis clustering algorithm was performed to estimate the patterns of pyroptosis regulation and classify the CC samples for further analysis. The stability and patterns of molecular clusters were adjusted by the consensus clustering algorithm (Wong, 1979). The “ConsensuClusterPlus” package was employed to cluster, and the process was performed 1,000 times (Wilkerson and Hayes, 2010).

To identify PR regulators genes, we need to estimate the expression level of different genes for studying the molecular feature among PR subtypes. We identified the DEGs with the empirical Bayesian approach in “limma” package, and we set the |log2-fold change| > 1 and false discovery rate (FDR) < 0.05 as the significance criteria.

To investigate the molecular feature among PR subtypes, we established gene set variation analysis (GSVA) enrichment analysis with “GSVA” R packages (Hänzelmann et al., 2013). The gene set of “c2. cp.kegg.v6.2. symbols” and “c5. all.v6.2. symbols.gmt” were gathered from the MSigDB database to be used in GSVA. H: Hallmark gene sets; C2: curated gene sets [including Kyoto Encyclopedia of Genes and Genomes (KEGG)] were downloaded from the MSigDB database to be used in gene set enrichment analysis (GSEA) with the software gsea 3.0. And we set the adjusted p < 0.05, nominal (NOM) p < 0.05, and FDR q < 0.05 as the statistically significance to identify the difference on biological process.

The Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data (ESTIMATE) was used to calculate the stromal score, immune-score, tumor-purity, and ESTIMATE-score for CC (Song et al., 2017). The enrichment levels of the 29 immune signatures were established based on the genes set from MSigDB database (Supplementary Table S1) with the single-sample GSEA (ssGSEA) (Ritchie et al., 2015; Bu et al., 2021). And the infection of 22 human immune cells in TME was established with cell-type identification by estimating relative subsets of RNA transcripts (CIBERSORT) web portal (https://cibersortx.stanford.edu/) and 1,000 permutations (Chai et al., 2019). The deconvolution algorithm output had a p-value <0.05 was set as successful and accurate deconvolution, which would be normalized to make their direct interpretation as cell fractions for comparison across different groups.

The Weighted Gene Co-expression Network Analysis (“WGCNA”) R package was employed to build the co-expression network of DEGs (van Houwelingen et al., 2006). The co-expression similarity matrix, Pearson’s correlation matrices, and average linkage method were involved in evaluating the correlations among the included genes. The Amn = |Cmn|β (Amn is theadjacency between gene m and gene n; Cmn, Pearson’s correlation between gene-m and gene-n; and β, soft thresholding parameter) could show that the strength of correlations contributes to the weighted adjacency matrix with a scale-free co-expression network. The topological overlap matrix (TOM) was used to identify the connectivity and dissimilarity of the co-expression network established with an appropriate β value.

The log-rank test and the Kaplan–Meier survival analysis were used to evaluate the difference in overall survival (OS) among different groups. We used the package “caret” to allocate all the CC patients in inner-training and inner-testing groups randomly through the 8:2 ratio, which contributed to enhance the generalization ability of model. The LASSO-penalized Cox regression model was used to evaluate the role of genes to identify signatures significantly associated with the patients’ OS. And the 10-fold cross validation was employed to prevent overfitting with the penalty parameter lambda.1se (Heagerty et al., 2000). The univariable and multivariate Cox regression analyses were used to identify the independent prognostic factors and to establish eight PR signatures and nomogram based on the forward and backward elimination methods. The area under the curve (AUC) and the time-dependent receiver operating characteristic (ROC) curve were used to evaluate the prognostic accuracy of the eight PR signatures model in inner-training and inner-testing groups with the package “survival ROC” (Pei et al., 2020). The PRM-scores were established based on the eight PR signatures model, and the median of PRM-scores was set as the cutoff value to the separate patients into high- and low-PRM-score groups. Bootstrap method was performed to validate the Cox model internally and externally. Bootstrap-corrected OS rates were calculated by averaging the Kaplan–Meier estimates based on 2,000 bootstrap samples.

A total of 27 PR regulators were identified in CC in this study with three eligible CC cohorts. We dissect the incidence of somatic mutations and molecular signatures of PR regulators in CC from TCGA-COAD (Figure 1A). The result showed that 109 of 590 CC samples experienced mutations of PR regulators, with frequency of 23.33%. It was found that the missense mutation exhibited the highest frequency variant classification. Both C>T ranked and SNPs were the most frequent alternatives in single-nucleotide variant (SNV) class and variant type. The NLRP7 exhibited the highest alteration frequency followed by SCAF11, while the PRKACA, CASP6, PYCARD, and TNF showed extremely low alteration frequency in CC samples (Figure 1B). To ascertain whether the above genetic variations influenced the expression of PR regulators in CC patients, we investigated the mRNA expression levels of regulators between normal and CC samples (Mann–Whitney U test; *p < 0.05; **p < 0.01; ***p < 0.001; p ≥ 0.05, not significant) (Figure 1C). The expression of CASP4, GASP8, GPX4, GSDMC, GZMB, IL1B, NOD1, NOD2, and PLCG1 was increased; while the expression of AIM2, CASP1, CASP3, CASP5, CASP6, CASP9, GSDMB, GZMB, IL18, NLRP1, and NLRP7 was decreased in CC samples compared with normal tissues. Correlation analysis was performed with genetic variation and expression variations of PR regulators in CC to further investigate the relationship among these regulators (left: genetic variation; right: expression variations) (Figure 1D). The correlation network containing all PR genes is presented in Figure 1E (red: positive correlations; blue: negative correlations).

To explore the potential biological molecular of PR regulators, we established a PR molecular subtype using consensus clustering analysis for CC patients. Three CC datasets with available clinical and follow-up information (GSE39582, GSE33113, and TCGA-COAD) were incorporated into one meta-cohort and clustered into three molecular subtypes (PR-A, PR-B, and PR-C) based on the expression of 27 PR regulators (Figure 2A). There are high intragroup correlations and low intergroup correlations in this classification pattern. There was also a significant difference in the survival among three subtypes (Figure 2B). The results of survival analysis proved that the OS of the PR-B and PR-C groups was significantly lower than that of the PR-A group according to the Kaplan–Meier curves of the CC cohorts (log-rank test, p < 0.01, Figure 2B). The expression of 27 PR regulators was different in three subtypes (ANOVA test, *p < 0.05; **p < 0.01; ***p < 0.001; p ≥ 0.05, not significant) (Figure 2C). In order to further portray the biological characteristics of these distinct molecular subtypes, we established GSVA enrichment analysis, including the KEGG and Gene Ontology (GO). The PR-A showed enrichment in terms of pathways associated with immune activation, including IL-17 production, T cell-mediated cytotoxicity, T cell-mediated, T-cell chemotaxis, and T-cell migration and differentiation. PR-B presented enrichment pathways including the proximal tubule bicarbonate reclamation, nitrogen metabolism, and tyrosine phosphorylation of STAT5 protein. While the enrichment pathways in PR-C were associated with immune suppression, including downregulation in natural killer (NK) cell activation involved in immune response, B-cell proliferation, and T-cell activation involved in immune response.

In addition, we tend to estimate the immune microenvironment among the PR molecular subtypes. The TME cell infiltration characteristics were calculated with ESTIMATE, including the tumor purity and immune-scores (ANOVA test, *p < 0.05; **p < 0.01; ***p < 0.001; p ≥ 0.05, not significant) (Figure 3A). The result showed that the immune-scores and ESTIMATE were the highest in PR-A among three subtypes, which suggested that the PR-A presented a high level of immune fully activation. The highest stromal-scores and tumor purity were in PR-C, and the lowest immune-scores were in PR-C, which suggested that the PR-C may characterized by the suppression of immunity. To investigate the proportions and differences of tumor infiltrating immune cell subsets among PR regulators subtypes, we employed a deconvolution algorithm with the CIBERSORT method (Figure 3B, Supplementary Figures 1A,B). The results noted that there were significant differences on the compositions of TME cell types among the three PR subtypes, which suggested that PR regulators may influence the types of TME infiltrating cell in CC. We found that the infiltration of activated immune cell in TME was abundant in PR-A, including the presence of CD8 T cells, activated NK cells, and B cells (ANOVA test, *p < 0.05; **p < 0.01; ***p < 0.001; p ≥ 0.05, not significant) (Figure 3C), which were same with the immune-scores from ESTIMATE. The high level of immunity may be related to the significant survival advantage (Bai et al., 2020). The PR-B was enriched with M1 macrophages, dendritic cells, plasma cells, and CD8 T cells. And the PR-C was enriched with M2 macrophages, naive B cell, CD4 T-cell memory resting, and T-cell regulatory cells (Tregs). The PR-C was reached with M2 macrophages, resting dendritic cells, and Tregs. And we quantify the enrichment levels of immunity related pathways and immune cells in CC via ssGSEA with a total of 29 immune-associated gene sets (Supplementary Figure S2). There was a significant difference in level of HLA genes among three subtypes. The checkpoint, CD8 T cells, HLA, MHC, and TILs were the highest in PR-A, which suggested the potentially ability for immune-inflamed. Based on the characterization of TME cell infiltration and biological molecular, PR-A was classified as immune-activated phenotype, with abundant immune cell infiltration and survival advantage; PR-B was classified as intermediate phenotype; and PR-C was classified as immune-excluded phenotype, characterized by the low immune response and high tumor purity. But the type of TME immune cells was the same among different subtypes, which showed that the PR regulators may regulate the level of immune cell infiltration and that they could not influence the types of cells in TME.

To further reveal the role of PR subtypes for prognosis and treatment of CC and apply the clusters to guide subsequent treatment, we established risk assessment tool-constructions based on the PR subtypes. All the genes were analyzed for co-expression network analysis using the WGCNA package (Figure 4A, Supplementary Figure S3). The association was built among the expression of gene and the PR clusters and clinical information based on the three eligible CC cohorts of gene-expression data (GSE39582, GSE33113, and TCGA-COAD). A total of 18 modules were identified; and the ME in the brown, yellow, red, and pink modules showed significantly higher association with PR regulators clusters than other modules in CC. From these modules, we identified 854 signature genes associated with the PR regulators (p < 0.05), which were selected for further analysis. Next, we estimated the independent prognostic signature of these genes using univariate Cox regression analysis, and the p-value <0.05 was considered to be the cutoff criteria. Patients from TCGA-COAD, GSE33113, and GSE39582 were randomly divided into inner-training and inner-testing groups through the 8:2 ratio. And we set GSE17538 as the independent validation cohort. Next, we established the LASSO-Cox regression model and cross validation to calculate the mean-squared error of genes with independent prognostic factors (Figure 4B). Eight genes, cytotoxic T-lymphocyte-associated protein 4 (CTLA4), chemokine (C-C motif) ligand 11 (CCL11), ninein (NIN), transmembrane protein 154 (TMEM154), kinesin family member 7 (KIF7), KIAA1671, ribonuclease P/MRP 14-kDa subunit (RPP14), and cadherin 19 (CDH19), were identified with the LASSO-Cox regression model and multivariate Cox regression analysis, which were used to establish the PRM (Figure 4B). All of these genes had significant independent prognostic factors in multivariate Cox regression analysis (Figure 4C). Besides these, eight genes expression were different in three PR subtypes (ANOVA test, *p < 0.05; **p < 0.01; ***p < 0.001; p ≥ 0.05, not significant) (Figure 4D). And CTLA4, CCL11, NIN, TMEM154, KIAA1671, RPP14, and CDH19 expressions were associated with immune-scores (Figure 4E). The prognostic index formula for CC was as follows: PRM-scores = [Status of CTLA4 * (−0.27274) + Status of CCL11 * (−0.05312) + Status of NIN * (0.30814) + Status of TMEM154 * (−0.21183) + Status of KIF7 * (0.55555) + Status of KIAA1671 * (−0.15928) + Status of RPP14 * (−0.34418) + Status of CDH19 * (0.45252)]. We divided colon patients into high- and low-PRM-score groups based on the median value, which was set as the cutoff value to divide the patient into high or low group in the validation cohorts.

The survival analysis suggested that the OS of the high-PRM-score group was significantly lower than that of the low-PRM-score group in inner-training cohort (log-rank test, p < 0.001, Figure 5A), as well as the Kaplan–Meier curves of the inner-testing cohort (log-rank test, p < 0.001, Figure 5E). The PRM-score distribution and the expression of eight PR significant genes in the inner-training and inner-testing cohorts are presented in Figures 5B,C,F–G. Then, ROC curves were used to estimate the validity of the eight PR risk assessment tool-constructions in CC cohorts. The AUCs were equal to 0.738 at 3 years and 0.782at 5 years in the inner-training group (Figure 5D, Supplementary Figures 4A, B). Similarly, the AUCs were equal to 0.708 at 3 years and 0.753 at 5 years in the inner-testing group (Figure 5H, Supplementary Figure 4C), which showed that the model could achieve satisfactory predictive accuracy in both the inner-training and inner-testing cohorts. We established the survival analysis and ROC curves in the independent validation cohort (GSE17538), which showed the significant difference in OS between high- and low-PRM-score groups (log-rank test, p < 0.001, Supplementary Figure 5A). The AUCs were equal to 0.644 at 3 years and 0.684 at 5 years in the independent validation group (Supplementary Figures 5B, C). And we established the model to predict the prognosis based on PR genes with “random Survival Forest” (Supplementary Figures 6A–C). And the 10-fold cross validation was employed to prevent overfitting.

We estimated the immune microenvironment between the eight genes related high- and low-PRM-score groups. Only the immune-scores were significant different between two groups, and the level of tumor-purity was the same in the two groups (Mann–Whitney U test, *p < 0.05; **p < 0.01; ***p < 0.001; p ≥ 0.05, not significant) (Figure 6A). We found that the infiltration of activated immune cell in TME was abundant in low-PRM-score groups, including the M1 macrophages, NK cells, CD4 T cells (Mann–Whitney U test, *p < 0.05; **p < 0.01; ***p < 0.001; p ≥ 0.05, not significant) (Figure 6B). To further evaluate the association between the expression of the tumor immune microenvironment and these eight genes, we analyzed the corrections between the 22 types of immune cell infiltration profiles and these eight genes (Supplementary Figure S7). GSEA was used to analyze potential biological characteristics of the PRM-score groups in CC patients. As shown in Figures 6C,D, according to the Hallmark and KEGG collection defined by MSigDB, the genes in the high-PRM-score group were mainly enriched in angiogenesis, KRAS signaling, and epithelial mesenchymal transition. And the genes in the low-PRM-score groups were mainly enriched in cell cycle, P53 signaling pathway, T-cell receptor signaling pathway, and PI3K/AKT/MTOR signaling.

The univariate and multivariable Cox regression models were applied to the inner-training cohort to evaluate the predictors of OS. Univariate analyses indicated that age, stage-N, stage-M, PRM-scores, and PR subtypes were associated with OS in CC patients (p < 0.05 in all cases, Table 1). Next, the multivariate Cox analyses found that age, stage-N, stage-M, PRM-scores, and PR subtypes were independent risk factors for OS based on forward and backward elimination methods (Table 1).

Because stage-N, age, stage-M, PRM-scores, and PR subtypes were predictive for OS in multivariate analysis, these variables were further included in the nomogram, which was for predicting the 1, 3, and 5-years OS for CC patients (Figure 7A). The weighted total score, calculated from these factors, was applied to predict the 1, 3, and 5-years OS of CC patients.

Besides, the model showed good accuracy for predicting the OS, and internal validation was performed using the inner-training cohort with a C-index of 0.739. Furthermore, the decision curve analysis (DCA) results of the nomograms also confirmed their clinical applicability for predicting the OS, with superior performance compared with PRM-scores and PR subtypes (Figure 7B). Calibration curves for the probability of OS at 3 and 5 years indicated satisfactory consistency between actual observation and nomogram-predicted OS probabilities in CC cohort (Figure 7C, Supplementary Figure 4B).