More divide by zero errors

[kautto]

I ended up with another (different) set of divide by zero errors:

/igm/apps/python/Python-2.7.12_install_static/lib/python2.7/site-packages/scipy-0.18.1-py2.7-linux-x86_64.egg/scipy/stats/_distn_infrastructure.py:408: RuntimeWarning: invalid value encountered in double_scalars
  return m3 / np.power(m2, 1.5)
src/distributions.py:40: RuntimeWarning: divide by zero encountered in double_scalars
  p = (va-av)/(va)
src/distributions.py:31: RuntimeWarning: divide by zero encountered in double_scalars
  f0 = sm/(r+sm)-p
src/distributions.py:32: RuntimeWarning: invalid value encountered in double_scalars
  f1 = np.sum(special.psi(data+r)) - n*special.psi(r) + n*np.log(r/(r+sm))
src/distributions.py:32: RuntimeWarning: divide by zero encountered in double_scalars
  f1 = np.sum(special.psi(data+r)) - n*special.psi(r) + n*np.log(r/(r+sm))
src/distributions.py:32: RuntimeWarning: invalid value encountered in log
  f1 = np.sum(special.psi(data+r)) - n*special.psi(r) + n*np.log(r/(r+sm))
/igm/apps/python/Python-2.7.12_install_static/lib/python2.7/site-packages/scipy-0.18.1-py2.7-linux-x86_64.egg/scipy/optimize/minpack.py:161: RuntimeWarning: The iteration is not making good progress, as measured by the 
  improvement from the last ten iterations.
  warnings.warn(msg, RuntimeWarning)

The interesting thing is that if I extract reads for just one chromosome, this doesn't happen. But if I run it on our full WES or WGS files, it breaks.

[relyanow]

Huh. I tested it on a whole genome file without errors. I will do some more testing to get to the bottom of this.

[kautto]

Could it be reference genome specific? I think the manuscript mentioned it being tested on hg19, but these samples are on hg38/GRCh38, so there's a whole different set of contigs there...

[relyanow]

That could be it. I'll do some testing with an hg38 aligned genome

[kautto]

I guess I was wrong, since I tried running through a handful of hg19/b37 samples, and they all give me that same error.

[relyanow]

Thanks, I was away from my computer yesterday, will get back to you about this after some more testing today.

[relyanow]

I tested the code on 3 datasets and was unable to replicate your specific error. However, I added some error handling that will hopefully prevent this error from occurring. Let me if the error persists. If so, could you send me a few lines from one of your BAM files so I can see if there may be some differences in the file format?

Thank you!

[kautto]

I've tried this on both b37 and GRCh38 files, and a decent amount (about a dozen or so) files in total, and it happens on all of them. I'm not sure that I can send any of the data since it's all patient samples. I'll try to troubleshoot the code on my end and keep you posted, since it'd be nice to get this working.