Learning DNA folding patterns with Recurrent Neural Networks 
Authors: Michal Rozenwald
Aleksandra Galitsyna
Ekaterina Khrameeva
Grigory Sapunov
Mikhail S. Gelfand
Venue:  (2020)
Keywords: Machine Learning
Recurrent Neural Networks
3D chromatin structure
topologically associating domains
computational biology.

Abstract: 
The recent expansion of machine learning applications to molecular biology proved to have a significant contribution to our understanding of biological systems, and genome functioning in particular. Technological advances enabled the collection of large epigenetic datasets, including information about various DNA binding factors (ChIP-Seq) and DNA spatial structure (Hi-C). Several studies have confirmed the correlation between DNA binding factors and Topologically Associating Domains (TADs) in DNA structure. However, the information about physical proximity represented by genomic coordinate was not yet used for the improvement of the prediction models.

In this research, we focus on Machine Learning methods for prediction of folding patterns of DNA in a classical model organism Drosophila melanogaster. The paper considers linear models with four types of regularization, Gradient Boosting and Recurrent Neural Networks for the prediction of chromatin folding patterns from epigenetic marks. The bidirectional LSTM RNN model outperformed all the models and gained the best prediction scores. This demonstrates the utilization of complex models and the importance of memory of sequential DNA states for the chromatin folding. We identify informative epigenetic features that lead to the further conclusion of their biological significance.

