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.

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/.

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.

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.

Epigenetic Regulations of Perineural Invasion in Head and Neck Squamous Cell Carcinoma.

Carcinomas of the oral cavity and oropharynx belong among the ten most common malignancies in the human population. The prognosis of head and neck squamous cell carcinoma (HNSCC) is determined by the degree of invasiveness of the primary tumor and by the extent of metastatic spread into regional and distant lymph nodes. Moreover, the level of the perineural invasion itself associates with tumor localization, invasion's extent, and the presence of nodal metastases. Here, we summarize the current knowledge about different aspects of epigenetic changes, which can be associated with HNSCC while focusing on perineural invasion (PNI). We review epigenetic modifications of the genes involved in the PNI process in HNSCC from the omics perspective and specific epigenetic modifications in OSCC or other neurotropic cancers associated with perineural invasion. Moreover, we summarize DNA methylation status of tumor-suppressor genes, methylation and demethylation enzymes and histone post-translational modifications associated with PNI. The influence of other epigenetic factors on the HNSCC incidence and perineural invasion such as tobacco, alcohol and oral microbiome is overviewed and HPV infection is discussed as an epigenetic factor associated with OSCC and related perineural invasion. Understanding epigenetic regulations of axon growth that lead to tumorous spread or uncovering the molecular control of axon interaction with cancer tissue can help to discover new therapeutic targets for these tumors.

Simultaneously targeting ErbB family kinases and PI3K in HPV-positive head and neck squamous cell carcinoma.

OBJECTIVES: To identify the most effective PI3K and EGFR inhibitors in HPV-positive head and neck squamous cell carcinoma (HNSCC) and investigate the efficacy of a combination of an ErbB family kinase inhibitor and a PI3K inhibitor to inhibit cell proliferation of HPV-positive HNSCC. MATERIALS AND METHOD: HPV-positive HNSCC cell lines were treated with the FDA approved ErbB kinase inhibitor, Afatinib or FDA-approved PI3K inhibitor, Copanlisib, alone or in combination, and phosphorylation and total protein levels of cells were assessed by Western blot analysis.Cell proliferation and apoptosis were examined by MTS assay, flow cytometry, and Western blots, respectively. RESULTS: Copanlisib more effectively inhibited cell proliferation in comparison to other PI3K inhibitors tested. HPV-positive HNSCC cells differentially responded to cisplatin, Afatinib, or Copanlisib. The combination of Afatinib and Copanlisib more effectively suppressed cell proliferation and induced apoptosis compared to either treatment alone. Mechanistically, the combination of Afatinib and Copanlisib completely blocked phosphorylation of EGFR, HER2, HER3, and Akt as well as significantly decreased the HPV E7 expression compared to either treatment alone. CONCLUSION: Afatinib and Copanlisib more effectively suppress cell proliferation and survival of HPV-positive HNSCC in comparison to either treatment alone.

Entinostat Decreases Immune Suppression to Promote Antitumor Responses in a HER2+ Breast Tumor Microenvironment.

Therapeutic combinations to alter immunosuppressive, solid tumor microenvironments (TME), such as in breast cancer, are essential to improve responses to immune checkpoint inhibitors (ICI). Entinostat, an oral histone deacetylase inhibitor, has been shown to improve responses to ICIs in various tumor models with immunosuppressive TMEs. The precise and comprehensive alterations to the TME induced by entinostat remain unknown. Here, we employed single-cell RNA sequencing on HER2-overexpressing breast tumors from mice treated with entinostat and ICIs to fully characterize changes across multiple cell types within the TME. This analysis demonstrates that treatment with entinostat induced a shift from a protumor to an antitumor TME signature, characterized predominantly by changes in myeloid cells. We confirmed myeloid-derived suppressor cells (MDSC) within entinostat-treated tumors associated with a less suppressive granulocytic (G)-MDSC phenotype and exhibited altered suppressive signaling that involved the NFκB and STAT3 pathways. In addition to MDSCs, tumor-associated macrophages were epigenetically reprogrammed from a protumor M2-like phenotype toward an antitumor M1-like phenotype, which may be contributing to a more sensitized TME. Overall, our in-depth analysis suggests that entinostat-induced changes on multiple myeloid cell types reduce immunosuppression and increase antitumor responses, which, in turn, improve sensitivity to ICIs. Sensitization of the TME by entinostat could ultimately broaden the population of patients with breast cancer who could benefit from ICIs.

Using DNA mechanics to predict in vitro nucleosome positions and formation energies.

