WEBVTT Kind: captions Language: en 00:00:00.960 --> 00:00:07.020 Awesome, yeah, thank you for the introduction.  Hopefully, folks can see my screen. Yeah,   00:00:07.020 --> 00:00:11.220 so hey everyone. As mentioned, my name is Niema  Moshiri. I’m an Assistant Teaching Professor in   00:00:11.220 --> 00:00:16.560 the Computer Science and Engineering department  at UC San Diego. My talk is going to be focused on   00:00:17.100 --> 00:00:23.700 some methods that my lab developed using the NSF  funding for speeding up viral genomic analysis.   00:00:24.540 --> 00:00:28.980 Specifically, today's talk is just going to  focus on how we've enabled massively scalable   00:00:28.980 --> 00:00:34.080 reference guided multiple sequence alignment  of complete viral genomes. We actually have   00:00:34.080 --> 00:00:39.360 done multiple other accelerations as well that  I just didn't have time to talk about today.   00:00:39.360 --> 00:00:43.440 I'll end with a link to my website if folks  are curious about how else you could speed   00:00:43.440 --> 00:00:49.320 up other aspects of this type of analysis. So let's let's get started. Just to give a   00:00:49.320 --> 00:00:54.360 little bit of context - here's a framework for  a standard viral phylogenetics workflow. And,   00:00:54.360 --> 00:00:57.840 you know, before I even talk about this, viral  phylogenetics is very important to be able to   00:00:57.840 --> 00:01:04.980 study how is the virus mutating over time. How is  it, kind of, branching off? How are the different   00:01:04.980 --> 00:01:09.900 samples that we collect all over the world  related? There's a ton of uses in the world   00:01:09.900 --> 00:01:14.220 of viral molecular epidemiology that are out of  the scope of this talk, but generally, having   00:01:14.220 --> 00:01:19.920 a phylogeny inferred from viral genomes is very  useful to have. Typically the workflow starts like   00:01:19.920 --> 00:01:24.900 this where you start with a bunch of unaligned  viral genome sequences, which I'm showing over   00:01:24.900 --> 00:01:30.840 here. The first step is usually multiple sequence  alignment where you try to kind of place these   00:01:30.840 --> 00:01:35.760 gaps in the different kind of positions of  each of the sequences to get them to line up   00:01:35.760 --> 00:01:41.940 better. This gives you some notion of sequence  homology after you do this. Then, given the   00:01:41.940 --> 00:01:47.160 multiple sequence alignment, we can then perform  phylogenetic inference to try to infer an unrooted   00:01:47.160 --> 00:01:53.760 evolutionary relationship between these sequences.  Then, typically after that, we do what's called   00:01:53.760 --> 00:01:59.640 rooting to determine what is the most likely  common ancestor of all of the sequences. That   00:01:59.640 --> 00:02:04.440 kind of then tells us what was the forward in  time evolutionary history of these sequences.   00:02:05.100 --> 00:02:09.720 Then, maybe you'll do some additional downstream  analyzes. Maybe you do transmission clustering.   00:02:10.560 --> 00:02:14.820 There's there's a lot of other analyzes that you  can do on the phylogeny and on the sequences.   00:02:15.960 --> 00:02:19.740 But this is kind of the the building blocks  of how you do all these other analyzes.   00:02:20.940 --> 00:02:26.520 So typically these steps over here are the  key computational bottlenecks. The multiple   00:02:26.520 --> 00:02:31.320 sequence alignment and then the phylogenetic  inference. In today's talk, I won't be talking   00:02:31.320 --> 00:02:35.340 about phylogenetic inference I'm just going to  be zooming in on multiple sequence alignment.  00:02:37.020 --> 00:02:40.500 So, some context - multiple sequence  alignment this is what's called an   00:02:40.500 --> 00:02:46.200 NP-Complete computational problem. What that  means - there's a very technical computer   00:02:46.200 --> 00:02:51.840 science term - but basically what this means  is there's no polynomial time exact solution.   00:02:51.840 --> 00:02:56.820 Basically it gave me a bunch of sequences  and asked me to come up with the optimal   00:02:56.820 --> 00:03:01.440 multiple sequence alignment. There is no way to  do this in polynomial time. It's very very slow.   00:03:03.000 --> 00:03:09.540 Heuristics have been developed to provide optim-  to provide approximate solutions. For example,   00:03:09.540 --> 00:03:13.440 you might have heard of ClustalOmega, MUSCLE,  and MAFFT. These are some kind of standard tools   00:03:13.440 --> 00:03:18.540 that are used in the space. However, even these  heuristics - they generally scale quadratically   00:03:18.540 --> 00:03:24.480 with respect to the number of sequences.  For context, the GISAID database, which is   00:03:24.480 --> 00:03:30.900 the database where most folks are storing their  full SARS-CoV-2 genomes, this database is growing   00:03:30.900 --> 00:03:38.400 extremely rapidly and as of today we have over 15  million SARS-CoV-2 sequences available from all   00:03:38.400 --> 00:03:43.680 over the world. The next epidemic is going to  be more and kind of sequencing genomes in real   00:03:43.680 --> 00:03:48.660 time. This is going to be a tool that hopefully we  continue to use in in the viral epidemics to come.   00:03:49.620 --> 00:03:56.280 We can expect this to be even more significant  of a big data problem. Currently with these tools   00:03:56.280 --> 00:04:01.140 like ClustalOmega, MAFFT, and MUSCLE, we're  looking at runtimes of decades to centuries,   00:04:01.140 --> 00:04:06.360 which, you know, for obvious reasons, if we're  trying to do real-time molecular analysis,   00:04:07.080 --> 00:04:13.260 decades or centuries is just a little bit too  slow. So how can we speed this up? Well, it turns   00:04:13.260 --> 00:04:17.880 out that the problem is actually a little bit  easier than what we're trying to solve. Multiple   00:04:17.880 --> 00:04:24.180 sequence alignment, in general, is kind of  assuming no homology of the sequences whatsoever.   00:04:24.180 --> 00:04:30.480 This is the time it takes to align completely  arbitrary sequences. But SARS-CoV-2 and with   00:04:30.480 --> 00:04:37.740 viruses in general we have a much simpler problem,  right? We have a lot of sequence homology. Even if   00:04:37.740 --> 00:04:44.460 the virus is mutating, you know, significantly  across the world, every single viral sequence   00:04:44.460 --> 00:04:49.560 that we obtain is going to be almost identical to  the reference gene. It's not going to be exactly   00:04:49.560 --> 00:04:54.660 identical, but it's going to be almost identical.  So we actually are facing a much simpler   00:04:54.660 --> 00:05:01.560 computational problem which is multiple sequence  alignment of highly similar sequences. So how can   00:05:01.560 --> 00:05:07.260 we use that feature to speed up this analysis? We can do what's called an align-to-reference   00:05:07.260 --> 00:05:11.580 approach. So instead of just trying to align  everything with each other all at once, what   00:05:11.580 --> 00:05:17.040 we could do is individual pair-wise alignments  against a reference unit. So in this figure,   00:05:17.880 --> 00:05:23.040 the thicker green bar at the top represents the  reference to our scope 2 genome, and each of these   00:05:23.040 --> 00:05:28.920 other colored genomes represents a sequence  that I collect from the real world. I want to   00:05:28.920 --> 00:05:36.180 align each of these to the reference genome. What  I could do is just one by one by one by one I can   00:05:36.180 --> 00:05:42.840 independently align each of these genome sequences  against the reference genome, which I could do   00:05:42.840 --> 00:05:48.540 each of these fairly quickly and I can do massive  parallelization because each of these pairwise   00:05:48.540 --> 00:05:53.580 alignments to the reference can be done completely  independently. I can parallelize however many   00:05:53.580 --> 00:05:59.280 cores my computer has, I can throw that many at  this problem. Then, once I've competed all those   00:05:59.280 --> 00:06:04.740 pairwise lines to the reference, I could use  the wrap wrench genome - I can kind of use its   00:06:04.740 --> 00:06:10.800 anchors, its positions as anchors to create the  columns of my multiple sequence line. For example,   00:06:10.800 --> 00:06:16.320 maybe I'll start with the first position of the  reference genome and I'll see, okay, well, in   00:06:16.320 --> 00:06:21.660 the red sequence this is the letter that aligned  to that position. In the orange sequence, this   00:06:21.660 --> 00:06:25.560 is the letter. In the pink sequence, this is the  letter. In the blue sequence, this is the letter.   00:06:25.560 --> 00:06:30.360 And I can kind of merge all of those letters into  one column of my multiple sequence alignment.   00:06:31.200 --> 00:06:34.380 And I could do the same thing for the  second position of my reference genome,   00:06:34.380 --> 00:06:39.180 the same thing for the third position, fourth  position, all the way across. And kind of position   00:06:39.180 --> 00:06:47.340 by position by position I can build my multiple  sequence alignment. This idea, this is really   00:06:47.340 --> 00:06:53.400 good because it is massively parallelizable and  it scales linearly with the number of sequences   00:06:53.400 --> 00:06:59.460 rather than quadratically. So it has much better  scalability as well. Do we have to implement this   00:06:59.460 --> 00:07:03.960 approach from scratch? We actually don't. It turns out that if you kind of step back   00:07:03.960 --> 00:07:10.680 and think about this problem, this is really  equivalent, in a sense, to the long read mapping   00:07:10.680 --> 00:07:15.900 problem. Let's just kind of step back and rethink  what is the problem that we're tackling. Our   00:07:15.900 --> 00:07:21.180 input is a reference genome and a bunch of long  sequences that are very similar to the reference   00:07:21.180 --> 00:07:26.220 genome. Our output is an alignment of each of  those sequences against the reference genome.   00:07:26.220 --> 00:07:33.420 This is exactly the same computational problem as  mapping long reads. Instead of having to reinvent   00:07:33.420 --> 00:07:37.740 the wheel, we could just kind of build off of  all of these really advanced techniques that   00:07:37.740 --> 00:07:43.980 folks have built to solve the long read mapping  problem and just kind of apply it to this context.  00:07:45.600 --> 00:07:51.960 To that aim, I developed a tool called ViralMSA  and what it does is it just wraps around existing   00:07:51.960 --> 00:07:56.460 long read mappers to perform this reference-guided  multiple sequence alignment. It kind of treats   00:07:56.460 --> 00:08:01.920 each of those genomes that I've collected as  long reads and it treats the reference genome   00:08:01.920 --> 00:08:08.040 as a reference genome. It just calls that read  mapper - I wrap against a few different read   00:08:08.040 --> 00:08:14.100 mappers just to demonstrate flexibility - but I  mainly suggest people use Minimap2 for both speed   00:08:14.100 --> 00:08:19.260 and accuracy. Then, given those read mapping  results, I can then - or given the mapping   00:08:19.260 --> 00:08:22.860 results - I can then just kind of compile  them into a single multiple sequence line.  00:08:25.320 --> 00:08:31.680 So all you have to do to run ViralMSA is you just  give ViralMSA a reference genome and a bunch of   00:08:31.680 --> 00:08:36.360 sequences to align. It'll automatically  handle indexing the reference genome,   00:08:36.960 --> 00:08:41.160 maybe if you give it an accessory number it'll  handle downloading and indexing the reference   00:08:41.160 --> 00:08:45.720 genome. It'll handle all of the pre-processing  and all the downstream stuff and it'll just   00:08:45.720 --> 00:08:50.580 output - it'll call the read mapper, it'll merge  the results into the multiple SQL alignment,   00:08:50.580 --> 00:08:55.440 and it'll just output a single standard file  that is your multiple sequence alignment.  00:08:57.300 --> 00:09:02.340 How does it do against existing tools? We we  did a benchmark experiment where we compared   00:09:02.340 --> 00:09:09.960 the runtime of ViralMSA wrapping around Minimap  2 compared against Virulign, which is an existing   00:09:09.960 --> 00:09:14.820 align to reference approach but that just kind  of implements its own from scratch aligned to   00:09:14.820 --> 00:09:20.940 reference. We also compared against MAFFT which is  typically considered one of the most commonly used   00:09:20.940 --> 00:09:25.860 multiple sequence segment tools. In this plot, on  the horizontal axis I have number of sequences.   00:09:25.860 --> 00:09:30.900 On the vertical axis I have total execution  time in seconds. This was done on complete   00:09:30.900 --> 00:09:36.660 SARS-CoV-2 genome sequences, so genome length  roughly 29,000. As we can see, the blue line,   00:09:36.660 --> 00:09:43.560 which is ViralMSA, is orders of magnitude faster  than the existing tools. Compared to VIRULIGN,   00:09:43.560 --> 00:09:50.580 which is also scaling linearly, we're getting  - and by the way, this this plot is a log scale   00:09:50.580 --> 00:09:57.600 plot - so compared to VIRULIGN, we're roughly  like a thousand times faster-ish. And with MAFFT,   00:09:57.600 --> 00:10:01.860 we're not quite as much faster, but you can  see that because MAFFT scales quadratically,   00:10:01.860 --> 00:10:08.160 our speed up with respect to MAFFT is actually  increasing as time progresses. Even at just a   00:10:08.160 --> 00:10:13.320 thousand sequences. we hit roughly a thousand  times faster and that gap kind of increases.  00:10:14.820 --> 00:10:19.080 Now, you might be wondering, okay, well, fast is  good but what's the point if it doesn't give me   00:10:19.080 --> 00:10:25.380 good alignments? We also compared the accuracy.  What we did was we took the multiple sequence   00:10:25.380 --> 00:10:30.480 alignment computed by MAFFT, took the multiple  sequence alignment computed by ViralMSA on a bunch   00:10:30.480 --> 00:10:40.500 of hand-curated alignments from HIV, Ebola, and  I'm blanking on the third virus, but basically we   00:10:40.500 --> 00:10:44.400 took viruses that we had curated alignments from  the Los Alamo- oh actually, no - this plot is from   00:10:44.400 --> 00:10:50.340 just HIV-1. From the Los Alamos National Lab we  took their curated multiple sequence alignments   00:10:50.340 --> 00:10:54.720 and we use that as ground truth. Then we saw  how does mapped in ViralMSA compare against   00:10:54.720 --> 00:11:01.980 the ground truth curated multiple sequence  alignment. If we compute pairwise distances of   00:11:01.980 --> 00:11:07.260 the sequences that we get in our alignment, and  then we do a mantle test for accuracy - we find   00:11:07.260 --> 00:11:12.840 the correlation between our pairwise distances,  and the pairwise is calculated directly from the   00:11:12.840 --> 00:11:19.080 true multiple sequence alignment, we see that  the correlation is negligibly different. Here,   00:11:19.620 --> 00:11:24.300 you know, we get a correlation coefficient  like 0.994 for viral MSA compared to 0.997   00:11:25.140 --> 00:11:31.080 for pairwise distance calculation. Actually,  when we calculated phylogenies, the phylogenies   00:11:31.080 --> 00:11:36.120 inferred using the ViralMSA multiple sequence  elements are actually slightly higher topological   00:11:36.120 --> 00:11:41.820 accuracy than those estimated from MAFFT. So,  negligibly so, but still what we're showing is   00:11:41.820 --> 00:11:46.800 that these are kind of essentially equivalent in  terms of accuracy for all intents and purposes.  00:11:48.720 --> 00:11:53.220 The conclusion - ViralMSA is a tool that enables  rapid multiple sequence alignment of ultra large   00:11:53.220 --> 00:11:57.720 viral datasets. It's open source, you  can find it on GitHub, and you know,   00:11:57.720 --> 00:12:04.080 please consider using it in your viral analyzes. And acknowledgments - I want to thank Heng Li,   00:12:04.080 --> 00:12:09.480 he's the developer of Minimap2 and it's really  his expertise in developing Minimap2 that enables   00:12:09.480 --> 00:12:13.680 ViralMSA's speed and performance. I want to  thank the NSF for the grant that supports this   00:12:13.680 --> 00:12:18.780 project. And the the research was also supported  using Google Cloud platform research credits.  00:12:19.980 --> 00:12:23.880 So I’ll save time for any questions  or I'm happy to end it here.