Minhash scaling code does not return hashes for many sequences

[etongit]

I have a very large list of sequences for which scaling using Minhash does not return any result. Even when scaled < len(sequence)

Example sequences -
ATGGTTGGGATCGACGGACCGTAAATATCGGCATCGAGAACCCCGACCTTTGCGCCTTCAGCCGCTAACGCCAGCGCCAGGTTTACCGCCGTGGACGATTTCCCCACCCC
ATTAGTAAAAAACATGAGCATGGCCTGGCAAAATGTACTGTATATCGTGGCCGCGATATTAGTAATCATGCTGTGCGTCTTTACGCTGATCATTCGCGGTAAAGCCAAAAGCGA

Minimum code reproduce -

 mh = MinHash(n=0, ksize=40, scaled=100)
 mh.add_sequence(sequence, force=False)
 print(mh.hashes.keys())

This prints -

KeysView({})

Where as

mh = MinHash(n=100, ksize=40)
works just fine.

[ctb]

oof, that's weird!

What happens with scaled=1? That will save all k-mers => hashes.

If you get nothing - then, for some reason, the file is not being read properly!

If you get something, then maybe further diagnosing needed 😅

[etongit]

scaled=1 works. It stops to work beyond some values (but only for some sequences). There is definitely a bug in the code.

[ctb]

(I'm not saying there's no bug, but I would be surprised if there were one this obvious, since we use sourmash quite heavily for lots of things ;).)

I would then also expect

sourmash sig downsample scaled1.sig --scaled 100 -o scaled100.sig

to produce empty output. Hmm.

So if it's not the sequence format, maybe it's the sequence content. Would you be able to attach a subsampled file here, or send some of the sequence to me at [email protected]?

You could also try intermediate values of scaled to see where it's losing hashes.

[etongit]

I linked two example sequences above, would that be enough?

# First sequence
ATGGTTGGGATCGACGGACCGTAAATATCGGCATCGAGAACCCCGACCTTTGCGCCTTCAGCCGCTAACGCCAGCGCCAGGTTTACCGCCGTGGACGATTTCCCCACCCC
# Second sequence
ATTAGTAAAAAACATGAGCATGGCCTGGCAAAATGTACTGTATATCGTGGCCGCGATATTAGTAATCATGCTGTGCGTCTTTACGCTGATCATTCGCGGTAAAGCCAAAAGCGA

[etongit]

So I ran this code to do a quick test -

from sourmash.minhash import MinHash

def compute_minhash(sequence, scale):
    #mh = MinHash(n=scale, ksize=40)
    mh = MinHash(n=0, ksize=40, scaled=scale)
    mh.add_sequence(sequence, force=False)
    print(f"{scale}: {len(mh.hashes.keys())}")

SEQUENCE1="ATGGTTGGGATCGACGGACCGTAAATATCGGCATCGAGAACCCCGACCTTTGCGCCTTCAGCCGCTAACGCCAGCGCCAGGTTTACCGCCGTGGACGATTTCCCCACCCC"
SEQUENCE2="ATTAGTAAAAAACATGAGCATGGCCTGGCAAAATGTACTGTATATCGTGGCCGCGATATTAGTAATCATGCTGTGCGTCTTTACGCTGATCATTCGCGGTAAAGCCAAAAGCGA"

for i in range(1, 110):
    compute_minhash(SEQUENCE1, i)

For SEQUENCE1, upto scale=31, there are values and after that it is 0. Scale=1 is 71 hashes.
For SEQUENCE2, upto scale = 49, there is are values and after that it is 0. Scale=1 is 75 hashes.