This article is mentioned in:

Introduction

Colorectal cancer (CRC) is one of the most common tumors worldwide, and its incidence and mortality rates are ranked third among all malignant tumors (1–3). Early-stage CRC is mainly treated with surgery (4), while advanced-stage or metastatic CRC is treated with chemotherapy combined with targeted therapy (5). Unfortunately, a number of patients with CRC have no obvious symptoms in the early stage of disease, leading to patients reaching advanced stage disease at the time of diagnosis, missing the best opportunity for surgery (6). Additionally, even after successful surgical resection, some patients still encounter local recurrence or distant metastasis (7). Chemotherapy and targeted therapy have improved the therapeutic effect to a certain extent, but there are also side effects and drug resistance issues (8). Therefore, it is of great significance to develop new biomarkers to detect early CRC, improve the cure rate of surgery and overcome the existing drug resistance mechanism to prolong the survival of patients. With the popularization of gene detection, CRC classification based on biomarker classification, such as RAS and BRAF mutation and microsatellite instability (MSI), contributes to treatment and prognosis prediction in CRC (9). Therefore, further exploration and improvement in biomarker classification is of great value to provide individualized and accurate treatment for patients with CRC.

Pyroptosis is a newly discovered type of programmed cell death (10). Contrary to apoptosis, pyroptosis is accompanied by changes in cell membrane permeability, water influx, cell rupture and the release of inflammatory factors (11), causing a strong inflammatory response. Pyroptosis can be induced via the classical inflammasome pathway (caspase-1-dependent) and non-classical inflammasome pathway (caspase-4/5/11-dependent) (12). Both pathways lead to the cleavage and activation of the gasdermin (GSD) family of proteins (GSDMA, GSDMB, GSDMA, GSDMD and GSDME) and finally result in membrane pore formation and cell death (13). A previous study has shown that pyroptosis is involved in the development of malignant tumors (14). In CRC, pyroptosis not only affects tumor angiogenesis (15,16) but also increases chemosensitivity (17). Therefore, pyroptosis-related genes (PRGs) may be biomarkers of CRC and improve treatment and prognosis.

In the present study, the genomic mutations and expression patterns of PRGs in CRC were analyzed using datasets from The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) databases. Distinct pyroptosis subgroups were identified and then the differential genes between the different subgroups were obtained to generate a PRG risk score. Further analysis was performed to determine whether the PRG risk score could assess prognosis, immune cell infiltration and drug sensitivity in CRC.

Materials and methods

Data sources

Somatic mutation data, mRNA expression data and the clinical information of patients with CRC were downloaded from TCGA (18) (https://portal.gdc.cancer.gov/repository) on February 20, 2022. There were 44 normal tissues and 568 CRC tissues in this dataset. The GSE39582 dataset (19) (platform GPL570; 579 patients with CRC) was downloaded from the GEO database (http://www.ncbi.nlm.nih.gov/geo/) on February 25, 2022. Some cases with missing clinical information were excluded from analysis. The CamBat function of the ‘SVA’ package in R (version 4.1.0; http://www.r-project.org) was used for merging and correcting the TCGA and GEO datasets.

Gene expression analysis and unsupervised clustering

A total of 52 PRGs were pooled from previous reviews (20–22) (Table SI). Gene differential analysis was performed using the ‘limma’ package in R (version 4.1.0), with P<0.05 as the cut-off. Unsupervised clustering analysis (23) was performed using the ‘ConsensusClusterPlus’ package (24) to identify pyroptosis-related subgroups. The repetition was set at 1,000 times to guarantee the stability of classification. Principal component analysis (PCA) was also conducted to investigate whether different subgroups had different characteristics.

Clinical samples, RNA extraction and reverse transcription-quantitative PCR (RT-qPCR)

The samples used for RT-qPCR were paired CRC tissues and adjacent normal tissues from 27 patients. These samples were collected from patients who underwent CRC surgery at Zhongshan People's Hospital (Zhongshan, China) but had not undergone chemotherapy, targeted therapy, immunotherapy or other internal medicine treatments. These tissue samples were collected from June 2022 to May 2023, and the age range of patients was 41–81 years old, with 14 females and 13 males. The inclusion criteria were as follows: i) All samples must have undergone pathological examination and be diagnosed with CRC; ii) fresh or flash frozen tissue samples were preferred to maintain tissue integrity and bioactivity; iii) each sample should have enough tissue to meet the experimental needs; iv) detailed clinical information available, including the patient's age, sex, tumor stage, grade and treatment history; and v) patients or their family members had signed the informed consent form for the use of their tissue samples in research. The exclusion criteria were as follows: i) Tissue samples with obvious infection, inflammation or severe necrosis; ii) samples containing a large number of non-tumor tissues (such as normal mucosa and adipose tissue) to reduce interference; and iii) samples from patients who had received radiotherapy or chemotherapy were excluded as these treatments may change the biological characteristics of tumors. The clinical information of the patients is listed in Table SII. The clinical staging system for patients with CRC used in the present study was the 8th edition of TNM staging system released by the American Joint Committee on cancer and the Union for international cancer control in 2017 (25).

