QChromosomeVisualizer: A new tool for 3D visualization of long simulations of polymer-like chromosome models.

Recent years have brought us great wealth of new types of experimental data on different aspects of chromatin state, from chromosome conformation assays, through super-resolution microscopic imaging to epigenetic modifications and lamina interaction assays. This rapid increase in data availability have motivated many novel approaches to 3D modeling of chromosomes, their conformations and dynamic behavior. Even though there are many tools already developed for molecular visualization in the field of structural bioinformatics, they are usually optimized for visualization of smaller molecules (like proteins) and much shorter trajectories. We have developed a novel approach to visualization of long trajectories of large polymers, typical in the field of chromatin modeling. Our software, called QChromosomeVisualizer (QCV), allows for quick visualization of long simulations containing thousands or even millions of frames and generating good looking still images and animations including spherical 360 videos that can be viewed in VR headsets. We believe that this kind of tools will be helpful for the broader community of researchers interested in modeling by allowing them to create new and clearer ways to communicate their results.

HERON: A Novel Tool Enables Identification of Long, Weakly Enriched Genomic Domains in ChIP-seq Data.

The explosive development of next-generation sequencing-based technologies has allowed us to take an unprecedented look at many molecular signatures of the non-coding genome. In particular, the ChIP-seq (Chromatin ImmunoPrecipitation followed by sequencing) technique is now very commonly used to assess the proteins associated with different non-coding DNA regions genome-wide. While the analysis of such data related to transcription factor binding is relatively straightforward, many modified histone variants, such as H3K27me3, are very important for the process of gene regulation but are very difficult to interpret. We propose a novel method, called HERON (HiddEn MaRkov mOdel based peak calliNg), for genome-wide data analysis that is able to detect DNA regions enriched for a certain feature, even in difficult settings of weakly enriched long DNA domains. We demonstrate the performance of our method both on simulated and experimental data.

K-mer Content Changes with Node Degree in Promoter-Enhancer Network of Mouse ES Cells.

Maps of Hi-C contacts between promoters and enhancers can be analyzed as networks, with cis-regulatory regions as nodes and their interactions as edges. We checked if in the published promoter-enhancer network of mouse embryonic stem (ES) cells the differences in the node type (promoter or enhancer) and the node degree (number of regions interacting with a given promoter or enhancer) are reflected by sequence composition or sequence similarity of the interacting nodes. We used counts of all k-mers (k = 4) to analyze the sequence composition and the Euclidean distance between the k-mer count vectors (k-mer distance) as the measure of sequence (dis)similarity. The results we obtained with 4-mers are interpretable in terms of dinucleotides. Promoters are GC-rich as compared to enhancers, which is known. Enhancers are enriched in scaffold/matrix attachment regions (S/MARs) patterns and depleted of CpGs. Furthermore, we show that promoters are more similar to their interacting enhancers than vice-versa. Most notably, in both promoters and enhancers, the GC content and the CpG count increase with the node degree. As a consequence, enhancers of higher node degree become more similar to promoters, whereas higher degree promoters become less similar to enhancers. We confirmed the key results also for human keratinocytes.

Automated inference of gene regulatory networks using explicit regulatory modules.

Gene regulatory networks are a popular tool for modelling important biological phenomena, such as cell differentiation or oncogenesis. Efficient identification of the causal connections between genes, their products and regulating transcription factors, is key to understanding how defects in their function may trigger diseases. Modelling approaches should keep up with the ever more detailed descriptions of the biological phenomena at play, as provided by new experimental findings and technical improvements. In recent years, we have seen great improvements in mapping of specific binding sites of many transcription factors to distinct regulatory regions. Recent gene regulatory network models use binding measurements; but usually only to define gene-to-gene interactions, ignoring regulatory module structure. Moreover, current huge amount of transcriptomic data, and exploration of all possible cis-regulatory arrangements which can lead to the same transcriptomic response, makes manual model building both tedious and time-consuming. In our paper, we propose a method to specify possible regulatory connections in a given Boolean network, based on transcription factor binding evidence. This is implemented by an algorithm which expands a regular Boolean network model into a "cis-regulatory" Boolean network model. This expanded model explicitly defines regulatory regions as additional nodes in the network, and adds new, valuable biological insights to the system dynamics. The expanded model can automatically be compared with expression data. And, for each node, a regulatory function, consistent with the experimental data, can be found. The resulting models are usually more constrained (by biologically-motivated metadata), and can then be inspected in in silico simulations. The fully automated method for model identification has been implemented in Python, and the expansion algorithm in R. The method resorts to the Z3 Satisfiability Modulo Theories (SMT) solver, and is similar to the RE:IN application (Yordanov et al., 2016). It is available on https://github.com/regulomics/expansion-network.

