The dataset consists of 361 samples (329 tumor tissues and 32 adjacent non-cancerous tissues from TCGA) (Blum et al., 2018). The RNA-seq data and corresponding clinical information were downloaded from TCGA, and RNA-seq data in TPM format standardized by the Toil process was downloaded from UCSC XENA (https://xenabrowser.net/datapages/) (Vivian et al., 2017). To further compare the mRNA expression level between tumor and non-tumor samples, the logarithm of 2-fold change (log2FC) was calculated. Clinical information of the OSCC patients included TMN stage, clinical stage, histologic grade, Overall Survival (OS) event, Disease Specific Survival (DSS) event, Progression Free Interval (PFI) event, etc. Samples with unclear or error information were excluded. For GGPS1 correlation analysis, the RNA-seq data from the TCGA (https://portal.gdc.cancer.gov/) in the HTSeq-FPKM format were converted to TPM format, and log2FC is converted. The significant correlation condition was |cor_spearman|≥0.4 and p _ spearman <0.05. Data without corresponding clinical information were excluded. Specifically, we included samples from specific oral cancer sites, including the alveolar ridge, base of tongue, buccal mucosa, floor of mouth, hard palate, oral cavity, and oral tongue, and excluded samples from non-oral cancer sites such as the hypopharynx, larynx, lip, oropharynx, and tonsil. STRING database (https://string-db.org/) aims to integrate all known and predicted associations between proteins and analyze their sources, thus, it was used to explore the protein-protein interaction (PPI) networks of GGPS1 (Szklarczyk et al., 2021). In addition, GeneMANIA (http://www.genemania.org), a web interface that identifies functionally similar genes using available genomics and proteomics data, was used for the analyses of the protein-protein interaction network of GGPS1 (Warde-Farley et al., 2010). For single gene difference analysis, the RNA-seq data from the TCGA (https://portal.gdc.cancer.gov/) in the HTSeq-Counts format were converted to log2FC (Love et al., 2014). The conditions of single gene difference analysis set as Log2FC > 2 and p-value <0.05. MSigDB (http://www.broadinstitute.org/msigdb) was used to perform gene set enrichment analysis(GSEA) (Liberzon et al., 2011; Liberzon et al., 2015).

The HPA database (https://www.proteinatlas.org) is an extensive project that aims to map the entire human proteome by utilizing a combination of antibody-based proteomics and other omics techniques. In this study, the protein expression of GGPS1 was explored using the HPA database, and immunohistochemical samples of oral mucosa from patient ID 3917 and head and neck squamous cell carcinoma (HNSC) from patient ID 3669 were extracted.

The TISDB (http://cis.hku.hk/TISIDB/) database is a web-based platform that integrates multiple resources related to tumor immunology. In our study, we utilized the TISDB database for analyzing drugs that target GGPS1. This database offers a user-friendly interface and provides access to diverse data resources related to tumor immunology.

Patients with OSCC were divided into high and low GGPS1 expression groups according to the median expression level. Kaplan-Meier (KM) survival analyses were used to evaluate the association between GGPS1 expression level and Overall Survival (OS), Disease Specific Survival (DSS), and Progression Free Survival (PFS) by the R package (survminer, version 0.4.9 and survival, version 3.2.10). Calculated hazard ratios (HR) were presented, and p < 0.05 was statistically significant.

The role of univariate and multivariate Cox regression was assessing whether gender, age, TNM stage, the expression of GGPS1and so on were independent prognostic factors for the survival of OSCC patients. This study calculated hazard ratios (HR) and 95% confidence intervals (CI), and the significance threshold was set at p < 0.05. R package (survival, version 3.2.10) was a procedure of data processing. HR > 1 indicates disadvantageous factors.

The nomogram was used to make a forecast for OSCC prognosis by predicting the 1-, 3- and 5-year survival probability. Calibration curves were drawn to show the prediction effect. The ROC curves with time dependence, by using the R package (pROC, version 1.17.0.1 and ggplot2, version 3.3.3), were obtained for assessing the diagnostic value of GGPS1 mRNA expression.

This study used the R package (DESeq2, version 1.26.0) to analyze Different Expressed Genes (Love et al., 2014). Based on the dataset of heatmaps, STRING database, and GeneMANIA, there were Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses to evaluate GGPS1-associated potential gene functions with the R package (org.Hs.e.g.,.db, version 3.10.0 and clusterProfiler, version 3.14.3) (Yu et al., 2012). Gene Set Enrichment Analysis (GSEA) for interpreting gene expression data assesses whether an a priori defined set of genes shows statistically significant, concordant differences between two biological states (Subramanian et al., 2005; Subramanian et al., 2007). Set the False Discovery Rate (FDR) < 0.25 and the p.adjust<0.05 condition to be significant enrichment. The computational method was carried out using the R package clusterProfiler (version 3.14.3) and the R package ggplot2 (version 3.3.3) for data visualization.

This study conducted statistical calculations and graphing by R software (version 3.6.3, R Core Team, Vienna, Austria) and Adobe Photoshop software (version 25.0.0.60, Adobe Inc., Mountain View, CA, United States). The correlation analysis between clinical information and gene expression was implemented by Cox regression. Same as in previous studies, p < 0.05 was set up as the cut-off criterion.

The clinical and gene expression data of 329 tumor tissues and 32 adjacent non-cancerous tissues were downloaded from TCGA. This part was grouped according to the clinical information of OSCC, such as TMN stage, clinical stage, gender, race, age, and so on. Samples with unclear or error information were excluded, p < 0.05 suggested that this clinical feature may be related to GGPS1 expression in terms of statistical significance (Table 1). The results showed that Primary therapy outcome, histologic grade, Overall Survival (OS) event, Disease Specific Survival (DSS) event, and Progression Free Interval (PFI) event may be related to GGPS1 expression.

To begin with, the useful data about expressions of GGPS1 were analyzed and the results showed that the expressions of GGPS1 in tumor tissues were significantly higher in normal tissues (Figure 1A). The data of tumor samples were divided into two parts with a line of demarcation that was derived from the median expression level of GGPS1. Next, KM survival analyses were performed to assess the connection between GGPS1 expression levels and Overall Survival (OS) events (Figure 1B), Disease Specific Survival (DSS) events (Figure 1C), and Progression Free Interval (PFI) events (Figure 1D). It is worth noting that the GGPS1 expressions were very closely linked to OS, DSS, and PFI of the patients of OSCC. In addition, we utilized the HPA database to investigate the expression of GGPS1 in normal oral mucosa and HNSC tissue. Our findings indicated that GGPS1 protein expression was not significantly detected in oral mucosal tissues, whereas it was detected in HNSC (Figures 1E, F). The above results taken together indicated that GGPS1 possibly impacts cancer and the high expression level not only shortened the Progression Free Interval but also presaged a worse prognosis.

Kruskal–Wallis Test concerning clinical characteristics of OSCC and expressions of GGPS1 were used to analyze the correlations between the two. First of all, data analyses of the expressions of GGPS1 in histologic grade showed that with the deterioration of histological stages (Figure 2A), the expression increased obviously, and it was statistically significant. Interestingly, through the statistics and analysis of the data, it was found that the expression levels of GGPS1 were irrelevant to age (Figure 2B), smoker (Figure 2C), and alcohol history (Figure 2D). Results showed that a higher expression level of GGPS1 predicted a worse prognosis.

The role of univariate and multivariate Cox regression was to assess whether the expression of GGPS1 and clinical information of OSCC such as gender, age, histologic grade, and lymphovascular invasion were independent prognostic factors for the survival of patients. There were results that N stage, lymphovascular invasion, and most importantly, the expression of GGPS1 were risk factors. According to the multivariate analysis, these three factors could be regarded as independent prognostic factors as well (Figure 3). The nomogram was used to make a forecast for OSCC prognosis by predicting the 1-, 3- and 5-year survival probability (Figure 4A). The calibration plots of the above nomogram-predicated survival probability were very close to the ideal reference line (Figure 4B). Furthermore, the ROC curves with time dependence, drafted to assess the diagnostic value of GGPS1 mRNA expression, showed that the above predictions are reliable (Figure 4C). The results revealed that GGPS1was an independent predictor of poor prognosis with excellent diagnostic value.

To further understand the gene functions of GGPS1, analyses of gene correlations were performed. According to the analysis results, there were 23,818 genes positively correlated with GGPS1, and in parallel, 1,121 genes negatively correlated (Supplementary Table 1). The top 20 genes positively or negatively correlated with GGPS1 were shown in heat maps (Figures 5A, B).

After obtaining representative genes correlated with GGPS1, the protein-protein interaction (PPI) network analysis was mapped. It was found that there were 10 genes associated with GGPS1 using the STRING database (Figure 6A), and 20 genes related to GGPS1 using the GeneMANIA database (Figure 6B). Results indicated that some genes such as IDI2, IDI1, FDPS, and FDFT1 related to cholesterol synthesis were most closely related to GGPS1 (Nakamura et al., 2015; Dong et al., 2020; Ershov et al., 2021). In addition, to explore the potential biological functions related to GGPS1 expression, GO/KEGG analysis was executed and outcomes showed that GGPS1 expressions are associated with a series of functions such as Terpenoid backbone biosynthesis, regulation of small molecule metabolic process, regulation of lipid metabolic process, cholesterol metabolic process, and so on (Supplementary Table 2, Figures 6C, D).

Single gene difference analysis of GGPS1 was performed and the corresponding volcano map was drawn (Supplementary Table 3, Figure 7A). The total number excluding the null value is 35,620. Selected threshold value is |log2(FC)|>2 & p.adj<0.05, and the number meeting this threshold value is 54, among which 27 are high expressions (log2FC is positive) and 27 are low expressions (log2FC is negative). Furthermore, Gene Set Enrichment Analysis (GSEA) analysis was carried out for GGPS1 gene function prediction (Supplementary Table 4). Outcomes of GSEA analyses showed that a high expression of GGPS1 was associated with the activation of ascorbate and aldarate metabolism (Figure 7B), pentose and glucuronate interconversions (Figure 7C), porphyrin and chlorophyll metabolism (Figure 7D), pathways in cancer (Figure 7E), chemokine signaling pathway (Figure 7F), oxidative phosphorylation (Figure 7G), and PPAR signaling pathway (Figure 7H). The above-mentioned results comprehen show that GGPS1 did play a role in the development of OSCC.

The development of drugs targeting GGPS1 is crucial for the effective treatment of cancer patients. In this study, we utilized the TISDB database to analyze the available drugs that target GGPS1, and identified a total of 14 small molecule drugs. Notably, farnesyl diphosphate (DB07780) and geranylgeranyl diphosphate (DB07841) exhibited the highest number of targets, with a total of 7 targets (Figure 8). Additional information can be found in Table 2.