[dp05yk]

Actually YEG, I did find a bug. Thanks for pointing this out to me. It was assigning the processor rank AFTER determining the filename. I guess every system behaves differently, so yours was assigning a rank of -1, hence the additional dash.

I hadn't caught this because I did most if not all of my testing with the sampe command as it seemed to be more popular.

I'll be uploading the latest revision to the sourceforge page today, thanks for the input!

[dp05yk]

This may be a small bug. I had to rename *.sai files for the above command to work. The files need to have an extra '-'. So [prefix]-0.sai needs to be named [prefix]--0.sai and so on for every file made with pBWA align.

Here's the pBWA align command I used :

mpirun -np 3 -hostfile hostfile pBWA aln -f a1 -t 24 ~/hg18/hg18 all.fq 29424134

That's... really strange. I just checked the code (for both revisions 21 and 30) and it seems like it should be functioning properly... both bwase and bwape take the entered prefix and concatenate "-%d.sai", where %d = processor rank.

[YEG]

I should probably just have showed you:

./pBWA samse -f out ~/hg18/hg18 a1 all.fq 29424134

This may be a small bug. I had to rename *.sai files for the above command to work. The files need to have an extra '-'. So [prefix]-0.sai needs to be named [prefix]--0.sai and so on for every file made with pBWA align.

Here's the pBWA align command I used :

mpirun -np 3 -hostfile hostfile pBWA aln -f a1 -t 24 ~/hg18/hg18 all.fq 29424134

[dp05yk]

Hi YEG,

You want to input the .sai prefix (a1)... Then pBWA will align every .sai file that you have with that prefix! Also, you just want to specify "out" as your -f parameter as pBWA will add the rank and .SAMs to the output files.

Also, you may want to use revision 30, always best to stay current! :-)

I should probably just have showed you:

./pBWA samse -f out ~/hg18/hg18 a1 all.fq 29424134


And that will align all of your .sai files at the same time!

Cheers,
Darren

[dp05yk]

Hi YEG,

You want to input the .sai prefix (a1)... Then pBWA will align every .sai file that you have with that prefix! Also, you just want to specify "out" as your -f parameter as pBWA will add the rank and .SAMs to the output files.

Also, you may want to use revision 30, always best to stay current! :-)

[YEG]

How does pBWI work with single-end reads? It's not really clear from the tutorial.
I aligned single-end reads using 3 processors. Now I have 3 *.sai files: a1-0.sai ... a1-2.sai. To get the sam file I tried:

pBWA samse -f out.sam ~/hg18/hg18 a1-0.sai all.fq 29424134

pBWA crashes and produces the following:

[bwa_sai2sam_se_core] fail to open file 'a1-0.sai--1.sai'. Abort!
[bart:29010] *** Process received signal ***
[bart:29010] Signal: Aborted (6)
[bart:29010] Signal code:  (-6)
[bart:29010] [ 0] /lib64/libpthread.so.0 [0x7f8f5393bc00]
[bart:29010] [ 1] /lib64/libc.so.6(gsignal+0x35) [0x7f8f528984e5]
[bart:29010] [ 2] /lib64/libc.so.6(abort+0x180) [0x7f8f528999b0]
[bart:29010] [ 3] pBWA [0x405309]
[bart:29010] [ 4] pBWA(bwa_sai2sam_se_core+0xca) [0x41597a]
[bart:29010] [ 5] pBWA(bwa_sai2sam_se+0x14a) [0x415e5a]
[bart:29010] [ 6] pBWA(main+0xe3) [0x427263]
[bart:29010] [ 7] /lib64/libc.so.6(__libc_start_main+0xfd) [0x7f8f52884a7d]
[bart:29010] [ 8] pBWA [0x404f69]
[bart:29010] *** End of error message ***
Aborted (core dumped)

I tried the same command with 'regular' bwa (minus the last argument) and it executed without a problem. What am I missing?

I am using 0.5.9-r21-MPI.

[bioinfosm]

We are seeing more and more that re-alignment is an amazing benefit but terribly slow.. it makes the regular bwa alignment seem so fast... I hope there are solutions in the works for re-alignment aspect, where one needs to take all reads in a window and cannot arbitrarily split and parallelize..

[seb567]

Hi Sebastien,



A mod value is the remainder after performing integer division. For instance, 7 % 3 = 1, since 3*2 = 7 - 1. So the way we handle sequence distribution for threading is:

Loop i = 1 to (num_seqs)