Endothelial cell differentiation is encompassed by changes in long range interactions between inactive chromatin regions.

Endothelial cells (ECs) differentiate from mesodermal progenitors during vasculogenesis. By comparing changes in chromatin interactions between human umbilical vein ECs, embryonic stem cells and mesendoderm cells, we identified regions exhibiting EC-specific compartmentalization and changes in the degree of connectivity within topologically associated domains (TADs). These regions were characterized by EC-specific transcription, binding of lineage-determining transcription factors and cohesin. In addition, we identified 1200 EC-specific long-range interactions (LRIs) between TADs. Most of the LRIs were connected between regions enriched for H3K9me3 involving pericentromeric regions, suggesting their involvement in establishing compartmentalization of heterochromatin during differentiation. Second, we provide evidence that EC-specific LRIs correlate with changes in the hierarchy of chromatin aggregation. Despite these rearrangements, the majority of chromatin domains fall within a pre-established hierarchy conserved throughout differentiation. Finally, we investigated the effect of hypoxia on chromatin organization. Although hypoxia altered the expression of hundreds of genes, minimal effect on chromatin organization was seen. Nevertheless, 70% of hypoxia-inducible genes situated within a TAD bound by HIF1α suggesting that transcriptional responses to hypoxia largely depend on pre-existing chromatin organization. Collectively our results show that large structural rearrangements establish chromatin architecture required for functional endothelium and this architecture remains largely unchanged in response to hypoxia.

Arabidopsis SWI/SNF chromatin remodeling complex binds both promoters and terminators to regulate gene expression.

ATP-dependent chromatin remodeling complexes are important regulators of gene expression in Eukaryotes. In plants, SWI/SNF-type complexes have been shown critical for transcriptional control of key developmental processes, growth and stress responses. To gain insight into mechanisms underlying these roles, we performed whole genome mapping of the SWI/SNF catalytic subunit BRM in Arabidopsis thaliana, combined with transcript profiling experiments. Our data show that BRM occupies thousands of sites in Arabidopsis genome, most of which located within or close to genes. Among identified direct BRM transcriptional targets almost equal numbers were up- and downregulated upon BRM depletion, suggesting that BRM can act as both activator and repressor of gene expression. Interestingly, in addition to genes showing canonical pattern of BRM enrichment near transcription start site, many other genes showed a transcription termination site-centred BRM occupancy profile. We found that BRM-bound 3΄ gene regions have promoter-like features, including presence of TATA boxes and high H3K4me3 levels, and possess high antisense transcriptional activity which is subjected to both activation and repression by SWI/SNF complex. Our data suggest that binding to gene terminators and controlling transcription of non-coding RNAs is another way through which SWI/SNF complex regulates expression of its targets.

Coordinated expression and genetic polymorphisms in Grainyhead-like genes in human non-melanoma skin cancers.

