VarDetect: a nucleotide sequence variation exploratory tool

Background Single nucleotide polymorphisms (SNPs) are the most commonly studied units of genetic variation. The discovery of such variation may help to identify causative gene mutations in monogenic diseases and SNPs associated with predisposing genes in complex diseases. Accurate detection of SNPs requires software that can correctly interpret chromatogram signals to nucleotides. Results We present VarDetect, a stand-alone nucleotide variation exploratory tool that automatically detects nucleotide variation from fluorescence based chromatogram traces. Accurate SNP base-calling is achieved using pre-calculated peak content ratios, and is enhanced by rules which account for common sequence reading artifacts. The proposed software tool is benchmarked against four other well-known SNP discovery software tools (PolyPhred, novoSNP, Genalys and Mutation Surveyor) using fluorescence based chromatograms from 15 human genes. These chromatograms were obtained from sequencing 16 two-pooled DNA samples; a total of 32 individual DNA samples. In this comparison of automatic SNP detection tools, VarDetect achieved the highest detection efficiency. Availability VarDetect is compatible with most major operating systems such as Microsoft Windows, Linux, and Mac OSX. The current version of VarDetect is freely available at .


Background
Following completion of the human genome project, detection and discovery of single nucleotide polymor-phisms (SNPs) is at the forefront of genomic research. The discovery of SNPs may help to identify causative gene mutations in monogenic diseases as well as SNPs associ-from Asia Pacific Bioinformatics Network (APBioNet) Seventh International Conference on Bioinformatics (InCoB2008) Taipei, Taiwan. 20-23 October 2008 ated with predisposing genes in complex diseases [1,2]. Most fluorescence based sequencers produce nucleotide signals (chromatograms) that must be base-called in order to detect the SNP or point mutation. The terms SNP and point mutation are considered synonymous for the algorithm described in this paper. The efficiency of nucleotide variation detection relies mainly on the accuracy of bioinformatic software used to base-call the chromatograms [3][4][5][6]. However, most base-calling tools developed for conventional sequencing may not be suitable for SNP detection because they usually misinterpret chromatogram traces at heterozygous base positions. The common sequencing artifacts, which cause most standard base-calling tools to miscall the chromatogram traces, include: 1) polymerase slipage, 2) loss of resolution, 3) contamination and 4) dye blob [7]. SNP discovery would be greatly accelerated if a reliable, automatic SNP discovery tool was available. A commercial automatic SNP detection program called Mutation Surveyor (SoftGenetics) was recently developed utilizing patented anti-correlation technology to increase the efficiency of SNP and mutation detection [8]. Non-commercial programs for SNP detection include PolyPhred (used together with Phred [6,9], Phrap and Consed programs [10]), novoSNP and Genalys. PolyPhred was designed in conjunction with the well-known Phred and Phrap programs to base-call and assemble input chromatograms prior to SNP detection and visualize the results using the Consed program. The current version of PolyPhred is 6.11 beta at the time of writing. The novoSNP program adopts three independent cumulative scores to identify SNPs, and is able to identify more true SNPs (lower false negative rate) than PolyPhred (version 3) [11]. Genalys software attempts to minimize the number of incorrectly assigned (false positive) SNPs by using peak base ratios and surrounding peak information to identify SNPs [12]. Despite the sophistication of the mathematical models widely used in these algorithms, they still report an unacceptably high number of false negative and/or false positive SNPs.
In this study, we present VarDetect, a sequence variation exploratory software to detect SNPs efficiently from fluorescence based chromatogram data. VarDetect supersedes existing automatic SNP detection tools through utilization of rules which account for the common sequence reading artifacts, combined with pre-calculated peak content base ratios. Furthermore, SNPs can be detected by this software using sequencing data obtained from single, or two-pooled DNA samples.

