Semi-nonparametric modeling of topological domain formation from epigenetic data

Related work

Previous work mainly focused on analyzing epigenetic data in an unsupervised way. Segway [16] and ChromHMM [17] take as input a collection of genomics datasets and learn chromatin states that exhibit similar epigenetic activity patterns which then have different interpretations such as transcriptionally active, Polycomb-repressed. Libbrecht et al. [18] improve Segway predictions by integrating Hi-C data which is not as abundant as histone data, whereas [19] jointly infers chromatin state maps in multiple genomes by a hierarchical model. However, none of these methods deal directly with TADs. Even though a subset of their chromatin states overlap with TADs, predicting TADs from them heuristically does not perform well. Additionally, they either ignore the histone densities, or make parametric distribution assumptions such as geometric or normal which are not always reflected in the true data. When modified to run in a supervised setting, they cannot capture the most informative subset of epigenetic elements.

The recent approach [20] proposes a supervised learning method based on random forests to predict TAD boundaries from histone modifications and chromatin proteins. In general, this approach is reported to perform quite accurately in predicting boundaries. However, it does not model interior TAD segments and it treats each segment independently ignoring the fact that TADs form as a result of the joint effects of multiple segments.

The likelihood function

Let V be the ordered set of genome restriction fragments (bins), where each bin v represents the interval \([vr-r+1\), vr], where r is the Hi-C resolution. Let M be the set of histone modifications (markers) over V. The marker data \(H = (h_{vm})\) is a \(|V| \times |M|\)-matrix where its (v, m)’th entry \(h_{vm}\) is the count of the occurrences of marker m inside segment v. Let \(d=[s, e]\) be a domain (interval) where s and e are its start and end boundaries respectively, \(\{s+1, \ldots , e-1\}\) are the segments inside d, and let \(D = \{[s_{1}, e_{1}], [s_{2}, e_{2}], \ldots , [s_{i}, e_{i}] \}\) be a partition of V where none of the domains overlap.

We propose a supervised, semi-nonparametric, high-dimensional model nTDP that uses H to model and predict D. Our model can be seen as a generalization of Conditional Random Field [21, 22] where we have continuous weights instead of binary features and where we model the marker effects semi-nonparametrically.

Specifically, we assume there are 3 types of segments in V that are relevant for modeling: those that are at the domain boundaries (\(V_{\mathbf {b}}\)), those that are in the interior of domains (\(V_{\mathbf {i}}\)), and those that are not part of a domain (\(V_{\mathbf {e}}\)), and we have \(V = V_{\mathbf {b}}\cup V_{\mathbf {i}} \cup V_{\mathbf {e}}\). For each marker type m, we have 3 types of effect functions, \(f^b_m(c, \mathbf {w_m^b}), f_m^i(c, \mathbf {w_m^i}), f_m^e(c, \mathbf {w_m^e})\), that will describe the relationship between marker count c and the fragment type (b, i, e) for marker type m. Here, \(\bf {w_m^b},\bf {w_m^i},\bf {w_m^e}\) are parameters that we will fit to determine the shape of the effect function. Thus, for example, \(f_m^i(c, \mathbf {w_m^i})\) will describe how a count of c for marker m influences whether the fragment is in the interior (i) of a domain.

Nonparametric form of the effect functions

Because the shape of the marker effect function is unknown, we choose the f functions from the nonparametric family of Bernstein basis polynomials. Bernstein polynomials can approximate any effect function and additionally can handle imposed shape constraints such as monotonicity and concavity.

Training

A nice feature of the objective (3) is that it is convex in its arguments, \(\{\mathbf {w^{b}_{m}}, \mathbf {w^{i}_{m}},\) \(\mathbf {w^{e}_{m}}\}_{m \in M}\), which follows from linearity, composition rules for convexity, and convexity of the negative logarithm. However, training involves several challenges: (a) computing the partition function \(Z^{q}_{|V^{q}|}\) in (3), and (b) estimating W so that the weights are sparse. We solve each of these challenges next.

Training extensions

We can also extend the problem to modeling multiple domain subclasses instead of a single class where domains are categorized into subclasses according to their gene-expression profiles such as TADs with highly-active genes, TADs with repressive genes, etc.

Inferring domains using the trained model

Experimental setup

We binned ChIP-Seq histone modification and DNase-seq data at \(40~\hbox {kb}\) resolution, estimate RPKM (Reads Per Kilobase per Million) measure for each bin, and transform values x in each bin by \(\log (x+1)\), which reduces the distorting effects of high values. In the case of 2 or more replicates, the RPKM-level for each bin is averaged to get a single histone modification file, in order to minimize batch-related differences. We convert any data mapped to hg19 (mm8) to hg18 (mm9) using UCSC liftOver tool. We define TADs over human IMR90, human embryonic stem (ES), and mouse ES cells Hi-C data [8] at 40 kb resolution after normalization by [27]. We use consensus domains from Armatus [28] as the true TAD partition by selecting threshold \(\gamma\) where maximum Armatus domain size is closest to the maximum Dixon et al. [8] domain size (\(\gamma =0.5\) for IMR90, \(\gamma =0.6\) for human ES, and \(\gamma =0.2\) for mouse ES cells).

