Promotech: a general tool for bacterial promoter recognition

Promoters are genomic regions where the transcription machinery binds to initiate the transcription of specific genes. Computational tools for identifying bacterial promoters have been around for decades. However, most of these tools were designed to recognize promoters in one or few bacterial species. Here, we present Promotech, a machine-learning-based method for promoter recognition in a wide range of bacterial species. We compare Promotech’s performance with the performance of five other promoter prediction methods. Promotech outperforms these other programs in terms of area under the precision-recall curve (AUPRC) or precision at the same level of recall. Promotech is available at https://github.com/BioinformaticsLabAtMUN/PromoTech. Supplementary Information The online version contains supplementary material available at (10.1186/s13059-021-02514-9).

BPROM [17], respectively, while bTSSfinder achieved higher accuracy than the other three assessed tools. These results are promising as they showed that it is possible to recognize promoters of several bacterial species even when the methods were designed for specific bacterial species. BPROM uses five relatively conserved motifs from E. coli to identify promoters, and bTSSfinder focuses on E. coli and three species of Cyanobacteria. Based on this, we hypothesized that predictive performance can be improved if a method is trained on data from a diverse group of bacterial species.
Here, we developed a general (species independent) bacterial promoter recognition method, Promotech, trained on a large data set of promoter sequences of nine distinct bacterial species belonging to five different phyla (namely, Actinobacteria, Chlamydiae, Firmicutes, Proteobacteria, and Spirochaetes). As promoters are typically located directly upstream of the transcription start site (TSS), we used published TSS global maps obtained using sequencing technology such as dRNA-seq [21] and Cappable-seq [22] to define promoter sequences. We trained and evaluated twelve random forests and recurnrent neural network models using these data to select Promotech's classification model (Fig. 1). Finally, we compared the performance of Promotech with that of five other bacterial promoter prediction methods on independent data from four bacterial species. Promotech outperformed all these five methods in terms of the area under the ROC curve (AUROC), the area under the precision-recall curve (AUPRC), and precision at a specific recall level.

Variety of training and validation data
We obtained a large amount of promoter sequences from published global TSS maps (listed in Table 10). On both the training and the validation data, we had bacterial species belonging to distinct phyla and having a wide range (from 30 to 72%) of GC content (Tables 1 and 2). In total, our training data contained 27,766 promoter sequences, and our validation data contained 11,615 promoter sequences (Supplementary Tables S2 and S3).

Model selection
To select Promotech's classification model, we built two and ten models using random forest (RF) [23,24] and recurrent neural networks (RNN) [25] (Fig. 1), respectively. The RF models consisted of one trained with hot-encoded features (RF-HOT) and another trained with tetra-nucleotide frequencies (RF-TETRA). To calculate the tetra-nucleotide frequency vector of a given sequence, the number of occurrences of each possible  (4-mer) in that sequence is divided by the total number of 4-mers in it. The optimal parameters of the RF models were selected using 10-fold grid search cross-validation on 75% of the training data. The RNN models consisted of five long short-term memory (LSTM) [26] and five Gated recurrent unit (GRU) [27] models having zero to four hidden layers and a word embedding [28] layer to obtain a numerical representation of the promoter sequences (Table S4 in Additional file 1). From now on, these RNN models are denoted as GRU-X or LSTM-X, where X indicates the number of hidden layers. All the models were trained using an unbalanced dataset with a 1:10 ratio of positive to negative instances to simulate the small number of promoters in a whole bacterial genome. Due to time constraints, we were unable to run grid search crossvalidation for the RNNs and assessed these models by randomly splitting the training data into 75% for training and 25% for testing. The best performing models per machine learning method in terms of AUPRC and AUROC were RF-HOT, GRU-1, and LSTM-4, as shown in Table 3. RF-HOT was the model with the highest AUPRC overall.

