BBNorm Guide
Written by Brian Bushnell
Last updated December 16, 2015
BBNorm是用来对高覆盖率的reads进行抽样的工具,这个过程可以极大的==加速组装==,==让无法组装的数据变得"可以组装"==,==同时一般也可以极大提高组装的质量==。此外还有:基于深度的-binning,kmer频度直方图,错误矫正,错误标记,基因组大小估计
BBNorm有四个显著的特点:
- BBNorm绝对不会run out of memory,但是如果kmers太多了会使得精度稍稍下降
- Multipass模式可以减少normalized的输出中的平均错误率,standard normalization可以富集包含错误的reads
- 速度超快,并且很易用
- 支持任何长度的kmer值
数据结构:
A Count-Min Sketch (CMS) is also called a "counting Bloom filter". 这是一种只存储值,不存keys并且忽略冲突(collisions)的hash table。为了避免collisions带来的负面影响,值被储存在多个位置,尽可能保障至少有一个不会和其他的冲突。在读取kmer数目时,所有位置都被读取,最小值被使用
- BBNorm 可以使用1, 2, 4, 8, 16, or 32 位每cell的CMSs结构. 当字节越大,最大计数值也可以更(up to 2^bits-1), 但是可供使用的cells更少; 例如,for example, 1GB RAM will accommodate (提供空间)4 billion 2-bit cells, 最大计数值可以为3, or 500 million 16-bit cells, 最大计数值为65535. 在我们数据预期覆盖度是200x时,设置32-bit的单个cell大小是完全没道理的。同时,存储一个kmer的计数值的位置可以被设置,从1到无限大都可以(默认值为3),越多存储位置会越准确(直到hash table变满),但是会更慢. 确定最佳hash数,请参考Bloom filters.
Memory and Capacity:
- BBNorm默认情况下会使用所有的运行内存,使用越多的运行内存结果越准确
- 当处理很大的数据集,却只有很小的运行内存可用时,会得到警告信息:
"Made hash table: hashes = 1 mem = 581.26 MB cells = 152.38M used = 93.540%
Warning: This table is extremely full, which may reduce accuracy. Ideal load is under 60% used. - 因此,可以看出,要获得最好的精度,load尽量小于60%
- 为了更加准确:
- 使用'prefilter'
- 在一个节点上使用更多的运行内存运行
- 质控数据或者对数据进行错误矫正
- 加大minprob的值,以减少虚假的kmers
- ==在实际运行中,即便是loads超过了90%,你还是可以得到比较好的normalization结果,但是histogram和statistics会被关闭==
请不要忽略上面的信息,在使用prefilter会把较少的kmers计数值存储在小的cells上(默认为2-bit),把大的计数值储存在大的cells上 (默认32-bit),这时候memory可以被更加高效的使用。 Prefilter在默认情况下是关闭的,但是当我们很需要==准确性==的时候,或者==tables实在是有点满了==(比如normalization时超过了50%,在error correction时我们希望loads更低),请务必打开这个参数。
- 也可以把primary cells 的值变小些, e.g. "bits=16"
- 或者先执行一次2-bit或4-bit的过滤,把那些非常低覆盖度的reads给去除掉,从而减少总的独特kmers数目。
- 我们也可以调整minprob参数,which ignores kmers with a probability of being error-free (based on quality scores) of below X。
- normalization在tables很满时问题也不大,但是错误矫正和histogram generation会变得没那么可靠
Shellscripts:
BBNorm有三个脚本(bbnorm.sh, ecc.sh, and khist.sh),他们都使用KmerNormalize,但配合不同的参数,在错误矫正和normalization的时候我们还是可以make kmer frequency histogram的。通过调整参数,这三个脚本的功能可以变得一样。
丢弃kmers,精确计数,错误矫正--Dumping Kmers, Exact Counts, and Error Correction:
BBNorm不能丢弃具体的kmers,因为BBNorm只存储了计数值。为了丢弃具体的kmer,可以使用kmercountexact,它可以生成具体kmer的频度分布直方图但,是这个时候不能再分析很大的数据了。由==由于Tadpole使用和kemrcountexact相同的数据结构,在错误矫正方面,tadpole优于BBNorm,因此如果我们有足够的内存,最好使用tadpole进行错误矫正。==
什么时候不要用BBNorm--When Not To Use BBNorm:
- 首先BBNorm不能把低丰度的数据丰度变高,它主要针对数据覆盖度太高或者分布抬不均匀(amplified single-cell, RNA-seq, viruses, metagenomes, etc),主要服务于组装。
- 在需要对样品进行定量的时候,用于mapping的数据是没有进行过矫正的。
- 同时如果数据的错误率太高了,也是不合适的,BBNorm主要用于Illumina和Ion Torrent数据。
- 如果要进行rare variants calling,也最好不要使用这个工具
临时文件管理和管道--Temp Files and Piping:
不可以使用管道
BBNorm needs to read input files multiple times (twice per pass), which means it is unable to accept piped input. In multipass mode, it also needs to generate temp files. The location of temp files can be specified with the "tmpdir" flag; it defaults to the environment variable $TMPDIR, which on Genepool points to local disk when available. Temp files will be cleaned up once BBNorm finishes.
核心数设置--Threads:
可以使用t参数设置核心数
BBNorm is fully multithreaded both when counting kmers and when doing other operations such as normalization. The counting is lock-free, using atomic counters. As a result, it will default to using all available hardware threads; this can be adjusted with the "t" flag.
==Usage Examples==
估计内存需求:
loglog.sh in=reads.fq
这个命令可以估计数据集中的独特kmers,从而估计kmer-counting programs such as BBNorm需要的内存量。这个过程是很快的,并且几乎不耗内存。如果loglog报告10亿kmers,那么有3hashes(每个计数值储存到三个地方)16bits(每个cell的大小)/kmer/8bits(每个字节有8位)/byte1,000,000,000kmers/0.5(占用50%的table)load=12 GB的内存便可以使得table消耗小于50%对reads数目进行normalization:
bbnorm.sh in=reads.fq out=normalized.fq target=100 min=5
这条命令会将覆盖度超过100的reads都下调到100,同时盖度低于5的reads都会被当做是错误从而丢弃
- To error-correct reads:
ecc.sh in=reads.fq out=corrected.fq
或者
bbnorm.sh in=reads.fq out=corrected.fq ecc=t keepall passes=1 bits=16 prefilter
这个可以进行错误矫正但是不丢弃任何reads
bits=16 prefilter
并不是必须的,但是可以通过存储kmer更加高效从而使得矫正更加准确
- 生成kmer分布直方图
khist.sh in=reads.fq khist=khist.txt peaks=peaks.txt
or equivalently
bbnorm.sh in=reads.fq khist=khist.txt peaks=peaks.txt passes=1 prefilter minprob=0 minqual=0 mindepth=0
- 直方图可以表示出在指定深度时,独特kmers的数目。比如,a point at "Depth=10, UniqueKmers=248028"指明248028个kmers的测序深度为10。这个直方图要在对数坐标上画,peak文件包含直方图中peak的位置,以及对基因组大小和几倍体的估计。==这个估计只有在数据是随机断裂DNA,倍型数在4以下(1或者2的时候最准确),并且污染很少时才可靠。如果在kmer直方图中没有明显的peaks,那么peaks文件中的结果是不可靠的。==
bbnorm.sh在这的其他参数 (minprob=0 minqual=0 mindepth=0
) 用于避免低丰度的kmers被丢掉
- 直方图可以表示出在指定深度时,独特kmers的数目。比如,a point at "Depth=10, UniqueKmers=248028"指明248028个kmers的测序深度为10。这个直方图要在对数坐标上画,peak文件包含直方图中peak的位置,以及对基因组大小和几倍体的估计。==这个估计只有在数据是随机断裂DNA,倍型数在4以下(1或者2的时候最准确),并且污染很少时才可靠。如果在kmer直方图中没有明显的peaks,那么peaks文件中的结果是不可靠的。==
- Normalization测序数据,并且进行错误矫正,在矫正前后都生成kmer直方图图
bbnorm.sh in=reads.fq out=normalized.fq target=100 min=5 prefilter ecc khist=khist_before.txt khistout=khist_after.txt
- To make a high-pass or low-pass filter:
bbnorm.sh in=reads.fq out=highpass.fq outt=lowpass.fq passes=1 target=999999999 min=10
这个命令可以把coverage大于10的输出到out中,小于10的输出到outt中
- 按照depth把数据分到3个区间中To split by depth into 3 bins:
bbnorm.sh in=reads.fq outlow=low.fq outmid=mid.fq outhigh=high.fq passes=1 lowbindepth=10 highbindepth=80
10x 以下放到 low.fq; 80x 以上放到 high.fq; 其他的放到 mid.fq。 ==特别的, 对于承兑数据,如果一端的depth低于低阈值,另一端depth高于高阈值,那么它们都去mid数据集==
- 用额外的序列信息辅助对kmer的统计:
bbnorm.sh in=reads.fq out=corrected.fq passes=1 ecc extra=genome.fa,more_reads.fq
这个指令会利用genome.fa和more_reads.fa中的kmer计数信息指导对in=
数据中reads的矫正。这个操作也可以在其他功能中使用,例如normalization。Extra只会被用于提供kmer丰度信息,不会再输出中明显出现
额外补充一点信息
- BBTools has another program functionally related to BBNorm, "kmercountexact.sh". It does NOT use probabilistic data structures, and uses locking rather than atomic counters, and as a result may not scale as well, and will run out of memory on large datasets. However, it is still extremely fast and memory-efficient - using ~15 bytes per kmer (with an optional count-min-sketch prefilter to remove low-count error kmers). It cannot normalize or error-correct, but it can generate the exact kmer count of a dataset as well as the exact kmer frequency histogram (and do rudimentary peak calling for genome size estimation). In practice, when I am interested in kmer frequency histograms, I use KmerCountExact for isolates, and BBNorm for metagenomes.