Transcription factor motif quality assessment requires systematic comparative analysis

Transcription factor (TF) binding site prediction remains a challenge in gene regulatory research due to degeneracy and potential variability in binding sites in the genome. Dozens of algorithms designed to learn binding models (motifs) have generated many motifs available in research papers with a subset making it to databases like JASPAR, UniPROBE and Transfac. The presence of many versions of motifs from the various databases for a single TF and the lack of a standardized assessment technique makes it difficult for biologists to make an appropriate choice of binding model and for algorithm developers to benchmark, test and improve on their models. In this study, we review and evaluate the approaches in use, highlight differences and demonstrate the difficulty of defining a standardized motif assessment approach. We review scoring functions, motif length, test data and the type of performance metrics used in prior studies as some of the factors that influence the outcome of a motif assessment. We show that the scoring functions and statistics used in motif assessment influence ranking of motifs in a TF-specific manner. We also show that TF binding specificity can vary by source of genomic binding data. We also demonstrate that information content of a motif is not in isolation a measure of motif quality but is influenced by TF binding behaviour. We conclude that there is a need for an easy-to-use tool that presents all available evidence for a comparative analysis.

We thank the reviewers for their comments. The paper has been updated in response to their comments as follows: 1. A sentence is added in the background section to update information on the performance of PWM models.
2. Most of the figures have been updated or modified, especially the captions. Figure 6 was added while Figure 5 was replaced.
3. We include analysis on the effect of negative sequences on motif ranking. We added a subsection, 'How choice of negative (background) sequence affects motif ranking' and Figure 6.
4. We include data on a re-run of the analysis using PBM data. A paragraph has been included in the Data subsection (2nd paragraph), a new subsection, 'Effect of PBM data on motif assessment' in the results section discussed in the discussion section.

5.
A subsection of the results section replaced with 'Effect of statistics on motif ranking'. This goes with the new Figure 5. In addition, Table 2, which offered supporting information to the old Figure 5, has also been removed.
6. The notations used in the formulas have been harmonized. 9. We have added a list of ENCODE ChIP-seq data used to the repository as Supplementary Table 3 (Table S3).
10. Finally, we have made other minor changes in response to the reviewers comments, as noted in our detailed response to the reviews.

Background
Understanding gene regulation remains a long-standing problem in biological research. The main players, transcription factors (TFs), are proteins that bind to short and potentially degenerate sequence patterns (motifs) at gene regulatory sites to promote or repress expression of target genes. The search for a code to predict binding sites and model binding affinity of TFs has led to several experimental techniques and motif discovery algorithms being developed ( Figure 1).
A position weight matrix (PWM) is the common form of representing TF binding specificity. For a motif of length L, the corresponding PWM is a 4 × L matrix of probabilities of observing a base b (A, C, G or T) at position i through L. Other forms of representing TF binding specificity have been introduced [1][2][3][4] , but Weirauch et al.
showed that a well-trained PWM performs comparably to some of the above well trained complex models 5 . However, recent studies 6-8 have reported significant improvement to the PWM by models that consider nucleotide inter-dependencies. The persistent popularity of PWM can be attributed to its simplicity and ease of use as well as the ease of visualizing a PWM using a sequence logo 9 . Motifs can be found using a variety of methods including algorithms that do de novo motif discovery from sequences containing binding sites [10][11][12] and in vitro methods such as protein binding microarrays (PBM) 13 and high-throughput systematic evolution of ligands by exponential enrichment (HT-SELEX) 14 .
Initially, the low resolution of the available experimental techniques for TF binding specificity detection was a hindrance to the quality of binding models. However, next generation sequencing and techniques like chromatin immunoprecipitation (ChIP) followed by deep sequencing (ChIP-seq) 15 and exonuclease cleavage in ChIP-exo 16 that measure TF in vivo occupancy, have improved the resolution to single-nucleotide level. In addition to providing high resolution data for motif discovery, they are a useful resource to test the quality of the available motifs since they are TF specific. However, no benchmark capable of assessing the growing range of published motifs is available and quality measures are largely subjective 17 .
Although it is possible that PWM models' ability to describe TF binding may be getting saturated, the lack of a robust approach to test the quality of a model and maximize the best-performing ones may also be hampering improvement in performance. How are the algorithms being developed, tested and improved? Furthermore, the number of motif finding algorithms from dissimilar data sets and subsequently the number of motif models for a single TF generated, continue to increase. There are at least 44 PWM motif models available in 14 different databases for Hnf4a alone. How does the end-user decide which motif to use? In this study, we review and test the approaches used to evaluate TF binding models.

Review of motif assessment approaches
The available motif assessment techniques can be divided into three categories: assess by binding site prediction, motif comparison or, by sequence scoring and classification.
Binding site prediction Early review and assessment of motif-finding algorithms tested tools on the ability to predict sites of motifs, known or inserted into the sequence. Tompa et al. tested motif discovery algorithms by their ability to predict sites of inserted motifs using statistical measures for site sensitivity and correlation coefficient 18 . In this first comprehensive study, they found that a motif assessment problem is complex and admitted that inserting random motifs fails to capture the biological condition of TF binding. Later, Hu et al. 19 used real RegulonDB binding data in a large-scale analysis of five motif-finding algorithms. The tools available at that time performed poorly -"15-25% accuracy at the nucleotide level and 25-35% at the binding site level for sequences of 400 nt long" -largely due to the poor quality of RegulonDB annotations 20 .
Sandve and colleagues 21-23 tested motif discovery algorithms using sequences with real and inserted binding sites as benchmarks; from Transfac, and the third-order Markov model respectively. Quest and colleagues 24 developed the Motif Tool Assessment Platform (MTAP) as an automated test of motif discovery tools. However, this was computationally expensive and was made obsolete by new experimental data and algorithms.
The most comprehensive assessment based on binding site prediction so far has been by the Regulatory Sequence Analysis Tools (RSAT) consortium. In their 'matrix quality' script, they use theoretical -information content (IC) and E-values -and empirical scores computed by predicting binding sites in RegulonDB, ChIP-chip and ChIP-seq positive and negative control sequences 20 .
Inadequate knowledge of TF binding sites has mainly hampered the ability to assess motifs and algorithms by binding site prediction. Predicting binding sites that are inserted or known in the sequences cannot accurately identify unknown true sites. Techniques that identify such sites may be penalized. Until TF binding sites are well annotated, this technique cannot be confidently utilized.