Model interpretation
To interpret the models created, we performed feature importance analysis to find out motifs recognized by the models. To do this, we obtained the feature importance ranking  from the RF models. First, the importance scores were calculated using permutationbased importance (also called mean decrease in accuracy) [29] and impurity-based importance [24] on RF-TETRA. The most important tetra-mers based on the impuritybased importance score from RF-TETRA were TATA, ATAA, TAAT, TTAT, AAAA, and TTTT. The test was repeated using the permutation-based importance score and permuting each feature five times. Both tests produced similar results having the same tetra-nucleotide sequences appearing at the top of the ranking, only varying their relative ranking position and score (Tables 4 and 5).
The feature importance analysis was repeated on the RF-HOT model. In RF-HOT, each feature represents the presence of one of the four possible nucleotides: adenine (A), thymine (T), guanine (G), and cytosine (C) for the current position in the range of − 39 to 0 relative to the TSS. Each nucleotide was represented as a 4-digit binary number, i.e., A (1000), G (0100), C (0010), and T (0001). The permutation and impurity-based feature importance ranking generated by RF-HOT provided the most important positions in the range of − 39 to 0 relative to the TSS and the nucleotide with the most relevance for each position. To have visual representations of these results, each nucleotide's importance score was plotted on a bar graph (Figs. 2 and 3). These results suggest that having adenine (A) and thymine (T) in the range of − 8 to − 12 relative to the TSS is highly important for promoter recognition. Additionally, these results suggest that the RF models learn to identify the Pribnow-Schaller box [30,31], which is a six-nucleotide consensus sequence (TATAAT), commonly located around 10 bp upstream from the TSS.

Genome-wide promoter prediction assessment
We designed our first two assessments to demonstrate that Promotech was able to make predictions for a whole bacterial genome. Five models were selected for this assessment: the best models per machine learning method (RF-HOT, GRU-1, and LSTM-4), GRU-0, and LSTM-3. GRU-0 and LSTM-3 were selected to evaluate the effect in performance if the model had one less hidden layer. As the models were trained on 40-nt-long sequences, to do whole-genome predictions, we needed to cut the genome in 40-nt-long sequences. Thus, we traversed each genome with a sliding window with a one nt step and a 40-nt window size. The sliced sequences were then pre-processed and fed to the model twice, first, using a forward strand configuration and then using a backwards strand configuration. These steps were repeated for each bacterium on the validation set, namely, M. smegmatis, L. phytofermentans, B. amyloliquefaciens, and R. capsulatus. This was a computationally demanding assessment, as, for example, the sliding window created 6,988,167 sequences of 40 nt when used on the M. smegmatis genome. Each sequence was given a second time with a backward strand configuration ending up with 13,976,336 sequences. Thus, each of the four selected models was executed roughly 14 million times for the M. smegmatis genome. The process took around 4 h to run per bacterium per model, including the sliding window, data pre-processing, and model's execution. Increasing the step size decreased the execution time but also decreased the model performance.
In the first assessment, predicted promoters were considered true positives if they have at least 10% sequence overlap with an actual promoter. Table 6 shows the average AUPRC  RF-HOT achieved the best overall AUPRC (0.14) and AUROC (0.82) ( Table 6). An AUPRC of 0.14 might seem low, but if one considers that there are millions of 40-nt-long sequences in a bacterial genome and only a few thousand of these sequences are actual promoters, then this performance is much better than random guessing. For example, M. smegmatis has four thousand actual promoters (Table 10) and 14 million 40-nt-long genomic sequences; thus, a random classifier has an average AUPRC of 0.0003. RF-HOT achieved an AUPRC of 0.27 in M. smegmatis genome which is roughly a thousand-fold improvement over random performance.
To gain insight into the behavior of the models, we visually inspected the location of the predicted promoters and observed that many predicted promoters were located nearby the actual promoters (Fig. 4). To account for this, we re-evaluated the models' performance to count as correct predictions those within 100 nt of an actual promoter. We called this task "the cluster promoter prediction. " Assessing the performance of the models using the cluster promoter prediction method increased AUPRC 2 to 6 times and AUROC by 1.5 times the values obtained in the first assessment (Figs. 5, 6, 7, and 8   Table 7). This suggests that our models predict promoters in the proximity of actual promoters but are unable to recognize the exact genomic location of the actual promoters.

