Segmentation and genome annotation algorithms

Maxwell W. Libbrecht (0000-0003-2502-0262)1*, Rachel C.W. Chan
(0000-0003-1009-6379)2,3, Michael M. Hoffman (0000-0002-4517-1562)2,3,4,5*

1 School of Computing Science, Simon Fraser University, Burnaby, BC, Canada

2 Department of Computer Science, University of Toronto, Toronto, ON, Canada

3 Princess Margaret Cancer Centre, Toronto, ON, Canada

4 Department of Medical Biophysics, University of Toronto, Toronto, ON, Canada

5 Vector Institute for Artificial Intelligence, Toronto, ON, Canada

These authors contributed equally to this work.

*<EMAIL_ADDRESS><EMAIL_ADDRESS>
## Abstract

SAGA algorithms are widely used to understand genome activity and gene
regulation. These algorithms take as input epigenomic datasets, such as
chromatin immunoprecipitation-sequencing (ChIP-seq) measurements of histone
modifications or transcription factor binding. They partition the genome and
assign a label to each segment such that positions with the same label exhibit
similar patterns of input data. SAGA algorithms discover categories of
activity such as promoters, enhancers, or parts of genes without prior
knowledge of known genomic elements. In this sense, they generally act in an
unsupervised fashion like clustering algorithms, but with the additional
simultaneous function of segmenting the genome. Here, we review the common
methodological framework that underlies these methods, review variants of and
improvements upon this basic framework, catalogue existing large-scale
reference annotations, and discuss the outlook for future work.

## Background and motivation

High-throughput sequencing technology has enabled numerous techniques for
genome-scale measurement of chemical and physical properties of chromatin and
associated molecules in individual cell types. Using sequencing assays, the
Encyclopedia of DNA Elements (ENCODE) Project, the Roadmap Epigenomics
Project, and myriad individual researchers have generated thousands of such
datasets. These datasets quantify various facets of gene regulation such as
genome-wide transcription-factor binding, histone modifications, open
chromatin, and RNA transcription. Each dataset measures a particular activity
at billions of positions, and the collection of datasets does so in hundreds
of samples across a variety of species and tissues. Transforming these
quantifications of diverse properties into a holistic understanding of each
part of the genome requires effective means for summarization. SAGA algorithms
(1) have emerged as the predominant way to summarize activity at each position
of the genome, distilling complex data into an interpretable précis of genomic
activity.

List of textboxfloats 1 Terminology

SAGA. We define a segmentation and genome annotation (SAGA) algorithm as a
procedure that:

  1. 1.

assigns to each position of a whole genome a label (“genome annotation”),

  2. 2.

from a set of multiple ($\geq 3$) classes,

  3. 3.

by

    1. (a)

integrating multiple independent observations at each position, and

    2. (b)

modeling dependence between adjacent positions (“segmentation”).

Previously, researchers have used several other terms to describe this task,
including “segmentation” [1], “chromatin state annotation” [2] and “semi-
automated genome annotation” [3].

Assay. An experiment that produces a measurement at each genomic position,
such ChIP-seq or assay for transposase-accessible chromatin-sequencing (ATAC-
seq).

Label. One of a finite set of classes assigned to each genomic segment that
shares similar activity. Other terms include “state” or “chromatin state”.

Sample. A population of cells on which one can perform an assay, such as a
primary tissue sample or a cell line. Other terms include “cell type”,
“epigenome”, or “biosample”.

Fig 1: Overview of segmentation and genome annotation (SAGA). First,
preprocessing transforms genomic assay sequencing reads into signal datasets.
Second, with signal datasets as input, a SAGA algorithm partitions the genome
and assigns an integer label to each segment, yielding an annotation. Third, a
researcher interprets the labels, assigning a biological interpretation to
each.

SAGA algorithms take as input a collection of genomic datasets, such as ChIP-
seq measurements of histone modifications or of transcription factor binding
(Figure 1). The SAGA task is to use the input datasets to partition the genome
into segments and assign a label to each segment. SAGA algorithms perform this
task in a way that leads to positions with the same label having similar
patterns in the input data.

Year • |  Name or description | References
---|---|---
2007 • |  HMMSeg | [1]
2010 • |  Chromatin colors | [4]
2010 • |  Chromatin states model | [5]
2012 • |  ChromHMM | [2, 6, 7, 8]
2012 • |  Segway | [9, 6, 10, 11]
2013 • |  TreeHMM | [12]
2015 • |  Spectacle | [13]
2015 • |  hiHMM | [14]
2015 • |  Ensembl Regulatory Build (with Segway, ChromHMM) | [15]
2015 • |  EpiCSeg | [16]
2015 • |  Segway+GBR | [3, 17]
2016 • |  IDEAS | [18, 19]
2017 • |  GenoSTAN | [20]
2017 • |  diHMM | [21]
2018 • |  iSeg | [22]
2018 • |  StatePaintR | [23]
2019 • |  RT States | [24]
2019 • |  ConsHMM | [25]
2020 • |  modHMM | [26]
2020 • |  SPIN | [27]
2020 • |  SegRNA | [28]
Table 1: Timeline of selected SAGA methods.

The first SAGA methods were developed in the 2000s, but have increased in
usage recently, thanks to the wide availability of genomic datasets (Table 1).
Large-scale genomic profiling projects such as ENCODE [29] and Roadmap
Epigenomics [7] produced SAGA annotations as a primary output. Researchers
have developed a large variety of SAGA strategies with the goal of improving
upon the basic SAGA framework.

In this review, we summarize the main strategies used by most SAGA methods.
Then, we discuss differences between methods, the challenges they face, and
the outlook for future work.

## Input data

### Experimental assays used for input data

SAGA methods typically use as input a number of different experimental
datasets, each describing some local property of the genome [30]. Such
properties might include chromatin accessibility or presence of some DNA-
binding protein. Although input data initially came from microarray methods
such as tiling arrays [31], they now usually come from sequence census assays
[32].

A standard collection of input datasets might measure histone modifications or
DNA-binding proteins (using assays like ChIP-seq [33] or cleavage under
targets and release using nuclease (CUT&RUN) [34]), and chromatin
accessibility (using assays like deoxyribonuclease-sequencing (DNase-seq) [35,
36] or ATAC-seq [37]). Supplying SAGA SAGA algorithm with datasets that
measure chromatin activity yields an output annotation that captures the
regulatory state of chromatin. Creating these chromatin activity annotations
has served as the predominant use of SAGA methods thus far.

Less frequently, researchers have gone beyond measurements of chromatin and
DNA-binding proteins and have used SAGA methods for other kinds of data. The
output annotation summarizes the input datasets, so the choice of input
greatly influences the annotation’s content and its subsequent interpretation.
SAGA methods can work for any sort of dense linear signal along the genome.
Individual studies have applied it DNA replication timing data [3, 24, 27],
interspecies comparative genomics data [25], and RNA-seq data [28]. Other
studies have even found ways to incorporate non-linear chromatin 3D genome
organization data into the SAGA framework [3, 27].

### Signal representation of genomic assays

Most genomic assay data so far has come from bulk samples of cells. These data
depict a noisy mixture of sampling an assayed property from the many cells
within the population. These cells may represent subpopulations of slightly
different types, or within different cell cycle stages. Thus, each
subpopulation might have different characteristics in the assayed properties.
In the mixture of cell subpopulations, only frequently sampled properties will
rise above background noise. By comparison, less frequently sampled properties
seen in a minority of cells, may remain indistinguishable from background
noise.

