National Energy Research Scientific Computing Center (NERSC) NERSC 2022 Early Career Achievement Award Seminars, 14 Dec 2022

Previous Meeting Next Meeting

⏯

youtube image

►

From YouTube: Chirag Jain: Sketch-based algorithms for large-scale whole-genome comparisons

Description

Sketch-based algorithms for large-scale whole-genome comparisons
Chirag Jain, Indian Institute of Science

December 14, 2022
9:00-9:30 Pacific Time

In the era of exponential data growth, sketching has become a standard algorithmic technique for rapid genome sequence comparison. I will describe our work on fast, lightweight approximate sequence mapping algorithm by using minimizer sampling and MinHash techniques. The proposed algorithm computes the positional origin of a query sequence in a given reference and estimates nucleotide-level identity under an assumed probabilistic model of mutations. We show an application of this algorithm in quantifying relatedness between two microbial genomes and its impact in tackling a long-standing biological question. Microbiologists are increasingly turning to whole-genome sequencing driven approaches to address fundamental questions associated with ecology. This involves quantifying similarity of two or more genomes, e.g., to check whether a newly sequenced genome is novel, or where else has it been seen before. We developed FastANI (Average Nucleotide Identity) software by using the proposed approximate sequence matching framework to quantify similarity of two or more genomes. Our algorithmic improvements, coupled with parallelizability allowed us to index entire database of 90,000 bacterial genomes, and compute pairwise ANI values among all pairs of genomes for the first time. This analysis sheds light on the extent to which microbes form discrete clusters (species).

A

So yeah and thank you all for joining this morning so good morning, everyone I'll speak about our uh work on sketch, based algorithms, for genome sequence, comparisons and then I'll show you their applications on large-scale whole genome comparisons.

A

So I think I would say that sequence alignment is a fundamental problem in bioinformatics and this is a problem where there is always a need for faster algorithms.

A

So uh I guess for those who have don't have background in bioinformatics. Here we are talking about query sequences, query, DNA sequences and for these sequences we are looking for approximate matches in a large set of reference sequences.

A

So in this toy example, here you see a query sequence and it has a Target match here with 80 identity, because 4 out of five characters are matching, but this is just a toy example in real situations we have long query, sequences and massive reference databases so as one example, the refseq database at ncbi, the number of genomes currently add up to 3.2 trillion, base pairs there and similarly, the rate of the data generation is also high in this field. So a single sequencing machine can produce one more than one terabits per data per experiment.

A

So, given this massive data, we always need faster algorithms for doing sequence. Analysis in this area.

A

So the problem statement looks like following so given a query sequence. We would like to identify all target mapping positions in the reference where it matches with greater than Upsilon percentage alignment identity where Epsilon is a user-specified cut off that they want for their application.

A

For example, Epsilon could be 90 or 80 percent, depending on the Divergence that the user wants to tolerate in their application.

A

Now again, given a query, we want to identify all target mapping positions and the exact algorithms for solving these problems have been known for long, but they need quality Grant time in terms of the query, size and the residence length.

A

So as a result of the quality complexity, these Exact Solutions cannot scale to the data sets we are talking about.

A

Now. There exists extreme heuristics, which are an alternate to the exact solutions. For example, blast is a popular software that implements till an extent realistics and the more recent software has minimap2, but these heuristics are inefficient for this task, and the second problem with these heuristics is that they lack theoretical guarantees on the output.

A

So here in this work, what we are looking for is we are looking for a approximate algorithm that can guarantee giving us this output with high probability, and this approximate solution will make use of Gamers now. Kmarts are nothing but short words of length k throughout this presentation, I'll be denoting using this circular. Dots I'll be denoting these scammers from both of these sequences.

A

Now we will use these schemers and we will use Jacquard similarity as a criteria to judge the alignment quality. So, for example, if I consider a particular substring from the reference B so suppose, I want to decide whether query a aligns well at position. I of reference b or not I'll.

A

