Deep learning and support vector machines for transcription start site identification

Generation of the dataset

We used the reference sequence of the human genome in its GRCh38.p13 version to have a reliable data source for our experiments. We chose this version because it is the most recent one available in Ensembl at the moment. However, the DNA sequence by itself is not enough, the specific TSS position of each transcript is needed. In this section, we explain the steps followed to generate the final dataset (Fig. 1). These steps are: raw data gathering, positive instances processing, negative instances generation, and data splitting by chromosomes.

First, we need an interface to download the raw data, which are made up of every transcript sequence in the human genome. We use Ensembl release 104 (Howe et al., 2020) and its utility BioMart (Smedley et al., 2009), which allows us to easily obtain large amounts of data. It also enables us to select a wide variety of interesting fields, including the transcription start and end sites. After filtering instances that present null values in any relevant field, this combination of the sequence and its flanks will form our raw dataset.

Once the sequences are available, we find the TSS position (given by Ensembl) and the two following bases to treat it as a codon. After that, 700 bases before this codon and 300 bases after it were concatenated, giving the final sequence of 1,003 nucleotides that is going to be used in our models. These specific window values have been used in Bhandari et al. (2021) and we have kept them because we find them interesting for comparison purposes.

One of the most sensitive parts of this dataset is the generation of negative instances. We cannot obtain this kind of data directly, so we need to generate them synthetically. To get examples of negative instances, i.e., sequences that do not represent a transcript start site, we select random DNA positions inside the transcripts that do not correspond to a TSS. Once we have selected the specific position, we obtain 700 bases ahead and 300 bases after it, as we did with the positive instances. Regarding the positive to negative ratio, in a similar problem, but studying TIS instead of TSS (Zhang et al., 2017), a ratio of 10 negative instances to each positive one was found optimal. Following this idea, we select 10 random positions from the transcript sequence of each positive codon and label them as negative instances.

After this process, we end up with 1,122,113 instances: 102,488 positive and 1,019,625 negative sequences. To validate and test our models, we need to split this dataset into three parts: training, validation, and testing. We have decided to make this differentiation by chromosomes, as in Perez-Rodriguez, Haro-Garcia & Garcia-Pedrajas (2020). Thus, we use chromosome 16 as validation because it is a good example of a chromosome with average characteristics. Then we selected samples from chromosomes 1, 3, 13, 19 and 21 to be part of the test set and used the rest of them to train our models.

Every step of this process can be replicated using the scripts available in https://github.com/JoseBarbero/EnsemblTSSPrediction.

Sequence encoding

Most machine learning methods cannot read DNA sequences directly. One of the most popular ways to encode nucleotides strings is one-hot encoding. Using this method, we represent every base using a binary array of four positions, each representing one possible base. Only one of the values of the array will be 1, while the others remain 0. Hence, we have a simple way to represent four different characters: A, C, G, and T. Although some approaches are using word embedding systems based on the natural language processing field (Wu et al., 2021), we found one-hot encoding adequate for the purposes of this work because of its simplicity. This is the encoding used for the deep learning methods. For SVM methods, we have used both this encoding (for SVM with RBF kernel) and directly the raw sequence (for SVM with string kernel).

Support vector machines (SVMs)

SVMs (Cortes & Vapnik, 1995) are one of the most used methods in machine learning and specially in binary classification. We find that SVMs are one of the few alternatives that can keep up with DNNs in complex problems such as TSS identification. These methods have some advantages over DNNs, like the possibility of introducing “biases” based on human expertise in our models through its kernel. Since this is essentially a string classification problem, we need to use a string kernel with our SVM. In Perez-Rodriguez, Haro-Garcia & Garcia-Pedrajas (2020), it was found that the weighted degrees kernel consistently gives the best results in the identification of TIS. Due to the similarity of the problem we are trying to solve, we decided to use a weighted degree kernel as well. This kernel is also the one that gave the best results in Sonnenburg et al. (2007), where they were compared with Markov chains and SVMs using other kernels.

SVM algorithms integrate a mathematical function named kernel. The purpose of this kernel is to implement the trick used by the SVM algorithms that allows instances to be represented as a relationship metric between two data points instead of being explicitly represented in a higher-dimensional space. In the case of string kernels, it is specially useful to get numerical values from the string sequences, which originally could not be used by the models.

