awk: cmd. line:1: fatal: division by zero attempted

Question

[Maryem]

hey, i'm struggling with this error. after i run this script ), i got this error

awk: cmd. line:1: fatal: division by zero attempted
Done.

I think the problem in AWK but i can't find it the output is 3 empty files : depth.feat.gz depth.feat.log depth.feat.stats. If any one can help Thank you

( I'm trying to run this tool in github : https://github.com/AshleyLab/scotch )

# SCRIPT : 
______________________
#!/bin/bash
#Calculates depth of coverage (includes or excludes SC bases 
#depending on whether called with unclipped bam or not from getFeatures.sh)
#and writes, gzipped, to outFeature

bam="$1"

bed="$2"

fastaRef="$3"

gatkJAR="$4"

tmpDir="$5"

mem="$6"

outFeature="$7"

outFeatureStats="$8"


outputBase=$(sed "s|\.gz$||" <<< "$outFeature")

log="${outputBase}.log"

#Use GATK to calculate depth across portion of BAM in BED

java -Djava.io.tmpdir="$tmpDir" -Xmx"$mem"g -jar "$gatkJAR" \

    -T DepthOfCoverage \

    -mbq 0 -mmq 0 \

    --omitIntervalStatistics \

    --omitLocusTable \

    --omitPerSampleStats \

    --includeRefNSites \

    -L "$bed" \

    -I "$bam" \

    -o /dev/stdout \

    -R "$fastaRef" 2> "$log" | \

grep -vi warn   `#0` | \

tail -n +2  `#1` | \

cut -f-2    `#2` | \

sed "s|:|\t|"   `#3` | \

tee >(awk -F"\t" '{SUM += $3} END {OFS=FS; print NR, SUM/NR}' > "$outFeatureStats") `#4` | \

gzip > "$outFeature"


#       1. Use tail to remove header
#       2. Use cut to extract chrom, pos, value columns
#       3. Use sed to parse into chrom\tpos\value format
#   4. Use awk to calculate the number of loci processed, and the mean feat value 
#added --includeRefNSites to make sure emits output for ambiguous bases (because read.feats will have)
echo Done.

---

[GenoMax]

Please post the complete command you are using to run this script.

[Maryem]

if you can check this page github (Run,Overview) : https://github.com/AshleyLab/scotch I used this command : python ~/scotch/scotch.py get-features-depth --project_dir=~/ABC123/ --chrom=1 --beds_dir=~/beds/ --fasta_ref=~/GRCh37.fa

[Maryem]

when i check the output i found this error message :

 A USER ERROR has occurred: Couldn't create output file, encountered exception: /dev/stdout.sample_interval_summary 


Set the system property GATK_STACKTRACE_ON_USER_EXCEPTION (--java-options '-DGATK_STACKTRACE_ON_USER_EXCEPTION=true') to print the stack trace.

Answer 1

[Pierre Lindenbaum]

https://gatk.broadinstitute.org/hc/en-us/articles/360041851491-DepthOfCoverage-BETA-#--output

Base file location to which to write coverage summary information, must be a path to a file Base file name about which to create the coverage information

this is the basename of all the files produced by DepthOfCoverage. There is no output on stdout, you can't pipe this. That's why you don't have any lines in your awk, that's why you got a division by 0 (NR==0)

[Maryem]

so i have to remplace stdout by something else ? i didn't get it

[Matthias Zepper]

Redirecting the output of a software to /dev/stdout is equivalent to output something to your current terminal, from where it can be consumed by other programs in a pipe. The current version of the GATK toolkit does no longer output anything on stdout, but exclusively writes the outputs into files. The awk step does a division by the number of lines received and because it receives no output, it tries to divide by 0, which obviously upsets the gods of math.

So to fix, this you either need to use an older version of the GATK toolkit, which was current when that Scotch software was last updated 3 years ago or can try to alter the script like so:

-o "$tmpDir" \

-R "$fastaRef" 2> "$log"


grep -vi warn "$tmpDir/*" | \

But there is no guarantee that splitting the command works, because the output may have changed in different ways. It is anyway quite a wild combination to have a Python script calling a Bash script wrapping a Java executable, so maybe there is a better maintained software available. I am somewhat biased because I know some developers, but I like Sarek for Indel calling ;-)

[Maryem]

Hi Matthias Zepper , Thank you for the explanation. I appreciate it. and thank's for the suggestion !! I will try it

[Matthias Zepper]

You're welcome. Good luck with your project!