Deciphering a proliferation-essential gene signature based on CRISPR-Cas9 screening to predict prognosis and characterize the immune microenvironment in HNSCC

Abstract

Background

Head and neck squamous cell carcinoma (HNSCC) is a highly aggressive malignancy with a poor prognosis. Identifying reliable prognostic biomarkers and therapeutic targets is crucial for improving patient outcomes. This study aimed to systematically identify proliferation-essential genes (PEGs) associated with HNSCC prognosis using CRISPR-Cas9 screening data.

Methods

CRISPR-Cas9 screening data from the DepMap database were used to identify PEGs in HNSCC cells. A prognostic PEGs signature was constructed using univariate Cox regression, least absolute shrinkage and selection operator (LASSO) Cox regression, and multivariate Cox regression analyses. The predictive accuracy of the signature was validated in internal and external datasets. Weighted gene co-expression network analysis (WGCNA), gene set enrichment analysis (GSEA), and immune infiltration analysis were used to investigate the underlying mechanism between high and low-risk patients. Random forest analysis and functional experiments were conducted to investigate the role of key proliferation essential genes in HNSCC progression.

Results

A total of 1511 PEGs were identified. A seven-gene prognostic PEGs signature (MRPL33, NAT10, PSMC1, PSMD11, RPN2, TAF7, and ZNF335) was developed and validated, demonstrating robust prognostic performance in stratifying HNSCC patients by survival risk. WGCNA and GSEA analyses revealed a marked downregulation of immune-related pathways in high-risk patients. Immune infiltration analysis validated those high-risk patients had reduced immune scores, stromal scores, and ESTIMATE scores, as well as decreased infiltration of multiple immune cell types. Among the identified genes, PSMC1 was highlighted as a pivotal regulator of HNSCC proliferation and migration, as confirmed by functional experiments.

Conclusions

This study identifies a novel PEGs signature that effectively predicts HNSCC prognosis and stratifies patients by survival risk. PSMC1 was identified as a key gene promoting malignant progression, offering potential as a therapeutic target for HNSCC.

Introduction

Head and neck squamous cell carcinoma (HNSCC) is a highly aggressive malignant tumor, ranking as the sixth most common cancer worldwide, with over 890,000 new cases and 450,000 deaths annually [1]. Despite advances in multimodal treatment strategies such as surgery, radiotherapy, chemotherapy, and immunotherapy, the prognosis of HNSCC remains poor, particularly for patients with advanced-stage disease [2, 3]. The disease's high recurrence rate and resistance to conventional treatments further worsen the prognosis, emphasizing the necessity of developing improved therapeutic strategies [4]. Moreover, these challenges underscore the urgent need for innovative approaches to overcome treatment resistance and improve patient outcomes [5, 6]. The major risk factors for HNSCC include tobacco use, alcohol consumption, and human papillomavirus infection, which collectively contribute to its complex etiology and genetic heterogeneity [7]. While immune checkpoint inhibitors (ICIs), including PD- 1/PD-L1 inhibitors, have demonstrated efficacy in enhancing survival, a significant fraction of patients respond inadequately, indicating limited effectiveness in certain subpopulations [8]. This disparity strongly underscores the need for new prognostic biomarkers and therapeutic targets. Due to the heterogeneity of HNSCC, discovering reliable molecular features for prognostic prediction and personalized therapeutic guidance remains a significant challenge [9].

In recent years, the development of genome-wide functional screening technologies, particularly CRISPR-Cas9, has enabled the systematic identification of genes essential for cancer cell proliferation and survival [10]. By inducing site-specific double-strand breaks in DNA, CRISPR-Cas9 enables precise gene knockout and serves as a potent tool for identifying cancer vulnerabilities [11]. Compared to RNA interference (RNAi) screening, CRISPR-Cas9 exhibits higher specificity and efficiency, significantly reducing off-target effects. Large-scale CRISPR-Cas9 screening projects, such as the Cancer Dependency Map (DepMap), have systematically mapped gene dependencies across hundreds of cancer cell lines, providing valuable resources for identifying key genes and potential therapeutic targets [12]. Despite these advancements, the application of CRISPR-Cas9 screening to HNSCC remains relatively underexplored. Using this database, researchers have identified proliferation-essential genes (PEGs) that play a critical role in tumor progression, but their roles in HNSCC remain insufficiently explored [13,14,15,16,17].

In this study, we systematically identified PEGs in HNSCC using the DepMap database and constructed a robust prognostic PEGs signature for predicting overall survival (OS). To develop this prognostic signature, we employed a combination of univariate Cox regression, least absolute shrinkage and selection operator (LASSO) regression and multivariate Cox regression analysis. The identified PEGs signature exhibited significant prognostic potential, enabling the effective classification of HNSCC patients into high-risk and low-risk groups based on OS. Importantly, this prognostic signature provides a novel framework for risk stratification, aiding clinical decision-making. To further elucidate the biological mechanisms underlying this PEGs signature, we performed weighted gene co-expression network analysis (WGCNA), gene ontology (GO) enrichment analysis, and gene set enrichment analysis (GSEA). These analyses revealed functional pathways associated with the identified PEGs. Notably, these pathways were predominantly enriched in cellular processes critical for tumor progression, including cell proliferation and immunosuppression. Furthermore, random forest analysis identified PSMC1 as a critical gene involved in HNSCC progression. Subsequently, we validated the oncogenic role of PSMC1 in HNSCC cells through cellular experiments, providing experimental evidence for its clinical relevance. These results highlight the potential of PSMC1 as a prognostic biomarker and therapeutic target, meriting further exploration.

Materials and methods

Data acquisition

