Comparison of machine learning and deep learning techniques in promoter prediction across diverse species

Selection of dataset

We selected yeast, A. thaliana and human genomes for our analysis. The genomes were selected using UCSC genome browser (Haeussler et al., 2019). We collected the datasets of approximately 41,671 sequences of A. thaliana (UCSC version araTha1), 61,546 sequences of humans (UCSC version hg38) each, and 6125 sequences of yeast (UCSC version sacCer3). However, after pre-processing and data cleaning, we randomly selected 35,000, 35,000 and 6000 sequences for A. thaliana, human and yeast respectively for experimentation. One thousand basepairs (bp) long putative promoter regions (−700 to +300 bp around TSS) were extracted in two steps as follows: As a first step, 700 bases upstream of the TSS were selected in the Table Browser using RefSeq Genes as an input to create a custom track. As a second step, 300 bases downstream were selected using custom track as an input to extract a 1,000 bases promoter sequence in a FASTA format. Equal number of background sequences were generated using a synthetic method of shuffling of promoter sequences for each organism (Caballero et al., 2014). The ratio of positive and negative sequences for each organism was 1:1. Each of these datasets were processed for the development of distinct models. The number of samples were divided into training and testing sets with the split of 90% and 10%, respectively. Out of 90% of training samples 10% of data was considered for validation and parameter tuning. We have utilised stratified random sampling for splitting of data, and therefore equal percentage of samples of each target class were used for training and testing of models. These datasets were further subjected to feature extraction, feature encoding and classification purpose (Fig. 1). Number of learning epochs used were 10 and the batch size was 128. Binary cross-entropy and sparce categorical cross-entropy were used as loss function for binary and multiclass classification, respectively. Adam optimizer was utilized with default learning rate. Early stopping callbacks was monitored through loss on the validation dataset. The efficiency of the models was also tested using human sequence data. 1,000 nucleotides long sequences were extracted from the human genome (hg38) using fasta-subsample script of MEME Suite (v 5.0.5) (Bailey et al., 2009). These sequences were used as negative dataset to test the efficiency of the model.

Data pre-processing

DNA sequences are 1D-channel of four nucleotides, adenine (A), thymine (T), guanine (G) and cytosine (C). Feature extraction process is used to generate input for ML/DL models from raw data of nucleotide sequences. For feature extraction we used k-merization process. In this process, a short sequence segment of k consecutive nucleotides (k-mers) with the stride of 1 were generated from parent sequences (Chor et al., 2009). We specifically studied the effect of feature-extraction parameter k-mer size on the prediction performance of the models. We used four different k-mer sizes for experimentation: 1-mer, 2-mer, 4-mer and 8-mer. Further, k-merized sequences were given as input to the feature encoding techniques.

Feature encoding

We used two different approaches for feature encoding viz. frequency-based tokenization and one-hot encoding. In frequency-based tokenization, each k-mer in the sequence is encoded into tokens based on its frequency and mapped to a unique index starting from 1. Zero is a reserved index not assigned to any k-mer. This helps transform each sequence of k-mers to a sequence of integers. Feature encoding was achieved using in-house script written in Python (version 3) with system specification as: 64bit operating system, x64 based processor, 8GB RAM, intel CORE i7 processor. We used a text tokenization utility class of Keras pre-processing and data augmentation module for feature encoding. Whereas for one-hot encoding each k-mer is represented by a new binary variable. For this we utilized the OneHotEncoder class of Scikit-learn preprocessing utility. The output of Frequency based tokenization and one-hot encoding were subjected to CNN based analysis. The output of frequency-based tokenization was also subjected to other ML and DL analyses viz. random forest and LSTM RNN.

Random forest classifier

