Skip to main content
Download PDF
Download ePub

Pyroptosis patterns influence the clinical outcome and immune microenvironment characterization in HPV-positive head and neck squamous cell carcinoma

Infectious Agents and Cancer volume 18, Article number: 30 (2023) Cite this article

Abstract

Background

Head and neck squamous cell carcinoma (HNSCC) is a heterogeneous tumor with diverse molecular pathological profiles. Recent studies have suggested the vital role of pyroptosis in tumor microenvironment. However, the expression patterns of pyroptosis in HPV-positive HNSCC are still unclear.

Methods

Unsupervised clustering analysis was used to identify the pyroptosis patterns based on the RNA-sequencing data of 27 pyroptosis-related genes (PRGs) in HPV-positive HNSCC samples. Random forest classifier and artificial neural network were performed to screen the signature genes associated with pyroptosis, which were verified in two independent external cohorts and qRT-PCR experiment. Principal component analysis was used to develop a scoring system, namely Pyroscore.

Results

The expression variations of 27 PRGs in HPV-positive HNSCC patients were analyzed from genomic and transcriptional domains. Two pyroptosis-related subtypes with distinct clinical outcomes, enrichment pathways and immune characteristics were identified. Next, six signature genes (GZMB, LAG3, NKG7, PRF1, GZMA and GZMH) associated with pyroptosis were selected for prognostic prediction. Further, a Pyroscore system was constructed to determine the level of pyroptosis in each patient. A low Pyroscore was featured by better survival time, increased immune cell infiltration, higher expression of immune checkpoint molecules and T cell-inflamed genes, as well as elevated mutational burden. The Pyroscore was also related to the sensitivity of chemotherapeutic agents.

Conclusions

The pyroptosis-related signature genes and Pyroscore system may be reliable predictors of prognosis and serve as mediators of immune microenvironment in patients with HPV-positive HNSCC.

Introduction

According to the statistics in 2018, head and neck squamous cell carcinoma (HNSCC) was the 6th leading cancer worldwide, which is a heterogeneous malignancy with poor survival rates [1]. Human papillomavirus (HPV) infection is an increasingly common risk factor of HNSCC. HNSCC could be divided into three molecular subtypes: classical (CL), basal (BA) and mesenchymal (MS), whereas HPV-positive HNSCC does not fall into any group and belong to CL and MS subtypes [2]. The multi-omics data also suggest that HPV-positive HNSCC is a unique clinical entity with a specific genetic landscape and propensity for cell death compared to HPV-negative HNSCC [3]. In clinical practice, tumor immunotherapy has revolutionized traditional methods and achieved some gratifying achievements. Food and Drug Administration has recommended immune checkpoint inhibitors (ICIs) such as pembrolizumab and nivolumab as the preferred choices for patients with recurrent or metastatic HNSCC [4]. Although HPV-positive patients have been proven to be more sensitive to immune therapy, only a limited subset of patients benefited from it. Therefore, there is an urgent need to explore the tumor microenvironment of HPV-positive HNSCC to improve clinical effectiveness and targeted therapy efficiency.

The form of cell death, including necroptosis, ferroptosis and pyroptosis, has attracted more and more attention in the context of tumors. Pyroptosis, an inflammatory programmed cell death, is characterized by the cleavage of gasdermin-family perforating proteins, followed by the swelling and rupture of cell membranes, releasing intracellular contents and thus, triggering a strong inflammatory response [5]. It has been reported that pyroptosis-related genes (PRGs) were aberrantly and differentially expressed in different cancers, contributing to tumor suppression or tumor progression [6]. GSDME has been shown to enhance the phagocytosis of macrophages, as well as the number and function of NK cells and CD8(+) T lymphocytes, thereby exerting tumor inhibition [7]. Awad et al. have confirmed that NLRP1 facilitated the progression of skin cancer by mediating the expression of IL-1b and IL-18 cytokines [8]. In HNSCC, Taxol treatment has been observed to induce pyroptosis in tumor cells by mediating Caspase-1/GSDMD [9]. Zhang et al. found that the reducing expression of CD38 could prevent pyroptosis in HNSCC [10]. Notably, a series of studies have shown the close association between pyroptosis and anti-tumor immunity [11]. Given the overlapping and counteracting effects of multiple pyroptotic components, understanding the overall profile of PRGs on specific tumors, rather than the individual regulation of each component, seemed to be a more effective strategy to elucidate the crosstalk of pyroptosis and immune microenvironment in HPV-positive HNSCC.

