Integrated analysis of oral tongue squamous cell carcinoma identifies key variants and pathways linked to risk habits, HPV, clinical parameters and tumor recurrence

Habits, clinical parameters and epidemiology

We collected tumor and matched control (adjacent normal and/or lymphocytes) samples from 50 patients diagnosed with OTSCC, with informed consent. Data from patient habits, epidemiology and clinical parameters are presented in Figure 1A and Additional file 1A. About two-thirds of the patients (N = 31) included in our study were in the younger age group (≤50yrs), with 20% female patients in the total pool. Approximately, 70% of the patients were positive for at least one risk habit, namely, smoking, alcohol consumption or chewing tobacco (33% of patients smoked tobacco, 40% consumed alcohol and 42% chewed tobacco). HPV infection status in the primary tumors was established with type-specific qPCR or HPV16 digital PCR. Thirty-three percent of the patients were deceased at the time of completing the analysis. About 60% of the tumors were moderately differentiated, 25% well differentiated and the rest were poorly differentiated. Among the patients recruited, 60% were node-positive, 70% had no recurrence, 9% had distant metastasis and 24% had loco-regional recurrence at the time of completing the analysis. The mean and median follow-up durations for patients were nearly 30 months and 21 months, respectively. About 27% of the tumors were early stage tumors (T1N0M0 and T2N0M0) and the rest 73% were late stage tumors (tumors belonging to the rest of the TNM stage).

Figure 1. Key variants in OTSCC and their relationship with habits, clinical and epidemiological parameters.

A. The OTSCC samples are represented in color-codes with their corresponding status on; node (P: positive, N: negative); stage (E: early, L: late), recurrence (Y: loco-regionally recurrent, N: non-recurrent and M: distant metastatic); grade (WD: well-differentiated, MD: moderately-differentiated and PD: poorly-differentiated); disease-free survival or DFS (L: low/≤12mo, M: mid/12–24mo and H: high/>24mo); HPV (P: positive and N: negative); and habits (chewing, alcohol and smoking, Y: yes and N: no). B. Somatic mutation frequency per megabase (MB) is represented as scatterplot with the median point as a fine dotted line. C. Genes with significant somatic variants. D. Frequency histogram of nineteen cancer-associated genes bearing somatic missense and nonsense variants (mutations and indels). E. Columns representing mutually exclusive sets of genes. F. Significant copy number insertions and deletions (CNVs), alongside the chromosome cytobands (the numbers of cancer-associated genes within each cytoband are listed on the right).

Discovery and validation of significant somatic variants and their relationship with other parameters

We re-discovered variants, as described previously¹² using whole-genome arrays, to validate the variant call accuracy as obtained from the exome sequencing data. We validated ~99% of the SNPs discovered from Illumina sequencing in both the tumor and matched control samples (Additional file 2). After filtering and annotation, we identified 19 cancer-associated genes bearing significantly altered somatic variants in OTSCC (Figure 1D). These were validated using Sanger sequencing in two sets of samples, one using the same tumor-control pairs used in the exome sequencing (the discovery set, Additional file 1A) and second, using an additional 36–60 primary tumors (validation set, Additional file 1B) for genes altered in ≥ 5% of the tumor samples. All the TP53 variants were validated in the discovery set. Three out of the four variants were validated for CASP8. The mutant alleles for the heterozygous variants in HLA-A, OBSCN, ING1, TTK and U2AF1 discovered by exome sequencing were difficult to interpret from the results of the validation using Sanger sequencing as they were present at a very low frequency (Additional file 3). Combining data from the validation set; the mutation frequencies for RASA1 and CDKN2A rose significantly to 10.71% and 16.47% in primary tumors respectively but those for TP53 and CASP8 remained largely unchanged (Additional file 3).

