Genome-wide localization of protein-DNA binding and histone modification by a Bayesian change-point method with ChIP-seq data

Multiplex-GAM: genome-wide identification of chromatin contacts yields insights overlooked by Hi-C

Technology for measuring 3D genome topology is increasingly important for studying gene regulation, for genome assembly and for mapping of genome rearrangements. Hi-C and other ligation-based methods have become routine but have specific biases. Here, we develop multiplex-GAM, a faster and more affordable version of genome architecture mapping (GAM), a ligation-free technique that maps chromatin contacts genome-wide. We perform a detailed comparison of multiplex-GAM and Hi-C using mouse embryonic stem cells. When examining the strongest contacts detected by either method, we find that only one-third of these are shared. The strongest contacts specifically found in GAM often involve 'active' regions, including many transcribed genes and super-enhancers, whereas in Hi-C they more often contain 'inactive' regions. Our work shows that active genomic regions are involved in extensive complex contacts that are currently underestimated in ligation-based approaches, and highlights the need for orthogonal advances in genome-wide contact mapping technologies.

Dynamic antagonism between key repressive pathways maintains the placental epigenome

DNA and Histone 3 Lysine 27 methylation typically function as repressive modifications and operate within distinct genomic compartments. In mammals, the majority of the genome is kept in a DNA methylated state, whereas the Polycomb repressive complexes regulate the unmethylated CpG-rich promoters of developmental genes. In contrast to this general framework, the extra-embryonic lineages display non-canonical, globally intermediate DNA methylation levels, including disruption of local Polycomb domains. Here, to better understand this unusual landscape's molecular properties, we genetically and chemically perturbed major epigenetic pathways in mouse trophoblast stem cells. We find that the extra-embryonic epigenome reflects ongoing and dynamic de novo methyltransferase recruitment, which is continuously antagonized by Polycomb to maintain intermediate, locally disordered methylation. Despite its disorganized molecular appearance, our data point to a highly controlled equilibrium between counteracting repressors within extra-embryonic cells, one that can seemingly persist indefinitely without bistable features typically seen for embryonic forms of epigenetic regulation.

Identification of Enriched Regions in ChIP-Seq Data via a Linear-Time Multi-Level Thresholding Algorithm

Chromatin immunoprecipitation (ChIP-Seq) has emerged as a superior alternative to microarray technology as it provides higher resolution, less noise, greater coverage and wider dynamic range. While ChIP-Seq enables probing of DNA-protein interaction over the entire genome, it requires the use of sophisticated tools to recognize hidden patterns and extract meaningful data. Over the years, various attempts have resulted in several algorithms making use of different heuristics to accurately determine individual peaks corresponding to unique DNA-protein. However, finding all the significant peaks with high accuracy in a reasonable time is still a challenge. In this work, we propose the use of Multi-level thresholding algorithm, which we call LinMLTBS, used to identify the enriched regions on ChIP-Seq data. Although various suboptimal heuristics have been proposed for multi-level thresholding, we emphasize on the use of an algorithm capable of obtaining an optimal solution, while maintaining linear-time complexity. Testing various algorithm on various ENCODE project datasets shows that our approach attains higher accuracy relative to previously proposed peak finders while retaining a reasonable processing speed.

Transcription factor protein interactomes reveal genetic determinants in heart disease

Congenital heart disease (CHD) is present in 1% of live births, yet identification of causal mutations remains challenging. We hypothesized that genetic determinants for CHDs may lie in the protein interactomes of transcription factors whose mutations cause CHDs. Defining the interactomes of two transcription factors haplo-insufficient in CHD, GATA4 and TBX5, within human cardiac progenitors, and integrating the results with nearly 9,000 exomes from proband-parent trios revealed an enrichment of de novo missense variants associated with CHD within the interactomes. Scoring variants of interactome members based on residue, gene, and proband features identified likely CHD-causing genes, including the epigenetic reader GLYR1. GLYR1 and GATA4 widely co-occupied and co-activated cardiac developmental genes, and the identified GLYR1 missense variant disrupted interaction with GATA4, impairing in vitro and in vivo function in mice. This integrative proteomic and genetic approach provides a framework for prioritizing and interrogating genetic variants in heart disease.

Brahma safeguards canalization of cardiac mesoderm differentiation