Random forest is an ensemble of decision trees constructed using a different sample from the original data (Breiman, 2001). First, each tree is built from a random bootstrapped sample of the training data. Second, at each split of the tree, RF model considers only a small subset of features for training. This helps to improve generalization by reducing variance. The final classification is obtained by combining results from the decision trees passed by votes. The bagging strategy of RF can effectively decrease the risk of overfitting when applied to large dimension data. Therefore, RF can handle a large data set with high dimensionality. Random forest has been widely used in the prediction of DNA-binding proteins (Wang, Yang & Jack, 2009), microarray data analysis (Yang et al., 2010), regulatory elements prediction (Li et al., 2017), etc. The encoded data are given as input to the Random Forest Classifier of scikit-learn ensemble module for classification of promoters and non-promoter sequences (Pedregosa et al., 2011). The value of n_estimator (number of trees in the forest) was 300, minimum number of samples required to split an internal node was 2 and bootstrap value was false.

LSTM-recurrent neural network

For sequential data, the flow of gradients for long durations can help in learning long-term dependencies. This can be achieved using Long Short-Term Memory (LSTM) RNN (Hochreiter & Schmidhuber, 1997). We fed the pre-processed encoded input to the LSTM layer. The output of the LSTM layer with 128 units was fed to a Dense layer and the output of the Dense layer with 64 units was fed to the Drop-out layer and used Relu activation. In the Drop-out layer half the units were dropped. The final state was mapped through a fully connected layer with Sigmoid activation. We used LSTM, Activation, Dense, Dropout and Embedding methods from layer module of Keras (version: 2.2.4) to develop this architecture. Performance evaluation parameters are further used to evaluate the performance of the LSTM model.

Convolutional neural network

Encoding helps to represent 1-D channel of DNA sequence in the form of a sequence of numerical values. This form of representation can be used as input to convolutional neural networks (Collobert & Weston, 2008; LeCun et al., 1998). To build CNN architecture, we utilized three 1D convolution layers, each followed by a max-pooling layer. Each convolutional layer uses filter of size 5. The pool size of the max-pooling layer was 4 with stride 1. The output of the max-pooling layer was then fed to three fully connected dense layers, consisting of 1,025, 512, 128 units and used the ReLu activation function with a 20% dropout. The final one is the classification layer and uses the sigmoid activation. We utilized Conv1D and MaxPooling1D methods from the layer module of Keras (version: 2.2.4) for experimentation. Performance evaluation parameters are further used to evaluate the performance of the CNN model.

Performance evaluation parameters

After the application of the prediction model on the dataset, we obtained true positive (TP) and true negative (TN) numbers of truly identified promoter and non-promoter sequences. We also obtained false positive (FP) and false-negative (FN) numbers of falsely identified promoter and non-promoter sequences. In order to evaluate the performance of classification models we computed accuracy (Acc), sensitivity (Sn), specificity (Sp) (Skaik, 2008) and Matthews correlation coefficient (MCC) (Matthews, 1975) using following equations:

${\displaystyle Accuracy=\frac{TP+TN}{TP+TN+FP+FN}}$ ${\displaystyle Sensitivity=\frac{TP}{TP+FN}}$ ${\displaystyle Specificity=\frac{TN}{TN+FP}}$ ${\displaystyle MCC=\frac{\left(\left(TP\ast TN\right)-\left(FP\ast FN\right)\right)}{\sqrt{\left(\left(TP+FP\right)\ast \left(TP+FN\right)\ast \left(TN+FP\right)\ast \left(TN+FN\right)\right)}}.}$

Selection of genomic datasets

We selected three eukaryotic genomes viz. yeast, A. thaliana, and humans, for classification of sequences into the promoter and non-promoter groups. Yeast being the smallest genome, resulted in a small dataset of ∼6,000 promoters. Large genomes of A. thaliana and humans generated datasets of ∼35,000 promoters each. Negative datasets were generated by reshuffling the promoter sequences. The data were divided into train and test sets with a split of 90% and 10%, respectively. The 10% of training data was used for validation and parameter tuning. The resulting datasets were used for comparison of feature encoding techniques and training ML/DL models. In order to evaluate the performance of the ML/DL models, a dataset of 600 random sequences from human genome was used.