Consider all KMS from query a and I'll consider k-mer's from this particular region of the reference and assume this region is of same length as the query I will consider the Jacquard similarity of the Kmart sets extracted from both of these regions and consider their intersection divided by the union. So the card similarity is nothing but intersection, divided by the Union and the goal. The reason for using jacquardes lady I'll come to it next in the next Slide.

A

The reason for using Jacquard similarity is that there are known techniques in literature to efficiently estimate the card. Similarity for example, one his Min hash technique from broader that was developed in 1997., and then we know minimizer techniques that are developed in database community that have been originally developed for scalable web document clustering.

A

So there are efficient techniques to estimate the cardinality and the second reason for using the card. Similarity is that if we assume, uh if you assume that mutations occur randomly in these sequences and suppose this mutations occur using some poisson distribution model, then it is possible to correlate sequence, identity and the card similarity.

A

In fact, we can show that the maximum likelihood estimate of sequence identity can be expressed in terms of Kmart length and jacquard's miracle, so this relationship helps us in going using the card similarity as a proxy or alignment identity, foreign.

A

Now the problem reduces to identifying all the positions in the reference database where this cutoff is satisfied.

A

So we are looking for all such positions where this Jacquard similarity exceeds the a certain cutoff indicated by Tau symbol here now there are two challenges here. First of all, we cannot evaluate all positions in the reference one by one, because that will be time computationally time probability and Second. Challenge here is that we do not want to store all payments in the reference database because, as I said, reference databases have more than trillion characters and therefore storing trillion Gamers from the reference will require massive amount of space.

A

So we want to optimize and address these two challenges here and in the next couple of slides. I want to show you how we address these two challenges, one by one.

A

So, first of all, I want to show you that we can estimate this. The cardinality using a new estimator that combines the known techniques. One is first, is the minimizer technique and second, is the mean hash technique now assume that gamers in a sequence are hashed into integers, so instead of Gamers I can use these numbers to indicate each gamer.

A

Minimizer sampling is nothing but a considering minimum gamer in conjugative Windows of fixed length. So, for example, if my window size is 5, I'll consider the Kmart with minimum hash value in the first window and I'll sample that Cayman from the first window. So here the camera with hash value, 2 is a minimum payment. I'll slide the window, wipe one position and I'll consider the next again. The same gamer is: has the minimum hash value in the next window, and so on? The method window has a new value which is three and so on.

A

So I'll continue to slide this window at one position at a time and I'll keep on sampling Kmarts with minimum hash value in each window.

A

So the the sample set of minimize sample set of Gamers is known as minimizers on top of the sample set of minimizers I'll apply this sketching technique where, for a particular fixed sketch size, let's say s equals three I'll further.

A

Consider the minimum s elements out of this set and using a combination of these two techniques, I will Define a new new Jacquard similarity estimator between two sequences, A and B.

A

Now the reason this estimator is helpful will become clear in the next slide. So we go back to the two challenges that we wanted to address. First, was we wanted to avoid storing all cables and the second was we don't want to evaluate all positions in the reference so for the first challenge, what we notice is that for the using the new estimator, we don't need to store all Caymans from the reference we. It suffices to store only the minimizers in memory and from the slimer's result.

A

We know that the expected set of minimizers that are sampled is the total sequence length divided by a factor of w by 2, where W is the window size that we used for minimize the sampling.

A

So this saves us Space by a factor of w by 2.. So this addresses the first challenge that we don't want to store all cameras in memory for the Second Challenge. We show that, if similarity exceeds cut off, then the number of minimizer matches have to be at least s times Tau, and what the reason why this helps is that we only if we consider and minimizer matches from the read to the reference. We only need to consider those positions where there are sufficient minimized effects.

A

This helps with the time issue that we don't need to evaluate all positions and we can only compute the card similarity and those positions in the reference where this cutoff is satisfied, and this helps us accelerate this process overall, using the new estimator we in the paper, we show that the proposed algorithm reports all desired mapping positions with high probability, so this algorithm also gives a guarantee on its output and which is a nice part. That is not only. It is fast, but it also offers theoretical guarantees.

A

