diff --git a/Bio/Bam/Evan.hs b/Bio/Bam/Evan.hs
--- a/Bio/Bam/Evan.hs
+++ b/Bio/Bam/Evan.hs
@@ -71,7 +71,24 @@
 
 
 -- | Removes syntactic warts from old read names or the read names used
--- in FastQ files.
+-- in FastQ files.  Supported conventions:
+--
+-- * A name suffix of @/1@ or @/2@ is turned into the first mate or second
+--   mate flag and the read is flagged as paired.
+--
+-- * Same for name prefixes of @F_@ or @R_@, respectively.
+--
+-- * A name prefix of @M_@ flags the sequence as unpaired and merged
+--
+-- * A name prefix of @T_@ flags the sequence as unpaired and trimmed
+--
+-- * 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.
+
 removeWarts :: BamRec -> BamRec
 removeWarts br = br { b_qname = name, b_flag = flags, b_exts = tags }
   where
diff --git a/Bio/Bam/Fastq.hs b/Bio/Bam/Fastq.hs
--- a/Bio/Bam/Fastq.hs
+++ b/Bio/Bam/Fastq.hs
@@ -1,121 +1,208 @@
--- | Parser for @FastA/FastQ@, 'ByteStream' 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.
+-- | Parser for @FastA/FastQ@, 'ByteStream' style, written such that it
+-- works well with module "Bio.Bam".
+--
+-- Input streams are broken into numbered lines, then into records.
+-- Records can start with empty lines, which are ignored, or random
+-- junk, which is ignored, but results in a warning, followed by a
+-- header indicating either a @FastA@ (begins with @\>@ or @;@) or
+-- @FastQ@ record (begins with @\@@).  More description lines begining
+-- with @;@ are allowed, and silently ignored.  All following lines not
+-- starting with @+@, @\>@, @;@ or @\@@ are sequence lines.  (Only) in a
+-- @FastQ@ record, this is followed by a separator line starting with a
+-- @+@, which is ignored, and exactly as many quality lines as there
+-- were sequence lines.  A missing separator results in a warning and
+-- the record being parsed without quality scores.
+--
+-- In sequence lines, IUPAC-IUB ambiguity codes are converted to
+-- 'Nucleotides', white space is skipped silently.  Any other character
+-- becomes an unknown base ('=' in SAM) and a warning is emitted.  Note
+-- that downstream tools are unlikely to handle the resulting unknown
+-- bases and/or empty records gracefully.  If the quality lines do not
+-- have the same total length as the sequence lines (this includes
+-- missing quality lines due to end-of-stream), a warning is emitted,
+-- and the record receives no quality scores (just as if it was a
+-- @FastA@ record).  Else, if the quality lines have a different layout
+-- than the sequence lines, a warning is emitted, but they are still
+-- used.
+--
+-- Quality scores must be stored as raw bytes with offset 33.  (Other
+-- variants, like 454's ASCII qualities and Solexa's raw bytes with
+-- offset 64 are difficult to detect, and extinct in the wild anyway.)
+-- If the second word of the header stores multiple fields, we try to
+-- extract Illumina's \"QC failed\" flag and either an index sequence or
+-- a read group name from it.
+--
+-- Other flags are commonly encoded into the sequence names.  We do not
+-- handle those here, but most of the conventions at MPI EVAN are dealt
+-- with by 'Bio.Bam.Evan.removeWarts'.
 
-module Bio.Bam.Fastq ( parseFastq, parseFastqWith, parseFastqCassava ) where
+module Bio.Bam.Fastq
+    ( parseFastq
+    , EmptyRecord(..)
+    , IncoherentQualities(..)
+    , IncongruentQualities(..)
+    , JunkFound(..)
+    , QualitiesMissing(..)
+    , SequenceHasGaps(..)
+    ) where
 
 import Bio.Bam.Header
 import Bio.Bam.Rec
-import Bio.Prelude hiding ( isSpace )
+import Bio.Prelude
 import Bio.Streaming
