{'fraction_non_alphanumeric': 0.03142177982469172, 'fraction_numerical': 0.03142177982469172, 'mean_word_length': 4.861123204179364, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Short Report on the Markov Property of DNA Sequences on 200-bp Genomic Units of Roadmap Genomics ChromHMM Annotations: A Computational Perspective

The non-coding DNA in eukaryotic genomes encodes a language that programs chromatin accessibility, transcription factor binding, and various other activities. The objective of this study was to determine the effect of the primary DNA sequence on the epigenomic landscape across a 200-base pair of genomic units by integrating 127 publicly available ChromHMM BED files from the Roadmap Genomics project. Nucleotide frequency profiles of 127 chromatin annotations stratified by chromatin variability were analyzed and integrative hidden Markov models were built to detect Markov properties of chromatin regions. Our aim was to identify the relationship between DNA sequence units and their chromatin variability based on integrated ChromHMM datasets of different cell and tissue types.


Introduction
Since large-scale epigenetic datasets such as Encyclopedia of DNA Elements (ENCODE) or Roadmap Genomics became publicly available [1,2], there has been a growing interest in predicting the function of non-coding DNA regions directly from their sequences [3][4][5][6]. The details of the dynamics of chromatin state conversions among different cell types reported that extensive signal variation exists in regulatory regions [7]. And recent studies based on ChromHMM datasets [8,9] provided novel insights into n-gram probabilistic language models for non-coding DNA regions stratified by chromatin variability.
As a follow-up to our preliminary study on ChromHMM datasets of ENCODE [10], we extend our discoveries and continue ongoing efforts to build comparative nucleotide frequency profiles stratified by the chromatin variability. We hope to detect Markov properties by analyzing datasets of the full range of 127 cell and tissue types provided by Roadmap Genomics.
We investigated whether some subsets of the annotated Roadmap Genomics 15-state model stratified by chromatin variability can be predicted by purely making n-gram models of DNA sequences. To do that, ChromHMM blocks of human genome were first dissected into nucleosome resolution of 200-bp units, which accounted for 1,965,764,166 units (127 × 15,478,458 of 200-bp units), and they were integrated into one BED file. Then each individual unit was assigned one dominant chromatin state, by analyzing the integrated BED file of ChromHMM. Next based on chromatin variability of each 200-bp unit-referred to as occurrence frequency-the number of different chromatin states in ChromHMM annotations we divided the units into two groups for each chromatin state: highly variable units and invariable units. By using highly variable 200-bp units as our control group, we were able to isolate some invariable chromatin units that showed strong Markov properties.

Methods
The process is explained in detail in the following sections and the basic steps include: combining 127 BED files into a single file, sorting 200-bp units by the frequency of chromatin variability, filtering out highly variable 200-bp units, building 5th order Markov Models, and evaluating prediction accuracy.

Combining 127 BED files into a single file
The Roadmap Epigenomics consortium has released a 15-state model of BED files from a joint analysis of 111 consolidated epigenomes with 16 additional epigenomes from the ENCODE project [2] in total 127 BED files for public download [11]. We downloaded 125 ChromHMM BED files from the Roadmap Genomics (E001.bed through E128.bed files). We excluded E60.bed and E64.bed files, as of August 20, 2017, to build comparative nucleotide frequency profiles (with human genome GRCh35/hg19) to detect their Markov properties.
When making 15-state ChromHMM BED files, the Roadmap Genomics consortium used a core set of five chromatin markers [2]. We investigated whether some subset of the annotated Roadmap Genomics 15-state model can be predicted by making pure n-gram models of DNA sequences, in reverse. To do that, ChromHMM blocks of human genome were initially dissected into different nucleosome resolution of 200-bp units and each individual unit was assigned one dominant chromatin state, by analyzing the 127 BED files of ChromHMM.
To elaborate on our method, BED format shown at the top of Fig. 1 shows a flexible way to define the data lines that are displayed in an annotation track. The four BED fields shown in each of the BED file represent chrom (the name of the chromosome), chromStart (the starting position of the feature in the chromosome), chromEnd (the ending position of the feature in the chromosome), and state (the 15 chromatin states of Roadmap Genomics, ranging from 1 to 15). For example, the chromatin state of E001 in Fig. 1, for the block from chr1: 9,800 to chr1: 10,600 is 9 (Het heterochromatin state), whereas the chromatin state of E002 in Fig. 1, for the block from chr1: 762,000 to chr1: 763,000 is 1 (TssA proximal promoter state).
For our study, it became critical to develop a functional annotation framework that could be generalized to different cell types. To build good predictive models in making the Markov models of human genomes, we modified the original BED files by dissecting ChromHMM blocks in each BED file into 200-bp units. For example, the original unit of E001 cell line in Fig. 1, ranging from chr1: 9,800 to chr1: 10,600 (a unit size of 800-bp) was dissected into four units of 200-bp blocks (from chr1: 9,800 to chr1: 10,000; from chr1: 10,000 to chr1: 10,200; from chr1: 10,200 to chr1: 10,400; and from chr1: 10,400 to chr1: 10,600), in a new BED file. Likewise, the original E002 unit in Fig. 1, ranging from chr1: 762,000

Filtering out highly variable 200-bp units
Since the details of the dynamics of chromatin state conversions among different cell types was reported it was noted that extensive signal variation exists in regulatory regions [7]. So, we needed a way to quantify signal variation in regulatory regions. Thus, we defined the unit variability count of chromatin states of a given 200-bp unit as the number of states where counts of occurrences were non-zero, to define and compare the observed consistency of each chromatin state at any given genomic position across all 127 epigenomes. Table 1 shows some randomly chosen highly variable and invariable 200-bp units of the integrated BED file, sorted by chromosome number. The degree of chromatin variability is marked as 'H' (high) or 'L' (low), as in the last column of the table. According to the table, the chromatin state variability count of unit chr14: 61,114,600-61,114,800 unit would be eleven, as there were eleven non-zero states (state 1, 2, 3, 5, 7, 10, 11, 12, 13, 14, and 15), and this unit is marked as 'H,' or highly variable.
We reasoned that we could use variability count statistics and the maximum likelihood decision rule to make optimal classifications for the Markov model, since uniform priors can be assumed if we only use 200-bp units with low chromatin state variability. In this way, highly variable 200-bp units where different chromatin states were frequently switched to other states across different tissues and cell types could be either eliminated, or used as controls in datasets, in practicing Markov models.  meaning that all 127 cell lines annotated the same state in these 200-bp units. The variability count of chromatin states of 1,808,431 units was two; meaning that all 127 cell lines were annotated as one of two states in these 200-bp units.
Furthermore, if a state with the number of occurrence count fewer than five was discarded in each of the units, the average variability count of chromatin states drops dramatically. Variability counts of less than 3 states accounted for 93.64% (38.59% ＋ 36.87% ＋ 18.18%) of overall 200-bp units. This means that most of these 200-bp units have a strong preference for a certain chromatin state.
Our model should not be based on a single cell line, therefore; it is critical to propose a new functional annotation framework that can be generalized to different cell types. This gave us good heuristic insight to design new Markov models for our study. A generalizable framework can be achieved through statistically-justifiable models. Based on the newly integrated BED file, we assigned a dominant chromatin state for each of the 200-bp unit, as was explained in the previous section. In this way, it was possible to assign just one or two dominant chromatin states for most of the 200bp units of the entire human genome.  (xi |xi−1, . . . , x1)].