Results and discussion
VarDetect main graphical user interface (GUI) contains four panes and one quick access toolbar ( Figure 1). By pressing the first button on the left of the toolbar (Figure  1a), the program wizard will respectively guide users to 1) load a reference sequence either in an eXtensible Markup Language (XML) or a FASTA format, 2) add input fluorescence based chromatogram traces (.ab1 or .scf extension), and 3) adjust the software parameters such as noise level, heterozygosity level, accepted trace quality and quality difference (δ). At the conclusion of these three steps, Var-Detect will perform the SNP detection resulting in the graphical traces (Figure 1b), list of predicted SNPs ( Figure  1c), detailed information of selected SNPs from the list (Figure 1d), and a whole-map view of the aligned input sequences against the reference sequence ( Figure 1e). The whole-map viewer allows users to oversee where the predicted SNPs are located. If the annotated reference sequence (XML format) is used, users can visualize the predicted SNPs to determine if they are novel. As an example, Figure 1e shows the nine SNPs which were previously reported on the reference sequence as short yellow vertical bars on a white horizontal line. The blue horizontal lines indicate the aligned sequences to the reference sequence. The short vertical bars running through these blue lines indicate the predicted SNPs. Furthermore, VarDetect provides an export function that can generate a SNP report using the NCBI SNP submission format.
We compared the features and accuracy of VarDetect with PolyPhred (version 6.11 beta) [10], Genalys (version 3.3.23a) [12], novoSNP (version 2.0.3) [11] and Mutation Surveyor (trial version 3.23) [8]. The efficiencies of these tools were measured in terms of false positive (FP) and false negative (FN) SNP counts. The term true positive count (TP) conveys the number of predicted SNPs which correspond to the actual, prior verified SNPs in the region of interest. The FP is defined as the number of SNPs predicted by the software which are in fact not present. The FN is the number of true SNPs which could not be identified by the software. Hence, the efficiency can be calculated in terms of precision (TP/(TP+FP)) and recall (TP/ (TP+FN)) and F-score (2 × (precision + recall)/(precision + recall)) [13]. The performance comparison was made by comparing the results obtained by running these tools against chromatogram traces used in our previous human SNP discovery of fifteen candidate genes. In that report, 171 SNPs were validated from a total of 77 reads [14]. We used these validated SNPs as our TP. Table 1 tabulates the experimental results on the fifteen candidate genes. To avoid biases caused by software parameter adjustment, automated detection was performed with the default parameter settings. The detail of these software default parameters are listed at the software website hosted at http://www.biotec.or.th/GI/tools/vardetect.
Of the five tools, VarDetect and novoSNP yielded the lowest false negative counts (23 and 11 respectively), which is of paramount importance in most SNP/mutation discovery projects. PolyPhred reported the fewest false positives (25), while novoSNP reported the most false positives (4539). VarDetect had the second lowest false positive count (149). The chromatograms and analysis results of this experiment can be obtained and visualized online in scalable vector graphics (SVG) format from the above Var-Detect website. Finally, we measured the software precision and recall ratios in order to compare them using Fscore [13]. These scores are presented in the last row of Table 1 where VarDetect had the highest F-score (63.25%), slightly greater than PolyPhred (62.94%). Despite the very similar efficiency, VarDetect reported considerably fewer FNs (about 3.5 times lower) than PolyPhred; hence, VarDetect is preferable for SNP discov-ery. This implies higher recall rate (higher sensitivity) which is the ability to detect SNPs even from low or ambiguous chromatogram signals. Nonetheless, higher recall rate has a tradeoff that is having lower precision. In other words, a lot more FPs would be predicted as the program sensitivity is improved. Mutation Surveyor, Genalys and novoSNP had low efficiencies due to much higher numbers of false positives (31.84%, 6.93%, and 6.56%, respectively). Overall, VarDetect is superior to other automatic SNP discovery tools because both FN and FP counts are minimized.
It is worth noting that most FPs picked up by VarDetect in this experiment came from two fundamental problems which are very difficult for the proposed heuristics to Illustration of VarDetect's graphical user interface Figure 1 Illustration of VarDetect's graphical user interface. The graphical user interface comprises four panes and a quick access toolbar: a) toolbar with a wizard button located on the left-most b) graphical view of input chromatogram traces c) list of predicted SNPs d) SNP information window and e) whole-map view.
resolve. First VarDetect utilizes the vicinity quality concept to implicate how much confidence we can identify a SNP since the quality of the observed base will be bad for both SNP and poor quality signal (see Figures 2 and 3). Because we have preset the software to have a higher recall rate and given that the vicinity quality is good, VarDetect would likely predict the observed base as SNP even with a small trace of turning point. Secondly, VarDetect does not explicitly trim both ends of the input chromatograms which usually contain low quality trace signals. Hence, with current high sensitivity settings, VarDetect tends to predict more FPs. To get around these problems, one must fine tune the parameters, e.g., noise threshold, vicinity/ observed quality acceptance etc., to alter the prediction behavior of VarDetect. Furthermore, we also compared different software features and tabulated them in the feature matrix table (Table  2). VarDetect, Genalys, novoSNP and Mutation Surveyor are more user-friendly since they are not dependent on other external software packages such as Consed, Phrap and Phred. They also offer an embedded graphical user interface while PolyPhred requires an external Consed tool installed to visualize the results. VarDetect and Poly-Phred can use batch SNP detection with command line interface (CLI), while the other softwares have only Graphical User Interface (GUI). In terms of processing allele frequencies from two-pooled DNA samples, VarDetect and Genalys are the only two softwares that offer such a feature. Furthermore, VarDetect can read in a reference sequence and corresponding validated SNP locations on the reference sequence to compare with the predicted SNPs, allowing indication of possible novel SNPs. Reference sequences and the corresponding previously reported SNPs can be downloaded directly from the ThaiSNP database [15] in XML format.

