The use of machine learning to discover regulatory networks controlling biological systems.

Biological systems are composed of a vast web of multiscale molecular interactors and interactions. High-throughput technologies, both bulk and single cell, now allow for investigation of the properties and quantities of these interactors. Computational algorithms and machine learning methods then provide the tools to derive meaningful insights from the resulting data sets. One such approach is graphical network modeling, which provides a computational framework to explicitly model the molecular interactions within and between the cells comprising biological systems. These graphical networks aim to describe a putative chain of cause and effect between interacting molecules. This feature allows for determination of key molecules in a biological process, accelerated generation of mechanistic hypotheses, and simulation of experimental outcomes. We review the computational concepts and applications of graphical network models across molecular scales for both intracellular and intercellular regulatory biology, examples of successful applications, and the future directions needed to overcome current limitations.

Insights into the mechanisms and structure of breakage-fusion-bridge cycles in cervical cancer using long-read sequencing.

Cervical cancer is caused by human papillomavirus (HPV) infection, has few approved targeted therapeutics, and is the most common cause of cancer death in low-resource countries. We characterized 19 cervical and four head and neck cancer cell lines using long-read DNA and RNA sequencing and identified the HPV types, HPV integration sites, chromosomal alterations, and cancer driver mutations. Structural variation analysis revealed telomeric deletions associated with DNA inversions resulting from breakage-fusion-bridge (BFB) cycles. BFB is a common mechanism of chromosomal alterations in cancer, and our study applies long-read sequencing to this important chromosomal rearrangement type. Analysis of the inversion sites revealed staggered ends consistent with exonuclease digestion of the DNA after breakage. Some BFB events are complex, involving inter- or intra-chromosomal insertions or rearrangements. None of the BFB breakpoints had telomere sequences added to resolve the dicentric chromosomes, and only one BFB breakpoint showed chromothripsis. Five cell lines have a chromosomal region 11q BFB event, with YAP1-BIRC3-BIRC2 amplification. Indeed, YAP1 amplification is associated with a 10-year-earlier age of diagnosis of cervical cancer and is three times more common in African American women. This suggests that individuals with cervical cancer and YAP1-BIRC3-BIRC2 amplification, especially those of African ancestry, might benefit from targeted therapy. In summary, we uncovered valuable insights into the mechanisms and consequences of BFB cycles in cervical cancer using long-read sequencing.

Matrix factorization and transfer learning uncover regulatory biology across multiple single-cell ATAC-seq data sets.

While the methods available for single-cell ATAC-seq analysis are well optimized for clustering cell types, the question of how to integrate multiple scATAC-seq data sets and/or sequencing modalities is still open. We present an analysis framework that enables such integration across scATAC-seq data sets by applying the CoGAPS Matrix Factorization algorithm and the projectR transfer learning program to identify common regulatory patterns across scATAC-seq data sets. We additionally integrate our analysis with scRNA-seq data to identify orthogonal evidence for transcriptional regulators predicted by scATAC-seq analysis. Using publicly available scATAC-seq data, we find patterns that accurately characterize cell types both within and across data sets. Furthermore, we demonstrate that these patterns are both consistent with current biological understanding and reflective of novel regulatory biology.

Integrated single-cell and bulk gene expression and ATAC-seq reveals heterogeneity and early changes in pathways associated with resistance to cetuximab in HNSCC-sensitive cell lines.

BACKGROUND: Identifying potential resistance mechanisms while tumour cells still respond to therapy is critical to delay acquired resistance. METHODS: We generated the first comprehensive multi-omics, bulk and single-cell data in sensitive head and neck squamous cell carcinoma (HNSCC) cells to identify immediate responses to cetuximab. Two pathways potentially associated with resistance were focus of the study: regulation of receptor tyrosine kinases by TFAP2A transcription factor, and epithelial-to-mesenchymal transition (EMT). RESULTS: Single-cell RNA-seq demonstrates heterogeneity, with cell-specific TFAP2A and VIM expression profiles in response to treatment and also with global changes to various signalling pathways. RNA-seq and ATAC-seq reveal global changes within 5 days of therapy, suggesting early onset of mechanisms of resistance; and corroborates cell line heterogeneity, with different TFAP2A targets or EMT markers affected by therapy. Lack of TFAP2A expression is associated with HNSCC decreased growth, with cetuximab and JQ1 increasing the inhibitory effect. Regarding the EMT process, short-term cetuximab therapy has the strongest effect on inhibiting migration. TFAP2A silencing does not affect cell migration, supporting an independent role for both mechanisms in resistance. CONCLUSION: Overall, we show that immediate adaptive transcriptional and epigenetic changes induced by cetuximab are heterogeneous and cell type dependent; and independent mechanisms of resistance arise while tumour cells are still sensitive to therapy.