As mentioned before, we used the weighted degrees (WD) kernel. This type of kernel allows us to represent the similarity of two sequences by counting the amount of k-mers of size k ∈ {1, …, d} that two specific sequences s_i and s_j of length L share. This function is defined by the Eq. (1), where 𝕀 returns 1 when the k-mer starting in position l is the same in sequences s_i and s_j and 0 if not. β_k is a weight factor defined by ${\beta }_k=2\frac{\left(d-k+1\right)}{d\left(d+1\right)}$. (1)${\displaystyle K\left({\mathbf{s}}_i,{\mathbf{s}}_j\right)=\sum _{k=1}^d{\beta }_k\sum _{l=1}^{L-k+1}\mathbb{I}\left(\right.{\mathbf{u}}_{k,l}\left({\mathbf{s}}_i\right)={\mathbf{u}}_{k,l}\left({\mathbf{s}}_j\right)\left)\right.}$

To have a similar model to evaluate the benefits of a string kernel SVM, we additionally used a numerical method, specifically an SVM with a radial basis function (RBF) kernel. This kernel is commonly used in machine learning with SVM classifiers. Experimenting with this kernel allows us to see the differences in performance between a string kernel and a numerical one. In this case, the kernel function is described by Eq. (2), which is defined by the squared Euclidean distance of two instances x and x′ and a γ value that determines the influence of the support vectors in the model. (2)${\displaystyle K\left(\mathbf{x},{\mathbf{x}}^{\prime }\right)=exp\left(-\gamma \parallel \mathbf{x}-{\mathbf{x}}^{\prime }{\parallel }^2\right).}$

The use of the RBF kernel requires the one-hot method to encode the instances, given that this is a numerical kernel and it is not possible to use the string instances as we did with the weighted degrees kernel.

All the parameters in Eqs. (1) and (2) as well as the regularization parameter in the support vector classifier model were defined by an extensive grid search cross-validation process.

Another advantage of SVMs is their interpretability, as opposed to the black-box nature of DNNs. Given the robust mathematical theory behind the method and the simplicity that characterizes the kernels, the models are much easier to understand for humans than DNNs. This simplicity also makes the model easier to tune, having fewer parameters to optimize, needing less memory to store the trained model and having a lighter and faster predictive model (Yoon et al., 2011). It is also a method without the problems of local minima, being, in principle, less prone to overfitting.

Several works can be found on SVMs applied to DNA sequences that demonstrate a good performance on similar problems. Specifically, great results were obtained in alternative splicing recognition (Rätsch, Sonnenburg & Schölkopf, 2005), promoter (Wang et al., 2018; Sato, 2018; Cassiano & Silva-Rocha, 2020) and promoter features prediction (Meng et al., 2017), transcription start identification in eukaryotes (Sonnenburg, Zien & Rätsch, 2006; Georgakilas, Perdikopanis & Hatzigeorgiou, 2020) and prokaryotes (Towsey, Gordon & Hogan, 2006) or translation initiation sites prediction (Perez-Rodriguez, Haro-Garcia & Garcia-Pedrajas, 2020). SVMs have also been combined with convolutional neural networks (Qian et al., 2018). This approach will be addressed in the next section.

The string kernel has been applied using Scikit-learn’s (Pedregosa et al., 2011) support vector classifier implementation, which allows using custom kernels. However, the kernel needs to be precomputed. This implies creating the gram matrix for the train and test datasets, being a computationally slow process. We tried to minimize this issue by parallelizing the process, although it is still noticeably slower than the deep learning methods that will be presented in the next section. The main parameter to choose is the value d, which represents the maximum length of k-mers to consider. We conducted a preliminary experiment using a reduced dataset, formed by a sample of 1% of the positive and negative instances, in order to test several values of d, the results of which can be seen in Table 1. Higher values generally achieve better results, but at a much higher running time. As we will address later, training times can be significant having to calculate the kernel gram matrix: a problem implying quadratic growth. For that reason, we determined that a d value of 10 is reasonable in this case. It improves the results significantly, being a bit less precise than 20 but at a considerably lower running time.

The most relevant inconvenience of support vectors machines are their memory and execution time requirements when used with very large data sets, precisely the data sets that appear in bioinformatics. For reference, in one of our preliminary experiments we needed to use a server with 1 TB of RAM, and the experiment took several months to complete, even parallelizing using 20 cores. These execution times were reduced in subsequent experiments using an approach similar to that suggested in Graf et al. (2004), the Cascade SVM. Specifically, we have used the implementation of the method available in https://github.com/fhebert/CascadeSVC for the RBF kernel, and a modified version to work with the WD kernel. The idea of this method is to get the support vectors by levels. In the first level, the original dataset is divided into several subsets (typically a number power of two), then an SVM is used on each subset to obtain a set of support vectors. In the next level of the cascade, two sets of support vectors from the previous level are combined and used on a new SVM, which in turn will give a set of support vectors that will be combined with others at the next level of the cascade. In the last level of the cascade, with only one SVM, the final set of support vectors that is obtained is added to each of the subsets of the first level of the cascade, and the process is repeated for a few additional iterations until convergence is reached. The application of this procedure allows reducing the duration of an experiment from several months to several hours. It also allows reducing the memory requirements.

