diff --git a/biohazard.cabal b/biohazard.cabal
--- a/biohazard.cabal
+++ b/biohazard.cabal
@@ -1,31 +1,26 @@
 Name:                biohazard
-Version:             0.6.16
+Version:             1.0.0
 Synopsis:            bioinformatics support library
 Description:         This is a collection of modules I separated from
                      various bioinformatics tools.
 Category:            Bioinformatics
 
-Homepage:            http://github.com/udo-stenzel/biohazard
+Homepage:            https://bitbucket.org/ustenzel/biohazard
 License:             BSD3
 License-File:        LICENSE
 
 Author:              Udo Stenzel
-Maintainer:          udo.stenzel@eva.mpg.de
+Maintainer:          u.stenzel@web.de
 Copyright:           (C) 2010-2017 Udo Stenzel
 
 Cabal-version:       >= 1.10
 Build-type:          Simple
-Tested-With:         GHC == 7.8.4, GHC == 7.10.1, GHC == 8.0.1
+Tested-With:         GHC == 7.8.4, GHC == 7.10.3, GHC == 8.0.2, GHC == 8.2.1
 
 source-repository head
   type:     git
   location: https://bitbucket.org/ustenzel/biohazard.git
 
-Flag debug
-  Description: enable additional sanity checks
-  Default:     False
-  Manual:      True
-
 Library
   Exposed-modules:     Bio.Adna,
                        Bio.Align,
@@ -61,11 +56,11 @@
 
   Build-depends:       async                    >= 2.0 && < 2.2,
                        attoparsec               >= 0.10 && < 0.14,
-                       base                     >= 4.6 && < 4.11,
-                       base-prelude             == 1.0.* || == 1.2.*,
+                       base                     >= 4.7 && < 4.11,
+                       base-prelude             == 1.2.*,
                        binary                   >= 0.7 && < 0.9,
                        bytestring               >= 0.10.2 && < 0.11,
-                       containers               >= 0.4.1 && < 0.6,
+                       containers               == 0.5.*,
                        directory                >= 1.2 && < 1.4,
                        exceptions               >= 0.6 && < 0.9,
                        filepath                 >= 1.3 && < 1.5,
@@ -91,6 +86,7 @@
                        DeriveDataTypeable,
                        FlexibleContexts,
                        FlexibleInstances,
+                       LambdaCase,
                        MultiParamTypeClasses,
                        NoImplicitPrelude,
                        OverloadedStrings,
@@ -113,12 +109,12 @@
                        UndecidableInstances
 
   Hs-source-dirs:      src
-  Install-Includes:    src/cbits/myers_align.h
+  Include-dirs:        src/cbits
+  Install-Includes:    myers_align.h
   C-sources:           src/cbits/loops.c,
                        src/cbits/mmap.c,
                        src/cbits/myers_align.c,
                        src/cbits/trim.c
   CC-options:          -fPIC
-
 
 -- :vim:tw=132:
diff --git a/src/Bio/Adna.hs b/src/Bio/Adna.hs
--- a/src/Bio/Adna.hs
+++ b/src/Bio/Adna.hs
@@ -1,4 +1,16 @@
-{-# LANGUAGE DeriveGeneric, CPP #-}
+{-# LANGUAGE DeriveGeneric #-}
+
+-- | Things specific to ancient DNA, e.g. damage models.
+--
+-- For aDNA, we need a substitution probability.  We have three options:
+-- use an empirically determined PSSM, use an arithmetically defined
+-- PSSM based on the /Johnson/ model, use a context sensitive PSSM based
+-- on the /Johnson/ model and an alignment.  Using /Dindel/, actual
+-- substitutions relative to a called haplotype would be taken into
+-- account.  Since we're not going to do that, taking alignments into
+-- account is difficult, somewhat approximate, and therefore not worth
+-- the hassle.
+
 module Bio.Adna (
     DmgStats(..),
     CompositionStats,
@@ -17,7 +29,11 @@
     Alignment(..),
     FragType(..),
     Subst(..),
+
     NPair,
+    npair,
+    fst_np,
+    snd_np,
 
     noDamage,
     univDamage,
@@ -41,20 +57,9 @@
 import qualified Data.Vector.Unboxed            as U
 import qualified Data.Vector.Unboxed.Mutable    as UM
 
--- ^ Things specific to ancient DNA, e.g. damage models.
---
--- For aDNA, we need a substitution probability.  We have three options:
--- use an empirically determined PSSM, use an arithmetically defined
--- PSSM based on the /Johnson/ model, use a context sensitive PSSM based
--- on the /Johnson/ model and an alignment.  Using /Dindel/, actual
--- substitutions relative to a called haplotype would be taken into
--- account.  Since we're not going to do that, taking alignments into
--- account is difficult, somewhat approximate, and therefore not worth
--- the hassle.
---
--- We represent substitution matrices by the type 'Mat44D'.  Internally,
+-- | We represent substitution matrices by the type 'Mat44D'.  Internally,
 -- this is a vector of packed vectors.  Conveniently, each of the packed
--- vectors represents all transition /into/ the given nucleotide.
+-- vectors represents all transitions /into/ the given nucleotide.
 
 newtype Mat44D = Mat44D (U.Vector Double) deriving (Show, Generic)
 newtype MMat44D = MMat44D (UM.IOVector Double)
@@ -250,8 +255,6 @@
 --   mind, 'substs5d5', 'substs5d3', 'substs5dd' are like 'substs5', but
 --   counting only reads where the 5' end is damaged, where the 3' end
 --   is damaged, and where both ends are damaged, respectively.
---
--- XXX  This got kind of ugly.  We'll see where this goes...
 
 data DmgStats a = DmgStats {
     basecompo5 :: CompositionStats,
@@ -274,14 +277,34 @@
 
 
 data FragType = Complete | Leading | Trailing deriving (Show, Eq)
-type NPair = ( Nucleotides, Nucleotides )
 
--- Alignment record, might have been gotten from practically anywhere
--- with varying completeness.  We record anything we can get, most is
--- optional.  Reference sequence is filled with Ns if missing.
+-- | Compact storage of a pair of ambiguous 'Nucleotides'.  Used to
+-- represent alignments in a way that is accessible even to assembly
+-- code.  The first and sencond field are stored in the low and high
+-- nybble, respectively.  See 'fst_np', 'snd_np', 'npair'.
+newtype NPair = NPair Word8 deriving (Eq, Ord)
+
+npair :: Nucleotides -> Nucleotides -> NPair
+npair (Ns r) (Ns q) = NPair $ shiftL q 4 .|. r .&. 0xF
+
+fst_np, snd_np :: NPair -> Nucleotides
+fst_np (NPair w) = Ns (w .&. 0xF)
+snd_np (NPair w) = Ns (shiftR w 4)
+
+instance Storable NPair where
+    sizeOf    _ = 1
+    alignment _ = 1
+    peek p = NPair <$> peek (castPtr p :: Ptr Word8)
+    poke p (NPair v) = poke (castPtr p :: Ptr Word8) v
+
+instance Show NPair where
+    showsPrec _ p = shows (fst_np p) . (:) '/' . shows (snd_np p)
+
+-- | Alignment record.  The reference sequence is filled with Ns if
+-- missing.
 data Alignment = ALN
-    { a_sequence :: !(U.Vector NPair)       -- the alignment proper
-    , a_fragment_type :: !FragType }    -- was the adapter trimmed?
+    { a_sequence :: !(VS.Vector NPair)      -- the alignment proper
+    , a_fragment_type :: !FragType }        -- was the adapter trimmed?
 
 addFragType :: BamMeta -> Enumeratee [BamRaw] [(BamRaw,FragType)] m b
 addFragType meta = mapStream $ \br -> (br, case unpackBam br of
@@ -353,15 +376,15 @@
         guard (not $ isUnmapped b)
         md <- getMd b
         let pps = alnFromMd b_seq b_cigar md
-            ref = U.map fromN $ U.filter ((/=) gap . fst) pps
+            ref = U.convert $ VS.map fromN $ VS.filter ((/=) gap) $ VS.map fst_np pps
         return (b, ft, ref, pps)) =$
     damagePatternsIter 0 rng it
   where
-    fromN (ns,_) | ns == nucsA = 2
-                 | ns == nucsC = 1
-                 | ns == nucsG = 3
-                 | ns == nucsT = 0
-                 | otherwise   = 4
+    fromN ns | ns == nucsA = 2
+             | ns == nucsC = 1
+             | ns == nucsG = 3
+             | ns == nucsT = 0
+             | otherwise   = 4
 
 -- | Common logic for statistics. The function 'get_ref_and_aln'
 -- reconstructs reference sequence and alignment from a Bam record.  It
@@ -370,7 +393,7 @@
 damagePatternsIter :: MonadIO m
                    => Int -> Int
                    -> Iteratee [Alignment] m b
-                   -> Iteratee [(BamRec, FragType, U.Vector Word8, U.Vector NPair)] m (DmgStats b)
+                   -> Iteratee [(BamRec, FragType, U.Vector Word8, VS.Vector NPair)] m (DmgStats b)
 damagePatternsIter ctx rng it = mapStream revcom_both =$ do
     let maxwidth = ctx + rng
     acc_bc <- liftIO $ UM.replicate (2 * 5 *    maxwidth) (0::Int)
@@ -378,13 +401,6 @@
     acc_cg <- liftIO $ UM.replicate (2 * 2 * 4 *     rng) (0::Int)
 
     it' <- flip mapStreamM it $ \(BamRec{..}, a_fragment_type, ref, a_sequence) -> liftIO $ do
-#ifdef DEBUG
-              when (U.any (<0) ref || U.any (>4) ref) . error $
-                    "Unexpected value in reference fragment: " ++ show ref
-#endif
-              let good_pairs     = U.indexed             a_sequence
-                  good_pairs_rev = U.indexed $ U.reverse a_sequence
-
               -- basecompositon near 5' end, near 3' end
               let (width5, width3) = case a_fragment_type of
                                             Leading -> (full_width, 0)
@@ -397,51 +413,52 @@
 
               -- For substitutions, decide what damage class we're in:
               -- 0 - no damage, 1 - damaged 5' end, 2 - damaged 3' end, 3 - both
-              let dmgbase = 2*4*4*rng * ( (if U.null a_sequence || U.head a_sequence /= (nucsC,nucsT) then 1 else 0)
-                                        + (if U.null a_sequence || U.last a_sequence /= (nucsC,nucsT) then 2 else 0) )
+              let dmgbase = 2*4*4*rng * ( (if VS.null a_sequence || VS.head a_sequence /= npair nucsC nucsT then 1 else 0)
+                                        + (if VS.null a_sequence || VS.last a_sequence /= npair nucsC nucsT then 2 else 0) )
 
               -- substitutions near 5' end
               let len_at_5 = case a_fragment_type of Leading  -> min rng (G.length b_seq)
                                                      Complete -> min rng (G.length b_seq `div` 2)
                                                      Trailing -> 0
-              U.forM_ (U.take len_at_5 good_pairs) $
-                    \(i,uv) -> withPair uv $ \j -> bump (j * rng + i + dmgbase) acc_st
+              flip G.imapM_ (VS.take len_at_5 a_sequence) $
+                    \i uv -> withPair uv $ \j -> bump (j * rng + i + dmgbase) acc_st
 
               -- substitutions at CpG sites near 5' end
-              U.zipWithM_
-                  (\(i,(u,v)) (_,(w,z)) ->
-                      when (u == nucsC && w == nucsG) $ do
-                          withNs v $ \y -> bump (  y   * rng +  i ) acc_cg
-                          withNs z $ \y -> bump ((y+4) * rng + i+1) acc_cg)
-                  (U.take len_at_5 good_pairs) (U.drop 1 good_pairs)
+              G.izipWithM_
+                  (\i uv wz ->
+                      when (fst_np uv == nucsC && fst_np wz == nucsG) $ do
+                          withNs (snd_np uv) $ \y -> bump (  y   * rng +  i ) acc_cg
+                          withNs (snd_np wz) $ \y -> bump ((y+4) * rng + i+1) acc_cg)
+                  (VS.take len_at_5 a_sequence) (VS.drop 1 a_sequence)
 
               -- substitutions near 3' end
               let len_at_3 = case a_fragment_type of Leading  -> 0
                                                      Complete -> min rng (G.length b_seq `div` 2)
                                                      Trailing -> min rng (G.length b_seq)
-              U.forM_ (U.take len_at_3 good_pairs_rev) $
-                    \(i,uv) -> withPair uv $ \j -> bump ((17+j) * rng -i -1 + dmgbase) acc_st
+              flip G.imapM_ (VS.take len_at_3 (VS.reverse a_sequence)) $
+                    \i uv -> withPair uv $ \j -> bump ((17+j) * rng -i -1 + dmgbase) acc_st
 
               -- substitutions at CpG sites near 3' end
-              U.zipWithM_
-                  (\(_,(u,v)) (i,(w,z)) ->
-                      when (u == nucsC && w == nucsG) $ do
-                          withNs v $ \y -> bump ((y+ 9) * rng - i-2) acc_cg
-                          withNs z $ \y -> bump ((y+13) * rng - i-1) acc_cg)
-                  (U.drop 1 good_pairs_rev) (U.take len_at_3 good_pairs_rev)
-
+              G.izipWithM_
+                  (\i wz uv ->
+                      when (fst_np uv == nucsC && fst_np wz == nucsG) $ do
+                          withNs (snd_np uv) $ \y -> bump ((y+ 9) * rng - i-2) acc_cg
+                          withNs (snd_np wz) $ \y -> bump ((y+13) * rng - i-1) acc_cg)
+                  (VS.take len_at_3 (VS.reverse a_sequence))
+                  (VS.drop 1 (VS.reverse a_sequence))
 
               return ALN{..}
 
+
     let nsubsts = 4*4*rng
         mk_substs off = sequence [ (,) (n1 :-> n2) <$> U.unsafeFreeze (UM.slice ((4*i+j)*rng + off*nsubsts) rng acc_st)
                                  | (i,n1) <- zip [0..] [nucA..nucT]
                                  , (j,n2) <- zip [0..] [nucA..nucT] ]
 
     accs <- liftIO $ DmgStats <$> sequence [ (,) nuc <$> U.unsafeFreeze (UM.slice (i*maxwidth) maxwidth acc_bc)
-                                           | (i,nuc) <- zip [2,1,3,0,4] [Just nucA, Just nucC, Just nucG, Just nucT, Nothing] ]
+                                           | (i,nuc) <- zip [2,1,3,0,4] [Just nucA,Just nucC,Just nucG,Just nucT,Nothing] ]
                               <*> sequence [ (,) nuc <$> U.unsafeFreeze (UM.slice (i*maxwidth) maxwidth acc_bc)
-                                           | (i,nuc) <- zip [7,6,8,5,9] [Just nucA, Just nucC, Just nucG, Just nucT, Nothing] ]
+                                           | (i,nuc) <- zip [7,6,8,5,9] [Just nucA,Just nucC,Just nucG,Just nucT,Nothing] ]
 
                               <*> mk_substs 0
                               <*> mk_substs 1
@@ -469,18 +486,16 @@
                    , substs3d3 = mconcat [ substs3d3 accs', substs3dd accs'] }
   where
     {-# INLINE withPair #-}
-    withPair (Ns u, Ns v) k = case pairTab `U.unsafeIndex` fromIntegral (16*u+v) of
-         j -> if j >= 0 then k j else return ()
+    withPair (NPair i) k =
+        case pairTab `U.unsafeIndex` fromIntegral i of
+            j -> when (j >= 0) (k j)
 
     !pairTab = U.replicate 256 (-1) U.//
-            [ (fromIntegral $ 16*u+v, x*4+y) | (Ns u,x) <- zip [nucsA, nucsC, nucsG, nucsT] [0,1,2,3]
-                                             , (Ns v,y) <- zip [nucsA, nucsC, nucsG, nucsT] [0,1,2,3] ]
+            [ (fromIntegral i, x*4+y) | (u,x) <- zip [nucsA, nucsC, nucsG, nucsT] [0,1,2,3]
+                                      , (v,y) <- zip [nucsA, nucsC, nucsG, nucsT] [0,1,2,3]
+                                      , let NPair i = npair u v ]
     {-# INLINE bump #-}
-#ifdef DEBUG
-    bump i v = UM.read v i >>= UM.write v i . succ
-#else
     bump i v = UM.unsafeRead v i >>= UM.unsafeWrite v i . succ
-#endif
 
     {-# INLINE withNs #-}
     withNs ns k | ns == nucsA = k 0
@@ -533,20 +548,20 @@
                                      | otherwise        = (x :-> y, U.zipWith (+) u v)
 
 
-revcom_both :: ( BamRec, FragType, U.Vector Word8, U.Vector (Nucleotides, Nucleotides) )
-            -> ( BamRec, FragType, U.Vector Word8, U.Vector (Nucleotides, Nucleotides) )
+revcom_both :: ( BamRec, FragType, U.Vector Word8, VS.Vector NPair )
+            -> ( BamRec, FragType, U.Vector Word8, VS.Vector NPair )
 revcom_both (b, ft, ref, pps)
     | isReversed b = ( b, ft, revcom_ref ref, revcom_pairs pps )
     | otherwise    = ( b, ft,            ref,              pps )
   where
-    revcom_ref   = U.reverse . U.map (\c -> if c > 3 then c else xor c 2)
-    revcom_pairs = U.reverse . U.map (compls *** compls)
+    revcom_ref   =  U.reverse .  U.map (\c -> if c > 3 then c else xor c 2)
+    revcom_pairs = VS.reverse . VS.map (\p -> npair (compls $ fst_np p) (compls $ snd_np p))
 
 
 -- | Reconstructs the alignment from reference, query, and cigar.  Only
 -- positions where the query is not gapped are produced.
-aln_from_ref :: U.Vector Word8 -> Vector_Nucs_half Nucleotides -> VS.Vector Cigar -> U.Vector NPair
-aln_from_ref ref0 qry0 cig0 = U.fromList $ step ref0 qry0 cig0
+aln_from_ref :: U.Vector Word8 -> Vector_Nucs_half Nucleotides -> VS.Vector Cigar -> VS.Vector NPair
+aln_from_ref ref0 qry0 cig0 = VS.fromList $ step ref0 qry0 cig0
   where
     step ref qry cig1
         | U.null ref || G.null qry || G.null cig1 = []
@@ -554,14 +569,14 @@
                       case G.unsafeTail cig1 of { cig ->
                       case op of {
 
-        Mat -> zipWith (\r q -> (nn r,q)) (G.toList (G.take n ref))
-                                          (G.toList (G.take n qry)) ++ step (G.drop n ref) (G.drop n qry) cig ;
-        Del ->                                                         step (G.drop n ref)           qry  cig ;
-        Ins ->    map (\q -> ( gap,  q )) (G.toList (G.take n qry)) ++ step           ref  (G.drop n qry) cig ;
-        SMa ->    map (\q -> ( gap,  q )) (G.toList (G.take n qry)) ++ step           ref  (G.drop n qry) cig ;
-        HMa ->   replicate n (gap, nucsN)                           ++ step           ref            qry  cig ;
-        Nop ->                                                         step           ref            qry  cig ;
-        Pad ->                                                         step           ref            qry  cig }}}
+        Mat -> zipWith (npair . nn) (G.toList (G.take n ref))
+                                    (G.toList (G.take n qry)) ++ step (G.drop n ref) (G.drop n qry) cig ;
+        Del ->                                                   step (G.drop n ref)           qry  cig ;
+        Ins -> map (npair gap) (G.toList (G.take n qry))      ++ step           ref  (G.drop n qry) cig ;
+        SMa -> map (npair gap) (G.toList (G.take n qry))      ++ step           ref  (G.drop n qry) cig ;
+        HMa -> replicate n (npair gap nucsN)                  ++ step           ref            qry  cig ;
+        Nop ->                                                   step           ref            qry  cig ;
+        Pad ->                                                   step           ref            qry  cig }}}
 
     nn 0 = nucsT
     nn 1 = nucsC
