Prognostic and chemotherapeutic implications of a novel four-gene pyroptosis model in head and neck squamous cell carcinoma

Data extraction and analysis

RNA-sequencing data (497 tumor samples and 44 normal samples, Transcripts Per Million value) were extracted from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/), specifically focusing on the expression levels of four PRG genes (CASP1, CASP3, GSDME, GSDMD). Further validation analysis was carried out using the GSE41613 cohort, which includes 97 HNSCC patients, deriving from the Gene Expression Omnibus (GEO). Any batch effect resulting from nonbiotech bias was corrected using the “combat” algorithm from the “sva” R package. Corresponding clinicopathologic data (including clinical stage, age, sex, grade, T-stage, N-stage, overall survival (OS) time, and status) were also collected.

Generation of pyroptosis-related gene score (PRGscore)

Four pyroptosis-related genes (CASP1, CASP3, GSDME, GSDMD), previously indicated to be involved in the initiation and development of HNSCC, were defined as the pyroptosis-related gene set (Wang et al., 2022; Rioja-Blanco et al., 2022; Zi et al., 2023; Zhang et al., 2020; Wang & Liu, 2023; Yue et al., 2019). A human gene transfer format (GTF) file was obtained from the Ensembl dataset (http://asia.ensembl.org) for gene ID annotation, which enabled the identification of specific mRNAs.

The expression levels of these pyroptosis-related genes were extracted and compared between HNSCC and corresponding normal tissue using the “limma” package. A gene set variation analysis (GSVA) index for each sample was generated by calculating a gene set enrichment score using the z-score augmenting method. The GSVA index was then considered as the PRGscore, representing each sample’s pyroptotic level. Pearson correlation analysis was conducted to compute the correlation between gene expression levels and the GSVA index. An association was considered significant if the correlation coefficient was >0.1 and P was < 0.05.

Prognostic analysis based on the PRGscore

The “surv_cutpoint” function from the “survminer” R package was used to determine the ideal cut-off values. These values were used to categorize samples into high or low PRGscore groups. Associations between the PRGscore and clinicopathological parameters were assessed using the chi-squared test (Table S1). A risk curve was plotted using R software to display the distinct survival status associated with high or low PRGscore for individuals from the TCGA or GEO cohorts. Kaplan–Meier survival and multivariable Cox proportional hazard regression analyses were performed to evaluate the PRGscore’s predictive value. To verify the COX proportional hazards assumption, the global proportional hazards assumption test using “cox.zph()” function, and graphed scaled Schoenfeld residuals of each covariate against the transformed time using “ggcoxzph()” function were applied. The prognostic value of the PRGscore was further validated using multivariable Cox regression and Kaplan–Meier survival analyses in the GSE41613 dataset.

IC50 estimation of routine chemotherapeutic drugs for HNSCC

The half-maximum inhibitory concentration (IC50) of six commonly used chemotherapeutic drugs for HNSCC (cisplatin, methotrexate, gemcitabine, 5-fluorouracil, paclitaxel, and Docetaxel) was estimated using the “pRRophetic” package (Geeleher, Cox & Huang, 2014). The “pRRophetic” R package is typically used to predict clinical chemotherapeutic responses using gene expression and drug sensitivity data from cell lines in the Cancer Genome Project (CGP). After extracting data from the CGP, a ridge regression model was established using certain genes as predictors and drug sensitivity (IC50) values as the outcome variable. Subsequently, the processed, standardized, filtered clinical tumor expression data (HNSCC data obtained from TCGA) was employed in this model to estimate drug sensitivity.

Immunotherapeutic efficacy predicting analysis

The study involved an analysis of the varying expression levels of immune checkpoint inhibitor -related genes, specifically Programmed Cell Death 1 Ligand 1(CD274), Cytotoxic T-Lymphocyte-Associated Protein 4(CTLA4), Programmed Cell Death 1(PDCD1), Indoleamine 2,3-Dioxygenase 1(IDO1), Programmed Cell Death 1 Ligand 2(PDCD1LG2), and Hepatitis A Virus Cellular Receptor 2(HAVCR2) between two cohorts. Additionally, to gauge the potential response to immunotherapy, the Immunopheno score (IPS) score from the TCGA-HNSC project was acquired through The Cancer Immunome Atlas (TCIA) database (https://tcia.at/home). This scoring was computed based on the gene expression data of immune-related genes.

Cell culture

Fourteen HNSCC cell lines, namely HSC-4, CAL-27, UM1, HN4, HN12, HN30, CAL-33, SCC-47, SCC-9, H314, H103, H357, H376, and SCC-25 were utilized for this study. Eight of these lines (HSC-4, CAL-27, UM1, HN4, HN12, HN30, CAL-33, and SCC-47) were maintained in DMEM/high-glucose (Cytiva, Logan, USA) supplemented with 10% fetal bovine serum (FBS, VivaCell, Shanghai, China) and 1% streptomycin–penicillin (Biosharp, Hefei, China). The remaining six lines, SCC-9, H314, H103, H357, H376, and SCC25, were grown in F12 medium (VivaCell, Shanghai, China) enriched with 10% FBS, 1% streptomycin–penicillin, and 0.5 µg/mL hydrocortisone. The cells were incubated in a humidified environment at 37 °C with 5% CO2. All these cells were carefully maintained in the central cell bank of State Key Laboratory of Oral Diseases, West China Hospital of Stomatology in Sichuan University. Cell lines were routinely tested to be negative for mycoplasma, and low passage cultures were applied in this experiment.

Quantitative polymerase chain reaction (q-PCR)

Total RNA from the HNSCC cells was isolated using TRIzol reagent (Invitrogen, San Diego, CA, USA). cDNA was synthesized from 2 µg of total RNA using the PrimeScript RT Reagent kit (TaKaRa, Shiga, Japan). The ChamQ Universal SYBR qPCR Master Mix reaction system (Vazyme, Nanjing, China) was used for q-PCR with specific primers. GAPDH served as an internal control. The relative mRNA expression levels were calculated using the 2-ΔΔCt method. The specific primers used in this study were provided by Sangon Biotech (Shanghai, China) and are listed in Table 1.

Chemosensitivity assay

The CCK8 assay was employed to determine the impact of cisplatin (MedChemExpress, Monmouth Junction, NJ, USA) and 5-fluorouracil (5-fluorouracil; MedChemExpress, Monmouth Junction, NY, USA) on the viability of HNSCC cells. Cell lines HSC-3, CAL-27, UM1, SCC9, H103, H357, H376, H400 and SCC25 were seeded in 96-well plates at a density of 2  × 10³ cells/well, while HSC-4, HN4, HN12, HN30, and H314 were seeded at 3  × 10³ cells/well. The cells were then treated with varying concentrations of cisplatin or 5-Fluorouracil, followed by incubation for 72 h. Cell viability was subsequently assessed by adding a 10% CCK8 solution (Yeasen Biotechnology, Shanghai, China) to each well, incubating for an additional 2 h, and measuring the absorbance at a wavelength of 450nm using Varioskan LUX (Thermo Scientific, Waltham, MA, USA). The half-maximal inhibitory concentration (IC50) values were then calculated based on the relative survival curves.

Statistical analysis

Differences in continuous variables between the HNSCC and normal tissues were evaluated using the Student’s t-test and Wilcoxon test. Survival curves were compared using the log-rank test in the Kaplan–Meier analysis. Cox regression analysis was utilized to calculate the hazard ratio (HR) and 95% confidence interval (CI) of pyroptosis-related genes and clinical parameters. All statistical analyses were conducted using R software (version 4.1.1; R Core Team, 2020) and GraphPad Prism 9.0.

Individual pyroptosis-related gene expression is associated with poor prognosis of HNSCC

The study design is illustrated in Fig. 1. The expression levels of four pyroptosis-related genes were analyzed in 497 HNSCC and 44 normal tissue samples obtained from the TCGA database. The expression of the four genes (CASP1, CASP3, GSDMD, and GSDME) was found to be significantly elevated in HNSCC tissues (Fig. 2A). Kaplan–Meier (KM) survival analysis of the corresponding genes suggested that higher expression levels of these genes were associated with worse prognosis (CASP3, P = 0.006; GSDME, P = 0.009, Fig. 2B). Moreover, independent prognostic analysis revealed that the odds ratio of these four genes were all greater than 1.00 (GSDME, P = 0.019, Fig. 2C). These findings suggest that the elevated expression of these four pyroptosis-related genes may be associated with worse prognosis of HNSCC.

High PRGsocre may predict unfavorable prognosis of HNSCC

A scoring system, the pyroptosis-related gene score, was created using the GSVA method with the z-score argument to quantify the overall level of pyroptosis. Expression levels of each gene is positively correlated with the PRGscore in both the TCGA and GEO cohorts (correlation coefficient > 0.3 and P < 0.05, Figs. 3A and S1A). Thus, we define this PRGscore as a signature to characterize pyroptosis activity. Although no significant correlation was found between the PRGscore and clinicopathological characteristics (P > 0.05, Table S1). Using the “surv_cutpoint” function of the “survminer” R package, samples were categorized into high and low PRGscore groups (Figs. 3B and S1B). KM survival analysis revealed that patients with a higher PRGscore had significantly worse prognosis (Fig. 3C, P = 0.014), which was verified in GEO cohort (Fig. S1C, P = 0.017). Multivariate Cox regression analysis demonstrated that PRGscore, as a categorical variable, was an independent predictor of inferior prognosis in both the TCGA and GEO cohorts (Fig. 3D, hazard ratio of PRGscore = 1.473, P = 0.007; Fig. S1D, hazard ratio of PRGscore = 1.840, P = 0.034). Meanwhile, the COX proportional hazards assumption was verified, with the result indicating that each covariate has a smooth fitting curve of Schoenfeld residuals changing over time that is almost parallel to the x-axis with p valu e > 0.05 (Fig. S4 and Table S2), thus confirming the applicability multivariate Cox regression analysis in this context. Collectively, these results suggest that a higher PRGscore is associated with worse prognosis in HNSCC.

High PRGscore serves as a protective factor in the cohort of patients receiving chemotherapy

To explore the impact of pyroptosis activity on response to chemotherapy treatment, we investigated whether high pyroptosis levels were a protective factor for prognosis in the chemotherapy cohort and the non-chemotherapy cohort within the TCGA HNSCC dataset. The KM survival analysis and multivariate Cox regression analysis indicated that a higher PRGscore was associated with worse prognosis in the patients not implementing chemotherapy (P = 0.018, Log-rank test, Fig. 4A; hazard ratio of PRGscore = 1.523, P = 0.006, Fig. 4B). However, patients within the higher PRGscore subgroup exhibited better prognosis in the chemotherapy cohort (P = 0.05, Log-rank test, Fig. 4C; hazard ratio of PRGscore = 0.175, P = 0.057, Fig. 4D). Given the observed association between higher levels of pyroptosis and better prognosis in the chemotherapy HNSCC patients, further investigation into the relationship between PRGscore and chemotherapeutic sensitivity is warranted.

High PRGscore may suggest better chemotherapeutic sensitivity in HNSCC

To initially investigate the potential association between PRGscore and chemotherapeutic sensitivity, the half-maximal inhibitory concentration (IC50) of commonly used chemotherapeutic agents for HNSCC was estimated using the pRRophetic algorithm. As hypothesized, our analysis revealed a negative correlation between PRGscore and the IC50 values of chemotherapeutic agents such as docetaxel, gemcitabine, paclitaxel, methotrexate, cisplatin, and 5-fluorouracil in HNSCC. This finding implies that patients with a high PRGscore may show increased sensitivity to chemotherapy (Figs. 5A–5B, P < 0.05).

To further verify the presumed relationship between overall pyroptosis levels and chemotherapy sensitivity, the IC50 values of cisplatin (Fig. 6) and 5-fluorouracil (Fig. 7) were evaluated in multiple HNSCC cell lines, along with the expression levels of the four pyroptosis-related genes (Fig. S2). Utilizing the same analytical method, the PRGscore of each cell line was calculated. As expected, Pearson correlation analysis revealed significantly negative relationship between IC50 of these HNSCC cells for cisplatin or 5-fluorouracil and pyroptosis activity of these cell lines, indicating that HNSCCs with higher PRGscore are more sensitive to chemotherapy (Fig. 5C, R = −0.65, P = 0.011 for cisplatin and Fig. 5D, R = −0.69, P = 0.04 for 5-fluorouracil).

High PRGsocre may be associated with better efficacy of immunotherapy

Furthermore, we conducted correlation analysis of gene expression and immunotherapy to investigate the potential role of PRGscore in predicting patients’ response to immunotherapy. The results of the immune checkpoint analysis demonstrated a positive correlation between PRGscore and several key immune checkpoints, including CD274, CTLA4, PDCD1, IDO1, PDCD1LG2, and HAVCR2. Notably, the high PRGscore group exhibited a tendency toward higher expression levels of these immune checkpoints significantly, as illustrated in Figs. S3A and S3B. Immunophenoscore (IPS) data for patients diagnosed with HNSCC were obtained from the Cancer Immuneome Database (TCIA). The analysis of IPS data revealed that individuals with elevated PRGscore or increased pyroptosis gene expression may benefit more from immunotherapy targeting PD1, CTLA4, or the combination of CTLA4 and PD1, as depicted in Fig. S3C to S3G. These suggested that individuals with higher PRGscore may benefit from specific immunotherapies.