Under review as a conference paper at ICLR 2020
LEARNING DNA FOLDING PATTERNS WITH RECURRENT NEURAL NETWORKS
Anonymous authors Paper under double-blind review
ABSTRACT
The recent expansion of Machine Learning applications to molecular biology proved to have a significant contribution to our understanding of biological sys- tems, and genome functioning in particular. Technological advances enabled the collection of large epigenetic datasets, including information about various DNA binding factors (ChIP-Seq) and DNA spatial structure (Hi-C). Several studies have confirmed the correlation between DNA binding factors and Topologically As- sociating Domains (TADs) in DNA structure. However, the information about physical proximity represented by genomic coordinate was not yet used for the improvement of the prediction models.
In this research, we focus on Machine Learning methods for prediction folding patterns of DNA in classical model organism Drosophila melanogaster. The paper considers linear models with four types of regularization, Gradient Boosting and Recurrent Neural Networks for the prediction of chromatin folding patterns from epigenetic marks. The bidirectional LSTM RNN model outperformed all the models and gained the best prediction scores. This demonstrates the utilization of complex models and the importance of memory of sequential DNA states for the chromatin folding. We identify informative epigenetic features that lead to the further conclusion of their biological significance.
INTRODUCTION
Machine Learning algorithms are used nowadays in multiple disciplines. In partieular, the utiliza- tion of these methods in molecular biology has significant impact understanding of cell processes (Eraslan et al. 2019). Investigating the large-scale DÑA structure, i.c. the spatial orga- nization of the genome, chromatin, is of the challenging tasks in the field. The relevance of this research is supported by multiple observations of interconnections between gene regulation, inheritance, disease and chromatin structure (Lupiáñez et al., 2016).
Although the chromatin structure is folded 10a 1 10s times, it maintains fundamental and vital pro- shown to act through the three-dimensional
cesses of the cell. Various regulation mechanisms structure formation. High-throughput experiments capturing contacting fragments of the genome, such as Hi-C, have unravelled many principles of chromosomal folding [Lieberman-Aided et al.) 2009). Although Hi-C-like techniques were developed ten years ago, the experiments of high qual- ity started to be published mainly during the last several years, and the protocol is still elaborate and expensive.
Hi-C has also revealed that chromosomes are subdivided into a set of self-interacting domains called Topologically Associating Domains (TADs) (Ulianov et 2016) that can be seen in Figure 1. TADS were shown to correlate with units of replication timing regulation in mammals (Pope et al. 2014] 4), as well with either active repressed epigenetic domains in Drosophila (Sexton al. 2012).
Various factors were shown to contribute to structure formation. ChIP-Seq is one of the high- throughput experiments dedicated to the detection of factors binding the DNA in vivo. The rapid growth of its data enables exploring the chromatin structure with more sophisticated and complex methods such Machine Learning. The datasets for various factors such as ChIP-Seg experiments for histone modifications become increasingly available in public databases (Ho et al. (2014). The relationship between TADs and epigenetics marks has been investigated recently (Ulianov et al.

---

Under review as a conference paper at ICLR 2020
2016 However, the mechanisms that underlie partitioning of the genome into TADs remain poorly understood. Moreover, there is no comprehensive work investigating all the factors that are publicly available yet.
Figure 1: Typical representation of Hi-C interaction map as genome-wide contact matrix, or heatmap. Bright triangles can be visible across the diagonal. These structures are called TADs (topologically associating domains) and interpreted as compact globules of interacting chromatin. Drosophila melanogaster S2-DRSC cells Chr 3R 22.8 -24.6 Mb
This study focuses bringing insights into the 3D chromatin structure using Machine Learning- The goal is to explore the principles of TAD folding and the role of epigenetics in this process. To that end, the analysis of Drosophila melanogaster chromatin was performed using Linear Regression models and Recurrent Neural Networks. Quality metrics were calculated, and informative features were investigated identify which of chromatin marks are most significant in predicting information about TADs.
In addition, the same techniques might be used to explore the 3D chromatin structure of mammals and humans in particular. Such reconstruction of the information about Hi-C map might be useful not only for understanding the chromatin structure formation but can also have various practical medical applications. For example, gliomagenesis and limb malformations in humans were demon- strated to be caused by chromosomal topology disruption (Krijger De Laat 2016).
LITERATURE REVIEW
Over the last decade, the volume of produced data has significantly increased and brought the op- portunity of applying complex and efficient methods. Several other studies were focused pre- dicting the 3D chromatin architecture using Machine Learning methods. One approach to this prob- lem is to the Hi-C map as input of the model, for example, Cristescu et al. (Cristescu et al.) 2018) presented the REcurrent Autoencoders for CHromatin 3D structure prediction (REACH-3D). REACH-3D reconstructs the chromatin structure, recovers several biological properties and have high correlation with microscopy measurements.
However, another approach is predict the information about the chromatin structure from other types of biological characteristics. In particular, Schreiber et al. (Schreiber et al. 2018) considered nucleotide sequence as input for deep Convolutional Neural Network. The objective of this chitecture was to estimate the Hi-C contacts. This Neural Network demonstrated that the predicted outcomes are related to histone modification, selected functional elements and replication timing which correlates with theoretical knowledge.
Moreover, another work that inspired this research was made by Ulianov al. (Uianov et al., 2016). They suggested that active chromatin and transcription play key role in chromosome par- titioning into TADs. It was shown that numerous transient interactions between nucleosomes of inactive chromatin lead to the formation of TADs that are potentially highly dynamic self-organized structures. On the other hand, nucleosomes, that tend to interact less often, influence the formation of inter-TADs and TAD boundaries. Ulianov et al. showed that active chromatin marks were prefer- ably present at TAD borders, and repressive histone modifications that reflect nucleosomes occu- pancy were depleted within inter-TADs, which reveals the correlation between TADs and chromatin
2

---

Under review as a conference paper at ICLR 2020
marks. Fortin et al. in (Fortin & Hansen 2015) succeeded in extracting knowledge from ChIP-See data of histone modification to analyze the chromatin structure. They constructed predictive model of the Hi-C that unrevealed the correlation with replication timing, which proves the hypothesis of the possibility of extracting information about the Hi-C contacts from Nucleotide Sequence and DNasel assay signal of Homo sapiens cell lines. A principal difference from all described works is that to our model explores the 3D chromatin characteristics of Drosophila melanogaster using set of ChIP-Seq data input. To the best of our knowledge, no other published work was conducted to predict Topologically Associated Domain characteristics from epigenetic marks.
DATA
3.1
INPUT DATA
Hi-C datasets for Drosophila melanogaster S2-DRSC cells were collected from Ulyanov et al. (Ulianov et al., 2016). Drosophila dm3 genome assembly subdivided into 5950 sequential genomic regions called bins, where each bin coresponds to 200 000 (20-Kb) DNA base pairs. Each bin be described by number of epigenetic features, estimated by ChIP-Seq We downloaded all available epigenetic datasets at the moment from the modENCODE database (Celniker et al. 2009) and processed identically (Ulianov et al., 2016).
Based on the current model of chromatin formation in Drosophila, we distinguish two ChIP-Seq sets. The first set has five biologically significant features: Chriz, CTCF, Su(Hw), H3K27me3, H3K27a. The second set contains Chriz, CTCF, Su(Hw), H3K27me3, H3K27a, BEAF-32, CP190, Smc3, GAF, H3K36me1, I, H3K36me3, H3K4mel, H3K9ac, H3K9mel, H3K9me2, H3K9me3, #14K16ac.
For normalization of the input data, each feature was centred by mean and scaled by variance. The example eight original ChIP-Seq features and their transformation is seen Appendix.
3.2 TARGET VALUE
Topologically Associating Domains (TADs) can be represented as the segmentation of the genome into discrete regions. However, this segmentation dependent on one or several parameters, cor- responding to the characteristic size of TADs. We sought for avoiding the problem of parameters selection in our approach. Thus we adopted the approach from (Ulanov et al., (016) and calculated the local characteristic of TAD formation of the genome, namely, gamma transitional.
The procedure of calculation is briefly described below. Armatus software (Filippova et al. 2014) used annotate Topologically Associating Domains (TADs) with scaling parameter gamma that determines the average size and the number of TADs. When gamma is fixed, each genomic bin is annotated as part of a TAD, miD or TAD boundary, as part of segmentation. We characterized each bin by the scaling parameter gamma at which this bin switches from being part of a TAD to being a part of an inter-TAD or TAD boundary. Given the higher the gamma value, the smaller the TADs are in the Armatus annotation. See the illustration in Figure 2.
Whole-genome Hi-C maps of chromatin folding set of S2-DRSC Drosophila cells were taken from and processed similarly to (Ulianov et al. 2016
METHODS
4.1 PROBLEM STATEMENT
To avoid ambiguity, let us clearly define our Machine Learning problem.
The objects are "bins" DNA sections of the length of 20,000 nucleotides with intersection of Drosophila melanogaster (see Introduction and Section 3.1 for more details).
-The features are ChIP- Seq epigenetic data on chromatin markers (Section 3.1).
-The target value transitional gamma- parameter transformation from TAD to inter-TAD, TAD
boundary (Section 3.2).
3

---

Under review as a conference paper at ICLR 2020
Heatmap of (Nam of reads)
Histegram of gamma values
3000
3400
7-
y=
22.8
23.2
23..
24
24.4
Mb
Chr
Drasophila mefamoguster, S2-DRSC c.l.
Tranritinaal gumma this point would be ie have only two fixod gomsa and turned to inter TAD or TAD boundary with the smaller gomma of
they were part of
Annotation of TADs different gamma parameter values is the left side. Higher gamma values correspond to smaller TADs. Transitional gamma is the value of gamma at which genomic bin switches from being a part of TAD to being part of an inter-TAD The histogram of the target value transitional gamma in the data is presented in the right part of this plot.
Figure 2:
TAD boundary.
The task to predict the characteristics of the 3D structure of chromatin transitional gamma. The aim is to identify which of chromatin marks (ChIP-Seq data) are most significant in predicting information about the Topologically Associating Domains (TADs).
4.2 LOss FUNCTION
As described in Section 3.2, the target object, transitional gamma, is continuous value from 0 to 10, which leads to solving a regression problem. The classical optimisation function in this type of problems is Mean Square Error (MSE). However, the distribution of the target is significantly unbalanced (Figure 2). The target value of most of the objects is in the diapason between 0 and 3. Nevertheless, the contribution of the error on objects with high also be high in the total score using Mean Square Error. Moreover, the biological nature of objects with high value of the transitional gamma is different from other objects. For DNA bins with transitional gamma value equal to 10, gamma value at which this bin passed from the TAD state to the Onter-TAD TAD boundary was found. To build model that accurately predicts the values of the transitional gamma for most objects, we have introduced our own custom loss function called modified weighted Mean Square Error (wMSE). might be reformulated as MSE multiplied by the weight (penalty) of the error, depending on the true value of the target variable.
value of the target will
n wMSE 72
(Verues Upred, )2 a
)²a Strues
&
where is the number of data points, Strue; is the true value for data point number i, Upred, is the predicted value for data point number i,
is equal to the maximum value of the ytrae values increased by 1 to avoid multiplying the error with 0.
As a result the model is penalized less for errors on objects with high value of the transitional gamma by using the weighting. The maximum values of the target value in the transitional gamma dataset is 10, thus a is equal to 11.
4.3 MODELS
To explore the relationships between the 3D chromatin structure and epigenetics data, we built Lin- Regression (LR) models, Gradient Boosting (GB) Regressors and Recurrent Neural Networks (RNN). The LR were applied with no regularization, either LI and L2 regularization or both of them. All the models were trained using the wMSE loss function. The Linear models chosen
4

---

Under review as a conference paper at ICLR 2020
to create benchmark for this problem other results of ML pipeline is publicly available for this dataset. also allows intuitive feature importance interpretation. It is worth mentioning that our input bins are sequentially ordered in the genome. Due to DNA connectivity and local properties of clustering, the target variable values might be vastly correlated. Thus, in order to increase the chance of learning this property of the biological data, we selected RNN models.
DNA is a long structured molecule formed out of nucleotides arranged in a linear sequence. DNA double-stranded which means each nucleotide has complementary pair, together called a base pair. DNA molecule might be several million base pairs (Mb) long and serves as the storage and the means of utilization of genetic information. The information content of DNA is equivalent read in forward and reverse direction, thus all local properties of its sequence should be independent of the selected direction. To use this property of DNA molecule, we selected bidirectional L.STM RNN architecture (Schuster & Paliwal 1997). The index of the middle bin is calculated as the floor division of the length of the input by? 2.
The variable parameters that we investigated in our LSTM model are:
sequence length of input RNN objects is a set of consecutive DNA bins with fixed input length called window size from to 10.
Number of LSTM Units: 1, 8, 16, 32, 64, 128, 256.
Number of training Epochs: 1, 4, 8, 16, 32, 64, 128, 256, 512. Early Stopping to automatically identify the optimal number of training epochs was used for the final models. Loss function: weighted Mean Square Error (wMSE), our custom evaluation function defined in Section 4.
Optimizer: Adam. The data was always randomly separated into three groups: train dataset 70%
of data, 20% test dataset and 10% for validation.
%
hp
IN Notiled
IV Natures)
Figure 3: Scheme of Bidirectional LSTM Recurrent Neural Networks with one output that was implemented. The values of .., are the DNA bins with input window size t, {hi,.., h2] are the hidden states of the RNN model, Ut/2 represents the corresponding target values Transitional gamma of bins the middle bin F2/2-
RESULTS
For each type of the models, we have performed training several times to get more consistent results. The results are presented for the observations of ten experiments. The weighted Mean Square Error (wMSE) that is defined in Section 4.2 was calculated for each experiment. The best score of the weighted Mean Square Error using Linear Regression with LI and L2 regularization (Elastic Net model) with parameter alpha equal to 0.2 was performed using a grid search. The wMSE of these experiments on train and test datasets was found and is presented in Table 1. The values of MSE, MAE and R² can be found in Table and Figure 10, where LR Linear Regression models, GB-X Grad Boosting models with X estimators. Feature importance can be analyzed by exploring the weight coefficients of the Linear models. The prediction is created based on the multiplication of each weight on the corresponding feature. Thus, larger absolute values of the feature result in the stronger influence of this particular feature on the prediction of the model. Thus we were able to extract the prioritization in terms of the influence of the features. After performing experiments on the first dataset with five ChIP- Seq characteristics, the resulting weights happen to be significantly stable as it is shown below in the table of feature coefficient of Linear Regression (Figure 4). As
5

---

Under review as a conference paper at ICLR 2020
result, we obtain that the most valuable in terms of the absolute value of the feature weight is Chriz,
then CTCF, H3K27ac and H3K27me3, when the weight of Su(Hw) is the smallest.
Chriz
CTCF
H3K27ac
H3K27me3
Su(Hw)
0
-0.48
-0.31
-0.21
0.19
0.05
-0.11
-0.23
2
-0.43
0032
-025
0.20
0.03
3 -0.45
-0.31
-0.2
0.20
0.04
4 -0.48
-0.31
-0.23
0.21
0.07
Figure 4: Weights of features of trained Linear Regression models (rows) on the dataset with 5
ChIP-Seq features (columns)
We adopted the same approach to the second dataset of ChIP-Seq characteristics. In comparison to the same application on dataset of five features the coefficient order by absolute weights values less stable (a table with sorting of the indexes of features by their weights can be seen in Sup- plementary Section). The numbers of occurrences of each of the feature indexes in the list of most influencing features were calculated. We have sorted the features based on this frequency number. Chriz was proved to be the most robustly reproduced influential factor. CTCF and CP190 were identified as the second degree on the scale of significant factors.
Another result worth mentioning is the selection of only one important feature Chriz out of both datasets while using the model of Linear Regression with L1 regularization (visualization can be found in the supplementary materials).
We implemented also Gradient Boosting (GB) for regression. The GB additive model has outper- formed the linear regression models. However, there was a strong tendency to over-fitting for wide range of variable parameters such as the number of estimators, learning rate, maximum depth of the individual regression estimators, minimum number of samples required split an internal node. The best results were observed while setting the 'n estimators': 100, 'max depth': 3 and n_estimators': 250, "max mexx_ddpph' 4, 'learning_rate" rate": 0.01 and they presented bellow in Table 1.
The main Neural Network that we were exploring is Bidirectional LSTM. As described above, the sequential relationship of the input objects in terms of the physical distance in the DNA justifies the usage of Recurrent Neural Networks. For each variation of parameters, experiments were conducted and evaluation metrics calculated (tables with results be found in Section 9).
To explore the dependencies of weighted Mean Square Error on the sizes of sequence length, Bidi- rectional RNN models were trained with different input window size and number of LSTM Units. The result is shown in Figure 5 where an optimal sequence size equals to input window size 6 and 64 LSTM Units was revealed.
This result has clear biological interpretation as the typical size of TADs from around 120 Kb, which corresponds to bins of 20,000 that turned out to have the strongest prediction scores.
As result, the Bidirectional LSTM Recurrent Neural Networks with 64 LSTM Units and sequence of 6 bins taken as input data were trained and achieved better evaluation scores than a constant prediction, Linear models and Gradient Boosting models (Table 1). The constant prediction was made using the mean value of the training dataset.
To explore the importance of each feature X from the input space, replaced the values of the corresponding column of the feature matrix with zeros. Further, we calculate the evaluation metries and check how significantly different they are from the results obtained while using the complete set of data (Figure 8).
The results of wMSE on the test set do not differ dramatically from using full dataset. When we drop out each of the five features, we get the same score of around 0.9 that is almost equal to using all of them together. This means that our RNN is able to achieve the same score with subset of these four features out of all five. The results of applying the same technique while omitting each feature one by using the second dataset of ChiP-Seq features allowed the evaluation of the biological

---

Under review as a conference paper at ICLR 2020
Weighted MSE SEQUENCE LENGTH change Train
Weighted MSE by NUMBER of LSTM UNITS change on Train
371
039
8985
222
3875
15
100
200
556
Inpul Bequaincal langts
Number of LSTM Units
Weighted MSE SEQUENCE LENGTH change on
Weighted MSE by NUMBER LSTM UNITS change C
UPDATE
weightrr Jergee
214
2
10
100
120
256
Lengih
Figure 5: Weighted Mean Square Error of trained Bidirectional LSTM. The upper of graphs present the results for the train dataset, the bottom row shows wMSE counted on the test objects. The left half shows the results of training RNN with 64 units for different sizes of sequence length. The rights half shows wMSE for training RNNs with an input sequence of bins for different number of LSTM Units. The green box is surrounding the optimal scores.
Table 1: Weighted MSE of all Models
5 features Test
18 features Train
Train
Test
Constant prediction
1.61
1.62
1.61
1.62
Linear Regression
1.19
1.19
1.13
1.13
Linear Regression +
1.16
1.16
1.10
1.10
Linear Regression + L2
1.19
1.18
1.13
1.14
Linear Regression + LI + L2
1.16
1.16
1.10
1.10
Grad Boosting 100 estimators
1.13
1.13
1.07
1.08
Grad Boosting 250 estimators
1.06
1.06
0.98
0.98
biL.TMM 64 units bins
0.85
0.85
0.72
0.72
impact of the features. These wMSE scores are presented in Figure 6 as well as the result of training the model on all features together.
The difference between the wMSE using all the features and omitting each one separately is pre- sented on Figure This provides the opportunity of identifying how valuable a particular biological characteristic is using RNN.
CONCLUSION
The ChIP-Seg data usage for chromatin folding patterns prediction was confirmed by training ML models with dignified evaluation scores. Moreover, the results were interpretable and biologically relevant.
Linear Regression models, Gradient Boosting Trees and Recurrent Neural Networks were for the first time applied new dataset of chromatin characteristics. All models have performed better than constant prediction with the mean value of the training dataset. The utilization of memory of previous states linearly ordered by DNA molecule improves the prediction significantly as the best

---

Under review as a conference paper at ICLR 2020
Table 2: Evaluation scores for all models
5 features
model type
MSF
MSE
MAE
MAE
R
Train
Test
Train
Test
Const
3.73
3.52
1.36
1.31
0
LR + L1
2.84
2.72
1.11
1.1
0.24
LR + L2
2.76
2.66
1.11
1.1
0.26
LR + +1 + L.2
2.79
2.68
1.1
1.09
0.25
GB-250
2.27
2.29
0.98
0.38
biLSTM RNN
2.36
2.9
0.92
1.01
0.33
18 features
LR + +1
2.72
2.63
1.07
1.07
0.27
LR + L2
2.61
2.55
1.07
1.09
0.3
LR + +1+
2.68
2.61
1.07
1.07
0.28
GB-250
2.27
2.29
0.98
0.38
biLSTM RNN
2.03
2.45
0.85
0.9
0.43
WMSE of 18 features but one
differefee using all but one features
1.0
100
100
065
00
860
889
Figure 6: On the left the weighted Mean Square Error Train and Test dataset while dropping out one of the input ChIP-Seq features and wMSE of training using all features. On the right the difference between the weighted Mean Square Error using all the features and dropping each one separately. Blue bars correspond to the wMSE a testing dataset, red training. Bidirectional LSTM with 64 units and six input bins was used.
results were obtained by bidirectional L.TTM RNN model. The optimal input window size was also equal to six which has biological meaning as it strongly aligns with the average TAD length.
Feature importance analysis of the input ChIP-Seq data was conducted. The Linear models weights provided biologically meaningful prioritization of the ChIP-Seq. Moreover, after training Linear Regression with L1 regularization detected one ChIP-Seq feature Chriz on both of the datasets as the most influencing. The results of applying Neural Network models allowed the evaluation of the biological impact of the features.
Exploration of the transferability of the models between different cell types and species might be an interesting development this work. More input features different biological nature, such as DNA sequence itself, is another direction of research.
The code is open sourced and the implemented pipeline can be easily adapted to any similar biolog- ical dataset of chromatin features.
8

---

Under review as a conference paper at ICLR 2020
REFERENCES
Susan E Celniker, Laura AL Dillon, Mark B Gerstein, Kristin C Gunsalus, Steven Henikoff, Gary H Karpen, Manolis Kellis, Eric C Lai, Jason D Lieb, David M MacAlpine, al. Unlocking the secrets of the genome. Nature, 459(7249):927, 2009.
Bianca-Cristina Cristescu, Zalán Borsos, John Lygeros, María Rodríguez Martínez, and Maria Anna Rapsomaniki. Inference of the anee-dimensional chromatin structure and its temporal behavior. 2018.
arXiv preprint arXiv:1811.0960
Gökcen Eraslan, Ziga Avsec, Julien Gagneur, and Fabian J Theis. Deep learning: new computational modelling techniques for genomics. Nature Reviews Genetics, pp. 1, 2019.
Darya Filippova, Rob Patro, Geet Duggal, and Carl Kingsford. Identification of alternative topolog- ical domains in chromatin. Algorithms for Molecular Biology, 9(1):14, 2014.
Jean-Philippe Fortin and Kasper Hansen. Reconstructing a/b compartments as revealed by hi-c using long-range correlations in epigenetic data. Genome biology, 16(1):180, 2015.
Joshua WK Ho. Youngsook L. Jung, Tao Liu, Burak H Alver, Soohyun Lee, Kohta Ikegami, Kyung- Ah Sohn, Aki Minoda, Michael Y Tolstorukov, Alex Appert, et al. Comparative analysis of metazoan chromatin organization. Nature, 512(7515):449, 2014.
Peter Hugo Lodewijk Krijger and Wouter De Laat. Regulation of aassase-associated gene expression in the 3d genome. Nature reviews Molecular cell biology, 17(12)771, 2016.
Erez Lieberman-Aiden, Nynke I. Van Berkum, Louise Williams, Maxim Imakaev, Tobias Ragoczy, Agnes Telling, Ido Amit, Bryan Lajoie, Peter J Sabo. Michael O Dorschner, et al. Comprehen- sive mapping of long-range interactions reveals folding principles of the human genome. science, 326(5950);289-29)) 2009.
Darío G Lupiáñez, Malte Spielmann, and Stefan Mundlos. Breaking tads: how alterations of chro- matin domains result in disease. Trends in Genetics, 32(4):225-237, 2016.
Benjamin D Pope, Tyrone Ryba, Vishnu Dileep, Feng Yue, Weisheng Wu, Olgert Denas, Daniel Vera, Yanli Wang. R Scott Hansen, Theresa K Canfield, al. Topologically associating domains are stable units of replication-timing regulation. Nature, 515(7527):402, 2014.
Jacob Schreiber, Maxwell Libbrecht, Jeffrey Bilmes, and William Noble. Nucleotide sequence and dnasei sensitivity predictive of 3d chromatin architecture. bioRxiv, pp. 103614, 2018.
Mike Schuster and Kuldip Paliwal. Bidirectional recurrent neural networks. IEEE Transactions on Signal Processing, 45(1 1):2673-2681, Nov 1997. doi: 10. 109/78.650093.
Tom Sexton, Fitan Yaffe, Ephraim Kenigsberg, Frédéric Bantignies, Benjamin Leblanc, Michael Hoichman, Hugues Parrinello, Amos Tanay, and Giacomo Cavalli. Three-dimensional folding and functional organization principles of the drosophila genome. Cell, 148(3):458 472, 2012.
Sergey V Ulianov, Ekaterina E Khramceva, Alexey A Gavrilov, Ilya M Flyamer, Pavel Kos. Elena A Mikhaleva, Aleksey A Penin, Maria D Logacheva, Maxim Imakaev, Alexander Chertovich, et al. Active chromatin and transcription play key role in chromosome partitioning into topo- logically associating domains. Genome research, 26(1):70-84, 2016.

