View Single Post
Old 01-12-2017, 10:14 AM   #45
Brian Bushnell
Super Moderator
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,695

I've described the algorithm in some detail in /bbmap/docs/guides/BBNormGuide.txt. I also wrote this a while back:


This program accepts input files of single or paired reads in fasta or fastq format, correcting substitution-type errors and normalizing the read depth to some desired target before outputting the result. All stages are multithreaded, allowing very high processing speed while still (optionally) maintaining strictly deterministic output.

Phase 1: Gather Kmer Frequencies

An input file of sequence data is read and processed. Each read is translated into a set of all constituent kmers of fixed k (default 31). Each kmer’s count is incremented in a shared table (a count-min sketch) whenever it is seen, so at the end of this phase, the frequencies of all kmers are known.

Phase 2: Correct and Normalize Reads

The input file is read a second time. Each read is again translated into an array of kmers, and each kmer’s count is read from the table. An error-free read is expected to have a relatively smooth profile of kmer counts, so each read is scanned for the presence of adjacent kmers with discrepant counts. For such a pair of adjacent kmers, the one with the high count is considered a “good kmer” (or genomic kmer) and the one with the low count is considered a possible “error kmer”. The single base covered by the error kmer but not the good kmer is considered the suspect “error base”. In addition to absolute cutoffs for counts of good and bad kmers, data with very high coverage is handled with a relative cutoff for the ratio of adjacent kmer counts.

All 4 possible replacements of the error base (A, C, G, T) are considered. For each replacement, the kmer covering the error base is regenerated, and its count read from the table. If exactly one of the four replacements has a count sufficiently high to be considered a good kmer and the other three are sufficiently low to be considered error kmers, then the error base is replaced accordingly, and the error is considered corrected. Otherwise, the error cannot be corrected; any prior corrections are rolled back, and the read is output unchanged.

If normalization is desired, the kmer counts from correction are re-used to determine whether a read should be discarded. If the median count is below a cutoff, the read is discarded as noise. Reads between the lower cutoff and the target depth are all retained. Otherwise, the median is above the target depth and the read is discarded with probability 1-(target/median). For paired reads, the greater median is compared to the cutoff, and the lesser median is compared to the target; the pair is either kept or discarded together. Normalization may be run using multiple passes for greater precision.
Note that I don't recommend BBNorm for error-correction anymore, though, since Tadpole does a much better job (which is possible because it uses exact kmer counts). So I just use BBNorm for normalization and depth partitioning.
Brian Bushnell is offline   Reply With Quote