|
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 October 14, 2020
|
|
jpayne@69
|
7
|
|
jpayne@69
|
8 Description: Counts the number of unique kmers in a file.
|
|
jpayne@69
|
9 Generates a kmer frequency histogram and genome size estimate (in peaks output),
|
|
jpayne@69
|
10 and prints a file containing all kmers and their counts.
|
|
jpayne@69
|
11 Supports K=1 to infinity, though not all values are allowed.
|
|
jpayne@69
|
12 SEE ALSO: bbnorm.sh/khist.sh, loglog.sh, and kmercountmulti.sh.
|
|
jpayne@69
|
13
|
|
jpayne@69
|
14 Usage: kmercountexact.sh in=<file> khist=<file> peaks=<file>
|
|
jpayne@69
|
15
|
|
jpayne@69
|
16 Input may be fasta or fastq, compressed or uncompressed.
|
|
jpayne@69
|
17 Output may be stdout or a file. out, khist, and peaks are optional.
|
|
jpayne@69
|
18
|
|
jpayne@69
|
19
|
|
jpayne@69
|
20 Input parameters:
|
|
jpayne@69
|
21 in=<file> Primary input file.
|
|
jpayne@69
|
22 in2=<file> Second input file for paired reads.
|
|
jpayne@69
|
23 amino=f Run in amino acid mode.
|
|
jpayne@69
|
24
|
|
jpayne@69
|
25 Output parameters:
|
|
jpayne@69
|
26 out=<file> Print kmers and their counts. This is produces a
|
|
jpayne@69
|
27 huge file, so skip it if you only need the histogram.
|
|
jpayne@69
|
28 fastadump=t Print kmers and counts as fasta versus 2-column tsv.
|
|
jpayne@69
|
29 mincount=1 Only print kmers with at least this depth.
|
|
jpayne@69
|
30 reads=-1 Only process this number of reads, then quit (-1 means all).
|
|
jpayne@69
|
31 dumpthreads=-1 Use this number of threads for dumping kmers (-1 means auto).
|
|
jpayne@69
|
32
|
|
jpayne@69
|
33 Hashing parameters:
|
|
jpayne@69
|
34 k=31 Kmer length (1-31 is fastest).
|
|
jpayne@69
|
35 prealloc=t Pre-allocate memory rather than dynamically growing; faster and more memory-efficient. A float fraction (0-1) may be specified, default 1.
|
|
jpayne@69
|
36 prefilter=0 If set to a positive integer, use a countmin sketch to ignore kmers with depth of that value or lower.
|
|
jpayne@69
|
37 prehashes=2 Number of hashes for prefilter.
|
|
jpayne@69
|
38 prefiltersize=0.2 Fraction of memory to use for prefilter.
|
|
jpayne@69
|
39 minq=6 Ignore kmers containing bases with quality below this. (TODO)
|
|
jpayne@69
|
40 minprob=0.0 Ignore kmers with overall probability of correctness below this.
|
|
jpayne@69
|
41 threads=X Spawn X hashing threads (default is number of logical processors).
|
|
jpayne@69
|
42 onepass=f If true, prefilter will be generated in same pass as kmer counts. Much faster but counts will be lower, by up to prefilter's depth limit.
|
|
jpayne@69
|
43 rcomp=t Store and count each kmer together and its reverse-complement.
|
|
jpayne@69
|
44
|
|
jpayne@69
|
45 Histogram parameters:
|
|
jpayne@69
|
46 khist=<file> Print kmer frequency histogram.
|
|
jpayne@69
|
47 histcolumns=2 2 columns: (depth, count). 3 columns: (depth, rawCount, count).
|
|
jpayne@69
|
48 histmax=100000 Maximum depth to print in histogram output.
|
|
jpayne@69
|
49 histheader=t Set true to print a header line.
|
|
jpayne@69
|
50 nzo=t (nonzeroonly) Only print lines for depths with a nonzero kmer count.
|
|
jpayne@69
|
51 gchist=f Add an extra histogram column with the average GC%.
|
|
jpayne@69
|
52
|
|
jpayne@69
|
53 Intersection parameters:
|
|
jpayne@69
|
54 ref=<file> An input assembly of the input reads.
|
|
jpayne@69
|
55 intersection=<file> Output file for a 2-D histogram of read and ref kmer counts.
|
|
jpayne@69
|
56 refmax=6 Maximum reference kmer depth to track; read depth is controlled by 'histmax'.
|
|
jpayne@69
|
57
|
|
jpayne@69
|
58 Smoothing parameters:
|
|
jpayne@69
|
59 smoothkhist=f Smooth the output kmer histogram.
|
|
jpayne@69
|
60 smoothpeaks=t Smooth the kmer histogram for peak-calling, but does not affect the output histogram.
|
|
jpayne@69
|
61 smoothradius=1 Initial radius of progressive smoothing function.
|
|
jpayne@69
|
62 maxradius=10 Maximum radius of progressive smoothing function.
|
|
jpayne@69
|
63 progressivemult=2 Increment radius each time depth increases by this factor.
|
|
jpayne@69
|
64 logscale=t Transform to log-scale prior to peak-calling.
|
|
jpayne@69
|
65 logwidth=0.1 The larger the number, the smoother.
|
|
jpayne@69
|
66
|
|
jpayne@69
|
67 Peak calling parameters:
|
|
jpayne@69
|
68 peaks=<file> Write the peaks to this file. Default is stdout.
|
|
jpayne@69
|
69 Also contains the genome size estimate in bp.
|
|
jpayne@69
|
70 minHeight=2 (h) Ignore peaks shorter than this.
|
|
jpayne@69
|
71 minVolume=5 (v) Ignore peaks with less area than this.
|
|
jpayne@69
|
72 minWidth=3 (w) Ignore peaks narrower than this.
|
|
jpayne@69
|
73 minPeak=2 (minp) Ignore peaks with an X-value below this.
|
|
jpayne@69
|
74 maxPeak=BIG (maxp) Ignore peaks with an X-value above this.
|
|
jpayne@69
|
75 maxPeakCount=12 (maxpc) Print up to this many peaks (prioritizing height).
|
|
jpayne@69
|
76 ploidy=-1 Specify ploidy; otherwise it will be autodetected.
|
|
jpayne@69
|
77
|
|
jpayne@69
|
78 Sketch parameters (for making a MinHashSketch):
|
|
jpayne@69
|
79 sketch=<file> Write a minhash sketch to this file.
|
|
jpayne@69
|
80 sketchlen=10000 Output the top 10000 kmers. Only kmers with at least mincount are included.
|
|
jpayne@69
|
81 sketchname= Name of output sketch.
|
|
jpayne@69
|
82 sketchid= taxID of output sketch.
|
|
jpayne@69
|
83
|
|
jpayne@69
|
84 Quality parameters:
|
|
jpayne@69
|
85 qtrim=f Trim read ends to remove bases with quality below minq.
|
|
jpayne@69
|
86 Values: t (trim both ends), f (neither end), r (right end only), l (left end only).
|
|
jpayne@69
|
87 trimq=4 Trim quality threshold.
|
|
jpayne@69
|
88 minavgquality=0 (maq) Reads with average quality (before trimming) below this will be discarded.
|
|
jpayne@69
|
89
|
|
jpayne@69
|
90 Overlap parameters (for overlapping paired-end reads only):
|
|
jpayne@69
|
91 merge=f Attempt to merge reads before counting kmers.
|
|
jpayne@69
|
92 ecco=f Error correct via overlap, but do not merge reads.
|
|
jpayne@69
|
93
|
|
jpayne@69
|
94 Java Parameters:
|
|
jpayne@69
|
95 -Xmx This will set Java's memory usage, overriding autodetection.
|
|
jpayne@69
|
96 -Xmx20g will specify 20 gigs of RAM, and -Xmx200m will specify 200 megs.
|
|
jpayne@69
|
97 The max is typically 85% of physical memory.
|
|
jpayne@69
|
98 -eoom This flag will cause the process to exit if an
|
|
jpayne@69
|
99 out-of-memory exception occurs. Requires Java 8u92+.
|
|
jpayne@69
|
100 -da Disable assertions.
|
|
jpayne@69
|
101 "
|
|
jpayne@69
|
102 }
|
|
jpayne@69
|
103
|
|
jpayne@69
|
104 #This block allows symlinked shellscripts to correctly set classpath.
|
|
jpayne@69
|
105 pushd . > /dev/null
|
|
jpayne@69
|
106 DIR="${BASH_SOURCE[0]}"
|
|
jpayne@69
|
107 while [ -h "$DIR" ]; do
|
|
jpayne@69
|
108 cd "$(dirname "$DIR")"
|
|
jpayne@69
|
109 DIR="$(readlink "$(basename "$DIR")")"
|
|
jpayne@69
|
110 done
|
|
jpayne@69
|
111 cd "$(dirname "$DIR")"
|
|
jpayne@69
|
112 DIR="$(pwd)/"
|
|
jpayne@69
|
113 popd > /dev/null
|
|
jpayne@69
|
114
|
|
jpayne@69
|
115 #DIR="$( cd "$( dirname "${BASH_SOURCE[0]}" )" && pwd )/"
|
|
jpayne@69
|
116 CP="$DIR""current/"
|
|
jpayne@69
|
117
|
|
jpayne@69
|
118 z="-Xmx1g"
|
|
jpayne@69
|
119 z2="-Xms1g"
|
|
jpayne@69
|
120 set=0
|
|
jpayne@69
|
121
|
|
jpayne@69
|
122 if [ -z "$1" ] || [[ $1 == -h ]] || [[ $1 == --help ]]; then
|
|
jpayne@69
|
123 usage
|
|
jpayne@69
|
124 exit
|
|
jpayne@69
|
125 fi
|
|
jpayne@69
|
126
|
|
jpayne@69
|
127 calcXmx () {
|
|
jpayne@69
|
128 source "$DIR""/calcmem.sh"
|
|
jpayne@69
|
129 setEnvironment
|
|
jpayne@69
|
130 parseXmx "$@"
|
|
jpayne@69
|
131 if [[ $set == 1 ]]; then
|
|
jpayne@69
|
132 return
|
|
jpayne@69
|
133 fi
|
|
jpayne@69
|
134 freeRam 3200m 84
|
|
jpayne@69
|
135 z="-Xmx${RAM}m"
|
|
jpayne@69
|
136 z2="-Xms${RAM}m"
|
|
jpayne@69
|
137 }
|
|
jpayne@69
|
138 calcXmx "$@"
|
|
jpayne@69
|
139
|
|
jpayne@69
|
140 kmercountexact() {
|
|
jpayne@69
|
141 local CMD="java $EA $EOOM $z $z2 -cp $CP jgi.KmerCountExact $@"
|
|
jpayne@69
|
142 echo $CMD >&2
|
|
jpayne@69
|
143 eval $CMD
|
|
jpayne@69
|
144 }
|
|
jpayne@69
|
145
|
|
jpayne@69
|
146 kmercountexact "$@"
|
