Overview of genetic and transcriptional variations of PRGs in HPV-positive HNSCC

In this study, protein-protein interaction networks functional enrichment analysis of the 27 pyroptosis-related molecules were analyzed by STRING platform, which showed widespread protein interactions of PRGs (Fig. 1A). Next, somatic copy number variation (CNV) of these PRGs in HPV-positive HNSCC were analyzed (Fig. 1B). Among them, the CNV of APIP, GSDMD, GSDMC and AIM2 significantly increased, while the CNV of CASP4, CASP5 and CASP1 remarkably decreased. Figure 1 C has shown the chromosomal localization information of these PRGs with CNV. In addition, we compared PRGs expression in normal and HPV-positive HNSCC tissues. As shown in Fig. 1D, the majority of PRGs were significantly elevated in tumor tissue, while GSDMA and CTSG were significantly decreased. Collectively, the above results have shown the genomic and transcription profiles of PRGs and their significant differences between normal and HPV-positive HNSCC tissues, suggesting the imbalance of PRGs and its potential role in the occurrence of HPV-positive HNSCC.

Fig. 1

Genomic and transcriptional variations of pyroptosis-related genes (PRGs) in HPV-positive HNSCC. (A) The potential protein-protein interaction networks among 27 pyroptosis-related proteins were drawn on the STRING website. (B) Frequency of copy number variations (CNV) gain and loss in 27 PRGs. (C) The chromosomal localization information of PRGs with CNV. (D) Expression levels of 27 PRGs in HPV-positive HNSCC and normal tissues. *p < 0.05, **p < 0.01, and ***p < 0.001

Identification of two different pyroptosis subtypes

To reveal the expression patterns of pyroptosis in HPV-positive HNSCC. The RNA expression of 138 HPV-positive HNSCC patients from TCGA-HNSCC and GSE65858 databases was extracted. With the unsupervised clustering algorithm in R language, we found that when the samples were clustered into two modules (k = 2), the clustering stability achieved satisfactory results (Fig. 2A and B and Additional file 1: Figure S2). The tumor samples were divided into two types, namely pyrocluster A (n = 66) and pyrocluster B (n = 72). K-M analysis showed that patients in pyrocluster A had a significantly better OS time than those in pyrocluster B (p = 0.037, Fig. 2C). Then, we examined the association between clinicopathological features and scorch molecule expressions in the two pyroclusters. As shown in Fig. 2D, pyrocluster B had significantly higher mortality and lower T staging (p < 0.05) compared to pyrocluster A. The majority of PRGs was significantly higher in pyrocluster A than that in pyrocluster B (Fig. 2E).

Fig. 2

Unsupervised learning to identify two pyroclusters based on the expression of 27 PRGs. (A) Consensus clustering matrix defining two patterns (k = 2) in TCGA-HNSCC and GSE65858 cohorts. (B) The cumulative distribution function of clustering. (C) The two pyroclusters have a substantial difference in overall survival time according to Kaplan-Meier survival analysis. (D) The heatmap shows the clinicopathological traits and unsupervised clustering of the 27 PRGs in TCGA-HNSCC and GSE65858 cohorts. Red denotes higher expression, whereas blue denotes lower expression. (E) Relative expression levels of 27 PRGs between the two pyroclusters

The two pyroptosis patterns associated with distinct biological processes and immune infiltration characteristics

Then, we analyzed the enrichment pathways of the two patterns to explore the underlying biological processes of pyroptosis in HPV-positive HNSCC. The GSVA enrichment analysis showed the differential enrichment pathways between the two pyroclusters. As the results suggested, pyrocluster A was enriched in immune-activated hallmarks such as Toll-like receptor signaling pathway, NOD-like receptor signaling pathway, T and B cell receptor signaling pathways (Fig. 3A and Additional file 2: Table S2). To further understand the influence of pyroptosis-related molecules on TME in HPV-positive HNSCC patients, we utilized the ESTIMATE algorithm to assess the immune and stromal scores as well as the relative tumor purity of each patient. We found that patients in pyrocluster A had remarkably higher immunescores and estimatescores, while patients in pyrocluster B had remarkably higher tumor purity (Fig. 3B C). Given the unique role of immune checkpoint molecules in tumor immune microenvironment (TIME), we explored the expression of PD-1, PD-L1, LAG3, GAL9 and CTLA4 and found that the expression of these molecules was significantly elevated in pyrocluster A than that in pyrocluster B (Fig. 3D). Moreover, we estimated the infiltration characteristics of 22 immune cells in pyroptosis patterns using the CIBERSORT algorithm. The results showed that the activated TIICs, especially B memory cells, CD8(+) T cells, CD4 memory-activated T cells, activated NK cells, and M1 macrophage cells infiltration were notably higher in pyrocluster A than those in pyrocluster B, while CD4 memory-resting T cells, regulatory T cells (Tregs), activated dendritic cells and activated mast cells infiltration were notably higher in pyrocluster B than those in pyrocluster A (Fig. 3E).

Fig. 3

