TransPrise: a novel machine learning approach for eukaryotic promoter prediction

Selection of genome annotation version

We selected rice chromosomes and Genome Annotation release 7 (MSUv7, http://rice.plantbiology.msu.edu). There are two commonly used annotations of rice: MSU (Kawahara et al., 2013) and FgeneSH (Zhang et al., 2008). The Fgenesh gene prediction set contains 18,389 high quality (5′ full, with mRNA support) gene models, while the MSU gene prediction set contains 20,367 high quality gene models (Tatarinova et al., 2016). We used Fgenesh mRNA-based gene prediction models, since Fgenesh-annotated promoters have a more pronounced nucleotide consensus as compared to the promoters annotated by MSU (Triska et al., 2017b). Fgenesh was successfully used to annotate several plant genomes (Chan et al., 2017a; Chan et al., 2017b; Davis et al., 2010; Ito et al., 2005; Jiang et al., 2015; Nasiri et al., 2013; Sanusi et al., 2018; Sheshadri et al., 2018; Yao et al., 2005). Therefore, we selected the Fgenesh annotation as the gold standard for our analysis. To obtain the highest quality dataset, pseudogenes, transposable elements, and genes with 5′ UTR shorter than 20 nt or longer than 1,000 nt have been excluded.

Training, validation and test sets

The procedure consists of two steps: classification (dividing the genome into “promoters” and “non-promoters”) and regression (finding the position of TSS inside the sequence identified as “promoter”). The genomic sequence of Oryza sativa had been divided into the testing and training sets. The training set contains three files:

1. Training “non-promoter” dataset contains sequences extracted from random genomic positions separated from experimentally validated transcription start sites by 2,000 nt. This dataset contains mostly intergenic regions. All sequences are 2,000 nt long.

2. Training “promoter” dataset contains sequences [TSS-1,000; TSS+999] from the all chromosomes with length 2,000 nt.

3. File with indicators of TSS positions, containing (2,000 × 1) matrices that correspond to positions of biologically validated TSS in every training sequence (“1” TSS, “0” not TSS position).

The same set of files was created for the testing dataset. Since the procedure has multiple steps (classification and regression), training and testing sets were selected at each step of the method.

The following procedure was used to assemble the dataset for the classification model:

1. “Non-promoters”: $\frac{1}{4}$ of the examples chosen from the training “non-promoter” dataset, randomly selecting 512 nt long sequences from 2,000 nt long regions.

2. “Promoters sans TSS”: $\frac{1}{4}$ of the examples were randomly selected from the training “promoter” dataset, making sure that the chosen 512 nt long fragment did not overlap the region [TSS-50, TSS+50].

3. “TSS vicinity”: $\frac{1}{2}$ of the examples extracted from the training “promoter” dataset, containing only one TSS in a random position within the 512 nt long sequence, with a restriction that it should be in the [250, 450] fragment.

The dataset for the regression model was assembled using sequences that contain one validated TSSs in a randomly selected position of the [250, 450] fragment. The datasets are represented as (512 × 4) nucleotide matrices M with 512 columns and 4 rows. The 1st row contains delta function δ(x_i = A)—it is equal to 1 if there is nucleotide “A” in the ith position of the sequence and 0 otherwise. Similarly, 2nd, 3rd, and 4th rows correspond to nucleotides C, G and T.

Model training

We implemented the Convolutional Neural Networks (CNN) using the Keras library for training (https://keras.io/).

Classification and Regression models training

The dataset for the classification contains equal numbers of positive and negative examples. The matrices (512 × 4) described above are input into the model. The CNN architecture (Fig. 1) started with four parallel convolutional layers (composed of 128 filters with 2, 4, 8 and 16 kernel sizes) ReLU, was used as activation function followed by concatenation. After concatenation layer we used convolution, batch normalization, max pooling layers twice. The first convolution had 128 filters and the second had 16. There was one kernel size and ReLU activation in both situations. To help regularize the model, we used the 0.5 Dropout technique. The signal is fed to two fully connected layers with ReLU activation functions consisting of 256 and 128 neurons, followed by batch normalization. The output layer had a sigmoid activation function.

Whole Genome Sequencing (WGS) processing

TransPrise also works with WGS data in fasta format. We prepared a Python script that with 250 min run time on a genome of 374 Mb long. It uses a sliding window for extraction of 512 nt long sequences with the step size 4 nt. If the classification step identifies a fragment as TSS-containing, this fragment will be passed to the regression step. The prediction vector is compared to the TSS model using a L1 norm. The default value for the similarity threshold is 600, but it can be modified by users. The only restriction of our algorithm is that alternative TSSs must be more than 100 nt of each other.

Model evaluation

For the purpose of the K-fold cross-validation, the dataset is randomly divided into K equal-size subsets. Of the K subsets, a single subset is retained as the validation data for testing the model, and the remaining (K-1) subsets are used for training. The cross-validation process is then repeated K times (the “folds”), with each of the K subsamples used exactly once as the validation data. Then, the results from K folds are averaged. The advantage of K-fold cross-validation is that all observations are used for both training and validation, and each observation is used for validation exactly once. We conducted 10-fold cross-validation (dataset was divided into training and validation sets in 9:1 ratio; validation set was used to avoid overfitting and find the optimal number of learning epochs). The ROC curves obtained in 10-fold cross-validation are presented in the “Results” section. We determined that the optimal number of learning epochs is five. After the model training, we tested our model using the test set and calculated Accuracy (Ac), Sensitivity (Se), Specificity (Sp), and the Matthews Correlation Coefficient (CC): ${\displaystyle \mathrm{Specificity}=\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FP}},}$ ${\displaystyle \mathrm{Sensitivity}=\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}},}$ ${\displaystyle \mathrm{Accuracy}=\frac{\mathrm{TP}+\mathrm{TN}}{\mathrm{TP}+\mathrm{TN}+\mathrm{FP}+\mathrm{FN}},}$ ${\displaystyle \mathrm{CC}=\frac{\mathrm{TP}\times \mathrm{TN}+\mathrm{FP}\times \mathrm{FN}}{\sqrt{\left(\mathrm{TP}+\mathrm{FP}\right)\left(\mathrm{TP}+\mathrm{FN}\right)\left(\mathrm{TN}+\mathrm{FP}\right)\left(\mathrm{TN}+\mathrm{FN}\right)}},}$ where TP—true positive, TN—true negative, FP—false positive, FN—false negative. The input of the regression model has the same format as the classification model input. TSS is assumed to be located at a random position between nucleotides 250 and 450. There is only one difference between classification and regression models: in the output layer the activation function is replaced by a linear function.