Differentiation proceeds along a continuum of increasingly fate-restricted intermediates, referred to as canalization^1,2. Canalization is essential for stabilizing cell fate, but the mechanisms that underlie robust canalization are unclear. Here we show that the BRG1/BRM-associated factor (BAF) chromatin-remodelling complex ATPase gene Brm safeguards cell identity during directed cardiogenesis of mouse embryonic stem cells. Despite the establishment of a well-differentiated precardiac mesoderm, Brm^-/- cells predominantly became neural precursors, violating germ layer assignment. Trajectory inference showed a sudden acquisition of a non-mesodermal identity in Brm^-/- cells. Mechanistically, the loss of Brm prevented de novo accessibility of primed cardiac enhancers while increasing the expression of neurogenic factor POU3F1, preventing the binding of the neural suppressor REST and shifting the composition of BRG1 complexes. The identity switch caused by the Brm mutation was overcome by increasing BMP4 levels during mesoderm induction. Mathematical modelling supports these observations and demonstrates that Brm deletion affects cell fate trajectory by modifying saddle-node bifurcations². In the mouse embryo, Brm deletion exacerbated mesoderm-deleted Brg1-mutant phenotypes, severely compromising cardiogenesis, and reveals an in vivo role for Brm. Our results show that Brm is a compensable safeguard of the fidelity of mesoderm chromatin states, and support a model in which developmental canalization is not a rigid irreversible path, but a highly plastic trajectory.

VSS: Variance-stabilized signals for sequencing-based genomic signals

MOTIVATION

A sequencing-based genomic assay such as ChIP-seq outputs a real-valued signal for each position in the genome that measures the strength of activity at that position. Most genomic signals lack the property of variance stabilization. That is, a difference between 0 and 100 reads usually has a very different statistical importance from a difference between 1,000 and 1,100 reads. A statistical model such as a negative binomial distribution can account for this pattern, but learning these models is computationally challenging. Therefore, many applications - including imputation and segmentation and genome annotation (SAGA) - instead use Gaussian models and use a transformation such as log or inverse hyperbolic sine (asinh) to stabilize variance.

RESULTS

We show here that existing transformations do not fully stabilize variance in genomic data sets. To solve this issue, we propose VSS, a method that produces variance-stabilized signals for sequencing-based genomic signals. VSS learns the empirical relationship between the mean and variance of a given signal data set and produces transformed signals that normalize for this dependence. We show that VSS successfully stabilizes variance and that doing so improves downstream applications such as SAGA. VSS will eliminate the need for downstream methods to implement complex mean-variance relationship models, and will enable genomic signals to be easily understood by eye.

AVAILABILITY

https://github.com/faezeh-bayat/VSS.

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

A flexible ChIP-sequencing simulation toolkit

Background

A major challenge in evaluating quantitative ChIP-seq analyses, such as peak calling and differential binding, is a lack of reliable ground truth data. Accurate simulation of ChIP-seq data can mitigate this challenge, but existing frameworks are either too cumbersome to apply genome-wide or unable to model a number of important experimental conditions in ChIP-seq.

Results

We present ChIPs, a toolkit for rapidly simulating ChIP-seq data using statistical models of key experimental steps. We demonstrate how ChIPs can be used for a range of applications, including benchmarking analysis tools and evaluating the impact of various experimental parameters. ChIPs is implemented as a standalone command-line program written in C++ and is available from https://github.com/gymreklab/chips.

Conclusions

ChIPs is an efficient ChIP-seq simulation framework that generates realistic datasets over a flexible range of experimental conditions. It can serve as an important component in various ChIP-seq analyses where ground truth data are needed.

Supplementary Information

The online version contains supplementary material available at 10.1186/s12859-021-04097-5.

Deciphering hierarchical organization of topologically associated domains through change-point testing

Background

The nucleus of eukaryotic cells spatially packages chromosomes into a hierarchical and distinct segregation that plays critical roles in maintaining transcription regulation. High-throughput methods of chromosome conformation capture, such as Hi-C, have revealed topologically associating domains (TADs) that are defined by biased chromatin interactions within them.

Results

We introduce a novel method, HiCKey, to decipher hierarchical TAD structures in Hi-C data and compare them across samples. We first derive a generalized likelihood-ratio (GLR) test for detecting change-points in an interaction matrix that follows a negative binomial distribution or general mixture distribution. We then employ several optimal search strategies to decipher hierarchical TADs with p values calculated by the GLR test. Large-scale validations of simulation data show that HiCKey has good precision in recalling known TADs and is robust against random collisions of chromatin interactions. By applying HiCKey to Hi-C data of seven human cell lines, we identified multiple layers of TAD organization among them, but the vast majority had no more than four layers. In particular, we found that TAD boundaries are significantly enriched in active chromosomal regions compared to repressed regions.

Conclusions

HiCKey is optimized for processing large matrices constructed from high-resolution Hi-C experiments. The method and theoretical result of the GLR test provide a general framework for significance testing of similar experimental chromatin interaction data that may not fully follow negative binomial distributions but rather more general mixture distributions.

Supplementary Information

The online version contains supplementary material available at 10.1186/s12859-021-04113-8.

RECAP reveals the true statistical significance of ChIP-seq peak calls

Motivation

Chromatin Immunopreciptation (ChIP)-seq is used extensively to identify sites of transcription factor binding or regions of epigenetic modifications to the genome. A key step in ChIP-seq analysis is peak calling, where genomic regions enriched for ChIP versus control reads are identified. Many programs have been designed to solve this task, but nearly all fall into the statistical trap of using the data twice—once to determine candidate enriched regions, and again to assess enrichment by classical statistical hypothesis testing. This double use of the data invalidates the statistical significance assigned to enriched regions, thus the true significance or reliability of peak calls remains unknown.