if (i % num_threads) = thread_id then process
else skip

End loop

Since we are performing the modulus on the loop counter with the number of threads, we are guaranteed to cycle through the numbers 0...num_threads for each consecutive sequence. This ensures that the sequences are evenly divided, and it also ensures that no threads will be competing, since thread i only processes sequences (i, i+num_threads, i+(2*num_threads), etc.

Does that make sense? Previously, threads would essentially fight over sequence distribution by locking and "reserving" sequences for processing. This is responsible for the 20% efficiency difference.

I know what a modulo is.

In Ray, a de novo assembler, we were using this approach to split sequences across MPI ranks. Now, we just do a simple partition on the sequences.

Given N sequences (which can be in many files, of course) and M MPI ranks, MPI rank 0 takes sequences 0 to (N/M)-1, MPI rank 1 takes sequences (N/M) to 2*(N/M)-1, and so on. Finally, the MPI rank M-1 (the last one) also takes the remaining N%M sequences.

The partition-wise approach has the advantage that each MPI rank knows where to start and where to end.

This function in bwaseqio.c performs the input reads file indexing. In order to distribute reads evenly over processes, each process receives a contiguous block of reads. However, with paired reads, we cannot assume that both reads files will be of exactly the same length (especially when dealing with SOLiD reads), so we need to index the files to find the start and end location of each contiguous block of reads. This is a one-processor job, and processor 0 essentially scans the file, marking the start and end locations of an evenly distributed reads block. Once it finds these positions, it sends them to processor i and this becomes processor i's block of input reads. The reason these are 8-byte numbers is because some input reads files are extremely large and are larger than 2^32 bytes.

Regardless, I think you could enhance your already-enhanced approach using message aggregation.

Example with 4 MPI ranks and 9 integers so send:

Without message aggregation

Rank 0 sends value 0 to Rank 1
Rank 0 sends value 1 to Rank 2
Rank 0 sends value 2 to Rank 3
Rank 0 sends value 3 to Rank 1
Rank 0 sends value 4 to Rank 2
Rank 0 sends value 5 to Rank 3
Rank 0 sends value 6 to Rank 1
Rank 0 sends value 7 to Rank 2
Rank 0 sends value 8 to Rank 3

(9 messages)

With message aggregation


Rank 0 sends values 0,3,6 to Rank 1
Rank 0 sends values 1,4,7 to Rank 2
Rank 0 sends values 2,5,8 to Rank 3

(3 messages)

In this toy example, agglomerated messages contains 3 values.

You can bundle 500 8-byte integers (4000 bytes) in a 4096-byte message, assuming the the envelope is at most 96 bytes.

So, in your case, agglomerated messages would contain 500 values and you would divide your number of sent messages by 500, which is good given that transiting a message between two MPI ranks that are not on the same computer is costly.

Sébastien http://Boisvert.info

[dp05yk]

Does threading improve over the "embarrassingly parallel" style of just splitting the input and running on num_available_cpus instances? I think any speedups in BWA would be with 1) aggressive optimization 2) re-workiing code to keep more data/code in close cache.

It's interesting you mention this, as I originally thought pBWA with x processors would be slower than BWA with x threads, however I found them to be competitive, and even found pBWA to be faster than BWA in most instances.

However, once I implemented my multithreading improvement into BWA, BWA with x threads outperforms pBWA with x threads. The advantage provided by pBWA is that you can run it on multiple processors, whereas multithreaded BWA only runs on as many cores as you possess on a single processor (or x2 if you have hyperthreading). Another advantage provided with pBWA is that you can combine parallelism and multithreading,

ie. run pBWA with 20 processes, each process on its own node, and each process spawning x threads, where each node is an x-core processor. I've been able to run pBWA with 2000+ cores and align 500million+ SOLiD reads in under an hour.

[dp05yk]

Hi Sebastien,

What exactly is a loop counter's modulo value ?

A mod value is the remainder after performing integer division. For instance, 7 % 3 = 1, since 3*2 = 7 - 1. So the way we handle sequence distribution for threading is:

Loop i = 1 to (num_seqs)

if (i % num_threads) = thread_id then process
else skip

End loop