Conclusion
We present the framework of a novel algorithm to interpret (base-call) fluorescence based chromatograms and efficiently detect the corresponding nucleotide variations in an automatic fashion. In this framework, three main heuristic procedures are employed: 1) Partitioning and Re-sampling (PnR) algorithm that may be used to basecall the bases with ambiguous signal, 2) calculation of the observed signal intensity ratio(Q o ) and vicinity intensity ratio (Q v ) and utilizing the differences between Q v and Q o (quality difference) to check whether the heterozygous Our testing set comprises 77 exonic contigs from 15 human genes. Of these 77 contigs, there are 171 SNPs experimentally validated by the experts [14]. novoSNP has lowest false negative result while generates heighest false positive result (4539). The number of detected SNPs and the false positive result identified by VarDetect are better than Genalys, novoSNP and Mutation Surveyor. Although VarDetect has worse performance on false positive, its SNP identification coverage (148 out of 171) is much higher than that of PolyPhred (90 out of 171) resulting in higher false positive and lower false negative result than PolyPhred' s. From these experiments, VarDetect had the best F-score (63.25%) that slightly greater than PolyPhred (62.94%). Mutation Surveyor, Genalys and novoSNP had low F-score that had 31.84%, 6.93% and 6.56%, respectively.
peaks are correctly called by the PnR algorithm, and 3) conversion of the chromatogram inputs to numeric code using CodeMap so that the variation can be correctly identified by computer.
The experimental results showed that VarDetect is more efficient than other existing tools, namely PolyPhred, nov-oSNP, Genalys, and Mutation Surveyor for detecting SNPs. VarDetect's heuristics minimize both false positive and negative errors reducing the effort needed to detect and validate SNPs, making it the tool of choice for automatic SNP detection. Furthermore, VarDetect offers the most features including the ability to detect SNPs from pooled DNA samples and the use of XML annotated reference sequence to cross check the SNP discovery results within the tool without using external applications. Var-Detect is platform independent since it was implemented in Java, allowing it to run on all major operating systems without recompiling the source codes.  Therefore, is 500/970 = 0.515. If the middle peak con-

