1. Establishment of cox regression model for colorectal cancer (CRC)
Using cox regression model, we can predict the risk of colorectal cancer (CRC) including colon tumor (TCGA_COAD) and rectal tumor (TCGA_READ). Transcriptome data from COAD and READ datasets were randomly classified into train group (70% samples, total 413 cases) and test group (30% samples, total 177 cases). By using multivariate linear regression analyses and stepwise regression analysis, 13 pyroptosis-related genes (PRGs) were selected and a harzard modeling was generated on the basis of these genes. The formula of the modeling for colon and rectal cancer is: risk_score = 0.116011 + 0.004995*AIM2-0.000642*CASP1 + 0.002337*CASP5 + 0.000557*CASP6-0.005383*CASP8-0.007034*CASP9-0.040302*ELANE + 0.000154*GPX4-0.000887*GSDMD-0.116512*NLRP7-0.004866*NOD2 + 0.050655*PJVK + 0.002196*PRKACA. We obtained 2 groups (high-hazard group and low-hazard group) by using the median hazard score. In training and test set, the group with high hazard score had poor survival prognosis than the group with low hazard score (Fig. 1A, B). AUC prediction value of ROC curves in 1-year, 2-year as well as 3-year in training set were 0.758,0.791 as well as 0.802, respectively, whereas AUC prediction value of ROC curves in 1-year,2-year as well as 3-year in test set were 0.630, 0.755 as well as 0.755, respectively. All the above data demonstrated that this harzard modeling on the basis of pyroptosis-related genes has great prediction capacity for prognosis of CRC patients.
2. Mutant profiles of PRGs in CRC
We overviewed the whole picture of somatic mutants as well as copy number mutants of 35 PRGs in CRC. As demonstrated in Fig. 2A, 147 of 526 (27.95%) CRC samples showed genetic mutants. Missense mutation was the most frequent mutant category (Fig. 2A). The results showed that NLRP7 has the top mutation frequencies, followed by NLRP3 and SCAF11, among 35 PRGs (Fig. 2A). Furthermore, wilcoxon rank-sum test was employed to evaluate effects of genetic mutations of PRGs on expression profile of PRGs in CRC. When compared to mutant samples, we found expression levels in 3 of 35 PRGs of unmutated samples were significantly differentially regulated (Fig. 2B, CASP3, NLRC4, PLCG1). Finally, we compared genetic mutations of 35 PRGs between groups with 2 levels of hazard score. It showed a difference in the mutant category for 35 PRGs between these 2 groups (Fig. 2C).
3. Pathways involved in the cox regression model for CRC
GSEA was applied to explore the pathways involved in the cox regression modeling for of CRC. In Fig. 3A-I, low risk scores were significantly related to "IL6_JAK_STAT3_SIGNALING","INFLAMMATORY_RESPONSE","ALLOGRAFT_REJECTION","TNFA_SIGNALING_VIA_NFKB","IL2_STAT5_SIGNALING","INTERFERON_GAMMA_RESPONSE","HEDGEHOG_SIGNALING","APICAL_JUNCTION","COMPLEMENT".
4. Determination of differentially expressed genes (DEGs) within pyroptosis subgroups
In order to clarify gene expression pattern of pyroptosis subgroups, we used the R package “Limma” to conduct analysis of differential expression between groups with 2 levels of hazard score. Based on the cut-off threshold (adj.p.value < 0.05 as well as |logFC|>1, we identified the entire 133 DEGs, including 2 upregulated and 131 downregulated (Fig. 4A). The differentially expressed DEGs were mostly significantly associated with “antigen binding”, “complement activation, classical pathway”, and “complement activation” among Go and Pathway terms (Fig. 4B).
5. Association between the pyroptosis-associated genes and inflammation- associated genes in CRC
We then evaluated whether pyroptosis-associated genes are correlated with inflammation-associated genes in CRC. We selected gene set variation analysis (GSVA) to caculate the spearman correlation between 100 inflammation-related genes and 35 pyroptosis-related genes. To our surprise, no relevance has been found between these 2 gene sets. However, when we used gene expression profile to analyze the spearman correlation between 100 inflammation-related genes and 35 pyroptosis-related genes, we found certain relationships existed. For instance, PRGs (GSDME, NLRC4, NLRP1, etc.) were shown to be significantly associated with several inflammation-related genes in the heatmap (Fig. 5A). To further characterize interaction relationships, we selected highly correlated pairs (correlation coefficients were more than 0.65) to set up protein-protein interaction (PPI) networks by appling STRING database (https://string-db.org/). Both positive activation interaction and negative inhibition interaction revealed relationships among these inflammation-related genes and pyroptosis-related genes (Fig. 5B).
6. Immune features evaluations of signature of pyroptosis-related genes in CRC
Then we assessed the infiltrates status of immune cells and found that plasma cell, regulatory T cell, resting NK cell and neutrophil were present in the low-hazard group, whereas the high-hazard group was enriched with resting dendritic cell(Figs. 6A). In the current results, the ESTIMATE score also displayed that the immune score of the high-hazard group was lower than that of the low-hazard group (Fig. 6B). In addition, stromal score, and ESTIMATE score were considerably lower in the high-hazard group than in the low-hazard group, while the tumor purity was conversely higher (Fig. 6C, D, E). We next used ssGSEA algorithm to detect infiltrating modes of immune cells within two divergent sections. Except for type 2 T helper cell, effector memory CD4 T cell, CD56brght natural killer cell, as well as memory B cell, we found the patients in the low-hazard group had considerably higher levels of immune-cell infiltration by other 24 kinds of immune-cells compared with the patients in the high-hazard group (Fig. 6F). Finally, we used gene expression profile to analyze the spearman correlation between 100 immune-related genes and 35 pyroptosis-related genes. We found certain relationships existed, in which we labeled highly correlated pairs (correlation coefficients were more than 0.6) with star. For instance, PRGs (GSDME, NLRC4, GZMA, etc.) were shown to be significantly associated with several immune-related genes in the heatmap (Fig. 6G). Our study showed that pyroptosis risk scores were associated with immune features and that increased immune activity in the low-hazard group may be involved in the anti-tumor immunity of CRC.
7. Pyroptosis risk score in the function of platinum chemotherapy in CRC
Platinum-based chemotherapy played an important role in CRC therapy. We asked whether pyroptosis factors are involved in the reaction of CRC patients to chemotherapy with platinum agents We screened TCGA clinical data to identify sick persons with/without platinum-containing chemotherapy. Among these patients, we caculated pyroptosis risk score respectively. Firstly, in the 129 patients with platinum-based chemotherapy, sick persons with a low pyroptosis points showed a trend of higher survival rate than did sick persons with a high pyroptosis points, all of whom were subjected to chemotherapy (p = 0.38; Fig. 7A). Secondly, in the 312 sick persons without platinum-based chemotherapy, sick persons with a high pyroptosis points showed significantly lower survival rate than did patients with a low pyroptosis points, all of whom were not subjected to chemotherapy (p < 0.0001; Fig. 7B). The global survival rate of sick persons with platinum-containing chemotherapy was higher than those without platinum-based chemotherapy. This result suggested that pyroptosis risk score did not relate with the responses to platinum chemotherapy in CRC.