In eukaryotic genomes, nucleosomes function to compact DNA and to regulate access to it both by simple physical occlusion and by providing the substrate for numerous covalent epigenetic tags. While competition with other DNA-binding factors and action of chromatin remodeling enzymes significantly affect nucleosome formation in vivo, nucleosome positions in vitro are determined by steric exclusion and sequence alone. We have developed a biophysical model, DNABEND, for the sequence dependence of DNA bending energies, and validated it against a collection of in vitro free energies of nucleosome formation and a set of in vitro nucleosome positions mapped at high resolution. We have also made a first ab initio prediction of nucleosomal DNA geometries, and checked its accuracy against the nucleosome crystal structure. We have used DNABEND to design both strong and weak histone- binding sequences, and measured the corresponding free energies of nucleosome formation. We find that DNABEND can successfully predict in vitro nucleosome positions and free energies, providing a physical explanation for the intrinsic sequence dependence of histone-DNA interactions.

Extracellular Vesicles Released by Tumor Endothelial Cells Spread Immunosuppressive and Transforming Signals Through Various Recipient Cells.

Head and neck squamous cell carcinoma (HNSCC) has a high recurrence and metastatic rate with an unknown mechanism of cancer spread. Tumor inflammation is the most critical processes of cancer onset, growth, and metastasis. We hypothesize that the release of extracellular vesicles (EVs) by tumor endothelial cells (TECs) induce reprogramming of immune cells as well as stromal cells to create an immunosuppressive microenvironment that favor tumor spread. We call this mechanism as non-metastatic contagious carcinogenesis. Extracellular vesicles were collected from primary HNSCC-derived endothelial cells (TEC-EV) and were used for stimulation of peripheral blood mononuclear cells (PBMCs) and primary adipose mesenchymal stem cells (ASCs). Regulation of ASC gene expression was investigated by RNA sequencing and protein array. PBMC, stimulated with TEC-EV, were analyzed by enzyme-linked immunosorbent assay and fluorescence-activated cell sorting. We validated in vitro the effects of TEC-EV on ASCs or PBMC by measuring invasion, adhesion, and proliferation. We found and confirmed that TEC-EV were able to change ASC inflammatory gene expression signature within 24-48 h. TEC-EV were also able to enhance the secretion of TGF-ß1 and IL-10 by PBMC and to increase T regulatory cell (Treg) expansion. TEC-EV carry specific proteins and RNAs that are responsible for Treg differentiation and immune suppression. ASCs and PBMC, treated with TEC-EV, enhanced proliferation, adhesion of tumor cells, and their invasion. These data indicate that TEC-EV exhibit a mechanism of non-metastatic contagious carcinogenesis that regulates tumor microenvironment and reprograms immune cells to sustain tumor growth and progression.

Expression of Programmed Cell Death Ligand 1 and Associated Lymphocyte Infiltration in Olfactory Neuroblastoma.

BACKGROUND: Programmed cell death ligand 1 (PD-L1) is a transmembrane glycoprotein that interacts with the receptor programmed cell death 1 (PD-1) to suppress T-cell activation, reduce adjacent tissue damage, and promote tolerance to self-antigens. Tumors may express PD-L1 as a mechanism to evade immune detection. Recent clinical trials have demonstrated the efficacy of PD-L1/PD-1 antagonists through activation of tumor-infiltrated CD8+ T cells. The aim of this study was to determine the expression pattern of PD-L1 and PD-1 in olfactory neuroblastoma (ONB) tumor cells and to determine the presence of PD-1+ and CD8+ lymphocytes in the ONB immune microenvironment. METHODS: Immunohistochemistry for expression of PD-L1, PD-1, and CD8 was performed on paraffin-embedded ONB tissue. RESULTS: Of the 10 primary site ONB samples, 4 demonstrated positive PD-L1 expression. Of PD-L1+ tumors, the 2 highest expressing samples were found to contain PD-1+ tumor cells. Of the 4 available metastatic samples, all of which arose from PD-L1- primary site ONB, 3 were positive for PD-L1 and contained PD-1+ tumor cells. PD-L1+ primary and metastatic tumors also demonstrated increased PD-1+ infiltrating lymphocytes in the tumor and stroma (11.6- and 4.62-fold increase) compared with PD-L1- samples (P < 0.05 and P = 0.068 respectively). PD-L1+ specimens demonstrated increased CD8+ lymphocytes in the tumor and stroma (7.46- and 2.14-fold increase) compared with PD-L1- tumors (P < 0.05 for both). CONCLUSIONS: These data demonstrate that a proportion of ONB primary and metastatic tumors express PD-L1 and possess an associated tumor and stromal infiltrate of PD-1+ and CD8+ lymphocytes.