BACKGROUND
the control of gene expression during development in drosophila and other metazoans is tightly directed by cis-acting regulatory sequences in the genome.
these dna sequences modulate expression of target genes by binding protein transcription factors   <cit> .
contact between a tf and dna sequence is mediated through the tf’s dna binding domain in a sequence dependent manner  <cit> .
each tf has one or more of a variety of different dna binding domains, including zinc fingers and homeoboxes  <cit> .
significant efforts have been undertaken to comprehend the organization of dna sequence at known binding regions and further understand how this influences the ability of a tf to bind.
our understanding of tf-dna interactions has been greatly aided by bioinformatic tools developed to analyze dna sequences obtained from experimental studies focused on identifying tf binding regions.
a key approach involves the construction of a position weight matrix   <cit> .
in pwm-based models, known binding regions for a given tf are first characterized by utilizing experimental data from dna footprinting assays, yeast one-hybrid assays, chromatin immunoprecipitation-sequencing  or protein binding microarrays   <cit> .
the binding regions are then aligned and trimmed to some minimal sequence length, l, and the frequency at which each nucleotide is observed at each position is recorded in a matrix of dimension  <dig> × l <cit>  .
once a pwm is constructed, these models aid in the discovery of de novo binding sites in silico, providing predictions for the location of additional binding regions in the genome, without the need for technically challenging in vitro binding assays  <cit> .
a host of computational tools have been developed that seek to streamline the discovery of de novo tf binding regions using pwms  <cit> .
however, a major limitation of pwms is their potential to lose information content during construction.
the lengths of pwms are often determined based on an optimal alignment between minimal sequences of varying length, potentially eliminating bordering regions crucial to determining a tf’s binding preference  <cit> .
extending pwms may therefore serve to increase their information content, and thus their predictive power  <cit> .
one limitation resulting from the experimental approaches to isolate tf-bound dna fragments  <cit> , is that there may be additional, but non-contiguous, bases that are fundamentally important to tf binding initiation  which are not represented in the experimental data and therefore not taken into account during traditional pwm construction.
a potential explanation for this lack of information content in canonical pwm construction is the omission of secondary binding by tfs with multiple dna binding domains  <cit> .
for example, in drosophila the hunchback tf has two distinct c2h2-type zinc-finger binding domains  <cit> .
if multiple dna binding domains contact sequences separately then each domain may contribute to the overall binding of the tf.
accordingly, in the case where there are two binding domains, one of the dna binding regions may be either: a) discarded because it fails to meet minimal fragment size requirements or b) incorporated into a combined alignment along with the sequences representing regions bound by the other binding domain.
either of these scenarios may lower the information content of the pwm.
the first scenario may result in a pwm that does not include all nucleotides necessary for in vivo binding .
the second, on the other hand, points to an even larger problem in pwm construction: the possibility that a tf may have two different modes of binding, and thus a single unique pwm is insufficient to predict all dna binding regions, which there has been strong evidence to support in the case of mammalian dna binding proteins  <cit> .
to address the limitations of pwms, we align and analyze predicted binding regions for four well-studied tfs in  <dig> cis-regulatory modules  that are essential to direct gene expression along the anterio-posterior axis in early drosophila development.
our analysis indicates that the current pwms for all four tfs examined exclude significant biases towards a given base, or bases, in specific positions in the neighboring sequences and that the information content of these pwms can be improved by including these additional sequences.
methods
cis-regulatory module and flanking sequences
we identify  <dig> genes of interest that display a pattern of differential expression along the anterio-posterior axis during drosophila development at or before stage  <dig>   <cit> .
in order to compile a database of crms, we utilize the redfly database  <cit>  to search each of the  <dig> genes and identify all those with in situ-verified crms .
since we are investigating the flanking region  of each pwm-predicted tf binding region, and tf binding sites may be predicted within the first or last 100 bps of a given crm, for each crm we obtained the entire sequence of the crm and 100 bp of flanking sequence in each direction.
the dna sequences of both the original crms and these extended crms are available in the additional file 1: dataset s <dig> and additional file 2: dataset s2).
to analyze sequence conservation in each crm we run evoprinterhd strict  <cit>  with default settings, specifying four drosophilid species for comparison: d. pseudoobscura, d. ananassae, d. erecta, and d. sechellia .
bioinformatic analysis
the pwms we use for caudal   <cit> , hunchback , knirps  and kruppel   <cit>  are as previously described.
for our analysis, we run patser  <cit>  with default settings .
to determine score ) cutoffs, we first observe the distribution of scores patser assigns, using each of the four individual pwms, to all the known binding regions used to construct each of the original respective pwms.
the cutoffs used are calculated by taking the 75th and 50th percentile cutoff of all these scores, and are referred to as 'strong’ and 'weak’ cutoffs respectively.
we then run patser on each of the original  <dig> crms  with each of the four pwms to predict binding regions.
only those regions scoring above the respective cutoff are used for further analysis and are referred to as core pwm predicted binding regions .
note that scores are all negative, so 'scoring above’ refers to patser outputting a ln less than or equal to the cutoff.
the strong cutoff is more stringent, only predicting binding regions that receive a score less than or equal to that obtained from the top 25% of known binding region scores , representing binding regions that are most similar to the consensus core binding region for the given tf.
the weak cutoff is less stringent, predicting binding regions that receive a score less than or equal to that obtained from the top 50% of known binding region scores , representing binding regions that are contained in a larger range of similarity to the consensus core binding region for the given tf.
one should note that these cutoff scores are tf-specific and are different for each of the four tfs analyzed.
in cases where overlapping binding regions were identified, both regions are included in all subsequent analyses.
lists of all the known binding regions for each tf, their corresponding ln obtained using patser, and whether they fall into the 75th percentile, 50th percentile, or neither are available in the additional file 3: table s <dig>
statistical analysis
for each tf we aggregate the core pwm-pbrs predicted by patser along with a fixed number, n, of bps of flanking sequence on each side, obtained from the  <dig> extended crms with the core pwm-pbr in the center and n bps on each side.
when binding sites are found on the reverse complement strand, we use the reverse complement sequence for analysis.
since each core pwm-pbr is exactly the same length  as the pwm associated with the given tf, each sequence of dna in the list is the same length.
this allows us to compute the frequency of each base  at each position .
these frequencies are shown for n =  <dig>  in the case of both strong and weak cutoff scores in figures  <dig> and  <dig>  and are listed in the additional file 4: table s <dig>
for each tf and cutoff score, given the list of pbrs , we run a chi-squared test on each position with the null hypothesis that at any given location, the distribution of a/c/g/t is exactly the same as the genome-wide distribution  <cit> , a:c:g:t. we note here that the overall nucleotide frequency in the  <dig> extended crms, a:c:g:t is not significantly different from the genome-wide distribution .
we analyze the results of these tests at three different confidence levels α =  <dig> , α =  <dig>  and α =  <dig> .
we choose more than one confidence level to control the familywise error rate for multiple comparisons.
a simple bonferroni correction leads to a corrected alpha value obtained by dividing alpha by the number of chi-squared tests .
thus, although the three alpha values stated can be interpreted without a bonferroni correction, α =  <dig>  can also be interpreted as a bonferonni corrected alpha corresponding to α =  <dig> .
the chi-squared values obtained are shown for n =  <dig> in figures  <dig> and  <dig>  and actual values for each nucleotide are in the additional file 4: table s <dig>  note that in figures  <dig> and  <dig>  the color-coding corresponds to the smallest alpha value of those analyzed in which the null hypothesis is rejected.
software availability
the web application that is used to run this analysis is freely available for non-commercial use at: http://drewell.sites.hmc.edu/projects/sequence_context_grapher.html.
jaspar database search
using the alignment produced by the bioinformatic analysis for the 'weak’  cutoff described above, a single pwm was constructed from a portion of the flanking sequence of hb showing statistically significant nucleotide bias  for use in a jaspar alignment search  <cit> .
the top  <dig> drosophila melanogaster tfs that give similarity scores to this pwm  are manually annotated for expression pattern in early  embryos  <cit> .
testing predictions for expanded hb pwms
when expanding the hb pwm, we use a 2-step process.
first, we choose an initial cutoff score and predict and align all binding sites using the original core pwm as described in 'bioinformatic analysis’.
the pwm-pbrs identified with the original core pwm using this initial cutoff score are then extended by including flanking regions of interest .
an extended pwm is constructed from the base frequencies of these extended pwm-pbrs.
next, patser is used with the extended pwm to determine the score for each of the extended pwm-pbrs constructed from the set of predicted hb binding sites with the initial cutoff score.
computing the percentiles of these scores, in the same way as is described for the original pwm analysis in 'bioinformatic analysis’, we obtain a secondary cutoff score.
lists of the extended pwm-pbrs, their corresponding ln obtained using patser with the corresponding extended hb pwm, and whether they fall into the 0th, 25th, 50th or 75th percentile are available in the additional file 5: table s <dig>  we run patser again with the extended pwm and secondary cutoff score on the iab7b crm, which contains one known functional hb binding site  <cit> .
this allows us to then compare the location of the predicted binding sites obtained using extended hb pwms of varying lengths to the known hb binding site to determine which pwms result in the lowest number of false positive and false negative predictions.
analyzing expanded pwms using chip-seq datasets
to analyze the predictive power of the extended pwms, we first choose initial cutoff scores representing the 25th and 50th percentile scores and generate pwm-pbrs for each of the four tfs as described in 'bioinformatic analysis’.
we then extend those core pwm-pbrs  to include the core and all highly significant  flanking sequence context-dependent biases , as well as the core with - <dig> to + <dig> flanking regions used to generate the pwm-pbr .
we use a secondary cutoff score corresponding to the 0th percentile as described in 'testing predictions for expanded hb pwms’.
chip-seq datasets for each tf from the bdgp  <cit>  are filtered to include only those peaks with more than 100 bp of sequence.
we run patser with each of the three different pwms for each tf on their respective tf peaks and score a true positive prediction when the pwm predicts at least one tf binding site.
for each chip-seq peak, we calculate the nucleotide distribution within the peak and create  <dig> 'scrambled peaks’, random dna sequences of the same length and nucleotide distribution.
we then run patser with each of the three different pwms for each tf on these scrambled peaks and score a false positive prediction when the pwm predicts at least one tf binding site.
both the true positive and false positive results for the 25th and 50th percentile initial cutoff scores at the 0th, 25th, 50th, and 75th percentile secondary cutoff scores are available in the additional file 6: table s <dig>
RESULTS
when considering the possibility of sequence context-dependence for tf binding, evidence has pointed toward the existence of nucleotide biases at positions in close proximity to a region experimentally verified or computationally predicted to bind a tf  <cit> .
to test this idea, we analyze  <dig> crms that are active during early drosophila development for predicted binding regions for four tfs using patser .
these four tfs: caudal , hunchback , knirps  and kruppel  are all critical for normal development and are present in spatially restricted patterns along the anterio-posterior axis in early embryogenesis  <cit> .
a number of in vivo confirmed minimal binding regions have been characterized for these tfs and the existing pwms for each of these factors range in size from  <dig> to 9 bp   <cit> .
of greatest importance for this study, their current canonical pwms have been proven to have greater predictive power for experimentally validated tf binding regions, when compared to other published pwms  <cit> .
if context-dependent biases are present in sequences near these characterized binding regions, we predict that these bases would be evolutionarily conserved.
by combining patser  <cit>  analysis with evoprinterhd  <cit>  analysis, we are able to identify several examples of extended regions of sequence conservation surrounding evolutionarily conserved tf binding regions, including portions of the even-skipped stripe  <dig> crm and the paired stripe  <dig> crm .
in all cases within the depicted portions of the even-skipped stripe  <dig> crm and paired stripe  <dig> crm, predicted tf binding regions are flanked by substantial extended sequence conservation on one or both sides.
this presents a testable hypothesis: that these regions of extended conservation contain functionally important flanking bases that are important for robust tf binding.
to address the hypothesis that there may be sequence context-dependent binding for the four tfs, we investigate the sequences 25 bp up- and downstream of defined core pwm-predicted binding regions  .
alignment of the core pwm-pbrs and their flanking regions for each individual tf does indeed reveal a statistically significant enrichment of certain bases outside of the core pwm-pbrs .
a very clear example of this enrichment with high statistical significance  is found at binding regions for hb.
using the weak cutoff value , beyond the hb core pwm-pbr there is context-dependent bias at the first two and the 7th nucleotide downstream  of the core pwm-pbr .
in addition, there are four clusters of context-dependent bias upstream of the hb core pwm-pbr at positions - <dig>  - <dig> to - <dig>  - <dig> to - <dig> and - <dig> .
kni seems to follow a similar pattern to hb, with nucleotide enrichment bias at  <dig> positions downstream and  <dig> upstream of the core pwm-pbr .
this enrichment bias is also seen for cad and kr , but is not as prevalent.
cad and kr display only short stretches of sequence with robust context-dependent bias, and in both cases these sequences are largely contiguous to the core pwm-pbrs .
for all four tfs the enrichment biases at positions neighboring the defined core pwm-pbrs could be incorporated in to expanded pwms.
to directly test the ability of expanded pwms to accurately predict tf binding sites, we examine the performance of multiple hb pwms of varying length on the iab7b crm from the bithorax complex.
this 828 bp crm has been extensively characterized  <cit>  and contains a single hb binding site in the highly-conserved 154 bp signature motif  <cit> , which is sufficient to account for the functional repression of the crm mediated by hb in the anterior half of the embryo.
using a relatively stringent initial cutoff , the 7 bp hb pwm corresponding to the core pwm-pbr fails to predict any binding sites in the enhancer at all secondary cutoff values .
extending the pwm by only 3 bp  allows the prediction of the true positive site .
as the pwm is expanded to include more flanking sequence it continues to perform better.
the most robust prediction of the true positive hb binding site is with the - <dig> to + <dig> pwm.
at this stringent initial cutoff value there are no false positive predictions with any of the hb pwms.
we also investigate the predictive power of the pwms under less stringent conditions, utilizing an initial cutoff equivalent to the 25th percentile .
using this less stringent initial cutoff there are many more predicted false positive hb binding sites .
the hb pwm corresponding to the core pwm-pbr gives the highest false negative rate; it fails to predict the single true positive hb site in the enhancer at any of the secondary cutoff values .
once more, as the pwm is extended it generally performs better.
the most robust true positive predictions are with the - <dig> to + <dig> hb pwm , suggesting that expanding the pwm to include all flanking sequences that show a statistically significant nucleotide bias improves the predictive power of the pwm.
with the less stringent initial cutoff value, in all cases, shorter hb pwms reduce the false positive rate but at the expense of increasing the false negative rate .
to further analyze the performance of the extended pwms, we generate three different pwms for each of the four tfs and test the predictive power on chip-seq peak sequences for each tf.
the three different pwms correspond to: the core pwm-pbr , a region encompassing the core pbr and all statistically significant flanking sequence context-dependent biases , and the core with 25 bp of flanking sequence on each side .
comparison of the ability of each of the three pwms to correctly identify at least one binding site in individual chip-seq peaks for each tf reveals that in each case, with the exception of the cad pwm at low stringency , the extended pwms significantly outperform the core pwms .
there are many potential reasons for the fact that all four tfs exhibit context-dependent biases within their extended binding regions, including specific interactions dependent on the physical constraints of the tf contacting dna, nucleosomes or co-factors.
to address the possibility that co-factors may in fact be recruited to these neighboring sequences we investigate a cluster of context-dependent sequence exhibiting strong statistical significance  for additional tf binding sites.
alignment of the pwm corresponding to the - <dig> to - <dig> dna region to all pwms available in the jaspar database identifies putative binding sites for other drosophila tfs .
of the top five binding sites, ranked according to their binding motif’s similarity to the pwm constructed from the - <dig> to - <dig> hb region, only two  and sd) correspond to tfs expressed in the early embryo .
none of the top five ranked tfs are known to have any functional interaction with hb, suggesting that secondary sequence context-dependence is not sufficiently explained by co-factor binding alone.
a further testable explanation for the fact that all four tfs exhibit context-dependent biases at their extended binding regions may be that these tfs have multiple dna binding domains, each of which contacts different nucleotides independently  <cit> , but together act to increase the tf’s binding affinity for target sequences.
if this is the case, these secondary dna binding domains may enhance, but not replace the function of the primary, canonical binding domain in the tf protein.
to assess this hypothesis, we compare the binding regions predicted by patser using both strong and weak cutoff scores .
while many of the nucleotide position biases persist at either cutoff score, we find a general decrease in the number and significance of context-dependent biases for all four tfs when we only consider strong sites .
for example, context-dependent biases within 15 bp of the kni pwm-pbr are greatly reduced when we only consider strong sites - of the  <dig> enriched biases found with the weak cutoff, only one remains .
the range of decrease across the four tfs is variable .
overall, these results suggest two separate binding regimes:  if a tf has secondary, non-contiguous binding, stronger core binding regions may overcome a paucity of context-dependent biases in the flanking sequences, whereas weaker core binding regions may depend more heavily on these biases;  if a tf exhibits contiguous biases, these biases may simply suggest a larger canonical core binding region.
however, in some cases, there are nucleotides that are found to be significant only when the strong cutoff score is used.
for example, two new biases are detected at + <dig> and - <dig> relative to the kni pwm-pbr .
this finding supports a hypothesis that strong and weak core binding regions may in fact have two different functional roles that allow the tf to bind in a different mode, again posing the question of whether it is valid or not to represent a tf binding with one unique pwm  <cit> .
to address the possibility that the presence of secondary sequence context-dependence in the case of a particular tf is due to the tf having multiple binding domains, as previously shown for a large number of mammalian tfs  <cit> , we investigate the binding domains of each of the four drosophila tfs.
hb has two groups of c2h2-type zinc fingers, separated by over  <dig> amino acids, while kr has only one group of fingers .
this may explain the stark difference between the significance profiles at flanking sequences neighboring defined pwm-pbrs for hb and kr .
in addition, it may reveal a relationship between the predicted binding strength and number of context-dependent biases for each tf.
one model is that the group of four zinc-fingers in the center of hb form the core dna-binding domain , while the other group of two auxiliary zinc-fingers form a secondary binding domain  with less contribution to the overall binding stability .
there is some experimental evidence to support this model, as the highly conserved multi-zinc finger tf ctcf has been shown to contact dna nucleotides at least 12 bp  apart at the human c-myc promoter  <cit> .
in comparison, the single group of zinc-fingers in kr form a single dna binding domain which can only bind the core pwm-pbr and directly adjacent sequences .
cad also has only one binding domain , potentially explaining cad's lower level of secondary sequence context-dependence, when compared to hb or kni.
of the four tfs, kni is the most puzzling, exhibiting secondary sequence context-dependency while having only a single annotated binding domain .
one potential explanation may be that, unlike in the case of kr where the two zinc fingers are only six amino acids apart, the two c4-type zinc fingers in kni are separated by  <dig> amino acids.
this increased separation may in fact allow the two zinc fingers to act as functionally distinct dna binding domains, as is the case for hb.
CONCLUSIONS
taken together these data suggest that current pwms may not be optimal to explain the complexity of tf binding.
although we only test four tfs in this study, we demonstrate that all four tfs exhibit context-dependent biases towards given nucleotides both contiguously with the defined minimal binding region and non-contiguously in flanking dna sequences, thus providing a foundation for this to be explored more broadly.
an additional intriguing question for future study will be to investigate if the context-dependent bias persists at predicted tf binding regions in other genomic regions that are not characterized as crms.
by taking these secondary context-dependencies into account, we propose that the information content of pwms can be expanded in many cases.
this expansion would not only provide better predictions of true tf binding regions in the genome, but may also help improve estimates of relative binding affinities at specific sites, allowing one to understand the molecular basis for the difference between weak and strong binding sites.
the ability to identify novel crms and decipher the sequence organization at crms relies heavily on a concrete understanding of tf binding preferences.
improving the information content of pwms and our comprehension of tf binding events will contribute to these continued efforts.
availability of supporting data
the data sets supporting the results of this article are included within the article
competing interests
the authors declare that they have no competing interests.
authors’ contributions
jls created the web application, ran all bioinformatic and statistical analyses, and tested expanded hb pwms.
asb compiled the  <dig> crms, ran the evoprinterhd analysis, performed the jaspar database search, and analyzed each tf’s binding domain.
the manuscript was written by jls, asb, rad and jmd, and the overall project was conceived and guided by rad and jmd.
all authors read and approved the final manuscript.
supplementary material
additional file 1
dataset s <dig>  dna sequences of the original crm .
the original crm folder  contains the  <dig> fasta files for the crms we used.
click here for file
additional file 2
dataset s <dig>  dna sequences of the extended crms .
the extended crm folder  contains those same  <dig> crms, extended by including the neighboring  <dig> bp both up- and downstream of the crm.
the extended crm fasta filenames are identical to the original, but with a ' + ’ appended to the filename .
click here for file
additional file 3: table s1
patser scores for predicted binding regions.
each sheet contains a list of the sequences in predicted binding regions for a particular tf.
sequences are ordered by their ln, as computed by patser .
the scores highlighted in yellow indicate those scoring at or above  less than or equal to) those included in the 75th percentile .
the scores highlighted in green indicate those scoring at or above  less than or equal to) those included in the 50th percentile .
recall that the 75th and 50th percentile cutoff scores are specific for each tf and thus vary from sheet to sheet.
click here for file
additional file 4: table s2
nucleotide frequencies in flanking sequences neighboring predicted binding regions.
each sheet corresponds to a particular tf and cutoff percentile .
the column labeled 'position’ assigns alphabetic values to the consensus core pwm-predicted binding regions  and numeric values to the neighboring flanking sequences, ranging from ‒ <dig> to ‒ <dig> upstream and from  <dig> to  <dig> downstream.
the next four columns, labeled 'a’, 'c’, 'g’, and 't’, contain the frequency observed for each nucleotide at the specified position when all pwm-pbrs are considered for that particular tf and cutoff percentile.
the column labeled 'chi‒squared’ contains the corresponding chi-squared value for each position.
click here for file
additional file 5: table s3
patser scores for extended pwm-pbrs.
each sheet contains a list of the sequences in the extended pwm-pbrs for a particular initial cutoff score and length of flanking region.
sequences are ordered by their ln, as computed by patser using the corresponding extended hb pwm .
the scores highlighted in yellow indicate those scoring at or above  less than or equal to) those included in the 75th percentile .
the scores highlighted in green indicate those scoring at or above  less than or equal to) those included in the 50th percentile .
the scores highlighted in orange indicate those scoring at or above  less than or equal to) those included in the 25th percentile .
the scores highlighted in blue indicate those scoring at or above  less than or equal to) those included in the 0th percentile .
this includes all sites detected with the initial cutoff score.
recall that the 0th, 25th, 50th and 75th percentile cutoff scores are specific for each initial cutoff score and length of flanking region and thus vary from sheet to sheet.
click here for file
additional file 6: table s4
true and false positive results for extended pwm-pbrs on chip-seq peaks.
each sheet corresponds to a particular tf.
the name of the tf, the total number of chip-seq peaks, and the number of scrambled peaks used for the analysis are listed in the first column.
in the other columns, each row corresponds to a patser run for the given tf, chip-seq peaks and scrambled peaks.
the primary and secondary cutoff percentiles used, as well as the number of nucleotides extended to the left  and right  of the core pwm-pbr are listed.
the results of each run are shown as the number of chip peaks containing at least one predicted binding site, which is used to calculate the true positive rate, and the number of scrambled chip peaks containing at least one predicted binding site, which is used to calculate the false positive rate.
the last column contains the ratio of true positive to false positive.
to highlight settings in which the sensitivity increased in at least one of the extended pwm-pbrs compared to the core, we have highlighted those true positive rates in green, to highlight settings in which the specificity increased in at least one of the extended pwm-pbrs compared to the core, we have highlighted those false positive rates in pink, and to highlight settings in which the ratio of true to false positives increased in at least one of the extended pwm-pbrs compared to the core, we have highlighted those ratios in yellow.
click here for file
acknowledgments
the research in this paper was supported by funding to r.a.d.
from the national institutes of health , the national science foundation  and howard hughes medical institute undergraduate science education program grants  to the biology department at harvey mudd college.
j.m.d.
was funded as a teaching and research postdoctoral fellow, supported in part by nsf grant dms- <dig>  and through research funds provided by amherst college.

 the control of gene expression during development in drosophila and other metazoans is tightly directed by cis-acting regulatory sequences in the genome. these dna sequences modulate expression of target genes by binding protein transcription factors. contact between a tf and dna sequence is mediated through the tf’s dna binding domain in a sequence dependent manner. each tf has one or more of a variety of different dna binding domains, including zinc fingers and homeoboxes. significant efforts have been undertaken to comprehend the organization of dna sequence at known binding regions and further understand how this influences the ability of a tf to bind. our understanding of tf-dna interactions has been greatly aided by bioinformatic tools developed to analyze dna sequences obtained from experimental studies focused on identifying tf binding regions. a key approach involves the construction of a position weight matrix. in pwm-based models, known binding regions for a given tf are first characterized by utilizing experimental data from dna footprinting assays, yeast one-hybrid assays, chromatin immunoprecipitation-sequencing or protein binding microarrays. the binding regions are then aligned and trimmed to some minimal sequence length, l, and the frequency at which each nucleotide is observed at each position is recorded in a matrix of dimension × l. once a pwm is constructed, these models aid in the discovery of de novo binding sites in silico, providing predictions for the location of additional binding regions in the genome, without the need for technically challenging in vitro binding assays. a host of computational tools have been developed that seek to streamline the discovery of de novo tf binding regions using pwms. however, a major limitation of pwms is their potential to lose information content during construction. the lengths of pwms are often determined based on an optimal alignment between minimal sequences of varying length, potentially eliminating bordering regions crucial to determining a tf’s binding preference. extending pwms may therefore serve to increase their information content, and thus their predictive power. one limitation resulting from the experimental approaches to isolate tf-bound dna fragments, is that there may be additional, but non-contiguous, bases that are fundamentally important to tf binding initiation which are not represented in the experimental data and therefore not taken into account during traditional pwm construction. a potential explanation for this lack of information content in canonical pwm construction is the omission of secondary binding by tfs with multiple dna binding domains. for example, in drosophila the hunchback tf has two distinct c2h2-type zinc-finger binding domains. if multiple dna binding domains contact sequences separately then each domain may contribute to the overall binding of the tf. accordingly, in the case where there are two binding domains, one of the dna binding regions may be either: a) discarded because it fails to meet minimal fragment size requirements or b) incorporated into a combined alignment along with the sequences representing regions bound by the other binding domain. either of these scenarios may lower the information content of the pwm. the first scenario may result in a pwm that does not include all nucleotides necessary for in vivo binding. the second, on the other hand, points to an even larger problem in pwm construction: the possibility that a tf may have two different modes of binding, and thus a single unique pwm is insufficient to predict all dna binding regions, which there has been strong evidence to support in the case of mammalian dna binding proteins. to address the limitations of pwms, we align and analyze predicted binding regions for four well-studied tfs in cis-regulatory modules that are essential to direct gene expression along the anterio-posterior axis in early drosophila development. our analysis indicates that the current pwms for all four tfs examined exclude significant biases towards a given base, or bases, in specific positions in the neighboring sequences and that the information content of these pwms can be improved by including these additional sequences. methods cis-regulatory module and flanking sequences we identify genes of interest that display a pattern of differential expression along the anterio-posterior axis during drosophila development at or before stage. in order to compile a database of crms, we utilize the redfly database to search each of the genes and identify all those with in situ-verified crms. since we are investigating the flanking region of each pwm-predicted tf binding region, and tf binding sites may be predicted within the first or last 100 bps of a given crm, for each crm we obtained the entire sequence of the crm and 100 bp of flanking sequence in each direction. the dna sequences of both the original crms and these extended crms are available in the additional file 1: dataset s and additional file 2: dataset s2). to analyze sequence conservation in each crm we run evoprinterhd strict with default settings, specifying four drosophilid species for comparison: d. pseudoobscura, d. ananassae, d. erecta, and d. sechellia. bioinformatic analysis the pwms we use for caudal, hunchback, knirps and kruppel are as previously described. for our analysis, we run patser with default settings. to determine score ) cutoffs, we first observe the distribution of scores patser assigns, using each of the four individual pwms, to all the known binding regions used to construct each of the original respective pwms. the cutoffs used are calculated by taking the 75th and 50th percentile cutoff of all these scores, and are referred to as'strong’ and'weak’ cutoffs respectively. we then run patser on each of the original crms with each of the four pwms to predict binding regions. only those regions scoring above the respective cutoff are used for further analysis and are referred to as core pwm predicted binding regions. note that scores are all negative, so'scoring above’ refers to patser outputting a ln less than or equal to the cutoff. the strong cutoff is more stringent, only predicting binding regions that receive a score less than or equal to that obtained from the top 25% of known binding region scores, representing binding regions that are most similar to the consensus core binding region for the given tf. the weak cutoff is less stringent, predicting binding regions that receive a score less than or equal to that obtained from the top 50% of known binding region scores, representing binding regions that are contained in a larger range of similarity to the consensus core binding region for the given tf. one should note that these cutoff scores are tf-specific and are different for each of the four tfs analyzed. in cases where overlapping binding regions were identified, both regions are included in all subsequent analyses. lists of all the known binding regions for each tf, their corresponding ln obtained using patser, and whether they fall into the 75th percentile, 50th percentile, or neither are available in the additional file 3: table s statistical analysis for each tf we aggregate the core pwm-pbrs predicted by patser along with a fixed number, n, of bps of flanking sequence on each side, obtained from the extended crms with the core pwm-pbr in the center and n bps on each side. when binding sites are found on the reverse complement strand, we use the reverse complement sequence for analysis. since each core pwm-pbr is exactly the same length as the pwm associated with the given tf, each sequence of dna in the list is the same length. this allows us to compute the frequency of each base at each position. these frequencies are shown for n =  in the case of both strong and weak cutoff scores in figures  and and are listed in the additional file 4: table s for each tf and cutoff score, given the list of pbrs, we run a chi-squared test on each position with the null hypothesis that at any given location, the distribution of a/c/g/t is exactly the same as the genome-wide distribution, a:c:g:t. we note here that the overall nucleotide frequency in the extended crms, a:c:g:t is not significantly different from the genome-wide distribution. we analyze the results of these tests at three different confidence levels α = , α =  and α = . we choose more than one confidence level to control the familywise error rate for multiple comparisons. a simple bonferroni correction leads to a corrected alpha value obtained by dividing alpha by the number of chi-squared tests. thus, although the three alpha values stated can be interpreted without a bonferroni correction, α =  can also be interpreted as a bonferonni corrected alpha corresponding to α = . the chi-squared values obtained are shown for n =  in figures  and and actual values for each nucleotide are in the additional file 4: table s note that in figures  and the color-coding corresponds to the smallest alpha value of those analyzed in which the null hypothesis is rejected. software availability the web application that is used to run this analysis is freely available for non-commercial use at: http://drewell.sites.hmc.edu/projects/sequence_context_grapher.html. jaspar database search using the alignment produced by the bioinformatic analysis for the'weak’ cutoff described above, a single pwm was constructed from a portion of the flanking sequence of hb showing statistically significant nucleotide bias for use in a jaspar alignment search. the top drosophila melanogaster tfs that give similarity scores to this pwm are manually annotated for expression pattern in early embryos. testing predictions for expanded hb pwms when expanding the hb pwm, we use a 2-step process. first, we choose an initial cutoff score and predict and align all binding sites using the original core pwm as described in'bioinformatic analysis’. the pwm-pbrs identified with the original core pwm using this initial cutoff score are then extended by including flanking regions of interest. an extended pwm is constructed from the base frequencies of these extended pwm-pbrs. next, patser is used with the extended pwm to determine the score for each of the extended pwm-pbrs constructed from the set of predicted hb binding sites with the initial cutoff score. computing the percentiles of these scores, in the same way as is described for the original pwm analysis in'bioinformatic analysis’, we obtain a secondary cutoff score. lists of the extended pwm-pbrs, their corresponding ln obtained using patser with the corresponding extended hb pwm, and whether they fall into the 0th, 25th, 50th or 75th percentile are available in the additional file 5: table s we run patser again with the extended pwm and secondary cutoff score on the iab7b crm, which contains one known functional hb binding site. this allows us to then compare the location of the predicted binding sites obtained using extended hb pwms of varying lengths to the known hb binding site to determine which pwms result in the lowest number of false positive and false negative predictions. analyzing expanded pwms using chip-seq datasets to analyze the predictive power of the extended pwms, we first choose initial cutoff scores representing the 25th and 50th percentile scores and generate pwm-pbrs for each of the four tfs as described in'bioinformatic analysis’. we then extend those core pwm-pbrs to include the core and all highly significant flanking sequence context-dependent biases, as well as the core with - to + flanking regions used to generate the pwm-pbr. we use a secondary cutoff score corresponding to the 0th percentile as described in'testing predictions for expanded hb pwms’. chip-seq datasets for each tf from the bdgp are filtered to include only those peaks with more than 100 bp of sequence. we run patser with each of the three different pwms for each tf on their respective tf peaks and score a true positive prediction when the pwm predicts at least one tf binding site. for each chip-seq peak, we calculate the nucleotide distribution within the peak and create'scrambled peaks’, random dna sequences of the same length and nucleotide distribution. we then run patser with each of the three different pwms for each tf on these scrambled peaks and score a false positive prediction when the pwm predicts at least one tf binding site. both the true positive and false positive results for the 25th and 50th percentile initial cutoff scores at the 0th, 25th, 50th, and 75th percentile secondary cutoff scores are available in the additional file 6: table s  when considering the possibility of sequence context-dependence for tf binding, evidence has pointed toward the existence of nucleotide biases at positions in close proximity to a region experimentally verified or computationally predicted to bind a tf. to test this idea, we analyze crms that are active during early drosophila development for predicted binding regions for four tfs using patser. these four tfs: caudal, hunchback, knirps and kruppel are all critical for normal development and are present in spatially restricted patterns along the anterio-posterior axis in early embryogenesis. a number of in vivo confirmed minimal binding regions have been characterized for these tfs and the existing pwms for each of these factors range in size from to 9 bp. of greatest importance for this study, their current canonical pwms have been proven to have greater predictive power for experimentally validated tf binding regions, when compared to other published pwms. if context-dependent biases are present in sequences near these characterized binding regions, we predict that these bases would be evolutionarily conserved. by combining patser analysis with evoprinterhd analysis, we are able to identify several examples of extended regions of sequence conservation surrounding evolutionarily conserved tf binding regions, including portions of the even-skipped stripe crm and the paired stripe crm. in all cases within the depicted portions of the even-skipped stripe crm and paired stripe crm, predicted tf binding regions are flanked by substantial extended sequence conservation on one or both sides. this presents a testable hypothesis: that these regions of extended conservation contain functionally important flanking bases that are important for robust tf binding. to address the hypothesis that there may be sequence context-dependent binding for the four tfs, we investigate the sequences 25 bp up- and downstream of defined core pwm-predicted binding regions. alignment of the core pwm-pbrs and their flanking regions for each individual tf does indeed reveal a statistically significant enrichment of certain bases outside of the core pwm-pbrs. a very clear example of this enrichment with high statistical significance is found at binding regions for hb. using the weak cutoff value, beyond the hb core pwm-pbr there is context-dependent bias at the first two and the 7th nucleotide downstream of the core pwm-pbr. in addition, there are four clusters of context-dependent bias upstream of the hb core pwm-pbr at positions - - to - - to - and -. kni seems to follow a similar pattern to hb, with nucleotide enrichment bias at positions downstream and upstream of the core pwm-pbr. this enrichment bias is also seen for cad and kr, but is not as prevalent. cad and kr display only short stretches of sequence with robust context-dependent bias, and in both cases these sequences are largely contiguous to the core pwm-pbrs. for all four tfs the enrichment biases at positions neighboring the defined core pwm-pbrs could be incorporated in to expanded pwms. to directly test the ability of expanded pwms to accurately predict tf binding sites, we examine the performance of multiple hb pwms of varying length on the iab7b crm from the bithorax complex. this 828 bp crm has been extensively characterized and contains a single hb binding site in the highly-conserved 154 bp signature motif, which is sufficient to account for the functional repression of the crm mediated by hb in the anterior half of the embryo. using a relatively stringent initial cutoff, the 7 bp hb pwm corresponding to the core pwm-pbr fails to predict any binding sites in the enhancer at all secondary cutoff values. extending the pwm by only 3 bp allows the prediction of the true positive site. as the pwm is expanded to include more flanking sequence it continues to perform better. the most robust prediction of the true positive hb binding site is with the - to + pwm. at this stringent initial cutoff value there are no false positive predictions with any of the hb pwms. we also investigate the predictive power of the pwms under less stringent conditions, utilizing an initial cutoff equivalent to the 25th percentile. using this less stringent initial cutoff there are many more predicted false positive hb binding sites. the hb pwm corresponding to the core pwm-pbr gives the highest false negative rate; it fails to predict the single true positive hb site in the enhancer at any of the secondary cutoff values. once more, as the pwm is extended it generally performs better. the most robust true positive predictions are with the - to + hb pwm, suggesting that expanding the pwm to include all flanking sequences that show a statistically significant nucleotide bias improves the predictive power of the pwm. with the less stringent initial cutoff value, in all cases, shorter hb pwms reduce the false positive rate but at the expense of increasing the false negative rate. to further analyze the performance of the extended pwms, we generate three different pwms for each of the four tfs and test the predictive power on chip-seq peak sequences for each tf. the three different pwms correspond to: the core pwm-pbr, a region encompassing the core pbr and all statistically significant flanking sequence context-dependent biases, and the core with 25 bp of flanking sequence on each side. comparison of the ability of each of the three pwms to correctly identify at least one binding site in individual chip-seq peaks for each tf reveals that in each case, with the exception of the cad pwm at low stringency, the extended pwms significantly outperform the core pwms. there are many potential reasons for the fact that all four tfs exhibit context-dependent biases within their extended binding regions, including specific interactions dependent on the physical constraints of the tf contacting dna, nucleosomes or co-factors. to address the possibility that co-factors may in fact be recruited to these neighboring sequences we investigate a cluster of context-dependent sequence exhibiting strong statistical significance for additional tf binding sites. alignment of the pwm corresponding to the - to - dna region to all pwms available in the jaspar database identifies putative binding sites for other drosophila tfs. of the top five binding sites, ranked according to their binding motif’s similarity to the pwm constructed from the - to - hb region, only two and sd) correspond to tfs expressed in the early embryo. none of the top five ranked tfs are known to have any functional interaction with hb, suggesting that secondary sequence context-dependence is not sufficiently explained by co-factor binding alone. a further testable explanation for the fact that all four tfs exhibit context-dependent biases at their extended binding regions may be that these tfs have multiple dna binding domains, each of which contacts different nucleotides independently, but together act to increase the tf’s binding affinity for target sequences. if this is the case, these secondary dna binding domains may enhance, but not replace the function of the primary, canonical binding domain in the tf protein. to assess this hypothesis, we compare the binding regions predicted by patser using both strong and weak cutoff scores. while many of the nucleotide position biases persist at either cutoff score, we find a general decrease in the number and significance of context-dependent biases for all four tfs when we only consider strong sites. for example, context-dependent biases within 15 bp of the kni pwm-pbr are greatly reduced when we only consider strong sites - of the enriched biases found with the weak cutoff, only one remains. the range of decrease across the four tfs is variable. overall, these results suggest two separate binding regimes: if a tf has secondary, non-contiguous binding, stronger core binding regions may overcome a paucity of context-dependent biases in the flanking sequences, whereas weaker core binding regions may depend more heavily on these biases; if a tf exhibits contiguous biases, these biases may simply suggest a larger canonical core binding region. however, in some cases, there are nucleotides that are found to be significant only when the strong cutoff score is used. for example, two new biases are detected at + and - relative to the kni pwm-pbr. this finding supports a hypothesis that strong and weak core binding regions may in fact have two different functional roles that allow the tf to bind in a different mode, again posing the question of whether it is valid or not to represent a tf binding with one unique pwm. to address the possibility that the presence of secondary sequence context-dependence in the case of a particular tf is due to the tf having multiple binding domains, as previously shown for a large number of mammalian tfs, we investigate the binding domains of each of the four drosophila tfs. hb has two groups of c2h2-type zinc fingers, separated by over amino acids, while kr has only one group of fingers. this may explain the stark difference between the significance profiles at flanking sequences neighboring defined pwm-pbrs for hb and kr. in addition, it may reveal a relationship between the predicted binding strength and number of context-dependent biases for each tf. one model is that the group of four zinc-fingers in the center of hb form the core dna-binding domain, while the other group of two auxiliary zinc-fingers form a secondary binding domain with less contribution to the overall binding stability. there is some experimental evidence to support this model, as the highly conserved multi-zinc finger tf ctcf has been shown to contact dna nucleotides at least 12 bp apart at the human c-myc promoter. in comparison, the single group of zinc-fingers in kr form a single dna binding domain which can only bind the core pwm-pbr and directly adjacent sequences. cad also has only one binding domain, potentially explaining cad's lower level of secondary sequence context-dependence, when compared to hb or kni. of the four tfs, kni is the most puzzling, exhibiting secondary sequence context-dependency while having only a single annotated binding domain. one potential explanation may be that, unlike in the case of kr where the two zinc fingers are only six amino acids apart, the two c4-type zinc fingers in kni are separated by amino acids. this increased separation may in fact allow the two zinc fingers to act as functionally distinct dna binding domains, as is the case for hb.  taken together these data suggest that current pwms may not be optimal to explain the complexity of tf binding. although we only test four tfs in this study, we demonstrate that all four tfs exhibit context-dependent biases towards given nucleotides both contiguously with the defined minimal binding region and non-contiguously in flanking dna sequences, thus providing a foundation for this to be explored more broadly. an additional intriguing question for future study will be to investigate if the context-dependent bias persists at predicted tf binding regions in other genomic regions that are not characterized as crms. by taking these secondary context-dependencies into account, we propose that the information content of pwms can be expanded in many cases. this expansion would not only provide better predictions of true tf binding regions in the genome, but may also help improve estimates of relative binding affinities at specific sites, allowing one to understand the molecular basis for the difference between weak and strong binding sites. the ability to identify novel crms and decipher the sequence organization at crms relies heavily on a concrete understanding of tf binding preferences. improving the information content of pwms and our comprehension of tf binding events will contribute to these continued efforts. availability of supporting data the data sets supporting the results of this article are included within the article competing interests the authors declare that they have no competing interests. authors’ contributions jls created the web application, ran all bioinformatic and statistical analyses, and tested expanded hb pwms. asb compiled the crms, ran the evoprinterhd analysis, performed the jaspar database search, and analyzed each tf’s binding domain. the manuscript was written by jls, asb, rad and jmd, and the overall project was conceived and guided by rad and jmd. all authors read and approved the final manuscript. supplementary material additional file 1 dataset s dna sequences of the original crm. the original crm folder contains the fasta files for the crms we used. click here for file additional file 2 dataset s dna sequences of the extended crms. the extended crm folder contains those same crms, extended by including the neighboring bp both up- and downstream of the crm. the extended crm fasta filenames are identical to the original, but with a' + ’ appended to the filename. click here for file additional file 3: table s1 patser scores for predicted binding regions. each sheet contains a list of the sequences in predicted binding regions for a particular tf. sequences are ordered by their ln, as computed by patser. the scores highlighted in yellow indicate those scoring at or above less than or equal to) those included in the 75th percentile. the scores highlighted in green indicate those scoring at or above less than or equal to) those included in the 50th percentile. recall that the 75th and 50th percentile cutoff scores are specific for each tf and thus vary from sheet to sheet. click here for file additional file 4: table s2 nucleotide frequencies in flanking sequences neighboring predicted binding regions. each sheet corresponds to a particular tf and cutoff percentile. the column labeled'position’ assigns alphabetic values to the consensus core pwm-predicted binding regions and numeric values to the neighboring flanking sequences, ranging from ‒ to ‒ upstream and from to downstream. the next four columns, labeled'a’,'c’,'g’, and't’, contain the frequency observed for each nucleotide at the specified position when all pwm-pbrs are considered for that particular tf and cutoff percentile. the column labeled'chi‒squared’ contains the corresponding chi-squared value for each position. click here for file additional file 5: table s3 patser scores for extended pwm-pbrs. each sheet contains a list of the sequences in the extended pwm-pbrs for a particular initial cutoff score and length of flanking region. sequences are ordered by their ln, as computed by patser using the corresponding extended hb pwm. the scores highlighted in yellow indicate those scoring at or above less than or equal to) those included in the 75th percentile. the scores highlighted in green indicate those scoring at or above less than or equal to) those included in the 50th percentile. the scores highlighted in orange indicate those scoring at or above less than or equal to) those included in the 25th percentile. the scores highlighted in blue indicate those scoring at or above less than or equal to) those included in the 0th percentile. this includes all sites detected with the initial cutoff score. recall that the 0th, 25th, 50th and 75th percentile cutoff scores are specific for each initial cutoff score and length of flanking region and thus vary from sheet to sheet. click here for file additional file 6: table s4 true and false positive results for extended pwm-pbrs on chip-seq peaks. each sheet corresponds to a particular tf. the name of the tf, the total number of chip-seq peaks, and the number of scrambled peaks used for the analysis are listed in the first column. in the other columns, each row corresponds to a patser run for the given tf, chip-seq peaks and scrambled peaks. the primary and secondary cutoff percentiles used, as well as the number of nucleotides extended to the left and right of the core pwm-pbr are listed. the results of each run are shown as the number of chip peaks containing at least one predicted binding site, which is used to calculate the true positive rate, and the number of scrambled chip peaks containing at least one predicted binding site, which is used to calculate the false positive rate. the last column contains the ratio of true positive to false positive. to highlight settings in which the sensitivity increased in at least one of the extended pwm-pbrs compared to the core, we have highlighted those true positive rates in green, to highlight settings in which the specificity increased in at least one of the extended pwm-pbrs compared to the core, we have highlighted those false positive rates in pink, and to highlight settings in which the ratio of true to false positives increased in at least one of the extended pwm-pbrs compared to the core, we have highlighted those ratios in yellow. click here for file acknowledgments the research in this paper was supported by funding to r.a.d. from the national institutes of health, the national science foundation and howard hughes medical institute undergraduate science education program grants to the biology department at harvey mudd college. j.m.d. was funded as a teaching and research postdoctoral fellow, supported in part by nsf grant dms- and through research funds provided by amherst college.
