Title: AN IMAGE REPRESENTATION BASED CONVOLUTIONAL NETWORK FOR DNA CLASSIFICATION

Abstract: The folding structure of the DNA molecule combined with helper molecules, also referred to as the chromatin, is highly relevant for the functional properties of DNA. The chromatin structure is largely determined by the underlying primary DNA sequence, though the interaction is not yet fully understood. In this paper we develop a convolutional neural network that takes an image-representation of primary DNA sequence as its input, and predicts key determinants of chromatin structure. The method is developed such that it is capable of detecting interactions between distal elements in the DNA sequence, which are known to be highly relevant. Our experiments show that the method outperforms several existing methods both in terms of prediction accuracy and training time.

1. INTRODUCTION:

DNA is perceived as a sequence over the letters {A,C,G,T }, the alphabet of nucleotides. This sequence constitutes the code that acts as a blueprint for all processes taking place in a cell. But beyond merely reflecting primary sequence, DNA is a molecule, which implies that DNA assumes spatial structure and shape. The spatial organization of DNA is achieved by integrating ("recruiting") other molecules, the histone proteins, that help to assume the correct spatial configuration. The combination of DNA and helper molecules is called chromatin; the spatial configuration of the chromatin, finally, defines the functional properties of local areas of the DNA (de Graaf & van Steensel, 2013) .

Chromatin can assume several function-defining epigenetic states, where states vary along the genome (Ernst et al., 2011) . The key determinant for spatial configuration is the underlying primary DNA sequence: sequential patterns are responsible for recruiting histone proteins and their chemical modifications, which in turn give rise to or even define the chromatin states. The exact configuration of the chromatin and its interplay with the underlying raw DNA sequence are under active research. Despite many enlightening recent findings (e.g. Brueckner et al., 2016; The EN-CODE Project Consortium, 2012; Ernst & Kellis, 2013) , comprehensive understanding has not yet been reached. Methods that predict chromatin related states from primary DNA sequence are thus of utmost interest. In machine learning, many prediction methods are available, of which deep neural networks have recently been shown to be promising in many applications (LeCun et al., 2015) . Also in biology deep neural networks have been shown to be valuable (see Angermueller et al. (2016) for a review).

Although DNA is primarily viewed as a sequence, treating genome sequence data as just a sequence neglects its inherent and biologically relevant spatial configuration and the resulting interaction between distal sequence elements. We hypothesize that a deep neural network designed to account for long-term interactions can improve performance. Additionally, the molecular spatial configuration of DNA suggests the relevance of a higher-dimensional spatial representation of DNA. However, due to the lack of comprehensive understanding with respect to the structure of the chromatin, sensible suggestions for such higher-dimensional representations of DNA do not exist.

One way to enable a neural net to identify long-term interactions is the use of fully connected layers. However, when the number of input nodes to the fully connected layer is large, this comes with a large number of parameters. We therefore use three other techniques to detect long-term interactions. First, most convolutional neural networks (CNNs) use small convolution filters. Using larger filters already at an early stage in the network allows for early detection of long-term interactions without the need of fully connected layers with a large input. Second, a deep network similar to the ResNet (He et al., 2015) or Inception (Szegedy et al., 2015) network design prevents features found in early layers from vanishing. Also, they reduce the size of the layers such that the final fully connected layers have a smaller input and don't require a huge number of parameters. Third, we propose a novel kind of DNA representation by mapping DNA sequences to higher-dimensional images using space-filling curves. Space-filling curves map a 1-dimensional line to a 2-dimensional space by mapping each element of the sequence to a pixel in the 2D image. By doing so, proximal elements of the sequence will stay in close proximity to one another, while the distance between distal elements is reduced.

The space-filling curve that will be used in this work is the Hilbert curve which has several advantages. (i): [Continuity] Hilbert curves optimally ensure that the pixels representing two sequence elements that are close within the sequence are also close within the image (Bader, 2016; Aftosmis et al., 2004) . (ii): [Clustering property] Cutting out rectangular subsets of pixels (which is what convolutional filters do) yields a minimum amount of disconnected subsequences (Moon et al., 2001) . (iii): If a rectangular subimage cuts out two subsequences that are disconnected in the original sequence, chances are maximal that the two different subsequences are relatively far apart (see our analysis in Appendix A).

The combination of these points arguably renders Hilbert curves an interesting choice for representing DNA sequence as two-dimensional images. (i) is a basic requirement for mapping short-term sequential relationships, which are ubiquitous in DNA (such as codons, motifs or intron-exon structure). (ii) relates to the structure of the chromatin, which -without all details being fully understood -is tightly packaged and organized in general. Results from Elgin (2012) indicate that when arranging DNA sequence based on Hilbert curves, contiguous areas belonging to identical chromatin states cover rectangular areas. In particular, the combination of (i) and (ii) motivate the application of convolutional layers on Hilbert curves derived from DNA sequence: rectangular subspaces, in other words, submatrices encoding the convolution operations, contain a minimum amount of disconnected pieces of DNA. (iii) finally is beneficial insofar as long-term interactions affecting DNA can also be mapped. This in particular applies to so-called enhancers and silencers, which exert positive (enhancer) or negative (silencer) effects on the activity of regions harboring genes, even though they may be far apart from those regions in terms of sequential distance.

