Assertion error with large dataset: sorted(tetra_z[org1].keys()) == sorted(tetra_z[org2].keys())

[ghbio]

Summary:

Assertion error sorted(tetra_z[org1].keys()) == sorted(tetra_z[org2].keys()) when using method TETRA in a large dataset (~3000 bacterial genomes)

Description:

When running average_nucleotide_identity.py with method TETRA in a dataset consisting on 3000 prokaryotic draft genomes it stops with the above mentioned AssertionError. The same command with a smaller dataset (500 genomes) in the same system works as expected.

Thanks for your help and developing and maintaining the software!

pyani Version:

pyani version 0.2.7

Python Version:

Python 2.7.16

Operating System:

CentOS 7

[widdowquinn]

That sounds odd, as if something in the other 2500 genomes is the issue.

Would you be in a position to try the other 2500 in batches of 500 to help try to isolate the problematic genome(s) so I can debug further?

[ghbio]

Sure... it'll take a while but I'll get back to you. Thanks again

[widdowquinn]

Thanks - smaller batches might be just as good (and will be quicker). It might be we miss the specific buggy interaction, though.

If they all seem to pass, I'll provide an alternative script that should save out much more state to help us track down the problem.

[ghbio]

I've run several batches o 500 genomes and through successive subsetting steps I've narrowed it down to a set of 10 genomes, which trigger the issue. They have nothing particular as far as I can see... do you have any suggestion or would you rather get the sequences to examine them yourself?

[widdowquinn]

Thanks for doing that! I'm afraid I can't really identify the issue without the sequences that cause the failure - and I can't test if I've fixed the problem without them, either. If they're not confidential, you could compress them and attach them here - if they are, then there are a number of more private transfer methods we could use.

[ghbio]

I'd rather send them to you privately. It's just a 1.5 Mb tarball.

[widdowquinn]

That's fine - my email address is leighton.pritchard@hutton.ac.uk. I'll treat them as confidential and delete when I'm done.

[widdowquinn]

Hi Pablo - that's the fix made.

The issue was that the code (which means I) was expecting all 256 4-mers to be present in each sequence. One of your sequences did not have one of the 4-mers present. The new code doesn't make this assumption (it's not a reasonable assumption…).

I've pushed the change as a release candidate, and I'll make a revised release on PyPI (and hopefully also on Bioconda) later today. In the meantime, cloning/installing the current HEAD from the repo (da53544) should get you the fix immediately.

Thanks again for using pyani and bringing this to my attention, and I'm glad the profile pic is appreciated :D

L.

[ghbio]

Hi Leighton,

thanks for your reply. I'll check the fix as soon as I can and I'll get back to you if I find any issue. Our dataset consists on bins from metagenomic reads. We have some fragmented genomes or even small viral genomes that may not fulfill the assumption you talk about. Again, thanks for your help.

Cheers,

Pablo

[widdowquinn]

Yes, it was a bad assumption for me to make!

I'll close this just now, but please do reopen the issue if you find a related problem.