Benchmarking tools for the alignment of functional noncoding DNA
© Pollard et al 2004
Received: 05 November 2003
Accepted: 21 January 2004
Published: 21 January 2004
Skip to main content
© Pollard et al 2004
Received: 05 November 2003
Accepted: 21 January 2004
Published: 21 January 2004
Numerous tools have been developed to align genomic sequences. However, their relative performance in specific applications remains poorly characterized. Alignments of protein-coding sequences typically have been benchmarked against "correct" alignments inferred from structural data. For noncoding sequences, where such independent validation is lacking, simulation provides an effective means to generate "correct" alignments with which to benchmark alignment tools.
Using rates of noncoding sequence evolution estimated from the genusDrosophila, we simulated alignments over a range of divergence times under varying models incorporating point substitution, insertion/deletion events, and short blocks of constrained sequences such as those found incis-regulatory regions. We then compared "correct" alignments generated by a modified version of the ROSE simulation platform to alignments of the simulated derived sequences produced by eight pairwise alignment tools (Avid, BlastZ, Chaos, ClustalW, DiAlign, Lagan, Needle, and WABA) to determine the off-the-shelf performance of each tool. As expected, the ability to align noncoding sequences accurately decreases with increasing divergence for all tools, and declines faster in the presence of insertion/deletion evolution. Global alignment tools (Avid, ClustalW, Lagan, and Needle) typically have higher sensitivity over entire noncoding sequences as well as in constrained sequences. Local tools (BlastZ, Chaos, and WABA) have lower overall sensitivity as a consequence of incomplete coverage, but have high specificity to detect constrained sequences as well as high sensitivity within the subset of sequences they align. Tools such as DiAlign, which generate both local and global outputs, produce alignments of constrained sequences with both high sensitivity and specificity for divergence distances in the range of 1.25–3.0 substitutions per site.
For species with genomic properties similar toDrosophila, we conclude that a single pair of optimally diverged species analyzed with a high performance alignment tool can yield accurate and specific alignments of functionally constrained noncoding sequences. Further algorithm development, optimization of alignment parameters, and benchmarking studies will be necessary to extract the maximal biological information from alignments of functional noncoding DNA.
The increasing availability of genome sequences of related organisms offers myriad opportunities to address questions in gene function, genome organization and evolution, but also presents new challenges for sequence analysis. Many classical tools for sequence analysis are obsolete, and there has been active effort in recent years to develop tools that work efficiently with whole genome data. Aligning long genomic sequences – the first step in many analyses – is substantially more complex and computational taxing than aligning short sequences, and many methods have been developed in recent years to address this challenge (reviewed in[1, 2]). Nevertheless, comparative genomic researchers are still faced with the task of making decisions such as which alignment tools to use and which genomes to compare for their particular application. Benchmarking studies that address both the selection of alignment methods and the choice of species can provide the needed framework for informed application of genomic alignment tools and biological discovery in the field of comparative genomics.
Research in alignment benchmarking has focused on the alignment of protein-coding sequences[3, 4], where independent evidence (either the three-dimensional structure of a protein sequence[5, 6] or cDNA sequence[7, 8]) is available to use as a "gold standard" to assess the relative performance of different alignment tools. In contrast, little is known about the relative performance of tools to align noncoding sequences, which comprise the vast majority of metazoan genomes and contain many functional sequences includingcis-regulatory elements that control gene regulation. For noncoding sequences, little external evidence is available to evaluate alignment tool performance. Benchmarking, however, can be achieved by simulating sequence divergencein silicowhere it is possible to generate sequences that are related by a known, "correct" alignment. Simulation experiments have been used extensively to assess the performance of different methods for phylogenetic reconstruction. Yet only a few studies to date have exploited simulated data to benchmark alignment tools[11–18], and currently none have done so explicitly for the purposes of functional noncoding sequence alignment.
Here we present results of a simulation-based benchmarking study designed to assess the performance of eight tools (Avid, BlastZ, Chaos, ClustalW, DiAlign, Lagan, Needle, and WABA) for the pairwise alignment of noncoding sequences. We have chosen to address the question of pairwise alignment since pairwise alignment methods often are used in the construction of multiple alignments, since the evaluation of pairwise alignment performance is more tractable than that of multiple alignment, and since pairwise alignment performance is an important part of a general assessment of noncoding alignment strategies. We have chosen to model noncoding sequence evolution in the genusDrosophilaas a biological system for methodological evaluation, because of the high quality sequence and annotations available forD. melanogaster[19, 20], and the recent availability of the genome sequence for the related species,D. pseudoobscura. In addition, because of the high rate of deletion as well as the relatively low density of repetitive DNA as compared with mammalian genomes[22–24],Drosophilanoncoding regions are likely to be enriched for sequences under functional constraint. Previous results indicate thatDrosophilanoncoding regions contain an abundance of short blocks of highly conserved sequences, but that the detection of these sequences is dependent on the alignment method used. Optimizing strategies for the accurate identification of functionally constrained noncoding sequences will play a critical role in the annotation ofcis-regulatory elements and other important noncoding sequences inDrosophilaas well as other metazoan genomes.
In this study, we use empirically-derived estimates to parameterize simulations of noncoding sequence evolution over a range of divergences that includes those between species commonly used in comparative genomics such asH. sapiens-M. musculus[26, 27],C. elegans-C. briggsae[28, 29] andD. melanogaster-D. pseudoobscura[30, 31]. Alignments of simulated descendent sequences produced by the tools under consideration were compared to correct alignments and various performance measures were calculated. In general, we find that global tools (Avid, ClustalW, DiAlign-G, Lagan, and Needle), which align the entirety of input sequences, tend to have the highest accuracy over entire sequences as well as within interspersed blocks of constrained sequences, but both measures were decreasing functions of divergence. Local tools (BlastZ, Chaos, DiAlign-L, and WABA), which align subsets of input sequences, tend to have the highest accuracy for the portion of the sequences they align, but the proportion of sequences included in their alignments decreased quickly with increasing divergence distance. For intermediate to high divergences, local tools also showed a high specificity for only aligning interspersed blocks of constrained sequences. Despite these general trends, we find that some tools can systematically out-perform others over a wide range of divergence distances. These results should prove useful for comparative genomics researchers and algorithm developers alike.
To link our simulations to species commonly used in comparative genomic analyses of noncoding DNA, we estimated silent site divergence (Ks) betweenH. sapiensvs.M. musculus,C. elegansvs.C. briggsae, andD. melanogastervs.D. pseudoobscura(see methods). Since estimates of Ksare highly dependent on methodology, we sought to generate estimates between these three species pairs using a single method. We estimate the mean (and median) of Ksmeasured in expected number of substitutions per silent site, for these species pairs to be:H. sapiensvs.M. musculus0.64 (0.56);C. elegansvs.C. briggsae, 1.39 (1.26); andD. melanogastervs.D. pseudoobscura, 2.40 (2.24). We note that these divergence estimates do not underlie our simulation, but rather are intended to frame the interpretation of our simulation results in a biological context.
Summary of parameters used in simulations of noncoding sequence evolution.
this work (Fig.1)
AT : GC
60 : 40
Transition / Transversion Bias
Point substitutions : Indels
10 : 1
Median constrained block length
D. melvs.D. vir
Mean density of constrained blocks
D. melvs.D. vir
Overall identity between derived pairs in the regimes without constrained blocks decreases to the random background of 0.26 (the sum of the squares of the mononucleotide frequencies) by 5.0 substitutions per site with and without indel evolution (Figure2A and2B, red crosses). In the regimes with constrained blocks, unconstrained sites have the same level of identity as entire sequences in the regimes without constrained blocks (Figure2C and2D, green diamonds), whereas the identity in the constrained blocks is much greater (Figure2C and2D, yellow x's). In the regimes with indel evolution, the fraction of the ancestral sequence remaining diminishes most quickly in the absence of constrained blocks (Figure2B, green triangles). In regime C (with constrained blocks and without indel evolution), the fraction of constrained sites in derived sequences matches the input parameter of 0.2 (Figure2C, blue checked-boxes). However, in regime D (with constrained blocks and indel evolution), the fraction of constrained sites in derived sequences decreases below the input parameter of 0.2 at large divergence distances (Figure2D, blue checked-boxes). This is because the derived sequences are on average longer than ancestral sequences in regime D, differing by 300–400 bp at 1 substitution per site, 400–500 bp at 2 substitutions per site and 700–800 bp at 5 substitutions per site. In our simulation there are equal input rates of insertion and deletion, however deletions are unable to extend into constrained blocks and are omitted, creating a net excess of insertions to deletions. This phenomenon was recently proposed as a possible explanation for differences in observed insertion:deletion ratios in unconstrained dead-on-arrival retrotransposon pseudogenes versus noncoding sequences flanking genes.
Unaligned pairs of derived sequences generated by ROSE were used as input to each of the eight genomic alignment tools (see Methods) and resulting alignments were compared to the simulated alignments produced by ROSE. Our objective was to test the off-the-shelf performance of these tools over a wide range of different divergences, so each tool was run using default parameter settings. In addition, BlastZ and Chaos were run using author suggested settings (BlastZ-A and Chaos-A), as described in the Methods. We note that the output of DiAlign can be treated as both a global alignment as well as a local alignment, so we analyzed both (DiAlign-G and DiAlign-L). Alignments produced by each tool were scored for the overall coverage and overall sensitivity for all regimes (A-D), and were also scored for constraint coverage, constraint sensitivity, constraint specificity, and local constraint sensitivity in the regimes with constrained blocks (C and D) (see Methods for details).
Constraint sensitivity of all tools is much better than overall sensitivity but, as with constraint coverage, the degree of improvement varies considerably across tools (Figure5C and5D). Similar to overall sensitivity, global tools tend to maintain the highest sensitivity out to large divergence distances in the presence of constrained sites. It is of note that in the presence of indel evolution (Figure5D), constraint sensitivity of the best performing global tools (as well as the local Dialign-L) closely parallels the decrease in identity of constrained sites (Figure2D), suggesting that they are attaining near-maximal constraint sensitivity. Most tools show only moderate decreases in constraint sensitivity in the presence of indel evolution but a few, like ClustalW, Chaos-A, and BlastZ have dramatic decreases in constraint sensitivity in the presence of indel evolution.
Local constraint sensitivity (Figure6) is equivalent to constraint sensitivity (Figure5) for the global tools, but for the local tools it differs in that it is a measure of their constraint sensitivity just within the subsequences they align. For BlastZ, BlastZ-A, Chaos, and DiAlign-L, local constraint sensitivity is nearly maximal (1.0) with and without indel evolution across all divergences studied. For Chaos-A and WABA, local constraint sensitivity varies with divergence distance and is less than the other local tools. Thus local tools can produce nearly perfect alignments within constraint blocks while maintaining relatively high constraint specificity, though it is important to note that this may not be meaningful if the coverage of a tool is extremely low (e.g. BlastZ, BlastZ-A, Chaos).
In this report we investigate the performance of eight pairwise genomic alignment tools to align functional noncoding DNA such as that found in metazoancis-regulatory regions. To do so, we have used a biologically-informed simulation approach to determine off-the-shelf performance over a range of divergence distances. This study provides important information regarding the ability of genomic alignment tools to identify and align constrained sequences in noncoding regions, which would not otherwise be possible. We argue that a simulation study is necessary to achieve our goal since large datasets of functionally annotated noncoding sequences are not available to use as "gold standards" of alignment accuracy. Likewise, datasets of large orthologous genomic regions spanning a range of divergence distances are only recently becoming available[31, 34]. As is common in alignment benchmarking[4, 17, 35], we have studied performance of alignment tools using default parameters since fundamental differences in objective functions, scoring matrices, the type and values of parameters, and algorithmic design prevent a systematic exploration of parameter space.
We have attempted to construct a realistic simulation of noncoding sequence evolution and test alignment performance for species with genomic properties similar toDrosophila. Noncoding alignment assessment for mammalian and other species with large, repeat-rich genomes would require modifications to our current simulation, such as the inclusion of ancestral repeats and lineage-specific transposition events. Moreover, as more becomes known about the substitution process in noncoding regions (especially those under weak primary sequence constraint), it will be important to implement more realistic models such as context-dependent substitution[36–38]. It would be also instructive to assess alignment performance based on a simulation that decouples suppression of indel rates from substitution rates, given the possibility that the spacing (but not the primary sequence) between conserved noncoding segments may be constrained. In addition, though we have attempted to be systematic in our evaluation of tools, we unfortunately cannot have included all available pairwise alignment tools. As new pairwise alignment tools emerge and old tools are modified or brought to our attention, we will update our results periodically on the web using the same set of simulated alignments presented here. Moreover, assessment of tools which take advantage of the phylogenetic information and higher signal-to-noise inherent in multiple alignments will be an essential extension to this work to provide a more general evaluation of strategies for noncoding alignment.
From the standpoint of the most biologically relevant simulation regime studied here (D, which includes indel evolution and interspersed blocks of constrained sequences), our results indicate that global alignment tools have the highest sensitivity in general to align orthologous sites accurately in noncoding sequences, as well as blocks of constrained sites (Figures4D,5D). We find that constraint sensitivity of the top global tools can be quite high (>75%) and limited only by sequence identity in constrained sites at intermediate divergence distances (1.25–3.0 substitutions per site), whereas overall sensitivity is relatively low beyond such intermediate divergence distances. The improved performance of global tools over local tools is largely a consequence of incomplete coverage of both constrained and unconstrained sites in alignments produced by local tools (Figure3). The subset of sequences aligned by the highest performing local tools, however, is accurately aligned and specifically corresponds to constrained sites (Figure6). In fact, most local tools can effectively discriminate between constrained and unconstrained sites to greater than 80% specificity at intermediate divergence distances while the constrained portions of their alignments are nearly perfectly aligned at large divergence distances. Finally, when compared with regime C (which excludes indel evolution but includes interspersed constrained blocks), it is clear that our model of indel evolution affects alignment coverage, sensitivity and specificity, but not enough to overturn these major trends.
These results have important implications for the analysis of functional noncoding sequences. First, if a researcher's goal is to alignall constrained sites in a noncoding region, then a global tool like Lagan will reliably produce the best results, but will require post-processing to identify constrained sequences[40, 41]. Conversely, if one's goal is to alignonly constrained blocks in a noncoding region, then a local tool like Chaos will reliably produce the best results, provided that complete recovery of all constrained sequences is not required. The distinct virtues of both global and local tools are currently incorporated in the output of only one alignment tool, DiAlign. For this reason, use of the global parse of DiAlign (DiAlign-G) can provide high coverage and sensitivity across entire noncoding regions, while use of the local parse of DiAlign (DiAlign-L) will specifically provide highly accurate alignments of blocks of constrained sites. In light of these results, we recommend the further development of global alignment tools that also output a local parse of high confidence local alignments contained within, which should be possible since local anchors are often used in the construction of the global alignment (e.g.[7, 8]).
Our results also indicate that for species with structural and evolutionary constraints on noncoding sequences such as those found inDrosophila, DiAlign can produce alignments with high coverage and sensitivity, as well as high specificity to detect constrained sites in the range of 1.25–3.0 substitutions per site. Since the divergence betweenD. melanogastervs.D. pseudoobscuraand betweenC. elegansvs.C. briggsaefalls within this range, we suggest that the use of DiAlign for detecting functionally constrained noncoding sequences will prove successful in these taxa on a genomic scale. In contrast, our results also indicate that species pairs such asH. sapiensand.M. musculusmay not be sufficiently diverged for a single pairwise comparison to provide the needed resolution to detect functionally constrained noncoding sequences, though differences in genome organization and evolution between flies and mammals require a more thorough evaluation of this claim. This conclusion, however, supports results based on Poisson modelling of point substitution that approximately 3 substitutions per site would be needed to detect functional constrained sites reliably in mammalian noncoding DNA.
Finally, the results presented here also imply that biological and technical conditions exist with which to study with accuracy the evolutionary events underlying the process ofcis-regulatory evolution in flies and worms. Current evolutionary models ofcis-regulatory sequence divergence posit the gain and loss of transcription factor binding sites, even under constant functional constraints[43, 44]. However, the absence of alignable binding sites in comparisons of divergent sequences may result from inaccuracies in alignment as well as thebona fideloss of transcription factor binding sites. We suggest that alignments of noncoding sequences using tools such as DiAlign in the range of 1.25–3.0 substitutions per site are of sufficient accuracy to measure binding site loss among divergent species pairs, such as the high levels recently reported in the genusDrosophila[45, 46].
Our study demonstrates that recently developed alignment tools have the potential to produce biologically meaningful alignments of functional noncoding DNA on a genome scale. Continued development of alignment algorithms in conjunction with parameter optimization and continued benchmarking will be necessary to provide the highest quality genomic alignments under the wide diversity of genomic and evolutionary scenarios to be studied.
To generate biologically relevant input sequences for our simulation, we estimated properties of noncoding sequences in the genome sequences of the fruitfly,D. melanogaster. First we extracted all noncoding regions from the Release 3D. melanogastergenomic sequences based on annotations in the Gadfly database[19–20 47]. This was accomplished by masking all DNA corresponding to coding exons, producing inter-coding-exon intervals. Subsequent to extracting noncoding regions, transposable elements were masked using annotations in Gadfly to create "pre-integration" noncoding sequences. In our analysis, we chose to treat all noncoding sequences (intergenic, intronic, untranslated region) together since many noncoding sequences cannot be unambiguously categorized because of alternative splicing or alternative promoter usage. Moreover, previous results revealed that similar evolutionary constraints act on intergenic and intronic sequences inDrosophila. Summary statistics of noncoding sequence lengths were calculated using the R statistical package (Figure1).
The probabilistic dependence of adjacent bases inD. melanogasternoncoding sequences was assessed by Markov chain analysis in order to create an accurate model of random noncoding sequences. TE-masked noncoding sequences were concatenated, and n-mers of size 1 to 10 were counted. Counts of reverse complementing n-mers were averaged, and used to estimate frequencies of each n-mer. Based on these counts and frequencies, we determined the likelihood of Markov chains of orders 1 through 9 describingDrosophilanoncoding sequences, and evaluated the likelihood of each Markov chain using the Bayesian information criterion[49, 51]. This analysis revealed thatD. melanogasternoncoding sequences are best modeled by a 7th-order Markov chain (data not shown). We therefore created the ancestral input sequences for our evolution simulations using a 7th-order Markov chain. We note that because our evolutionary simulation models bases independently (see below), the higher order structure of these ancestral input sequences was not maintained in the more divergent derived output sequences. Nevertheless, sequences generated by a 0th-order Markov chain gave qualitatively and quantitatively similar simulation and alignment results, with correlation among performance measures for the 0th-order and 7th-order generated sequences exceeding an r2of 0.97 (data not shown).
Estimates of silent site divergence (Ks) betweenH. sapiensvs.M. musculus,C. elegansvs.C. briggsae, andD. melanogastervs.D. pseudoobscurawere obtained using the yn00 method in PAML (version 3.13)[52, 53]. The mean and median of Kswere calculated for 29 fly, 193 worm, and 153 mammalian coding sequence alignments taken from references[31, 28] and, respectively.
Noncoding sequence evolution was simulated using a modified version of the sequence simulation program ROSE. In general, in the absence of large datasets of noncoding sequences from closely relatedDrosophilaspecies, we have taken estimates of noncoding evolution from previous results reported in the literature. Beginning with ancestral sequences, evolution occurred on two descendent branches of equal length under the HKY model of point substitution, with a transition/transversion bias of 2 to reflect the nucleotide and transition biases observed inDrosophilanoncoding sequences[25, 55–56]. The substitution rate was set to 0.01 such that a branch length unit was on average 0.01 substitutions per site. Total branch lengths spanned a range of divergence times from 0.25 to 5.0 substitutions per site. Insertion/deletion evolution was based on the length distribution of polymorphic indels estimated in, and occurred at a 10-fold lower rate than point substitution, approximating relative rates estimated in[22, 23].
To model the evolution of constrained blocks in noncoding sequences a modification of the ROSE sequence simulation program was developed to map constraints on ancestral sequences onto derived sequences (available for download as ROSE version 1.3 from). Constraints on noncoding sequences were modelled as short blocks of highly conserved sequences typical ofcis-regulatory sequences, and follow a lognormal distribution with parameters estimated in. On average, interspersed blocks of constrained sites accounted for 20% of the sites in ancestral sequences, a conservative estimate of constraint inDrosophilanoncoding DNA. Parameters used in our simulations are summarized in Table1.
Estimation of evolutionary distance for simulated alignments was performed using the F84 model of sequence evolution in the DnaDist program of the PHYLIP package with a transition:transversion ratio of 1.0 (note that a transition:transversion ratio of 1.0 in PHYLIP is equivalent to a transition/transversion bias of 2 in ROSE, see discussion in). Summary statistics for the simulations were calculated using the R statistical package (Figure2).
The alignment tools tested in this study were chosen based on the criteria that they are (1) publicly available, (2) run in batch mode from the command line and are able to produce (3) strictly co-linear, (4) error-free, pairwise genomic alignments of sequences (5) up to 10 Kb in length. Tools like BBA (5), Bl2seq (3), DBA (4), MUMmer (3), Owen (2) and SSEARCH (3) were not evaluated since they do not satisfy one of these criteria. We now briefly describe the tools that we tested.
Avid is a pairwise global alignment tool whose general strategy for aligning two sequences is to anchor and align iteratively. A set of maximal (but not necessarily unique) matches between the sequences is constructed using a suffix tree. Dynamic programming is used to order and orient the longest matches, which are then fixed. For each subsequence remaining between the fixed matches, the process is repeated until every base is aligned. When sequences are short and the matches make up less than half of the total sequence, the program defaults to the Needleman-Wunsch algorithm.
The Chaos/Lagan suite of tools consists of a pairwise local alignment tool, Chaos, and a global alignment tool, Lagan. Chaos starts by finding all words between the two sequences of a specified length and a specified maximum number of mismatches. These words are then chained together if they are close together in both sequences. These maximal chains are then scored and all chains that are above a specified threshold are returned. Lagan starts by running Chaos with conservative parameter settings and then finds the optimal path through the maximal chains using dynamic programming. Lagan then recursively calls Chaos with increasingly more permissive parameters on the regions between each maximal chain in the optimal path. When the recursion has created a dense map of maximal chains that have been ordered with dynamic programming, Lagan runs the Needleman-Wunsch algorithm on the whole length of both sequences but puts close bounds around the maximal chains to provide the final global alignment. Chaos was run on default parameters as well as using parameters suggested by the authors: word length = 7, number of degeneracies = 1, score cut-off = 20 and extension mode on.
BlastZ is a pairwise local alignment tool that is based on the gapped BLAST algorithm that has been redesigned for the alignment of long genomic sequences. BlastZ first removes lineage-specific interspersed repeats from each sequence, then searches for short near-perfect matches between the two sequences. Each match is extended first using gap-free dynamic programming and if it scores above a specified threshold it will be extended using dynamic programming with gaps; extended matches that score above a specified threshold are then kept. Part of the unique implementation of BlastZ is that it can be forced to return alignments that are both unique within each sequence as well as collinear with respect to each other. To satisfy our strict collinear requirement, we ran BlastZ with both of these options. Blastz was also run using the author's suggestion of lowering the score cut-off (k) to 2000 (BlastZ-A).
DiAlign (v. 2.1) is a segment-to-segment alignment algorithm. Like the BLAST algorithms, DiAlign looks for short ungapped segments that have a similarity that deviates from what would be expected by random chance, keeping segments with a score above a certain threshold. These high scoring segments are then aligned into a collinear global alignment using a dynamic programming algorithm. DiAlign produces a global alignment but distinguishes high confidence columns of an alignment from low confidence columns. We used DiAlign as both a global (DiAlign-G) and a local (DiAlign-L) alignment tool.
ClustalW (v. 1.8) was used on default settings. ClustalW is a progressive multiple alignment tool that reduces to the Needleman-Wunsch algorithm in the pair-wise case with default parameters of a match score of 1.9, mismatch penalty of 0, a gap open penalty of 10 and a gap extension penalty of 0.1.
The second implementation of the Needleman-Wunsch algorithm used in this study is the needle program in the EMBOSS suite of tools. needle was used with default parameter settings of a match score of 5, a mismatch penalty of 4, a gap open penalty of 10 and a gap extension penalty of 0.5.
The final tool tested, WABA, is a three-tier alignment algorithm. The first tier partitions the first sequence into overlapping windows of 2 Kb and then defines a synteny map of high scoring 2 Kb windows of the first sequence onto the second sequence. The second tier then carefully aligns syntenic regions using a seven-state, pair Hidden Markov Model that includes separate query and database insertion/deletion states, high and low noncoding conservation states, as well as three coding states (one for each position in a codon). The final tier then attempts to assemble individual alignments together into a more global alignment.
The performance of alignment tools was assessed using six basic measures: overall coverage, overall sensitivity, constraint coverage, constraint sensitivity, constraint specificity and local constraint sensitivity. Overall coverage and overall sensitivity were measured for all four evolutionary regimes (A-D) while the constraint measures were only measured in the two regimes that included constrained blocks (C, D). Alignments produced by each alignment tool were parsed to generate the statistics, which were then used to calculate each performance measure.
Each site in an alignment produced by a tool (a site being a base in one strand of a column of an alignment) can have two simulated alignment states, two constraint states, three tool alignment states, and two conditional tool alignment states. The two simulated alignment states are "homolog" (h), ungapped sites in the simulated alignments, and "no homolog" (nh), gapped sites in the simulated alignments. Simulations without indel evolution have only homolog sites since there are no gaps in the simulated alignments. The two constraint states are "constrained" (c), sites in constraint blocks, and "unconstrained" (u), sites not in constrained blocks. The three tool alignment states are "aligned" (a), sites aligned in the tool alignment, "gapped" (g), sites gapped in the tool alignment, and "not aligned" (na), sites not included in a local tool alignment. The two conditional tool alignment states are "aligned correctly" (ac), sites aligned to the same site in both the tool and simulated alignments, and "aligned incorrectly" (ai), sites aligned to different sites in the tool and simulated alignments. There are fourteen possible combinations of these states (e.g. homolog constrained aligned correctly, h_c_ac), giving us fourteen statistics to calculate for each estimated alignment. Counts for each statistic were used to calculate the following measures:
Overall coverage is the fraction of ungapped sites in a simulated alignment that are included in a tool alignment. Overall Coverage = (h_c_ac + h_c_ai + h_c_g + h_u_ac + h_u_ai + h_u_g) / (h_c_ac + h_c_ai + h_c_g + h_c_na + h_u_ac + h_u_ai + h_u_g + h_u_na)
Overall sensitivity is the fraction of ungapped sites in a simulated alignment that are aligned to the correct base in a tool alignment. Overall Sensitivity = (h_c_ac + h_u_ac) / (h_c_ac + h_c_ai + h_c_g + h_c_na + h_u_ac + h_u_ai + h_u_g + h_u_na)
Constraint coverage is the fraction of ungapped constrained sites in a simulated alignment that are included in a tool alignment. Constraint Coverage = (h_c_ac + h_c_ai + h_c_g) / (h_c_ac + h_c_ai + h_c_g + h_c_na)
Constraint sensitivity is the fraction of ungapped constrained sites in a simulated alignment that are aligned to the correct base in a tool alignment. Constraint Sensitivity = (h_c_ac) / (h_c_ac + h_c_ai + h_c_g + h_c_na)
Constraint specificity is the fraction of unconstrained sites in a simulated alignment that are gapped or not included in a tool alignment. Constraint Specificity = (h_u_g + h_u_na + nh_u_g + nh_u_na) / (h_u_ac + h_u_ai + h_u_g + h_u_na + nh_u_a + nh_u_g + nh_u_na)
Local constraint sensitivity is the fraction of sites that are both, contained in a tool alignment and are ungapped constrained sites in a simulated alignment, that are aligned to the correct base in the tool alignment. Local Constraint Sensitivity = (h_c_ac) / (h_c_ac + h_c_ai + h_c_g)
For each of these six measures, a mean and standard error of the mean were calculated for up to 1000 replicates (local tools do not always return an alignment and replicates which produced no alignment were not counted toward the mean) using R.
We thank Mark Yandell and Chris Mungall for providing scripts to extract noncoding sequences, Erwin Frise for assistance with cluster computing, and Cristian Castillo-Davis and Anton Nekrutenko for kindly providing worm and mammalian coding sequence alignments. We thank Venky Nandagopal for discussions on alignment specificity, Hunter Fraser, Emily Hare, Alan Moses, and Monty Slatkin for critical reading of the manuscript, and Alex Kondrashov and one anonymous reviewer for helpful comments on the manuscript. This work was supported by NIH grant HG00750 to G. Rubin. DAP is supported by NIH training grant T32 HG00047 to D. Rokhsar and J. Rine. CMB is supported by NIH training grant T32 HL07279 to E. Rubin and by a Royal Society USA Research Fellowship. MBE is a Pew Scholar in the Biomedical Sciences. Research was conducted at the Lawrence Berkeley National Laboratory under Department of Energy contract DE-AC0376SF00098.
This article is published under license to BioMed Central Ltd. This is an Open Access article: verbatim copying and redistribution of this article are permitted in all media for any purpose, provided this notice is preserved along with the article's original URL.