The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/) included the gene expression profiling, somatic mutation, and clinical phenotype data of HNSCC. The RNA-seq data were then log2 transformed and converted to TPM format. All patients were assured to have a survival time longer than 0 days, and samples with missing survival time and survival status were eliminated when processing the TCGA-HNSCC data. Screening produced 499 HNSCC samples and 44 control samples. Furthermore, the GSE41613 data was collected from Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) database. A sum of 97 tumor samples were acquired by selecting the GEO cohort with survival time, converting the probes to Symbol according to the annotation file, and excluding samples lacking clinical follow-up information and overall generation rate data. The set of NAD⁺-related genes was then acquired from the MSigDB database (https://www.gsea-msigdb.org/gsea/msigdb/index.jsp), which contained the reactome database (R-HSA-196807) and KEGG pathway database (pathway: hsa00760).

Using ssGSEA in the “GSVA” package, we determined the NMRGs correlation scores for every TCGA-HNSCC sample (18). We next used the “WGCNA” package to create weighted gene co-expression network to find co-expression networks and select genes from various clusters (19, 20). All samples and missing genes are first clustered. Second, the optimal soft threshold power (β=16) is found using the “pickSoftThreshold” R function to identify significant correlations between modules more effectively. Then, using the requirement of at least 60 genes per module (minModuleSize = 60) to identify gene modules, we conducted a hierarchical cluster analysis. Last but not least, we employed the R package "Heatmap" (21) to extract various module signature genes according to the first principal component of module expression. Then, we assessed the relationship between the module genes and the diagnosis of clinical signature to test the association between module and signature scores. The genes contained in the modules were extracted by filtering the modules with the highest correlations (22).

In the TCGA cohort, the “limma” package was utilized to find DEGs between HNSCC and control samples (23). Using p-adj < 0.05 and |log 2(FC)|>log2(1.5), the gene expression profile was professionally summarized, quartile normalized, and background adjusted in order to screen for significant DEGs. Following their intersection with the DEGs, the midnightblue and green modular genes were found and examined using the R package “clusterProfiler” (24) to examine the module genes’ KEGG function and gene ontology (GO) (the screening criteria was p-value < 0.05). To assess the modular gene enrichment pathways and biological processes, respectively, we created bubble diagrams by charting the top 10 functions enriched to the three terms of the GO enrichment results and the top 10 enriched KEGG pathway results. We used the R package “clusterProfiler” to compute the GSEA of the high- and low-risk groups of TCGA-HNSCC to look into the pathways of various biological processes in various subgroups. The KEGG database was used the reference for enriched pathways during analysis (25).

To find genes significantly linked with prognosis in TCGA-HNSCC patients (p<0.05), the R package “survival” (26) was used to conduct univariate Cox proportional risk regression on intersecting genes. To enhance the model generalization, 10-fold cross-validation was employed and LASSO Cox regression analysis of the R package "glmnet" (27) to compress the genes in order to maximize gene number in the risk model. Furthermore, multifactorial stepwise regression analysis was employed to check for important genes and correlation coefficients that were independently linked to the prognosis of HNSCC, and risk scores were computed for every patient. The following is the formula: Riskscore=Σβi×Expi, where Expi is the expression of each gene gathered, i is the gene expression level, and β is the associated gene’s Cox regression coefficient. Following zscore normalization, the Riskscore was used to assign the TCGA-HNSCC patients into high- and low-risk groups by the Riskscore’s optimal critical value. The R package “survminer” (28) was then utilized to conduct survival analysis between the low- and high-risk groups. Kaplan-Meier (KM) survival curves were then displayed for prognosis analysis, followed by using log-rank test to evaluate significant differences. Further, we examined the prediction of the model by displaying time-dependent receiver operating characteristic (ROC) curves using the R package “timeROC” (29) and calculated 1-, 2-, 3-, 4- and 5-year area under the curve (AUC). Finally, we validated the GSE41613 dataset using the same methodology to better validate the stability and reliability of our constructed clinical prognostic model based on risk-related gene signatures.

Since genomic mutations are closely associated with disease onset and progression (30), we analyzed each sample in the TCGA-HNSCC cohort for gene mutations. The mutation dataset of HNSCC samples was processed using MuTect 2 software (31), and the mutation profiles of the top 10 mutated genes in the low- and high-risk groups were plotted separately.

The association between Riskscore and immune function in HNSCC was evaluated by analyzing the immune infiltration of the TCGA dataset samples using the R package “estimate” (32) and expressed as their respective scores (StromalScore, ImmuneScoreh and ESTIMATEScore). The association between the Riskscore of the TCGA dataset and the 10 immune cell scores was calculated using the “MCPcounter” package (33). According to the transcriptomic expression profiles of the samples, we computed the scores of 28 tumor-infiltrating immune cells (34) in the TCGA-HNSCC cohort with the ssGSEA function of the R package "GSVA".

MEM medium (Gibco, USA) and DMEM medium (Gibco, Grand Island, NY, USA) were used to culture human pharyngeal squamous carcinoma cell line (FaDu) and human normal squamous epithelial cell line (NOK) ordered from Procell Life Science and Technology Co. Ltd (Wuhan, China), respectively. All cultures were added with 1% penicillin-streptomycin (Solarbio, Beijing, China) and 10% fetal bovine serum (Clark, Richmond, VA, USA). All the cells were cultured in an incubator with 5% CO₂ at 37°C. Utilizing Lipofectamine 3000 Transfection Reagent (Thermo Fisher Scientific, Waltham, MA, USA), si-PSME1 and negative control si-RNA were transiently transfected. si-PSME1 sequences were as follows, sense. UGGAUUUGUACCAUUCUUCUG, antisense: GAAGAAUGGUACAAAUCCAAG.

Total RNA from NOK and FaDu cells was separated applying RNA Extraction Kit (TRIzol, Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s protocols. The purity and concentration of the total RNA were assessed, and cDNA templates were generated using the HiScript II kit (Vazyme, Nanjing, China). Quantitative real-time PCR (qRT-PCR) was conducted using specific primers and the KAPA SYBR^® FAST kit (Sigma Aldrich, St Louis, MO, USA). GAPDH was an internal control, and the 2^-ΔΔCT method was used for data analysis (35). Table 1 shows the primer sequences of the specific genes.

Scratch and transwell assays were subsequently carried out to examine the effect of PSME1 expression on FaDu cell migration and invasion. For migration assays, collective cell migration was detected in a wound healing assay. Transfected cells were inoculated into 6-well plates (5 × 10⁵/ml). 2 ml of cell suspension was inoculated into 6-well plates and incubated with 5% CO₂ in an incubator at 37°C. When the cell density was approximately 80%, the monolayer was scraped with a 10 μL plastic pipette tip to create a uniform wound. PBS was used to wash the monolayers, which were then incubated in a non-FBS medium. The wound edge distances between two edges of migrating cell sheet were imaged at 0 h and 48 h, respectively. All the experiments were conducted three times.

For the invasion assay, 1 × 10⁵ cells were inoculated into the upper chamber covered with 10% Matrigel (Corning, Inc., Corning, NY, USA) for 24-h incubation. After the incubation, cells in the upper chamber were eliminated by swabbing, while those on the lower chamber were then fixed by 4% paraformaldehyde and dyed by 0.1% crystal violet solution. These migrated or invaded cells in the lower chamber were counted under a microscope using 6 different fields of view (36).

All statistical analyses were performed using Prism 8 (GraphPad Software, San Diego, CA, USA) and R software version 3.6.0 ((R Foundation, Vienna, Austria)). Wilcoxon rank-sum test was utilized to calculate the difference between the two groups of continuous variables. Correlations were calculated using the spearman method, and the log-rank test was employed to compare the survival between patients in each subgroup. p<0.05 was defined to be statistically different.

Next, we used the ssGSEA method to determine each sample’s NMRG score in the TCGA dataset. NMRG-related gene modules were identified using the R package “WGCNA”. To satisfy the scale-free topology of the network, we selected a soft threshold power of 16 to construct the topological network ( Figure 1A ). 9 co-expression modules were ultimately produced when the module correlation was computed and the module contained a minimum of 60 genes ( Figures 1B, C ). Out of the 9 modules, the grey module had a comparatively large number of genes, followed by the salmon module, as seen in Figure 1D . The gene that was unable to aggregate to other modules was known as the grey module. To select clinically important modules, we calculated the correlation of each module with NMRG scores and plotted a heat map of module-shape correlation. Among the nine modules, significant strong positive correlations were found between midnightblue and green modules and NMRG scores (midnightblue: cor = 0.7, p = 7.87e-76; green: cor = 0.49, p = 1.25e-31 and Figure 1E ).

Figure 2A shows the volcano map of DEGs. Next, we found DEGs between HNSCC cases and control samples in the TCGA dataset. In total, we found 1206 strongly down-regulated genes and 4954 significantly up-regulated genes. After that, we found 359 intersecting genes by taking the genes of the midnightblue and green modules and DEGs ( Figure 2B ). We used GO and KEGG to enrich the intersecting genes in order to investigate their regulatory involvement in the pathophysiology of HNSCC. According to KEGG analysis, pathways such as epstein-barr virus infection and cytokine-cytokine receptor interaction were significantly enriched with the intersecting genes ( Figure 2C ). GO enrichment analysis showed that the BPs in which the intersecting genes were largely involved were defense pathways, for instance, T cell activation and defense response to another organism ( Figure 2D ). The CCs that were mainly localized were side of the membrane, secretory granule membrane, endocytic vesicle and other structures ( Figure 2E ). The MFs most significantly enriched for intersecting genes were pathways such as cytokine receptor activity and chemokine activity ( Figure 2F ).

To create the risk model, we split the TCGA-HNSCC samples into training and test sets in a 5:5 ratio. For LASSO, its effectiveness hinges on the sparsity assumption, which may not hold in the complex genetic landscape of TCGA-HNSC, potentially causing misidentification of prognostic genes and coefficient underestimation (37). Moreover, the data-dependent variable selection process increases the risk of overfitting. As for stepwise regression, the arbitrary selection criteria based on statistical significance, the order-dependence of variable entry/removal, and the susceptibility to overfitting when the gene to sample size ratio is high can all lead to suboptimal model performance and biased results (38). In order to remove redundant confounding genes and identify the genes that have the biggest influence on patients’ prognosis, the R package “survival” was employed to run univariate Cox proportional risk regression on the aforementioned intersecting genes using the training set data. To minimize the gene range in the risk model, we compressed these genes using the “glmnet” package’s LASSO Cox regression approach. To enhance the model generalization, we further performed 10-fold cross-validation ( Figure 3A ). We then used multifactorial stepwise regression analysis to determine which six distinctive genes (ICOS, PSME1, SERPINA1, SH3KBP1, SP100 and ZAP70) were independently linked to prognosis ( Figure 3B ). Characteristics indicative of the prognostic outcomes in the TCGA-HNSCC training set were developed based on the expressions of the characterized genes and the regression coefficients as described below: Riskscore = (-0.685*ICOS) +0.749*PSME1 + 0.25*SERPINA1 + 0.41*SH3KBP1 + 0.24*SP100+(-0.368*ZAP70). Based on the best critical value of the Riskscore, TCGA training set patients were classified into low-risk and high-risk groups. KM curves showed that compared with the high-risk group, patients in the TCGA-HNSCC training set (p< 0.0001), validation set (p< 0.0019), and the low-risk group of the TCGA cohort (p< 0.0001) had better overall survival (OS) ( Figures 3C–E ). To investigate the diagnostic accuracy of the prognostic risk model, ROC analyses for 1-, 2-, 3-, 4- and 5-year prognostic predictions were conducted using the “timeROC” R package. The results indicated that the TCGA-HNSCC training set, test set, and TCGA cohort displayed high AUC values at 1-, 2-, 3-, 4- and 5-years (training set: 0.77, 0.78, 0.78, 0.75 and 0.7; test set: 0.59, 0.65, 0.59, 0.62 and 0.52; TCGA cohort: 0.68, 0.71, 0.69, 0.68 and 0.62, Figures 3F–H ), demonstrating good classification accuracy for prognosis features. Subsequently, we analyzed the expressions of the characterized genes between patients in the low- and high-risk groups of the TCGA cohort. ICOS and ZAP70 were low-expressed in the high-risk group than the low-risk group, whereas PSME1, SERPINA1, SH3KBP1 and SP100 had markedly higher expression in the high-risk group than the low-risk group ( Figure 3I ).

We utilized the GSE41613 dataset to assess the model robustness using comparable models and equivalence coefficients to those utilized in the training set to confirm the stability and dependability of our developed model using NMRG-related signature genes. The training set finding that the prognostic outcomes of high-risk HNSCC patients were more unfavorable (p<0.0001, Figure 3J ) was supported by the validation data. For the 1-, 3- and 5-year periods, the GSE41613 validation set’s AUC values were 0.72, 0.73 and 0.74, respectively Figure 3K ). For patients with HNSCC, it showed that the prognostic model had good prognostic prediction.

We further analyzed the gene mutations in the low- and high-risk groups. In the TCGA-HNSCC samples, we observed that in the high-risk group, 141 (94%) HNSCC patients out of 150 samples showed high-frequency mutations top 10 genes, of which the top 3 mutated genes were TP53 (73%), TTN (39%) and FAT1 (26%) ( Figure 4A ). In contrast, the top 10 mutated genes in HNSCC patients showed mutations in 315 out of 344 samples (91.57%) in the low-risk group, with TP53 as the gene with the most mutations of 65% in the samples, followed by TTN and FAT1 in 39% and 21% respectively ( Figure 4B ). Subsequently, we revealed the mutational co-occurrence or mutual exclusion patterns between the top 10 mutant gene pairs in different risk subgroups of HNSCC, respectively. Statistical results showed frequent co-occurrence between the TNN gene with MUC16 gene, CDKN2A gene with TP53 gene and FAT1 genes in the high-risk group (p < 0.05, Figure 4C ). Significant co-occurrence was also shown in the low-risk group, especially between TTN and SYNE1, LRP1B, MUC16 and CSMD3 (p < 0.05, Figure 4D ). This may imply that the co-occurring genes have a synergistic role in the associated pathological processes.

To analyze the association between the Riskscore and immune microenvironment of HNSCC, we calculated immune cell infiltration using different methods. We first used ESTIMATE algorithm to assess immune cell infiltration and found that the HNSCC high-risk group had lower immune infiltration ( Figure 5A ). Based on MCPcounter, the immune cell score was calculated for the TCGA dataset. The results showed that the myeloid dendritic cells, immune cell scores of T cells, B lineage, CD8 T cells, cytotoxic lymphocytes were all lower in the high-risk group ( Figure 5B ). Using the ssGSEA function of the R package “GSVA”, the scores of 28 types of immune cells in HNSCC were analyzed, and their correlations with Riskscore and signature genes were computed. The data showed that Riskscore was closely negatively linked to the scores of the majority of the immune cells in HNSCC, including MDSC, immature B cells, activated CD8 T cells, mast cells, activated CD4 T cells, activated B cells. On the other hand, the signature genes were positively linked to the scores of the majority of the immune cells in HNSCC, while Riskscore was closely negatively linked to the scores of most of the cells ( Figure 5C ).

To investigate the differences in biological pathways in different risk groups, we performed a KEGG pathway enrichment analysis of DEG between high- and low-risk groups of HNSCC patients. Based on the enrichment results, it was found that the high-risk group was significantly enriched in the pathways of the ribosome, proteasome, and spliceosome ( Figure 6A ). Low-risk group was notably enriched in metabolism-related pathways, including linoleic acid metabolism and arachidonic acid metabolism ( Figure 6B ).

The mRNA expressions of 6 key genes in NOK and FaDu cells using qPCR. The corresponding results demonstrated that in comparison to NOK cells, the levels of ICOS, PSME1, SERPINA1, and SH3KBP1 genes were notably upregulated in FaDu cells, while the expressions of PSME1 and SERPINA1 genes were higher ( Figure 7A ). It has been reported that PSME1 has been shown to serve as a therapeutic target in a variety of tumors (39, 40). Therefore, we explored the role of PSME1 in HNSCC progression. Here, this study analyzed the effect of the PSME1 gene on FaDu cell invasion and metastasis by wound healing assay and transwell assay. The results showed that the reduction of PSME1 significantly inhibited the migration and metastasis of FaDu cells ( Figures 7B, C ).