Pyroptosis-related genes are closely related to the tumor microenvironment
We explored the relationship between 50 PRGs and the immune microenvironment in the IMvigor210 cohort. First of all, we used ssGSEA to calculate the enrichment scores of 16 immune cells and 13 critical immune pathways in the transcriptome data of 348 mUC patients. The relationship heat map showed that most genes showed immune activation states, such as GZMB, GZMA, TNF, NOD2, NLRP6, NLRP3, NLRC4, IRF1, IL16, CASP5, CASP4, CASP1, and AIM2; a small number of genes showed immunosuppressive state, such as TP63 and CHMP4C (Figure 1A). GZMA and GZMB were highly positively correlated with CD8+ T cell and cytolytic activity scores, and cor > 0.8. There are four kinds of tumor therapeutic effects: PD (progressive disease), SD (stable disease), PR (partial response), and CR (complete response). PD/SD is the ICB non-response group, PR/CR is the ICB response group. We divided the patients into two groups according to PD/SD and PR/CR and analyzed the differences. The results of thermography showed that the expression of IRF1, GZMA, GZMB, HMGB1, CYCS, CASP6, and CHMP4B increased in the PR/CR group, while the expression of GSDMC, NLRP1, NLRP3, and GPX4 decreased in the PR/CR group, which are potential markers of ICB therapy (Figure 1B). To explore the interaction network of PRG in mUC, we calculated the correlation between 50 PRG pairs and constructed a co-expression network with Cytoscape. The hub gene was found out by the MCODE algorithm, which is the pink part of the network map (Figure 1C). These genes are highly correlated with the immune-related gene set of ssGSEA. We intersected this hub and ICB differential genes to obtain GZMA, GZMB, IRF1, and NLRP3 and observed their expression differences in three immunophenotypes (desert, excluded, and inflamed). It was found that GZMB, IRF1, and NLRP3 were all low expressed in desert typing, while high expression was found in inflamed typing (Figure 1D).
Construction of pyroptosis related gene prognostic index（PRGPI）
To explore the prognostic ability of 50 genes in immunotherapy, we performed KM analysis in IMvigor210 and GSE176307 cohorts. IMvigor210 has 19 genes that can divide 348 mUC patients into high and low-risk groups (Supplementary figure 1), while GSE176307 has 14 genes that can classify 88 mUC patients into risk stratification (Supplementary figure 2). We divided the genes with prognostic ability into HR < 1 and HR > 1 for display in Figure 2A. It was found that GSDMD, GZMB, IRF1, NLRP2, NLRP7, and GZMA showed HR < 1 in both cohorts, while TP63 showed HR > 1 in both cohorts. We use the LASSO algorithm to construct the prognostic model in the IMvigor210 cohort. The prognostic model constructed by this method has the characteristics of a few variables as possible and robust stability. With the increase of lambda value, the coefficient of the gene becomes 0, which indicates that the gene has little effect on the model and can be eliminated (Figure 2B). The results of 10X cross-validation showed that when the lambda value was 3, the partial likelihood deviance of the gene was the smallest (Figure 2C). Finally, the formula of the model is as follows:
PRGPI = GZMB * (-0.00185) + IRF1* (-0.00403) + TP63 * (0.00067)
Through this formula, we can calculate the PRGPI of each patient. The best cut-off value was found by using the "res. cut" function. According to the cut-off value, the patients were divided into two groups for KM analysis, and the results were statistically significant (p < 0.001). The prognosis of the high-PRGPI group was significantly worse than that of the low-PRGPI group (Figure 2D). ROC curve analysis showed that the AUC values of 0.5 years, 1 year, and 1.5 years were 0.562, 0.604, and 0.597. AUC values greater than 0.5 indicate that PRGPI has a certain prognostic ability. Figure 2F is the heat map of GZMB, IRF1, TP63 expression in high-PRGPI and low-PRGPI groups. Figure 2G is the heat map of the distribution of PRGPI in patients. We used univariate cox to analyze the prognostic ability of stage, gender and PRGPI, and it was found that only PRGPI had statistical significance. PRGPI had the ability of risk stratification in TCGA_stageⅠ, TCGA_stageⅢ, and male of mUC. We then explored the expression of GZMB, IRF1, and TP63 in the tumor microenvironment using bladder cancer single-cell data set GSE145281. The results showed that GZMB was significantly overexpressed in NK cells and CD8 Teff cells (Supplementary figure 3A), IRF1 was expressed in all immune cells (Supplementary figure 3B). At the same time, TP63 was not detected in all immune cells (Supplementary figure 3C).
Analysis of the relationship between PRGPI and clinical features
PRGPI is closely related to clinical factors and can guide the treatment of ICB. The Lund classification of bladder cancer can be divided into "Genomically unstable", "Basal/SCC-like", "UroA", "Infiltrated" and "UroB", and the immunophenotype is "inflamed", "excluded" and "desert". We drew the Sankey Diagram of PRGPI grouping, Lund typing, immunophenotyping, and ICB response (Figure 3A). We found an obvious enrichment trend: high-PRGPI, to UroA/Genomically unstable, to excluded/desert, to SD/PD. In PRGPI immunophenotyping, the score of desert typing was the highest, and that of inflamed typing was the lowest (Figure 3B). Infiltrated typing was the lowest in Lund typing, and UroA typing was the highest (Figure 3C). There was no difference in the expression of PRGPI between TMB ≥ 10 and TMB < 10. In the analysis of PRGPI and stage, the PRGPI of stage 1 was higher than that of the stage in 2 / 3 / 4. There is no relationship between PRGPI and smoking and gender (Figure 3D). We conducted a chi-square test for immunotherapy response and PRGPI in the IMvigor210 cohort. In the CR/PR group, the low-PRGPI group (35%) proportion was significantly higher than that of the high-PRGPI group (18%); in the SD/PD group, the proportion of the low-PRGPI group (65%) was significantly lower than that of the high-PRGPI (82%), with a p-value of 0.002, indicating that PRGPI can distinguish the ICB treatment response (Figure 3E). We used the ROC curve to compare the predictive ability of PRGPI, TMB, PDCD1, and CD274 for ICB treatment. The results showed that the AUC value of TMB was 0.725, PRGPI was 0.6, PDCD1 was 0.54, CD274 was 0.564 (Figure 3F). The prediction ability of PRGPI is stronger than that of PDCD1 and CD274 but weaker than that of TMB. Therefore, we combine TMB and PRGPI to predict the treatment and prognosis of ICB. It was found that 8% of SD/PD patients with TMB-H/PRGPI-L accounted for 34% of CR / PR patients, which could significantly benefit from ICB treatment (Figure 3G). The overall survival rate of the TMB-H/PRGPI-L group was significantly higher than that of other groups (Figure 3H).
PRGPI can also have a successful prognosis in GSE176307
We used the ICB treatment cohort GSE176307 to evaluate the ability to predict PRGPI. GSE176307 included 34 cases of Atezolizumab treatment, 1 case of Avelumab treatment, 2 cases of Durvalumab treatment, 5 cases of Nivolumab treatment, and 46 cases of Pembrolizumab treatment. KM analysis of PRGPI in GSE176307 showed that the prognosis of H-PRGPI was significantly worse than that of L-PRGPI (Figure 4A). The ROC curve shows that the AUC values of 0.5, 1, and 1.5 years are 0.667, 0.694, and 0.629, respectively. The results showed that PRGPI had a good ability of prognostic classification in GSE176307 (Figure 4B). Figure 4C shows the patient distribution, gene expression heat map, and survival status under PRGPI grouping. We used PRGPI to predict the progression-free survival of patients and found that PRGPI still has a good predictive ability. Both KM analysis and ROC curve show good prediction efficiency (Figure 4D-4E). The PRGPI and ICB treatment response chi-square test showed that PRGPI had a classification ability (p = 0.046, Figure 4F).
PRGPI can also predict survival in the cohort of bladder cancer dominated by chemotherapy
TCGA-BLCA and GSE31684 are two bladder cancer cohorts based on chemotherapy, containing 403 and 93 samples, respectively. The KM analysis of PRGPI in TCGA-BLCA was statistically significant (p = 0.005 and figure 5A). The results of the ROC curve show that the AUC values of 3 years, 5 years, and 7 years are all greater than 0.5. The picture on the right of Figure 5A shows the sample distribution, the heat map of gene expression, and the survival status of H-PRGPI and L-PRGPI in TCGA-BLCA. Similarly, we also verified the prognostic ability of PRGPI at GSE31684 (Figure 5B). These results suggest that PRGPI has a very generalization ability.
Construction of clinical predictive nomogram
To make it easier for doctors to make clinical decisions, we integrated the nomogram of the mUC immunotherapy cohort of the stage, gender, and PRGPI. Figure 6A is the nomogram of the IMvigor210 cohort. The calibration curve shows that the curves of 0.5, 1, and 1.5 years are all close to the middle diagonal line (Figure 6B), and the results of ROC curves show that the AUC values of 0.5, 1, and 1.5 years are 0.595, 0.606 and 0.613 (Figure 6C).
Analysis of characteristics related to PRGPI and immunity
We extracted the expression of 24 HLA molecules from the IMvigor210 cohort. The expression of most HLA molecules in H-PRGPI was lower than that in L-PRGPI, and there was significant statistical significance (Figure 7A). Then, we observed the difference between the immune cell score and the key immune gene set score between H-PRGPI and L-PRGPI. The immune score of the H-PRGPI group was generally decreased (Figure 7B-7C), indicating that the H-PRGPI group showed a state of immunosuppression. Then we explore the relationship between PDCD1, CD274, and PRGPI. The results show that PDCD1 and CD274 are overexpressed in L-PRGPI, and the correlation is more significant than 0.7 (Figure 7D-7G). Finally, we conducted a joint KM analysis of PDCD1 and PRGPI. The prognosis of patients with the H-PDCD1/L-PRGPI subgroup was the best (Figure 7H). The combined analysis of CD274 and PRGPI also showed that patients with H-CD274/L-PRGPI had a better prognosis than those in other groups (Figure 7I). We also verified these results in the GSE176307 cohort. The results of GSE176307 were very similar to those of IMvigor210 (Supplementary figure 4).
Gene set enrichment analysis
We carried out gene set analysis in IMvigor210 and GSE176307 to explore the related pathways of H-PRGPI and L-PRGPI. The results showed that only the related pathway (Figure 8A-8B) enriched to L-PRGPI was found in both cohorts. In IMvigor210, the first six pathways enriched by L-PRGPI are "ANTIGEN PROCESSING AND PRESENTATION", "AUTOIMMUNE THYROID DISEASE", "CELL ADHESION MOLECULES CAMS", "CYTOKINE-CYTOKINE RECEPTOR INTERACTION", "JAK STAT SIGNALING PATHWAY" and "NATURAL KILLER CELL MEDIATED CYTOTOXICITY". In GSE176307, the first six pathways enriched by L-PRGPI are "ANTIGEN PROCESSING AND PRESENTATION", "CELL ADHESION MOLECULES CAMS", "CYTOKINE-CYTOKINE RECEPTOR INTERACTION", "LEISHMANIA INFECTION", "SYSTEMIC LUPUS ERYTHEMATOSUS" and "VIRAL MYOCARDITIS".