So next, given this algorithm, I want to apply this algorithm in real world setting where we want to do over genome comparisons at a massive scale.

A

Now, why do we want to do this whole genome comparisons, and this is for a particular application for species delineation in microbiology now whether bacteria forms species or discrete clusters has been debated for decades in microbiology literature and this particular the particular reason for this debate has been that the microbes are capable of exchanging genes through mechanism of horizontal Gene transfer.

A

So due to this genetic transfer possible among them, it's not clear whether they maintain these discrete genetic boundaries or not.

A

So our in this work, we wanted to take a data driven approach where we wanted to do large scale, slusting analysis of all publicly available genomes and see whether the results indicate discrete clusters or not.

A

We'll be using this measure called average nucleotide, Identity or Ani. In short, now Ani is commonly used in my in genomic literature to first species classification. So it is helpful to quantify similarity between two genomes and computationally. It is calculated as mean alignment identity of conserved genes between two genomes.

A

So if there are, uh if a genome has thousand genes, then we look for the conserved pair of genes between two genomes and we average their identity to compute. This A N I measure which, if and the higher the Ani, is more similar to more similar. The two genomes are.

A

If we look at the literature for a I solvers, the first Ani solver was implemented using blast tool, but this last based AI solver, was too slow and over the years people replaced blast with other alignment. Software, for example, number and you search, but still the maximum speed up that we got was 4.7 x over the blast based implementation.

A

Recently there are other approximate methods, but they were found to be inapplicable for draft assemblies.

A

So we used our Jacquard similarity based sequence, matching framework and implemented this software called Fast a I to show that this faster software is accurate and it is two to three orders of magnitude faster than blast based a I implementation.

A

So the implementation works simply by looking for reciprocal best matching genes which tools identify the conserved set of genes between two genomes, but instead of using alignment based methods. Now we are using our sketch based algorithm for sequence, comparisons and note that we can estimate the alignment identities directly from Jacquard similarity. We don't need to explicitly do sequence alignment for this task.

A

If we look at the results on how this works in practice, we used five reference databases with different composition.

A

So if you see there are five databases with different set of genomes and in each case we have a single query genome that we test foreign each database. So in each case we are querying, we are Computing pairwise, a and I of the query genome against all in the reference database.

A

First, we show the runtime comparison with blast based a I method. So here we see that if we ignore the one time indexing cost for the database, the compute time against the blast based a I method is two to three orders of magnitude faster for all. Five data sets.

A

This is single thread single process execution, and the reason for speed up here is that we are avoiding the need for alignment step and all the alignment. The Ani estimation is happening directly: Through The, jacquard's Narrative of gamers.

A

Now the next question is: are we compromising accuracy if we are faster to do this? We check the correlation of estimated a I with the blast based A9 method, assuming that the blast based a I method is accurate and we also compare our method with the another approximate tool which is match. So Mash is also a popular tool in this area, which does alignment 3, Sequence comparison, and here we show that password gives those two ideal correlation for all five data sets.

A

So in all five cases the black line is indicating x equals to y, and if you see all the A9 estimates are closely following the Y equals X plot. In all five data sets.

A

So this also shows that Fastener is accurate in uh variable settings so now that this method is accurate and significantly faster than previous methods, we wanted to use this method to do a large scale, clustering of all publicly available genomes, and for that we did a all to all pairwise comparison among all 90 000 genomes that are available in genbank, and this analysis was used to detect whether there are discrete clusters that are formed in bacteria or not.

A

Now, in the results on the right hand, side I am showing you the distribution of pairwise Ni values that we see from this all to all pairwise comparisons.

A

What we see here is that most of the pairs that corresponding to same species had a Ani value of more than 95 and most of the pairs corresponding to different species had a n. I value less than 83., with a few exceptions in the middle, so most of the pairs corresponding to different name species have dni value less than 83, leaving a big gap in the middle.

A

So this bio model distribution indicated us that indeed, these genomes form discrete clusters and this whole experiment, even though fast is orders of magnitude faster than previous methods. This whole experiment took us 77, 000 CPU, core hours, which was paralyzed across multiple threads and nodes to achieve reasonable time to completion.