BACKGROUND: The Grainyhead-like (GRHL) transcription factors have been linked to many different types of cancer. However, no previous study has attempted to investigate potential correlations in expression of different GRHL genes in this context. Furthermore, there is very little information concerning damaging mutations and/or single nucleotide polymorphisms in GRHL genes that may be linked to cancer. METHODS: DNA and RNA were extracted from human non-melanoma skin cancers (NMSC) and adjacent normal tissues (n = 33 pairs of samples). The expression of GRHL genes was measured by quantitative real time PCR. Regulation of GRHL expression by miRNA was studied using cell transfection methods and dual-luciferase reporter system. Targeted deep sequencing of GRHL genes in tumor samples and control tissues were employed to search for mutations and single nucleotide polymorphisms. Single marker rs141193530 was genotyped with pyrosequencing in additional NMSC replication cohort (n = 176). Appropriate statistical and bioinformatic methods were used to analyze and interpret results. RESULTS: We discovered that the expression of two genes - GRHL1 and GRHL3 - is reduced in a coordinated manner in tumor samples, in comparison to the control healthy skin samples obtained from the same individuals. It is possible that both GRHL1 and GRHL3 are regulated, at least to some extent, by different strands of the same oncogenic microRNA - miR-21, what would at least partially explain observed correlation. No de novo mutations in the GRHL genes were detected in the examined tumor samples. However, some single nucleotide polymorphisms in the GRHL genes occur at significantly altered frequencies in the examined group of NMSC patients. CONCLUSIONS: Non-melanoma skin cancer growth is accompanied by coordinated reduced expression of epidermal differentiation genes: GRHL1 and GRHL3, which may be regulated by miR-21-3p and -5p, respectively. Some potentially damaging single nucleotide polymorphisms in GRHL genes occur with altered frequencies in NMSC patients, and they may in particular impair the expression of GRHL3 gene or functioning of encoded protein. The presence of these polymorphisms may indicate an increased risk of NMSC development in affected people.

Coordinated repression and activation of two transcriptional programs stabilizes cell fate during myogenesis.

Molecular models of cell fate specification typically focus on the activation of specific lineage programs. However, the concurrent repression of unwanted transcriptional networks is also essential to stabilize certain cellular identities, as shown in a number of diverse systems and phyla. Here, we demonstrate that this dual requirement also holds true in the context of Drosophila myogenesis. By integrating genetics and genomics, we identified a new role for the pleiotropic transcriptional repressor Tramtrack69 in myoblast specification. Drosophila muscles are formed through the fusion of two discrete cell types: founder cells (FCs) and fusion-competent myoblasts (FCMs). When tramtrack69 is removed, FCMs appear to adopt an alternative muscle FC-like fate. Conversely, ectopic expression of this repressor phenocopies muscle defects seen in loss-of-function lame duck mutants, a transcription factor specific to FCMs. This occurs through Tramtrack69-mediated repression in FCMs, whereas Lame duck activates a largely distinct transcriptional program in the same cells. Lineage-specific factors are therefore not sufficient to maintain FCM identity. Instead, their identity appears more plastic, requiring the combination of instructive repressive and activating programs to stabilize cell fate.

Potential protective role of Grainyhead-like genes in the development of clear cell renal cell carcinoma.

The involvement of Grainyhead-like (GRHL) transcription factors in various cancers is well documented. However, little is known about their role in clear cell renal cell carcinoma (ccRCC). We discovered that the expression of two of these factors-GRHL1 and GRHL2-are downregulated in ccRCC samples, and their expression is correlated with the expression of VHL gene. This suggests a functional link between the GRHL transcription factors and one of the best known tumor suppressors. Although the GRHL genes are not mutated in ccRCC, some of the single nucleotide polymorphisms in these genes may indicate an increased risk of ccRCC development and/or may allow to assess patients' prognoses and predict their responses to various forms of therapy. Silencing of GRHL2 expression in non-tumorigenic kidney cell line results in increased cell proliferation, increased resistance to apoptosis, as well as changes in the levels of selected proteins involved in the pathogenesis of ccRCC. These changes support the potential role for GRHL2 as a suppressor of ccRCC.

A Specialized Histone H1 Variant Is Required for Adaptive Responses to Complex Abiotic Stress and Related DNA Methylation in Arabidopsis.

