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.
•
written
by
Maryem •
1 answer
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)
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 ;-)
Hi Matthias Zepper , Thank you for the explanation. I appreciate it. and thank's for the suggestion !! I will try it
You're welcome. Good luck with your project!
Log in to answer this question.
Please post the complete command you are using to
runthis script.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
when i check the output i found this error message :