Often, the property examined by an epigenomic assay is exhibited or not
exhibited by some position of a single chromosome in a single cell, with no
gradations between the extremes. For example, at some nucleotide of one
chromosome in a single cell, an interrogated histone modification is either
present or it is not. A single diploid cell has two copies of the chromosome.
Thus, at that position, each eudiploid cell can have only 0, 1, or 2 instances
of the histone modification.

Summing or averaging discrete counts over a population of cells leads to a
representation of the assay data called “signal”. Signal appears as a
continuous-scale measurement. Signal arises, however, only from the
aggregation of position-specific properties, which in each cell may have only
a small number of potential ordinal value at the moment of observation.

Unlike epigenomic assays, transcriptomic assays can measure any number of
transcript copies of one position per cell. Despite similar data
representations, one must avoid the temptation to interpret epigenomic signal
intensity as one might interpret transcriptomic signal intensity. For a
transcriptomic assay, greater signal intensity for a transcriptomic assay
might reflect a greater “level” of some transcriptional property within each
cell. For an epigenomic assay, greater signal intensity indicates primarily
that a higher number of cells within a sample have the property of interest.

In both the epigenomic and transcriptomic cases, it remains difficult or
impossible to untangle the contribution to higher signal intensity that arises
from frequency of molecular activity within each cell of a subpopulation from
that from the composition of subpopulations within a whole bulk population.
Improvements in single-cell assays, however, may enable SAGA algorithms on
data from single cells in the near future (see “Outlook for future work”).

### Preprocessing of input data

SAGA methods generally use a signal representation of the input data. This
signal representation originates from raw experimental data, such as
sequencing reads, by way of a preprocessing procedure. For simplicity, we
describe the steps of preprocessing as if a human analyst conducted them all
individually, although some SAGA software packages might perform some steps
without manual intervention:

  * Required preprocessing for all SAGA methods:

    1. 1.

The analyst transforms the experimental data into raw numeric signal data.

       * •

For sequencing data, the analyst:

         1. 1.

Aligns each sequencing read to the reference genome,

         2. 2.

May choose to extend each read to an estimated length of the DNA fragment it
begins, and

         3. 3.

Computes the number of reads per base or extended reads per base for each
genomic position [9, 10].

       * •

For microarray data, the analyst:

         1. 1.

Acquires microarray signal intensity for the experimental sample and for a
control sample, and

         2. 2.

Computes the ratio of experimental intensity to control intensity.

    2. 2.

The analyst chooses units to represent the strength of activity at each
position and may perform further transformation of the raw numeric signal data
into these units.

       * •

For sequencing data, the analyst commonly uses one of:

         * •

Read count (no transformation),

         * •

Fold enrichment of observed data relative to a control [6], or

         * •

$-\log_{10}$ Poisson p-values indicating the likelihood of statistically
significant peaks relative to control [19].

The latter two units can mitigate experimental artifacts because they compare
to a control experiment such as a ChIP input control.

       * •

For microarray data, the analyst commonly performs $\log_{2}$ transformation
of the intensity ratios [38, 4, 39].

  * Optional preprocessing or preprocessing required only for specific SAGA methods:

    1. 3.

The analyst may normalize data to harmonize signal across cell types [40].
Normalization proves especially important when annotating multiple cell types
(see “Annotating multiple cell types”).

    2. 4.

To prevent large outlier signal values from dominating the results, the
analyst may transform signals using one of three variance-stabilizing
transformations of each signal value $x$:

       * •

$\operatorname{asinh}x$ [9],

       * •

$\log_{2}(x+\text{pseudocount})$ [19], or

       * •

an empirical variance-stabilizing transformation [41].

    3. 5.

The analyst may downsample 1-$\mathrm{bp}$ resolution signal into bins (see
“Spatial resolution”). This involves computing one of:

       * •

Average read count,

       * •

Reads per million mapped reads fold enrichment [42],

       * •

Total count of reads [16, 43, 44], or

       * •

Maximum count of reads of each bin [6, 18].

Binning greatly decreases the computational cost of the SAGA algorithm and can
improve the data’s statistical properties.

    4. 6.

The analyst may binarize numeric signal data into presence/absence values [2,
45, 12, 21, 46]. Binarizing signal simplifies analysis by avoiding issues
related to the choice of units, but eliminates all but
$1\text{\,}\mathrm{bit}$ of information about signal intensity per bin.

### Missing data

Genomic assays almost always cannot produce signal for every region of the
whole genome. Regions where an assay cannot provide reliable information about
the interrogated property constitute “missing data” for that assay. Missing
data in sequencing assays may arise due to unmappable sequences, which occur
when repetitive reads do not uniquely map to a region [47, 48]. Missing data
in microarray assays comes from regions covered by no microarray probes. There
are three main ways to treat regions of missing data: (1) by treating missing
data as 0-valued data, (2) by decreasing the model resolution, averaging over
available data so that the missing data has limited impact, or (3) statistical
marginalization over the missing data [49, 9].

When analyzing coordinated assays across multiple cell types, researchers may
have to contend with having no data on some properties within a subset of cell
types. This represents another kind of missing data: one with an entire
dataset missing rather than only data at specific positions. Researchers can
impute [21] entire missing datasets through tools such as ChromImpute [50],
PREDICTD [51] or Avocado [52]. Alternatively they can use SAGA SAGA model with
built-in capability for handling the missing datasets [12].

##  HMM formulation

Many SAGA methods rely on HMM hidden Markov model (HMM), a probabilistic model
of the relationships between sequences of observed events and the unobservable
hidden states which generate the observed events. The structure of HMMs, and
similar models such as dynamic Bayesian networks [53], naturally reflect the
SAGA task of clustering observed data generated by processes that act on
sequences of genomic positions.

### Simple HMM example

As an illustration of a simple HMM, consider a dog, Rover, and his owner,
Thomas. Thomas is 5 years old and too short to see out of the windows in his
home. Rover can leave the house through his dog door and loves taking walks,
playing indoors, and napping. Every morning, he will either wait by the door
for Thomas, play with his squeaky toys, or sleep in. Whichever action he takes
depends on the weather he sees outdoors. For example, on rainy days Rover will
more likely nap or play with his toys indoors.

Thomas must infer the state of the weather outside, hidden to him, based on
the behavior he observes from Rover. Thomas knows the general weather patterns
outside near his home—for example, that rainy weather likely continues across
multiple days.

This scenario fits well into the HMM HMM framework. It has a sequence of
observations (Rover’s behavior) generated by hidden, non-independent
unobservables (the weather outside). One would like to infer the sequence of
hidden unobservables based on the sequence of observations.

### Mathematical formulation

A

$Q_{1}$Hidden$Q_{t}$$Q_{T}$$X_{1}$Observed$X_{t}$$X_{T}$……

B

rrainy$A(\textrm{r}|\textrm{r})=0.7$ $\neg{\textrm{r}}$not
rainy$A(\neg{\textrm{r}}|\neg{\textrm{r}})=0.2$$A(\neg{\textrm{r}}|\textrm{r})=0.3$$A(\textrm{r}|\neg{\textrm{r}})=0.8$
Fig 2: Two representations of HMM HMM. (a) Conditional dependence diagram
representation of an unrolled HMM with sequence of hidden states
${\\{Q_{t}\\}}_{t=1}^{T}$ and sequence of observations
${\\{X_{t}\\}}_{t=1}^{T}$. In this representation, each node represents a
hidden discrete (white rectangle) or observed continuous (grey circle) random
variable. For every index $t$, each hidden random variable $Q_{t}$ takes on
some value $q_{t}$; similarly, each observed variable $X_{t}$ takes on some
value $x_{t}$. $X_{t}$ may represent either scalar or vector observations.
Solid arcs represent conditional dependence relationships between random
variables. (b) State transition diagram representation of Rover and Thomas’s
weather example. In this representation, each node represents a potential
value of the hidden variable $Q_{t}$. The hidden variable takes on values r
(rainy) or $\neg{\textrm{r}}$ (not rainy) on any given day $t$. Solid arcs
represent transitions between hidden states, which have transition
probabilities $A$.