The gene expression profiles and clinical data of HNSCC patients were retrieved from The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/). The external datasets GSE41613 and GSE65858 were sourced from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/). All datasets were preprocessed using standardized methods, including normalization and log2 transformation via the DESeq2 R package (https://bioconductor.org/packages/DESeq2/), ensuring data consistency and quality. The study workflow is illustrated in Fig. 1.

Fig. 1

The flowchart of the whole study

Identification of PEGs in HNSCC cells

To identify proliferation essential genes PEGs in HNSCC, we utilized gene dependency probability scores from the DepMap (https://depmap.org), which estimate the likelihood that a gene behaves as a common essential gene in a given cell line. These scores incorporate corrections for screen quality and experimental variability. Genes with an average dependency probability > 0.5 across HNSCC cell lines were considered PEGs and selected for further prognostic evaluation.

GSEA

To evaluate the functional impact of specific gene knockouts, we used CERES gene effect scores from the DepMap database (https://depmap.org), which quantify the magnitude of gene knockout effects while correcting for confounding factors such as copy number variation and screen quality. Lower CERES scores indicate stronger negative effects on cell viability. All genes were ranked in ascending order based on their average CERES scores in HNSCC cell lines, and GSEA was performed to identify hallmark pathways associated with cell proliferation. Additionally, GSEA was used to compare pathway enrichment between high-risk and low-risk groups based on the PEGs signature. Analyses were conducted using the clusterProfiler R [18] package, with pathways considered significantly enriched at an adjusted p-value < 0.05.

Functional enrichment analyses

We conducted functional enrichment analysis, encompassing GO analysis, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis, and hallmark pathway enrichment analysis, utilizing the R packages clusterProfiler, org.Hs.eg.db (https://bioconductor.org/packages/org.Hs.eg.db/), and enrichplot (https://bioconductor.org/packages/enrichplot/).

Construction and validation of a prognostic signature based on PEGs

We randomly divided 493 patients from the TCGA-HNSC dataset into training and testing cohorts in a 1:1 ratio. In the training cohort, univariate Cox regression analysis was first performed to evaluate the prognostic value of PEGs in HNSCC. Subsequently, LASSO regression (implemented via the glmnet R package) and stepwise multivariate Cox regression analysis were used to identify prognostically relevant genes. The risk score for each patient was calculated using the following formula: Risk Score = (exp gene1 × coef1) + (exp gene2 × coef2) + … + (exp gene n × coefn). In this formula,"exp"represents the mRNA expression level of each gene, while"coef"corresponds to the regression coefficient obtained from multivariate Cox regression analysis. HNSCC patients were classified into high-risk and low-risk groups based on the median risk score. Kaplan–Meier survival analysis was used to assess survival differences between the two groups, and statistical significance was evaluated using the log-rank test. Receiver operating characteristic (ROC) curve analysis was performed using the pROC R package [19], with 1-year, 3-year, and 5-year survival as endpoints to evaluate predictive performance. The area under the ROC curve (AUC), sensitivity, and specificity were used to assess predictive accuracy. To validate the robustness of the PEGs signature, independent validation was performed using two internal datasets (TCGA-test and TCGA-entire) and two external datasets (GSE41613 and GSE65858) to ensure its generalizability and stability.

Association between PEGs signature and clinicopathological characteristics

Visualization of the association between the PEGs signature and different clinicopathological features was performed using the R packages ggplot2 (https://CRAN.R-project.org/package=ggplot2), ggpubr (https://CRAN.R-project.org/package=ggpubr), and pheatmap (https://CRAN.R-project.org/package=pheatmap). Stratified OS analysis was performed using the survival (https://CRAN.R-project.org/package=survival) and survminer (https://CRAN.R-project.org/package=survminer) R packages to assess subgroup-specific effects. Additionally, univariate and multivariate Cox regression analyses were performed using the survival R package to identify independent prognostic factors.

Genomic alterations and tumor mutation burden analysis

Tumor mutation burden (TMB) and mutation profiles in the high-risk and low-risk groups were evaluated using the R packages maftools [20], ggplot2, and forestPlot (https://CRAN.R-project.org/package=forestplot).

Nomogram construction

The prognostic nomogram was developed using the rms R package (https://CRAN.R-project.org/package=rms), incorporating prognostic variables such as clinical N stage, clinical M stage, lymphovascular invasion, perineural invasion, and risk score. The reliability of the nomogram was validated using calibration curves, which assess the correlation between predicted probabilities and actual occurrences to determine predictive accuracy. Predictive performance was further quantified by calculating the concordance index (C-index) and AUC, evaluated using the R packages riskRegression (https://CRAN.R-project.org/package=riskRegression) and timeROC [21] to quantitatively measure model discrimination.

WGCNA

The WGCNA R package [22] was used to construct a weighted gene co-expression network to identify gene modules associated with the PEGs signature. A soft threshold of 12 was applied in the TCGA-HNSC dataset to maintain the scale-free topology of the network. Genes exhibiting similar expression patterns were grouped into respective modules, with the module most correlated with the PEGs signature chosen for further functional enrichment analysis.

Immune infiltration analysis

The ESTIMATE R package [23] was employed to assess immune infiltration and the tumor immune microenvironment by computing immune scores, stromal scores, and ESTIMATE composite scores. Furthermore, single-sample gene set enrichment analysis (ssGSEA) was conducted with the GSVA R package [24] to quantify the infiltration levels of 28 tumor-infiltrating immune cell subsets. Spearman and Pearson correlation analyses were applied to evaluate the association between immune cell infiltration levels and risk scores.

Cell culture

The HNSCC cell lines SCC- 9 and SAS were obtained from the Cell Bank of the Chinese Academy of Sciences (Shanghai, China). SCC- 9 cells were cultured in DMEM/F12 (Gibco, USA), and SAS cells were cultured in DMEM (Gibco, USA). Both media were supplemented with 10% fetal bovine serum (FBS, Gibco, USA) and 1% penicillin–streptomycin solution (P1400, Solarbio, China) at 37 °C in a humidified atmosphere with 5% CO₂.

PSMC1 knockdown

To knock down PSMC1 expression, custom-synthesized small interfering RNAs (siRNAs) sequences were ordered from GenePharma (Suzhou, China). The following siRNA sequences were used: si-NC sense: 5′- UUC UCC GAA CGU GUC ACG UTT − 3′, si-NC antisense: 5′- ACG UGA CAC GUU CGG AGA ATT − 3′, si-PSMC1 - 1 sense: 5′-GAA CCU UGG AAG AGA UCA UTT − 3′, si-PSMC1 - 1 antisense: 5′- AUG AUC UCU UCC AAG GUU CTT − 3′, si-PSMC1 - 2 sense: 5′- GCU GAU GGA UGA CAC GGA UTT − 3′, si-PSMC1 - 2 antisense: 5′- AUC CGU GUC AUC CAU CAG CTT − 3′. For the transfection of PSMC1 siRNA, SCC- 9 and SAS cells were seeded at a concentration of 2 × 10⁵ cells per well in 6-well plates. When the cells reached 70–80% confluence, the siRNAs (si-NC and si-PSMC1) were transfected into SCC- 9 and SAS cells using Lipofectamine 3000 (L3000015, Thermo Fisher, USA) according to the manufacturer's protocol.

Real-time qPCR

To knock down the expression of PSMC1 in SCC- 9 and SAS cells, we used the method described above. Total RNA was then extracted from HNSCC cells using the SevenFast® Total RNA Extraction Kit for Cells (SM130, SevenBio, China) following the manufacturer's instructions. The RNA concentration was measured using the NanoDrop One spectrophotometer (Thermo Fisher, USA). Complementary DNA (cDNA) synthesis was performed using the TransScript® Uni All-in-One First-Strand cDNA Synthesis SuperMix for qPCR (One-Step gDNA Removal) (AU341, Transgen, China) following the manufacturer's guidelines. Real-time qPCR was performed on the Quant Studio™ 6 Real-Time PCR System (Thermo Fisher, USA) using MicroAmp Fast 0.2 mL 96-well Reaction Plates (Thermo Fisher, USA) and SYBR Green Mix (A25742, Thermo Fisher, USA) for the qPCR reactions. GAPDH was used as a housekeeping gene, and PSMC1 expression fold-change relative to the control sample was calculated using the 2^(-∆∆Ct) method. Primer pairs for PCR were obtained from Sangon Biotech. The following primers were used: PSMC1 forward: 5′- CAC ACT CAG TGC CGG TTA AAA − 3′, PSMC1 reverse: 5′- GTA GAC ACG ATG GCA TGA TTG T − 3′, GAPDH forward: 5′- TTG CCA TCA ATG ACC CCT TCA − 3′, GAPDH reverse: 5′- CGC CCC ACT TGA TTT TGG A − 3′.

CCK- 8 assay

Cell viability was evaluated using the CCK- 8 assay. SCC- 9 and SAS cells transfected with si-NC or si-PSMC1 were resuspended in full culture medium and plated into 96-well plates at a density of 2,000 cells per well. The cells were cultured for the specified periods. Subsequently, 10 µl of CCK- 8 solution (BS350 A, Biosharp, China) was added to each well, and the plates were incubated for 2 h at 37 °C in a humidified incubator. The number of viable cells was determined by measuring absorbance at 450 nm using a spectrophotometer. Each experiment was performed with five biological replicates (n = 5) to ensure reproducibility.

Colony formation assay

SCC- 9 and SAS cells transfected with si-control or si-PSMC1 were resuspended in a complete culture medium and seeded into 12-well plates at a density of 2000 cells per well. Cells were cultured for 7 days in a 37 °C incubator with 5% CO₂ and saturated humidity. Afterward, the cells were fixed using a 4% paraformaldehyde solution (P1110, Solarbio, China), followed by staining with 0.1% crystal violet (G1062, Solarbio, China). The total area covered by the colonies was captured using a camera, and the number of colonies was counted using ImageJ software. Each experiment was performed in triplicate (n = 3) to ensure reproducibility.

Migration assays

SCC- 9 and SAS cells transfected with si-control or si-PSMC1 were resuspended in a complete culture medium at a concentration of 500,000 cells/mL, and 70 µL of the cell suspension was pipetted into each chamber of the cell culture insert (80209, Ibidi, Germany). Once the cells adhered, the Ibidi chamber was removed. Cell migration was observed at 0 h, 6 h, and 9 h under an inverted microscope, and images were taken at 10 × magnification. Wound healing was measured using ImageJ software, and the healing rate was calculated as follows: Wound healing rate = (0 h area—n h area)/(0 h area) × 100%.

Statistical analysis

Statistical analyses were conducted using R software (version 4.2.2). For experimental data (e.g., qPCR, CCK- 8, colony formation, and migration assays), GraphPad Prism (version 9.5) was used. Each experimental group included at least three biological replicates, and differences between groups were assessed using two-tailed Student’s t-tests. Statistical significance was denoted as follows: *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001. For survival analysis, univariate Cox regression was performed to evaluate the prognostic value of PEGs, and false discovery rate (FDR) correction was applied using the Benjamini–Hochberg method. Genes with FDR-adjusted p-values < 0.05 were considered statistically significant. For immune infiltration analysis, Spearman correlation was used to assess the association between the risk score and immune cell abundance estimated by multiple deconvolution algorithms. FDR correction was also applied to correlation p-values, and adjusted p < 0.05 was considered significant.

Results

Identification of PEGs in HNSCC cells

We utilized genome-wide CRISPR-mediated loss-of-function screening data from the DepMap database to identify candidate PEGs associated with HNSCC. First, we used CERES scores combined with GSEA analysis to explore key hallmark pathways influencing the proliferation of HNSCC cell lines. As shown in Fig. 2A, hallmark pathways promoting HNSCC cell proliferation include MYC targets, oxidative phosphorylation, E2 F targets, G2M checkpoint, and DNA repair. The top 10 enriched pathways, ranked by normalized enrichment score (NES), are presented in Fig. 2B.

Fig. 2

Identification of PEGs in HNSCC cells. GSEA of hallmark pathways associated with HNSCC cell proliferation based on CERES scores from the DepMap database (a). Enrichment plots of the top 10 hallmark pathways ranked by NES (b). GO enrichment analysis of the identified 1,511 PEGs (c). KEGG pathway enrichment analysis of the identified PEGs (d). Abbreviations: PEGs, proliferation-essential genes; HNSCC, head and neck squamous cell carcinoma; GSEA, gene set enrichment analysis; DepMap, cancer dependency map; NES, normalized enrichment score; GO, gene ontology; KEGG, kyoto encyclopedia of genes and genomes

Next, to determine which genes are essential for HNSCC cell growth, we analyzed gene dependency scores for all HNSCC cell lines (including 72 HNSCC cell lines) in the DepMap database. Genes with an average dependency score greater than 0.5 across all HNSCC cell lines were defined as PEGs. A total of 1,511 PEGs were identified in HNSCC cells (Supplementary Table S1), along with their corresponding CRISPRGeneEffect and CRISPRGeneDependency scores. Furthermore, GO enrichment analysis revealed that these genes were primarily enriched in ribonucleoprotein complex biogenesis, chromosomal region, and catalytic RNA-related activities (Fig. 2C). KEGG enrichment analysis further revealed that the spliceosome pathway was the most significantly enriched (Fig. 2D).

Construction and validation of the PEGs signature

First, all cases in the TCGA-HNSC dataset were randomly divided into a training cohort (n = 247) and a testing cohort (n = 246). Subsequently, key genes were selected using univariate Cox regression analysis and LASSO Cox regression analysis, resulting in 14 genes: ALG2, GAPDH, HSPA5, MRPL33, NAT10, PGK1, PSMC1, PSMD11, PSMD2, RPN2, TAF7, TXNRD1, VBP1, and ZNF335. Further optimization using multivariate Cox regression analysis led to the final selection of seven genes: MRPL33, NAT10, PSMC1, PSMD11, RPN2, TAF7, and ZNF335 (Fig. 3A-C). A risk score model was developed based on these seven genes, and its calculation formula is as follows: Riskscore = 0.44719 × MRPL33 + 0.59750 × NAT10 + 0.49953 × PSMC1 + 0.50911 × PSMD11 + 0.59609 × RPN2 + 0.44624 × TAF7 + (− 0.79222) × ZNF335 (Fig. 3D). Patients were stratified into high-risk and low-risk groups according to the median risk score. The distribution of risk score, survival time, and patient classification in the two groups are depicted in Fig. 3E, whereas Fig. 3F shows the expression patterns of the seven genes across the risk groups. As shown in Fig. 3G, survival analysis revealed that high-risk patients exhibited significantly poorer survival outcomes compared to the low-risk group, highlighting the model's prognostic utility. The robustness of the model was further confirmed by ROC curve analysis, which yielded AUC values of 0.718, 0.794, and 0.766 for 1-year, 3-year, and 5-year OS, respectively (Fig. 3H), demonstrating its high predictive accuracy. The chromosomal positions of the seven genes are illustrated in Fig. 3I.

Fig. 3

Construction of the PEGs signature. LASSO Cox regression analysis for feature selection of prognostic PEGs (a-b). Multivariate Cox regression analysis of the identified prognostic PEGs (c). Risk coefficients of the prognostic PEGs (d). Risk score distribution (top) and survival status of HNSCC patients (bottom) (e). Heatmap displaying the expression profiles of the prognostic PEGs in high-risk and low-risk groups (f). Kaplan–Meier survival analysis comparing overall survival (OS) between high-risk and low-risk patients (g). Receiver operating characteristic (ROC) curve analysis for 1-year, 3-year, and 5-year OS (h). Chromosomal locations of the seven genes in the PEGs signature (i). Abbreviations: PEGs, proliferation-essential genes; LASSO, least absolute shrinkage and selection operator; HNSCC, head and neck squamous cell carcinoma; OS, overall survival; ROC, receiver operating characteristic; AUC, area under the curve. Statistical significance: *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001, ns: not significant

To evaluate the effectiveness and reliability of the PEGs signature, survival analysis, and ROC analysis were conducted in two internal cohorts (TCGA-test and TCGA-entire) and two external cohorts (GSE41613 and GSE65858). The distribution of risk score and survival status of HNSCC patients are shown in Fig. 4A. Kaplan–Meier survival analysis (Fig. 4B) demonstrated that in all internal and external cohorts, patients in the low-risk group had significantly higher survival rates than those in the high-risk group, consistent with the findings from the training cohort. Furthermore, ROC analysis indicated that the AUC values remained high across different cohorts, further emphasizing the robustness of this signature (Fig. 4C).

Fig. 4

Validation of the PEGs signature. Risk score distribution (top) and survival status (bottom) of HNSCC patients in the TCGA-test, TCGA-entire, GSE41613, and GSE65858 cohorts (a). Kaplan–Meier survival analysis comparing OS between high- and low-risk patients in the four cohorts (b). ROC curve analysis for 1-year, 3-year, and 5-year OS in the four cohorts (c). Abbreviations: PEGs, proliferation-essential genes; HNSCC, head and neck squamous cell carcinoma; TCGA, The Cancer Genome Atlas; OS, overall survival; ROC, receiver operating characteristic; AUC, area under the curve. Statistical significance: *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001, ns: not significant

Association between PEGs signature and clinicopathological characteristics

To evaluate the clinical applicability of the PEGs signature, we conducted a comprehensive clinical correlation analysis. Heatmaps and box plots (Fig. 5A-B) demonstrated significant associations between the PEGs signature and various clinicopathological features. The results indicated that patients with advanced clinical T stage (T3/T4), positive perineural invasion, and deceased status had significantly higher risk scores. The stratified analysis further revealed that in different subgroups (stratified by age, gender, clinical stage, lymphovascular invasion, and perineural invasion), patients with high-risk scores consistently exhibited poorer OS (Fig. 6A-P). Importantly, the prognostic value of this signature remained consistent across diverse clinical subgroups, reinforcing its stability.

Fig. 5

Association between the PEGs signature and clinicopathological characteristics. Heatmap displaying the relationship between the PEGs signature and various clinicopathological features, including age, gender, clinical stage (T, N, M), lymphovascular invasion, perineural invasion, and survival status (a). Box plots illustrating the distribution of risk scores across different clinical subgroups (b). Abbreviations: PEGs, proliferation-essential genes

Fig. 6

Stratified survival analysis of the PEGs signature across different clinical subgroups. Kaplan–Meier survival analysis for patients stratified by age (a-b), gender (c-d), clinical T stage (e–f), clinical N stage (g-h), clinical M stage (i-j), clinical stage (k-l), lymphovascular invasion status (m–n), perineural invasion status (o-p). Abbreviations: PEGs, proliferation-essential genes

Genomic alterations and TMB analysis

Analysis of genomic alterations demonstrated that the mutation rate was markedly lower in the low-risk group compared to the high-risk group (89.75% vs. 97.12%) (Fig. 7A-B). Furthermore, the mutation landscape differed significantly between the high- and low-risk groups (Fig. 7C). Among high-risk patients, TP53, SH3PXD2 A, RTN4, FBN1, MIA3, and PLEKHG1 mutations were significantly more frequent (p < 0.01), whereas HRAS, CYLD, HFM1, NIPBL, ROBO2, and RIMS2 mutations were more prevalent in low-risk patients (p < 0.05). Beyond mutation profiling, we also assessed differences in TMB between the high- and low-risk groups (Fig. 7D-F). Notably, TMB was significantly elevated in HNSCC patients within the high-risk group (Fig. 7D). Kaplan–Meier survival analysis indicated that patients with high TMB had lower survival rates (Fig. 7E). Moreover, HNSCC patients characterized by both high TMB and high-risk scores had the poorest survival outcomes (Fig. 7F).

Fig. 7

Genomic alterations and TMB analysis in high- and low-risk HNSCC patients. Mutation landscape of the low-risk (a) and high-risk (b) groups, showing the frequency and distribution of somatic mutations in HNSCC patients. Comparison of significantly mutated genes between the low-risk and high-risk groups (c). Violin plot comparing TMB levels between the low-risk and high-risk groups (d). Kaplan–Meier survival analysis of HNSCC patients stratified by TMB (e). Kaplan–Meier survival analysis of HNSCC patients stratified by both TMB status and risk score (f). Abbreviations: TMB, tumor mutation burden; HNSCC, head and neck squamous cell carcinoma; OR, odds ratio; CI, confidence interval. Statistical significance: *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001, ns: not significant

Nomogram construction

To assess whether our signature is an independent prognostic factor for HNSCC patients, we performed univariate and multivariate Cox regression analyses. The results indicated that both risk score and perineural invasion were independent prognostic factors (Fig. 8A-B), suggesting that the risk score can predict HNSCC prognosis independently of other clinical variables. To enhance the clinical applicability of our study, we constructed a nomogram model based on significant prognostic factors, including clinical N stage, clinical M stage, lymphovascular invasion, perineural invasion, and risk score (Fig. 8C). This nomogram provides an intuitive and practical statistical tool for predicting survival in HNSCC patients, facilitating the clinical application of the risk score. The results demonstrated that this nomogram exhibited strong predictive performance for 1-year, 3-year, and 5-year OS in HNSCC patients (Fig. 8D). The C-index of the nomogram was significantly higher than that of individual predictors (Fig. 8E), highlighting the superiority of the integrated model. Additionally, Kaplan–Meier survival analysis showed that patients with higher nomogram scores had significantly worse OS (Fig. 8F). The ROC AUC values for 1-year, 3-year, and 5-year OS predictions using the nomogram were 0.714, 0.778, and 0.724, respectively (Fig. 8G), further emphasizing the robustness and clinical significance of this prognostic nomogram.

Fig. 8

Nomogram construction and validation. Univariate (a) and multivariate (b) Cox regression analyses of prognostic factors, including clinical and the PEGs signature. Nomogram incorporating significant prognostic factors for predicting 1-year, 3-year, and 5-year OS (c). Calibration curves assessing the predictive accuracy of the nomogram for 1-year, 3-year, and 5-year OS (d). Time-dependent C-index comparison between the nomogram and individual prognostic factors (e). Kaplan–Meier survival analysis based on nomogram-derived risk stratification (f). ROC curve analysis for 1-year, 3-year, and 5-year OS predictions using the nomogram (g). Abbreviations: OS, overall survival; Concordance index, C-index; HNSCC, head and neck squamous cell carcinoma; ROC, receiver operating characteristic; AUC, area under the curve

WGCNA and GSEA

We constructed a WGCNA network to explore regulatory relationships associated with the PEGs signature. In the TCGA-HNSC dataset, a soft threshold of β = 12 was applied, and 12 distinct gene modules were identified based on the topological overlap matrix (TOM) (Fig. 9A-B). Among these, the yellow module exhibited the strongest correlation with the risk score (Fig. 9B). Further analysis revealed that the yellow module had the strongest negative correlation with the risk score (correlation coefficient = − 0.53, p = 6.9e- 106) (Fig. 9C). GO enrichment analysis indicated that the genes in this module were primarily associated with leukocyte-mediated immunity, the external side of the plasma membrane, and immune receptor activity (Fig. 9D). Additionally, KEGG enrichment analysis revealed significant enrichment in cytokine-receptor interactions, chemokine signaling pathways, and cell adhesion molecule pathways (Fig. 9E). Hallmark pathway enrichment analysis further revealed that these genes were predominantly involved in immune-related processes, including allograft rejection, inflammatory response, interferon-gamma response, complement activation, IL2-STAT5 signaling, and IL6-JAK-STAT3 signaling (Fig. 9F).

Fig. 9

WGCNA analyses of the PEGs signature. Cluster dendrogram of genes based on WGCNA (a). Module-trait relationship heatmap showing the correlation between gene modules and risk groups (b). Scatter plot illustrating the correlation between gene significance and module membership in the yellow module (c). GO enrichment analysis of genes in the yellow module (d). KEGG pathway enrichment analysis of genes in the yellow module (e). Hallmark pathway enrichment analysis of genes in the yellow module (f). Abbreviations: WGCNA, weighted gene co-expression network analysis; PEGs, proliferation-essential genes; GO, gene ontology; KEGG, Kyoto encyclopedia of genes and genomes

To further elucidate the potential biological mechanisms underlying survival differences between high- and low-risk groups driven by the PEGs signature, we conducted a GSEA analysis. Hallmark pathway enrichment analysis indicated that in the TCGA dataset, the PEGs signature was associated with signaling pathways such as oxidative phosphorylation, MYC targets, mTORC1 signaling, glycolysis, allograft rejection, IL6-JAK-STAT3 signaling, interferon response, IL2-STAT5 signaling, and inflammatory response (Fig. 10A-C). Notably, these pathways significantly overlapped with those enriched in the yellow module, particularly immune-related pathways. Furthermore, immune-related pathways—including allograft rejection, IL6-JAK-STAT3 signaling, interferon-gamma response, interferon-alpha response, IL2-STAT5 signaling, inflammatory response, and complement pathways—were significantly downregulated in the high-risk group (Fig. 10C).

Fig. 10

GSEA of hallmark pathways associated with the PEGs signature. Bubble plot showing hallmark pathways significantly enriched in the high-risk and low-risk groups based on GSEA (a). Enrichment plots of key hallmark pathways activated in the high-risk group (b). Enrichment plots of immune-related hallmark pathways suppressed in the high-risk group (c). Abbreviations: GSEA, gene set enrichment analysis; PEGs, proliferation-essential genes

Immune infiltration analysis

To explore immune infiltration and the tumor immune microenvironment in HNSCC, we initially assessed immune scores, stromal scores, and ESTIMATE scores in high- and low-risk patient groups. Our findings indicated that immune, stromal, and ESTIMATE scores were markedly reduced in high-risk HNSCC patients (Fig. 11A). Furthermore, ssGSEA was employed to examine the correlation between risk score and tumor-infiltrating immune cells. The analysis revealed that most immune cell infiltration levels were markedly decreased in high-risk patients compared to the low-risk group (Fig. 11B). Specifically, 19 types of immune cells exhibited significantly lower infiltration levels in the high-risk group (Fig. 11C). Additional correlation analysis revealed a negative association between the risk score and the majority of immune cell types (Fig. 12A). A more detailed correlation analysis was performed for immune cell types that exhibited a significant negative correlation with the risk score (cor < − 0.1, p < 0.05) (Fig. 12B-P).

Fig. 11

Immune infiltration analysis in high- and low-risk HNSCC patients. Box plots comparing ESTIMATE score, immune score, and stromal score between high- and low-risk groups (a). Heatmap displaying the infiltration levels of tumor-infiltrating immune cells in high- and low-risk groups based on ssGSEA (b). Violin plots showing the distribution of immune cell infiltration levels between high- and low-risk groups (c). Abbreviations: HNSCC, head and neck squamous cell carcinoma; ssGSEA, single-sample gene set enrichment analysis. Statistical significance: *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001, ns: not significant

Fig. 12

Correlation analysis between the PEGs signature and immune cell infiltration. Bubble plot showing the correlation between risk score and various tumor-infiltrating immune cells (a). Scatter plots depicting the detailed correlation between risk score and selected immune cell types with significant negative associations (b). Abbreviations: PEGs, proliferation-essential genes

To further validate the immune landscape differences, we performed multi-algorithm immune cell deconvolution using CIBERSORT, TIMER, XCELL, EPIC, MCPCOUNTER, and QUANTISEQ. The results consistently revealed that high-risk patients exhibit decreased infiltration of anti-tumor immune cells and increased infiltration of immunosuppressive cells, supporting an immunosuppressive tumor microenvironment (Fig.S1).

Identification of PSMC1 as the most important gene in PEGs signature

Next, random forest analysis was performed to determine the key gene within the PEGs signature for further study. The findings demonstrated that PSMC1 was the most pivotal gene among the seven PEGs (Fig. 13A-C). Pan-cancer analysis of the TCGA database showed that PSMC1 exhibited elevated expression across multiple tumor types compared to normal tissues (Fig. 13D). Importantly, PSMC1 expression was markedly elevated in HNSCC tissues (p < 0.0001, Fig. 13D-E). ROC curve analysis indicated that PSMC1 expression effectively differentiated HNSCC from normal tissues, with an AUC of 0.812 (95% confidence interval: 0.763–0.846), implying that PSMC1 could be a potential diagnostic biomarker for HNSCC (Fig. 13F). Kaplan–Meier survival analysis demonstrated that high PSMC1 expression was associated with significantly poorer OS in HNSCC patients (Fig. 13G-H). Furthermore, CRISPR-based gene effect analysis (CERES score) indicated that PSMC1 depletion markedly suppressed the proliferation of all HNSCC cell lines (Fig. 13G). These results suggest that PSMC1 is a key player in HNSCC progression and could be a promising therapeutic target.

Fig. 13

Identification of PSMC1 as the key gene in the PEGs signature. Random forest analysis ranking the importance of the seven PEGs, identifying PSMC1 as the most pivotal gene (a-c). Pan-cancer analysis of PSMC1 expression across multiple tumor types in the TCGA database (d). Box plot comparing PSMC1 expression between tumor and para-tumor tissues in HNSCC (e). ROC curve analysis evaluating the diagnostic performance of PSMC1 in distinguishing HNSCC from para-tumor tissues (f). Kaplan–Meier survival analysis showing the association between PSMC1 expression and OS in the TCGA-HNSC (g) and GSE41613 (h) cohorts. CERES scores of PSMC1 in all HNSCC cell lines (i). Abbreviations: PEGs, proliferation-essential genes; HNSCC, head and neck squamous cell carcinoma; TCGA, The Cancer Genome Atlas; ROC, receiver operating characteristic; OS, overall survival. Statistical significance: *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001, ns: not significant

Furthermore, to explore the potential mechanisms underlying the oncogenic and immunosuppressive roles of PSMC1, we performed KEGG pathway enrichment analysis based on genes associated with PSMC1 expression. The results revealed that high PSMC1 expression was negatively associated with multiple immune-related pathways, including primary immunodeficiency, allograft rejection, antigen processing and presentation, T cell receptor signaling pathway, and intestinal immune network for IgA production (Fig. S2). These findings suggest that PSMC1 may contribute to the formation of an immunosuppressive tumor microenvironment in HNSCC, further supporting its role in cancer progression.

Knockdown PSMC1 inhibits the proliferation and migration of HNSCC cells

To investigate the functional role of PSMC1, specific PSMC1-siRNA was transfected into HNSCC cells (SAS and SCC- 9), effectively silencing PSMC1 expression (Fig. 14A). The CCK- 8 assay demonstrated that silencing PSMC1 markedly suppressed HNSCC cell proliferation (Fig. 14B-C). The colony formation assay further verified that PSMC1 knockdown substantially impaired the clonogenic potential of HNSCC cells (Fig. 14D-E). Additionally, the scratch wound assay demonstrated that PSMC1 silencing significantly suppressed the migration capability of HNSCC cells (Fig. 14F-G). In summary, these results suggest that PSMC1 is a key regulator of HNSCC cell proliferation and migration, and targeting PSMC1 may offer a promising therapeutic approach for HNSCC.

Fig. 14

Knockdown of PSMC1 inhibits the proliferation and migration of HNSCC cells. qRT-PCR analysis confirming the knockdown efficiency of PSMC1-siRNA in SAS and SCC- 9 cells (a). CCK- 8 assay evaluating the proliferation of SAS (b) and SCC- 9 (c) cells after PSMC1 knockdown. Colony formation assay assessing the clonogenic potential of SAS (d) and SCC- 9 (e) cells following PSMC1 knockdown. Scratch wound healing assay measuring the migration rate of SAS (f) and SCC- 9 (g) cells after PSMC1 knockdown. Abbreviations: HNSCC, head and neck squamous cell carcinoma; qRT-PCR, quantitative real-time polymerase chain reaction; CCK- 8, cell counting kit- 8; siRNA, small interfering RNA. Statistical significance: *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001, ns: not significant

Discussion

HNSCC remains a highly aggressive malignancy with a poor prognosis, highlighting the urgent need for novel prognostic biomarkers and therapeutic targets [25]. Despite advances in multimodal treatment strategies, the survival rate of HNSCC patients remains unsatisfactory, particularly for those diagnosed at advanced stages [2, 3]. In recent years, CRISPR-Cas9 screening has emerged as a transformative tool for genome-wide functional studies, enabling the systematic identification of essential genes and pathways. This technology has proven particularly powerful in uncovering therapeutic vulnerabilities in cancers, offering new opportunities for the development of targeted therapies and precision oncology approaches [12]. Recent studies have highlighted the critical role of PEGs in tumor progression [13,14,15,16,17]; however, their specific contributions to HNSCC remain largely unexplored. Therefore, this study aimed to systematically identify PEGs associated with HNSCC prognosis and investigate their biological significance.

Based on CRISPR-Cas9 screening data from the DepMap database, we identified 1,511 PEGs in HNSCC cells. By employing a stringent selection strategy, incorporating univariate Cox regression, LASSO Cox regression, and multivariate Cox regression analyses, we developed a prognostic PEGs signature consisting of MRPL33, NAT10, PSMC1, PSMD11, RPN2, TAF7, and ZNF335. The PEGs signature successfully classified HNSCC patients into low-risk and high-risk groups, with the high-risk group exhibiting significantly poorer survival outcomes. Moreover, the PEGs signature was independently validated across internal and external datasets, reinforcing its potential clinical utility.

The clinical relevance analysis demonstrated that an increased risk score was strongly correlated with poor clinical characteristics, such as advanced T stage (T3/T4), positive perineural invasion, and elevated mortality. Additionally, genomic alteration analysis showed that the high-risk group exhibited significantly increased TMB and mutation rates compared to the low-risk group. Importantly, patients characterized by both high TMB and high-risk scores exhibited the poorest survival outcomes. Moreover, subgroup analysis validated that this signature maintained its prognostic significance in various clinical subgroups, reinforcing its stability. To enhance clinical applicability, we developed a nomogram by incorporating this PEGs signature with other prognostic factors, which exhibited high predictive accuracy for 1-year, 3-year, and 5-year OS in HNSCC patients. From a translational perspective, the PEGs-based prognostic model has the potential to be implemented in clinical settings using a standardized qRT-PCR panel. Given the limited number of genes involved, qRT-PCR offers a feasible and cost-effective platform for evaluating PEGs expression in tumor biopsy samples. This approach may enable personalized risk stratification and inform treatment planning in HNSCC patients, thereby bridging the gap between molecular profiling and clinical decision-making.

Several prognostic models have been reported for HNSCC, including those based on TP53 mutation status, ferroptosis-related gene signatures, and long non-coding RNAs (lncRNAs) [26,27,28]. For example, a TP53 mutation-associated model constructed from ten differentially expressed genes was shown to stratify patients into high- and low-risk groups with significantly different survival outcomes and predicted responses to immunotherapy and chemotherapy [26]. However, the prognostic value of TP53-related models may be limited by the heterogeneity and complexity of TP53 mutations across HNSCC subtypes. Another study developed a ferroptosis-related prognostic score by identifying molecular subtypes and constructing an 8-gene signature. Although the FPRS model showed good predictive performance and was linked to immune infiltration and immune escape, it was based solely on transcriptomic correlations and lacked direct functional validation [27]. In addition, a nomogram based on an 8-lncRNA signature was developed and validated using TCGA data and qRT-PCR in 102 clinical specimens. This model achieved a 3-year survival AUC of 0.74 and was proposed as a potential diagnostic and prognostic tool. However, the biological roles of many lncRNAs in tumor proliferation remain poorly understood, which may limit the mechanistic interpretability of such models [28]. In contrast, our PEGs-based signature was derived from genome-scale CRISPR-Cas9 functional screening data, which directly reflects gene essentiality for tumor cell fitness. This functionally grounded approach provides robust biological relevance beyond correlation-based models. Furthermore, our model demonstrated consistent and high prognostic accuracy across multiple independent cohorts, suggesting strong generalizability. Therefore, the PEGs signature may serve as a complementary and potentially more robust prognostic tool for HNSCC, with both mechanistic and translational implications.

To elucidate the underlying mechanisms contributing to the prognostic disparity between high- and low-risk patients, we performed WGCNA and GSEA analyses, which identified a marked downregulation of immune-related pathways in high-risk patients. Further immune infiltration analysis validated those high-risk patients had reduced immune scores, stromal scores, and ESTIMATE scores, as well as decreased infiltration of multiple immune cell types. These findings suggest that immune suppression may contribute to the poorer prognosis observed in high-risk patients, highlighting the potential interplay between the PEGs signature and the tumor immune microenvironment.

PSMC1, identified through random forest analysis, emerged as the most critical gene among the seven PEGs in our prognostic signature. As a key ATPase subunit of the 26S proteasome, PSMC1 plays a crucial role in protein degradation and cellular homeostasis [29, 30]. In various cancers, PSMC1 is frequently upregulated, promoting tumor progression by enhancing cell proliferation, immune evasion, and drug resistance [31,32,33,34]. Due to its association with poor prognosis, PSMC1 is considered a potential therapeutic target and prognostic biomarker in cancer research. In HNSCC, our findings suggest that PSMC1 is not only essential for tumor cell proliferation but also strongly associated with poor patient outcomes, as demonstrated by Kaplan–Meier survival analysis. The significant correlation between high PSMC1 expression and shorter OS highlights its potential as both a prognostic biomarker and a therapeutic target. Further in vitro functional experiments demonstrated that PSMC1 knockdown significantly inhibited HNSCC cell proliferation and migration, reinforcing its role in tumor progression. Collectively, these findings position PSMC1 as a promising candidate for further investigation in HNSCC biology and therapy. In addition to its role in promoting tumor proliferation and poor prognosis, PSMC1 may also play a part in modulating the tumor immune microenvironment. KEGG pathway enrichment analysis based on PSMC1 expression revealed significant negative associations with several immune-related pathways, such as antigen processing and presentation, T cell receptor signaling, and intestinal immune responses. These results imply that PSMC1 may contribute to immune evasion in HNSCC, potentially through the suppression of key immune signaling pathways. This immunosuppressive function, combined with its tumor-promoting effects, highlights the multifaceted role of PSMC1 in HNSCC and supports its potential as a therapeutic target worthy of further investigation.

Despite the promising findings, this study has several limitations. First, although the oncogenic role of PSMC1 was validated in vitro, further in vivo studies are needed to confirm its tumorigenic potential. Second, while our prognostic signature was validated using retrospective datasets, prospective clinical validation with well-annotated patient cohorts, including clinical and follow-up data, is required to confirm its clinical utility. We plan to collaborate with clinical centers to address this in future studies. Third, although we identified an association between the PEGs signature and immune infiltration, the underlying molecular mechanisms remain to be elucidated. Future studies should explore the precise role of PSMC1 in HNSCC, integrate multi-omics approaches, and conduct prospective clinical trials to further validate the prognostic significance of the PEGs signature in precision oncology.

Conclusions

In conclusion, this study systematically identified a PEGs signature (MRPL33, NAT10, PSMC1, PSMD11, RPN2, TAF7, and ZNF335) in HNSCC using CRISPR-Cas9 screening data from the DepMap database. The PEGs signature demonstrated strong prognostic value, effectively stratifying HNSCC patients into high- and low-risk groups and providing a novel tool for risk assessment. Furthermore, PSMC1 was identified as a key oncogenic driver in HNSCC, playing a crucial role in tumor progression. Functional validation experiments confirmed that PSMC1 knockdown significantly inhibited HNSCC cell proliferation and migration, highlighting its potential as both a prognostic biomarker and a therapeutic target. These findings offer new insights into the molecular mechanisms underlying HNSCC progression and provide potential avenues for personalized treatment strategies.