Transcriptome data, genomic mutation data, and corresponding clinical information from HNSCC patients, which contained 44 normal and 495 primary tumor tissues, were downloaded from The Cancer Genome Atlas (TCGA). GSE65858, containing 270 tumor samples with survival data, was obtained from the Gene Expression Omnibus (GEO) database to verify the prognostic signature. GSE30784 and GSE37991 were used to differentiate the expression levels of the genes in the signature between normal and HNSCC tissues. Moreover, we downloaded the expression data of one normal cell line (HaCaT) and six HNSCC cell lines (93VU147T, SCC61, SCC047, SCC090, SCC25, and SQ20B) from GSE62027. A total of 169 glycosyltransferase-related genes (GTRGs) were derived from a previous study (Mohamed Abd-El-Halim et al., 2021).

A volcano plot and heatmap were used to visualize the DE-GTRGs, which were defined as those with p < 0.05 and |log2 (foldchange)| > 1. Then, a protein–protein interaction (PPI) network of DE-GTRGs was constructed by the GeneMANIA database to identify the interactions of glycosyltransferases. Moreover, the correlation between the expression of various DE-GTRGs was also analyzed.

Univariate and multivariate Cox regression analyses were used to screen GTRGs to construct the prognostic signature. The signature reflects both the expression levels of the selected genes and their relative regression coefficient weights calculated from the multivariate Cox analysis. Patients in the training set were classified into high-risk and low-risk groups based on the median risk score. The Kaplan–Meier (KM) method was used to describe the differences in overall survival (OS) between the two groups. Receiver operating characteristic (ROC) curve analysis was further used to demonstrate the specificity and sensitivity of the signature. Principal component analysis (PCA) and T-distributed stochastic neighbor embedding (tSNE) methods were used to conduct dimension reduction analysis for all patients to evaluate the effect of the signature on the ability to distinguish between them. Finally, we incorporated the risk score and clinical parameters into the univariate and multivariate Cox regression analyses to prove that the risk score was an independent risk factor for prognosis.

Clinical parameters including age, sex, tumor grade, tumor stage, T stage, and lymph node metastasis were used to perform a stratified analysis of OS. In addition, we analyzed differences in the risk score between subgroups based on the aforementioned clinical parameters.

A nomogram was established based on clinical parameters and the risk score. The calibration curves were drawn to assess the consistency between the actual and predicted probabilities of 1-year, 2-year, and 3-year survival. Furthermore, nomograms and calibration curves were also constructed according to the GSE65858 dataset.

GSEA was performed to identify the pathways enriched in the high-risk group or in the low-risk group according to the reference gene set Hallmark and KEGG analysis. An |NSE| > 1, a p value <0.05 and a false discovery rate (FDR) < 0.25 were set as the cut-off values.

The R package “ESTIMATE” was used to analyze the immune, stromal, and estimate scores and tumor purity between the high- and low-risk groups. TIMER, CIBERSORT, QUANTISEQ, MCPCOUNTER, and EPIC were used to compare the differences in immune cells between the two groups using the Wilcoxon test. Moreover, we analyzed the differences in immune cells and immune function between the two groups by using the ssGSEA method. Immune checkpoints, obtained from the TISIDB, were used to evaluate the differences between the two groups. Finally, to assess the value of the signature in patients receiving immunotherapy, multiple datasets were obtained to evaluate whether the signature could be used in immunotherapy patient cohorts. The cohorts included Van Allen’s cohort (CTLA-4 blockade in 40 metastatic melanoma patients) (Van Allen et al., 2015), Braun DA’s cohort (PD1 blockade in 311 advanced clear cell renal cell carcinoma patients) (Braun et al., 2020), Riaz N’s cohort (nivolumab in 50 advanced melanoma patients) (Riaz et al., 2017) and David Liu’s cohort (PD1 blockade in 121 metastatic melanoma patients) (Liu et al., 2019a).

TMB was compared between the high- and low-risk groups and was used to analyze its prognostic value in HNSCC patients through KM analysis. In addition, survival analysis was performed according to TMB and the risk score. Gene mutation frequency was also analyzed between the two groups based on the mutation data from the TCGA, and the top 30 mutated genes are shown on a waterfall plot. We also obtained genes with significant mutation differences between the two groups.