Formally, we can define HMM HMM over time $t\in\\{1,\dots,T\\}$ as follows
[54, 55]. Let the sequence of observed events
$\boldsymbol{X}={\\{X_{t}\\}}_{t=1}^{T}$ consist of each observed event
$X_{t}$ at every time $t$. Let the sequence of hidden states
$\boldsymbol{Q}={\\{Q_{t}\\}}_{t=1}^{T}$ consist of each hidden state $Q_{t}$
at every time $t$. Each $Q_{t}$ takes on a value $q_{t}$ from a set of $m$
possible hidden state values (Figure 2a).

Under the Markov assumption, the probability of realizing state value
$q_{t+1}$ at the next time step $t+1$ depends only on the current state value
$q_{t}$:

$P(Q_{t+1}=q_{t+1}|Q_{t}=q_{t},Q_{t-1}=q_{t-1},\dots,Q_{1}=q_{1})=P(Q_{t+1}=q_{t+1}|Q_{t}=q_{t}).$

We define the transition probability
$A(q_{t+1}|q_{t})=P(Q_{t+1}=q_{t+1}|Q_{t}=q_{t})$, which reflects the
frequency of moving from state $q_{t}$ to state $q_{t+1}$.

We define the emission probability $B(x_{t}|q_{t})=P(X_{t}=x_{t}|Q_{t}=q_{t})$
as the probability that the observable $X_{t}$ is $x_{t}$ if the present
hidden state $Q_{t}=q_{t}$. Specifically, we assume that $B(x_{t}|q_{t})$
depends only on $Q_{t}=q_{t}$, such that

$P(X_{t}=x_{t}|Q_{t}=q_{t},Q_{t-1}=q_{t-1},\dots,Q_{1}=q_{1})=P(X_{t}=x_{t}|Q_{t}=q_{t}).$

Finally, we define the hidden state probability at the first time step as
$\pi_{0}(q_{0})=P(Q_{0}=q_{0})$. We can fully define HMM HMM $M=(A,B,\pi_{0})$
by specifying all of $A$, $B$ and $\pi_{0}$.

In the case of Rover and Thomas, we have $m=2$ possible hidden states (rainy,
not-rainy) and $3$ possible observations (Rover is napping, playing indoors,
or waiting by the door). To Thomas, the hidden variable $Q_{t}$ captures the
weather outside, while the observed variable $X$ captures Rover’s behavior.
The probability of the state of the weather outside changing on a day-to-day
basis is defined by the transition probabilities $A$ (Figure 2b). The
probability of Rover’s behavior, given the weather, is defined by the emission
probabilities $B$.

### Algorithms for inference, decoding, and training

#### Inference

The main task one uses HMMs for is to quantify how well some predicted
sequence of hidden states fits the observed data. Other common tasks like
decoding or training serve as variations of, or build on, this inference task.

In HMMs inference, we can compute the likelihood of any sequence of hidden
states $\boldsymbol{Q}$. We use the sequence of observed events
$\boldsymbol{X}$ and the model probabilities $M$ to compute the likelihood
function $P(\boldsymbol{X}|\boldsymbol{Q},M)$. The likelihood function is the
probability that our predicted sequence of hidden states produced our observed
sequence of observed states. We often compute the likelihood function using
the forward-backward algorithm [56, 57].

#### Viterbi decoding

Given a sequence of observed events $\boldsymbol{X}$, we often wish to know
the maximum likelihood sequence of corresponding hidden states
$\boldsymbol{Q}$. For example, if Thomas observes that in the past 3 mornings,
Rover slept, played, and then slept again, what weather sequence outside is
most likely?

To answer this question, we decode the optimal sequence of hidden states
$\boldsymbol{q}^{*}$ such that
$\boldsymbol{q}^{*}=\operatorname*{arg\,max}_{\boldsymbol{Q}}P(\boldsymbol{Q}|\boldsymbol{X},M)$.
The Viterbi algorithm [58] provides an efficient solution for this problem,
making it unnecessary to compare the likelihood for every possible sequence of
hidden states.

#### Training

Usually, we do not know the model parameters $(A,B,\pi_{0})$ and must learn
them from data. We define training as the process of learning these
parameters, and training data as the sequence of observations upon which we
learn. An efficient algorithm that finds the global optimum parameter values
for some training data does not exist. Instead, researchers commonly train
HMMs using expectation-maximization (EM) [59] algorithms such as the Baum-
Welch algorithm [60], which find a local optimum. Other reviews [54] describe
inference and training methods in more detail.

###  HMMs for SAGA

We can readily apply the HMM formalization to genomic data for use in SAGA
methods. Instead of time, we define the dynamic axis $t$ in terms of physical
position along a chromosome. Each position $t$ refers to a single base pair
or, in the case of lower-resolution models, a fixed-size region (see “Spatial
resolution”). The observation at each genomic position usually represents
genomic signal (see “Input data”). Each position’s hidden state represents its
label (see “Understanding labels”). As a result, decoding the most probable
sequence of hidden states reveals the most probable sequence of labels across
the genome. We call this resulting sequence of labels an annotation.

Many SAGA methods use HMM HMM structure [39, 2, 9, 12, 42, 21, 16, 24], or
generalizations thereof. For example, DBNs are generalizations of HMMs that
can model connections between variables over adjacent time steps. Methods such
as Segway [9] use DBN DBN model in their approach to the SAGA problem. This
can make it easier to extend the model to tasks such as semi-supervised,
instead of unsupervised, annotation [61].

## Understanding labels

SAGA methods are unsupervised. The labels they produce usually begin with
integer designations without any essential meaning. Ideally each label
corresponds to a particular category of genomic element. To make this
correspondence explicit, we must assign a biological interpretation, such as
“Enhancer” or “Transcribed gene”, to each label.

Usually, one makes assignments of labels to biological interpretations in a
post-processing step. In post-processing, a researcher compares each label to
known biological phenomena and assigns an interpretation that matches the
researcher’s understanding of molecular biology. For example, a label
characterized by the histone modification H3K36me3 (associated with
transcription) and enriched in annotated gene bodies might have the
interpretation “Transcribed”. A label characterized by H3K27ac and H3K4me1,
both histone modifications canonically associated with active enhancers, might
have the interpretation “Enhancer” [29].

The interpretation process provides an opportunity to discover new categories
of genomic elements. For example, one SAGA study found that their model
consistently produces a label corresponding to transcription termination
sites. Previously, none had described a distinctive epigenetic signature for
transcription termination [6].

Manual interpretation proves time-consuming for human analysis. Applying SAGA
to multiple cell types independently exacerbates this problem (see “Annotating
multiple cell types”).

Two existing methods automate the label interpretation process: expert rules
and machine learning. In both cases, an interpretation program considers the
information that a researcher would use for interpretation. This includes
examining the relationship between labels and individual input data
properties. It also includes reviewing colocalization of labels with features
in previously created annotations. These annotations may have come from SAGA
approaches or other manual or automated methods.

In the expert rule approach, an analyst designs rules about what properties a
given label must have to receive a particular interpretation. The analyst then
applies these rules to assign interpretations to labels from all models [15].

In the machine learning approach, one trains a classifier on previous manual
annotations. The classifier then learns a model that assigns interpretations
to labels given their properties [11]. One analysis [11] found that automatic
interpretation agreed with manual for 77% of labels, compared to 19% expected
by chance.

