Currently, all operations are not optimized at all and use huge amount of FFI calls. I believe it can be done better.

Not a bug, but a bit confusing.
I see flagstat uses -n for threads, but sort uses -t

Every alignment record should have properties for easy access to flag bits, with human-readable names, and they are to be taken into account in algorithms, at least whether the read is mapped or not.

Testing Value code...

object.Exception@utils/range.d(21): Amount of elements to prefetch must be positive

./run_unittests(pure @safe bool std.exception.enforce!(bool, "utils/range.d", 21).enforce(bool, lazy const(char)[])+0x2b) [0x81642df]

At the moment, type conversions are very restrictive. If Value holds integer, and I want to get integer representation as a string, I have to use first to!int, and then to!string. Moreover, the user should be able to convert value holding byte to short, int, long, string, float, and double. And so on. The user really should not care about the actual type as long as he converts it to a type it can be converted to.

That would make working with tags way easier. I want to be able to do that, for instance:

Value v = 25;
float f = to!float(v); // determine, that value holds int, then convert this int to float
string s = to!string(v); // the same, but convert to string

to!float(to!int(v)) is not what the user wants to type, that's for sure.

I realise this is a dev tree, even so, I would like to be able to run the benchmark :)

verona:~/tmp/BAMread/benchmarks$ make
dmd ../bamfile.d ../utils/range.d ../chunkinputstream.d ../bgzfrange.d ../samheader.d ../reference.d ../alignment.d ../tagstorage.d ../tagvalue.d ../utils/switchendianness.d ../virtualoffset.d ../randomaccessmanager.d ../bai/read.d ../bai/bin.d ../bai/chunk.d ../bai/utils/algo.d ../utils/algo.d benchmarks.d -ofrun_benchmarks -O -release -inline -version=serial
../randomaccessmanager.d(14): Error: module memoize is in file 'utils/memoize.d' which cannot be read
import path[0] = /usr/include/i386-linux-gnu/dmd/phobos
import path[1] = /usr/include/i386-linux-gnu/dmd/druntime/import
make: *** [serial] Error 1

Hi,

On occasion I see this error "sambamba-view: not enough data in stream".
It does not happen always and if I re_run the command it usually works, but it seems it gets this on a regular basis...

We run this on an isilon system that gets heavily used, so maybe there is some timeout setting or something so that if it does not get the data as fast as it wants, it barfs?

Thanks

Thon

This can be done relatively easily once there's support for constructing/modifying alignments.

It is possible to make a simple query parser using Pegged so that users can write queries to get needed reads.

E.g., alignments.select("refseq:(chr1 OR chr2) AND region:[10000 TO 15000] AND quality > 50")

The idea is that it can be useful for web apps, and also can allow easy interfacing with dynamic languages, without passing sophisticated objects through FFI.

outdated quite a bit

I got stuck with this error:
export LD_LIBRARY_PATH=/home/avilella/src/dmd2/linux/lib64
export PATH=$PATH:/home/avilella/src/dmd2/linux/bin64
cd ~/src/sambamba/latest/sambamba
make clean
make
mkdir -p build/
rdmd --compiler=dmd --build-only -IBioD -version=development -ofbuild/sambamba main.d
sambamba/merge.d(211): Error: non-shared method bio.bam.utils.samheadermerger.SamHeaderMerger.this is not callable using a shared object
sambamba/merge.d(211): Error: no constructor for SamHeaderMerger
make: *** [all] Error 1

Many bioinformaticians have never used D, and it can be problematic for them to assess features of the new library. The script should work on all platforms.

With one of files from 1000 Genomes, using index generated by Samtools, when I specify '20' as the region in sambamba view, it throws an exception, trying to read some non-existent BGZF block. However, if '20:1-100000000' is specified in command-line, everything is OK.

(Perhaps, has to do with some erroneous type conversions between size_t and uint... will look into that tomorrow)

There should be info about Biogems & Bioruby as well as description how to install and use the library. Having some code examples would be nice.

Big-endian system support is absolutely untested. QEMU can emulate big-endianness when using PowerPC architecture.

