s2fieldsofstudy:
Computer Science
Biology

text:
Quantifying the Displacement of Mismatches in Multiple Sequence Alignment Benchmarks

Multiple Sequence Alignment (MSA) methods are typically benchmarked on sets of reference alignments. The quality of the alignment can then be represented by the sum-of-pairs (SP) or column (CS) scores, which measure the agreement between a reference and corresponding query alignment. Both the SP and CS scores treat mismatches between a query and reference alignment as equally bad, and do not take the separation into account between two amino acids in the query alignment, that should have been matched according to the reference alignment. This is significant since the magnitude of alignment shifts is often of relevance in biological analyses, including homology modeling and MSA refinement/manual alignment editing. In this study we develop a new alignment benchmark scoring scheme, SPdist, that takes the degree of discordance of mismatches into account by measuring the sequence distance between mismatched residue pairs in the query alignment. Using this new score along with the standard SP score, we investigate the discriminatory behavior of the new score by assessing how well six different MSA methods perform with respect to BAliBASE reference alignments. The SP score and the SPdist score yield very similar outcomes when the reference and query alignments are close. However, for more divergent reference alignments the SPdist score is able to distinguish between methods that keep alignments approximately close to the reference and those exhibiting larger shifts. We observed that by using SPdist together with SP scoring we were able to better delineate the alignment quality difference between alternative MSA methods. With a case study we exemplify why it is important, from a biological perspective, to consider the separation of mismatches. The SPdist scoring scheme has been implemented in the VerAlign web server (http://www.ibi.vu.nl/programs/veralignwww/). The code for calculating SPdist score is also available upon request.

Introduction
Since its inception in the 1970s-early 1980s [1,2], Multiple Sequence Alignment (MSA) has been one of the most prominent computational techniques in modern molecular biology, and its significance and popularity has only increased since the advent of next-generation sequencing (NGS) techniques. MSA plays a crucial part in various types of biological analyses, e.g. inference of phylogenetic relationships between organisms, conserved sequence elements, functional regions, correlated mutations, and is a vital prerequisite for protein homology modeling [3][4][5][6][7][8][9]. MSA methods are used to align multiple homologous sequences, i.e. sequences that are assumed to share a common evolutionary origin; in the resulting alignment the matched residues should indicate that they derive from a single residue in an ancestral sequence, implying that they are structurally superposable and/or involved in a common function [10].
There are two main approaches which are used to assess the quality of a MSA based on sequence information. One is a reference-independent alignment score, that can be calculated solely from sequence information of a single MSA, which is very useful in case curated or structural information on proteins of interest is unavailable. An example of such a score is the sumof-pairs scoring scheme (which we will refer to as standalone_SP score) [11]. This score is calculated for a stand-alone MSA by exploiting an evolutionary model from which probabilities of pairwise residue conservations and mutations are derived [12,13]. Such a score is typically used as the objective score for MSA methods, but may also be used to judge the quality of a resulting MSA. The standalone_SP score, calculated using general residue substitution matrices, is rarely used to benchmark MSA methods as the more evolutionary divergent sequence sets will attain lower standalone_SP scores than less divergent sets, making comparison and normalization difficult. Furthermore, an alignment with a high standalone_SP score does not necessarily correspond to a biologically correct alignment as this score does not take structural similarity and conservation of functional residues into account. Therefore its use as benchmark score is not insightful.
The most widely used approach to benchmark alignment quality is referred to as referencedependent benchmarking, which aims to compute the similarity between the MSA in question with a gold standard reference alignment of the same sequence set. Here, there is no a-priori tendency for more divergent alignments to score lower, although attaining high query-reference alignment similarity is known to be more difficult for evolutionary divergent cases [14].
To address these complications, more recent techniques for assessing alignment quality abounded, which make use of additional information such as protein tertiary structure information. However, since these methods do not rely on sets of reference alignments, they are referred to as reference-independent benchmark techniques [15]. Armougom et al. [16] exploited the structural distance between pairs of residues to score pairwise sequence alignment quality without the need for structural superpositioning of the associated tertiary structures. This type of alignment assessment methods are especially useful in pairwise target-template alignment quality assessment, for example in homology modelling and threading algorithms [17,18].
The main goal of MSA benchmarking is to measure the ability of different MSA methods to reproduce trusted reference alignments [19]. As already mentioned, these reference alignments are normally obtained from a manually (or semi-automatically) curated alignment database, such as BAliBASE [20], OXBENCH [21], PREFAB [22], or SABmark [23]. It has been shown that in addition to the divergence of the sequences, the performance of an MSA method depends also on several other factors such as sequence lengths, the presence or absence of large insertions and N/C terminal extensions, or over-representation of one or more protein subfamilies [14]. Therefore, a thorough comparative analysis of MSA methods requires accurate reference alignments representing various alignment problems to objectively assess the quality of MSA methods. The ultimate benefit of MSA benchmarking to the user is knowing what is the best method for a specific alignment problem. Therefore, the availability of scoring schemes that bring out different strengths and weaknesses of an alignment method is important.
The similarity between a reference alignment and the corresponding query alignment is commonly calculated using the column score (CS) and the sum-of-pairs (SP) score. The CS score calculates the query-reference alignment similarity as the fraction of columns in the reference alignment that is exactly reproduced in the query alignment, meaning that a single mistake in a MSA column (i.e. a mismatch in a single sequence) will disqualify the entire column ( Fig 1A). The SP score determines the similarity by the fraction of aligned residue pairs in the reference alignment that is correctly reproduced in the query alignment (Fig 1A), i.e. the proportion of mistakes over all residue pairs in a column is reflected in the score. As a result, the SP score is softer and less strict than the CS score. Note that this type of SP score is not the same as the standalone_SP score as it takes the reference alignment and not the amino acid substitution matrix into account. Fig 1 illustrates the difference between CS and SP scores and their limitations. In the toy example ( Fig 1B) one can observe that the CS score plummets to 0 due to minute shifts in the sequences. This indicates that the CS score is highly influenced by the number of the sequences and the gap locations compared to the SP score. The toy example also shows that query alignment A is more similar to the reference alignment than query alignment B, hence query alignment A should have a better score than query alignment B. However, this is not the case if the similarity between the query and the reference alignments are to be measured with the CS and SP score. Since both scoring schemes only take matches and mismatches into account, they are not able to discriminate whether one mismatch, i.e. a misplaced residue, is more shifted than the other.
Due to their nature, SP and CS scores may not be discriminative for MSAs that contain long stretches of gaps. For example, one MSA method may detect a long insertion/deletion appropriately, but get the boundaries slightly wrong, resulting in a small shift, while another method may miss the insertion/deletion altogether. The former case is arguably less incorrect, but may yield similar SP and CS scores as the latter (see the case study in the "Results and Discussion" section for a real-life example of a similar situation in MSA benchmark). The key problem is that the SP and CS scores do not take the severity of displacement of mismatches into account. In practice the length of misalignments often plays a major role in biological analyses where the accuracy of an MSA is of paramount importance, for example in analyses where manual alignment editing is required to achieve a biologically sound alignment. In such cases, misalignment of (stretches of) identical or similar residue pairs with relatively small shifts will be much easier detected by the human eye than misplaced residues positioned further away from their correct position in the alignment.
The purpose of this study is to investigate whether one can obtain a more fine-grained indication of how much a query alignment actually differs from the reference by taking distances between mismatched residue pairs into account in MSA comparison, which is particularly meaningful in the case of alignments with long stretches of insertions/deletions. In order to do this we developed a variant of the SP score, SPdist, which takes into account the distance of mismatched residue pairs instead of the binary classification into match and mismatch state as in the case of the standard SP score. Cline et al. (2001) devised a distance-based scoring scheme for pairwise alignments called shift score which also aims to differentiate the treatment between : Toy alignments to illustrate the difference punishment for alignment shifts in CS, SP, and SPdist (Eq 8) scores. All sequences in this contain a small sequence motif H-x-D-E which has to be aligned correctly. Both query A and B contain the same number of mismatches (indicated by the same SP score for both queries), the sole difference between these two queries is that sequence 2 in query B is shifted five positions further. In this situation one might argue that although both alignment are wrong, i.e. the sequence motifs are not aligned properly, query A is a slightly better alignment than B as the motifs in query A is positioned closer to each other than in query B. major and minor shifts in alignment comparison [24]. In this study we will also compare the behaviour of our SPdist score with the shift score.

MSA benchmark data set
BAliBASE is one of the most widely used alignment database to assess the quality of MSA methods. It contains various MSAs obtained from manually curated 3D structural superpositions [20]. The alignments in BAliBASE are organized in reference sets, each representing a specific alignment problem: Reference 1 (RV11 and RV12) contains alignments of equidistant sequences; Reference 2 (RV20) holds alignments of protein families with with one or more socalled orphan sequences(< 25% identical); Reference 3 (RV30) comprises alignments of divergent protein subfamilies with < 25% sequence similarity between subfamilies; Reference 4 (RV40) has alignments of sequences with large N/C terminal extensions; Reference 5 (RV50) alignments with large internal insertions; while Reference 6-8 contain highly specialized alignment problems such as, transmembrane proteins, proteins with multiple repeats, and inverted domains. Each reference set in BAliBASE contains two types of alignments: a full-length alignment and a truncated one containing only highly conserved regions. We selected the full-length alignments for this study as they are more akin to real alignment situations. Out of the above eight reference sets in BAliBASE, we only selected sets 1 to 5 for our analysis since the remaining sets represent more specific alignment situations that are less amenable to the generic MSA methods assessed in the current study. The five selected benchmark sets contain 218 reference alignments in total. Although the aforementioned reference alignment databases PREFAB [22] and SABmark [23] are much larger than BAliBASE, both these datasets only provide pairwise alignments and no true MSA as reference for each alignment case, which is undesirable for the SPdist score introduced in this paper. Similarly, the smaller OXBENCH reference set [21] comprises multiple sequence alignments containing regions that are annotated as unreliably aligned, which also makes this database less amenable to evaluating the SPdist score.

Generation of multiple sequence alignments
Various multiple sequence alignment programs were selected for this study, namely ClustalW2 (ver. 2.0.10) [25], Clustal Omega (ver. 1.2.0) [26], MAFFT (ver. 7.215) [27], Muscle (ver. 3.7) [22], ProbCons (ver 1.12) [28], and Praline [29]. As the aim of this study is to analyse the behavior of the SP and SPdist scoring schemes, and not to select the best alignment program, we made no attempt to optimize these alignments programs and ran them using their default settings for all test sets. This holds particularly for the Praline method, where, out of a number of different alignment strategies, the suboptimal standard progressive alignment technique was used as a control. An exception was made for MAFFT, which was run using both the default and optimized parameter settings.

SP and CS scoring schemes
The column score (CS) is calculated as the percentage of correctly aligned columns for a reference and query alignment, having M sequences and N columns in the reference alignment: where Q(k i ) is a lookup function that translates the column of a residue k i of sequence i in the reference alignment into the corresponding column in the query alignment (Fig 2). The compare function checks if the positions of the residues in a column k of the query alignment correspond to a single column in the reference alignment. If all residues compared are in the same column, the result is 1, and 0 otherwise. For all columns in the reference alignment the positions of the residues in the query alignment are compared. The outcomes over all columns are then summed and divided by the total number of columns in the reference alignment. The sum-of-pairs score (SP) is calculated by taking all aligned pairs of residues in the reference and Illustration of the mapping of aligned residue pairs k i and k j in reference alignment onto the corresponding pairs in the query alignment (Q(k i ) and Q(k j )). The distance for residue pair k i and k j is calculated as the absolute difference in alignment columns between Q(k i ) and Q(k j ). doi:10.1371/journal.pone.0127431.g002 Quantifying the Displacement of Mismatches in MSA Benchmarks the corresponding pairs in the query alignment into account: Here, each aligned pair k i and k j in the reference alignment is compared to the corresponding pair Q(k i ) and Q(k j ) in the query alignment. If the pair is indeed also aligned in the query alignment a 1 is added; a 0 otherwise. The total score over all pairwise comparisons is divided by the total number of matched residue pairs 1 2 NM M À 1 ð Þin the reference alignment.

Distance-based SP score (SPdist)
In this work we develop a new score, the SPdist score, that is similar to the SP score but incorporates distance information between misaligned residue pairs rather than a score of 1 for matches and 0 for mismatches. For each aligned pair of residues in the reference alignment, now the distance between the corresponding two residues in the query alignment is determined. As the final SPdist score, the average mismatch distance over all aligned pairs in the reference alignment is calculated. The way mismatch distances are incorporated in the SPdist score is inspired by similar measures used in the field of structural alignment, such as the distance-based scoring scheme employed in the ProSup structural alignment method [30] and the Global distance test (GDT_TS) score which is mainly employed in assessing structural similarity between alternative 3D structures of a single protein, where one structure represents the template and the other the model [31]. The GDT_TS scoring scheme determines the number of template-model Cα atom pairs (G(v)) that are within the range of a predetermined distance cut-off, v. The score is then calculated as the average percentage of residues whose template-model Cα atom pairs fit within four predetermined distance cut-offs: 1, 2, 4, or 8 Å [32]; as shown in the following formula: where t is number of all Cα atoms in the template structure. In SPdist we similarly employ binning as a pre-processing method to be able to weight the effects of various alignment shift distances; for example, large shifts would otherwise tend to dominate the final score. SPdist scoring works as follows: First, we calculate the distance between all query residue pairs that are matched in the reference alignment. Then, these residue pair distances are grouped into separate bins, based on distance thresholds t, followed by SPdist(t) calculation over all the thresholds (Eqs 6 and 7) SPdistðtÞ ¼ where D t is number of residue pairs whose distance is less than t. Note that within a bin, a mismatch shift that is smaller than the threshold specified for that bin is counted positively, such that a higher SPdist score for that particular bin implies less mismatch displacement and hence higher similarity between reference and query. Since the scores are cumulative over the bins, i.e. each next bin (for larger alignment shifts) will always contain the SPdist score of the preceding bin, the bin scores grow uniformly. Both the SPdist score for the individual thresholds and the overall average score are scaled between 0 and 1, and can thus be interpreted similarly to the classical SP score. A single weighted average of SPdist score over all selected distance thresholds can be calculated as follows: Where W t is the weight of threshold t. The SPdist scoring scheme allows the user to specify different weights for the different thresholds; for example, the weights of the distance thresholds can be set based on expected frequencies gleaned from a specific collection of alignments. In this study we employed five thresholds (t = 1, 4, 10, 22, 46) with equal weights (W t = 1) to calculate the avg_SPdist score. The thresholds selected are loosely based upon the lengths of secondary and super-secondary structures [33][34][35], so that the severity of mismatch punishment is calibrated according to misalignment of these structural features. Note that SPdist(1) SP.

Shift score
The shift score [24] is somewhat akin to the SPdist score as both aim to take the displacement of mismatches into account in alignment comparison. However, these scoring schemes differ greatly in the way the score is normalized. For a pair of sequences X and Y aligned in query alignment Q and in reference alignment R, the shift score is calculated as follows: where jRj is number of aligned residue pairs in alignment R, CS R (Q i ) is the score for position i in alignment R, S(a i ) is the score for residue a i with respect to the reference alignment, and shift (a i ) is the measure of how much residue a i is shifted. if residue a i is not aligned to another residue in either query or reference alignment then shift(a i ) is undefined. The shift score ranges from − to 1.0, where is a small number (usually set by default to 0.2) used to differentiate minor shifts from major ones. If two alignments are completely identical, their shift score is 1.0. In this study we employed the shift score as it is implemented in qscore package (http://www.drive5.com/qscore/), because the original web server for calculating the shift score is unfortunately no longer available.

The behaviour of SP and SPdist scores in MSA benchmark
In order to understand the effect of the severity of mismatch shifts on the alignment score, we benchmarked six different MSA methods against BAliBASE reference alignments using the standard SP score, shift score and the new SPdist score. As our goal is solely to observe the behavior of SP, SPdist, and shift scores, all of the MSA methods that we used were ran using their default settings.
We compiled the scores produced by the three scoring schemes to analyse whether the performance rankings produced by these scoring schemes are similar to each other. We calculated the SPdist score (Eq 6) for a large range of distance thresholds (i.e. 1-100), to probe the behaviour of alignment methods as the distance threshold is varied (Fig 3). From this figure it becomes apparent that the MSA methods which obtain the best SPdist scores over small thresholds are not guaranteed to have the best SPdist scores over large thresholds. Hence it is important to consider mismatches using different distance thresholds to evaluate the performance of different MSA methods. Fig 4 shows the distribution of the benchmark scores of the different alignment methods; note that the SP score is identical to SPdist(1). The figure shows that the SP and shift scores produce very similar rankings while the SPdist scoring scheme shows fluctuations in the ranking of the MSA methods across different distance thresholds. The SPdist scores for the various mismatch distance thresholds reveal that the different MSA methods differ in their tendency to insert smaller or larger mismatch shifts.
For example, one can observe in Figs 3 and 4 that ProbCons alignments contain few small mismatches as it has the best SP/SPdist(1) score; yet at higher distance thresholds the scores drop relative to the other methods. The distribution of the number of gaps in the alignments (Fig 5) shows clearly that overall less accurate methods (ClustalW2 and Praline) tend to generate alignments with longer mismatches.
Moreover, it is striking that for quite a large number of reference alignment sets the MSAs produced by the methods never converge at large distances (Fig 4) which is apparent from the spread of individual scores at larger distance thresholds. This shows that for all methods considered here, there is room for improvement by making these more accurate over longer distances.
The MSA benchmark scores described in this study are not true metrics because they violate the principle of symmetry as the similarity between the reference and query MSAs are based on Here the reference sets 1-5 of the BAliBASE dataset are used to calculate the scores. The ranking for the different methods changes at larger distance thresholds, compared to both the SP and the shift score; note that SP score = SPdist (1). Note that the variance of individual alignment score remain high at large distance thresholds. The average values associated with the median values displayed in this figure are reproduced in S1 Table. doi:10.1371/journal.pone.0127431.g004 Quantifying the Displacement of Mismatches in MSA Benchmarks the fraction of matched residues in the reference MSA that are correctly reproduced in the query alignment; therefore d(a, b) 6 ¼ d(b, a). Consequently, the scores cannot be used to directly compare benchmark results obtained from different datasets [36].
Lastly, we wanted to analyse the relationship between SP and SPdist scoring schemes in relation to the difficulty reproducing the reference alignments. Here, we relate difficulty to the extent of gap insertion in the reference alignments, since alignments with many gaps will more easily accommodate mismatches with larger shifts. We accordingly define "gapness" as a relative measure of the number of gaps in a MSA with respect to the minimum and maximum number of possible gaps in the pairwise sequence alignments underlying the MSA (Eq 12). The rationale behind the gapness measure is to normalize the severity of misplacement according to what can be expected in each of the alignments. Gapness is defined as follows: where M is the number of sequences, Xij is the number of gaps observed in sequence i and j, max ij is the maximum number of gaps possible in the pairwise alignment between sequence i and j (max ij = length seq. i + length seq. j). While min ij is the minimum number of gaps possible in the pairwise alignment between sequence i and j (min ij = jlength seq. i-length seq. jj).  Fig 6 shows the comparison of SP and the avg_SPdist in terms of "gapness". We analysed the correlation between avg_SPdist and SP score with the "gapness" score of the reference alignments by calculating the Spearman's rank correlation. Both scoring schemes show negative correlation with respect to the "gapness" score, with the avg_SPdist score having a slightly stronger (Spearman's rho = -0.72) anti-correlation than the SP score (Spearman's rho = -0.69). Although the gapness score is not measuring amino acid relatedness, it is related to the difficulty of reproducing reference alignments as reflected in the negative correlation with the SP score. As expected, the SPdist scoring scheme is more sensitive to reference alignments containing a larger proportion of gaps than the SP scoring scheme.

Web server
VerAlign is a web tool for calculating the similarity between two MSAs according to the SP, CS, and SPdist scores. This web tool is accessible on the IBIVU website at VU University Amsterdam (http://www.ibi.vu.nl/programs/veralignwww/). VerAlign provides the user with a compact and intuitive user interface. The alignments to be compared must be in fasta or MSF format and can either be entered manually in the text boxes provided or uploaded as files (Fig 7). The user can specify the number of thresholds as well as the associated threshold values and weights via the user interface. Once a VerAlign job is submitted, the user is presented with a holding page that refreshes automatically and displays the results page when the job is completed. The output page displays SP, CS and the SPdist scores of the query alignment. The SPdist scores of the various distance thresholds are displayed alongside the avg_SPdist score. This allows the user to tune an alignment method on the behaviour of mismatch shifts over specified distance thresholds. Fig 8 shows an example of the VerAlign result page, showing the reference and query alignments and their similarity scores. The columns in the reference alignment that are regarded as highly conserved (i.e. containing 20% or less gaps) are highlighted using a rainbow colour scheme. The colouring scheme of the reference alignment is mapped directly onto the query alignment so that the user can easily track the difference between reference and query alignments.

Case study
To illustrate the importance of aligning residues approximately correct in terms of misplacement, as opposed to completely wrong, we show different alignments for reference set BB30009 as a case study (Fig 9). In this figure one can observe that this set consists of related protein subfamilies with well-defined core block regions, which are indicated by rainbow- Quantifying the Displacement of Mismatches in MSA Benchmarks colored columns. Furthermore, the subfamilies in this set also appear to be quite different in terms of sequence length.
In this case study we compared the query alignments generated by Clustal Omega and Praline because these two methods generated alignments with the largest contrast in terms of SP and SPdist scores for this set. Both query alignments obtained very similar SP scores: 0.48 for Clustal Omega and 0.46 for Praline. However, the difference in terms of avg_SPdist score between these query alignments is much more discernible, namely 0.59 for Clustal Omega and 0.78 for Praline. This outcome is highly interesting since one would not expect SP and avg_SPdist scoring to produce such different benchmark rankings. It is fair to stress, however, that this case study goes against the general trend in the observed benchmark tests, as Clustal Omega is shown to be the overall best performing MSA method in our benchmark.
A closer examination of reference and query alignments (Fig 9) reveals that Clustal Omega managed to get the intra-subfamilies alignment slightly better than Praline, but overall the Praline alignment contains more but less displaced mismatches. Since the SP scoring scheme treats every mismatch as equal, it is not surprising that the SP score for Clustal Omega alignment is better than the SP score for Praline alignment. However we can also observe from Fig 9 that Clustal Omega totally misaligns the subfamilies, whereas Praline aligns these more correctly. The relatively long inter-subfamily mismatches in the Clustal Omega alignment are appropriately punished in the SPdist scoring scheme resulting in the Praline alignment having a better avg_SPdist score.
This case study exemplifies that even though SP and SPdist scores are conceptually similar, the different treatment of mismatches may result in very different scores, which in turn may lead to different rankings in MSA benchmarking. Therefore, by additionally taking SPdist into account, one can get a more comprehensive assessment of the similarity between query and reference alignments. A thorough similarity comparison between two multiple sequence alignments is not only crucial in an MSA benchmark, but also plays a role in, for example, protein homology modeling: If one were to use the alignments of this case study to make a homology model, the Clustal Omega alignment would yield a considerably worse model as compared to a model obtained from the Praline alignment, since the latter alignment positions most of the conserved regions correctly or within relatively close proximity to each other.
Even though both the SPdist and shift scores take displacements of mismatches into account, the difference in the normalization procedures of these scores results in different behaviour between both scoring schemes. SPdist allows multiple mismatch shift thresholds, whereas the shift score only has a single constant () which is used to differentiate between minor and major shift errors. The multiple mismatch thresholds in SPdist allow the user to determine the graveness of mismatch shifts in a more detailed way. S1 Fig shows that the SPdist score appears more expressive in the case of alignments of distantly related proteins, while the shift score shows more variance in the case of alignments of closely related sequences. In the benchmark case study of BAliBASE reference set BB30009 (Fig 9), where large mismatch shifts occur, the shift score appears unable to distinguish between the query alignments generated by Clustal Omega and Praline as both queries obtained shift score of 0.28.
SPdist analysis in alignment method optimization. We compared the benchmark results of Clustal Omega and MAFFT, the latter of which was ran ran using its default (high-speed) and high-accuracy settings in order to illustrate how the SPdist scoring scheme can be employed to analyse the influence of parameter selection on alignment behaviour. The default setting of MAFFT is primarily aimed at aligning large datasets (where speed is often of the essence), whereas the "auto" setting aims to improve the accuracy by selecting the appropriate alignment optimization strategy which best suits the size of the input dataset. From Fig 10 one can observe a marked improvement on the quality of MAFFT alignments generated using the "auto" setting with respect to the alignments generated using the default setting. Interestingly, although MAFFT "auto" performs better at short distance thresholds, the convergence to a fully correct alignment at long distance thresholds appears to be slightly better for Clustal Omega. However the difference in terms of standard error is not significant. This shows that taking mismatch distance into account helps to obtain a richer performance comparison between the methods.

Conclusion
Based on our observations, we argue that the severity of alignment shifts is often of relevance in numerous biological analyses, e.g. in homology modeling and manual alignment editing. It is therefore important to take the magnitude of mismatch displacement into account in MSA benchmark scoring procedures. Furthermore we showed that by using the SP score together with the SPdist score, one can discriminate the alignment quality better for alignments of distantly related sequences and/or difficult alignments (i.e. alignments that contains many gaps). Whilst classical scoring schemes such as SP scoring are useful for determining the overall similarity between reference and query alignments, scrutinizing the severity of mismatches as implemented in the SPdist score may provide a more salient similarity measure between the reference and the query alignments, particularly from the biological user perspective. The ability of observing detailed alignment quality is further enhanced by monitoring the behaviour of mismatch shifts over different distance thresholds (Figs 3, 4, and 10). For example, if an alignment column shows conserved Cysteine residues, a shift of a few positions can readily be corrected by a user, whereas longer shifts will lead the user to miss the overall conservation of the Cysteine. Therefore we recommend to use SP together with SPdist scoring in order to be able to measure to what extent a multiple sequence alignment actually is displaced with respect to the corresponding reference alignment. Understanding the behaviour of individual MSA methods in this respect provides a step towards optimized manual or semi-automatic MSA correction. Since MSA methods typically contain many adjustable parameters [37], this work provides a novel objective function that can be used when optimizing MSA methods on a reference alignment set.  Table. Average SP, SPdist, and shift scores. The highest score for each category is written in bold. Notice that although SP and SPdist scores are conceptually similar to each other, they may lead to different rankings in an MSA benchmark. By observing the difference in rankings between bins with various cutoffs, one can see that some MSA methods are less prone to create large misalignment shifts than the others. (XLS)


extfieldsofstudy: Physics Medicine

text:
Cell Tracking Accuracy Measurement Based on Comparison of Acyclic Oriented Graphs

Tracking motile cells in time-lapse series is challenging and is required in many biomedical applications. Cell tracks can be mathematically represented as acyclic oriented graphs. Their vertices describe the spatio-temporal locations of individual cells, whereas the edges represent temporal relationships between them. Such a representation maintains the knowledge of all important cellular events within a captured field of view, such as migration, division, death, and transit through the field of view. The increasing number of cell tracking algorithms calls for comparison of their performance. However, the lack of a standardized cell tracking accuracy measure makes the comparison impracticable. This paper defines and evaluates an accuracy measure for objective and systematic benchmarking of cell tracking algorithms. The measure assumes the existence of a ground-truth reference, and assesses how difficult it is to transform a computed graph into the reference one. The difficulty is measured as a weighted sum of the lowest number of graph operations, such as split, delete, and add a vertex and delete, add, and alter the semantics of an edge, needed to make the graphs identical. The measure behavior is extensively analyzed based on the tracking results provided by the participants of the first Cell Tracking Challenge hosted by the 2013 IEEE International Symposium on Biomedical Imaging. We demonstrate the robustness and stability of the measure against small changes in the choice of weights for diverse cell tracking algorithms and fluorescence microscopy datasets. As the measure penalizes all possible errors in the tracking results and is easy to compute, it may especially help developers and analysts to tune their algorithms according to their needs.

Introduction
The cornerstone of many modern live-cell imaging experiments is the ability to automatically track and analyze the motility of cells in time-lapse microscopy images [1,2]. Cell tracking is an essential step in understanding a large variety of complex biological processes such as the immune response, embryonic development, or tumorigenesis [3].
Automated cell tracking can be formulated as a problem of identifying and segmenting all desired cell occurrences and describing their temporal relationships in the time-lapse series. Because cells can migrate, undergo division or cell death, collide, or enter and leave the field of view, a cell tracking algorithm suitable for daily practice must reliably address all these events and provide a data structure that thoroughly characterizes the behavior of tracked objects, which could be either whole cells or cell nuclei depending on the application.
State-of-the-art cell tracking approaches can be broadly classified into two categories [4]: tracking by detection [5][6][7][8][9] and tracking by model evolution [4,[10][11][12][13]. The former paradigm generally involves two steps. First, a cell or cell nucleus segmentation algorithm identifies all target objects in the entire time-lapse series separately for each frame. Second, the detected objects are associated between successive frames, typically by optimizing a probabilistic objective function. In contrast, the latter paradigm solves both steps simultaneously, usually using either parametric or implicit active contour models.
Regardless of the particular algorithm used, its tracking results can be mathematically represented using an acyclic oriented graph. The vertices of such a graph correspond to the detected objects while its edges coincide with the temporal relationships between them. Non-dividing objects have one successor at most, whereas those that undergo division have two or even more successors in the case of abnormal division. Cell lineage tracking results represented by an acyclic oriented graph form a forest of trees in the graph theory terminology.
With the increasing number of cell tracking algorithms, there is a natural demand for objective comparisons of their performance. In general, there are two aspects of cell tracking algorithms, which are worth being evaluated: segmentation accuracy and tracking accuracy. The former one characterizes the ability of an algorithm to precisely identify pixels (or voxels) occupied by the objects in the images. It usually leads to the comparison of reference and computed regions based on their overlap or distances between their contours [14,15]. The tracking accuracy evaluates the ability of an algorithm to correctly detect individual objects of interest and follow them in time.
There are two popular approaches to measuring the tracking accuracy. One approach is based on the ratio of completely reconstructed tracks to the total number of ground-truth tracks [4,16]. The second computes the ratio of correct temporal relationships within reconstructed tracks to the total number of temporal relationships within ground-truth tracks [16,17]. Both approaches quantify, at different scales, how well the cell tracking algorithms are able to reconstruct a particular ground-truth reference. However, they neither penalize detecting spurious tracks nor account for division events, which are often evaluated separately [4,16].
A comprehensive framework for evaluating the performance of the detection and tracking algorithms was established in the field of computer vision [18]. Nevertheless, it targets only topologically stable objects, such as human faces, text boxes, and vehicles. Therefore, it cannot be applied to cell tracking applications because the tracked objects can divide over time or disappear after undergoing cell death. Similarly, another evaluation framework [19], established for comparing the performance of particle tracking methods, does not consider division events, ruling out its ability to evaluate correct cell lineage reconstruction.
In this paper, we propose a tracking accuracy measure that penalizes all possible errors in tracking results and aggregates them into a single value. The measure assesses the difficulty of transforming a computed acyclic oriented graph into a given ground-truth reference. Such difficulty is measured as a weighted sum of the lowest number of graph operations required to make the graphs identical.
The proposed measure can serve not only algorithm developers, but also analysts in order to choose the most suitable algorithm and tune its parameters with respect to all tracking events by optimizing a single criterion. A typical scenario is to create ground truth and evaluate the prospective algorithms on a part of image data, and let the most suitable algorithm run on the rest of it. An alternative way of comparing the performance of algorithms without the need for ground truth has been proposed recently in [20,21]. However, this approach creates the ranking based on a pairwise comparison of the algorithms and therefore the absolute performance of the algorithms remains unknown.
A prototype of the proposed measure has been continuously used in the individual editions of the Cell Tracking Challenge (http://www.codesolorzano.com/celltrackingchallenge/) being an open and ongoing competition focused on objective and systematic comparison of state-ofthe-art cell tracking algorithms [22]. In comparison to [22] where the measure prototype was only sketched out and primarily used as a black-box tool with fixed weights for ranking the algorithms, the contribution of this paper is twofold. First, we provide a rigorous mathematical description of the proposed measure and an extensive study of its behavior, in particular with respect to the choice of weights, based on diverse fluorescence microscopy datasets. It is shown that the proposed measure is robust and stable against small changes in the choice of weights and that the weights can be set to compile rankings strongly correlated with human expert appraisal while reflecting the importance of a particular type of error. Second, a slight modification in the definition of edge-related graph operations, which has no practical impact on the ranking compilation, allows us to formulate the necessary condition for the choice of weights, which guarantees the measure value to be not only the weighted sum of the lowest number of graph operations but also the minimum weighted cost of transforming a computed graph into a given ground-truth reference.

Proposed cell tracking accuracy measure informally
The main purpose of the proposed measure is to evaluate the ability of cell tracking algorithms to detect all desired objects and follow them in time. Although it does not directly evaluate the accuracy of segmented regions, reliable object detection is a very important factor in this measure as well.
In fact, the measure counts the number of all detection as well as linking errors committed by the algorithm. It counts the number of missed objects (FN-false negatives), the number of extra detected objects (FP-false positives), and the number of missed splits (required to correctly segment clusters, NS). Having those three numbers of errors, we can aggregate them into one number as the weighted sum w NS NS + w FN FN + w FP FP with non-negative weights w NS , w FN , and w FP . The ability of the algorithm to correctly identify temporal relationships between the objects is evaluated by counting the number of errors in object linking. Namely, it counts the number of missing links (EA), the number of redundant links (ED), and the number of links with the wrong semantics (EC). These numbers are again aggregated into one number as the weighted sum w EA EA + w ED ED + w EC EC with non-negative weights w EA , w ED , and w EC . The weights are penalties for individual types of errors and can, for example, reflect the manual effort needed to correct the errors in some particular software.
The number of committed errors can be calculated by counting differences between the ground-truth reference and the computed result where each can be mathematically represented by an acyclic oriented graph. From the computational point of view, a critical part of the proposed measure is the existence of a unique way of pairing reference and computed objects (i.e., graph vertices). To this end, we pair a reference object with a computed one if and only if the latter covers the majority of the former, which guarantees the uniqueness of established pairing and thus straightforward computation without any optimization. Interestingly, such a simple test does not exist for the particle tracking problem, where particles are considered volumeless, making similar evaluation procedure computationally unfeasible.
In the rest of this section, the basic terminology and notation necessary to establish a connection between cell lineage tracking results and acyclic oriented graphs is introduced, and the proposed measure is formally defined.

Basic terminology and notation
Let N be the number of frames in a time-lapse series, L be a set of object labels shared among all frames, and ? = 2 L be a background label. Let a marker M t i be a set of pixels or voxels of a unique label i 2 L in the t-th frame, t 2 {0, . . ., N − 1}, related to a particular object. Such markers can be manually created by experts, automatically computed by cell tracking algorithms, established by combining both in an edit-based framework [23], or inherently generated using a simulation toolkit as digital phantoms [24,25]. A track θ i is defined as the longest temporal series of markers M t init i ; Á Á Á ; M t end i , 0 t init t end N − 1, without temporal gaps. When a particular object temporarily disappears from a frame, the corresponding track terminates and a new one with a unique label is established once the object reappears in another frame. Analogously, when a particular object undergoes division, the track of the mother object terminates and new daughter tracks with unique labels are initiated. In either case, the new established tracks are descendants of the terminated tracks. We denote a set of all tracks for the particular time-lapse series with the symbol Θ. To keep the information about relationships between individual tracks, let P : Y ! L [ f?g be a parent function defined as : For simplicity, we also define two other functions, I : Y ! f0; Á Á Á ; N À 1g and T : Y ! f0; Á Á Á ; N À 1g, returning zero-based indices of the initial and terminal frames, t init and t end , for each track. Finally, tracking results for the particular time-lapse series can be expressed as a quadruple ðY; P; I ; T Þ. Any quadruple ðY; P; I ; T Þ can be directly transformed into an acyclic oriented graph G = (V, E) where a set of vertices V is composed of all markers present in tracks θ i 2 Θ and a set of oriented edges E & V × V represents temporal relationships between the markers. More precisely, a pair ðM t 1 i ; M t 2 j Þ is an edge of the graph G if and only if either i = j^t 2 = t 1 + 1 or In the former case, the edge connects two successive markers within a single track, whereas the terminal marker of θ i is linked to the initial one of θ j in the latter case. Hereinafter, we refer to the former edge as track link and to the latter one as parent link. Over the set of edges E, we define a function S : E ! fT; Pg describing the semantics of an edge e 2 E as SðeÞ ¼ T for track links and SðeÞ ¼ P for parent links. Note that the orientation of the edges follows the ascending temporal ordering of markers within as well as between tracks, which ensures acyclicity of the graph G. An example of the graph G is depicted in Fig 1. and brevity, we denote the reference graph vertices as R t i , and the computed graph vertices as C t j . To determine whether a reference marker R t i was found, we exploit a simple binary test checking whether a computed marker C t j covers the majority of the reference marker R t i . Note that each reference marker can be assigned to one computed marker at most, whereas one computed marker can have multiple reference markers assigned using this detection test. For instance, the latter happens when a cell division was not detected by the algorithm, leaving the daughter cells clustered in a single marker C t j . In the case of positive detection test, we write R t i ⋐ C t j and say that R t i is assigned to C t j . Analogously, we consider edges e R ¼ ðR t 1 i ; R t 2 j Þ and e C ¼ ðC t 1 k ; C t 2 l Þ to match if and only if the corresponding vertices have positive detection tests (i.e., e R ⋐ e C $ R t 1 i ⋐ C t 1 k^R t 2 j ⋐ C t 2 l ). Based on the detection test, we classify the reference vertices V R as either true positive V TP R or false negative V FN R in the following way: • True positives: the correctly detected objects (i.e., reference markers assigned to a computed marker): We denote the number of true positive vertices as TP ¼ jV TP R j.
• False negatives: the missed objects (i.e., reference markers not assigned to any computed marker): We denote the number of false negative vertices as FN ¼ jV FN R j.
Note that each vertex is included in exactly one set and their union contains all reference vertices (i.e., Similarly, we classify the computed vertices V C as true positive V TP C , false positive V FP C , or non-split vertices V VS in the following way: • False positives: the extra detected objects (i.e., computed markers without any reference marker assigned): We denote the number of false positive vertices as FP ¼ jV FP C j.
• Non-split vertices: the computed markers with more than one reference marker assigned: We denote the number of non-split vertices as VS ¼ jV VS C j. When m reference markers (m > 1) are assigned to a single computed marker, m − 1 split vertex operations need to be performed to locally equalize the number of vertices of the reference and computed graphs. These operations decompose the computed marker into m non-empty, disjoint markers, such that each of the m reference markers is assigned to exactly one of them. The total number of split vertex operations can be easily obtained as the difference between the number of true positive reference vertices and the number of the computed graph vertices with a reference marker assigned: • True positives: the computed markers with exactly one reference marker assigned: Note that each computed vertex is included in exactly one set and their union contains all computed vertices (i.e., . Knowing vertex classification, we define edge-related errors by comparing the reference edges with those in an induced subgraphĜ C ¼ ðV C ;Ê C Þ of the computed graph G C by a vertex setV C ¼ V TP C , which is formed of the uniquely matching computed vertices (i.e., those with exactly one reference vertex assigned) and all their incident edges (i.e., . First, we define the set of redundant edges in the computed graph. These are the induced subgraph edges without any counterpart in the reference graph: We denote the number of redundant edges as ED ¼ jE FP C j. Note that the computed edges attached to a false positive vertex or non-split vertex are not included because they are inherently removed together with deleting the false positive vertices or with splitting the non-split vertices.
Analogously, we define the set of missing edges in the computed graph. These are the reference graph edges without any counterpart in the induced subgraph: We denote the number of missing edges as EA ¼ jE FN C j.
Finally, we define the set of edges with wrong semantics in the computed graph. These are the matching edges between the reference graph and the induced subgraph, which differ in semantics: We denote the number of such edges as EC ¼ jE CS C j. The transformation of any computed graph into the reference one involves the following procedure. First, one can build a binary matrix with |V R | columns and |V C | rows containing the results of the detection test given by Eq (2) for each pair of vertices ðR t i ; C t j Þ, R t i 2 V R , C t j 2 V C . Due to the majority overlap criterion, each reference vertex is assigned to at most one computed vertex. Therefore, either none or only one match can appear in each column, resulting in either false negative or true positive classification of the respective reference vertices. Similarly, the matrix rows without any match reveal false positive vertices, and those with multiple matches correspond to the computed vertices that need to be split. Note that the number of matches in each such row decremented by one gives the number of splits that must be executed on the particular non-split vertex to locally equalize the number of vertices in both graphs. Because the reference as well as computed markers are spatially localized, the vertex matching is unique. In total, NS split vertex operations, FN add vertex operations, and FP delete vertex operations need to be performed to have the vertex sets of both graphs matching. Subsequently, we remove redundant edges, add missing ones, and finally correct those with wrong semantics. These operations are also unique because any of them cannot be replaced by a reasonable combination of the others. This requires ED delete edge operations, EA add edge operations, and EC alter the edge semantics operations, respectively.
The weighted sum of the executed operations is considered as the cost of transforming the computed graph into the reference one (AOGM measure): With the assumption of non-negative weights along with at least one weight being positive, the AOGM measure is bounded below by zero. Its value is equal to zero when a computed graph is identical to the ground-truth reference or when it contains only errors penalized with a zero weight. The AOGM value increases, theoretically to infinity, with the increasing complexity of the transformation that converts the computed graph into the reference one, where complexity is judged with respect to operations penalized with non-zero weights. The higher the AOGM value is, the worse output an algorithm has provided and the worse its ranking is. An example of calculating the AOGM measure is illustrated in Fig 2. The ability of an algorithm to detect all important objects can be measured using the AOGM measure when keeping only vertex-related weights positive (i.e., w NS , w FN , w FP > 0; w ED = w EA = w EC = 0). We further refer to such variant of the AOGM measure as AOGM-D. Analogously, when keeping only the edge-related weights positive (i.e., w NS = w FN = w FP = 0; w ED , w EA , w EC > 0), the AOGM measure evaluates the ability of an algorithm to follow objects in time (i.e., its association skills). We further refer to such variant of the AOGM measure as AOGM-A.

Results and Discussion
In this section, we first discuss the minimality of the proposed measure, describe testing data used for the experimental evaluation, and present their properties in terms of tracking errors. Next, we study how sensitive the AOGM measure is to the choice of weights, and how its behavior coincides with human expert appraisal. Finally, we discuss the size of necessary ground truth.

Measure minimality
The computation of the AOGM measure is a deterministic procedure, which always returns a single value derived from the number of differences between the computed and reference graphs. Because we can, thanks to the detection test given by Eq (2), uniquely match the graphs, the set of errors is also unique and can be easily determined. The allowed graph operations to correct the errors are add, delete, and split a vertex and add, delete, and alter the semantics of an edge. All but the split vertex operation are independent in a sense that they cannot be substituted by a reasonable combination of the others and each operation directly corrects one error. A non-split vertex with m reference vertices assigned could also be corrected using the following sequence of operations: delete the non-split vertex and add m new vertices. If the cost of this sequence is higher than the cost of m − 1 split vertex operations, the AOGM value is equal to the minimum cost necessary to transform the computed graph into the reference one. We formulate the necessary condition, which guarantees the AOGM value to be not only the weighted sum of the lowest number of graph operations but also the minimum weighted sum, as where m ? is the maximum number of reference vertices assigned to a single non-split vertex among all non-split vertices. For practical purposes, the minimality condition can be weakened to w NS w FN to be independent of the input graphs. Note that the measure prototype used in the first Cell Tracking Challenge calculated the edge-related errors over the entire computed graph and therefore, some of the erroneous edges attached to non-split vertices were inherently corrected along with the vertex splitting while the others were not. This complicated the measure minimality reasoning because the cost of correcting a single non-split vertex may have also involved the penalties for a variable number of delete edge and alter the edge semantics operations on the left side of Eq (13) and for a variable number of add edge operations on the right side of Eq (13), and therefore we slightly modified the measure definition. By calculating the edge-related errors over the induced subgraph only, the correction of vertex-related errors is clearly separated from the edge-related ones, allowing minimality condition to be formulated only using the vertex-related weights. In general, the proposed measure penalizes the splitting errors slightly higher than the measure used in the first Cell Tracking Challenge, in particular due to a slight increase in the number of add edge operations, being in average approximately 1.88 per non-split vertex, caused by inherent deleting of all edges attached to non-split vertices during their splittings. However, such increase does not have any practical impact on the compiled rankings and they remain unchanged for all algorithms and datasets used in the first Cell Tracking Challenge regardless of the way how the edges attached to non-split vertices are handled.

Testing data
We experimentally validated the proposed measure using the results of four cell tracking algorithms (COM-US, HEID-GE, KTH-SE, and PRAG-CZ as named in [22]) that participated in the first Cell Tracking Challenge. These algorithms provided complete tracking results for the entire competition dataset repository [22], which allowed us to study the AOGM measure behavior under distinct scenarios involving, in particular, various cell phenotypes (i.e., shape, density, and motion model) and acquisition configurations (i.e., imaging system, image data dimensionality, and time step). The basic properties of the competition datasets are listed in Table 1. As our primary interest lies in analyzing the AOGM measure behavior with respect to the choice of weights in Eq (12), rather than in compiling a ranking of the competing algorithms and in discussing their strengths and weaknesses, details about the algorithms are omitted. We further refer to them as A1, A2, A3, and A4 in a random but fixed order.
In this work, two types of ground-truth references were used, depending on the origin of a particular dataset, referred as real and synthetic datasets. They were adopted without any modification from [22]. For the real datasets, three experts were requested to place a quintessential marker inside every object and establish temporal relationships between the markers to provide an acyclic oriented graph for each time-lapse series. Because it is widely documented in the literature that humans commit errors when performing manual tracking [26,27], the reference graphs were created based on majority voting to reduce human errors. For the synthetic datasets, the reference graphs were inherently generated by a simulation toolkit [24,25] along with the time-lapse series to be analyzed.

Distribution of tracking errors
We focused on the statistical analysis of errors in tracking results produced by the four tested algorithms to reveal the distribution of errors under distinct scenarios displayed in the competition datasets ( Table 1). The average number and standard deviation of particular errors, normalized either per ground-truth vertex or per ground-truth edge to facilitate the comparison between datasets with different numbers of objects, are listed in Table 2. The distribution of errors that the algorithms made in particular datasets is shown in Fig 3. The measured values revealed three key observations. First, from the computed standard deviations of errors, it can be observed that the tested algorithms made different types of errors in each particular dataset. Second, there is a preponderance of the add edge operations over the add vertex ones because correcting a false negative vertex requires adding at least one edge to integrate the vertex into a particular track, except when a whole track is formed of a single vertex. Furthermore, missing edges may need to be added also due to inherent deletion of all edges attached to non-split vertices. Third, the predominance of a specific error seems to be dataset-dependent: the false negative detection predominates in Datasets 5 and 8, the incorrectly clustered objects in Dataset 3, and the false positive detection in the remaining datasets.  Table 2. The average number and standard deviation of errors in tracking results normalized either per ground-truth vertex or per ground-truth edge depending on a particular error type. The second and third columns list the aggregated numbers of vertices and edges in the ground-truth references per dataset.

Sensitivity to the choice of weights
We investigated the behavior of the AOGM measure with respect to the choice of weights in Eq (12). As a reference weight configuration, we adopted the weights used in the first Cell Tracking Challenge, which reflected the effort needed to correct a particular error manually: w NS = 5 for vertex splitting, w FN = 10 for vertex adding, w FP = 1 for vertex deleting, w ED = 1 for edge Cell Tracking Accuracy Measurement Using Acyclic Oriented Graphs deleting, w EA = 1.5 for edge adding, and w EC = 1 for altering the edge semantics. Note that such weight configuration satisfies the minimality condition given by Eq (13) for any m ? . We wanted to know how a change in the setting of weights would have affected the reference rankings of the four tested algorithms in terms of tracking ability. To this end, we computed sectors in weight space leading to the same rankings. The analysis was carried out with respect to the pair of the two most influential weights, w NS and w FN . The other weights were fixed at their original values, serving as normalization factors. By solving systems of linear inequalities in two unknown weights, w Ã NS and w Ã FN , within the rectangular domain h0, 10i × h0, 20i, we studied how many transpositions in the obtained rankings occurred, compared to those compiled for the reference configuration ðw Ã NS ; w Ã FN Þ ¼ ð5; 10Þ. The domain was chosen to range from 0 up to a double magnitude of the reference weight in each axis, being centered at the reference weight configuration. Due to the linearity of Eq (12), the solutions of these systems form polygonal sectors, each consisting of the configurations leading to the same ranking. The results for all datasets are shown in Fig 4. The borderlines between the sectors correspond to the configurations, for which at least two algorithms gained the same AOGM value. Almost all tested configurations led to no more than one transposition in the reference rankings, and even in the half of the datasets, we observed practically no changes in the reference rankings. We also carried out a similar study for another two pairs of weights, namely (w FP , w FN ) and (w EA , w FN ), which turned out to be the weights of most frequent errors across all the datasets. The number of sectors and their relative areas are listed in Table 3. Note that in all of the cases the largest sector was formed of the configurations leading to the reference ranking. The sectors with no more than one transposition occupied the whole studied domain in all but three cases. However, in these three cases, they covered more than 99% of the domain.
Based on the information depicted in Fig 4 and Table 3, it can be concluded that rankings of the tested algorithms compiled using the AOGM measure are highly robust against small changes in the choice of weights. This demonstrates that the reference configuration ðw Ã NS ; w Ã FN Þ ¼ ð5; 10Þ compiles the predominant rankings of the tested algorithms for diverse datasets (Table 1) because it belongs to the sectors of the largest area (the pink sectors in Fig 4, and the column 0 in Table 3).

Comparison with human expert appraisal
To validate the AOGM measure against human expert appraisal, we asked three independent human experts to rank the algorithms according to two criteria, which were formulated as (1) Table 3. The number of sectors with different rankings and their relative areas within the rectangular domain centered at the reference weight configuration for individual datasets. The results for three different pairs of weights, (w NS , w FN ), (w FP , w FN ), and (w EA , w FN ), are presented. The columns 0, 1, and 2 contain the relative areas of the sectors with the particular number of transpositions given by the column label. "the ability of the algorithm to detect all important objects", shortly detection performance, and (2) "the ability of the algorithm to follow objects in time", shortly association performance. We did not provide any special instructions on how the experts were supposed to rank the algorithms and we let their decision completely on their opinion. In particular, we did not provide them with the ground truth and did not require any quantification, because it would tend toward the manual counting of errors in the results, which is basically what the AOGM measure does automatically, provided the ground truth exists. The experts allocated points to the algorithms from 0 to 3 (3 means the best) for detection and association separately and the total score in each category was computed by summing up the points. The overall score was obtained by pure addition of the total score in each category. The rankings compiled using AOGM-D and AOGM-A variants of the AOGM measure were compared with human expert rankings on all the six real 2-D sequences. The results are summarized in Tables 4, 5 and 6. They indicate high correlation, with the Kendall rank correlation coefficient τ ranging from 0.78 up to 0.83, between the rankings on the basis of the proposed method and human expert opinions. The perfect agreement occurred for three analyzed sequences (Dataset 1/Sequence 1, Dataset 1/Sequence 2, and Dataset 4/Sequence 2).
The difference in the rankings for Dataset 4/Sequence 1 exists because the human experts preponderated the existence of non-split objects and false positive objects over the existence of false negative objects during the evaluation process. Tables 7 and 8 show that the AOGM measure could have compiled the same ranking as that compiled by the human experts if their Table 4. The detection ranking comparison of the four tested algorithms (A1-A4) with human expert appraisal on real 2-D dataset sequences using AOGM-D with weights w NS = 5, w FN = 10, w FP = 1; w ED = w EA = w EC = 0. The black bullets mark the matches in the rankings.  Table 5. The association ranking comparison of the four tested algorithms (A1-A4) with human expert appraisal on real 2-D dataset sequences using AOGM-A with weights w NS = w FN = w FP = 0; w ED = 1, w EA = 1.5, w EC = 1. The black bullets mark the matches in the rankings. preference on the type of errors was reflected in the weights. The same observation holds for Dataset 5/Sequence 1. Similarly, the swap in AOGM-A for Dataset 4/Sequence 1 exists because AOGM-A puts more emphasis on the missing edges rather than on their semantics. If the weights were accordingly altered, we could have obtained the human expert ranking. In Dataset 5/Sequence 2, A3 was slightly worse than A4 with respect to all types of errors. However, the human experts ranked this algorithm higher. We suppose it is because the human expert comparison was sometimes very difficult and subjective, especially if the algorithms behaved Table 6. The overall ranking comparison of the four tested algorithms (A1-A4) with human expert appraisal in terms of detection on real 2-D dataset sequences using AOGM with weights w NS = 5, w FN = 10, w FP = 1; w ED = 1, w EA = 1.5, w EC = 1. The overall AOGM values and expert points were obtained by pure addition of the respective quantities for the detection (Table 4) and association (  Table 8. The overall ranking comparison of the four tested algorithms (A1-A4) with human expert appraisal on real 2-D dataset sequences using AOGM with weights w NS = 10, w FN = 1, w FP = 10; w ED = 1, w EA = 1.5, w EC = 1. The overall AOGM values and expert points were obtained by pure addition of the respective quantities for the detection (Table 7) and association ( similarly, and it was impracticable to manually compute the precise number of errors due to image data complexity (e.g., Dataset 5 contains hundreds of cells per frame). The experiments demonstrate that AOGM measure can compile expected ranking, being strongly correlated with human expert one, provided the weights reflect the importance of a particular type of error.

Size of ground truth
Typical cell tracking experiments produce hundreds to thousands of 2-D or 3-D images. A decision on what algorithm to use and how to optimally set its parameters is often not an easy task.
A widely adopted approach is to run several algorithms with different parameter settings and visually compare their results. The size of image data often makes such evaluation practicable on its limited subset only, although being tricky, especially for 3-D experiments, and subjective due to inter-operator variability and a graphical user interface used.
An alternative approach is to create ground truth for a part of image data and use the AOGM measure for objective and accurate evaluation of the algorithms. In general, the creation of ground truth is laborious, although the amount of work can be reduced by adopting an edit-based framework [23]. Nevertheless, such effort pays off soon when multiple algorithms need to be evaluated and their parameters optimally tuned.
However, any general recommendation on the size of ground truth to get unbiased evaluation of the algorithm performance can be hardly stated. It depends on the complexity of image data and particular application. Instead, we computed the temporal evolution of the AOGM measure for each algorithm and each real dataset (Fig 5). It can be observed that the number of frames after which the rankings of the tested algorithms stabilized varied across the datasets. It can also be observed often a linear increase in the AOGM value, indicating the tested algorithms committed approximately the same number of errors over time.

Conclusions
In this paper, the AOGM measure, an accuracy measure for objective and systematic comparison of cell tracking algorithms that are capable of providing segmentation of individual cells rather than simplifying them as single-point objects, has been defined and analyzed. Treating tracking results as an acyclic oriented graph, the proposed measure assesses how difficult it is to transform a computed graph into a given ground-truth graph. The cost of such a transformation is defined as the weighted sum of the lowest number of graph operations needed to make the graphs identical. The behavior of the AOGM measure was analyzed on tracking results provided by four state-of-the-art cell tracking algorithms that participated in the first Cell Tracking Challenge. The performed analyses verified its robustness and stability against small changes in the choice of weights for diverse fluorescence microscopy datasets.
As the weights chosen in Eq (12) are not biologically motivated, they reflect the effort needed for performing a particular graph operation, the AOGM measure is application-independent. A different choice of weights can serve purposes other than objectively and systematically comparing the performance of multiple cell tracking algorithms as presented in this paper. For example, the AOGM measure might also be useful to tune the parameters of each single module involved in a cell tracking algorithm under development for a specific dataset. Its AOGM-D variant allows one to determine the optimal parameters of a cell detection module. Furthermore, setting the weight w NS to 1 and the others to 0 allows one to optimize a cluster separation module. Cell Tracking Accuracy Measurement Using Acyclic Oriented Graphs In addition to being robust, stable, and flexible, the AOGM measure is also universal. The measure penalizes all possible errors in tracking results, not concentrating on only a single, often application-limited aspect of cell tracking as the existing approaches do [4,16,17]. It can be applied to datasets with various characteristics, even nearly degenerated cases such as those showing no division or containing only a single cell. Furthermore, it can be used for evaluating the performance of any cell tracking algorithm irrespective of its nature because it evaluates its final output. Therefore, developers can use the AOGM measure to easily compare the performance of a cell tracking algorithm under development to that of existing algorithms, whereas analysts can use it to determine the optimal parameters of a chosen algorithm for a dataset to be analyzed. The software for computing the AOGM measure is made publicly available at http://cbia.fi.muni.cz/aogm/ or as a supplementary material S1 Software, free of charge for noncommercial and research purposes.
Supporting Information S1 Software. AOGMMeasure. This package contains a routine for computing the AOGM measure. (ZIP) S1 File. Supporting Information File. This supporting information file contains all relevant data used for generating the results described in the manuscript. (XLSX) Author Contributions