Motif comparison
Novel motifs can be assessed by comparison to 'reference motifs' using the sum of square deviation, Euclidean distance and other statistics that measure divergence between two PWMs 25,26 . Thomas-Chollier et al. proposed a motif comparison approach for their RSAT algorithm where they combine multiple metrics, including Pearson's correlation, width normalized correlation, logo dot product, correlation of IC, normalized Sandelin-Wasserman, sum of squared distances and normalized Euclidean similarity for each matrix pair 27 . They then unified all of these scores to ranks whereby the mean of the ranks is considered the overall score.
Assessing motifs by comparison, as currently implemented, only tests similarity to the available motifs with little information on quality and ranks of the motifs. It assumes accuracy of 'reference motifs', with no way of assessing novel ones. In addition, the definition of 'reference motifs' remains largely subjective.

Assessment by scoring
Motif assessment has since shifted towards scoring positive sequences known to contain binding sites and negative background sequences without binding sites, driven by high-throughput sequencing techniques 5,28-30 . This avoids the need to identify binding sites a priori by focusing on the ability to classify the two sets of sequences. The differences in the assessments arise from the choice of sequences to use as positive and negative, the thresholds used to identify binding sites, the length of the sequences in both sets, the scoring function and the statistic used to quantify the performance of the tool.
For ChIP-seq data, the main difference is that the length of sequences (250bp 28 ; 600bp 30 , 100bp 5 or 60bp 31 ) and the choice of negative sets (300bp downstream; 28,30 ; random sequences, 5000bp from a transcription start site (TSS) or random genomic sequences 5 , or flanking sequences 31 ) differ greatly in sequence scoring. In addition Agius et al. 31 , test PWMs and support vector regression (SVR) models in the 36bp sliding window of the test sequences, a deviation from the rest of the techniques. All these differences, in addition to the scoring functions and statistics used, can lead to variations in the results of comparative analyses. Users and algorithm developers therefore have to frequently re-invent the wheel to test their tools. Figure 1 shows the evolution of experimental motif discovery assessment techniques. We have not focused on the experimental techniques or motif discovery algorithms as excellent reviews are already available 17,32 . Rather, we focus on TF binding models represented as a PWM and aim to determine how the choice and length of benchmark sequences, scoring functions, and the statistics influence motif assessment. We hope that this study will highlight some of the pitfalls in the previous motif assessments and provide a starting point for a standard in motif assessment that will ensure comparability and reuse of results. We used motifs from a number of databases and publications listed in Table 1. The TFs used in this analysis were selected based on availability of ChIP-seq data with motifs in at least 10 motif databases. We converted these motifs from their various formats into MEME format and scored the positive and negative sequences with GOMER, occupancy, energy and log-odds scoring functions.

Data
We quantify how each motif performs using AUC, MNCP, Spearman's and Pearson correlation ( Figure 2). This was implemented in a Python module which is available free from https:// github.com/kipkurui/Kibet-F1000Research. This repository also contains raw data and Ipython notebooks that document how to reproduce the analysis we describe in this paper.  Scoring functions When testing motifs by scoring ChIP-seq or PBM data, multiple scoring functions are available, which may affect the outcome. In the section that follows, we describe the scoring functions tested, as well as provide a review of how they have been previously applied.

GOMER scoring
The GOMER scoring framework was introduced by Granek et al. 49 but adapted for PBM sequence scoring 37,38 . It seeks to compute the probability g(S,Θ) = exp(f (S,Θ)) that a TF, given PWM Θ, will bind to at least one of the sub-sequences of S. This assumes that each site can be bound independently where L is the length of sequence S, and S t:t+k is the sub-sequence of S from position t to t+k inclusive. See Chen et al. 37 for more details.

Occupancy score
The occupancy score calculates the occupancy of a PWM Θ for sub-sequence S i of length k as the product of the probabilities of each base in S i using equation 2.
For a sequence S of length L, the sum of the occupancies of all subsequences S i (sum occupancy) 28,50 , the maximum score (maximum occupancy) 30 , or the average occupancy (average motif affinity -AMA) have been used. Sum occupancy is defined in equation 3:

BEEML-PBM energy scoring
The energy scoring framework of binding energy estimation by maximum likelihood for protein binding microarrays (BEEML-PBM) 4 computes the logarithm of base frequencies with the idea that this is proportional to the energy contributions of the bases. The binding energy at each location is computed; the lower the binding energy, the higher the binding affinity. For each sequence, the sub-sequence with the lowest binding energy represents the score of the sequence. It has mainly been used to score PBM data 5,30 .
The probability that sub-sequence S i is bound is given by equation 4, where, for a sub-sequence S i , E(S i ) is given equation 5, 1 ( ) ( , ) ( , ), for binding site of length L, ∈(b, t) is the energy contribution of base b while S i (b, t) is an indicator function of site t within S i (1 with base b, 0 otherwise).

Log-odds scoring
In log-odds scoring, used by a majority of the MEME Suite tools 51 , the score for a given site is the sum of the log-odds ratios of a PWM at the match site. For a sub-sequence S i of length L scored using PWM Θ, the log-odds score is given by equation 6, , 1 ( , , ) ( , ) .
where p b the background probability (uniform background probability of 0.25 is used) and S i (b, t) is an indicator function of site t as in equation 5.
The score for a given sequence can then be derived by summing individual scores or by finding the maximum score. Sum log-odds scoring has generally been used by MEME Suite tools while maximum log-odds scoring has also been used to compare motifs represented differently (PWM, k-mer and SVM models) against one another 30,31 . Each of these approaches has inherent advantages but may produce inconsistent results.