Currently, when constructing alignments, memory is copied from decompressed blocks. That shouldn't be the case, slices should be used instead. That will save about 10% of iteration time, also reducing GC contention and making multithreaded version faster. (By the way, using slices instead of memcpy can really become a killer feature. That's simply impossible in 'traditional' languages.)

Creating BAI is a must, maybe also consider alternative indexing schemes (creating nested containment lists, for example).

I ran this on an Ubuntu 64-bit VM with 4-cores (running on a Windows 7 64-bit laptop). It seems like running on more cores actually slows down the view command in this instance.

I repeated the experiment at a cluster farm with isilon filesystem on a smaller file, and I also get slower results when using higher number of cores.

I also tried it on a node that has a solid state drive partition by copying the bam and .bai files there, and I get similar results both in time taken and taking slightly more time on higher number of cores than lower.

Any ideas?

avilella@ubuntu64:~/src/sambamba/latest/sambamba$ for t in seq 4 -1 1; do echo $t===; time ~/src/sambamba/latest/sambamba/build/sambamba view -t $t NA12878.bam 1>/dev/null; done
4===
[info] file size: 4366637089
[info] file size: 4366637089
[info][BGZF range] EOF, current offset is 4366637089
[info][chunkinputstream] len == 0, start offset = 4366637061, end offset = 4366637089
[info][chunkinputstream] len == 0, start offset = 4366637089, end offset = 4366637089
[info][bamRead range] empty, read 0 bytes, expected 4

real 26m28.151s
user 20m1.075s
sys 0m31.534s
3===
[info] file size: 4366637089
[info] file size: 4366637089
[info][BGZF range] EOF, current offset is 4366637089
[info][chunkinputstream] len == 0, start offset = 4366637061, end offset = 4366637089
[info][chunkinputstream] len == 0, start offset = 4366637089, end offset = 4366637089
[info][bamRead range] empty, read 0 bytes, expected 4

real 14m37.343s
user 17m41.930s
sys 0m27.646s
2===
[info] file size: 4366637089
[info] file size: 4366637089
[info][BGZF range] EOF, current offset is 4366637089
[info][chunkinputstream] len == 0, start offset = 4366637061, end offset = 4366637089
[info][chunkinputstream] len == 0, start offset = 4366637089, end offset = 4366637089
[info][bamRead range] empty, read 0 bytes, expected 4

real 13m49.845s
user 16m7.328s
sys 0m21.113s
1===
[info] file size: 4366637089
[info] file size: 4366637089
[info][BGZF range] EOF, current offset is 4366637089
[info][chunkinputstream] len == 0, start offset = 4366637061, end offset = 4366637089
[info][chunkinputstream] len == 0, start offset = 4366637089, end offset = 4366637089
[info][bamRead range] empty, read 0 bytes, expected 4

real 10m40.287s
user 12m16.670s
sys 0m15.185s

I repeated the experiment on a smaller file at a cluster farm with isilon filesystem, and this is what I get:
8===
[info] file size: 1647602665
[info] file size: 1647602665
[info][BGZF range] EOF, current offset is 1647602665
[info][chunkinputstream] len == 0, start offset = 1647602637, end offset = 1647602665
[info][chunkinputstream] len == 0, start offset = 1647602665, end offset = 1647602665
[info][bamRead range] empty, read 0 bytes, expected 4

real 1m34.347s
user 1m56.074s
sys 0m4.066s
7===
[info] file size: 1647602665
[info] file size: 1647602665
[info][BGZF range] EOF, current offset is 1647602665
[info][chunkinputstream] len == 0, start offset = 1647602637, end offset = 1647602665
[info][chunkinputstream] len == 0, start offset = 1647602665, end offset = 1647602665
[info][bamRead range] empty, read 0 bytes, expected 4

real 1m26.204s
user 1m57.226s
sys 0m1.818s
6===
[info] file size: 1647602665
[info] file size: 1647602665
[info][BGZF range] EOF, current offset is 1647602665
[info][chunkinputstream] len == 0, start offset = 1647602637, end offset = 1647602665
[info][chunkinputstream] len == 0, start offset = 1647602665, end offset = 1647602665
[info][bamRead range] empty, read 0 bytes, expected 4

