|
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 April 30, 2020
|
|
jpayne@69
|
7
|
|
jpayne@69
|
8 Description: Calculates per-scaffold or per-base coverage information from an unsorted sam or bam file.
|
|
jpayne@69
|
9 Supports SAM/BAM format for reads and FASTA for reference.
|
|
jpayne@69
|
10 Sorting is not needed, so output may be streamed directly from a mapping program.
|
|
jpayne@69
|
11 Requires a minimum of 1 bit per reference base plus 100 bytes per scaffold (even if no reference is specified).
|
|
jpayne@69
|
12 If per-base coverage is needed (including for stdev and median), at least 4 bytes per base is needed.
|
|
jpayne@69
|
13
|
|
jpayne@69
|
14 Usage: pileup.sh in=<input> out=<output>
|
|
jpayne@69
|
15
|
|
jpayne@69
|
16 Input Parameters:
|
|
jpayne@69
|
17 in=<file> The input sam file; this is the only required parameter.
|
|
jpayne@69
|
18 ref=<file> Scans a reference fasta for per-scaffold GC counts, not otherwise needed.
|
|
jpayne@69
|
19 fastaorf=<file> An optional fasta file with ORF header information in PRODIGAL's output format. Must also specify 'outorf'.
|
|
jpayne@69
|
20 unpigz=t Decompress with pigz for faster decompression.
|
|
jpayne@69
|
21 addfromref=t Allow ref scaffolds not present in sam header to be added from the reference.
|
|
jpayne@69
|
22 addfromreads=f Allow ref scaffolds not present in sam header to be added from the reads.
|
|
jpayne@69
|
23 Note that in this case the ref scaffold lengths will be inaccurate.
|
|
jpayne@69
|
24
|
|
jpayne@69
|
25 Output Parameters:
|
|
jpayne@69
|
26 out=<file> (covstats) Per-scaffold coverage info.
|
|
jpayne@69
|
27 rpkm=<file> Per-scaffold RPKM/FPKM counts.
|
|
jpayne@69
|
28 twocolumn=f Change to true to print only ID and Avg_fold instead of all 6 columns.
|
|
jpayne@69
|
29 countgc=t Enable/disable counting of read GC content.
|
|
jpayne@69
|
30 outorf=<file> Per-orf coverage info to this file (only if 'fastaorf' is specified).
|
|
jpayne@69
|
31 outsam=<file> Print the input sam stream to this file (or stdout). Useful for piping data.
|
|
jpayne@69
|
32 hist=<file> Histogram of # occurrences of each depth level.
|
|
jpayne@69
|
33 basecov=<file> Coverage per base location.
|
|
jpayne@69
|
34 bincov=<file> Binned coverage per location (one line per X bases).
|
|
jpayne@69
|
35 binsize=1000 Binsize for binned coverage output.
|
|
jpayne@69
|
36 keepshortbins=t (ksb) Keep residual bins shorter than binsize.
|
|
jpayne@69
|
37 normcov=<file> Normalized coverage by normalized location (X lines per scaffold).
|
|
jpayne@69
|
38 normcovo=<file> Overall normalized coverage by normalized location (X lines for the entire assembly).
|
|
jpayne@69
|
39 normb=-1 If positive, use a fixed number of bins per scaffold; affects 'normcov' and 'normcovo'.
|
|
jpayne@69
|
40 normc=f Normalize coverage to fraction of max per scaffold; affects 'normcov' and 'normcovo'.
|
|
jpayne@69
|
41 delta=f Only print base coverage lines when the coverage differs from the previous base.
|
|
jpayne@69
|
42 nzo=f Only print scaffolds with nonzero coverage.
|
|
jpayne@69
|
43 concise=f Write 'basecov' in a more concise format.
|
|
jpayne@69
|
44 header=t (hdr) Include headers in output files.
|
|
jpayne@69
|
45 headerpound=t (#) Prepend header lines with '#' symbol.
|
|
jpayne@69
|
46 stdev=t Calculate coverage standard deviation.
|
|
jpayne@69
|
47 covminscaf=0 (minscaf) Don't print coverage for scaffolds shorter than this.
|
|
jpayne@69
|
48 covwindow=0 Calculate how many bases are in windows of this size with
|
|
jpayne@69
|
49 low average coverage. Produces an extra stats column.
|
|
jpayne@69
|
50 covwindowavg=5 Average coverage below this will be classified as low.
|
|
jpayne@69
|
51 k=0 If positive, calculate kmer coverage statstics for this kmer length.
|
|
jpayne@69
|
52 keyvalue=f Output statistics to screen as key=value pairs.
|
|
jpayne@69
|
53 mincov=1 When calculating percent covered, ignore bases under this depth.
|
|
jpayne@69
|
54
|
|
jpayne@69
|
55 Processing Parameters:
|
|
jpayne@69
|
56 strandedcov=f Track coverage for plus and minus strand independently.
|
|
jpayne@69
|
57 startcov=f Only track start positions of reads.
|
|
jpayne@69
|
58 stopcov=f Only track stop positions of reads.
|
|
jpayne@69
|
59 secondary=t Use secondary alignments, if present.
|
|
jpayne@69
|
60 softclip=f Include soft-clipped bases in coverage.
|
|
jpayne@69
|
61 minmapq=0 (minq) Ignore alignments with mapq below this.
|
|
jpayne@69
|
62 physical=f (physcov) Calculate physical coverage for paired reads. This includes the unsequenced bases.
|
|
jpayne@69
|
63 tlen=t Track physical coverage from the tlen field rather than recalculating it.
|
|
jpayne@69
|
64 arrays=auto Set to t/f to manually force the use of coverage arrays. Arrays and bitsets are mutually exclusive.
|
|
jpayne@69
|
65 bitsets=auto Set to t/f to manually force the use of coverage bitsets.
|
|
jpayne@69
|
66 32bit=f Set to true if you need per-base coverage over 64k; does not affect per-scaffold coverage precision.
|
|
jpayne@69
|
67 This option will double RAM usage (when calculating per-base coverage).
|
|
jpayne@69
|
68 delcoverage=t (delcov) Count bases covered by deletions or introns as covered.
|
|
jpayne@69
|
69 True is faster than false.
|
|
jpayne@69
|
70 dupecoverage=t (dupes) Include reads flagged as duplicates in coverage.
|
|
jpayne@69
|
71 samstreamer=t (ss) Load reads multithreaded to increase speed.
|
|
jpayne@69
|
72
|
|
jpayne@69
|
73 Trimming Parameters:
|
|
jpayne@69
|
74 ** NOTE: These are applied before adding coverage, to allow mimicking **
|
|
jpayne@69
|
75 ** tools like CallVariants, which uses 'qtrim=r trimq=10 border=5' **
|
|
jpayne@69
|
76 qtrim=f Quality-trim. May be set to:
|
|
jpayne@69
|
77 f (false): Don't quality-trim.
|
|
jpayne@69
|
78 r (right): Trim right (3') end only.
|
|
jpayne@69
|
79 l (left): Trim right (5') end only.
|
|
jpayne@69
|
80 rl (both): Trim both ends.
|
|
jpayne@69
|
81 trimq=-1 If positive, quality-trim to this threshold.
|
|
jpayne@69
|
82 border=0 Ignore this many bases on the left and right end.
|
|
jpayne@69
|
83
|
|
jpayne@69
|
84 Output format (tab-delimited):
|
|
jpayne@69
|
85 ID, Avg_fold, Length, Ref_GC, Covered_percent, Covered_bases, Plus_reads, Minus_reads, Read_GC, Median_fold, Std_Dev
|
|
jpayne@69
|
86
|
|
jpayne@69
|
87 ID: Scaffold ID
|
|
jpayne@69
|
88 Length: Scaffold length
|
|
jpayne@69
|
89 Ref_GC: GC ratio of reference
|
|
jpayne@69
|
90 Avg_fold: Average fold coverage of this scaffold
|
|
jpayne@69
|
91 Covered_percent: Percent of scaffold with any coverage (only if arrays or bitsets are used)
|
|
jpayne@69
|
92 Covered_bases: Number of bases with any coverage (only if arrays or bitsets are used)
|
|
jpayne@69
|
93 Plus_reads: Number of reads mapped to plus strand
|
|
jpayne@69
|
94 Minus_reads: Number of reads mapped to minus strand
|
|
jpayne@69
|
95 Read_GC: Average GC ratio of reads mapped to this scaffold
|
|
jpayne@69
|
96 Median_fold: Median fold coverage of this scaffold (only if arrays are used)
|
|
jpayne@69
|
97 Std_Dev: Standard deviation of coverage (only if arrays are used)
|
|
jpayne@69
|
98
|
|
jpayne@69
|
99 Java Parameters:
|
|
jpayne@69
|
100
|
|
jpayne@69
|
101 -Xmx This will set Java's memory usage, overriding
|
|
jpayne@69
|
102 autodetection. -Xmx20g will
|
|
jpayne@69
|
103 specify 20 gigs of RAM, and -Xmx200m will specify 200 megs.
|
|
jpayne@69
|
104 The max is typically 85% of physical memory.
|
|
jpayne@69
|
105 -eoom This flag will cause the process to exit if an out-of-memory
|
|
jpayne@69
|
106 exception occurs. Requires Java 8u92+.
|
|
jpayne@69
|
107 -da Disable assertions.
|
|
jpayne@69
|
108
|
|
jpayne@69
|
109 Please contact Brian Bushnell at bbushnell@lbl.gov if you encounter any problems.
|
|
jpayne@69
|
110 "
|
|
jpayne@69
|
111 }
|
|
jpayne@69
|
112
|
|
jpayne@69
|
113 #This block allows symlinked shellscripts to correctly set classpath.
|
|
jpayne@69
|
114 pushd . > /dev/null
|
|
jpayne@69
|
115 DIR="${BASH_SOURCE[0]}"
|
|
jpayne@69
|
116 while [ -h "$DIR" ]; do
|
|
jpayne@69
|
117 cd "$(dirname "$DIR")"
|
|
jpayne@69
|
118 DIR="$(readlink "$(basename "$DIR")")"
|
|
jpayne@69
|
119 done
|
|
jpayne@69
|
120 cd "$(dirname "$DIR")"
|
|
jpayne@69
|
121 DIR="$(pwd)/"
|
|
jpayne@69
|
122 popd > /dev/null
|
|
jpayne@69
|
123
|
|
jpayne@69
|
124 #DIR="$( cd "$( dirname "${BASH_SOURCE[0]}" )" && pwd )/"
|
|
jpayne@69
|
125 CP="$DIR""current/"
|
|
jpayne@69
|
126
|
|
jpayne@69
|
127 z="-Xmx1g"
|
|
jpayne@69
|
128 z2="-Xms1g"
|
|
jpayne@69
|
129 set=0
|
|
jpayne@69
|
130
|
|
jpayne@69
|
131 if [ -z "$1" ] || [[ $1 == -h ]] || [[ $1 == --help ]]; then
|
|
jpayne@69
|
132 usage
|
|
jpayne@69
|
133 exit
|
|
jpayne@69
|
134 fi
|
|
jpayne@69
|
135
|
|
jpayne@69
|
136 calcXmx () {
|
|
jpayne@69
|
137 source "$DIR""/calcmem.sh"
|
|
jpayne@69
|
138 setEnvironment
|
|
jpayne@69
|
139 parseXmx "$@"
|
|
jpayne@69
|
140 if [[ $set == 1 ]]; then
|
|
jpayne@69
|
141 return
|
|
jpayne@69
|
142 fi
|
|
jpayne@69
|
143 freeRam 3200m 84
|
|
jpayne@69
|
144 z="-Xmx${RAM}m"
|
|
jpayne@69
|
145 z2="-Xms${RAM}m"
|
|
jpayne@69
|
146 }
|
|
jpayne@69
|
147 calcXmx "$@"
|
|
jpayne@69
|
148
|
|
jpayne@69
|
149 pileup() {
|
|
jpayne@69
|
150 local CMD="java $EA $EOOM $z -cp $CP jgi.CoveragePileup $@"
|
|
jpayne@69
|
151 echo $CMD >&2
|
|
jpayne@69
|
152 eval $CMD
|
|
jpayne@69
|
153 }
|
|
jpayne@69
|
154
|
|
jpayne@69
|
155 pileup "$@"
|
