jpayne@69: #!/bin/bash jpayne@69: jpayne@69: usage(){ jpayne@69: echo " jpayne@69: Written by Brian Bushnell jpayne@69: Last modified March 21, 2019 jpayne@69: jpayne@69: Description: Calculates observed quality scores from mapped sam/bam files. jpayne@69: Generates matrices for use in recalibrating quality scores. By default, jpayne@69: the matrices are written to /ref/qual/ in the current directory. jpayne@69: jpayne@69: If you have multiple sam/bam files demultiplexed from a single sequencing run, jpayne@69: it is recommended to use all of them as input for increased statistical power. jpayne@69: Once the matrices are generated, recalibration can be done on mapped or jpayne@69: unmapped reads; you may get better results by recalibrating the fastq and jpayne@69: remapping the calibrated reads. jpayne@69: jpayne@69: Note! Diploid organisms with a high heterozygousity rate will induce jpayne@69: inaccurate recalibration at the high end of the quality scale unless SNP jpayne@69: locations are masked or variations are called. For example, recalibrating jpayne@69: human reads mapped to an unmasked human reference would generate an jpayne@69: expected maximal Q-score of roughly 30 due to the human 1/1000 SNP rate. jpayne@69: Variations can be ignored by using the callvars flag or providing jpayne@69: a file of variations. jpayne@69: jpayne@69: Usage: jpayne@69: jpayne@69: Step 1. Generate matrices (from mapped sam or bam files): jpayne@69: calctruequality.sh in= path= jpayne@69: jpayne@69: Step 2. Recalibrate reads (any kind of files): jpayne@69: bbduk.sh in= out= recalibrate jpayne@69: jpayne@69: jpayne@69: Parameters (and their defaults) jpayne@69: jpayne@69: Input parameters: jpayne@69: in= Sam file or comma-delimited list of files. Alignments jpayne@69: must use = and X cigar symbols, or have MD tags, or jpayne@69: ref must be specified. jpayne@69: reads=-1 Stop after processing this many reads (if positive). jpayne@69: samstreamer=t (ss) Load reads multithreaded to increase speed. jpayne@69: unpigz=t Use pigz to decompress. jpayne@69: jpayne@69: Output parameters: jpayne@69: overwrite=t (ow) Set to true to allow overwriting of existing files. jpayne@69: path=. Directory to write quality matrices (within /ref subdir). jpayne@69: write=t Write matrices. jpayne@69: showstats=t Print a summary. jpayne@69: pigz=f Use pigz to compress. jpayne@69: jpayne@69: Other parameters: jpayne@69: t=auto Number of worker threads. jpayne@69: passes=2 Recalibration passes, 1 or 2. 2 is slower but gives more jpayne@69: accurate quality scores. jpayne@69: recalqmax=42 Adjust max quality scores tracked. The actual highest jpayne@69: quality score allowed is recalqmax-1. jpayne@69: trackall=f Track all available quality metrics and produce all jpayne@69: matrices, including the ones that are not selected for jpayne@69: quality adjustment. Reduces speed, but allows testing the jpayne@69: effects of different recalibration matrices. jpayne@69: indels=t Include indels in quality calculations. jpayne@69: jpayne@69: Variation calling: jpayne@69: varfile= Use the variants in this var file, instead of calling jpayne@69: variants. The format can be produced by CallVariants. jpayne@69: vcf= Use the variants in this VCF file, instead of jpayne@69: calling variants. jpayne@69: callvars=f Call SNPs, and do not count them as errors. jpayne@69: ploidy=1 Set the organism's ploidy. jpayne@69: ref= Required for variation-calling. jpayne@69: jpayne@69: *** 'Variant-Calling Cutoffs' flags in callvariants.sh are also supported *** jpayne@69: jpayne@69: Selecting matrices: jpayne@69: loadq102= For each recalibration matrix, enable or disable that matrix with t/f. jpayne@69: You can specify pass1 or pass2 like this: loadq102_p1=f loadq102_p2=t. jpayne@69: The default is loadqbp_p1=t loadqbp_p2=t loadqb123_p=t. jpayne@69: clearmatrices=f If true, clear all the existing matrix selections. For example: jpayne@69: 'clearmatrices loadqbp_p1' jpayne@69: This would ignore defaults and select only qbp for the first pass. jpayne@69: jpayne@69: Available matrices: jpayne@69: q102 Quality, leading quality, trailing quality. jpayne@69: qap Quality, average quality, position. jpayne@69: qbp Quality, current base, position. jpayne@69: q10 Quality, leading quality. jpayne@69: q12 Quality, trailing quality. jpayne@69: qb12 Quality, leading base, current base. jpayne@69: qb012 Quality, two leading bases, current base. jpayne@69: qb123 Quality, leading base, current base, trailing base. jpayne@69: qb234 Quality, current base, two trailing bases. jpayne@69: q12b12 Quality, trailing quality, leading base, current base. jpayne@69: qp Quality, position. jpayne@69: q Current quality score only. jpayne@69: jpayne@69: jpayne@69: Java Parameters: jpayne@69: -Xmx This will set Java's memory usage, overriding autodetection. jpayne@69: -Xmx20g will 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 jpayne@69: out-of-memory 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="-Xmx2g" jpayne@69: z2="-Xms2g" 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: calctruequality() { jpayne@69: local CMD="java $EA $EOOM $z $z2 -cp $CP jgi.CalcTrueQuality $@" jpayne@69: echo $CMD >&2 jpayne@69: eval $CMD jpayne@69: } jpayne@69: jpayne@69: calctruequality "$@"