Statistical measures of performance
With the scores of each motif for the sequences acquired, binding prediction can be evaluated by various statistics. The area under the receiver operating characteristic curve (AUC) 52 has been widely used, especially with the advent of PBM 5,28,37 . In addition to popularizing AUC, Clarke et al. 52 also introduced a novel metric, mean normalized conditional probability (MNCP), for quantifying the correlation between DNA features and gene regulation. This statistic has been applied for motif assessment in GIMME motifs 53 and is said to be less affected by the presence of false positives compared with AUC since it places emphasis on true positives. MNCP is a rank-based statistic that determines if mean occurrence of a motif in test sequences is higher than the mean occurrence in a random set. Each set of sequences is ranked based on the mean occurrence, and the MNCP calculated by finding the mean of the normalized ratio of the two sets of ranks. We use MNCP to test how it contributes to better prediction in an effort to encourage its use.
Pearson and Spearman's rank correlation are still widely used as a measure of motif performance. Spearman's rank correlation has been used for PBM and ChIP-seq sequences 28 while Pearson's correlation was used by Weirauch et al. 5 . However, Weirauch et al. cautioned on the use of Spearman's correlation for PBM data citing its inability to exclude low intensity probes. We wish to check the usefulness of correlation statistics in motif assessment.
In addition to comparing the scoring approaches, we use CentriMo version 4.10.0 in differential mode 54 -an option that tests differences in motif enrichment between two sequence sets -in a novel way for motif assessment. We set differential mode parameters for local rather than central enrichment of all the input motifs in the positive (primary) and negative (control) set, as described in the Data section, by using a very large threshold. The negative log of the E-value is used as the measure of a motif's enrichment and rank. Motif enrichment analysis has previously been performed 36 using the FIMO algorithm 55 to scan for motif matches in sequences and calculate an enrichment value.

Results
Length of sequences has a little effect on motif performance The size of the putative binding region -length of the sequences in each data set -is to some extent a proxy for how accurate the ChIP-seq experiment was. If the result was accurate a narrow region should contain the true site. For the three variants of sequence length, we did not observe a significant effect (p=0.113, for 50 and 100; p=0.0545, 50 and 250; p=0.678, 100 and 250bp -Wilcoxon rank-sum test) on the scoring of the sequences ( Figure 3). The scores assigned for each sequence length, however, seems to indicate how the TFs bind. Motifs with higher scores at lower sequence length (50 or 250bp) are generally enriched at the ChIP-seq peak, which is also a strong indicator of direct binding 56 . This is consistent with a previous observation that a successful ChIP-seq experiment localizes binding within about 100bp of the true site 57 . Others with significantly better AUC values at 250bp sequence length like Elf1 (p=0.017, Wilcoxon rank-sum test) and Sp1 (p=0.013, Wilcoxon rank-sum test) 58 , are known to bind cooperatively.
Tissue or cell line of the data could affect enrichment Transcription factors bind to their possible sites in a sequencespecific manner. Some actually have alternative binding motifs depending on the tissue or cell line. Unless the purpose of motif assessment is to identify tissue-specific binding, if data is available from more than one cell line, an average of the scores should be used. For example, in Figure 4, we show that the rank correlation of the motif scores in different cell lines can be as low as 0.8 for GOMER scoring (or as low as 0.65 when energy scoring is used), and not 1 or very close to 1 as would be expected if the cell line had no effect. In addition, FOXA1_1.GUERTIN motif is differentially enriched only in the A549 cell line (although this could be an outlier).
In light of this possible effect, the results displayed throughout this paper are based on the mean score of all the available ChIP-seq data sets to avoid a bias towards cell line-specific motifs.
How choice of negative (background) sequences affect motif ranking In motif discovery, the choice of background sequences has significant effects on the motifs identified. We sought, therefore, to test whether motif scoring would be affected in a similar way. In addition to downstream sequences, we used a dinucleotide shuffled set from the positive sequences. The scores obtained using dinucleotide shuffled positive sequences were always lower than those for downstream sequences. We then computed and plotted the rank correlation of scores normalized by maximum score for each TF, from which we find that it affects the ranks of the motifs ( Figure 5) in a TF-specific manner. However, the scores from the two sets of negative sequences used are not significantly different (p=0.484, Wilcoxon rank-sum test). For Myb, the low correlation could be attributed to how it binds, indirectly 59 .

Effect of statistic on motif ranking
The statistic used, whether it measures scores correlation or ability to classify the two sets of sequences, will definitely have an effect on how we interpret the results of the analysis. Generally, the motifs ranks based on AUC and MNCP statistics' are not significantly different (p=0.52, Wilcoxon rank-sum test), but the ranks based on  Pearson and Spearman's differ significantly from MNCP or AUC scores (p=0.006 and 0.002 respectively, Wilcoxon rank-sum test). The large standard deviations of the correlation statistics' scores, as shown by the error bars in the Figure 6, shows how unreliable the use of correlation statistics to rank the motifs can be. The correlation scores are also quite low.

Effect of scoring function is transcription factor specific
We tested the ability of PWM models to discriminate positive (top 500 peaks of width 100bp centred on the peak) and negative (500 peaks 100bp wide located 500bp downstream from the positive) sequence sets using five scoring functions. Maximum and sum log-odds scoring had low discriminative power for most of the motifs when AUC (Figure 7) and MNCP (Figure 8) statistical measures are used. However, sum log-odds scoring had some good performance (over 0.55 AUC scores) for some TF motifs like Max, Nrf1, Tcf3 and Pax5.
There is no significant difference in performance when GOMER, energy or occupancy scores (sum, maximum and AMA) are used for scoring ( Figure 7B) with AUC statistic (see Table_S1 for details of statistical significance). Also, we did not observe any significant   difference (p=0.85, Wilcoxon rank-sum test) between sum occupancy and maximum (Table 2), contrary to a claim by Orenstein et al. 28 . When using MNCP, there is a higher rank correlation among the scores assigned by the different scoring functions except logodds scoring ( Figure 8B). When using AUC or MNCP statistic, Ctcf, Egr1 and Hnf4a score significantly higher in energy while other TFs like Pou2f2 and Esrra, the preference is reversed. These observations are reflective of the inherent features of the scoring functions or the statistics used.

Motif length and information content
Motif length has little bearing on the quality of motif, independent of other factors. However, specific motifs with very high IC such as those from POUR have a lower performance (Figure 9). In the same light, those motifs with low IC also have a lower performance in vivo.
The heat map in Figure 9 shows how the motif scores from the four discriminative functions correlate with motif length, full-length  IC and average IC. The examples have no consistent correlation between the IC and the scores ( Figure 9A). However, there is a negative correlation between both the total IC and motif length, and the scores except for sum log-odds scoring, which has no significant correlation (p=0.34, correlation p-value). This shows that motif length, rather than the IC of the motifs, negatively influences the scores assigned by these functions. This is not a general rule. Some TFs exemplify a different scenario. For example, Egr1 ( Figure 9B) has a positive correlation between IC and scores and a negative correlation with motif length (except for sum log-odds scoring), showing that it has a highly specific binding site 60 . Mef2a, on the other hand, has a positive correlation between motif length and scores showing preference for longer low information motifs ( Figure 9C). This could also reflect variability in binding sites.
Ctcf has the highest negative correlation for average IC, with a neutral to positive correlation for motif length ( Figure 9D), which may indicate preference for longer low IC motifs.

