genrich fails with floating point error (136) on too much input

[malcook]

I am getting floating point error while using genrich on

314,153,343 TOTAL READS

coming from seven name sorted bams of replicate atac-seq data with

read count 21,001,420 31,212,472 56,249,110 43,494,211 53,324,666 42,859,376 66,012,088

any advice or request for data to reproduce?

my workaround is to drop some replicates.

[jsh58]

At what stage does the error occur? You should get some indication of this from the verbose output.

If there is a subset of the data on which I can reproduce the data, please provide it and I will have a look.

[malcook]

running

Genrich --version Genrich, version 0.6 Copyright (C) 2018 John M. Gaspar (jsh58@wildcats.unh.edu)

run1 quit with floating point error output ends right after it shows it has begun processing the 2nd ("#1") file, as below. run2 I swapped the order of the 1st two input files. Surprise! It completed: "Peaks identified: 18324 (6097652bp)" run3 I swapped them back (superstition) and same error as run1. run4 I ran with just 1st 2 input file and same error same error as run1

Genrich --verbose -S -j -f AC_LL_H_000.narrowPeak.pq.bedgraph -r -e chrM -g 300 -a 5.0 -y -t 'AC_LL_H_000.10.n.bam AC_LL_H_000.3.n.bam' -q .05 -o AC_LL_H_000.narrowPeak &> AC_LL_H_000.narrowPeak.stderrout

run1 log file: https://research.stowers.org/mec/genrich72/AC_LL_H_000.narrowPeak.stderrout

2 input files: https://research.stowers.org/mec/genrich72/

[jsh58]

I am unable to download AC_LL_H_000.10.n.bam. I get only a truncated 62KB file.

[malcook]

apologies - can you please find them instead in https://ftp.stowers.org/#/public/mec/genrich72/

[jsh58]

Thanks for the files. I was able to reproduce the error, diagnose the bug, and fix it. It occurred because the file AC_LL_H_000.3.n.bam had no unpaired alignments, and Genrich failed to recognize this and threw this error when trying to remove PCR duplicates from this nonexistent set.

I reran the updated Genrich on the input files, twice:

$ ~/Genrich/Genrich -t AC_LL_H_000.10.n.bam,AC_LL_H_000.3.n.bam -echrM -g300 -a5 -q0.05 -o run1.peak -jyrv
$ ~/Genrich/Genrich -t AC_LL_H_000.3.n.bam,AC_LL_H_000.10.n.bam -echrM -g300 -a5 -q0.05 -o run2.peak -jyrv

They produced identical outputs, as expected. Here is the verbose output:

Peak-calling parameters:
  Genome length: 1371702787bp
  Significance threshold: -log(q) > 1.301
  Min. AUC: 5.000
  Max. gap between sites: 300bp
Peaks identified: 1450 (681706bp)

Sorry for the difficulty. Thanks for letting me know and providing the files!

[malcook]

Fantastic!

Will you make another versioned release I can refer to in my METHODS?

Cheers,

Malcolm

[jsh58]

Yes, I just made v0.6.1.