real 1m21.599s
user 1m56.810s
sys 0m1.728s

On the SSD partition:
avilella_scratch]$ for t in seq 8 -2 1; do echo $t===; time /home/avilella/src/sambamba/latest/sambamba/build/sambamba view -t $t TruSeqNano-350.chr19.bwamem.sorted.bam 1>/dev/null; done
8===
[info] file size: 1395101757
[info] file size: 1395101757
[info][BGZF range] EOF, current offset is 1395101757
[info][chunkinputstream] len == 0, start offset = 1395101729, end offset = 1395101757
[info][chunkinputstream] len == 0, start offset = 1395101757, end offset = 1395101757
[info][bamRead range] empty, read 0 bytes, expected 4

real 1m12.507s
user 1m43.656s
sys 0m2.190s
6===
[info] file size: 1395101757
[info] file size: 1395101757
[info][BGZF range] EOF, current offset is 1395101757
[info][chunkinputstream] len == 0, start offset = 1395101729, end offset = 1395101757
[info][chunkinputstream] len == 0, start offset = 1395101757, end offset = 1395101757
[info][bamRead range] empty, read 0 bytes, expected 4

real 1m8.823s
user 1m41.281s
sys 0m1.879s

Even after optimizations, it's about 1.5x slower than samtools.

The function is too long and cries out to be split into smaller functions.

It would be very useful to have an effective implementation of such thing.

For 'QNAME' sorting order it's trivial since paired reads come in adjacent lines, but the typical case is that of coordinate sorted BAM file, and another sorting takes a lot of time (and disk space as well), though paired reads are usually (relatively) near each other since reads are short. But not all the time. I've wrote a simple utility, and for a sample file it showed that average distance between paired reads is about 40, but the distribution has extreeeemely looooong taaaaail, it well might be that first read is at the beginning of the file while the second one is at the end. But the average suggests that additional sorting is not the best way to handle it.

Picard uses such kind of thing for validating mate existence and fields, and it stores temporary data on disk (picard/sam/SamFileValidator.java). Basically, when we find one read from a pair, store its read name; when we find its mate — remove that read from the map.

Amount of lines in Java implementation of disk-stored map is overwhelming... BUT seems to me those guys just like reinventing the wheel. Kyoto Cabinet comes to mind.

The help for 0.3.1 says this:

sambamba-flagstat [options] <input.bam>
-n, --nthreads=NTHREADS
-p, --show-progress

But I can't seem to get -p and -n working together:

These work:

sambamba flagsat input.bam 
sambamba flagsat -p input.bam 
sambamba flagsat -n=32 input.bam 

These don't:

sambamba flagsat -n 32 input.bam   # typical short option in Unix
sambamba flagsat -n=32 -p input.bam   
sambamba flagsat -n=32 --show-progress input.bam   

I won't go through all the combinations, you get my point.

The I/O performance data I gathered also indicates that using a large buffer for writes is helpful. Currently, sambamba view uses the default buffer size on stdout; I obtained a significant improvement by using setvbuf(3) like so:

import std.c.stdio : stdout, setvbuf, _IOFBF;
import std.c.stdlib : malloc;
[...]
char *buf = cast(char*)malloc(1048576);
setvbuf(stdout, buf, _IOFBF, 1048576);

This 1 MB buffer resulted in ~2 seconds being spent in write(2) calls, rather than ~4 with the default buffer size.

Similar results were obtained when using writeBAM with a BufferedFile; the default 8K buffer size resulted in 7-8 seconds spent in write(2) calls, vs. about 2.2 seconds when using a 128K buffer. (No significant improvement was observed with a 1M buffer.)

I'm not sure what the best way to standardize this across the codebase would be, given my limited D experience, so I'm not attempting to submit a pull request for this one at this point.

rdmd ./examples/snpcaller.d NA12043.chrom20.LS454.ssaha2.CEU.low_coverage.20101123.bam.v1