Results

Using simulated and real ChIP-seq data, we show that three well-known peak callers, MACS, SICER and diffReps, output biased P-values and false discovery rate estimates that can be many orders of magnitude too optimistic. We propose a wrapper algorithm, RECAP, that uses resampling of ChIP-seq and control data to estimate a monotone transform correcting for biases built into peak calling algorithms. When applied to null hypothesis data, where there is no enrichment between ChIP-seq and control, P-values recalibrated by RECAP are approximately uniformly distributed. On data where there is genuine enrichment, RECAP P-values give a better estimate of the true statistical significance of candidate peaks and better false discovery rate estimates, which correlate better with empirical reproducibility. RECAP is a powerful new tool for assessing the true statistical significance of ChIP-seq peak calls.

Availability and implementation

The RECAP software is available through www.perkinslab.ca or on github at https://github.com/theodorejperkins/RECAP.

Supplementary information

Supplementary data are available at Bioinformatics online.

Seq-ing answers: Current data integration approaches to uncover mechanisms of transcriptional regulation

Advancements in the field of next generation sequencing lead to the generation of ever-more data, with the challenge often being how to combine and reconcile results from different OMICs studies such as genome, epigenome and transcriptome. Here we provide an overview of the standard processing pipelines for ChIP-seq and RNA-seq as well as common downstream analyses. We describe popular multi-omics data integration approaches used to identify target genes and co-factors, and we discuss how machine learning techniques may predict transcriptional regulators and gene expression.

Improved detection of epigenomic marks with mixed-effects hidden Markov models

Chromatin immunoprecipitation followed by next-generation sequencing (ChIP-seq) is a technique to detect genomic regions containing protein-DNA interaction, such as transcription factor binding sites or regions containing histone modifications. One goal of the analysis of ChIP-seq experiments is to identify genomic loci enriched for sequencing reads pertaining to DNA bound to the factor of interest. The accurate identification of such regions aids in the understanding of epigenomic marks and gene regulatory mechanisms. Given the reduction of massively parallel sequencing costs, methods to detect consensus regions of enrichment across multiple samples are of interest. Here, we present a statistical model to detect broad consensus regions of enrichment from ChIP-seq technical or biological replicates through a class of zero-inflated mixed-effects hidden Markov models. We show that the proposed model outperforms existing methods for consensus peak calling in common epigenomic marks by accounting for the excess zeros and sample-specific biases. We apply our method to data from the Encyclopedia of DNA Elements and Roadmap Epigenomics projects and also from an extensive simulation study.

A map of direct TF-DNA interactions in the human genome