@@ -572,8 +587,8 @@
 
 -- | Reconstructs the alignment from query, cigar, and md.  Only
 -- positions where the query is not gapped are produced.
-alnFromMd :: Vector_Nucs_half Nucleotides -> VS.Vector Cigar -> [MdOp] -> U.Vector NPair
-alnFromMd qry0 cig0 md0 = U.fromList $ step qry0 cig0 md0
+alnFromMd :: Vector_Nucs_half Nucleotides -> VS.Vector Cigar -> [MdOp] -> VS.Vector NPair
+alnFromMd qry0 cig0 md0 = VS.fromList $ step qry0 cig0 md0
   where
     step qry cig1 md
         | G.null qry || G.null cig1 || null md = []
@@ -584,21 +599,23 @@
     step' qry op n cig (MdDel [] : md) = step' qry op n cig md
 
     step' qry Mat n cig (MdNum m : md)
-            | n <  m = map (\q -> (q,q)) (G.toList (G.take n qry)) ++ step  (G.drop n qry)           cig (MdNum (m-n) : md)
-            | n >  m = map (\q -> (q,q)) (G.toList (G.take m qry)) ++ step' (G.drop m qry) Mat (n-m) cig                md
-            | n == m = map (\q -> (q,q)) (G.toList (G.take n qry)) ++ step  (G.drop n qry)           cig                md
-    step' qry Mat n cig (MdRep c : md) =         ( c, G.head qry )  : step' (G.tail   qry) Mat (n-1) cig                md
+            | n <  m =    map twin (G.toList (G.take n qry)) ++ step  (G.drop n qry)           cig (MdNum (m-n) : md)
+            | n >  m =    map twin (G.toList (G.take m qry)) ++ step' (G.drop m qry) Mat (n-m) cig                md
+            | n == m =    map twin (G.toList (G.take n qry)) ++ step  (G.drop n qry)           cig                md
+    step' qry Mat n cig (MdRep c : md) = npair c (G.head qry) : step' (G.tail   qry) Mat (n-1) cig                md
     step'   _ Mat _   _          _     = []
 
     step' qry Del n cig (MdDel (_:ss) : md) = step' qry Del (n-1) cig (MdDel ss : md)
     step'   _ Del _   _               _     = []
 
-    step' qry Ins n cig                 md  = map ((,) gap) (G.toList (G.take n qry)) ++ step (G.drop n qry) cig md
-    step' qry SMa n cig                 md  = map ((,) gap) (G.toList (G.take n qry)) ++ step (G.drop n qry) cig md
-    step' qry HMa n cig                 md  =                replicate n (gap, nucsN) ++ step           qry  cig md
-    step' qry Nop _ cig                 md  =                                            step           qry  cig md
-    step' qry Pad _ cig                 md  =                                            step           qry  cig md
+    step' qry Ins n cig                 md  = map (npair gap) (G.toList (G.take n qry)) ++ step (G.drop n qry) cig md
+    step' qry SMa n cig                 md  = map (npair gap) (G.toList (G.take n qry)) ++ step (G.drop n qry) cig md
+    step' qry HMa n cig                 md  =             replicate n (npair gap nucsN) ++ step           qry  cig md
+    step' qry Nop _ cig                 md  =                                              step           qry  cig md
+    step' qry Pad _ cig                 md  =                                              step           qry  cig md
 
+    twin q = npair q q
+
 -- | Number of mismatches allowed by BWA.
 -- @bwa_cal_maxdiff thresh len@ returns the number of mismatches
 -- @bwa aln -n $tresh@ would allow in a read of length @len@.  For
@@ -620,18 +637,6 @@
 --      return 2;
 --   }
 -- @
---      double sum, y = 1.0;
---      int k, x = 1;
---      for (k = 1, sum = elambda; k < 1000; ++k) {
---          y *= l * err;
---          x *= k;
---          sum += elambda * y / x;
---          if (1.0 - sum < thres) return k;
---      }
---      return 2;
---   }
--- @
---
 
 bwa_cal_maxdiff :: Double -> Int -> Int
 bwa_cal_maxdiff thresh len = k_fin-1
