|
jpayne@69
|
1 #!/bin/bash
|
|
jpayne@69
|
2
|
|
jpayne@69
|
3 usage(){
|
|
jpayne@69
|
4 echo "
|
|
jpayne@69
|
5 Written by Brian Bushnell
|
|
jpayne@69
|
6 Last modified September 20, 2022
|
|
jpayne@69
|
7
|
|
jpayne@69
|
8 Description: Error corrects reads and/or filters by depth, storing
|
|
jpayne@69
|
9 kmer counts in a count-min sketch (a Bloom filter variant).
|
|
jpayne@69
|
10 This uses a fixed amount of memory. The error-correction algorithm is taken
|
|
jpayne@69
|
11 from Tadpole; with plenty of memory, the behavior is almost identical to
|
|
jpayne@69
|
12 Tadpole. As the number of unique kmers in a dataset increases, the accuracy
|
|
jpayne@69
|
13 decreases such that it will make fewer corrections. It is still capable
|
|
jpayne@69
|
14 of making useful corrections far past the point where Tadpole would crash
|
|
jpayne@69
|
15 by running out of memory, even with the prefilter flag. But if there is
|
|
jpayne@69
|
16 sufficient memory to use Tadpole, then Tadpole is more desirable.
|
|
jpayne@69
|
17
|
|
jpayne@69
|
18 Because accuracy declines with an increasing number of unique kmers, it can
|
|
jpayne@69
|
19 be useful with very large datasets to run this in 2 passes, with the first
|
|
jpayne@69
|
20 pass for filtering only using a 2-bit filter with the flags tossjunk=t and
|
|
jpayne@69
|
21 ecc=f (and possibly mincount=2 and hcf=0.4), and the second pass using a
|
|
jpayne@69
|
22 4-bit filter for the actual error correction.
|
|
jpayne@69
|
23
|
|
jpayne@69
|
24 Usage: bbcms.sh in=<input file> out=<output> outb=<reads failing filters>
|
|
jpayne@69
|
25
|
|
jpayne@69
|
26 Example of use in error correction:
|
|
jpayne@69
|
27 bbcms.sh in=reads.fq out=ecc.fq bits=4 hashes=3 k=31 merge
|
|
jpayne@69
|
28
|
|
jpayne@69
|
29 Example of use in depth filtering:
|
|
jpayne@69
|
30 bbcms.sh in=reads.fq out=high.fq outb=low.fq k=31 mincount=2 ecc=f hcf=0.4
|
|
jpayne@69
|
31
|
|
jpayne@69
|
32 Error correction and depth filtering can be done simultaneously.
|
|
jpayne@69
|
33
|
|
jpayne@69
|
34 File parameters:
|
|
jpayne@69
|
35 in=<file> Primary input, or read 1 input.
|
|
jpayne@69
|
36 in2=<file> Read 2 input if reads are in two files.
|
|
jpayne@69
|
37 out=<file> Primary read output.
|
|
jpayne@69
|
38 out2=<file> Read 2 output if reads are in two files.
|
|
jpayne@69
|
39 outb=<file> (outbad/outlow) Output for reads failing mincount.
|
|
jpayne@69
|
40 outb2=<file> (outbad2/outlow2) Read 2 output if reads are in two files.
|
|
jpayne@69
|
41 extra=<file> Additional comma-delimited files for generating kmer counts.
|
|
jpayne@69
|
42 ref=<file> If ref is set, then only files in the ref list will be used
|
|
jpayne@69
|
43 for kmer counts, and the input files will NOT be used for
|
|
jpayne@69
|
44 counts; they will just be filtered or corrected.
|
|
jpayne@69
|
45 overwrite=t (ow) Set to false to force the program to abort rather than
|
|
jpayne@69
|
46 overwrite an existing file.
|
|
jpayne@69
|
47
|
|
jpayne@69
|
48 Hashing parameters:
|
|
jpayne@69
|
49 k=31 Kmer length, currently 1-31.
|
|
jpayne@69
|
50 hashes=3 Number of hashes per kmer. Higher generally reduces
|
|
jpayne@69
|
51 false positives at the expense of speed; rapidly
|
|
jpayne@69
|
52 diminishing returns above 4.
|
|
jpayne@69
|
53 ksmall= Optional sub-kmer length; setting to slightly lower than k
|
|
jpayne@69
|
54 can improve memory efficiency by reducing the number of hashes
|
|
jpayne@69
|
55 needed. e.g. 'k=31 ksmall=29 hashes=2' has better speed and
|
|
jpayne@69
|
56 accuracy than 'k=31 hashes=3' when the filter is very full.
|
|
jpayne@69
|
57 minprob=0.5 Ignore kmers with probability of being correct below this.
|
|
jpayne@69
|
58 memmult=1.0 Fraction of free memory to use for Bloom filter. 1.0 should
|
|
jpayne@69
|
59 generally work; if the program crashes with an out of memory
|
|
jpayne@69
|
60 error, set this lower. You may be able to increase accuracy
|
|
jpayne@69
|
61 by setting it slightly higher.
|
|
jpayne@69
|
62 cells= Option to set the number of cells manually. By default this
|
|
jpayne@69
|
63 will be autoset to use all available memory. The only reason
|
|
jpayne@69
|
64 to set this is to ensure deterministic output.
|
|
jpayne@69
|
65 seed=0 This will change the hash function used. Useful if running
|
|
jpayne@69
|
66 iteratively with a very full table. -1 uses a random seed.
|
|
jpayne@69
|
67 symmetricwrite=t (sw) Increases counting accuracy for a slight speed penalty.
|
|
jpayne@69
|
68 Could be slow on very low-complexity sequence.
|
|
jpayne@69
|
69
|
|
jpayne@69
|
70 Depth filtering parameters:
|
|
jpayne@69
|
71 mincount=0 If positive, reads with kmer counts below mincount will
|
|
jpayne@69
|
72 be discarded (sent to outb).
|
|
jpayne@69
|
73 hcf=1.0 (highcountfraction) Fraction of kmers that must be at least
|
|
jpayne@69
|
74 mincount to pass.
|
|
jpayne@69
|
75 requireboth=t Require both reads in a pair to pass in order to go to out.
|
|
jpayne@69
|
76 When true, if either read has a count below mincount, both
|
|
jpayne@69
|
77 reads in the pair will go to outb. When false, reads will
|
|
jpayne@69
|
78 only go to outb if both fail.
|
|
jpayne@69
|
79 tossjunk=f Remove reads or pairs with outermost kmer depth below 2.
|
|
jpayne@69
|
80 (Suggested params for huge metagenomes: mincount=2 hcf=0.4 tossjunk=t)
|
|
jpayne@69
|
81
|
|
jpayne@69
|
82 Error correction parameters:
|
|
jpayne@69
|
83 ecc=t Perform error correction.
|
|
jpayne@69
|
84 bits= Bits used to store kmer counts; max count is 2^bits-1.
|
|
jpayne@69
|
85 Supports 2, 4, 8, 16, or 32. 16 is best for high-depth data;
|
|
jpayne@69
|
86 2 or 4 are for huge, low-depth metagenomes that saturate the
|
|
jpayne@69
|
87 bloom filter otherwise. Generally 4 bits is recommended for
|
|
jpayne@69
|
88 error-correction and 2 bits is recommended for filtering only.
|
|
jpayne@69
|
89 ecco=f Error-correct paired reads by overlap prior to kmer-counting.
|
|
jpayne@69
|
90 merge=t Merge paired reads by overlap prior to kmer-counting, and
|
|
jpayne@69
|
91 again prior to correction. Output will still be unmerged.
|
|
jpayne@69
|
92 smooth=3 Remove spikes from kmer counts due to hash collisions.
|
|
jpayne@69
|
93 The number is the max width of peaks to be smoothed; range is
|
|
jpayne@69
|
94 0-3 (3 is most aggressive; 0 disables smoothing).
|
|
jpayne@69
|
95 This also affects tossjunk.
|
|
jpayne@69
|
96
|
|
jpayne@69
|
97
|
|
jpayne@69
|
98 Java Parameters:
|
|
jpayne@69
|
99 -Xmx This will set Java's memory usage, overriding autodetection.
|
|
jpayne@69
|
100 -Xmx20g will specify 20 gigs of RAM, and -Xmx200m will
|
|
jpayne@69
|
101 specify 200 megs. The max is typically 85% of physical memory.
|
|
jpayne@69
|
102 -eoom This flag will cause the process to exit if an out-of-memory
|
|
jpayne@69
|
103 exception occurs. Requires Java 8u92+.
|
|
jpayne@69
|
104 -da Disable assertions.
|
|
jpayne@69
|
105
|
|
jpayne@69
|
106 Please contact Brian Bushnell at bbushnell@lbl.gov if you encounter any problems.
|
|
jpayne@69
|
107 "
|
|
jpayne@69
|
108 }
|
|
jpayne@69
|
109
|
|
jpayne@69
|
110 #This block allows symlinked shellscripts to correctly set classpath.
|
|
jpayne@69
|
111 pushd . > /dev/null
|
|
jpayne@69
|
112 DIR="${BASH_SOURCE[0]}"
|
|
jpayne@69
|
113 while [ -h "$DIR" ]; do
|
|
jpayne@69
|
114 cd "$(dirname "$DIR")"
|
|
jpayne@69
|
115 DIR="$(readlink "$(basename "$DIR")")"
|
|
jpayne@69
|
116 done
|
|
jpayne@69
|
117 cd "$(dirname "$DIR")"
|
|
jpayne@69
|
118 DIR="$(pwd)/"
|
|
jpayne@69
|
119 popd > /dev/null
|
|
jpayne@69
|
120
|
|
jpayne@69
|
121 #DIR="$( cd "$( dirname "${BASH_SOURCE[0]}" )" && pwd )/"
|
|
jpayne@69
|
122 CP="$DIR""current/"
|
|
jpayne@69
|
123
|
|
jpayne@69
|
124 z="-Xmx4g"
|
|
jpayne@69
|
125 z2="-Xms4g"
|
|
jpayne@69
|
126 set=0
|
|
jpayne@69
|
127
|
|
jpayne@69
|
128 if [ -z "$1" ] || [[ $1 == -h ]] || [[ $1 == --help ]]; then
|
|
jpayne@69
|
129 usage
|
|
jpayne@69
|
130 exit
|
|
jpayne@69
|
131 fi
|
|
jpayne@69
|
132
|
|
jpayne@69
|
133 calcXmx () {
|
|
jpayne@69
|
134 source "$DIR""/calcmem.sh"
|
|
jpayne@69
|
135 setEnvironment
|
|
jpayne@69
|
136 parseXmx "$@"
|
|
jpayne@69
|
137 if [[ $set == 1 ]]; then
|
|
jpayne@69
|
138 return
|
|
jpayne@69
|
139 fi
|
|
jpayne@69
|
140 freeRam 4000m 84
|
|
jpayne@69
|
141 z="-Xmx${RAM}m"
|
|
jpayne@69
|
142 z2="-Xms${RAM}m"
|
|
jpayne@69
|
143 }
|
|
jpayne@69
|
144 calcXmx "$@"
|
|
jpayne@69
|
145
|
|
jpayne@69
|
146 bloomfilter() {
|
|
jpayne@69
|
147 local CMD="java $EA $EOOM $z $z2 -cp $CP bloom.BloomFilterCorrectorWrapper $@"
|
|
jpayne@69
|
148 echo $CMD >&2
|
|
jpayne@69
|
149 eval $CMD
|
|
jpayne@69
|
150 }
|
|
jpayne@69
|
151
|
|
jpayne@69
|
152 bloomfilter "$@"
|