Immunohistochemistry images of PYGL, ALG3, FUT2, and KDELC1 were obtained from the HPA database. UALCAN was utilized to compare the protein levels of PYGL, EXT2, FUT2, and KDELC1 between normal and tumor tissues. Finally, a KM plotter was used to analyze the prognostic value of the five genes in HNSCC patients.

A flow chart of the current study is shown in Figure 1A. We first obtained the expression values of 169 glycosyltransferase genes in normal and tumor tissues from the TCGA, and differential analysis was performed according to the abovementioned criteria. Twenty-nine DE-GTRGs were identified, among which 10 were downregulated and 19 were upregulated (Figures 1B,C). Based on the DE-GTRGs, a PPI network was constructed through the GeneMANIA database, and significant interactions were observed among these glycosyltransferases (Figure 1D). In addition, a generally positive or negative correlation was found among these DE-GTRGs (Figure 1E).

Univariate and multivariate Cox regression analyses were used to screen genes associated with patient prognoses and construct a prognostic signature. The results of the univariate analysis showed that five DE-GTRGs were associated with the prognosis of HNSCC patients, of which PYGL, ALG3, KDELC1, and EXT2 were risk factors and FUT2 was a protective factor (Figure 2A). Then, a prognostic signature based on five GTRGs was constructed utilizing multivariate Cox analysis, and the coefficients of the five genes are shown in Figure 2B. Risk score = (0.043622*KDELC1) + (0.009262*ALG3) + (0.00617*PYGL) + (0.003267*EXT2) - (0.00738*FUT2). The median risk score was applied to divide the HNSCC patients into a high-risk group and a low-risk group. Patients in the low-risk group had a better prognosis than those in the high-risk group (Figure 2C). As the risk score increased, there were more deaths and shorter survival times (Figure 2D). The areas under the curve (AUCs) of the signature at 1, 2, and 3 years were 0.619, 0.656, and 0.675, respectively (Figure 2E). The AUC of the risk score was 0.630, which was higher than that of the other six clinical parameters (Figure 2F). PCA and tSNE analyses were used to reduce dimensionality in all patients, and we found that patients with different risk scores could be distinguished significantly (Figures 2G,H). To analyze the prognostic value of the signature, univariate and multivariate analyses were used, and the results showed that the risk score was an independent risk factor for prognosis (Figures 2I,J).

To validate the signature constructed through the TCGA, we obtained the expression profile and clinical parameters of the GSE65858 dataset. Similar to the training model, KM analysis showed poor prognosis in the high-risk group (Figure 3A). An increasing number of deaths were observed as the risk score increased (Figure 3B). ROC curve analysis revealed that the signature was of great value for patient prognosis. The AUCs of the signature at 1, 2, and 3 years were high, and the AUC of the risk score was 0.625, similar to the training model and higher than the other clinical parameters (Figures 3C,D). The results of PCA and tSNE analyses suggested that the prognostic signature can effectively distinguish high- and low-risk patients (Figures 3E,F). Again, the signature was an independent risk factor for patients in the external dataset (Figures 3G,H).

To determine the value of the signature in different clinical subgroups, we performed a stratified analysis. We divided the patients into age ≤ 65 years, age >65 years, female, male, grades 1–2, grades 3–4, stages I–II, stages III–IV, T1–2, T3–4, N0, and N1–3 groups based on various clinical parameters (age, sex, tumor grade, tumor stage, T stage, and lymph node metastasis). The signature showed good efficacy in differentiating patient outcomes across all subgroups (Figure 4), indicating that the signature can be applied to all patients regardless of there being multiple clinical variables.