Methods
Calculation of vicinity peak intensity ratio of the base-call position (arrowed)   tains only an adenine signal, the ratio would be 500/500 = 1. This calculated ratio can be used to filter out parts of chromatogram that contain low quality peaks. The observed peak content base ratio is then used to calculate an average peak content base ratio of the  Figure 3). The peak content base ratio of the i th base from the above definition reflects the changes that occur when peak intensity is altered by having two or more different signals coincident at the i th position. Therefore, we can tentatively identify the heterozygous state at the i th base by observing the difference between the observed and vicinity peak content base ratios.
However, the problem of sampling the values at the i th position is that the intensity ratio may be interpreted as being a polymorphic site or just an unreadable nucleotide with noisy signal data. Figure 4 demonstrates that even though the intensity ratios from both nucleotide positions (boxes 1 and 2) are identical, the base-call for box 1 can be made more accurately than box 2. The various signals at box 2 most likely reflect sequence reading artifacts, rather than true sequence polymorphism. VarDetect considers peak shapes in addition to the intensity ratio, which can help to reduce these reading artifacts. Regardless of the nature of the base-calling ambiguity, the concept of intensity ratios for the base under observation and adjacent bases provides a good foundation for the SNP discovery process. The sub-procedures used for automatic SNP discovery and verification comprise the following three steps: 1. Reading nucleotides (base-calling) from chromatogram traces 2. Alignment of input sequences to the reference sequence

Detection of SNPs and insertions/deletions (indels)
Each step also comprises a two-fold process, namely rough and fine data processing to ensure the accuracy of the resulting data. Processed data collected from base-calling and trace-alignment are analyzed using the aforementioned intensity ratio concept.
The peak intensity ratio approach may not correctly base-call different peak patterns Figure 4 The peak intensity ratio approach may not correctly base-call different peak patterns. The Q o value from both boxes 1 and 2 are identical (0.551); however the black peak in box 2 is misinterpreted as a primary peak, since it clearly overshoots from the adjacent base position.

Reading nucleotides (base-calling) from chromatogram traces
Currently, most chromatogram trace data come in two formats, namely .ab1 and .scf. Although the .ab1 extension format is proprietary, there are numerous bioinformatics tools, e.g., Phred [6,9], PolyPhred [10], 4Peaks [16], FinchTV [17], Genalys [12], novoSNP [11], Mutation Surveyor [8] and software libraries including BioPerl [18,19], BioJava [19], and BioPython [19] that can read this file format. We used Java and BioJava to develop Var-Detect because of its high portability. From Equations 1 and 2, the intensity ratios can be pre-computed while base-calling of an input chromatogram is being processed. Each nucleotide position contains intensity values of A, T, G, and C as required by the definitions. The algorithm to calculate such ratios strictly follows these definitions.  Figure  5). The base-call array stores the sampling interval where information from the nucleotide content arrays should be read. Figure 5 shows sample array contents which can be read as TGGAC by observing the highest amplitude content, 5 in this case, specified by the positions held in the B [] array. If more than one nucleotide content array has significant amplitude, the base-call at this position could be interpreted as heterozygous. Base-calling configuration settings such as the noise level and heterozygosity ratio are preset in base-calling programs. These settings are however not optimal for identifying SNPs for all fluorescence based sequencers. Therefore, VarDetect allows users to adjust these values to fit their experiments. Figure 6 depicts the base-call parameters used in VarDetect. The highest signal is assigned as the primary peak, and the lower signal is assigned as the secondary peak. Signal content below the noise level is ignored. The heterozygosity level in this setting roughly estimates the typically observed nucleotide mixture ratio for biallelic DNA. Many base-calling programs misinterpret heterozygous nucleotides as "N", or calling error, resulting in false negative SNP detection. VarDetect can accommodate chromatogram input data from mixed DNA samples, which is discussed in detail below.

Base-calling
Decaying signals from adjacent bases (a frequent reading artifact) can lead to base miscalling. Consider the scenarios in Figure 7. At the base-call position in Figure 7a, both C and T nucleotide signals are significant in the context of the primary peak, indicative of the true heterozygosity. On the other hand, if part of the primary peak (T) from the adjacent base overlaps with the primary peak at the basecall position (C) (Figure 7b), the base-calling program may misinterpret this as a heterozygous C/T base instead of a homozygous C base, with T disregarded as background noise. Getting accurate information from input traces is thus crucial for SNP detection. The following section discusses a novel base re-sampling and calling technique, which is an improvement over current base-calling programs.

Re-sampling and calling
In this step, we used the Partitioning and Re-sampling (PnR) technique to improve the base-calling procedure.