diff --git a/src/Bio/Bam.hs b/src/Bio/Bam.hs
--- a/src/Bio/Bam.hs
+++ b/src/Bio/Bam.hs
@@ -1,3 +1,5 @@
+-- | Umbrella module for most of what's under 'Bio.Bam'.
+
 module Bio.Bam (
     module Bio.Bam.Fastq,
     module Bio.Bam.Filter,
@@ -19,6 +21,4 @@
 import Bio.Bam.Trim
 import Bio.Bam.Writer
 import Bio.Iteratee
-
--- ^ Umbrella module for most of what's under 'Bio.Bam'.
 
diff --git a/src/Bio/Bam/Evan.hs b/src/Bio/Bam/Evan.hs
--- a/src/Bio/Bam/Evan.hs
+++ b/src/Bio/Bam/Evan.hs
@@ -14,7 +14,7 @@
 
 -- | Fixes abuse of flags valued 0x800 and 0x1000.  We used them for
 -- low quality and low complexity, but they have since been redefined.
--- If set, we clear them and store them into the ZD field.  Also fixes
+-- If set, we clear them and store them into the ZQ field.  Also fixes
 -- abuse of the combination of the paired, 1st mate and 2nd mate flags
 -- used to indicate merging or trimming.  These are canonicalized and
 -- stored into the FF field.  This function is unsafe on BAM files of
@@ -38,9 +38,9 @@
         is_trimmed = flags' .&. (flagPaired .|. flagFirstMate .|. flagSecondMate) == flagSecondMate
         newflags = (if is_merged then eflagMerged else 0) .|. (if is_trimmed then eflagTrimmed else 0)
 
-        -- Extended flags, renamed to avoid collision with BWA Goes like this:  if FF is there, use
-        -- it.  Else check if XF is there _and_is_numeric_.  If so, use it and remove it, and set FF
-        -- instead.  Else use 0 and leave it alone.  Note that this solves the collision with BWA,
+        -- Extended flags, renamed to avoid collision with BWA.  Goes like this:  if FF is there, use
+        -- it.  Else check if XF is there __and is numeric__.  If so, use it, remove it, and set FF
+        -- instead.  Else use 0 and leave it alone.  Note that this resolves the collision with BWA,
         -- since BWA puts a character there, not an int.
         cleaned_exts = case (lookup "FF" (b_exts b), lookup "XF" (b_exts b)) of
                 ( Just (Int i), _ ) -> updateE "FF" (Int (i .|. newflags))                (b_exts b)
@@ -51,7 +51,7 @@
 
 -- | Fixes typical inconsistencies produced by Bwa: sometimes, 'mate unmapped' should be set, and we
 -- can see it, because we match the mate's coordinates.  Sometimes 'properly paired' should not be
--- set, because one mate in unmapped.  This function is generally safe, but needs to be called only
+-- set, because one mate is unmapped.  This function is generally safe, but needs to be called only
 -- on the output of affected (older?) versions of Bwa.
 fixupBwaFlags :: BamRec -> BamRec
 fixupBwaFlags b = b { b_flag = fixPP $ b_flag b .|. if mu then flagMateUnmapped else 0 }
@@ -65,6 +65,7 @@
         -- If either mate is unmapped, remove "properly paired".
         fixPP f | f .&. (flagUnmapped .|. flagMateUnmapped) == 0 = f
                 | otherwise = f .&. complement flagProperlyPaired
+
 
 -- | Removes syntactic warts from old read names or the read names used
 -- in FastQ files.
diff --git a/src/Bio/Bam/Fastq.hs b/src/Bio/Bam/Fastq.hs
--- a/src/Bio/Bam/Fastq.hs
+++ b/src/Bio/Bam/Fastq.hs
@@ -1,3 +1,8 @@
+-- | Parser for @FastA/FastQ@, 'Iteratee' style, based on
+-- "Data.Attoparsec", and written such that it is compatible with module
+-- 'Bio.Bam'.  This gives import of @FastA/FastQ@ while respecting some
+-- local (to MPI EVAN) conventions.
+
 module Bio.Bam.Fastq ( parseFastq, parseFastq', parseFastqCassava ) where
 
 import Bio.Bam.Header
@@ -11,19 +16,13 @@
 import qualified Data.ByteString.Char8              as S
 import qualified Data.Vector.Generic                as V
 
--- ^ Parser for @FastA/FastQ@, 'Iteratee' style, based on
--- "Data.Attoparsec", and written such that it is compatible with module
--- 'Bio.Bam'.  This gives import of @FastA/FastQ@ while respecting some
--- local conventions.
-
 -- | Reader for DNA (not protein) sequences in FastA and FastQ.  We read
 -- everything vaguely looking like FastA or FastQ, then shoehorn it into
 -- a BAM record.  We strive to extract information following more or
--- less established conventions from the header, but we won't support
--- everything under the sun.  The recognized syntactical warts are
--- converted into appropriate flags and removed.  Only the canonical
--- variant of FastQ is supported (qualities stored as raw bytes with
--- base 33).
+-- less established conventions from the header, but don't aim for
+-- completeness.  The recognized syntactical warts are converted into
+-- appropriate flags and removed.  Only the canonical variant of FastQ
+-- is supported (qualities stored as raw bytes with offset 33).
 --
 -- Supported additional conventions:
 --
@@ -36,18 +35,13 @@
 --
 -- * A name prefix of @T_@ flags the sequence as unpaired and trimmed
 --
--- * A name prefix of @C_@, either before or after any of the other
+-- * A name prefix of @C_@, optionally before or after any of the other
 --   prefixes, is turned into the extra flag @XP:i:-1@ (result of
 --   duplicate removal with unknown duplicate count).
 --
 -- * A collection of tags separated from the name by an octothorpe is
 --   removed and put into the fields @XI@ and @XJ@ as text.
 --
--- * In 'parseFastqCassava' only, if the first word of the description
---   has at least four colon separated subfields, the first if used to
---   flag first/second mate, the second is the \"QC failed\" flag, and
---   the fourth is the index sequence.
---
 -- Everything before the first sequence header is ignored.  Headers can
 -- start with @\>@ or @\@@, we treat both equally.  The first word of
 -- the header becomes the read name, the remainder of the header is
@@ -63,6 +57,13 @@
 parseFastq :: Monad m => Enumeratee Bytes [ BamRec ] m a
 parseFastq = parseFastq' (const id)
 
+-- | Like 'parseFastq', but also
+--
+-- * If the first word of the description has at least four colon
+--   separated subfields, the first is used to flag first/second mate,
+--   the second is the \"QC failed\" flag, and the fourth is the index
+--   sequence.
+
 parseFastqCassava :: Monad m => Enumeratee Bytes [ BamRec ] m a
 parseFastqCassava = parseFastq' (pdesc . S.split ':' . S.takeWhile (' ' /=))
   where
@@ -78,7 +79,6 @@
 -- which can modify the parsed record.  Note that the quality field can
 -- end up empty.
 
-{-# WARNING parseFastq' "parseFastq' no longer removes syntactic warts!" #-}
 parseFastq' :: Monad m => ( Bytes -> BamRec -> BamRec ) -> Enumeratee Bytes [ BamRec ] m a
 parseFastq' descr it = do skipJunk ; convStream (parserToIteratee $ (:[]) <$> pRec) it
   where
diff --git a/src/Bio/Bam/Filter.hs b/src/Bio/Bam/Filter.hs
--- a/src/Bio/Bam/Filter.hs
+++ b/src/Bio/Bam/Filter.hs
@@ -1,3 +1,5 @@
+-- | Quality filters adapted from prehistoric pipeline.
+
 module Bio.Bam.Filter (
     filterPairs, QualFilter,
     complexSimple, complexEntropy,
@@ -13,8 +15,6 @@
 import Prelude
 
 import qualified Data.Vector.Generic as V
-
--- ^ Quality filters adapted from old pipeline.
 
 -- | A filter/transformation applied to pairs of reads.  We supply a
 -- predicate to be applied to single reads and one to be applied to
diff --git a/src/Bio/Bam/Header.hs b/src/Bio/Bam/Header.hs
--- a/src/Bio/Bam/Header.hs
+++ b/src/Bio/Bam/Header.hs
@@ -317,13 +317,14 @@
 
 
 -- | Compares two sequence names the way samtools does.
--- samtools sorts by "strnum_cmp":
--- . if both strings start with a digit, parse the initial
+-- samtools sorts by \"strnum_cmp\":
+--
+-- * if both strings start with a digit, parse the initial
 --   sequence of digits and compare numerically, if equal,
 --   continue behind the numbers
--- . else compare the first characters (possibly NUL), if equal
+-- * else compare the first characters (possibly NUL), if equal
 --   continue behind them
--- . else both strings ended and the shorter one counts as
+-- * else both strings ended and the shorter one counts as
 --   smaller (and that part is stupid)
 
 compareNames :: Seqid -> Seqid -> Ordering
diff --git a/src/Bio/Bam/Index.hs b/src/Bio/Bam/Index.hs
--- a/src/Bio/Bam/Index.hs
+++ b/src/Bio/Bam/Index.hs
@@ -9,8 +9,7 @@
     eneeBamRefseq,
     eneeBamSubseq,
     eneeBamRegions,
-    eneeBamUnaligned,
-    subsampleBam
+    eneeBamUnaligned
 ) where
 
 import Bio.Bam.Header
@@ -18,15 +17,12 @@
 import Bio.Bam.Rec
 import Bio.Bam.Regions              ( Region(..), Subsequence(..) )
 import Bio.Iteratee
-import Bio.Iteratee.Bgzf
 import Bio.Prelude
 import System.Directory             ( doesFileExist )
 import System.FilePath              ( dropExtension, takeExtension, (<.>) )
-import System.Random                ( randomRIO )
 
 import qualified Bio.Bam.Regions                as R
-import qualified Control.Exception              as E
-import qualified Data.IntMap                    as M
+import qualified Data.IntMap.Strict             as M
 import qualified Data.ByteString                as B
 import qualified Data.Vector                    as V
 import qualified Data.Vector.Mutable            as W
@@ -115,7 +111,7 @@
     , let boff' = max boff cpt
     , boff' < eoff ]
   where
-    !cpt = maybe 0 snd $ lookupLE beg cpts
+    !cpt = maybe 0 snd $ M.lookupLE beg cpts
 
 -- list of bins for given range of coordinates, from Heng's horrible code
 binList :: BamIndex a -> Int -> Int -> [Int]
@@ -333,51 +329,3 @@
 eneeBamRegions :: Monad m => BamIndex b -> [R.Region] -> Enumeratee [BamRaw] [BamRaw] m a
 eneeBamRegions bi = foldr ((>=>) . uncurry (eneeBamSubseq bi)) return . R.toList . R.fromList
 
-
-lookupLE :: M.Key -> M.IntMap a -> Maybe (M.Key, a)
-lookupLE k m = case ma of
-    Just a              -> Just (k,a)
-    Nothing | M.null m1 -> Nothing
-            | otherwise -> Just $ M.findMax m1
-  where (m1,ma,_) = M.splitLookup k m
-
-
--- | Subsample randomly from a BAM file.  If an index exists, this
--- produces an infinite stream taken from random locations in the file.
---
--- XXX It would be cool if we could subsample from multiple BAM files.
--- It's a bit annoying to code: we'd probably read the indices up front,
--- estimate how many reads we'd find in each file, then open them
--- recursively to form a monad stack where the merging function has to
--- select randomly where to read from.  Hm.
-
-subsampleBam :: (MonadIO m, MonadMask m) => FilePath -> Enumerator' BamMeta [BamRaw] m b
-subsampleBam fp o = liftIO (E.try (readBamIndex fp)) >>= subsam
-  where
-    -- no index, so just stream
-    subsam (Left e) = enumFile defaultBufSize fp >=> run $
-                      joinI $ decompressBgzfBlocks $
-                      joinI $ decodeBam $ \hdr ->
-                      takeWhileE (isValidRefseq . b_rname . unpackBam) (o hdr)
-                            `const` (e::E.SomeException)
-
-    -- with index: chose random bins and read from them
-    subsam (Right bix) = withFileFd fp $ \fd -> do
-                         hdr <- enumFdRandom defaultBufSize fd >=> run $
-                                joinI $ decompressBgzfBlocks' 1 $
-                                joinI $ decodeBam return
-                         go fd (o hdr)
-      where
-        !ckpts = U.fromList . V.foldr ((++) . M.elems) [] $ refseq_ckpoints bix
-
-        go fd o1 = enumCheckIfDone o1 >>= go' fd
-
-        go'  _ (True,  o2) = return o2
-        go' fd (False, o2) = do i <- liftIO $ randomRIO (0, U.length ckpts -1)
-                                enum fd i o2 >>= go fd
-
-        enum fd i = enumFdRandom defaultBufSize fd               $=
-                    decompressBgzfBlocks' 1                      $=
-                    (\it -> do seek . fromIntegral $ ckpts U.! i
-                               convStream getBamRaw it)          $=
-                    takeStream 512
diff --git a/src/Bio/Bam/Pileup.hs b/src/Bio/Bam/Pileup.hs
--- a/src/Bio/Bam/Pileup.hs
+++ b/src/Bio/Bam/Pileup.hs
@@ -1,4 +1,15 @@
 {-# LANGUAGE Rank2Types, DeriveGeneric #-}
+
+-- | 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.
+
 module Bio.Bam.Pileup where
 
 import Bio.Bam.Header
@@ -10,63 +21,12 @@
 import qualified Data.Vector.Generic    as V
 import qualified Data.Vector.Unboxed    as U
 
--- ^ Genotype Calling:  like Samtools(?), but for aDNA
---
--- The goal for this module is to call haploid and diploid single
--- nucleotide variants the best way we can, including support for aDNA.
--- Indel calling is out of scope, we only do it "on the side".
---
--- The cleanest way to call genotypes under all circumstances is
--- probably the /Dindel/ approach:  define candidate haplotypes, align
--- each read to each haplotype, then call the likely haplotypes with a
--- quality derived from the quality scores.  This approach neatly
--- integrates indel calling with ancient DNA and makes a separate indel
--- realigner redundant.  However, it's rather expensive in that it
--- requires inclusion of an aligner, and we'd need an aligner that is
--- compatible with the chosen error model, which might be hard.
---
--- Here we'll take a short cut:  We do not really call indels.  Instead,
--- these variants are collected and are assigned an affine score.  This
--- works best if indels are 'left-aligned' first.  In theory, one indel
--- variant could be another indel variant with a sequencing error---we
--- ignore that possibility for the most part.  Once indels are taken
--- care off, SNVs are treated separately as independent columns of the
--- pileup.
---
--- Regarding the error model, there's a choice between /samtools/ or the
--- naive model everybody else (GATK, Rasmus Nielsen, etc.) uses.  Naive
--- is easy to marry to aDNA, samtools is (probably) better.  Either way,
--- we introduce a number of parameters (@eta@ and @kappa@ for
--- /samtools/, @lambda@, @delta@, @delta_ss@ for /Johnson/).  Running a
--- maximum likehood fit for those may be valuable.  It would be cool, if
--- we could do that without rerunning the complete genotype caller, but
--- it's not a priority.
---
--- So, outline of the genotype caller:  We read BAM (minimally
--- filtering; general filtering is somebody else's problem, but we might
--- want to split by read group).  We will scan each read's CIGAR line in
--- concert with the sequence and effective quality.  Effective quality
--- is the lowest available quality score of QUAL, MAPQ, and BQ.  For
--- aDNA calling, the base is transformed into four likelihoods based on
--- the aDNA substitution matrix.
---
--- So, either way, we need something like "pileup", where indel variants
--- are collected as they are (any length), while matches are piled up.
---
--- Regarding output, we certainly don't want to write VCF or BCF.  (No
--- VCF because it's ugly, no BCF, because the tool support is
--- non-existent.)  It will definitely be something binary.  For the GL
--- values, small floating point formats may make sense: half-precision
--- floating point's representable range would be 6.1E-5 to 6.5E+5, 0.4.4
--- minifloat from Bio.Util goes from 0 to 63488.
-
-
 -- | 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.
-data PrimChunks = Seek Int PrimBase                                   -- ^ skip to position (at start or after N operation)
+data PrimChunks = Seek {-# UNPACK #-} !Int PrimBase                   -- ^ skip to position (at start or after N operation)
                 | Indel [Nucleotides] [DamagedBase] PrimBase          -- ^ observed deletion and insertion between two bases
                 | EndOfRead                                           -- ^ nothing anymore
   deriving Show
@@ -109,8 +69,8 @@
 -- | 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
--- appropriate.  We also do apply a substitution matrix to each base,
--- which must be supplied along with the read.
+-- needed.  We also apply a substitution matrix to each base, which must
+-- be supplied along with the read.
 {-# INLINE decompose #-}
 decompose :: DmgToken -> BamRaw -> [PosPrimChunks]
 decompose dtok br =
@@ -118,7 +78,7 @@
     then [] else [(b_rname, b_pos, isReversed b, pchunks)]
   where
     b@BamRec{..} = unpackBam br
-    pchunks = firstBase b_pos 0 0 (maybe [] id $ getMd b)
+    pchunks = firstBase b_pos 0 0 (fromMaybe [] $ getMd b)
 
     !max_cig = V.length b_cigar
     !max_seq = V.length b_seq
@@ -129,27 +89,14 @@
     -- and BAQ.  If QUAL is invalid, we replace it (arbitrarily) with
     -- 23 (assuming a rather conservative error rate of ~0.5%), BAQ is
     -- added to QUAL, and MAPQ is an upper limit for effective quality.
-    get_seq :: Int -> Nucleotides -> DamagedBase
-    get_seq i = case b_seq `V.unsafeIndex` i of                                 -- nucleotide
-            n | n == nucsA -> DB nucA qe dtok dmg
-              | n == nucsC -> DB nucC qe dtok dmg
-              | n == nucsG -> DB nucG qe dtok dmg
-              | n == nucsT -> DB nucT qe dtok dmg
-              | otherwise  -> DB nucA (Q 0) dtok dmg
-      where
-        !q   = case b_qual `V.unsafeIndex` i of Q 0xff -> Q 30 ; x -> x         -- quality; invalid (0xff) becomes 30
-        !q'  | i >= B.length baq = q                                            -- no BAQ available
-             | otherwise = Q (unQ q + (B.index baq i - 64))                     -- else correct for BAQ
-        !qe  = min q' b_mapq                                                    -- use MAPQ as upper limit
-        !dmg = if i+i > max_seq then i-max_seq else i
 
-    get_seq' :: Int -> DamagedBase
-    get_seq' i = case b_seq `V.unsafeIndex` i of                                -- nucleotide
-            n | n == nucsA -> DB nucA qe dtok dmg nucsA
-              | n == nucsC -> DB nucC qe dtok dmg nucsC
-              | n == nucsG -> DB nucG qe dtok dmg nucsG
-              | n == nucsT -> DB nucT qe dtok dmg nucsT
-              | otherwise  -> DB nucA (Q 0) dtok dmg n
+    get_seq :: Int -> (Nucleotides -> Nucleotides) -> DamagedBase
+    get_seq i f = case b_seq `V.unsafeIndex` i of                                -- nucleotide
+            n | n == nucsA -> DB nucA qe dtok dmg (f n)
+              | n == nucsC -> DB nucC qe dtok dmg (f n)
+              | n == nucsG -> DB nucG qe dtok dmg (f n)
+              | n == nucsT -> DB nucT qe dtok dmg (f n)
+              | otherwise  -> DB nucA (Q 0) dtok dmg (f n)
       where
         !q   = case b_qual `V.unsafeIndex` i of Q 0xff -> Q 30 ; x -> x         -- quality; invalid (0xff) becomes 30
         !q'  | i >= B.length baq = q                                            -- no BAQ available
@@ -176,7 +123,7 @@
       where
         -- We have to treat (MdNum 0), because samtools actually
         -- generates(!) it all over the place and if not handled as a
-        -- special case, it looks like an incinsistend MD field.
+        -- special case, it looks like an inconsistent MD field.
         drop_del n (MdDel ns : mds')
             | n < length ns = MdDel (drop n ns) : mds'
             | n > length ns = drop_del (n - length ns) mds'
@@ -194,11 +141,11 @@
     nextBase !wt !pos !is !ic !io mds = case mds of
         MdNum   0 : mds' -> nextBase wt pos is ic io mds'
         MdDel  [] : mds' -> nextBase wt pos is ic io mds'
-        MdNum   1 : mds' -> nextBase' (get_seq' is      ) mds'
-        MdNum   n : mds' -> nextBase' (get_seq' is      ) (MdNum (n-1) : mds')
-        MdRep ref : mds' -> nextBase' (get_seq  is ref  ) mds'
-        MdDel   _ : _    -> nextBase' (get_seq  is nucsN) mds
-        [              ] -> nextBase' (get_seq  is nucsN) [ ]
+        MdNum   1 : mds' -> nextBase' (get_seq is id   ) mds'
+        MdNum   n : mds' -> nextBase' (get_seq is id   ) (MdNum (n-1) : mds')
+        MdRep ref : mds' -> nextBase' (get_seq is $ const ref  ) mds'
+        MdDel   _ : _    -> nextBase' (get_seq is $ const nucsN) mds
+        [              ] -> nextBase' (get_seq is $ const nucsN) [ ]
       where
         nextBase' ref mds' = Base wt ref b_mapq $ nextIndel  [] [] (pos+1) (is+1) ic (io+1) mds'
 
@@ -223,11 +170,9 @@
             Mat :* cl | io == cl  -> nextIndel  ins     del   pos     is  (ic+1) 0 mds
                       | otherwise -> indel del out $ nextBase (length del) pos is ic io mds -- ends up generating a 'Base'
       where
-        indel d o k = rlist o `seq` Indel d o k
+        indel d o k = foldr seq (Indel d o k) o
         out    = concat $ reverse ins
-        isq cl = [ get_seq i gap | i <- [is..is+cl-1] ] : ins
-        rlist [    ] = ()
-        rlist (a:as) = a `seq` rlist as
+        isq cl = [ get_seq i $ const gap | i <- [is..is+cl-1] ] : ins
 
         -- We have to treat (MdNum 0), because samtools actually
         -- generates(!) it all over the place and if not handled as a
@@ -243,10 +188,11 @@
 -- fitlering (so not very useful), but we keep them because it's easy to
 -- track them.
 
-data CallStats = CallStats { read_depth       :: {-# UNPACK #-} !Int       -- number of contributing reads
-                           , reads_mapq0      :: {-# UNPACK #-} !Int       -- number of (non-)contributing reads with MAPQ==0
-                           , sum_mapq         :: {-# UNPACK #-} !Int       -- sum of map qualities of contributing reads
-                           , sum_mapq_squared :: {-# UNPACK #-} !Int }     -- sum of squared map qualities of contributing reads
+data CallStats = CallStats
+    { read_depth       :: {-# UNPACK #-} !Int       -- number of contributing reads
+    , reads_mapq0      :: {-# UNPACK #-} !Int       -- number of (non-)contributing reads with MAPQ==0
+    , sum_mapq         :: {-# UNPACK #-} !Int       -- sum of map qualities of contributing reads
+    , sum_mapq_squared :: {-# UNPACK #-} !Int }     -- sum of squared map qualities of contributing reads
   deriving (Show, Eq, Generic)
 
 instance Monoid CallStats where
@@ -259,22 +205,6 @@
                             , sum_mapq         = sum_mapq x + sum_mapq y
                             , sum_mapq_squared = sum_mapq_squared x + sum_mapq_squared y }
 
--- | Genotype likelihood values.  A variant call consists of a position,
--- some measure of qualities, genotype likelihood values, and a
--- representation of variants.  A note about the GL values:  @VCF@ would
--- normalize them so that the smallest one becomes zero.  We do not do
--- that here, since we might want to compare raw values for a model
--- test.  We also store them in a 'Double' to make arithmetics easier.
--- Normalization is appropriate when converting to @VCF@.
---
--- If GL is given, we follow the same order used in VCF:
--- \"the ordering of genotypes for the likelihoods is given by:
--- F(j/k) = (k*(k+1)/2)+j.  In other words, for biallelic sites the
--- ordering is: AA,AB,BB; for triallelic sites the ordering is:
--- AA,AB,BB,AC,BC,CC, etc.\"
-
-type GL = U.Vector Prob
-
 newtype V_Nuc  = V_Nuc  (U.Vector Nucleotide)  deriving (Eq, Ord, Show)
 newtype V_Nucs = V_Nucs (U.Vector Nucleotides) deriving (Eq, Ord, Show)
 
@@ -292,10 +222,10 @@
 type IndelPile = [( Qual, ([Nucleotides], [DamagedBase]) )]   -- a list of indel variants
 
 -- | Running pileup results in a series of piles.  A 'Pile' has the
--- basic statistics of a 'VarCall', but no GL values and a pristine list
--- of variants instead of a proper call.  We emit one pile with two
--- 'BasePile's (one for each strand) and one 'IndelPile' (the one
--- immediately following) at a time.
+-- 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 'BasePile's (one for each strand) and one 'IndelPile' (the
+-- one immediately following) at a time.
 
 data Pile' a b = Pile { p_refseq     :: {-# UNPACK #-} !Refseq
                       , p_pos        :: {-# UNPACK #-} !Int
@@ -308,27 +238,6 @@
 -- | Raw pile.  Bases and indels are piled separately on forward and
 -- backward strands.
 type Pile = Pile' (BasePile, BasePile) (IndelPile, IndelPile)
-
--- | Simple single population model.  'prob_div' is the fraction of
--- homozygous divergent sites, 'prob_het' is the fraction of
--- heterozygous variant sites among sites that are not homozygous
--- divergent.
-data SinglePop = SinglePop { prob_div :: !Double, prob_het :: !Double }
-
--- | Computes posterior  genotype probabilities from likelihoods under
--- the 'SinglePop' model.
-{-# INLINE single_pop_posterior #-}
-single_pop_posterior :: ( U.Unbox a, Ord a, Floating a )
-                     => SinglePop -> Int -> U.Vector (Prob' a) -> U.Vector (Prob' a)
-single_pop_posterior SinglePop{..} refix lks = U.map (/ U.sum v) v
-  where
-    priors = U.replicate (U.length lks)
-                         ((1/3) * prob_het  * (1-prob_div))                                 -- hets
-                U.// [ (refix, (1-prob_het) * (1-prob_div)) ]                               -- ref
-                U.// [ (i,               (1/3) * prob_div ) | i <- hixes, i /= refix ]      -- homs
-
-    v = U.zipWith (\l p -> l * toProb (realToFrac p)) lks priors
-    hixes = takeWhile (< U.length lks) $ scanl (+) 0 [2..]
 
 -- | The pileup enumeratee takes 'BamRaw's, decomposes them, interleaves
 -- the pieces appropriately, and generates 'Pile's.  The output will
diff --git a/src/Bio/Bam/Reader.hs b/src/Bio/Bam/Reader.hs
--- a/src/Bio/Bam/Reader.hs
+++ b/src/Bio/Bam/Reader.hs
@@ -1,3 +1,12 @@
+-- | Parsers for BAM and SAM.
+--
+-- TONOTDO:
+--
+-- * Reader for gzipped\/bzipped\/bgzf'ed SAM.  Storing SAM is a bad idea,
+--   so why would anyone ever want to compress, much less index it?
+-- * Proper support for the "=" symbol.  It's completely alien to the
+--   ususal representation of sequences.
+
 module Bio.Bam.Reader (
     Block(..),
     decompressBgzfBlocks,
@@ -52,21 +61,13 @@
 import qualified Data.Vector.Storable               as VS
 import qualified Data.Vector.Unboxed                as U
 
--- ^ Parsers for BAM and SAM.  We employ an @Iteratee@ interface, and we
--- strive to support everything possible in BAM.  The implementation of
--- nucleotides is somewhat lacking:  the "=" symbol is not understood.
---
--- TONOTDO:
--- - Reader for gzipped/bzipped/bgzf'ed SAM.  Storing SAM is a bad idea,
---   so why would anyone ever want to compress, much less index it?
-
 type BamrawEnumeratee m b = Enumeratee' BamMeta Bytes [BamRaw] m b
 type BamEnumeratee m b = Enumeratee' BamMeta Bytes [BamRec] m b
 
 isBamOrSam :: MonadIO m => Iteratee Bytes m (BamEnumeratee m a)
 isBamOrSam = maybe decodeSam wrap `liftM` isBam
   where
-    wrap enee it' = enee (\hdr -> mapStream unpackBam (it' hdr)) >>= lift . run
+    wrap enee it' = enee (mapStream unpackBam . it') >>= lift . run
 
 
 -- | Checks if a file contains BAM in any of the common forms,
@@ -94,16 +95,17 @@
     decodeSamLoop (meta_refs meta) (inner meta)
 
 decodeSamLoop :: Monad m => Refs -> Enumeratee [Bytes] [BamRec] m a
-decodeSamLoop refs inner = convStream (liftI parse_record) inner
-  where !refs' = M.fromList $ zip [ nm | BamSQ { sq_name = nm } <- F.toList refs ] [toEnum 0..]
-        ref x = M.lookupDefault invalidRefseq x refs'
+decodeSamLoop refs = convStream (liftI parse_record)
+  where
+    !refs' = M.fromList $ zip [ nm | BamSQ { sq_name = nm } <- F.toList refs ] [toEnum 0..]
+    ref x = M.lookupDefault invalidRefseq x refs'
 
-        parse_record (EOF x) = icont parse_record x
-        parse_record (Chunk []) = liftI parse_record
-        parse_record (Chunk (l:ls)) | "@" `S.isPrefixOf` l = parse_record (Chunk ls)
-        parse_record (Chunk (l:ls)) = case P.parseOnly (parseSamRec ref) l of
-            Right  r -> idone [r] (Chunk ls)
-            Left err -> icont parse_record (Just $ iterStrExc $ err ++ ", " ++ show l)
+    parse_record (EOF x) = icont parse_record x
+    parse_record (Chunk []) = liftI parse_record
+    parse_record (Chunk (l:ls)) | "@" `S.isPrefixOf` l = parse_record (Chunk ls)
+    parse_record (Chunk (l:ls)) = case P.parseOnly (parseSamRec ref) l of
+        Right  r -> idone [r] (Chunk ls)
+        Left err -> icont parse_record (Just $ iterStrExc $ err ++ ", " ++ show l)
 
 -- | Parser for SAM that doesn't look for a header.  Has the advantage
 -- that it doesn't stall on a pipe that never delivers data.  Has the
@@ -119,7 +121,7 @@
                   <*> snum <*> sequ <*> quals <*> exts <*> pure 0
   where
     sep      = P.endOfInput <|> () <$ P.char '\t'
-    word     = P.takeTill ((==) '\t') <* sep
+    word     = P.takeTill ('\t' ==) <* sep
     num      = P.decimal <* sep
     num'     = P.decimal <* sep
     snum     = P.signed P.decimal <* sep
@@ -143,7 +145,7 @@
 
     value    = P.char 'A' *> P.char ':' *> (Char <$>               anyWord8) <|>
                P.char 'i' *> P.char ':' *> (Int  <$>     P.signed P.decimal) <|>
-               P.char 'Z' *> P.char ':' *> (Text <$> P.takeTill ((==) '\t')) <|>
+               P.char 'Z' *> P.char ':' *> (Text <$>   P.takeTill ('\t' ==)) <|>
                P.char 'H' *> P.char ':' *> (Bin  <$>               hexarray) <|>
                P.char 'f' *> P.char ':' *> (Float . realToFrac <$> P.double) <|>
                P.char 'B' *> P.char ':' *> (
@@ -186,7 +188,7 @@
 -- @Iteratee@.)
 isBgzfBam  = do b <- isBgzf
                 k <- if b then joinI $ enumInflate GZip defaultDecompressParams isPlainBam else return Nothing
-                return $ (\_ -> (joinI . decompressBgzfBlocks . decodeBam)) `fmap` k
+                return $ const (joinI . decompressBgzfBlocks . decodeBam) `fmap` k
 
 isGzipBam  = do b <- isGzip
                 k <- if b then joinI $ enumInflate GZip defaultDecompressParams isPlainBam else return Nothing
@@ -212,14 +214,12 @@
 concatDefaultInputs it0 = liftIO getArgs >>= \fs -> concatInputs fs it0
 
 concatInputs :: (MonadIO m, MonadMask m) => [FilePath] -> Enumerator' BamMeta [BamRaw] m a
-concatInputs [        ] = \k -> enumFd defaultBufSize stdInput (decodeAnyBam k) >>= run
-concatInputs (fp0:fps0) = \k -> enum1 fp0 k >>= go fps0
+concatInputs [       ] = enumFd defaultBufSize stdInput . decodeAnyBam >=> run
+concatInputs (fp0:fps) = enum1 fp0 >=> go fps
   where
-    enum1 "-" k1 = enumFd   defaultBufSize stdInput (decodeAnyBam k1) >>= run
-    enum1  fp k1 = enumFile defaultBufSize fp       (decodeAnyBam k1) >>= run
-
-    go [       ] = return
-    go (fp1:fps) = enum1 fp1 . const >=> go fps
+    go = foldr (\fp1 -> (>=>) (enum1 fp1 . const)) return
+    enum1 "-" = enumFd   defaultBufSize stdInput . decodeAnyBam >=> run
+    enum1  fp = enumFile defaultBufSize fp       . decodeAnyBam >=> run
 
 mergeDefaultInputs :: (MonadIO m, MonadMask m)
     => (BamMeta -> Enumeratee [BamRaw] [BamRaw] (Iteratee [BamRaw] m) a)
@@ -270,7 +270,7 @@
                                    nm <- endianRead4 LSB >>= iGetString . fromIntegral
                                    ln <- endianRead4 LSB
                                    return $! acc |> BamSQ (S.init nm) (fromIntegral ln) []
-                             ) Z.empty $ [1..nref]
+                             ) Z.empty [1..nref]
 
     -- Need to merge information from header into actual reference list.
     -- The latter is the authoritative source for the *order* of the
@@ -282,7 +282,7 @@
         in meta { meta_refs = fmap (\s -> maybe s (mmerge' s) (M.lookup (sq_name s) tbl)) refs }
 
     mmerge' l r | sq_length l == sq_length r = l { sq_other_shit = sq_other_shit l ++ sq_other_shit r }
-               | otherwise                  = l -- contradiction in header, but we'll just ignore it
+                | otherwise                  = l -- contradiction in header, but we'll just ignore it
 
 
 {-# INLINE getBamRaw #-}
diff --git a/src/Bio/Bam/Rec.hs b/src/Bio/Bam/Rec.hs
--- a/src/Bio/Bam/Rec.hs
+++ b/src/Bio/Bam/Rec.hs
@@ -6,13 +6,6 @@
 -- do not have support for ambiguity codes, and the "=" symbol is not
 -- understood.
 
--- TODO:
--- - Automatic creation of some kind of index.  If possible, this should
---   be the standard index for sorted BAM and/or the newer CSI format.
---   Optionally, a block index for slicing of large files, even unsorted
---   ones.  Maybe an index by name and an index for group-sorted files.
---   Sensible indices should be generated whenever a file is written.
-
 module Bio.Bam.Rec (
     BamRaw,
     bamRaw,
@@ -103,7 +96,7 @@
         where
             w = fromEnum op .|. shiftL num 4
 
--- | extracts the aligned length from a cigar line
+-- | Extracts the aligned length from a cigar line.
 -- This gives the length of an alignment as measured on the reference,
 -- which is different from the length on the query or the length of the
 -- alignment.
diff --git a/src/Bio/Bam/Regions.hs b/src/Bio/Bam/Regions.hs
--- a/src/Bio/Bam/Regions.hs
+++ b/src/Bio/Bam/Regions.hs
@@ -1,10 +1,10 @@
 module Bio.Bam.Regions where
 
 import Bio.Bam.Header ( Refseq(..) )
-import Data.List ( foldl' )
-import qualified Data.IntMap as IM
-import Prelude
+import Bio.Prelude
 
+import qualified Data.IntMap.Strict as IM
+
 data Region = Region { refseq :: !Refseq, start :: !Int, end :: !Int }
   deriving (Eq, Ord, Show)
 
@@ -32,7 +32,7 @@
 addInt :: Int -> Int -> Subsequence -> Subsequence
 addInt b e (Subsequence m0) = Subsequence $ merge_into b e m0
   where
-    merge_into x y m = case lookupLT y m of
+    merge_into x y m = case IM.lookupLT y m of
         Just (u,v) | x < u && y <= v -> merge_into x v $ IM.delete u m    -- extend to the left
                    | x < u           -> merge_into x y $ IM.delete u m    -- subsume
                    | y <= v          -> m                                 -- subsumed
@@ -40,11 +40,6 @@
         _                            -> IM.insert  x y m                  -- no overlap
 
 overlaps :: Int -> Int -> Subsequence -> Bool
-overlaps b e (Subsequence m) = case lookupLT e m of
+overlaps b e (Subsequence m) = case IM.lookupLT e m of
         Just (_,v) -> b < v
         Nothing    -> False
-
-lookupLT :: IM.Key -> IM.IntMap a -> Maybe (IM.Key, a)
-lookupLT k m | IM.null m1 = Nothing
-             | otherwise  = Just $ IM.findMax m1
-  where (m1,_) = IM.split k m
diff --git a/src/Bio/Bam/Rmdup.hs b/src/Bio/Bam/Rmdup.hs
--- a/src/Bio/Bam/Rmdup.hs
+++ b/src/Bio/Bam/Rmdup.hs
@@ -12,7 +12,7 @@
 
 import qualified Data.ByteString        as B
 import qualified Data.ByteString.Char8  as T
-import qualified Data.Map               as M
+import qualified Data.Map.Strict        as M
 import qualified Data.Vector.Generic    as V
 import qualified Data.Vector.Storable   as VS
 import qualified Data.Vector.Unboxed    as U
@@ -115,10 +115,10 @@
     same_pos u v = b_cpos u == b_cpos v
     b_cpos u = (b_rname u, b_pos u)
 
-    nice_sort x = sortBy (comparing (V.length . b_seq . snd)) x
+    nice_sort = sortBy $ comparing (V.length . b_seq . snd)
 
     mapGroups f o = tryHead >>= maybe (return o) (\a -> eneeCheckIfDone (mg1 f a []) o)
-    mg1 f a acc k = tryHead >>= \mb -> case mb of
+    mg1 f a acc k = tryHead >>= \case
                         Nothing -> return . k . Chunk . f $ a:acc
                         Just b | same_pos a b -> mg1 f a (b:acc) k
                                | otherwise -> eneeCheckIfDone (mg1 f b []) . k . Chunk . f $ a:acc
@@ -209,7 +209,7 @@
         (half_unaligned, half_aligned) = partition isUnmapped raw_half_pairs
 
         mkMap :: Ord a => (BamRec -> a) -> [BamRec] -> (M.Map a (Int,Decision), [BamRec])
-        mkMap f x = let m1 = M.map (\xs -> (length xs, collapse xs)) $ accumMap f id x
+        mkMap f x = let m1 = M.map (length &&& collapse) $ accumMap f id x
                     in (M.map (second fst) m1, concatMap (snd.snd) $ M.elems m1)
 
         (pairs',r1)   = mkMap (\b -> (b_mate_pos b,   b_strand b, b_mate b)) pairs
@@ -257,7 +257,7 @@
 
         in case M.lookup str singles of
             Nothing                    -> (              (m,v) : r,     unal )
-            Just (n, Consensus      w) -> ( (n, add_xp_of w v) : r,     unal )      -- XXX do we need this w?!
+            Just (n, Consensus      w) -> ( (n, add_xp_of w v) : r,     unal )
             Just (n, Representative w) -> ( (n, add_xp_of w v) : r, w : unal )
 
     -- No more pairs, delegate the problem
@@ -412,11 +412,11 @@
     xa' = nub' [ T.split ';' xas | Just (Text xas) <- map (lookup "XA" . b_exts) brs ]
 
     modify_extensions es = foldr ($!) es $
-        [ let vs = [ v | Just v <- map (lookup k . b_exts) brs ]
+        [ let vs = mapMaybe (lookup k . b_exts) brs
           in if null vs then id else updateE k $! maj vs | k <- do_maj ] ++
         [ let vs = [ v | Just (Int v) <- map (lookup k . b_exts) brs ]
           in if null vs then id else updateE k $! Int (rmsq vs) | k <- do_rmsq ] ++
-        [ deleteE k | k <- useless ] ++
+        map deleteE useless ++
         [ updateE "NM" $! Int nm'
         , updateE "XP" $! Int (foldl' (\a b -> a `oplus` extAsInt 1 "XP" b) 0 brs)
         , if null xa' then id else updateE "XA" $! (Text $ T.intercalate (T.singleton ';') xa')
@@ -528,7 +528,7 @@
   where
     overhangs = not (isUnmapped b) && b_pos b < l && l < b_pos b + alignedLength (b_cigar b)
 
-    do_wrap = case split_ecig (l - b_pos b) $ toECig (b_cigar b) (maybe [] id $ getMd b) of
+    do_wrap = case split_ecig (l - b_pos b) $ toECig (b_cigar b) (fromMaybe [] $ getMd b) of
                   (left,right) -> [ Right $ b { b_cigar = toCigar  left }            `setMD` left
                                   , Left  $ b { b_cigar = toCigar right, b_pos = 0 } `setMD` right ]
 
@@ -569,12 +569,10 @@
 mask_all (Ins' n ec) = SMa' n $ mask_all ec
 mask_all (SMa' n ec) = SMa' n $ mask_all ec
 
--- | Argh, this business with the CIGAR operations is a mess, it gets
--- worse when combined with MD.  Okay, we will support CIGAR (no "=" and
--- "X" operations) and MD.  If we have MD on input, we generate it on
--- output, too.  And in between, we break everything into /very small/
--- operations.  (Yes, the two terminating constructors are a weird
--- hack.)
+-- | Extended CIGAR.  This subsumes both the CIGAR string and the
+-- optional MD field.  If we have MD on input, we generate it on output,
+-- too.  And in between, we break everything into /very small/
+-- operations.
 
 data ECig = WithMD                      -- terminate, do generate MD field
           | WithoutMD                   -- terminate, don't bother with MD
@@ -589,7 +587,7 @@
 
 
 toECig :: VS.Vector Cigar -> [MdOp] -> ECig
-toECig cig md = go (VS.toList cig) md
+toECig = go . VS.toList
   where
     go        cs  (MdNum  0:mds) = go cs mds
     go        cs  (MdDel []:mds) = go cs mds
diff --git a/src/Bio/Bam/Trim.hs b/src/Bio/Bam/Trim.hs
--- a/src/Bio/Bam/Trim.hs
+++ b/src/Bio/Bam/Trim.hs
@@ -30,7 +30,7 @@
 import qualified Data.Vector.Storable                   as W
 
 -- | Trims from the 3' end of a sequence.
--- @trim_3\' p b@ trims the 3' end of the sequence in @b@ at the
+-- @trim_3' p b@ trims the 3' end of the sequence in @b@ at the
 -- earliest position such that @p@ evaluates to true on every suffix
 -- that was trimmed off.  Note that the 3' end may be the beginning of
 -- the sequence if it happens to be stored in reverse-complemented form.
@@ -229,7 +229,7 @@
 merged_seq :: (V.Vector v Nucleotides, V.Vector v Qual)
            => Int -> v Nucleotides -> v Qual -> v Nucleotides -> v Qual -> v Nucleotides
 merged_seq l b_seq_r1 b_qual_r1 b_seq_r2 b_qual_r2 = V.concat
-        [ V.take (l - len_r2) (b_seq_r1)
+        [ V.take (l - len_r2) b_seq_r1
         , V.zipWith4 zz          (V.take l $ V.drop (l - len_r2) b_seq_r1)
                                  (V.take l $ V.drop (l - len_r2) b_qual_r1)
                      (V.reverse $ V.take l $ V.drop (l - len_r1) b_seq_r2)
@@ -246,7 +246,7 @@
 merged_qual :: (V.Vector v Nucleotides, V.Vector v Qual)
             => Word8 -> Int -> v Nucleotides -> v Qual -> v Nucleotides -> v Qual -> v Qual
 merged_qual qmax l b_seq_r1 b_qual_r1 b_seq_r2 b_qual_r2 = V.concat
-        [ V.take (l - len_r2) (b_qual_r1)
+        [ V.take (l - len_r2) b_qual_r1
         , V.zipWith4 zz           (V.take l $ V.drop (l - len_r2) b_seq_r1)
                                   (V.take l $ V.drop (l - len_r2) b_qual_r1)
                       (V.reverse $ V.take l $ V.drop (l - len_r1) b_seq_r2)
diff --git a/src/Bio/Bam/Writer.hs b/src/Bio/Bam/Writer.hs
--- a/src/Bio/Bam/Writer.hs
+++ b/src/Bio/Bam/Writer.hs
@@ -1,3 +1,6 @@
+-- | Printers for BAM and SAM.  BAM is properly supported, SAM can be
+-- piped to standard output.
+
 module Bio.Bam.Writer (
     IsBamRec(..),
     encodeBamWith,
@@ -29,12 +32,6 @@
 import qualified Data.Vector.Unboxed                as U
 import qualified Data.Sequence                      as Z
 
--- ^ Printers for BAM.  We employ an @Iteratee@ interface, and we strive
--- to keep BAM records in their encoded form.  This is most compact and
--- often faster, since it saves the time for repeated decoding and
--- encoding, if that's not strictly needed.
-
-
 -- | write in SAM format to stdout
 -- This is useful for piping to other tools (say, AWK scripts) or for
 -- debugging.  No convenience function to send SAM to a file exists,
@@ -118,11 +115,7 @@
 pushBamRaw :: BamRaw -> BgzfTokens -> BgzfTokens
 pushBamRaw = TkLnString . raw_data
 
--- | writes BAM encoded stuff to a file
--- XXX This should(!) write indexes on the side---a simple block index
--- for MapReduce style slicing, a standard BAM index or a name index
--- would be possible.  When writing to a file, this makes even more
--- sense than when writing to a @Handle@.
+-- | Writes BAM encoded stuff to a file.
 writeBamFile :: IsBamRec r => FilePath -> BamMeta -> Iteratee [r] IO ()
 writeBamFile fp meta =
     C.bracket (liftIO $ openBinaryFile fp WriteMode)
@@ -135,12 +128,7 @@
 pipeBamOutput :: IsBamRec r => BamMeta -> Iteratee [r] IO ()
 pipeBamOutput meta = encodeBamWith 0 meta =$ mapChunksM_ (liftIO . S.hPut stdout)
 
--- | writes BAM encoded stuff to a @Handle@
--- We generate BAM with dynamic blocks, then stream them out to the file.
---
--- XXX This could write indexes on the side---a simple block index
--- for MapReduce style slicing, a standard BAM index or a name index
--- would be possible.
+-- | Writes BAM encoded stuff to a 'Handle'.
 writeBamHandle :: (MonadIO m, IsBamRec r) => Handle -> BamMeta -> Iteratee [r] m ()
 writeBamHandle hdl meta = encodeBamWith 6 meta =$ mapChunksM_ (liftIO . S.hPut hdl)
 
diff --git a/src/Bio/Base.hs b/src/Bio/Base.hs
--- a/src/Bio/Base.hs
+++ b/src/Bio/Base.hs
@@ -52,11 +52,6 @@
 import qualified Data.ByteString.Char8 as S
 import qualified Data.Vector.Unboxed   as U
 
-#if __GLASGOW_HASKELL__ == 704
-import Data.Vector.Generic          ( Vector(..) )
-import Data.Vector.Generic.Mutable  ( MVector(..) )
-#endif
-
 -- | A nucleotide base.  We only represent A,C,G,T.  The contained
 -- 'Word8' ist guaranteed to be 0..3.
 newtype Nucleotide = N { unN :: Word8 } deriving ( Eq, Ord, Enum, Ix, Storable )
@@ -299,7 +294,7 @@
 
     readList ('-':cs) = readList cs
     readList (c:cs) | isSpace c = readList cs
-                    | otherwise = case readsPrec 0 (c:cs) of
+                    | otherwise = case reads (c:cs) of
                             [] -> [ ([],c:cs) ]
                             xs -> [ (n:ns,r2) | (n,r1) <- xs, (ns,r2) <- readList r1 ]
     readList [] = [([],[])]
diff --git a/src/Bio/Iteratee.hs b/src/Bio/Iteratee.hs
--- a/src/Bio/Iteratee.hs
+++ b/src/Bio/Iteratee.hs
@@ -1,6 +1,3 @@
--- | Basically a reexport of "Data.Iteratee" less the names that clash
--- with "Prelude" plus a handful of utilities.
-
 module Bio.Iteratee (
     iGetString,
     iterLoop,
@@ -102,19 +99,16 @@
                         else it a >>= iterLoop it
 infixl 1 $==
 {-# INLINE ($==) #-}
--- | Compose an 'Enumerator\'' with an 'Enumeratee', giving a new
--- 'Enumerator\''.
+-- | Compose an 'Enumerator'' with an 'Enumeratee', giving a new
+-- 'Enumerator''.
 ($==) :: Monad m => Enumerator' hdr input m (Iteratee output m result)
                  -> Enumeratee      input             output m result
                  -> Enumerator' hdr                   output m result
 ($==) enum enee iter = run =<< enum (enee . iter)
 
--- | Merge two 'Enumerator\''s into one.  The header provided by the
--- inner 'Enumerator\'' is passed to the output iterator, the header
--- provided by the outer 'Enumerator\'' is passed to the merging iteratee
---
--- XXX  Something about those headers is unsatisfactory... there should
---      be an unobtrusive way to combine headers.
+-- | Merge two 'Enumerator''s into one.  The header provided by the
+-- inner 'Enumerator'' is passed to the output iterator, the header
+-- provided by the outer 'Enumerator'' is passed to the merging iteratee
 
 {-# INLINE mergeEnums' #-}
 mergeEnums' :: (Nullable s2, Nullable s1, Monad m)
@@ -259,9 +253,9 @@
 pushQ a (QQ l  f b) = QQ (l+1) f (a:b)
 
 popQ :: QQ a -> Maybe (a, QQ a)
-popQ (QQ l (a:[]) b) = Just (a, QQ (l-1) (reverse b) [])
-popQ (QQ l (a:fs) b) = Just (a, QQ (l-1) fs b)
 popQ (QQ _ [    ] _) = Nothing
+popQ (QQ l [ a  ] b) = Just (a, QQ (l-1) (reverse b) [])
+popQ (QQ l (a:fs) b) = Just (a, QQ (l-1) fs b)
 
 cancelAll :: MonadIO m => QQ (Async a) -> m ()
 cancelAll (QQ _ ff bb) = liftIO $ mapM_ cancel (ff ++ bb)
@@ -302,7 +296,7 @@
     go mv i
         | i == n    = return n
         | otherwise =
-            tryHead >>= \x -> case x of
+            tryHead >>= \case
                 Nothing -> return i
                 Just  a -> liftIO (VM.write mv i a) >> go mv (i+1)
 
@@ -310,7 +304,7 @@
 stream2vector :: (MonadIO m, VG.Vector v a) => Iteratee [a] m (v a)
 stream2vector = liftIO (VM.new 1024) >>= go 0
   where
-    go !i !mv = tryHead >>= \x -> case x of
+    go !i !mv = tryHead >>= \case
                   Nothing -> liftIO $ VG.unsafeFreeze $ VM.take i mv
                   Just  a -> do mv' <- if VM.length mv == i then liftIO (VM.grow mv (VM.length mv)) else return mv
                                 when (i `rem` 0x10000 == 0) $ liftIO performGC
@@ -318,7 +312,7 @@
                                 go (i+1) mv'
 
 withFileFd :: (MonadIO m, MonadMask m) => FilePath -> (Fd -> m a) -> m a
-withFileFd filepath iter = CMC.bracket
+withFileFd filepath = CMC.bracket
     (liftIO $ openFd filepath ReadOnly Nothing defaultFileFlags)
-    (liftIO . closeFd) iter
+    (liftIO . closeFd)
 
diff --git a/src/Bio/Iteratee/Base.hs b/src/Bio/Iteratee/Base.hs
--- a/src/Bio/Iteratee/Base.hs
+++ b/src/Bio/Iteratee/Base.hs
@@ -15,7 +15,6 @@
   -- ** Control functions
   ,run
   ,tryRun
-  ,mapIteratee
   ,ilift
   ,ifold
   -- ** Creating Iteratees
@@ -168,9 +167,6 @@
 instance NullPoint s => MonadTrans (Iteratee s) where
   lift m = Iteratee $ \onDone _ -> m >>= flip onDone (Chunk emptyP)
 
--- instance (MonadBase b m, Nullable s, NullPoint s) => MonadBase b (Iteratee s m) where
-  -- liftBase = lift . liftBase
-
 instance (MonadIO m, Nullable s, NullPoint s) => MonadIO (Iteratee s m) where
   liftIO = lift . liftIO
 
@@ -216,14 +212,6 @@
     onCont' _ Nothing  = return $ maybeExc (toException EofException)
     onCont' _ (Just e) = return $ maybeExc e
     maybeExc e = maybe (Left (E.throw e)) Left (fromException e)
-
--- |Transform a computation inside an @Iteratee@.
-mapIteratee :: (NullPoint s, Monad n, Monad m) =>
-  (m a -> n b)
-  -> Iteratee s m a
-  -> Iteratee s n b
-mapIteratee f = lift . f . run
-{-# DEPRECATED mapIteratee "This function will be removed, compare to 'ilift'" #-}
 
 -- | Lift a computation in the inner monad of an iteratee.
 --
diff --git a/src/Bio/Iteratee/Bgzf.hsc b/src/Bio/Iteratee/Bgzf.hsc
--- a/src/Bio/Iteratee/Bgzf.hsc
+++ b/src/Bio/Iteratee/Bgzf.hsc
@@ -17,8 +17,6 @@
 import Foreign.C.String                     ( withCAString )
 import Foreign.C.Types                      ( CInt(..), CChar(..), CUInt(..), CULong(..) )
 import Foreign.Marshal.Alloc                ( mallocBytes, free, allocaBytes )
-import Foreign.Ptr                          ( nullPtr, castPtr, Ptr, plusPtr, minusPtr )
-import Foreign.Storable                     ( peekByteOff, pokeByteOff )
 
 import qualified Data.ByteString            as S
 import qualified Data.ByteString.Unsafe     as S
diff --git a/src/Bio/Iteratee/Builder.hs b/src/Bio/Iteratee/Builder.hs
--- a/src/Bio/Iteratee/Builder.hs
+++ b/src/Bio/Iteratee/Builder.hs
@@ -1,12 +1,8 @@
 {-# LANGUAGE ForeignFunctionInterface #-}
--- | Buffer builder to assemble Bgzf blocks.  The plan is to serialize
+-- | Buffer builder to assemble Bgzf blocks.  The idea is to serialize
 -- stuff (BAM and BCF) into a buffer, then Bgzf chunks from the buffer.
 -- We use a large buffer, and we always make sure there is plenty of
--- space in it (to avoid redundant checks).  Whenever a block is ready
--- to be compressed, we stick it into a MVar.  When we run out of space,
--- we simply use a new buffer.  Multiple threads grab pieces from the
--- MVar, compress them, pass them downstream through another MVar.  A
--- final thread restores the order and writes the blocks.
+-- space in it (to avoid redundant checks).
 
 module Bio.Iteratee.Builder (
     BB(..),
@@ -94,12 +90,12 @@
     withForeignPtr arr1 $ \d ->
         withForeignPtr (buffer b) $ \s ->
              copyBytes d (plusPtr s (off b)) (used b - off b)
-    return $ BB { buffer = arr1
-                , size   = sz'
-                , off    = 0
-                , used   = used b - off b
-                , mark   = if mark  b == maxBound then maxBound else mark  b - off b
-                , mark2  = if mark2 b == maxBound then maxBound else mark2 b - off b }
+    return BB{ buffer = arr1
+             , size   = sz'
+             , off    = 0
+             , used   = used b - off b
+             , mark   = if mark  b == maxBound then maxBound else mark  b - off b
+             , mark2  = if mark2 b == maxBound then maxBound else mark2 b - off b }
 
 compressChunk' :: Int -> ForeignPtr Word8 -> Int -> Int -> IO B.ByteString
 compressChunk' lv fptr off len =
@@ -134,7 +130,7 @@
             case tk of TkEnd -> liftI $ go bb k
                        _     -> go1 bb k tk
 
-        | otherwise = do
+        | otherwise =
             eneeCheckIfDone (flush_blocks tk bb { off = off bb + maxBlockSize }) $
                 k $ Chunk [compressChunk' lv (buffer bb) (off bb) maxBlockSize]
 
@@ -149,7 +145,7 @@
 fillBuffer :: BB -> BgzfTokens -> IO (BB, BgzfTokens)
 fillBuffer bb0 tk = withForeignPtr (buffer bb0) (\p -> go_slowish p bb0 tk)
   where
-    go_slowish p bb tk1 = go_fast p bb (used bb) tk1
+    go_slowish p bb = go_fast p bb (used bb)
 
     go_fast p bb use tk1 = case tk1 of
         -- no space?  not our job.
@@ -178,8 +174,7 @@
 
         TkString   s tk'
             -- Too big, can't handle.  We will get progressively bigger
-            -- buffers and eventually handle it; for very large strings,
-            -- it works, but isn't ideal.  XXX
+            -- buffers and eventually handle it.
             | B.length s > size bb - use -> return (bb { used = use },tk')
 
             | otherwise  -> do let ln = B.length s
@@ -192,8 +187,7 @@
 
         TkLnString s tk'
             -- Too big, can't handle.  We will get progressively bigger
-            -- buffers and eventually handle it; for very large strings,
-            -- it works, but isn't ideal.  XXX
+            -- buffers and eventually handle it.
             | B.length s > size bb - use - 4 -> return (bb { used = use },tk')
 
             | otherwise  -> do let ln = B.length s
diff --git a/src/Bio/Iteratee/Bytes.hs b/src/Bio/Iteratee/Bytes.hs
--- a/src/Bio/Iteratee/Bytes.hs
+++ b/src/Bio/Iteratee/Bytes.hs
@@ -226,7 +226,7 @@
 -- This is provided as a higher-performance alternative to enumWords, and
 -- is equivalent to treating the stream as a Data.ByteString.Char8.ByteString.
 enumWordsBS :: Monad m => Enumeratee Bytes [Bytes] m a
-enumWordsBS iter = convStream getter iter
+enumWordsBS = convStream getter
   where
     getter = liftI step
     lChar = isSpace . C.last
diff --git a/src/Bio/Iteratee/Exception.hs b/src/Bio/Iteratee/Exception.hs
--- a/src/Bio/Iteratee/Exception.hs
+++ b/src/Bio/Iteratee/Exception.hs
@@ -118,7 +118,7 @@
   fromException = enumExceptionFromException
 
 -- |Create an enumerator exception from a @String@.
-data EnumStringException = EnumStringException String
+newtype EnumStringException = EnumStringException String
   deriving (Show, Typeable)
 
 instance Exception EnumStringException where
@@ -130,7 +130,7 @@
 enStrExc = EnumException . EnumStringException
 
 -- |The enumerator received an 'IterException' it could not handle.
-data EnumUnhandledIterException = EnumUnhandledIterException IterException
+newtype EnumUnhandledIterException = EnumUnhandledIterException IterException
   deriving (Show, Typeable)
 
 instance Exception EnumUnhandledIterException where
@@ -177,7 +177,7 @@
   fromIterException = Just
 
 -- |A seek request within an @Iteratee@.
-data SeekException = SeekException FileOffset
+newtype SeekException = SeekException FileOffset
   deriving (Typeable, Show)
 
 instance Exception SeekException where
@@ -197,7 +197,7 @@
 instance IException EofException where
 
 -- |An @Iteratee exception@ specified by a @String@.
-data IterStringException = IterStringException String deriving (Typeable, Show)
+newtype IterStringException = IterStringException String deriving (Typeable, Show)
 
 instance Exception IterStringException where
   toException   = iterExceptionToException
diff --git a/src/Bio/Iteratee/IO.hs b/src/Bio/Iteratee/IO.hs
--- a/src/Bio/Iteratee/IO.hs
+++ b/src/Bio/Iteratee/IO.hs
@@ -56,10 +56,10 @@
 -- |The enumerator of a POSIX File Descriptor: a variation of @enumFd@ that
 -- supports RandomIO (seek requests).
 enumFdRandom :: MonadIO m => Int -> Fd -> Enumerator Bytes m a
-enumFdRandom bs fd iter = enumFdCatch bs fd handler iter
+enumFdRandom bs fd = enumFdCatch bs fd handler
   where
     handler (SeekException off) =
-        Nothing <$ (liftIO . fdSeek fd AbsoluteSeek $ fromIntegral off)
+        liftIO $ fdSeek fd AbsoluteSeek (fromIntegral off) >> return Nothing
 
 enumFile' :: (MonadIO m, MonadMask m) =>
   (Int -> Fd -> Enumerator s m a)
diff --git a/src/Bio/Iteratee/Iteratee.hs b/src/Bio/Iteratee/Iteratee.hs
--- a/src/Bio/Iteratee/Iteratee.hs
+++ b/src/Bio/Iteratee/Iteratee.hs
@@ -2,7 +2,6 @@
             ,RankNTypes
             ,FlexibleContexts
             ,ScopedTypeVariables
-            ,BangPatterns
             ,DeriveDataTypeable #-}
 
 -- |Monadic and General Iteratees:
@@ -516,7 +515,7 @@
 -- | Enumerate chunks from a list
 --
 enumList :: (Monad m) => [s] -> Enumerator s m a
-enumList chunks = go chunks
+enumList = go
  where
   go [] i = return i
   go xs' i = runIter i idoneM (onCont xs')
@@ -544,8 +543,8 @@
   (st -> m (Either SomeException ((Bool, st), s)))
   -> st
   -> Enumerator s m a
-enumFromCallback c st =
-  enumFromCallbackCatch c (\NotAnException -> return Nothing) st
+enumFromCallback c =
+  enumFromCallbackCatch c (\NotAnException -> return Nothing)
 
 -- Dummy exception to catch in enumFromCallback
 -- This never gets thrown, but it lets us
diff --git a/src/Bio/Iteratee/List.hs b/src/Bio/Iteratee/List.hs
--- a/src/Bio/Iteratee/List.hs
+++ b/src/Bio/Iteratee/List.hs
@@ -770,7 +770,7 @@
             onDone a str'@(Chunk c') =
                 od (a, newLen - fromIntegral (length c')) str'
             onDone a str'@(EOF _) = od (a, n) str'
-            onCont f' mExc = oc (go newLen f') mExc
+            onCont f' = oc (go newLen f')
 {-# INLINE countConsumed #-}
 
 -- ------------------------------------------------------------------------
@@ -808,8 +808,5 @@
 
 -- | Folds a monadic function over an 'Iteratee'.
 foldStreamM :: Monad m => (b -> a -> m b) -> b -> Iteratee [a] m b
-foldStreamM k = foldChunksM go
-  where
-    go b [   ] = return b
-    go b (h:t) = k b h >>= \b' -> go b' t
+foldStreamM = foldChunksM . foldM
 {-# INLINE foldStreamM #-}
diff --git a/src/Bio/Prelude.hs b/src/Bio/Prelude.hs
--- a/src/Bio/Prelude.hs
+++ b/src/Bio/Prelude.hs
@@ -6,12 +6,6 @@
     module System.Posix.Files,
     module System.Posix.IO,
     module System.Posix.Types,
-#if !MIN_VERSION_base_prelude(1,1,0)
-    module Foreign.Storable,
-    module Foreign.Ptr,
-    module Foreign.ForeignPtr,
-    module Foreign.StablePtr,
-#endif
 
     Bytes, LazyBytes,
     HashMap,
@@ -26,11 +20,6 @@
 #endif
 #endif
 
-#if !MIN_VERSION_base(4,7,0)
-    isLeft,
-    isRight,
-#endif
-
 #if !MIN_VERSION_base(4,8,0)
     first,
     second,
@@ -49,10 +38,8 @@
 import BasePrelude
 #if MIN_VERSION_base(4,9,0)
                     hiding ( EOF, log1p, log1pexp, log1mexp, expm1 )
-#elif MIN_VERSION_base(4,7,0)
-                    hiding ( EOF )
 #else
-                    hiding ( EOF, block )
+                    hiding ( EOF )
 #endif
 
 #if !MIN_VERSION_base(4,8,0)
@@ -75,13 +62,6 @@
 import System.Posix.IO
 import System.Posix.Types
 
-#if !MIN_VERSION_base_prelude(1,1,0)
-import Foreign.Storable
-import Foreign.Ptr
-import Foreign.ForeignPtr
-import Foreign.StablePtr
-#endif
-
 import qualified Data.ByteString.Unsafe as B
 import qualified Data.ByteString.Lazy   as BL
 import qualified Data.ByteString.Char8  as S
@@ -112,12 +92,6 @@
 instance Unpack Text       where unpack = T.unpack
 instance Unpack String     where unpack = id
 
-#if !MIN_VERSION_base(4,7,0)
-isLeft, isRight :: Either a b -> Bool
-isLeft = either (const False) (const True)
-isRight = either (const True) (const False)
-#endif
-
 fdPut :: Fd -> Bytes -> IO ()
 fdPut fd s = B.unsafeUseAsCStringLen s $ \(p,l) ->
              throwErrnoIf_ (/= fromIntegral l) "fdPut" $
@@ -128,7 +102,7 @@
 
 withFd :: FilePath -> OpenMode -> Maybe FileMode -> OpenFileFlags
        -> (Fd -> IO a) -> IO a
-withFd fp om fm ff k = bracket (openFd fp om fm ff) closeFd k
+withFd fp om fm ff = bracket (openFd fp om fm ff) closeFd
 
 -- | Converts 'Bytes' into 'Text'.  This uses UTF8, but if there is an
 -- error, it pretends it was Latin1.  Evil as this is, it tends to Just
diff --git a/src/Bio/TwoBit.hs b/src/Bio/TwoBit.hs
--- a/src/Bio/TwoBit.hs
+++ b/src/Bio/TwoBit.hs
@@ -1,3 +1,14 @@
+-- | Would you believe it?  The 2bit format stores blocks of Ns in a table at
+-- the beginning of a sequence, then packs four bases into a byte.  So it
+-- is neither possible nor necessary to store Ns in the main sequence, and
+-- you would think they aren't stored there, right?  And they aren't.
+-- Instead Ts are stored which the reader has to replace with Ns.
+--
+-- The sensible way to treat these is probably to just say there are two
+-- kinds of implied annotation (repeats and large gaps for a typical
+-- genome), which can be interpreted in whatever way fits.  And that's why
+-- we have 'Mask' and 'getSubseqWith'.
+
 module Bio.TwoBit (
         TwoBitFile(..),
         TwoBitSequence(..),
@@ -27,24 +38,11 @@
 import           Data.Binary.Get
 import qualified Data.ByteString                as B
 import qualified Data.ByteString.Lazy           as L
-import qualified Data.IntMap                    as I
+import qualified Data.IntMap.Strict             as I
 import qualified Data.HashMap.Lazy              as M
 import qualified Data.Vector.Unboxed            as U
 import           System.Random
 
--- ^ Would you believe it?  The 2bit format stores blocks of Ns in a table at
--- the beginning of a sequence, then packs four bases into a byte.  So it
--- is neither possible nor necessary to store Ns in the main sequence, and
--- you would think they aren't stored there, right?  And they aren't.
--- Instead Ts are stored which the reader has to replace with Ns.
---
--- The sensible way to treat these is probably to just say there are two
--- kinds of implied annotation (repeats and large gaps for a typical
--- genome), which can be interpreted in whatever way fits.  And that's why
--- we have 'Mask' and 'getSubseqWith'.
---
--- TODO:  use binary search for the Int->Int mappings on the raw data?
-
 data TwoBitFile = TBF {
     tbf_raw :: B.ByteString,
     -- This map is intentionally lazy.  May or may not be important.
@@ -91,7 +89,7 @@
 
     readBlockList = getWord32 >>= \n -> liftM2 zip (repM n getWord32) (repM n getWord32)
 
--- | Repeat monadic action 'n' times.  Returns result in reverse(!)
+-- | Repeat monadic action @n@ times.  Returns result in reverse(!)
 -- order, but doesn't build a huge list of thunks in memory.
 repM :: Monad m => Int -> m a -> m [a]
 repM n0 m = go [] n0
@@ -150,8 +148,8 @@
 -- lowercasing.  Here, we take a user supplied function to apply
 -- masking.
 getSubseqWith :: (Nucleotide -> Mask -> a) -> TwoBitFile -> Range -> [a]
-getSubseqWith maskf tbf (Range { r_pos = Pos { p_seq = chr, p_start = start }, r_length = len }) = do
-    let sq1 = maybe (error $ unpack chr ++ " doesn't exist") id $ M.lookup chr (tbf_seqs tbf)
+getSubseqWith maskf tbf Range{ r_pos = Pos { p_seq = chr, p_start = start }, r_length = len } = do
+    let sq1 = fromMaybe (error $ unpack chr ++ " doesn't exist") $ M.lookup chr (tbf_seqs tbf)
     let go = getFwdSubseqWith tbf sq1
     if start < 0
         then reverse $ take len $ go (maskf . cmp_nt) (-start-len)
@@ -162,8 +160,8 @@
 
 -- | Works only in forward direction.
 getLazySubseq :: TwoBitFile -> Position -> [Nucleotide]
-getLazySubseq tbf (Pos { p_seq = chr, p_start = start }) = do
-    let sq1 = maybe (error $ unpack chr ++ " doesn't exist") id $ M.lookup chr (tbf_seqs tbf)
+getLazySubseq tbf Pos{ p_seq = chr, p_start = start } = do
+    let sq1 = fromMaybe (error $ unpack chr ++ " doesn't exist") $ M.lookup chr (tbf_seqs tbf)
     let go  = getFwdSubseqWith tbf sq1
     if start < 0
         then error "sorry, can't go backwards"
diff --git a/src/Bio/Util/MMap.hs b/src/Bio/Util/MMap.hs
--- a/src/Bio/Util/MMap.hs
+++ b/src/Bio/Util/MMap.hs
@@ -1,7 +1,7 @@
 {-# LANGUAGE ForeignFunctionInterface #-}
 module Bio.Util.MMap where
 
-import Bio.Prelude hiding ( left, right, chr )
+import Bio.Prelude
 import Data.ByteString.Internal ( fromForeignPtr )
 import Foreign.C.Types
 
diff --git a/src/Bio/Util/Numeric.hs b/src/Bio/Util/Numeric.hs
--- a/src/Bio/Util/Numeric.hs
+++ b/src/Bio/Util/Numeric.hs
@@ -1,19 +1,18 @@
+-- | Random useful stuff I didn't know where to put.
+
 module Bio.Util.Numeric (
     wilson, invnormcdf, choose,
     estimateComplexity, showNum, showOOM,
     log1p, expm1, (<#>),
     log1mexp, log1pexp,
-    lsum, llerp,
-    sigmoid2, isigmoid2
+    lsum, llerp
                 ) where
 
 import Data.List ( foldl1' )
 import Data.Char ( intToDigit )
 import Prelude
 
--- ^ Random useful stuff I didn't know where to put.
-
--- | calculates the Wilson Score interval.
+-- | Calculates the Wilson Score interval.
 -- If @(l,m,h) = wilson c x n@, then @m@ is the binary proportion and
 -- @(l,h)@ it's @c@-confidence interval for @x@ positive examples out of
 -- @n@ observations.  @c@ is typically something like 0.05.
@@ -33,9 +32,9 @@
 showNum = triplets [] . reverse . show
   where
     triplets acc [] = acc
-    triplets acc (a:[]) = a:acc
-    triplets acc (a:b:[]) = b:a:acc
-    triplets acc (a:b:c:[]) = c:b:a:acc
+    triplets acc [a] = a:acc
+    triplets acc [a,b] = b:a:acc
+    triplets acc [a,b,c] = c:b:a:acc
     triplets acc (a:b:c:s) = triplets (',':c:b:a:acc) s
 
 showOOM :: Double -> String
@@ -51,8 +50,8 @@
                                             '0' : if s == '.' then [] else [s]
                         | otherwise = findSuffix (y*0.001) ss
 
--- Stolen from Lennart Augustsson's erf package, who in turn took it rom
--- http://home.online.no/~pjacklam/notes/invnorm/ Accurate to about 1e-9.
+-- Stolen from Lennart Augustsson's erf package, who in turn took it from
+-- <http://home.online.no/~pjacklam/notes/invnorm/> Accurate to about 1e-9.
 invnormcdf :: (Ord a, Floating a) => a -> a
 invnormcdf p =
     let a1 = -3.969683028665376e+01
@@ -108,25 +107,31 @@
 -- How many different things are there?
 --
 -- Let the total number be @m@.  The copy number follows a Poisson
--- distribution with paramter @\lambda@.  Let @z := e^{\lambda}@, then
+-- distribution with paramter @\lambda@.  Let \( z := e^{\lambda} \), then
 -- we have:
 --
---   P( 0 ) = e^{-\lambda} = 1/z
---   P( 1 ) = \lambda e^{-\lambda} = ln z / z
---   P(>=1) = 1 - e^{-\lambda} = 1 - 1/z
---
---   singles = m ln z / z
---   total   = m (1 - 1/z)
---
---   D := total/singles = (1 - 1/z) * z / ln z
---   f := z - 1 - D ln z = 0
+-- \[
+--   P( 0 ) = e^{-\lambda} = \frac{1}{z}                    \\
+--   P( 1 ) = \lambda e^{-\lambda} = \frac{\ln z}{z}                   \\
+--   P(\ge 1) = 1 - e^{-\lambda} = 1 - \frac{1}{z}                    \\
+-- \]
+-- \[
+--   \mbox{singles} = m \frac{\ln z}{z}                   \\
+--   \mbox{total}   = m \left( 1 - \frac{1}{z} \right)                  \\
+-- \]
+-- \[
+--   D := \frac{\mbox{total}}{\mbox{singles}} = (1 - \frac{1}{z}) * \frac{z}{\ln z}                  \\
+--   f := z - 1 - D \ln z = 0
+-- \]
 --
 -- To get @z@, we solve using Newton iteration and then substitute to
 -- get @m@:
 --
---   df/dz = 1 - D/z
---   z' := z - z (z - 1 - D ln z) / (z - D)
---   m = singles * z /log z
+-- \[
+--   df/dz = 1 - D/z                                    \\
+--   z' = z - \frac{ z (z - 1 - D \ln z) }{ z - D }     \\
+--   m = \mbox{singles} * \frac{z}{\ln z}
+-- \]
 --
 -- It converges as long as the initial @z@ is large enough, and @10D@
 -- (in the line for @zz@ below) appears to work well.
@@ -144,7 +149,7 @@
     m = fromIntegral singles * zz / log zz
 
 
--- | Computes @log (exp x + exp y)@ without leaving the log domain and
+-- | Computes \( \ln \left( e^x + e^y \right) \) without leaving the log domain and
 -- hence without losing precision.
 infixl 5 <#>
 {-# INLINE (<#>) #-}
@@ -152,7 +157,7 @@
 x <#> y = if x >= y then x + log1pexp (y-x) else y + log1pexp (x-y)
 
 -- | Computes @log (1+x)@ to a relative precision of @10^-8@ even for
--- very small @x@.  Stolen from http://www.johndcook.com/cpp_log_one_plus_x.html
+-- very small @x@.  Stolen from <http://www.johndcook.com/cpp_log_one_plus_x.html>
 {-# INLINE log1p #-}
 log1p :: (Floating a, Ord a) => a -> a
 log1p x | x < -1 = error "log1p: argument must be greater than -1"
@@ -163,20 +168,20 @@
         | otherwise = (1 - 0.5*x) * x
 
 
--- | Computes @exp x - 1@ to a relative precision of @10^-10@ even for
--- very small @x@.  Stolen from http://www.johndcook.com/cpp_expm1.html
+-- | Computes \( e^x - 1 \) to a relative precision of @10^-10@ even for
+-- very small @x@.  Stolen from <http://www.johndcook.com/cpp_expm1.html>
 {-# INLINE expm1 #-}
 expm1 :: (Floating a, Ord a) => a -> a
 expm1 x | x > -0.00001 && x < 0.00001 = (1 + 0.5 * x) * x       -- Taylor approx
         | otherwise                   = exp x - 1               -- direct eval
 
--- | Computes @log (1 - exp x)@, following Martin Mächler.
+-- | Computes \( \ln (1 - e^x) \), following Martin Mächler.
 {-# INLINE log1mexp #-}
 log1mexp :: (Floating a, Ord a) => a -> a
 log1mexp x | x > - log 2 = log (- expm1 x)
            | otherwise   = log1p (- exp x)
 
--- | Computes @log (1 + exp x)@, following Martin Mächler.
+-- | Computes \( \ln (1 + e^x) \), following Martin Mächler.
 {-# INLINE log1pexp #-}
 log1pexp :: (Floating a, Ord a) => a -> a
 log1pexp x | x <=  -37 = exp x
@@ -185,14 +190,14 @@
            | otherwise = x
 
 
--- | Computes \( \log ( \sum_i e^{x_i} ) \) sensibly.  The list must be
+-- | Computes \( \ln ( \sum_i e^{x_i} ) \) sensibly.  The list must be
 -- sorted in descending(!) order.
 {-# INLINE lsum #-}
 lsum :: (Floating a, Ord a) => [a] -> a
-lsum xs = foldl1' (\x y -> if x >= y then x + log1pexp (y-x) else err) xs
-    where err = error $ "lsum: argument list must be in descending order"
+lsum = foldl1' (\x y -> if x >= y then x + log1pexp (y-x) else err)
+    where err = error "lsum: argument list must be in descending order"
 
--- | Computes \( \log \left( c e^x + (1-c) e^y \right) \).
+-- | Computes \( \ln \left( c e^x + (1-c) e^y \right) \).
 {-# INLINE llerp #-}
 llerp :: (Floating a, Ord a) => a -> a -> a -> a
 llerp c x y | c <= 0.0  = y
@@ -200,20 +205,9 @@
             | x >= y    = log     c  + x + log1p ( (1-c)/c * exp (y-x) )        -- Hmm.
             | otherwise = log1p (-c) + y + log1p ( c/(1-c) * exp (x-y) )        -- Hmm.
 
--- | Binomial coefficient:  @n `choose` k == n! / ((n-k)! k!)@
+-- | Binomial coefficient: \( \mbox{choose n k} = \frac{n!}{(n-k)! k!} \)
 {-# INLINE choose #-}
 choose :: Integral a => a -> a -> a
-n `choose` k = product [n-k+1 .. n] `div` product [2..k]
-
-
--- | Kind-of sigmoid function that maps the reals to the interval
--- @[0,1)@.  Good to compute a probability without introducing boundary
--- conditions.
-sigmoid2 :: (Fractional a, Floating a) => a -> a
-sigmoid2 x = y*y where y = (exp x - 1) / (exp x + 1)
-
--- | Inverse of 'sigmoid2'.
-isigmoid2 :: (Fractional a, Floating a) => a -> a
-isigmoid2 y = log $ (1 + sqrt y) / (1 - sqrt y)
+choose n k = product [n-k+1 .. n] `div` product [2..k]
 
 
diff --git a/src/Bio/Util/Zlib.hs b/src/Bio/Util/Zlib.hs
--- a/src/Bio/Util/Zlib.hs
+++ b/src/Bio/Util/Zlib.hs
@@ -10,10 +10,9 @@
 
 -- | Decompresses Gzip or Bgzf and passes everything else on.  In
 -- reality, it simply decompresses Gzip, and when done, looks for
--- another Gzip stream.  Trailing garbage is returned as is, therefore,
--- uncompressed files are passed through.  Since there is a small chance
--- to attempt compression of an uncompressed stream, the original data
--- is returned in case of an error.
+-- another Gzip stream.  Since there is a small chance to attempt
+-- decompression of an uncompressed stream, the original data is
+-- returned in case of an error.
 decompressGzip :: L.ByteString -> L.ByteString
 decompressGzip s = case L.uncons s of
     Just (31, s') -> case L.uncons s' of
diff --git a/src/cbits/myers_align.h b/src/cbits/myers_align.h
--- a/src/cbits/myers_align.h
+++ b/src/cbits/myers_align.h
@@ -1,10 +1,10 @@
 #ifndef INCLUDED_MYERS_ALIGN
 #define INCLUDED_MYERS_ALIGN
 
-enum myers_align_mode { 
-	myers_align_globally = 0,
-	myers_align_is_prefix = 1,
-	myers_align_has_prefix = 2 } ;
+enum myers_align_mode {
+    myers_align_globally = 0,
+    myers_align_is_prefix = 1,
+    myers_align_has_prefix = 2 } ;
 
 //! \brief aligns two sequences in O(nd) time
 //! This alignment algorithm following Eugene W. Myers: "An O(ND)