Linker (H1) histones play critical roles in chromatin compaction in higher eukaryotes. They are also the most variable of the histones, with numerous nonallelic variants cooccurring in the same cell. Plants contain a distinct subclass of minor H1 variants that are induced by drought and abscisic acid and have been implicated in mediating adaptive responses to stress. However, how these variants facilitate adaptation remains poorly understood. Here, we show that the single Arabidopsis (Arabidopsis thaliana) stress-inducible variant H1.3 occurs in plants in two separate and most likely autonomous pools: a constitutive guard cell-specific pool and a facultative environmentally controlled pool localized in other tissues. Physiological and transcriptomic analyses of h1.3 null mutants demonstrate that H1.3 is required for both proper stomatal functioning under normal growth conditions and adaptive developmental responses to combined light and water deficiency. Using fluorescence recovery after photobleaching analysis, we show that H1.3 has superfast chromatin dynamics, and in contrast to the main Arabidopsis H1 variants H1.1 and H1.2, it has no stable bound fraction. The results of global occupancy studies demonstrate that, while H1.3 has the same overall binding properties as the main H1 variants, including predominant heterochromatin localization, it differs from them in its preferences for chromatin regions with epigenetic signatures of active and repressed transcription. We also show that H1.3 is required for a substantial part of DNA methylation associated with environmental stress, suggesting that the likely mechanism underlying H1.3 function may be the facilitation of chromatin accessibility by direct competition with the main H1 variants.

Optimally choosing PWM motif databases and sequence scanning approaches based on ChIP-seq data.

BACKGROUND: For many years now, binding preferences of Transcription Factors have been described by so called motifs, usually mathematically defined by position weight matrices or similar models, for the purpose of predicting potential binding sites. However, despite the availability of thousands of motif models in public and commercial databases, a researcher who wants to use them is left with many competing methods of identifying potential binding sites in a genome of interest and there is little published information regarding the optimality of different choices. Thanks to the availability of large number of different motif models as well as a number of experimental datasets describing actual binding of TFs in hundreds of TF-ChIP-seq pairs, we set out to perform a comprehensive analysis of this matter. RESULTS: We focus on the task of identifying potential transcription factor binding sites in the human genome. Firstly, we provide a comprehensive comparison of the coverage and quality of models available in different databases, showing that the public databases have comparable TFs coverage and better motif performance than commercial databases. Secondly, we compare different motif scanners showing that, regardless of the database used, the tools developed by the scientific community outperform the commercial tools. Thirdly, we calculate for each motif a detection threshold optimizing the accuracy of prediction. Finally, we provide an in-depth comparison of different methods of choosing thresholds for all motifs a priori. Surprisingly, we show that selecting a common false-positive rate gives results that are the least biased by the information content of the motif and therefore most uniformly accurate. CONCLUSION: We provide a guide for researchers working with transcription factor motifs. It is supplemented with detailed results of the analysis and the benchmark datasets at http://bioputer.mimuw.edu.pl/papers/motifs/ .

BNFinder2: Faster Bayesian network learning and Bayesian classification.

SUMMARY: Bayesian Networks (BNs) are versatile probabilistic models applicable to many different biological phenomena. In biological applications the structure of the network is usually unknown and needs to be inferred from experimental data. BNFinder is a fast software implementation of an exact algorithm for finding the optimal structure of the network given a number of experimental observations. Its second version, presented in this article, represents a major improvement over the previous version. The improvements include (i) a parallelized learning algorithm leading to an order of magnitude speed-ups in BN structure learning time; (ii) inclusion of an additional scoring function based on mutual information criteria; (iii) possibility of choosing the resulting network specificity based on statistical criteria and (iv) a new module for classification by BNs, including cross-validation scheme and classifier quality measurements with receiver operator characteristic scores. AVAILABILITY AND IMPLEMENTATION: BNFinder2 is implemented in python and freely available under the GNU general public license at the project Web site https://launchpad.net/bnfinder, together with a user's manual, introductory tutorial and supplementary methods.

Effective modeling of the chromatin structure by coarse-grained methods.