expanding pwms to include sequence context-dependence will increase the information content in pwms and facilitate a more efficient functional identification and dissection of crms. in this proof-of-principle study, we analyze crms and focus on four tfs that control transcription of target genes along the anterio-posterior axis of the embryo early in development. gene expression in the drosophila embryo is controlled by functional interactions between a large network of protein transcription factors and specific sequences in dna cis-regulatory modules.
 gene expression in the drosophila embryo is controlled by functional interactions between a large network of protein transcription factors and specific sequences in dna cis-regulatory modules. the binding site sequences for any tf can be experimentally determined and represented in a position weight matrix. pwms can then be used to predict the location of tf binding sites in other regions of the genome, although there are limitations to this approach as currently implemented.  in this proof-of-principle study, we analyze crms and focus on four tfs that control transcription of target genes along the anterio-posterior axis of the embryo early in development. for all four of these tfs, there is some degree of conserved flanking sequence that extends beyond the predicted binding regions. a potential role for these conserved flanking sequences may be to enhance the specificity of tf binding, as the abundance of these sequences is greatly diminished when we examine only predicted high-affinity binding sites.  expanding pwms to include sequence context-dependence will increase the information content in pwms and facilitate a more efficient functional identification and dissection of crms. transcription factorbinding siteposition weight matrixenhancercis-regulatory moduledrosophila 