Chromatin immunoprecipitation followed by sequencing (ChIP-seq) is the most popular assay to identify genomic regions, called ChIP-seq peaks, that are bound in vivo by transcription factors (TFs). These regions are derived from direct TF-DNA interactions, indirect binding of the TF to the DNA (through a co-binding partner), nonspecific binding to the DNA, and noise/bias/artifacts. Delineating the bona fide direct TF-DNA interactions within the ChIP-seq peaks remains challenging. We developed a dedicated software, ChIP-eat, that combines computational TF binding models and ChIP-seq peaks to automatically predict direct TF-DNA interactions. Our work culminated with predicted interactions covering >4% of the human genome, obtained by uniformly processing 1983 ChIP-seq peak data sets from the ReMap database for 232 unique TFs. The predictions were a posteriori assessed using protein binding microarray and ChIP-exo data, and were predominantly found in high quality ChIP-seq peaks. The set of predicted direct TF-DNA interactions suggested that high-occupancy target regions are likely not derived from direct binding of the TFs to the DNA. Our predictions derived co-binding TFs supported by protein-protein interaction data and defined cis-regulatory modules enriched for disease- and trait-associated SNPs. We provide this collection of direct TF-DNA interactions and cis-regulatory modules through the UniBind web-interface (http://unibind.uio.no).

Cell-autonomous and redundant roles of Hey1 and HeyL in muscle stem cells: HeyL requires Hes1 to bind diverse DNA sites

The undifferentiated state of muscle stem (satellite) cells (MuSCs) is maintained by the canonical Notch pathway. Although three bHLH transcriptional factors, Hey1, HeyL, and Hes1, are considered to be potential effectors of the Notch pathway exerting anti-myogenic effects, neither HeyL nor Hes1 inhibits myogenic differentiation of myogenic cell lines. Furthermore, whether these factors work redundantly or cooperatively is unknown. Here, we showed cell-autonomous functions of Hey1 and HeyL in MuSCs using conditional and genetic null mice. Analysis of cultured MuSCs revealed anti-myogenic activity of both HeyL and Hes1. We found that HeyL forms heterodimeric complexes with Hes1 in living cells. Moreover, our ChIP-Seq experiments demonstrated that, compared with HeyL alone, HeyL-Hes1 heterodimer bound with high affinity to specific sites in the chromatin including the cis-element of Hey1. Finally, the analyses of myogenin promoter activity showed that HeyL and Hes1 acted synergistically to suppress myogenic differentiation. Collectively, those results suggest that HeyL and Hey1 function redundantly in MuSCs, and that HeyL requires Hes1 for effective DNA binding and biological activity.

Clustering-local-unique-enriched-signals (CLUES) promotes identification of novel regulators of ES cell self-renewal and pluripotency

BACKGROUND

Key regulators of developmental processes can be prioritized through integrated analysis of ChIP-Seq data of master transcriptional factors (TFs) such as Nanog and Oct4, active histone modifications (HMs) such as H3K4me3 and H3K27ac, and repressive HMs such as H3K27me3. Recent studies show that broad enrichment signals such as super-enhancers and broad H3K4me3 enrichment signals play more dominant roles than short enrichment signals of the master TFs and H3K4me3 in epigenetic regulatory mechanism. Besides the broad enrichment signals, up to ten thousands of short enrichment signals of these TFs and HMs exist in genome. Prioritization of these broad enrichment signals from ChIP-Seq data is a prerequisite for such integrated analysis.

RESULTS

Here, we present a method named Clustering-Local-Unique-Enriched-Signals (CLUES), which uses an adaptive-size-windows strategy to identify enriched regions (ERs) and cluster them into broad enrichment signals. Tested on 62 ENCODE ChIP-Seq datasets of Ctcf and Nrsf, CLUES performs equally well as MACS2 regarding prioritization of ERs with the TF's motif. Tested on 165 ENCODE ChIP-Seq datasets of H3K4me3, H3K27me3, and H3K36me3, CLUES performs better than existing algorithms on prioritizing broad enrichment signals implicating cell functions influenced by epigenetic regulatory mechanism in cells. Most importantly, CLUES helps to confirm several novel regulators of mouse ES cell self-renewal and pluripotency through integrated analysis of prioritized broad enrichment signals of H3K4me3, H3K27me3, Nanog and Oct4 with the support of a CRISPR/Cas9 negative selection genetic screen.

CONCLUSIONS

CLUES holds promise for prioritizing broad enrichment signals from ChIP-Seq data. The download site for CLUES is https://github.com/Wuchao1984/CLUESv1.

ChromTime: modeling spatio-temporal dynamics of chromatin marks

To model spatial changes of chromatin mark peaks over time we develop and apply ChromTime, a computational method that predicts peaks to be either expanding, contracting, or holding steady between time points. Predicted expanding and contracting peaks can mark regulatory regions associated with transcription factor binding and gene expression changes. Spatial dynamics of peaks provide information about gene expression changes beyond localized signal density changes. ChromTime detects asymmetric expansions and contractions, which for some marks associate with the direction of transcription. ChromTime facilitates the analysis of time course chromatin data in a range of biological systems.

Methionine metabolism influences genomic architecture and gene expression through H3K4me3 peak width

Nutrition and metabolism are known to influence chromatin biology and epigenetics through post-translational modifications, yet how this interaction influences genomic architecture and connects to gene expression is unknown. Here we consider, as a model, the metabolically-driven dynamics of H3K4me3, a histone methylation mark that is known to encode information about active transcription, cell identity, and tumor suppression. We analyze the genome-wide changes in H3K4me3 and gene expression in response to alterations in methionine availability in both normal mouse physiology and human cancer cells. Surprisingly, we find that the location of H3K4me3 peaks is largely preserved under methionine restriction, while the response of H3K4me3 peak width encodes almost all aspects of H3K4me3 biology including changes in expression levels, and the presence of cell identity and cancer-associated genes. These findings may reveal general principles for how nutrient availability modulates specific aspects of chromatin dynamics to mediate biological function.

Genome-wide analysis of PDX1 target genes in human pancreatic progenitors

Objective

Homozygous loss-of-function mutations in the gene coding for the homeobox transcription factor (TF) PDX1 leads to pancreatic agenesis, whereas heterozygous mutations can cause Maturity-Onset Diabetes of the Young 4 (MODY4). Although the function of Pdx1 is well studied in pre-clinical models during insulin-producing β-cell development and homeostasis, it remains elusive how this TF controls human pancreas development by regulating a downstream transcriptional program. Also, comparative studies of PDX1 binding patterns in pancreatic progenitors and adult β-cells have not been conducted so far. Furthermore, many studies reported the association between single nucleotide polymorphisms (SNPs) and T2DM, and it has been shown that islet enhancers are enriched in T2DM-associated SNPs. Whether regions, harboring T2DM-associated SNPs are PDX1 bound and active at the pancreatic progenitor stage has not been reported so far.

Methods

In this study, we have generated a novel induced pluripotent stem cell (iPSC) line that efficiently differentiates into human pancreatic progenitors (PPs). Furthermore, PDX1 and H3K27ac chromatin immunoprecipitation sequencing (ChIP-seq) was used to identify PDX1 transcriptional targets and active enhancer and promoter regions. To address potential differences in the function of PDX1 during development and adulthood, we compared PDX1 binding profiles from PPs and adult islets. Moreover, combining ChIP-seq and GWAS meta-analysis data we identified T2DM-associated SNPs in PDX1 binding sites and active chromatin regions.

Results

ChIP-seq for PDX1 revealed a total of 8088 PDX1-bound regions that map to 5664 genes in iPSC-derived PPs. The PDX1 target regions include important pancreatic TFs, such as PDX1 itself, RFX6, HNF1B, and MEIS1, which were activated during the differentiation process as revealed by the active chromatin mark H3K27ac and mRNA expression profiling, suggesting that auto-regulatory feedback regulation maintains PDX1 expression and initiates a pancreatic TF program. Remarkably, we identified several PDX1 target genes that have not been reported in the literature in human so far, including RFX3, required for ciliogenesis and endocrine differentiation in mouse, and the ligand of the Notch receptor DLL1, which is important for endocrine induction and tip-trunk patterning. The comparison of PDX1 profiles from PPs and adult human islets identified sets of stage-specific target genes, associated with early pancreas development and adult β-cell function, respectively. Furthermore, we found an enrichment of T2DM-associated SNPs in active chromatin regions from iPSC-derived PPs. Two of these SNPs fall into PDX1 occupied sites that are located in the intronic regions of TCF7L2 and HNF1B. Both of these genes are key transcriptional regulators of endocrine induction and mutations in cis-regulatory regions predispose to diabetes.

Conclusions

Our data provide stage-specific target genes of PDX1 during in vitro differentiation of stem cells into pancreatic progenitors that could be useful to identify pathways and molecular targets that predispose for diabetes. In addition, we show that T2DM-associated SNPs are enriched in active chromatin regions at the pancreatic progenitor stage, suggesting that the susceptibility to T2DM might originate from imperfect execution of a β-cell developmental program.

•

PDX1 ChIP-seq analysis revealed 5664 target genes in human pancreatic progenitors, including unreported target genes.

•

Comparison of PDX1 profiles from PPs and adult human islets identified stage-specific PDX1 target gene sets.

•

T2DM-associated SNPs are enriched in active chromatin regions from iPSC-derived PPs.

•

Three SNPs fall into PDX1 occupied sites, located in intronic regions of the developmental regulatory TFs TCF7L2 and HNF1B.

Recent advances in ChIP-seq analysis: from quality management to whole-genome annotation

Chromatin immunoprecipitation followed by sequencing (ChIP-seq) analysis can detect protein/DNA-binding and histone-modification sites across an entire genome. Recent advances in sequencing technologies and analyses enable us to compare hundreds of samples simultaneously; such large-scale analysis has potential to reveal the high-dimensional interrelationship level for regulatory elements and annotate novel functional genomic regions de novo. Because many experimental considerations are relevant to the choice of a method in a ChIP-seq analysis, the overall design and quality management of the experiment are of critical importance. This review offers guiding principles of computation and sample preparation for ChIP-seq analyses, highlighting the validity and limitations of the state-of-the-art procedures at each step. We also discuss the latest challenges of single-cell analysis that will encourage a new era in this field.

RNA polymerase II primes Polycomb-repressed developmental genes throughout terminal neuronal differentiation

Polycomb repression in mouse embryonic stem cells (ESCs) is tightly associated with promoter co‐occupancy of RNA polymerase II (RNAPII) which is thought to prime genes for activation during early development. However, it is unknown whether RNAPII poising is a general feature of Polycomb repression, or is lost during differentiation. Here, we map the genome‐wide occupancy of RNAPII and Polycomb from pluripotent ESCs to non‐dividing functional dopaminergic neurons. We find that poised RNAPII complexes are ubiquitously present at Polycomb‐repressed genes at all stages of neuronal differentiation. We observe both loss and acquisition of RNAPII and Polycomb at specific groups of genes reflecting their silencing or activation. Strikingly, RNAPII remains poised at transcription factor genes which are silenced in neurons through Polycomb repression, and have major roles in specifying other, non‐neuronal lineages. We conclude that RNAPII poising is intrinsically associated with Polycomb repression throughout differentiation. Our work suggests that the tight interplay between RNAPII poising and Polycomb repression not only instructs promoter state transitions, but also may enable promoter plasticity in differentiated cells.

The Short Isoform of BRD4 Promotes HIV-1 Latency by Engaging Repressive SWI/SNF Chromatin-Remodeling Complexes

BET proteins commonly activate cellular gene expression, yet inhibiting their recruitment paradoxically reactivates latent HIV-1 transcription. Here we identify the short isoform of BET family member BRD4 (BRD4S) as a corepressor of HIV-1 transcription. We found that BRD4S was enriched in chromatin fractions of latently infected T cells, and it was more rapidly displaced from chromatin upon BET inhibition than the long isoform. BET inhibition induced marked nucleosome remodeling at the latent HIV-1 promoter, which was dependent on the activity of BRG1-associated factors (BAF), an SWI/SNF chromatin-remodeling complex with known repressive functions in HIV-1 transcription. BRD4S directly bound BRG1, a catalytic subunit of BAF, via its bromodomain and extraterminal (ET) domain, and this isoform was necessary for BRG1 recruitment to latent HIV-1 chromatin. Using chromatin immunoprecipitation sequencing (ChIP-seq) combined with assay for transposase-accessible chromatin coupled to high-throughput sequencing (ATAC-seq) data, we found that the latent HIV-1 promoter phenotypically resembles endogenous long terminal repeat (LTR) sequences, pointing to a select role of BRD4S-BRG1 complexes in genomic silencing of invasive retroelements.

Mixture model with multiple allocations for clustering spatially correlated observations in the analysis of ChIP-Seq data

Model-based clustering is a technique widely used to group a collection of units into mutually exclusive groups. There are, however, situations in which an observation could in principle belong to more than one cluster. In the context of next-generation sequencing (NGS) experiments, for example, the signal observed in the data might be produced by two (or more) different biological processes operating together and a gene could participate in both (or all) of them. We propose a novel approach to cluster NGS discrete data, coming from a ChIP-Seq experiment, with a mixture model, allowing each unit to belong potentially to more than one group: these multiple allocation clusters can be flexibly defined via a function combining the features of the original groups without introducing new parameters. The formulation naturally gives rise to a 'zero-inflation group' in which values close to zero can be allocated, acting as a correction for the abundance of zeros that manifest in this type of data. We take into account the spatial dependency between observations, which is described through a latent conditional autoregressive process that can reflect different dependency patterns. We assess the performance of our model within a simulation environment and then we apply it to ChIP-seq real data.

Flipping between Polycomb repressed and active transcriptional states introduces noise in gene expression

Polycomb repressive complexes (PRCs) are important histone modifiers, which silence gene expression; yet, there exists a subset of PRC-bound genes actively transcribed by RNA polymerase II (RNAPII). It is likely that the role of Polycomb repressive complex is to dampen expression of these PRC-active genes. However, it is unclear how this flipping between chromatin states alters the kinetics of transcription. Here, we integrate histone modifications and RNAPII states derived from bulk ChIP-seq data with single-cell RNA-sequencing data. We find that Polycomb repressive complex-active genes have greater cell-to-cell variation in expression than active genes, and these results are validated by knockout experiments. We also show that PRC-active genes are clustered on chromosomes in both two and three dimensions, and interactions with active enhancers promote a stabilization of gene expression noise. These findings provide new insights into how chromatin regulation modulates stochastic gene expression and transcriptional bursting, with implications for regulation of pluripotency and development.Polycomb repressive complexes modify histones but it is unclear how changes in chromatin states alter kinetics of transcription. Here, the authors use single-cell RNAseq and ChIPseq to find that actively transcribed genes with Polycomb marks have greater cell-to-cell variation in expression.

Targeted Degradation of CTCF Decouples Local Insulation of Chromosome Domains from Genomic Compartmentalization

The molecular mechanisms underlying folding of mammalian chromosomes remain poorly understood. The transcription factor CTCF is a candidate regulator of chromosomal structure. Using the auxin-inducible degron system in mouse embryonic stem cells, we show that CTCF is absolutely and dose-dependently required for looping between CTCF target sites and insulation of topologically associating domains (TADs). Restoring CTCF reinstates proper architecture on altered chromosomes, indicating a powerful instructive function for CTCF in chromatin folding. CTCF remains essential for TAD organization in non-dividing cells. Surprisingly, active and inactive genome compartments remain properly segregated upon CTCF depletion, revealing that compartmentalization of mammalian chromosomes emerges independently of proper insulation of TADs. Furthermore, our data support that CTCF mediates transcriptional insulator function through enhancer blocking but not as a direct barrier to heterochromatin spreading. Beyond defining the functions of CTCF in chromosome folding, these results provide new fundamental insights into the rules governing mammalian genome organization.

Features that define the best ChIP-seq peak calling algorithms

Chromatin immunoprecipitation followed by sequencing (ChIP-seq) is an important tool for studying gene regulatory proteins, such as transcription factors and histones. Peak calling is one of the first steps in the analysis of these data. Peak calling consists of two sub-problems: identifying candidate peaks and testing candidate peaks for statistical significance. We surveyed 30 methods and identified 12 features of the two sub-problems that distinguish methods from each other. We picked six methods GEM, MACS2, MUSIC, BCP, Threshold-based method (TM) and ZINBA] that span this feature space and used a combination of 300 simulated ChIP-seq data sets, 3 real data sets and mathematical analyses to identify features of methods that allow some to perform better than the others. We prove that methods that explicitly combine the signals from ChIP and input samples are less powerful than methods that do not. Methods that use windows of different sizes are more powerful than the ones that do not. For statistical testing of candidate peaks, methods that use a Poisson test to rank their candidate peaks are more powerful than those that use a Binomial test. BCP and MACS2 have the best operating characteristics on simulated transcription factor binding data. GEM has the highest fraction of the top 500 peaks containing the binding motif of the immunoprecipitated factor, with 50% of its peaks within 10 base pairs of a motif. BCP and MUSIC perform best on histone data. These findings provide guidance and rationale for selecting the best peak caller for a given application.