Building two-state fifth order Markov models for each of the 15 chromatin states
Usually, it consists of states corresponding to a biological meaning (e.g., chromatin states) and allows transitions between these states in a biologically meaningful way. The model can define a probability distribution on DNA sequences together with chromatin states.
Our rationale for using HMMs was that invariable units contained in ChromHMM chromatin states can follow a regulatory grammar. They do not form a short fragment of motifs or DNA signatures, but form as a continuous, longer stretch of sequences.
Thus, the 200-bp units were subdivided into 15 chromatin states, based on dominant states. Then, based on chromatin variability of each 200-bp unit, we further divided the 200-bp units into two groups for each of the chromatin state: highly

Results
We directly investigated whether our HMMs based on invariable units have the discriminating power against control datasets of highly variable units. As these HMMs could be used as a binary Naive Bayes classifier, we calculated the sequence of each 200-bp unit that maximized our Markov models. For each chromatin state, the given HMM tries to capture the statistical differences in the two hidden states of 'H' (High) and 'L' (low). As was shown in Fig. 3; based on nucleotide frequency profiles and given a random sequence x1,x2,…,xn, we calculated sequences π1, π2,…,πn of chromatin states that maximized the probability between highly-variable and invariable states, for the 15 × 2 Hidden Markov models. Table 2 shows the prediction accuracy for the 15 chromatin states. We defined a correctly predicted unit as one when predicted results matched the dominant chromatin state of annotations of each of the testing units against control datasets or highly-variable units.
According to Table 2, HMM trained from invariable 1_TssA-dominant units shows the highest prediction accuracy of 85.3%. Investigating the dynamics of chromatin state conversions of highly variable units, we found that most significant state switches are between active states. Among eight active chromatin states, the order of prediction accuracy can be sorted as follows.
Among seven inactive chromatin states, 10_TssBiv and 11_BivFlnk, and 12_EnhBiv showed a high precision rate above 90%.
These results and additional properties of the model suggest that n-grams related to invariant chromatin regions are an inherent and biologically-informative feature of the genome. The framework enables us to infer about coordinated differences in marks by studying chromatin state variability of 200-bp units.

Discussion
We extended our previous study of conditional character-  Table 1). This explains the reason why recall of state 15 was relatively high (91.0% in Table 2). The average frequencies for the state three and six were 20.0 and 36.7, respectively. These states are not usually in a dominant state in most of the 200-bp units, and thus are not appropriate for HMM models.
One limitation of our study is that when building HMM models from 200-bp units of DNA sequences, we did not consider whether the unit should be read forward or in reverse, at this time. We just averaged the n-gram counts of both cases when building emission or transition tables.
Another limitation of our study is that the unit length is set at 200 base pairs without exception. The states of HMMs models based on chromatin states can have an explicit length distribution of the sequence emitted in chromatin state. Future HMM models should consider the DNA sequences of states together with emission lengths. Due to the current lack of data availability this is beyond the scope of this paper.
In conclusion, identifying functional regions in the human genome is a major goal in human genetics. N-gram models, including the most notable n-gram model the Hidden Markov model, have been extensively studied in the field of bioinformatics. However, the Markov property of nucleotide sequences associated with chromatin states at whole human genome scale has rarely been reported in the literature. Consequently, little research has been carried out to explore n-gram models associated with whole genome-wide chromatin maps.
It is important to note that we only used DNA sequences contained in epigenetic datasets in modelling Markov chains. Our study showed that some subsets of the invariable chromatin states possessed strong Markov properties. Though our study is preliminary, it is significant in that it has potential to be used to construct statistical models necessary for developing algorithms to predict function directly from sequence when combined with SNPs, motifs, and other resources in future studies. ORCID: Hyun-Seok Park: https://orcid.org/0000-0002-1

Conflicts of Interest
No potential conflicts of interest relevant to this article was reported.