Comparison of feature encoding techniques

We tested efficiency and computational requirement of one-hot encoding and frequency-based tokenization (FBT) techniques. For this purpose, we used yeast dataset and tested the feature encoding methods on k-mer sizes 1 and 2 (Table 1). We found that at both k-mer sizes the training time was higher for one-hot encoding without significant improvement in sensitivity and specificity. For one-hot encoding, the testing accuracy of CNN model was 95% for 1-mer and 96% for 2-mer whereas for FBT accuracy achieved was 97% for 1-mer and 96% for 2-mer. We also trained higher k-mer sizes (more than 2-mer) however the training time required was very high at the described configuration. Based on these results and the fact that we wanted to use k-mer sizes 2, 4 and 8; we decided to use FBT for further work.

Comparison of ML/DL algorithms in binary classification

The first set of experiments was conducted to perform binary classification of sequences into the promoter and non-promoter groups. The performance of CNN, LSTM and RF are shown in Table 2 for the yeast, A. thaliana, and human datasets. The average performance of the CNN model on the test set of sequences was: Sn-95.33%, Sp-97.5%, Acc-96%, Mcc-92.5%. It was found that the computed CNN model performed better than the LSTM and RF for all k-mer sizes. For k-mer size 2, RF performed better than LSTM. However, LSTM tends to work significantly better for k-mer sizes 4 and 8 (Figs. 2A–2C). Then we analysed the efficiency of these models on each of the three organisms. In all three species, CNN outperformed LSTM and RF in all evaluation metrics (Figs. 2D–2F). Both CNN and RF do not show any pattern in change in the accuracy with k-mer size. Accuracy of CNN remains high throughout for all k-mer sizes. Accuracy of LSTM on the other hand improves with increase in k-mer size. However, we have observed that the computational requirement for LSTM increases with increase in k-mer size.

Cross species evaluation using ML/DL models

We also performed a cross-species evaluation using CNN and RF. The purpose of the cross-species evaluation was to test whether sequence structures that underlie promoter regions are conserved across species (Lai et al., 2019). In cross-species evaluation, we trained the CNN and RF models using one organism’s data and evaluated on the rest of the organisms for all k-mer sizes considered. We can see that the average accuracy performance of CNN dropped to 72.77% and that of RF dropped to 55.66%. This result suggests that there is a large variation in sequence structures of promoters across species (Fig. 3).

Testing of ML/DL model with random genome sequences

Finally, in order to evaluate the performance of these ML/DL models, we also used random sequences from human genome as negative data with testing data. The prediction performance of CNN model was better than LSTM and RF. The average accuracy achieved by CNN model is 91%, along with 88% sensitivity, 99% specificity and 85% Mcc (Table 3).

Comparison of ML/DL algorithms in multiclass-classification

From cross-species evaluation we found that promoters from the different species have different sequence structures. Therefore, we performed the multispecies classification of sequences to a promoter category. We created a sample dataset of ∼66,000 promoter sequences from yeast, Arabidopsis thaliana, and human. Multiclass classification performance of CNN, LSTM, and RF models for various k-mer sizes are shown in Table 4. The average prediction accuracy score of the CNN model is 98% whereas the average prediction accuracy score for LSTM is 96% and RF is 80%. CNN achieved the highest sensitivity and specificity among all techniques. The CNN and LSTM architectures outperformed the RF classifier for all k-mer sizes (CNN model: Sn-97%, Sp-99% MCC-97%; LSTM model: Sn-92%, Sp-97%, MCC-94%; RF: Sn-58%, Sp-87%, MCC-65%). For multiclass classification of sequences into promoter groups, the performance of RF in terms of accuracy and sensitivity has decreased with an increase in k-mer size. LSTM achieved the highest sensitivity and specificity at k-mer 4. The deep learning techniques CNN and LSTM outperformed RF (Fig. 4A), whereas change in the k-mer size does not show a significant difference in the prediction performance (Fig. 4B).