Performance comparison with existing methods
To compare Promotech's performance with that of other existing methods, we used four independent test datasets containing promoters found by global TSS mapping using sequencing technologies. As in the previous assessments, the independent datasets included B. amyloliquefaciens, L. phytofermentans, M. smegmatis, and R. capsulatus. On this assessment, the datasets have a 1:1 ratio of positive to negative instances (Supplementary Table S3). As these datasets contained thousands instead of millions of sequences, we were able to include the RF-TETRA model that failed to run on a whole genome (due to memory issues). As Promotech's goal is to be applicable to a wide range of bacterial species, we compared Promotech models with other multi-species methods such as bTSSFinder [11] and G4PromFinder [10]. Additionally, we included BPROM [17] in the comparative assessment, as it is the most commonly used promoter prediction program (as per Google Scholar, BPROM's manuscript has been cited 547 times). Finally, we also compared Promotech's performance with two recent E. coli-specific methods: MULTiPly [4], designed for various sigma factors, and iPro70-FMWin [6], designed for sigma 70. In total, this benchmark included Promotech's six models and five other bacterial promoter prediction tools.
Promotech's random forests models (RF-TETRA and RF-HOT) consistently achieved the highest AUPRC and AUROC across the four bacterial species (Tables 8 and 9). RF-TETRA achieved the highest average AUPRC and AUROC among all the methods.  Among the five other bacterial promoter prediction tools, iPro70-FMWin showed the best predictive performance but still substantially lower than Promotech's. Based on these results, we selected RF-HOT as Promotech's predictive model for genome-wide promoter prediction. For recognizing promoters on 40-nt-long genomic sequences in data sets containing up to thousands of sequences, we recommend both RF-HOT and RF-TETRA models.
Additionally, we compared the performance of RF-TETRA and RF-HOT for identifying E. coli promoters against that of E. coli-specific tools. To do this, we obtained Promotech's predictions on a balanced data set with 2860 experimentally validated E. coli promoters collected from RegulonDB [32]. This data set has been used to evaluate the performance of several E. coli promoter prediction tools [3,19,33]. The average 5-fold cross-validation MCC and accuracy reported on this data set [3,19,33]

Conclusions
Based on our results, we recommend (1) to use E. coli-specific tools to predict E. coli promoters as they can identify E. coli promoters more accurately than a general bacterial promoter identification method such as Promotech and (2) to use Promotech to identify promoters in bacterial species other than E. coli, as we have shown Promotech   (Tables 8 and 9).
In sum, Promotech is a promoter recognition tool suited for general (species independent) bacterial promoter detection that is able to perform promoter recognition on a whole bacterial genome. Promotech is available under the GNU General Public License v3.0 at [34,35].

Methods
The goal of this study was to develop a general tool to recognize bacterial promoters. To do this, we assembled a large data set of promoter sequences from various bacterial species, generated twelve machine learning models, and selected the best models based on AUROC and AUPRC. Our best models were compared with five existing tools (BPROM [17], bTSSFinder [11], MULTIPly [4], iPro70-FMWin [6], and G4PromFinder [10]) using a validation data set, not used for training, of four bacterial species.

Collecting data
Bacterial TSS detected by next-generation sequencing (NGS) approaches, namely, dRNAseq [21] and Cappable-seq [22], were collected from the literature (Table 10). We obtained promoter genomic coordinates and the corresponding sequences using BEDTools [36]. Eσ covers DNA from roughly 55 bp upstream to 15 bp downstream of the TSS [1]. As the promoter region is not located downstream of the TSS and Eσ covers 15 bp downstream  of the TSS, we excluded 15 bp from both sides of the region covered by Eσ and took as the promoter sequence the 40 bp upstream of the TSS. The bacterial species included in this study are listed in Table 10.

Generating positive and negative instance sets
A Nextflow [37] pipeline was designed to obtain the promoter (positive) sequences from the TSS coordinates by taking the genome FASTA file and the TSS coordinates as input and obtaining the promoter coordinates as 40 bp upstream from the TSS to the TSS using BEDTools' slop command. The BEDTools' subtract command was used to delete duplicates, and the getfasta command was used to obtain the FASTA sequences from the promoter coordinates.
To obtain non-promoter (negative) sequences, we used BEDTools' random command to obtain random genomic coordinates and the getfasta command to obtain the corresponding genomic sequences. Negative sequences overlapping positive sequences were excluded from the training data. Note that some of these negative instances might in fact be actual promoters, and thus, predictive performance is conservatively assessed.
The training data sets created have a 1:10 ratio of positive to negative instances (unbalanced). For the validation data set, we created a data set with a 1:1 ratio of positive to negative instances (balanced). The total number of positive and negative instances per bacterium is shown in Supplementary Tables S2 and S3 in Additional file 1. AUPRC is roughly the weighted average precision across all recall levels. A perfect classifier has an AUPRC of 1, while a random classifier has an AUPRC of 0.5 in a balanced data set. These results were obtained in balanced data sets (i.e., with a 1:1 ratio of positive to negative instances). The numbers in bold indicate the model with the highest AUPRC. For BPROM, bTSSFinder, and G4PromFinder, the numbers between brackets indicate precision and recall achieved as these tools did not provide a probability associated to each instance in the data set AUROC is roughly the likelihood that a positive instance will get a higher probability of being a promoter sequence than a negative instance. These results were obtained in data sets (not seen during training) with a 1:1 ratio of positive to negative instances. The numbers in bold indicate the model with the highest AUROC. For BPROM, bTSSFinder, and G4PromFinder, the numbers between brackets indicate true-positive rate and false-positive rate obtained as these tools did not provide a probability associated to each instance in the data set

Machine learning models
We used two machine learning methods: recurrent neural networks (RNNs) [25] and random forest (RF) [23,24]. Both methods have been successfully used before to classify genomic sequences. Random forest is a popular machine learning method for its ability to identify feature importance and handles many data types (continuous, categorical, and binary). It is well-suited for high-dimensional data and avoids over-fitting by its voting scheme among the ensemble of trees within it [38]. RNNs are also well-suited for genomic sequence analysis due to their ability to handle variable-length inputs, detecting sequential patterns and retaining information through time.

Model selection
Due to the lengthy training time of the RNNs, we were unable to run CV for the RNN models. Thus, to select the best model, we split our training data in 75% for training and 25% for testing. Models were trained with 75% of the data and then compared to each other on their performance in the 25% left-out data. After selecting the best models, these models were retrained using all of the training data and the resulting models used for whole-genome promoter prediction and comparative assessment with the other tools.

Random forest
Two RF models were generated; the first was trained using hot-encoded features; this meant that the nucleotides (A, G, C, T) were transformed into binary vector representations [1000], [0100], [0010], and [0001], respectively. This model is henceforth referred to as RF-HOT. The second model was trained using tetra-nucleotide frequencies calculated using the scikit-bio library [39] and denoted as RF-TETRA. The models were created using the Sklearn's RandomForestClassifier [40] combined with a 10-fold grid search CV to handle the hyper-parameter optimization. The hyper-parameter search space was max_features (m): [None, "sqrt", "log2"] and n_estimators (n): [1000, 2000, 3000]; both models were trained using an unbalanced data set with a 1:10 ratio of positive to negative instances. The best hyper-parameters found by grid search CV for both RF models were m = "log2" and n = 2000 with class weights values of {0 : 0.53, 1 : 10.28}. In the last column, a T or V indicates whether the bacterium is reserved for training or validation, respectively. Additional information is included such as the number of TSS per bacterium, the genome's length, the next-generation sequencing technology used to obtain the TSSs, and the literature sources' PubMed ID ( if PubMed ID is missing, then at the time of this publication, the source manuscript was still in preparation

Recurrent neural networks
Two types of RNNs were trained, long short-term memory unit (LSTM) [26] and gated recurrent unit (GRU) [27], using word embeddings representation of the tetra-nucleotide sequences calculated using the Keras' Tokenizer class [41]. The models were designated as GRU-X or LSTM-X where X indicates the number of hidden layers. All models were manually tuned with an architecture consisting of one embedding layer, one GRU or LSTM layer, zero to four dense layers with dropout to reduce overfitting, one binary output, Adam optimizer function, and binary cross-entropy loss function (Table S4 in Additional file 1).

Computer infrastructure
All RF and RNN models were trained on the Compute Canada's Beluga Cluster [42] configured with four NVidia V100SXM2 16GB GPUs, eight Intel Gold 6148 Skylake @ 2.4 GHz CPUs, and managed using SLURM commands.

Model assessment
Three assessments were performed to evaluate the models' performance. The first consisted in scanning each bacterial genome using a 40-nt sliding window. In total, the number of generated sequences ranged from 4 to 7 million depending on the genome size. Models were given as input each of the sliding window sequences. Models then outputted per 40-nt sequence the probability of having a promoter within that sequence. To be counted as a true positive, the predicted promoter sequence had to have at least 10% sequence overlap with an actual promoter sequence. All other predicted promoters were considered false positives. To determine whether a predicted promoter overlapped with an actual promoter, we used the BEDTools intersect command with the parameters -s and -f 0.1.
In the second assessment, we also scanned each bacterial genome using a 40-nt sliding window. However, in this assessment, we considered a predicted promoter a true positive if it was within 100 nt of an actual promoter. In this setting, we use BEDtools' closest command to find the five predicted promoters closest to an actual promoter. Then, those closest predicted promoters less than 100 nt away, upstream or downstream, from an actual promoter were counted as true positives. All other predicted promoters were considered false positives.
In the third assessment, we used the validation balanced data set obtained as described above. In this assessment, we included BPROM, bTSSfinder, G4Promfinder, iPro70-FMWin, and MULTiPLy. MULTiPLy and G4PromFinder accept 40-nt-long sequences, but BPROM and bTSSFinder require sequences 250-nt-long and iPro70-FMWin requires 81nt-long sequences. We used BEDTools' slopBed [38] command to extend the sequences in the data set from 40 to the required length. As G4Promfinder is written in Python, we integrated into our own pipeline and fed the sequences directly to G4Promfinder. To run, BPROM and bTSSfinder, we wrote each sequence to a file and then ran the programs through a shell script called from our own pipeline. MULTiPLy was tested separately, as it was developed in Matlab, so a short script was written to feed it each bacterium's data set. We ran BPROM and MULTiPLy with their default values. bTSSFinder was run with the parameters -c 1 -t e and -h 2 indicating to search for the highest ranking promoter regardless of promoter class, use E. coli mode, and search on both strands. All other bTSSFinder parameters were left at their default values. G4Promfinder only requires as input sequences in FASTA format. iPro70-FMWin's results were obtained from its website inputting the sequences in FASTA format.
Performance metrics used were (1) the area under the precision-recall curve (AUPRC), where precision is the number of true positives divided by the total number of predicted positives and recall (also called sensitivity or true-positive rate) is the number of true positives divided by the total number of actual positives and (2) the area under the ROC curve (AUROC), where true-positive rate is the same as recall and false-positive rate is the number of false positives divided by the total number of predicted positives.
Additional file 1: The Additional file 1 is a PDF file that includes the following tables: Table S1 -A summary of promoter prediction approaches in the last twelve years.

Acknowledgments
This research was enabled in part by support provided by ACENET (www.ace-net.ca/) and Compute Canada (www. computecanada.ca).

Peer review information
Tim Sands was the primary editor of this article and managed its editorial process and peer review in collaboration with the rest of the editorial team.

Review history
The review history is available as Additional file 2.  Table 10 lists the studies used during data collection. Promotech is available under the GNU General Public License v3.0 at [34]. All data, scripts, and pipelines used are available at [35].

Declarations
Ethics approval and consent to participate Not applicable.

Consent for publication
Not applicable.