Comparison of motif databases
We have shown that the effect of scoring algorithms is TF-specific.
We also test to see how, overall, the different databases (DBs) are ranked against each other. For TFs with more than one motif in a given DB, we use the best performing one to represent the DB. We also use motif enrichment-based assessment using CentriMo version 4.10.0 to provide more evidence to scoring based techniques' results. Motif enrichment analysis compares how various motifs in foreground sequences are enriched compared with background sequences. In comparing how two or more motifs for the same TF perform, the level of enrichment of the motif in sequences known to contain possible binding sites of the TF compared to some background should provide a measure of the quality of the motif.
Figure 10 provides a summary of ranking of the various databases for the given TFs. We observed that the performance of a majority of the motif databases did not differ much, except for TF2DNA motifs, but the ranking or the performance of individual motifs differs. This further supports the observation of TF-specific performance of scoring function, algorithms and DBs. It also shows that no single database currently outperforms the others for all TFs. There is agreement in ranking of the best (ZLAB and HOCOMOCO) and worst performing (TF2DNA and SWISSREGULON) DBs. We observe that, compared with GOMER ( Figure 10A), the ranks for most DBs remain the same when using energy ( Figure 10B) except for POUR and JOLMA. This shows that motifs from these DBs, or at least the best performing ones, may be favoured by energy scoring. It is also noteworthy that POUR and GUERTIN DB motifs, discovered and tested on ENCODE ChIP-seq data, generally perform poorly.

Effect of PBM data on motif assessment
To test whether the conclusions of the paper are only linked to ChIPseq data, we re-ran the whole analysis using PBM data from the UniPROBE database 13 . It is important to note that we only found 9 TFs that had PBM data from the set used in ChIP-seq analysis. Since this may bias the comparison, we compared with a similar set in ChIP-seq and found the observations below were not affected by the difference in number of TFs used. These observations include: 1. A much higher energy score in PBM ( Figure S1 and Figure S2) compared with ChIP-seq (Figure 7 and Figure 8). We also observe a much lower correlation between the energy and the occupancy scores.
2. A stronger negative correlation between the occupancy scores and motif length -0.47 compared with -0.28 of energy scoring ( Figure S3), an observation not made when using ChIP-seq data (Figure 9). This may actually explain observation 1.
3. Motifs generated using the PBM technique perform best when using occupancy scores with MNCP or energy scores with AUC or MNCP, except when occupancy scoring and AUC are used ( Figure S4). Poor ranking of UniPROBE PBM-derived motifs by GOMER-AUC may be linked to the fact that they penalize long motifs -UniPROBE motifs are know to be long (mostly over 14bp).
4. Energy scoring with either MNCP or AUC, or occupancy scoring with MNCP display similar behaviour: a preference for specific motifs, which may be longer or have a higher IC. This supported by the high negative correlation between motif length and occupancy scores with AUC ( Figure S3).
5. TF2DNA motifs perform better when PBM data is used ( Figure S4A) compared with ChIP-seq data ( Figure 10A), and especially so when GOMER scoring is used together with AUC statistic. It is not immediately clear what the cause of the difference of performance of TF2DNA motifs in PBM and ChIP-seq data is, but the short length (7bp) of TF2DNA motifs and the fact that they were generated in vitro could provide some explanation given that PBM data are generated in 8-mers and PBM is also an in vitro technique.

