Mercurial > repos > rliterman > csp2

comparison CSP2/CSP2_env/env-d9b9114564458d9d-741b3de822f2aaca6c6caa4325c4afce/opt/bbmap-39.01-1/kmercountexact.sh @ 69:33d812a61356

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