biohazard 2.0 → 2.1

raw patch · 29 files changed

+1286/−707 lines, 29 filesdep +monad-controldep +optparse-applicativedep +transformers-basedep −bifunctorsdep ~attoparsecdep ~basedep ~base-preludebuild-type:Customsetup-changedPVP ok

version bump matches the API change (PVP)

Dependencies added: monad-control, optparse-applicative, transformers-base

Dependencies removed: bifunctors

Dependency ranges changed: attoparsec, base, base-prelude, bytestring, directory, exceptions, hashable, primitive, semigroups, text, transformers, unix, unordered-containers, vector-algorithms

API changes (from Hackage documentation)

Files

• Bio/Bam/Evan.hs +18/−1
• Bio/Bam/Fastq.hs +184/−97
• Bio/Bam/Filter.hs +9/−7
• Bio/Bam/Header.hs +26/−28
• Bio/Bam/Index.hs +8/−7
• Bio/Bam/Pileup.hs +2/−11
• Bio/Bam/Reader.hs +39/−40
• Bio/Bam/Rec.hs +37/−13
• Bio/Bam/Rmdup.hs +116/−103
• Bio/Bam/Trim.hs +163/−188
• Bio/Bam/Writer.hs +7/−7
• Bio/Base.hs +41/−39
• Bio/Prelude.hs +6/−4
• Bio/Streaming.hs +28/−16
• Bio/Streaming/Bgzf.hs +35/−8
• Bio/Streaming/Bytes.hs +17/−7
• Bio/Streaming/Parse.hs +28/−22
• Bio/TwoBit.hs +26/−21
• Bio/Util/MMap.hs +0/−28
• Bio/Util/MMap.hsc +43/−0
• Bio/Util/Numeric.hs +21/−10
• Bio/Util/Text.hs +2/−0
• CHANGELOG.md +10/−1
• Control/Monad/Log.hs +247/−0
• Setup.hs +50/−2
• biohazard.cabal +37/−26
• cbits/mmap.c +0/−10
• cbits/trim.c +78/−11
• cbits/util.c +8/−0

@@ -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

@@ -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 

@@ -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 }  

@@ -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.

@@ -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. 

@@ -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 ()

@@ -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

@@ -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

@@ -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 }

@@ -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 #-} 

@@ -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')

@@ -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 #-}

@@ -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(..) ) 

@@ -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

@@ -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

@@ -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' #-}  -- --------------------------------------------------------------------------

@@ -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') 

@@ -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 

@@ -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 ())--

@@ -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 ()+

@@ -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.

@@ -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

@@ -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.

@@ -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)+

@@ -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+

@@ -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:

@@ -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 ) ; -}

@@ -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 ; } 

@@ -0,0 +1,8 @@+#include <sys/mman.h>++void my_munmap(void *len, unsigned char *p)+{+    munmap( p, (size_t)len ) ;+}++