The somatic mutation frequency per megabase (MB) ranged from 10–45 with a median around 25 (Figure 1B). The median value for transition to transversion (ti/tv) ratio for both the tumor and its matched control samples was ~2.5 (Additional file 4). Overall, T->C changes were most frequent, followed by G->A and then T->G. Habits (smoking and alcohol consumption), nodal status, HPV infection, tumor grade and stage had no significant impact on the distribution of these nucleotides (Additional file 5). We used the workflow described in the Methods section to identify somatic mutations and indels in tumor samples following which we used three functional tools, IntOGen¹⁹, MutSigCV²¹ and MuSiC2²² for variant interpretations (Additional file 6). In order to identify genes harboring significant variants, we used the intersection of these tools, following the criteria that the somatic variants be callable in the matched control sample and present in a single sequencing read in the control sample. This resulted in a final list of 19 cancer-associated genes (Figure 1C), which were divided into three categories with varying mutation frequencies (Figure 1D). The three frequency tiers were ≥ 30% (TP53), 6–30% (RASA1, CASP8 and CDKN2A) and 2–5% (NOTCH1, NOTCH2, DMD and PIK3CA were prominent among them).

Next, we looked for mutual exclusivity of finding somatic variants in the genes and found that many of these genes harbor variants in a mutually exclusive manner across samples (Figure 1E), suggesting the possibility that there might be some common pathway(s) involved in the development of OTSCC. We observed mutual exclusivity among somatic variants in NOTCH1 and NOTCH2 genes, and expanded this finding to identifying 15 such mutually exclusive sets (Figure 1E). Among them, CDKN2A, HLA-A and TTK form a mutually exclusive set with TP53; RASA1, OBSCN, HLA-A, AJUBA and TTK are mutually exclusive with either NOTCH1 alone, or NOTCH2 and ANK3 together; NOTCH1, NOTCH2, HLA-A, AJUBA, ANK3, TTK, MLL2, ING1 or KEAP1, are mutually exclusive with CASP8 alone, or FAT1 and DMD together; FAT1, HLA-A, AJUBA, ANK3, TTK, MLL2, ING1 or KEAP1, are mutually exclusive with PIK3CA or DMD or NOTCH1 and OBSCN, or CDKN2A and OBSCN; U2AF1, MLL2 and TTK form a small mutually exclusive set. We juxtaposed the positions of the somatic variants from final list of all 19 genes (Additional file 7) detected in OTSCC against those found in the TCGA data using the cBioPortal. We found that the somatic variants in OTSCC were in the same domains where mutations were observed earlier in many of the genes (Additional file 7).

CNV analyses using data from the whole-genome single nucleotide polymorphism (SNP) genotyping arrays revealed a large chunk of chromosome 9, bearing cancer-associated genes like CDKN2A, NF1 and MRPL4, to be affected in about 17% of the tumors (Figure 1F and Additional file 8). We found several CNVs of short stretches (in low kb range) within chromosomes 6–8, 11, 17 and X in many tumors.

Linking habits, HPV infection, nodal status, tumor grade and recurrence, with genes harboring somatic variants and the associated pathways

We further classified the 19 cancer-associated genes from the previous analyses and linked those with habits, clinical parameters and HPV infection. Among the genes harboring significant somatic variants, we found CDKN2A to be mutated only in the never-smokers and past smokers, PIK3CA to be mutated only in the smokers, and TP53 to be mutated at a 20% greater frequency in the smokers, CASP8 mutation has a 12% greater frequency in those that consumed alcohol or chewed tobacco. RASA1 was exclusively mutated only in the non-chewers (Figure 2A). HPV-negative patients harbored somatic variants in DMD and PIK3CA, while HPV-positive patients alone had somatic variants in RASA1. Only the moderate- and well-differentiated tumor samples harbored variants in CASP8, while NOTCH1 was mutated largely in the poorly-differentiated tumors. Node-positive tumors had a 19% greater occurrence of TP53 variants. Somatic variants in RASA1 occurred exclusively in the late stage tumors (Figure 2A). We further studied the association of affected cancer-related signaling pathways with habits and clinical parameters, and found that recurrence and HPV infection had the highest impact (Figure 2B). The Procaspase-8 activation, Notch, p53 and Wnt signaling pathways were linked most with many of the clinical parameters, HPV infection and habits (Figure 2B).

Figure 2. Relationship between genes harboring somatic variants with clinical-, epidemiological parameters and signaling pathways.