We trained the classification and regression models using the Keras Library and performed 10-fold cross validation procedures for them. The algorithm for TSS prediction is presented in Fig. 1. For every fold, we carried out five learning epochs. The complete learning time was 35 s on average. We performed 10-fold cross-validation and calculated the average value of mean absolute error (MAE) to estimate the accuracy of TSS position prediction, where y_i—position of TSS in test set (assumed to be accurate), and x_i—predicted position of TSS. ${\displaystyle \mathrm{MAE}=\frac{\sum _{i=1}^n\left|y_i-x_i\right|}{n}.}$

Figure 2 shows the flowchart for model building and validation.

Classification

To reduce the influence of how the data is split on the resulting testing statistics, we carried out a 10-fold cross-validation (CV) procedure for “internal” validation and supplemented it by the “external” validation of the rice chromosome 2 (excluded from the training set). We have randomly divided the dataset into 10 partitions and used each partition as the testing data while training the model on the remaining partitions. The ROC (Receiver Operating Characteristic) curve represents dependence of sensitivity on the specificity. It is a graph showing the performance of a classification model at all classification thresholds. AUC-ROC curves for classifier obtained in 10-fold “internal” cross-validation are presented in Fig. 3. Accuracy = 0.88, Se = 0.84, Sp = 0.92, CC = 0.79, AUC = 0.94.

Then we have selected the best model and evaluated its performance on rice chromosome 2 (“external” validation dataset) and compared it with TSSPlant. The dataset for “external” validation contained 2,000 nucleotide sequences with length 512 nt, with 1,000 examples—“non-TSS” sequences and 1,000—“TSS” sequences. For the “external” validation we calculated the Matthews correlation coefficient (MCC), Accuracy (Ac), Sensitivity (Se), Specificity (Sp), and Area Under the ROC Curve (AUC-ROC) for classification models. The results of the “external” validation are presented in Table 1.

We have clearly shown that the classification model (Stage 1 of the pipeline) already achieves high accuracy (0.88). Regression (Stage 2) will further refine the prediction.

Regression

Figure 4 shows the error density curves obtained in the 10-fold cross validation procedure for regression models. The mean absolute error (MAE) for regression model was 29.19 nt. The mean difference between predicted and true TSS, if the guess is random, is 66 nt.

We compared the accuracy of promoter prediction for TSSPlant and TransPrise. For this purpose, we selected rice chromosomes 1 and 2. The testing was performed for entire chromosomes, without filtering for possible TSS containing regions, located upstream of the ATG.

These two chromosomes contain 5,298 high quality experimentally validated TSS, defined as follows:

1. Locus does not correspond to transposable element

2. Locus has experimental support (full-length mRNA)

3. If multiple isoforms are predicted, the “representative” is used

4. Size of the 5′ UTR is at least 20 nt.

For these two chromosomes, TSSPlant predicted 153,009 TSSs, while TransPrise has found 13,765 sites. Of TSSPlant predictions, 10,721 (∼7%) were located within 1,000 nt from validated TSS. Of TransPrise predictions, 3,989 (∼29%) were located within 1,000 nt from validated TSS. Additionally, TransPrise predictions tend to be closer to validated TSS than TSSPlant predictions (Fig. 5).

Validation of TransPrise in Homo sapiens

We have used Cap Analysis of Gene Expression (CAGE) data from the DBTSS database (https://dbtss.hgc.jp/), using the HG19 version of the human genome (Suzuki et al., 2018). The model was trained on all chromosomes except chromosome 8 and tested on chromosome 8. Accuracy of the classification model was 0.778, sensitivity 0.816, specificity 0.74, and Matthews correlation coefficient 0.57. The mean absolute error of the regression model was 47.986 nt.