Reusable, extensible, and modifiable R scripts and Kepler workflows for comprehensive single set ChIP-seq analysis

Background

There has been an enormous expansion of use of chromatin immunoprecipitation followed by sequencing (ChIP-seq) technologies. Analysis of large-scale ChIP-seq datasets involves a complex series of steps and production of several specialized graphical outputs. A number of systems have emphasized custom development of ChIP-seq pipelines. These systems are primarily based on custom programming of a single, complex pipeline or supply libraries of modules and do not produce the full range of outputs commonly produced for ChIP-seq datasets. It is desirable to have more comprehensive pipelines, in particular ones addressing common metadata tasks, such as pathway analysis, and pipelines producing standard complex graphical outputs. It is advantageous if these are highly modular systems, available as both turnkey pipelines and individual modules, that are easily comprehensible, modifiable and extensible to allow rapid alteration in response to new analysis developments in this growing area. Furthermore, it is advantageous if these pipelines allow data provenance tracking.

Results

We present a set of 20 ChIP-seq analysis software modules implemented in the Kepler workflow system; most (18/20) were also implemented as standalone, fully functional R scripts. The set consists of four full turnkey pipelines and 16 component modules. The turnkey pipelines in Kepler allow data provenance tracking. Implementation emphasized use of common R packages and widely-used external tools (e.g., MACS for peak finding), along with custom programming. This software presents comprehensive solutions and easily repurposed code blocks for ChIP-seq analysis and pipeline creation. Tasks include mapping raw reads, peakfinding via MACS, summary statistics, peak location statistics, summary plots centered on the transcription start site (TSS), gene ontology, pathway analysis, and de novo motif finding, among others.