Correction: Integrated single-cell and bulk gene expression and ATAC-seq reveals heterogeneity and early changes in pathways associated with resistance to cetuximab in HNSCC-sensitive cell lines.

An amendment to this paper has been published and can be accessed via a link at the top of the paper.

Splice Expression Variation Analysis (SEVA) for inter-tumor heterogeneity of gene isoform usage in cancer.

Motivation: Current bioinformatics methods to detect changes in gene isoform usage in distinct phenotypes compare the relative expected isoform usage in phenotypes. These statistics model differences in isoform usage in normal tissues, which have stable regulation of gene splicing. Pathological conditions, such as cancer, can have broken regulation of splicing that increases the heterogeneity of the expression of splice variants. Inferring events with such differential heterogeneity in gene isoform usage requires new statistical approaches. Results: We introduce Splice Expression Variability Analysis (SEVA) to model increased heterogeneity of splice variant usage between conditions (e.g. tumor and normal samples). SEVA uses a rank-based multivariate statistic that compares the variability of junction expression profiles within one condition to the variability within another. Simulated data show that SEVA is unique in modeling heterogeneity of gene isoform usage, and benchmark SEVA's performance against EBSeq, DiffSplice and rMATS that model differential isoform usage instead of heterogeneity. We confirm the accuracy of SEVA in identifying known splice variants in head and neck cancer and perform cross-study validation of novel splice variants. A novel comparison of splice variant heterogeneity between subtypes of head and neck cancer demonstrated unanticipated similarity between the heterogeneity of gene isoform usage in HPV-positive and HPV-negative subtypes and anticipated increased heterogeneity among HPV-negative samples with mutations in genes that regulate the splice variant machinery. These results show that SEVA accurately models differential heterogeneity of gene isoform usage from RNA-seq data. Availability and implementation: SEVA is implemented in the R/Bioconductor package GSReg. Contact: bahman@jhu.edu or favorov@sensi.org or ejfertig@jhmi.edu. Supplementary information: Supplementary data are available at Bioinformatics online.

CancerInSilico: An R/Bioconductor package for combining mathematical and statistical modeling to simulate time course bulk and single cell gene expression data in cancer.

Bioinformatics techniques to analyze time course bulk and single cell omics data are advancing. The absence of a known ground truth of the dynamics of molecular changes challenges benchmarking their performance on real data. Realistic simulated time-course datasets are essential to assess the performance of time course bioinformatics algorithms. We develop an R/Bioconductor package, CancerInSilico, to simulate bulk and single cell transcriptional data from a known ground truth obtained from mathematical models of cellular systems. This package contains a general R infrastructure for running cell-based models and simulating gene expression data based on the model states. We show how to use this package to simulate a gene expression data set and consequently benchmark analysis methods on this data set with a known ground truth. The package is freely available via Bioconductor: http://bioconductor.org/packages/CancerInSilico/.

Discovery and development of differentially methylated regions in human papillomavirus-related oropharyngeal squamous cell carcinoma.

Human papillomavirus (HPV)-related oropharyngeal squamous cell carcinoma (OPSCC) exhibits a different composition of epigenetic alterations. In this study, we identified differentially methylated regions (DMRs) with potential utility in screening for HPV-positive OPSCC. Genome wide DNA methylation was measured using methyl-CpG binding domain protein-enriched genome sequencing (MBD-seq) in 50 HPV-positive OPSCC tissues and 25 normal tissues. Fifty-one DMRs were defined with maximal methylation specificity to cancer samples. The Cancer Genome Atlas (TCGA) methylation array data was used to evaluate the performance of the proposed candidates. Supervised hierarchical clustering of 51 DMRs found that HPV-positive OPSCC had significantly higher DNA methylation levels compared to normal samples, and non-HPV-related head and neck squamous cell carcinoma (HNSCC). The methylation levels of all top 20 DNA methylation biomarkers in HPV-positive OPSCC were significantly higher than those in normal samples. Further confirmation using quantitative methylation specific PCR (QMSP) in an independent set of 24 HPV-related OPSCCs and 22 controls showed that 16 of the 20 candidates had significant higher methylation levels in HPV-positive OPSCC samples compared with controls. One candidate, OR6S1, had a sensitivity of 100%, while 17 candidates (KCNA3, EMBP1, CCDC181, DPP4, ITGA4, BEND4, ELMO1, SFMBT2, C1QL3, MIR129-2, NID2, HOXB4, ZNF439, ZNF93, VSTM2B, ZNF137P and ZNF773) had specificities of 100%. The prediction accuracy of the 20 candidates rang from 56.2% to 99.8% by receiver operating characteristic analysis. We have defined 20 highly specific DMRs in HPV-related OPSCC, which can potentially be applied to molecular-based detection tests and improve disease management.

