Multi-threaded (faster) SAMtools

[nilshomer]

I have been working on speeding up reading and writing within SAMtools by creating a multi-threaded block gzip reader and writer. I have an alpha version working. I would appreciate some feedback and testing, just don't use it for production systems yet. Thank-you!

http://github.com/nh13/samtools/tree/pbgzip

NB: I would be happy describe the implementation, and collaborate to get this into Picard too.

[Richard Finney]

Any benchmarks?

[nilshomer]

Copied here from http://sourceforge.net/mailarchive/m...sg_id=28915492

I am working on benchmarking the samtools commands today, and will post back.

A 4GB SAM file was used on a dual-hex-core (12 cores) computer. I
benchmarked compression then decompression, making sure the resulting files
were the same. Decompression seems to be limited by IO.

Name Compression Time Decompression Time
bgzip 485.64 39.93
pbgzip -n 1 481.57 40.02
pbgzip -n 2 240.85 41.03
pbgzip -n 4 122.05 41.79
pbgzip -n 8 63.17 41.17
pbgzip -n 12 43.12 41.65
pbgzip -n 16 39.59 41.48
pbgzip -n 20 37.03 42.41
pbgzip -n 24 34.90 47.24

[nilshomer]

Updated #s on a few commands:
Command samtools psamtools
view BAM 29.45 19.2
view -b BAM 207.51 19.36
view -S SAM 44.89 44.43
view -Sb SAM 222.64 32.62
sort 206.32 25.17
mpileup 6574.2 7252.08
depth 17.64 7.47
index 11.96 1.93
flagstat 11.73 1.73
calmd -b 209.25 22.86
rmdup -s 154.88 22.08
reheader 0.76 0.74
cat 1.54 1.37

[Richard Finney]

Looks good!
question
1) why is mpileup slower?

[nilshomer]

Working on it. I am doing this in my free time, so having one perform worse isn't that bad so far.

[krobison]

Really cool!!

Do you have benchmarks for retrieving specific reads for a region? For mpileup of a specific region or a list of targets?

Any idea if this will work with the Bio:B::Sam perl module (which must be linked in to samtools)

What are the prospects for merging this with the main samtools development?

[nilshomer]

The seeks are just as fast, so no speedup/slowdown on seeking, but then there should speedup reading from that point on, assuming there are at least a basal number of reads in the region (otherwise there is no work to be done). For mpileup, it doesn't process the regions in parallel, if that is what you were implying.

I posted to the samtools list with response, so I have no hypothesis as to the inclusion of this (of course it needs more testing first). It generally is difficult to get things included there. I have more hope for Picard.

Pysam and the SAM perl module should not notice the difference in the API, though there is no good mechanism yet for determining the # of threads to use (it autodetects the # of cores).

[krobison]

I also see the sort command now gives an option to pick an algorithm. What a blast to the past?

Any heuristics on what algorithm might perform better in what setting?

And why no bubble sort option :-)

[nilshomer]

I also see the sort command now gives an option to pick an algorithm. What a blast to the past?

Any heuristics on what algorithm might perform better in what setting?

And why no bubble sort option :-)

I just used Heng's ksort.h library. I like introsort, but mergesort is by default in the original samtools.

I have also been toying with multi-threaded sort, which sort of works in the new version, except I didn't take time to do a proper multi-way merge (one implementation requires the calculation of evenly spaced pivots). Maybe wait a few more weekends.

[adaptivegenome]

I have been working on speeding up reading and writing within SAMtools by creating a multi-threaded block gzip reader and writer. I have an alpha version working. I would appreciate some feedback and testing, just don't use it for production systems yet. Thank-you!

http://github.com/nh13/samtools/tree/pbgzip

NB: I would be happy describe the implementation, and collaborate to get this into Picard too.

So are you saying you made a parallelized version of BZIP2? We have also been playing around with this. We parallelized the compression and decompression steps in the read/write functions of samtools for a local realignment tool we built.

I would love to learn more about what you are doing as I would hate to duplicate anything you are going to already do!

[colindaven]

As far as parallel (g)zip goes pigz works wonders : http://zlib.net/pigz/

[adaptivegenome]

PIGZ is very very fast however it produces file sizes that are much larger than BZIP2. Is this your experience as well?

It would be really nice to be able to simply parallelize BZIP2. We have tried to do this a little bit but certainly don't have completed product yet.

[lh3]

Firstly I greatly appreciate and strongly support Nils effort in multithreading samtools. The change is likely to be merged to samtools.

Re sorting algorithm: samtools sort does stable sorting (i.e. preserving the relative order of records having the same coordinate). In some rare/non-typical use cases, this feature is useful. Merge sort is stable. Introsort is not.

Re pigz: someone told me on biostar that pigz is not very scalable with many cores. If this is true (I have not tried), this must be because the gzip format has long range dependencies. bzip2 and bgzip are much easier to parallelize and probably more scalable. In addition, bzip2 has a parallel version pbzip2 which the same person told me scales very well with the number of CPU cores.

Re bzip2: I have argued a couple times here (years ago) and also on the samtools list that the key reason samtools uses gzip instead of bzip2 is because gzip is 5-10X faster on decompression. With bzip2, most samtools command will be 2-10 times slower. I think for huge data sets that need to be read frequently, gzip is always preferred over bzip2.

[adaptivegenome]

I think it is worth figuring out the best way to compress/decompress. Our nodes have 64 cores so I will do some tests and see how BZIP2 and GZIP scale. I'll post what I find on this thread.

In the meantime a quick internet search turned up this:

http://nerdbynature.de/s9y/?251

Firstly I greatly appreciate and strongly support Nils effort in multithreading samtools. The change is likely to be merged to samtools.

Re sorting algorithm: samtools sort does stable sorting (i.e. preserving the relative order of records having the same coordinate). In some rare/non-typical use cases, this feature is useful. Merge sort is stable. Introsort is not.

Re pigz: someone told me on biostar that pigz is not very scalable with many cores. If this is true (I have not tried), this must be because the gzip format has long range dependencies. bzip2 and bgzip are much easier to parallelize and probably more scalable. In addition, bzip2 has a parallel version pbzip2 which the same person told me scales very well with the number of CPU cores.

Re bzip2: I have argued a couple times here (years ago) and also on the samtools list that the key reason samtools uses gzip instead of bzip2 is because gzip is 5-10X faster on decompression. With bzip2, most samtools command will be 2-10 times slower. I think for huge data sets that need to be read frequently, gzip is always preferred over bzip2.