Since we are performing the modulus on the loop counter with the number of threads, we are guaranteed to cycle through the numbers 0...num_threads for each consecutive sequence. This ensures that the sequences are evenly divided, and it also ensures that no threads will be competing, since thread i only processes sequences (i, i+num_threads, i+(2*num_threads), etc.

Does that make sense? Previously, threads would essentially fight over sequence distribution by locking and "reserving" sequences for processing. This is responsible for the 20% efficiency difference.

Correct me if I am wrong, but you use MPI_Send/MPI_Recv in bwaseqio.c because you don't want rank 0 to receive anything, but you you want all the other ranks in MPI_COMM_WORLD to receive data prepared by rank 0.

What exactly are these 8-byte numbers that you are sending there ?

For any MPI rank greater than 0, the MPI rank 0 will send numerous 8-byte values individually. I can see there that you could manually group them and send them in agglomerates. Messages are usually sent eagerly when they are at most 4096 bytes in MPICH2 and Open-MPI.

Let us assume that the envelope is at most 96 bytes. Then, the data size will be at least 4000 bytes -- or 500 8-byte integers. So, since you will be manually grouping your 8-byte communications into a few larger 4000-byte communications, the input/output will be at most 500 times faster considering the overhead of sending individual small messages.

This function in bwaseqio.c performs the input reads file indexing. In order to distribute reads evenly over processes, each process receives a contiguous block of reads. However, with paired reads, we cannot assume that both reads files will be of exactly the same length (especially when dealing with SOLiD reads), so we need to index the files to find the start and end location of each contiguous block of reads. This is a one-processor job, and processor 0 essentially scans the file, marking the start and end locations of an evenly distributed reads block. Once it finds these positions, it sends them to processor i and this becomes processor i's block of input reads. The reason these are 8-byte numbers is because some input reads files are extremely large and are larger than 2^32 bytes.

[Richard Finney]

Does threading improve over the "embarrassingly parallel" style of just splitting the input and running on num_available_cpus instances? I think any speedups in BWA would be with 1) aggressive optimization 2) re-workiing code to keep more data/code in close cache.

[seb567]

Hello dp05yk,

I am looking forward for your published paper about this work.

As for the improvements to multithreading, I noticed that BWA introduces too much thread competition, which isn't a problem for small amounts of threads, but when I was going up to 24 threads I noticed that pBWA was actually faster for 24 parallel processes than BWA was for 24 threads which surprised me. I changed the way BWA handles multithreading by basing it purely on a loop counter's mod value and removed all the sequence locking, etc, and multithreading improved by ~20% for higher amounts of threads.

Quite interesting !

What exactly is a loop counter's modulo value ?

Do you suspect that removing the locks was THE game changer that allowed a significant 20% improvement ?

I think multi-threading samse/sampe is a huge contribution.

I also think that it is a huge contribution !



Correct me if I am wrong, but you use MPI_Send/MPI_Recv in bwaseqio.c because you don't want rank 0 to receive anything, but you you want all the other ranks in MPI_COMM_WORLD to receive data prepared by rank 0.

What exactly are these 8-byte numbers that you are sending there ?

For any MPI rank greater than 0, the MPI rank 0 will send numerous 8-byte values individually. I can see there that you could manually group them and send them in agglomerates. Messages are usually sent eagerly when they are at most 4096 bytes in MPICH2 and Open-MPI.

Let us assume that the envelope is at most 96 bytes. Then, the data size will be at least 4000 bytes -- or 500 8-byte integers. So, since you will be manually grouping your 8-byte communications into a few larger 4000-byte communications, the input/output will be at most 500 times faster considering the overhead of sending individual small messages.






Awesome work !

Cordially,

Sébastien

[seb567]

We started out with openMPI and found many shortcoming and switched to MPICH2 which works much better based on our benchmarking tests.

I am also curious to hear what shortcomings were encountered.

MPICH2 and Open-MPI are both compliant with the MPI 2.2 standard !

http://www.mpi-forum.org/docs/docs.html

[dp05yk]

I think we had novoalignMPI out before pBWA so it wouldnt be the first parallel alignment tool. We started out with openMPI and found many shortcoming and switched to MPICH2 which works much better based on our benchmarking tests.

Thanks for the info! I'll edit that out of my original post.

By the way, what shortcomings did you find with OpenMPI? I'm curious as to any problems you ran into, as I may have had the same.

[zee]

I think we had novoalignMPI out before pBWA so it wouldnt be the first parallel alignment tool. We started out with openMPI and found many shortcoming and switched to MPICH2 which works much better based on our benchmarking tests.

BTW nice work and it would be great to see more projects like this for SNP calling, Indel detection, SNV detection, coverage calculation, etc.