Discussion
We have described a comparative analysis on the effect of scoring functions, ChIP-seq test data processing and statistics on motif assessment. We chose to focus on TF binding models represented as a PWM, since it is most commonly used. The review reveals the complexity of the motif assessment problem, showing no appropriate solution is available so far. The available techniques focus on testing motif algorithms or the experimental techniques used, but little has been done to compare the available motifs for a given TF. There is a need for a tool, accessible and easy to use by end-users, to aid in choosing motifs.
The use of 100 or 250bp sequence length provides necessary discrimination for the TFs tested (Figure 3). The performance was also found to be TF specific, an observation that could reflect inherent binding behaviour; direct, indirect or cooperative binding of the TF. This supports the observation that direct binding can be inferred from ChIP-seq peaks 56 . We also confirm that 100bp provides acceptable specificity in motif assessment given that most TF binding sites are less than 30bp 57 .
Since TF binding is cell line specific 61 , users should be aware of bias caused by use of one cell line in an assessment. We observe that the use of more than one cell line reduces the bias towards cell line specific motifs (Figure 4).
The MNCP rank-order metric is similar to AUC but derived by plotting true positive hits against all sequences' scores. This places emphasis on true positives and therefore is less affected by false positives. Most of the observations from the PBM-based analysis support the conclusion that energy scoring prefers specific motifs (long or with a high IC). We also observe an agreement when energy scoring is used with AUC and MNCP, or occupancy scoring. In MNCP, the preference for specific motifs is expected because the MNCP score provides a rank-based measure of the ratio of mean occurrence of a motif in test sequences and a random set. These observations are not conclusive and further research may be required. Although there is no clear winner among the scoring function, occupancy-based (GOMER, AMA, sum and max) and energy scoring functions are preferred. We recommend, based on the presented evidence, using occupancy scoring with the MNCP statistic or energy scoring with normal AUC or the MNCP statistic.
There is no significant correlation (p=0.513, correlation p-value) between the IC and the motif scores ( Figure 9). This compares with the conflicting observations that the best-quality motifs may have low IC motifs 5 , or high IC motifs 62 . Weirauch et al. did not normalize for motif length, which results in high IC motifs that are generally longer but not necessarily more specific 5 . A shorter motif with higher IC per position will be more specific but have lower total IC. We argue that the effect of IC on motif quality is dependent on the TF binding behaviour. TFs with short and specific binding sites will favour high IC while those with long and variable binding sites will be more accurately modelled with low IC. Furthermore, it has been shown the low IC flanking motif sites contribute to specificity of binding in vivo 62 . We have also shown that the techniques used in motif assessment have a direct effect on motif discovery. We observe how motifs generated from similar data using the same techniques could have highly variable performance in POUR, ZLAB and GUERTIN motifs (Figure 10). This difference in quality can be explained by motif discovery algorithms used, data processing as well as the assessment techniques used in each motif discovery pipeline. POUR motifs are learned from full-length sequences of the top 250 peaks using five motif finding algorithms (MEME, MDscan, Trawler, AlignAce and Weeder) 41 , the ZLAB group used 100bp of the top 500 sequences centred on the ChIP-seq peaks using MEME-ChIP 63 , while GUERTIN reports the top 5 motifs for each technique generated using MEME. For quality assessment, POUR 41 used a TFM-PVALUE 64 to match motifs against the testing ChIPseq data set and the most enriched motifs against a background composed of intergenic non-repetitive regions. ZLAB group used FIMO 55 , which uses a log likelihood score for motif scanning.
The worst performing motifs, from TF2DNA, are generated from 3D models of TF from experimental or homology-modelled PDB files. When generating these models, they determined the accuracy of the models based on similarity to UniPROBE and JASPAR motifs at a given threshold. They claimed their technique successfully identifies true motifs 41-81% of the time depending on the quality of templates used in modelling 3D structures. We speculate that part of the reason for this low performance could be use of motif comparison against 'reference motifs' as a measure of motif quality, in addition to being in vitro derived. Better performance of TF2DNA motifs in PBM data ( Figure S4) further supports this observation. Although JASPAR and UniPROBE are widely used, reliance on reference motifs is problematic, especially with the advent of motif databases like HOCOMOCO and CIS-BP that have motifs with better prediction quality. As a general principle, it is not reasonable to use historical data as a benchmark for assessing potentially better new methods.
We also show that the choice of data used in motif assessment has a direct effect on the ranks of the motifs. It goes without saying that PBM-derived motifs will perform better when tested with PBM data or for ChIP-seq based motifs tested on ChIP-seq data. The main criteria for choosing the test data should be based on the intended use of the motifs. In addition, we confirm the effect of negative (background) sequences in motif assessment, an effect well known in motif discovery.
We have confirmed that motif assessment has transcription-specific variability, an observation previously made 65 . Assessments should be less focused on how a particular motif database or algorithm performs but on how different motifs, for a particular TF, perform compared to the other potential motifs. For the end user, no single database can provide the sole measure of quality of new data. This raises the need for collation of the different motifs tested using a variety of motif assessments to provide information to the end user on their ranks.

Conclusions
We have demonstrated that the scoring techniques used in motif assessment influence ranking of motifs in a TF-specific manner. Motif assessments and tools developed to date have focused on comparing algorithms, experimental techniques or databases. This does not help the user choose which motif to use for a given TF. Some TFs reviewed here have at least 44 PWM motifs available, raising the need for a tool that can be utilized to rank these motifs. We have also shown that data processing as well as motif assessment can have a significant influence on the quality of motifs derived. Therefore, the choice of an assessment approach should be given as much thought as that of the motif discovery algorithm to use. We have also shown the effect of IC on motif quality is influenced by TF binding behaviour.
In short, a single measure of motif quality is likely to remain elusive, pointing to the need for tools and methods for comparative analysis using multiple methods. Lessons learned from this analysis will be useful in a number of ways. Firstly, we are working on a web-based application that can allow users to compare motifs available in different databases for a specific TF. Secondly, we intend to extend the motif by comparison approach to avoid 'reference motifs' bias. Thirdly, we have shown the effect of motif scoring on motif discovery. We intend to use the robust motif assessment techniques we introduce to improve motif finding.

Data and software availability
Data, software, supplementary files and documentation for 'Transcription factor motif quality assessment requires systematic comparative analysis' are available from Github: https://github. com/kipkurui/Kibet-F1000Research.
Archived files at the time of publication are available from Zenodo: doi: 10.5281/zenodo.46440 66 Author contributions CK designed and performed the analysis and wrote the first draft. PM supervised the work and contributed to subsequent drafts. All authors read and approved the final manuscript.

Competing interests
No competing interests were disclosed.

Grant information
The financial assistance of the South African National Research Foundation (NRF) towards this research is hereby acknowledged. Opinions expressed and conclusions arrived at, are those of the authors and are not necessarily to be attributed to the NRF. PM funding: NRF/IFR Grant 85362; CK: DST Innovation Doctoral Scholarship.

Supplementary information
This section provides supporting figures for the paper. Figure S1. Effect of scoring function on motif ranking using AUC statistic for PBM data. A. For each transcription factor (TF), the mean AUC score is used to represent it for each scoring functions used. In B, we show how the ranks assigned to various motifs for a given TF by each scoring function are correlated. It displays the pairwise rank correlation for all TFs in A. Sumlog: Sum log-odds function, Sumoc: sum occupancy score and Maxoc: maximum occupancy. Figure S1 for details. Figure S3. Effect of motif length and IC on scoring functions using PBM data. In this figure, we show the correlation of motif length, full length information content (IC) and the assessment scores, to determine how performance of scoring functions are influenced by motif characteristics. For each motif, the information content is calculated based on information theory for the whole length and also normalized for length. The results for average motif affinity (AMA) and maximum occupancy are similar to sum occupancy, and are not included. Figure S4. Ranking of motif databases when based on PBM data. We compare the motif databases by using the best ranking for each motif using GOMER and energy AUC and MNCP values, and CentriMo enrichment values. For each scoring function, the scores for each TF are normalized by dividing each value with the maximum, which are then averaged to rank the different databases.

Figure S2. Effect of scoring function on motif ranking based on MNCP statistic in PBM data. See caption in
46. Contreras-Moreira B: 3D-footprint: a database for the structural analysis of F1000Research 1.

2.
3. I am happy to state that the authors thoroughly addressed most of my concerns regarding the previous version of the manuscript. However, I still have a few minor comments, which the authors might consider, as listed in the following.