Conclusions

These pipelines range from those performing a single task to those performing full analyses of ChIP-seq data. The pipelines are supplied as both Kepler workflows, which allow data provenance tracking, and, in the majority of cases, as standalone R scripts. These pipelines are designed for ease of modification and repurposing.

ChIP-BIT: Bayesian inference of target genes using a novel joint probabilistic model of ChIP-seq profiles

Chromatin immunoprecipitation with massively parallel DNA sequencing (ChIP-seq) has greatly improved the reliability with which transcription factor binding sites (TFBSs) can be identified from genome-wide profiling studies. Many computational tools are developed to detect binding events or peaks, however the robust detection of weak binding events remains a challenge for current peak calling tools. We have developed a novel Bayesian approach (ChIP-BIT) to reliably detect TFBSs and their target genes by jointly modeling binding signal intensities and binding locations of TFBSs. Specifically, a Gaussian mixture model is used to capture both binding and background signals in sample data. As a unique feature of ChIP-BIT, background signals are modeled by a local Gaussian distribution that is accurately estimated from the input data. Extensive simulation studies showed a significantly improved performance of ChIP-BIT in target gene prediction, particularly for detecting weak binding signals at gene promoter regions. We applied ChIP-BIT to find target genes from NOTCH3 and PBX1 ChIP-seq data acquired from MCF-7 breast cancer cells. TF knockdown experiments have initially validated about 30% of co-regulated target genes identified by ChIP-BIT as being differentially expressed in MCF-7 cells. Functional analysis on these genes further revealed the existence of crosstalk between Notch and Wnt signaling pathways.

Detecting broad domains and narrow peaks in ChIP-seq data with hiddenDomains

Background

Correctly identifying genomic regions enriched with histone modifications and transcription factors is key to understanding their regulatory and developmental roles. Conceptually, these regions are divided into two categories, narrow peaks and broad domains, and different algorithms are used to identify each one. Datasets that span these two categories are often analyzed with a single program for peak calling combined with an ad hoc method for domains.

Results

We developed hiddenDomains, which identifies both peaks and domains, and compare it to the leading algorithms using H3K27me3, H3K36me3, GABP, ESR1 and FOXA ChIP-seq datasets. The output from the programs was compared to qPCR-validated enriched and depleted sites, predicted transcription factor binding sites, and highly-transcribed gene bodies. With every method, hiddenDomains, performed as well as, if not better than algorithms dedicated to a specific type of analysis.

Conclusions

hiddenDomains performs as well as the best domain and peak calling algorithms, making it ideal for analyzing ChIP-seq datasets, especially those that contain a mixture of peaks and domains.

Electronic supplementary material

The online version of this article (doi:10.1186/s12859-016-0991-z) contains supplementary material, which is available to authorized users.

Methylation of RNA polymerase II non-consensus Lysine residues marks early transcription in mammalian cells

Dynamic post-translational modification of RNA polymerase II (RNAPII) coordinates the co-transcriptional recruitment of enzymatic complexes that regulate chromatin states and processing of nascent RNA. Extensive phosphorylation of serine residues at the largest RNAPII subunit occurs at its structurally-disordered C-terminal domain (CTD), which is composed of multiple heptapeptide repeats with consensus sequence Y1-S2-P3-T4-S5-P6-S7. Serine-5 and Serine-7 phosphorylation mark transcription initiation, whereas Serine-2 phosphorylation coincides with productive elongation. In vertebrates, the CTD has eight non-canonical substitutions of Serine-7 into Lysine-7, which can be acetylated (K7ac). Here, we describe mono- and di-methylation of CTD Lysine-7 residues (K7me1 and K7me2). K7me1 and K7me2 are observed during the earliest transcription stages and precede or accompany Serine-5 and Serine-7 phosphorylation. In contrast, K7ac is associated with RNAPII elongation, Serine-2 phosphorylation and mRNA expression. We identify an unexpected balance between RNAPII K7 methylation and acetylation at gene promoters, which fine-tunes gene expression levels.

