50 patients with OSCC were enrolled in this study after providing informed consent. This study was approved by the institutional review board of MacKay Memorial Hospital (approval numbers: 12MMHIS178 and 15MMHIS104). Tumor specimens were collected from patients during OSCC surgery. Laser capture microdissection was performed to isolate relatively pure tumor cells for DNA extraction according to previously established protocols (19). 10 mL of whole blood were collected in Vacutainer tubes containing ethylenediaminetetraacetic acid as anticoagulant (Becton Dickinson, Franklin Lakes, NJ). Genomic DNA was extracted from blood or tumor specimens using the QIAamp DNA Blood Mini Kit according to the manufacturer’s instructions (QIAGEN, Hilden, Germany).

Demographic data, including age, sex, clinical stage, perineural invasion, and lymphovascular invasion were retrospectively obtained from the patients’ medical records. Clinical staging was performed according to the American Joint Committee on Cancer (AJCC 7th edition) guidelines for tumor, node, and metastasis TNM classification (20). None of the patients which were enrolled in this study had received adjuvant chemotherapy or radiotherapy before surgery. WES data from the TCGA-HNSCC dataset was collected and downloaded from the Genomic Data Commons portal (https://portal.gdc.cancer.gov/). Using the TCGA-HNSCC dataset, all cases where the primary site was located at the tongue, lip, mouth floor, tonsil, gums, palate, or oropharynx, were included. This dataset was called “TCGA-OSCC” dataset. In total, 387 OSCC patients with somatic mutation data were analyzed herein. Only mutations in the coding regions and in splicing sites were retained, whereas mutations in the introns, intergenic regions, and in untranslated regions (UTR) were filtered out.

WES was performed with the SureSelect Human All Exon v6 + UTR Enrichment Kit (Agilent, Santa Clara, CA), followed by sequencing on a NextSeq500 DNA sequencer (Illumina, San Diego, CA). Downstream analysis was performed as previously described (21). The software used for WES analysis is listed in Table S1 . Somatic mutation was called using our pipelines were shown in Figure S1 . Sequencing data was aligned to the human genome (NCBI build GRCh38/UCSC hg38) using BWA-MEM. SAMtools was used for the file format conversion from SAM to BAM. Six different programs were used to call somatic mutations, including two callers (Muse and SomaticSniper) for single nucleotide variants (SNVs) (22, 23) and four callers (Mutect2, Strelka2, VarScan2, and VarDict) for both SNVs and short insertion and deletion variants (indels) (9, 23–26). These variant callers were run with default parameters and further filtering of the data was based on the following criteria: (1) Mutations called by Muse, SomaticSniper, and Mutect2 were labeled as “PASS” in the FILTER column, (2) Strelka2 algorithm with WES default parameters was used to identify somatic mutations (27). Only variants labeled “PASS” were considered high-quality variants (28). (3) Mutations called by VarScan2 were labeled as “Somatic” in the somatic status column, and (4) mutations called by VarDict were labeled as “StrongSomatic” in the INFO column. Thereafter, mutations outside the targeted region were removed. The filtered mutations were considered somatic mutations. The command line calls somatic mutations as described in Table S2 . All WES data in this manuscript were submitted to Short Reads Archive under the BioProject accession PRJNA749133 and SRA Run Selector project (https://www.ncbi.nlm.nih.gov/Traces/study/?acc=PRJNA749133&o=acc_s%3Aa).

The somatic mutations were annotated using Ensembl Variant Effect Predictor (version 102, https://asia.ensembl.org/Homo_sapiens/Tools/VEP) (21). Potential mutational driver genes in OSCC were identified and annotated using the InToGene platform (https://www.intogen.org/search) and Bailey et al. datasets (29, 30).

To reduce false-positive somatic mutations which might originate from germline mutations or might have been accidentally generated during sample preparation, DNA amplification, sequencing, and ambiguous mapping (31), the following filter flags were used to annotate and filter somatic mutations: (1) removing common polymorphisms (SNPs): minor allele frequency in the 1000 Genomes Project or The Genome Aggregation Database (gnomAD) > 1%; (2) removing Panel-of-Normal (PoN): A normal panel was created from 50 normal samples using GATK (Genome Analysis Toolkit) tool CreateSomaticPanelOfNormals; (3) removing oxodG artifacts: oxidation-damaged base 8-oxodG (8-oxoguanine) was identified and filtered using GATK tools CollectSequencingArtifactMetrics and FilterByOrientationBias; (4) removing strand bias, multiallelic site, and clustered events: strand bias is a type of sequencing bias wherein one DNA strand is favored over another by the variant. A multiallelic site is a genomic locus that contains two or more alternative (Alt) alleles. Clustered events are several variants that are clustered in a region of the genome. These artifacts were estimated and filtered using the GATK tools CalculateContamination, GetPileupSummaries, and FilterMutectCalls; and (v) Alt allele count filter: we filtered out mutations with Alt alleles in tumor derived date (T_Alt) < 4 and mutations with Alt alleles in normal data (N_Alt) ≥ 4 ( Figure 2A ) (32).

Somatic mutations were called using multiple callers and stored as VCF (Variant Call Format) files. Following variant calling and filtering, VCF files from six tools were merged according to each sample ID and genomic position (e.g., ID-chr1-123). The number of variants hits detected for each mutation in the merge files was then counted. Thereafter, mutations that were not identified by two or more variant callers were removed. The tumor mutation burden (TMB) was calculated using the number of non-synonymous mutations per mega-base (Mb) in the target region of SureSelect Human All Exon v6 + UTR (91.08 Megabases). The target region BED file is available online at the SureDesign website (https://earray.chem.agilent.com/suredesign/).

Sanger sequencing and IGV was performed to validate somatic mutations (33, 34). For Sanger sequencing, individual primer sets designed by Primer3 (version 0.4.0) are listed in Table S3 . Polymerase chain reactions (PCRs) were performed using the KAPA LongRange HotStart PCR Kit (KAPA Biosystems, Wilmington, MA, USA). Amplicons were sequenced on an ABI 3730xl DNA Analyzer (Applied Biosystems, Foster City, CA, USA) with the BigDye Terminator Cycle Sequencing Kit (Applied Biosystems).

The top 20 most frequently mutated genes [mutated in at least 7/50 (14%) patients] were selected and examined for false-positive rates using IGV. Mutations were considered “true-positive” based on the following criteria: (1) number of Alt alleles < 3 in normal cells and ≥ 3 in tumors; (2) both forward and reverse strands have at least one mutant allele; (3) number of mismatches within a 40 bp window ≤3; and (4) allelic configurations of the mutation are multiallelic variants (21).

The 200 most frequently mutated genes [mutated in at least 4/50 (8%) patients] in our data and in the TCGA-OSCC dataset were selected for the creation of a circular plot using Circos-0.69 (http://circos.ca/software/). Track 1 (inner circle) visualizes the mutation frequency of the genes in the TCGA-OSCC dataset ( Figure 4 ), whereas Track 2 (outer cycle) illustrates the genomic profile from our WES data.

The 200 most frequently mutated genes in our study were imported into the analysis pipeline of the Ingenuity Pathway Analysis (IPA, QIAGEN, CA, USA; http://www.ingenuity.com/products/pathways_analysis.html). IPA was used to examine which canonical pathways were enriched within our candidate genes. In IPA, Fisher’s exact test was used to determine whether a canonical pathway was enriched within a data set. A −log [P value] > 1.3 (corresponding to a P value of <0.05) was set as the threshold for statistical significance. The United States Food and Drug Administration (FDA) Table of Pharmacogenomic Biomarkers in Drug Labels was used to match our candidate genes against FDA-approved drugs (https://www.fda.gov/).

Data are presented as mean ± standard error of mean. The Chi-squared test, Fisher’s exact test, and Mann–Whitney U test were used for statistical analysis. The receiver operating characteristic curve (ROC) analysis was used to identify the optimal thresholds of read-depth (DP) and minor allele frequency (MAF) in mutation calling. Overall survival (OS) was defined as the duration from the first date of diagnosis to death or last date of follow-up. Kaplan–Meier analysis was used to compare the OS between two groups. P values less than 0.05 were considered to be statistically significant.

In this study, we have collected and analyzed tumor specimens and matched blood samples from 50 patients with OSCC. 48 patients were male and 2 were female, with an average age of 59.6 years (range: 40–89 years). All patients were confirmed to have squamous cell carcinoma. The most common primary sites were the bucca (28%, 14/50) and the gingiva (24%, 12/50). The detailed clinical characteristics of our study subjects are described in Table 1 .

Six different variant calling tools were used to identify somatic mutations in the 50 paired tumor/normal samples: Muse, Mutect2, SomaticSniper, Sterlka2, VarScan2, and VarDict ( Figure 1 ). All six callers were used to identify somatic SNVs. Moreover, somatic indels were identified using Mutect2, Strelka2, VarScan2, and VarDict. A total of 163,069 somatic mutations were found by the six callers using the default parameters in the target region (coding and splicing region). Muse, Mutect2, SomaticSniper, Sterlka2, VarScan2, and VarDict detected 10,019, 79,933, 9037, 58,070, 3656, and 37,240 somatic mutations, respectively ( Figure 1 ). To reduce false-positive mutations from variant calling, filter flags were used to annotate and assess somatic mutations ( Figure 2A and Table S4 ). Somatic mutations marked with filter flags, including “common SNP”, “oxodG” oxidative damage, “StrandBias”, “Multiallelic site”, “Clustered events”, and “PoN” variants were then filtered out. Furthermore, somatic mutations that were lacking sufficient evidence to be called a somatic mutation, such as those with low Alt alleles in the tumor sample (T_Alt < 4 alleles) and high Alt alleles in normal samples (N_Alt ≥4 alleles), were filtered out. As result, the filtered (PASS only) data contained 8,730, 23,691, 3,574, 54,846, 2,683, and 23,545 somatic mutations detected by Muse, Mutect2, SomaticSniper, Sterlka2, VarScan2, and VarDict, respectively ( Figure 1 and Table S4 ). As depicted in Figure S2A , somatic mutations marked “PASS” filter flag were successfully verified by IGV and Sanger sequencing. However, the somatic mutations which have been labeled “N_Alt ≥4,” “T_Alt < 4,” “Multiallelic site,” and “Clustered event” were not confirmed by Sanger sequencing ( Figures S2B–E ). Furthermore, our results show that all variant callers except Mutect2 and VarDict mis-detected dinucleotide mutations as SNVs. In total, 53 dinucleotide mutations were identified in our study ( Table S5 ), all of which were confirmed by Sanger sequencing ( Figure S3 ). Finally, the VCF files from each caller were merged ( Figure 1 ). Somatic mutations that were only identified by a single caller were excluded. Overall, the ≥2 caller agreeing data set only recovered 8% of the variants in the unfiltered data set. Our study identified 13,730 somatic mutations in 7,729 unique genes, including 3,057 synonymous, 8,541 missense, 17 start loss, 854 stop gain, 9 stop loss, 3 stop retain, 969 splicing site, 53 dinucleotide, 183 frameshift, and 44 in frame mutations ( Figure 1 ).

Figure 2 and Table S6 summarize the process of filtering mutations in the six callers. These variant callers detected a different amount of SNVs and indels ( Figures 2B, D ). For SNVs, Mutect2, Strelka2, and VarDict produced the largest number of unfiltered SNVs (70,532, 58,033, and 36,060, respectively), while Muse, SomaticSniper, and VarScan2 called the smallest number of SNVs (10,019, 9,037, and 3,430, respectively) ( Figure 2B and Table S6 ). Moreover, Mutect2 and SomaticSniper had high rates of SNVs that did not satisfy the filtering criteria (non-PASS SNVs) (68.2% and 60.5%, respectively), whereas Muse, Strelka2, VarScan2, and VarDict had low rates of non-PASS SNVs (12.9%, 5.5%, 25.7%, and 37.1%, respectively). These non-PASS SNVs mainly consisted of “PoN” and “T_Alt < 4” in Mutect2 (39.9% and 15.1%, respectively) and “PoN” and “N_Alt >4” in SomaticSniper (46.0% and 13.2%, respectively) ( Figure 2C ). For indels, Mutect2 found more unfiltered indels compared to Strelka2, VarScan2, and VarDict (9,401, 37, 226, and 1,180, respectively) ( Figure 2D and Table S6 ). After filtering, high rates of non-PASS indels were found in Mutect2, whereas Strelka2, VarScan2, and VarDict had low rates of non-PASS indels (86.6%, 37.8%, 40.7% and 33.9%, respectively). These non-PASS indels mainly consisted of “PoN” and “T_Alt < 4” in Mutect2 (32.6% and 30.4%, respectively), “Clustered events” and “T_Alt < 4” in Strelka2 (13.5% and 10.8%, respectively), “PoN” in VarScan2 (24.8%), and “PoN” and “T_Alt < 4” in VarDict (17.4% and 14.7%, respectively) ( Figure 2E ). As such, approximately 23% (37,287/163,069) of unfiltered mutations consisted of PoN variants. Removing mutations marked with the “PoN” filter flag was a crucial step in reducing false-positive rates during WES analysis.

Muse and VarScan2 removed the lowest number of SNVs after filtering with filter flags and removing variants found by a single caller only ( Table S6 ). Strelka2 removed the lowest number of indels after variant filtering. Mutect2 removed the largest number of both SNVs and indels after filtering with filter flags and removing variants called by a single caller only. Although SomaticSniper removed a large number of the SNVs selection after filtering with filter flags, it removed a smaller number of SNVs after removing variants called by a single caller. Strelka2 recovered 94.5% of the SNVs after filtering with filter flags but only recovered 19.5% of the same after removing variants called by a single caller. Moreover, Strelka2 found a very small number of indels only.

In order to evaluate the performance of variant calling, IGV was used to visualize and examine somatic mutation in the six callers. The 20 most frequently mutated genes [mutated in at least 7/50 (14%) patients] were utilized for this analysis. Table S7 shows that the most frequent non-PASS genes during IGV examination were MUC16, MUC19, KMT2D, and TTN. To assess the reliability of the IGV examination, a number of mutations for each top 4 IGV non-PASS genes was selected for confirmation by Sanger sequencing. Figure S4 shows that both IGV-passed and -non-passed mutations could be confirmed by Sanger sequencing. The results of the IGV evaluation were also consistent with those of Sanger sequencing ( Table S8 ). In Table S7 , MUC16 and MUC19 genes had a large number of false-positive mutations [2100/2121 (99.0%) and 1900/1906 (99.7%), respectively] that did not satisfy the IGV filtering criteria. The IGV screenshot ( Figure S5 ) demonstrates that recurrent false-positive variants were observed at the MUC16 and MUC19 loci in our WES data. Therefore, MUC16 and MUC19 mutations were removed to reduce false-positive calls. A total of 2276 mutations were retained after removing mutations in MUC16 and MUC19, which then were used to evaluate the validation statistics of the mutation callers.

Table S9 shows that filtered mutations had significantly higher rates of IGV-PASS compared to unfiltered mutations (P < 0.001). The filtered mutations were also associated with a significantly higher rate of IGV-PASS in Mutect2 and SomaticSniper (both P < 0.001 and P < 0.001). An increased number of callers agreeing on the mutations was associated with an increase in IGV-PASS rates. However, in ≥2, ≥3, ≥4, ≥5 and 6 caller agreeing groups, filtered mutations significantly increased IGV-PASS rates compared to unfiltered mutations. Both filtering with filter flags and selection with mutation calling times were important procedures to reduce false-positive calls.

Furthermore, this study examined whether DP and MAF were associated with false-positive variant calls. Figure S6A shows the MAF and DP for unfiltered mutations in different mutation calling times. Accordingly, our results show that mutations with high MAF (MAF ≥ 0.8) and high DP (DP ≥ 1000) were mainly distributed in the one caller agreeing group. After filtering with filter flags, the high MAF and high DP mutations were significantly reduced in the one caller group ( Figure S6B ). IGV examination data shows that mutations with high DP (DP ≥ 1000) were mainly distributed in the IGV non-PASS group ( Figure S6C ). The aforementioned data suggest that read depths of 1000 and a MAF of 0.8 can be selected to eliminate false-positives mutations.

Receiver operating characteristic (ROC) curve was used to evaluate the optimal cut-off value of DP and MAF in distinguishing IGV-PASS mutations ( Figures S6E, G ). The AUCs for DP and MAF were 0.478 (95% CI: 0.459–0.496) and 0.607 (95% CI: 0.585–0.628), respectively. At a threshold of 175 reads for DP, the false positive rate and false negative rate were 51% and 45.1% in separating IGV-PASS mutations patients from normal IGV non-PASS mutations. At a threshold of 0.182 for MAF, the false positive rate and false negative rate were 38.5 and 86.3%, respectively. The Figures S6E, G showed the relationship between false positive rate and false negative rate with DP or MAF, respectively.

TMB was calculated as the number of somatic mutations in the coding region per Mb. In the TCGA-OSCC dataset, 71,890 non-synonymous mutations were identified in 387 patients (185.76 mutations per patient). In our unfiltered data, a total of 163,069 mutations were identified in 50 patients (3261.38 mutations per patient). Our findings show that the mean TMB was 35.81 mutations/Mb per patient. After filtering procedures, 13,730 mutations were retained (274.6 mutations per patient), with the TMB decreasing to 3.01 mutations/Mb per patient in filtered mutations. As depicted in Table 1 , high TMB was significantly associated with old age (P = 0.034) and advanced clinical stage (P = 0.021). The median TMB (0.96 mutations/Mb per patient) was considered the cut-off point for assessing outcomes in OSCC. Patients with higher TMB exhibited a poorer outcome compared to those with lower TMB (P = 0.041, Figure 3A ). In TCGA-OSCC dataset, patients with higher TMB also had a poorer outcome compared to those with lower TMB (P = 0.026, Figure 3B ).

CIRCOS plots were used to illustrate the mutational landscape in our study and in the TCGA-OSCC dataset ( Figure 4 ). The outer cycle illustrates the genomic profile from our WES data, with the most frequently mutated genes being TP53 (62%), FAT1 (40%), NOTCH1 (28%), and TTN (26%) ( Figure 4 ). The inner circle illustrates the mutation frequency of the genes in the TCGA-OSCC dataset, with the most mutations observed in TP53 (68%), TTN (42%) FAT1 (26%), and CDKN2A (22%). The mutational landscape in our study appears to be similar to that in the TCGA-OSCC dataset.

We found five novel mutated genes in at least six different patients (≥10%) included herein, which were not detected in the TCGA-OSCC dataset ( Table S10 ). The mutation frequencies of SKOR2, CPLANE1, CCDC168, STARD9, and NSD2 were 14%, 12%, 12%, 12%, and 10%, respectively. The InToGene platform and Bailey et al. data sets were used to predict potential mutational driver genes in OSCC. A total of 53 recognized mutational driver genes were found therein, among whom TP53 and FAT1 had high mutation rates ( Table S11 ). Table 2 shows the differences in the mutation frequency distribution between our most frequent genes (MAF ≥ 10%) and those in the TCGA-OSCC dataset. Accordingly, among our most frequent genes, 73 genes were identified to have high MAF, whereas only 1 had low MAF. Several of these genes we found were tumor suppressor genes, including FAT1, EPHA2, ASXL1, PTPRT, USP9X, IGF2R, SPTBN1, and PLCB3 (39–43, 45, 47, 48), whereas ARHGEF10L and NSD2 function as oncogenes (44, 55). TRRAP has been found to be involved in the regulation of stemness in ovarian cancer stem cells (35), while DNMT1 and EPHB1 can regulate tumor progression (38, 46). Other genes have been reported to been involved in the control of proliferation, invasion, and apoptosis in cancer cells (36, 37, 49–54, 56–59).

Furthermore, this study examined the relationship between non-synonymous mutation status and clinical parameters in the 20 most frequently mutated genes ( Figure 5 and Table S12 ). Accordingly, somatic mutation in ASXL1 was significantly associated with older age (P = 0.010) ( Table S12 ). In addition, FAT1 mutations were significantly associated with advanced clinical stage (P = 0.033) and marginally significantly associated with T stage (P = 0.050). Histological grade was significantly associated with PKD1L1 mutations (P = 0.044) and was marginally significantly associated with FAT3 mutations (P = 0.087). The presence of TTN mutations was significantly associated with perineural invasion (P = 0.038). Mutation of FMN2 was marginally significantly associated with T stage (P = 0.093) and histological grade mutations (P = 0.087). Survival analysis revealed that patients with CASP8 (log-rank P = 0.004), CUBN (log-rank P = 0.025), and USP9X mutations (log-rank P = 0.018) had significantly reduced OS rates ( Figures 6A–C ). Figure S7 showed that distribution of the 20 most frequently mutated genes and clinical features in the TCGA-OSCC patients. In the TCGA-HNSCC dataset, USP9X mutation was found to be associated with poor OS (log-rank P = 0.010, Figure S8C ). CASP8 mutation was marginal significantly associated with OS (log-rank P = 0.070, Figure S8C ). However, there were no significant differences in OS between CUBN mutation group and non-CUBN mutation group (log-rank P =0.885, Figure S8B ).

Pathway enrichment analysis of the 200 most frequently mutated genes was performed using the IPA. Accordingly, upstream regulator analysis in OSCC ( Figure S9 ) showed that AKT, TP53, and ERK were the most predicted upstream regulators controlling different gene clusters in OSCC. Table 3 shows that four canonical pathways had a P value < 0.05 using IPA, including the clathrin-mediated endocytosis signaling, NFκB signaling, PEDF signaling, and the calcium signaling pathways. Furthermore, we evaluated the association between OS and mutations of genes in canonical pathways. Among the patients included herein, 27 (54%) had a mutation in the NFκB signaling-related gene set, including CARD10, CASP8, EP300, FGFR1, IGF2R, and PIK3CA ( Figure 7A ). Mutations in CARD10, CASP8, EP300, FGFR1, IGF2R, and PIK3CA were observed in 5 (10%), 11 (22%), 6 (12%), 4 (8%), 6(12%), and 7 (14%) patients, respectively. Patients with mutations in the NFκB signaling-related gene set had poorer prognosis compared to those without such mutations ( Figure 6D ). However, no significant correlation between mutation status in NFκB signaling-related gene set and OS was found in TCGA-OSCC dataset ( Figure S8D ).

Mutations in the calcium signaling-related gene set were observed in 27 (54%) patients ( Figure 7B ). Mutations in the calcium signaling-related gene set, including ADCY2, PLCB1, PLCB3, ITPR1, ITPR2, and PIK3CA, were observed in 4 (8%), 4 (8%), 5 (10%), 4 (8%), 6 (12%), and 7 (14%) patients, respectively. Mutations in the calcium signaling-related gene set were associated with poor outcomes in OSCC ( Figure 6E ). In TCGA-OSCC dataset, there was no significant correlation of mutation status in calcium signaling-related gene set with OS ( Figure S8E ).

Thereafter, identified FDA-approved drugs associated with our candidate genes. Accordingly, seven candidate genes were found which may be targeted by FDA-approved drugs. They are involved in the regulation of the NFκB signaling-related and calcium signaling-related pathways ( Figure 7C ). These genes include PIK3CA and six receptor tyrosine kinases, namely FGFR1, FGFR2, FGFR3, EGFR, ALK, and ROS1 ( Figure 7C ).