Biological processes and immune characteristics of the two pyroclusters. (A) The heatmap shows the biological processes in the two pyroclusters analyzed by GSVA. Red denotes activated pathways and blue denotes inhibited pathways. (B-C) Tumor immune microenvironment scores (B) and tumor purity (C) between the pyroclusters. (D) The expression of immune checkpoint molecules between the two pyroclusters. (E) The abundance of 22 tumor immune infiltrating cells between the two pyroclusters

Two scorch gene subtypes were recognized using prognostic DEGs

By conducting the “limma” package, the DEGs between the two pyroptosis patterns were identified. To explore the biological mechanisms involved in pyroptosis-related molecules, functional enrichment analysis of DEGs was performed. GO analysis exhibited the top 5 cellular components, molecular functions and biological processes respectively, revealing that immune-related processes, especially major histocompatibility complex (MHC) proteins and T-cell were activated (Fig. 4A and Additional file 2: Table S3). The antigen processing and presentation process by MHC II affected T cell recognition and tumor cell killing [27]. KEGG analysis indicated that the DEGs were mainly enriched in immune and inflammation related pathways, such as Th1 and Th2 cell differentiation, Th17 cell differentiation and inflammatory bowel disease (Fig. 4B and Additional file 2: Table S3).

Next, the univariate Cox regression analysis was performed to screen out the prognosis related DEGs. With p < 0.001 as the screening criterion, a total of 66 prognostic DEGs between the two scorch patterns were identified (Additional file 2: Table S4). Consequently, 138 HNSCC HPV-positive patients were classified into 2 genomic subtypes using the consensus clustering algorithm based on the expression of the 66 DEGs (Fig. 4C and Additional file 1: Figure S3). K-M analysis method showed the significant difference of OS time between the gene subtypes (p = 0.023) (Fig. 4D). The heatmap showed the expression of the 66 genes in the geneclusters and the pyroptosis patterns as well as the relationship between these genes and clinical parameters (Fig. 4E). In addition, the two geneclusters exhibited significant differences in PRGs expression (Additional file 1: Figure S4).

Fig. 4

Generation of three pyroptosis-related geneclusters based on prognostic differentially expressed genes (DEGs). (A-B) Gene Ontology (A) and Kyoto Encyclopedia of Genes and Genomes (B) enrichment analysis of DEGs between the two pyroclusters. (C) Consensus clustering matrix defining two genomic subtypes (k = 2). (D) Kaplan-Meier survival analysis shows significant differences in overall survival time between genecluster A and genecluster B. (E) The heatmap shows the correlation between clinicopathological traits and the geneclusters. Red denotes high expression and blue denotes low expression

Random forest and neural network algorithms to explore the signature genes

To find the signature genes associated with the pyroptosis, we utilized the random forest algorithm to assess the importance of 66 prognostic DEGs (Fig. 5A). The variable importance was measured by the Gini coefficient method and we defined genes with an importance score greater than 3 and ranked the top 6 as the signature genes, including GZMB, LAG3, NKG7, PRF1, GZMA and GZMH (Fig. 5B). Consequently, the artificial neural network model depicting the distribution of the signature genes in pyroclusters was shown in Fig. 5C. The weights of the signature genes were shown in Additional file 2: Table S5. The heatmap showed the signature genes could distinguish the two pyroclusters (Fig. 5D). It could be seen that all signature genes were low expressed in pyrocluster A, whereas high expressed in pyrocluster B. Then, the classification efficiency of the model was evaluated by receiver operating characteristic (ROC) curves in the training group (TCGA-HNSCC and GSE65858 cohorts). As shown in Fig. 5E, the area under the curve (AUC) of the training group was up to 0.997. The performance was confirmed in the external GEO datasets (GSE3292 and GSE6792 cohorts). As shown in Fig. 5F, the AUC of the test group was up to 0.892. In addition, the mRNA expression levels of the signature genes in HPV-positive HNSCC tissues were verified to be higher than that in adjacent normal tissues detected by qRT-PCR (Fig. 5G).

Fig. 5

Exploration of the signature genes by random forest classifier and neural network model. (A) The relationship plot between the error rate and the number of decision trees. When the number of decision trees is about 400, the error rate is relatively stable. (B) The importance score of the signature genes using the Gini coefficient method. The X-axis denotes the signature genes, and the Y-axis denotes the importance index. (C) The neural network diagram demonstrates the relationship between the signature genes and the two pyroclusters. (D) The signature genes were differentially expressed in the two pyroclusters. (E-F) The receiver operating characteristic curve analysis of the classification efficiency of the signature genes in training sets (TCGA-HNSCC and GSE65858 cohorts) and external testing sets (GSE3292 and GSE6792 cohorts). (G) Expression levels of the 6 signature genes in HPV-positive HNSCC and their adjacent normal tissues using qRT-PCR.

Construction and validation of the prognostic pyroscore