## Spatial resolution

Baroque music often employs a musical architecture known as “ternary form”.
Specifically, pieces of this structure follow a general “ABA” pattern,
whereupon the second “A” section recapitulates the first with some variation.
Each section contains multiple musical “sentences”, which may repeat or vary.
Just like linguistic sentences, each musical sentence contains clusters of
notes, or motifs, between “breaths” in the musical articulation. Finer
examination of the motifs shows they contain a few notes and chords each.
Finer examination of the notes themselves shows they behave just like isolated
phonemes in speech, with little meaning on their own.

The genome resembles a musical composition in that one observes different
behaviors at different scales. The scale of genomic behavior one wishes to
observe influences the choice of SAGA method and parameters chosen for the
method. To observe nucleosome-scale behavior such as genes, promoters and
enhancers, one desires $\sim 1\text{\times}{10}^{3}\text{\,}\mathrm{bp}$
segments. To describe behavior on the scale of topological domains [62], one
desires segments of length
$1\text{\times}{10}^{5}\text{\,}\mathrm{bp}1\text{\times}{10}^{6}\text{\,}\mathrm{bp}$
[17, 3, 1].

The most important parameter influencing segment length is the underlying
resolution of the SAGA method. As noted above (see “Input data”), most SAGA
methods downsample data into bins. To observe nucleosome-scale segment lengths
($\sim 1\text{\times}{10}^{3}\text{\,}\mathrm{bp}$), one should use one should
use $100\text{\,}\mathrm{bp}200\text{\,}\mathrm{bp}$ resolution [2, 9, 18]. To
observe domain-scale segment lengths ($\sim
1\text{\times}{10}^{5}\text{\,}\mathrm{bp}$), one should use $\sim
1\text{\times}{10}^{4}\text{\,}\mathrm{bp}$ resolution [4, 3, 27]. Segway [9]
and RoboCOP [63] provide some of few SAGA methods optimized for single-base
resolution inference, and can identify behavior on a $1\text{\,}\mathrm{bp}$
scale. While most existing SAGA methods handle data at just one genomic scale,
there exist methods capable of learning from data at multiple genomic scales
[21].

Limitations of of the experimental data itself influence the choice of SAGA
model resolution. Spatial imprecision in ChIP-seq data gives it an inherent
resolution of about $10\text{\,}\mathrm{bp}$. More precise assays such as
ChIP-exo [64] and ChIP-nexus [65] can approach $1\text{\,}\mathrm{bp}$
resolution. Conversely, assays like DNA adenine methyltransferase
identification (DamID) and Repli-seq have a coarser resolution of $\geq
100\text{\,}\mathrm{bp}$.

The desired scale may also influence the choice of input data. When aiming to
annotate at the domain scale, one should include data with activity at this
scale, such as replication time data and Hi-C data [4, 3, 27, 24]. The
inclusion of long-range contact information from Hi-C data poses a challenge
because standard algorithms for HMMs cannot be used for a probabilistic model
that includes long-range dependencies. Therefore, one must instead use
alternative approaches such as graph-based regularization [3] or approximate
inference [27].

SAGA methods model segment length through their transition parameters. HMM
models assume a geometric distribution in determination of a segment’s length
[66]. Related DBN methods can include constraints to tune segment length
further. Constraints include the enforcement of a minimum or maximum segment
length [9]. Enforcement of a minimum segment length ensures that one does not
obtain segments shorter than the effective resolution of the underlying data
or biological phenomena. Probabilistic models often additionally use a prior
distribution on the transition parameters during training to encourage them to
produce shorter or longer segment lengths.

## Choosing the number of labels

Most SAGA methods require the user to define the number of labels. Using more
labels increases the granularity of the resulting annotation at the cost of
added complexity. Typically, the number of labels ranges from 5–20, with more
recent work favoring 10–15 labels.

One might think to make the choice of number of labels automatically with a
statistical approach. The Akaike information criterion (AIC), Bayes
information criterion (BIC), and factorized information criterion (FIC) [67]
measure the statistical support a particular number of labels has. Instead of
a fixed number of labels, one may give the model flexibility to choose the
number of labels during training and include a hyperparameter that encourages
it to choose a higher or lower number [14]. Or one might define labels
according to local minima in an optimization based on a network model of
assays [46]. One could even exhaustively assign a separate label to every
observed presence/absence pattern in binary data [43].

In practice, however, researchers rarely use these statistical approaches for
determining the number of labels. Optimizing an information criterion does not
necessarily yield the most interpretable annotation. Interpretability reigns
supreme in most SAGA applications. End users find annotations most useful when
they have about 5–20 labels for two reasons. First, most can only articulate
that many known distinctions between classes of genomic elements. Second, even
if one could find meaningful distinctions between a large number of labels,
few using the resulting annotations could keep fine distinctions between such
a large number of patterns in their working memory. [68] Even if a statistical
approach supported the use of 50 labels, the complexity of such an annotation
would make it impractical for most users.

## Annotating multiple cell types

Fig 3: Annotating multiple cell types. (a) Datasets generated by the ENCODE
and Roadmap Epigenomics consortia as of 2019. The black cells represent the
datasets actually generated out of a larger number of potential combinations
of cell type and assay type. (b) Annotating 6 datasets from 3 different
samples: 3 from cell type A, 2 from cell type B, and 1 from cell type C.
Colored letters over signal data indicate data associated with those samples.
One can use three different SAGA strategies with this collection of datasets:
_Independent_ : Performing training and inference completely independently on
each sample. This yields a different annotation for each sample.
_Concatenated_ (horizontal sharing): Training a single model across all cell
types. This yields one annotation per sample with a shared label set. Each
sample must have the same datasets, necessitating imputation of any missing
datasets. _Stacked_ (vertical sharing): Performing training and inference on
datasets from all samples. This yields a single pan-cell-type annotation.

There now exist epigenomics datasets describing hundreds of biological samples
(Figure 3a). Researchers have correspondingly adapted SAGA methods to work
with many samples simultaneously.

We use the term “sample” to refer to some population of cells on which one can
perform an epigenomic assay. A sample could correspond to a primary tissue
sample, a cell line, cells affected by some perturbation such as drug or
disease, or even cells from different species.

The simplest approach for annotating multiple samples involves simply training
a separate model on each sample [11] (Figure 3b). The large number of models
produced by this approach necessitates using an automated label interpretation
process (see “Label interpretation”).

Two categories of approaches aim to share information across samples. The
first, “horizontal sharing” or “concatenated” approaches, share information
between samples to inform the label-training process. The second, “vertical
sharing” or “stacked” approaches, share position-specific information to
inform the label assignment of each position.

### Horizontal sharing: emphasizing similarities across samples for learning
labels

The simplest way to remove the need for interpreting multiple models is to
apply a single model across many samples. To do this, one can treat each
sample as referring to separate copies of a longer genome added horizontally
after the first one in a “concatenated” approach (Figure 3b). One performs
concatenated training and inference little differently than if the data from
different samples pertained to different chromosomes in the same genome.
Because all samples share a single concatenated model, researchers need only
perform post-processing interpretation once.

The concatenated approach has wide usage [69, 6, 7], but has two downsides.
First, concatenated SAGA requires that every sample has data from the same
assays. In practice, this criterion often does not hold true. This means
that—unless these assays are imputed or treated as missing (see “Vertical
sharing: emphasizing similarities across samples in positional
information”)—one must exclude data for an assay conducted in even all but one
samples. In a simple concatenated approach, one cannot annotate a sample which
lacks even one dataset present in the others.

