BACKGROUND
the availability of genome sequences of closely related species  has provided opportunities to solve several key biological problems such as the inference of phylogenetic trees, reconstruction of ancestral genomes, estimation of evolutionary rates, identification of conserved and non-conserved regions, and more generally the study of genome structure and evolution.
the alignment of multiple sequences, highlighting regions of homology among the sequences and predicting nucleotide level relationships among them, plays a critical role in such analyses.
numerous attempts have been made to develop accurate and efficient methods to solve the multiple sequence alignment problem , offering us much flexibility, as well as difficulty, in choosing the most appropriate tool for the task.
another important task related to multiple alignment is the annotation of insertions and deletions  in the alignment, a task that has received some attention in recent years  <cit>  in light of the realization that indels may be responsible for genomic variation as much as nucleotide substitutions are  <cit> , and that indels may affect regional mutation rates  <cit> .
given the availability of multiple tools to perform either of these two tasks, a researcher faces two important questions: "which of the tools should i use for my task?"
and "how accurate will the tool be on my data?"
answers to these come from studies that use data sets  where the true answers are known, to evaluate and compare different tools.
the design of benchmarks therefore directly affects the reliability of bioinformatics analyses that use those tools.
the two most widely used benchmarking approaches for alignment tools are  to make use of biological sequences and their manually curated alignments from databases such as homstrad  <cit> , balibase  <cit> , and sabmark  <cit> , or  to simulate the evolution of biological sequences by using specialized tools such as dawg  <cit> , rose  <cit>  and indelible  <cit> .
the main advantage of the former approach is the use of real biological sequences and alignments that are produced by using protein structure information.
this approach does not apply to non-coding dna sequences, whose alignments form the basis of regulatory comparative genomics.
therefore, simulation-based benchmarks have been widely adopted in this context  <cit> .
the simulation approach, however, is highly dependent on its parameters that reflect the underlying evolutionary processes and their rates.
it is not clear how to choose "correct" settings for these parameters and how to assess if the simulated sequences mimic real data well enough for claims about alignment accuracy, both in relative terms  and in the absolute, to generalize from the benchmarks to the real world setting.
we address these questions in this work, whose main contributions are the following.
1) we present a new simulation-based benchmarking method that is based on the entire spectrum of values of its parameters as inferred from real data.
this is in contrast to existing approaches that rely on the average observed values of the parameters.
2) we quantify the difficulty of aligning a data set by leveraging recent developments  <cit>  on estimating alignment accuracy without requiring the "true" alignments.
we reason that if the synthetic data sets truly mimic real orthologous sequences, the difficulty of aligning them ought to match that for the real data.
this is the key insight used to determine how realistic a particular benchmark  is, and we use this idea to show that the novel simulation method produces far more realistic benchmarks than the existing approach.
3) using our new benchmarks, we evaluate and compare the accuracy of six multiple alignment tools  on drosophila non-coding sequences.
the specific alignment task we consider is that of global alignment of ~1- <dig> kbp long sequences, and our conclusions may not apply to the task of local alignment, which was studied in  <cit> .
we are able to estimate the accuracy of alignment for specific sets of drosophila genomes, and find these to be very different from previously reported values.
we also evaluate two schemes for annotating insertions and deletions specifically, and find their accuracy to be comparable, and close to optimal.
4) we find that data sets with an excess of deletions over insertions are more amenable to accurate alignment than those with an excess of insertions, suggesting an implicit bias  with respect to their treatment of indels, even though none of the evaluated tools explicitly makes a distinction between insertions and deletions.
RESULTS
simulation of non-coding sequences by a traditional method
modeling of dna sequence evolution has been studied extensively in the past, and state-of-the-art simulation programs  <cit>  draw on various aspects of such models.
simulation of non-coding sequences  <cit>  incorporates current understanding of the architecture of such sequences in terms of regions of evolutionary constraint, for example by stipulating the presence of short  subsequences that evolve at a much slower rate than the rest of the sequence.
we refer the reader to  <cit>  for a comprehensive description of these approaches, which form the foundation of our own work reported here.
these simulation programs rely crucially on the values of their parameters .
the parameters serve to fully specify the stochastic processes from which evolutionary events  will be sampled, and prescribe the expected frequency of those events in the generated data sets.
variation in the frequency of these events, which underlie the difficulty of alignment tasks, results from the inherent randomness of the simulation process, i.e., the differences in random choices made from one "run" of the process to another.
it is natural to ask if the resulting variability across data sets in a synthetic benchmark is comparable to the corresponding variability observed in real orthologous sequences.
the question is particularly relevant due to the heterogeneity of non-coding sequences with respect to the density of functional elements and also motivated by the known variation in evolutionary rates across loci  <cit> .
we began by implementing the above-mentioned simulation paradigm, which we call the "traditional" paradigm, by incorporating the "constraint blocks" idea of pollard et al.
<cit>  into the dawg simulation program  <cit> .
parameters, including phylogeny, branch lengths, indel frequency, and various parameters related to conserved blocks were set based on previously published values from the literature  <cit>  or estimated by us from published multiple alignments of drosophila non-coding sequences .
a key difference in our implementation was that branch lengths  were estimated from non-coding sequences themselves, instead of synonymous substitution rates from coding sequences, as has been done previously.
we elaborate on this important issue later in this section.
we considered the alignments of real drosophila sequences from eight species , computed the sum of branch lengths of the phylogenetic tree estimated from ~ <dig> kbp segments of alignment, and found the distribution of this statistic to have a large variance across the genome .
the same distribution, when computed from  <dig> synthetic data sets generated using the traditional simulator described above, and the same alignment program, shows a very sharp peak around the mean .
we note that the means of the two distributions are similar , since the benchmark was parameterized by the average substitution rates observed in real data.
this is the first clear evidence that existing simulators fall short of representing the range of conservation levels in real data.
since substitution rates are generally correlated with indel rates, a large variance in the former implies a corresponding variance in indel frequencies, which of course lie at the root of the alignment problem.
this suggests that if we could measure the "difficulty of alignment" in any region of the genome , we ought to see a large variability in this measure across the genome.
moreover, if the observed distribution of the alignment difficulty measure is comparable to that in a benchmark, we would be confident in making claims about performance of alignment tools based on that benchmark.
the problem is that measuring alignment difficulty on real data requires knowledge of their true alignment, which is unavailable.
recent work by landan and graur  <cit>  showed that a reasonable surrogate for the accuracy of an alignment program on a data set can be computed even without the true alignment.
they reasoned that good alignments should be invariant to the orientation of the input sequences, and therefore defined the "heads or tails " alignment quality score as the agreement between two alignments, one generated from original sequences and the other from their reversed versions.
hall  <cit>  later showed that there is a clear positive correlation between hot alignment quality scores and the real alignment accuracy measured by comparison with the true alignment.
this remarkable finding inspired us to formulate the following strategy for quantifying the spectrum of alignment difficulty in data sets.
we computed the hot alignment quality score on the computed alignment of a data set, and used this score as a surrogate for the alignment difficulty of the data set.
low values of the alignment quality score indicate that the data set is particularly hard to align, and high values are suggestive of an "easy" data set.
as shown in figure 2a, the distributions of the score were significantly different between synthetic and real data sets.
alignment quality scores for 83% of the synthetic sequences are above  <dig>  whereas close to 50% of real sequences had scores below this range.
this strongly suggests that by and large the synthetic sequences simulated by the traditional approach are easier to align than real sequences, even though the former were generated with evolutionary parameters mirroring their real data counterparts.
in particular, the variance of alignment quality  is much smaller in synthetic data sets.
simulation based on a mixture model of parameters
we hypothesized that the above observation about synthetic data sets was due to the use of a single setting of the branch lengths, and the relatively low variability resulting from the randomness of the process itself .
if this is true, then one way to alleviate the problem would be to allow for multiple phylogenies for simulation of different data sets, with the variability of branch lengths across phylogenies introducing an additional source of data set variability.
we therefore considered a set of k =  <dig> phylogenies {ϕ <dig>  ϕ <dig>  ..., ϕk} that are scaled versions of the original phylogeny ϕ <dig>  i.e., every branch length in phylogeny ϕi is a constant factor τi times the corresponding branch length in ϕ <dig>   we modified the simulator to first sample at random one of the k phylogenies, and simulate according to this setting of branch lengths, with all other parameters being fixed as before.
in other words, the distribution of alignment quality scores from the new simulation process is a mixture distribution, with components parameterized by different phylogenies and the probability of sampling any particular phylogeny being the mixture weight.
we estimated an upper bound on the agreement between this mixture distribution and the observed distribution of alignment quality scores, by maximum likelihood training of mixture weights, through expectation-maximization algorithm  <cit> .
this "best fit" mixture distribution is shown in figure 2b, along with the real data distribution, and reveals a much stronger agreement between the two distributions, as compared to figure 2a.
the same trend was seen when allowing for a set of values of the "substitution to indel ratio" parameter , keeping all other parameters, including the phylogeny, fixed .
these results strongly suggested that the use of a range of parameter values instead of a single value has great impact on the variability of alignment difficulty in synthetic data sets, and has the potential to lead to the generation of realistic sequences.
simulation based on parameter sampling
the above results, while encouraging in terms of better reproducing the genomic variability of alignment difficulty, were obtained by fitting parameters of the simulation process so as to best match real data.
we next asked if we could achieve the same or better agreement between the synthetic and real data distributions without having seen the real distribution of alignment quality scores.
this would then allow us to use the observed agreement as a relatively unbiased assessment of how realistic the benchmark is.
developing the mixture model idea from the previous section, we now computed for each parameter the entire distribution of values observed in real data alignments, just as the traditional approach estimates the average of these values.
the simulation process was now made to sample each parameter independently from its empirical distribution, and then generate a data set based on the sampled parameter values.
the benchmark thus constructed  was examined for its distribution of alignment quality scores, and as seen in figure 2d, this distribution was remarkably close to that observed in real sequences.
in other words, the newly constructed benchmark meets our pre-specified criterion for a "realistic" benchmark.
the above analysis was performed using the sum-of-pairs score , which is the simplest of the scores defined in the hot approach  <cit> .
we repeated all analyses with another score, called the hot column score , and observed the same trends , although the agreement between synthetic and real data distributions was not as strong now as with the sps  .
assessment of multiple alignment tools
accuracy of multiple alignments
we used our new benchmark to evaluate and compare six leading multiple alignment tools that are publicly available and can align dna sequences.
these are clustalw  <dig> .
<dig>  <cit> , dialign-tx  <dig> .
<dig>  <cit> , mafft  <dig>   <cit> , mavid  <dig>  build  <dig>  <cit> , mlagan  <dig>   <cit> , and pecan  <dig>   <cit> .
we performed the assessment with varying numbers of species, k =  <dig>  ...,  <dig>  for each choice of k,  <dig> sets of sequences corresponding to k different drosophila species were simulated and the above alignment tools were run with default parameters or with the best setting recommended by their authors.
we then compared the resulting alignments to the "true" alignments reported by the simulation program, using the following three commonly used evaluation measures  <cit> :  alignment agreement, which is the fraction of aligned base pairs  in the predicted alignment that agree with the true alignment,  alignment sensitivity, which is the fraction of aligned base pairs of the true alignment that agree with the predicted alignment, and  alignment specificity, which is the fraction of aligned base pairs of the predicted alignment that agree with the true alignment.
whereas the alignment agreement score considers aligned base pairs as well as bases aligned to gaps, the sensitivity and specificity scores are calculated only from aligned base pairs.
the results of our evaluations are shown in figure  <dig> and additional files  <dig> and  <dig>  .
the pecan alignment program was found to be superior by all three measures, across all values of k. its performance degrades more slowly  than the other tools, as a result of which the gap between pecan and the other tools became larger more species were included in the tests.
the average alignment agreement in five species alignments produced by pecan  was close to 80%, but degraded to ~67% when aligning all eight species.
we performed the same evaluations by limiting ourselves to those data sets  that had an excess of insertions over deletions, and separately to those data sets with an excess of deletions .
surprisingly, we saw a clear difference between these two classes of data sets, with most tools performing significantly worse when there was an excess of insertions in the data set.
for example, on data sets with k =  <dig>  clustalw showed an alignment agreement of 36% or 46% depending on whether there was an excess of insertions or deletions .
the same trend was seen in terms of the alignment sensitivity and specificity measures.
noticably, pecan was largely unaffected by this dichotomy of data sets.
the evaluation measures used above consider all pairs of species in the k-species alignment and sum the accuracy values obtained from all pairs, without regard to the varying divergences of different pairs.
in an attempt to address this issue, we separately measured the alignment accuracy of different pairs of species , limiting ourselves to the eight-species data sets.
all trends reported above were also seen in this alternative view of the results .
the alignment agreement, using pecan, for d. melanogaster with d. yakuba, d. anannassae, d. pseudoobscura and d. willistoni was found to be 96%, 77%, 71% and 60% respectively.
disagreement with estimates based on existing benchmark
we found a substantial disagreement between our performance estimates and those previously reported by pollard et al.
<cit>  using their own benchmark.
for instance, the alignment sensitivity for the d. melanogaster - d. pseudoobscura pair comes out to be ~70% in our assessment and ~40% by their estimates, using the mlagan alignment tool.
we observe such gaps  also for alignment specificity, and for other species pairs and alignment programs as well .
while this discordance could be in part due to the fact that our benchmark employs a spectrum of parameter values to achieve greater and more realistic variability, we believe the major difference here is that even the average substitution rate, a key parameter in both simulation programs, is widely different between their study and ours.
the estimate used by pollard et al.
<cit>   is based on silent positions in codons, while our estimate  reflects the average subsitution frequency  seen in non-coding sequences.
in light of the results of figure 2d, where we show that our benchmark accurately mirrors the range of alignment difficulty in real data, the use of non-coding sequences in estimating this key parameter seems better justified.
we investigated this issue with additional tests.
we collected data sets representing the d. melanogaster - d. pseudoobscura pair from pollard et al.
<cit> , as well as from our benchmark and the real genomes.
the alignment quality score  distributions were computed for each type of benchmark, and are shown in figure  <dig>  we observed a close agreement between our data sets and the real orthologous sequences, while the pollard et al.
<cit>  data sets were harder to align on average, consistent with the greater substitution rate used there.
as noted in methods, the overall substitution frequency observed in non-coding sequences may be viewed as an average of the corresponding frequency in conserved blocks and the much higher frequency outside conserved blocks.
this average is determined by two key parameters α, the fraction of sequence length that falls into conserved blocks, and β, the ratio of the evolutionary rate of conserved blocks to that outside blocks.
given that the divergence estimate used by pollard et al.
<cit>  for these two species is ~ <dig>   substitutions per site, if we are to treat this value as the neutral rate  in non-coding sequences, what values of α and β would lead to the observed overall substitution frequency of  <dig> ?
we determined that if β =  <dig> , as was used by pollard et al.
<cit>  , α has to be ~ <dig> , i.e., about 92% of non-coding sequences have to be conserved blocks, which is far higher than most current estimates of this parameter  <cit> .
similarly, if we are to trust the values of α =  <dig>  and β =  <dig> , as was used by pollard et al.
<cit>  , then the overall divergence, after averaging between conserved blocks and non-blocks, would be ~ <dig>  substitutions per site, far greater than what is observed .
we therefore concluded that the use of synonymous substitution rates as the neutral rate for non-coding sequence is likely to lead to benchmarks with overly "diverged" sequences that are more difficult to align than real sequences from those species.
assessment of indel annotation schemes
traditional alignment programs mark the predicted locations of insertions and deletions as "gaps", and do not proceed to annotate these gaps as being insertions or deletions.
this latter task has received some attention recently with at least two "indel annotation schemes" being published, based on maximum-parsimony  and probabilistic-models  respectively.
we examined the accuracy of these two alignment-related tools on our new benchmark.
we noted that the best alignment agreement score  is ~70% for d. melanogaster - d. pseudoobscura, and decreases to ~60% when a more diverged species  is added.
reasoning that phylogenies for which computed alignments are largely inaccurate would not be suitable for insertion/deletion annotation in any case, we chose to limit our assessment to the following five drosophila species: d. melanogaster, d. simulans, d. yakuba, d. ananassae, and d. pseudoobscura .
the "true" alignment  was provided to the two indel annotation tools and the insertion/deletion annotations on each of the five terminal branches  of the phylogeny were compared to the "true" annotations.
the following three measures were used for assessment, borrowed from  <cit> :  indel count agreement, which is the agreement of indel counts between true and predicted annotations,  indel ratio agreement, which is the agreement of the ratio of the number of insertions to the total number of indels between the two annotations, and  indel annotation coverage, which is the fraction of indel positions on which the two annotations agree .
as summarized in table  <dig>  indel count agreement scores of the two tools were very similar to each other and close to optimal  for most species except d. pseudoobscura, the species with the longest terminal branch in the phylogeny.
indel ratio agreement scores of both tools were close to optimal  in all five species.
while the sensitivity scores of indel annotation coverage of the two tools were above 90% across all five species, the specificity scores were above 90% only for the four species except d. pseudoobscura.
the loss of accuracy on the d. pseudoobscura branch is presumably due to the fact that there is no "outgroup" species to aid disambiguation of insertions and deletions on this branch.
we further discuss the implications of these observations in the next section.
we also repeated our assessment for sequences with an excess of insertions or of deletions, as above, but no significant differences was observed between these two categories .
aindel count agreement
bindel ratio agreement
cindel annotation coverage
discussion
choosing the most suitable tool for aligning orthologous sequences is essential to studies in comparative genomics and in molecular evolution, making it critical to develop accurate benchmarking methodology.
in this study, we propose a novel simulation-based approach to generate realistic data sets mimicking orthologous non-coding sequences from multiple drosophila species.
this new simulation method exploits the spectrum of values of evolutionary statistics  seen across a genome.
we take advantage of an objective "alignment quality" measure to show that the synthetic sequences produced agree with real sequences not only in terms of evolutionary statistics, but are also as easy or hard to align as real data sets.
in this sense, our evaluation results are more likely to reflect the actual accuracy values of alignment-related tools on data from drosophila species.
we note that our strategy of sampling parameters  from their empirical distribution has parallels with traditional bayesian inference where one integrates over  a prior distribution on parameters, rather than using a single point estimate.
a key step in our benchmark construction was the ability to assess the quality of an alignment without access to the corresponding true alignment.
this ability has been the result of several recent publications by other authors.
prakash and tompa  <cit>  developed statistical methods to assess if a multiple sequence alignment appears contaminated with one or more unrelated sequences, based on which they identified regions of whole genome alignments as being suspect.
the development of the "hot" method by landan and graur  <cit>  then came as a breakthrough to assess the reliability of multiple sequence alignments.
later on, landan and graur  <cit>  extended the hot method to take advantage of co-optimal alternative alignments generated by progressive alignment tools.
however, the implementation of this method is too dependent on the specific procedures of a progressive alignment method, making the original hot score  <cit>  a natural choice for our purpose.
while our benchmark is shown to be very close to real sequences in terms of the distribution of hot sps, we are cautioned by the discrepancy observed between simulated and real sequences in terms of the hot cs, an alternative alignment quality score from the same authors .
this is likely the product of properties of non-coding sequences that are not adequately represented in our simulation process.
for example, modeling the functional constraints embedded in non-coding sequences through short conserved blocks  is surely an oversimplification of the complexity of genomic architecture.
important progress has been made on this front, in the form of specialized evolutionary simulators that model transcription factor binding site evolution in realistic ways  <cit> .
each of these simulators makes specific assumptions about cis-regulatory architecture, vis-a-vis the density and evolution of binding sites.
however, it is not yet clear which, if any, of these different assumed models of regulatory sequence evolution is most suited to represent the variability in constraint patterns across different regions of the genome.
our simplistic "conserved block" model  seems to be a good approximation that captures the most prominent patterns in orthologous non-coding sequences, in terms of alignment difficulty.
we expect that future research on more realistic models of cis-regulatory architecture will lead us to replace the alternating arrangement of conserved blocks and faster evolving segments with a pattern more in line with reality.
future work may also include careful modeling of genomic repeats and repeat generating evolutionary events, since repeat-rich genomes may present additional challenges for the alignment task.
our proposed framework of sampling evolutionary parameters before running the simulation process will remain equally important in future benchmarks that implement such sophisticated models.
some clarification is in order with respect to our manner of choosing substitution rates for the simulation process, since it marks a significant departure from traditional thinking.
the latter, as embodied in the work of pollard et al.
<cit> , prescribes that the "unconstrained" parts of the sequence evolve with nucleotide substitution rate equal to that infered from synonymous mutations in the nearby gene .
this rate  is widely different from the value observed in real non-coding sequence alignments .
one could argue that this gap may be offset if we set an appropriate frequency of conserved positions , resulting in an average substitution rate that is close to the empirically observed value.
however, this turned out not be the case for any realistic setting of the frequency of conserved positions .
we therefore chose to be guided by existing estimates of the frequency and length distribution of conserved blocks, with substitution rates that are some constant β  times the "neutral" rate outside of the blocks, and set this neutral rate so that the average rate for the entire sequence matches observed values.
our choice reflects the philosophy that simulated data sets ought to match real data in terms of various evolutionary statistics and net alignment difficulty, and the discordance of the used neutral rate from synonymous substitution rates is ignored for the sake of practicality.
to our knowledge, no previous benchmarking study has evaluated the effect of insertions and deletions on the performance of alignment tools.
some studies  <cit>  have used equal frequencies for insertions and deletions and focused on the collective effects of indels.
here, we attempted to elucidate the differing effects of insertions and deletions by separately summarizing results for the two extreme cases where the number of insertions is at least two times the frequency of deletions and vice versa.
the results were surprising, and indicated that most multiple alignment tools find it harder to accurately align data sets with an excess of insertions than those with more deletions .
löytynoja and goldman  <cit>  offered valuable insight into a possible source of this asymmetry, pointing out that progressive alignment methods  "end up penalizing single insertion events multiple times".
we speculate therefore, as they did, that claims about insertion/deletion frequencies along the genome should be preceded by an examination of the alignment method's accuracy in regimes of high insertion frequency.
finally, a note about our findings on insertion/deletion annotation.
indelign  <cit>  is a probabilistic tool that annotates insertions and deletions by maximum likelihood training of an evolutionary model.
sbinfer  <cit>  is a greedy algorithm that reconstructs ancestral sequences based on the maximum parsimony principle, and therefore allows us to infer insertion/deletion annotations.
to assess these two tools without being confounded by errors of an alignment program, we examined their performance on the true alignments.
we found the two programs to have comparable accuracy on our benchmark for the five drosophila species.
while the accuracy was close to optimal on four of the five terminal branches, we observed that both tools over-estimate insertions as well as deletions on the longest branch , while accurately predicting the ratio of insertions to deletions.
we note that the d. pseudoobscura branch in the phylogenetic tree originates from the root of the tree, and we would expect to have better annotation results for this branch if an appropriate outgroup species was used.
for studies that intend to use insertion to deletion ratio profiling to identify loci with unusual evolutionary patterns  it may be safe to examine all five terminal branches of this tree; however, for the more common requirement of accurately annotating insertion and deletion events, e.g., to study gain and loss patterns of specific classes of transcription factor binding sites  <cit> , we do not recommend using events on the d. pseudoobscura branch.
CONCLUSIONS
we have presented a novel method for generating benchmarks of non-coding sequence alignments, that relies on a spectrum of parameter values reflecting the genome-wide variation of those parameters.
we have shown our benchmarks to accurately match the difficulty of aligning real data, by taking advantage of recent developments in measurement of alignment quality.
benchmark evaluations on drosophila non-coding sequences suggest a greater accuracy of multiple alignment tools  than previously reported, and points to a clear asymmetry in the handling of insertions versus deletions by most alignment tools.

 the availability of genome sequences of closely related species has provided opportunities to solve several key biological problems such as the inference of phylogenetic trees, reconstruction of ancestral genomes, estimation of evolutionary rates, identification of conserved and non-conserved regions, and more generally the study of genome structure and evolution. the alignment of multiple sequences, highlighting regions of homology among the sequences and predicting nucleotide level relationships among them, plays a critical role in such analyses. numerous attempts have been made to develop accurate and efficient methods to solve the multiple sequence alignment problem, offering us much flexibility, as well as difficulty, in choosing the most appropriate tool for the task. another important task related to multiple alignment is the annotation of insertions and deletions in the alignment, a task that has received some attention in recent years in light of the realization that indels may be responsible for genomic variation as much as nucleotide substitutions are, and that indels may affect regional mutation rates. given the availability of multiple tools to perform either of these two tasks, a researcher faces two important questions: "which of the tools should i use for my task?" and "how accurate will the tool be on my data?" answers to these come from studies that use data sets where the true answers are known, to evaluate and compare different tools. the design of benchmarks therefore directly affects the reliability of bioinformatics analyses that use those tools. the two most widely used benchmarking approaches for alignment tools are to make use of biological sequences and their manually curated alignments from databases such as homstrad, balibase, and sabmark, or to simulate the evolution of biological sequences by using specialized tools such as dawg, rose and indelible. the main advantage of the former approach is the use of real biological sequences and alignments that are produced by using protein structure information. this approach does not apply to non-coding dna sequences, whose alignments form the basis of regulatory comparative genomics. therefore, simulation-based benchmarks have been widely adopted in this context. the simulation approach, however, is highly dependent on its parameters that reflect the underlying evolutionary processes and their rates. it is not clear how to choose "correct" settings for these parameters and how to assess if the simulated sequences mimic real data well enough for claims about alignment accuracy, both in relative terms and in the absolute, to generalize from the benchmarks to the real world setting. we address these questions in this work, whose main contributions are the following. 1) we present a new simulation-based benchmarking method that is based on the entire spectrum of values of its parameters as inferred from real data. this is in contrast to existing approaches that rely on the average observed values of the parameters. 2) we quantify the difficulty of aligning a data set by leveraging recent developments on estimating alignment accuracy without requiring the "true" alignments. we reason that if the synthetic data sets truly mimic real orthologous sequences, the difficulty of aligning them ought to match that for the real data. this is the key insight used to determine how realistic a particular benchmark is, and we use this idea to show that the novel simulation method produces far more realistic benchmarks than the existing approach. 3) using our new benchmarks, we evaluate and compare the accuracy of six multiple alignment tools on drosophila non-coding sequences. the specific alignment task we consider is that of global alignment of ~1- kbp long sequences, and our conclusions may not apply to the task of local alignment, which was studied in. we are able to estimate the accuracy of alignment for specific sets of drosophila genomes, and find these to be very different from previously reported values. we also evaluate two schemes for annotating insertions and deletions specifically, and find their accuracy to be comparable, and close to optimal. 4) we find that data sets with an excess of deletions over insertions are more amenable to accurate alignment than those with an excess of insertions, suggesting an implicit bias with respect to their treatment of indels, even though none of the evaluated tools explicitly makes a distinction between insertions and deletions.  simulation of non-coding sequences by a traditional method modeling of dna sequence evolution has been studied extensively in the past, and state-of-the-art simulation programs draw on various aspects of such models. simulation of non-coding sequences incorporates current understanding of the architecture of such sequences in terms of regions of evolutionary constraint, for example by stipulating the presence of short subsequences that evolve at a much slower rate than the rest of the sequence. we refer the reader to for a comprehensive description of these approaches, which form the foundation of our own work reported here. these simulation programs rely crucially on the values of their parameters. the parameters serve to fully specify the stochastic processes from which evolutionary events will be sampled, and prescribe the expected frequency of those events in the generated data sets. variation in the frequency of these events, which underlie the difficulty of alignment tasks, results from the inherent randomness of the simulation process, i.e., the differences in random choices made from one "run" of the process to another. it is natural to ask if the resulting variability across data sets in a synthetic benchmark is comparable to the corresponding variability observed in real orthologous sequences. the question is particularly relevant due to the heterogeneity of non-coding sequences with respect to the density of functional elements and also motivated by the known variation in evolutionary rates across loci. we began by implementing the above-mentioned simulation paradigm, which we call the "traditional" paradigm, by incorporating the "constraint blocks" idea of pollard et al. into the dawg simulation program. parameters, including phylogeny, branch lengths, indel frequency, and various parameters related to conserved blocks were set based on previously published values from the literature or estimated by us from published multiple alignments of drosophila non-coding sequences. a key difference in our implementation was that branch lengths were estimated from non-coding sequences themselves, instead of synonymous substitution rates from coding sequences, as has been done previously. we elaborate on this important issue later in this section. we considered the alignments of real drosophila sequences from eight species, computed the sum of branch lengths of the phylogenetic tree estimated from ~ kbp segments of alignment, and found the distribution of this statistic to have a large variance across the genome. the same distribution, when computed from synthetic data sets generated using the traditional simulator described above, and the same alignment program, shows a very sharp peak around the mean. we note that the means of the two distributions are similar, since the benchmark was parameterized by the average substitution rates observed in real data. this is the first clear evidence that existing simulators fall short of representing the range of conservation levels in real data. since substitution rates are generally correlated with indel rates, a large variance in the former implies a corresponding variance in indel frequencies, which of course lie at the root of the alignment problem. this suggests that if we could measure the "difficulty of alignment" in any region of the genome, we ought to see a large variability in this measure across the genome. moreover, if the observed distribution of the alignment difficulty measure is comparable to that in a benchmark, we would be confident in making claims about performance of alignment tools based on that benchmark. the problem is that measuring alignment difficulty on real data requires knowledge of their true alignment, which is unavailable. recent work by landan and graur showed that a reasonable surrogate for the accuracy of an alignment program on a data set can be computed even without the true alignment. they reasoned that good alignments should be invariant to the orientation of the input sequences, and therefore defined the "heads or tails " alignment quality score as the agreement between two alignments, one generated from original sequences and the other from their reversed versions. hall later showed that there is a clear positive correlation between hot alignment quality scores and the real alignment accuracy measured by comparison with the true alignment. this remarkable finding inspired us to formulate the following strategy for quantifying the spectrum of alignment difficulty in data sets. we computed the hot alignment quality score on the computed alignment of a data set, and used this score as a surrogate for the alignment difficulty of the data set. low values of the alignment quality score indicate that the data set is particularly hard to align, and high values are suggestive of an "easy" data set. as shown in figure 2a, the distributions of the score were significantly different between synthetic and real data sets. alignment quality scores for 83% of the synthetic sequences are above whereas close to 50% of real sequences had scores below this range. this strongly suggests that by and large the synthetic sequences simulated by the traditional approach are easier to align than real sequences, even though the former were generated with evolutionary parameters mirroring their real data counterparts. in particular, the variance of alignment quality is much smaller in synthetic data sets. simulation based on a mixture model of parameters we hypothesized that the above observation about synthetic data sets was due to the use of a single setting of the branch lengths, and the relatively low variability resulting from the randomness of the process itself. if this is true, then one way to alleviate the problem would be to allow for multiple phylogenies for simulation of different data sets, with the variability of branch lengths across phylogenies introducing an additional source of data set variability. we therefore considered a set of k = phylogenies {ϕ ϕ..., ϕk} that are scaled versions of the original phylogeny ϕ i.e., every branch length in phylogeny ϕi is a constant factor τi times the corresponding branch length in ϕ we modified the simulator to first sample at random one of the k phylogenies, and simulate according to this setting of branch lengths, with all other parameters being fixed as before. in other words, the distribution of alignment quality scores from the new simulation process is a mixture distribution, with components parameterized by different phylogenies and the probability of sampling any particular phylogeny being the mixture weight. we estimated an upper bound on the agreement between this mixture distribution and the observed distribution of alignment quality scores, by maximum likelihood training of mixture weights, through expectation-maximization algorithm. this "best fit" mixture distribution is shown in figure 2b, along with the real data distribution, and reveals a much stronger agreement between the two distributions, as compared to figure 2a. the same trend was seen when allowing for a set of values of the "substitution to indel ratio" parameter, keeping all other parameters, including the phylogeny, fixed. these results strongly suggested that the use of a range of parameter values instead of a single value has great impact on the variability of alignment difficulty in synthetic data sets, and has the potential to lead to the generation of realistic sequences. simulation based on parameter sampling the above results, while encouraging in terms of better reproducing the genomic variability of alignment difficulty, were obtained by fitting parameters of the simulation process so as to best match real data. we next asked if we could achieve the same or better agreement between the synthetic and real data distributions without having seen the real distribution of alignment quality scores. this would then allow us to use the observed agreement as a relatively unbiased assessment of how realistic the benchmark is. developing the mixture model idea from the previous section, we now computed for each parameter the entire distribution of values observed in real data alignments, just as the traditional approach estimates the average of these values. the simulation process was now made to sample each parameter independently from its empirical distribution, and then generate a data set based on the sampled parameter values. the benchmark thus constructed was examined for its distribution of alignment quality scores, and as seen in figure 2d, this distribution was remarkably close to that observed in real sequences. in other words, the newly constructed benchmark meets our pre-specified criterion for a "realistic" benchmark. the above analysis was performed using the sum-of-pairs score, which is the simplest of the scores defined in the hot approach. we repeated all analyses with another score, called the hot column score, and observed the same trends, although the agreement between synthetic and real data distributions was not as strong now as with the sps. assessment of multiple alignment tools accuracy of multiple alignments we used our new benchmark to evaluate and compare six leading multiple alignment tools that are publicly available and can align dna sequences. these are clustalw., dialign-tx., mafft, mavid build, mlagan, and pecan. we performed the assessment with varying numbers of species, k =..., for each choice of k, sets of sequences corresponding to k different drosophila species were simulated and the above alignment tools were run with default parameters or with the best setting recommended by their authors. we then compared the resulting alignments to the "true" alignments reported by the simulation program, using the following three commonly used evaluation measures: alignment agreement, which is the fraction of aligned base pairs in the predicted alignment that agree with the true alignment, alignment sensitivity, which is the fraction of aligned base pairs of the true alignment that agree with the predicted alignment, and alignment specificity, which is the fraction of aligned base pairs of the predicted alignment that agree with the true alignment. whereas the alignment agreement score considers aligned base pairs as well as bases aligned to gaps, the sensitivity and specificity scores are calculated only from aligned base pairs. the results of our evaluations are shown in figure and additional files and. the pecan alignment program was found to be superior by all three measures, across all values of k. its performance degrades more slowly than the other tools, as a result of which the gap between pecan and the other tools became larger more species were included in the tests. the average alignment agreement in five species alignments produced by pecan was close to 80%, but degraded to ~67% when aligning all eight species. we performed the same evaluations by limiting ourselves to those data sets that had an excess of insertions over deletions, and separately to those data sets with an excess of deletions. surprisingly, we saw a clear difference between these two classes of data sets, with most tools performing significantly worse when there was an excess of insertions in the data set. for example, on data sets with k = clustalw showed an alignment agreement of 36% or 46% depending on whether there was an excess of insertions or deletions. the same trend was seen in terms of the alignment sensitivity and specificity measures. noticably, pecan was largely unaffected by this dichotomy of data sets. the evaluation measures used above consider all pairs of species in the k-species alignment and sum the accuracy values obtained from all pairs, without regard to the varying divergences of different pairs. in an attempt to address this issue, we separately measured the alignment accuracy of different pairs of species, limiting ourselves to the eight-species data sets. all trends reported above were also seen in this alternative view of the results. the alignment agreement, using pecan, for d. melanogaster with d. yakuba, d. anannassae, d. pseudoobscura and d. willistoni was found to be 96%, 77%, 71% and 60% respectively. disagreement with estimates based on existing benchmark we found a substantial disagreement between our performance estimates and those previously reported by pollard et al. using their own benchmark. for instance, the alignment sensitivity for the d. melanogaster - d. pseudoobscura pair comes out to be ~70% in our assessment and ~40% by their estimates, using the mlagan alignment tool. we observe such gaps also for alignment specificity, and for other species pairs and alignment programs as well. while this discordance could be in part due to the fact that our benchmark employs a spectrum of parameter values to achieve greater and more realistic variability, we believe the major difference here is that even the average substitution rate, a key parameter in both simulation programs, is widely different between their study and ours. the estimate used by pollard et al. is based on silent positions in codons, while our estimate reflects the average subsitution frequency seen in non-coding sequences. in light of the results of figure 2d, where we show that our benchmark accurately mirrors the range of alignment difficulty in real data, the use of non-coding sequences in estimating this key parameter seems better justified. we investigated this issue with additional tests. we collected data sets representing the d. melanogaster - d. pseudoobscura pair from pollard et al., as well as from our benchmark and the real genomes. the alignment quality score distributions were computed for each type of benchmark, and are shown in figure we observed a close agreement between our data sets and the real orthologous sequences, while the pollard et al. data sets were harder to align on average, consistent with the greater substitution rate used there. as noted in methods, the overall substitution frequency observed in non-coding sequences may be viewed as an average of the corresponding frequency in conserved blocks and the much higher frequency outside conserved blocks. this average is determined by two key parameters α, the fraction of sequence length that falls into conserved blocks, and β, the ratio of the evolutionary rate of conserved blocks to that outside blocks. given that the divergence estimate used by pollard et al. for these two species is ~ substitutions per site, if we are to treat this value as the neutral rate in non-coding sequences, what values of α and β would lead to the observed overall substitution frequency of? we determined that if β =, as was used by pollard et al., α has to be ~, i.e., about 92% of non-coding sequences have to be conserved blocks, which is far higher than most current estimates of this parameter. similarly, if we are to trust the values of α = and β =, as was used by pollard et al., then the overall divergence, after averaging between conserved blocks and non-blocks, would be ~ substitutions per site, far greater than what is observed. we therefore concluded that the use of synonymous substitution rates as the neutral rate for non-coding sequence is likely to lead to benchmarks with overly "diverged" sequences that are more difficult to align than real sequences from those species. assessment of indel annotation schemes traditional alignment programs mark the predicted locations of insertions and deletions as "gaps", and do not proceed to annotate these gaps as being insertions or deletions. this latter task has received some attention recently with at least two "indel annotation schemes" being published, based on maximum-parsimony and probabilistic-models respectively. we examined the accuracy of these two alignment-related tools on our new benchmark. we noted that the best alignment agreement score is ~70% for d. melanogaster - d. pseudoobscura, and decreases to ~60% when a more diverged species is added. reasoning that phylogenies for which computed alignments are largely inaccurate would not be suitable for insertion/deletion annotation in any case, we chose to limit our assessment to the following five drosophila species: d. melanogaster, d. simulans, d. yakuba, d. ananassae, and d. pseudoobscura. the "true" alignment was provided to the two indel annotation tools and the insertion/deletion annotations on each of the five terminal branches of the phylogeny were compared to the "true" annotations. the following three measures were used for assessment, borrowed from: indel count agreement, which is the agreement of indel counts between true and predicted annotations, indel ratio agreement, which is the agreement of the ratio of the number of insertions to the total number of indels between the two annotations, and indel annotation coverage, which is the fraction of indel positions on which the two annotations agree. as summarized in table indel count agreement scores of the two tools were very similar to each other and close to optimal for most species except d. pseudoobscura, the species with the longest terminal branch in the phylogeny. indel ratio agreement scores of both tools were close to optimal in all five species. while the sensitivity scores of indel annotation coverage of the two tools were above 90% across all five species, the specificity scores were above 90% only for the four species except d. pseudoobscura. the loss of accuracy on the d. pseudoobscura branch is presumably due to the fact that there is no "outgroup" species to aid disambiguation of insertions and deletions on this branch. we further discuss the implications of these observations in the next section. we also repeated our assessment for sequences with an excess of insertions or of deletions, as above, but no significant differences was observed between these two categories. aindel count agreement bindel ratio agreement cindel annotation coverage discussion choosing the most suitable tool for aligning orthologous sequences is essential to studies in comparative genomics and in molecular evolution, making it critical to develop accurate benchmarking methodology. in this study, we propose a novel simulation-based approach to generate realistic data sets mimicking orthologous non-coding sequences from multiple drosophila species. this new simulation method exploits the spectrum of values of evolutionary statistics seen across a genome. we take advantage of an objective "alignment quality" measure to show that the synthetic sequences produced agree with real sequences not only in terms of evolutionary statistics, but are also as easy or hard to align as real data sets. in this sense, our evaluation results are more likely to reflect the actual accuracy values of alignment-related tools on data from drosophila species. we note that our strategy of sampling parameters from their empirical distribution has parallels with traditional bayesian inference where one integrates over a prior distribution on parameters, rather than using a single point estimate. a key step in our benchmark construction was the ability to assess the quality of an alignment without access to the corresponding true alignment. this ability has been the result of several recent publications by other authors. prakash and tompa developed statistical methods to assess if a multiple sequence alignment appears contaminated with one or more unrelated sequences, based on which they identified regions of whole genome alignments as being suspect. the development of the "hot" method by landan and graur then came as a breakthrough to assess the reliability of multiple sequence alignments. later on, landan and graur extended the hot method to take advantage of co-optimal alternative alignments generated by progressive alignment tools. however, the implementation of this method is too dependent on the specific procedures of a progressive alignment method, making the original hot score a natural choice for our purpose. while our benchmark is shown to be very close to real sequences in terms of the distribution of hot sps, we are cautioned by the discrepancy observed between simulated and real sequences in terms of the hot cs, an alternative alignment quality score from the same authors. this is likely the product of properties of non-coding sequences that are not adequately represented in our simulation process. for example, modeling the functional constraints embedded in non-coding sequences through short conserved blocks is surely an oversimplification of the complexity of genomic architecture. important progress has been made on this front, in the form of specialized evolutionary simulators that model transcription factor binding site evolution in realistic ways. each of these simulators makes specific assumptions about cis-regulatory architecture, vis-a-vis the density and evolution of binding sites. however, it is not yet clear which, if any, of these different assumed models of regulatory sequence evolution is most suited to represent the variability in constraint patterns across different regions of the genome. our simplistic "conserved block" model seems to be a good approximation that captures the most prominent patterns in orthologous non-coding sequences, in terms of alignment difficulty. we expect that future research on more realistic models of cis-regulatory architecture will lead us to replace the alternating arrangement of conserved blocks and faster evolving segments with a pattern more in line with reality. future work may also include careful modeling of genomic repeats and repeat generating evolutionary events, since repeat-rich genomes may present additional challenges for the alignment task. our proposed framework of sampling evolutionary parameters before running the simulation process will remain equally important in future benchmarks that implement such sophisticated models. some clarification is in order with respect to our manner of choosing substitution rates for the simulation process, since it marks a significant departure from traditional thinking. the latter, as embodied in the work of pollard et al., prescribes that the "unconstrained" parts of the sequence evolve with nucleotide substitution rate equal to that infered from synonymous mutations in the nearby gene. this rate is widely different from the value observed in real non-coding sequence alignments. one could argue that this gap may be offset if we set an appropriate frequency of conserved positions, resulting in an average substitution rate that is close to the empirically observed value. however, this turned out not be the case for any realistic setting of the frequency of conserved positions. we therefore chose to be guided by existing estimates of the frequency and length distribution of conserved blocks, with substitution rates that are some constant β times the "neutral" rate outside of the blocks, and set this neutral rate so that the average rate for the entire sequence matches observed values. our choice reflects the philosophy that simulated data sets ought to match real data in terms of various evolutionary statistics and net alignment difficulty, and the discordance of the used neutral rate from synonymous substitution rates is ignored for the sake of practicality. to our knowledge, no previous benchmarking study has evaluated the effect of insertions and deletions on the performance of alignment tools. some studies have used equal frequencies for insertions and deletions and focused on the collective effects of indels. here, we attempted to elucidate the differing effects of insertions and deletions by separately summarizing results for the two extreme cases where the number of insertions is at least two times the frequency of deletions and vice versa. the results were surprising, and indicated that most multiple alignment tools find it harder to accurately align data sets with an excess of insertions than those with more deletions. löytynoja and goldman offered valuable insight into a possible source of this asymmetry, pointing out that progressive alignment methods "end up penalizing single insertion events multiple times". we speculate therefore, as they did, that claims about insertion/deletion frequencies along the genome should be preceded by an examination of the alignment method's accuracy in regimes of high insertion frequency. finally, a note about our findings on insertion/deletion annotation. indelign is a probabilistic tool that annotates insertions and deletions by maximum likelihood training of an evolutionary model. sbinfer is a greedy algorithm that reconstructs ancestral sequences based on the maximum parsimony principle, and therefore allows us to infer insertion/deletion annotations. to assess these two tools without being confounded by errors of an alignment program, we examined their performance on the true alignments. we found the two programs to have comparable accuracy on our benchmark for the five drosophila species. while the accuracy was close to optimal on four of the five terminal branches, we observed that both tools over-estimate insertions as well as deletions on the longest branch, while accurately predicting the ratio of insertions to deletions. we note that the d. pseudoobscura branch in the phylogenetic tree originates from the root of the tree, and we would expect to have better annotation results for this branch if an appropriate outgroup species was used. for studies that intend to use insertion to deletion ratio profiling to identify loci with unusual evolutionary patterns it may be safe to examine all five terminal branches of this tree; however, for the more common requirement of accurately annotating insertion and deletion events, e.g., to study gain and loss patterns of specific classes of transcription factor binding sites, we do not recommend using events on the d. pseudoobscura branch.  we have presented a novel method for generating benchmarks of non-coding sequence alignments, that relies on a spectrum of parameter values reflecting the genome-wide variation of those parameters. we have shown our benchmarks to accurately match the difficulty of aligning real data, by taking advantage of recent developments in measurement of alignment quality. benchmark evaluations on drosophila non-coding sequences suggest a greater accuracy of multiple alignment tools than previously reported, and points to a clear asymmetry in the handling of insertions versus deletions by most alignment tools. 