1.1. RELATED WORK:

Since Watson and Crick first discovered the double-helix model of DNA structure in 1953 (Watson & Crick, 1953) , researchers have attempted to interpret biological characteristics from DNA. DNA sequence classification is the task of determining whether a sequence S belongs to an existing class C, and this is one of the fundamental tasks in bio-informatics research for biometric data analysis (Z. Xing, 2010). Many methods have been used, ranging from statistical learning (Vapnik, 1998) to machine learning methods (Michalski et al., 2013) . Deep neural networks (LeCun et al., 2015) form the most recent class of methods used for DNA sequence classification (R.R. Bhat, 2016; Salimans et al., 2016; Zhou & Troyanskaya, 2015; Angermueller et al., 2016) .

Both Pahm et al. (2005) and Higashihara et al. (2008) use support vector machines (SVM) to predict chromatin state from DNA sequence features. While Pahm et al. (2005) use the entire set of features as input to the SVM, Higashihara et al. (2008) use random forests to pre-select a subset of features that are expected to be highly relevant for prediction of chromatin state to use as input to the SVM. Only Nguyen et al. (2016) use a CCN as we do. There are two major differences between their approach and ours. First and foremost, the model architecture is different: the network in Nguyen et al. (2016) consists of two convolution layers followed by pooling layers, a fully connected layer and a sigmoid layer, while our model architecture is deeper, uses residual connections to reuse the learned features, has larger convolution filters and has small layers preceding the fully connected layers (see Methods). Second, while we use a space-filling curve to transform the sequence data into an image-like tensor, Nguyen et al. (2016) keep the sequential form of the input data.

Apart from Elgin (2012) , the only example we are aware of where Hilbert curves were used to map DNA sequence into two-dimensional space is from Anders (2009) , who demonstrated the power of Hilbert curves for visualizing DNA. Beyond our theoretical considerations, these last two studies suggest there are practical benefits of mapping DNA using Hilbert curves.

1.2. CONTRIBUTION:

Our contributions are twofold. First, we predict chromatin state using a CNN that, in terms of architecture, resembles conventional CNNs for image classification and is designed for detecting distal relations. Second, we propose a method to transform DNA sequence patches into two-dimensional image-like arrays to enhance the strengths of CNNs using space-filling curves, in particular the Hilbert curve. Our experiments demonstrate the benefits of our approach: the developed CNN decisively outperforms all existing approaches for predicting the chromatin state in terms of prediction performance measures as well as runtime, an improvement which is further enhanced by the convolution of DNA sequence to a 2D image. In summary, we present a novel, powerful way to harness the power of CNNs in image classification for predicting biologically relevant features from primary DNA sequence.

2. METHODS:



2.1. DNA SEQUENCE REPRESENTATION:

We transform DNA sequences into images through three steps. First, we represent a sequence as a list of k-mers. Next, we transform each k-mer into a one-hot vector, which results in the sequence being represented as a list of one-hot vectors. Finally, we create an image-like tensor by assigning each element of the list of k-mers to a pixel in the image using Hilbert curves. Each of the steps is explained in further detail below.

From a molecular biology point of view, the nucleotides that constitute a DNA sequence do not mean much individually. Instead, nucleotide motifs play a key role in protein synthesis. In bioinformatics it is common to consider a sequence's k-mers, defined as the k-letter words from the alphabet {A,C,G,T} that together make up the sequence. In computer science the term q-gram is more frequently used, and is often applied in text mining (Tomovic et al., 2006) . As an example, the sequence TGACGAC can be transformed into the list of 3-mers {TGA, GAC, ACG, CGA, GAC} (note that these are overlapping). The first step in our approach is thus to transform the DNA sequence into a list of k-mers. Previous work has shown that 3-mers and 4-mers are useful for predicting epigenetic state (Pahm et al., 2005; Higashihara et al., 2008) . Through preliminary experiments, we found that k = 4 yields the best performance: lower values for k result in reduced accuracy, while higher values yield a high risk of overfitting. Only for the Splice dataset (see experiments) we used k = 1 to prevent overfitting, as this is a small dataset.