(hang)

A lot of new functionality was added recently, but all of it was tested only with DMD. Got segfault when tried to run unittests compiled with GDC. Needs investigation.

Sambamba tool lacks progressbar. Most BAM files are quite large, and it's desirable to see progressbar so as to estimate processing time. It should look like wget one.

Now the tool works with local files only, and I can add method, returning current file offset, to BamFile, and call it, say, every 1000 or 10000 alignments to update progressbar.

Some more issue tags would be useful

in progress
in testing
critical
important
newbie
small fix

The index tool should throw an error if the bin field is set incorrectly, because that would make the resulting index incorrect.

Turns out that switching to parallel decompression introduced a bug, files are truncated.

Seems that heuristic that it relies upon is wrong, or maybe there's a bug in algorithm. By design, number of detected SNPs shouldn't depend on chunk size.

However, there is a dependency:

  chunk size, bytes      |       #SNPs
  -------------------------------------------------
         10_000                    1309
         20_000                     635
         30_000                     964
         40_000                    1037
         50_000                     618
        100_000                     618
        150_000                     436
        200_000                     436

It should:

• be able to read from config file a list of errors to validate together with their priorities
• have an option to check only top N errors from a list sorted by priority
• have an option to stop after M errors with priority >= P
• have several verbosity modes (e.g. only summary, invalid read names + summary, invalid read names + error messages + summary, invalid reads in SAM format + error messages + summary)
• validate mate existence and paired read fields as Picard does (with the aid of KyotoCabinet on-disk hashtable)
• effectively use multithreading

As of now, there's a lot of templated code that doesn't restrict compile-time parameters in any way. Concepts should be introduced that specify the required behaviour, e.g. ReadRange, PileupColumnRange, PileupObject, etc.

Basically, the code in Sambamba tries to avoid overhead from D classes. However, object-oriented programming is still widely presented, as it makes it much easier to understand and write the code. As such, typically used programming constructs are

• structs - instead of classes
• alias this - instead of inheritance
• compile-time polymorphism - instead of run-time polymorphism
• (currently lacking) concepts - instead of interfaces

The purpose of concepts is two-fold -- making it easier to read the code of this library, and also introduce somewhat stricter typing.

External sort algorithms are to be studied and evaluated.

Currently, Ruby interpreter just crashes when D runtime throws an exception. That should be fixed.

Ideally, including tags. At least modifying alignments should be easy: change properties of primitive types in place, and for arrays use std.array.replaceSlice (with updating offsets before that).
As for tags, the same idea can be applied, and as I don't currently store any offsets, it's going to be even easier a bit.

The benefits of this approach are memory reuse and easy output to bam file (just do endian conversion if needed, and write the whole chunk).

Need to find a workaround, and submit a patch to Phobos if possible.

Currently, there are no convenient methods to edit header text, such as addComment(), addRgLine(), and so on. All this stuff should be encapsulated, and access should be through properties. This becomes important now that the library has BAM output capabilities.

The library should support validation of every alignment record, and SAM header. Each entity should have 'is_valid' property, and then functions like filterValid() and createValidationReport() should be available. Ruby code should have access to those high-level functions.

The library should support working with standard streams, though stdin clearly has limitations on random access.

When changing the chunk size you get the same SNPs, but probabilities may differ (probably on border SNPs). Example:

2474c2474

< 13695679 C T 696.606

13695679 C T 699.458
4656c4656

< 29623222 T T,C 248.637

29623222 T T,C 270.62
5438c5438

< 32990049 G T,G 490.482

32990049 G T,G 496.058
5576c5576

< 33509522 C T,C 492.262

33509522 C T,C 492.342

Add ability to clear all tags, and to append tag without checking if it already exists (for better performance).

Due to Amdahl's law, that's the biggest bottleneck at the moment. Every opportunity should be considered, including using ASM in some places.

Currently, there're two such places: reading the file contents, and reading alignments from the stream composed of decompressed BGZF blocks.

That would be useful for pipelining with other tools.