we also examine the accuracy of two existing tools for annotating insertions versus deletions, and find their performance to be close to optimal in drosophila non-coding sequences if provided with the true alignments. we find that the conventional simulation approach, which relies on empirically estimated values for various parameters such as substitution rate or insertion/deletion rates, is unable to generate synthetic sequences reflecting the broad genomic variation in conservation levels. apart from helping to select the most effective tools, these benchmarks will help practitioners of comparative genomics deal with the effects of alignment errors, by providing accurate estimates of the extent of these errors. we tackle this problem with a new method for simulating non-coding sequence evolution, by relying on genome-wide distributions of evolutionary parameters rather than their averages. we evaluate six widely used multiple alignment tools in the context of drosophila non-coding sequences, and find the accuracy to be significantly different from previously reported values. this allows us to make significant progress towards estimating the alignment accuracy of current tools in an absolute sense, going beyond only a relative assessment of different tools. however, it is not clear if such benchmarks truly represent real sequence data from any given group of species, in terms of the difficulty of alignment tasks. we then generate synthetic data sets to mimic orthologous sequences from the drosophila group of species, and show that these data sets truly represent the variability observed in genomic data in terms of the difficulty of the alignment task. with the continued development of new computational tools for multiple sequence alignment, it is necessary today to develop benchmarks that aid the selection of the most effective tools. we have developed a method to generate benchmarks for multiple alignments of drosophila non-coding sequences, and shown it to be more realistic than traditional benchmarks.
 with the continued development of new computational tools for multiple sequence alignment, it is necessary today to develop benchmarks that aid the selection of the most effective tools. simulation-based benchmarks have been proposed to meet this necessity, especially for non-coding sequences. however, it is not clear if such benchmarks truly represent real sequence data from any given group of species, in terms of the difficulty of alignment tasks.  we find that the conventional simulation approach, which relies on empirically estimated values for various parameters such as substitution rate or insertion/deletion rates, is unable to generate synthetic sequences reflecting the broad genomic variation in conservation levels. we tackle this problem with a new method for simulating non-coding sequence evolution, by relying on genome-wide distributions of evolutionary parameters rather than their averages. we then generate synthetic data sets to mimic orthologous sequences from the drosophila group of species, and show that these data sets truly represent the variability observed in genomic data in terms of the difficulty of the alignment task. this allows us to make significant progress towards estimating the alignment accuracy of current tools in an absolute sense, going beyond only a relative assessment of different tools. we evaluate six widely used multiple alignment tools in the context of drosophila non-coding sequences, and find the accuracy to be significantly different from previously reported values. interestingly, the performance of most tools degrades more rapidly when there are more insertions than deletions in the data set, suggesting an asymmetric handling of insertions and deletions, even though none of the evaluated tools explicitly distinguishes these two types of events. we also examine the accuracy of two existing tools for annotating insertions versus deletions, and find their performance to be close to optimal in drosophila non-coding sequences if provided with the true alignments.  we have developed a method to generate benchmarks for multiple alignments of drosophila non-coding sequences, and shown it to be more realistic than traditional benchmarks. apart from helping to select the most effective tools, these benchmarks will help practitioners of comparative genomics deal with the effects of alignment errors, by providing accurate estimates of the extent of these errors. 