.travis.yml

 language: d
-script:
-- rdmd -unittest -d test/unittests.d
+
+os:
+  - linux
+  - osx
+
+d:
+  - dmd
+  - ldc

README.md

-### BioD
+# BioD [![Build Status](https://travis-ci.org/biod/BioD.svg?branch=master)](https://travis-ci.org/biod/BioD) [![DUB Package](https://img.shields.io/badge/dub-v0.1.0-red.svg)](https://code.dlang.org/packages/biod)
 
 [BioD](https://github.com/biod/BioD) is a fast and memory efficient bioinformatics library written in the [D programming language](http://dlang.org).
+
 BioD aims to:
 
 * Provide a platform for writing high-performance bioinformatics applications in D. BioD achieves this by:
-     - automatic parallelization of tasks where possible for example reading and writing BAM files.
+  - automatic parallelization of tasks where possible for example reading and writing BAM files
   - reducing the GC overhead by avoiding unnecessary memory allocations
-* Offer support for manipulating common biological data formats.
-* Write clear documented and maintainable codebase.
+* Offer support for manipulating common biological data formats
 
+## Why D?
 
-### Usage
-See the [examples directory](https://github.com/biod/BioD/tree/master/examples) for examples and usage.
+D is a language that suits parallel programming because of guarantees
+the compiler provides. D is both a low-level language and a high-level
+hybrid OOP/FP language. There is no other programming language that
+matches those features. Also, D templating/generics is far easier that
+that of C++ or, say, Scala.
 
+That is not to say that D is an easy language. A powerful toolbox will
+be complicated. If you want to do everything with a hammer, maybe
+better choose Java instead ;).
 
-#### BioD core
-* [Artem Tarasov](https://github.com/lomereiter)
-* [Pjotr Prins](https://github.com/pjotrp)
+For more information about D find Andrei Alexandrecu's D book. It is a
+classic. Ali Çehreli's book also is recommended.
+
+## Current development
+
+Our current focus is to provide a bamreader and bamwriter that is
+really fast and easy to use. We believe the BAM format is here to stay
+for the foreseeable future in pipelines. With D we have an good way to
+write performance parsers, particularly with three typical scenarios:
+
+1. Go through a BAM file a read at a time
+2. Go through a BAM file a nucleotide at a time (pileup)
+3. Go through a BAM file with a sliding window
+
+The sliding window is a derivation of the first - a read at a time or
+a nucleotide at a time.
+
+At this point this functionality is mostly in BioD, but not in an
+intuitive way. We are building up this functionality and will give
+examples (WIP).
+
+# Install
+
+The current default is to provide the path to the checked out repo to the D-compiler. For example
+in sambamba we use
 
+    DFLAGS = -wi -I. -IBioD -g
+
+# Usage
+
+See the [examples directory](https://github.com/biod/BioD/tree/master/examples)
+for examples and usage.
+
+BioD is also a crucial part of the [sambamba](https://github.com/biod/sambamba) tool.
+
+# Contributing
+
+Simply clone the repository on github and put in a pull request.
+
+# BioD contributors and support
+
+See
+[contributors](https://github.com/biod/BioD/graphs/contributors). For
+support use the [issue tracker](https://github.com/biod/BioD/issues) or contact
+
+* [Pjotr Prins](https://github.com/pjotrp)
+* [Artem Tarasov](https://github.com/lomereiter)
+* [George Githinji](https://github.com/George-Githinji)
 
-##### License
-This library is licensed under the MIT license. Please see the license file for more details
+# License
 
-##### Citation
+BioD is licensed under the liberal MIT (expat) [license](./LICENSE).

RELEASE-NOTES.md

+## ChangeLog v0.2.1 (20181004)
+
++ Fix bunch of deprecation warnings
++ Fixes sort removing tags from @RG header https://github.com/biod/sambamba/issues/356
++ Fixes https://github.com/biod/BioD/issues/37
+
+## ChangeLog v0.2.0 (20180915)
+
+64-bit compilation should be fine on ldc 1.10
+
++ Added FastQ parser (thanks George Githinji @george-githinji)
++ Added rewritten Bamreader f4a6c1c55c8903e948b793adbb150704d1e267b2 (thanks Pjotr Prins @pjotrp)
++ Improved big-endian support and other bug fixes (thanks Artem Tarasov @lomereiter)
++ Include SAM/BAM AH tag (thanks Indraniel Das @indraniel)
++ D undeaD compilation fixes (thanks John Colvin @John-Colvin)
++ Meson build definition added (thanks Matthias Klumpp @ximion)
++ Update outputstream.d (thanks Brett T. Hannigan @godotgildor)
++ Moved Cigar into its own module 7994406592f277bc6950739aaf75c5f948cb7928
++ Added 'asserte' method which throws an exception on assert (failure)
++ Bug fixes:
+  * #31 Bug in bio/core/utils/outbuffer.d: _heap.length is not correctly set affects reading long SAM records
+  * https://github.com/biod/sambamba/issues/210 correct pseudo-bin calculation (thanks @jfarek)
+  * JSON header output fixes #331

VERSION

+0.2.1

bio/bam/bai/indexing.d

 /*
     This file is part of BioD.
-    Copyright (C) 2012-2013    Artem Tarasov <lomereiter@gmail.com>
+    Copyright (C) 2012-2017    Artem Tarasov <lomereiter@gmail.com>
 
     Permission is hereby granted, free of charge, to any person obtaining a
     copy of this software and associated documentation files (the "Software"),
@@ -30,12 +30,12 @@ import bio.bam.constants;
 import bio.bam.bai.bin;
 import bio.core.bgzf.chunk;
 
-import std.stream;
+import undead.stream;
 import std.array;
 import std.algorithm;
 import std.system;
 import std.exception;
-import std.c.string;
+import core.stdc.string;
 
 // Suppose we have an alignment which covers bases on a reference,
 // starting from one position and ending at another position.
@@ -117,7 +117,8 @@ struct IndexBuilder {
         void updateMetadata(ref BamReadBlock read) {
             if (read.ref_id == -1) {
                 ++_no_coord;
-            } else if (read.is_unmapped) {
+            } else {
+                if (read.is_unmapped) {
                     ++_unmapped;
                 } else {
                     ++_mapped;
@@ -127,6 +128,7 @@ struct IndexBuilder {
                     _beg_vo = cast(ulong)read.start_virtual_offset;
                 _end_vo = cast(ulong)read.end_virtual_offset;
             }
+        }
 
         void updateLinearIndex() {
             assert(_prev_read.ref_id >= 0);

bio/bam/baifile.d

@@ -28,7 +28,7 @@ public import bio.bam.bai.bin;
 import bio.core.bgzf.virtualoffset;
 import bio.bam.constants;
 
-import std.stream;
+import undead.stream;
 import std.system;
 import std.exception;
 import std.algorithm;

bio/bam/baseinfo.d

@@ -30,6 +30,7 @@ import bio.bam.read;
 import bio.bam.tagvalue;
 import bio.bam.iontorrent.flowcall;
 import bio.bam.md.core;
+import bio.bam.cigar;
 
 import std.range;
 import std.conv;

bio/bam/cigar.d

+/*
+    This file is part of BioD.
+    Copyright (C) 2012-2016    Artem Tarasov <lomereiter@gmail.com>
+
+    Permission is hereby granted, free of charge, to any person obtaining a
+    copy of this software and associated documentation files (the "Software"),
+    to deal in the Software without restriction, including without limitation
+    the rights to use, copy, modify, merge, publish, distribute, sublicense,
+    and/or sell copies of the Software, and to permit persons to whom the
+    Software is furnished to do so, subject to the following conditions:
+
+    The above copyright notice and this permission notice shall be included in
+    all copies or substantial portions of the Software.
+
+    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+    IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+    FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+    AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+    LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+    DEALINGS IN THE SOFTWARE.
+
+*/
+
+module bio.bam.cigar;
+
+import std.algorithm;
+import std.range;
+import std.conv;
+import std.format;
+import std.exception;
+import std.system;
+import std.traits;
+import std.array;
+import std.bitmanip;
+import core.stdc.stdlib;
+
+import bio.core.base;
+import bio.core.utils.format;
+
+import bio.bam.abstractreader;
+import bio.bam.cigar;
+import bio.bam.writer;
+import bio.bam.tagvalue;
+import bio.bam.bai.bin;
+
+import bio.bam.md.core;
+
+import bio.bam.utils.array;
+import bio.bam.utils.value;
+import bio.core.utils.switchendianness;
+
+import bio.bam.thirdparty.msgpack : Packer, unpack;
+
+/**
+  Represents single CIGAR operation
+ */
+struct CigarOperation {
+    static assert(CigarOperation.sizeof == uint.sizeof);
+    /*
+        WARNING!
+
+      It is very essential that the size of
+      this struct is EXACTLY equal to uint.sizeof!
+
+      The reason is to avoid copying of arrays during alignment parsing.
+
+      Namely, when some_pointer points to raw cigar data,
+      we can just do a cast. This allows to access those data
+      directly, not doing any memory allocations.
+    */
+
+    private uint raw; // raw data from BAM
+
+    private static ubyte char2op(char c) {
+        switch(c) {
+            case 'M': return 0;
+            case 'I': return 1;
+            case 'D': return 2;
+            case 'N': return 3;
+            case 'S': return 4;
+            case 'H': return 5;
+            case 'P': return 6;
+            case '=': return 7;
+            case 'X': return 8;
+            default:  return 15; // 15 is used as invalid value
+        }
+    }
+
+    /// Length must be strictly less than 2^28.
+    /// $(BR)
+    /// Operation type must be one of M, I, D, N, S, H, P, =, X.
+    this(uint length, char operation_type) {
+        enforce(length < (1<<28), "Too big length of CIGAR operation");
+        raw = (length << 4) | char2op(operation_type);
+    }
+
+    this(uint _raw) {
+        raw = _raw;
+    }
+
+    /// Operation length
+    uint length() @property const nothrow @nogc {
+        return raw >> 4;
+    }
+
+    /// CIGAR operation as one of MIDNSHP=X.
+    /// Absent or invalid operation is represented by '?'
+    char type() @property const nothrow @nogc {
+        return "MIDNSHP=X????????"[raw & 0xF];
+    }
+
+    // Each pair of bits has first bit set iff the operation is query consuming,
+    // and second bit set iff it is reference consuming.
+    //                                            X  =  P  H  S  N  D  I  M
+    private static immutable uint CIGAR_TYPE = 0b11_11_00_00_01_10_10_01_11;
+
+    /// True iff operation is one of M, =, X, I, S
+    bool is_query_consuming() @property const nothrow @nogc {
+        return ((CIGAR_TYPE >> ((raw & 0xF) * 2)) & 1) != 0;
+    }
+
+    /// True iff operation is one of M, =, X, D, N
+    bool is_reference_consuming() @property const nothrow @nogc {
+        return ((CIGAR_TYPE >> ((raw & 0xF) * 2)) & 2) != 0;
+    }
+
+    /// True iff operation is one of M, =, X
+    bool is_match_or_mismatch() @property const nothrow @nogc {
+        return ((CIGAR_TYPE >> ((raw & 0xF) * 2)) & 3) == 3;
+    }
+
+    /// True iff operation is one of 'S', 'H'
+    bool is_clipping() @property const nothrow @nogc {
+        return ((raw & 0xF) >> 1) == 2; // 4 or 5
+    }
+
+    void toSam(Sink)(auto ref Sink sink) const
+        if (isSomeSink!Sink)
+    {
+        sink.write(length);
+        sink.write(type);
+    }
+
+    void toString(scope void delegate(const(char)[]) sink) const {
+        toSam(sink);
+    }
+}
+
+alias CigarOperation[] CigarOperations;
+
+bool is_unavailable(CigarOperations cigar) @property nothrow @nogc {
+  return (cigar.length == 1 && cigar[0].raw == '*');
+}
+
+/// Forward range of extended CIGAR operations, with =/X instead of M
+/// Useful for, e.g., detecting positions of mismatches.
+struct ExtendedCigarRange(CigarOpRange, MdOpRange) {
+    static assert(isInputRange!CigarOpRange && is(Unqual!(ElementType!CigarOpRange) == CigarOperation));
+    static assert(isInputRange!MdOpRange && is(Unqual!(ElementType!MdOpRange) == MdOperation));
+
+    private {
+        CigarOpRange _cigar;
+        MdOpRange _md_ops;
+        CigarOperation _front_cigar_op;
+        MdOperation _front_md_op;
+        uint _n_mismatches;
+        bool _empty;
+    }
+
+    ///
+    this(CigarOpRange cigar, MdOpRange md_ops) {
+        _cigar = cigar;
+        _md_ops = md_ops;
+        fetchNextCigarOp();
+        fetchNextMdOp();
+    }
+
+    /// Forward range primitives
+    bool empty() @property const {
+        return _empty;
+    }
+
+    /// ditto
+    CigarOperation front() @property {
+        debug {
+            import std.stdio;
+            writeln(_front_cigar_op, " - ", _front_md_op);
+        }
+
+        if (_front_cigar_op.type != 'M')
+            return _front_cigar_op;
+
+        if (_n_mismatches == 0) {
+            assert(_front_md_op.is_match);
+            uint len = min(_front_md_op.match, _front_cigar_op.length);
+            return CigarOperation(len, '=');
+        }
+
+        assert(_front_md_op.is_mismatch);
+        return CigarOperation(min(_n_mismatches, _front_cigar_op.length), 'X');
+    }
+
+    /// ditto
+    ExtendedCigarRange save() @property {
+        typeof(return) r = void;
+        r._cigar = _cigar.save;
+        r._md_ops = _md_ops.save;
+        r._front_cigar_op = _front_cigar_op;
+        r._front_md_op = _front_md_op;
+        r._n_mismatches = _n_mismatches;
+        r._empty = _empty;
+        return r;
+    }
+
+    /// ditto
+    void popFront() {
+        if (!_front_cigar_op.is_match_or_mismatch) {
+            if (_front_cigar_op.is_reference_consuming)
+                fetchNextMdOp();
+            fetchNextCigarOp();
+            return;
+        }
+
+        auto len = _front_cigar_op.length;
+        if (_n_mismatches > 0) {
+            enforce(_front_md_op.is_mismatch);
+
+            if (len > _n_mismatches) {
+                _front_cigar_op = CigarOperation(len - _n_mismatches, 'M');
+                _n_mismatches = 0;
+                fetchNextMdOp();
+            } else if (len < _n_mismatches) {
+                _n_mismatches -= len;
+                fetchNextCigarOp();
+            } else {
+                fetchNextCigarOp();
+                fetchNextMdOp();
+            }
+        } else {
+            enforce(_front_md_op.is_match);
+            auto n_matches = _front_md_op.match;
+
+            if (len > n_matches) {
+                _front_cigar_op = CigarOperation(len - n_matches, 'M');
+                fetchNextMdOp();
+            } else if (len < n_matches) {
+                _front_md_op.match -= len;
+                fetchNextCigarOp();
+            } else {
+                fetchNextCigarOp();
+                fetchNextMdOp();
+            }
+        }
+    }
+
+    private {
+        void fetchNextCigarOp() {
+            if (_cigar.empty) {
+                _empty = true;
+                return;
+            }
+
+            _front_cigar_op = _cigar.front;
+            _cigar.popFront();
+        }
+
+        void fetchNextMdOp() {
+            if (_md_ops.empty)
+                return;
+
+            _n_mismatches = 0;
+
+            _front_md_op = _md_ops.front;
+            _md_ops.popFront();
+
+            if (_front_md_op.is_mismatch) {
+                _n_mismatches = 1;
+                while (!_md_ops.empty && _md_ops.front.is_mismatch) {
+                    _md_ops.popFront();
+                    _n_mismatches += 1;
+                }
+            }
+        }
+    }
+}
+
+auto makeExtendedCigar(CigarOpRange, MdOpRange)(CigarOpRange cigar, MdOpRange md_ops) {
+    return ExtendedCigarRange!(CigarOpRange, MdOpRange)(cigar, md_ops);
+}

bio/bam/md/reconstruct.d

@@ -23,6 +23,7 @@
 */
 module bio.bam.md.reconstruct;
 
+import bio.bam.cigar;
 import bio.bam.read;
 import bio.bam.md.core;
 
@@ -215,7 +216,7 @@ auto dna(T)(T read)
 unittest {
 
     import std.stdio;
-    writeln("Testing reconstruction of reference from MD tags and CIGAR");
+    // stderr.writeln("Testing reconstruction of reference from MD tags and CIGAR");
 
     // Test reference reconstruction from MD and CIGAR.
     // (Tests are taken from http://davetang.org/muse/2011/01/28/perl-and-sam/)

bio/bam/multireader.d

@@ -66,7 +66,7 @@ bool compare(T)(auto ref T r1, auto ref T r2) {
         sorting_order = r1[1].merged_header.sorting_order;
 
     if (sorting_order == SortingOrder.coordinate)
-        return compareCoordinates(r1[0], r2[0]);
+        return compareCoordinatesAndStrand(r1[0], r2[0]);
     else if (sorting_order == SortingOrder.queryname)
         return compareReadNames(r1[0], r2[0]);
     else

bio/bam/pileup.d

@@ -69,6 +69,7 @@
 /// ---------------------------------------------------------
 module bio.bam.pileup;
 
+import bio.bam.cigar;
 import bio.bam.read;
 import bio.bam.md.reconstruct;
 import bio.bam.splitter;
@@ -760,7 +761,7 @@ unittest {
     auto reference = to!string(dna(reads));
 
     import std.stdio;
-    writeln("Testing pileup (low-level aspects)...");
+    // stderr.writeln("Testing pileup (low-level aspects)...");
 
     auto pileup = makePileup(reads, true, 796, 849, false);
     auto pileup2 = makePileup(reads, true, 0, ulong.max, false);
@@ -864,6 +865,11 @@ struct PileupChunkRange(C) {
     private bool _use_md_tag;
     private ulong _start_from;
     private ulong _end_at;
+    private int _chunk_right_end;
+
+    private int computeRightEnd(ref ElementType!C chunk) {
+        return chunk.map!(r => r.position + r.basesCovered()).reduce!max;
+    }
 
     this(C chunks, bool use_md_tag, ulong start_from, ulong end_at) {
         _chunks = chunks;
@@ -885,13 +891,12 @@ struct PileupChunkRange(C) {
                     break;
                 }
 
-                auto last_read = _current_chunk[$-1];
-                if (last_read.position + last_read.basesCovered() > start_from) {
+                _chunk_right_end = computeRightEnd(_current_chunk);
+                if (_chunk_right_end > start_from)
                     break;
             }
         }
     }
-    }
 
     bool empty() @property {
         return _empty;
@@ -900,7 +905,7 @@ struct PileupChunkRange(C) {
     auto front() @property {
         auto end_pos = _current_chunk[$-1].position;
         if (_chunks.empty || _chunks.front[0].ref_id != _current_chunk[$-1].ref_id)
-            end_pos += _current_chunk[$-1].basesCovered();
+            end_pos = _chunk_right_end;
 
         return makePileup(chain(_prev_chunk, _current_chunk),
                           _use_md_tag,
@@ -921,6 +926,8 @@ struct PileupChunkRange(C) {
             if (_current_chunk[0].ref_id >= 0) break;
         }
 
+        _chunk_right_end = computeRightEnd(_current_chunk);
+
         // if we changed reference, nullify prev_chunk
         if (_prev_chunk.length > 0 &&
             _prev_chunk[$ - 1].ref_id == _current_chunk[0].ref_id)

bio/bam/randomaccessmanager.d

@@ -50,6 +50,7 @@ import std.traits;
 import std.exception;
 import std.container;
 import std.parallelism;
+static import std.file;
 
 debug {
     import std.stdio;

bio/bam/read.d

@@ -51,6 +51,7 @@ import bio.core.base;
 import bio.core.utils.format;
 
 import bio.bam.abstractreader;
+import bio.bam.cigar;
 import bio.bam.writer;
 import bio.bam.tagvalue;
 import bio.bam.bai.bin;
@@ -71,234 +72,8 @@ import std.exception;
 import std.system;
 import std.traits;
 import std.array;
-import std.c.stdlib;
-
-/**
-  Represents single CIGAR operation
- */
-struct CigarOperation {
-    static assert(CigarOperation.sizeof == uint.sizeof);
-    /*
-        WARNING!
-
-      It is very essential that the size of 
-      this struct is EXACTLY equal to uint.sizeof!
-
-      The reason is to avoid copying of arrays during alignment parsing.
-
-      Namely, when some_pointer points to raw cigar data,
-      we can just do a cast. This allows to access those data
-      directly, not doing any memory allocations. 
-    */
-
-    private uint raw; // raw data from BAM
-
-    private static ubyte char2op(char c) {
-        switch(c) {
-            case 'M': return 0;
-            case 'I': return 1;
-            case 'D': return 2;
-            case 'N': return 3;
-            case 'S': return 4;
-            case 'H': return 5;
-            case 'P': return 6;
-            case '=': return 7;
-            case 'X': return 8;
-            default:  return 15; // 15 is used as invalid value
-        }
-    }
-
-    /// Length must be strictly less than 2^28.
-    /// $(BR)
-    /// Operation type must be one of M, I, D, N, S, H, P, =, X.
-    this(uint length, char operation_type) {
-        enforce(length < (1<<28), "Too big length of CIGAR operation");
-        raw = (length << 4) | char2op(operation_type);
-    }
-
-    /// Operation length
-    uint length() @property const nothrow {
-        return raw >> 4;
-    }
-  
-    /// CIGAR operation as one of MIDNSHP=X.
-    /// Absent or invalid operation is represented by '?'
-    char type() @property const nothrow {
-        return "MIDNSHP=X????????"[raw & 0xF];
-    }
-
-    // Each pair of bits has first bit set iff the operation is query consuming,
-    // and second bit set iff it is reference consuming.
-    //                                            X  =  P  H  S  N  D  I  M
-    private static immutable uint CIGAR_TYPE = 0b11_11_00_00_01_10_10_01_11;
-
-    /// True iff operation is one of M, =, X, I, S
-    bool is_query_consuming() @property const {
-        return ((CIGAR_TYPE >> ((raw & 0xF) * 2)) & 1) != 0;
-    }
-
-    /// True iff operation is one of M, =, X, D, N
-    bool is_reference_consuming() @property const {
-        return ((CIGAR_TYPE >> ((raw & 0xF) * 2)) & 2) != 0;
-    }
-
-    /// True iff operation is one of M, =, X
-    bool is_match_or_mismatch() @property const {
-        return ((CIGAR_TYPE >> ((raw & 0xF) * 2)) & 3) == 3;
-    }
-
-    /// True iff operation is one of 'S', 'H'
-    bool is_clipping() @property const {
-        return ((raw & 0xF) >> 1) == 2; // 4 or 5
-    }
-
-    private void toSam(Sink)(auto ref Sink sink) const 
-        if (isSomeSink!Sink)
-    {
-        sink.write(length);
-        sink.write(type);
-    }
-
-    void toString(scope void delegate(const(char)[]) sink) const {
-        toSam(sink);
-    }
-}
-
-/// Forward range of extended CIGAR operations, with =/X instead of M
-/// Useful for, e.g., detecting positions of mismatches.
-struct ExtendedCigarRange(CigarOpRange, MdOpRange) {
-    static assert(isInputRange!CigarOpRange && is(Unqual!(ElementType!CigarOpRange) == CigarOperation));
-    static assert(isInputRange!MdOpRange && is(Unqual!(ElementType!MdOpRange) == MdOperation));
-  
-    private {
-        CigarOpRange _cigar;
-        MdOpRange _md_ops;
-        CigarOperation _front_cigar_op;
-        MdOperation _front_md_op;
-        uint _n_mismatches;
-        bool _empty;
-    }
-
-    ///
-    this(CigarOpRange cigar, MdOpRange md_ops) {
-        _cigar = cigar;
-        _md_ops = md_ops;
-        fetchNextCigarOp();
-        fetchNextMdOp();
-    }
-
-    /// Forward range primitives
-    bool empty() @property const {
-        return _empty;
-    }
-
-    /// ditto
-    CigarOperation front() @property {
-        debug {
-            import std.stdio;
-            writeln(_front_cigar_op, " - ", _front_md_op);
-        }
-        
-        if (_front_cigar_op.type != 'M')
-            return _front_cigar_op;
-
-        if (_n_mismatches == 0) {
-            assert(_front_md_op.is_match);
-            uint len = min(_front_md_op.match, _front_cigar_op.length);
-            return CigarOperation(len, '=');
-        }
-    
-        assert(_front_md_op.is_mismatch);
-        return CigarOperation(min(_n_mismatches, _front_cigar_op.length), 'X');
-    }
-
-    /// ditto
-    ExtendedCigarRange save() @property {
-        typeof(return) r = void;
-        r._cigar = _cigar.save;
-        r._md_ops = _md_ops.save;
-        r._front_cigar_op = _front_cigar_op;
-        r._front_md_op = _front_md_op;
-        r._n_mismatches = _n_mismatches;
-        r._empty = _empty;
-        return r;
-    }
-
-    /// ditto
-    void popFront() {
-        if (!_front_cigar_op.is_match_or_mismatch) {
-            if (_front_cigar_op.is_reference_consuming)
-                fetchNextMdOp();
-            fetchNextCigarOp();
-            return;
-        }
-
-        auto len = _front_cigar_op.length;
-        if (_n_mismatches > 0) {
-            enforce(_front_md_op.is_mismatch);
-      
-            if (len > _n_mismatches) {
-                _front_cigar_op = CigarOperation(len - _n_mismatches, 'M');
-                _n_mismatches = 0;
-                fetchNextMdOp();
-            } else if (len < _n_mismatches) {
-                _n_mismatches -= len;
-                fetchNextCigarOp();
-            } else {
-                fetchNextCigarOp();
-                fetchNextMdOp();
-            }
-        } else {
-            enforce(_front_md_op.is_match);
-            auto n_matches = _front_md_op.match;
-      
-            if (len > n_matches) {
-                _front_cigar_op = CigarOperation(len - n_matches, 'M');
-                fetchNextMdOp();
-            } else if (len < n_matches) {
-                _front_md_op.match -= len;
-                fetchNextCigarOp();
-            } else {
-                fetchNextCigarOp();
-                fetchNextMdOp();
-            }
-        }
-    }
-  
-    private {
-        void fetchNextCigarOp() {
-            if (_cigar.empty) {
-                _empty = true;
-                return;
-            }
-        
-            _front_cigar_op = _cigar.front;
-            _cigar.popFront();
-        }
-
-        void fetchNextMdOp() {
-            if (_md_ops.empty)
-                return;
-
-            _n_mismatches = 0;
-      
-            _front_md_op = _md_ops.front;
-            _md_ops.popFront();
-
-            if (_front_md_op.is_mismatch) {
-                _n_mismatches = 1;
-                while (!_md_ops.empty && _md_ops.front.is_mismatch) {
-                    _md_ops.popFront();
-                    _n_mismatches += 1;
-                }
-            }
-        }
-    }
-}
-
-auto makeExtendedCigar(CigarOpRange, MdOpRange)(CigarOpRange cigar, MdOpRange md_ops) {
-    return ExtendedCigarRange!(CigarOpRange, MdOpRange)(cigar, md_ops);
-}
+import std.bitmanip;
+import core.stdc.stdlib;
 
 /**
   BAM record representation.
@@ -704,7 +479,7 @@ struct BamRead {
     /*
       Constructs the struct from memory chunk
       */
-    this(ubyte[] chunk) {
+    this(ubyte[] chunk, bool fix_byte_order=true) {
 
         // Switching endianness lazily is not a good idea:
         //
@@ -720,7 +495,7 @@ struct BamRead {
         _chunk = chunk;
         this._is_slice = true;
 
-        if (std.system.endian != Endian.littleEndian) {
+        if (fix_byte_order && std.system.endian != Endian.littleEndian) {
             switchChunkEndianness();
 
             // Dealing with tags is the responsibility of TagStorage.
@@ -1232,7 +1007,7 @@ private:
         assert(n < 16);
         // http://graphics.stanford.edu/~seander/bithacks.html#ConditionalSetOrClearBitsWithoutBranching
         ushort mask = cast(ushort)(1 << n);
-        _flag = (_flag & ~mask) | ((-cast(int)b) & mask);
+        _flag = (_flag & ~cast(int)(mask)) | ((-cast(int)b) & mask);
     }
 
     // If _chunk is still a slice, not an array, duplicate it.
@@ -1507,7 +1282,7 @@ unittest {
     import std.stdio;
     import std.math;
 
-    writeln("Testing BamRead behaviour...");
+    // stderr.writeln("Testing BamRead behaviour...");
     auto read = BamRead("readname",
                         "AGCTGACTACGTAATAGCCCTA",
                         [CigarOperation(22, 'M')]);
@@ -1578,6 +1353,7 @@ unittest {
 
     {
     import std.typecons;
+    static import bio.bam.thirdparty.msgpack;
     auto packer = bio.bam.thirdparty.msgpack.packer(Appender!(ubyte[])());
     read.toMsgpack(packer);
     auto data = packer.stream.data;
@@ -1748,7 +1524,9 @@ unittest {
 ///   $(UL
 ///     $(LI two reads (in this case, reference IDs are also taken into account))
 ///     $(LI read and integer in any order)))
-bool compareCoordinates(R1, R2)(const auto ref R1 a1, const auto ref R2 a2)
+
+/// This function takes read direction into account (used for original samtools style sorting)
+bool compareCoordinatesAndStrand(R1, R2)(const auto ref R1 a1, const auto ref R2 a2)
     if (isBamRead!R1 && isBamRead!R2)
 {
     if (a1.ref_id == -1) return false; // unmapped reads should be last
@@ -1760,6 +1538,16 @@ bool compareCoordinates(R1, R2)(const auto ref R1 a1, const auto ref R2 a2)
     return !a1.is_reverse_strand && a2.is_reverse_strand;
 }
 
+bool compareCoordinates(R1, R2)(const auto ref R1 a1, const auto ref R2 a2)
+    if (isBamRead!R1 && isBamRead!R2)
+{
+    if (a1.ref_id == -1) return false; // unmapped reads should be last
+    if (a2.ref_id == -1) return true;
+    if (a1.ref_id < a2.ref_id) return true;
+    if (a1.ref_id > a2.ref_id) return false;
+    return (a1.position < a2.position);
+}
+
 bool compareCoordinates(R1, R2)(const auto ref R1 a1, const auto ref R2 a2)
     if (isBamRead!R1 && isIntegral!R2)
 {

bio/bam/reader.d

@@ -94,8 +94,8 @@ class BamReader : IBamSamReader {
       }
       -------------------------------------------
      */
-    this(std.stream.Stream stream, 
-         std.parallelism.TaskPool task_pool = std.parallelism.taskPool) {
+    this(undead.stream.Stream stream,
+         TaskPool task_pool = taskPool) {
         _source_stream = new EndianStream(stream, Endian.littleEndian);
         _task_pool = task_pool;
 

bio/bam/readrange.d

 /*
     This file is part of BioD.
-    Copyright (C) 2012-2014    Artem Tarasov <lomereiter@gmail.com>
+    Copyright (C) 2012-2016    Artem Tarasov <lomereiter@gmail.com>
 
     Permission is hereby granted, free of charge, to any person obtaining a
     copy of this software and associated documentation files (the "Software"),
@@ -28,11 +28,11 @@ import bio.bam.abstractreader;
 import bio.bam.reader;
 import bio.core.bgzf.inputstream;
 import bio.core.bgzf.virtualoffset;
-import bio.core.utils.switchendianness;
 
-import std.stream;
+import undead.stream;
 import std.algorithm;
 import std.system;
+import std.bitmanip;
 
 /// Read + its start/end virtual offsets
 struct BamReadBlock {
@@ -133,14 +133,11 @@ private:
         }
         beforeNextBamReadLoad();
 
-        // (FIXME: that won't work on Big Endian systems!)
-        
         // Here's where _empty is really set!
-        int block_size = void;
-        ubyte* ptr = cast(ubyte*)(&block_size);
+        ubyte[int.sizeof] tmp = void;
         auto _read = 0;
         while (_read < int.sizeof) {
-            auto _actually_read = _endian_stream.readBlock(ptr, int.sizeof - _read);
+            auto _actually_read = _endian_stream.readBlock(tmp.ptr + _read, int.sizeof - _read);
             if (_actually_read == 0) {
                 version(development) {
                     import std.stdio;
@@ -150,18 +147,15 @@ private:
                 return;
             }
             _read += _actually_read;
-            ptr += _actually_read;
         }
 
+        int block_size = littleEndianToNative!int(tmp);
+
         version(extraVerbose) {
             import std.stdio;
             stderr.writeln("[uncompressed] record size: ", block_size);
         }
 
-        if (std.system.endian != Endian.littleEndian) {
-            switchEndianness(&block_size, int.sizeof);
-        }
-
         ubyte[] data = void;
         if (_reader !is null && _reader._seqprocmode) {
             if (block_size > _buffer.length)

bio/bam/reference.d

@@ -47,7 +47,7 @@ import bio.bam.region;
 import bio.bam.randomaccessmanager;
 import bio.core.bgzf.virtualoffset;
 
-import std.stream;
+import undead.stream;
 import std.exception;
 import std.array;
 

bio/bam/referenceinfo.d

@@ -23,7 +23,7 @@
 */
 module bio.bam.referenceinfo;
 
-import std.stream;
+import undead.stream;
 import std.exception;
 import std.array;
 

bio/bam/snpcallers/maq.d

@@ -5,7 +5,7 @@ module bio.bam.snpcallers.maq;
  */
 
 import core.stdc.math;
-import std.math : LN2, LN10, isnan;
+import std.math : LN2, LN10, isNaN;
 import std.traits;
 import std.range;
 import std.algorithm;

bio/bam/tagvalue.d

@@ -296,7 +296,7 @@ string injectOpCast() {
     }
 
     cs ~= `(is(T == string)) {` ~
-          `  if (is_string) {`
+          `  if (is_string) {` ~
           `    return bam_typeid == 'Z' ? u.Z : u.H;`~
           `  } else if (is_integer || is_float || is_character) {`~
           `    `~injectSwitchPrimitive("string")~

bio/bam/thirdparty/msgpack.d

@@ -59,7 +59,7 @@ else
 // for Converting Endian using ntohs and ntohl;
 version(Windows)
 {
-    import std.c.windows.winsock;
+    import core.stdc.windows.winsock;
 }
 else
 {
@@ -82,7 +82,7 @@ static if (real.sizeof == double.sizeof) {
     import std.numeric;
 }
 
-version(unittest) import std.file, std.c.string;
+version(unittest) import std.file, core.stdc.string;
 
 
 @trusted:
@@ -344,7 +344,7 @@ struct PackerImpl(Stream) if (isOutputRange!(Stream, ubyte) && isOutputRange!(St
      * Params:
      *  withFieldName = serialize class / struct with field name
      */
-    this(bool withFieldName = false)
+    this(bool withFieldName)
     {
         withFieldName_ = withFieldName;
     }
@@ -3145,7 +3145,7 @@ struct Value
      *  type  = the type of value.
      */
     @safe
-    this(Type type = Type.nil)
+    this(Type type)
     {
         this.type = type;
     }
@@ -3617,7 +3617,7 @@ struct Value
                         ret ^= value.toHash();
                     }
                     return ret;
-                } catch assert(0);
+                } catch (Exception) assert(0);
         }
     }
 }
@@ -4656,7 +4656,7 @@ mixin template MessagePackable(Members...)
             if (withFieldName) {
                 packer.beginMap(this.tupleof.length);
                 foreach (i, member; this.tupleof) {
-                    pack(getFieldName!(typeof(this), i));
+                    packer.pack(getFieldName!(typeof(this), i));
                     packer.pack(member);
                 }
             } else {