The interphase chromatin structure is extremely complex, precise and dynamic. Experimental methods can only show the frequency of interaction of the various parts of the chromatin. Therefore, it is extremely important to develop theoretical methods to predict the chromatin structure. In this publication, we implemented an extended version of the SBS model described by Barbieri et al. and created the ChroMC program that is easy to use and freely available (https://github.com/regulomics/chroMC) to other users. We also describe the necessary factors for the effective modeling of the chromatin structure in Drosophila melanogaster. We compared results of chromatin structure predictions using two methods: Monte Carlo and Molecular Dynamic. Our simulations suggest that incorporating black, non-reactive chromatin is necessary for successful prediction of chromatin structure, while the loop extrusion model with a long range attraction potential or Lennard-Jones (with local attraction force) as well as using Hi-C data as input are not essential for the basic structure reconstruction. We also proposed a new way to calculate the similarity of the properties of contact maps including the calculation of local similarity.Communicated by Ramaswamy H. Sarma.

Predicting spatial and temporal gene expression using an integrative model of transcription factor occupancy and chromatin state.

Precise patterns of spatial and temporal gene expression are central to metazoan complexity and act as a driving force for embryonic development. While there has been substantial progress in dissecting and predicting cis-regulatory activity, our understanding of how information from multiple enhancer elements converge to regulate a gene's expression remains elusive. This is in large part due to the number of different biological processes involved in mediating regulation as well as limited availability of experimental measurements for many of them. Here, we used a Bayesian approach to model diverse experimental regulatory data, leading to accurate predictions of both spatial and temporal aspects of gene expression. We integrated whole-embryo information on transcription factor recruitment to multiple cis-regulatory modules, insulator binding and histone modification status in the vicinity of individual gene loci, at a genome-wide scale during Drosophila development. The model uses Bayesian networks to represent the relation between transcription factor occupancy and enhancer activity in specific tissues and stages. All parameters are optimized in an Expectation Maximization procedure providing a model capable of predicting tissue- and stage-specific activity of new, previously unassayed genes. Performing the optimization with subsets of input data demonstrated that neither enhancer occupancy nor chromatin state alone can explain all gene expression patterns, but taken together allow for accurate predictions of spatio-temporal activity. Model predictions were validated using the expression patterns of more than 600 genes recently made available by the BDGP consortium, demonstrating an average 15-fold enrichment of genes expressed in the predicted tissue over a naïve model. We further validated the model by experimentally testing the expression of 20 predicted target genes of unknown expression, resulting in an accuracy of 95% for temporal predictions and 50% for spatial. While this is, to our knowledge, the first genome-wide approach to predict tissue-specific gene expression in metazoan development, our results suggest that integrative models of this type will become more prevalent in the future.

Challenges for modeling global gene regulatory networks during development: insights from Drosophila.

Development is regulated by dynamic patterns of gene expression, which are orchestrated through the action of complex gene regulatory networks (GRNs). Substantial progress has been made in modeling transcriptional regulation in recent years, including qualitative "coarse-grain" models operating at the gene level to very "fine-grain" quantitative models operating at the biophysical "transcription factor-DNA level". Recent advances in genome-wide studies have revealed an enormous increase in the size and complexity or GRNs. Even relatively simple developmental processes can involve hundreds of regulatory molecules, with extensive interconnectivity and cooperative regulation. This leads to an explosion in the number of regulatory functions, effectively impeding Boolean-based qualitative modeling approaches. At the same time, the lack of information on the biophysical properties for the majority of transcription factors within a global network restricts quantitative approaches. In this review, we explore the current challenges in moving from modeling medium scale well-characterized networks to more poorly characterized global networks. We suggest to integrate coarse- and find-grain approaches to model gene regulatory networks in cis. We focus on two very well-studied examples from Drosophila, which likely represent typical developmental regulatory modules across metazoans.

Dynamic CRM occupancy reflects a temporal map of developmental progression.

Development is driven by tightly coordinated spatio-temporal patterns of gene expression, which are initiated through the action of transcription factors (TFs) binding to cis-regulatory modules (CRMs). Although many studies have investigated how spatial patterns arise, precise temporal control of gene expression is less well understood. Here, we show that dynamic changes in the timing of CRM occupancy is a prevalent feature common to all TFs examined in a developmental ChIP time course to date. CRMs exhibit complex binding patterns that cannot be explained by the sequence motifs or expression of the TFs themselves. The temporal changes in TF binding are highly correlated with dynamic patterns of target gene expression, which in turn reflect transitions in cellular function during different stages of development. Thus, it is not only the timing of a TF's expression, but also its temporal occupancy in refined time windows, which determines temporal gene expression. Systematic measurement of dynamic CRM occupancy may therefore serve as a powerful method to decode dynamic changes in gene expression driving developmental progression.

Careful feature selection is key in classification of Alzheimer's disease patients based on whole-genome sequencing data.

Despite great increase of the amount of data from genome-wide association studies (GWAS) and whole-genome sequencing (WGS), the genetic background of a partially heritable Alzheimer's disease (AD) is not fully understood yet. Machine learning methods are expected to help researchers in the analysis of the large number of SNPs possibly associated with the disease onset. To date, a number of such approaches were applied to genotype-based classification of AD patients and healthy controls using GWAS data and reported accuracy of 0.65-0.975. However, since the estimated influence of genotype on sporadic AD occurrence is lower than that, these very high classification accuracies may potentially be a result of overfitting. We have explored the possibilities of applying feature selection and classification using random forests to WGS and GWAS data from two datasets. Our results suggest that this approach is prone to overfitting if feature selection is performed before division of data into the training and testing set. Therefore, we recommend avoiding selection of features used to build the model based on data included in the testing set. We suggest that for currently available dataset sizes the expected classifier performance is between 0.55 and 0.7 (AUC) and higher accuracies reported in literature are likely a result of overfitting.

HiCEnterprise: identifying long range chromosomal contacts in Hi-C data.

MOTIVATION: Computational analysis of chromosomal contact data is currently gaining popularity with the rapid advance in experimental techniques providing access to a growing body of data. An important problem in this area is the identification of long range contacts between distinct chromatin regions. Such loops were shown to exist at different scales, either mediating relatively short range interactions between enhancers and promoters or providing interactions between much larger, distant chromosome domains. A proper statistical analysis as well as availability to a wide research community are crucial in a tool for this task. RESULTS: We present HiCEnterprise, a first freely available software tool for identification of long range chromatin contacts not only between small regions, but also between chromosomal domains. It implements four different statistical tests for identification of significant contacts for user defined regions or domains as well as necessary functions for input, output and visualization of chromosome contacts. AVAILABILITY: The software and the corresponding documentation are available at: github.com/regulomics/HiCEnterprise. SUPPLEMENTARY INFORMATION: Supplemental data are available in the online version of the article and at the website regulomics.mimuw.edu.pl/wp/hicenterprise.

The role of epigenetic modifications, long-range contacts, enhancers and topologically associating domains in the regulation of glioma grade-specific genes.

Genome-wide studies have uncovered specific genetic alterations, transcriptomic patterns and epigenetic profiles associated with different glioma types. We have recently created a unique atlas encompassing genome-wide profiles of open chromatin, histone H3K27ac and H3Kme3 modifications, DNA methylation and transcriptomes of 33 glioma samples of different grades. Here, we intersected genome-wide atlas data with topologically associating domains (TADs) and demonstrated that the chromatin organization and epigenetic landscape of enhancers have a strong impact on genes differentially expressed in WHO low grade versus high grade gliomas. We identified TADs enriched in glioma grade-specific genes and/or epigenetic marks. We found the set of transcription factors, including REST, E2F1 and NFKB1, that are most likely to regulate gene expression in multiple TADs, containing specific glioma-related genes. Moreover, many genes associated with the cell-matrix adhesion Gene Ontology group, in particular 14 PROTOCADHERINs, were found to be regulated by long-range contacts with enhancers. Presented results demonstrate the existence of epigenetic differences associated with chromatin organization driving differential gene expression in gliomas of different malignancy.

Mapping chromatin accessibility and active regulatory elements reveals pathological mechanisms in human gliomas.

Chromatin structure and accessibility, and combinatorial binding of transcription factors to regulatory elements in genomic DNA control transcription. Genetic variations in genes encoding histones, epigenetics-related enzymes or modifiers affect chromatin structure/dynamics and result in alterations in gene expression contributing to cancer development or progression. Gliomas are brain tumors frequently associated with epigenetics-related gene deregulation. We perform whole-genome mapping of chromatin accessibility, histone modifications, DNA methylation patterns and transcriptome analysis simultaneously in multiple tumor samples to unravel epigenetic dysfunctions driving gliomagenesis. Based on the results of the integrative analysis of the acquired profiles, we create an atlas of active enhancers and promoters in benign and malignant gliomas. We explore these elements and intersect with Hi-C data to uncover molecular mechanisms instructing gene expression in gliomas.