Deep learning

It could be said that the emergence of deep learning was a turning point for artificial intelligence and, recently, also in bioinformatics. DNNs have set a new state-of-the-art performance in tasks like variant calling (Poplin et al., 2018) or protein folding prediction (Jumper et al., 2021), and have been found very valuable in promoter identification (Bhandari et al., 2021). Since the first work using neural networks in promoter prediction (Demeler & Zhou, 1991), models such as multilayer perceptrons (MLP) (Shahmuradov, Umarov & Solovyev, 2017), convolutional neural networks (Umarov et al., 2019; Pachganov et al., 2019; Wei et al., 2021) or a combination of CNN and recurrent neural networks (RNNs) (Oubounyt et al., 2019; Zhang et al., 2017) have appeared in problems related to promoter identification. Meanwhile, it is much harder to find deep learning approaches for TSS identification than other related problems, and the ones available (Mahdi & Rouchka, 2009; Zheng, Li & Hu, 2020), do not offer a comparison that takes into account SVMs as the most established method in this field.

Unlike SVMs, DNNs’ interpretability is a challenging task, but on the other hand, they do not require as much expertise in the field of application as with SVMs kernel design. DNNs can learn the particularities of the problem by themselves. Even though DNNs do not need as much knowledge in the molecular biology field as SVMs do to create appropriate kernels, they do need that expertise in the deep learning field to be able to develop and fine-tune the model. Furthermore, DNNs need significantly more data to be trained in complex problems such as TSS prediction, and can fall in local minima and tend to overfit more than SVMs (Salman & Liu, 2019).

In Khan et al. (2020), in the context of promoter identification, the researchers concluded that MLPs can achieve better results than SVMs. We believe that extending this research to TSS identification could be promising. Furthermore, string kernels were not taken into account, only RBF and sigmoid kernels were tested, and we think that the inclusion of CNNs and RNNs will enrich the study.

We selected three of the most representative types of DNNs that can work with one-hot encoded sequences: a CNN, a LSTM (Hochreiter & Schmidhuber, 1997) and a combination of both. We chose these architectures because they are well established in the deep learning field and are good examples of the kind of models we can find used in related problems research (Wei et al., 2021; Zhang et al., 2017; Zuallaert et al., 2018).

In the experiments, the deep learning models have been optimized using the validation dataset mentioned in Section 2.1. The values of the hyperparameters have been found using a grid search cross-validation process. A batch size of 32, a learning rate of 0.001, and the Adam optimizer (Kingma & Ba, 2017) were used in the experiments for every deep learning model. The number of epochs was set to 100, although an early stopping mechanism would stop the training process after 10 epochs without any improvement. All DNNs have been trained by applying binary cross entropy as their loss function.

Later, in the experiments, the architectures used for the neural networks are faithful to those used in the papers in which these architectures are presented, with a number of parameters ranging from 100,000 to 600,000. Modifications that increase the number of parameters have been avoided, because this could lead to generalization problems in some of the models.

Convolutional neural networks (CNNs)

CNNs are a kind of artificial neural network based on the use of filtering operations called convolutions, which are usually used as a complement to classic perceptrons and pooling layers. These convolutions allow the network to represent the input information with different feature maps. CNNs have excelled in tasks related to image classification (Krizhevsky, Sutskever & Hinton, 2012), to the point that it has even led scientists to look for ways to encode all kinds of data as if they were images (Sharma et al., 2019) and turn different problems, for example, variant calling (Poplin et al., 2018), into image classification tasks, in which CNNs can achieve excellent results (Nguyen et al., 2016). In this case, we are working with a similar representation of the data, coming from DNA sequences of 1003 bases one-hot encoded in arrays of length 4. Thus, having N samples, we start with an input of N × 1003 × 4, which can be seen as a set of 1D images with four channels.

The structure of the CNN (Fig. 2A) is based on Bhandari et al. (2021) considering that it has been proven effective and that we want to keep similar structures with comparison purposes. Therefore, the network is composed by three 1D convolution layers of size 5, stride 1 and 32 filters. Each convolution layer is followed by a max-pooling layer with size 4. Following these steps, there are three fully connected layers with 1,024, 512 and 128 units, each one followed by a 20% dropout step (Salman & Liu, 2019). Finally, there is a classification layer using the sigmoid activation function.

