jpayne@69: #!/bin/bash jpayne@69: jpayne@69: usage(){ jpayne@69: echo " jpayne@69: Written by Brian Bushnell jpayne@69: Last modified April 30, 2020 jpayne@69: jpayne@69: Description: Calculates per-scaffold or per-base coverage information from an unsorted sam or bam file. jpayne@69: Supports SAM/BAM format for reads and FASTA for reference. jpayne@69: Sorting is not needed, so output may be streamed directly from a mapping program. jpayne@69: Requires a minimum of 1 bit per reference base plus 100 bytes per scaffold (even if no reference is specified). jpayne@69: If per-base coverage is needed (including for stdev and median), at least 4 bytes per base is needed. jpayne@69: jpayne@69: Usage: pileup.sh in= out= jpayne@69: jpayne@69: Input Parameters: jpayne@69: in= The input sam file; this is the only required parameter. jpayne@69: ref= Scans a reference fasta for per-scaffold GC counts, not otherwise needed. jpayne@69: fastaorf= An optional fasta file with ORF header information in PRODIGAL's output format. Must also specify 'outorf'. jpayne@69: unpigz=t Decompress with pigz for faster decompression. jpayne@69: addfromref=t Allow ref scaffolds not present in sam header to be added from the reference. jpayne@69: addfromreads=f Allow ref scaffolds not present in sam header to be added from the reads. jpayne@69: Note that in this case the ref scaffold lengths will be inaccurate. jpayne@69: jpayne@69: Output Parameters: jpayne@69: out= (covstats) Per-scaffold coverage info. jpayne@69: rpkm= Per-scaffold RPKM/FPKM counts. jpayne@69: twocolumn=f Change to true to print only ID and Avg_fold instead of all 6 columns. jpayne@69: countgc=t Enable/disable counting of read GC content. jpayne@69: outorf= Per-orf coverage info to this file (only if 'fastaorf' is specified). jpayne@69: outsam= Print the input sam stream to this file (or stdout). Useful for piping data. jpayne@69: hist= Histogram of # occurrences of each depth level. jpayne@69: basecov= Coverage per base location. jpayne@69: bincov= Binned coverage per location (one line per X bases). jpayne@69: binsize=1000 Binsize for binned coverage output. jpayne@69: keepshortbins=t (ksb) Keep residual bins shorter than binsize. jpayne@69: normcov= Normalized coverage by normalized location (X lines per scaffold). jpayne@69: normcovo= Overall normalized coverage by normalized location (X lines for the entire assembly). jpayne@69: normb=-1 If positive, use a fixed number of bins per scaffold; affects 'normcov' and 'normcovo'. jpayne@69: normc=f Normalize coverage to fraction of max per scaffold; affects 'normcov' and 'normcovo'. jpayne@69: delta=f Only print base coverage lines when the coverage differs from the previous base. jpayne@69: nzo=f Only print scaffolds with nonzero coverage. jpayne@69: concise=f Write 'basecov' in a more concise format. jpayne@69: header=t (hdr) Include headers in output files. jpayne@69: headerpound=t (#) Prepend header lines with '#' symbol. jpayne@69: stdev=t Calculate coverage standard deviation. jpayne@69: covminscaf=0 (minscaf) Don't print coverage for scaffolds shorter than this. jpayne@69: covwindow=0 Calculate how many bases are in windows of this size with jpayne@69: low average coverage. Produces an extra stats column. jpayne@69: covwindowavg=5 Average coverage below this will be classified as low. jpayne@69: k=0 If positive, calculate kmer coverage statstics for this kmer length. jpayne@69: keyvalue=f Output statistics to screen as key=value pairs. jpayne@69: mincov=1 When calculating percent covered, ignore bases under this depth. jpayne@69: jpayne@69: Processing Parameters: jpayne@69: strandedcov=f Track coverage for plus and minus strand independently. jpayne@69: startcov=f Only track start positions of reads. jpayne@69: stopcov=f Only track stop positions of reads. jpayne@69: secondary=t Use secondary alignments, if present. jpayne@69: softclip=f Include soft-clipped bases in coverage. jpayne@69: minmapq=0 (minq) Ignore alignments with mapq below this. jpayne@69: physical=f (physcov) Calculate physical coverage for paired reads. This includes the unsequenced bases. jpayne@69: tlen=t Track physical coverage from the tlen field rather than recalculating it. jpayne@69: arrays=auto Set to t/f to manually force the use of coverage arrays. Arrays and bitsets are mutually exclusive. jpayne@69: bitsets=auto Set to t/f to manually force the use of coverage bitsets. jpayne@69: 32bit=f Set to true if you need per-base coverage over 64k; does not affect per-scaffold coverage precision. jpayne@69: This option will double RAM usage (when calculating per-base coverage). jpayne@69: delcoverage=t (delcov) Count bases covered by deletions or introns as covered. jpayne@69: True is faster than false. jpayne@69: dupecoverage=t (dupes) Include reads flagged as duplicates in coverage. jpayne@69: samstreamer=t (ss) Load reads multithreaded to increase speed. jpayne@69: jpayne@69: Trimming Parameters: jpayne@69: ** NOTE: These are applied before adding coverage, to allow mimicking ** jpayne@69: ** tools like CallVariants, which uses 'qtrim=r trimq=10 border=5' ** jpayne@69: qtrim=f Quality-trim. May be set to: jpayne@69: f (false): Don't quality-trim. jpayne@69: r (right): Trim right (3') end only. jpayne@69: l (left): Trim right (5') end only. jpayne@69: rl (both): Trim both ends. jpayne@69: trimq=-1 If positive, quality-trim to this threshold. jpayne@69: border=0 Ignore this many bases on the left and right end. jpayne@69: jpayne@69: Output format (tab-delimited): jpayne@69: ID, Avg_fold, Length, Ref_GC, Covered_percent, Covered_bases, Plus_reads, Minus_reads, Read_GC, Median_fold, Std_Dev jpayne@69: jpayne@69: ID: Scaffold ID jpayne@69: Length: Scaffold length jpayne@69: Ref_GC: GC ratio of reference jpayne@69: Avg_fold: Average fold coverage of this scaffold jpayne@69: Covered_percent: Percent of scaffold with any coverage (only if arrays or bitsets are used) jpayne@69: Covered_bases: Number of bases with any coverage (only if arrays or bitsets are used) jpayne@69: Plus_reads: Number of reads mapped to plus strand jpayne@69: Minus_reads: Number of reads mapped to minus strand jpayne@69: Read_GC: Average GC ratio of reads mapped to this scaffold jpayne@69: Median_fold: Median fold coverage of this scaffold (only if arrays are used) jpayne@69: Std_Dev: Standard deviation of coverage (only if arrays are used) jpayne@69: jpayne@69: Java Parameters: jpayne@69: jpayne@69: -Xmx This will set Java's memory usage, overriding jpayne@69: autodetection. -Xmx20g will jpayne@69: specify 20 gigs of RAM, and -Xmx200m will specify 200 megs. jpayne@69: 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="-Xmx1g" jpayne@69: z2="-Xms1g" 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 3200m 84 jpayne@69: z="-Xmx${RAM}m" jpayne@69: z2="-Xms${RAM}m" jpayne@69: } jpayne@69: calcXmx "$@" jpayne@69: jpayne@69: pileup() { jpayne@69: local CMD="java $EA $EOOM $z -cp $CP jgi.CoveragePileup $@" jpayne@69: echo $CMD >&2 jpayne@69: eval $CMD jpayne@69: } jpayne@69: jpayne@69: pileup "$@"