New feature request: ANI distance

[cumbof]

Hi @rsharris
this is not actually an an issue, but a request to add a new feature in HowDeSBT if possible.

I’m using your software as a dependency of MetaSBT where I’m building SBTs starting from genomes, microbial genomes in particular.
It would be very useful if HowDeSBT could provide a way to compute the Average Nucleotide Identity (ANI) measure between two bloom filters with the bfdistance subcommand.

It would be very useful if it could also compute the ANI between filters while running a query. Let’s say that I have an SBT built with genomes. Now, I have a new genome (its bloom filter representation) and I want to establish which is the closest one in the tree according to the ANI distance.

Computing the ANI usually means performing alignment, but we could use the following formula as a very good estimation of ANI based on the number of active bits in the union and intersection between two bloom filters:

1 - (1 + (1/kmer_size) * log((2*jaccard_index) / (1+jaccard_index)))

Where kmer_size is obviously the size of the kmers used to build the bloom filters, and the jaccard_index is simply the number of active bits in the intersection of two bloom filters over the number of active bits at the union of the same two bloom filters intersection/union.

I’m currently computing this measure within MetaSBT by running the bfdistance subcommand twice with --show:intersect and --show:union. It works as expected, but of course it’s very inefficient. It would be much faster if HowDeSBT could provide this feature natively.

Let me know if it makes sense and if there is any chance we could see this feature implemented in the near future.
Thanks!

[rsharris]

Having taken a glance at bfdistance.c, I'll add jaccard as a distance metric for bfdistance. That will be pretty easy to do given the other metrics already implemented.

I'm less sure about the formula for estimating ANI. log in that formula is probably log-base-e, or maybe log-base-2, right? I kind of expect there should be an e^something wrapping part of that expression. Do you have a citation for that formula? If I implement it I want to include a citation as a comment in the code.

[rsharris]

I've added --show:jaccard to bfdistance and pushed it as version 2.00.17.

AVVERTIMENTO!! I haven't tested that at all, beyond compiling it. Please test it and let me know if it works.

I would test it, but it would take me a while to familiarize myself with the package and how I've tested stuff in the past. I haven't used this much in the last five or more years, so all that is tucked in a back closet in my brain. And, at the risk of repeating these famous last words: "this function was so easy to implement it couldn't possibly be wrong."

[cumbof]

Thanks Bob! I’ll test it out and I’ll let you know ASAP!
For what concerns the citation: https://doi.org/10.1186/s13059-016-0997-x Formula (4)

[cumbof]

In any case, I think the Jaccard is enough. I can easily compute that formula within Python.
But yeah, this reduces the number of calls to HowDeSBT from 2 to only 1, and this is definitely going to make my pipeline way faster than it is now

[rsharris]

Thanks for the reference. With Paul Medvedev, I worked on a couple of papers related to that. I didn't immediately recognize it in the form you have it, but I should have. In any case, "log" in that format is natural log (base e). One thing to be aware of with that estimator (mash distance) is that it's great when sequences are 'close', but the further the get the further the estimator is from the truth (from the standpoint of expectation).

I'm still not completely on board with the formula you have though. It seems like ANI should be 1-D. D is defined in the paper as -(1/k)ln(2j/(1+j)) where j is jaccard index. So 1-D would be 1+(1/k)ln(2j/(1+j)). Your formula is subtracting that from 1, so it seems like your formula is 1-(1-D) which is D.

[cumbof]

Yes, that’s natural log, correct, and yes, I agree with you that the further two sequences get, the further the estimator is from the truth.

For what concern the last point, I should have a look at the formula again but yes, your interpretation seems correct to me.

I haven’t been able to test —show:jaccard yet. I’ll probably go through that tomorrow, but I’ll let you know for sure if there is any issue.