PatternMarkers & GWCoGAPS for novel data-driven biomarkers via whole transcriptome NMF.

SUMMARY: Non-negative Matrix Factorization (NMF) algorithms associate gene expression with biological processes (e.g. time-course dynamics or disease subtypes). Compared with univariate associations, the relative weights of NMF solutions can obscure biomarkers. Therefore, we developed a novel patternMarkers statistic to extract genes for biological validation and enhanced visualization of NMF results. Finding novel and unbiased gene markers with patternMarkers requires whole-genome data. Therefore, we also developed Genome-Wide CoGAPS Analysis in Parallel Sets (GWCoGAPS), the first robust whole genome Bayesian NMF using the sparse, MCMC algorithm, CoGAPS. Additionally, a manual version of the GWCoGAPS algorithm contains analytic and visualization tools including patternMatcher, a Shiny web application. The decomposition in the manual pipeline can be replaced with any NMF algorithm, for further generalization of the software. Using these tools, we find granular brain-region and cell-type specific signatures with corresponding biomarkers in GTEx data, illustrating GWCoGAPS and patternMarkers ascertainment of data-driven biomarkers from whole-genome data. AVAILABILITY AND IMPLEMENTATION: PatternMarkers & GWCoGAPS are in the CoGAPS Bioconductor package (3.5) under the GPL license. CONTACT: gsteinobrien@jhmi.edu or ccolantu@jhmi.edu or ejfertig@jhmi.edu. SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

Structural analysis of nucleosomal barrier to transcription.

Thousands of human and Drosophila genes are regulated at the level of transcript elongation and nucleosomes are likely targets for this regulation. However, the molecular mechanisms of formation of the nucleosomal barrier to transcribing RNA polymerase II (Pol II) and nucleosome survival during/after transcription remain unknown. Here we show that both DNA-histone interactions and Pol II backtracking contribute to formation of the barrier and that nucleosome survival during transcription likely occurs through allosterically stabilized histone-histone interactions. Structural analysis indicates that after Pol II encounters the barrier, the enzyme backtracks and nucleosomal DNA recoils on the octamer, locking Pol II in the arrested state. DNA is displaced from one of the H2A/H2B dimers that remains associated with the octamer. The data reveal the importance of intranucleosomal DNA-protein and protein-protein interactions during conformational changes in the nucleosome structure on transcription. Mechanisms of nucleosomal barrier formation and nucleosome survival during transcription are proposed.

Characterization of functionally active gene fusions in human papillomavirus related oropharyngeal squamous cell carcinoma.

The Cancer Genome Atlas (TCGA) sequencing analysis of head and neck squamous cell carcinoma (HNSCC) recently reported on gene fusions, however, few human papillomavirus (HPV) positive samples were included, and the functional relevance of identified fusions was not explored. We therefore performed an independent analysis of gene fusions in HPV-positive oropharyngeal SCC (OPSCC). RNA sequencing was performed on 47 HPV-positive OPSCC primary tumors and 25 normal mucosal samples from cancer unaffected controls on an Illumina TruSeq platform. MapSplice2 was used for alignment and identification of fusion candidates. Putative fusions with less than five spanning reads, detected in normal tissues, or that mapped to the same gene were filtered out. Selected fusions were validated by RT-PCR and Sanger sequencing. Within 47 HPV-positive OPSCC tumors, 282 gene fusions were identified. Most fusions (85.1%) occurred in a single tumor, and the remaining fusions recurred in 2-16 tumors. Gene fusions were associated with significant up regulation of 16 genes (including EGFR and ERBB4) and down regulation of four genes (PTPRT, ZNF750, DLG2, SLCO5A1). Expression of these genes followed similar patterns of up regulation and down regulation in tumors without these fusions compared to normal tissue. Five of six gene fusions selected for validation were confirmed through RT-PCR and sequencing. This integrative analysis provides a method of prioritizing functionally relevant gene fusions that may be expanded to other tumor types. These results demonstrate that gene fusions may be one mechanism by which functionally relevant genes are altered in HPV-positive OPSCC.