A. Histograms showing relationship between genes with significant somatic variants and various clinical and epidemiological parameters. For genes solely mutated in one of the clinical or epidemiological categories, or those mutated at a >= 5% frequency between two categories. B. Stack net charts of relative patient fraction (%) for each of the eight cancer-associated signaling pathways and their relationship with various clinical and epidemiological parameters.

Differentially expressed genes in OTSCC

Significant (q val ≤ 0.05) differentially expressed genes with a fold change of at least 1.5 revealed a consistent pattern of differential expression across the tumor samples (21 up- and 23 down-regulated genes, Figure 3A and Additional file 9). Genes involved in peroxisome proliferator-activated receptor (PPAR) signaling- (e.g., MMP1) and ECM-receptor interaction pathways (LAMC2 and SPP1) were up-regulated and CRNN, APOD, SCARA5 and RERGL were down-regulated in a majority of tumors (Figure 3A). Next, we studied the pathways involving genes with aberrant expression and their link with HPV infection and other clinical parameters. Genes in the arachidonic acid metabolism and Toll-like receptors were differentially expressed in patients with no smoking history (never smokers or past smokers) and alcohol habits (Figure 3B). SERPINE1 (a gene in HIF-1 signaling pathway) was differentially expressed in patients that are habits-negative. The NF-κ-B signaling pathway was differentially expressed only in metastasized tumors.

Figure 3. Differentially expressed genes, affected pathways and their relationship with clinical and epidemiological parameters.

A. Expression changes (green – up- and red – down-regulation) representing significantly differentially expressed genes in tumors. B. Stacked histograms representing relative patient fraction (%) for each of the 19 cancer-associated pathways and their relationship with clinical and epidemiological parameters.

Functional studies with CASP8 in OTSCC cell lines

CASP8 is mutated in a significant number of oral tongue tumors [this study, 5, 7]. Caspase-8 is an important and versatile protein that plays a role in both apoptotic (extrinsic or death receptor-mediated) and non-apoptotic processes^13,14. We studied the functional consequences of CASP8 knockdown through a siRNA-mediated method in an HPV-positive UM:SCC-47¹⁵ and an HPV-negative UPCI:SCC040¹⁶ OTSCC cell lines. Prior to the functional assay, the concentration of siRNA required for silencing, extent of CASP8 knockdown and cisplatin sensitivity (IC₅₀) in both these cell lines was tested (Additional file 10). The invasion of cells was greater in both UM:SCC-47 and UPCI:SCC040 cell lines when CASP8 was knocked down (Figure 4A). To analyze the effect of caspase-8 on the migration property of cells, scratches were made on the confluent monolayer of cells and the wound closure area was measured at different time points (0hr, 15hr, 23hr & 42hr, Figure 4B). The wound closure was faster in CASP8 knockdown HPV-negative cells compared to the HPV-positive cells. At 15hr, 23hr and 48hrs, about 65%, 90% and 100% of the wound got closed respectively in the HPV-negative cell line compared to 50%, 70% and 85% respectively during the same time period in the HPV-positive cell lines (Figure 4B). siRNA knockdown of CASP8 rescued the chemo-sensitivity caused by cisplatin treatment as evident by the MTT (3-(4,5-dimethylthiazol-2-yl)-2,5-diphenyltetrazolium bromide) survival assay (Figure 4C). Interestingly, we found that the extent of rescue is greater in the HPV-negative cell line compared to the HPV16-positive one.

Figure 4. Role of CASP8 in HPV-positive and HPV-negative OTSCC cell lines.

Results from A. Matrigel cell invasion assay (plotted with respect to the control cells), B. Wound healing assay, and C. MTT cell survival assay (plotted with respect to the control cells) in UPCI:SCC040 (HPV-negative) and UM:SCC-47 (HPV-positive) cell lines.

Tumor recurrence prediction using random forests