-import Bio.Streaming.Parse ( parse, atto )
-import Data.Attoparsec.ByteString.Char8
-        ( char, skipSpace, satisfy, inClass, skipWhile, takeTill
-        , scan, isSpace, isSpace_w8, (<?>) )
+import Bio.Streaming.Bytes ( lines' )
 
-import qualified Data.Attoparsec.ByteString.Char8   as P
 import qualified Data.ByteString                    as B
 import qualified Data.ByteString.Char8              as C
-import qualified Data.Vector.Generic                as V
+import qualified Data.ByteString.Unsafe             as B
+import qualified Data.Vector.Storable               as W
 import qualified Streaming.Prelude                  as Q
-import qualified Bio.Streaming.Bytes                as S
 
--- | 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 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:
---
--- * A name suffix of @/1@ or @/2@ is turned into the first mate or second
---   mate flag and the read is flagged as paired.
---
--- * Same for name prefixes of @F_@ or @R_@, respectively.
---
--- * A name prefix of @M_@ flags the sequence as unpaired and merged
---
--- * A name prefix of @T_@ flags the sequence as unpaired and trimmed
---
--- * 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.
---
--- 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
--- ignored.  The sequence can be split across multiple lines;
--- whitespace, dashes and dots are ignored, IUPAC-IUB ambiguity codes
--- are accepted as bases, anything else causes an error.  The sequence
--- ends at a line that is either a header or starts with @\+@, in the
--- latter case, that line is ignored and must be followed by quality
--- scores.  There must be exactly as many Q-scores as there are bases,
--- followed immediately by a header or end-of-file.  Whitespace is
--- ignored.
+-- | Emitted when random text is found instead of a header.
+data JunkFound = JunkFound !Int !Bytes deriving (Typeable, Show)
 
-parseFastq :: Monad m => ByteStream m r -> Stream (Of BamRec) m (Either SomeException r)
-parseFastq = parseFastqWith (const id)
+instance Exception JunkFound where
+    displayException (JunkFound n s) = printf "junk found at line %d: %s" n (unpack s)
 
--- | 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.
+-- | Emitted when a quality separator was expected, but not found.
+data QualitiesMissing = QualitiesMissing !Int !Bytes deriving (Typeable, Show)
 
-parseFastqCassava :: Monad m => ByteStream m r -> Stream (Of BamRec) m (Either SomeException r)
-parseFastqCassava = parseFastqWith (pdesc . C.split ':' . C.takeWhile (' ' /=))
+instance Exception QualitiesMissing where
+    displayException (QualitiesMissing 0 _) = printf "expected '+' symbol at end of file"
+    displayException (QualitiesMissing n s) = printf "expected '+' symbol at line %d, but found %s" n (unpack s)
+
+-- | Emitted when a quality record does not fit the sequence record.
+data IncoherentQualities = IncoherentQualities !Int !Bytes deriving (Typeable, Show)
+
+instance Exception IncoherentQualities where
+    displayException (IncoherentQualities n s) = printf "quality record of incorrect length ignored at line %d (%s)" n (unpack s)
+
+-- | Emitted when a quality record has different layout than the
+-- sequence.
+data IncongruentQualities = IncongruentQualities !Int !Bytes deriving (Typeable, Show)
+
+instance Exception IncongruentQualities where
+    displayException (IncongruentQualities n s) = printf "quality and sequence have different layouts at line %d (%s)" n (unpack s)
+
+-- | Emitted when a sequence record contains strange characters
+data SequenceHasGaps = SequenceHasGaps !Int !Bytes deriving (Typeable, Show)
+
+instance Exception SequenceHasGaps where
+    displayException (SequenceHasGaps n cs) = printf "undefined characters %s stored as unknown bases at line %d" (show cs) n
+
+data EmptyRecord = EmptyRecord !Int !Bytes deriving (Typeable, Show)
+
+instance Exception EmptyRecord where
+    displayException (EmptyRecord n s) = printf "(effectively) empty record at line %d (%s)" n (unpack s)
+
+
+{-# INLINE parseFastq #-}
+parseFastq :: MonadLog m => ByteStream m r -> Stream (Of BamRec) m r
+parseFastq = Q.unfoldr go . Q.zip (Q.enumFrom (1::Int)) . lines'
   where
-    pdesc (num:flg:_:idx:_) br = br { b_flag = sum [ if num == "1" then flagFirstMate .|. flagPaired else 0
-                                                   , if num == "2" then flagSecondMate .|. flagPaired else 0
-                                                   , if flg == "Y" then flagFailsQC else 0
-                                                   , b_flag br .&. complement (flagFailsQC .|. flagSecondMate .|. flagPaired) ]
-                                    , b_exts = if C.all (`C.elem` "ACGTN") idx then insertE "XI" (Text idx) (b_exts br) else b_exts br }
-    pdesc _ br = br
+    go = inspect >=> \case
+        Left r                                    ->  return (Left r)
+        Right ((i,h) :> ls)
+            | B.null h                            ->  go ls
+            | C.head h == '>' || C.head h == ';'  ->  goFasta (i, B.tail h) (Q.dropWhile isDescr ls)
+            | C.head h == '@'                     ->  goFastq (i, B.tail h) (Q.dropWhile isDescr ls)
+            | otherwise                           ->  logMsg Warning (JunkFound i h) >> go ls
 
--- | Same as 'parseFastq', but a custom function can be applied to the
--- description string (the part of the header after the sequence name),
--- which can modify the parsed record.  Note that the quality field can
--- end up empty.
+    isDescr  (_,s) = not (B.null s) &&  C.head s == ';'
+    isHeader (_,s) = not (B.null s) && (C.head s == '>' || C.head s == '@' || C.head s == ';')
+    isSep    (_,s) = not (B.null s) &&  C.head s == '+'
+    isSeq x = not (isHeader x || isSep x)
 
-parseFastqWith :: Monad m => (Bytes -> BamRec -> BamRec) -> ByteStream m r -> Stream (Of BamRec) m (Either SomeException r)
-parseFastqWith descr = Q.unfoldr (liftM twiddle . parse (const $ atto pRec)) . skipJunk
+    goFasta h ls = do
+        sq :> ls' <- Q.toList $ Q.span isSeq ls
+        make_record h sq Nothing ls'
+
+    goFastq h ls = do
+        sq :> ls1  <- Q.toList $ Q.span isSeq ls
+        Q.next ls1 >>= \case
+            Right (sep, ls2)
+                | isSep sep -> do
+                    qs :> ls3 <- Q.toList $ Q.splitAt (length sq) ls2
+                    if sum (map (B.length . snd) qs) /= sum (map (B.length . snd) sq) then do
+                        logMsg Error $ uncurry IncoherentQualities  h
+                        make_record h sq Nothing ls3
+                      else do
+                        when (map (B.length . snd) qs /= map (B.length . snd) sq) $
+                            logMsg Warning $ uncurry IncongruentQualities h
+                        make_record h sq (Just qs) ls3
+                | otherwise -> do
+                    logMsg Error $ uncurry QualitiesMissing sep
+                    make_record h sq Nothing (Q.cons sep ls2)
+            Left x -> do
+                    logMsg Error $ QualitiesMissing 0 C.empty
+                    make_record h sq Nothing (pure x)
+
+    make_record h sq qs k = do
+        when (ngaps > 0) . logMsg Warning $
+            SequenceHasGaps (fst $ head sq) (B.filter (isGap . toNucleotides) $ B.concat $ map snd sq)
+        when (l == 0) . logMsg Warning $ uncurry EmptyRecord h
+        return $ Right (r,k)
+      where
+        !l     = sum $ map (B.length . snd) sq
+        (!nseq, !ngaps) = mkSeq l sq
+        !qual  = case qs of Nothing -> Nothing
+                            Just [] -> Nothing
+                            Just  q -> Just $! mkQual l q
+
+        (!qname, !descr) = B.break (== 32) (snd h)
+        !fflag = B.drop 1 . C.dropWhile (/= ':') $ descr
+
+        !r = if B.length fflag < 2 || C.index fflag 1 /= ':' || (C.head fflag /= 'Y' && C.head fflag /= 'N')
+             then nullBamRec { b_qname = qname, b_seq = nseq, b_qual = qual }
+             else let !flag | C.head fflag /= 'Y' = b_flag nullBamRec .|. flagFailsQC
+                            | otherwise           = b_flag nullBamRec
+
+                      !sample = B.takeWhile (/=32). B.drop 1. C.dropWhile (/=':'). B.drop 1. C.dropWhile (/=':') $ fflag
+                      !exts | B.null sample                    =  [                   ]
+                            | C.all (`C.elem` "ACGTN") sample  =  [("XI", Text sample)]
+                            | otherwise                        =  [("RG", Text sample)]
+                  in nullBamRec { b_qname = qname, b_seq = nseq, b_qual = qual, b_flag = flag, b_exts = exts }
+
+
+mkSeq :: Int -> [(Int,Bytes)] -> (Vector_Nucs_half Nucleotides, Int)
+mkSeq ltot xs0 = unsafePerformIO $ do
+    fp <- mallocForeignPtrBytes (shiftR (ltot+1) 1)
+    g  <- withForeignPtr fp $ \p -> go_even p 0 xs0
+    return ( Vector_Nucs_half 0 ltot fp, g )
   where
-    twiddle        (Left  e)  = Left (Left  e)
-    twiddle (Right (Left  r)) = Left (Right r)
-    twiddle (Right (Right a)) = Right a
+    go_even !p !g (s:ss)  =  B.unsafeUseAsCStringLen (snd s) $ \(q,l) -> go1_even ss p q g l
+    go_even  _ !g [    ]  =  return g
 
-    isCBase   = inClass "ACGTUBDHVSWMKRYNacgtubdhvswmkryn"
-    canSkip c = isSpace c || c == '.' || c == '-'
-    isHdr   c = c == '@' || c == '>'
+    go_odd !p !g !a (s:ss)  =  B.unsafeUseAsCStringLen (snd s) $ \(q,l) -> go1_odd a ss p q g l
+    go_odd !p !g !a [    ]  =  poke p a >> return g
 
-    pRec   = (satisfy isHdr <?> "start marker") *> (makeRecord <$> pName <*> (descr <$> P.takeWhile ('\n' /=)) <*> (pSeq >>= pQual))
-    pName  = takeTill isSpace <* skipWhile (\c -> c /= '\n' && isSpace c)  <?> "read name"
-    pSeq   =     (:) <$> satisfy isCBase <*> pSeq
-             <|> satisfy canSkip *> pSeq
-             <|> pure []                                                   <?> "sequence"
+    go1_odd a ss !p !q !g !l
+        | l > 0      =  do !b <- unNs . toNucleotides <$> peekByteOff q 0
+                           poke p $ a `shiftL` 4 .|. b
+                           go1_even ss (plusPtr p 1) (plusPtr q 1) (g + fromEnum (b == 0)) (l-1)
+        | otherwise  =  go_odd p g a ss
 
-    pQual sq = (,) sq <$> (char '+' *> skipWhile ('\n' /=) *> pQual' (length sq) <* skipSpace <|> return C.empty)  <?> "qualities"
-    pQual' n = B.filter (not . isSpace_w8) <$> scan n step
-    step 0 _ = Nothing
-    step i c | isSpace c = Just i
-             | otherwise = Just (i-1)
+    go1_even ss !p !q !g !l
+        | l > 1      =  do !a <- unNs . toNucleotides <$> peekByteOff q 0
+                           !b <- unNs . toNucleotides <$> peekByteOff q 1
+                           poke p $ a `shiftL` 4 .|. b
+                           go1_even ss (plusPtr p 1) (plusPtr q 2) (g + fromEnum (a == 0) + fromEnum (b == 0)) (l-2)
+        | l > 0      =  do !a <- unNs . toNucleotides <$> peekByteOff q 0
+                           go_odd p (g + fromEnum (a == 0)) (a `shiftL` 4) ss
+        | otherwise  =  go_even p g ss
 
-skipJunk :: Monad m => ByteStream m r -> ByteStream m r
-skipJunk = lift . S.nextByte >=> check
+
+mkQual :: Int -> [(Int, Bytes)] -> W.Vector Qual
+mkQual ltot qs0 = unsafePerformIO $ do
+    fp <- mallocForeignPtrBytes ltot
+    withForeignPtr fp $ \p -> go (castPtr p) qs0
+    return $! W.unsafeFromForeignPtr0 fp ltot
   where
-    check (Right (c,s)) | bad c     = skipJunk . S.drop 1 . S.dropWhile (c2w '\n' /=) $ s
-                        | otherwise = S.cons c s
-    check (Left r)                  = pure r
-    bad c = c /= c2w '>' && c /= c2w '@'
+    go !p (s:ss)  =  B.unsafeUseAsCStringLen (snd s) $ \(q,l) -> go1 ss p (castPtr q) l
+    go  _ [    ]  =  return ()
 
-makeRecord :: Bytes -> (BamRec->BamRec) -> (String, Bytes) -> BamRec
-makeRecord name extra (sq,qual) = extra $ nullBamRec
-        { b_qname = name, b_seq = V.fromList $ read sq, b_qual = V.fromList $ map (Q . subtract 33) $ B.unpack qual }
+    go1 ss !p !q !l
+        | l > 0      =  do peek q >>= poke p . subtract (33::Word8)
+                           go1 ss (plusPtr p 1) (plusPtr q 1) (pred l)
+        | otherwise  =  go p ss
 
diff --git a/Bio/Bam/Filter.hs b/Bio/Bam/Filter.hs
--- a/Bio/Bam/Filter.hs
+++ b/Bio/Bam/Filter.hs
@@ -82,30 +82,32 @@
     total = fromIntegral $ V.length $ b_seq b
     ent   = sum [ fromIntegral c * log (total / fromIntegral c) | c <- counts, c /= 0 ] / log 2
 
--- | Filter on average quality.  Reads without quality string pass.
+-- | Filter on average quality.  Reads without qualities pass.
 {-# INLINE qualityAverage #-}
 qualityAverage :: Int -> QualFilter
 qualityAverage q b = if p then b else b'
   where
     b' = setQualFlag 'Q' $ b { b_flag = b_flag b .|. flagFailsQC }
-    p  = let total = V.foldl' (\a x -> a + fromIntegral (unQ x)) 0 $ b_qual b
-         in total >= q * V.length (b_qual b)
+    p  = case b_qual b of Nothing -> True
+                          Just qs -> let total = V.foldl' (\a x -> a + fromIntegral (unQ x)) 0 qs
+                                     in total >= q * V.length qs
 
 -- | Filter on minimum quality.  In @qualityMinimum n q@, a read passes
 -- if it has no more than @n@ bases with quality less than @q@.  Reads
--- without quality string pass.
+-- without qualities pass.
 {-# INLINE qualityMinimum #-}
 qualityMinimum :: Int -> Qual -> QualFilter
 qualityMinimum n (Q q) b = if p then b else b'
   where
     b' = setQualFlag 'Q' $ b { b_flag = b_flag b .|. flagFailsQC }
-    p  = V.length (V.filter (< Q q) (b_qual b)) <= n
+    p  = case b_qual b of Nothing -> True
+                          Just qs -> V.length (V.filter (< Q q) qs) <= n
 
 
 -- | Convert quality scores from old Illumina scale (different formula
 -- and offset 64 in FastQ).
 qualityFromOldIllumina :: BamRec -> BamRec
-qualityFromOldIllumina b = b { b_qual = V.map conv $ b_qual b }
+qualityFromOldIllumina b = b { b_qual = V.map conv <$> b_qual b }
   where
     conv (Q s) = let s' :: Double
                      s' = exp $ log 10 * (fromIntegral s - 31) / (-10)
@@ -116,6 +118,6 @@
 -- | Convert quality scores from new Illumina scale (standard formula
 -- but offset 64 in FastQ).
 qualityFromNewIllumina :: BamRec -> BamRec
-qualityFromNewIllumina b = b { b_qual = V.map (Q . subtract 31 . unQ) $ b_qual b }
+qualityFromNewIllumina b = b { b_qual = V.map (Q . subtract 31 . unQ) <$> b_qual b }
 
 
diff --git a/Bio/Bam/Header.hs b/Bio/Bam/Header.hs
--- a/Bio/Bam/Header.hs
+++ b/Bio/Bam/Header.hs
@@ -56,6 +56,7 @@
 import Data.ByteString.Builder      ( Builder, byteString, char7, intDec, word16LE )
 
 import qualified Data.Attoparsec.ByteString.Char8   as P
+import qualified Data.ByteString                    as B
 import qualified Data.ByteString.Char8              as S
 import qualified Data.HashMap.Strict                as H
 import qualified Data.Vector                        as V
@@ -88,8 +89,8 @@
 -- (arbitrarily) prepended to the first existing chain, or forms a new
 -- singleton chain if none exists.
 
-addPG :: Maybe Version -> IO (BamMeta -> BamMeta)
-addPG vn = do
+addPG :: MonadIO m => Maybe Version -> m (BamMeta -> BamMeta)
+addPG vn = liftIO $ do
     args <- getArgs
     pn   <- getProgName
 
@@ -168,9 +169,6 @@
 
 instance Hashable BamSQ
 
-bad_seq :: BamSQ
-bad_seq = BamSQ (error "no SN field") (error "no LN field") []
-
 data BamPG pp = BamPG {
         pg_pref_name :: Bytes,
         pg_prev_pg :: Maybe pp,
@@ -185,13 +183,9 @@
 instance Show (f (Fix f)) => Show (Fix f) where
     showsPrec p (Fix f) = showsPrec p f
 
-instance Hashable1 f => Hashable (Fix f) where
-    hashWithSalt salt (Fix f) = liftHashWithSalt hashWithSalt salt f
-
-instance Hashable1 BamPG
-
-bad_pg :: BamPG pp
-bad_pg = BamPG (error "no ID field") Nothing []
+instance Hashable (Fix BamPG) where
+    hashWithSalt s (Fix (BamPG n Nothing  o)) = hashWithSalt               (hashWithSalt s n)    o
+    hashWithSalt s (Fix (BamPG n (Just p) o)) = hashWithSalt (hashWithSalt (hashWithSalt s n) p) o
 
 
 -- | Possible sorting orders from bam header.  Thanks to samtools, which
@@ -219,7 +213,7 @@
 
 parseBamMeta :: P.Parser BamMeta
 parseBamMeta = fixupMeta . foldl' (flip ($)) emptyHeader
-               <$> many (parseBamMetaLine <* P.skipWhile (=='\t') <* P.char '\n')
+               <$> many (parseBamMetaLine <* P.skipWhile (=='\t') <* P.char '\n') <* P.endOfInput
 
 -- Bam header in the process of being parsed.  Better suited for
 -- collecting lines than 'BamMeta'.
@@ -234,13 +228,13 @@
 emptyHeader = PreBamMeta mempty [] H.empty [] []
 
 
--- | Fixes a bam header after parsing.  It turns accumulated lists in to
--- vectors, and it handles the program lines.  Program lines come in as
--- an arbitrary graph.  It chould be a linear chain, but this
+-- | Fixes a bam header after parsing.  It turns accumulated lists into
+-- vectors, throws errors for mandatory fields that weren't parsed
+-- correctly, and it handles the program (PG) lines.  Program lines come
+-- in as an arbitrary graph.  It should be a linear chain, but this
 -- isn't guaranteed in practice.  We decompose the graph into chains by
 -- tracing from nodes with no predecessor, or from an arbitrary node if
--- all nodes have predecessors.  Tracing stops once it would form a
--- cycle.
+-- all nodes have predecessors.  Tracing stops if it would form a cycle.
 fixupMeta :: PreBamMeta -> BamMeta
 fixupMeta PreBamMeta{..} = BamMeta
     { meta_hdr        = pmeta_hdr
@@ -287,14 +281,18 @@
              (\fns meta -> meta { pmeta_hdr = foldr ($) (pmeta_hdr meta) fns })
                <$> P.sepBy1 (P.choice [hdvn, hdso, hdother]) tabs
 
-    sqLine = P.string "SQ\t" >>
-             (\fns meta -> meta { pmeta_refs = foldr ($) bad_seq fns : pmeta_refs meta })
-               <$> P.sepBy1 (P.choice [sqnm, sqln, sqother]) tabs
+    sqLine = do _ <- P.string "SQ\t"
+                fns <- P.sepBy1 (P.choice [sqnm, sqln, sqother]) tabs
+                let sq = foldr ($) (BamSQ "" (-1) []) fns
+                guard (not . B.null $ sq_name sq) P.<?> "SQ:NM field"
+                guard (sq_length sq >= 0) P.<?> "SQ:LN field"
+                pure $ \meta -> meta { pmeta_refs = sq : pmeta_refs meta }
 
-    pgLine = P.string "PG\t" >>
-             (\fns meta -> let pg = foldr ($) bad_pg fns
-                           in meta { pmeta_pgs = H.insert (pg_pref_name pg) pg (pmeta_pgs meta) })
-               <$> P.sepBy1 (P.choice [pgid, pgpp, pgother]) tabs
+    pgLine = do _ <- P.string "PG\t"
+                fns <- P.sepBy1 (P.choice [pgid, pgpp, pgother]) tabs
+                let pg = foldr ($) (BamPG "" Nothing []) fns
+                guard (not . B.null $ pg_pref_name pg) P.<?> "PG:ID field"
+                pure $ \meta -> meta { pmeta_pgs = H.insert (pg_pref_name pg) pg (pmeta_pgs meta) }
 
     hdvn = P.string "VN:" >>
            (\a b hdr -> hdr { hdr_version = (a,b) })
@@ -415,7 +413,7 @@
 -- | Reference sequence in Bam
 -- Bam enumerates the reference sequences and then sorts by index.  We
 -- need to track that index if we want to reproduce the sorting order.
-newtype Refseq = Refseq { unRefseq :: Word32 } deriving (Eq, Ord, Ix, Bounded)
+newtype Refseq = Refseq { unRefseq :: Word32 } deriving (Eq, Ord, Ix, Bounded, Hashable)
 
 instance Show Refseq where
     showsPrec p (Refseq r) = showsPrec p r
@@ -539,8 +537,8 @@
          | isDigit (S.head s) = do (n,t) <- S.readInt s
                                    (MdNum n :) <$> readMd t
          | S.head s == '^'    = let (a,b) = S.break isDigit (S.tail s)
-                                in (MdDel (map toNucleotides $ S.unpack a) :) <$> readMd b
-         | otherwise          = (MdRep (toNucleotides $ S.head s) :) <$> readMd (S.tail s)
+                                in (MdDel (map toNucleotides $ B.unpack a) :) <$> readMd b
+         | otherwise          = (MdRep (toNucleotides $ B.head s) :) <$> readMd (S.tail s)
 
 -- | Normalizes a series of 'MdOp's and encodes them in the way BAM and
 -- SAM expect it.
diff --git a/Bio/Bam/Index.hs b/Bio/Bam/Index.hs
--- a/Bio/Bam/Index.hs
+++ b/Bio/Bam/Index.hs
@@ -4,6 +4,7 @@
     readBamIndex,
     readBaiIndex,
     readTabix,
+    IndexFormatError(..),
 
     Region(..),
     Subsequence(..),
@@ -143,11 +144,10 @@
 xs ~~ [] = xs
 
 
-data IndexFormatError = IndexFormatError Bytes deriving Typeable
+data IndexFormatError = IndexFormatError Bytes deriving (Typeable, Show)
 
-instance Exception IndexFormatError
-instance Show IndexFormatError where
-    show (IndexFormatError m) = "index signature " ++ show m ++ " not recognized"
+instance Exception IndexFormatError where
+    displayException (IndexFormatError m) = "index signature " ++ show m ++ " not recognized"
 
 {- | Reads any index we can find for a file.
 
@@ -204,7 +204,7 @@
                    | otherwise = llim bin (dp-3) (sf+3)
             where ix = (1 `shiftL` dp - 1) `div` 7
 
-withIndexedBam :: (MonadIO m, MonadMask m) => FilePath -> (BamMeta -> BamIndex () -> Handle -> m r) -> m r
+withIndexedBam :: (MonadIO m, MonadLog m, MonadMask m) => FilePath -> (BamMeta -> BamIndex () -> Handle -> m r) -> m r
 withIndexedBam f k = do
     idx <- liftIO $ readBamIndex f
     bracket (liftIO $ openBinaryFile f ReadMode) (liftIO . hClose) $ \hdl -> do
@@ -322,15 +322,16 @@
 'BamRaw' records of the correct reference sequence only, and produces an
 empty stream if the sequence isn't found.
 -}
-streamBamRefseq :: MonadIO m => BamIndex b -> Handle -> Refseq -> Stream (Of BamRaw) m ()
+streamBamRefseq :: (MonadIO m, MonadLog m) => BamIndex b -> Handle -> Refseq -> Stream (Of BamRaw) m ()
 streamBamRefseq BamIndex{..} hdl (Refseq r)
     | Just ckpts <- refseq_ckpoints V.!? fromIntegral r
     , Just (voff, _) <- M.minView ckpts
     , voff /= 0 = void $
                   Q.takeWhile ((Refseq r ==) . b_rname . unpackBam) $
-                  Q.unfoldr (P.parseIO getBamRaw) $
+                  Q.unfoldr (P.parseLog Error getBamRaw) $
                   streamBgzf hdl voff Nothing
     | otherwise = pure ()
+
 
 {- | Reads from a Bam file the part with unaligned reads.
 
diff --git a/Bio/Bam/Pileup.hs b/Bio/Bam/Pileup.hs
--- a/Bio/Bam/Pileup.hs
+++ b/Bio/Bam/Pileup.hs
@@ -115,7 +115,7 @@
               | 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   = maybe (Q 23) (`V.unsafeIndex` i) b_qual                          -- 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
@@ -270,7 +270,7 @@
 -- encountered or a t end of file.
 
 {-# INLINE pileup #-}
-pileup :: Stream (Of PosPrimChunks) IO b -> Stream (Of Pile) IO b
+pileup :: Monad m => Stream (Of PosPrimChunks) m b -> Stream (Of Pile) m b
 pileup = runPileM pileup' finish (Refseq 0) 0 ([],[]) (Empty,Empty)
   where
     finish () _ _ ([],[]) (Empty,Empty) inp = lift (Q.effects inp)
@@ -318,15 +318,6 @@
     return a = PileM $ \k -> k a
     {-# INLINE (>>=) #-}
     m >>=  k = PileM $ \k' -> runPileM m (\a -> runPileM (k a) k')
-
---  -- XXX
--- {-# INLINE get_refseq #-}
--- get_refseq :: PileM m Refseq
--- get_refseq = PileM $ \k r -> k r r
-
--- {-# INLINE get_pos #-}
--- get_pos :: PileM m Int
--- get_pos = PileM $ \k r p -> k p r p
 
 {-# INLINE upd_pos #-}
 upd_pos :: (Int -> Int) -> PileM m ()
diff --git a/Bio/Bam/Reader.hs b/Bio/Bam/Reader.hs
--- a/Bio/Bam/Reader.hs
+++ b/Bio/Bam/Reader.hs
@@ -4,6 +4,7 @@
     decodeBam,
     decodeBamFile,
     decodeBamFiles,
+    IncompatibleRefs(..),
 
     decodePlainBam,
     decodePlainSam,
@@ -43,7 +44,7 @@
 is reliably recognized, anything else is treated as SAM.  The offsets
 stored in BAM records make sense only for uncompressed or bgzf'd BAM.
 -}
-decodeBam :: MonadIO m
+decodeBam :: (MonadIO m, MonadLog m)
           => S.ByteStream m r
           -> m (BamMeta, Stream (Of BamRaw) m r)
 decodeBam = getBgzfHdr >=> S.splitAt' 4 . pgunzip >=> unbam
@@ -55,7 +56,7 @@
     pgunzip (Just _,  hdr, s) =  bgunzip (S.consChunk hdr s)
 {-# INLINE decodeBam #-}
 
-decodeBamFile :: (MonadIO m, MonadMask m) => FilePath -> (BamMeta -> Stream (Of BamRaw) m () -> m r) -> m r
+decodeBamFile :: (MonadIO m, MonadLog m, MonadMask m) => FilePath -> (BamMeta -> Stream (Of BamRaw) m () -> m r) -> m r
 decodeBamFile f k = streamFile f $ decodeBam >=> uncurry k
 {-# INLINE decodeBamFile #-}
 
@@ -66,17 +67,18 @@
 unrelated bam files.  All files are opened at the same time, which might
 run into the file descriptor limit given some ridiculous workflows.
 -}
-decodeBamFiles :: (MonadMask m, MonadIO m) => [FilePath] -> ([(BamMeta, Stream (Of BamRaw) m ())] -> m r) -> m r
+decodeBamFiles :: (MonadMask m, MonadLog m, MonadIO m) => [FilePath] -> ([(BamMeta, Stream (Of BamRaw) m ())] -> m r) -> m r
 decodeBamFiles [      ] k = k []
 decodeBamFiles ("-":fs) k = decodeBam (streamHandle stdin)   >>= \b -> decodeBamFiles fs $ \bs -> k (b:bs)
 decodeBamFiles ( f :fs) k = streamFile f $ \s -> decodeBam s >>= \b -> decodeBamFiles fs $ \bs -> k (b:bs)
 {-# INLINE decodeBamFiles #-}
 
-decodePlainBam :: MonadIO m => S.ByteStream m r -> m (BamMeta, Stream (Of BamRaw) m r)
+decodePlainBam :: MonadLog m => S.ByteStream m r -> m (BamMeta, Stream (Of BamRaw) m r)
 decodePlainBam =
-    S.parseIO (const getBamMeta) >=>
-    either (const . liftIO $ throwM S.EofException)
-           (return . fmap (Q.unfoldr (S.parseIO getBamRaw)))
+    S.parse (const getBamMeta) >=> \case
+        Left (exception, rest) -> logMsg Error exception >> S.effects rest >>= \r -> pure (mempty, pure r)
+        Right (Left r)         -> logMsg Error S.EofException >> pure (mempty, pure r)
+        Right (Right (h,s))    -> return (h, Q.unfoldr (S.parseLog Error getBamRaw) s)
 
 getBamMeta :: Monad m => S.Parser r m BamMeta
 getBamMeta = liftM2 mmerge get_bam_header get_ref_array
@@ -103,42 +105,41 @@
                 | otherwise                  = l -- contradiction in header, but we'll just ignore it
 {-# INLINABLE getBamMeta #-}
 
-
-data ShortRecord = ShortRecord deriving Typeable
-
-instance Exception ShortRecord
-instance Show ShortRecord where show _ = "short BAM record"
-
 getBamRaw :: Monad m => Int64 -> S.Parser r m BamRaw
 getBamRaw o = do
         bsize <- fromIntegral `liftM` S.getWord32
-        when (bsize < 32) $ throwM ShortRecord
         s <- S.getString bsize
         unless (B.length s == bsize) S.abortParse
-        return $ bamRaw o s
+        bamRaw o s
 {-# INLINABLE getBamRaw #-}
 
 {- | Streaming parser for SAM files.
 
-It parses plain uncompressed SAM
-and returns a result compatible with 'decodePlainBam'.  Since it is
-supposed to work the same way as the BAM parser, it requires a
-symbol table for the reference names.  This is extracted from the @SQ
-lines in the header.  Note that reading SAM tends to be inefficient;
-if you care about performance at all, use BAM.
--}
-decodePlainSam :: MonadIO m => S.ByteStream m r -> m (BamMeta, Stream (Of BamRaw) m r)
+It parses plain uncompressed SAM and returns a result compatible with
+'decodePlainBam'.  Since it is supposed to work the same way as the BAM
+parser, it requires a symbol table for the reference names.  This is
+extracted from the @SQ lines in the header.  Note that reading SAM tends
+to be inefficient; if you care about performance at all, use BAM.  -}
+
+decodePlainSam :: (MonadLog m, MonadIO m) => S.ByteStream m r -> m (BamMeta, Stream (Of BamRaw) m r)
 decodePlainSam s = do
     (hdr,rest) <- either (\r -> (mempty, pure r)) id `liftM` S.parseIO (const $ S.atto parseBamMeta) s
-    let !refs = M.fromList $ zip [ nm | BamSQ { sq_name = nm } <- V.toList (unRefs (meta_refs hdr))] [toEnum 0..]
-        ref x = M.lookupDefault invalidRefseq x refs
-    return (hdr, Q.mapM (liftIO . either throwM packBam . getSamRec ref) (S.lines' rest))
+    let !refs  = M.fromList $ zip [ nm | BamSQ { sq_name = nm } <- V.toList (unRefs (meta_refs hdr))] [toEnum 0..]
+        ref  x = M.lookupDefault invalidRefseq x refs
+        report = fmap (const Nothing) . logMsg Error
+        use    = fmap Just . liftIO . packBam
+        strm   = Q.concat . Q.mapM (either report use <=< getSamRec ref) $ S.lines' rest
+    return (hdr, strm)
 
 
-getSamRec :: (Bytes -> Refseq) -> Bytes -> Either S.ParseError BamRec
+getSamRec :: MonadLog m => (Bytes -> Refseq) -> Bytes -> m (Either S.ParseError BamRec)
 getSamRec ref s = case P.parseOnly record s of
-    Left  e -> Left $ S.ParseError [unpack s] e
-    Right b -> Right b
+    Left  e                                         -> pure . Left $ S.ParseError [unpack s] e
+    Right b -> case b_qual b of
+        Nothing                                     -> pure $ Right b
+        Just qs | W.length qs == V.length (b_seq b) -> pure $ Right b
+                | otherwise                         -> do logMsg Warning $ LengthMismatch (b_qname b)
+                                                          pure . Right $ b { b_qual = Nothing }
   where
     record = do b_qname <- word
                 b_flag  <- num
@@ -150,7 +151,7 @@
                 b_mpos  <- subtract 1 <$> num
                 b_isize <- snum
                 b_seq   <- sequ
-                b_qual  <- quals <*> pure b_seq
+                b_qual  <- quals
                 b_exts  <- exts
                 let b_virtual_offset = 0
                 return BamRec{..}
@@ -163,12 +164,11 @@
 
     rnext    = id <$ P.char '=' <* sep <|> const . ref <$> word
     sequ     = (V.empty <$ P.char '*' <|>
-               V.fromList . map toNucleotides . C.unpack <$> P.takeWhile is_nuc) <* sep
+               V.fromList . map toNucleotides . B.unpack <$> P.takeWhile is_nuc) <* sep
 
-    quals    = defaultQs <$ P.char '*' <* sep <|> bsToVec <$> word
+    quals    = Nothing <$ P.char '*' <* sep <|> bsToVec <$> word
         where
-            defaultQs sq = W.replicate (V.length sq) (Q 0xff)
-            bsToVec qs _ = W.fromList . map (Q . subtract 33) $ B.unpack qs
+            bsToVec = Just . W.fromList . map (Q . subtract 33) . B.unpack
 
     cigar    = [] <$ P.char '*' <* sep <|>
                P.manyTill (flip (:*) <$> P.decimal <*> cigop) sep
@@ -193,11 +193,10 @@
     is_nuc = P.inClass "acgtswkmrybdhvnACGTSWKMRYBDHVN"
 
 
-data IncompatibleRefs = IncompatibleRefs FilePath FilePath deriving Typeable
+data IncompatibleRefs = IncompatibleRefs FilePath FilePath deriving (Typeable, Show)
 
-instance Exception IncompatibleRefs
-instance Show IncompatibleRefs where
-    show (IncompatibleRefs a b) = "references in " ++ a ++ " and " ++ b ++ " are incompatible"
+instance Exception IncompatibleRefs where
+    displayException (IncompatibleRefs a b) = "references in " ++ a ++ " and " ++ b ++ " are incompatible"
 
 guardRefCompat :: MonadThrow m => (FilePath,BamMeta) -> (FilePath,BamMeta) -> m ()
 guardRefCompat (f0,hdr0) (f1,hdr1) =
@@ -216,7 +215,7 @@
 for the result, and an exception is thrown if one of the subsequent
 headers is incompatible with the first one.
 -}
-concatInputs :: (MonadIO m, MonadMask m) => [FilePath] -> (BamMeta -> Stream (Of BamRaw) m () -> m r) -> m r
+concatInputs :: (MonadIO m, MonadLog m, MonadMask m) => [FilePath] -> (BamMeta -> Stream (Of BamRaw) m () -> m r) -> m r
 concatInputs fs0 k = streamInputs fs0 (go1 $ fs0 ++ repeat "-")
   where
     go1 fs = inspect >=> \case
@@ -241,7 +240,7 @@
 exception unless every input is compatible with the effective reference
 list.
 -}
-mergeInputsOn :: (Ord x, MonadIO m, MonadMask m)
+mergeInputsOn :: (Ord x, MonadIO m, MonadLog m, MonadMask m)
               => (BamRaw -> x) -> [FilePath]
               -> (BamMeta -> Stream (Of BamRaw) m () -> m r) -> m r
 mergeInputsOn _ [] k = decodeBam (streamHandle stdin) >>= uncurry k
diff --git a/Bio/Bam/Rec.hs b/Bio/Bam/Rec.hs
--- a/Bio/Bam/Rec.hs
+++ b/Bio/Bam/Rec.hs
@@ -6,17 +6,20 @@
     bamRaw,
     virt_offset,
     raw_data,
+    br_copy,
 
     BamRec(..),
     unpackBam,
     nullBamRec,
     getMd,
+    LengthMismatch(..),
+    BrokenRecord(..),
 
     Cigar(..),
     CigOp(..),
     alignedLength,
 
-    Nucleotides(..), Vector_Nucs_half,
+    Nucleotides(..), Vector_Nucs_half(..),
     Extensions, Ext(..),
     extAsInt, extAsString, setQualFlag,
     deleteE, insertE, updateE, adjustE,
@@ -29,7 +32,6 @@
     isMateReversed,
     isFirstMate,
     isSecondMate,
-    isAuxillary,
     isSecondary,
     isFailsQC,
     isDuplicate,
@@ -91,7 +93,9 @@
   where l (op :* n) = if op == Mat || op == Del || op == Nop then n else 0
 
 
--- | internal representation of a BAM record
+-- | More convenient representation of a BAM record.
+-- Invariant: Either @b_qual == Nothing@ or
+-- @fmap V.length b_qual == Just (V.length b_seq)@.
 data BamRec = BamRec {
         b_qname :: Bytes,
         b_flag  :: Int,
@@ -103,7 +107,7 @@
         b_mpos  :: Int,
         b_isize :: Int,
         b_seq   :: Vector_Nucs_half Nucleotides,
-        b_qual  :: VS.Vector Qual,
+        b_qual  :: Maybe (VS.Vector Qual),
         b_exts  :: Extensions,
         b_virtual_offset :: Int64 -- ^ virtual offset for indexing purposes
     } deriving Show
@@ -120,7 +124,7 @@
         b_mpos  = invalidPos,
         b_isize = 0,
         b_seq   = V.empty,
-        b_qual  = VS.empty,
+        b_qual  = Nothing,
         b_exts  = [],
         b_virtual_offset = 0
     }
@@ -131,6 +135,14 @@
     Just (Char mdfield) -> readMd $ B.singleton mdfield
     _                   -> Nothing
 
+data LengthMismatch = LengthMismatch !Bytes deriving (Typeable, Show)
+instance Exception LengthMismatch where
+    displayException (LengthMismatch nm) =
+        "length of nucleotide and quality sequences" ++
+        (if B.null nm then "" else " in " ++ unpack nm) ++
+        " must be equal"
+
+
 -- | A vector that packs two 'Nucleotides' into one byte, just like Bam does.
 data Vector_Nucs_half a = Vector_Nucs_half !Int !Int !(ForeignPtr Word8)
 
@@ -204,13 +216,22 @@
 data BamRaw = BamRaw { virt_offset :: {-# UNPACK #-} !Int64
                      , raw_data    :: {-# UNPACK #-} !Bytes }
 
+br_copy :: BamRaw -> BamRaw
+br_copy (BamRaw o r) = BamRaw o $! B.copy r
+
+data BrokenRecord = BrokenRecord !Int [Int] !Bytes deriving (Typeable, Show)
+instance Exception BrokenRecord where
+    displayException (BrokenRecord ln m raw) =
+        "broken BAM record of length " ++ shows ln " , need offsets " ++ shows m ", begins with " ++ show (B.unpack raw)
+
 -- | Smart constructor.  Makes sure we got a at least a full record.
 {-# INLINE bamRaw #-}
-bamRaw :: Int64 -> Bytes -> BamRaw
-bamRaw o s = if good then BamRaw o s else error $ "broken BAM record " ++ shows (S.length s, m) " " ++ show (S.unpack (S.take 10 s))
+bamRaw :: MonadThrow m => Int64 -> Bytes -> m BamRaw
+bamRaw o s
+    | S.length s <    32  =  throwM $ BrokenRecord (S.length s) [32] (S.take 16 s)
+    | S.length s < sum m  =  throwM $ BrokenRecord (S.length s)   m  (S.take 16 s)
+    | otherwise           =  pure $ BamRaw o s
   where
-    good | S.length s < 32 = False
-         | otherwise       = S.length s >= sum m
     m = [ 32, l_rnm, l_seq, (l_seq+1) `div` 2, l_cig * 4 ]
     l_rnm = fromIntegral (B.unsafeIndex s  8) - 1
     l_cig = fromIntegral (B.unsafeIndex s 12)             .|. fromIntegral (B.unsafeIndex s 13) `shiftL`  8
@@ -231,8 +252,7 @@
         b_qname = B.unsafeTake l_read_name $ B.unsafeDrop 32 $ raw_data br,
         b_cigar = VS.unsafeCast $ VS.unsafeFromForeignPtr fp (off0+off_c) (4*l_cigar),
         b_seq   = Vector_Nucs_half (2 * (off_s+off0)) l_seq fp,
-        b_qual  = VS.unsafeCast $ VS.unsafeFromForeignPtr fp (off0+off_q) l_seq,
-
+        b_qual  = mk_qual,
         b_exts  = unpackExtensions $ S.drop off_e $ raw_data br,
         b_virtual_offset = virt_offset br }
   where
@@ -246,6 +266,10 @@
         l_seq       = getWord32 16
         l_cigar     = getInt16  12
 
+        mk_qual | l_seq == 0                                 =  Just VS.empty
+                | B.unsafeIndex (raw_data br) off_q == 0xff  =  Nothing
+                | otherwise  =  Just $ VS.unsafeFromForeignPtr (castForeignPtr fp) (off0+off_q) l_seq
+
         getInt8 :: Num a => Int -> a
         getInt8  o = fromIntegral (B.unsafeIndex (raw_data br) o)
 
@@ -261,6 +285,7 @@
         getInt32 :: Num a => Int -> a
         getInt32 o = fromIntegral (getWord32 o :: Int32)
 
+
 -- | A collection of extension fields.  A 'BamKey' is actually two ASCII
 -- characters.
 type Extensions = [( BamKey, Ext )]
@@ -335,7 +360,7 @@
 
 
 isPaired, isProperlyPaired, isUnmapped, isMateUnmapped, isReversed,
-    isMateReversed, isFirstMate, isSecondMate, isAuxillary, isSecondary,
+    isMateReversed, isFirstMate, isSecondMate, isSecondary,
     isFailsQC, isDuplicate, isSupplementary,
     isTrimmed, isMerged, isAlternative, isExactIndex :: BamRec -> Bool
 
@@ -347,7 +372,6 @@
 isMateReversed   = flip testBit  5 . b_flag
 isFirstMate      = flip testBit  6 . b_flag
 isSecondMate     = flip testBit  7 . b_flag
-isAuxillary      = flip testBit  8 . b_flag
 isSecondary      = flip testBit  8 . b_flag
 isFailsQC        = flip testBit  9 . b_flag
 isDuplicate      = flip testBit 10 . b_flag
diff --git a/Bio/Bam/Rmdup.hs b/Bio/Bam/Rmdup.hs
--- a/Bio/Bam/Rmdup.hs
+++ b/Bio/Bam/Rmdup.hs
@@ -2,7 +2,8 @@
             rmdup, Collapse, cons_collapse, cheap_collapse,
             cons_collapse_keep, cheap_collapse_keep,
             check_sort, normalizeTo, wrapTo,
-            ECig(..), toECig, setMD, toCigar
+            ECig(..), toECig, setMD, toCigar,
+            RmdupException(..), BrokenMD(..)
     ) where
 
 import Bio.Bam.Header
@@ -10,31 +11,33 @@
 import Bio.Prelude hiding ( left, right )
 import Bio.Streaming
 
-import qualified Data.ByteString        as B
-import qualified Data.ByteString.Char8  as T
-import qualified Data.Map.Strict        as M
-import qualified Data.Vector.Generic    as V
-import qualified Data.Vector.Storable   as W
-import qualified Data.Vector.Unboxed    as U
-import qualified Streaming.Prelude      as Q
+import qualified Data.ByteString                as B
+import qualified Data.ByteString.Char8          as T
+import qualified Data.HashMap.Strict            as M
+import qualified Data.Vector                    as VV ( fromList )
+import qualified Data.Vector.Algorithms.Intro   as VV ( sortBy )
+import qualified Data.Vector.Generic            as V
+import qualified Data.Vector.Storable           as W
+import qualified Data.Vector.Unboxed            as U
+import qualified Streaming.Prelude              as Q
 
-data Collapse = Collapse {
-        collapse  :: [BamRec] -> (Decision,[BamRec]),    -- cluster to consensus and stuff or representative and stuff
+data Collapse m = Collapse {
+        collapse  :: [BamRec] -> m (Decision,[BamRec]),  -- cluster to consensus and stuff or representative and stuff
         originals :: [BamRec] -> [BamRec] }              -- treatment of the redundant original reads
 
 data Decision = Consensus      { fromDecision :: BamRec }
               | Representative { fromDecision :: BamRec }
 
-cons_collapse :: Qual -> Collapse
+cons_collapse :: MonadLog m => Qual -> Collapse m
 cons_collapse maxq = Collapse (do_collapse maxq) (const [])
 
-cons_collapse_keep :: Qual -> Collapse
+cons_collapse_keep :: MonadLog m => Qual -> Collapse m
 cons_collapse_keep maxq = Collapse (do_collapse maxq) (map (\b -> b { b_flag = b_flag b .|. flagDuplicate }))
 
-cheap_collapse :: Collapse
+cheap_collapse :: Monad m => Collapse m
 cheap_collapse = Collapse do_cheap_collapse (const [])
 
-cheap_collapse_keep :: Collapse
+cheap_collapse_keep :: Monad m => Collapse m
 cheap_collapse_keep = Collapse do_cheap_collapse (map (\b -> b { b_flag = b_flag b .|. flagDuplicate }))
 
 
@@ -103,17 +106,18 @@
 -- duplicates of each other.  The typical label function would extract
 -- read groups, libraries or samples.
 
-rmdup :: (Monad m, Ord l) => (BamRec -> l) -> Bool -> Collapse
-      -> Q.Stream (Of BamRec) m r -> Q.Stream (Of (Int,BamRec)) m r
-rmdup label strand_preserved collapse_cfg =
+rmdup :: (MonadLog m, MonadThrow m, Hashable l, Eq l)
+      => Bool -> Collapse m -> Q.Stream (Of (l,BamRec)) m r -> Q.Stream (Of ((l,Int),BamRec)) m r
+rmdup strand_preserved collapse_cfg =
     -- Easiest way to go about this:  We simply collect everything that
     -- starts at some specific coordinate and group it appropriately.
     -- Treat the groups separately, output, go on.
-    check_sort (b_cpos . snd) "internal error, output isn't sorted anymore" .
-    mapGroups ( sortBy (comparing $ V.length . b_seq . snd)
-              . concatMap (do_rmdup strand_preserved collapse_cfg)
-              . M.elems . accumMap label id ) .
-    check_sort b_cpos "input must be sorted for rmdup to work"
+    check_sort (b_cpos . snd) InternalError .
+    mapGroups ( fmap (V.modify (VV.sortBy (comparing $ V.length . b_seq . snd)) . VV.fromList)
+              . fmap (M.foldrWithKey (\l mrs -> (++) [ ((l,n),r) | (n,r) <- mrs ]) [])
+              . mapM (do_rmdup strand_preserved collapse_cfg)
+              . accumMap fst snd ) .
+    check_sort (b_cpos . snd) UnsortedInput
   where
     b_cpos = b_rname &&& b_pos
 
@@ -121,17 +125,22 @@
 
     mg1 f acc a s =
         lift (inspect s) >>= \case
-            Left r  -> Q.each (f (a:acc)) >> pure r
+            Left r                                  ->  lift (f (a:acc)) >>= Q.each >> pure r
             Right (b :> s')
-                | b_cpos a == b_cpos b -> mg1 f (b:acc) a s'
-                | otherwise            -> Q.each (f (a:acc)) >> mg1 f [] b s'
+                | b_cpos (snd a) == b_cpos (snd b)  ->  mg1 f (b:acc) a s'
+                | otherwise                         ->  lift (f (a:acc)) >>= Q.each >> mg1 f [] b s'
 
-check_sort :: (Monad m, Ord b) => (a -> b) -> String -> Stream (Of a) m r -> Stream (Of a) m r
+data RmdupException = UnsortedInput | InternalError deriving (Typeable, Show)
+instance Exception RmdupException where
+    displayException UnsortedInput = "input must be sorted for rmdup to work"
+    displayException InternalError = "internal error, output isn't sorted anymore"
+
+check_sort :: (MonadThrow m, Ord b) => (a -> b) -> RmdupException -> Stream (Of a) m r -> Stream (Of a) m r
 check_sort f msg = lift . inspect >=> either pure step
   where
     step (a :> s) = lift (inspect s) >>= either (Q.cons a . pure) (step' a)
     step' a (b :> s)
-        | f a > f b = fail $ "rmdup: " ++ msg
+        | f a > f b = lift $ throwM msg
         | otherwise = Q.cons a $ step (b :> s)
 {-# INLINE check_sort #-}
 
@@ -196,33 +205,36 @@
    (4) See 'merge_singles' for how it's actually done.
 -}
 
-do_rmdup :: Bool -> Collapse -> [BamRec] -> [(Int,BamRec)]
-do_rmdup strand_preserved Collapse{..} rds =
-        results ++ map ((,) 0) (originals (leftovers ++ r1 ++ r2 ++ r3))
-    where
-        (results, leftovers) = merge_singles singles' unaligned' $
-                [ (str, second fromDecision b) | ((_,str  ),b) <- M.toList merged' ] ++
-                [ (str, second fromDecision b) | ((_,str,_),b) <- M.toList pairs' ]
-
-        (raw_pairs, raw_singles)       = partition isPaired rds
+do_rmdup :: Monad m => Bool -> Collapse m -> [BamRec] -> m [(Int,BamRec)]
+do_rmdup strand_preserved Collapse{..} rds = do
+    let (raw_pairs, raw_singles)       = partition isPaired rds
         (merged, true_singles)         = partition (liftA2 (||) isMerged isTrimmed) raw_singles
 
         (pairs, raw_half_pairs)        = partition b_totally_aligned raw_pairs
         (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 (length &&& collapse) $ accumMap f id x
-                    in (M.map (second fst) m1, concatMap (snd.snd) $ M.elems m1)
+        unaligned'                     = accumMap b_strand id half_unaligned
 
-        (pairs',r1)   = mkMap (\b -> (b_mate_pos b,   b_strand b, b_mate b)) pairs
-        (merged',r2)  = mkMap (\b -> (alignedLength (b_cigar b), b_strand b))           merged
-        (singles',r3) = mkMap                         b_strand (true_singles++half_aligned)
-        unaligned'    = accumMap b_strand id half_unaligned
+    (pairs',r1)   <- mkMap (\b -> (b_mate_pos b, b_strand b, b_mate b))    pairs
+    (merged',r2)  <- mkMap (\b -> (alignedLength (b_cigar b), b_strand b)) merged
+    (singles',r3) <- mkMap (\b ->  b_strand b)                            (true_singles++half_aligned)
 
-        b_strand b = strand_preserved && isReversed  b
-        b_mate   b = strand_preserved && isFirstMate b
+    let (results, leftovers) = merge_singles singles' unaligned' $
+                [ (str, second fromDecision b) | ((_,str  ),b) <- M.toList merged' ] ++
+                [ (str, second fromDecision b) | ((_,str,_),b) <- M.toList pairs' ]
+    pure $ results ++ map ((,) 0) (originals (leftovers ++ r1 ++ r2 ++ r3))
+  where
+    mkMap f x = do m1 <- traverse (\xs -> (,) (length xs) <$> collapse xs) $ accumMap f id x
+                   pure (M.map (second fst) m1, concatMap (snd.snd) $ M.elems m1)
 
+    b_strand b = strand_preserved && isReversed  b
+    b_mate   b = strand_preserved && isFirstMate b
+    b_mate_pos = b_mrnm &&& b_mpos &&& isUnmapped &&& isMateUnmapped
+    b_totally_aligned  = (&&) <$> not . isUnmapped <*> not . isMateUnmapped
 
+
+
+
 -- | Merging information about true singles, merged singles,
 -- half-aligned pairs, actually aligned pairs.
 --
@@ -240,8 +252,8 @@
 -- everything(?).  Then we don't have a mate for the consensus... though
 -- we could decide to duplicate one mate read to get it.
 
-merge_singles :: M.Map Bool (Int,Decision)              -- strand --> true singles & half aligned
-              -> M.Map Bool [BamRec]                    -- strand --> half unaligned
+merge_singles :: M.HashMap Bool (Int,Decision)          -- strand --> true singles & half aligned
+              -> M.HashMap Bool [BamRec]                -- strand --> half unaligned
               -> [ (Bool, (Int, BamRec)) ]              -- strand --> paireds & mergeds
               -> ([(Int,BamRec)],[BamRec])              -- results, leftovers
 
@@ -255,7 +267,7 @@
     -- strand, it goes to the leftovers.
     go ( (str, (m,v)) : paireds) =
         let (r,l) = merge_singles (M.delete str singles) (M.delete str unaligneds) paireds
-            unal  = M.findWithDefault [] str unaligneds ++ l
+            unal  = M.lookupDefault [] str unaligneds ++ l
 
         in case M.lookup str singles of
             Nothing                    -> (              (m,v) : r,     unal )
@@ -282,7 +294,7 @@
 -- singleton.  Duplication is ugly, so in this case, we force it to
 -- singleton.
 
-merge_halves :: M.Map Bool [BamRec]                     -- strand --> half unaligned
+merge_halves :: M.HashMap Bool [BamRec]                 -- strand --> half unaligned
              -> [(Bool, (Int,Decision))]                -- strand --> true singles & half aligned
              -> ([(Int,BamRec)],[BamRec])               -- results, leftovers
 
@@ -290,7 +302,8 @@
 -- we may still need it for something else to be emitted.  (While that
 -- would be strange, making sure the BAM file stays completely valid is
 -- probably better.)
-merge_halves unaligneds ((_, (n, Consensus v)) : singles) = ( (n, v { b_flag = b_flag v .&. complement pflags }) : r, l )
+merge_halves unaligneds ((_, (n, Consensus v)) : singles) =
+    ( (n, v { b_flag = b_flag v .&. complement pflags }) : r, l )
   where
     (r,l)  = merge_halves unaligneds singles
     pflags = flagPaired .|. flagProperlyPaired .|. flagMateUnmapped .|. flagMateReversed .|. flagFirstMate .|. flagSecondMate
@@ -301,32 +314,23 @@
 -- result.  Everything else goes to leftovers.  If the representative
 -- happens to be unpaired, no mate is found and that case therefore is
 -- handled smoothly.
-merge_halves unaligneds ((str, (n, Representative v)) : singles) = ((n,v) : map ((,)1) (take 1 same) ++ r, drop 1 same ++ diff ++ l)
+merge_halves unaligneds ((str, (n, Representative v)) : singles) =
+    ((n,v) : map ((,)1) (take 1 same) ++ r, drop 1 same ++ diff ++ l)
   where
     (r,l)          = merge_halves (M.delete str unaligneds) singles
-    (same,diff)    = partition (is_mate_of v) $ M.findWithDefault [] str unaligneds
+    (same,diff)    = partition (is_mate_of v) $ M.lookupDefault [] str unaligneds
     is_mate_of a b = b_qname a == b_qname b && isPaired a && isPaired b && isFirstMate a == isSecondMate b
 
 -- No more singles, all unaligneds are leftovers.
-merge_halves unaligneds [] = ( [], concat $ M.elems unaligneds )
-
-
-
-
-type MPos = (Refseq, Int, Bool, Bool)
-
-b_mate_pos :: BamRec -> MPos
-b_mate_pos b = (b_mrnm b, b_mpos b, isUnmapped b, isMateUnmapped b)
-
-b_totally_aligned :: BamRec -> Bool
-b_totally_aligned b = not (isUnmapped b || isMateUnmapped b)
+merge_halves unaligneds [] =
+    ( [], concat $ M.elems unaligneds )
 
 
-accumMap :: Ord k => (a -> k) -> (a -> v) -> [a] -> M.Map k [v]
+accumMap :: (Hashable k, Eq k) => (a -> k) -> (a -> v) -> [a] -> M.HashMap k [v]
 accumMap f g = go M.empty
   where
     go m [    ] = m
-    go m (a:as) = let ws = M.findWithDefault [] (f a) m ; g' = g a
+    go m (a:as) = let ws = M.lookupDefault [] (f a) m ; g' = g a
                   in g' `seq` go (M.insert (f a) (g':ws) m) as
 
 
@@ -359,22 +363,25 @@
    X0, X1, XT, XS, XF, XE, BC, LB, RG, XI, YI, XJ, YJ
          majority vote -}
 
-do_collapse :: Qual -> [BamRec] -> (Decision, [BamRec])
-do_collapse maxq [br] | V.all (<= maxq) (b_qual br) = ( Representative br, [  ] )     -- no modifcation, pass through
-                      | otherwise                   = ( Consensus   lq_br, [br] )     -- qualities reduced, must keep original
+do_collapse :: MonadLog m => Qual -> [BamRec] -> m (Decision, [BamRec])
+do_collapse maxq [br]
+    | maybe True (V.all (<= maxq)) (b_qual br) = pure ( Representative br, [  ] ) -- no modification, pass through
+    | otherwise                                = pure ( Consensus   lq_br, [br] ) -- qualities reduced, must keep original
   where
-    lq_br = br { b_qual  = V.map (min maxq) $ b_qual br
+    lq_br = br { b_qual  = V.map (min maxq) <$> b_qual br
                , b_virtual_offset = 0
                , b_qname = b_qname br `B.snoc` c2w 'c' }
 
-do_collapse maxq  brs = ( Consensus b0 { b_exts  = modify_extensions $ b_exts b0
-                                       , b_flag  = failflag .&. complement flagDuplicate
-                                       , b_mapq  = Q $ rmsq $ map (unQ . b_mapq) $ good brs
-                                       , b_cigar = cigar'
-                                       , b_seq   = V.fromList $ map fst cons_seq_qual
-                                       , b_qual  = V.fromList $ map snd cons_seq_qual
-                                       , b_qname = b_qname b0 `B.snoc` 99
-                                       , b_virtual_offset = 0 }, brs )              -- many modifications, must keep everything
+do_collapse maxq  brs = do
+    forM_ ws $ logMsg Warning
+    pure ( Consensus b0 { b_exts  = modify_extensions $ b_exts b0
+                        , b_flag  = failflag .&. complement flagDuplicate
+                        , b_mapq  = Q $ rmsq $ map (unQ . b_mapq) $ good brs
+                        , b_cigar = cigar'
+                        , b_seq   = V.fromList $ map fst cons_seq_qual
+                        , b_qual  = Just $ V.fromList $ map snd cons_seq_qual
+                        , b_qname = b_qname b0 `B.snoc` 99
+                        , b_virtual_offset = 0 }, brs )        -- many modifications, must keep everything
   where
     !b0 = minimumBy (comparing b_qname) brs
     !most_fail = 2 * length (filter isFailsQC brs) > length brs
@@ -392,23 +399,17 @@
     good = filter ((==) cigar' . b_cigar)
 
     cons_seq_qual = [ consensus maxq [ (V.unsafeIndex (b_seq b) i, q)
-                                     | b <- good brs, let q = if V.null (b_qual b) then Q 23 else b_qual b V.! i ]
+                                     | b <- good brs
+                                     , let q = maybe (Q 23) (V.! i) (b_qual b) ]
                     | i <- [0 .. len - 1] ]
         where !len = V.length . b_seq . head $ good brs
 
-    md' = case [ (b_seq b,md,b) | b <- good brs, Just md <- [ getMd b ] ] of
-                [               ] -> []
+    (md',ws) = case [ (b_seq b,md,b) | b <- good brs, Just md <- [ getMd b ] ] of
+                [               ] -> ([],Nothing)
                 (seq1, md1,b) : _ -> case mk_new_md' [] (V.toList cigar') md1 (V.toList seq1) (map fst cons_seq_qual) of
-                    Right x -> x
-                    Left (MdFail cigs ms osq nsq) -> error $ unlines
-                                    [ "Broken MD field when trying to construct new MD!"
-                                    , "QNAME: " ++ show (b_qname b)
-                                    , "POS:   " ++ shows (unRefseq (b_rname b)) ":" ++ show (b_pos b)
-                                    , "CIGAR: " ++ show cigs
-                                    , "MD:    " ++ show ms
-                                    , "refseq:  " ++ show osq
-                                    , "readseq: " ++ show nsq ]
-
+                    Right x -> (x,Nothing)
+                    Left (MdFail cigs ms osq nsq) -> ([],Just e)
+                        where e = BrokenMD (b_qname b) (b_rname b) (b_pos b) cigs ms osq nsq
 
     nm' = sum $ [ n | Ins :* n <- W.toList cigar' ] ++ [ n | Del :* n <- W.toList cigar' ] ++ [ 1 | MdRep _ <- md' ]
     xa' = nub' [ T.split ';' xas | Just (Text xas) <- map (lookup "XA" . b_exts) brs ]
@@ -428,9 +429,21 @@
     do_rmsq = map fromString $ words "AM AS MQ PQ SM UQ"
     do_maj  = map fromString $ words "X0 X1 XT XS XF XE BC LB RG XI XJ YI YJ"
 
-minViewBy :: (a -> a -> Ordering) -> [a] -> (a,[a])
-minViewBy  _  [    ] = error "minViewBy on empty list"
-minViewBy cmp (x:xs) = go x [] xs
+
+data BrokenMD = BrokenMD Bytes Refseq Int [Cigar] [MdOp] [Nucleotides] [Nucleotides] deriving (Typeable, Show)
+instance Exception BrokenMD where
+    displayException (BrokenMD qname rname pos cigs ms osq nsq)
+        = unlines [ "Broken MD field when trying to construct new MD!"
+                  , "QNAME: " ++ show qname
+                  , "POS:   " ++ shows (unRefseq rname) ":" ++ show pos
+                  , "CIGAR: " ++ show cigs
+                  , "MD:    " ++ show ms
+                  , "refseq:  " ++ show osq
+                  , "readseq: " ++ show nsq ]
+
+
+minViewBy :: (a -> a -> Ordering) -> a -> [a] -> (a,[a])
+minViewBy cmp x xs = go x [] xs
   where
     go m acc [    ] = (m,acc)
     go m acc (a:as) = case m `cmp` a of GT -> go a (m:acc) as
@@ -478,19 +491,19 @@
     let accs :: U.Vector Int
         accs = U.accum (+) (U.replicate 16 0) [ (fromIntegral n, fromIntegral q) | (Ns n,Q q) <- nqs ]
     in case sortBy (flip $ comparing snd) $ zip [Ns 0 ..] $ U.toList accs of
-        (n0,q0) : (_,q1) : _ ->
-            let qr = fromIntegral $ (q0-q1) `min` fromIntegral maxq
-            in if qr > 3 then (n0, Q qr) else (nucsN, Q 0)
-        _ -> error "can't happen"
+        (n0,q0) : (_,q1) : _ | qr > 3 -> (n0, Q qr)
+            where qr = fromIntegral $ (q0-q1) `min` fromIntegral maxq
+        _                             -> (nucsN, Q 0)
 
 
 -- Cheap version: simply takes the lexically first record, adds XP field
-do_cheap_collapse :: [BamRec] -> ( Decision, [BamRec] )
-do_cheap_collapse [b] = ( Representative                     b, [] )
-do_cheap_collapse  bs = ( Representative $ replaceXP new_xp b0, bx )
+do_cheap_collapse :: Monad m => [BamRec] -> m ( Decision, [BamRec] )
+do_cheap_collapse [    ] = pure ( Representative            nullBamRec, [] )
+do_cheap_collapse [  b ] = pure ( Representative                     b, [] )
+do_cheap_collapse (b:bs) = pure ( Representative $ replaceXP new_xp b0, bx )
   where
-    (b0, bx) = minViewBy (comparing b_qname) bs
-    new_xp   = foldl' (\a b -> a `oplus` extAsInt 1 "XP" b) 0 bs
+    (b0, bx) = minViewBy (comparing b_qname) b bs
+    new_xp   = foldl' (\a -> oplus a . extAsInt 1 "XP") (extAsInt 1 "XP" b) bs
 
 replaceXP :: Int -> BamRec -> BamRec
 replaceXP new_xp b0 = b0 { b_exts = updateE "XP" (Int new_xp) $ b_exts b0 }
diff --git a/Bio/Bam/Trim.hs b/Bio/Bam/Trim.hs
--- a/Bio/Bam/Trim.hs
+++ b/Bio/Bam/Trim.hs
@@ -1,28 +1,34 @@
--- | Trimming of reads as found in BAM files.  Implements trimming low
--- quality sequence from the 3' end.
+-- | Trimming and fusing of reads as found in BAM files.
+--
+-- This API is remarkably ugly because the core loop is implemented in
+-- C.  This requires the adapters to be in storable vectors, and since
+-- they shouldn't be constantly copied around, the ugly 'withADSeqs'
+-- function is needed.  The performance gain seems to be worth it,
+-- though.
 
 module Bio.Bam.Trim
         ( trim_3
         , trim_3'
         , trim_low_quality
+        , AD_Seqs
+        , withADSeqs
         , default_fwd_adapters
         , default_rev_adapters
         , find_merge
         , mergeBam
         , find_trim
         , trimBam
-        , mergeTrimBam
-        , twoMins
         , merged_seq
         , merged_qual
         ) where
 
 import Bio.Bam.Header
 import Bio.Bam.Rec
-import Bio.Bam.Rmdup        ( ECig(..), setMD, toECig )
+import Bio.Bam.Rmdup               ( ECig(..), setMD, toECig )
 import Bio.Prelude
-import Bio.Streaming
-import Foreign.C.Types      ( CInt(..) )
+import Control.Monad.Trans.Control ( MonadBaseControl, control )
+import Foreign.C.Types             ( CInt(..) )
+import Foreign.Marshal.Array       ( allocaArray )
 
 import qualified Data.ByteString                        as B
 import qualified Data.Vector.Generic                    as V
@@ -36,31 +42,33 @@
 -- Also note that trimming from the 3' end may not make sense for reads
 -- that were constructed by merging paired end data (but we cannot take
 -- care of that here).  Further note that trimming may break dependent
--- information, notably the "mate" information of the mate and many
--- optional fields.
+-- information, notably the "mate" information and many optional fields.
+-- Since the intention is to trim based on quality scores, reads without
+-- qualities are passed along unchanged.
 
 trim_3' :: ([Nucleotides] -> [Qual] -> Bool) -> BamRec -> BamRec
-trim_3' p b | b_flag b `testBit` 4 = trim_rev
-            | otherwise            = trim_fwd
-  where
-    trim_fwd = let l = subtract 1 . length . takeWhile (uncurry p) $
-                            zip (inits . reverse . V.toList $ b_seq b)
-                                (inits . reverse . V.toList $ b_qual b)
-               in trim_3 l b
+trim_3' p b = case b_qual b of
+    Nothing                         ->  b
+    Just qs | b_flag b `testBit` 4  ->  trim_3 len_rev b
+            | otherwise             ->  trim_3 len_fwd b
+      where
+        len_fwd = subtract 1 . length . takeWhile (uncurry p) $
+                      zip (inits . reverse . V.toList $ b_seq b)
+                          (inits . reverse . V.toList $ qs)
 
-    trim_rev = let l = subtract 1 . length . takeWhile (uncurry p) $
-                            zip (inits . V.toList $ b_seq  b)
-                                (inits . V.toList $ b_qual b)
-               in trim_3 l b
+        len_rev = subtract 1 . length . takeWhile (uncurry p) $
+                      zip (inits . V.toList $ b_seq  b)
+                          (inits . V.toList $ qs)
 
+
 trim_3 :: Int -> BamRec -> BamRec
 trim_3 l b | b_flag b `testBit` 4 = trim_rev
            | otherwise            = trim_fwd
   where
     trim_fwd = let (_, cigar') = trim_back_cigar (b_cigar b) l
                    c = modMd (takeECig (V.length (b_seq  b) - l)) b
-               in c { b_seq   = V.take (V.length (b_seq  c) - l) (b_seq  c)
-                    , b_qual  = V.take (V.length (b_qual c) - l) (b_qual c)
+               in c { b_seq   = V.take (V.length (b_seq c) - l)  $  b_seq  c
+                    , b_qual  = V.take (V.length (b_seq c) - l) <$> b_qual c
                     , b_cigar = cigar'
                     , b_exts  = map (\(k,e) -> case e of
                                         Text t | k `elem` trim_set
@@ -70,8 +78,8 @@
 
     trim_rev = let (off, cigar') = trim_fwd_cigar (b_cigar b) l
                    c = modMd (dropECig l) b
-               in c { b_seq   = V.drop l (b_seq  c)
-                    , b_qual  = V.drop l (b_qual c)
+               in c { b_seq   = V.drop l  $  b_seq  c
+                    , b_qual  = V.drop l <$> b_qual c
                     , b_pos   = b_pos c + off
                     , b_cigar = cigar'
                     , b_exts  = map (\(k,e) -> case e of
@@ -160,15 +168,15 @@
 -- | Finds the merge point.  Input is list of forward adapters, list of
 -- reverse adapters, sequence1, quality1, sequence2, quality2; output is
 -- merge point and two qualities (YM, YN).
-find_merge :: [W.Vector Nucleotides] -> [W.Vector Nucleotides]
+find_merge :: AD_Seqs -> AD_Seqs
            -> W.Vector Nucleotides -> W.Vector Qual
            -> W.Vector Nucleotides -> W.Vector Qual
-           -> (Int, Int, Int)
-find_merge ads1 ads2 r1 q1 r2 q2 = (mlen, score2 - score1, plain_score - score1)
-  where
-    plain_score = 6 * (V.length r1 + V.length r2)
-    (score1, mlen, score2) = twoMins plain_score (V.length r1 + V.length r2) $
-                             merge_score ads1 ads2 r1 q1 r2 q2
+           -> IO (Int, Int, Int)
+find_merge pads1 pads2 r1 q1 r2 q2 =
+        with_fw_seq r1 q1 $ \pr1 ->
+        with_fw_seq r2 q2 $ \pr2 ->
+        with_rc_seq r2 q2 $ \prv2 -> do
+            min_merge_score pads1 pads2 pr1 pr2 prv2
 
 -- | Overlap-merging of read pairs.  We shall compute the likelihood
 -- for every possible overlap, then select the most likely one (unless it
@@ -189,41 +197,45 @@
 -- would further limit the returned quality!  (In practice, map quality
 -- later imposes a limit anyway, so no worries...)
 
-mergeBam :: Int -> Int -> [W.Vector Nucleotides] -> [W.Vector Nucleotides] -> BamRec -> BamRec -> [BamRec]
-mergeBam lowq highq ads1 ads2 r1 r2
-    | V.null (b_seq r1) && V.null (b_seq r2) = [              ]
-    | qual1 < lowq || mlen < 0               = [ r1', r2'     ]
-    | qual1 >= highq && mlen == 0            = [              ]
-    | qual1 >= highq                         = [           rm ]
-    | mlen < len_r1-20 || mlen < len_r2-20   = [           rm ]
-    | otherwise         = map flag_alternative [ r1', r2', rm ]
-  where
-    len_r1    = V.length  $ b_seq  r1
-    len_r2    = V.length  $ b_seq  r2
+mergeBam :: Int -> Int
+         -> AD_Seqs -> AD_Seqs
+         -> BamRec -> BamRec -> IO [BamRec]
+mergeBam lowq highq ads1 ads2 r1 r2 = do
+    let len_r1    = V.length  $ b_seq  r1
+        len_r2    = V.length  $ b_seq  r2
 
-    b_seq_r1  = V.convert $ b_seq  r1
-    b_seq_r2  = V.convert $ b_seq  r2
-    b_qual_r1 = V.convert $ b_qual r1
-    b_qual_r2 = V.convert $ b_qual r2
+        b_seq_r1  = V.convert $ b_seq  r1
+        b_seq_r2  = V.convert $ b_seq  r2
+        b_qual_r1 = fromMaybe (V.map (const (Q 23)) b_seq_r1) (b_qual r1)
+        b_qual_r2 = fromMaybe (V.map (const (Q 23)) b_seq_r2) (b_qual r2)
 
-    (mlen, qual1, qual2) = find_merge ads1 ads2 b_seq_r1 b_qual_r1 b_seq_r2 b_qual_r2
+    (mlen, qual1, qual2) <- find_merge ads1 ads2 b_seq_r1 b_qual_r1 b_seq_r2 b_qual_r2
 
-    flag_alternative br = br { b_exts = updateE "FF" (Int $ extAsInt 0 "FF" br .|. eflagAlternative) $ b_exts br }
-    store_quals      br = br { b_exts = updateE "YM" (Int qual1) $ updateE "YN" (Int qual2) $ b_exts br }
-    pair_flags = flagPaired.|.flagProperlyPaired.|.flagMateUnmapped.|.flagMateReversed.|.flagFirstMate.|.flagSecondMate
+    let flag_alternative br = br { b_exts = updateE "FF" (Int $ extAsInt 0 "FF" br .|. eflagAlternative) $ b_exts br }
+        store_quals      br = br { b_exts = updateE "YM" (Int qual1) $ updateE "YN" (Int qual2) $ b_exts br }
+        pair_flags = flagPaired.|.flagProperlyPaired.|.flagMateUnmapped.|.flagMateReversed.|.flagFirstMate.|.flagSecondMate
 
-    r1' = store_quals r1
-    r2' = store_quals r2
-    rm  = store_quals $ merged_read mlen (fromIntegral $ min 63 qual1)
+        r1' = store_quals r1
+        r2' = store_quals r2
+        rm  = store_quals $ merged_read mlen (fromIntegral $ min 63 qual1)
 
-    merged_read l qmax = nullBamRec {
-                b_qname = b_qname r1,
-                b_flag  = flagUnmapped .|. complement pair_flags .&. b_flag r1,
-                b_seq   = V.convert $ merged_seq l b_seq_r1 b_qual_r1 b_seq_r2 b_qual_r2,
-                b_qual  = V.convert $ merged_qual qmax l b_seq_r1 b_qual_r1 b_seq_r2 b_qual_r2,
-                b_exts  = let ff = if l < len_r1 then eflagTrimmed else 0
-                          in updateE "FF" (Int $ extAsInt 0 "FF" r1 .|. eflagMerged .|. ff) $ b_exts r1 }
+        merged_read l qmax =
+            nullBamRec
+                { b_qname = b_qname r1
+                , b_flag  = flagUnmapped .|. complement pair_flags .&. b_flag r1
+                , b_seq   = V.convert $  merged_seq l b_seq_r1 b_qual_r1 b_seq_r2 b_qual_r2
+                , b_qual  = Just $ merged_qual qmax l b_seq_r1 b_qual_r1 b_seq_r2 b_qual_r2
+                , b_exts  = let ff = if l < len_r1 then eflagTrimmed else 0
+                            in updateE "FF" (Int $ extAsInt 0 "FF" r1 .|. eflagMerged .|. ff) $ b_exts r1 }
 
+    return $ case () of
+        _ | V.null (b_seq r1) && V.null (b_seq r2) -> [              ]
+          | qual1 < lowq || mlen < 0               -> [ r1', r2'     ]
+          | qual1 >= highq && mlen == 0            -> [              ]
+          | qual1 >= highq                         -> [           rm ]
+          | mlen < len_r1-20 || mlen < len_r2-20   -> [           rm ]
+          | otherwise         -> map flag_alternative [ r1', r2', rm ]
+
 {-# INLINE merged_seq #-}
 merged_seq :: (V.Vector v Nucleotides, V.Vector v Qual)
            => Int -> v Nucleotides -> v Qual -> v Nucleotides -> v Qual -> v Nucleotides
@@ -262,47 +274,44 @@
 
 -- | Finds the trimming point.  Input is list of forward adapters,
 -- sequence, quality; output is trim point and two qualities (YM, YN).
-find_trim :: [W.Vector Nucleotides]
+find_trim :: AD_Seqs
           -> W.Vector Nucleotides -> W.Vector Qual
-          -> (Int, Int, Int)
-find_trim ads1 r1 q1 = (mlen, score2 - score1, plain_score - score1)
-  where
-    plain_score = 6 * V.length r1
-    (score1, mlen, score2) = twoMins plain_score (V.length r1) $
-                             merge_score ads1 [V.empty] r1 q1 V.empty V.empty
+          -> IO (Int, Int, Int)
+find_trim pads1 r1 q1 =
+        withADSeqs [W.empty]              $ \pads2 ->
+        with_fw_seq r1 q1                   $ \pr1 ->
+        min_merge_score pads1 pads2 pr1 (FW_Seq nullPtr nullPtr 0) (RC_Seq nullPtr nullPtr 0)
 
 -- | Trimming for a single read:  we need one adapter only (the one coming
 -- /after/ the read), here provided as a list of options, and then we
 -- merge with an empty second read.  Results in up to two reads (the
 -- original, possibly flagged, and the trimmed one, definitely flagged,
 -- and two qualities).
-trimBam :: Int -> Int -> [W.Vector Nucleotides] -> BamRec -> [BamRec]
-trimBam lowq highq ads1 r1
-    | V.null (b_seq r1)              = [          ]
-    | mlen == 0 && qual1 >= highq    = [          ]
-    | qual1 < lowq || mlen < 0       = [ r1'      ]
-    | qual1 >= highq                 = [      r1t ]
-    | otherwise = map flag_alternative [ r1', r1t ]
-  where
-    -- the "merge" score if there is no trimming
-
-    b_seq_r1 = V.convert $ b_seq r1
-    b_qual_r1 = V.convert $ b_qual r1
+trimBam :: Int -> Int -> AD_Seqs -> BamRec -> IO [BamRec]
+trimBam lowq highq ads1 r1 = do
+    let b_seq_r1 = V.convert $ b_seq r1
+    (mlen, qual1, qual2) <- find_trim ads1 b_seq_r1 $
+                            fromMaybe (V.map (const (Q 23)) b_seq_r1) (b_qual r1)
 
-    (mlen, qual1, qual2) = find_trim ads1 b_seq_r1 b_qual_r1
+    let flag_alternative br = br { b_exts = updateE "FF" (Int $ extAsInt 0 "FF" br .|. eflagAlternative) $ b_exts br }
+        store_quals      br = br { b_exts = updateE "YM" (Int qual1) $ updateE "YN" (Int qual2) $ b_exts br }
 
-    flag_alternative br = br { b_exts = updateE "FF" (Int $ extAsInt 0 "FF" br .|. eflagAlternative) $ b_exts br }
-    store_quals      br = br { b_exts = updateE "YM" (Int qual1) $ updateE "YN" (Int qual2) $ b_exts br }
+        r1'  = store_quals r1
+        r1t  = store_quals $ trimmed_read mlen
 
-    r1'  = store_quals r1
-    r1t  = store_quals $ trimmed_read mlen
+        trimmed_read l = nullBamRec {
+                b_qname = b_qname r1,
+                b_flag  = flagUnmapped .|. b_flag r1,
+                b_seq   = V.take l  $  b_seq  r1,
+                b_qual  = V.take l <$> b_qual r1,
+                b_exts  = updateE "FF" (Int $ extAsInt 0 "FF" r1 .|. eflagTrimmed) $ b_exts r1 }
 
-    trimmed_read l = nullBamRec {
-            b_qname = b_qname r1,
-            b_flag  = flagUnmapped .|. b_flag r1,
-            b_seq   = V.take l $ b_seq  r1,
-            b_qual  = V.take l $ b_qual r1,
-            b_exts  = updateE "FF" (Int $ extAsInt 0 "FF" r1 .|. eflagTrimmed) $ b_exts r1 }
+    return $ case () of
+        _ | V.null (b_seq r1)              -> [          ]
+          | mlen == 0 && qual1 >= highq    -> [          ]
+          | qual1 < lowq || mlen < 0       -> [ r1'      ]
+          | qual1 >= highq                 -> [      r1t ]
+          | otherwise -> map flag_alternative [ r1', r1t ]
 
 
 -- | For merging, we don't need the complete adapters (length around 70!),
@@ -313,8 +322,8 @@
 -- (defined here in the direction they would be sequenced in):  Genomic
 -- R2, Multiplex R2, Fraft P7.
 
-default_fwd_adapters :: [ W.Vector Nucleotides ]
-default_fwd_adapters = map (W.fromList. map toNucleotides)
+default_fwd_adapters :: [W.Vector Nucleotides]
+default_fwd_adapters = map (W.fromList . map toNucleotides . map c2w)
          [ {- Genomic R2   -}  "AGATCGGAAGAGCGGTTCAG"
          , {- Multiplex R2 -}  "AGATCGGAAGAGCACACGTC"
          , {- Graft P7     -}  "AGATCGGAAGAGCTCGTATG" ]
@@ -323,8 +332,8 @@
 -- the reverse read (defined in the direction they would be sequenced in
 -- as part of the second read):  Genomic R1, CL 72.
 
-default_rev_adapters :: [ W.Vector Nucleotides ]
-default_rev_adapters = map (W.fromList. map toNucleotides)
+default_rev_adapters :: [W.Vector Nucleotides]
+default_rev_adapters = map (W.fromList . map toNucleotides . map c2w)
          [ {- Genomic_R1   -}  "AGATCGGAAGAGCGTCGTGT"
          , {- CL72         -}  "GGAAGAGCGTCGTGTAGGGA" ]
 
@@ -337,106 +346,72 @@
 -- position elsewhere.  (Yes, this ignores base composition.  It doesn't
 -- matter enough.)
 
-merge_score
-    :: [ W.Vector Nucleotides ]                 -- 3' adapters as they appear in the first read
-    -> [ W.Vector Nucleotides ]                 -- 5' adapters as they appear in the second read
-    -> W.Vector Nucleotides -> W.Vector Qual    -- first read, qual
-    -> W.Vector Nucleotides -> W.Vector Qual    -- second read, qual
-    -> Int                                      -- assumed insert length
-    -> Int                                      -- score (roughly sum of qualities at mismatches)
-merge_score fwd_adapters rev_adapters !read1 !qual1 !read2 !qual2 !l
-    =   6 * (l `min` V.length read1)                                        -- read1, part before adapter
-      + 6 * (max 0 (l - V.length read1))                                    -- read2, part before overlap
-
-      + foldl' (\acc fwd_ad -> min acc
-                    (match_adapter l read1 qual1 fwd_ad +                   -- read1, match with forward adapter
-                     6 * (max 0 (V.length read1 - V.length fwd_ad - l)))    -- read1, part after (known) adapter
-               ) maxBound fwd_adapters
-
-      + foldl' (\acc rev_ad -> min acc
-                    (match_adapter l read2 qual2 rev_ad +                   -- read2, match with reverse adapter
-                     6 * (max 0 (V.length read2 - V.length rev_ad - l)))    -- read2, part after (known) adapter
-               ) maxBound rev_adapters
-
-      + match_reads l read1 qual1 read2 qual2
-
-{-# INLINE match_adapter #-}
-match_adapter :: Int -> W.Vector Nucleotides -> W.Vector Qual -> W.Vector Nucleotides -> Int
-match_adapter !off !rd !qs !ad
-    | V.length rd /= V.length qs = error "read/qual length mismatch"
-    | efflength <= 0             = 0
-    | otherwise
-        = fromIntegral . unsafePerformIO $
-          W.unsafeWith rd $ \p_rd ->
-          W.unsafeWith qs $ \p_qs ->
-          W.unsafeWith ad $ \p_ad ->
-          prim_match_ad (fromIntegral off)
-                        (fromIntegral efflength)
-                        p_rd p_qs p_ad
-  where
-    !efflength =  (V.length rd - off) `min` V.length ad
-
-foreign import ccall unsafe "prim_match_ad"
-    prim_match_ad :: CInt -> CInt
-                  -> Ptr Nucleotides -> Ptr Qual
-                  -> Ptr Nucleotides -> IO CInt
-
-
--- | Computes overlap score for two reads (with qualities) assuming an
--- insert length.
-{-# INLINE match_reads #-}
-match_reads :: Int -> W.Vector Nucleotides -> W.Vector Qual -> W.Vector Nucleotides -> W.Vector Qual -> Int
-match_reads !l !rd1 !qs1 !rd2 !qs2
-    | V.length rd1 /= V.length qs1 || V.length rd2 /= V.length qs2 = error "read/qual length mismatch"
-    | efflength <= 0                                               = 0
-    | otherwise
-        = fromIntegral . unsafePerformIO $
-          W.unsafeWith rd1 $ \p_rd1 ->
-          W.unsafeWith qs1 $ \p_qs1 ->
-          W.unsafeWith rd2 $ \p_rd2 ->
-          W.unsafeWith qs2 $ \p_qs2 ->
-          prim_match_reads (fromIntegral minidx1)
-                           (fromIntegral maxidx2)
-                           (fromIntegral efflength)
-                           p_rd1 p_qs1 p_rd2 p_qs2
-  where
-    -- vec1, forward
-    !minidx1 = (l - V.length rd2) `max` 0
-    -- vec2, backward
-    !maxidx2 = l `min` V.length rd2
-    -- effective length
-    !efflength = ((V.length rd1 + V.length rd2 - l) `min` l) `max` 0
-
+min_merge_score
+    :: AD_Seqs              -- 3' adapters as they appear in the first read
+    -> AD_Seqs              -- 5' adapters as they appear in the second read
+    -> FW_Seq               -- first read, prepped
+    -> FW_Seq               -- second read, qual, prepped
+    -> RC_Seq               -- second read, qual, reversed and prepped
+    -> IO (Int,Int,Int)     -- best length, min score, 2nd min score
+min_merge_score (AD_Seqs !p_fwd_ads !p_fwd_lns !n_fwd_ads) (AD_Seqs !p_rev_ads !p_rev_lns !n_rev_ads)
+            (FW_Seq !p_rd1 !p_qs1 !l1) (FW_Seq !p_rd2 !p_qs2 !l2) (RC_Seq !p_rrd2 !p_rqs2 _) =
+    allocaArray 2 $ \pmins ->
+        liftM3 (,,)
+               (fromIntegral <$>
+                prim_merge_score p_fwd_ads p_fwd_lns (fromIntegral n_fwd_ads)
+                                 p_rev_ads p_rev_lns (fromIntegral n_rev_ads)
+                                 p_rd1 p_qs1 (fromIntegral l1)
+                                 p_rd2 p_qs2 (fromIntegral l2)
+                                 p_rrd2 p_rqs2 pmins)
+               (fromIntegral <$> peekElemOff pmins 0)
+               (fromIntegral <$> peekElemOff pmins 1)
 
-foreign import ccall unsafe "prim_match_reads"
-    prim_match_reads :: CInt -> CInt -> CInt
+foreign import ccall unsafe "prim_merge_score"
+    prim_merge_score :: Ptr (Ptr Nucleotides) -> Ptr CInt -> CInt
+                     -> Ptr (Ptr Nucleotides) -> Ptr CInt -> CInt
+                     -> Ptr Nucleotides -> Ptr Qual -> CInt
+                     -> Ptr Nucleotides -> Ptr Qual -> CInt
                      -> Ptr Nucleotides -> Ptr Qual
-                     -> Ptr Nucleotides -> Ptr Qual -> IO CInt
+                     -> Ptr CInt -> IO CInt
 
 
-{-# INLINE twoMins #-}
-twoMins :: (Bounded a, Ord a) => a -> Int -> (Int -> a) -> (a,Int,a)
-twoMins a0 imax f = go a0 (-1) maxBound 0 0
-  where
-    go !m1 !i1 !m2 !i2 !i
-        | i == imax = (m1,i1,m2)
-        | otherwise =
-            case f i of
-                x | x < m1    -> go  x  i m1 i1 (i+1)
-                  | x < m2    -> go m1 i1  x  i (i+1)
-                  | otherwise -> go m1 i1 m2 i2 (i+1)
 
+data AD_Seqs = AD_Seqs !(Ptr (Ptr Nucleotides)) !(Ptr CInt) !Int
+data FW_Seq  = FW_Seq !(Ptr Nucleotides) !(Ptr Qual) !Int
+data RC_Seq  = RC_Seq !(Ptr Nucleotides) !(Ptr Qual) !Int
 
-mergeTrimBam :: Monad m
-             => Int -> Int -> [W.Vector Nucleotides] -> [W.Vector Nucleotides]
-             -> Stream (Of BamRec) m r -> Stream (Of BamRec) m r
-mergeTrimBam lowq highq fwd_ads rev_ads = go
-  where
-    go = lift . inspect >=> either pure go1
+-- Maybe pad with something suitable?
+withADSeqs :: MonadBaseControl IO m => [W.Vector Nucleotides] -> (AD_Seqs -> m r) -> m r
+withADSeqs ads0 k =
+    control                                                 $ \run_io ->
+    allocaArray (length ads0)                               $ \pps ->
+    allocaArray (length ads0)                               $ \pls ->
+    let go !n [    ] = run_io (k $! AD_Seqs pps pls n)
+        go !n (v:vs) = W.unsafeWith v $ \pa -> do
+                       pokeElemOff pps n pa
+                       pokeElemOff pls n (fromIntegral (W.length v))
+                       go (succ n) vs
+    in go 0 ads0
 
-    go1 (r1 :> s) | isPaired r1 = lift (inspect s) >>= go2 r1
-                  | otherwise   = each (trimBam lowq highq fwd_ads r1) >> go s
+-- Maybe pad with something suitable?
+with_fw_seq :: W.Vector Nucleotides -> W.Vector Qual -> (FW_Seq -> IO r) -> IO r
+with_fw_seq ns qs k
+    | W.length ns == W.length qs
+        = W.unsafeWith ns    $ \p_ns ->
+          W.unsafeWith qs    $ \p_qs ->
+          k (FW_Seq p_ns p_qs $ W.length ns)
+    | otherwise
+        = throwIO $ LengthMismatch "forward adapter"
+{-# INLINE with_fw_seq #-}
 
-    go2 r1 (Left          _) = error $ "Lone mate found: " ++ show (b_qname r1)
-    go2 r1 (Right (r2 :> s)) = each (mergeBam lowq highq fwd_ads rev_ads r1 r2) >> go s
+-- Maybe pad with something suitable?
+with_rc_seq :: W.Vector Nucleotides -> W.Vector Qual -> (RC_Seq -> IO r) -> IO r
+with_rc_seq ns qs k
+    | W.length ns == W.length qs
+        = W.unsafeWith (W.reverse $ W.map compls ns)    $ \p_rns ->
+          W.unsafeWith (W.reverse qs)                   $ \p_rqs ->
+          k (RC_Seq p_rns p_rqs $ W.length ns)
+    | otherwise
+        = throwIO $ LengthMismatch "reverse adapter"
+{-# INLINE with_rc_seq #-}
 
diff --git a/Bio/Bam/Writer.hs b/Bio/Bam/Writer.hs
--- a/Bio/Bam/Writer.hs
+++ b/Bio/Bam/Writer.hs
@@ -19,7 +19,7 @@
 import Bio.Streaming.Bgzf
 
 import Data.ByteString.Builder.Prim ( (>*<) )
-import Data.ByteString.Internal     ( ByteString(..) )
+import Data.ByteString.Internal     ( fromForeignPtr )
 import Data.ByteString.Lazy         ( foldrChunks )
 import Foreign.Marshal.Alloc        ( alloca )
 
@@ -69,11 +69,11 @@
         | otherwise                            =  B.byteString (sq_name $ getRef refs mrnm)
 
     buildSeq  = E.primUnfoldrFixed E.word8 (vuncons $ \(Ns x) -> B.index "-ACMGRSVTWYHKDBN" $ fromIntegral $ x .&. 15)
-    buildQual = E.primUnfoldrFixed E.word8 (vuncons $ \(Q  q) -> q + 33)
+    buildQual = maybe (B.char7 '*') (E.primUnfoldrFixed E.word8 (vuncons $ \(Q q) -> q + 33))
 
     buildExt (BamKey k,v) = B.char7 '\t' <>
-                            B.word8 (fromIntegral k .&. 0xff) <>
-                            B.word8 (shiftR (fromIntegral k) 8 .&. 0xff) <>
+                            B.word8 (fromIntegral         k   ) <>
+                            B.word8 (fromIntegral (shiftR k 8)) <>
                             B.char7 ':' <>
                             buildExtVal v
 
@@ -182,9 +182,9 @@
     . TkWord32 (fromIntegral b_isize)
     . TkString b_qname
     . TkWord8 0
-    . W.foldr ((.) . TkWord8) id (W.unsafeCast b_cigar :: W.Vector Word8)
+    . TkMemCopy (W.unsafeCast b_cigar)
     . pushSeq b_seq
-    . W.foldr ((.) . TkWord8 . unQ) id b_qual
+    . maybe (TkMemFill (V.length b_seq) 0xff) (TkMemCopy . W.unsafeCast) b_qual
     . foldr ((.) . pushExt) id b_exts
     . TkEndRecord
   where
@@ -238,7 +238,7 @@
 packBam :: BamRec -> IO BamRaw
 packBam br = do bb <- newBuffer 1000
                 (bb', TkEnd) <- store_loop bb (pushBamRec br TkEnd)
-                return . bamRaw 0 $ PS (buffer bb') 4 (used bb' - 4)
+                bamRaw 0 $ fromForeignPtr (buffer bb') 4 (used bb' - 4)
   where
     store_loop bb tk = do (bb',tk') <- fillBuffer bb tk
                           case tk' of TkEnd -> return (bb',tk')
diff --git a/Bio/Base.hs b/Bio/Base.hs
--- a/Bio/Base.hs
+++ b/Bio/Base.hs
@@ -34,10 +34,14 @@
 
 import BasePrelude
 #if MIN_VERSION_base(4,9,0)
-                             hiding ( log1pexp, log1mexp )
+                             hiding ( log1pexp, log1mexp, (<>) )
+#else
+                             hiding ( (<>) )
 #endif
 import Bio.Util.Numeric             ( log1pexp, log1mexp )
 
+import Data.ByteString.Internal     ( c2w )
+import Data.Semigroup               ( Semigroup(..) )
 import qualified Data.ByteString.Char8 as C
 import qualified Data.Vector.Unboxed   as U
 
@@ -65,6 +69,9 @@
     minBound = Ns  0
     maxBound = Ns 15
 
+instance Semigroup Nucleotides where Ns a <> Ns b = Ns (a .|. b)
+instance Monoid Nucleotides where mempty = Ns 0 ; mappend = (<>)
+
 nucToNucs :: Nucleotide -> Nucleotides
 nucToNucs (N x) = Ns $ 1 `shiftL` fromIntegral (x .&. 3)
 
@@ -72,7 +79,7 @@
 -- represent a value @p@, we store @-10 * log_10 p@.  Operations work
 -- directly on the \"Phred\" value, as the name suggests.  The same goes
 -- for the 'Ord' instance:  greater quality means higher \"Phred\"
--- score, meand lower error probability.
+-- score, which means lower error probability.
 
 newtype Qual = Q { unQ :: Word8 } deriving ( Eq, Ord, Storable, Bounded )
 
@@ -185,43 +192,38 @@
     } deriving (Show, Eq, Ord)
 
 
--- | Converts a character into a 'Nucleotides'.
--- The usual codes for A,C,G,T and U are understood, '-' and '.' become
--- gaps and everything else is an N.
-toNucleotide :: Char -> Nucleotide
-toNucleotide c = if ord c < 128 then N (ar `U.unsafeIndex` ord c) else N 0
-  where
-    ar = U.replicate 128 0 U.//
-          ( [ (ord          x,  n) | (x, N n) <- pairs ] ++
-            [ (ord (toUpper x), n) | (x, N n) <- pairs ] )
-
-    pairs = [ ('a', nucA), ('c', nucC), ('g', nucG), ('t', nucT) ]
-
+-- | Converts a character into a 'Nucleotide'.
+-- The codes for A,C,G understood, everything else is a T.  (Error
+-- detection is the caller's responsibility.
+toNucleotide :: Word8 -> Nucleotide
+toNucleotide x | x .|. 32 == c2w 'a'  =  nucA
+               | x .|. 32 == c2w 'c'  =  nucC
+               | x .|. 32 == c2w 'g'  =  nucG
+               | otherwise            =  nucT
+{-# INLINABLE toNucleotide #-}
 
 -- | Converts a character into a 'Nucleotides'.
--- The usual codes for A,C,G,T and U are understood, '-' and '.' become
--- gaps and everything else is an N.
-toNucleotides :: Char -> Nucleotides
-toNucleotides c = if ord c < 128 then Ns (ar `U.unsafeIndex` ord c) else nucsN
-  where
-    ar = U.replicate 128 (unNs nucsN) U.//
-          ( [ (ord          x,  n) | (x, Ns n) <- pairs ] ++
-            [ (ord (toUpper x), n) | (x, Ns n) <- pairs ] )
-
-    Ns a `plus` Ns b = Ns (a .|. b)
-
-    pairs = [ ('a', nucsA), ('c', nucsC), ('g', nucsG), ('t', nucsT),
-              ('u', nucsT), ('-', gap),  ('.', gap),
-              ('b', nucsC `plus` nucsG `plus` nucsT),
-              ('d', nucsA `plus` nucsG `plus` nucsT),
-              ('h', nucsA `plus` nucsC `plus` nucsT),
-              ('v', nucsA `plus` nucsC `plus` nucsG),
-              ('k', nucsG `plus` nucsT),
-              ('m', nucsA `plus` nucsC),
-              ('s', nucsC `plus` nucsG),
-              ('w', nucsA `plus` nucsT),
-              ('r', nucsA `plus` nucsG),
-              ('y', nucsC `plus` nucsT) ]
+-- The usual codes for A,C,G,T and U are understood along with the IUPAC
+-- ambiguity codes, '-' and '.' become gaps and everything else is an N.
+toNucleotides :: Word8 -> Nucleotides
+toNucleotides x | x .|. 32 == c2w 'a'  =  nucsA
+                | x .|. 32 == c2w 'c'  =  nucsC
+                | x .|. 32 == c2w 'g'  =  nucsG
+                | x .|. 32 == c2w 't'  =  nucsT
+                | x .|. 32 == c2w 'n'  =  nucsN
+                | x .|. 32 == c2w 'u'  =  nucsT
+                | x .|. 32 == c2w 'b'  =  nucsC <> nucsG <> nucsT
+                | x .|. 32 == c2w 'd'  =  nucsA <> nucsG <> nucsT
+                | x .|. 32 == c2w 'h'  =  nucsA <> nucsC <> nucsT
+                | x .|. 32 == c2w 'v'  =  nucsA <> nucsC <> nucsG
+                | x .|. 32 == c2w 'w'  =  nucsA <> nucsT
+                | x .|. 32 == c2w 's'  =  nucsC <> nucsG
+                | x .|. 32 == c2w 'm'  =  nucsA <> nucsC
+                | x .|. 32 == c2w 'k'  =  nucsG <> nucsT
+                | x .|. 32 == c2w 'y'  =  nucsC <> nucsT
+                | x .|. 32 == c2w 'r'  =  nucsA <> nucsG
+                | otherwise            =  gap
+{-# INLINABLE toNucleotides #-}
 
 -- | Tests if a 'Nucleotides' is a base.
 -- Returns 'True' for everything but gaps.
@@ -279,10 +281,10 @@
     showList l = (map showNucleotides l ++)
 
 instance Read Nucleotides where
-    readsPrec _ (c:cs) = [(toNucleotides c, cs)]
+    readsPrec _ (c:cs) = [(toNucleotides (c2w c), cs)]
     readsPrec _ [    ] = []
     readList s = let (hd,tl) = span (\c -> isAlpha c || isSpace c || '-' == c) s
-                 in [(map toNucleotides $ filter (not . isSpace) hd, tl)]
+                 in [(map (toNucleotides . c2w) $ filter (not . isSpace) hd, tl)]
 
 -- | Complements a Nucleotides.
 {-# INLINE compl #-}
diff --git a/Bio/Prelude.hs b/Bio/Prelude.hs
--- a/Bio/Prelude.hs
+++ b/Bio/Prelude.hs
@@ -5,6 +5,7 @@
     module Bio.Util.Text,
     module Control.Monad.Catch,
     module Control.Monad.IO.Class,
+    module Control.Monad.Log,
     module Control.Monad.Trans.Class,
     module Data.Bifunctor,
     module Data.List.NonEmpty,
@@ -14,13 +15,12 @@
     Bytes, LazyBytes,
     Generic1(..),
     Hashable(..),
-    Hashable1(..),
-    Hashable2(..),
     HashMap,
     HashSet,
     IntMap,
     IntSet,
     NonEmpty(..),
+    PrimMonad(..),
     Semigroup(..),
     Text, LazyText
                    ) where
@@ -31,6 +31,7 @@
                           , Handler, handle, handleJust
                           , finally, try, tryJust, onException
                           , mask, mask_, uninterruptibleMask, uninterruptibleMask_
+                          , exitSuccess, exitFailure, exitWith
 #if MIN_VERSION_base(4,9,0)
                           , log1p, log1pexp, log1mexp, expm1
 #endif
@@ -40,6 +41,8 @@
 import Bio.Util.Text
 import Control.Monad.Catch
 import Control.Monad.IO.Class
+import Control.Monad.Log            ( LIO, execWithParser_, MonadLog(..), Level(..) )
+import Control.Monad.Primitive      ( PrimMonad(..) )
 import Control.Monad.Trans.Class
 import Data.Bifunctor
 import Data.ByteString              ( ByteString )
@@ -47,13 +50,12 @@
 import Data.Semigroup               ( Semigroup(..) )
 import Data.Text                    ( Text )
 import Data.Hashable                ( Hashable(..) )
-import Data.Hashable.Lifted         ( Hashable1(..), Hashable2(..) )
 import Data.HashMap.Strict          ( HashMap )
 import Data.HashSet                 ( HashSet )
 import Data.IntMap                  ( IntMap )
 import Data.IntSet                  ( IntSet )
 import GHC.Generics                 ( Generic1(..) )
-import System.IO                    ( hPrint, hPutStr, hPutStrLn, stderr, stdout, stdin
+import System.IO                    ( stdin, stdout, stderr
                                     , openBinaryFile, withBinaryFile, IOMode(..)
                                     , hFlush, hSeek, hClose, SeekMode(..) )
 
diff --git a/Bio/Streaming.hs b/Bio/Streaming.hs
--- a/Bio/Streaming.hs
+++ b/Bio/Streaming.hs
@@ -1,3 +1,9 @@
+{-# LANGUAGE CPP, TypeFamilies #-}
+#if __GLASGOW_HASKELL__ >= 800
+{-# OPTIONS_GHC -Wno-orphans #-}
+#else
+{-# OPTIONS_GHC -fno-warn-orphans #-}
+#endif
 module Bio.Streaming
     ( MonadIO(..)
     , MonadMask
@@ -9,6 +15,7 @@
     , streamInputs
     , withOutputFile
 
+    , UnwantedTerminal(..)
     , protectTerm
     , psequence
     , progressGen
@@ -30,10 +37,14 @@
 import Streaming                     hiding ( (<>) )
 import Streaming.Internal                   ( Stream(..) )
 import Streaming.Prelude                    ( each )
-import System.IO
+import System.IO                            ( hIsTerminalDevice )
 
 import qualified Streaming.Prelude      as Q
 
+instance (Functor f, PrimMonad m) => PrimMonad (Stream f m) where
+    type PrimState (Stream f m) = PrimState m
+    primitive = lift . primitive
+
 {- | Default buffer size in elements.
 
 Since we often want to merge many files, a read should take more time
@@ -74,6 +85,10 @@
                          liftIO $ hClose h
 {-# INLINE streamInputs #-}
 
+data UnwantedTerminal = UnwantedTerminal deriving (Typeable, Show)
+instance Exception UnwantedTerminal where
+    displayException _ = "cowardly refusing to write binary data to terminal"
+
 {- | Protects the terminal from binary junk.
 
 If @s@ is a 'Stream', then @protectTerm s@ throws an error if 'stdout'
@@ -84,7 +99,7 @@
 protectTerm :: (Functor f, MonadIO m) => Stream f m r -> Stream f m r
 protectTerm str = do
     t <- liftIO $ hIsTerminalDevice stdout
-    when t . liftIO . throwM $ ErrorCall "cowardly refusing to write binary data to terminal"
+    when t . liftIO . throwM $ UnwantedTerminal
     str
 {-# INLINE protectTerm #-}
 
@@ -167,27 +182,24 @@
           GT -> Step (b :> go str0 rest1)
 {-# INLINABLE mergeStreamsBy #-}
 
--- | A general progress indicator that prints some message after a set
+-- | A general progress indicator that logs some message after a set
 -- number of records have passed through.
-progressGen :: MonadIO m
-            => (Int -> a -> String) -> Int -> (String -> IO ())
-            -> Q.Stream (Q.Of a) m r -> Q.Stream (Q.Of a) m r
-progressGen msg sz put = go 0
+progressGen :: MonadLog m => (Int -> a -> String) -> Int -> Q.Stream (Q.Of a) m r -> Q.Stream (Q.Of a) m r
+progressGen msg sz = go 0
   where
-    go   !n = lift . Q.next >=> either pure (step $ succ n)
-    step !n (a,s) = do when (n `mod` sz == 0) . liftIO . put $ "\27[K" ++ msg n a ++ "\r"
+    go   !n = lift . Q.next >=> either fin (step $ succ n)
+    step !n (a,s) = do when (n `mod` sz == 0) . lift . logString_ $ msg n a
                        Q.cons a (go n s)
+    fin r = r <$ lift (logString_ "")
 
--- | A simple progress indicator that prints the number of records.
-progressNum :: MonadIO m
-            => String -> Int -> (String -> IO ())
-            -> Q.Stream (Q.Of a) m r -> Q.Stream (Q.Of a) m r
+-- | A simple progress indicator that logs the number of records.
+progressNum :: MonadLog m => String -> Int -> Q.Stream (Q.Of a) m r -> Q.Stream (Q.Of a) m r
 progressNum msg = progressGen (\n _ -> msg ++ " " ++ showNum n)
 
--- | A simple progress indicator that prints a position every set number
+-- | A simple progress indicator that logs a position every set number
 -- of passed records.
-progressPos :: MonadIO m
-            => (a -> (Refseq, Int)) -> String -> Refs -> Int -> (String -> IO ())
+progressPos :: MonadLog m
+            => (a -> (Refseq, Int)) -> String -> Refs -> Int
             -> Q.Stream (Q.Of a) m r -> Q.Stream (Q.Of a) m r
 progressPos f msg refs =
     progressGen $ \_ a -> let (!rs1, !po1) = f a
diff --git a/Bio/Streaming/Bgzf.hs b/Bio/Streaming/Bgzf.hs
--- a/Bio/Streaming/Bgzf.hs
+++ b/Bio/Streaming/Bgzf.hs
@@ -16,13 +16,15 @@
     BclArgs(..),
     BclSpecialType(..),
     loop_dec_int,
-    loop_bcl_special
+    loop_bcl_special,
+    CompressionError(..),
+    DecompressionError(..)
                             ) where
 
 import Bio.Prelude
 import Bio.Streaming
 import Foreign.C.Types                     ( CInt(..) )
-import Foreign.Marshal.Utils               ( copyBytes, with )
+import Foreign.Marshal.Utils               ( copyBytes, fillBytes, with )
 
 import qualified Bio.Streaming.Bytes        as S
 import qualified Data.ByteString            as B
@@ -89,6 +91,9 @@
                 | TkString   {-# UNPACK #-} !B.ByteString BgzfTokens -- a raw string
                 | TkDecimal  {-# UNPACK #-} !Int          BgzfTokens -- roughly ':%d'
 
+                | TkMemFill {-# UNPACK #-} !Int {-# UNPACK #-} !Word8   BgzfTokens
+                | TkMemCopy {-# UNPACK #-} !(V.Vector Word8)            BgzfTokens
+
                 | TkSetMark                               BgzfTokens -- sets the first mark
                 | TkEndRecord                             BgzfTokens -- completes a BAM record
                 | TkEndRecordPart1                        BgzfTokens -- completes part 1 of a BCF record
@@ -129,13 +134,21 @@
              , mark   = if mark  b == maxBound then maxBound else mark  b - off b
              , mark2  = if mark2 b == maxBound then maxBound else mark2 b - off b }
 
+data CompressionError = CompressionError !CInt deriving (Typeable,Show)
+instance Exception CompressionError where
+    displayException (CompressionError rc) = "compress_chunk failed: " ++ show rc
+
+data DecompressionError = DecompressionError !CInt deriving (Typeable,Show)
+instance Exception DecompressionError where
+    displayException (DecompressionError rc) = "decompress_chunk failed: " ++ show rc
+
 compressChunk :: Int -> ForeignPtr Word8 -> Int -> Int -> IO B.ByteString
 compressChunk lv fptr off slen =
     withForeignPtr fptr                 $ \ptr ->
     B.createAndTrim 65536               $ \buf ->
     with 65536                          $ \p_len -> do
         rc <- compress_chunk buf p_len (plusPtr ptr off) (fromIntegral slen) (fromIntegral lv)
-        when (rc /= 0 && rc /= 1) . error $ "compress_chunk failed: " ++ show rc
+        when (rc /= 0 && rc /= 1) $ throwIO $ CompressionError rc
         fromIntegral <$> peek p_len
 
 decompressChunk :: B.ByteString -> IO B.ByteString
@@ -144,7 +157,7 @@
     peekByteOff psrc (B.length ck - 4)            >>= \dlen ->
     B.create (fromIntegral (dlen::Word32))          $ \pdest -> do
         rc <- decompress_chunk pdest (fromIntegral dlen) (castPtr psrc) (fromIntegral $ B.length ck)
-        when (rc /= 0) . error $ "decompress_chunk failed: " ++ show rc
+        when (rc /= 0) $ throwIO $ DecompressionError rc
 
 
 -- | Expand a chain of tokens into a buffer, sending finished pieces
@@ -220,10 +233,10 @@
         TkDouble   x tk' -> do pokeByteOff p use x
                                go_fast p bb (use + 8) tk'
 
+        -- The next three may be too big to handle.  By returning with
+        -- unfinished business, we will get progressively bigger buffers
+        -- and eventually handle it just fine.
         TkString   s tk'
-            -- Too big, can't handle.  By returning with unfinished
-            -- business, we will get progressively bigger buffers and
-            -- eventually handle it.
             | B.length s > size bb - use -> return (bb { used = use },tk1)
 
             | otherwise  -> do let ln = B.length s
@@ -231,6 +244,21 @@
                                     copyBytes (p `plusPtr` use) q ln
                                go_fast p bb (use + ln) tk'
 
+        TkMemFill ln c tk'
+            | ln > size bb - use -> return (bb { used = use },tk1)
+
+            | otherwise  -> do fillBytes (p `plusPtr` use) c ln
+                               go_fast p bb (use + ln) tk'
+
+        TkMemCopy v tk'
+            | V.length v > size bb - use -> return (bb { used = use },tk1)
+
+            | otherwise  -> do let ln = V.length v
+                               V.unsafeWith v $ \q ->
+                                    copyBytes (p `plusPtr` use) q ln
+                               go_fast p bb (use + ln) tk'
+
+
         TkDecimal  x tk' -> do ln <- int_loop (p `plusPtr` use) (fromIntegral x)
                                go_fast p bb (use + fromIntegral ln) tk'
 
@@ -247,7 +275,6 @@
         TkEndRecordPart2 tk' -> do let !l = use - mark2 bb
                                    pokeByteOff p (mark bb) (fromIntegral l :: Word32)
                                    go_slowish p bb { used = use, mark = maxBound } tk'
-
 
         TkBclSpecial special_args tk' -> do
             l <- loop_bcl_special (p `plusPtr` use) special_args
diff --git a/Bio/Streaming/Bytes.hs b/Bio/Streaming/Bytes.hs
--- a/Bio/Streaming/Bytes.hs
+++ b/Bio/Streaming/Bytes.hs
@@ -124,7 +124,7 @@
         ,AllocationStrategy,ChunkIOStream(..),buildStepToCIOS
         ,byteStringFromBuffer,safeStrategy,defaultChunkSize)
 import GHC.Exts                                (SpecConstrAnnotation(..))
-import Streaming                               (Of(..),Identity(..),destroy)
+import Streaming                               (MFunctor(..),Of(..),Identity(..),destroy)
 import Streaming.Internal                      (Stream (..))
 import System.Directory                        (renameFile)
 
@@ -143,8 +143,8 @@
 -- also contains an offset, which will be needed to track the virtual
 -- offsets in the BGZF decode.
 
-data ByteStream m r =
-  Empty r
+data ByteStream m r
+  = Empty r
   | Chunk {-# UNPACK #-} !Bytes {-# UNPACK #-} !Int64 (ByteStream m r)
   | Go (m (ByteStream m r))
 
@@ -187,6 +187,13 @@
   lift ma = Go $ liftM Empty ma
   {-# INLINE lift #-}
 
+instance MFunctor ByteStream where
+  hoist f = loop where
+    loop (Empty     r) = Empty r
+    loop (Chunk c o s) = Chunk c o (loop s)
+    loop (Go        m) = Go (f (liftM loop m))
+  {-# INLINEABLE hoist #-}
+
 instance (r ~ ()) => IsString (ByteStream m r) where
   fromString = chunk . fromString
   {-# INLINE fromString #-}
@@ -682,8 +689,7 @@
 
 {- | Turns a 'ByteStream' into a stream of strict 'Bytes' that divide at
      newline characters. The resulting strings do not contain newlines.
-     This will cost memory if the lines are very long, and it does not
-     recognize DOS line endings. -}
+     This will cost memory if the lines are very long.  -}
 
 lines' :: Monad m => ByteStream m r -> Stream (Of Bytes) m r
 lines' = loop1 []
@@ -691,15 +697,19 @@
     loop1 :: Monad m => [Bytes] -> ByteStream m r -> Stream (Of Bytes) m r
     loop1 acc text =
       case text of
-        Empty r -> Return r
+        Empty r -> r <$ unless (null acc) (Q.yield (checkCR $ B.concat (reverse acc)))
         Go m    -> Effect $ liftM (loop1 acc) m
         Chunk c o cs
           | B.null c  -> loop1 acc cs
           | otherwise ->
               case B.elemIndex 10 c of
-                Just  i -> Q.cons (if null acc then B.take i c else B.concat (reverse (B.take i c : acc)))
+                Just  i -> Q.cons (checkCR $ if null acc then B.take i c else B.concat (reverse (B.take i c : acc)))
                                   (loop1 [] (Chunk (B.drop (i+1) c) (o+1 + fromIntegral i) cs))
                 Nothing -> loop1 (c:acc) cs
+    checkCR s
+        | B.null s        =  s
+        | B.last s == 13  =  B.init s
+        | otherwise       =  s
 {-# INLINABLE lines' #-}
 
 -- --------------------------------------------------------------------------
diff --git a/Bio/Streaming/Parse.hs b/Bio/Streaming/Parse.hs
--- a/Bio/Streaming/Parse.hs
+++ b/Bio/Streaming/Parse.hs
@@ -1,10 +1,12 @@
 {-# LANGUAGE Rank2Types #-}
+-- | Parsers for use with 'ByteStream's.
 module Bio.Streaming.Parse
     ( Parser
     , ParseError(..)
     , EofException(..)
     , parse
     , parseIO
+    , parseLog
     , parseM
     , abortParse
     , isFinished
@@ -18,8 +20,6 @@
     , atto
     ) where
 
--- ^ Parsers for use with 'ByteStream's.
-
 import Bio.Prelude                       hiding ( drop )
 import Bio.Streaming.Bytes                      ( ByteStream )
 
@@ -30,10 +30,10 @@
 
 newtype Parser r m a = P {
     runP :: forall x .
-            (a -> ByteStream m r -> m x)
-         -> (r -> m x)
-         -> (SomeException -> m x)
-         -> ByteStream m r -> m x }
+            (a -> ByteStream m r -> m x)                -- successful parse
+         -> (r -> m x)                                  -- end of input stream
+         -> (SomeException -> ByteStream m r -> m x)    -- exception and remaining input
+         -> ByteStream m r -> m x }                     -- input, result
 
 instance Functor (Parser r m) where
     fmap f p = P $ \sk -> runP p (sk . f)
@@ -53,12 +53,13 @@
     lift m = P $ \sk _rk _ek s -> m >>= \a -> sk a s
 
 instance MonadThrow (Parser r m) where
-    throwM e = P $ \_sk _rk ek _s -> ek (toException e)
+    throwM e = P $ \_sk _rk ek -> ek (toException e)
 
 modify :: (ByteStream m r -> ByteStream m r) -> Parser r m ()
 modify f = P $ \sk _rk _ek -> sk () . f
 
-parse :: Monad m => (Int64 -> Parser r m a) -> ByteStream m r -> m (Either SomeException (Either r (a, ByteStream m r)))
+parse :: Monad m => (Int64 -> Parser r m a) -> ByteStream m r
+      -> m (Either (SomeException, ByteStream m r) (Either r (a, ByteStream m r)))
 parse p = go
   where
     go    (S.Empty     r)             = return $ Right $ Left r
@@ -66,14 +67,18 @@
     go ck@(S.Chunk c o s) | B.null  c = go s
                           | otherwise = runP (p o) (\a t -> return . Right $ Right (a,t))
                                                    (return . Right . Left)
-                                                   (return . Left)
+                                                   (curry $ return . Left)
                                                    ck
 
 parseIO :: MonadIO m => (Int64 -> Parser r m a) -> ByteStream m r -> m (Either r (a, ByteStream m r))
-parseIO p = parse p >=> either (liftIO . throwM) return
+parseIO p = parse p >=> either (liftIO . throwM . fst) return
 
+parseLog :: MonadLog m => Level -> (Int64 -> Parser r m a) -> ByteStream m r -> m (Either r (a, ByteStream m r))
+parseLog lv p = parse p >=> either throw_it pure
+  where throw_it (ex,rest) = logMsg lv ex >> Left <$> S.effects rest
+
 parseM :: MonadThrow m => (Int64 -> Parser r m a) -> ByteStream m r -> m (Either r (a, ByteStream m r))
-parseM p = parse p >=> either throwM return
+parseM p = parse p >=> either (throwM . fst) return
 
 abortParse :: Monad m => Parser r m a
 abortParse = P $ \_sk rk _ek -> S.effects >=> rk
@@ -96,7 +101,7 @@
 dropLine = modify $ S.drop 1 . S.dropWhile (/= 10)
 
 getByte :: Monad m => Parser r m Word8
-getByte = P $ \sk _rk ek -> S.nextByte >=> either (const $ ek (toException EofException)) (uncurry sk)
+getByte = P $ \sk _rk ek -> S.nextByte >=> either (ek (toException EofException) . pure) (uncurry sk)
 
 getString :: Monad m => Int -> Parser r m B.ByteString
 getString l = liftFun $ liftM Q.lazily . S.splitAt' l
@@ -110,16 +115,17 @@
 isolate :: Monad m => Int -> Parser (ByteStream m r) m a -> Parser r m a
 isolate l p = P $ \sk rk ek -> runP p (\a -> S.effects >=> sk a)
                                       (S.effects >=> rk)
-                                      ek . S.splitAt (fromIntegral l)
+                                      (\e rest -> ek e (join rest)) .
+                               S.splitAt (fromIntegral l)
 
-data ParseError = ParseError {errorContexts :: [String], errorMessage :: String}
-    deriving (Show, Typeable)
 
-data EofException = EofException
-    deriving (Show, Typeable)
+data EofException = EofException deriving (Show, Typeable)
+instance Exception EofException where displayException _ = "end-of-file"
 
-instance Exception ParseError
-instance Exception EofException
+data ParseError = ParseError {errorContexts :: [String], errorMessage :: String} deriving (Show, Typeable)
+instance Exception ParseError where
+    displayException (ParseError ctx msg)
+        = "Parse error at " ++ intercalate ", " ctx ++ ": " ++ msg
 
 atto :: Monad m => A.Parser a -> Parser r m a
 atto = go . A.parse
@@ -127,13 +133,13 @@
     go k = P $ \sk rk ek ->
         S.nextChunk >=> \case
             Left r -> case k B.empty of
-                      A.Fail _ err dsc -> ek $ toException (ParseError err dsc)
-                      A.Partial _      -> ek $ toException EofException
+                      A.Fail _ err dsc -> ek (toException (ParseError err dsc)) (pure r)
+                      A.Partial _      -> ek (toException EofException) (pure r)
                       A.Done rest v    -> sk v (S.consChunk rest (pure r))
             Right (c,s')
                 | B.null c -> runP (go k) sk rk ek s'
                 | otherwise -> case k c of
-                      A.Fail _ err dsc -> ek $ toException (ParseError err dsc)
+                      A.Fail _ err dsc -> ek (toException (ParseError err dsc)) s'
                       A.Partial k'     -> runP (go k') sk rk ek s'
                       A.Done rest v    -> sk v (S.consChunk rest s')
 
diff --git a/Bio/TwoBit.hs b/Bio/TwoBit.hs
--- a/Bio/TwoBit.hs
+++ b/Bio/TwoBit.hs
@@ -33,19 +33,19 @@
         Mask(..)
     ) where
 
-import           Bio.Prelude hiding ( left, right, chr )
+import           Bio.Prelude         hiding ( left, right, chr )
 import           Bio.Util.MMap
 import           Bio.Util.Storable
-import           Control.Monad.Trans.State ( StateT(..), get, evalStateT )
+import           Control.Monad.Trans.State  ( StateT(..), get, evalStateT )
 import qualified Data.ByteString                as B
-import qualified Data.ByteString.Unsafe         as B
 import qualified Data.IntMap.Strict             as I
 import qualified Data.HashMap.Lazy              as M
 import qualified Data.Vector.Unboxed            as U
-import           Foreign.C.Types ( CChar )
+import           Foreign.C.Types            ( CChar )
+import           Foreign.ForeignPtr.Unsafe  ( unsafeForeignPtrToPtr )
 
 data TwoBitFile = TBF {
-    tbf_raw :: B.ByteString,
+    tbf_raw :: ForeignPtr CChar,
     -- This map is intentionally lazy.  May or may not be important.
     tbf_seqs :: !(M.HashMap Bytes TwoBitSequence)
 }
@@ -60,10 +60,9 @@
 -- the file is modified in any way.
 openTwoBit :: FilePath -> IO TwoBitFile
 openTwoBit fp = do
-        raw <- unsafeMMapFile fp
-        B.unsafeUseAsCString raw $ \praw ->
-        -- return $ flip runGet (L.fromChunks [raw]) $ do
-            flip evalStateT praw $ do
+        (_sz, raw) <- mmapFile fp
+        withForeignPtr raw $
+              evalStateT $ do
                     sig <- getWord32be
                     getWord32 <- case sig :: Word32 of
                             0x1A412743 -> return getWord32be
@@ -91,10 +90,10 @@
 getByteString :: Int -> Get Bytes
 getByteString l = StateT $ \p -> B.packCStringLen (p,l) >>= \s -> return (s, plusPtr p l)
 
-mkBlockIndex :: B.ByteString -> Get Int -> Int -> TwoBitSequence
+mkBlockIndex :: ForeignPtr CChar -> Get Int -> Int -> TwoBitSequence
 mkBlockIndex raw getWord32 ofs =
-    unsafePerformIO $
-    B.unsafeUseAsCString raw $ \praw ->
+    unsafeDupablePerformIO $
+    withForeignPtr raw $ \praw ->
     evalStateT getBlock (plusPtr praw ofs)
   where
     getBlock = do p0 <- get
@@ -132,10 +131,10 @@
 getFwdSubseqWith TBF{..} TBS{..} nt start =
     do_mask (takeOverlap start tbs_n_blocks `mergeBlocks` takeOverlap start tbs_m_blocks) start .
     drop (start .&. 3) .
-    B.foldr toDNA [] .
-    B.drop (tbs_dna_offset + (start `shiftR` 2)) $ tbf_raw
+    toDNA $ plusPtr (unsafeForeignPtrToPtr tbf_raw) (tbs_dna_offset + (start `shiftR` 2))
   where
-    toDNA b = (++) [ 3 .&. (b `shiftR` x) | x <- [6,4,2,0] ]
+    toDNA p = let !b = unsafeDupablePerformIO $ do peek p <* touchForeignPtr tbf_raw
+              in [ 3 .&. (b `shiftR` x) | x <- [6,4,2,0] ] ++ toDNA (plusPtr p 1)
 
     do_mask            _ _ [] = []
     do_mask [          ] _ ws = map (`nt` None) ws
@@ -160,10 +159,10 @@
 mergeBlocks [     ] [     ] = []
 
 
--- | Extract a subsequence and apply masking.  TwoBit file can represent
--- two kinds of masking (hard and soft), where hard masking is usually
--- realized by replacing everything by Ns and soft masking is done by
--- lowercasing.  Here, we take a user supplied function to apply
+-- | Extract a subsequence and apply masking.  TwoBit files can
+-- represent two kinds of masking (hard and soft), where hard masking is
+-- usually realized by replacing everything by Ns and soft masking is
+-- done by 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
@@ -175,6 +174,7 @@
   where
     fwd_nt = (!!) [nucT, nucC, nucA, nucG] . fromIntegral
     cmp_nt = (!!) [nucA, nucG, nucT, nucC] . fromIntegral
+{-# INLINE getSubseqWith #-}
 
 -- | Works only in forward direction.
 getLazySubseq :: TwoBitFile -> Position -> [Nucleotide]
@@ -183,15 +183,16 @@
     let go  = getFwdSubseqWith tbf sq1
     if start < 0
         then error "sorry, can't go backwards"
-        -- then reverse $ take len $ go (maskf . cmp_nt) (-start-len)
         else go fwd_nt start
   where
     fwd_nt n _ = [nucT, nucC, nucA, nucG] !! fromIntegral n
+{-# INLINE getLazySubseq #-}
 
 
 -- | Extract a subsequence without masking.
 getSubseq :: TwoBitFile -> Range -> [Nucleotide]
 getSubseq = getSubseqWith const
+{-# INLINE getSubseq #-}
 
 -- | Extract a subsequence with typical masking:  soft masking is
 -- ignored, hard masked regions are replaced with Ns.
@@ -202,6 +203,7 @@
     mymask n Soft = nucToNucs n
     mymask _ Hard = nucsN
     mymask _ Both = nucsN
+{-# INLINE getSubseqMasked #-}
 
 -- | Extract a subsequence with masking for biologists:  soft masking is
 -- done by lowercasing, hard masking by printing an N.
@@ -212,6 +214,7 @@
     mymask n Soft = toLower (showNucleotide n)
     mymask _ Hard = 'N'
     mymask _ Both = 'N'
+{-# INLINE getSubseqAscii #-}
 
 
 getSeqnames :: TwoBitFile -> [Bytes]
@@ -291,6 +294,8 @@
                        | off < s+l -> Just (4, (succ off, nbs))
                        | otherwise -> Just (y, (succ off, nbs'))
       where
-        x = B.index tbf_raw (tbs_dna_offset + off `shiftR` 2)
+        x = unsafeDupablePerformIO $
+            withForeignPtr tbf_raw $ \p ->
+            peekByteOff p (tbs_dna_offset + off `shiftR` 2)
         y = x `shiftR` (6 - 2 * (off .&. 3)) .&. 3     -- T,C,A,G
 
diff --git a/Bio/Util/MMap.hs b/Bio/Util/MMap.hs
deleted file mode 100644
--- a/Bio/Util/MMap.hs
+++ /dev/null
@@ -1,28 +0,0 @@
-{-# LANGUAGE ForeignFunctionInterface #-}
-module Bio.Util.MMap ( unsafeMMapFile ) where
-
-import BasePrelude
-import Data.ByteString.Internal ( fromForeignPtr, ByteString )
-import Foreign.C.Types
-import System.Posix.Files
-import System.Posix.IO
-
-unsafeMMapFile :: FilePath -> IO ByteString
-unsafeMMapFile fp =
-    bracket (openFd fp ReadOnly Nothing defaultFileFlags) closeFd $ \fd -> do
-        stat <- getFdStatus fd
-        let size = fromIntegral (fileSize stat)
-        if size <= 0
-            then return mempty
-            else do
-                ptr <- c_mmap size (fromIntegral fd)
-                if ptr == nullPtr
-                    then error "unable to mmap file"
-                    else do
-                          fptr <- newForeignPtrEnv c_munmap (intPtrToPtr $ fromIntegral size) ptr
-                          return $ fromForeignPtr fptr 0 (fromIntegral size)
-
-foreign import ccall unsafe  "my_mmap"   c_mmap   :: CSize -> CInt -> IO (Ptr Word8)
-foreign import ccall unsafe "&my_munmap" c_munmap :: FunPtr (Ptr () -> Ptr Word8 -> IO ())
-
-
diff --git a/Bio/Util/MMap.hsc b/Bio/Util/MMap.hsc
new file mode 100644
--- /dev/null
+++ b/Bio/Util/MMap.hsc
@@ -0,0 +1,43 @@
+module Bio.Util.MMap ( mmapFile, createMmapFile ) where
+
+import BasePrelude
+import Foreign.C.Error      ( getErrno, errnoToIOError )
+import Foreign.C.Types
+import System.Posix.Files   ( fileSize, getFdStatus, setFdSize )
+import System.Posix.IO      ( openFd, closeFd, defaultFileFlags, OpenMode(ReadOnly,ReadWrite) )
+import System.Posix.Types   ( Fd(..), COff(..) )
+
+#include <sys/mman.h>
+
+-- | Maps a whole file into memory, returns the size in bytes and a
+-- 'ForeignPtr' to the contents.
+mmapFile :: FilePath -> IO (Int, ForeignPtr a)
+mmapFile fp =
+    bracket (openFd fp ReadOnly Nothing defaultFileFlags) closeFd $ \fd -> do
+        size <- fileSize <$> getFdStatus fd
+        if size <= 0
+            then (,) 0 <$> newForeignPtr_ nullPtr
+            else do
+                ptr <- mmap nullPtr (fromIntegral size) (#const PROT_READ) (#const MAP_SHARED) fd 0
+                if ptrToIntPtr ptr == #const MAP_FAILED
+                    then do errno <- getErrno
+                            ioError $ errnoToIOError "mmapFile" errno Nothing (Just fp)
+                    else (,) (fromIntegral size) <$> newForeignPtrEnv my_munmap (intPtrToPtr $ fromIntegral size) ptr
+
+-- | Creates a new file of a desired initial size, maps it into memory,
+-- and calls a function to fill it.  That function returns a pointer to
+-- the first unused byte in the file, and it is truncated accordingly.
+createMmapFile :: FilePath -> CSize -> (Ptr a -> IO (Ptr a, b)) -> IO b
+createMmapFile fp sz k =
+    bracket (openFd fp ReadWrite (Just 0x1b6) defaultFileFlags) closeFd $ \fd -> do
+        setFdSize fd (fromIntegral sz)
+        bracket (mmap nullPtr sz (#const PROT_READ | PROT_WRITE) (#const MAP_SHARED) fd 0)
+                (flip munmap sz) $ \p -> do
+            (p',r) <- k p
+            setFdSize fd (fromIntegral $ minusPtr p' p)
+            return r
+
+foreign import ccall unsafe "&my_munmap"        my_munmap :: FunPtr (Ptr () -> Ptr a -> IO ())
+foreign import ccall unsafe "sys/mman.h mmap"   mmap      :: Ptr a -> CSize -> CInt -> CInt -> Fd -> COff -> IO (Ptr a)
+foreign import ccall unsafe "sys/mman.h munmap" munmap    :: Ptr a -> CSize -> IO ()
+
diff --git a/Bio/Util/Numeric.hs b/Bio/Util/Numeric.hs
--- a/Bio/Util/Numeric.hs
+++ b/Bio/Util/Numeric.hs
@@ -1,8 +1,7 @@
 -- | Random useful stuff I didn't know where to put.
-
 module Bio.Util.Numeric (
     wilson, invnormcdf, choose,
-    estimateComplexity, showNum, showOOM,
+    estimateComplexity, showNum, showOOM, readOOM,
     log1p, expm1, (<#>),
     log1mexp, log1pexp,
     lsum, llerp
@@ -11,6 +10,7 @@
 import Prelude
 import Data.Char ( intToDigit )
 import Data.List ( foldl1' )
+import Options.Applicative.Builder ( eitherReader, ReadM )
 
 -- | Calculates the Wilson Score interval.
 -- If @(l,m,h) = wilson c x n@, then @m@ is the binary proportion and
@@ -37,18 +37,29 @@
     triplets acc [a,b,c] = c:b:a:acc
     triplets acc (a:b:c:s) = triplets (',':c:b:a:acc) s
 
-showOOM :: Double -> String
-showOOM x | x < 0 = '-' : showOOM (negate x)
-          | otherwise = findSuffix (x*10) ".kMGTPEZY"
+
+showOOM :: (Enum a, Num a, Ord a) => a -> String
+showOOM x | x < 0     = '-' : showOOM (negate x)
+          | otherwise = findSuffix (fromEnum (x * 100 + 5) `div` 10) ".kMGTPEZY"
   where
-    findSuffix _ [] = "many"
-    findSuffix y (s:ss) | y < 100  = intToDigit (round y `div` 10) : case (round y `mod` 10, s) of
+    findSuffix _ [    ] = "many"
+    findSuffix y (s:ss) | y < 100  = intToDigit (div y 10) : case (mod y 10, s) of
                                             (0,'.') -> [] ; (0,_) -> [s] ; (d,_) -> [s, intToDigit d]
-                        | y < 1000 = intToDigit (round y `div` 100) : intToDigit ((round y `mod` 100) `div` 10) :
+                        | y < 1000 = intToDigit (div y 100) : intToDigit (mod y 100 `div` 10) :
                                             if s == '.' then [] else [s]
-                        | y < 10000 = intToDigit (round y `div` 1000) : intToDigit ((round y `mod` 1000) `div` 100) :
+                        | y < 10000 = intToDigit (div y 1000) : intToDigit (mod y 1000 `div` 100) :
                                             '0' : if s == '.' then [] else [s]
-                        | otherwise = findSuffix (y*0.001) ss
+                        | otherwise = findSuffix (div y 1000) ss
+
+readOOM :: (Read a, Num a) => ReadM a
+readOOM = eitherReader $ \s -> case reads s of
+    [(n,[ ])] -> Right n
+    [(n,"k")] -> Right $ n * 1000
+    [(n,"M")] -> Right $ n * 1000000
+    [(n,"G")] -> Right $ n * 1000000000
+    [(n,"T")] -> Right $ n * 1000000000000
+    _         -> Left $ "unable to parse: " ++ show s
+
 
 -- 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.
diff --git a/Bio/Util/Text.hs b/Bio/Util/Text.hs
--- a/Bio/Util/Text.hs
+++ b/Bio/Util/Text.hs
@@ -13,6 +13,7 @@
 import qualified Codec.Compression.Zlib.Internal as Z
 import qualified Data.ByteString.Char8           as S
 import qualified Data.ByteString.Lazy            as L
+import qualified Data.ByteString.Lazy.Char8      as C ( unpack )
 import qualified Data.ByteString.Lazy.Internal   as L ( ByteString(..) )
 import qualified Data.Text                       as T
 
@@ -20,6 +21,7 @@
 -- opposite of 'IsString'.
 class Unpack s where unpack :: s -> String
 
+instance Unpack L.ByteString where unpack = C.unpack
 instance Unpack S.ByteString where unpack = S.unpack
 instance Unpack T.Text       where unpack = T.unpack
 instance Unpack String       where unpack = id
diff --git a/CHANGELOG.md b/CHANGELOG.md
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,4 +1,13 @@
-# 2.0.0
+# 2.1 (2020-01-06)
+
+ * Proper error handling.  Instead of calling error, now either a
+   warning is logged or a proper exception is thrown.
+
+ * Faster fastq parser, faster fusion of read pairs.
+
+ * Cleanup of messy code.  Ghc >= 7.10 is now required.
+
+# 2.0 (2019-05-28)
 
  * Switched from "iteratee" to "streaming".  The code looks much cleaner
    and is easier to understand.
diff --git a/Control/Monad/Log.hs b/Control/Monad/Log.hs
new file mode 100644
--- /dev/null
+++ b/Control/Monad/Log.hs
@@ -0,0 +1,247 @@
+{-# LANGUAGE UndecidableInstances, TypeFamilies #-}
+module Control.Monad.Log
+    ( MonadLog(..)
+    , Level(..)
+    , LoggingConf(..)
+    , Logged(..)
+    , LIO
+    , withLogging
+    , withLogging_
+    , logOptions
+    , execWithParser
+    , execWithParser_
+    , PanicCall(..)
+    , panic
+    ) where
+
+import BasePrelude           hiding ( try, catchIOError )
+import Control.Monad.Base           ( MonadBase(..) )
+import Control.Monad.Catch
+import Control.Monad.Primitive
+import Control.Monad.Trans.Class
+import Control.Monad.Trans.Control
+import Control.Monad.Trans.Reader
+import Control.Monad.Trans.RWS.Strict ( RWST )
+import GitVersion                   ( gitFullVersion )
+import Options.Applicative
+import Paths_biohazard              ( version )
+import Streaming
+import System.IO                    ( hPutStr, hPutStrLn, hFlush, stderr, openFile, IOMode(..) )
+
+import qualified Data.Vector                 as V
+
+-- | Severity levels for logging.
+data Level = Debug      -- ^ Message only useful for debugging.  Typically ignored.
+           | Info       -- ^ Purely informative message, e.g. progress reports.  Sometimes printed.
+           | Notice     -- ^ Something remarkable, but harmless.  Sometimes printed, but not collected.
+           | Warning    -- ^ Something unexpected, but usually not a problem.  Typically printed, but not collected.
+           | Error      -- ^ Recoverable error, will normally result in `ExitFailure 1`.  Printed and collected.
+    deriving ( Show, Eq, Ord, Enum, Bounded, Ix )
+
+color_coded :: Level -> String -> String
+color_coded Debug   s = "\27[90m"   ++ s ++ "\27[0m"        -- gray
+color_coded Info    s = "\27[34m"   ++ s ++ "\27[0m"        -- blue
+color_coded Notice  s = "\27[32;1m" ++ s ++ "\27[0m"        -- bold green
+color_coded Warning s = "\27[33m"   ++ s ++ "\27[0m"        -- yellow
+color_coded Error   s = "\27[31;1m" ++ s ++ "\27[0m"        -- bold red
+
+-- | Monads in which messages can be logged.  Any 'Exception' can be
+-- logged; it is reported and/or collected, but does not abort any
+-- computation.
+class Monad m => MonadLog m where
+    -- | Logs a message at a given level.  Depending on settings, the
+    -- message may be printed and/or stored.
+    logMsg :: Exception e => Level -> e -> m ()
+
+    -- | Updates the progress indicator.  The message should not contain
+    -- line feeds, as it is intended to fit on one line and be
+    -- overwritten repeatedly.
+    logString_ :: String -> m ()
+
+    -- | Prints a progress indication.  The message should persist on
+    -- the user's terminal.
+    logStringLn :: String -> m ()
+
+instance (MonadLog m, Monoid w) => MonadLog (RWST r w s m) where
+    logMsg    l e = lift (logMsg    l e)
+    logString_  e = lift (logString_  e)
+    logStringLn e = lift (logStringLn e)
+
+
+-- | Adds logging to any 'MonadIO' type.  Warnings are printed
+-- to stderr immediately, but we remember whether any were emitted.  If
+-- so, we exit with an error code.  The advantage over @WarningT IO@ is
+-- that the warnings are tracked even if the computation exits with an
+-- exception.  Progress indicators are sent to the controlling terminal,
+-- and dicarded if none exists.
+newtype Logged m a = Logged { runLogged :: ReaderT (LoggingConf, Journal) m a }
+  deriving ( Functor, Applicative, Alternative, Monad, MonadTrans, MonadIO, MonadThrow, MonadCatch, MonadMask, MFunctor )
+
+instance MonadTransControl Logged where
+    type StT Logged a = StT (ReaderT (LoggingConf, Journal)) a
+    liftWith = defaultLiftWith Logged runLogged
+    restoreT = defaultRestoreT Logged
+
+instance MonadBase b m => MonadBase b (Logged m) where
+    liftBase = lift . liftBase
+
+instance MonadBaseControl b m => MonadBaseControl b (Logged m) where
+    type StM (Logged m) a = StM (ReaderT (LoggingConf, Journal) m) a
+    liftBaseWith f        = defaultLiftBaseWith f
+    restoreM              = defaultRestoreM
+
+instance PrimMonad m => PrimMonad (Logged m) where
+    type PrimState (Logged m) = PrimState m
+    primitive                 = lift . primitive
+
+type LIO = Logged IO
+
+data LoggingConf = LoggingConf
+    { reporting_level :: Level      -- ^ minimum 'Level' to print a message
+    , logging_level   :: Level      -- ^ minimum 'Level' to remember a message
+    , error_level     :: Level      -- ^ minimum 'Level' that results in a call to 'exitFailure'
+    , max_log_size    :: Int        -- ^ number of messages to keep at any given level
+    , want_progress   :: Bool }
+  deriving Show
+
+data Journal = Journal
+    { logged_messages :: V.Vector (IORef [SomeException])     -- ^ collected messages per level
+    , num_messages    :: V.Vector (IORef Int)                 -- ^ number of collected messages per level
+    , error_exit      :: IORef Bool
+    , cterminal       :: Maybe Handle
+    , spinner         :: IORef String }
+
+instance MonadIO m => MonadLog (Logged m) where
+    logMsg lv e = Logged $ ReaderT $ \(LoggingConf{..},Journal{..}) -> do
+        when (lv >= reporting_level) $ liftIO $ do
+            -- clear spinner
+            forM_ cterminal $ \h -> tryIO $ hPutStr h "\r\27[K" >> hFlush h
+            pn <- getProgName
+            hPutStrLn stderr $ color_coded lv $ printf "%s: [%s] %s" pn (show lv) (displayException e)
+            hFlush stderr
+            -- restore spinner
+            forM_ cterminal $ \h -> readIORef spinner >>= \s ->
+                hPutStr h ("\27[?7l" ++ s ++ "\27[?7h") >> hFlush h
+        when (lv >= logging_level) $ liftIO $
+            atomicModifyIORef' (num_messages V.! fromEnum lv)
+                (\num -> if num < max_log_size then (succ num, True) else (num, False)) >>=
+            flip when (atomicModifyIORef (logged_messages V.! fromEnum lv)
+                (\es -> (toException e : es, ())))
+        when (lv >= error_level) $ liftIO $
+            atomicWriteIORef error_exit True
+
+    logString_ m = Logged $ ReaderT $ \(LoggingConf{..},Journal{..}) ->
+        liftIO $ forM_ cterminal $ \h -> do
+            pn <- getProgName
+            let s = if null m then m else pn ++ ": " ++ m
+            writeIORef spinner s
+            tryIO $ hPutStr h ("\r\27[K\27[?7l" ++ s ++ "\27[?7h") >> hFlush h
+
+    logStringLn m = Logged $ ReaderT $ \(LoggingConf{..},Journal{..}) ->
+        liftIO $ forM_ cterminal $ \h -> do
+            s <- readIORef spinner
+            tryIO $ hPutStr h ("\r\27[K" ++ m ++ "\n\27[?7l" ++ s ++ "\27[?7h") >> hFlush h
+
+
+withLogging_ :: (MonadIO m, MonadMask m) => LoggingConf -> Logged m a -> m a
+withLogging_ conf = withLogging conf >=> either (liftIO . exitWith) pure
+
+withLogging :: (MonadIO m, MonadMask m) => LoggingConf -> Logged m a -> m (Either ExitCode a)
+withLogging conf (Logged k) = do
+    journal <- let n = fromEnum (maxBound :: Level) - fromEnum (minBound :: Level) + 1
+               in liftIO $ Journal <$> V.replicateM n (newIORef [])
+                                   <*> V.replicateM n (newIORef 0)
+                                   <*> newIORef False
+                                   <*> bool (pure Nothing) (tryIO $ openFile "/dev/tty" WriteMode) (want_progress conf)
+                                   <*> newIORef []
+
+    r  <- try $ runReaderT k (conf,journal)
+    liftIO $ do
+        ws  <- V.mapM readIORef (logged_messages journal)
+        nws <- V.mapM readIORef (num_messages journal)
+        pn  <- getProgName
+        forM_ (cterminal journal) $ \h -> do
+            s <- readIORef (spinner journal)
+            tryIO $ unless (null s) (hPutStrLn h []) >> hClose h
+
+        do let eff_warnings  =     [ (l,e) | l <- [minBound ..], l < error_level conf,     e <- ws V.! fromEnum l ]
+               neff_warnings = sum [   n   | l <- [minBound ..], l < error_level conf, let n = nws V.! fromEnum l ]
+           unless (neff_warnings == 0) $ do
+               hPrintf stderr "%s: there were %d warnings\n" pn neff_warnings
+               forM_ eff_warnings $ \(l,e) -> hPutStrLn stderr . color_coded l $ displayException e
+               unless (neff_warnings - length eff_warnings <= 0 || null eff_warnings) $
+                   hPrintf stderr "(and %d more)\n" (neff_warnings - length eff_warnings)
+
+        do let eff_errors    =     [ (l,e) | l <- [error_level conf ..],                   e <- ws V.! fromEnum l ]
+               neff_errors   = sum [     n | l <- [error_level conf ..],               let n = nws V.! fromEnum l ]
+           unless (null eff_errors) $ do
+               hPrintf stderr "%s: there were %d (non-catastrophic) errors\n" pn neff_errors
+               forM_ eff_errors $ \(l,e) -> hPutStrLn stderr . color_coded l $ displayException e
+               unless (neff_errors - length eff_errors <= 0 || null eff_errors) $
+                   hPrintf stderr "(and %d more)\n" (neff_errors - length eff_errors)
+
+        case r of
+          Left  e -> do case fromException e of
+                            Just UserInterrupt -> hPutStrLn stderr $ pn ++ ": Interrupted"
+                            _                  -> hPutStrLn stderr $ pn ++ ": catastrophic error: " ++ displayException e
+                        return . Left $ ExitFailure 2
+
+          Right x -> bool (Right x) (Left $ ExitFailure 1) <$> readIORef (error_exit journal)
+
+
+-- | General wrapper around main.  Runs a command line parser with added
+-- standard options (logging and usage related), runs the actual main
+-- function, prints collected warnings and caught exceptions, and exits
+-- appropriately:  `exitWith (ExitFailure 2)` if an exception was
+-- caught, `exitFailure` if there were warnings of sufficient severity,
+-- and `exitSuccess` otherwise.
+
+execWithParser_ :: Parser a -> Maybe Version -> Maybe String -> InfoMod (a,LoggingConf) -> (a -> LIO b) -> IO b
+execWithParser_ opts prog_ver prog_git_ver inf =
+    execWithParser opts prog_ver prog_git_ver inf >=> either exitWith pure
+
+execWithParser :: Parser a -> Maybe Version -> Maybe String -> InfoMod (a,LoggingConf)
+               -> (a -> LIO b) -> IO (Either ExitCode b)
+execWithParser opts prog_ver prog_git_ver inf k = do
+    pn <- getProgName
+    let verStr = printf "%s%s (%s) using biohazard-%s (%s)" pn
+                        (maybe "" (('-':) . showVersion) prog_ver) (fromMaybe "release" prog_git_ver)
+                        (showVersion version) (fromMaybe "release" gitFullVersion)
+        verOpt = infoOption verStr (short 'V' <> long "version" <> help "Print version number and exit")
+    (a,cf) <- execParser $ info ((,) <$> opts <*> logOptions <* verOpt <* helper) inf
+    withLogging cf (k a)
+
+logOptions :: Parser LoggingConf
+logOptions =
+    LoggingConf
+    <$> (foldl (&) Notice <$> many
+            (flag' more (long "quiet" <> help "Print only important messages") <|>
+             flag' less (long "verbose" <> help "Print also trivial messages")))
+
+    <*> (foldl (&) Warning <$> many
+            (flag' more (long "drop-errors" <> help "Remember only critical messages") <|>
+             flag' less (long "keep-warnings" <> help "Remember also minor messages")))
+
+    <*> (foldl (&) Error <$> many
+            (flag' more (long "warn-ignore" <> help "Fail only after critical errors") <|>
+             flag' less (long "warn-error" <> help "Fail also after warnings")))
+
+    <*> option auto (long "journal-size" <> metavar "NUM" <> help "Hold up to NUM errors in memory" <> value 20)
+    <*> switch (long "progress" <> help "Print progress reports to the terminal")
+  where
+    more, less :: (Enum a, Bounded a, Eq a) => a -> a
+    more a = if a == maxBound then a else succ a
+    less a = if a == minBound then a else pred a
+
+
+-- | An exception than can be thrown when it doesn't seem warranted to
+-- define a custom exception.  Transports a message.
+data PanicCall = PanicCall String deriving (Typeable, Show)
+instance Exception PanicCall where displayException (PanicCall msg) = msg
+
+panic :: MonadIO m => String -> m a
+panic = liftIO . throwIO . PanicCall
+
+tryIO :: IO k -> IO (Maybe k)
+tryIO k = catchIOError (Just <$> k) (\_ -> pure Nothing)
+
diff --git a/Setup.hs b/Setup.hs
--- a/Setup.hs
+++ b/Setup.hs
@@ -1,2 +1,50 @@
-import Distribution.Simple
-main = defaultMain
+import Data.List                            ( isPrefixOf )
+import Distribution.Compat.Exception        ( tryIO )
+import Distribution.PackageDescription      ( PackageDescription(..) )
+import Distribution.Simple                  ( UserHooks(..), defaultMainWithHooks, simpleUserHooks )
+import Distribution.Simple.LocalBuildInfo   ( LocalBuildInfo(..) )
+import Distribution.Simple.Program.Run      ( getProgramInvocationOutput, simpleProgramInvocation )
+import Distribution.Simple.Utils            ( createDirectoryIfMissingVerbose )
+import Distribution.Verbosity               ( Verbosity, silent )
+
+main :: IO ()
+main = defaultMainWithHooks $ simpleUserHooks
+    { buildHook = \pkgDesc lbi hooks flags -> do
+         git_version <- gitDescribe
+         git_branch <- gitBranch
+         let git_full_version = fmap (maybe "" (++":") git_branch ++) git_version
+         createDirectoryIfMissingVerbose silent False (buildDir lbi ++ "/autogen")
+         writeFileIfChanged (buildDir lbi ++ "/autogen/GitVersion.hs") $ unlines
+            [ "module GitVersion (gitVersion,gitBranch,gitFullVersion) where"
+            , "import Prelude"
+            , "gitVersion :: Maybe String"
+            , "gitVersion = " ++ show git_version
+            , "{-# NOINLINE gitVersion #-}"
+            , "gitBranch :: Maybe String"
+            , "gitBranch = " ++ show git_branch
+            , "{-# NOINLINE gitBranch #-}"
+            , "gitFullVersion :: Maybe String"
+            , "gitFullVersion = " ++ show git_full_version
+            , "{-# NOINLINE gitFullVersion #-}" ]
+         buildHook simpleUserHooks pkgDesc lbi hooks flags
+    }
+
+writeFileIfChanged :: FilePath -> String -> IO ()
+writeFileIfChanged fp new = do
+    old <- tryIO $ do s <- readFile fp ; length s `seq` return s
+    if old == Right new then return () else writeFile fp new
+
+gitDescribe :: IO (Maybe String)
+gitDescribe =
+    fmap (either (const Nothing) (Just . unwords . lines)) $ tryIO $
+    getProgramInvocationOutput silent $ simpleProgramInvocation "git" ["describe", "--always", "--long"]
+
+gitBranch :: IO (Maybe String)
+gitBranch =
+    fmap (either (const Nothing) get_branch) $ tryIO $
+    getProgramInvocationOutput silent $ simpleProgramInvocation "git" ["branch"]
+  where
+    get_branch s = case map words . filter ("*" `isPrefixOf`) $ lines s of
+        [_:b:_] -> Just b
+        _       -> Nothing
+
diff --git a/biohazard.cabal b/biohazard.cabal
--- a/biohazard.cabal
+++ b/biohazard.cabal
@@ -1,5 +1,6 @@
+Cabal-version:       2.0
 Name:                biohazard
-Version:             2.0  
+Version:             2.1
 Synopsis:            bioinformatics support library
 Description:         This is a collection of modules I separated from
                      various bioinformatics tools.
@@ -12,16 +13,19 @@
 
 Author:              Udo Stenzel
 Maintainer:          u.stenzel@web.de
-Copyright:           (C) 2010-2017 Udo Stenzel
+Copyright:           (C) 2010-2019 Udo Stenzel
 
-Cabal-version:       >= 1.10
-Build-type:          Simple
-Tested-with:         GHC == 7.8.4, GHC == 7.10.3, GHC == 8.0.2, GHC == 8.2.2, GHC == 8.4.3, GHC == 8.6.1
+Build-type:          Custom
+Tested-with:         GHC == 7.10.3, GHC == 8.0.2, GHC == 8.2.2, GHC == 8.4.3, GHC == 8.6.1, GHC == 8.8.1
 
 source-repository head
   type:     git
   location: https://bitbucket.org/ustenzel/biohazard.git
 
+custom-setup
+  setup-depends:       base                     >= 4.8 && < 4.14,
+                       Cabal                    >= 2.0 && < 3.1
+
 Library
   Exposed-modules:     Bio.Adna,
                        Bio.Align,
@@ -51,31 +55,36 @@
                        Bio.Util.Nub,
                        Bio.Util.Numeric,
                        Bio.Util.Storable,
-                       Bio.Util.Text
+                       Bio.Util.Text,
+                       Control.Monad.Log
 
-  Build-depends:       attoparsec               >= 0.10 && < 0.14,
-                       base                     >= 4.7 && < 4.13,
-                       base-prelude             == 1.2.0.*,
-                       bytestring               >= 0.10.2 && < 0.11,
+  Other-modules:       Paths_biohazard, GitVersion
+  Autogen-modules:     Paths_biohazard, GitVersion
+
+  Build-depends:       attoparsec               >= 0.13 && < 0.14,
+                       base                     >= 4.8 && < 4.14,
+                       base-prelude             >= 1.2 && < 1.2.1,
+                       bytestring               >= 0.10.6 && < 0.11,
                        containers               >= 0.5 && < 0.7,
-                       directory                >= 1.2.1 && < 1.4,
-                       exceptions               >= 0.6 && < 0.11,
-                       hashable                 >= 1.0 && < 1.3,
-                       primitive                >= 0.5 && < 0.7,
+                       directory                >= 1.2.2 && < 1.4,
+                       exceptions               >= 0.6.1 && < 0.11,
+                       hashable                 >= 1.2.3.2 && < 1.4,
+                       monad-control            == 1.0.*,
+                       optparse-applicative     >= 0.13 && < 0.16,
+                       primitive                >= 0.6.1 && < 0.8,
                        stm                      >= 2.4 && < 2.6,
                        streaming                >= 0.1.4.2 && < 0.3,
-                       text                     >= 1.0 && < 1.3,
-                       transformers             >= 0.4.1 && < 0.6,
-                       unix                     >= 2.5 && < 2.8,
-                       unordered-containers     >= 0.2.3 && < 0.3,
+                       text                     >= 1.2.0.2 && < 1.3,
+                       transformers             >= 0.4.2 && < 0.6,
+                       transformers-base        == 0.4.*,
+                       unix                     >= 2.7.1 && < 2.8,
+                       unordered-containers     >= 0.2.5.1 && < 0.3,
                        vector                   >= 0.11 && < 0.13,
-                       vector-algorithms        >= 0.3 && < 0.8,
+                       vector-algorithms        >= 0.8 && < 0.9,
                        zlib                     == 0.6.*
 
-  if !impl(ghc >= 7.10)
-    build-depends: bifunctors == 5.*
   if !impl(ghc >= 8.0)
-    build-depends: semigroups == 0.18.*
+    build-depends: semigroups >= 0.18 && < 0.20
 
   Ghc-options:         -Wall
   if impl(ghc >= 8.0)
@@ -103,16 +112,18 @@
                        ExistentialQuantification,
                        ForeignFunctionInterface,
                        Rank2Types,
-                       TypeFamilies
+                       StandaloneDeriving,
+                       TypeFamilies,
+                       UndecidableInstances
 
+  Build-tool-depends:  hsc2hs:hsc2hs
   Hs-source-dirs:      .
   Include-dirs:        cbits
-  Install-Includes:    myers_align.h
+  Install-includes:    myers_align.h
   C-sources:           cbits/loops.c,
-                       cbits/mmap.c,
                        cbits/myers_align.c,
                        cbits/trim.c
+                       cbits/util.c
                        cbits/zlib.c
-  CC-options:          -fPIC
 
 -- :vim:tw=132:
diff --git a/cbits/mmap.c b/cbits/mmap.c
deleted file mode 100644
--- a/cbits/mmap.c
+++ /dev/null
@@ -1,10 +0,0 @@
-#include <sys/mman.h>
-
-unsigned char *my_mmap(size_t len, int fd) {
-        void *result = mmap(0, len, PROT_READ, MAP_SHARED, fd, 0);
-        return (unsigned char*)( result == MAP_FAILED ? 0 : result );
-}
-
-void my_munmap(void *len, unsigned char *p) {
-        munmap( p, (size_t)len ) ; 
-}
diff --git a/cbits/trim.c b/cbits/trim.c
--- a/cbits/trim.c
+++ b/cbits/trim.c
@@ -1,28 +1,32 @@
 #include <stdint.h>
+#include <limits.h>
 
-static const uint8_t compls[] =
-    { 0, 8, 4, 12, 2, 10, 6, 14, 1, 9, 5, 13, 3, 11, 7, 15 } ;
+#define MIN(a,b) ((a)<(b)?(a):(b))
+#define MAX(a,b) ((a)>(b)?(a):(b))
 
-int prim_match_reads( int i1
-                    , int i2
-                    , int r
+/* This version receives the second sequence reverse-complemented and
+ * the second quality vector reversed.  Allows for a forward loop,
+ * doesn't need the lookup table. */
+int prim_match_reads( int i1       // matches from i2
+                    , int i2       // ???
+                    , int r        // length to match
                     , const uint8_t *rd1
                     , const uint8_t *qs1
-                    , const uint8_t *rd2
-                    , const uint8_t *qs2 )
+                    , const uint8_t *rc_rd2
+                    , const uint8_t *rv_qs2 )
 {
     int acc = 0 ;
     while( r != 0 )
     {
-        --i2 ;
         uint8_t n1 = rd1[ i1 ] ;
-        uint8_t n2 = rd2[ i2 ] ;
+        uint8_t n2 = rc_rd2[ i2 ] ;
         uint8_t q1 = qs1[ i1 ] ;
-        uint8_t q2 = qs2[ i2 ] ;
+        uint8_t q2 = rv_qs2[ i2 ] ;
 
-        acc += (n1 & 0xF) == compls[ n2 & 0xF ] ? 0 : 5 + (q1 < q2 ? q1 : q2) ;
+        acc += (n1 & 0xF) == (n2 & 0xF) ? 0 : 5 + (q1 < q2 ? q1 : q2) ;
 
         ++i1 ;
+        ++i2 ;
         --r ;
     }
     return acc ;
@@ -42,5 +46,68 @@
                (qs[ i+off ] < 25 ? qs[ i+off ] : 25) ;
     }
     return acc ;
+}
+
+// return the length that minimizes the score
+int prim_merge_score (
+        uint8_t **p_fwd_ads,
+        int *p_fwd_lns,
+        int n_fwd_adapters,
+
+        uint8_t **p_rev_ads,
+        int *p_rev_lns,
+        int n_rev_adapters,
+
+        uint8_t *p_rd1,
+        uint8_t *p_qs1,
+        int l1,
+        uint8_t *p_rd2,
+        uint8_t *p_qs2,
+        int l2,
+        uint8_t *p_rrd2,
+        uint8_t *p_rqs2,
+
+        int *scores )
+{
+    int minl = -1 ;
+    int min0 = 6 * (l1+l2) ;
+    int min1= INT_MAX ;
+
+    for( int l = 0 ; l <= l1+l2 ; ++l )
+    {
+        int ff = INT_MAX ;
+        for( int i = 0 ; i != n_fwd_adapters ; ++i )
+        {
+            int efflength = MIN(l1 - l, p_fwd_lns[i]) ;
+            int x = prim_match_ad( l, efflength, p_rd1, p_qs1, p_fwd_ads[i] ) ;
+            ff = MIN(ff, x + 6 * MAX( 0, l1 - p_fwd_lns[i] - l)) ;
+        }
+
+        int rr = INT_MAX ;
+        for( int i = 0 ; i != n_rev_adapters ; ++i )
+        {
+            int efflength = MIN(l2 - l,p_rev_lns[i]) ;
+            int x = prim_match_ad( l, efflength, p_rd2, p_qs2, p_rev_ads[i] ) ;
+            rr = MIN(rr, x + 6 * MAX( 0, l1 - p_rev_lns[i] - l)) ;
+        }
+
+        int minidx1 = MAX(l - l2, 0) ; // vec1, forward
+        int minidx2 = MAX(l2 - l, 0) ; // vec2, forward-on-reverse
+        int efflength = MIN(l1 + l2 - l, l) ; // effective length
+
+        int mm = efflength <= 0 ? 0 : prim_match_reads( minidx1, minidx2, efflength, p_rd1, p_qs1, p_rrd2, p_rqs2 ) ;
+        int z = 6 * l + ff + rr + mm ;
+
+        if( z < min0 ) {
+            min1 = min0 ;
+            min0 = z ;
+            minl = l ;
+        } else if( z < min1 ) {
+            min1 = z ;
+        }
+    }
+    scores[0] = min1 - min0 ;
+    scores[1] = 6*(l1+l2) - min0 ;
+    return minl ;
 }
 
diff --git a/cbits/util.c b/cbits/util.c
new file mode 100644
--- /dev/null
+++ b/cbits/util.c
@@ -0,0 +1,8 @@
+#include <sys/mman.h>
+
+void my_munmap(void *len, unsigned char *p)
+{
+    munmap( p, (size_t)len ) ;
+}
+
+