NF-κB and stat3 transcription factor signatures differentiate HPV-positive and HPV-negative head and neck squamous cell carcinoma.

Using high-throughput analyses and the TRANSFAC database, we characterized TF signatures of head and neck squamous cell carcinoma (HNSCC) subgroups by inferential analysis of target gene expression, correcting for the effects of DNA methylation and copy number. Using this discovery pipeline, we determined that human papillomavirus-related (HPV+) and HPV- HNSCC differed significantly based on the activity levels of key TFs including AP1, STATs, NF-κB and p53. Immunohistochemical analysis confirmed that HPV- HNSCC is characterized by co-activated STAT3 and NF-κB pathways and functional studies demonstrate that this phenotype can be effectively targeted with combined anti-NF-κB and anti-STAT therapies. These discoveries correlate strongly with previous findings connecting STATs, NF-κB and AP1 in HNSCC. We identified five top-scoring pair biomarkers from STATs, NF-κB and AP1 pathways that distinguish HPV+ from HPV- HNSCC based on TF activity and validated these biomarkers on TCGA and on independent validation cohorts. We conclude that a novel approach to TF pathway analysis can provide insight into therapeutic targeting of patient subgroup for heterogeneous disease such as HNSCC.

Correcting transcription factor gene sets for copy number and promoter methylation variations.

Gene set analysis provides a method to generate statistical inferences across sets of linked genes, primarily using high-throughput expression data. Common gene sets include biological pathways, operons, and targets of transcriptional regulators. In higher eukaryotes, especially when dealing with diseases with strong genetic and epigenetic components such as cancer, copy number loss and gene silencing through promoter methylation can eliminate the possibility that a gene is transcribed. This, in turn, can adversely affect the estimation of transcription factor or pathway activity from a set of target genes, as some of the targets may not be responsive to transcriptional regulation. Here we introduce a simple filtering approach that removes genes from consideration if they show copy number loss or promoter methylation, and demonstrate the improvement in inference of transcription factor activity in a simulated dataset based on the background expression observed in normal head and neck tissue.

Radiomic Biomarkers of Locoregional Recurrence: Prognostic Insights from Oral Cavity Squamous Cell Carcinoma preoperative CT scans.

This study aimed to identify CT-based imaging biomarkers for locoregional recurrence (LR) in Oral Cavity Squamous Cell Carcinoma (OSCC) patients. Our study involved a retrospective review of 78 patients with OSCC who underwent surgical treatment at a single medical center. An approach involving feature selection and statistical model diagnostics was utilized to identify biomarkers. Two radiomics biomarkers, Large Dependence Emphasis (LDE) of the Gray Level Dependence Matrix (GLDM) and Long Run Emphasis (LRE) of the Gray Level Run Length Matrix (GLRLM) of the 3D Laplacian of Gaussian (LoG σ = 3), have demonstrated the capability to preoperatively distinguish patients with and without LR, exhibiting exceptional testing specificity (1.00) and sensitivity (0.82). The group with LRE > 2.99 showed a 3-year recurrence-free survival rate of 0.81, in contrast to 0.49 for the group with LRE ≤ 2.99. Similarly, the group with LDE > 120 showed a rate of 0.82, compared to 0.49 for the group with LDE ≤ 120. These biomarkers broaden our understanding of using radiomics to predict OSCC progression, enabling personalized treatment plans to enhance patient survival.

Clinical and genomic characterization of chemoradiation-resistant HPV-positive oropharyngeal squamous cell carcinoma.