The main concept of PnR is the detection of secondary
Base-call parameter setting in VarDetect Figure 6 Base-call parameter setting in VarDetect. The highest signal is determined as its primary peak, the lower signal is determined as the secondary peak. The signal contents below the noise level are ignored. The heterozygosity level in this setting roughly estimates nucleotide mixture ratio when dealing with pooled DNA.
Computer representation (array) of chromatogram traces Figure 5 Computer representation (array) of chromatogram traces.
peak turning point patterns. VarDetect divides the chromatogram peak at the base-call position into four equal parts (partitions) and focuses on the two central partitions (shaded boxes in Figure 8a). The two vectors and are created by connecting the points that the curve segment of the secondary peak crosses over the two partitions ( Figure  8b). Let ⊥ be a perpendicular vector of u by rotating it 90° counter-clockwise (Figure 8c). Then the secondary peak curve has a turning point if the dot product of ⊥ · produces a negative value. In other words, if the angle θ between ⊥ and is obtuse, this secondary peak can be interpreted as a heterozygous peak pattern (Figure 8d). In the case of artifactual secondary peak patterns resulting from primary peak overlap from an adjacent base (described above in Figure 7b), the angle between ⊥ and is acute and PnR would determine this secondary peak amplitude as a decaying event and ignore this signal (Figure 9).

Alignment of input sequences to the reference sequence
After re-sampling and base-calling, the input sequences are then aligned against the reference sequence using a local alignment method. There are two steps in this alignment process: 1) pre-alignment and 2) alignment enhancement.

Pre-alignment
Since this tool uses the direct method to search for SNPs, alignment of input sequences to a reference sequence is required. The reference sequence in FASTA format can be obtained from the NCBI public database. VarDetect simplifies the pre-alignment task by linearly searching a local match of m contiguous bases greater than or equal to p percent. In other words, each individual sequence is aligned to the reference sequence by sliding a window of m adjacent bases (W m ) along the reference sequence until a match of p percent or greater is found. Since the noisy parts of input chromatograms can be filltered out using the intensity peak ratio concept explained previously, the candidate window W m is selected from good intensity areas of the input chromatogram. From this observation, we investigate the vicinity peak ratio such that [k] ≥ 90%, starting from the position i ≥ (0.2 × N), where N is the total number of peaks (or 20% of the total number of bases).
The 90% quality value is used to guarantee that the selected regions are readable and good enough for automatic SNP detection. The W m window is formed by inclusively extending the next m bases from the accepted observed base i. The W m window is chosen for each chromatogram trace based on this selection scheme. Each trace is aligned with the reference sequence using W m as its representative in matching the pattern whose percent similarity is greater than a given value. This algorithm, called "Quick Alignment using Sliding Window", is applied repeatedly to more than one W m window to optimize the alignment. The alignment process is performed on both forward and reverse orientations and applied iteratively to each chromatogram. Effect of signal intensity decay on base-calling. Correct (a) and Incorrect (b) base-calling interpretation due to signal intensity decay.

Alignment enhancement
VarDetect uses a sliding window W m to get several local alignment candidate fragments. To complete the overall alignment task, VarDetect increases the window size of each candidate fragment by extending the base-matching from both ends of the window. The longest matching result is then selected. The algorithm is described in Figure  10   and ⊥ is obtuse, this secondary peak could be interpreted as being heterozygous peak (d).

SNP selection
u v u u u v v u one individual can accelerate and reduce the cost for SNP discovery. However, there is the limitation that different DNA samples will be sequenced with different sequencing reaction efficiencies, owing to variable DNA quality and concentration, and variable affinity of DNA polymerase for different nucleotides [12,14]. Despite this limitation, VarDetect can still accurately calculate allele frequencies and detect SNPs from chromatogram traces derived from pooled DNA samples.
When DNA samples from two individuals are pooled, there are five possible biallelic combinations, each producing different chromatogram patterns (Table 3). At a given base i, if both individuals are homozygous (biallelic combinations XXXX and YYYY) or heterozygous (XXYY), the primary and secondary peak intensities would not be different from the data obtained from single samples analyzed separately. On the other hand, if one individual is heterozygous, while the other homozygous (XXXY and XYYY), a weak secondary peak would be observed which must be accounted for by the SNP detection algorithm, or otherwise that SNP would be missed. For this situation, the δ term in VarDetect approximates the amplitude of the secondary peak i, since both Q v and Q o are obtained from the ratio of primary intensity signals normalized to the sum of all signals.

