3.1 Identify DEGs
The samples in GSE52088 were screened, DOK cells were selected as the oral precancerous lesion group. SCC-15, SCC-25 and CAL-27 cells were employed as the OSCC group. Principal component analysis (PCA) was performed on four samples from the two groups. As shown in Figure 1A, there was a significant difference in gene expression between the two groups. Then the expression matrix of two groups of samples was obtained. A total of 439 DEGs were obtained after screening the missing and blank values and filtrating according to the P < 0.001 and | log2 (Fold Change) | > 2 standard. Compared with the oral precancerous lesion group, 215 down-regulated genes and 224 up-regulated genes were found in OSCC group (Figure 1B, 1C).
3.2 Identification of potential key genes associated with both OLK and HNSCC
4584 OLK related genes and 15759 HNSCC related genes were downloaded from CTD, which were used to explore the potential key genes of OLK and HNSCC. A total of 123 genes were obtained by taking the intersection of the OLK related genes, HNSCC related genes and 439 DEGs, then a Venn diagram was drawn (Figure 2A). After that the association between diseases and genes was scored by CTD, the top 50 genes in OLK and HNSCC were obtained, respectively (Figure 2B, 2C). So the intersection of two groups were taken finally. The results showed that there were 12 DEGs (RDH10, TCEA2, SWAP70, BCAP31, LAPTM4B, VTI1A, NR3C1, COX6A1, SOX2, HMOX1, TFAP2A and BIRC5) had high correlation scores with both OLK and HNSCC (Figure 2D).
3.3 PPI network construction and identification of 12 potential key gene networks
The PPI network of 439 DEGs was constructed with STRING database, including a total of 420 nodes and 1826 edges (Figure 3A). Subsequently, the PPI of 12 potential key genes was isolated from the PPI of 439 DEGs (Figure 3B). Except for RDH10, the majority of potential key genes had many interactions, in which BIRC5, SOX2, NR3C1 and COX6A1 had more interactions with each other, indicating that they were in a relatively core position in the network.
3.4 Functional enrichment analysis of DEGs
The DEGs contained in Figure 3B was conducted with enrichment analysis of GO and KEGG pathways. There were 104 GO analysis results, including 51 BP, 33 CC and 20 MF (P < 0.05). The top 20 results were selected according to the p-value and the circle graph was drawed (Figure 4A). BP analysis showed that DEGs mainly enriched in mitotic cell cycle, cell division and chromosome segregation. In CC, DEGs were mainly distributed in the necleoplasm, nucleus and cytosol. MF analysis showed that DEGs were significantly enriched in protein binding, microtubule binding, ubiquitin-cytochrome-C reductase activity and ATP-dependent microtubule motor activity. In addition, KEGG pathway analysis showed that DEGs mainly focused on oxidative phosphorylation, cell cycle, thermogenesis, chemical carcinogene-reactive oxygen species, prion diseases, phagosome, pathways of neurodegeneration-multiple diseases pathways, etc. (Figure 4B). Detailed results of the GO enrichment analysis were provided in Table 2.
3.5 Validation of potential key genes in the external dataset
The expression levels of the 12 potential key genes were extracted from each sample in the GSE26549 dataset, and the expression hot map was drawn (Figure 5A). According to the sample description, all patients were divided into two groups: OLK without oral cancer development group and OLK with oral cancer development group. Then the expression levels of 12 potential key genes were analyzed statistically. As a result, the expressions of LAPTM4B and COX6A1 in OLK with oral cancer development group were significantly higher than those in OLK without Oral cancer development group (P < 0.01), while expression of NR3C1 was contrary to the expressions of LAPTM4B and COX6A1 (P < 0.05, Figure 5B). The expression levels of other potential key genes were not significantly different between the two groups. ROC analysis showed that the area under the curve (AUC) value of LAPTM4B, NR3C1 and COX6A1 to distinguish OLK from OSCC was 0.679, 0.633 and 0.675, respectively, which were higher than other clinical parameters (Age: 0.636; Gender: 0.517; Drinking habit: 0.507; Smoking habit: 0.497; Treatment measures: 0.501; Pathological grade: 0.572) (Figure 5C). When the three indicators were combined to distinguish OLK from OSCC, the AUC value was 0.753, indicating that the combination of the three treatments had a high diagnostic accuracy for OLK cancer (Figure 5D). Then the differences of the three key genes in different pathological grades of OLK were compared. The expression of LAPTM4B was significantly different between hyperplasia and dysplasia groups (P < 0.001), while the expression of NR3C1 and COX6A1 was not different between the two groups (Figure 5E). The expression of LAPTM4B was significantly different between each of the hyperplasia, mild dysplasia and moderate dysplasia three groups (P < 0.05). The expression of NR3C1 was significantly different between mild dysplasia and moderate dysplasia groups (P < 0.05), and there was no significant difference in other genes between the three groups (Figure 5F). Unfortunately, there were only 2 cases of severe dysplasia in the dataset, which could not be further statistically analyzed.
3.6 Immune infiltration analysis of potential key genes
Immune infiltration analysis showed that the correlation between three potential key genes and immune cells was different. LAPTM4B was weakly correlated with all the immune cells (correlation coefficient | r | < 0.4), of which CD56dim cells have strongest correlation (r = 0.291) (Figure 6A-C). NR3C1 was strongly positively correlated with Tcm and T Helper cells with correlation coefficients of 0.418 and 0.422, respectively (Figure 6D-F). COX6A1 was strongly negatively correlated with Eosinophils and Tcm, with correlation coefficients of -0.436 and -0.422, respectively (Figure 6G-I).
3.7 Expression levels of key genes in cells
The expression levels of three target genes in four groups of cells were analyzed.
The mean relative expression levels of LAPTM4B in DOK, SCC15, SCC25 and CAL27 cells were 1.0005, 1.5395, 1.2281 and 1.6368, respectively. The mean relative expression levels of NR3C1 in the above cells were 1.0013, 0.5693, 0.4163 and 0.2141, respectively. The mean relative expression levels of COX6A1 in the above cells were 1.0013, 2.2940, 2.2959 and 2.0781, respectively. These results showed that the expression levels of LAPTM4B (Figure 7A) and COX6A1 (Figure 7B) were positively correlated with OLK canceration, while NR3C1 (Figure 7C) was negatively correlated with OLK canceration (P < 0.05). Real-time PCR results of the three potential key genes were consistent with the results of bioinformatics analysis, which further confirmed the possibility of the three genes could serve as key genes for predicting OLK canceration.