In natural language processing, it is common to use word embeddings as GloVe or word2vec or onehot vectors (Goldberg, 2017) . The latter approach is most suitable for our method. Each element of such a vector corresponds to a word, and a vector of length N can thus be used to represent N different words. A one-hot vector has a one in the position that corresponds to the word the position is representing, and a zero in all other positions. In order to represent all k-mers in a DNA sequence, we need a vector of length 4 k , as this is the number of words of length k that can be constructed from the alphabet {A,C,G,T}. For example, if we wish to represent all 1-mers, we can do so using a one-hot vector of length 4, where A corresponds to [1 0 0 0], C to [0 1 0 0], G to [0 0 1 0] and T to [0 0 0 1]. In our case, the DNA sequence is represented as a list of 4-mers, which can be converted to a list of one-hot vectors each of length 4 4 = 256.

Our next step is to transform the list of one-hot vectors into an image. For this purpose, we aim to assign each one-hot vector to a pixel. This gives us a 3-dimensional tensor, which is similar in shape to the tensor that serves as an input to image classification networks: the color of a pixel in an RGB-colored image is represented by a vector of length 3, while in our approach each pixel is represented by a one-hot vector of length 256.

What remains now is to assign each of the one-hot vectors in the list to a pixel in the image. For this purpose, we can make use of space-filling curves, as they can map 1-dimensional sequences to a 2-dimensional surface preserving continuity of the sequence (Bader, 2016; Aftosmis et al., 2004) . Various types of space-filling curves are available. We have compared the performance of several such curves, and concluded that Hilbert curves yield the best performance (Appendix A). This corresponds with our intuition: the Hilbert curve has several properties that are advantageous in the case of DNA sequences, as discussed in the introduction section.

The Hilbert curve is a well-known space-filling curve that is constructed in a recursive manner: in the first iteration, the curve is divided into four parts, which are mapped to the four quadrants of a square (Fig. 1a ). In the next iteration, each quadrant is divided into four sub-quadrants, which, in a similar way, each hold 1/16 th of the curve (Fig. 1b ). The quadrants of these sub-quadrants each hold 1/64 th of the curve, etc (Figs. 1c and d ).

By construction, the Hilbert curve yields a square image of size 2 n × 2 n , where n is the order of the curve (see Fig. 1 ). However, a DNA sequence does not necessarily have 2 n * 2 n k-mers. In order to fit all k-mers into the image, we need to choose n such that 2 n * 2 n is at least the number of k-mers in the sequence, and since we do not wish to make the image too large, we pick the smallest such n. In many cases, a large fraction of the pixels then remains unused, as there are fewer k-mers than pixels in the image. By construction, the used pixels are located in upper half of the image. Cropping the picture by removing the unused part of the image yields rectangular images, and increases the fraction of the image that is used (Figure 1e ).

In most of our experiments we used sequences with a length of 500 base pairs, which we convert to a sequence of 500 -4 + 1 = 497 4-mers. We thus need a Hilbert curve of order 5, resulting in an image of dimensions 2 5 × 2 5 × 256 = 32 × 32 × 256 (recall that each pixel is assigned a one-hot vector of length 256). Almost half of the resulting 1024 pixels are filled, leaving the other half of the image empty which requires memory. We therefore remove the empty half of the image and end up with an image of size 16 × 32 × 256.

The data now has the appropriate form to input in our model. 

2.2. NETWORK ARCHITECTURE:

Modern CNNs or other image classification systems mainly focus on gray-scale images and standard RGB images, resulting in channels of length 1 or 3, respectively, for each pixel. In our approach, each pixel in the generated image is assigned a one-hot vector representing a k-mer. For increasing k, the length of the vector and thus the image dimension increases. Here, we use k = 4 resulting in 256 channels, which implies that each channel contains very sparse information. Due to the curse of dimensionality standard network architectures applied to such images are prone to severe overfitting.

Here, we design a specific CNN for the kind of high dimensional image that is generated from a DNA sequence. The architecture is inspired by ResNet (He et al., 2015) and Inception (Szegedy et al., 2015) . The network has L layers and each layer implements a non-linear function F l (x l ) where l is the index of the hidden layer with output x l+1 . The function F l (x l ) consists of various layers such as convolution (denoted by c), batch normalization (bn), pooling (p) and non-linear activation function (af ).

