A pyroptosis-related signature in colorectal cancer: exploring its prognostic value and immunological characteristics

Dataset collection and pre-processing

RNA sequencing (RNA-seq) data (count values) and simple nucleotide variation in patients with CRC and the relevant clinical characteristics were obtained from TCGA dataset (https://portal.gdc.cancer.gov). The transcriptomic data of healthy colon tissues were retrieved from the GTEx database (http://xena.ucsc.edu/). To enable comparison with the GTEx dataset, we normalized the count values from both datasets using log2 (X + 1) transformation. For the accuracy of the functional analysis, we transformed the count value to tpm (transcripts per million) in the subsequent analyses. We excluded the cases lacking the survival information and cases with vague survival information. The somatic mutation data of patients with CRC was used to analyze the mutation profiles using the “maftools” package. MSI data were acquired from the cBioPortal website (https://www.cbioportal.org/).

The GSE39582 dataset from GEO (http://www.ncbi.nlm.nih.gov/geo) was chose as the external validation cohort. The samples in the validation set have survival information and their tumor stages are similar to those in the training set (Table S1). The single‑cell transcriptomic datasets GSE14677, GSE13639 were acquired from the Tumor Immune Single-cell Hub (TISCH).

The Institutional Review Boards of Zhejiang Cancer Hospital granted approval (No. IRB-2021-291) for this study, which was conducted in accordance with the ethical guidelines outlined in the Declaration of Helsinki. A total of 35 patients’ tissues of rectal cancer and their clinical infromation were collected from Zhejiang Cancer Hospital. All participants signed informed consent according to the Institutional Review Boards of Zhejiang Cancer Hospital. All tissue were stored at −80 °C once dissected from patients until sequencing. Total RNA was extracted from the CRC tissues. After assessing the quality and integrity of the RN, the enriched mRNA was fragmented into small pieces and used as a template for cDNA synthesis with random primers. The second-strand cDNA was synthesized, followed by end repair, A-tailing, and adapter ligation. The ligated products were amplified by PCR to construct the cDNA libraries. The libraries were quantified using qPCR and sequenced on an Illumina HiSeq platform. The read counts for each gene were obtained using HTSeq (Illumina Corporation, Foster City, CA, USA).

Classification of patients based on pyroptosis patterns

Thirty-three PRGs were extracted from prior studies (Duan et al., 2016) are shown in Table S1. On the basis of PRGs, we implemented consensus clustering to identify distinct PR patterns using the k-means method. We performed it with the “ConsensuClusterPlus” package and set the parameter “repetition” as 1,000 for stability (Wilkerson & Hayes, 2010).

Functional enrichment analysis of the DEGs

We adopted the “DESeq” package to identify DEGs between different PR patterns, setting the significance criteria for differential expression at P < 0.05 and an absolute Log2 fold change (FC) >1. We applied the “clusterProfiler” package (Yu et al., 2012) for functional enrichment analysis of DEGs. We retrieved immune-related genes from the Gene Ontology database (http://geneontology.org/) and performed ssGSEA to assess differences between clusters in immune-related pathways. The hallmarker genesets were obtained from the Molecular Signatures Database (MSigDB) (https://www.gsea-msigdb.org/gsea/msigdb/index.jsp) for Gene Set Variation Analysis (GSVA).

Establishment of the pyroptosis-related prognostic model

To evaluate the survival significance of the DEGs, we utilized univariate and multivariate Cox regression analyses. Thirty-four genes (P < 0.05) were selected for the next analysis. To screen optimal genes, the ‘glmnet’ package was applied to carry out the LASSO Cox regression analysis (Duan et al., 2016). Then, we used the ten selected genes and their corresponding coefficients to construct a prognostic signature. The score based on the pyroptosis signature was calculated using the subsequent equation: PR score = 0.112 ×Expr_FCRL1 + 0.252×Expr_WNT16 + 0.130×Expr_GRIK2 + 0.070×Expr _ZMAT1 – 0.005×Expr_ZG16 + 0.103×Expr_DRD4 + 0.044×Expr_MAPK12 + 0.220×Expr_OR51B5 + 0.004×Expr_PRSS21 + 0.004×Expr_MAGEA3, where Expr_i is the gene expression level. Utilizing the median risk score as a threshold, we categorized patients into the high-score group and low-score group. The interaction network was mapped using Cytoscape (version 3.8.0) (Shannon et al., 2003).

Assessment of the immunological characteristics

We evaluated the immunological profile of CRC patients, including inhibitory immune checkpoints, cancer immunity cycle activity, and immune cells infiltration. The genes regulating the clinical response to ICBs are characterized from 20 solid cancers and established The Cancer Immunome Atlas (TCIA) (Auslander et al., 2018).

Comprising seven steps, the cancer immunity cycle reflects the intricate interplay between the processing of neoantigens and the prevention of autoimmunity (Chen & Mellman, 2013). To present TIICs precisely, we assessed the infiltration status by applying several independent algorithms comprehensively, including CIBERSORT, TIP, and xCell (Chen & Mellman, 2013; Aran, Hu & Butte, 2017). In addition, we performed single-cell analysis to explore the correlation of key genes and TME with TISCH (Sun et al., 2021).

Analysis of chemotherapy drugs and nomogram construction

DrugBank is a publicly available database that integrates bioinformatics and medicinal chemistry to provide information on drugs and their targets (Uhlén et al., 2015).

For applying the signature in clinical settings, we used the key genes of the signature to construct the nomogram using ‘rms’, ‘nomogramEx’, and ‘regplot’ packages.

The levels of key proteins in clinical specimens

The Human Protein Atlas (HPA, version: 18.1) (https://www.proteinatlas.org/) aims to provide a public platform for researchers to explore the gene expression landscape of 24,000 human proteins (Uhlen et al., 2017). We explored the the protein expression of key genes applying the HPA database.

Statistical analysis

Correlation coefficients were calculated using Pearson chi-square test.

We used the Kaplan-Meier method from the “survival” package to generate survival curves of which statistical significance was estimated with log-rank tests. The time-dependent ROC curves with AUC values were created using “pROC” and the “time ROTC” packages. The Wilcoxon test and Kruskal-Wallis test were utilized to examine the differences between two groups and multiple groups, respectively. All statistical analyzes were performed with R (R Core Team, 2023). The statistical threshold for significance was P < 0.05.

Expression of PRGs in CRC

The expression of 33 PRGs were compared between the tumor and normal tissues after pre-processing the data from TCGA and GTEx. The results demonstrated that the disparities in the expression of 27 PRGs were significant (Fig. 2A). To explore the mutation profiles of PRGs in CRC, we analyzed somatic mutation data using the VarScan2 Annotation. The mutation information of the PRGs showed that 114 patients had pyroptosis-related regulator mutations (Fig. 2B). The gene with the highest frequency of mutations was NLRP7, and the mutation rate of the five genes was up to 5%. C-T was the predominant single-nucleotide variation in CRC (Fig. 2C). The high proportion of differentially expressed PRGs and the mutation profile implied that pyroptosis played a vital role in CRC.

Identification of a classification pattern mediated by PRGs in CRC

We performed a consensus clustering analysis of CRC patients based on PRGs in TCGA cohort. The results showed that two different regulation patterns were identified, including 106 cases in cluster 1 (C1) and 318 cases in cluster 2 (C2) (Fig. 3A). To explore the biological differences between the two pyroptosis-related clusters, we performed a series of analyzes. First, the differential analysis was performed to identify 489 DEGs. The differential expression of DEGs was visualized using a volcano plot (Fig. 3B).

Enrichment analyses for GO and KEGG demonstrated that the DEGs were enriched in hydrogen peroxide metabolism, spliceosome snRNP complex, and channel activity (Fig. 3C). To profile the characteristics of the two clusters, we performed GSVA analysis using the hallmark genesets. The analysis revealed significant differences in various aspects, including metabolism, stress response, and immune between the two clusters. For instance, disparities were observed in pathways such as glucose metabolism, reactive oxygen species pathway, and the complement system (Fig. 3D and Fig. S1A). Considering the association of pyroptosis and immunity, we used GSEA with the annotation of the immune system (GO:0002376) to investigate the differences in the immunity of the two pyroptosis-related clusters. The enrichment of various pathways presented significant differences, including natural killer cell activation, T cell-mediated immunity, antigen processing and presentation, etc., (Fig. 3E and Fig. S1B).

Construction and verification of the pyroptosis-related prognostic signature

Using the differences in the two pyroptosis-related patterns, we built a prognostic model and validated its stability and veracity as follows. First, DEGs were screened using the univariate and multivariate Cox regression analyses for prognostic prediction. Next, we narrowed down the candidate genes by applying the LASSO Cox regression model with a minimum penalty parameter (λ) (Fig. 4A). Ultimately, we obtained 10 key genes with which we developed a pyroptosis-related prognostic model. With this signature, each patient could obtain a corresponding score named the pyroptosis-related (PR) score. We then performed multivariate Cox analysis on key genes to reveal their impact on prognosis. (Fig. S1). In Fig. 4B, we presented the interaction network between key genes,which was constructed based on gene co-expression patterns.

To assess the prognostic value of PR score, univariate and multivariate Cox regression analyses were performed between it and clinical characteristics. The hazard ratio (HR) of PR score was 2.49 (95% CI [1.958–3.16]; P < 0.001), indicating that PR score could serve as an independent prognostic factor in CRC (Fig. 4C).

With an increase in the PR score, the survival status of patients and expression of key genes exhibited significant differences. (Fig. 5A). We performed Kaplan–Meier and time-receiver operating characteristic curves (ROC) analyses to assess the discriminatory performance of the prognostic signature. The analysis suggested that patients in the high-score group had worse prognosis (P < 0.001) (Fig. 5B). The AUC for 1-, 2-, and 3-year in the ROC curve analysis was 0.72, 0.75, and 0.72, respectively, indicating excellent discriminatory ability of the prognostic signature (Fig. 5C). We validated the signature by applying it to GSE39582 dataset and performing Kaplan-Meier and time-ROC analyses. The results verified the robustness of the signature (Figs. 5D–5F).

Tumor immunity-associated characteristics

GSEA of the DEGs revealed the differences in the enrichment of multiple immune-related pathways. To explore the utility of the signature in the tumor immune microenvironment, we performed several analyses to evaluate the relationship between the PR score and the characteristics of tumor immunity.

We conducted Pearson correlation analysis to investigate the correlation of the PR score and the immune checkpoints. As presented in Fig. 6A, PR score was closely associated with various immune checkpoints, including PD1, PD-L1, and CTLA4, which is valuable in immunotherapy. Additionally, our study suggested that PR score may serve as an indicator of cancer immunogenicity in CRC, as demonstrated by the higher microsatellite instability (MSI) scores in the high-score group based on the results of Wilcoxon analysis (Fig. 6B). The above results indicated the role of PR score in immunotherapy. We combined the data of TCIA to analyze the correlation between PR scores and immune cells. The infiltration of multiple immune cells, including B cells, several subtypes of T cells, macrophage M2, neutrophils, and regulatory T cells (Treg) cells, was found to be significantly different between the two groups (Fig. 6C). Moreover, the abundance of B cells, CD8 T cells, dendritic cells, M2 macrophages, and neutrophils were strongly related to the PR score (Fig. S2). The immunity cycle sorts the process of antitumor immunity into seven steps, composed of 23 processes. Activities of multiple processes in the cycle were changed, including the priming and activation (Step 3) and infiltration of immune cells to tumors (Step 4) (CD8 T, CD4 T, macrophage, Th2, NK, and Treg cell recruitment) (Fig. S3). By further analyzing the correlation, we found that the PR score was related to the recruitment of various immune cells (Fig. 6D).

To assess the infiltration of TIICs with different algorithms for verification, we downloaded the data of immune cells from xCell based on the ssGSEA algorithm to investigate the variation as the PR score changes. The heatmap showed that the infiltration of multiple above-mentioned immune cells was associated with the PR score, which is consistent with the findings from TCIA. (Fig. 6E). The Sankey diagram suggested good coherence between the PR score and immune score, matrix score, and microenvironment score (Fig. 6F).

Single-cell analysis for the key genes of the signature

We screened two scRNA-seq atlases of CRC from TISCH database to explore the immunological characteristics of the signature. The 10 samples of GSE14677 are tumor tissues from patients without therapy. And the five samples of GSE13639 are tumor tissues from patients who accepted adoptive cell therapy. The percentage of CD4Tconv cells was the biggest in untreated patients and the number of CD8Tcm cells was predominated in patients after immunotherapy (Figs. 7A–7D). To analysis the association of the signature and immune cells, we analysed the expression of genes with positive coefficient in cell clusters. The intense expression of these genes in B cells indicated positive relevance between PR score and the infiltration of B cell (Fig. 7E). That was consistent with the results of TCIA and xCell (Figs. 6C and 6E). By contrast, the expression of these genes presented different distribution features in GSE13639. The level of these genes was highest in CD8Tcm cells (Fig. 7F). The difference should be caused by adoptive cell therapy and it indicated the interaction between immunotherapy and PR score.

Potential clinical application of the signature

Given the role of pyroptosis in oncotherapy, we extracted the target genes in CRC from the Drugbank database and studied the association between the PR score and the target genes. Figure 8A showed the differences in the expression of many genes targeted by drugs, such as oxaliplatin, irinotecan, bevacizumab, nivolumab, and ipilimumab, suggesting different responses for two groups to chemotherapy, molecular targeted therapy, and immunotherapy.

Neo-adjuvant therapy (nCRT) palyed an irreplaceable role for CRC patients, particularly rectal cancer. For the exploration of nCRT with the signature, we analysis the the response of nCRT and the PR score. The result indicated that patients with higher PR scores were more likely to be sensitive to nCRT, while those with low scores were more likely to be resistant (Fig. 8B). The further analysis of score and the response of nCRT revealed that score of patients resistant to nCRT was significantly lower than the counterparts (Fig. 8C, P < 0.05). The evaluation of immune cells with xCell was consistant with the training set, namely, patients belonged to the high-score group had better infiltration of immune cells (Figs. 8D–8F).

Because of the value of the PR score in predicting the prognosis, we generated a nomograph featuring 10 key genes for the clinical management of patients with CRC (Fig. 8G). It would help clinical doctors in making personalized treatment efficient.

The protein expression of key genes in human tissue

To validate the reliability of the prognostic signature at protein level, we investigated the expression of most key genes in human tissue samples with HPA database except for WNT16, DRD4 and MAGEA3, which were unavailable. As shown in Fig. 9, the protein expression of FCRL1, GRIK2, MAPK12, and OR51B5 showed differential level between normal colon tissue and tumor tissue (Fig. 9).