Total RNA was extracted from patient tissue samples by TRIzol reagent (Invitrogen; Thermo Fisher Scientific, Inc.) according to the manufacturer's protocol. Subsequently, the extracted RNA was reverse transcribed using the PrimeScript RT reagent Kit with a gDNA Eraser (Takara Bio, Inc.), according to the manufacturer's instructions. The cDNAs were subjected to SYBR Green-based (Thermo Fisher Scientific, Inc.) qPCR analysis. qPCR reactions were carried out on an ABI 7500 Real-Time PCR System (Applied Biosystems; Thermo Fisher Scientific, Inc.) with the following thermocycling conditions: Initial denaturation at 95°C for 3 min; denaturation at 95°C for 10 sec, annealing at 60°C for 30 sec and extension at 72°C for 30 sec, which was repeated for 40 cycles; melting curve analysis at 95°C for 15 sec, 60°C for 1 min, 95°C for 15 sec. β-actin was used as the reference gene to normalize the cDNA input. The relative expression levels of target genes were quantified using the ΔΔCq method (26). The primers used in the qPCR are listed in (Table SIII).

Gene Set Variation Analysis (GSVA)

GSVA (27) was performed using the ‘GSVA’ packages in R (version 4.1.0) to analyze the differences in biological signaling pathways between the different subgroups, with adjusted P<0.05 as the cut-off. ‘c2.cp. Kegg.v7.4.symbols’ was downloaded from the MSigDB database (https://www.gsea-msigdb.org/gsea/msigdb) to conduct GSVA.

CIBERSORT analysis

Initially, the gene expression matrix was extracted from the raw data, which included the gene expression levels of all samples. Subsequently, the CIBERSORT algorithm in R (version 4.1.0) was employed to estimate the abundance of 23 types of immune cells in each sample. The inputs for the algorithm comprised the gene expression matrix and the cell type-specific gene sets. The output results provided the relative abundance of each immune cell type in every sample.

Estimation of immune cell infiltration in the tumor microenvironment (TME)

A single-sample gene set enrichment analysis (ssGSEA) was conducted to compare the immune cell infiltration profiles between the different subgroups. In total, 23 types of TME-infiltrating immune cells, such as regulatory T cells (Tregs), activated CD8+ T cells and macrophages, were obtained from the study by Charoentong et al (28). R (version 4.1.0; http://www.r-project.org) was used to perform ssGSEA analysis to calculate the abundance of each type of infiltrating cells in the CRC TME. The enrichment scores were used to indicate the degree of infiltration. ES >0.5 or ES <-0.5 were selected as the threshold for significant enrichment or depletion. Adjusted P<0.05 indicated that the enrichment score was statistically significant.

Identification of differentially expressed genes (DEGs) in the different subgroups

The ‘Venn diagram’ package in R (version 4.1.0) was used to identify the overlapping differential genes between the different subgroups. Then, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed to determine the molecular mechanisms of the overlapping differential genes using the ‘clusterProfiler’ package in R (version 4.1.0).

Kaplan-Meier analyses

Kaplan-Meier survival analyses were performed using R (version 4.1.0) to evaluate the association between gene expression levels and patient survival. The median expression level of the target gene was used as the cut-off to classify patients into the high-expression and low-expression groups. Patients with gene expression levels higher than or equal to the median were classified as the high-expression group, while those with expression levels below the median were classified as the low-expression group.

Identification of prognostic genes and the generation of a risk score model

A univariate Cox regression analysis was performed to identify genes that were significantly associated (P<0.01) with the survival of patients with CRC. TCGA and GEO data were integrated and randomly divided into the training and test cohorts. The risk score model was constructed using LASSO analysis and the ‘glmnet’ R package. The risk score was calculated as follows: Risk score= Σi n=1 vi × ci, where n indicates the quantity of genes, vi indicates the expression of genes and ci indicates the regression coefficient of genes i. According to the median risk score of the training cohort, both cohorts were divided into the high-risk and low-risk groups. Then, receiver operating characteristic (ROC) curved were used to assess the prognostic practicability of the risk score. Using the ‘regplot’ and ‘survival’ packages in R, a nomogram was constructed based on the clinical characteristics to predict the patient survival time (1-, 3- and 5-year survival). Calibration curves were produced using the ‘rms’ package in R to verify the nomogram.

Relationship between the risk score and drug sensitivity

Genomics of Drug Sensitivity in Cancer (GDSC) (https://www.cancerrxgene.org/) is a database providing the relationship between the antitumor drug response and genomic features. According to the list of drugs, the ‘PRRophetic’ package in R was used to construct the ridge regression model. By analyzing the expression profiles based on the drug IC50, the relationship between the risk score and the sensitivity to antitumor drugs was shown.

ESTIMATE analysis

All analyses were conducted in the R (version 4.1.0; R Core Team) environment. The main R packages used include estimate and ggplot2. The estimate package was utilized to calculate the ImmuneScore, StromalScore and ESTIMATEScore, while the ggplot2 package was employed for data visualization.

Statistical analysis

All statistical analyses were performed using R (version 4.1.0). The ‘maftools’ package was used to analyze the mutation frequency of genes. The ‘RCircos’ package was used to analyze the copy number variation landscape of PRGs in 23 pairs of chromosomes. In the qPCR experiment, the comparison between two groups was performed using a paired two-sample t-test, while the comparison between two groups in other analyses was conducted using an unpaired two-sample t-test. One-way ANOVA was used to compare three or more groups, followed by Bonferroni correction as the post hoc test. All P-values were two-sided, and P<0.05 was considered to indicate a statistically significant difference.

Results

Landscape of the genetic variation of the PRGs in CRC

A flow chart of the present study is presented in Fig. S1. Using the mutation data from TCGA, the somatic mutations and copy number variations (CNVs) of 52 PRGs in CRC were analyzed. Among 399 patients, 297 exhibited PRG mutations, with a frequency of 74.44% (Fig. 1A). Among them, TP53 (55%) had the highest mutation frequency, followed by NLRP7 (6%). In total, 8 PRGs, including CASP6, PRKACA, PYCARD, TNF, CHMP4A, CYCS, GSDME and PJVK, did not show any mutations in the CRC samples. The location of CNVs of PRGs on chromosomes is shown in Fig. 1B. The CNV frequency of PRGs was also analyzed (Fig. 1C), and it was found that the copy number of most genes of the GSD family, such as GSDMA, GSDMB, GSDMC and GSDMD, was amplified, while CASP9, CASP3, IRF2 and CHMP7 were mainly deleted. Next, the mRNA expression levels of 52 PRGs were investigated (Fig. 1D). Compared with normal tissues, most PRGs with increased CNV, such as GSDMA and GSDMC, showed significantly higher expression in CRC samples, and an opposite pattern was observed for genes with decreased CNV, such as CASP9 and IRF2. However, for certain PRGs, such as GSDMB and TP53, the alterations of CNV may not be the main factor affecting their expression, as these PRGs with increased CNV had significantly decreased expression in tumor tissues. These findings suggest an uncertain association between the gene mutational intensity and expression level.

Figure 1. Genetic mutational and expression landscape of PRGs in CRC. (A) The waterfall plot demonstrates the mutation profiles of 52 PRGs in the TCGA COAD-READ. (B) The location of CNVs of 52 PRGs on chromosomes in the TCGA-COAD. (C) CNV frequency of PRGs in the TCGA COAD-READ. (D) The mRNA expression levels of 52 PRGs between normal and CRC samples in the TCGA COAD-READ and GSE39582 dataset. **P<0.01, ***P<0.001. CRC, colorectal cancer; CNV, copy number variation; PRGs, pyroptosis-related genes; TCGA, The Cancer Genome Atlas; COAD, colonic adenocarcinoma; READ, rectal adenocarcinoma.

Generation of pyroptosis subgroups in CRC

The data from 1,089 patients with CRC in TCGA and GSE39582 were integrated to investigate the potential biological function of PRGs (Table SIV). Univariate Cox and Kaplan-Meier analyses were conducted to reveal the prognostic values of the PRGs (Table SV). Then, according to the expression level and prognostic value of the PRGs, a correlation network was constructed (Fig. 2A and Table SVI). There were 7 PRGs (GZMB, CYCS, CASP3, CASP1, CASP6, IRF1 and NLRP1) with statistical significance in the univariate and Kaplan-Meier analyses (all P<0.05; Fig. S2 and Table SVII). Generally, if the upregulation of gene expression leads to a worse prognosis, the gene is termed a risk factor. Conversely, if the upregulation of a gene expression leads to an improved prognosis, the gene is termed a favorable factor. NLRP1 was the only risk factor, and increased expression of NLRP1 worsened the prognosis of CRC. By contrast, increased expression of GZMB, CYCS, CASP3, CASP1, CASP6 and IRF1 predicted a more favorable prognosis. The mRNA expression levels of the 7 genes were also verified by RT-qPCR using clinical samples obtained from Zhongshan People's Hospital (paired CRC tissues and adjacent normal tissues; n=30) (Fig. 2B). It was observed that, compared with the adjacent normal tissue, GZMB expression was elevated, while the expression of the other 6 genes (NLRP1, CYCS, CASP3, CASP1, CASP6 and IRF1) was decreased. Based on the expression patterns of the 7 PRGs, a consensus clustering and PCA analysis were performed. Notably, patients were clustered into two subtypes (clusters A and B; Fig. 2C and D). Furthermore, there were statistically significant differences in the prognosis of patients between clusters A and B (Fig. 2E). The differences in clinical features between the two groups of patients were also analyzed (Fig. 2F). It was found that there were no significant differences in the clinical characteristics between the two groups of patients, but the PRGs in cluster B were more inclined towards high expression, especially GZMB, IRF1, CASP1, CASP3 and CASP6.

Figure 2. Generation of PRG subgroups in CRC. (A) The correlation network of PRG interactions in The Cancer Genome Atlas cohort. (B) Analysis of the mRNA relative expression of the 7 genes by reverse transcription quantitative PCR (n=30). (C) The consensus matrix heatmap of consensus clustering analysis (k=2). (D) PC analysis for the expression of PRGs indicating a notable difference in clusters A and B. (E) Kaplan-Meier analysis of clusters A and B. (F) The expression of the 7 PRGs and the clinicopathological variables of clusters A and B. *P<0.05, **P<0.01, ***P<0.001. PCA, principle component; PRGs, pyroptosis-related genes.

Difference in the TME between the PRG subgroups

A GSVA was conducted to investigate the differences in the biological signaling pathways of the PRG subgroups. Notably, cluster B was enriched in immune-associated pathways, such as ‘NOD-LIKE RECEPTOR SIGNALING PATHWAY’, ‘T CELL RECEPTOR SIGNALING PATHWAY’, ‘CHEMOKINE SIGNALING PATHWAY’ and ‘NATURAL KILLER CELL-MEDIATED CYTOTOXICITY’ (Fig. 3A). Therefore, we speculated that cluster B may have a more favorable prognosis than cluster A due to the rich immune cell infiltration in its TME. The CIBERSORT algorithm was used to estimate the infiltration levels of 23 types of immune cells and verify our hypothesis. As expected the abundance of almost all immune cells was significantly increased in cluster B (Fig. 3B).

Figure 3. Tumor microenvironment differences in the PRG subgroups. (A) The heatmap of Gene Set Variation Analysis result, showing the active biological pathways of clusters A and B. (B) The infiltrating levels of 23 immune cells in different subpopulations. *P<0.05, **P<0.01, ***P<0.001. PRG, pyroptosis-related gene; TCGA, The Cancer Genome Atlas.

Identification of DEGs

Considering that phenotypic changes are closely related to differential expression of genes, the expression levels of all genes were analyzed using the ‘limma’ package and 76 DEGs between clusters A and B were obtained (Table SVIII). Next, functional enrichment analysis was conducted to explore the underlying biological function of these DEGs. Both GO and KEGG analyses showed that DEGs were mainly enriched in immunomodulation-associated pathways (Fig. 4A and B). Then, univariate Cox analysis was performed and 9 DEGs (P<0.01) with significant effects on survival were identified (Table SIX). A consensus clustering analysis was performed and patients were divided into three gene clusters (clusters A-C) to investigate the specific prognostic value of the 9 DEGs (Fig. S3). It was found that patients in cluster A had the shortest OS (Fig. 4C). The expression of most PRGs was significantly different in the three gene clusters (Fig. 4D). The relationship between the three gene clusters and clinicopathologic features is shown in Fig. 4E.

Figure 4. Screening out DEGs and revealing their regulatory functions. (A) Gene Ontology and (B) Kyoto Encyclopedia of Genes and Genomes analyses demonstrated the abundance of immunomodulation-associated pathways. (C) Kaplan-Meier analysis of gene clusters A, B and C in the TCGA COAD-READ and GSE39582 cohort. (D) PRGs have significant expression differences in the three different clusters. (E) The relationship between the three gene clusters and clinicopathological features. *P<0.05, **P<0.01, ***P<0.001. PRG, pyroptosis-related gene; TCGA, The Cancer Genome Atlas; COAD, colonic adenocarcinoma, READ, rectal adenocarcinoma.

Construction and validation of the PRG score

LASSO and multivariate Cox analysis were conducted for the 9 pyroptosis cluster-associated prognostic DEGs to construct an optimal PRG score. In total, 5 DEGs, including GZMB, CASP1, LINC00261, MMP3 and CKMT2, were used (Fig. S4). These 5 genes were not only statistically significant, but also contributed the most in the model, which may more accurately reflect the survival risk of the patients. The reason why 9 genes were not used is that reducing the number of genes can reduce the complexity of the model, avoid overfitting and making the model more powerful in predicting new data. Additionally, 5 genes may mean that the model is easier to understand and apply, which is convenient for practical operation in clinical practice. The formula used for constructing the PRGs score was as follows: Risk score=(expression of GZMB × −0.0990) + (expression of CASP1 × −0.1354) + (expression of LINC00261 × −0.2266) + (expression of MMP3 × −0.0849) + (expression of CKMT2 × −0.1529). The distribution of patients with CRC in the two pyroptosis clusters, three gene clusters and two PRG score groups is shown in Fig. 5A. It was found that there were significant differences in the PRG scores between the different PRG clusters and gene clusters. Consistent with the results of the survival analysis, the PRG risk scores of PRG cluster A and gene cluster A were the highest (Fig. 5B and C), suggesting that patients in these clusters had a higher risk of death (Fig. 5D). Furthermore, Kaplan-Meier analysis of the training cohort also confirmed that patients with high PRG risk scores had a worse prognosis (Fig. 5D), and the area under the curve (AUCs) for 1-, 3- and 5-year OS were 0.680, 0.721, and 0.698, respectively (Fig. 5E), which indicates that the model has good prediction performance. A risk plot of the PRG score was also generated. With an increase in PRG risk score, the OS of patients with CRC decreased and mortality increased (Fig. 5F). The same method was used to calculate the PRG risk score of the entire cohort (Fig. S5A-C) and the test cohort (Fig. S5D and E) to verify the robustness of the PRG risk score. Based on the median score of the training cohort, patients in the test cohort and the entire cohort were assigned to the high-risk and low-risk subgroups. Similar results were obtained in the training group, with good AUC values (Fig. S5B and E). These findings indicate that the PRG risk score has a great prognostic value for patients with CRC.

Figure 5. Generation of the PRG risk score and further validation. (A) The relationship between the two pyroptosis clusters, three gene clusters and the two PRG risk score groups. (B) PRG risk scores of PRG clusters A and B, as well as (C) gene clusters A, B and C. (D) Kaplan-Meier analysis of the training cohort. (E) The receiver operating characteristic curves showing the AUCs of 1–3-, and 5-year OS. (F) Risk plot of the PRG risk score. (G) A nomogram for predicting the prognosis of patients with colorectal cancer. (H) The calibration curve of the nomogram. (I) Decision Curve Analyses to evaluate the clinical value of the nomogram model. AUC, area under the curve; OS, overall survival; PRG, pyroptosis-related gene; T, tumor; N, node; M, metastasis.

Constructing a nomogram for patients with CRC

As the PRG risk score was closely related to the OS of patients with CRC, by integrating the PRG risk score and clinical parameters, a nomogram was constructed to predict the 1-, 3- and 5-year OS of patients with CRC (Fig. 5G). The calibration curve of the nomogram is shown in Fig. 5H, implying great accuracy between actual observations and predicted values. Decision curve analysis indicated that the nomogram curve was higher than the other curves, suggesting that within the high-risk threshold range (0 to 1), the prediction results of the nomogram model may better guide clinical decision-making (Fig. 5I).

Relationship between PRG risk score and the TME

The most important components of the TME are stromal cells and immune cells, and the immune score calculated using the ESTIMATE algorithm is an important index to assess the TME (29,30). The correlation between the PRG risk score and the immune score was analyzed. It was noted that the PRG risk score was inversely associated with the immune score (Fig. 6A), suggesting that the PRG risk score could be used to evaluate the abundance of immune cells in the TME of CRC. In addition, it was demonstrated that the PRG risk score was positively correlated with the abundance of M0 macrophages Tregs, M2 macrophages, memory B cells and follicular helper T cells, whereas it was negatively correlated with the abundance of naive B cells, activated dendritic cells, resting natural killer cells, plasma cells, resting memory CD4+ T cells and activated memory CD4+ T cells (Fig. 6B). This suggests that CRC tumors with high PRG risk scores may have a tumor immune microenvironment more prone to immune evasion. The correlation between GZMB, CASP1, LINC00261, MMP3 and CKMT2 and immune cell infiltration was also explored. Most immune cells were closely related to the 5 genes (Fig. 6C). Among these 5 genes, the expression levels of CASP1, GZMB and MMP3 were significantly correlated with almost all immune cells, suggesting that they are most critical in CRC immune cell infiltration. In particular, GZMB showed a strong positive correlation with the activated memory CD4+ T cells.

Figure 6. Assessment of the TME in distinct groups. (A) The relationship between the PRG risk score and the immune score. The immune score was calculated with the ESTIMATE algorithm. (B) The correlation of the PRG risk score with the immune cell types. (C) The correlation between the 5 genes in the PRG risk score construction and the infiltration of immune cells. *P<0.05, **P<0.01, ***P<0.001. PRG, pyroptosis-related gene; TME, tumor microenvironment; NK, natural killer.

Relationship between the PRG risk score and MSI

A number of studies have demonstrated that MSI is associated with the efficacy of immune checkpoint inhibitors. Patients with MSI-high can benefit from immune checkpoint inhibitors (31–34). The results suggest that patients with high PRG risk scores are more likely to have microsatellite stability (Fig. 7A and B). These findings indicate that patients with low PRG risk scores are more likely to benefit from immune checkpoint inhibitors. Consistent with the results of the MSI analysis, correlation analysis indicated that the PRG risk score is negatively associated with tumor mutation burden (Fig. 7C), although the P-value was not statistically significant. The somatic mutation distribution of different patterns of PRG risk scores in TCGA dataset were also investigated. The results showed that APC, TP53, KRAS and PIK3CA mutations were more common in the high PRG risk score subgroup (Fig. 7D and E). Previous studies (35–37) have shown that these 4 genes are the key genes in the development of CRC.

Figure 7. Relationship between the PRG risk score and MSI/TMB. (A and B) Relationships between the MSI pattern and the PRG risk score in the TCGA COAD-READ. (C) Relationships between the TMB pattern and the PRG risk score in the TCGA-COAD. (D and E) Waterfall plots showing somatic mutation of the high and low PRG risk score subgroups in the TCGA COAD-READ. PRG, pyroptosis-related gene; TMB, tumor mutational burden; MSI, microsatellite instability; -L, low; -H, high; MSS, microsatellite stability; TCGA, The Cancer Genome Atlas; COAD, colonic adenocarcinoma; READ, rectal adenocarcinoma.

Drug sensitivity analysis

The IC50 values of 138 antitumor drugs in TCGA database were measured to determine whether the PRG risk score could predict the response to antitumor drugs. It was found that patients with low PRG risk scores were sensitive to cisplatin, paclitaxel, gemcitabine, sorafenib, camptothecin and Epothilone.B. (Fig. 8A-F), while patients with high PRG risk scores were sensitive to imatinib, dasatinib, shikonin and CGP.60474 (VEGFR-2 inhibitor) (Fig. 8G-J).

Figure 8. Drug sensitivity analysis. (A-F) Sensitive drugs for patients with low PRG risk score CRC. (G-J) Sensitive drugs for patients with high PRG risk score CRC.CRC, colorectal cancer; PRG, pyroptosis-related gene.

Discussion

Pyroptosis is an important type of cancer cell death (31). Different from apoptosis and autophagy, pyroptosis is often followed by strong immune and inflammatory responses, implying that PRGs may markedly modulate the tumor immune microenvironment (6,23). A previous study reported that pyroptosis has a complex effect on cancer, which varies in different tissues and genetic backgrounds (38). In liver cancer, 17 β-estradiol was found to have antitumor effects by activating the NLRP3 inflammasome and pyroptosis (39). In gastric cancer, a previous study showed that PRGs can regulate tumor-related signaling pathways and modulate the TME (40). Tang et al (16) reported that pyroptosis plays a pivotal role in tumor cell growth and metastasis in CRC. Miguchi et al (40) found that TGFBR2 can upregulate GSDME expression, which contributes to tumor cell proliferation and tumorigenesis. Therefore, accumulative evidence has demonstrated that the expression levels of PRGs significantly affect cancer progression. However, the biological function of most PRGs in CRC is still unknown. Therefore, it is necessary to clarify the mutation/expression profiles and functional characteristics of PRGs in CRC.

In the present study, the transcriptional alterations and expression patterns of 52 PRGs in CRC were examined using data from TCGA and GEO datasets. Although a certain correlation between the mutational intensity and expression level was not found, most PRGs were abnormally expressed in patients with CRC and GZMB, CYCS, CASP3, CASP1, CASP6, IRF1 and NLRP1 were related to prognosis. Using the unsupervised clustering method, patients with CRC were divided into two clusters (clusters A and B). A significant difference was observed in the clinical outcome, immune cell infiltration and cell signaling pathways between the two clusters. Based on the DEGs of PRG clusters A and B, three gene clusters (gene clusters A, B, and C) with different clinical features were obtained. Furthermore, the PRG risk score was generated to differentiate between pyroptosis subgroups. Both cluster A and gene cluster A, with the highest PRG risk score, had the poorest clinical outcome, suggesting that a high PRG risk score may indicate a poor prognosis in CRC. The findings of the present study also confirmed that the PRG risk score was closely linked to the clinicopathological features of CRC. The predictive value of the PRG risk score was validated by ROC for 1-, 3- and 5-year OS. A nomogram to estimate the 1-, 3- and 5-year OS of patients with CRC was also established by integrating the PRG risk score and clinical parameters. The calibration curve showed that this nomogram had great accuracy.

Similar to the findings of the present study, a recent study has indicated that pyroptosis can lead to cell rupture, proinflammatory cytokine release and immune cell infiltration (41). In the present study, cluster B was a population with upregulated PRGs and a more favorable survival. According to the results of the GSVA, cluster B was enriched in signaling pathways involved in inflammation and immune response, such as natural killer cell-mediated cytotoxicity, nod-like receptor signaling pathway, T cell receptor signaling pathway and chemokine signaling pathway. Consistently, the abundance of almost all immune cells was significantly increased in cluster B, indicating that upregulation of PRGs can activate inflammatory signaling pathways and promote immune cell infiltration. GO and KEGG analyses of gene cluster-related DEGs reported similar results. The relationship between the PRG risk score and the immune score was also analyzed to elucidate the correlation between pyroptosis and the immune response and it was observed that patients with CRC and low PRG risk scores had higher immune scores. Since genes related to PRG risk score were all favorable factors, a low PRG risk score predicted high PRG expression. The correlation between the PRG risk score and the tissue infiltration of each type of immune cell was also analyzed. It was found that a high PRG risk score was positively associated with the abundance of follicular helper T cells, M0 macrophages, M2 macrophages, memory B cells and Tregs. M2 macrophages are a type of anti-inflammatory macrophage that primarily inhibit immune responses by secreting anti-inflammatory factors. In the TME, M2 macrophages can promote tumor growth and angiogenesis while suppressing antitumor immune responses (42). Tregs are a type of T cell with immunosuppressive functions, primarily maintaining immune balance by inhibiting the function of other immune cells. In the TME, Tregs can suppress antitumor immune responses, thereby promoting tumor growth (43). A previous study has reported that a high abundance of Tregs can inhibit the antitumor immune response, leading to a worse prognosis (44). These findings are consistent with the findings of the present study and explain why high PRG risk scores are associated with poor clinical outcomes. Pyroptosis is closely correlated with immune cell infiltration and the PRG risk score may predict the tumor immune microenvironment in CRC.

Antitumor drug resistance is a major cause of progression and mortality in patients with CRC (45,46). Current antineoplastic drugs have limited efficacy; therefore, identifying patients who are sensitive to antitumor drugs can improve the efficacy of treatment and reduce resistance. In the present study, by integrating the PRG risk score with drug-related data from the GDSC database, it was found that patients with low PRG risk scores had an improved response to cisplatin, paclitaxel, gemcitabine, sorafenib, camptothecin and Epothilone.B. By contrast, patients with a high PRG risk score may have an improved response to imatinib, dasatinib, shikonin and CGP.60474 (a VEGFR-2 inhibitor).

There have been studies that have reported pyroptosis-related prognostic models for CRC; however, only one model also included the MMP3 gene, while the other genes were different or similar. These previous studies generally have limitations such as a small sample size, few included genes and no experimental validation (47–49). The present study solves the aforementioned limitations and provides a more accurate prognostic model for CRC treatment. However, the present study had some limitations. First, the present study employed retrospective data from public databases and the accuracy of the results needs to be verified by prospective studies. Second, further experiments are needed to explore the relationship between the PRG risk score and immune cell infiltration in CRC. Third, application of the PRG risk score has certain limitations. CRC has a high degree of heterogeneity and there may be significant differences in gene expression patterns between different individuals, which may lead to large differences in the performance of the model in different patients. The changes in gene expression levels over time are very dynamic, and a single detection may not capture these dynamic changes, which can affect the accuracy of model predictions. Different genetic testing platforms and technologies may produce different results. For instance, there may be differences in the results of RNA-sequencing and microarray chips, which can affect the consistency and reproducibility of the model. Additionally, the relatively high cost of genetic testing may limit its widespread application in resource-limited areas. In terms of potential clinical application value, the PRG risk score can help doctors to identify high-risk patients and develop personalized treatment plans. Based on the correlation between the PRG risk score and the tumor immune microenvironment, the state of the tumor immune microenvironment could also be preliminarily determined by the PRG risk score to guide the decision-making in immunotherapy. In addition, the PRG risk score can also be used for disease monitoring, regularly monitoring the changes in the expression of pyroptosis genes, so as to understand the progress of the disease in a timely manner and provide a basis for adjusting the treatment plan.

In conclusion, the mutation and expression characteristics of PRGs in CRC were analyzed and a prognostic PRG signature was constructed. This signature may help estimate the immune cell infiltration and therapeutic response in CRC. Thus, this signature may advance the treatment and prognosis evaluation of CRC.

Supplementary Material

Supporting Data

Supporting Data

Supporting Data

Supporting Data

Supporting Data

Supporting Data

Supporting Data

Supporting Data

Supporting Data

Supporting Data

Acknowledgements

Not applicable.

Funding

This study was funded by the Science and Technology Projects Funding of Zhongshan Science and Technology Bureau (grant no. 2023B1009).

Availability of data and materials

The data generated in the present study may be requested from the corresponding author.

Authors' contributions

CH and WD contributed to conceptualization, data analysis, funding acquisition and writing the original draft; YL and YQ contributed to data analysis; SZ conducted the RT-qPCR experiment; XJ contributed to conception, design, methodology, and data curation which involved management activities to annotate (produce metadata), scrub data and maintain research data (including software code, where it is necessary for interpreting the data itself; JX contributed to validation, supervision, conception and design. CH and WD confirm the authenticity of all the raw data. All authors read and approved the final version of the manuscript.

Ethics approval and consent to participate

The samples used for RT-qPCR were collected from patients who underwent colorectal cancer surgery at Zhongshan People's Hospital (Zhongshan, China). The studies involving human participants were reviewed and approved by the Committee of the Zhongshan People's Hospital. The participants provided their written informed consent to participate in this study.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Glossary

Abbreviations

Abbreviations:

References