After cataloging the significantly altered genes in OTSCC, we wanted to see whether there is a relationship between the altered genes and loco-regional recurrence of tumors and metastasis. In order to do this, we used an ensemble machine-learning method implemented by variable elimination using random forests¹¹ (Figure 5). We used multiple testing correction and the 0.632 bootstrapping method¹⁷ to estimate false positives. We discovered a 38-gene minimal signature that discriminated between the non-recurring, loco-regionally recurring and distant metastatic tumors (Figure 5). The .632+ bootstrap errors, indicative of prediction specificity, varied across non-recurrent, recurrent and distant metastatic tumors. The median error was low (0.03) and intermediate (0.3) for the non-recurrent and the loco-regionally recurrent categories respectively but was relatively higher (1.0) for the metastatic tumors. The errors were proportional to the number of representative samples within each category.

Figure 5. A minimal gene signature for tumor recurrence.

Genes harboring somatic variants (in color) that are a part of the minimal signature set for tumor recurrence derived from random forest analyses are used.

Major signaling pathways implicated in OTSCC

We looked at significant pathways altered in OTSCC, taking into account all the molecular changes in tumors and found apoptosis, HIF, Notch, mTOR, p53, PI3K/Akt, Wnt and Ras to be some of the key signaling pathways affected in OTSCC (Figure 6). In addition, histone methylation, cell cycle/immunity and mRNA splicing processes were also affected. The complete list of pathways is provided in Additional file 11.

Figure 6. Significantly affected pathways in OTSCC.

Genes harboring significant somatic variants and with expression changes in tumors were used in Cytoscape to derive a set of important signaling pathways implicated in OTSCC.

Informed consent, ethics approval and patient samples used in the study

Informed consent was obtained voluntarily from each patient enrolled in the study. Ethical approval (NHH/MEC-CL/2014/197) was obtained from the Institutional Ethics Committees of the Mazumdar Shaw Medical Centre. Matched control (blood and/or adjacent normal tissue) and tumor specimens were collected and used in the study. Patients diagnosed and treated at the cancer clinic of the Mazumdar Shaw Medical Centre for oral tongue tumors were subjected to a screening procedure before being enrolled in the study. Only those patients, where the histological sections confirmed the presence of squamous cell carcinoma with at least 70% tumor cells in the specimen, were used in the current study. At the time of admission, patients were asked about the habits (chewing, smoking and/or alcohol consumption). Fifty treatment-naïve patients who underwent staging according to AJCC criteria, and curative intent treatment as per NCCN guideline involving surgery with or without post-operative adjuvant radiation or chemo-radiation at the Mazumdar Shaw Medical Centre were accrued for the study (Additional file 1). Post-treatment surveillance was carried out by clinical and radiographic examinations as per the NCCN guidelines.

HPV detection

HPV was detected by using q-PCR (Applied Biosystems 9700) using HPV16- and HPV18-specific TaqMan probes and primers, and digital PCR (BioRad QX100) using TaqMan probes and primers to detect HPV in primary tumor samples. The primers, probes and cycling conditions for q-PCR and ddPCR were as follows. For q-PCR: 5’ GCA CAG AGC TGC AAA CAA CT 3’; 3’ GCA TAA ATC CCG AAA AGC AA 5’; probe-ATTAGAATGTGTGTACTGCAAGCA-FAM-BHQ and 5’ TGA CAC TGT GCC TCA ATC CT 3’; 3’ AGA GCC ACT TGG AGA GGG AG 5’; Probe-TGCCTGCTTCACCTGGCAGC-VIC-BHQ for HPV16 and HPV18 respectively. The cycling conditions for q-PCR were: 95°C : 3 min, 95°C : 30 sec, 55°C for HPV16 and 60°C for HPV18 : 30 sec, 72°C : 30 sec for 40 cycles. For ddPCR: 5’ ACT GTC AAA AGC CAC TGT GT 3’; 3’ GCT GGG TTT CTC TAC GTG TT 5’ and Probe-AGGGGTCGGTGGACCGGTCGATGT-FAM-BHQ for HPV16. The cycling conditions for ddPCR were: 95°C: 10 min, 95°C : 315 sec, 55°C : 20 sec for 40 cycles.

Exome sequencing, read QC, alignment, variant discovery and post-processing filters