Metrics

To evaluate the models, we have employed different metrics, each adapted to a different context. The first one is the accuracy (Acc) (Eq. (3)) and it is acceptable as a first approach to measure model performance. (3)${\displaystyle Acc=\frac{Numberofcorrectpredictions}{Totalnumberofpredictionsmade}.}$

However, accuracy is not the best option in situations involving imbalanced data. Because of that, we decided to use more appropriate metrics, like binary cross-entropy (BC) (Eq. (4)) which is a standard metric in binary classification (4)${\displaystyle BC=-\frac{1}{N}\sum y_i\cdot log\left(p\left(y_i\right)\right)+\left(1-y_i\right)\cdot log\left(1-p\left(y_i\right)\right)}$ where y is the actual class and p(y_i) is the probability that the model assigns to the instance i being positive. It will be 1 if the instance belongs to the positive class and 0 otherwise. A good model gives low values of BC.

Another metric that works well in the context of imbalanced data is the F1-score (Eq. (7)) which takes the harmonic mean of the precision and recall metrics, that are defined in terms of the true positives (TP), false positives (FP), true negatives (TN) and false negatives (FN) all obtained from the confusion matrix. (5)${\displaystyle Precision=\frac{TP}{TP+FP}}$ (6)${\displaystyle Recall=\frac{TP}{TP+FN}}$ (7)${\displaystyle F1=\frac{2\cdot Precision\cdot Recall}{Precision+Recall}.}$

Finally, we used the area under the ROC curve (AUC), a thorough metric for imbalanced datasets, which is the area under the curve when plotting the sensitivity (Sn) against the specificity (Sp): (8)${\displaystyle Sn=\frac{TP}{TP+FN}}$ (9)${\displaystyle Sp=\frac{TN}{FP+TN}.}$ The AUC values are between 0 and 1 with values closer to 1 associated with better predictions.

Each metric in this study was calculated using a five-fold validation process. This process has a minor peculiarity since the training, validation and test sets are defined by chromosomes, which prevents the use of a classical five-fold cross-validation. In the case of the deep leaning models, we retrained the models five times using different seeds. In the cascade SVM experiments, we used random splits of the dataset in the first step of the cascade process, which leads to five independent results.

Evaluation of the positive to negative ratio in the dataset generation

In this article, we also wanted to test whether the use of a ratio of one to ten between positive and negative instances is the correct choice, as suggested in Zhang et al. (2017). Different executions have been run with each method using a one to one and a one to ten ratio between positive and negative instances. Although for the training process a ratio of one to one has been used, for the test process, a one to ten ratio has been kept, given this is the approximated ratio we could find in real data.

It consists of generating only one negative instance for each positive one instead of ten, which could be considered a process of undersampling on the majority class.

We evaluated two methods to obtain a one to one ratio of positive and negative instances. The first is very simple, since the negative instances are generated synthetically with the dataset. It consists of generating only one negative instance for each positive one, instead of ten, which could be considered a process of undersampling on the majority class. For the second method, we used SMOTE (Synthetic Minority Oversampling Technique), which is based on generating synthetic instances to oversample the minority class, leading to a ratio of one positive to one negative instance. The downside of SMOTE is that it needs to be applied to numerical values, and the weighted degree kernels work with sequences. This means that SMOTE can not be applied to our string kernel method.

As we can see in Table 2, most of the methods get the best results when they are trained with a positive to negative ratio of one to ten, instead of one to one. Also, using SMOTE does not seem to improve the performance either. This seems to be in line with the results reported in Blagus & Lusa (2013) where the authors claim that SMOTE does not work well with large numbers of features. This is precisely what is happening here, where in order to apply SMOTE we are using one-hot encoding, what increases to 4012 the number of characteristics of each sequence.

The conclusion that can be drawn from these results is that trying to balance the data set, for example, applying a subsampling to the majority class, is not a good heuristic to improve the results in this type of problem. Nor does any improvement seem to be obtained by generating artificial instances with SMOTE. For this reason, the rest of the experiments will be executed applying a positive to negative ratio of one to ten. Finally, from these results we can deduce that the use of more negative instances is helping to obtain better models and that the method chosen to generate these instances is adequate.

Comparing SVM and deep learning results

The first magnitude that we take into account is the training time of each model. Given that the model only has to be trained once, it could be reasonable to use a model whose training is slower if the predictions are significantly better. In this case, we find notable differences between the methods in Fig. 4, where the training times for each model can be seen.