A

Now one limitation of this experiment is that we are using public, this all the publicly available genomes, and we know that some of the species have been extensively sequenced compared to others.

A

So there is a there's a there could be a problem of sampling bias in database because some species are over represented and some species are underrepresented in the database.

A

But recently what we found out there was a similar distribution reported by a study from UC Berkeley, where they tested pasty and I on 5000 genomes that were derived from metagenomics.

A

So these genomes were extracted directly from natural habitats such as ocean soil, gut microbiome and so on, and in all cases we are seeing by model distributions like we observed in our in our case.

A

So these results sort of shed light on the fact that, indeed, these bacterial genomes form discrete clusters and 90. This 95, a I cut off, can be used as a reasonable measure for species classification.

A

So, to summarize, we discussed a fast sketch based algorithm to map sequences to large genome reference databases, and this was enabled by a noble estimator for Jacquard similarity.

A

This tool provides positional as well as identity estimates for each query sequence.

A

We saw that it is applicable for computing or genome pairwise similarities at scale, and when we run such analysis at a massive scale, we can always gather some new insights that were not possible before so. For in this case, we saw that we saw got new insights into bacterial species population genetics using fast in act based pastry, which indicated that in indeed these species form discrete clusters in the environment, so I'll close by concluding all the collaborators in these projects.

A

So this projects were done in collaboration with my colleagues and mentors at Georgia, Tech and NIH, and of course, I have to thank nurse because all the compute we redirect nurse- and uh this and many other projects are for ongoing projects are being supported by nurse resources and I'm. Very thankful for that.

A

Thank you and I'm happy to take questions.

B

Yes, thank you very much. That was really interesting. uh Are there any questions.

B

Just speak up: if you have any I'll start with one, um do you have a particular uh problem or one you know, of that um you're applying this to next or what's what's the next? What's the next big thing that this technique is going to be applied to.

A

Actually uh we have um so this is the general General technique and I would say. Whole genome comparison is just one of the applications that it can be used for, for example, we ourselves have used it for multiple applications, uh for example metagenomic profiling alignment of long reads: uh people are also using it for uh pan genome analysis, panginum graph, construction and so on.

A

So in general, we think that, in contrast to existing approaches, this has a unique ability of providing a guarantee on identifying all target mappings above a user specified cut off which, and such guarantee is not offered by existing tools. So in all the applications where that's a guarantee helps. uh This is useful.

B

Okay, great and um how is the code implemented and can it take advantage of gpus now.

A

No, so it's a CPU only right now we use openmp for uh multi-threading and um and actually I would say it's a mostly embarrassingly parallel problem, because if you see all this clustering based analysis where we are doing uh all these all to all pairwise comparisons, all these are independent computations. That can happen and and even if the there's, a single query genome still, we can split the reference database across threads to do the compute. So the implementation is in C, plus plus with openmp support.

A

uh I would say that mostly in for sequence, analysis, tasks and especially sequence mapping, because the computation is highly irregular, uh mostly I've, seen that CPU based implementations currently offer better performance for such tasks, at least. If I look at the state of the art, implementations.

B

Okay is: is this a? Is this a technique or an algorithm that um can benefit from more memory on the nodes.

A

Yes, more memory can be helpful. For example, uh like usually the memory consumption of this algorithm is dependent on the sampling rate. Actually, so, as I said, the sketching methods, work by sampling of cameras and the sampling rate is decided by the algorithm, depending on the Epsilon cut off like this, and if the silent cutoff is uh relaxed like, for example, if the user is wanting low identity cut off, then this algorithm would have to sample at a higher density.

A

So so it improves the specificity of the algorithm if we sample more so yes, so in, in a sense that, yes, the algorithm will be more accurate if we reduce the sorry if we increase the sampling. Okay, all right.

B

Okay, any other questions. Let's see, um if not I want to thank you again, and um it was really great to hear from you and again. This recording will be available soon, um so have a good rest of the day. Everybody- and thank you.

A

Thank you, bye.