Open Peer Review
In the introduction, the authors state that Chip-seq has improved "resolution to single-nucleotide level", whereas at the beginning of the Results section, they state that "a successful ChIP-seq experiment localizes binding within about 100bp of the true site", which still is quite accurate but not exactly single-nucleotide accuracy. Please clarify.
In section "Background", sub-section "Motif comparison", the authors list as possible measures of motif divergence "sum of square deviation, Euclidean distance". As the former is just the square of the latter, I would suggest to give another example (for instance, correlation-based meaures) in the introduction of that sub-section.
In section "Methods", sub-section "Data", the authors now describe the selection of negative sequences, which were "extracted 500bp downstream from the highest coordinate". a) However, as ChIP-seq peaks lack an orientation, I wonder which direction "downstream" refers to. Does this mean that the authors considered the forward strand of the available genomic sequence, or do they also take, e.g., the orientation of closely located genes into account?
b) In the response, the authors explain that "our focus was to get negative sequences [...] which maintains the nucleotide composition". As I mentioned in my previous review, sequences located 500bp downstream of a transcription factor binding site might already be located in coding sequence, which affects nucleotide composition. Could the authors please comment on this issue?
In section "Methods", sub-section "Data", the authors explain that they "only found nine TFs that had comparable data in ChIP-seq and PWM". However, there are additional, similar data sets for at least Foxo1, Zfx, Tbx5, and Nr5a2 available (see Grau ., 2013 for details). Please clarify.

et al
In equation 1, parentheses might help to spot the argument of the product.
In equation 4, the parameter μ is not explained.
In Figure 3, the authors discuss the influence of sequence length on AUC. However, I would consider the influence of sequence length on the ranking of motifs more relevant for the topic of the manuscript.
In section "Results", sub-section "Effect of PBM data on motif assessment", 5th item of the list, the authors state that "PBM data are generated in 8-mers". As far as I know, PBMs are typically designed such that they cover all 10-mers (see, e.g., ). http://the_brain.bwh.harvard.edu/pbm.html In addition, the authors might clarify which data they used for the assessment (the probe sequences and intensities, I assume?). The manuscript "Transcription factor motif quality assessment requires systematic comparative analysis" by Kibet and Machanick addresses the assessment of transcription factor binding motifs. This question is especially important for selecting appropriate motifs for computational predictions given the large number of different motifs for the same transcription factor available from databases. Kibet and Machanick specifically consider different measures for motif scoring and assessment, and investigate different factors that might influence the assessment and, hence, the chosen motif.
F1000Research factors that might influence the assessment and, hence, the chosen motif.
The topic is of great relevance in any research dealing with sequence motifs and a systematic analysis of the factors influencing their assessment may help to develop a standardized framework for motif assessment. However, I have several reservations regarding the current version of the manuscript as outlined below.
As a general comment (that does not necessarily require a response by the authors), I found it slightly disappointing that the present manuscript does pose many important questions and potential obstacles in motif assessment, but does not provide a solution, be it guidelines for reasonable motif assessment or be it even a platform for performing such analyses. DATA: 1. I wonder why the authors decided to only consider ChIP-seq data but no data (PBMs or in-vitro SELEX). While binding may of greater relevance for many applications, problems like cell in-vivo type-specificity of motifs (also addressed by the authors) would have a minor influence. In addition, competitive or interaction effects with other transcription factors might be ruled out. Finally, some of the motif sources considered derive their PWMs from data. In summary, an analysis also using in-vitro in-vitro data might affect the conclusions of the paper.
2. For all ChIP-seq data sets under consideration, the authors extract (only) the top 500 highest scoring sequences for the assessment. This may have largely differing effects for different transcription factors, where, for instance, one transcription factor might have several hundred ChIP-seq positive regions, whereas another transcription factor might have tens of thousands of ChIP-seq positive regions. Hence, in one case also lowly occupied sequences are collected whereas in the other case, the positive data set may only contain the most stringent binding sequences. This may affect all downstream analyses and, for instance, could be one of the reasons why the authors observe transcription factor-specific effects for some factors. Hence, I would strongly suggest to conduct the analysis with a transcription factor-specific selection of sequences (where the simplest idea might be to use just a percentile).
3. a) The authors state that they construct a "similar" negative set. Here, the authors should clearly define what "similar" means, how sequences are selected, and how many negative sequences are in the set. b) In addition, the specific selection of negative sequences described by the authors (500bp "downstream", where "downstream" is also unclear as ChIP-seq regions lack an orientation), might introduce a specific bias, because under the assumption that transcription factor binding sites are often located close to the transcription start, which might mean that the negative sequences may already be coding and, hence, per se different from promoter sequences. c) Finally, from my experience, the choice of the negative data set strongly affects the performance assessment of motifs. Hence, the authors might consider to test an additional set of negative sequences (e.g., di-nucleotide shuffled positive sequences) in their analysis.

METHODS:
4. The "Methods" section, especially formulas needs substantial revision: a) In general, notation should be harmonized between the different formulas. For instance, the sequence S appears with different indexes with different meanings; the indicator function is denoted by S_i(b,m) in eqn. (5) and by I(S_{i,b}) in equation (6).
In addition:

F1000Research
In addition: b) In eqn. (1), parentheses are missing around (1 -P(...)). In addition, the notation S_{t+1:t+k}^{i} is not explained c) In eqn. (2), it is unclear if i and [S_{t=i}] are indexes or if this should denote a product of \theta, i and S_{t=i} (which I consider unlikely). In addition, the variable t (in the index) is neither bound nor explained. d) In eqn. (3), the upper limit of the sum is |s|, where it should be |S|, I assume. In addition, there seems to be something missing (a \theta?) in the product. e) Before eqn. (5), the authors refer to T as the length of the sequence. However, considering the formula, the length should be L, and the first sum from A to T refers to the alphabet. In addition, eqn. (6) again denotes the sum over the alphabet differently. f) In eqn. (5), the text refers to sequence S but the formula to sequence S_i g) In eq. (6), the variable P_b is not defined (the authors later only refer to p, which might have the same meaning). In addition, the authors to not explain, which background distribution they use in the assessment, which will be relevant, e.g., for the results presented in Fig. 6. 5. The energy scoring framework (eqn. 4 and 5) and the LogOdds scoring framework are formally defined only for sub-sequences and it remains unclear how these are applied to longer sequences from ChIP-seq. Are those subjected to the occupancy definitions (maximum and/or average) as well?
6. LogOdds scoring is referred to as "Log likelihood scoring" in the section's title (page 6, left column), which is not fully correct.
7. On page 6, right column, second paragraph, the authors state that they "wish to check the usefulness of correlation in motif assessment" (which I would find interesting), but I did not find any results regarding correlation as performance measure in the results. 9. On page 6, penultimate paragraph, the authors state that "the Foxa motif from the POUR data set is significantly differentially enriched only in the A549 cell line", which I could not read from Fig. 4. Please clarify.
10. On page 8, right column, the authors state that "MNCP prefers specific motifs, which will have more true positives". Could the authors elaborate on these findings and also possibly give an (mathematical) explanation?
11. In Fig. 6, the authors show AUC values for different motifs and scoring functions. a) First, it remains unclear which data sets have been used in this analysis for the different transcription factors. Is it just the average over all motifs and data sets for each factor? b) Second, I did, unfortunately, not get the general idea of this analysis. If I understood it correctly, the main question of this manuscript is to study the effects of different factors on motif assessment with the goal of selecting the most appropriate motif for a given transcription factor. However, here it seems to be that exactly this information is averaged out. Wouldn't be the more interesting question how the scoring functions affect the ranking (by AUC) of the different motifs for each transcription factor? utility of the motifs learned is predicting binding of the motifs. That said, we agree that data in vivo used in testing the motifs does have an effect on the ranking of the motifs. This is an observation we confirm with our re-analysis using PBM data. We have included a section on how assessment in PBM and ChIP-seq are influenced differently (Effect of PBM data on motif assessment).
2. Our choice for top 500 sequences was informed by our understanding of previous research. However we did not make it clear that such prior work supports this point, and we have now cited a reference. As advised, we decided to test if this would affect the results from this analysis. The ChIP-seq peaks we use have a median of 14000 peaks, the highest having 92,258 peaks and a minimum of 101 peaks. Where the number of the available peaks was less than 500, we used all the peaks. Given the median number of peaks, we found 5% of the peaks to be appropriate and we used this when 5% the of the total was more than 500, else we used top 500 peaks (or all of them, for data sets smaller than this). We also tested for 10% of the peaks. In all this, we found that the size of the peaks used had no significant effect on the results obtained. We, therefore, eliminate that as being one of the reasons for cell line specific binding. This may, however, have an effect on cell line specific ranking behaviour even if we did not observe that in our examples, given that the number of peaks differs for a given TF in different cell lines. We will definitely consider this suggestion when developing our tool to avoid any potential bias caused by this.
3. a) The manuscript has been updated. By similar we mean in size and sequence length. b) In our analysis, downstream is based on the coordinates of the peaks. We extracted sequences located 500bp from the highest coordinate (highest coordinate + 500). Our focus was to get negative sequences which are not expected to contain binding sites for the given TF but which maintains the nucleotide composition. That distance, whether it falls in a promoter region or not, should be appropriate in for our specific analysis. c) On other negative sequences that can be chosen, we agree that this can have some influence in the analysis. The scores obtained when a negative set generated using dinucleotide shuffled positive sequences were always lower than those from downstream sequences. However, the ranking of the motifs did not change in any significant manner. We expect random negative sets to have a significant influence on the ranks of the motifs and their probable difference in GC content from the binding region makes their appropriateness questionable.
4. The notations used in the formulas have been updated for uniformity. Thank you! 5. For energy scoring, the subsequence with the lowest energy is used to represent the sequence while for logOdds scoring, the score can either be obtained by getting the sum or the maximum score for all the sub-sequences. Clarified in the paper, thank you. 6. Corrections have been made in response to 6, 8 and 9. 7. We have updated the Figure 5 to include details on the usefulness of correlation statistics. We find them to produce significantly different ranks from MNCP or AUC or even between Pearson and Spearman correlations. 8. Done 9. Done F1000Research 9. Done 10. We have updated the paper in to add a line giving more information about MNCP. Simply put, the MNCP is a rank-based statistic that determines if the mean occurrence of a motif in test sequences is higher that the mean occurrence in a random set. Each set of sequences is ranked based on the mean occurrence, and the MNCP is calculated by finding the mean of the normalized ratio of the two ranks. 11. We have updated Figure 6 (now Figure 7) to address the comments. Our earlier figure actually averaged the information on the effect of scoring functions on the ranks of the motifs. We have updated by using the rank correlation of the motifs for various TFs to show how it affects ranking. 12. Table 3 (now Table 2) has been updated to include maximum occupancy.
13. On Egr1 motifs correlation of motif IC and scores, we have updated the statement to be in accordance with the data.
14. In Figure 8 (now Figure 9) we have updated the figure to only retain relevant columns. We have also corrected the error that led to identical entries' correlation being more than 1.
15. On why we chose to use the best motif's score to represent a database, we argue that since the focus of this analysis is to test our ability to choose for the best motif, irrespective of the database, we find using the best motif score to represent the DB to be sufficient. Besides, using median will still lead to biased results since DBs with many motifs of low quality and a few of high quality will be poorly ranked.
16. We only introduce CentriMo at a later stage of our analysis as an alternative method of scoring techniques to motif assessment. The focus of the paper was to systematically assess the factors that do influence motif assessment, so we wanted to maintain that focus. 17. We have taken your suggestion on Figure 9 (now 10) to normalize the scores. Thank you.
18. On the performance of TF2DNA, we agree that the low performance would be expected. We also believe that a different approach to motif assessment during motif discovery may have produced better motifs. In addition, testing using PBM data produced a much better performance. This may be a consequence of the motifs being short and only generated using methods. in vitro 19. The background section has been updated to include to making the observations balanced and including recent citations. 20. We accept that our statement on the lack of significant improvement of the motifs may have been misleading and unsupported. We have updated it to reflect current evidence. This type of meta-analysis will be of interest to a wide audience. However, the current manuscript needs considerable revision. In particular, the connections between the data presented and the conclusions reached need to be strengthened and clarified. Furthermore, a lot more clarification about what is being shown in the figures is needed to properly evaluate the conclusions. Below I have outlined specific examples through to Figure 8. The figure legends could definitely use more detail to help clarify what is being shown, and there needs to be more explicit and careful connection between the data and the conclusions (see examples below). I think that the type of analyses contained in this manuscript will be of interest to a wide audience; however, the manuscript needs to be substantially revised.
For each peak file, the 500 highest scored sequences were identified "after eliminating Methods/Data. repeat masked sites". It is a little unclear what this statement means. Does that mean that no peak was selected if there was any repeat masked seqeuence within the 50, 100 and 250bp windows? Or was the repeat masked sequence just masked and the genomic window extended to attain the 50, 100 and 250 bp cutoffs? Also, for the negative set, does 'similar' mean length-matched? It was exactly clear how this negative set was constructed.
It was not clear why only a subset of 15 of the Encode ChIP-seq datasets were used Figure 3 /results. and shown here, and how many datasets were used in the 'Average'? Also, the figure caption notes that 'all the motifs for the 15 TFs' were used, but it's not clear how many that was and whether the reported AUC values were averages over their individual AUC values? I little more clarification would be helpful.
The authors write, "Unless the interest is tissue-specific binding, if more than one set of data is Page 6. available, an average should be used". Used for what? For motif discovery? Why was 'energy scoring' used for this enrichment analysis, while GOMER scoring was used in   Figure 3? Are the results dependent on these scoring differences? If not, then for consistency sake, it would be helpful to limit the enrichment analyses to a single scoring scheme.
The authors conclude, "the Foxa motif from the POUR data set is significantly Page 6/ Figure 4. differentially enriched only in the A549 cell line and not so much in the other cell lines". I have no idea to what the authors are referring here, and this is the only conclusion from Figure 4. There are 5 different FOXA_discX.POUR motifs, all of which seem to score about the same on the different ChIP datasets. There is a FOXA1_2.GUERTIN that seems to be quite different, but this seems like an outlier within the dataset. I do not see how the data supports the contention that there are specific FOXA motifs that are better suited to particular ChIP datasets, it seems that for the most part they agree. Much more clarity is needed here.
Page 8 / Figure 5. "However, in some situations like Hnf4a and Ctcf, they are not ( Figure 5)". I only see Ctcf data represented in Figure 5, this should be clarified.
"The motifs ranked higher only by MNCP are generally long or with high IC (Table 2)". It would be much easier to see this if they were indicated somehow in Figure 5, perhaps with arrows are stars or something. Second, these conclusions don't seem to follow from the data at all. The CTCF_disc1.POUR seems also to score high with Energy_AUC, so it's not clear that the MNCP is the only factor of relevance here. The CTCF.1_5.ZLAB seems to be most affected by the Energy vs GOMER scoring, and not the MNCP approach. Even if these issues were resolved, it is impossible to know whether these motifs are 'generally long or with high IC' from Table 2, because the other motifs aren't shown. It would be much clearer if the mean and variance of the length & IC for all motifs were also provided for context, or even better correlate the relative score AUC to MNCP differences by length or IC, to truly see if a trend exists.