We solved the training optimization problem by L-BFGS [25]. We use the public implementation of Armatus [28], and obtain histone modifications from NIH Roadmap Epigenomics [29] and UCSC Encode [30]. Code and datasets can be found at http://www.cs.cmu.edu/~ckingsf/research/ntdp. nTDP is reasonably fast: we train on all human IMR90 chromosomes in less than 3 h on a MacBook Pro with 2.5Ghz processor and \(8\hbox {Gb}\) Ram. The iterative procedure in general converges in fewer than 10 iterations.

We prevent overfitting by following a two-step nested cross-validation which has inner and outer steps. The outer K-fold cross-validation, for example, trains on all autosomal human chromosomes except the one to be predicted. Within each loop of outer cross-validation, we perform \((K-1)\)-fold inner cross-validation to estimate the regularization parameters.

nTDP finds a small subset of modifications predictive of TADs

We identified a minimal set of histone marks that can model TADs as follows: we run nTDP independently on each chromosome of human IMR90 to obtain 21 sets of marks. These sets overlap significantly across all chromosome pairs (hypergeometric \(p < 0.05\) for all pair-wise comparisons), and a total of 16 modifications cover all chromosomes. Despite the regularization, the weights of several of these marks are still very close to 0, so we identify a non-redundant subset of the modifications by Bayesian information criterion (BIC) [21] which penalizes model complexity more strongly.

As we increase the number of included modifications from 1 to 16, the BIC decrease nearly stops after 4 modifications, with some additional small reduction up to 6 modifications. The sets of 4 and 6 modifications that were most informative are: {H3K36me3, H3K4me1, H3K4me3, H3K9me3} and {H3K4me3, H3K79me2, H3K27ac, H3K9me3, H3K36me3, H4K20me1}. These non-redundant set of elements are preserved when we repeat this procedure between species. We find that only these 4–6 modifications are needed to accurately predict TADs.

These marks are common in good models

The 4 modifications {H3K36me3, H3K4me1, H3K4me3, H3K9me3} are also enriched among a collection of high quality training solutions. We measure the agreement between estimated and true partitions by normalized variation of information \(NVI=\frac{VI}{\log |V|}\) [31] where VI measures the similarity between two partitions and lower score means better performance. We analyze the fraction of models with 4 histone modifications for which NVI score is at least \(95\%\) of optimum NVI score obtained by running nTDP over all modifications as in Fig. 1. We find 161, 139 solutions satisfying this criteria among 1820 candidates for human IMR90 and human ES histone modifications respectively. We find the 4 histone modifications above to be significantly overrepresented in the set of models for both human IMR90 and ES cells (hypergeometric \(p < 0.0001\)). As a boundary case, restricting the effect function to be a linear function of model parameters in human IMR90 cells does not significantly change the results as in Fig. 2. In another species, mouse ES cells, these 4 histone modifications are also the most informative predictors of TADs as in Fig. 3. These significance values combined with the results above suggest the importance of the identified modifications in TADs.

Predicting TADs from histone marks in human

nTDP is able to predict domain boundaries accurately using 4 histone marks alone in both human IMR90 and human ES cells. We compare TAD prediction performance of nTDP with the chromatin state partition predicted by Segway [16] in terms of NVI. Even though Segway does not predict TADs directly, its chromatin state partition can still be used as a baseline. Training with all 28 histone modifications instead of with the identified 4 modifications does not lead to a major performance increase as shown in Fig. 4a even though it increases the training time approximately 4 times for human IMR90 cells. Restricting the effect function to be monotonic and concave only slightly increases the performance. Chromatin states inferred by Segway do not directly correspond to TADs which leads to a lower TAD prediction performance even though they have other meaningful interpretations.

We find combinatorial effects of histone modifications to be important for accurate domain prediction since none of the modifications can achieve NVI score better than 0.2 when considered independently. To verify that there are not inherent structures in the data that can lead to an easy prediction, we randomly shuffle domains in the training set by preserving their lengths without shuffling modifications, which NVI score is never better than 0.3 in all chromosomes showing the importance of histone modification distributions in TADs.

nTDP also predicts TADs accurately across different species as well as across different cell types as in Fig. 4b–d. For example, if we train on human IMR90 data, the model we obtain is still able to recover domains in human ES cells (Fig. 4b). Using the identified 4 histone modifications achieves NVI score of 0.13 in human ES whereas using all modifications achieves slightly lower 0.11 on average. This performance difference is not significant except chromosomes 7 and 21 in human ES. This holds true across species as well: training on human ES data, for example, produces a model that can work well on mouse ES data. Accurate prediction of TADs by training with the identified 4 histone modifications across different species and cell types suggests the consistency of the identified modifications across species and cell types.

Multiscale analysis of the predicted TADs

We find that our predicted TADs match true TADs more accurately at different scales defined by different Armatus \(\gamma\)’s as in Fig. 5. We observe a slight performance improvement if we define true TAD partition at lower Armatus \(\gamma\) values in human IMR90 which correspond to longer TADs. This figure suggests that some of our wrong TAD predictions may actually correspond to longer TAD blocks which we erroneously interpret as incorrect due to a scale mismatch.