MUSIC: identification of enriched regions in ChIP-Seq experiments using a mappability-corrected multiscale signal processing framework

We present MUSIC, a signal processing approach for identification of enriched regions in ChIP-Seq data, available at music.gersteinlab.org. MUSIC first filters the ChIP-Seq read-depth signal for systematic noise from non-uniform mappability, which fragments enriched regions. Then it performs a multiscale decomposition, using median filtering, identifying enriched regions at multiple length scales. This is useful given the wide range of scales probed in ChIP-Seq assays. MUSIC performs favorably in terms of accuracy and reproducibility compared with other methods. In particular, analysis of RNA polymerase II data reveals a clear distinction between the stalled and elongating forms of the polymerase.

Saturation analysis of ChIP-seq data for reproducible identification of binding peaks

Chromatin immunoprecipitation coupled with next-generation sequencing (ChIP-seq) is a powerful technology to identify the genome-wide locations of transcription factors and other DNA binding proteins. Computational ChIP-seq peak calling infers the location of protein–DNA interactions based on various measures of enrichment of sequence reads. In this work, we introduce an algorithm, Q, that uses an assessment of the quadratic enrichment of reads to center candidate peaks followed by statistical analysis of saturation of candidate peaks by 5′ ends of reads. We show that our method not only is substantially faster than several competing methods but also demonstrates statistically significant advantages with respect to reproducibility of results and in its ability to identify peaks with reproducible binding site motifs. We show that Q has superior performance in the delineation of double RNAPII and H3K4me3 peaks surrounding transcription start sites related to a better ability to resolve individual peaks. The method is implemented in C+l+ and is freely available under an open source license.

Genome Wide Mapping of Foxo1 Binding-sites in Murine T Lymphocytes

The Forkhead box O (Foxo) family of transcription factors has a critical role in controlling the development, differentiation, and function of T cells. However, the direct target genes of Foxo transcription factors in T cells have not been well characterized. In this study, we focused on mapping the genome wide Foxo1-binding sites in naïve CD4⁺ T cells, CD8⁺ T cells, and Foxp3⁺ regulatory T (Treg) cells. By using chromatin immunoprecipitation coupled with deep sequencing (ChIP-Seq), we identified Foxo1 binding sites that were shared among or specific to the three T cell populations. Here we describe the experiments, quality controls, as well as the deep sequencing data. Part of the data analysis has been published by Ouyang W et al. in Nature 2012[1] and Kim MV et al. in Immunity 2013[2], and the associated data set were uploaded to NCBI Gene Expression Omnibus.

JAMM: a peak finder for joint analysis of NGS replicates

MOTIVATION

Although peak finding in next-generation sequencing (NGS) datasets has been addressed extensively, there is no consensus on how to analyze and process biological replicates. Furthermore, most peak finders do not focus on accurate determination of enrichment site widths and are not widely applicable to different types of datasets.

RESULTS

We developed JAMM (Joint Analysis of NGS replicates via Mixture Model clustering): a peak finder that can integrate information from biological replicates, determine enrichment site widths accurately and resolve neighboring narrow peaks. JAMM is a universal peak finder that is applicable to different types of datasets. We show that JAMM is among the best performing peak finders in terms of site detection accuracy and in terms of accurate determination of enrichment sites widths. In addition, JAMM's replicate integration improves peak spatial resolution, sorting and peak finding accuracy.

AVAILABILITY AND IMPLEMENTATION

JAMM is available for free and can run on Linux machines through the command line: http://code.google.com/p/jamm-peak-finder.