Figure 6.
Please clarify in the figure legend whether these values are for averages over multiple ChIP datasets (as was discussed above), and if so how these averages are determined.
"Maximum and sum log-odds scoring had low discriminative power for most of the motifs when all three statistical measures are used ( Figure 6)". What are the three statistical measures you're referring to, and where's the data? I only see data for AUC. Please clarify.
Please be explicit in the figure legend about what the 'Mean' and 'Median' refer to (i.e., mean Table 3. and median AUC values calculated over X single motif analyses described in Figure 6) "The variation in the scores is particularly reduced when MNCP statistic is used Figure 7/ page 10. (Figure 7)". How am I supposed to see this? What is a significant difference in MNCP and how does it compare to a difference in AUC. Based on the coloring scheme presented the results in Figure 6 and Figure 7 look very similar-it is not clear at all that there is any qualitative difference between these two figures except for the different measures used (i.e., an appropriate normalization might make them near equivalent).
It is not clear (nor mentioned) what is being shown in this figure. I assume -but I could be Figure 8. wrong -that we're looking at AUC values for each factor (i.e., Mef2a etc) averaged over some ChIP-seq datasets, but how are these being compared to each other? Further, how is Motif_IC which is a function just of the PWM being compared to a scoring function. I can't speak to the conclusions being reached as I don't currently know what data is being shown. Much more clarification is warranted in the text and figure caption.
I have read this submission. I believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard, however I have significant reservations, as outlined above.
No competing interests were disclosed.

Competing Interests:
Author Response 10 Feb 2016 , Rhodes University, South Africa Caleb Kipkurui Thank you very much for taking the time to review our paper and provide recommendations.Your comments have been very helpful in improving the paper.

Methods/Data
On repeat masked sequences, we have updated the paper to clarify the we did not include any sequences in test or negative sequences that contained masked positions. By a similar set of negative sequences, we mean matched in length and number of sequences.

Figure 3/Results
The figure caption is updated to clarify. The number of the TFs used was decided based on the availability of ChIP-seq data as well as having motifs in more than 10 of the databases used. A list of the motifs used is provided in the data repository as well as the specific ChIP-seq data used. See also Methods / Data paragraph 3.

Page 6:
On cell line specific binding, an average of the scores of all the available cell lines should be used in motif assessment. We have updated the statement for clarity.

Figure 4
We have changed to using results from GOMER scoring since they are similar; the effect described is only pronounced in Energy scoring.

Page 6/Figure 4
We had incorrectly mentioned the wrong motif to be significantly enriched. We have corrected this and also provided further evidence to the effect that the cell line used in the assessment does actually have an effect on the ranking of the motifs. The conclusions remain valid.
We acknowledge that the figure we had used did not present the intended information correctly. We changed the figure to present the general information on the effect of statistics on the ranking of the motifs. We observe that, when normalized, the MNCP and AUC scores do not differ, except for slight difference in some TFs like Hnf4a, Ctcf, Gata3. However, the Pearson and Spearman's correlation scores vary greatly. The plot of the standard deviation of scores as represented by error bars in Figure 6 demonstrates why we consider correlation scores to be F1000Research and Spearman's correlation scores vary greatly. The plot of the standard deviation of scores as represented by error bars in Figure 6 demonstrates why we consider correlation scores to be reliable than the other scores. We have added clarification of this point to the paper. Thank you again for pointing out the problem. The caption has been updated for clarity.  Table 2  Our apologies for the lack of detail. The figure caption has been updated for clarity. We correlate the scores for the various motifs (for each scoring function) to the length and information content of the motifs to determine whether the scores obtained are in any way influenced by the motif characteristics.
No competing interests were disclosed. Competing Interests: