matUtils summary non-deterministic

[aofarrel]

From @majagarbulinska with title matutils summary random seed problem? Different results with the same command.

Hi! Thank you very much for the package!

matUtils summary -t test1.tsv -i public-latest.all.masked.pb -g ncbiGenes_.gtf -f NC_045512v2.fa.gz
matUtils summary -t test2.tsv -i public-latest.all.masked.pb -g ncbiGenes_.gtf -f NC_045512v2.fa.gz

diff test2.tsv test1.tsv | grep "^>" | wc -l
gives a number >0

Exactly the same command with exactly the same files give us different results. Would there be a possibility of including a random seed parameter to get the same results every time we run the same command for reproducible results?

Thank you!

This was originally posted in the usher wiki repo in Sep 2021: bpt26/usher_wiki#5
I have not yet attempted to replicate this issue with the current version of matUtils.

[aofarrel]

I spent some time trying to replicate this, and in summary: I'm getting slightly different results than OP, but the underlying problem (matUtils summary is non-deterministic) is definitely real as of matUtils 0.7.0.

I don't know which files OP is using so I used this pangolin tree, this fasta which I --fast gzipped, and this gtf in the pathogengenomics/usher:latest Docker image:

matUtils summary -t test1.tsv -i dvol/covid/public-2021-06-09.all.masked.nextclade.pangolin.pb -g dvol/covid/GCF_009858895.2_ASM985889v3_genomic.gtf -f dvol/covid/NC_045512v2.fasta.gz
matUtils summary -t test2.tsv -i dvol/covid/public-2021-06-09.all.masked.nextclade.pangolin.pb -g dvol/covid/GCF_009858895.2_ASM985889v3_genomic.gtf -f dvol/covid/NC_045512v2.fasta.gz

... and while unlike OP, I'm not getting any differences in diff (I tried both gnu diff and Apple diff)...

> diff test2.tsv test1.tsv | grep "^>" | wc -l
0

...it turns out the output is indeed non-identical.

>md5sum test1.tsv
db768762023908fd62ea54a062487511  test1.tsv
>md5sum test2.tsv 
8ca12a40669162e28253fa5e842dfe41  test2.tsv
>cmp -l test1.tsv test2.tsv | head
     144  65   0
     145  61   0
     148  65   0
     149  61   0
     294 140 241
     298 140 241
     480  57 121
     484  57 121
     591   0 117
     593   0 124

These TSVs do not include the "loaded in xxx msecs" tags, so it's not that. There's something wrong with the codon_changes column. Note how test1.tsv got G51>T51 but test2.tsv got G>T.

>head test1.tsv -c 150
node_id	aa_mutations	nt_mutations	codon_changes	leaves_sharing_mutations
CHN/YN-0306-466/2020|MT396241.1|2020-03-06	GU280_gp01:X5216X	G15910T	G51>T51	
>head test2.tsv -c 150
node_id	aa_mutations	nt_mutations	codon_changes	leaves_sharing_mutations
CHN/YN-0306-466/2020|MT396241.1|2020-03-06	GU280_gp01:X5216X	G15910T	G>T

Upon repeating a few more times...

>md5sum test*.tsv       
db768762023908fd62ea54a062487511  test1.tsv
8ca12a40669162e28253fa5e842dfe41  test2.tsv
1edd937dce765b95de872dcd63700b2b  test3.tsv
1239471498b087e9267efb6b231465b7  test4.tsv
86fc6fc7dc64aee6022a8e57332a8343  test5.tsv
71952c5b3ab2fe83ea1e20b5fea9be06  test6.tsv

...it looks like you can get G_l>T_l too. Sometimes.

head -n 3 test4.tsv
node_id	aa_mutations	nt_mutations	codon_changes	leaves_sharing_mutations
CHN/YN-0306-466/2020|MT396241.1|2020-03-06	GU280_gp01:X5216X	G15910T	G_l>T_l	1
node_2	GU280_gp02:X1117X	A24913T	A>T	2

[aofarrel]

@AngieHinrichs looked close into this and discovered this doesn't seem to happen when the fasta isn't gzipped. matUtils summary expects a non-gzipped reference sequence, and when it doesn't get one, it behaves unexpectedly.

[AngieHinrichs]

Thanks @aofarrel -- yes, there are at least two places where this error could have been detected instead of happily producing garbage output:

• translate.cpp build_reference makes no attempt to verify that the input file is plain fasta as expected. It should look for '>' at the start of the file and make sure that non-header lines contain only alpha or whitespace characters.
• translate.cpp build_codon_map does not perform bounds checking on the reference sequence when building codon strings using coordinates from the GTF input.

In addition, it would be nice to detect and support .gz input, but at a minimum those errors should cause matUtils to exit with a helpful error message.