In this paper, we comprehensively analyzed the expression profile of 27 PRGs in patients with HPV-positive HNSCC on the basis of genomic and transcriptional data. The study aimed to identify the pyroptosis-related subgroups, screen out the candidate genes and develop a pyroptosis scoring method to investigate the survival duration and immunological landscape of HPV-positive HNSCC.

Materials and methods

Data acquisition and processing

The analysis process of the article was presented in Additional file 1: Figure S1. Two HNSCC cohorts (TCGA-HNSCC and GSE65858) were used as the training set. The RNA-sequencing data and related clinical information of patients with HPV-positive HNSCC were retrieved from The Cancer Genome Atlas (TCGA) repository (https://portal.gdc.cancer.gov/repository) and the Gene Expression Omnibus (GEO) dataset (https://www.ncbi.nlm.nih.gov/geo/). The results for fragments per kilobase million were standardized to transcripts per kilobase million. The two datasets were combined, and the “ComBat” algorithm was used to correct batch effects. A total of 138 HPV-positive HNSCC samples and 44 normal samples with the corresponding clinical data, including overall survival (OS) time, age, gender, stage and TNM staging were enrolled. Besides, GSE3292 and GSE6792 databases from the GEO platform were combined and used as the external validation cohort. The raw data were obtained from publicly available databases and no ethical events were involved.

Unsupervised clustering analysis of PRGs

We summarized data on 27 PRGs from previous studies [12,13,14]. The pyroclusters and pyroptosis-related geneclusters were identified by consensus clustering algorithm of R language ConsensuClusterPlus [15]. The selection criteria for the value of K: the modules were relatively uniformly separated and the cumulative distribution function decreased slowly [16].

Gene set variation analysis (GSVA)

The “GSVA” package in R and a reference set (c2.cp.kegg.v7.4.symbols.gmt) retrieved from the MSigDB database were used to run the GSVA assay [17].

Protein-protein interactions (PPI) analysis

Interactions between pyroptosis-related molecules were identified on the STRING website (https://string-db.org/).

Identification and functional enrichment analysis of differentially expressed genes (DEGs)

The DEGs between the two pyroptosis patterns were determined by “limma” R package using the empirical Bayesian method [18]. Adjusted p < 0.001 and logFC = 1 were set as screening criteria. The functional enrichment analysis of DEGs, including the gene ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, were performed by “clusterProfiler” R package [19].

Random forest and artificial neural network screening for signature genes

First, univariate Cox regression analysis was performed to extract significant prognostic DEGs with p < 0.001 as the screening criterion. Next, the prognostic DEGs were input into randomForest software package. The variable importance value was obtained by the Gini coefficient method and decreasing mean square error [20]. Genes with an importance score greater than 3 and ranked the top 6 were selected as signature genes for subsequent analysis. Then, a neural network model of the signature genes was constructed by neuralnet package [21].

Analysis of tumor immune infiltration characteristics

The ESTIMATE algorithm was utilized to assess the immune and stromal scores and the relative purity of the tumor for each patient [22]. The CIBERSORT algorithm was utilized to calculate the proportion of 22 human immune cell subpopulations of each patient [23]. Correlations between the expression of tumor-infiltrating immune cells (TIICs) and Pyroscore were calculated using six algorithms (XCELL, TIMER, QUANTISEQ, MCPCOUNTER, EPIC and CIBERSORT-ABS algorithms) on the TIMER 2.0 platform [24]. Only differential TIICs (p < 0.05) were incorporated for the subsequent correlation analysis between TIICs and Pyroscore.

Quantitative real-time polymerase chain reaction (qRT-PCR)

Five pairs of HPV-positive HNSCC and the adjacent normal tissues were collected from Stomatological Hospital of Xi’an Jiaotong University. The study was approved by the Ethical Review Committee of the Stomatological Hospital of Xi’an Jiaotong University. All patients received written informed consent. All the specimens were subjected to qRT-PCR. Total RNA was extracted from tissues according to the instructions of RNAiso Plus (Takara, Tokyo, Japan) kit. According to PrimeScript™ RT-PCR (Takara) instructions, RNA was reversely transcribed into cDNA in 20-µL volumes. The cDNA was treated by SYBR®Premix Ex Taq™II (Takara). The primers used in this experiment were shown in Additional file 2: Table S1.

Generation of the pyroscore

We constructed a scoring system to assess the pyroptosis expression level of each patient. The procedures for establishment of Pyroscore were as follows:

First, we extracted the overlapping DEGs between the pyroptosis patterns. Then, the HPV-positive HNSCC patients were divided into several groups using an unsupervised clustering method. The consensus clustering algorithm was performed to determine the number of geneclusters as well as their stability. Next, we performed the prognostic analysis for each overlapping DEG using univariate Cox regression analysis and the genes with the significant prognosis were extracted for further analysis. Finally, the principal component analysis (PCA) was conducted to construct pyroptosis relevant score system. Both principal components 1 and 2 were selected to act as signature scores. After obtaining the prognostic value of each gene signature score, we applied a method similar to the gene expression grade index to assess the Pyroscore of each patient [25].

$${Pyroscore}_{i}=\sum (PC{1}_{i}+PC{2}_{i})$$

.

Where i is the expression of overlapping DEGs with a significant prognosis between pyroptosis patterns.

Correlation with drug sensitivity

Based on the Genomics of Drug Sensitivity in Cancer database (https://www.cancerrxgene.org/), we used the R package “pRRophetic” to predict the chemosensitivity of each sample. The estimated semi-inhibitory concentration (IC50) for each specific chemotherapeutic agent was obtained by the “linearRidge” function.

Statistical analysis

Data processing was completed by using the PERL language program (version 5.32.1). Statistical analyses were performed utilizing R software and related calculate packages (version 4.1.1). For clusters of two normally distributed variables, unpaired Student’s t-tests were used for analysis. For clusters of two non-normally distributed variables, the Wilcoxon rank-sum test was used to analyze. For multiple comparisons, Kruskal-Wallis and one-way ANOVA tests were used. Correlations between the variables were analyzed using Spearman’s correlation coefficient. The “Surv-cut point” function in the R package ‘Survminer’ was used to evaluate the cutoff point of each group. The “Surv-cut point” function was based on maximally selected rank statistics, which enabled the fast assessment of all potential cut points, to find the separating partition and determine the best cut-off point of continuous variables [26]. Survival curves were generated using the Kaplan-Meier (K-M) method and the log-rank test was used to determine differences between groups. If not special explanation above, p < 0.05 was deemed as statistically significant and all tests were two-sides.

Results

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
figure 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

Full size image

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
figure 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

Full size image

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
figure 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

Full size image

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
figure 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

Full size image

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
figure 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.

Full size image

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
figure 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

Full size image

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
figure 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

Full size image

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
figure 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

Full size image

Discussion

The intense inflammatory response elicited by pyroptosis was confirmed to affect TIME and tumor progression [30, 31]. Although immunotherapy has made great strides in HNSCC, tumor heterogeneity has yielded unpredictable clinical responses, which was a serious challenge for antitumor therapy. HPV infection has been reported to encompass much of the heterogeneity of HNSCC [2]. Therefore, we focused on the pyroptosis patterns in HPV-positive HNSCC to better understand and manage this disease.

In this study, we firstly explored genetic and transcriptional variations of 27 PRGs in HPV-positive HNSCC and found that the most expression of PRGs was significantly increased in tumor samples than that in normal samples. Then, two pyroptosis subtypes were identified based on the expression of PRGs. Surprisingly, the two pyroptosis patterns exhibited distinct prognosis, biological processes and TIME characteristics. Compared with pyrocluster B, pyrocluster A had significantly better survival time. Besides, pyrocluster A was significantly associated with more enriched immune pathways, higher immune cell abundance, higher expression of immune checkpoint molecules and activated TIICs, which indicated that pyrocluster A was an immune “hot” tumor compared to pyrocluster B. A “hot” tumor indicated a better treatment response to immunotherapy. Studies also have reported that pyroptosis inducers combined with PD-1 inhibitors could turn a “cold” tumor into a “hot” tumor that the immune system could recognize, which demonstrated the enormous therapeutic application value of pyroptosis [32].

Next, we explored the pyroptosis subtypes related genes and found that the inflammatory and immune-related biological processes were involved, especially the activation of MHC proteins and T-cells. Based on these pyroptosis-related genes, six signature genes associated with pyroptosis (GZMB, LAG3, NKG7, PRF1, GZMA and GZMH) were identified and their prognostic value was verified. Given the importance of pyroptosis in prognosis and TIME, we developed a novel scoring system, Pyroscore, to quantify pyroptosis levels in patients with HPV-positive HNSCC. Herein, patients with low Pyroscores had significantly better OS time than those with high Pyroscores. Besides, Pyroscore was negatively associated with the expression of most PRGs, especially GZMA, GZMB and PRF1. Previous reports have shown that when GZMA cleaved GSDMB or GZMB cleaved GSDME, NK cells and CD8(+) T cells directly triggered pyroptosis of tumor cells [7, 32]. The killing of tumor cells by NK cells and cytotoxic T lymphocytes was the ultimate pivotal event of antitumor immunity and was previously thought to be non-inflammatory. In GSDM-expressing tumors, the killing effect of immune cells could be transformed into inflammatory pyroptosis when gasdermin was directly cleaved and activated by granzymes [33]. It was proved that GSDME-mediated pyroptosis induced by chemotherapy drugs played a role in anti-tumor response of oral cancer [34]. Recently, PRF1 has been reported to be a prognostic marker with implications on the immune infiltration of HNSCC [35]. Our findings have suggested the feasibility of targeted pyroptosis to aid the clinical treatment for patients with HNSCC.

TIME refers to the milieu in which tumor cells live, which includes blood vessels, extracellular matrix, immune cells, and signaling molecules, etc [36]. The occurrence and development of tumor is the result of the interaction between tumor cells and their microenvironment. TIICs are closely associated with immunotherapy sensitivity and understanding the composition of immune cells in TIME can reveal the heterogeneity of tumors [37]. In this study, pyrocluster A associated with low Pyroscore was an actived immune signature. Numerous evidence has reported that effector and memory T cells exert tumor-killing functions, whereas Tregs weaken anti-tumor immune responses and exert immunosuppressive effects [38]. Herein, pyrocluster A was infiltrated with CD8(+) T cells and CD4 memory-activated T cells, while pyrocluster B showed higher infiltration of CD4 memory-resting T cells and Tregs. Besides, Pyroscore was also strongly correlated with CD8(+) T cells. M1-type macrophages exerted pro-inflammatory and anti-tumor effects, but tumor-associated macrophages in TIME were M2-type and promoted angiogenesis and tumor invasion by secreting Th2 cytokines [39]. Pyrocluster A had higher M1 macrophage cells infiltration than pyrocluster B. Recent studies revealed that B cells were the strongest prognostic element and involved in the immune response [40]. In this study, pyrocluster A with better survival time had higher infiltration of B cells, which was consistent with the previous studies.

Tumor cell resistance to the immunosuppressive microenvironment, such as through PD-1, is a major dilemma in antitumor therapy. Blocking immune checkpoint proteins through ICIs can reawaken the immune system. However, once activated, T cells must be sufficient to distinguish tumors from normal cells. The presence of immunogenic neoantigens on the surface of cancer cells helps the immune system to recognize cancer cells in the context of MHC [41]. Due to the increased neoantigens produced by somatic tumor mutations, the higher the number of mutations (higher the TMB), the higher the likelihood that certain neoantigens presented by MHC proteins are immunogenic and thus can help T cells to recognize and eradicate cancer cells [29]. Besides, studies have demonstrated that T cell-inflamed GEP was a unique feature of T-cell activation and could manifest clinical antitumor efficacy [28]. The high T cell-inflamed GEP was associated with a significant pan-cancer survival benefit and favored anti-PD-1 immunotherapy [42]. Patients in pyrocluster A with low Pyroscore had higher TMB, elevated expression of immune checkpoint molecules (CTAL-4, GAL9, LAG3, PD-1 and PD-L1) and T cell–inflamed genes, which indicates a better response to immunotherapy, especially to ICIs treatment. In addition, patients in high and low Pyroscore groups had different sensitivities to commonly used chemotherapeutic agents, suggesting that Pyroscore could assist in guiding clinical drug administration.

However, there are some limitations in this study. First, we did not find clinical data on HPV-positive HNSCC patients treated with ICIs to further validate our results. Second, the original data used in this study is retrospective from the cross-queue, and prospective studies are necessary to avoid the inherent errors of retrospective research in the future.

Conclusions

The need for personalized treatment of patients has been proposed because of the heterogeneity of tumors. We comprehensively analyzed the pyroptosis patterns of HPV-positive HNSCC and revealed the underlying biological processes and immune infiltration characteristics affected by the PRGs. Besides, the pyroptosis-related signature genes and the Pyroscore system were developed to explore the feasibility of targeted pyroptosis assisted immunotherapy for HPV-positive HNSCC. The landscape of pyroptosis facilitates our understanding of the TIME and lays the foundation to improve the immunotherapeutic efficacy of HNSCC.

Data Availability

The raw data in this study could be downloaded from the following databases: (https://portal.gdc.cancer.gov/repository and https://www.ncbi.nlm.nih.gov/geo/). The relevant codes are available from the corresponding author on reasonable requests.

Abbreviations

AUC:

Area under the curve

CNV:

Copy number variation

DEGs:

Differentially expressed genes

GEO:

Gene Expression Omnibus

GEP:

Gene expression profile

GO:

Gene ontology

GSVA:

Gene set variation analysis

HNSCC:

Head and neck squamous cell carcinoma

HPV:

Human papillomavirus

ICIs:

Immune checkpoint inhibitors

K-M:

Kaplan-Meier

KEGG:

Kyoto Encyclopedia of Genes and Genomes

MHC:

Major histocompatibility complex

OS:

Overall survival

PPI:

Protein-protein interactions

PRGs:

Pyroptosis-related genes

qRT-PCR:

Quantitative real-time polymerase chain reaction

TCGA:

The Cancer Genome Atlas

TIICs:

Tumor-infiltrating immune cells

TIME:

Tumor immune microenvironment

TMB:

Tumor mutational burden

References

  1. Shield KD, Ferlay J, Jemal A, Sankaranarayanan R, Chaturvedi AK, Bray F, et al. The global incidence of lip, oral cavity, and pharyngeal cancers by subsite in 2012. CA Cancer J Clin. 2017;67(1):51–64.

    Article  PubMed  Google Scholar 

  2. Raudenska M, Balvan J, Masarik M. Cell death in head and neck cancer pathogenesis and treatment. Cell Death Dis. 2021;12(2):192.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  3. Zhou C, Parsons JL. The radiobiology of HPV-positive and HPV-negative head and neck squamous cell carcinoma. Expert Rev Mol Med. 2020;22:e3.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  4. Long Q, Huang C, Meng Q, Peng J, Yao F, Du D, et al. TNF patterns and Tumor Microenvironment characterization in Head and Neck squamous cell carcinoma. Front Immunol. 2021;12:754818.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  5. Liu X, Lieberman J. A mechanistic understanding of pyroptosis: the fiery death triggered by invasive infection. Adv Immunol. 2017;135:81–117.

    Article  CAS  PubMed  Google Scholar 

  6. Loveless R, Bloomquist R, Teng Y. Pyroptosis at the forefront of anticancer immunity. J Exp Clin Cancer Res. 2021;40(1):264.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  7. Zhang Z, Zhang Y, Xia S, Kong Q, Li S, Liu X, et al. Gasdermin E suppresses tumour growth by activating anti-tumour immunity. Nature. 2020;579(7799):415–20.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Awad F, Assrawi E, Louvrier C, Jumeau C, Giurgea I, Amselem S, et al. Photoaging and skin cancer: is the inflammasome the missing link? Mech Ageing Dev. 2018;172:131–7.

    Article  CAS  PubMed  Google Scholar 

  9. Wang X, Li H, Li W, Xie J, Wang F, Peng X, et al. The role of Caspase-1/GSDMD-mediated pyroptosis in Taxol-induced cell death and a taxol-resistant phenotype in nasopharyngeal carcinoma regulated by autophagy. Cell Biol Toxicol. 2020;36(5):437–57.

    Article  PubMed  Google Scholar 

  10. Zhang MJ, Gao W, Liu S, Siu SP, Yin M, Ng JC, et al. CD38 triggers inflammasome-mediated pyroptotic cell death in head and neck squamous cell carcinoma. Am J Cancer Res. 2020;10(9):2895–908.

    CAS  PubMed  PubMed Central  Google Scholar 

  11. Tang R, Xu J, Zhang B, Liu J, Liang C, Hua J, et al. Ferroptosis, necroptosis, and pyroptosis in anticancer immunity. J Hematol Oncol. 2020;13(1):110.

    Article  PubMed  PubMed Central  Google Scholar 

  12. Van Opdenbosch N, Lamkanfi M. Caspases in cell death, inflammation, and Disease. Immunity. 2019;50(6):1352–64.

    Article  PubMed  PubMed Central  Google Scholar 

  13. Broz P, Pelegrin P, Shao F. The gasdermins, a protein family executing cell death and inflammation. Nat Rev Immunol. 2020;20(3):143–57.

    Article  CAS  PubMed  Google Scholar 

  14. Chen X, Chen H, Yao H, Zhao K, Zhang Y, He D, et al. Turning up the heat on non-immunoreactive tumors: pyroptosis influences the tumor immune microenvironment in bladder cancer. Oncogene. 2021;40(45):6381–93.

    Article  CAS  PubMed  Google Scholar 

  15. Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics. 2010;26(12):1572–3.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  16. Qin X, Li J, Hu W, Yang J. Machine Learning K-Means Clustering Algorithm for Interpolative Separable Density Fitting to accelerate Hybrid Functional Calculations with Numerical Atomic Orbitals. J Phys Chem A. 2020;124(48):10066–74.

    Article  CAS  PubMed  Google Scholar 

  17. Hanzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013;14:7.

    Article  PubMed  PubMed Central  Google Scholar 

  18. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43(7):e47.

    Article  PubMed  PubMed Central  Google Scholar 

  19. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16(5):284–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Tian Y, Yang J, Lan M, Zou T. Construction and analysis of a joint diagnosis model of random forest and artificial neural network for heart failure. Aging. 2020;12(24):26221–35.

    Article  PubMed  PubMed Central  Google Scholar 

  21. Beck MW. NeuralNetTools: visualization and analysis tools for neural networks. J Stat Softw. 2018;85(11):1–20.

    Article  PubMed  PubMed Central  Google Scholar 

  22. Yoshihara K, Shahmoradgoli M, Martinez E, Vegesna R, Kim H, Torres-Garcia W, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013;4:2612.

    Article  PubMed  Google Scholar 

  23. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. 2015;12(5):453–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Li T, Fu J, Zeng Z, Cohen D, Li J, Chen Q, et al. TIMER2.0 for analysis of tumor-infiltrating immune cells. Nucleic Acids Res. 2020;48(W1):W509–W14.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Sotiriou C, Wirapati P, Loi S, Harris A, Fox S, Smeds J, et al. Gene expression profiling in breast cancer: understanding the molecular basis of histologic grade to improve prognosis. J Natl Cancer Inst. 2006;98(4):262–72.

    Article  CAS  PubMed  Google Scholar 

  26. Hothorn T, Zeileis A. Generalized maximally selected statistics. Biometrics. 2008;64(4):1263–9.

    Article  PubMed  Google Scholar 

  27. Roche PA, Furuta K. The ins and outs of MHC class II-mediated antigen processing and presentation. Nat Rev Immunol. 2015;15(4):203–16.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Ott PA, Bang YJ, Piha-Paul SA, Razak ARA, Bennouna J, Soria JC, et al. T-Cell-inflamed gene-expression Profile, programmed death ligand 1 expression, and Tumor Mutational Burden Predict Efficacy in patients treated with Pembrolizumab Across 20 cancers: KEYNOTE-028. J Clin Oncol. 2019;37(4):318–27.

    Article  PubMed  Google Scholar 

  29. Alspach E, Lussier DM, Miceli AP, Kizhvatov I, DuPage M, Luoma AM, et al. MHC-II neoantigens shape tumour immunity and response to immunotherapy. Nature. 2019;574(7780):696–701.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Hou J, Hsu JM, Hung MC. Molecular mechanisms and functions of pyroptosis in inflammation and antitumor immunity. Mol Cell. 2021;81(22):4579–90.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Hsu SK, Li CY, Lin IL, Syue WJ, Chen YF, Cheng KC, et al. Inflammation-related pyroptosis, a novel programmed cell death pathway, and its crosstalk with immune therapy in cancer treatment. Theranostics. 2021;11(18):8813–35.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  32. Zhou Z, He H, Wang K, Shi X, Wang Y, Su Y et al. Granzyme A from cytotoxic lymphocytes cleaves GSDMB to trigger pyroptosis in target cells. Science. 2020;368(6494).

  33. Zhang Z, Zhang Y, Lieberman J. Lighting a fire: can we harness pyroptosis to Ignite Antitumor Immunity? Cancer Immunol Res. 2021;9(1):2–7.

    Article  PubMed  PubMed Central  Google Scholar 

  34. Wang S, Zhang MJ, Wu ZZ, Zhu SW, Wan SC, Zhang BX, et al. GSDME is related to prognosis and response to chemotherapy in oral Cancer. J Dent Res. 2022;101(7):848–58.

    Article  CAS  PubMed  Google Scholar 

  35. Fan C, Hu H, Shen Y, Wang Q, Mao Y, Ye B, et al. PRF1 is a prognostic marker and correlated with immune infiltration in head and neck squamous cell carcinoma. Transl Oncol. 2021;14(4):101042.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Turley SJ, Cremasco V, Astarita JL. Immunological hallmarks of stromal cells in the tumour microenvironment. Nat Rev Immunol. 2015;15(11):669–82.

    Article  CAS  PubMed  Google Scholar 

  37. Hinshaw DC, Shevde LA. The Tumor Microenvironment innately modulates Cancer Progression. Cancer Res. 2019;79(18):4557–66.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Oh DY, Kwek SS, Raju SS, Li T, McCarthy E, Chow E, et al. Intratumoral CD4(+) T cells mediate anti-tumor cytotoxicity in human bladder Cancer. Cell. 2020;181(7):1612–25. e13.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Pan Y, Yu Y, Wang X, Zhang T. Tumor-Associated Macrophages in Tumor Immunity. Front Immunol. 2020;11:583084.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Petitprez F, de Reynies A, Keung EZ, Chen TW, Sun CM, Calderaro J, et al. B cells are associated with survival and immunotherapy response in sarcoma. Nature. 2020;577(7791):556–60.

    Article  CAS  PubMed  Google Scholar 

  41. Kreiter S, Vormehr M, van de Roemer N, Diken M, Lower M, Diekmann J, et al. Mutant MHC class II epitopes drive therapeutic immune responses to cancer. Nature. 2015;520(7549):692–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Trujillo JA, Sweis RF, Bao R, Luke JJ. T cell-inflamed versus Non-T cell-inflamed tumors: a conceptual Framework for Cancer Immunotherapy Drug Development and Combination Therapy Selection. Cancer Immunol Res. 2018;6(9):990–1000.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

We sincerely acknowledge the contributions of open projects.

Funding

This study was funding by the National Nature Science Foundation of China (grant number: 61627811).

Author information

Authors and Affiliations

  1. Key Laboratory of Shaanxi Province for Craniofacial Precision Medicine Research, College of Stomatology, Xi’an Jiaotong University, 98# XiWu Road, Xi’an, 710004, Shaanxi, P.R. China

    Doudou Li & Yuxia Hou

  2. Department of Orthodontics, College of Stomatology, Xi’an Jiaotong University, 98# XiWu Road, Xi’an, 710004, Shaanxi, P.R. China

    Doudou Li & Yuxia Hou

  3. Department of Oral and Maxillofacial Surgery, College of Stomatology, Xi’an Jiaotong University, 98# XiWu Road, Xi’an, 710004, Shaanxi, P.R. China

    Dong Ma

Authors
  1. Doudou Li
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Dong Ma
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Yuxia Hou
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

Conception and design: L.D.D. and H.Y.X.; Acquisition of data: L.D.D. and M.D.; Analysis and interpretation of data: L.D.D.; Drafting the article: L.D.D.; Revising the article: L.D.D. and H.Y.X. All authors read and approved the final paper.

Corresponding author

Correspondence to Yuxia Hou.

Ethics declarations

Ethics approval and consent to participate

The study was approved by the Ethical Review Committee of the Xi’an Jiaotong University (NO.2020370). All patients received written informed consent. The research adhered to the principles of the Declaration of Helsinki.

Consent for publication

Not applicable.

Competing interests

The authors declare no conflict of interest.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Electronic supplementary material

Below is the link to the electronic supplementary material.

Additional file 1

:

Figure S1 The workflow of the study.

Figure S2 Clustering of pyroclusters and consensus matrix heatmaps for k = 3–9.

Figure S3 Unsupervised clustering of pyroptosis-related geneclusters and consensus matrix heatmaps for k = 3–9.

Figure S4 Relative expression levels of 27 PRGs between the two geneclusters. *p < 0.05, **p < 0.01, and ***p < 0.001.

Additional file 2

:

Table S1 The primer sequences for qRT-PCR.

Table S2 The activation states of biological pathways in pyroptosis subtypes by GSVA enrichment analysis.

Table S3 Functional annotation of the differentially expressed genes between the two pyroptosis subtypes.

Table S4 Prognostic analysis of 66 pyroptosis-related genes using a univariate Cox regression analysis.

Table S5 The output results of the neural network model.

Table S6 Correlation between tumor-infiltrating immune cells and Pyroscore.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Li, D., Ma, D. & Hou, Y. Pyroptosis patterns influence the clinical outcome and immune microenvironment characterization in HPV-positive head and neck squamous cell carcinoma. Infect Agents Cancer 18, 30 (2023). https://doi.org/10.1186/s13027-023-00507-w

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s13027-023-00507-w

Keywords

Infectious Agents and Cancer

ISSN: 1750-9378

Contact us