Second, data from different samples can have artifactual differences or batch
effects. Applying the same model across multiple cell types assumes that all
datasets from the same assay type have similar statistical properties. This
can result in label distributions to vary wildly across samples and
biologically implausible sample-specific labels. Data normalization can help
abate the problem of different statistical properties between samples, but
usually does so incompletely. This problem is particularly significant when
using continuous signal. In contrast, binarizing the data (see “Input data”)
can cover up some experimental biases.

One might expect that concatenated annotation would benefit training by
increasing the amount of training data. As it turns out, multiplying the
amount of training data does not significantly aid the training process, as
the types of labels vary little across samples. Most complex eukaryotic
organisms studied with SAGA have very large genomes, and just one sample
provides plenty of training data. In fact, for computational efficiency,
researchers often train on just a fraction of the available samples [7], a
fraction of the genome from a given sample [9] or both.

### Vertical sharing: emphasizing similarities across samples in positional
information

Another class of multi-sample SAGA methods shares position-specific
information across samples as part of the annotation process. These methods
take advantage of the non-independence of biological activity across samples
at a genomic position. For example, if a given position has an active enhancer
label in many samples, it is more likely to act as an active enhancer in a new
sample.

The simplest type of vertical sharing approach learns a model on data from all
samples jointly (Figure 3b). One can implement this “stacked” approach by
adding datasets from all samples vertically into a single combined model. A
stacked model captures patterns of activity across multiple cell types. For
example, a stacked model, unlike an independent model, can find a label for an
enhancer active in cell type A and cell type B but inactive in cell type C.

Although conceptually simple, the stacked approach tends not to work well for
more than several cell types. Stacking fails with larger number of cell types
because each pattern of activity requires its own label. Therefore, the number
of labels must grow exponentially in the number of samples. A simple stacked
model that treats all assays as independent also ignores the relationship
between assays on the same cell type, or the same assay type on different cell
types.

A second approach uses a concatenated model that additionally learns a
position-specific preference over the labels for each position. Through this
preference, data from one sample can influence inference on another. Two
implementations have applied variants of this hybrid horizontal-vertical
sharing approach. First, TreeHMM [12] uses a cellular lineage tree as part of
its input. For each genomic position, TreeHMM models statistical dependency
between the label of a parent cell type and that of a child cell type. Second,
IDEAS [18] uses a similar approach to TreeHMM, but dynamically identifies
groups of related samples rather than using a fixed developmental structure.
The IDEAS model allows these sample groups to vary across positions, which
allows for different relationships between samples in different genomic
regions.

A third approach for vertical sharing uses a pairwise prior to transfer
position-specific information between cell types [17, 3]. In other words, if
position $i$ and position $j$ received the same label in many other samples,
then they should be more likely to receive the same label in an additional
sample. In contrast to the other methods in this section, the pairwise prior
approach does not require the use of concatenated annotation. Therefore, the
pairwise approach has the advantage of not requiring the same available data
in all cell types.

A fourth approach imputes missing datasets in the target cell type, then
applies any of the above annotation methods to the imputed data [50].
Imputation provides three advantages. First, it ensures that all target cell
types have the same set of datasets. Second, one can conduct imputation
entirely as a preprocessing step, allowing the use of any SAGA method
afterward. Third, the imputation process can normalize some artifactual
differences between datasets, making concatenated annotation more appropriate.

Vertical sharing approaches have the downside that one cannot understand the
annotation of each sample in isolation. This arises from the influence on
label assignments in one sample by data from other samples. In particular,
vertical sharing tends to mask differences between samples. For example, if
some position has an enhancer label in many samples, vertical sharing
approaches will annotate that position as an enhancer in a target cell type
too, even with no enhancer-related data in the target cell type.

## Using and visualizing SAGA annotations

Name |  Method | Labels | Samples | Datasets per sample | Reference
---|---|---|---|---|---
Roadmap (5 mark) |  ChromHMM concatenated | 15 | 127 | 5 | [7]
Roadmap (6 mark) |  ChromHMM concatenated | 18 | 98 | 6 | [7]
Roadmap (imputed) |  Impute+ChromHMM | 25 | 127 | 5–12 | [7]
IDEAS |  IDEAS | 20 | 127 | 5 | [19]
Segway Encyclopedia |  Segway independent | 14–33 | 164 | 3–126 | [11]
Segway domains |  Segway+GBR | 5 | 8 | 12 | [3]
Ensembl Regulatory Build |  ChromHMM independent | 25 | 89 | 6–10 | [15]
Table 2: Existing large-scale human reference annotations.

A number of resources can aid in the application of SAGA algorithms and
annotations. Reference annotations exist for many cell types (Table 2). These
obviate the need for a user of the annotation to actually run SAGA SAGA
method. Alternatively, if the user must run SAGA SAGA algorithm on their own
data, standardized protocols describe how to perform this process [70, 8].

Fig 4:  Visualizations of SAGA annotations. (a) Genome browser display showing
164 cell type annotations for a $20\text{\,}\mathrm{kbp}$ region on human
chromosome 15 (GRCh37/hg19) [71]. Each annotation has 8 labels: Promoter
(red), Enhancer (orange), Transcribed (green), Permissive regulatory (yellow),
Bivalent (purple), Facultative heterochromatin (light blue), Constitutive
heterochromatin (black), Quiescent (grey), and Low Confidence (light grey).
(b) Importance score (CAAS) for the same region. Total height at each position
indicates the position’s estimated importance. Height of a given color band
denotes the contribution towards importance of the associated label. (c)
Genome-wide visualization of the SAGA annotation for 164 samples aggregated
over GENCODE [72] protein-coding gene components. Rows: the 9 labels of the
annotation. Columns: gene components, including $10\text{\,}\mathrm{kbp}$
flanking regions upstream and downstream. Each cell shows the enrichment of
the row’s label with a position along the column’s component. Figures derived
from [11].

Most users of SAGA annotations view them through one of three visualization
strategies. The first, and most common, strategy displays individual
annotations as individual rows or “tracks” on a genome browser (Figure 4a). In
each row, the browser displays the segments of that annotation for a region of
one chromosome, usually indicating the label by color. Popular genome browsers
for displaying segmentations include the University of California, Santa Cruz
(UCSC) Genome Browser [73], the Washington University in St. Louis (WashU)
Epigenome Browser [74], and Ensembl [75].