Introduction: Most patients with HPV-positive oropharyngeal squamous cell carcinoma (OPSCC) have an excellent response to chemoradiation, and trials are now investigating de-escalated treatment. However, up to 25% of patients with HPV-positive OPSCC will experience recurrence, and up to 5% will even progress through primary treatment. Currently, there are no molecular markers to identify patients with poor prognosis who would be harmed by de-escalation. Herein we report the clinical and genomic characteristics of persistent HPV-positive OPSCC after definitive platinum-based chemoradiation therapy. Methods: Patients with HPV-positive OPSCC treated with curative intent platinum-based chemoradiation between 2007 and 2017 at two institutions and with a persistent locoregional disease were included. We evaluated clinical characteristics, including smoking status, age, stage, treatment, and overall survival. A subset of five patients had tissue available for targeted exome DNA sequencing and RNA sequencing. Genomic analysis was compared to a previously published cohort of 47 treatment-responsive HPV+ OPSCC tumors after batch correction. Mutational landscape, pathway activation, and OncoGPS tumor states were employed to characterize these tumors. Results: Ten patients met the inclusion criteria. The tumor and nodal stages ranged from T1 to T4 and N1 to N2 by AJCC 8th edition staging. All patients were p16-positive by immunohistochemistry, and eight with available in situ hybridization were confirmed to be HPV-positive. The 1-year overall survival from the time of diagnosis was 57%, and the 2-year overall survival was 17%. TP53 mutations were present in three of five (60%) persistent tumors compared to 2% (one of 47) of treatment-responsive HPV-positive tumors (p = 0.008). Other genes with recurrent mutations in persistent HPV-positive OPSCC tumors were NF1, KMT2D, PIK3C2B, and TFGBR2. Compared to treatment-responsive HPV-positive tumors, persistent tumors demonstrated activation of DNA Repair and p53, EMT, MYC, SRC, and TGF-beta signaling pathways, with post-treatment samples demonstrating significant activation of the PI3K-EMT-Stem pathways compared to pretreatment samples. Conclusion: Chemoradiation-resistant HPV-positive OPSCC occurs infrequently but portends a poor prognosis. These tumors demonstrate higher rates of p53 mutation and activation of MYC, SRC, and TGF-beta pathways. A comparison of tumors before and after treatment demonstrates PI3K-EMT-Stem pathways post-treatment in HPV-positive tumors with persistent disease after platinum-based chemoradiation.

Radiomic biomarkers of locoregional recurrence: prognostic insights from oral cavity squamous cell carcinoma preoperative CT scans.

Introduction: This study aimed to identify CT-based imaging biomarkers for locoregional recurrence (LR) in Oral Cavity Squamous Cell Carcinoma (OSCC) patients. Methods: Computed tomography scans were collected from 78 patients with OSCC who underwent surgical treatment at a single medical center. We extracted 1,092 radiomic features from gross tumor volume in each patient's pre-treatment CT. Clinical characteristics were also obtained, including race, sex, age, tobacco and alcohol use, tumor staging, and treatment modality. A feature selection algorithm was used to eliminate the most redundant features, followed by a selection of the best subset of the Logistic regression model (LRM). The best LRM model was determined based on the best prediction accuracy in terms of the area under Receiver operating characteristic curve. Finally, significant radiomic features in the final LRM model were identified as imaging biomarkers. Results and discussion: Two radiomics biomarkers, Large Dependence Emphasis (LDE) of the Gray Level Dependence Matrix (GLDM) and Long Run Emphasis (LRE) of the Gray Level Run Length Matrix (GLRLM) of the 3D Laplacian of Gaussian (LoG σ=3), have demonstrated the capability to preoperatively distinguish patients with and without LR, exhibiting exceptional testing specificity (1.00) and sensitivity (0.82). The group with LRE > 2.99 showed a 3-year recurrence-free survival rate of 0.81, in contrast to 0.49 for the group with LRE ≤ 2.99. Similarly, the group with LDE > 120 showed a rate of 0.82, compared to 0.49 for the group with LDE ≤ 120. These biomarkers broaden our understanding of using radiomics to predict OSCC progression, enabling personalized treatment plans to enhance patient survival.

A polar barrier to transcription can be circumvented by remodeler-induced nucleosome translocation.

Many eukaryotic genes are regulated at the level of transcript elongation. Nucleosomes are likely targets for this regulation. Previously, we have shown that nucleosomes formed on very strong positioning sequences (601 and 603), present a high, orientation-dependent barrier to transcription by RNA polymerase II in vitro. The existence of this polar barrier correlates with the interaction of a 16-bp polar barrier signal (PBS) with the promoter-distal histone H3-H4 dimer. Here, we show that the polar barrier is relieved by ISW2, an ATP-dependent chromatin remodeler, which translocates the nucleosome over a short distance, such that the PBS no longer interacts with the distal H3-H4 dimer, although it remains within the nucleosome. In vivo, insertion of the 603 positioning sequence into the yeast CUP1 gene results in a modest reduction in transcription, but this reduction is orientation-independent, indicating that the polar barrier can be circumvented. However, the 603-nucleosome is present at the expected position in only a small fraction of cells. Thus, the polar barrier is probably non-functional in vivo because the nucleosome is not positioned appropriately, presumably due to nucleosome sliding activities. We suggest that interactions between PBSs and chromatin remodelers might have significant regulatory potential.