The first part of the network has the objective to reduce the sparseness of the input image (Figure 2 ), and consists of the consecutive layers [c, bn, c, bn, af, p] . The main body of the network enhances the ability of the CNN to extract the relevant features from DNA space-filling curves. For this purpose, we designed a specific Computational Block inspired by the ResNet residual blocks (He et al., 2015) . The last part of the network consists of 3 fully-connected layers, and softmax is used to obtain the output classification label. The complete model is presented in Table 1 , and code is available on Github (https://github.com/Bojian/Hilbert-CNN/tree/master). A simplified version of our network with two Computational Blocks is illustrated in Figure 2 .

Computation Block. In the Computation Block first the outputs of two Residual Blocks and one identity mapping are summed, followed by a bn and an af layer (Figure 2 ). In total, the computational block has 4 convolutional layers, two in each Residual Block (see Figure 3 ). The Residual Block first computes the composite function of five consecutive layers, namely [c, bn, af, c, bn], followed by the concatenation of the output with the input tensor. The residual block concludes with an af .

The Residual Block can be viewed as a new kind of non-linear layer denoted by:

Residual l [k j ,k j+1 ,d link ,d out ],

where k j and k j+1 are the respective filter sizes of the two convolutional layers. d link and d out are the dimensions of the outputs of the first convolutional layer and the Residual Block, respectively, where d link < d out ; this condition simplifies the network architecture and reduces the computational cost. The Computational Block can be denoted as C[k1,k2,k3,k4] with the two residual blocks defined as Residual

1 [k 1 ,k 2 ,d link ,d out ] and Residual 2 [k 3 ,k 4 ,d link ,d out ].

Note that here we chose the same d link and d out for both Residual Blocks in a Computational Block.

Implementation details. Most convolutional layers use small squared filters of size 2, 3 and 4, except for the layers in the first part of the network, where large filters are applied to capture long range features. We use Exponential Linear Units (ELU, Clevert et al. (2015) ) as our activation function af to reduce the effect of gradient vanishing: preliminary experiments showed that ELU preformed significantly better than other activation functions (data not shown). For the pooling layers p we used Average pooling. Average pooling outperformed Max pooling in terms of prediction accuracy by more than 2% in general, as it reduces the high variance of the sparse generated images. Cross entropy was used as the loss function. 

3. EXPERIMENTS:

We test the performance of our approach using ten publicly available datasets from Pokholok et al. (2005) . The datasets contain DNA subsequences with a length of 500 base pairs. Each sequence is labeled either as "positive" or "negative", indicating whether or not the subsequence contains regions that are wrapped around a histone protein. The ten datasets each consider a different type of histone protein, indicated by the name of the dataset. Details can be found in Table 2 .

A randomly chosen 90% of the dataset is used for training the network, 5% is used for validation and early stopping, and the remaining (5%) is used for evaluation. We train the network using the AdamOptimizer (Kingma & Ba, 2017) 1 . The learning rate is set to 0.003, the batch-size was set to 300 samples and the maximum number of epochs is 10. After each epoch the level of generalization is measured as the accuracy obtained on the validation set. We use early stopping to prevent overfitting. To ensure the model stops at the correct time, we combine the GL α measurement (Prechelt, 1998) of generalization capability and the No-Improvement-In-N-Steps(Nii-N) method (Prechelt, 1998) . For instance, Nii-2 means that the training process is terminated when generalization capability is not improved in two consecutive epochs.

We compare the performance of our approach, referred to as HCNN, to existing models. One of these is the support vector machine (SVM) model by Higashihara et al. (2008) , for which results are available in their paper. Second, in tight communication with the authors, we reconstructed the Seq-CNN model presented in Nguyen et al. (2016) (the original software was no longer available), see Appendix C for detailed settings. Third, we constructed the commonly used LSTM, where the so-called of the sequence is used as input. A 4-mer profile is a list containing the number of occurrences of all 256 4-mers of the alphabet {A,C,G,T} in a sequence. Preliminary tests showed that using all 256 4-mers resulted in overfitting, and including only the 100 most frequent 4-mers is sufficient. Details of the LSTM architecture can be found in Table 8 in Appendix C.

In order to assess the effect of using a 2D representation of the DNA sequence in isolation, we compare HCNN to a neural network using a sequential representation as input. We refer to this model as seq-HCNN. As in HCNN, the DNA sequence is converted into a list of kmer representing one-hot vectors, though the mapping of the sequence into a 2D image is omitted. The network architecture is a "flattened" version of the one used in HCNN: for example, a 7×7 convolution filter in HCNN is transformed to a 49×1 convolution filter in the 1D-sequence model. As a summary of model size, the Seq-CNN model contains 1.1M parameters, while both HCNN and seq-HCNN have 961K parameters, and the LSTM has 455K parameters.

In order to test whether our method is also applicable to DNA sequence classification tasks other than chromatin state prediction only, we performed additional tests on the splice-junction gene sequences dataset from Lichman (2013) . Most of the DNA sequence is unused, and splice-junctions refer to positions in the genetic sequence where the transition from an unused subsequence (intron) to a used subsequence (exon) or vice versa takes place. The dataset consists of DNA subsequences of length 61, and each of the sequences is known to be an intron-to-exon splice-junction, an exon-to-intron splice junction or neither. As the dataset is relatively small, we used 1-mers instead of 4-mers.

Note that the sequences are much shorter than for the other datasets, resulting in smaller images (dimensions 8 × 8 × 4).

4. RESULTS:

The results show that SVM and Seq-CNN were both outperformed by HCNN and seq-HCNN; LSTM shows poor performance. HCNN and seq-HCNN show similar performance in terms of prediction accuracy, though HCNN shows more consistent results over the ten folds indicating that using a 2D representation of the sequence improves robustness. Furthermore, HCNN yields better performance than seq-HCNN in terms of precision, recall, AP and AUC (Table 5 ). It thus enables to reliably vary the tradeoff between recall and false discoveries. HCNN outperforms all methods in terms of training time (Table 4 ).

The good performance of HCNN observed above may either be attributable to the conversion from DNA sequence to image, or to the use of the Hilbert curve. In order to address this question, we adapted our approach by replacing the Hilbert curve with other space-filling curves and compared their prediction accuracy. Besides the Hilbert curve, other space-filling curves exist (Moon et al., 2001) (see Appendix A ). In Figure 4 , we compare the performance of our model with different mapping strategies in various datasets as displayed. We find that the images generated by the spacefilling Hilbert curve provide the best accuracy on most datasets and the 1-d sequence performs worst.

5. DISCUSSION:

In this paper we developed a CNN that outperforms the state-of-the-art for prediction of epigenetic states from primary DNA sequence. Indeed, our methods show improved prediction accuracy and training time compared to the currently available chromatin state prediction methods from Pahm In HCNN on the other hand the number of nodes is strongly reduced before introducing a fully connected layer. Third, the use of a two-dimensional input further enhances the model's capabilities of incorporating long-term interactions.

We showed that seq-HCNN and HCNN are not only capable of predicting chromatin state, but can also predict the presence or absence of splice-junctions in DNA subsequences. This suggests that our approach could be useful for DNA sequence classification problems in general.

Hilbert curves have several properties that are desirable for DNA sequence classification. The intuitive motivation for the use of Hilbert curves is supported by good results when comparing Hilbert curves to other space-filling curves. Additionally, Hilbert curves have previously been shown to be useful for visualization of DNA sequences (Anders, 2009) .

The main limitation of Hilbert curves is their fixed length, which implies that the generated image contains some empty spaces. These spaces consume computation resources; nevertheless, the 2D representation still yields reduced training times compared to the 1D-sequence representation, presumably due to the high degree of optimization for 2D inputs present in standard CNN frameworks.

Given that a substantial part of the improvements in performance rates are due to our novel architecture, we plan on investigating the details of how components of the architecture are intertwined with improvements in prediction performance in more detail. We also plan to further investigate why Hilbert curves yield the particular advantages in terms of robustness and false discovery control we have observed here.

A COMPARISON OF SPACE-FILLING CURVES WITH REGARD TO LONG-TERM

INTERACTIONS:

As noted before, long-term interactions are highly relevant in DNA sequences. In this section consider these long-term interactions in four space-filling curves: the reshape curve, the snake curve, the Hilbert curve and the diag-snake curve. See Fig. 5 for an illustration. As can be seen in Fig. 5 , mapping a sequence to an image reduces the distance between two elements that are far from one another in the sequence, while the distance between nearby elements does not increase. Each of the curves does have a different effect on the distance between faraway elements. In order to assess these differences, we use a measure that is based on the distance between two sequence elements as can be observed in the image. We denote this distance by • L reshape (A,P ) = 3 √ 2 for the reshape curve;

• L snakesnake (A,P ) = 3 for the snake curve;

• L diag-snake (A,P ) = 3 √ 2 for the diagonal snake curve.

• L hilbert (A,P ) = 3 for the Hilbert curve;

We now introduce the following measure:

Γ(C) = mean(∆(C)) max(∆(C)) ,

where the ∆(C) is the set of the weighted distances between all pairs of the elements in the sequence.

Here, ∆(C) is a set containing the distance between any two sequence elements, weighted by their distance in the sequence:

∆(C) = {L seq (x,y) • L C (x,y) | x, y ∈ C, x = y}.

Note that a low max(∆(C)) relative to mean(∆(C)) implies that long-term interactions are strongly accounted for, so a high Γ(C) is desirable.

Γ(C) is evaluated for the four space-filling curves as well as the sequence representation for sequences of varying lengths.The results show that the Hilbert curve yields the highest values for Γ(C) 6 and thus performs best in terms of retaining long-distance interactions.

:

The LSTM model was implemented in Keras(Chollet et al.,

2015), all other models were implemented in Tensorflow(Abadi et al., 2015).

Figures and Tables:

Figure 1: Hilbert Curve and cropped image

Figure2: A simplified version of HCNN with two computational blocks. BN is short for "Batch Normalization".

Figure 3: Residual Block

Figure 4: HCNN with different mapping strategies

Figure 5: Space-filling curves

L C (x,y) where x, y ∈ S, with S the sequence and C the curve under consideration. Then for the sequence S = {A,B,C,D,E, • • • ,P } we obtain • L seq (A,P ) = 15 for the sequence;

Network Architecture. The output size is given as height, width, channels.

Table content: <table><tr><td>Layers</td><td>Description</td><td>Output size</td></tr><tr><td>Input</td><td/><td>16,32, 256</td></tr><tr><td>Conv1</td><td>7 × 7 conv, BN</td><td>16,32, 64</td></tr><tr><td>Conv2</td><td>5 × 5 conv, BN</td><td>16,32, 64</td></tr><tr><td>Activation</td><td/><td/></tr><tr><td>Average-Pool</td><td colspan="2">2 × 2 filter, stride 2 8, 16, 64</td></tr><tr><td colspan="3">ComputationBlock [8,4,4,3], d link = 4 8, 16, 32</td></tr><tr><td colspan="3">ComputationBlock [3,3,3,3], d link = 4 8, 16, 32</td></tr><tr><td>Average-Pool</td><td colspan="2">2 × 2 filter, stride 2 4, 8, 32</td></tr><tr><td colspan="3">ComputationBlock [2,4,4,3], d link = 4 4, 8, 32</td></tr><tr><td colspan="3">ComputationBlock [2,2,2,2], d link = 4 4, 8, 32</td></tr><tr><td colspan="3">ComputationBlock [3,2,2,3], d link = 4 4, 8, 32</td></tr><tr><td>Average-Pooling</td><td colspan="2">2 × 2 filter, stride 2 2, 4, 32</td></tr><tr><td>BN,Activation</td><td/><td/></tr><tr><td>Average-Pool</td><td colspan="2">2 × 2 filter, stride 2 1, 2, 32</td></tr><tr><td>Classification layer1</td><td>64D FC layer activation dropout with 0.5</td><td/></tr><tr><td>Classification layer2</td><td>100D FC layer activation dropout with 0.5</td><td/></tr><tr><td>Classification</td><td>50D FC layer,</td><td/></tr><tr><td>layer3</td><td>softmax</td><td/></tr></table>

DNA Datasets

Table content: <table><tr><td>Name</td><td colspan="2">Samples Description</td></tr><tr><td>H3</td><td>14965</td><td>H3 occupancy</td></tr><tr><td>H4</td><td>14601</td><td>H4 occupancy</td></tr><tr><td>H3K9ac</td><td>27782</td><td>H3K9 acetylation relative to H3</td></tr><tr><td>H3K14ac</td><td>33048</td><td>H3K14 acetylation relative to H3</td></tr><tr><td>H4ac</td><td>34095</td><td>H4 acetylation relative to H3</td></tr><tr><td>H3K4me1</td><td>31677</td><td>H3K4 monomethy-lation relative to H3</td></tr><tr><td>H3K4me2</td><td>30683</td><td>H3K4 dimethylation relative to H3</td></tr><tr><td>H3K4me3</td><td>36799</td><td>H3K4 trimethylation relative to H3</td></tr><tr><td colspan="2">H3K36me3 34880</td><td>H3K36 trimethylation relative to H3</td></tr><tr><td colspan="2">H3K79me3 28837</td><td>H3K79 trimethylation relative to H3</td></tr><tr><td>Splice</td><td>3190</td><td>Splice-junction Gene Sequences</td></tr></table>

Prediction accuracy obtained with an SVM-based method, Seq-CNN fromNguyen et al. (2016), LSTM, seq-HCNN and HCNN. The results for SVM are taken from Table12inHigashihara et al. (2008). In the splice dataset, Seq-CNN performed best when using 4-mers, while for HCNN and seq-HCNN 1-mers yielded the best performance.

Table content: <table><tr><td>Dataset</td><td>SVM</td><td colspan="2">LSTM Seq-CNN</td><td>seq-HCNN</td><td>HCNN</td></tr><tr><td colspan="3">H3 86.47% 64.13%</td><td colspan="3">79.25% 86.86 ± 1.563% 87.34 ±0.263%</td></tr><tr><td colspan="3">H4 87.82% 63.82%</td><td colspan="3">81.86% 87.31 ± 0.952% 87.33±0.264%</td></tr><tr><td colspan="3">H3K9ac 75.08% 63.07%</td><td colspan="3">68.76% 78.47 ± 0.699% 79.19±0.239%</td></tr><tr><td colspan="3">H3K14ac 73.28% 68.31%</td><td colspan="3">68.31% 75.06 ±0.987% 74.79±0.226%</td></tr><tr><td colspan="3">H4ac 72.06% 60.63%</td><td colspan="3">64.80% 77.04 ± 1.256% 77.06±0.233%</td></tr><tr><td colspan="3">H3K4me1 69.71% 60.43%</td><td colspan="3">62.60% 73.47 ±0.789% 73.21±0.221%</td></tr><tr><td colspan="3">H3K4me2 68.97% 61.45%</td><td colspan="3">62.38% 73.91 ± 0.631% 74.27±0.224%</td></tr><tr><td colspan="3">H3K4me3 68.57% 58.03%</td><td colspan="3">62.33% 74.54 ±0.865% 74.45±0.225%</td></tr><tr><td colspan="3">H3K36me1 75.19% 60.78%</td><td colspan="3">72.20% 77.18 ±0.973% 77.03±0.232%</td></tr><tr><td colspan="3">H3K79me1 80.58% 63.84%</td><td colspan="3">75.07% 81.66 ±1.264% 81.63±0.246%</td></tr><tr><td colspan="3">Splice 94.70% 96.23%</td><td colspan="3">91.82% 93.21 ± 1.645% 94.11±0.284%</td></tr></table>

Training times, presented as min:sec.We believe that the improved performance over the CNN developed byNguyen et al. (2016) (Seq-CNN) is a result of three factors. First, our network uses larger convolutional filters, allowing the model to detect long-distance interactions. Second, despite HCNN being deeper, both HCNN and seq-HCNN have a smaller number of parameters, allowing for faster optimization. This is due to the size of the layer preceding the fully connected layer, which is large in the method proposed

Table content: <table><tr><td>Dataset</td><td colspan="4">LSTM seq-CNN seq-HCNN HCNN</td></tr><tr><td>H3</td><td>35:43</td><td>95:23</td><td>6:47</td><td>3:40</td></tr><tr><td>H4</td><td>45:32</td><td>95:53</td><td>5:12</td><td>3:12</td></tr><tr><td>H3K9ac</td><td>76:06</td><td>173:18</td><td>17:24</td><td>7:40</td></tr><tr><td>H3K14ac</td><td>81:21</td><td>180:56</td><td>17:42</td><td>13:24</td></tr><tr><td>H4ac</td><td>93:32</td><td>181:33</td><td>24:48</td><td>17:32</td></tr><tr><td>H3K4me1</td><td>93:44</td><td>192:20</td><td>18:30</td><td>10:38</td></tr><tr><td>H3K4me2</td><td>94:22</td><td>188:13</td><td>18:23</td><td>14:38</td></tr><tr><td>H3K4me3</td><td>96:03</td><td>162:32</td><td>20:40</td><td>11:33</td></tr><tr><td colspan="2">H3K36me3 93:48</td><td>161:12</td><td>21:52</td><td>16:37</td></tr><tr><td colspan="2">H3K79me3 64:28</td><td>158:34</td><td>14:25</td><td>10:13</td></tr><tr><td>Splice</td><td>6:42</td><td>35:12</td><td>3:42</td><td>1:30</td></tr><tr><td colspan="5">et al. (2005), Higashihara et al. (2008) and Nguyen et al. (2016) as well as an LSTM model. Addi-</td></tr><tr><td colspan="5">tionally, we showed that representing DNA-sequences with 2D images using Hilbert curves further</td></tr><tr><td colspan="5">improves precision and recall as well as training time as compared to a 1D-sequence representation.</td></tr></table>

Recall, Precision, area under precision-recall curve (AP) and area under ROC curve (AUC) for seq-HCNN and HCNN. The reported values are the means over ten folds.

Table content: <table><tr><td>Dataset</td><td colspan="2">Recall seq-HCNN HCNN</td><td colspan="2">Precision seq-HCNN HCNN</td><td colspan="2">AP seq-HCNN HCNN</td><td>AUC seq-HCNN HCNN</td></tr><tr><td>H3</td><td>85.67%</td><td colspan="2">87.33% 85.67%</td><td colspan="2">87.33% 90.33%</td><td>93.33% 91.00%</td><td>93.67%</td></tr><tr><td>H4</td><td>87.00%</td><td/><td/><td/><td/></tr></table>

Peer review (iclr2018):
 Dear editors, the authors addressed all of my comments and clearly improved their manuscript over multiple iterations. I therefore increased my rating from ‘6: Marginally above acceptance threshold’ to ‘7: Good paper, accept’. Please note that the authors made important edits to their manuscript after the ICLR deadline and could hence not upload their most current version, which you can from https://file.io/WIiEw9. If you decided to publish the manuscript, I hence highly suggest using this (https://file.io/WIiEw9) version. Best, -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- The authors present Hilbert-CNN, a convolutional neural network for DNA sequence classification. Unlike existing methods, their model does not use the raw one-dimensional (1D) DNA sequence as input, but two-dimensional (2D) images obtained by mapping sequences to images using spacing-filling Hilbert-Curves. They further present a model (Hilbert-CNN) that is explicitly designed for Hilbert-transformed DNA sequences. The authors show that their approach can increase classification accuracy and decrease training time when applied to predicting histone-modification marks and splice junctions. Major comments ============= 1. The motivation of transforming sequences into images is unclear and claimed benefits are not sufficiently supported by experiments. The essence of deep neural networks is to learn a hierarchy of features from the raw data instead of engineering features manually. Using space filling methods such as Hilbert-curves to transform (DNA) sequences into images can be considered as unnecessary feature-engineering. The authors claim that ‘CNNs have proven to be most powerful when operating on multi-dimensional input, such as in image classification’, which is wrong. Sequence-based convolutional and recurrent models have been successfully applied for modeling natural languages (translation, sentiment classification, …), acoustic signals (speech recognition, audio generation), or biological sequences (e.g. predicting various epigenetic marks from DNA as reviewed in Angermueller et al). They further claim that their method can ‘better take the spatial features of DNA sequences into account’ and can better model ‘long-term interactions’ between distant regions. This is not obvious since Hilbert-curves map adjacent sequence characters to pixels that are close to each other as described by the authors, but distant characters to distant pixels. Hence, 2D CNN must be deep enough for modeling interactions between distant image features, in the same way as a 1D CNN. Transforming sequences to images has several drawbacks. 1) Since the resulting images have a small width and height but many channels, existing 2D CNNs such as ResNet or Inception can not be applied, which also required the authors to design a specific model (Hilbert-CNN). 2) Hilbert-CNN requires more memory due to empty image regions. 3) Due to the high number of channels, convolutional filters have more parameters. 4) The sequence-to-image transformation makes model-interpretability hard, which is in particular important in biology. For example, motifs of the first convolutional layers can not be interpreted as sequence motifs (as described in Angermueller et al) and it is unclear how to analyze the influence of sequence characters using attention or gradient-based methods. The authors should more clearly motivate their model in the introduction, tone-down the benefit of sequence-to-image transformations, and discuss drawbacks of their model. This requires major changes of introduction and discussion. 2. The authors should more clearly describe which and how they optimized hyper-parameters. The authors should optimize the most important hyper-parameters of their model (learning rate, batch size, weight decay, max vs. average pooling, ELU vs. ReLU, …) and baseline models on a holdout validation set. The authors should also report the validation accuracy for different sequence lengths, k-mer sizes, and space filling functions. Can their model be applied to longer sequences (>= 1kbp) which had been shown to improve performance (e.g. 10.1101/gr.200535.115)? Does Figure 4 show the performance on the training, validation, or test set? 3. It is unclear if the performance gain is due the proposed sequence-to-image transformation, or due to the proposed network architecture (Hilbert-CNN). It is also unclear if Hilbert-CNNs are applicable to DNA sequence classification tasks beyond predicting chromatin states and splice junctions. To address these points, the authors should compare Hilbert-CNN to models of the same capacity (number of parameters) and optimize hyper-parameters (k-mer size, convolutional filter size, learning rate, …) in the same way as they did for Hilbert-CNN. The authors should report the number of parameters of all models (Hilbert-CNN, Seq-CNN, 1D-sequence-CNN (Table 5), and LSTM (Table 6), …) in an additional table. The authors should also compare Hilbert-CNN to the DanQ architecture on predicting epigenetic markers using the same dataset as reported in the DanQ publication (DOI: 10.1093/nar/gkw226). The authors should also compare Hilbert-CNNs to gapped-kmer SVM, a shallow model that had been successfully applied for genomic prediction tasks. 4. The authors should report the AUC and area under precision-recall curve (APR) in additional to accuracy (ACC) in Table 3. 5. It is unclear how training time was measured for baseline models (Seq-CNN, LSTM, …). The authors should use the same early stopping criterion as they used for training Hilber-CNNs. The authors should also report the training time of SVM and gkm-SVM (see comment 3) in Table 3. Minor comments ============= 1. The authors should avoid uninformative adjectives and clutter throughout the manuscript, for example ‘DNA is often perceived’, ‘Chromatin can assume’, ‘enlightening’, ‘very’, ‘we first have to realize’, ‘do not mean much individually’, ‘very much like the tensor’, ‘full swing’, ‘in tight communication’, ‘two methods available in the literature’. The authors should point out in section two that k-mers can be overlapping. 2. Section 2.1: One-hot vectors is not the only way for embedding words. The authors should also mention Glove and word2vec. Similar approaches had been applied to protein sequences (DOI: 10.1371/journal.pone.0141287) 3. The authors should more clearly describe how Hilbert-curves map sequences to images and how images are cropped. What does ‘that is constructed in a recursive manner’ mean? Simply cropping the upper half of Figure 1c would lead to two disjoint sequences. What is the order of Figure 1e? 4. The authors should consistently use ‘channels’ instead of ‘full vector of length’ to denote the dimensionality of image pixels. 5. The authors should use ‘Batch norm’ instead of ‘BN’ in Figure 2 for clarification. 6. Hilber-CNN is similar to ResNet (DOI: 10.1371/journal.pone.0141287), which consists of multiple ‘residual blocks’, where each block is a sequence of ‘residual units’. A ‘computational block’ in Hilbert-CNN contains two parallel ‘residual blocks’ (Figure 3) instead of a sequence of ‘residual units’. The authors should use ‘residual block’ instead of ‘computational block’, and ‘residual units’ as in the original ResNet publication. The authors should also motivate why two residual units/blocks are applied parallely instead of sequentially. 7. Caption table 1: the authors should clarify if ‘Output size’ is ‘height, width, channels’, and explain the notation in ‘Description’ (or refer to the text.)