A second visualization strategy integrates annotations of all samples (Figure
4b). This visualization stacks all labels for a given position on top of one
another and scales the vertical axis by an estimate of functional importance
of that position. Two methods can estimate this importance: Epilogos
(https://epilogos.altius.org/), which emphasizes rare activity, and the CAAS,
which measures activity that is correlated with evolutionary conservation
[11].

A third visualization strategy aggregates information about each label across
the entire genome. This shows the enrichment of each label at positions of
known significance, such as gene components (Figure 4c) or curated enhancers.
Tools such as Segtools [76] and deepTools [77] can create these
visualizations.

SAGA annotations can provide valuable reference datasets to other analyses and
tools. The assignment of one and only one label from a small set to every
mappable part of the genome greatly eases downstream analyses. SAGA
annotations summarize genomic activity in a much simpler way than the
individual input datasets, and even than processed versions of the input
datasets such as peak calls.

Most SAGA annotations are in the tab-delimited browser extensible data (BED)
format (https://genome.ucsc.edu/FAQ/FAQformat.html#format1). This makes it
easy to remix SAGA annotations with other datasets using powerful software
such as BEDTools [78]. SAGA annotations form building blocks for methods for
integrative analysis of genomic data such as CADD [79].

## Conclusions and outlook for future work

SAGA methods provide a powerful and flexible tool for analyzing genomic
datasets. These methods promise to continue to play an important role as
researchers generate more datasets. Despite the large existing literature,
future work could still address many challenges.

##### Alternate scales and data types.

Nucleosome-scale annotations
($100\text{\,}\mathrm{bp}1000\text{\,}\mathrm{bp}$ segments) of histone
modifications have wide usage. While annotations of different data types or at
different length scales exist, they are less widely used. Currently, there
exist reference domain annotations with segments of length
$1\text{\times}{10}^{5}\text{\,}\mathrm{bp}1\text{\times}{10}^{6}\text{\,}\mathrm{bp}$
for only a small number of samples [4, 42, 80, 3], and few or no annotations
at other scales.

##### Data preprocessing

. Genome annotations would improve with better processing and normalization of
input datasets. Representations such as fold enrichment used by existing
methods seem primitive compared to more rigorous quantification schemes used
in RNA-seq analysis such as transcripts per million (TPM). One could also
improve SAGA preprocessing by more frequently incorporating information from
multi-mapping reads [81].

##### Confidence estimates.

Most methods do not report any measure of confidence in their predictions. Two
types of confidence would prove useful. First, one would often like to know
the level of confidence that a position in some sample has label X and not
label Y. Second, in many cases one would like to have confidence in a
differential labeling between two samples—that cell type A and cell type B
have different labels. Two methods work towards a solution for the second
problem [82, 83], but there remains much room for further work.

##### Continuous representations.

Existing SAGA methods output a discrete annotation, assigning a single label
to each position. In this discrete approach, annotations cannot easily
represent varying strength in activity of genomic elements or elements that
simultaneously multiple types of activity. A continuous annotation approach
analogous to the topic models used for text document classification might
address this limitation [84].

##### Single cell data.

Existing SAGA methods use data from bulk samples of cells. Increasing
availability of data from single-cell assays necessitates the development of
methods that can leverage this additional information.

##### Pan-cell-type annotation.

The semantics of genome annotations correspond poorly to the way most
molecular biologists conceptualize genomic elements. Most existing annotations
are cell-type-specific—the annotation states that a given locus acts as an
active enhancer in cell type A. In contrast, researchers often state that a
given locus “is an enhancer”.

In contrast, other annotations—such of those of protein-coding genes—serve as
a pan-cell-type characterization. Each gene has a fixed location, and only its
expression varies across samples.

There exists a need for pan-cell-type epigenome annotations. Such an
annotation would define fixed intervals for regulatory elements such as
promoters, enhancers, and insulators, and it would specify in which samples
each element is active. Specifically targeting this task in the SAGA model
could improve results over pan-cell-type annotations assembled from multiple
cell-type-specific SAGA models [11].

## Author contributions

Conceptualization, M.W.L. and M.M.H.; Funding acquisition, M.M.H.;
Investigation, M.W.L., R.C.W.C., and M.M.H.; Visualization, M.W.L. and
R.C.W.C.; Writing — original draft, M.W.L., R.C.W.C., and M.M.H.; Writing —
review & editing, M.W.L. and M.M.H.

## Acknowledgments

This work was supported by the Natural Sciences and Engineering Research
Council of Canada (RGPIN-2015-03948 to M.M.H.)

## Competing interests

The authors declare no competing interests.

## Acronyms

AIC

    Akaike information criterion
ATAC-seq

    assay for transposase-accessible chromatin-sequencing
BED

    browser extensible data
BIC

    Bayes information criterion
CAAS

    conservation-associated activity score
ChIP-seq

    chromatin immunoprecipitation-sequencing
CUT&RUN

    cleavage under targets and release using nuclease
DamID

    DNA adenine methyltransferase identification
DNase-seq

    deoxyribonuclease-sequencing
DBN

    dynamic Bayesian network
EM

    expectation-maximization
ENCODE

    Encyclopedia of DNA Elements
FIC

    factorized information criterion
HMM

    hidden Markov model
SAGA

    segmentation and genome annotation
TPM

    transcripts per million
UCSC

    University of California, Santa Cruz
WashU

    Washington University in St. Louis

## References

  * 1. Day N, Hemmaplardh A, Thurman RE, Stamatoyannopoulos JA, Noble WS.  Unsupervised segmentation of continuous genomic data.  Bioinformatics. 2007;23(11):1424–1426.
  * 2. Ernst J, Kellis M.  ChromHMM: automating chromatin-state discovery and characterization.  Nature Methods. 2012;9(3):215–216.
  * 3. Libbrecht MW, Ay F, Hoffman MM, Gilbert DM, Bilmes JA, Noble WS.  Joint annotation of chromatin state and chromatin conformation reveals relationships among domain types and identifies domains of cell-type-specific expression.  Genome Research. 2015;25(4):544–557.
  * 4. Filion GJ, van Bemmel JG, Braunschweig U, Talhout W, Kind J, Ward LD, et al.  Systematic protein location mapping reveals five principal chromatin types in _Drosophila_ cells.  Cell. 2010;143(2):212–224.
  * 5. Ernst J, Kellis M.  Discovery and characterization of chromatin states for systematic annotation of the human genome.  Nature Biotechnology. 2010;28(8):817–825.
  * 6. Hoffman MM, Ernst J, Wilder SP, Kundaje A, Harris RS, Libbrecht M, et al.  Integrative annotation of chromatin elements from ENCODE data.  Nucleic Acids Research. 2012;41(2):827–841.
  * 7. Kundaje A, Meuleman W, Ernst J, Bilenky M, Yen A, Heravi-Moussavi A, et al.  Integrative analysis of 111 reference human epigenomes.  Nature. 2015;518(7539):317–330.
  * 8. Ernst J, Kellis M.  Chromatin-state discovery and genome annotation with ChromHMM.  Nature Protocols. 2017;12(12):2478–2492.
  * 9. Hoffman MM, Buske OJ, Wang J, Weng Z, Bilmes JA, Noble WS.  Unsupervised pattern discovery in human chromatin structure through genomic segmentation.  Nature Methods. 2012;9(5):473–476.
  * 10. Chan RC, Libbrecht MW, Roberts EG, Bilmes JA, Noble WS, Hoffman MM.  Segway 2.0: Gaussian mixture models and minibatch training.  Bioinformatics. 2018;34(4):669–671.
  * 11. Libbrecht MW, Rodriguez OL, Weng Z, Bilmes JA, Hoffman MM, Noble WS.  A unified encyclopedia of human functional DNA elements through fully automated annotation of 164 human cell types.  Genome Biology. 2019;20(1):180.
  * 12. Biesinger J, Wang Y, Xie X.  Discovering and mapping chromatin states using a tree hidden Markov model.  BMC Bioinformatics. 2013;14(Suppl 5):S4.
  * 13. Song J, Chen KC.  Spectacle: fast chromatin state annotation using spectral learning.  Genome Biology. 2015;16(1):33.
  * 14. Sohn KA, Ho JW, Djordjevic D, Jeong Hh, Park PJ, Kim JH.  hiHMM: Bayesian non-parametric joint inference of chromatin state maps.  Bioinformatics. 2015;31(13):2066–2074.
  * 15. Zerbino DR, Wilder SP, Johnson N, Juettemann T, Flicek PR.  The Ensembl regulatory build.  Genome Biology. 2015;16(1):56.
  * 16. Mammana A, Chung HR.  Chromatin segmentation based on a probabilistic model for read counts explains a large portion of the epigenome.  Genome Biology. 2015;16(1):151.
  * 17. Libbrecht MW, Hoffman MM, Bilmes JA, Noble WS.  Entropic graph-based posterior regularization.  In: Proceedings of the International Conference on Machine Learning; 2015\. p. 1992–2000.
  * 18. Zhang Y, An L, Yue F, Hardison RC.  Jointly characterizing epigenetic dynamics across multiple human cell types.  Nucleic Acids Research. 2016;44(14):6721–6731.
  * 19. Zhang Y, Hardison RC.  Accurate and reproducible functional maps in 127 human cell types via 2D genome segmentation.  Nucleic Acids Research. 2017;45(17):9823–9836.
  * 20. Zacher B, Michel M, Schwalb B, Cramer P, Tresch A, Gagneur J.  Accurate promoter and enhancer identification in 127 ENCODE and roadmap epigenomics cell types and tissues by GenoSTAN.  PLoS One. 2017;12(1):e0169249.
  * 21. Marco E, Meuleman W, Huang J, Glass K, Pinello L, Wang J, et al.  Multi-scale chromatin state annotation using a hierarchical hidden Markov model.  Nature Communications. 2017;8:15011.
  * 22. Girimurugan SB, Liu Y, Lung PY, Vera DL, Dennis JH, Bass HW, et al.  iSeg: an efficient algorithm for segmentation of genomic and epigenomic data.  BMC Bioinformatics. 2018;19(1):131.
  * 23. Coetzee SG, Ramjan Z, Dinh HQ, Berman BP, Hazelett DJ.  StateHub-StatePaintR: rapid and reproducible chromatin state evaluation for custom genome annotation.  F1000Research. 2020;7(214):214.
  * 24. Poulet A, Li B, Dubos T, Rivera-Mulia JC, Gilbert DM, Qin ZS.  RT States: systematic annotation of the human genome using cell type-specific replication timing programs.  Bioinformatics. 2019;35(13):2167–2176.
  * 25. Arneson A, Ernst J.  Systematic discovery of conservation states for single-nucleotide annotation of the human genome.  Communications Biology. 2019;2(1):248.
  * 26. Benner P, Vingron M.  ModHMM: A modular supra-Bayesian genome segmentation method.  Journal of Computational Biology. 2020;27(4):442–457.
  * 27. Wang Y, Zhang Y, Zhang R, van Schaik T, Zhang L, Sasaki T, et al.  SPIN reveals genome-wide landscape of nuclear compartmentalization.  bioRxiv. 2020;2020.03.09.982967.
  * 28. Mendez M, FANTOM Consortium Main Contributors, Scott MS, Hoffman MM.  Unsupervised analysis of multi-experiment transcriptomic patterns with SegRNA identifies unannotated transcripts.  bioRxiv. 2020;2020.07.28.225193.
  * 29. ENCODE Project Consortium.  An integrated encyclopedia of DNA elements in the human genome.  Nature. 2012;489(7414):57–74.
  * 30. Zitnik M, Nguyen F, Wang B, Leskovec J, Goldenberg A, Hoffman MM.  Machine learning for integrating data in biology and medicine: Principles, practice, and opportunities.  Information Fusion. 2019;50:71–91.
  * 31. ENCODE Project Consortium, et al.  Identification and analysis of functional elements in 1% of the human genome by the ENCODE pilot project.  Nature. 2007;447(7146):799–816.
  * 32. Wold B, Myers RM.  Sequence census methods for functional genomics.  Nature Methods. 2007;5(1):19–21.
  * 33. Barski A, Cuddapah S, Cui K, Roh TY, Schones DE, Wang Z, et al.  High-resolution profiling of histone methylations in the human genome.  Cell. 2007;129(4):823–837.
  * 34. Skene PJ, Henikoff S.  An efficient targeted nuclease strategy for high-resolution mapping of DNA binding sites.  eLife. 2017;6:e21856.
  * 35. Boyle AP, Davis S, Shulha HP, Meltzer P, Margulies EH, Weng Z, et al.  High-resolution mapping and characterization of open chromatin across the genome.  Cell. 2008;132(2):311–322.
  * 36. Hesselberth JR, Chen X, Zhang Z, Sabo PJ, Sandstrom R, Reynolds AP, et al.  Global mapping of protein-DNA interactions in vivo by digital genomic footprinting.  Nature Methods. 2009;6(4):283–289.
  * 37. Buenrostro JD, Giresi PG, Zaba LC, Chang HY, Greenleaf WJ.  Transposition of native chromatin for fast and sensitive epigenomic profiling of open chromatin, DNA-binding proteins and nucleosome position.  Nature Methods. 2013;10(12):1213–1218.
  * 38. Schuettengruber B, Ganapathi M, Leblanc B, Portoso M, Jaschek R, Tolhuis B, et al.  Functional anatomy of polycomb and trithorax chromatin landscapes in _Drosophila_ embryos.  PLoS Biology. 2009;7(1):e1000013.
  * 39. Kharchenko PV, Alekseyenko AA, Schwartz YB, Minoda A, Riddle NC, Ernst J, et al.  Comprehensive analysis of the chromatin landscape in _Drosophila melanogaster_.  Nature. 2011;471(7339):480–485.
  * 40. Xiang G, Keller CA, Giardine B, An L, Li Q, Zhang Y, et al.  S3norm: simultaneous normalization of sequencing depth and signal-to-noise ratio in epigenomic data.  Nucleic Acids Research. 2020;48(8):e43.
  * 41. Bayat F, Libbrecht M.  Variance-stabilized units for sequencing-based genomic signals.  bioRxiv. 2020;2020.01.31.929174.
  * 42. Larson JL, Huttenhower C, Quackenbush J, Yuan GC.  A tiered hidden Markov model characterizes multi-scale chromatin states.  Genomics. 2013;102(1):1–7.
  * 43. Taudt A, Nguyen MA, Heinig M, Johannes F, Colome-Tatche M.  chromstaR: Tracking combinatorial chromatin state dynamics in space and time.  bioRxiv. 2016;038612.
  * 44. Zehnder T, Benner P, Vingron M.  Predicting enhancers in mammalian genomes using supervised hidden Markov models.  BMC Bioinformatics. 2019;20(1):157.
  * 45. Hamada M, Ono Y, Fujimaki R, Asai K.  Learning chromatin states with factorized information criteria.  Bioinformatics. 2015;31(15):2426–2433.
  * 46. Zhou J, Troyanskaya OG.  Probabilistic modelling of chromatin code landscape reveals functional diversity of enhancer-like chromatin states.  Nature Communications. 2016;7:10528.
  * 47. Derrien T, Estellé J, Sola SM, Knowles DG, Raineri E, Guigó R, et al.  Fast computation and applications of genome mappability.  PLoS One. 2012;7(1):e30377.
  * 48. Karimzadeh M, Ernst C, Kundaje A, Hoffman MM.  Umap and Bismap: quantifying genome and methylome mappability.  Nucleic Acids Research. 2018;46(20):e120.
  * 49. Lian H, Thompson WA, Thurman R, Stamatoyannopoulos JA, Noble WS, Lawrence CE.  Automated mapping of large-scale chromatin structure in ENCODE.  Bioinformatics. 2008;24(17):1911–1916.
  * 50. Ernst J, Kellis M.  Large-scale imputation of epigenomic datasets for systematic annotation of diverse human tissues.  Nature Biotechnology. 2015;33(4):364–376.
  * 51. Durham TJ, Libbrecht MW, Howbert JJ, Bilmes J, Noble WS.  PREDICTD parallel epigenomics data imputation with cloud-based tensor decomposition.  Nature Communications. 2018;9(1):1402.
  * 52. Schreiber J, Durham T, Bilmes J, Noble WS.  Avocado: a multi-scale deep tensor factorization method learns a latent representation of the human epigenome.  Genome Biol. 2020;21(1):81.
  * 53. Dean T, Kanazawa K.  A model for reasoning about persistence and causation.  Computational Intelligence. 1989;5(2):142–150.
  * 54. Bilmes JA.  What HMMs can do.  IEICE Transactions on Information and Systems. 2006;89(3):869–891.
  * 55. Yoon BJ.  Hidden Markov models and their applications in biological sequence analysis.  Current Genomics. 2009;10(6):402–415.
  * 56. Ferguson JD.  Variable duration models for speech.  Proceedings of Symposium on the Application of Hidden Markov Models to Text and Speech. 1980; p. 143–179.
  * 57. Levinson SE.  Continuously variable duration hidden Markov models for automatic speech recognition.  Computer Speech & Language. 1986;1(1):29–45.
  * 58. Viterbi A.  Error bounds for convolutional codes and an asymptotically optimum decoding algorithm.  IEEE Transactions on Information Theory. 1967;13(2):260–269.
  * 59. Dempster AP, Laird NM, Rubin DB.  Maximum likelihood from incomplete data via the EM algorithm.  Journal of the Royal Statistical Society: Series B (Methodological). 1977;39(1):1–22.
  * 60. Baum LE, Petrie T, Soules G, Weiss N.  A maximization technique occurring in the statistical analysis of probabilistic functions of Markov chains.  Annals of Mathematical Statistics. 1970;41(1):164–171.
  * 61. Chan RC, McNeil M, Roberts EG, Mendez M, Libbrecht MW, Hoffman MM.  Semi-supervised segmentation and genome annotation.  bioRxiv. 2020;2020.01.30.926923.
  * 62. Dixon JR, Selvaraj S, Yue F, Kim A, Li Y, Shen Y, et al.  Topological domains in mammalian genomes identified by analysis of chromatin interactions.  Nature. 2012;485(7398):376–380.
  * 63. Mitra S, Zhong J, MacAlpine D, Hartemink AJ.  RoboCOP: Jointly computing chromatin occupancy profiles for numerous factors from chromatin accessibility data.  bioRxiv. 2020;2020.06.03.132001.
  * 64. Rhee HS, Pugh BF.  Comprehensive genome-wide protein-DNA interactions detected at single-nucleotide resolution.  Cell. 2011;147(6):1408–1419.
  * 65. He Q, Johnston J, Zeitlinger J.  ChIP-nexus enables improved detection of _in vivo_ transcription factor binding footprints.  Nature Biotechnology. 2015;33(4):395.
  * 66. Codogno M, Fissore L.  Duration modelling in finite state automata for speech recognition and fast speaker adaptation.  In: ICASSP’87. IEEE International Conference on Acoustics, Speech, and Signal Processing. vol. 12. IEEE; 1987. p. 1269–1272.
  * 67. Fujimaki R, Morinaga S.  Factorized Asymptotic Bayesian Inference for Mixture Modeling.  In: Lawrence ND, Girolami M, editors. Proceedings of the Fifteenth International Conference on Artificial Intelligence and Statistics. vol. 22 of Proceedings of Machine Learning Research. La Palma, Canary Islands; 2012. p. 400–408.
  * 68. Cowan N.  The magical number 4 in short-term memory: A reconsideration of mental storage capacity.  Behavioral and Brain Sciences. 2001;24(1):87–114.
  * 69. Ernst J, Kheradpour P, Mikkelsen TS, Shoresh N, Ward LD, Epstein CB, et al.  Mapping and analysis of chromatin state dynamics in nine human cell types.  Nature. 2011;473(7345):43–49.
  * 70. Roberts EG, Mendez M, Viner C, Karimzadeh M, Chan RC, Ancar R, et al.  Semi-automated genome annotation using epigenomic data and Segway.  bioRxiv. 2016;080382.
  * 71. Church DM, Schneider VA, Graves T, Auger K, Cunningham F, Bouk N, et al.  Modernizing reference genome assemblies.  PLoS Biol. 2011;9(7):e1001091.
  * 72. Frankish A, Diekhans M, Ferreira AM, Johnson R, Jungreis I, Loveland J, et al.  GENCODE reference annotation for the human and mouse genomes.  Nucleic Acids Research. 2019;47(D1):D766–D773.
  * 73. Kent WJ, Sugnet CW, Furey TS, Roskin KM, Pringle TH, Zahler AM, et al.  The human genome browser at UCSC.  Genome Research. 2002;12(6):996–1006.
  * 74. Zhou X, Maricque B, Xie M, Li D, Sundaram V, Martin EA, et al.  The human epigenome browser at Washington University.  Nature Methods. 2011;8(12):989–990.
  * 75. Zerbino DR, Achuthan P, Akanni W, Amode MR, Barrell D, Bhai J, et al.  Ensembl 2018.  Nucleic Acids Research. 2017;46(D1):D754–D761.
  * 76. Buske OJ, Hoffman MM, Ponts N, Le Roch KG, Noble WS.  Exploratory analysis of genomic segmentations with Segtools.  BMC Bioinformatics. 2011;12(1):415.
  * 77. Ramírez F, Dündar F, Diehl S, Grüning BA, Manke T.  deepTools: a flexible platform for exploring deep-sequencing data.  Nucleic Acids Research. 2014;42(W1):W187–W191.
  * 78. Quinlan AR, Hall IM.  BEDTools: a flexible suite of utilities for comparing genomic features.  Bioinformatics. 2010;26(6):841–842.
  * 79. Kircher M, Witten DM, Jain P, O’Roak BJ, Cooper GM, Shendure J.  A general framework for estimating the relative pathogenicity of human genetic variants.  Nature Genetics. 2014;46(3):310–315.
  * 80. Rao SS, Huntley MH, Durand NC, Stamenova EK, Bochkov ID, Robinson JT, et al.  A 3D map of the human genome at kilobase resolution reveals principles of chromatin looping.  Cell. 2014;159(7):1665–1680.
  * 81. Zeng X, Li B, Welch R, Rojo C, Zheng Y, Dewey CN, et al.  Perm-seq: mapping protein-DNA interactions in segmental duplication and highly repetitive regions of genomes with prior-enhanced read mapping.  PLoS Comput Biol. 2015;11(10):e1004491.
  * 82. Yen A, Kellis M.  Systematic chromatin state comparison of epigenomes associated with diverse properties including sex and tissue type.  Nature Communications. 2015;6:7973.
  * 83. Ebert P, Schulz MH.  Fast detection of differential chromatin domains with SCIDDO.  Bioinformatics. 2020; p. btaa960.
  * 84. Chen B, Kenari NS, Libbrecht MW.  Continuous chromatin state feature annotation of the human epigenome.  bioRxiv. 2018;473017.

  *[ChIP-seq]: chromatin immunoprecipitation-sequencing
  *[ENCODE]: Encyclopedia of DNA Elements
  *[SAGA]: segmentation and genome annotation
  *[ATAC-seq]: assay for transposase-accessible chromatin-sequencing
  *[CUT&RUN]: cleavage under targets and release using nuclease
  *[DNase-seq]: deoxyribonuclease-sequencing
  *[HMM]: hidden Markov model
  *[HMMs]: hidden Markov model
  *[EM]: expectation-maximization
  *[DBNs]: dynamic Bayesian network
  *[DBN]: dynamic Bayesian network
  *[DamID]: DNA adenine methyltransferase identification
  *[AIC]: Akaike information criterion
  *[BIC]: Bayes information criterion
  *[FIC]: factorized information criterion
  *[CAAS]: conservation-associated activity score
  *[UCSC]: University of California, Santa Cruz
  *[WashU]: Washington University in St. Louis
  *[BED]: browser extensible data
  *[TPM]: transcripts per million