Error processsing large PAF files from minimap2

[brettChapman]

Hi Erik

I've been running Seqwish with the output from minimap2 and edyeet for comparison, across multiple chromosomes and also testing with multiple whole genomes.

I've found the PAF files from minimap2 to be around twice the size of the PAF files from edyeet. I've tried processing the PAF file from minimap2 which is a 19GB PAF file (from an alignment of chromosome 5H from 4 genomes). However, it ends with a bad_alloc error due to inadequate memory. I'm going to try and run on a larger machine to see if I can get it to complete. The PAF file from Edyeet on the other can be processed with no problems.

I have also found that running edyeet with s50K across 4 whole genomes results in an error (likely a memory issue, as the job was outright killed), but when running with s100K it seems to continue the alignment with no problems, thus far (it's still running after 17 days). Regardless, I fear the resulting PAF file will still be too large to covert into a GFA file using the nodes from my cluster.

In regards to processing large PAF files, would it be conceivable to perform the alignment in parts (like we have discussed before), but instead of merging each PAF file and running the concatenated PAF file through Seqwish, we instead run Seqwish on each smaller PAF file, produce several GFA files, which we then each run through say Smoothxg, then run through a tool like GFAKluge (https://github.com/edawson/gfakluge), to merge all the GFA together, and then run Smoothxg on the final GFA at the end. Would breaking up the process into parts right up to producing the final GFA be a good way of reducing memory overhead? I know Smoothxg should probably be run on the whole GFA, but would iteratively running Smoothxg on several GFAs, and then running it on the final merged smoothed GFAs, ultimately reduce the resource overhead for each individual job? Or would this approach add more issues than it resolves? I see that you're one of the developers of GFAKluge, but it hasn't been updated in a while. Is GFAKluge something worth looking into to add to my pangenome workflow, to help reduce memory overhead? Thanks.

[ekg]

Hi Brett, What -n and how many cores did you give to edyeet? (I typically run -n less than the number of input genomes.) That seems longer than I would have expected, although I would expect it to take a while there. I need to add output that estimates time to completion. How much memory do you have on your system? I've never run out of memory for PAF reading and I have worked with 25 human assemblies at once. As for breaking the problem apart, that's a very good idea. There are several ways, including graph to graph alignment. I'll be exploring them.

…

On Mon, Oct 5, 2020, 04:53 Brett Chapman ***@***.***> wrote: Hi Erik I've been running Seqwish with the output from minimap2 and edyeet for comparison, across multiple chromosomes and also testing with multiple whole genomes. I've found the PAF files from minimap2 to be around twice the size of the PAF files from edyeet. I've tried processing the PAF file from minimap2 which is a 19GB PAF file (from an alignment of chromosome 5H from 4 genomes). However, it ends with a bad_alloc error due to inadequate memory. I'm going to try and run on a larger machine to see if I can get it to complete. The PAF file from Edyeet on the other can be processed with no problems. I have also found that running edyeet with s50K across 4 whole genomes results in an error (likely a memory issue, as the job was outright killed), but when running with s100K it seems to continue the alignment with no problems, thus far (it's still running after 17 days). Regardless, I fear the resulting PAF file will still be too large to covert into a GFA file using the nodes from my cluster. In regards to processing large PAF files, would it be conceivable to perform the alignment in parts (like we have discussed before), but instead of merging each PAF file and running the concatenated PAF file through Seqwish, we instead run Seqwish on each smaller PAF file, produce several GFA files, which we then each run through say Smoothxg, then run through a tool like GFAKluge (https://github.com/edawson/gfakluge), to merge all the GFA together, and then run Smoothxg on the final GFA at the end. Would breaking up the process into parts right up to producing the final GFA be a good way of reducing memory overhead? I know Smoothxg should probably be run on the whole GFA, but would iteratively running Smoothxg on several GFAs, and then running it on the final merged smoothed GFAs, ultimately reduce the resource overhead for each individual job? Or would this approach add more issues than it resolves? I see that you're one of the developers of GFAKluge, but it hasn't been updated in a while. Is GFAKluge something worth looking into to add to my pangenome workflow, to help reduce memory overhead? Thanks. — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#67>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AABDQEPEW7KEO32JY3Q5TLDSJEYJZANCNFSM4SEEA7KA> .

[brettChapman]

Hi Erik

I'm using -n 16 which is all the cores on one of my nodes. Do you suggest I change to -n 2 or -n 4? I'm currently testing with 4 genomes. I do plan to go all the way up to 20 genomes, however I'm doubting it will scale well on a single node, hence why I'm look at different ways to break the problem down and do as much processing across multiple nodes as I possibly can.

On each node I have 64GB RAM.

If you think it would be a good idea, I'll try working with GFAKluge and see how the final output looks. I'll try and rerun edyeet/minimap2 and seqwish with only a few cores and see if I can get them to complete that way first.

A graph to graph alignment sounds very interesting!

[ekg]

Hi Brett,

…

-t 16 would be 16 threads of processing. -n is the number of secondary alignments to keep. By default, edyeet will keep one best mapping (as in mashmap2 mapping), but setting -n tells it to keep up to that many. If you have a four genomes, you could set -n to 2, which would be 3 mappings for each sequence. If everything was collinear and similar, then this would mean that you basically map each sequence to the others. Leaving -t at 1 and running -n 16 might by why you've seen such a long processing time. Is that possible? As for merging GFAs, gfakluge won't do that. We need a process to map GFAs together. Thus far, there is nothing available, but there is interest in making one. I need to review what you wrote to try to understand it better. I suggest working to get the scalability of the mapping step up by reducing the number of secondary mappings to a lower level.

On Mon, Oct 5, 2020 at 8:42 AM Brett Chapman ***@***.***> wrote: Hi Erik I'm using -n 16 which is all the cores on one of my nodes. Do you suggest I change to -n 2 or -n 4? I'm currently testing with 4 genomes. I do plan to go all the way up to 20 genomes, however I'm doubting it will scale well on a single node, hence why I'm look at different ways to break the problem down and do as much processing across multiple nodes as I possibly can. On each node I have 64GB RAM. If you think it would be a good idea, I'll try working with GFAKluge and see how the final output looks. I'll try and rerun edyeet/minimap2 and seqwish with only a few cores and see if I can get them to complete that way first. A graph to graph alignment sounds very interesting! — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#67 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AABDQEN7PE56Z2MAOMCHE63SJFTGRANCNFSM4SEEA7KA> .

[brettChapman]

Hi Erik

Sorry, I said -n 16. I meant to say -t 16.

For the job which is taking over 2 weeks, I ran with these parameters:

edyeet -X -s 100000 -p 75 -n 20 -a 75 -k 11 -t 16 all_genomes.fasta all_genomes.fasta > pangenome.paf

Would you suggest I use -t 2 and -n 3 instead? Perhaps the -n value could be a function of the number of genomes (number of genomes minus 1), and -t the same but minus 2, and take the choice of values for these away from the user, to apply an optimal number to use?

In regards to GFAKluge, I saw one of its functions on it's Github page was to merge GFA, unless the GFA is a different version or format to the GFAs Seqwish generates?

[ekg]

I think the -k 11 is a potential problem. I would definitely leave that at the default (16). Also, -n 20 seems far too high for just 4 genomes. I would keep that in the order of the number of genomes. To understand the way things scale, you could run -n 0 for each pair of genomes, getting only the best mapping. Ideally, you'd work your way up from a single pair of chromosomes. Setting -s higher tends to reduce runtime in the mashmap2 step. There will be fewer potential mappings. The cost of scanning over the target doesn't change much with larger segment sizes. Did this edyeet job start to write its final output? I definitely need to add logging with expected time to completion. It is relatively easy to do.

…

On Mon, Oct 5, 2020 at 9:25 AM Brett Chapman ***@***.***> wrote: Hi Erik Sorry, I said -n 16. I meant to say -t 16. For the job which is taking over 2 weeks, I ran with these parameters: edyeet -X -s 100000 -p 75 -n 20 -a 75 -k 11 -t 16 all_genomes.fasta all_genomes.fasta > pangenome.paf Would you suggest I use -t 2 and -n 3 instead? Perhaps the -n value could be a function of the number of genomes (number of genomes minus 1), and -t the same but minus 2, and take the choice of values for these away from the user, to apply an optimal number to use? In regards to GFAKluge, I saw one of its functions on it's Github page was to merge GFA, unless the GFA is a different version or format to the GFAs Seqwish generates? — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#67 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AABDQELIA6DO46YQSN5ELP3SJFYHXANCNFSM4SEEA7KA> .

[brettChapman]

Thanks for the explanation. I'll keep -n to the total genome count, and keep -k to 16 and see how it runs again. I'll leave the current running edyeet job going just to see how long it takes.

I'll also try a larger -s value later to see how the run time changes.

In regards to the -s value with aligning smaller regions, what would be an appropriate value for multiple gene sequence regions? I tried running edyeet on an extracted gene region +/- 1000bp either side, and found -s larger than the longest sequence length to kill the job. I then ran -s with the length of the shortest sequence among all the gene sequences, and it ran fine. However there were a few gaps in the alignment compared to a minimap2 job I ran as well. Perhaps I could run -s at half the length of the shortest sequence? As a s100K value for whole genomes/chromosomes would be only a fraction of the roughly 600Mbps chromosome length, so a much shorter -s value (perhaps in the range of 500 to 1000) for aligning gene regions may be more appropriate (for my example gene I ran -s with a value of around 5000)

The currently running edyeet job still hasn't started outputting any final output.

A log with expected time to completion would definitely help a lot with planning and testing out parameters quickly.

[ekg]

You can use edyeet -M and -s set to your smallest input sequence size (or some number lower than that) to get merged mappings for your genes. If you're aligning shorter sequences, you might also want to look at wfmash. That has some memory limitations for sequences >50kb (things grow quadratically wrt the score of the alignment) but might give you better alignment quality because it's affine gapped.

…

On Mon, Oct 5, 2020 at 9:53 AM Brett Chapman ***@***.***> wrote: Thanks for the explanation. I'll keep -n to the total genome count, and keep -k to 16 and see how it runs again. I'll leave the current running edyeet job going just to see how long it takes. I'll also try a larger -s value later to see how the run time changes. In regards to the -s value with aligning smaller regions, what would be an appropriate value for multiple gene sequence regions? I tried running edyeet on an extracted gene region +/- 1000bp either side, and found -s larger than the longest sequence length to kill the job. I then ran -s with the length of the shortest sequence among all the gene sequences, and it ran fine. However there were a few gaps in the alignment compared to a minimap2 job I ran as well. Perhaps I could run -s at half the length of the shortest sequence? As a s100K value for whole genomes/chromosomes would be only a fraction of the roughly 600Mbps chromosome length, so a much shorter -s value (perhaps in the range of 500 to 1000) for aligning gene regions may be more appropriate (for my example gene I ran -s with a value of around 5000) The currently running edyeet job still hasn't started outputting any final output. A log with expected time to completion would definitely help a lot with planning and testing out parameters quickly. — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#67 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AABDQENWFOB6HTLT7JFBZO3SJF3PDANCNFSM4SEEA7KA> .

[brettChapman]

Thanks. I'll try edyeet with -M and also give wfmash a go and see how the results compare.

[ekg]

They both now write estimate time to completion during mapping and alignment.

…

On Tue, Oct 6, 2020, 06:38 Brett Chapman ***@***.***> wrote: Thanks. I'll try edyeet with -M and also give wfmash a go and see how the results compare. — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#67 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AABDQEIKBSGIJKIFFSFGZ33SJKNKXANCNFSM4SEEA7KA> .

[brettChapman]

Thanks! I'll try out the new builds.

In regards to running edyeet with a lower -n value and sticking with -k 16, I just ran edyeet overnight and it failed with 4 genomes:

I ran this:

edyeet -X -s 100000 -p 75 -n 3 -a 75 -k 16 -t 16 all_genomes.fasta all_genomes.fasta > pangenome.paf

The job received this error:

srun: error: node-5: task 0: Killed
srun: Force Terminated job step 2119.1

The job I started up over 2 weeks ago using a -n 20 and -k 11, is still running. I imagine it will eventually get killed off as well.

I'm guessing it's a memory problem. Could I try any other parameters, perhaps a much larger -s value, or do you think at this stage by adding more genomes the problem is simply inadequate memory, and the only way forward is by splitting up the alignments to 2 genomes per job and merging the PAF files is the only way to go? I've concluded as much with my minimap2 runs. I'm currently testing multiple alignment jobs and concatenation of PAF files with minimap2 now. Each genome is roughly 4Gbp.

[ekg]

It looks like a memory problem. Try to break it down as you've suggested. First look at memory consumption for a pair of chromosomes that you expect to be homologous. If that's doable, you can write a script to get the PAF for all homologous chromosome pairs. That'd be fine as input unless you want to see translocations.

…

On Wed, Oct 7, 2020, 03:44 Brett Chapman ***@***.***> wrote: Thanks! I'll try out the new builds. In regards to running edyeet with a lower -n value and sticking with -k 16, I just ran edyeet overnight and it failed with 4 genomes: I ran this: edyeet -X -s 100000 -p 75 -n 3 -a 75 -k 16 -t 16 all_genomes.fasta all_genomes.fasta > pangenome.paf The job received this error: srun: error: node-5: task 0: Killed srun: Force Terminated job step 2119.1 The job I started up over 2 weeks ago using a -n 20 and -k 11, is still running. I imagine it will eventually get killed off as well. I'm guessing it's a memory problem. Could I try any other parameters, perhaps a much larger -s value, or do you think at this stage by adding more genomes the problem is simply inadequate memory, and the only way forward is by splitting up the alignments to 2 genomes per job and merging the PAF files is the only way to go? I've concluded as much with my minimap2 runs. I'm currently testing multiple alignment jobs and concatenation of PAF files with minimap2 now. Each genome is roughly 4Gbp. — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#67 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AABDQEK6TUIY42B6QHQNMPDSJPBWJANCNFSM4SEEA7KA> .

[ekg]

I would also start backing off on the parameters. What level of divergence do you expect between the sequences? Set -p and -a a little below that. Using a shorter segment length is probably not a bad idea.

[brettChapman]

Ok, thanks I'll try with lower parameter values. Perhaps -a 50 -p 50 and -s 50000 or -s 75000.

The divergence between varieties is quite variable. Some varieties are recently diverged being crossed by breeders over the last decade or two, or perhaps more recently, while others are wild-type varieties and may have diverged <10K years ago.

[ekg]

50% divergence is too high. It actually maxes out at 70%. But that's the limit for what we can meaningfully align with this strategy. If these are all varieties of the same species then I would actually keep -p and -a lower. Try setting them close to your expected pairwise diversity. In this case 90 or 95 might be good. This will make the alignment much faster and dramatically less memory intensive.

…

On Thu, Oct 8, 2020, 09:58 Brett Chapman ***@***.***> wrote: Ok, thanks I'll try with lower parameter values. Perhaps -a 50 -p 50 and -s 50000 or -s 75000. The divergence between varieties is quite variable. Some varieties are recently diverged being crossed by breeders over the last decade or two, or perhaps more recently, while others are wild-type varieties and may have diverged <10K years ago. — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#67 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AABDQEPFZGHNGBD247WA67TSJVWKNANCNFSM4SEEA7KA> .

[brettChapman]

Ok, so perhaps -a 90 -p 90 -s 50000 ? I'll submit another job with these values.