The above studies demonstrated a robust association among pyroptosis, tumor immunity and survival in HNSCC HPV-positive patients, so we introduced a unique scoring scheme, the Pyroscore, to assess the pyroptosis expression level of each patient. By employing principal component analysis, the sum of the top two contributing components was used as the Pyroscore for each patient. By conducting the R package “Survminer”, the optimal cut-off was calculated and then patients were divided into high Pyroscore and low Pyroscore groups. As shown in Fig. 6A, the distributions of genecluster, pyrocluster, Pyroscore, and survival status were plotted as a Sankey diagram. Of note, Pyroscore was confirmed to be effective in evaluating the pyroclusters and geneclusters. Comparing the level of scorch death among the subtypes, it was found that pyrocluster A and genecluster A had notably lower Pyroscore, while pyrocluster B and genecluster B had notably higher Pyroscore (Fig. 6B C). In addition, in the training and two validation cohorts, K-M analysis showed significantly higher OS time in low Pyroscore groups than those in high groups (p = 0.006; p = 0.024; p = 0.035;Fig. 6D - F). Moreover, Pyroscore was found to be negatively associated with the expression of most PRGs, especially GZMA, GZMB and PRF1 (Fig. 6G). Overall, Pyrocluster A and genecluster A with lower Pyroscore showed better OS time, which was consistent with our previous results.

Fig. 6

Construction and validation of Pyroscore. (A) Alluvial diagram showing the distributions of pyroclusters, geneclusters, Pyroscores and survival outcomes. (B) Pyroscore is differentially expressed in pyrocluster A and B. (C) Pyroscore is differentially expressed in genecluster A and B. (D) Kaplan-Meier survival analysis shows a significant difference in overall survival time between high and low Pyroscore groups in training sets (TCGA-HNSCC and GSE65858 cohorts). (E) Kaplan-Meier analysis of overall survival time between high and low Pyroscore groups in internal testing set (TCGA-HNSCC cohort). (F) Kaplan-Meier analysis of overall survival time between high and low Pyroscore groups in internal testing set (GSE65858 cohort). (G) The correlations among Pyroscore and PRGs. Red dots denote positive correlations and blue dots denote negative correlations. Larger and deeper dots represent stronger correlations

Pyroscore associated with immune cell infiltration and TME

The abundance of immune cells was thoroughly examined using XCELL, TIMER, QUANTISEQ, MCPCOUNTER, EPIC and CIBERSORT-ABS algorithms. Surprisingly, Pyroscore was negatively correlated with the level of infiltration of most immune cells, except for uncharacterized cells (EPIC) and NK-resting cells (CIBERSORT-ABS) (Fig. 7A and Additional file 2: Table S6). Of note, CD 8(+) T cell was strongly correlated with Pyroscore through the detection by different platforms. Then, we used the ESTIMATE algorithm to evaluate the immune and matrix components in the TME. Patients in the low Pyroscore group had higher stromalscore, immunescore and ESTIMATEscore as well as lower tumor purity compared to the high Pyroscore group (Fig. 7B C). In addition, T cell–inflamed gene expression profile (GEP) is a symbol of CD 8(+) T cell activation and closely reflects the tumor’s sensitivity to immune response [28]. Of the 18 T cell-inflamed genes, 17 genes were considerably higher in the low Pyroscore groups than that in the high Pyroscore groups (Fig. 7D).

Fig. 7

The immune characteristics of different Pyroscore subgroups. (A) Bubble plot depicting the correlation between Pyroscore and immune infiltrating cells. (B-C) Tumor immune microenvironment scores (B) and tumor purity (C) in different Pyroscore subgroups. (D) The expression profile of 18 T cell-inflamed genes in different Pyroscore subgroups

Clinical implications of pyroscore in antitumor therapy

Considering the critical value of immune checkpoints in immunotherapy, we compared the expression levels of CTAL-4, GAL9, LAG3, PD-1 and PD-L1 in Pyroscore groups. We found that the low Pyroscore group had higher expression of CTAL-4, GAL9, LAG3, PD-1and PD-L1, suggesting that patients with low Pyroscore may be more sensitive to immunosuppressive therapies (Fig. 8A). In addition, studies have shown that tumor mutational burden (TMB) may assist in predicting clinical response, which may be associated with neoantigens generated by the mutation, thereby enhancing the immune system’s ability to recognize cancer cells [29]. As shown in Fig. 8B, we found that patients with low Pyroscore showed a significantly higher TMB than patients with high Pyroscore. Further, compared with the high TMB group, patients with low TMB have a better OS time (Fig. 8C). We next selected the commonly used chemotherapeutic agents for the treatment of HNSCC in order to examine drug sensitivity variations between patients with high and low Pyroscores. To compute IC50 drug sensitivity values, we utilized the R package ‘pRRophetic.‘ The IC50 values of paclitaxel and docetaxel were considerably lower in the high Pyroscore group, whereas the values of axitinib, methotrexate, and gemcitabine were significantly higher in the high Pyroscore group compared to the low Pyroscore group (Fig. 8D-H).

Fig. 8

The role of Pyroscore in anti-tumor therapy. (A) The expression of immune checkpoint molecules in Pyroscore groups. (B) The tumor mutational burden in Pyroscore groups. (C) Kaplan-Meier survival analysis shows significant differences in overall survival time between high and low Pyroscore groups. (D-G) The sensitivity of commonly used drugs is significantly different in Pyroscore groups