jpayne@69: #!/bin/bash jpayne@69: jpayne@69: usage(){ jpayne@69: echo " jpayne@69: Written by Brian Bushnell jpayne@69: Last modified September 20, 2022 jpayne@69: jpayne@69: Description: Error corrects reads and/or filters by depth, storing jpayne@69: kmer counts in a count-min sketch (a Bloom filter variant). jpayne@69: This uses a fixed amount of memory. The error-correction algorithm is taken jpayne@69: from Tadpole; with plenty of memory, the behavior is almost identical to jpayne@69: Tadpole. As the number of unique kmers in a dataset increases, the accuracy jpayne@69: decreases such that it will make fewer corrections. It is still capable jpayne@69: of making useful corrections far past the point where Tadpole would crash jpayne@69: by running out of memory, even with the prefilter flag. But if there is jpayne@69: sufficient memory to use Tadpole, then Tadpole is more desirable. jpayne@69: jpayne@69: Because accuracy declines with an increasing number of unique kmers, it can jpayne@69: be useful with very large datasets to run this in 2 passes, with the first jpayne@69: pass for filtering only using a 2-bit filter with the flags tossjunk=t and jpayne@69: ecc=f (and possibly mincount=2 and hcf=0.4), and the second pass using a jpayne@69: 4-bit filter for the actual error correction. jpayne@69: jpayne@69: Usage: bbcms.sh in= out= outb= jpayne@69: jpayne@69: Example of use in error correction: jpayne@69: bbcms.sh in=reads.fq out=ecc.fq bits=4 hashes=3 k=31 merge jpayne@69: jpayne@69: Example of use in depth filtering: jpayne@69: bbcms.sh in=reads.fq out=high.fq outb=low.fq k=31 mincount=2 ecc=f hcf=0.4 jpayne@69: jpayne@69: Error correction and depth filtering can be done simultaneously. jpayne@69: jpayne@69: File parameters: jpayne@69: in= Primary input, or read 1 input. jpayne@69: in2= Read 2 input if reads are in two files. jpayne@69: out= Primary read output. jpayne@69: out2= Read 2 output if reads are in two files. jpayne@69: outb= (outbad/outlow) Output for reads failing mincount. jpayne@69: outb2= (outbad2/outlow2) Read 2 output if reads are in two files. jpayne@69: extra= Additional comma-delimited files for generating kmer counts. jpayne@69: ref= If ref is set, then only files in the ref list will be used jpayne@69: for kmer counts, and the input files will NOT be used for jpayne@69: counts; they will just be filtered or corrected. jpayne@69: overwrite=t (ow) Set to false to force the program to abort rather than jpayne@69: overwrite an existing file. jpayne@69: jpayne@69: Hashing parameters: jpayne@69: k=31 Kmer length, currently 1-31. jpayne@69: hashes=3 Number of hashes per kmer. Higher generally reduces jpayne@69: false positives at the expense of speed; rapidly jpayne@69: diminishing returns above 4. jpayne@69: ksmall= Optional sub-kmer length; setting to slightly lower than k jpayne@69: can improve memory efficiency by reducing the number of hashes jpayne@69: needed. e.g. 'k=31 ksmall=29 hashes=2' has better speed and jpayne@69: accuracy than 'k=31 hashes=3' when the filter is very full. jpayne@69: minprob=0.5 Ignore kmers with probability of being correct below this. jpayne@69: memmult=1.0 Fraction of free memory to use for Bloom filter. 1.0 should jpayne@69: generally work; if the program crashes with an out of memory jpayne@69: error, set this lower. You may be able to increase accuracy jpayne@69: by setting it slightly higher. jpayne@69: cells= Option to set the number of cells manually. By default this jpayne@69: will be autoset to use all available memory. The only reason jpayne@69: to set this is to ensure deterministic output. jpayne@69: seed=0 This will change the hash function used. Useful if running jpayne@69: iteratively with a very full table. -1 uses a random seed. jpayne@69: symmetricwrite=t (sw) Increases counting accuracy for a slight speed penalty. jpayne@69: Could be slow on very low-complexity sequence. jpayne@69: jpayne@69: Depth filtering parameters: jpayne@69: mincount=0 If positive, reads with kmer counts below mincount will jpayne@69: be discarded (sent to outb). jpayne@69: hcf=1.0 (highcountfraction) Fraction of kmers that must be at least jpayne@69: mincount to pass. jpayne@69: requireboth=t Require both reads in a pair to pass in order to go to out. jpayne@69: When true, if either read has a count below mincount, both jpayne@69: reads in the pair will go to outb. When false, reads will jpayne@69: only go to outb if both fail. jpayne@69: tossjunk=f Remove reads or pairs with outermost kmer depth below 2. jpayne@69: (Suggested params for huge metagenomes: mincount=2 hcf=0.4 tossjunk=t) jpayne@69: jpayne@69: Error correction parameters: jpayne@69: ecc=t Perform error correction. jpayne@69: bits= Bits used to store kmer counts; max count is 2^bits-1. jpayne@69: Supports 2, 4, 8, 16, or 32. 16 is best for high-depth data; jpayne@69: 2 or 4 are for huge, low-depth metagenomes that saturate the jpayne@69: bloom filter otherwise. Generally 4 bits is recommended for jpayne@69: error-correction and 2 bits is recommended for filtering only. jpayne@69: ecco=f Error-correct paired reads by overlap prior to kmer-counting. jpayne@69: merge=t Merge paired reads by overlap prior to kmer-counting, and jpayne@69: again prior to correction. Output will still be unmerged. jpayne@69: smooth=3 Remove spikes from kmer counts due to hash collisions. jpayne@69: The number is the max width of peaks to be smoothed; range is jpayne@69: 0-3 (3 is most aggressive; 0 disables smoothing). jpayne@69: This also affects tossjunk. jpayne@69: jpayne@69: jpayne@69: Java Parameters: jpayne@69: -Xmx This will set Java's memory usage, overriding autodetection. jpayne@69: -Xmx20g will specify 20 gigs of RAM, and -Xmx200m will jpayne@69: specify 200 megs. The max is typically 85% of physical memory. jpayne@69: -eoom This flag will cause the process to exit if an out-of-memory jpayne@69: exception occurs. Requires Java 8u92+. jpayne@69: -da Disable assertions. jpayne@69: jpayne@69: Please contact Brian Bushnell at bbushnell@lbl.gov if you encounter any problems. jpayne@69: " jpayne@69: } jpayne@69: jpayne@69: #This block allows symlinked shellscripts to correctly set classpath. jpayne@69: pushd . > /dev/null jpayne@69: DIR="${BASH_SOURCE[0]}" jpayne@69: while [ -h "$DIR" ]; do jpayne@69: cd "$(dirname "$DIR")" jpayne@69: DIR="$(readlink "$(basename "$DIR")")" jpayne@69: done jpayne@69: cd "$(dirname "$DIR")" jpayne@69: DIR="$(pwd)/" jpayne@69: popd > /dev/null jpayne@69: jpayne@69: #DIR="$( cd "$( dirname "${BASH_SOURCE[0]}" )" && pwd )/" jpayne@69: CP="$DIR""current/" jpayne@69: jpayne@69: z="-Xmx4g" jpayne@69: z2="-Xms4g" jpayne@69: set=0 jpayne@69: jpayne@69: if [ -z "$1" ] || [[ $1 == -h ]] || [[ $1 == --help ]]; then jpayne@69: usage jpayne@69: exit jpayne@69: fi jpayne@69: jpayne@69: calcXmx () { jpayne@69: source "$DIR""/calcmem.sh" jpayne@69: setEnvironment jpayne@69: parseXmx "$@" jpayne@69: if [[ $set == 1 ]]; then jpayne@69: return jpayne@69: fi jpayne@69: freeRam 4000m 84 jpayne@69: z="-Xmx${RAM}m" jpayne@69: z2="-Xms${RAM}m" jpayne@69: } jpayne@69: calcXmx "$@" jpayne@69: jpayne@69: bloomfilter() { jpayne@69: local CMD="java $EA $EOOM $z $z2 -cp $CP bloom.BloomFilterCorrectorWrapper $@" jpayne@69: echo $CMD >&2 jpayne@69: eval $CMD jpayne@69: } jpayne@69: jpayne@69: bloomfilter "$@"