The expression matrix was downloaded from Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/). To improve the accuracy of the study, we screened three eligible data sets. These data sets included tongue cancer and normal tissues. GSE13601, GSE34105, and GSE34106 contain 57, 78, and 43 samples, respectively. The SVA software package in R language was used to combine the samples of these three data sets to form a combined data set with 178 samples. Then, we downloaded the expression data of head and neck tumors from the TCGA database and screened the tongue squamous cell carcinoma samples, a total of 147 samples, including 140 tumor samples and 7 normal samples.

The limma package and DESeq package in R language were used for differential expression analysis, and the standard cutoff was set to false discovery rate (FDR) < 0.05 and |log2 (FC)| ≥ 1.5 for both GEO and TCGA, respectively (8, 9). The volcano package and heatmap package in R language were used to visualize DEGs. Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis, Gene Ontology (GO) analysis, and Gene Set Enrichment Analysis (GSEA) were performed for DEGs. Metascape website (https://metascape.org/) was used for pathway enrichment analysis of DEGs in TCGA. The t-test was used to compare the expression in the two groups, and ANOVA was used to compare that in multiple groups. We further verified the expression of upregulated DEGs using the GEPIA (http://gepia.cancer-pku.cn/).

The interaction between proteins corresponding to each gene was analyzed using String website (https://string-db.org/), and then the results were visualized using CytoScope. Genemania website (http://genemania.org/) was used for functional enrichment analysis of DEGs.

The DEGs obtained from combined dataset were crossed with those obtained from TCGA, including 5 up-regulated genes and 15 down-regulated genes.

Kaplan–Meier Plotter website (https://kmplot.com/) (citation) was used for survival analysis of the screened HUB genes.

We downloaded three data sets GSE13601, GSE34105, and GSE34106 containing tongue squamous cell carcinoma (TSCC) in the GEO database. Each data set contains tumor samples and normal tissue samples. We used the SVA package to merge these data sets and obtain a combined data set. To identify differences between TSCC samples and normal tongue samples, we performed principal component analysis (PCA) and found that the two groups could be separated in the combined data set (Figure 1A). We used the limma package [|log2 (FC)| ≥ 1.5, p < 0.05] in R to analyze the 178 samples and screened for DEGs. We identified 170 DEGs, including 104 upregulated genes and 66 downregulated genes, which were shown in the volcano diagram, where the red on the right represents the upregulated differentially expressed genes and the blue one represents the downregulated differentially expressed genes (Figure 1B). In the heatmap, groups were made according to data sets and sample characteristics to explore the expression of different genes. The top annotation represents the different data sets. DEGs are marked on the right side of the figure (Figure 1C). To further explore the functions of the screened DEGs, we used R language for further functional enrichment analysis of these genes in the KEGG database and enriched some up-regulated pathways such as cytokine–cytokine receptor interaction, viral protein interaction with cytokine and cytokine receptor, and primary immunodeficiency (Figure 1D). Then, GO analysis was performed on these genes, and the pathways enriched in upregulated genes were mainly cartilage development, connective tissue development, skeletal system development, and chondrocyte differentiation in the biological process (BP) module (Figure 1E). In the cellular component (CC) module, the enrichment pathways were mainly the collagen-containing extracellular matrix, basement membrane, apical part of the cell, and laminin complex. In the molecular function (MF) module, the pathways were mainly enriched in receptor ligand activity, signaling receptor activator activity, and metalloendopeptidase activity. GSEA analysis was also performed on the first six upregulated pathways, and the results were consistent with the previous analysis (Figure 2). The KEGG correspondence is in Supplementary Table S2.

DEGs of TSCC vs. normal samples in the TCGA database was performed using the DESeq package in R language [|log2 (FC)| ≥ 1.5, p < 0.05]. A total of 644 DEGs upregulated and 945 downregulated DEGs were found, as shown in the volcano diagram (Figure 3A). We used Metascape website to analyze these upregulated DEGs and found that the mainly enriched GO pathways were the aminoglycoside antibiotic metabolic process, proximal/distal pattern formation, and negative regulation of blood coagulation (Figure 3B). KEGG enrichment analysis of upregulated DEGs showed that the upregulated pathways were mainly folate biosynthesis, steroid hormone biosynthesis, complement and coagulation cascades, Staphylococcus aureus infection, and cytokine–cytokine receptor interaction. The main downregulated pathways were salivary secretion, calcium signaling pathway, and pancreatic secretion (Figure 3C).

Differentially expressed genes (DEGs) acquired in GEO were crossed with those in TCGA, and 5 common upregulated genes (CCL20, SCG5, SPP1, FOLR3, and KRT75) and 15 common downregulated genes (MAOB, HLFSELENBP1, CRISP3, TMPRSS2, TLX1, TGFBR3, APOD, PPP1R3C, SLITRK5, MFAP4, STATH, GSTM5, RNASE4, COX7A1, FHL1, SLC25A4, EFHD1, MYH11, MGP, SCGB2A1, LPIN1, ELF5, ADH1B, CKMT2) were screened (Figure 4A). The protein–protein interaction (PPI) network analysis of these five upregulated genes indicated that proteins expressed by other genes were all functionally related, except FOLR3 (Figure 4B). Then, we conducted interaction network analysis of these genes individually (Figure 4C).

We further verified the expression of upregulated DEGs using the GSE13601 data set. This data set includes tumor samples and their paired normal samples. After removing the individual samples from the database, we obtained 20 pairs of paired samples and then compared the expression levels of CCL20, SCG5, SPP1, FOLR3, and KRT75 genes in tumor tissues vs. normal tissues. It was shown that the expression of these five genes was higher in TSCC than that in normal tissue both in the whole sample and the most paired samples (Figure 5A). The same results can be obtained using the TCGA data set (Figure 5B).

To explore the relationship between the expression of these Hub genes and clinical stage, we used the TCGA database to explore and found that the expression level of CCL20 was significantly increased in patients with stages III and IV compared with stage I (p < 0.05) and the expression of SPP1 and SCG5 in stage IV was significantly higher than that in stage I patients (p < 0.05) (Figure 5C).

To further explore the clinical data of TCGA, multivariate Cox regression was performed using clinical data from the TCGA database (Figure 6). We find that risk increased when T stage > 2. Similarly, years of smoking also has an impact on the survival of patients; when the smoking time is more than 29 years, the risk is significantly increased. In addition, we found that when the TSCC stage is higher than grade III, age > 59, N stage > 1, and the number of cigarettes smoked per day > 2.5, the risk is increased; however, the P value of multivariate Cox analysis did not meet the statistical standards. Kaplan–Meier analysis was performed on clinical data using SPSS software, and we found that patients with a stage higher than III had a significantly worse prognosis than patients with a stage lower than III. AJCC_pathologic_T had a poor prognosis when it was higher than T2. Similarly, the AJCC_pathologic_N stage higher than N1 and smoking time more than 29 years remained risk factors (Figure 7).

These upregulated hub genes were further verified by GEPIA, and the results showed that the expression of these five upregulated genes in head and neck squamous cell carcinoma (HNSCC) was also higher than that in normal tissues (Figure 8). Kaplan–Meier Plotter website was used to analyze the relationship between the expression of these genes and prognosis, and the results showed that high expression of these genes was strongly associated with poor prognosis in HNSCC (p < 0.05) (Figure 9).