The times for the deep learning methods vary from nearly 2 h of the CNN to 12 h of the bidirectional LSTM, while for the Cascade SVM method we got execution times of almost 4 days using the string kernel. On the other hand, the SVM with an RBF kernel can be trained in 2 h using the cascade approach and the Intel(R) Extension (https://github.com/intel/scikit-learn-intelex) for Scikit-learn (which, unfortunately, cannot be used when the SVM uses a string kernel). Without the cascade approach, the SVM with a string kernel could take several months before being able to complete the training process.

This has to do with the weighted degree kernel. This kind of kernel can achieve notably better results working with strings, as seen in Perez-Rodriguez, Haro-Garcia & Garcia-Pedrajas (2020), but it comes at a great computational cost. In this case, the kernel has been developed in Python using the Scikit-learn (Pedregosa et al., 2011) tool to use personalized kernels, which implies precomputing the entire kernel matrix. This kernel computation can be optimized using faster languages like C, so we implemented the gram matrix calculation in C to reduce computational time. To minimize these problems, the matrix computing step has also been programmed to be parallelized. In our case, the SVM experiments have been run on an Intel(R) Xeon(R) Gold 6150 CPU @ 2.70 GHz using 20 cores. Despite that, its training process is remarkably slower than that of the other methods.

In order to study the effects of the amount of training data on the performance of each method we performed experiments using the 10% (Table 3) and 50% (Table 4) of the data in order to have a fair comparison between SVM and the rest of the algorithms, given that SVM methods could perform better than deep neural networks in situations with fewer data available. Now, the reduction in training times achieved by working with smaller datasets has allowed us to include regular SVM methods to analyze the effects of cascading optimization on the results.

In Table 3 we can see how, even using only 10% of the data, SVM gets the worst results due to a strong overfitting (note that the regularization parameter was the best found after testing several values using a grid search), except in the case of the RBF kernel using cascade SVM, which makes the model generalize better but notably reduces its prediction performance. These results seem to indicate that, although SVMs are typically more robust to overfitting than deep networks, for this particular problem the measures against overfitting of deep networks, the dropout layers, are more effective than the regularization of SVM. When using half of the data, in Table 4 we can see how, even though the results experience a slight improvement, the SVM models keep overfitting with significantly lower results than the rest of the methods. As we can see, SVMs consistently gets the worst results in every metric but accuracy, which is related to the unbalanced nature of the data. In metrics more suited to unbalanced data, SVM performs significantly worse than the deep learning methods.

Once SVM has been discarded as an adequate model for this problem, our purpose is to analyze which of the five different architectures based on DNNs gets the best results, now using the entire dataset. In Fig. 5 and Table 5, we can see how the CNN model is unable to keep a similar performance to the rest of the methods. Apparently, the model is too simple to detect complex patterns in the sequence. On the other hand, the addition of LSTMs makes the models achieve better results, even higher than the algorithms including a convolutional step in their architecture. This CNN addition makes the model a bit faster by reducing the sequence length given to the LSTM, thus improving GPU parallelization, but its results are slightly inferior. Something that actually helps the model is the bidirectional structure of the LSTM, being the best model for every metric analyzed. This makes its training twice as slow as the standard LSTM, from an average training time of 6.69 to 11.4 h using an NVIDIA Titan Xp GPU, that is, a 70% increase in time. Because of the significant difference in results between the RBF kernel and the rest of the methods showed in Table 5, the RBF SVM is not included in Fig. 5 to provide a clearer visualization of the other algorithms.

Generalization performance

In addition to the performance analysis of the models using the human genome dataset, to test if the obtained models could be used on other datasets without retraining, we created a new dataset from the mouse genome data in Ensembl. The process followed to generate this dataset was the same as already described in Section 2.1, and the dataset is also available from https://github.com/JoseBarbero/EnsemblTSSPrediction.

The results of this last experiment are shown in Table 6. The performance degrades noticeably for all models, something that should not be a surprise, as we are using them on sequences of a different species than the one used to train them. But within this worse performance, it is interesting to see how the SVM models achieve slightly better results than the rest of the models in accuracy and binary cross-entropy, specifically the string kernel. Although, if we focus on the metrics best suited for unbalanced data, we can see how the models that stand out are actually those using deep learning. In particular, it is surprising how, according to the AUC-ROC metric, the LSTM has a generalization far above the rest of the methods. The LSTM model has the simplest architecture of all the deep neural networks evaluated, which could explain a lower overfit to the training dataset leading to a better generalization ability, extending even to sequences of other species.