---

Under review as a conference paper at ICLR 2020
A
APPENDIX
Before scaling
After scaling
herrs
300
Gin 1
Geored
hetres
Rearhh
Figure 7: Histograms of (A) original ChIP-Seq data and (B) normalized data of first 8 ChIP-Seq features of all input bins.
Results of input all features bul one Soq san.s wat Exll Staping:
missing_
npopoc
Wmsse
Wmso trai
ature
n_tries
hs
Chrizz
mean
21
0.9 0.01
0.87 0.01
CTCF
mean
30
0.89 0.89&00.01 0.01
0.04
Su(Hw)
32
0.89 0.02
0.63 0.03
H3K27me3
mean
32
0.89
0.01
0.00
H3K27ac
mean
32
0.89 0.01
0.85
FULL
27
0.91
0.00
Figure 8: Weighted Mean Square Error Train and Test dataset while dropping out one of the input ChIP-Seq features and using this full dataset. Bidirectional LSTM with units with six input bins with Early Stopping was used.
10

---

Under review as a conference paper at ICLR 2020
CPR8
IOHEN!
13R2nd)
SEAL-
327
1.5
125
322
1310
6.9
-2.5
014
613
E12
O/E
10000
337
+25
423
-25
-1.0
1.15
-23
418
E13
513
012
300
430
4143
120
435
439
125
120
3.15
018
E14
019
01
038
027
-3X
10001N
137
125
124
120
E19
488
08
619
612
022
028
088
0.8
613
410
Figure 9: Weights of features trained Linear Regression models (corresponding to the rows) of the second dataset of ChIP-Seg characteristics (columns).
Chip-Seq
featuree
model type
mse_train
nse_test test
nae_train
nee_test
P"2
WMSE_tre train
wHSE_test
const
3.73
3.52
1.36
1.31
1.61
1.61
8.25
best
8.24
1.16
2.76
2.66
0.25
1.19
1.19
LR
L1 L2
3.13
2.97
1.19
1.16
8.16
1.26
1.26
LR
L1
L2 best
2.79
2,68
1.1
1.89
0.25
1.16
1.16
GB-189
2.78
2.84
1.12
1.11
0.25
1.13
1.19
08-250
2.45
2.48
1.63
0.34
18 Chip-Seq features
model
train
r*2
wMSE, train
WHSE, test
const
3.73
3.52
1.36
1.31
1.61
1.61
LR
2.61
2.55
1.09
1.87
0.3
1.13
1.13
3.61
3.42
1.32
1.29
8.03
1.54
1.54
best
2.72
1.87
1.67
0.27
1.1
1.1
L2 best
2.68
2.61
1.87
0.28
1.1
1.1
0B-109
2.69
2.73
1.00
1.88
0.27
1.98
1.98
GB-259
2.27
2.29
8.98
0.38
0.98
8.98
b1LSTH ANN
2.83
2.45
6.85
0.9
8.43
8.72
8.72
Figure 10: MSE, MAE, R², weighted MSE metrics for various ML models experiments. Here "LR" stands for Linear Regression models, "GB-X" Grad Boosting models with X estimators, "* best" means that the presented scores for the best of models of type *
11


