biohazard-1.0.0: bioinformatics support library

Safe HaskellNone



Pileup, similar to Samtools

Pileup turns a sorted sequence of reads into a sequence of "piles", one for each site where a genetic variant might be called. We will scan each read's CIGAR line and MD field in concert with the sequence and effective quality. Effective quality is the lowest available quality score of QUAL, MAPQ, and BQ. For aDNA calling, a base is represented as four probabilities, derived from a position dependent damage model.



data PrimChunks Source #

The primitive pieces for genotype calling: A position, a base represented as four likelihoods, an inserted sequence, and the length of a deleted sequence. The logic is that we look at a base followed by some indel, and all those indels are combined into a single insertion and a single deletion.


Seek !Int PrimBase

skip to position (at start or after N operation)

Indel [Nucleotides] [DamagedBase] PrimBase

observed deletion and insertion between two bases


nothing anymore

data PrimBase Source #



more chunks


data DamagedBase Source #

Represents our knowledge about a certain base, which consists of the base itself (A,C,G,T, encoded as 0..3; no Ns), the quality score (anything that isn't A,C,G,T becomes A with quality 0), and a substitution matrix representing post-mortem but pre-sequencing substitutions.

Unfortunately, none of this can be rolled into something more simple, because damage and sequencing error behave so differently.

Damage information is polymorphic. We might run with a simple version (a matrix) for calling, but we need more (a matrix and a mutable matrix, I think) for estimation.



reference base from MD field


newtype DmgToken Source #




decompose :: DmgToken -> BamRaw -> [PosPrimChunks] Source #

Decomposes a BAM record into chunks suitable for piling up. We pick apart the CIGAR and MD fields, and combine them with sequence and quality as appropriate. Clipped bases are removed/skipped as needed. We also apply a substitution matrix to each base, which must be supplied along with the read.

data CallStats Source #

Statistics about a genotype call. Probably only useful for fitlering (so not very useful), but we keep them because it's easy to track them.




newtype V_Nuc Source #


V_Nuc (Vector Nucleotide) 


Eq V_Nuc Source # 


(==) :: V_Nuc -> V_Nuc -> Bool #

(/=) :: V_Nuc -> V_Nuc -> Bool #

Ord V_Nuc Source # 


compare :: V_Nuc -> V_Nuc -> Ordering #

(<) :: V_Nuc -> V_Nuc -> Bool #

(<=) :: V_Nuc -> V_Nuc -> Bool #

(>) :: V_Nuc -> V_Nuc -> Bool #

(>=) :: V_Nuc -> V_Nuc -> Bool #

max :: V_Nuc -> V_Nuc -> V_Nuc #

min :: V_Nuc -> V_Nuc -> V_Nuc #

Show V_Nuc Source # 


showsPrec :: Int -> V_Nuc -> ShowS #

show :: V_Nuc -> String #

showList :: [V_Nuc] -> ShowS #

type BasePile = [DamagedBase] Source #

Map quality and a list of encountered bases, with damage information and reference base if known.

type IndelPile = [(Qual, ([Nucleotides], [DamagedBase]))] Source #

Map quality and a list of encountered indel variants. The deletion has the reference sequence, if known, an insertion has the inserted sequence with damage information.

data Pile' a b Source #

Running pileup results in a series of piles. A Pile has the basic statistics of a VarCall, but no likelihood values and a pristine list of variants instead of a proper call. We emit one pile with two BasePiles (one for each strand) and one IndelPile (the one immediately following) at a time.




(Show b, Show a) => Show (Pile' a b) Source # 


showsPrec :: Int -> Pile' a b -> ShowS #

show :: Pile' a b -> String #

showList :: [Pile' a b] -> ShowS #

type Pile = Pile' (BasePile, BasePile) (IndelPile, IndelPile) Source #

Raw pile. Bases and indels are piled separately on forward and backward strands.

pileup :: Enumeratee [PosPrimChunks] [Pile] IO b Source #

The pileup enumeratee takes BamRaws, decomposes them, interleaves the pieces appropriately, and generates Piles. The output will contain at most one BasePile and one IndelPile for each position, piles are sorted by position.

This top level driver receives BamRaws. Unaligned reads and duplicates are skipped (but not those merely failing quality checks). Processing stops when the first read with invalid br_rname is encountered or a t end of file.

newtype PileM m a Source #

The pileup logic keeps a current coordinate (just two integers) and two running queues: one of active PrimBases that contribute to current genotype calling and on of waiting PrimBases that will contribute at a later point.

Oppan continuation passing style! Not only is the CPS version of the state monad (we have five distinct pieces of state) somewhat faster, we also need CPS to interact with the mechanisms of Iteratee. It makes implementing yield, peek, and bump straight forward.





Monad (PileM m) Source # 


(>>=) :: PileM m a -> (a -> PileM m b) -> PileM m b #

(>>) :: PileM m a -> PileM m b -> PileM m b #

return :: a -> PileM m a #

fail :: String -> PileM m a #

Functor (PileM m) Source # 


fmap :: (a -> b) -> PileM m a -> PileM m b #

(<$) :: a -> PileM m b -> PileM m a #

Applicative (PileM m) Source # 


pure :: a -> PileM m a #

(<*>) :: PileM m (a -> b) -> PileM m a -> PileM m b #

liftA2 :: (a -> b -> c) -> PileM m a -> PileM m b -> PileM m c #

(*>) :: PileM m a -> PileM m b -> PileM m b #

(<*) :: PileM m a -> PileM m b -> PileM m a #

type PileF m r = Refseq -> Int -> ([PrimBase], [PrimBase]) -> (Heap, Heap) -> (Stream [Pile] -> Iteratee [Pile] m r) -> Stream [PosPrimChunks] -> Iteratee [PosPrimChunks] m (Iteratee [Pile] m r) Source #

The things we drag along in PileM. Notes: * The active queue is a simple stack. We add at the front when we encounter reads, which reverses them. When traversing it, we traverse reads backwards, but since we accumulate the BasePile, it gets reversed back. The new active queue, however, is no longer reversed (as it should be). So after the traversal, we reverse it again. (Yes, it is harder to understand than using a proper deque type, but it is cheaper. There may not be much point in the reversing, though.)

upd_pos :: (Int -> Int) -> PileM m () Source #

yieldPile :: CallStats -> BasePile -> BasePile -> CallStats -> IndelPile -> IndelPile -> PileM m () Source #

Sends one piece of output downstream. You are not expected to understand how this works, but inlining eneeCheckIfDone plugged an annoying memory leak.

pileup' :: PileM m () Source #

The actual pileup algorithm. If active contains something, continue here. Else find the coordinate to continue from, which is the minimum of the next waiting coordinate and the next coordinate in input; if found, continue there, else we're all done.

p'feed_input :: PileM m () Source #

Feeds input as long as it starts at the current position

p'check_waiting :: PileM m () Source #

Checks waiting queue. If there is anything waiting for the current position, moves it to active queue.

p'scan_active :: PileM m ((CallStats, BasePile), (CallStats, BasePile), (CallStats, IndelPile), (CallStats, IndelPile)) Source #

Separately scans the two active queues and makes one BasePile from each. Also sees what's next in the PrimChunks: Indels contribute to two separate IndelPiles, Seeks are pushed back to the waiting queue, EndOfReads are removed, and everything else is added to two fresh active queues.

data Heap Source #

We need a simple priority queue. Here's a skew heap (specialized to strict Int priorities and PrimBase values).


Node !Int PrimBase Heap Heap