Minimap2: parallelization with -sr

[brentp]

with short reads, the default parallelization isn't keeping more than 3 of 4 cores busy most of the time.

It seems to improve quite a bit by bumping the batch size. I tried doing this from the command-line with K, but I think maybe it gets overridden by the -x sr opts?

If I modify the source as below, it keeps more cores (tested on 12 cpus) busy for a higher percentage of time.

diff --git a/map.c b/map.c
index ac66e81..2273d32 100644
--- a/map.c
+++ b/map.c
@@ -98,7 +98,7 @@ int mm_set_opt(const char *preset, mm_idxopt_t *io, mm_mapopt_t *mo)
                mo->best_n = 20;
                mo->mid_occ = 1000;
                mo->max_occ = 5000;
-               mo->mini_batch_size = 50000000;
+               mo->mini_batch_size = 250000000;
        } else if (strcmp(preset, "splice") == 0 || strcmp(preset, "cdna") == 0) {
                io->is_hpc = 0, io->k = 15, io->w = 5;
                mo->flag |= MM_F_SPLICE | MM_F_SPLICE_FOR | MM_F_SPLICE_REV;

[lh3]

Which version are you using? What is the exact command line?

[brentp]

command is:

minimap2 -ax sr -t $cpus -K 13900M -R $rg $mmi <(bzcat $r1) <(bzcat $r2)

(the -K has no effect) $mmi is g1k_v37_decoy.fa

The bzcat is slow, but not the bottleneck.

I was using commit ce06188

[lh3]

Non-functioning -K is a bug. I will fix.

There is also another way to improve parallelization efficiency at the cost of higher peak RAM. I will experiment a bit. Thanks!

[lh3]

I have fixed the -K issue via 2801ed9.

In theory, a large -K only helps when one read in a minibatch takes longer to map than the rest of reads (this happens when you map very long contigs). For short reads, each read takes little time to map. A large -K shouldn't help. On my machine, this is the case. A large -K increases the time on loading the first and writing the last minimatch. The total wall clock time increases although during mapping, minimap2 appears to run at full capacity.

7c555f9 implements another way to improve multi-threading performance. The change to minimap2.1 explains what this new -2 option is doing. This option is now applied by default when you use -x sr.

Also, I guess you might be piping the minimap2 output to gzip. In my test, running minimap2 on 8 threads will nearly saturate the bandwidth of gzip -3. This means when you pipe SAM stream to gzip for on-the-fly compression, you won't reduce wall-clock time further unless you choose level-1 or level-2 gzip compression or use pigz to keep up with the throughput of minimap2.

[brentp]

This helps quite a lot. It is still not fully using 12 cpus. I am piping through samblaster into samtools sort. Neither samblaster nor samtools sort are using much CPU--though samblaster seems to be at about 15% which is higher than I've ever seen it with bwa mem, even with 64 cpus.

Also, after your recent change, the bzcat processes are at ~65% CPU so it's likely approaching IO limits.

thanks!

[lh3]

Further parallel efficiency can be achieved with SNAP-like strategy: integrating markduping and sorting into the mapper to reduce unnecessary interprocess I/O and to fine schedule all threads. That would be another project.

I am closing this for now. Thanks for the report!

[brentp]

by the way, for about every 20 or so reports to STDERR of [M::worker_pipeline ...] mapped 495050 sequences

I see: [W::sam_parse1] mapped query cannot have zero coordinate; treated as unmapped

... ( when piping to samtools sort )

[brentp]

seems to exclusively happen with hs37d5 chromosome

[lh3]

This should be a bug. I have not run enough reads to catch it. Did samtools say which reads had zero coordinate? If not, I will change the code later to print violating reads to stderr. I need examples to debug.

[brentp]

it doesn't. but I changed it. getting reads now. a lot of grepping...

[brentp]

Here's on from phix:
r1:

@HWI-ST807:461:C2P0JACXX:6:1314:5605:18097 1:N:0:TGACCA
TTCGATAAAAATGATTGGCGTATACAACCTGCAGAGTTTTATCGCTTCCATGACGCAGAAGTTAACACTTTCGGATATTTCTGATGAGTCGAAAAATTATC
+
@@@DDDDA>BF>FF9@4EFD?F@+AFCG;CDGG*1?*:?D499;?DDDBF@FHIIIIG:3?=?)77;)7.7?=B?=,(;>;>AD#################

r2

@HWI-ST807:461:C2P0JACXX:6:1314:5605:18097 2:N:0:TGACCA
TGAAAATCCCAACGAAACACTCTCTACCATTTTCAAATTGTAACTCATATATAAACACCTCCTTGACAATCGTGCTAAACACATTAACACACTTCTACTCA
+
#####################################################################################################

aligned to g1k_v37_decoy.fa. outputs:

HWI-ST807:461:C2P0JACXX:6:1314:5605:18097       73      phix    1       44      33S68M  =       1       0       TTCGATAAAAATGATTGGCGTATACAACCTGCAGAGTTTTATCGCTTCCATGACGCAGAAGTTAACACTTTCGGATATTTCTGATGAGTCGAAAAATTATC    @@@DDDDA>BF>FF9@4EFD?F@+AFCG;CDGG*1?*:?D499;?DDDBF@FHIIIIG:3?=?)77;)7.7?=B?=,(;>;>AD#################    tp:A:P  cm:i:7  s1:i:53 s2:i:0  NM:i:0  ms:i:136        AS:i:136        nn:i:0
HWI-ST807:461:C2P0JACXX:6:1314:5605:18097       133     phix    0       0       *       =       1       0       TGAAAATCCCAACGAAACACTCTCTACCATTTTCAAATTGTAACTCATATATAAACACCTCCTTGACAATCGTGCTAAACACATTAACACACTTCTACTCA    #####################################################################################################

[lh3]

Thanks a lot. Fixed via dea3b60.

[brentp]

thank you.