VarDetect: a nucleotide sequence variation exploratory tool

Following completion of the human genome project, detection and discovery of single nucleotide polymorphisms (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 associated 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–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.

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, VarDetect 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.

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 VarDetect website. Finally, we measured the software precision and recall ratios in order to compare them using F-score [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 discovery. 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 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.

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 base-call 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 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, novoSNP, 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. VarDetect is platform independent since it was implemented in Java, allowing it to run on all major operating systems without recompiling the source codes.

where $Q_o^i$ is the signal intensity of a nucleotide at the i^th position (the base-call position), and is the ratio of the highest signal intensity c ⁱ (b), where b ∈ {A, T, G, C}, to the sum of the signal intensities of adenine c ⁱ (A), thymine c ⁱ (T), cytosine c ⁱ (C), and guanine c ⁱ (G), respectively.

This term is the arithmetic mean of the signal intensities which flank $Q_o^i$ to the left for k bases and to the right for k bases. In other words, it is the summation of observed intensities of k bases toward the left and right of $Q_o^i$ divided by 2k (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:

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.

Base-calling

Base-calling translates the chromatogram traces to one of the four nucleotide characters A, T, G, or C along the sequence. When two or more different bases are detected at a calling position, the International Union of Pure and Applied Chemistry (IUPAC) ambiguous nucleotide codes are assigned to that position. At the end of this step, a corresponding sequence is obtained as the base-calling result. By using BioJava, the data from a chromatogram trace is represented by five arrays: base-call array or B [], and nucleotide content arrays: A [], T [], C [], and G [], for adenine, thymine, cytosine, and guanine, respectively (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.

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 base-call 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 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 $\overrightarrow{u}$ and $\overrightarrow{v}$ are created by connecting the points that the curve segment of the secondary peak crosses over the two partitions (Figure 8b). Let $\overrightarrow{u}$^⊥ 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 $\overrightarrow{u}$^⊥·$\overrightarrow{v}$ produces a negative value. In other words, if the angle θ between $\overrightarrow{u}$^⊥ and $\overrightarrow{v}$ 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 $\overrightarrow{u}$^⊥ and $\overrightarrow{v}$ is acute and PnR would determine this secondary peak amplitude as a decaying event and ignore this signal (Figure 9).

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

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 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). CodeMap individually observes variations associated with A, C, G or T and performs the conversion from the biallelic values at each base position to 0, 1, and 2 according to the following rule. Two prominent values, primary and secondary peak intensities, at each base i^th are symbolized as $b_{\text{pri}}^i$ and $b_{\text{sec}}^i$ respectively.

Four code sequences Θ[A], Θ[C], Θ[G], and Θ[T] collectively termed Θ[N], where N represents either A, C, G or T, can be constructed by:

$\begin{array}{cc}\Theta [N],& \text{base }{b}^i=\{\begin{array}{ll}2\hfill & \text{if }N\text{ is }{b}_{\text{pri}}^i\hfill \\ 1\hfill & \text{if }N\text{ is }{b}_{\sec }^i\hfill \\ 0\hfill & \text{if }N\text{ is not }{b}^i\hfill \end{array}\end{array}$

(3)

Using Equation 3, CodeMap converts the chromatograms in Figures 11 to numeric arrays. The homozygous base is converted into 0 and 2 codes (Figure 11a) while heterozygous base in Figure 11b (Θ[T]) is converted to 1. VarDetect can make use of this code in conjunction with the aforementioned δ value to automatically detect SNPs.

In addition to identifying nucleotide substitutions, VarDetect 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 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.

	Pattern Counts
Pattern	Θ[A]	Θ[C]	Θ[G]	Θ[T]	Total	Comment
2(1/2)	4	4	2	3	13	1 insertion
2?[1](1/2)	1	2	0	0	3	2 insertion
2?[2](1/2)	0	1	0	0	1	3 insertion
(1/2)2	1	2	0	1	4	1 deletion
(1/2)?[1]2	0	1	0	0	1	2 deletion
(1/2)?[2]2	1	1	0	3	5	3 deletion

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).