CodeMap
After the re-sampling and calling procedure is completed, accurate base-calling can be done on the chromatogram while noise which may interfere with SNP detection is Illustration of PnR analysis Figure 9 Illustration of PnR analysis. Partitioning (a) and Re-sampling (b) of chromatogram with rising (red) peak. ⊥ is a perpendicular vector of by rotating it 90° counter-clockwise (c). The secondary peak curve has no turning point since the dot product of ⊥ · . produces a positive value (d). Therefore, this peak is interpreted as a homozygous peak. u u u v removed. This "processed" chromatogram is then converted into a numeric code via our CodeMap technique prior to SNP and indel detection in VarDetect. This section describes the generation and interpretation of the code. Figure 11 shows examples of a homozygous chromatogram trace (Figure 11a), a chromatogram with one C/T SNP at the fifth base (Figure 11b), and a chromatogram with T-insertion at the first base (Figure 11c Illustration of traces with indels and their CodeMap analysis In addition to identifying nucleotide substitutions, Var-Detect also automatically detects indels through CodeMap (Figure 11c and its corresponding numeric arrays). When one base (T) is inserted, the following bases are shifted by one frame to the right. Such indels cause misinterpretation errors in most base-calling approaches.
To overcome this problem, the indel chromatograms have to be manually edited by skilled operators.
In the CodeMap view, the code sequence which one would expect to often see from this phenomenon is 2(1/ 2), the code 2 following with either 1 or 2. With an observation on any list in base array (Θ[N]), CodeMap generates this insertion pattern by isolating the correlation between the shifted bases with each possible nucleotide type that is eventually identified by VarDetect. For the adenine array (Θ[A]) in Figure 11c, there are two patterns of "2 1" or "2 2", with four counts, "2 1", "2 1", "2 2", and "2 1". Pattern "2 1" is derived from an observed base A (represented by code 2) being shifted by one position whose base content is not an adenine; hence, it appears in the next position as a heterozygous peak (represented by code 1). Pattern "2 2" is derived from an observed base A (represented by code 2), which is shifted by one position whose content is an adenine. Thus, the next position becomes homozygous A (represented by code 2). These patterns also occur in the cytosine (Θ[C]), guanine (Θ[G]), and thymine (Θ[T]) arrays. For one base deletion, the pattern to be detected is reversed to (1/2)2, the code 1 or 2 following with code 2.
The CodeMap principle may be applied to detect more than one base indel. One would expect to see the pattern 2?[k](1/2) code 2 following with k number of 0, 1 or 2 (we use "?[k]" to represent them) and then finishing with either code 1 or 2. For example, when bases CG are inserted in front of the sequence, the pattern 2? [1](1/2) can be observed. In this data set, other patterns are also possible: 2(1/2) (one base insertion), (1/2)2 (one base deletion), 2? [2](1/2) (three bases insertion) and others (Table 4). In general, N-base indel patterns are generated by a N -1 combination of 0, 1, and 2 codes. Therefore, there are multiple patterns possible in each base array. VarDetect simply counts the number of occurrences of each pattern and uses the patterns with the highest frequencies to determine indel types. Table 4 illustrates the Illustration of VNTR with ATG deletion and its CodeMap analysis Figure 12 Illustration of VNTR with ATG deletion and its CodeMap analysis. CodeMap converts chromatogram of trinucleotide repeats (r 1 , r 2 , r 3 ) to the corresponding numeric arrays (2 0 0 2 0 0) on the right (a). When a set of trinucleotide repeats is deleted, CodeMap reveals specific numeric patterns (underlined) on the right (b), which match with (1/2)? [2]2 pattern shown in Table  4.
pattern counting from Figure 11c. Here the dominant pattern with 13 occurrences is 2(1/2). This CodeMap technique along with pattern detection can be used to efficiently detect SNPs and indels automatically.
Furthermore, for simple indels found in variable number of tandem repeats (VNTR), where the repeat unit composed of more than one nucleotide can increase or decrease in number within a tandem repeat array in the sample sequence, we can theoretically apply Codemap to detect such a case ( Figure 12).