Reviews:
[{"review_id": "ryg8ROfaYH", "reviewer": "AnonReviewer1", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The authors utilise DNA spatial structures called Topologically-associative domains (TADs) from Hi-C data (of the Drosophila fly) and epigenetic marks (such as binding factors) from Chiq-seq data, thereby making use of physical proximity, to predict DNA folding patterns. The authors use a bidirectional LSTM RNN further emphasising that memory of the DNA states contributes to chromatin folding structures.\n\nThe paper offers for a good read. \n\nBelow are comments for improvement and clarification.\n\na) There is only one equation in the paper and this is also not given clearly. What is the K being summed over? \n\nb) wMSE is an old concept and depending on the objective function, the equation can vary. Therefore, change the sentence to read that the authors have a ‘modified' wMSE instead. \n\nc) Section 3.2, Page 3 last line: what is [5]?\n\nd) Section 4.3: 3rd paragraph, last line. Consider formalising the sentence, it will make it more readable. \n\ne) Figure 5: Right panel (top and bottom): correct the x-axis label\n\nf) There is no discussion explaining why regression is better or worse than a neural network, in this application setting. \n\ng) Figure 6 and associated experiment is very interesting and important. The authors should elaborate why, in certain cases, there were huge negative errors in training whereas the test error was positive.", "extra_scores": {"experience_assessment": "I have read many papers in this area.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper at least twice and used my best judgement in assessing the paper.", "review_assessment:_checking_correctness_of_experiments": "I assessed the sensibility of the experiments.", "review_assessment:_checking_correctness_of_derivations_and_theory": "N/A"}, "extra_text": {}}, {"review_id": "rJe0SjlScS", "reviewer": "AnonReviewer3", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The paper predicts DNA folding using different machine learning methods. The authors show that LSTM out performs other methods. They attribute this to the memory of sequential DNA and LSTM model structure. They also propose a weighted mean squared error that improves the performance of the proposed model.\n\nThe authors compare the LSTM model with other classical approaches showing better performance based on predictive and quality metrics, applied to Hi-C data for drosophila, for predicting TADs.\n\nMy major concern is that it is not clear if the improvement is a by-product of LSTM without the proposed new metric. A fair comparison would also consider similar loss function designs for other approaches or at least comparing to a vanilla LSTM model. \n\nAlso the approach lacks illustration of generalizability. The definition of the loss function is also very specific (why 11?) and I wonder if this is generalizable to other Hi-C datasets or predictions based on other epigenetic features beyond ChIP-seq, e.g. ATAC-seq. \n", "extra_scores": {"experience_assessment": "I have published one or two papers in this area.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper at least twice and used my best judgement in assessing the paper.", "review_assessment:_checking_correctness_of_experiments": "I assessed the sensibility of the experiments.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I assessed the sensibility of the derivations and theory."}, "extra_text": {}}, {"review_id": "HylmAQEc_B", "reviewer": "AnonReviewer2", "summary": "", "questions": "", "limitations": "", "rating": 1, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The authors compared a bidirectional LSTM to traditional machine learning models for genomic contact prediction. Although the paper presents an interesting application of ML, I vote for rejection since i)  the paper is similar to previously published works and lacks methodological novelty, ii) the description about the data and methods is not written clearly enough, and iii) the evaluation needs to be strengthened by additional baseline models and evaluation metrics.\n\nMajor comments\n=============\n1. As pointed out by Maria Anna Rapsomaniki, the paper is similar to previously published papers, which are not cited.\n\n2. It is unclear which features were used as inputs. How were the features that are described in section 3.1 represented? Are these binary values or the mean of the Hi-C signal over the genomic region? What does 20kb mean? How were genomic segments of different chromosomes treated?\n\n3. The motivation of using a biLSTM in section 4.3 is unclear. What has the choice of a biLSTM to do with the fact that DNA is double stranded? The DNA is not used as input to the model and genomic contacts can be formed between non-adjacent segments. Please compare recurrent models to non-recurrent models such as fully connected or convolutional networks.\n\n4. The authors used the weighted MSE as evaluation metric, which was used as training loss of the biLSTM, but it is unclear if the same loss was used for linear and gradient boosting regression. To understand if the performance differences are due to the loss function or model architecture, the authors should use the same loss function for training all models, and use additional evaluation metrics such as the MSE and R^2 score.\n\n5. The authors used linear regression weights to quantify feature importance (figure 4). It is unclear if the biLSTM assigned the same importance to features. To quantify the importance of features learned by the biLSTM, the authors could consider correlating the activations of LSTM units with input features, and also analyze if importances depend on the position of features.\n\n6. Which hyper-parameters of the biLSTM and baseline method did the authors tune and how?\n\n7. The authors should compare the described weighted MSE to the mean squared logarithmic error loss, which is commonly used for penalizing large errors less.\n\n8. Why did the authors use the activation of the center LSTM hidden state instead of concatenating the last hidden state of the forward and reverse LSTM? By using the center hidden state, half of the forward and reverse LSTM activations are ignored.", "extra_scores": {"experience_assessment": "I have published in this field for several years.", "review_assessment:_checking_correctness_of_derivations_and_theory": "N/A", "review_assessment:_checking_correctness_of_experiments": "I assessed the sensibility of the experiments.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper thoroughly."}, "extra_text": {}}]

Author rebuttal:

Meta-review: The authors consider the problem of predicting DNA folding patterns.                                                               
They use a range of simple, linear models and find that a bi-LSTM architecture                                                     
yielded best performance.                                                                                                          
                                                                                                                                   
This paper is below acceptance.                                                                                                    
Reviewers pointed out strong similarity to previously published work.                                                              
Furthermore the manuscript lacked in clarity, leaving uncertain eg details about                                                   
experimental details. 
Decision: Reject