Exome libraries were prepared using Agilent SureSelect, Illumina TruSeq and Nextera exome capture kits (Additional file 14) following manufacturers’ specifications. Paired end sequencing was performed using HiSeq 2500 or GAIIx and raw reads were generated using standard Illumina base caller (HCS 2.0). Read pairs were filtered using in house scripts (Additional file 15 and Additional file 16) and only those reads having ≥75% bases with ≥ 20 phred score and ≤ 15 Ns were used for sequence alignment against human hg19 reference genome using NovoAlign (v3.00.05)²⁷. The aligned files (*.sam) were processed using Samtools (v0.1.12a)²⁸ and only uniquely mapped reads from NovoAlign were considered for variant calling. The alignments were pre-processed using GATK (v1.2-62)²⁹ in three steps before variant calling. First, the indels were realigned using the known indels from 1000G (phase1) data. Second, duplicates were removed using Picard (v1.39). Third, base quality recalibration was done using CountCovariates and TableRecalibration from GATK (v1.2-62), taking into account known SNPs and indels from dbSNP (build 138). Finally, UnifiedGenotyper from GATK (v2.5-2) was used for variant calling, using known SNPs and indels from dbSNP (build 138). Raw variants from GATK were filtered to only include the PASS variants (standard call confidence ≥ 50) within the merged exomic bait boundaries. Two out of 50 tumor samples did not confirm to the QC standards, therefore excluded from all further analyses. Therefore, all the downstream analyses were restricted to 48 primary tumors. The variants were further flagged as novel or present in either dbSNP138 or COSMIC (v67) databases, based on their overlap. In addition to GATK, we also used Dindel³⁰ to call indels. Both GATK and Dindel calls were filtered for microsatellite repeats (flagged as STR). The raw variant calls were used to estimate frequencies of nucleotide changes and transition:transversion (ti/tv) ratios. Exome-filtered PASS variants specific to the tumor samples, with respect to both location and actual call, were retained as somatic variants, which were further filtered to exclude variants where the region bearing the variant was not callable in the matched control sample, and those where the matched control sample had even one read covering the variant allele.

Scripts used to perform various filtering steps are provided in Additional file 16. The numbers of raw reads, after QC, alignment statistics, numbers of variants pre- and post-filters are provided in Additional file 2.

Detection of cross-contamination and identification of significant somatic variants

We estimated cross-contamination using ContEst (June 2013)³¹ in the tumor samples (Additional file 16). Locus-wise and gene-wise driver scores were estimated by CRAVAT³² using the head and neck cancer database with the CHASM³³ analysis option. Genes with a CHASM score of at least 0.35 were considered significant for comparison with other functional analyses (Additional file 16). Somatic mutations were normalized with respect to the exome bait size (MB) to calculate the somatic mutation frequency per MB.

Annotation and functional analyses of variants

Annotation and functional analyses of somatic variants was performed using IntoGen (web version 2.4)^34,35, MutSigCV (v1.3.01)^36,37 and MuSiC2 (v0.1)³⁷. Somatic variants, filtered to contain only those callable in the matched normal but not covered by any read in the control samples (VCF), were used for IntoGen with the 'cohort analyses' option. We also ran MutsigCV1.3 with these variants using coverage from un-filtered variants of all tumor samples (Additional file 16). Pooled alignments for all normal and tumor samples (BAM), each, along with pooled variants for all normal samples (MAF) were analyzed using MuSiC2 to calculate the background mutation rates (bmrs) for all genes, and identify a list of significantly mutated genes (p-value of convolution test ≤ 0.05; Additional file 16). A condensed list of 19 genes, common between at least two analyses was compiled (Figure 1D).

SNP genotyping and validation using Illumina whole-genome Omni LCG arrays