Identification of CT-based non-invasive Radiographic Biomarkers for Overall Survival Stratification in Oral Cavity Squamous Cell Carcinoma.

This study addresses the limited non-invasive tools for Oral Cavity Squamous Cell Carcinoma OSCC survival prediction by identifying Computed Tomography (CT)-based biomarkers for improved prognosis. A retrospective analysis was conducted on data from 149 OSCC patients, including radiomics and clinical. An ensemble approach involving correlation analysis, score screening, and the Sparse-L1 algorithm was used to select functional features, which were then used to build Cox Proportional Hazards models (CPH). Our CPH achieved a 0.70 concordance index in testing. The model identified two CT-based radiomics features, Gradient-Neighboring-Gray-Tone-Difference-Matrix-Strength (GNS) and normalized-Wavelet-LLL-Gray-Level-Dependence-Matrix-Large-Dependence-High-Gray-Level-Emphasis (HLE), as well as smoking and alcohol usage, as survival biomarkers. The GNS group with values above 14 showed a hazard ratio of 0.12 and a 3-year survival rate of about 90%. Conversely, the GNS group with values less than or equal to 14 had a 49% survival rate. For normalized HLE, the high-end group (HLE > -0.415) had a hazard ratio of 2.41, resulting in a 3-year survival rate of 70%, while the low-end group (HLE <= -0.415) had a 36% survival rate. These findings contribute to our knowledge of how radiomics can be used to anticipate the outcome and tailor treatment plans from people with OSCC.

Identification of CT-based non-invasive radiomic biomarkers for overall survival prediction in oral cavity squamous cell carcinoma.

This study addresses the limited non-invasive tools for Oral Cavity Squamous Cell Carcinoma (OSCC) survival prediction by identifying Computed Tomography (CT)-based biomarkers to improve prognosis prediction. A retrospective analysis was conducted on data from 149 OSCC patients, including CT radiomics and clinical information. An ensemble approach involving correlation analysis, score screening, and the Sparse-L1 algorithm was used to select functional features, which were then used to build Cox Proportional Hazards models (CPH). Our CPH achieved a 0.70 concordance index in testing. The model identified two CT-based radiomics features, Gradient-Neighboring-Gray-Tone-Difference-Matrix-Strength (GNS) and normalized-Wavelet-LLL-Gray-Level-Dependence-Matrix-Large-Dependence-High-Gray-Level-Emphasis (HLE), as well as stage and alcohol usage, as survival biomarkers. The GNS group with values above 14 showed a hazard ratio of 0.12 and a 3-year survival rate of about 90%. Conversely, the GNS group with values less than or equal to 14 had a 49% survival rate. For normalized HLE, the high-end group (HLE > - 0.415) had a hazard ratio of 2.41, resulting in a 3-year survival rate of 70%, while the low-end group (HLE ≤ - 0.415) had a 36% survival rate. These findings contribute to our knowledge of how radiomics can be used to predict the outcome so that treatment plans can be tailored for patients people with OSCC to improve their survival.

Enhancing oral squamous cell carcinoma prediction: the prognostic power of the worst pattern of invasion and the limited impact of molecular resection margins.

Objective: Oral squamous cell carcinoma (OSCC) originates from the mucosal lining of the oral cavity. Almost half of newly diagnosed cases are classified as advanced stage IV disease, which makes resection difficult. In this study, we investigated the pathological features and mutation profiles of tumor margins in OSCC. Methods: We performed hierarchical clustering of principal components to identify distinct patterns of tumor growth and their association with patient prognosis. We also used next-generation sequencing to analyze somatic mutations in tumor and marginal tissue samples. Results: Our analyses uncovered that the grade of worst pattern of invasion (WPOI) is strongly associated with depth of invasion and patient survival in multivariable analysis. Mutations were primarily detected in the DNA isolated from tumors, but several mutations were also identified in marginal tissue. In total, we uncovered 29 mutated genes, mainly tumor suppressor genes involved in DNA repair including BRCA genes; however none of these mutations significantly correlated with a higher chance of relapse in our medium-size cohort. Some resection margins that appeared histologically normal harbored tumorigenic mutations in TP53 and CDKN2A genes. Conclusion: Even histologically normal margins may contain molecular alterations that are not detectable by conventional histopathological methods, but NCCN classification system still outperforms other methods in the prediction of the probability of disease relapse.