pileup_example.d

import bio.bam.reader;
import bio.bam.pileup;

import bio.core.base;
import bio.core.tinymap;

import std.stdio, std.parallelism, std.algorithm, std.array, std.range;

void main(string[] args) {
    defaultPoolThreads = 4;

    auto bam = new BamReader(args[1]);
    auto pileup = bam.reads().makePileup();

    foreach (column; pileup) {
        auto occurrences = TinyMap!(Base5, uint, useDefaultValue)(0);

        auto groups = column.reads().array()
            .sort!((r1, r2) => r1.strand < r2.strand, SwapStrategy.stable)
            .group!((r1, r2) => r1.strand == r2.strand &&
                                r1.sequence.equal(r2.sequence));

        foreach (read; groups.map!(rn => rn[0]))
            if (read.current_base != '-')
                occurrences[Base5(read.current_base)] += 1;

        auto supported_variants = TinyMap!(Base5, bool, useDefaultValue)(false); 
        immutable threshold = 4;
        foreach (base, n; occurrences)
            if (n >= threshold)
                supported_variants[base] = true;

        // ...
        writeln(bam.reference_sequences[column.ref_id].name, ": ",
                column.position, "\t",
                "ACGT".filter!(nuc => nuc.Base5() in supported_variants));
    }
}