Due to the important value of the signature in different clinical subgroups, we then analyzed its correlation with clinical parameters. We found that the signature correlated significantly with tumor stage and tumor size; that is, in advanced tumors (stages III–IV) or when the tumor size was large (T3-4), the risk score was higher (Figures 5A,B). We also integrated the signature and clinical features to construct a nomogram. The nomogram was constructed to predict the 1-year, 2-year, and 3-year survival probabilities in the TCGA cohort (Figure 5C) and in the GSE65858 cohort (Figure 5D). The calibration curve showed good agreement between the predicted and actual 1-year, 2-year, and 3-year survival rates in both the training and validation cohorts (Figures 5E,F). GSEA showed that the high-risk group was associated with a variety of pathways related to tumor development, such as angiogenesis, epithelial–mesenchymal transition (EMT), glycolysis, hypoxia, and the TGF beta signaling pathway (Figure 5G), while the low-risk group was mainly enriched in IL2-STAT5 signaling, KRAS signaling, metabolism-related pathways, and others (Figure 5H).

We first evaluated the immune score, stromal score, ESTIMATE score, and tumor purity. The immune and ESTIMATE scores were significantly lower in the high-risk group, while tumor purity was higher in the high-risk group (Figure 6A). We also compared the differences in immune cell infiltration between the two groups according to various algorithms. Based on the TIMER, CIBERSORT, QUANTISEQ, MCPCOUNTER, and EPIC, we found that B cells, CD8+T cells, dendritic cells, and Tregs were lower in the high-risk group and that cancer-associated fibroblasts (CAFs) were higher in the high-risk group (Figures 6B,C). In addition, we estimated the difference in immune cells and immune function between the two groups by using the ssGSEA algorithm. The results showed that immune cell infiltration was similar to that of the previous algorithm (Figure 7A), while multiple immune function scores were lower in the high-risk group, including checkpoints (Figure 7B). We, therefore, analyzed the differences between immune checkpoints in the two groups. Most checkpoints were significantly different between the two groups, and the expression levels of immune checkpoints were lower in the high-risk group, a pattern similar to that using the ssGSEA algorithm (Figure 7C).

Immunotherapy has been considered to play an important role in a variety of malignant tumors and can markedly improve the prognosis. Therefore, it is necessary to identify people who may benefit from immunotherapy. We evaluated the value of the prognostic signature in immunotherapy for malignant tumors using four cohorts. Patients in the high-risk group had a shorter survival time (OS and progression-free survival) than those in the low-risk group in the four cohorts (Figures 7D–I).

TMB was calculated based on the somatic mutation obtained from the TCGA and was compared between the two groups. TMB in the high-risk group was higher than that in the low-risk group (Figure 8A), and the higher TMB group was associated with shorter survival time (Figure 8B). In addition, a stratified analysis based on TMB and the risk score showed significantly worse outcomes in the high-TMB + high-risk group than in the low-TMB + low-risk group (Figure 8C). Gene mutations visualized as a waterfall plot showed that the five most frequent somatic mutations in the high-risk group were those in TP53, TTN, FAT1, CDKN2A, and MUC16 (Figure 8D), while the five most frequent somatic mutations in the low-risk group were those in TP53, TTN, FAT1, NOTCH1, and PIK3CA (Figure 8E). Differential analysis between the two groups showed that the mutation frequencies of TP53, CSMD1, NPAP1, AJUBA, CDKN2A, MUC17, and PRDM9 were higher in the high-risk group (Figure 8F).

The mRNA expression levels of the five genes in GSE30784 revealed that PYGL, ALG3, EXT2, and KDELC1 were upregulated while FUT2 was downregulated in HNSCC (Figure 9A). The differential expression analysis of the genes in the 40 paired normal and tumor tissues in GSE37991 showed similar results (Figure 9B). Furthermore, the expression level detected in one normal and six tumor cell lines was essentially consistent with that in tissues (Figure 9C). The immunohistochemistry images of PYGL, ALG3, FUT2, and KDELC1 in normal oral mucosa and HNSCC tissues are shown in Figure 9D. Meanwhile, the protein levels in the UALCAN database suggested that PYGL, EXT2, and KDELC1 were higher while FUT2 was lower in HNSCC tissues (Figure 9E). Finally, we performed survival analysis using KM curves in the KM plotter database. The results showed that the prognosis of patients with high expression of PYGL, ALG3, EXT2, and KDELC1 was worse, while high expression of FUT2 indicated longer survival time (Figure 9F).