{'fraction_non_alphanumeric': 0.031588686227861484, 'fraction_numerical': 0.011736716891356067, 'mean_word_length': 5.091948470209339, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
iPro-WAEL: a comprehensive and robust framework for identifying promoters in multiple species

Abstract Promoters are consensus DNA sequences located near the transcription start sites and they play an important role in transcription initiation. Due to their importance in biological processes, the identification of promoters is significantly important for characterizing the expression of the genes. Numerous computational methods have been proposed to predict promoters. However, it is difficult for these methods to achieve satisfactory performance in multiple species. In this study, we propose a novel weighted average ensemble learning model, termed iPro-WAEL, for identifying promoters in multiple species, including Human, Mouse, E.coli, Arabidopsis, B.amyloliquefaciens, B.subtilis and R.capsulatus. Extensive benchmarking experiments illustrate that iPro-WAEL has optimal performance and is superior to the current methods in promoter prediction. The experimental results also demonstrate a satisfactory prediction ability of iPro-WAEL on cross-cell lines, promoters annotated by other methods and distinguishing between promoters and enhancers. Moreover, we identify the most important transcription factor binding site (TFBS) motif in promoter regions to facilitate the study of identifying important motifs in the promoter regions. The source code of iPro-WAEL is freely available at https://github.com/HaoWuLab-Bioinformatics/iPro-WAEL.


INTRODUCTION
The transcription process mainly involves initiation, elongation and termination, in which initiation is the most complicated (1). Promoters located proximal to the transcription start sites (TSS) are non-coding DNA regions, which are essential for initiating the transcription of a particular gene by cooperating with RNA polymerase (RNAP) (2,3) and are critical for gene expression in different species. In prokaryotes, promoters are involved in many biological processes, such as the transcription of most genes (4), heatshock response, nitrogen fixation (5), the expression of flagella (6) and so forth. In eukaryotes, the initiator region or the downstream promoter element of promoters controls the exact position where transcription starts (7). Additionally, promoters often cooperate with their distal regulatory elements via chromatin loops and tend to be involved in developmental diseases, tumorigenesis and spatiotemporal gene expression (8)(9)(10). Therefore, the identification of promoters is crucial to investigating the regulation of gene expression.
Despite advances in computational methods for predicting promoters, several shortcomings in these methods remain to be addressed. Firstly, the feature encoding schemes in their studies and the architecture of their models is relatively simple, resulting in unsatisfactory performance in predicting promoters. Secondly, most of their models have validated the performance of predicting promoters on only one or few species. Therefore, it is not clear whether these methods can be implemented for the identification of promoters in multiple species. Although CNNProm (18) and DeePromoter (19) reported their performance on a few datasets, it is time-consuming to retune the parameters of the models for each species. Finally, the model interpretation is critical for exploring transcription factor binding site (TFBS) motifs in the sequences involved in promoters, but unfortunately, previous studies have not touched upon this aspect.
Therefore, we propose a comprehensive and robust framework, named iPro-WAEL (a Weighted Average Ensemble Learning-based model for identifying promoters), to identify promoters in multiple species by integrating an RF model and a CNN model. Based on 13 datasets of seven species, we compare it with the performance of eight promoter prediction models on without retuning parameters. Besides, we explore the TFBS motifs that have a significant impact on human promoter regions. Our main contributions are as follows: (i) We generate four new human datasets using stricter criteria. (ii) We propose a novel and robust weighted average ensemble learning-based model (iPro-WAEL) to identify the promoters in multiple species. (iii) We demonstrate the optimal performance of iPro-WAEL on multiple species and it is superior to stateof-the-art predictors. (iv) We identify the most important TFBS motifs in promoter regions which are consistent with the previous studies and identify many potentially important but previous underexplored motifs. The overall framework of this study is shown in Figure 1.

Data collection and processing
To develop an effective and robust model, it is necessary to establish the dataset using a strict criterion. Whalen et al. (20) identified human promoters using ENCODE Segway (21) and ChromHMM (22) annotations for GM12878, K562, HeLa-S3 and HUVEC cell lines. We obtain the sequences using BEDTools (23). Among these promoters, the length of most promoters is <3000-nt. Therefore, we remove promoters >3000-nt in length and remove the redundant sequences with a similarity of more than 80% in each cell line using the CD-HIT program (24). Finally, we treat them as positive samples of the benchmark dataset.
Different from previous studies in which negative samples are obtained from random genomic coordinates of nonpromoter regions (15)(16)(17)(18), this study uses the sequences with the highest similarity to the positive samples from the nonpromoter regions as negative samples. Specifically, for a positive sample, assuming its length is L, the following steps are adopted to generate negative samples: (i) Utilize a window of length L to slide 1000 times upstream and downstream of the positive sequence with a stride of 1-nt, and each time a sequence whose length equals to the length of the positive sample is obtained. (ii) If the sliding range covers the genomic coordinates of other positive samples, this sample will be removed from the positive dataset to avoid the generated negative samples overlapping with other positive samples, otherwise, it is difficult to define the true category of the generated samples. (iii) If the sliding range does not cover the genome coordinates of other positive samples, the similarity between these 2000 sequence segments and the positive sequence is calculated, and the sequence segment with the highest similarity is treated as the negative sample. This ensures that the generated negative samples and positive samples have a certain similarity so that the trained model is more robust. (iv) Remove sequences containing 'N' in positive samples and negative samples ('N' means that the base at this position is uncertain. Removing positive samples containing 'N' also removes the corresponding negative samples, and removing negative samples containing 'N' also removes the corresponding positive samples). Finally, to objectively evaluate the performance of the model, we use 20% of the sequences in each cell line as the independent test set. To ensure the reliability of the data, a positive sample and its corresponding negative sample are divided into a training set or test set simultaneously. The details of the human datasets are shown in Figure 2.

Feature encoding schemes
In this study, we employ six sequence encoding schemes to encode sequences, including reverse complement k-mer (RCKmer), mismatch k-mer (Mismatch), the composition of k-spaced nucleic acid pairs (CKSNAP), trinucleotide physicochemical properties (TPCP), pseudo trinucleotide composition (PseTNC) and Word2vec. Multiple feature encoding schemes comprehensively extract sequence information, and fusing features can effectively improve the representation ability of features. These encoding schemes are elaborated in the following sections.
Reverse complement k-mer. The RCKmer is a variant of k-mer, in which the k-mer is not strand-specific, so reverse complements are collapsed into a single feature (25,26). For instance, 'AGGT' is the reverse complement with 'ACCT'; 'TAG' is the reverse complement with 'CTA', and 'TCAGA' is the reverse complement with 'TCTGA'. In this study, we set k = 5 and thus the dimension of the RCKmer-based feature is 512.
Mismatch k-mer. The Mismatch is also a variant of k-mer, in which an error is allowed (26). For instance, when we calculate the occurrence number of 'GCA', we place 'GCA', 'ACA', 'CCA', 'TCA', 'GAA', 'GGA', 'GTA', 'GCC', 'GCG' and 'GCT'. The sum of the occurrence number of these 3-tuples is treated as the occurrence number of 'GCA'. In this study, we set k = 5, and thus the dimension of the Mismatch-based feature is 1024.
Composition of k-spaced nucleic acid pairs. The CKSNAP is used to calculate the frequency of all possible nucleotide pairs separated by any k space or less (26,27). For instance, in the sequence of 'ANNNTNNC', 'AT' is a 3-space nucleotide pair and 'TC' is a 2-space nucleotide pair. The CKSNAP-based feature is defined as follows: where N k (TT) represents the occurrence number of k-space nucleotide pair 'TT' and L is the length of the sequence.
In this study, we set k = 5 and thus the dimension of the CKSNAP-based feature is 96.
Trinucleotide physicochemical properties. The TPCP encoding has been successfully applied to DNA N4-Methylcytosine site prediction (28), which is defined as follows: where u i is the ith physicochemical property of the trinucleotide and f NNN denotes the normalized frequency of  trinucleotide NNN. In this study, we utilize 12 physicochemical properties (Supplementary Table S1), and thus the dimension of the TPCP-based feature is 768.
Pseudo trinucleotide composition. The PseTNC considers both the local sequence-order information and long-range sequence-order effects (26,29), which are defined as follows: where represents the number of the total counted ranks of the correlations along the sequence, f u represents the normalized frequency of trinucleotide, is the weighted factor and j is calculated as follows: where L represents the length of the sequence, P k (R i R i+1 ) is the value of kth physicochemical indices of R i R i+1 at position i, and μ is the number of physicochemical indices. In this study, we set λ = 2 and utilize six physicochemical indices to predict promoters, including angular parameters (twist, tilt and roll) and translational parameters (shift, slide and rise). Therefore, the dimension of the PseTNC-based feature is 66.
Word2vec. Word embedding techniques have achieved a great success in natural language processing (NLP) applications. Recently, word embedding techniques have been widely used in the bioinformatics community to tackle the limitation that the kmer-based features of different sequences may be very similar despite their orders being re-Nucleic Acids Research, 2022, Vol. 50, No. 18 10281 versed (30)(31)(32). In this study, we divide the sequence into 'word' to keep the order information of sequence. Specifically, a DNA sequence described as follows: where N 1 denotes the nucleotide at the first position and L is the length of the sequence. We take k consecutive nucleotides as a 'word'. Thus, the ith 'word' in the sequence is described as: After the above process, the DNA sequence is divided into sentences containing multiple words. Then we utilize the DNA sequences in the promoter dataset to form the corpus (each DNA sequence is the sentence in the corpus) and the 4 k types of k-mer segments to form vocabulary (each k-mer segment is a word in the vocabulary). To guarantee complete independence of the independent test set, we only utilize sequences in the training set to form the corpus and vocabulary. We train the language model by using the continuous skip-gram (Skip-gram) model in word2vec and obtain the feature vector for each 'word'. Then we concatenate the feature vectors of all 'word' in the sequence and treat them as the features of the sequence. Suppose each 'word' is embedded as a feature vector of dimension D, the feature dimension of each sequence is D × (Lk + 1). Therefore, if the feature dimensions obtained by sequences of different lengths in our dataset are not equal, the features cannot be trained in a deep learning model. To solve this problem, we transform the features of each sequence into equal dimensions using an adaptive pooling layer. In this study, each word has a length of 5 and is embedded as an 8-dimension feature vector. However, the sequence lengths in human datasets are unequal, which makes it impossible to train in deep learning models. Thus we additionally use adaptive pooling operations to further perform feature extraction. Given that 1000 is a centered number relative to the number of words in human promoter sequences, we set the output dimension of the adaptive pooling layer to 1000 × 8 ('1000' indicates 1000 words and '8' indicates that the feature dimension of each word is 8). Therefore, the dimension of the Word2vec-based feature is 8000.

iPro-WAEL architecture
iPro-WAEL is a weighted average ensemble learning model that integrates a random forests (RF) model and a CNN model, in which the RF model and CNN model use different features. The framework of iPro-WAEL is shown in Figure 3, and the RF model, the CNN model and the weighted average algorithm are depicted in the following section.
Random forest. Random forest (RF) is a combination of tree predictors, which is widely used in the bioinformatics community and shows excellent performance using traditional sequence-based features for prediction (15,17,33). Therefore, in iPro-WAEL, the RF model is trained using sequence-based features, including RCKmer, Mismatch, CKSNAP, TPCP and PseTNC.
Convolutional neural network-based model. CNN has been widely used in sequence-based prediction by combining word embedding techniques (34)(35)(36)(37). The convolutional layer performs convolution calculation using the convolution kernel to extract different features of the input, and more complex features can be extracted by stacking multiple layers. Besides, the convolutional layer is usually followed by a pooling layer to reduce the number of parameters in the network. In iPro-WAEL, the CNN model is trained using the Word2vec-based features and starts with three combinations of a convolutional layer and a maxpooling layer. We utilize 'relu' as the activation function in each CNN layer to enhance the nonlinear characteristics of the neural network. Then we add a dropout layer to avoid over-fitting, followed by two dense layers, which contain 64 nodes and 1 node, respectively, and utilize 'relu' and 'sigmoid' as activation functions, respectively. Finally, if the predicted value exceeds 0.5, the predicted result will be a positive example; otherwise, it will be a negative example. Besides, to obtain the model with better generalization, we divide 1/7 of the training set into the validation set and introduce an 'early stopping' mechanism, which is used to stop training early to avoid overfitting when the loss value of the model in the validation set does not decrease for five consecutive epochs.
Weighted averaging algorithm. The weighted averaging method obtains the combined output by averaging the outputs of individual models with different weights implying different importance (38). Compared to the simple averaging method, it can give more weight to the model with good performance. Specifically, the final output of iPro-WAEL is calculated as follows: where w i is the weight of the ith model, h i (x) is the output of the ith model and the weights are assumed to be constrained by To obtain the most suitable weights and complete independence of the independent test sets, we further divide the original training set into the training set and the dataset used for obtaining weights (weighted set) according to the ratio of 7:1. Therefore, the ratio of the training set, weighted set and test set is 7:1:2 in the RF model, and the ratio of the training set, weighted set, validation set and test set is 6:1:1:2 in the CNN model. Then we obtain the optimal weights by minimizing the loss values in the weighted set. Specifically, the loss value in the weighted set is calculated by the cross-entropy loss of the predicted values and the true labels, which is defined as follows: where y represents the true labels, and p represents the predicted values. To minimize the loss value, we utilize Sequential Least Squares Programming (SLSQP) algorithm from  the 'scipy' library, which uses the quasi-Newton method with a BFGS update of the B-matrix and an L1-test function in the step-length algorithm. Finally, the weights obtained with the minimum loss value are the optimal weights.

Performance evaluation of different features used for RF model
We consider more than 50 sequence features provided by IlearnPlus (26), BioSeq-Analysis2.0 (39) and BioSeq-BLM (40), but many of these feature extraction methods are not applicable for the unequal-length sequences. Therefore, we extract the applicable 29 sequence-based features of the training set in the GM12878 cell line and perform ten-fold cross-validation to evaluate the performance of these features using Area Under the Curve (AUC), Accuracy and Matthews Correlation Coefficient (MCC) (Supplementary Method). The results are shown in Supplementary Table S2. It is worth mentioning that the original one-hot encoding, DBE, TF-IDF and LDA-TFIDF cannot be used for encoding unequal-length sequences, but the fixed-length option provided by BioSeq-BLM makes it possible. Unexpectedly, however, although the good performance was achieved with fixed-length in the previous study (34), the performance of features with fixed-length in this study is poor. The first reason is that the lengths of the sequences in this study vary greatly (from tens of bp to 3000 bp), resulting in incalculable information loss caused by truncation and padding. Another reason is that the negative samples in this study have a certain similarity with the positive samples, and if the length is fixed, the similarity between the positive samples and the negative samples will be further enlarged, resulting in poor performance.
Among other 25 methods that can directly extract features of unequal-length sequences, many features have competitive performance in predicting promoters, especially Mismatch, RCKmer, CKSNAP, TPCP and PseTNC, showing that they are informative in predicting promoters. To further demonstrate the feature representation ability to fuse these five features, we visualize the distribution of the samples encoded by fused features in four cell lines using tdistributed Stochastic Neighbor Embedding (t-SNE) (41). We find that the positive samples and negative samples in the four cell lines are distributed in two clusters (Figure 4), proving that the fused features have excellent performance in distinguishing promoters and non-promoters. Therefore, we fuse and normalize these five features and use them as the input of the RF model.

Parameter optimization
To obtain the best-performing model, we evaluate the performance of the RF model and CNN model with different parameters, including the number of trees, the learning rate, the number of kernels and the size of kernels. We determine optimal parameters using the grid search algorithm on the dataset of the GM12878 cell line. To ensure that the performance of the model is only affected by changes in parameters, we set a fixed random seed for the RF model and the CNN model. The performance of the RF model and CNN model with different parameters are shown in Supplementary Tables S3 and S4, respectively. It can be seen that the performance of models is affected by the parameters. The RF model with a tree number of 300 performs best, and the CNN model with a learning rate of 0.001, a kernel number of 32 and a kernel size of 11 performs best. Therefore, we utilize these parameters to construct the RF model and CNN model.

Ten-fold cross-validation on the training datasets
Previous studies have shown that random division of the test set may lead to a large difference in performance between the training and test sets (42,43). To comprehensively evaluate the performance of iPro-WAEL, we perform tenfold cross-validation on the training set of four cell lines. The results are shown in Supplementary Table S5. It can be seen that iPro-WAEL achieves extremely high performance in predicting human promoters. It is worth mentioning that although the parameters are determined by param- eter optimization on the GM12878 cell line, the model with the same parameters also has excellent performance on the other three cell lines. These results indicate that iPro-WAEL has splendid performance and robustness in predicting promoters.

Performance evaluation on the independent test set and comparison with the state-of-the-art predictors
To comprehensively evaluate the performance of the models, we evaluate and compare the performance of iPro-WAEL with other eight computational methods, including iPro70-PseZNC, vw Z-curve, iPro54-PseKNC, iPromoter-2L, CNNProm, Promotech, DeePromoter and iPSW(2L)-PseKNC. We apply them to seven species (13 datasets), including humans (four cell lines), four prokaryotes (E. coli, B. amyloliquefaciens, B. subtilis and R. capsulatus). According to the evaluation method in the previous study (3), E. coli contains E. coli-general and E. coli-sigma70 and two eukaryotes (Mouse and Arabidopsis, each eukaryote contains a TATA dataset and a non-TATA dataset). Given that some of these datasets are imbalanced, we randomly subsample from the class with more data to construct a balanced dataset. The details of other species datasets are shown in Supplementary Table S6. All datasets are divided into the training set and test set according to the ratio of the human promoter dataset, and all methods utilize the same training set to train these models and evaluate their performance on the same test set. Due to the limitations of some methods, these methods cannot be applied to some datasets (Supplementary Table S7), and the reasons are shown in Supplementary Information. Also of note, to demonstrate and compare the generalization of the models, all methods including iPro-WAEL keep the parameters unchanged on the datasets of different species.
We evaluate the performance of nine models on 13 datasets and calculate the average performance of the models on the species with multiple types of datasets. The results are shown in Figure 5 and the detailed results are shown in Supplementary Tables S8-S14. It can be seen that the excellent performance of some models in one species is not transferable. For instance, as far as ACC is concerned, Z- curve slightly underperforms (about 1%) our model on the human dataset, but its performance on the mouse dataset is extremely different from our model, about 10% lower than that of our model. These results indicate that some methods are limited in their ability to identify promoters in multiple species. Fortunately, the iPro-WAEL proposed in this study has satisfactory performance on multiple species. The line representing iPro-WAEL in Figure 5 is always in the highest position, showing the effectiveness of iPro-WAEL on multiple species. Specifically, iPro-WAEL achieves optimal performance on six of seven species and is <0.1% lower than the best model on the remaining one (Arabidopsis). Additionally, although the second best performing model varies to a large extent on different datasets, ipro-WAEL outperforms these methods by a comparable or significant advantage.
To intuitively compare the performance of iPro-WAEL and the methods proposed by previous studies, we compare iPro-WAEL and these methods individually based on the maximum intersection of the dataset to which these models could be applied. For instance, only Human, Arabidopsis and Mouse are applicable to both DeePromoter and iPro-WAEL. Therefore, we compare the average performance of the two models on these three species. For another example, all species are applicable to iPromoter-2L, CNNProm, Promotech, iPSW(2L)-PseKNC, Z-curve and iPro-WAEL, so we compare the average performance of these models on all species. Based on this, we divide all datasets into three subsets according to their applicable conditions and compare iPro-WAEL and the methods proposed in previous studies.
It can be seen from Supplementary Table S15 that the performance of iPro-WAEL is significantly better than other methods. As far as ACC is concerned, iPro-WAEL obtains 1.7-18.6% higher than the second best performing model on the three subsets. However, it cannot be ignored that DeePromoter, which has the smallest difference from iPro-WAEL in performance, has very strict conditions of use. Due to its model architecture, DeePromoter has a very high limit on sequence length. Therefore, DeePromoter is difficult to be applied to prokaryotic promoter detection (the length of datasets in most prokaryotic is 81 bp or less, but DeePromoter can only predict sequences over 140 bp), and thus they also only report the performance on eukaryotic in their paper. Among other promoter predictors that can apply to both prokaryotes and eukaryotes, iPromoter-2L and iPro-WAEL have the smallest gap. But even so, iPro-WAEL still obtains 2.3% higher than iPromoter-2L. Overall, these results demonstrate that iPro-WAEL is an effective tool in predicting promoters and can be widely used to predict promoters in different species with satisfactory performance.

Cross cell lines validation
To explore the potential relationship between promoters in different cell lines, we further employ iPro-WAEL in the cross-cell lines and compare its performance with Z-curve, which is the best-performing method on Human datasets in the previous studies. We train the model on the training set of one cell line and predict the test set of the other three cell lines. The results are shown in Supplementary Tables S16 and S17. It can be seen that the performance of iPro-WAEL on cross-cell lines validation outperforms the bestperforming method (Z-curve) in previous studies and is not significantly different from that on the same cell line. It is noted that the data used to obtain the weights is from the training set, which ensures the independence of the test set when conducting cross-cell lines validation. These results indicate that iPro-WAEL is effective and robust for effectively predicting promoters in different cell lines and indicate potential similarity in sequence structure between promoters in different cell lines.

Identification of promoters and enhancers
Recently, some studies have shown that sequence architectures of enhancers and promoters are remarkably similar (44)(45)(46)(47)(48). Therefore, it is natural to ask whether iPro-WAEL can distinguish between promoters and enhancers. Whalen et al. (20) identified enhancers using the same method as promoters, most of which are less than 1000-nt in length. Therefore, we obtain the sequence using BEDTools (23) and remove enhancers over 1000-nt in length. Then we remove the redundant sequences with >80% similarity in each cell line using the CD-hit program (24). Given that the number of enhancers is much larger than that of promoters, we randomly subsample from the enhancer sequences to obtain the same number of sequences as the promoter sequences to construct a balanced dataset. The enhancer sequences are used as negative samples, and the promoter sequences are used as positive samples. The details are shown in Supplementary Table S18. Similarly, we use 20% of the sequences in each cell line as the independent test set. We keep the parameters of iPro-WAEL unchanged and the results of identifying promoters and enhancers are shown in Supplementary Table S19. It can be seen that iPro-WAEL achieves exceptionally high performance in distinguishing between promoters and enhancers, even though sequence architectures of enhancers and promoters are remarkably similar. These results demonstrate the superior robustness of iPro-WAEL and some sequence-level differences in the sequence architectures of enhancers and promoters.

Prediction ability on the promoters annotated in another method
Identification of promoters annotated with different methods may vary. Therefore, we further verify whether iPro-WAEL can predict promoters annotated by another method. Whalen et al. (20) identified promoters using Roadmap Epigenomics ChromHMM annotations for IMR90 and NHEK cell lines. We utilize the same method (see Method) to obtain negative sequences for these two cell lines and utilize 20% of the sequences in each cell line as the independent test set. Finally, the number of training sets for IMR90 and NHEK cell lines is 8192 and 8170, respectively, and the number of test sets for IMR90 and NHEK cell lines is 2048 and 2044 respectively. Similarly, we keep the parameters of iPro-WAEL unchanged. Then we evaluate the performance of iPro-WAEL on these two cell lines. It can be seen from Supplementary Table S20 that the performance of iPro-WAEL in predicting promoters annotated by different methods slightly decreased but is still satisfactory, which demonstrates that iPro-WAEL can capture potential features of promoters annotated by different methods and is of high capacity of generalization.

Estimate of important TFBS motifs in human promoter regions
Transcription factors play an important role in gene transcription through direct binding to their motifs in the genome and are involved in a large number of human diseases (49,50). To identify the important TFBS motifs that tend to bind in promoter regions, we propose a method to calculate the importance scores of the TFBS motifs. Specifically, given that the RCKmer and Mismatch used in the RF model can effectively reflect the importance of motifs, we integrate three methods to calculate the importance scores of the subsequence segment by interpreting these two features, including permutation-based importance, Gini importance and SHAP value. To integrate the three methods, we normalize and sum the importance scores of these two features obtained by each method. Then we perform a summation of importance scores representing the same subsequence segments, which is calculated as follows: score segment = score Mi smatch segment + score RC Kmer segment (11) where score segment represents the importance score of the segment, score Mi smatch segment represents the importance score of the Mismatch-based feature corresponding to the 'segment' and score RC Kmer segment represents the importance score of the RCKmer-based feature corresponding to the 'segment'. For instance, the importance score of segment 'AAAAA' is calculated by the sum of the importance scores of the segment 'AAAAA' in the Mismatch-based feature and RCKmerbased feature. It is noted that reverse complements are collapsed into a single feature. Therefore, the RCKmerbased feature importance of a segment corresponds to the RCKmer-based feature importance of its reverse complementary segment. For example, the RCKmer-based importance score of 'TTTTT' corresponds to the RCKmer-based importance score of 'AAAAA'.
Given that the segments are relatively short and thus may match multiple motifs at the same time, for each segment, we obtain the three motifs with the highest matching score with the motif position weight matrices (PWMs) from the HOCOMOCO Human v11 database (51) and assign them the same importance score. If multiple subsequences match the same motif, the importance scores of these subsequences are summed as the importance score of the motif. Then we categorize the motifs into different importance levels based on quantiles spaced at 20% of the importance score. We display the most important motifs (top 20%) in Figure 6. It can be seen from Figure 6 that the distribution of motifs' importance score in different cell lines is remarkably consistent, which can explain the reason for the high performance of iPro-WAEL on cross-cell lines validation and applicability of iPro-WAEL in predicting the promoters annotated in another method.
Besides, we find that many estimated important motifs are highly consistent with the previous studies. Among the most important motifs identified by iPro-WAEL, both the CTCF motif and CTCFL motif play a significant role in mediating chromatin loops (52,53); YY1 contributes to enhancer-promoter structural interactions and forms dimers to promote promoter-involved DNA interactions (54,55); ZNF143 is directly involved in chromatin looping (20) and provides sequence specificity to cohesin-associated chromatin loops at promoters (56), where ZNF143 directly binds to promoters and facilitates chromatin inter-actions that link promoters to distal regulatory elements (57); CUX1 is significantly enriched at promoters and contributes to predicting promoter-enhancer interactions (20); SP1 regulates chromatin looping between enhancer and distal promoter, and restores DNA looping and transcriptional activity by inserting at the promoter (58,59); REST is a transcriptional repressor with high predictive importance, as pointed out in TargetFinder (20); TFAP2C is a pioneer factor and is identified by TargetFinder as the most pre-dictive feature for predicting enhancer-promoter interactions (20). In addition, some next most important (top20%-top40%) motifs are consistent with some critical but underestimated motifs in promoter-involved interactions, as pointed out in TargetFinder (20), including SRF, MAX, TBP, CEBPB and RUNX. These results indicate that the formation of promoters is a complicated procedure and tends to be affected by a variety of factors.

DISCUSSION AND CONCLUSION
Identification of promoters can be achieved by some annotation methods such as genome-wide mapping of histone modifications (11), ENCODE Segway (21) and ChromHMM (22), which are time-consuming. In addition, it remains unclear what the potentially important TFBS motifs exist in human promoters. In this study, we propose a comprehensive and robust framework, named iPro-WAEL, to answer this question. We integrate two different models in iPro-WAEL, including the RF model, which utilizes traditional sequence-based features, and the CNN model, which utilizes word embedding techniques to extract the sequence features. From the analysis of a series of examinations on multiple datasets, it is concluded that iPro-WAEL has optimal performance and robustness, and it is superior to the existing methods. Besides, some previous methods have limitations in applying certain datasets. Fortunately, the iPro-WAEL proposed in this study does not suffer from this limitation. It is worth mentioning that iPromoter-2L noticed this problem, thus their webserver splits the sequence into multiple 81-nt segments when predicting sequences with length over 81-nt, but still cannot make predictions when the sequence length is less than 81-nt. Overall, iPro-WAEL is a competitive and robust model for predicting promoters in multiple species and can be widely applied to datasets with different sequence lengths.
In addition, we estimate potentially important TFBS motifs in promoter regions using the computational method we designed. We identify many critical sequence motifs for predicting promoters, including ZNF143, YY1, CTCF, CUX1, SP1, REST, TFAP2C, SRF, MAX, TBP, CEBPB and RUNX, which are consistent with the previous study. However, to the best of our knowledge, several potentially important motifs identified by iPro-WAEL have received little or no attention. Therefore, these motifs may be involved in the formation of promoters and gene regulation through underappreciated or potentially new biological interactions with certain proteins, and may be the key to analyzing the difference and association of promoters in different cell lines of humans.
There are also some areas in which our method can be improved. For instance, we identify the importance of motifs by the importance of features reflecting the frequency of 5mer segments, but indeed, the Word2vec-based features also reflect the infor-mation of the 5-mer segment because its window size is 5. However, we cannot interpret Word2vecbased features due to the difficulty of interpreting the embedding space intuitively. Therefore, we do not integrate word2vec information when calculating the 5-mer segment importance scores. Indeed, the field of NLP also faces this dilemma, although word embedding techniques are widely used in the NLP field. Although the RF model has better performance and higher weights than the CNN model on human cells, which indicates that interpreting the features used by the RF model may have higher confidence, our experimental results also show that integrating the RF model and CNN model can further improve the performance in predicting promoters. Therefore, if the importance scores of the Word2vec-based features can be calculated reasonably, the reliability of the method for analyzing the importance of motifs may be further improved when being combined with the importance scores of the features used in this study. Additionally, we have demonstrated a satisfactory performance when using iPro-WAEL trained on one cell line to predict promoters in another cell line, indicating that promoters from different cell lines may be similar at the motif level and related in some way. However, the performance of using models trained on prokaryotes to predict the promoters of eukaryotes is unsatisfactory, indicating that promoters of different species may be associated with different motifs. Overall, it will be interesting to explore the differences and associations between promoters of different cell lines and different species by analyzing the importance of motifs, which will provide a novel way to analyze cell specificity and species specificity.