High quality DNA (200ng), quantified by Qubit 2.0 (Invitrogen), was used as the starting material for whole-genome genotyping experiments following the manufacturer’s specifications. Briefly, the genomic DNA was denatured at room temperature (RT) for 10 mins using 0.1N NaOH, neutralized and used for whole genome amplification (WGA) under isothermal conditions, at 37°C for 20 hrs. Post-WGA, the DNA was enzymatically fragmented at 37°C for 1hr. The fragmented DNA was precipitated with isopropanol at 4°C and resuspended in hybridization buffer. The samples were then denatured at 95°C for 20 mins, cooled at RT for 30 mins and 35µl of each sample was loaded onto the Illumina HumanOmni 2.5-8 beadchip for hybridization (20hrs at 48°C) in a hybridization chamber. The unhybridized probes were washed away and the Chips (HumanOmni 2.5-8 v1.0 and v1.1, Additional file 2) were prepared for staining, single base extension and scanning using Illumina’s HiScan system.

We filtered the SNP locations to retain only those, called without any error, contained within the exome boundaries as per the sequencing baits, and which were callable (covered by at least five sequencing reads). At these locations, we estimated the overlap for individual SNP calls, i.e., chr/pos/ref/alt and for no calls; i.e., chr/pos/ref/ref; between sequencing and array platforms (Additional file 16).

Discovering Copy number Variations (CNVs) and Loss of Heterozygosity (LOH)

CNVs and LOHs were identified using cnvPartition 3.1.6 plugin in Illumina GenomeStudio v2011.1, with default settings except for a minimum coverage of at least 10 probes per CNV/LOH with a confidence score threshhold of at least 100 (Additional file 17). Somatic CNVs and LOHs were extracted by filtering out any region common to CNVs and LOHs detected in its matched control. Somatic CNVs and LOHs were further filtered with respect to common and disease-related CNVs and LOHs using CNVAnnotator³⁸. Overlaps with common CNVs and LOHs were discarded, reporting only the overlaps with disease-related, and novel CNVs and LOHs. We categorized the CNVs and LOHs within each cytoband and reported those with an occurrence in at least 10% of the patient samples.

Gene expression assay

Gene expression profiling was carried out using Illumina HumanHT-12 v4 expression BeadChip (Illumina, San Diego, CA) in tumor and matched normal tissues (Additional file 9) following manufacturer’s specifications. Total RNA was extracted from 20mg of tissue using PureLink RNA (Invitrogen) and RNeasy (Qiagen) Mini kits. RNA quality was checked using Agilent Bioanalyzer 2100 using RNA Nano6000 chip. Samples with poor RNA integrity numbers (RIN) (<7), indicating partial degradation of RNA, were processed using Illumina WGDASL assay as per manufacturer's recommendations. The RNA samples with no degradation were labelled using Illumina TotalPrep RNA Amplification kit (Ambion) and processed according to the array manufacturer’s recommendations. Gene expression data was collected using Illumina’s HiScan and analyzed with the GenomeStudio (v2011.1 Gene Expression module 1.9.0) and all assay controls were checked to ensure quality of the assay and chip scanning. Raw signal intensities were exported from GenomeStudio for pre-processing and analyzed using R further.

Gene-wise expression intensities for tumor and matched control samples from GenomeStudio were transformed and normalized using VST (Variance Stabilizing Transformation) and LOESS methods, respectively, using the R package lumi³⁹. The data was further batch-corrected using ComBat⁴⁰ (Additional file 16). The pre-processed intensities for tumor and matched control samples were subjected to differential expression analyses using the R package, limma⁴¹ (Additional file 16). Genes with significant expression changes (adjusted P value <= 0.05) and fold change of at least 1.5 were followed up with further functional analyses.

Recurrence prediction using random forests

We used presence or absence of somatic mutations/indels data in the entire set of genes for all the OTSCC patients, along with their recurrence patterns as training set for the random forests¹¹ analyses using the varSelRF package in R. This method performs both backward elimination of variables and selection based on their importance spectrum, and predicts recurrence patterns in the same set by iteratively eliminating 2% of the least important predictive variables until the current OOB (out-of-bag) error rate becomes larger than the initial or previous OOB error rates. In order to understand the specificity of the best minimalistic predictors of tumor recurrence, we estimated the 0.632+ error rate¹⁷ over 50 bootstrap replicates. We used the varSelRFBoot function from the varSelRF Bioconductor package to perform bootstrapping. The .632+ method is described by the following formula:

$$Er{r}^{{0.632}^{\prime }}=Er{r}^{0.632}+\frac{\left(Er{r}^1-err\right)\left(.368\cdot .632\cdot R^{\prime }\right)}{1-.368\cdot R^{\prime }}$$

where Err^(.632'), Err^(.632), Err⁽¹⁾ and err are errors estimated by the .632+ method, the original .632 method, leave-one-out bootstrap method and err represents the error. R’ represents a value between 0 and 1. Another popular error correction method used is leave-one-out bootstrap method. The .632+ method was designed to correct the upward bias in the leave-one-out and the downward bias in the original .632 bootstrap methods.

For all iterations of all random forest analyses, we confirmed that the variable importance remained the same before and after correcting for multiple hypotheses comparisons using pre- and post- Benjamin-Hochberg FDR-corrected P values. R commands for variable elimination using random forests, 0.632+ bootstrapping and re-computing importance values after multiple comparisons testing are provided in Additional file 16.

Pathway analyses

Consensus list of genes from analysis, filtering and annotation of variant calls and from differential expression analysis using whole genome micro-arrays, were mapped to pathways using the web version of Graphite Web⁴² employing KEGG and Reactome databases. The network of interactions between genes was drawn originally using CytoScape (v3.1.1)⁴³ using the .sif file created by Graphite Web (Additional file 16).

Data visualization

We used Circos (v0.66)⁴⁴ (Additional file 18 and Additional file 19) for multi-dimensional data visualization. Additionally, we used the cbioportal protal (http://www.cbioportal.org/) to visualize variants within the 19 genes harboring significant variants. All of the mandatory fields accepted by Mutation Mapper were provided for select genes from our study to create structural representations for each gene including domains. Such diagrams from our study, the HNSCC study and all cancer studies from TCGA were collated using the image-editing tool, GIMP (v2.8.0) (www.gimp.org). SNPs and indels were visualized for each individual tumor sample using IGV (v1.5.54)⁴⁵, along with the reads supporting variants (Additional file 20).

Validation of somatic variants using Sanger sequencing

Primers were designed using the NCBI primer designing tool (http://www.ncbi.nlm.nih.gov/tools/primer-blast/index.cgi?LINK_LOC=BlastHome) and used in Sanger sequencing for validation. The sequences of all primers (IDT) used for validation is provided in Additional file 21. We tested the specificity of the designed primers using UCSC's tool, In Silico PCR. The variant-bearing region was amplified by using specific primers and used in Sanger sequencing (Additional file 14). The somatic variants were confirmed by sequencing in the entire tumor and matched control DNA set used for the exome sequencing followed by further validation in 60 additional tumor samples (Additional file 1B).

Cell culture and knockdown of CASP8 gene

The human OTSCC cell lines UPCI:SCC040 (gift from Dr. Susan Gollin, University of Pittsburgh, PA, USA)¹⁶ and UM-SCC47 (gift from Dr. Thomas Carey, University of Michigan, MI, USA)¹⁵ were used in the study. All the cells were maintained in Dulbecco’s Modified Eagles’ Media (DMEM) supplemented with 10% FBS, 1% MEM nonessential amino acids solution & 1% penicillin/streptomycin mixture (Gibco) at 37°C with 5% CO₂ incubator.

We performed the siRNA-based knockdown using UPCI:SCC040 and UM:SCC47 cell lines for CASP8 gene. The expression of Caspase-8 was transiently knocked down using ON-TARGETplus Human CASP8 smart pool siRNA (L-003466-00-0010; Dharmacon) along with an ON-TARGETplus Non-targeting siRNA (D-001810-01-20; Dharmacon). The transfection efficiency for the two cell lines (UPCI:SCC040 and UM:SCC47) were optimized using siGLO Red Transfection Indicator (D-001630; Dharmacon). The siRNA duplexes were transfected using Lipofectamine-2000 according to the manufacturer’s instructions (Invitrogen). The siRNA-oligo complexes medium was changed 8 hrs post transfection. The efficiency of transfection along with the mRNA expression was analyzed at 24 and 48 hrs post transfection by qRT-PCR. The specific down-regulation of CASP8 was confirmed by three independent experiments.

RNA isolation and quantitative real-time PCR

RNA was extracted from cell pellets and tissues using RNeasy Mini kit spin columns (Qiagen) following manufacturer’s protocol. Genomic DNA contamination was removed by RNase-Free DNase Set (Qiagen) and the total RNA was eluted in nuclease free water (Ambion). The RNA samples were estimated using Qubit 2.0 fluorometer (Invitrogen) and the integrity was checked by gel electrophoresis. The RNA samples were stored at -80°C until further used. The cDNA was synthesized with 400ng total RNA, using a SuperScript-III first strand cDNA synthesis kit, and following the manufacturer’s instructions (Invitrogen). The cDNA was then subjected for quantitative real-time PCR (q-RT-PCR) using KAPA SYBR FAST qPCR Master Mix (KK4601, KAPA). The primer pairs used for testing the expression of caspase-8 in q-RT-PCR were, forward 5'-ATGATGACATGAACCTGCTGGA-3' and reverse 5’-CAGGCTCTTGTTGATTTGGGC-3'. The amplification was done on Stratagene MX300P real time machine. The cycling conditions were: step-1 95°C-3min, step-2 95°C-3sec, step-3 60°C-60sec then repeat steps 2–3 for 40 cycles following dissociation curve at 60°C-60sec, 95°C-1min, 60°C-60sec.

To normalize inter-sample variation in RNA input, the expression values were normalized with GAPDH. All amplification reactions were done in triplicates, using nuclease free water as negative controls. The differential gene expression was calculated by using the comparative C_T method of relative quantification⁴⁶.

Assessment of cell viability (confirm)

MTT cell proliferation assay was performed as per manufacturer’s instructions (Sigma) to assess cell viability. Briefly, cells were seeded on 96-well plates containing DMEM with 10% FBS & incubated overnight. After treatment with 0.1% DMSO (vehicle control), or Cisplatin for 48 hrs, medium was changed and 100 µl of MTT solution (1mg/ml) was added to each well. The cells were further incubated for 4hrs at 37°C. The formazan crystals in viable cells were dissolved by adding 100µl of dimethyl sulfoxide (DMSO) (Merck). The absorbance was recorded at 540 nm using reference wavelength of 690 nm on micro plate reader (Tecan Systems). Data were normalized to vehicle treatment, and the cell viability was calculated using GraphPad Prism software (version 4.03; La Jolla, CA). All the experiments were performed in triplicates.

Wound healing assay

Cells were cultured up to 80% confluency in 12 well plates; serum-starved for 24 hrs and then wounded using a 200µl pipette tip. The wound was washed with 1× PBS and the cells were grown in DMEM containing 10% FBS. Cells were imaged at 10× magnification at 0 hr, 15 hrs, 23 hrs and 42 hrs. For each well, three wounds were made and the migration distance was photographed and measured using Carl Zeiss software (Zeiss). Each experiment was performed in triplicates.

Matrigel invasion assay

The ECM gel (E1270, Sigma) was thawed overnight at 4°C and plated at requisite concentrations (for UPCI:SCC040: 1.5mg/ml and UMSCC047: 2mg/ml) onto the transwell inserts and incubated overnight in the CO₂ incubator at 37°C with 5% CO₂. Cells were serum-starved for overnight, harvested, counted and seeded (UPCI:SCC040: 50,000 cells and UMSCC047: 20,000 cells per well) on top of the matrigel transwell-inserts (2 mg/ml) in serum-free medium as per manufacturer’s specifications (Sigma). D-MEM containing 10% FBS and 1% NEAA was added to the lower chamber. The 24-well plates containing matrigel inserts with cells were incubated in 37°C incubator for 48 hrs. At the end of incubation time, cells in the upper chamber were removed with cotton swabs and cells that invaded the Matrigel to the lower surface of the insert were fixed with 4% paraformaldehyde (Merk Milipore), permeabilized with 100% methanol, stained with Giemsa (Sigma), mounted on glass slides with DPX mounting agent and counted under a light microscope (Zeiss). Each experiment was performed in triplicates.