A total of 24 VMRGs were obtained from earlier reviews and are provided in Supplementary Table 1 (23), along with RNA-Seq data (HTSeq-FPKM), clinical information, and lifespan data for patients with HNSCC downloaded from UCSC Xena (http://xena.ucsc.edu/). In total, 494 HNSCC and 44 normal tissue samples were obtained. The GSE65858 (n = 270 samples) dataset was obtained from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/) for gene regulation profiling and longevity analyses of the corresponding patients for validation using outside sources. The ENSEMBL gene IDs were converted to gene symbol IDs to maintain consistency. Genes with expression levels < 50% in the samples were omitted from the analysis. The human HNSCC scRNA sequencing dataset GSE139324 was sourced from the GEO database. To reduce potential batch effects arising from integrating data derived from different public databases, we applied batch correction to the expression matrix. All raw expression values were first log2-transformed and normalized (log2(x+1)). Subsequently, the ComBat function in the “sva” R package was used to model batch variables and adjust for systematic biases related to batch and platform differences. Then, the batch-corrected data were used for downstream analyses, including differential expression analysis, prognostic model construction, and immune infiltration assessment. To evaluate the effectiveness of the correction, we performed principal component analysis (PCA), which demonstrated that the samples were no longer clustered by batch after adjustment, indicating that batch-related variability was substantially mitigated.

We performed an unbiased cluster analysis of gene expression utilizing the “ConsensusClusterPlus” program to find various VM-related clusters. The “K-Means” algorithm was employed, utilizing “Euclidean” distance as the metric, along with resampling of 80% of the items and conducting 1000 replications. The ideal k value was determined based on the ratio of unclear clusters. To investigate differences in clinical features between clusters, we generated heatmaps using the “pheatmap” R package. The PCA algorithm, along with the “ggplot2” R package, was employed to assess and visualize the sample distribution across different clusters.

The Wilcoxon rank-sum test was used to assess differences in the expression status of VMRGs between healthy and HNSCC cells. Subsequently, we plotted the copy number variation (CNV) of different VMRGs using “ggplot” and analyzed patients with different subtypes using the “gistic2” module of GenePattern. Waterfall plots were generated using “maftools” to detect somatic mutations of VMRGs and different subtypes in HNSCC. Additionally, the tumor mutational burden of each sample was calculated to determine the correlation between subtypes and tumor mutational burden.

For the VMRGs, the STRING database was used to construct a protein–protein interaction network. Meanwhile, the “limma” R package was used to examine differentially expressed genes (DEGs) between clustered subtypes, with |log2 FC| > 1 and FDR < 0.05 as cut-offs. Based on VM-related genes and DEGs between clustered subtypes, bioconcentration entries were obtained using Gene Ontology (GO) functional and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses (p < 0.05). Gene set variation analysis (GSVA) of the clustered subtypes was conducted using the “GSVA” R program. The threshold for significantly enriched functions was set at p < 0.05. GSEA was performed using the clusterprofile package (p < 0.05).

Using the ESTIMATE approach and the CIBERSORT algorithm, we determined that different molecular clusters connected to the VM had significantly different immune microenvironments and immune infiltrating cells. Furthermore, we compared the differences in immune checkpoints and HLA-related genes to evaluate the potential therapeutic efficacy of immunotherapy in different molecular clusters. Hypoxia scores were obtained from the cBioportal.

VM-related prognostic genes were identified from the clustered subtype-related DEGs using univariate Cox proportional hazards analysis with a cutoff of p < 0.05. To identify optimal prognostic genes, we conducted LASSO using the glmnet package and utilized stepwise selection based on the Akaike Information Criterion (stepAIC) in conjunction with the MASS program to identify the most significant subtype-related genes contributing to colon cancer prognosis. A multivariate Cox proportional hazards regression model was used to incorporate the prognostic genes with the highest predictive value. To create the VM score system, the coefficients of the multivariate Cox proportional hazards regression model were multiplied by the normalized prognostic gene expression levels:

Two categories of patients with HNSCC were formed using the median score from the VM scoring method: those who scored high on the risk scale and those who scored low. Using Kaplan–Meier survival analysis and log-rank statistical techniques, we identified the disparate overall survival (OS) rates between the two risk groups. Moreover, we conducted external validation of the cohort.

Nomograms are widely employed in clinical practice to visualize prognostic models. In our study, a significance criterion of p < 0.05 was used in the univariate and multivariate Cox regression analyses to identify important risk factors. Then, these identified risk factors were utilized as inputs for constructing a nomogram using the “rms” package. This nomogram provides a graphical representation in which cumulative points derived from the input variables enable the estimation of 1-year, 3-year, and 5-year mortality rates. Using the “ggplot2” software, the clinical usefulness of the nomogram relative to other markers was evaluated via decision curve analysis (DCA). This approach allows for a thorough assessment of the nomogram’s performance by considering the trade-offs between the potential benefits and drawbacks across different decision thresholds.

The “Seurat” R software was used to convert the 10× scRNA-seq data to a Seurat object. Groups of cells containing fewer than three mitochondrial genes were expressed at a rate greater than 10%, and cells with fewer than 50 genes were eliminated. The top 1500 most variable genes were used for PCA. Based on the top 15 main components, a cell clustering analysis was conducted using the “FindNeighbors” and “FindClusters” functions. The marker genes of different cell clusters were located using the “FindAllMarkers” tool, with the criteria being a |log2FC| > 1 and an FDR< 0.01 threshold. Additionally, cluster annotation was performed using “CellMarker 2.0”, to identify various cell types. The activity of a particular gene set in each cell was measured using the “ssGSEA” tool included in the Seurat package.

Each patient in the Cancer Genome Atlas (TCGA)-HNSCC dataset had a TIDE score determined using the tumor Page 4immune system breakdown and exclusions (TIDE; http://tide.dfci.atherard.edu/) to assess the immunotherapy response between the high- and low-risk groups. Additionally, the correlation between each patient’s immune-infiltrating cell levels and the prognostic model was computed for the TCGA-HNSCC cohort. The “oncoPredict” R package was used to determine the half-maximal inhibitory concentration (IC50) of chemotherapeutic drugs with data retrieved from the Genomics of Drug Sensitivity in Cancer (https://www.cancerrxgene.org/) database. The Wilcoxon signed-rank test was used to investigate the difference in IC50 values between the high- and low-risk groups and analyze the correlation between the risk score and drug sensitivity.

Human HNSCC cells (FADU, HSC-3 and SCC-25) were grown in MEM and DMEM/F12 supplemented with 10% fetal bovine serum, penicillin, and streptomycin. HaCat cells were grown in DMEM as a control. All cells were cultivated in an incubator at 37°C with 5% CO₂. Total RNA was extracted from the cells using a TRIzol RNA extraction kit, and reverse transcription was performed to convert the extracted RNA into cDNA. RT-qPCR was used to measure the expression quantity of VM-related genes, and the 2^-ΔΔCt method was used to quantify the suppression effects. All samples were analyzed in duplicates.

All statistical analyses were performed using R software (v4.3.1). The Wilcoxon test was used for pairwise comparisons between two groups, whereas the Kruskal–Wallis test was used for multiple group comparisons (* p < 0.05, ** p < 0.01, *** p < 0.001, **** p < 0.0001). The Kaplan–Meier method and log-rank test were used to analyze survival. Statistical significance was set at p < 0.05.

Figure 1 presents a flow diagram of the study. Our analysis of TCGA data revealed that most VM genes, including LAMC2, LOXL2, MAPK1, MMP2, MMP9, and PIK3CA, showed differential expression between HNSCC and normal tissues when compared using a boxplot ( Figure 2A ). To clarify the complex significance of VM-associated proteins, we built a network of protein–protein interactions. We found that five putative hub genes, MMP2, TGFB1, KDR, MMP9, and SNAI1, may play crucial roles in the pathological process of HNSCC ( Figure 2B ). The transcriptome relationships were investigated, and close correlations were observed among VM-related genes ( Figure 2C ); SNAI2 expression showed a strong positive correlation with LAMC2 expression, while LOXL2 expression showed a strong positive correlation with MMP2 expression, indicating that they function together. Moreover, we investigated the molecular alteration landscape of VMRGs in HNSCC and found that the most frequent variants were nonsense mutations ( Figure 2D ). The top five mutated genes were PIKK3CA, NOTCH1, EPHA2, ROCK1, and KDR. Next, we investigated the prevalence of CNV mutations and found that genes linked to VM exhibited notable changes in their CNV ( Figure 2E ). Enrichment analysis revealed that VMRGs were associated with proteoglycans in cancer, epithelial cell migration, endothelial cell migration, tissue migration, cell-cell junction organization, and epithelial cell proliferation ( Figures 2F, G ).

Subsequently, by increasing the clustering matrix (k) from 2 to 10, we analyzed VMRG expression in 494 samples using an unsupervised clustering method to examine the HNSCC categorization. When k = 2, consensus clustering demonstrated optimal performance, and two molecular subtypes were found, designated Clusters A and B ( Figures 3A, B ). Cluster A contained 283 samples, whereas Cluster B contained 211. The relationship between these molecular subtypes and clinical variables, including stage, grade, age, and survival status, are shown in a heatmap ( Figures 3C, D ). We conducted GSEA and determined the most significantly enhanced signaling pathways to investigate possible variations in biological activities among VM-related molecular subtypes. Although Cluster A was richer in pathways related to cancer-related proteoglycans and actin cytoskeleton regulation, Cluster B was mostly linked to drug metabolism and DNA replication-related pathways ( Figures 3E, F ). Additionally, GSVA was used to generate a heatmap of differentially enriched biological processes. Compared to Cluster B, Cluster A exhibited higher enrichment of extracellular matrix and cancer-related proteoglycans ( Figure 3G ).

Additionally, we evaluated the variations in immune-infiltrating cells, stromal scores, and immunological scores between the two molecular subtypes. The immune infiltration analysis showed that the Cluster A group exhibited a significant increase in the levels of 28 different types of immune cells ( Figure 4A ). Specifically, the abundances of CD4 memory resting T cells, regulatory T cells, monocytes, M1 macrophages, dendritic cells, and neutrophils were remarkably higher in Cluster A, whereas those of B cell memory, plasma cells, CD8 T cells, CD4 memory activated T cells, resting NK cells, M2 macrophages, and dendritic cells were significantly higher in Cluster B ( Figure 4B ). Additionally, Cluster A had considerably higher immunological and stromal scores than Cluster B ( Figure 4C ). Previous studies have shown that increased immune checkpoint expression correlates with a better response to immune checkpoint inhibitor therapy. Therefore, we investigated the expression levels of immune checkpoints across different molecular subtypes. As shown in Figure 4D , major immune checkpoints, including CD40, CD276, CTLA-4, PDCD1LG2, and NRP1, were overexpressed in the Cluster A subtype. Immune evasion was positively correlated with the TIDE prediction score, indicating immunotherapy resistance. In the TCGA-HNSCC cohort, Cluster A had a significantly higher TIDE score than that of Cluster B ( Figure 4E ). Additionally, hypoxia-responsive gene expression analysis revealed that Cluster A had a higher hypoxia score ( Figure 4F ). These findings suggest that the immunological microenvironments of the two molecular subtypes are significantly different.

We further investigated the link between the VM molecular subtypes and genomic changes (including CNV alterations and mutations). A higher non-synonymous tumor mutational burden was found in the protein-coding regions of the genomes of Cluster B patients ( Figure 5A ). The 20 genes with the highest mutation frequency in both subtypes are presented in Figure 5B . Notably, the mutation frequency of TP53 was higher in Cluster A (75%) compared with that in Cluster B (58%) ( Figure 5C ). CNV analysis revealed a significant difference in CNV patterns between VM-related subtypes ( Figure 5D ).

Based on the clustered subtype-related DEGs, 39 VM-related prognostic genes were identified using univariate Cox regression analysis ( Supplementary Table 2 ). Nine potential genes were further filtered out using LASSO Cox regression analysis to reduce the risk of overfitting and the range of candidate genes ( Figure 6A ). Finally, stepAIC Cox regression analysis was used to select five genes (CALML5, FMOD, PLAU, DEFB1, and CKM) for the VM-related prognostic model ( Figure 6B ). The expression levels of these five genes are shown in Supplementary Figure 1 . Each patient’s risk score was computed using the formula below: The risk score was calculated by adding the expressions of CALML5*(-0.0740), FMOD*(-0.2117), PLAU*(0.1514), DEFB1*(-0.0583), and CKM*(0.0581) ( Figure 6C ). Subsequently, patients were divided into low- and high-risk groups based on the median risk score. In both the TCGA-HNSCC ( Figure 6D ; p < 0.0001) and GSE65858 cohorts, the OS rate was higher for patients in the low-risk category compared with that in the high-risk group ( Figure 6G ; p = 0.015). Time-sensitive ROC curves were generated to assess the predictive ability of the prognostic model. The 1-, 3-, and 5-year AUCs of the TCGA-HNSCC cohort were 0.675, 0.705, and 0.63, respectively ( Figure 6E ). Similarly, in the GSE65858 cohort, the 1-, 3-, and 5-year AUCs were 0.687, 0.671, and 0.599, respectively ( Figure 6H ). Figures 6F, I illustrate the danger score distribution and life expectancy for the TCGA-HNSCC and GSE65858 cohorts, respectively. These findings support the strong performance of the VM-related prognostic model in predicting the prognosis of patients with HNSCC across several datasets.

Compared with other standard clinical features, risk scores emerged as independent predictive indicators in patients with HNSCC, as demonstrated by univariate and multivariate Cox analyses ( Figures 7A, B ). Univariate and multivariate Cox analyses of the TCGA-HNSCC cohort demonstrated that age, radiation exposure, and risk score were independent predictors of disease outcome. To determine the clinical significance of VM risk models, we developed a nomogram incorporating radiation and age to predict OS in patients with HNSCC based on the TCGA-HNSCC dataset ( Figure 7C ). The nomogram model demonstrated superior prognostic performance relative to the gene signature model. The prognoses of the high- and low-risk groups significantly differed (p < 0.0001; Figure 7D ). The 1-, 3-, and 5-year survival rates had AUC values of 0.772, 0.726, and 0.682, respectively ( Figure 7E ). The calibration curves demonstrated the accuracy of the model in forecasting survival rates at 1, 3, and 5 years ( Figure 7F ). Based on the DCA results the nomogram model was the most effective predictor ( Figure 7G ). Our nomogram exhibited a robust predictive capability and clinical relevance in evaluating the prognosis of patients with HNSCC based on these significant clinical characteristics.

To assess the immune infiltration landscape in HNSCC, we used the CIBERSORT algorithm to quantify the number of immune cells that infiltrated each sample ( Figure 8A ). The abundances of monocytes, macrophages, resting NK cells, and resting mast cells were higher in the high-risk group, whereas those of plasma cells, regulatory T cells, gamma delta T cells, activated NK cells, and activated mast cells were higher in the low-risk group. Additionally, five genes in the prognostic model were strongly correlated with the abundances of immune cells that infiltrate tumors, specifically macrophages. M1 macrophage abundance was negatively associated with DEFB1 and CALML5 and positively correlated with PLAU and FMOD, whereas resting dendritic cell abundance exhibited a positive correlation with both CKM and FMOD expression ( Figure 8B ). Next, we examined the detailed distribution of model genes in patients with HNSCC using single-cell RNA transcriptome data (GSE139324). The single-cell annotation results are provided in Supplementary Figure 2 . The dot plot revealed that most model genes were significantly expressed in the mast cells ( Figures 8C–E ). Additionally, analysis of the TIDE scores of all patients with HNSCC showed that the high-risk group had substantially higher scores, as evidenced by a positive correlation ( Figure 8F ). These findings suggest that immunotherapy may not be advantageous for patients with high-risk HNSCC. Moreover, we determined whether there was a distinct correlation between risk values and immune-related genes. The bar plot suggests that low-risk patients may have more robust immune activities ( Figure 8G ).

By comparing the sensitivity of the two risk groups to those of commonly used chemotherapeutic agents, we aimed to identify effective pharmaceuticals for patients with HNSCC and use this information to guide precision treatment. The half-lives of paclitaxel, docetaxel, AZD, buparlisib, cisplatin, cyclophosphamide, and docetaxel were lower in high-risk patients ( Figures 9A–I ), suggesting enhanced sensitivity to these treatments.

We confirmed the expression of five genes involved in the development of the prognostic model using RT-qPCR in four different types of HNSCC cells and haplotypes. The results revealed that FMOD and PLAU levels were significantly increased in HNSCC cells, whereas CKM and DEFB1 levels were significantly decreased ( Figures 10A–E ).