biohazard 0.6.16 → 1.0.0

raw patch · 31 files changed

+351/−587 lines, 31 filesdep ~basedep ~base-preludedep ~containersPVP ok

version bump matches the API change (PVP)

Dependency ranges changed: base, base-prelude, containers

API changes (from Hackage documentation)

Files

• biohazard.cabal +10/−14
• src/Bio/Adna.hs +107/−102
• src/Bio/Bam.hs +2/−2
• src/Bio/Bam/Evan.hs +6/−5
• src/Bio/Bam/Fastq.hs +17/−17
• src/Bio/Bam/Filter.hs +2/−2
• src/Bio/Bam/Header.hs +5/−4
• src/Bio/Bam/Index.hs +3/−55
• src/Bio/Bam/Pileup.hs +39/−130
• src/Bio/Bam/Reader.hs +30/−30
• src/Bio/Bam/Rec.hs +1/−8
• src/Bio/Bam/Regions.hs +5/−10
• src/Bio/Bam/Rmdup.hs +13/−15
• src/Bio/Bam/Trim.hs +3/−3
• src/Bio/Bam/Writer.hs +5/−17
• src/Bio/Base.hs +1/−6
• src/Bio/Iteratee.hs +11/−17
• src/Bio/Iteratee/Base.hs +0/−12
• src/Bio/Iteratee/Bgzf.hsc +0/−2
• src/Bio/Iteratee/Builder.hs +12/−18
• src/Bio/Iteratee/Bytes.hs +1/−1
• src/Bio/Iteratee/Exception.hs +4/−4
• src/Bio/Iteratee/IO.hs +2/−2
• src/Bio/Iteratee/Iteratee.hs +3/−4
• src/Bio/Iteratee/List.hs +2/−5
• src/Bio/Prelude.hs +2/−28
• src/Bio/TwoBit.hs +17/−19
• src/Bio/Util/MMap.hs +1/−1
• src/Bio/Util/Numeric.hs +40/−46
• src/Bio/Util/Zlib.hs +3/−4
• src/cbits/myers_align.h +4/−4

@@ -1,31 +1,26 @@ Name:                biohazard-Version:             0.6.16+Version:             1.0.0 Synopsis:            bioinformatics support library Description:         This is a collection of modules I separated from                      various bioinformatics tools. Category:            Bioinformatics -Homepage:            http://github.com/udo-stenzel/biohazard+Homepage:            https://bitbucket.org/ustenzel/biohazard License:             BSD3 License-File:        LICENSE  Author:              Udo Stenzel-Maintainer:          udo.stenzel@eva.mpg.de+Maintainer:          u.stenzel@web.de Copyright:           (C) 2010-2017 Udo Stenzel  Cabal-version:       >= 1.10 Build-type:          Simple-Tested-With:         GHC == 7.8.4, GHC == 7.10.1, GHC == 8.0.1+Tested-With:         GHC == 7.8.4, GHC == 7.10.3, GHC == 8.0.2, GHC == 8.2.1  source-repository head   type:     git   location: https://bitbucket.org/ustenzel/biohazard.git -Flag debug-  Description: enable additional sanity checks-  Default:     False-  Manual:      True- Library   Exposed-modules:     Bio.Adna,                        Bio.Align,@@ -61,11 +56,11 @@    Build-depends:       async                    >= 2.0 && < 2.2,                        attoparsec               >= 0.10 && < 0.14,-                       base                     >= 4.6 && < 4.11,-                       base-prelude             == 1.0.* || == 1.2.*,+                       base                     >= 4.7 && < 4.11,+                       base-prelude             == 1.2.*,                        binary                   >= 0.7 && < 0.9,                        bytestring               >= 0.10.2 && < 0.11,-                       containers               >= 0.4.1 && < 0.6,+                       containers               == 0.5.*,                        directory                >= 1.2 && < 1.4,                        exceptions               >= 0.6 && < 0.9,                        filepath                 >= 1.3 && < 1.5,@@ -91,6 +86,7 @@                        DeriveDataTypeable,                        FlexibleContexts,                        FlexibleInstances,+                       LambdaCase,                        MultiParamTypeClasses,                        NoImplicitPrelude,                        OverloadedStrings,@@ -113,12 +109,12 @@                        UndecidableInstances    Hs-source-dirs:      src-  Install-Includes:    src/cbits/myers_align.h+  Include-dirs:        src/cbits+  Install-Includes:    myers_align.h   C-sources:           src/cbits/loops.c,                        src/cbits/mmap.c,                        src/cbits/myers_align.c,                        src/cbits/trim.c   CC-options:          -fPIC-  -- :vim:tw=132:

@@ -1,4 +1,16 @@-{-# LANGUAGE DeriveGeneric, CPP #-}+{-# LANGUAGE DeriveGeneric #-}++-- | Things specific to ancient DNA, e.g. damage models.+--+-- For aDNA, we need a substitution probability.  We have three options:+-- use an empirically determined PSSM, use an arithmetically defined+-- PSSM based on the /Johnson/ model, use a context sensitive PSSM based+-- on the /Johnson/ model and an alignment.  Using /Dindel/, actual+-- substitutions relative to a called haplotype would be taken into+-- account.  Since we're not going to do that, taking alignments into+-- account is difficult, somewhat approximate, and therefore not worth+-- the hassle.+ module Bio.Adna (     DmgStats(..),     CompositionStats,@@ -17,7 +29,11 @@     Alignment(..),     FragType(..),     Subst(..),+     NPair,+    npair,+    fst_np,+    snd_np,      noDamage,     univDamage,@@ -41,20 +57,9 @@ import qualified Data.Vector.Unboxed            as U import qualified Data.Vector.Unboxed.Mutable    as UM --- ^ Things specific to ancient DNA, e.g. damage models.------ For aDNA, we need a substitution probability.  We have three options:--- use an empirically determined PSSM, use an arithmetically defined--- PSSM based on the /Johnson/ model, use a context sensitive PSSM based--- on the /Johnson/ model and an alignment.  Using /Dindel/, actual--- substitutions relative to a called haplotype would be taken into--- account.  Since we're not going to do that, taking alignments into--- account is difficult, somewhat approximate, and therefore not worth--- the hassle.------ We represent substitution matrices by the type 'Mat44D'.  Internally,+-- | We represent substitution matrices by the type 'Mat44D'.  Internally, -- this is a vector of packed vectors.  Conveniently, each of the packed--- vectors represents all transition /into/ the given nucleotide.+-- vectors represents all transitions /into/ the given nucleotide.  newtype Mat44D = Mat44D (U.Vector Double) deriving (Show, Generic) newtype MMat44D = MMat44D (UM.IOVector Double)@@ -250,8 +255,6 @@ --   mind, 'substs5d5', 'substs5d3', 'substs5dd' are like 'substs5', but --   counting only reads where the 5' end is damaged, where the 3' end --   is damaged, and where both ends are damaged, respectively.------ XXX  This got kind of ugly.  We'll see where this goes...  data DmgStats a = DmgStats {     basecompo5 :: CompositionStats,@@ -274,14 +277,34 @@   data FragType = Complete | Leading | Trailing deriving (Show, Eq)-type NPair = ( Nucleotides, Nucleotides ) --- Alignment record, might have been gotten from practically anywhere--- with varying completeness.  We record anything we can get, most is--- optional.  Reference sequence is filled with Ns if missing.+-- | Compact storage of a pair of ambiguous 'Nucleotides'.  Used to+-- represent alignments in a way that is accessible even to assembly+-- code.  The first and sencond field are stored in the low and high+-- nybble, respectively.  See 'fst_np', 'snd_np', 'npair'.+newtype NPair = NPair Word8 deriving (Eq, Ord)++npair :: Nucleotides -> Nucleotides -> NPair+npair (Ns r) (Ns q) = NPair $ shiftL q 4 .|. r .&. 0xF++fst_np, snd_np :: NPair -> Nucleotides+fst_np (NPair w) = Ns (w .&. 0xF)+snd_np (NPair w) = Ns (shiftR w 4)++instance Storable NPair where+    sizeOf    _ = 1+    alignment _ = 1+    peek p = NPair <$> peek (castPtr p :: Ptr Word8)+    poke p (NPair v) = poke (castPtr p :: Ptr Word8) v++instance Show NPair where+    showsPrec _ p = shows (fst_np p) . (:) '/' . shows (snd_np p)++-- | Alignment record.  The reference sequence is filled with Ns if+-- missing. data Alignment = ALN-    { a_sequence :: !(U.Vector NPair)       -- the alignment proper-    , a_fragment_type :: !FragType }    -- was the adapter trimmed?+    { a_sequence :: !(VS.Vector NPair)      -- the alignment proper+    , a_fragment_type :: !FragType }        -- was the adapter trimmed?  addFragType :: BamMeta -> Enumeratee [BamRaw] [(BamRaw,FragType)] m b addFragType meta = mapStream $ \br -> (br, case unpackBam br of@@ -353,15 +376,15 @@         guard (not $ isUnmapped b)         md <- getMd b         let pps = alnFromMd b_seq b_cigar md-            ref = U.map fromN $ U.filter ((/=) gap . fst) pps+            ref = U.convert $ VS.map fromN $ VS.filter ((/=) gap) $ VS.map fst_np pps         return (b, ft, ref, pps)) =$     damagePatternsIter 0 rng it   where-    fromN (ns,_) | ns == nucsA = 2-                 | ns == nucsC = 1-                 | ns == nucsG = 3-                 | ns == nucsT = 0-                 | otherwise   = 4+    fromN ns | ns == nucsA = 2+             | ns == nucsC = 1+             | ns == nucsG = 3+             | ns == nucsT = 0+             | otherwise   = 4  -- | Common logic for statistics. The function 'get_ref_and_aln' -- reconstructs reference sequence and alignment from a Bam record.  It@@ -370,7 +393,7 @@ damagePatternsIter :: MonadIO m                    => Int -> Int                    -> Iteratee [Alignment] m b-                   -> Iteratee [(BamRec, FragType, U.Vector Word8, U.Vector NPair)] m (DmgStats b)+                   -> Iteratee [(BamRec, FragType, U.Vector Word8, VS.Vector NPair)] m (DmgStats b) damagePatternsIter ctx rng it = mapStream revcom_both =$ do     let maxwidth = ctx + rng     acc_bc <- liftIO $ UM.replicate (2 * 5 *    maxwidth) (0::Int)@@ -378,13 +401,6 @@     acc_cg <- liftIO $ UM.replicate (2 * 2 * 4 *     rng) (0::Int)      it' <- flip mapStreamM it $ \(BamRec{..}, a_fragment_type, ref, a_sequence) -> liftIO $ do-#ifdef DEBUG-              when (U.any (<0) ref || U.any (>4) ref) . error $-                    "Unexpected value in reference fragment: " ++ show ref-#endif-              let good_pairs     = U.indexed             a_sequence-                  good_pairs_rev = U.indexed $ U.reverse a_sequence-               -- basecompositon near 5' end, near 3' end               let (width5, width3) = case a_fragment_type of                                             Leading -> (full_width, 0)@@ -397,51 +413,52 @@                -- For substitutions, decide what damage class we're in:               -- 0 - no damage, 1 - damaged 5' end, 2 - damaged 3' end, 3 - both-              let dmgbase = 2*4*4*rng * ( (if U.null a_sequence || U.head a_sequence /= (nucsC,nucsT) then 1 else 0)-                                        + (if U.null a_sequence || U.last a_sequence /= (nucsC,nucsT) then 2 else 0) )+              let dmgbase = 2*4*4*rng * ( (if VS.null a_sequence || VS.head a_sequence /= npair nucsC nucsT then 1 else 0)+                                        + (if VS.null a_sequence || VS.last a_sequence /= npair nucsC nucsT then 2 else 0) )                -- substitutions near 5' end               let len_at_5 = case a_fragment_type of Leading  -> min rng (G.length b_seq)                                                      Complete -> min rng (G.length b_seq `div` 2)                                                      Trailing -> 0-              U.forM_ (U.take len_at_5 good_pairs) $-                    \(i,uv) -> withPair uv $ \j -> bump (j * rng + i + dmgbase) acc_st+              flip G.imapM_ (VS.take len_at_5 a_sequence) $+                    \i uv -> withPair uv $ \j -> bump (j * rng + i + dmgbase) acc_st                -- substitutions at CpG sites near 5' end-              U.zipWithM_-                  (\(i,(u,v)) (_,(w,z)) ->-                      when (u == nucsC && w == nucsG) $ do-                          withNs v $ \y -> bump (  y   * rng +  i ) acc_cg-                          withNs z $ \y -> bump ((y+4) * rng + i+1) acc_cg)-                  (U.take len_at_5 good_pairs) (U.drop 1 good_pairs)+              G.izipWithM_+                  (\i uv wz ->+                      when (fst_np uv == nucsC && fst_np wz == nucsG) $ do+                          withNs (snd_np uv) $ \y -> bump (  y   * rng +  i ) acc_cg+                          withNs (snd_np wz) $ \y -> bump ((y+4) * rng + i+1) acc_cg)+                  (VS.take len_at_5 a_sequence) (VS.drop 1 a_sequence)                -- substitutions near 3' end               let len_at_3 = case a_fragment_type of Leading  -> 0                                                      Complete -> min rng (G.length b_seq `div` 2)                                                      Trailing -> min rng (G.length b_seq)-              U.forM_ (U.take len_at_3 good_pairs_rev) $-                    \(i,uv) -> withPair uv $ \j -> bump ((17+j) * rng -i -1 + dmgbase) acc_st+              flip G.imapM_ (VS.take len_at_3 (VS.reverse a_sequence)) $+                    \i uv -> withPair uv $ \j -> bump ((17+j) * rng -i -1 + dmgbase) acc_st                -- substitutions at CpG sites near 3' end-              U.zipWithM_-                  (\(_,(u,v)) (i,(w,z)) ->-                      when (u == nucsC && w == nucsG) $ do-                          withNs v $ \y -> bump ((y+ 9) * rng - i-2) acc_cg-                          withNs z $ \y -> bump ((y+13) * rng - i-1) acc_cg)-                  (U.drop 1 good_pairs_rev) (U.take len_at_3 good_pairs_rev)-+              G.izipWithM_+                  (\i wz uv ->+                      when (fst_np uv == nucsC && fst_np wz == nucsG) $ do+                          withNs (snd_np uv) $ \y -> bump ((y+ 9) * rng - i-2) acc_cg+                          withNs (snd_np wz) $ \y -> bump ((y+13) * rng - i-1) acc_cg)+                  (VS.take len_at_3 (VS.reverse a_sequence))+                  (VS.drop 1 (VS.reverse a_sequence))                return ALN{..} +     let nsubsts = 4*4*rng         mk_substs off = sequence [ (,) (n1 :-> n2) <$> U.unsafeFreeze (UM.slice ((4*i+j)*rng + off*nsubsts) rng acc_st)                                  | (i,n1) <- zip [0..] [nucA..nucT]                                  , (j,n2) <- zip [0..] [nucA..nucT] ]      accs <- liftIO $ DmgStats <$> sequence [ (,) nuc <$> U.unsafeFreeze (UM.slice (i*maxwidth) maxwidth acc_bc)-                                           | (i,nuc) <- zip [2,1,3,0,4] [Just nucA, Just nucC, Just nucG, Just nucT, Nothing] ]+                                           | (i,nuc) <- zip [2,1,3,0,4] [Just nucA,Just nucC,Just nucG,Just nucT,Nothing] ]                               <*> sequence [ (,) nuc <$> U.unsafeFreeze (UM.slice (i*maxwidth) maxwidth acc_bc)-                                           | (i,nuc) <- zip [7,6,8,5,9] [Just nucA, Just nucC, Just nucG, Just nucT, Nothing] ]+                                           | (i,nuc) <- zip [7,6,8,5,9] [Just nucA,Just nucC,Just nucG,Just nucT,Nothing] ]                                <*> mk_substs 0                               <*> mk_substs 1@@ -469,18 +486,16 @@                    , substs3d3 = mconcat [ substs3d3 accs', substs3dd accs'] }   where     {-# INLINE withPair #-}-    withPair (Ns u, Ns v) k = case pairTab `U.unsafeIndex` fromIntegral (16*u+v) of-         j -> if j >= 0 then k j else return ()+    withPair (NPair i) k =+        case pairTab `U.unsafeIndex` fromIntegral i of+            j -> when (j >= 0) (k j)      !pairTab = U.replicate 256 (-1) U.//-            [ (fromIntegral $ 16*u+v, x*4+y) | (Ns u,x) <- zip [nucsA, nucsC, nucsG, nucsT] [0,1,2,3]-                                             , (Ns v,y) <- zip [nucsA, nucsC, nucsG, nucsT] [0,1,2,3] ]+            [ (fromIntegral i, x*4+y) | (u,x) <- zip [nucsA, nucsC, nucsG, nucsT] [0,1,2,3]+                                      , (v,y) <- zip [nucsA, nucsC, nucsG, nucsT] [0,1,2,3]+                                      , let NPair i = npair u v ]     {-# INLINE bump #-}-#ifdef DEBUG-    bump i v = UM.read v i >>= UM.write v i . succ-#else     bump i v = UM.unsafeRead v i >>= UM.unsafeWrite v i . succ-#endif      {-# INLINE withNs #-}     withNs ns k | ns == nucsA = k 0@@ -533,20 +548,20 @@                                      | otherwise        = (x :-> y, U.zipWith (+) u v)  -revcom_both :: ( BamRec, FragType, U.Vector Word8, U.Vector (Nucleotides, Nucleotides) )-            -> ( BamRec, FragType, U.Vector Word8, U.Vector (Nucleotides, Nucleotides) )+revcom_both :: ( BamRec, FragType, U.Vector Word8, VS.Vector NPair )+            -> ( BamRec, FragType, U.Vector Word8, VS.Vector NPair ) revcom_both (b, ft, ref, pps)     | isReversed b = ( b, ft, revcom_ref ref, revcom_pairs pps )     | otherwise    = ( b, ft,            ref,              pps )   where-    revcom_ref   = U.reverse . U.map (\c -> if c > 3 then c else xor c 2)-    revcom_pairs = U.reverse . U.map (compls *** compls)+    revcom_ref   =  U.reverse .  U.map (\c -> if c > 3 then c else xor c 2)+    revcom_pairs = VS.reverse . VS.map (\p -> npair (compls $ fst_np p) (compls $ snd_np p))   -- | Reconstructs the alignment from reference, query, and cigar.  Only -- positions where the query is not gapped are produced.-aln_from_ref :: U.Vector Word8 -> Vector_Nucs_half Nucleotides -> VS.Vector Cigar -> U.Vector NPair-aln_from_ref ref0 qry0 cig0 = U.fromList $ step ref0 qry0 cig0+aln_from_ref :: U.Vector Word8 -> Vector_Nucs_half Nucleotides -> VS.Vector Cigar -> VS.Vector NPair+aln_from_ref ref0 qry0 cig0 = VS.fromList $ step ref0 qry0 cig0   where     step ref qry cig1         | U.null ref || G.null qry || G.null cig1 = []@@ -554,14 +569,14 @@                       case G.unsafeTail cig1 of { cig ->                       case op of { -        Mat -> zipWith (\r q -> (nn r,q)) (G.toList (G.take n ref))-                                          (G.toList (G.take n qry)) ++ step (G.drop n ref) (G.drop n qry) cig ;-        Del ->                                                         step (G.drop n ref)           qry  cig ;-        Ins ->    map (\q -> ( gap,  q )) (G.toList (G.take n qry)) ++ step           ref  (G.drop n qry) cig ;-        SMa ->    map (\q -> ( gap,  q )) (G.toList (G.take n qry)) ++ step           ref  (G.drop n qry) cig ;-        HMa ->   replicate n (gap, nucsN)                           ++ step           ref            qry  cig ;-        Nop ->                                                         step           ref            qry  cig ;-        Pad ->                                                         step           ref            qry  cig }}}+        Mat -> zipWith (npair . nn) (G.toList (G.take n ref))+                                    (G.toList (G.take n qry)) ++ step (G.drop n ref) (G.drop n qry) cig ;+        Del ->                                                   step (G.drop n ref)           qry  cig ;+        Ins -> map (npair gap) (G.toList (G.take n qry))      ++ step           ref  (G.drop n qry) cig ;+        SMa -> map (npair gap) (G.toList (G.take n qry))      ++ step           ref  (G.drop n qry) cig ;+        HMa -> replicate n (npair gap nucsN)                  ++ step           ref            qry  cig ;+        Nop ->                                                   step           ref            qry  cig ;+        Pad ->                                                   step           ref            qry  cig }}}      nn 0 = nucsT     nn 1 = nucsC@@ -572,8 +587,8 @@  -- | Reconstructs the alignment from query, cigar, and md.  Only -- positions where the query is not gapped are produced.-alnFromMd :: Vector_Nucs_half Nucleotides -> VS.Vector Cigar -> [MdOp] -> U.Vector NPair-alnFromMd qry0 cig0 md0 = U.fromList $ step qry0 cig0 md0+alnFromMd :: Vector_Nucs_half Nucleotides -> VS.Vector Cigar -> [MdOp] -> VS.Vector NPair+alnFromMd qry0 cig0 md0 = VS.fromList $ step qry0 cig0 md0   where     step qry cig1 md         | G.null qry || G.null cig1 || null md = []@@ -584,21 +599,23 @@     step' qry op n cig (MdDel [] : md) = step' qry op n cig md      step' qry Mat n cig (MdNum m : md)-            | n <  m = map (\q -> (q,q)) (G.toList (G.take n qry)) ++ step  (G.drop n qry)           cig (MdNum (m-n) : md)-            | n >  m = map (\q -> (q,q)) (G.toList (G.take m qry)) ++ step' (G.drop m qry) Mat (n-m) cig                md-            | n == m = map (\q -> (q,q)) (G.toList (G.take n qry)) ++ step  (G.drop n qry)           cig                md-    step' qry Mat n cig (MdRep c : md) =         ( c, G.head qry )  : step' (G.tail   qry) Mat (n-1) cig                md+            | n <  m =    map twin (G.toList (G.take n qry)) ++ step  (G.drop n qry)           cig (MdNum (m-n) : md)+            | n >  m =    map twin (G.toList (G.take m qry)) ++ step' (G.drop m qry) Mat (n-m) cig                md+            | n == m =    map twin (G.toList (G.take n qry)) ++ step  (G.drop n qry)           cig                md+    step' qry Mat n cig (MdRep c : md) = npair c (G.head qry) : step' (G.tail   qry) Mat (n-1) cig                md     step'   _ Mat _   _          _     = []      step' qry Del n cig (MdDel (_:ss) : md) = step' qry Del (n-1) cig (MdDel ss : md)     step'   _ Del _   _               _     = [] -    step' qry Ins n cig                 md  = map ((,) gap) (G.toList (G.take n qry)) ++ step (G.drop n qry) cig md-    step' qry SMa n cig                 md  = map ((,) gap) (G.toList (G.take n qry)) ++ step (G.drop n qry) cig md-    step' qry HMa n cig                 md  =                replicate n (gap, nucsN) ++ step           qry  cig md-    step' qry Nop _ cig                 md  =                                            step           qry  cig md-    step' qry Pad _ cig                 md  =                                            step           qry  cig md+    step' qry Ins n cig                 md  = map (npair gap) (G.toList (G.take n qry)) ++ step (G.drop n qry) cig md+    step' qry SMa n cig                 md  = map (npair gap) (G.toList (G.take n qry)) ++ step (G.drop n qry) cig md+    step' qry HMa n cig                 md  =             replicate n (npair gap nucsN) ++ step           qry  cig md+    step' qry Nop _ cig                 md  =                                              step           qry  cig md+    step' qry Pad _ cig                 md  =                                              step           qry  cig md +    twin q = npair q q+ -- | Number of mismatches allowed by BWA. -- @bwa_cal_maxdiff thresh len@ returns the number of mismatches -- @bwa aln -n $tresh@ would allow in a read of length @len@.  For@@ -620,18 +637,6 @@ --      return 2; --   } -- @---      double sum, y = 1.0;---      int k, x = 1;---      for (k = 1, sum = elambda; k < 1000; ++k) {---          y *= l * err;---          x *= k;---          sum += elambda * y / x;---          if (1.0 - sum < thres) return k;---      }---      return 2;---   }--- @---  bwa_cal_maxdiff :: Double -> Int -> Int bwa_cal_maxdiff thresh len = k_fin-1

@@ -1,3 +1,5 @@+-- | Umbrella module for most of what's under 'Bio.Bam'.+ module Bio.Bam (     module Bio.Bam.Fastq,     module Bio.Bam.Filter,@@ -19,6 +21,4 @@ import Bio.Bam.Trim import Bio.Bam.Writer import Bio.Iteratee---- ^ Umbrella module for most of what's under 'Bio.Bam'. 

@@ -14,7 +14,7 @@  -- | Fixes abuse of flags valued 0x800 and 0x1000.  We used them for -- low quality and low complexity, but they have since been redefined.--- If set, we clear them and store them into the ZD field.  Also fixes+-- If set, we clear them and store them into the ZQ field.  Also fixes -- abuse of the combination of the paired, 1st mate and 2nd mate flags -- used to indicate merging or trimming.  These are canonicalized and -- stored into the FF field.  This function is unsafe on BAM files of@@ -38,9 +38,9 @@         is_trimmed = flags' .&. (flagPaired .|. flagFirstMate .|. flagSecondMate) == flagSecondMate         newflags = (if is_merged then eflagMerged else 0) .|. (if is_trimmed then eflagTrimmed else 0) -        -- Extended flags, renamed to avoid collision with BWA Goes like this:  if FF is there, use-        -- it.  Else check if XF is there _and_is_numeric_.  If so, use it and remove it, and set FF-        -- instead.  Else use 0 and leave it alone.  Note that this solves the collision with BWA,+        -- Extended flags, renamed to avoid collision with BWA.  Goes like this:  if FF is there, use+        -- it.  Else check if XF is there __and is numeric__.  If so, use it, remove it, and set FF+        -- instead.  Else use 0 and leave it alone.  Note that this resolves the collision with BWA,         -- since BWA puts a character there, not an int.         cleaned_exts = case (lookup "FF" (b_exts b), lookup "XF" (b_exts b)) of                 ( Just (Int i), _ ) -> updateE "FF" (Int (i .|. newflags))                (b_exts b)@@ -51,7 +51,7 @@  -- | Fixes typical inconsistencies produced by Bwa: sometimes, 'mate unmapped' should be set, and we -- can see it, because we match the mate's coordinates.  Sometimes 'properly paired' should not be--- set, because one mate in unmapped.  This function is generally safe, but needs to be called only+-- set, because one mate is unmapped.  This function is generally safe, but needs to be called only -- on the output of affected (older?) versions of Bwa. fixupBwaFlags :: BamRec -> BamRec fixupBwaFlags b = b { b_flag = fixPP $ b_flag b .|. if mu then flagMateUnmapped else 0 }@@ -65,6 +65,7 @@         -- If either mate is unmapped, remove "properly paired".         fixPP f | f .&. (flagUnmapped .|. flagMateUnmapped) == 0 = f                 | otherwise = f .&. complement flagProperlyPaired+  -- | Removes syntactic warts from old read names or the read names used -- in FastQ files.

@@ -1,3 +1,8 @@+-- | Parser for @FastA/FastQ@, 'Iteratee' style, based on+-- "Data.Attoparsec", and written such that it is compatible with module+-- 'Bio.Bam'.  This gives import of @FastA/FastQ@ while respecting some+-- local (to MPI EVAN) conventions.+ module Bio.Bam.Fastq ( parseFastq, parseFastq', parseFastqCassava ) where  import Bio.Bam.Header@@ -11,19 +16,13 @@ import qualified Data.ByteString.Char8              as S import qualified Data.Vector.Generic                as V --- ^ Parser for @FastA/FastQ@, 'Iteratee' style, based on--- "Data.Attoparsec", and written such that it is compatible with module--- 'Bio.Bam'.  This gives import of @FastA/FastQ@ while respecting some--- local conventions.- -- | Reader for DNA (not protein) sequences in FastA and FastQ.  We read -- everything vaguely looking like FastA or FastQ, then shoehorn it into -- a BAM record.  We strive to extract information following more or--- less established conventions from the header, but we won't support--- everything under the sun.  The recognized syntactical warts are--- converted into appropriate flags and removed.  Only the canonical--- variant of FastQ is supported (qualities stored as raw bytes with--- base 33).+-- less established conventions from the header, but don't aim for+-- completeness.  The recognized syntactical warts are converted into+-- appropriate flags and removed.  Only the canonical variant of FastQ+-- is supported (qualities stored as raw bytes with offset 33). -- -- Supported additional conventions: --@@ -36,18 +35,13 @@ -- -- * A name prefix of @T_@ flags the sequence as unpaired and trimmed ----- * A name prefix of @C_@, either before or after any of the other+-- * A name prefix of @C_@, optionally before or after any of the other --   prefixes, is turned into the extra flag @XP:i:-1@ (result of --   duplicate removal with unknown duplicate count). -- -- * A collection of tags separated from the name by an octothorpe is --   removed and put into the fields @XI@ and @XJ@ as text. ----- * In 'parseFastqCassava' only, if the first word of the description---   has at least four colon separated subfields, the first if used to---   flag first/second mate, the second is the \"QC failed\" flag, and---   the fourth is the index sequence.--- -- Everything before the first sequence header is ignored.  Headers can -- start with @\>@ or @\@@, we treat both equally.  The first word of -- the header becomes the read name, the remainder of the header is@@ -63,6 +57,13 @@ parseFastq :: Monad m => Enumeratee Bytes [ BamRec ] m a parseFastq = parseFastq' (const id) +-- | Like 'parseFastq', but also+--+-- * If the first word of the description has at least four colon+--   separated subfields, the first is used to flag first/second mate,+--   the second is the \"QC failed\" flag, and the fourth is the index+--   sequence.+ parseFastqCassava :: Monad m => Enumeratee Bytes [ BamRec ] m a parseFastqCassava = parseFastq' (pdesc . S.split ':' . S.takeWhile (' ' /=))   where@@ -78,7 +79,6 @@ -- which can modify the parsed record.  Note that the quality field can -- end up empty. -{-# WARNING parseFastq' "parseFastq' no longer removes syntactic warts!" #-} parseFastq' :: Monad m => ( Bytes -> BamRec -> BamRec ) -> Enumeratee Bytes [ BamRec ] m a parseFastq' descr it = do skipJunk ; convStream (parserToIteratee $ (:[]) <$> pRec) it   where

@@ -1,3 +1,5 @@+-- | Quality filters adapted from prehistoric pipeline.+ module Bio.Bam.Filter (     filterPairs, QualFilter,     complexSimple, complexEntropy,@@ -13,8 +15,6 @@ import Prelude  import qualified Data.Vector.Generic as V---- ^ Quality filters adapted from old pipeline.  -- | A filter/transformation applied to pairs of reads.  We supply a -- predicate to be applied to single reads and one to be applied to

@@ -317,13 +317,14 @@   -- | Compares two sequence names the way samtools does.--- samtools sorts by "strnum_cmp":--- . if both strings start with a digit, parse the initial+-- samtools sorts by \"strnum_cmp\":+--+-- * if both strings start with a digit, parse the initial --   sequence of digits and compare numerically, if equal, --   continue behind the numbers--- . else compare the first characters (possibly NUL), if equal+-- * else compare the first characters (possibly NUL), if equal --   continue behind them--- . else both strings ended and the shorter one counts as+-- * else both strings ended and the shorter one counts as --   smaller (and that part is stupid)  compareNames :: Seqid -> Seqid -> Ordering

@@ -9,8 +9,7 @@     eneeBamRefseq,     eneeBamSubseq,     eneeBamRegions,-    eneeBamUnaligned,-    subsampleBam+    eneeBamUnaligned ) where  import Bio.Bam.Header@@ -18,15 +17,12 @@ import Bio.Bam.Rec import Bio.Bam.Regions              ( Region(..), Subsequence(..) ) import Bio.Iteratee-import Bio.Iteratee.Bgzf import Bio.Prelude import System.Directory             ( doesFileExist ) import System.FilePath              ( dropExtension, takeExtension, (<.>) )-import System.Random                ( randomRIO )  import qualified Bio.Bam.Regions                as R-import qualified Control.Exception              as E-import qualified Data.IntMap                    as M+import qualified Data.IntMap.Strict             as M import qualified Data.ByteString                as B import qualified Data.Vector                    as V import qualified Data.Vector.Mutable            as W@@ -115,7 +111,7 @@     , let boff' = max boff cpt     , boff' < eoff ]   where-    !cpt = maybe 0 snd $ lookupLE beg cpts+    !cpt = maybe 0 snd $ M.lookupLE beg cpts  -- list of bins for given range of coordinates, from Heng's horrible code binList :: BamIndex a -> Int -> Int -> [Int]@@ -333,51 +329,3 @@ eneeBamRegions :: Monad m => BamIndex b -> [R.Region] -> Enumeratee [BamRaw] [BamRaw] m a eneeBamRegions bi = foldr ((>=>) . uncurry (eneeBamSubseq bi)) return . R.toList . R.fromList --lookupLE :: M.Key -> M.IntMap a -> Maybe (M.Key, a)-lookupLE k m = case ma of-    Just a              -> Just (k,a)-    Nothing | M.null m1 -> Nothing-            | otherwise -> Just $ M.findMax m1-  where (m1,ma,_) = M.splitLookup k m----- | Subsample randomly from a BAM file.  If an index exists, this--- produces an infinite stream taken from random locations in the file.------ XXX It would be cool if we could subsample from multiple BAM files.--- It's a bit annoying to code: we'd probably read the indices up front,--- estimate how many reads we'd find in each file, then open them--- recursively to form a monad stack where the merging function has to--- select randomly where to read from.  Hm.--subsampleBam :: (MonadIO m, MonadMask m) => FilePath -> Enumerator' BamMeta [BamRaw] m b-subsampleBam fp o = liftIO (E.try (readBamIndex fp)) >>= subsam-  where-    -- no index, so just stream-    subsam (Left e) = enumFile defaultBufSize fp >=> run $-                      joinI $ decompressBgzfBlocks $-                      joinI $ decodeBam $ \hdr ->-                      takeWhileE (isValidRefseq . b_rname . unpackBam) (o hdr)-                            `const` (e::E.SomeException)--    -- with index: chose random bins and read from them-    subsam (Right bix) = withFileFd fp $ \fd -> do-                         hdr <- enumFdRandom defaultBufSize fd >=> run $-                                joinI $ decompressBgzfBlocks' 1 $-                                joinI $ decodeBam return-                         go fd (o hdr)-      where-        !ckpts = U.fromList . V.foldr ((++) . M.elems) [] $ refseq_ckpoints bix--        go fd o1 = enumCheckIfDone o1 >>= go' fd--        go'  _ (True,  o2) = return o2-        go' fd (False, o2) = do i <- liftIO $ randomRIO (0, U.length ckpts -1)-                                enum fd i o2 >>= go fd--        enum fd i = enumFdRandom defaultBufSize fd               $=-                    decompressBgzfBlocks' 1                      $=-                    (\it -> do seek . fromIntegral $ ckpts U.! i-                               convStream getBamRaw it)          $=-                    takeStream 512

@@ -1,4 +1,15 @@ {-# LANGUAGE Rank2Types, DeriveGeneric #-}++-- | Pileup, similar to Samtools+--+-- Pileup turns a sorted sequence of reads into a sequence of \"piles\",+-- one for each site where a genetic variant might be called.  We will+-- scan each read's CIGAR line and MD field in concert with the sequence+-- and effective quality.  Effective quality is the lowest available+-- quality score of QUAL, MAPQ, and BQ.  For aDNA calling, a base is+-- represented as four probabilities, derived from a position dependent+-- damage model.+ module Bio.Bam.Pileup where  import Bio.Bam.Header@@ -10,63 +21,12 @@ import qualified Data.Vector.Generic    as V import qualified Data.Vector.Unboxed    as U --- ^ Genotype Calling:  like Samtools(?), but for aDNA------ The goal for this module is to call haploid and diploid single--- nucleotide variants the best way we can, including support for aDNA.--- Indel calling is out of scope, we only do it "on the side".------ The cleanest way to call genotypes under all circumstances is--- probably the /Dindel/ approach:  define candidate haplotypes, align--- each read to each haplotype, then call the likely haplotypes with a--- quality derived from the quality scores.  This approach neatly--- integrates indel calling with ancient DNA and makes a separate indel--- realigner redundant.  However, it's rather expensive in that it--- requires inclusion of an aligner, and we'd need an aligner that is--- compatible with the chosen error model, which might be hard.------ Here we'll take a short cut:  We do not really call indels.  Instead,--- these variants are collected and are assigned an affine score.  This--- works best if indels are 'left-aligned' first.  In theory, one indel--- variant could be another indel variant with a sequencing error---we--- ignore that possibility for the most part.  Once indels are taken--- care off, SNVs are treated separately as independent columns of the--- pileup.------ Regarding the error model, there's a choice between /samtools/ or the--- naive model everybody else (GATK, Rasmus Nielsen, etc.) uses.  Naive--- is easy to marry to aDNA, samtools is (probably) better.  Either way,--- we introduce a number of parameters (@eta@ and @kappa@ for--- /samtools/, @lambda@, @delta@, @delta_ss@ for /Johnson/).  Running a--- maximum likehood fit for those may be valuable.  It would be cool, if--- we could do that without rerunning the complete genotype caller, but--- it's not a priority.------ So, outline of the genotype caller:  We read BAM (minimally--- filtering; general filtering is somebody else's problem, but we might--- want to split by read group).  We will scan each read's CIGAR line in--- concert with the sequence and effective quality.  Effective quality--- is the lowest available quality score of QUAL, MAPQ, and BQ.  For--- aDNA calling, the base is transformed into four likelihoods based on--- the aDNA substitution matrix.------ So, either way, we need something like "pileup", where indel variants--- are collected as they are (any length), while matches are piled up.------ Regarding output, we certainly don't want to write VCF or BCF.  (No--- VCF because it's ugly, no BCF, because the tool support is--- non-existent.)  It will definitely be something binary.  For the GL--- values, small floating point formats may make sense: half-precision--- floating point's representable range would be 6.1E-5 to 6.5E+5, 0.4.4--- minifloat from Bio.Util goes from 0 to 63488.-- -- | The primitive pieces for genotype calling:  A position, a base -- represented as four likelihoods, an inserted sequence, and the -- length of a deleted sequence.  The logic is that we look at a base -- followed by some indel, and all those indels are combined into a -- single insertion and a single deletion.-data PrimChunks = Seek Int PrimBase                                   -- ^ skip to position (at start or after N operation)+data PrimChunks = Seek {-# UNPACK #-} !Int PrimBase                   -- ^ skip to position (at start or after N operation)                 | Indel [Nucleotides] [DamagedBase] PrimBase          -- ^ observed deletion and insertion between two bases                 | EndOfRead                                           -- ^ nothing anymore   deriving Show@@ -109,8 +69,8 @@ -- | Decomposes a BAM record into chunks suitable for piling up.  We -- pick apart the CIGAR and MD fields, and combine them with sequence -- and quality as appropriate.  Clipped bases are removed/skipped as--- appropriate.  We also do apply a substitution matrix to each base,--- which must be supplied along with the read.+-- needed.  We also apply a substitution matrix to each base, which must+-- be supplied along with the read. {-# INLINE decompose #-} decompose :: DmgToken -> BamRaw -> [PosPrimChunks] decompose dtok br =@@ -118,7 +78,7 @@     then [] else [(b_rname, b_pos, isReversed b, pchunks)]   where     b@BamRec{..} = unpackBam br-    pchunks = firstBase b_pos 0 0 (maybe [] id $ getMd b)+    pchunks = firstBase b_pos 0 0 (fromMaybe [] $ getMd b)      !max_cig = V.length b_cigar     !max_seq = V.length b_seq@@ -129,27 +89,14 @@     -- and BAQ.  If QUAL is invalid, we replace it (arbitrarily) with     -- 23 (assuming a rather conservative error rate of ~0.5%), BAQ is     -- added to QUAL, and MAPQ is an upper limit for effective quality.-    get_seq :: Int -> Nucleotides -> DamagedBase-    get_seq i = case b_seq `V.unsafeIndex` i of                                 -- nucleotide-            n | n == nucsA -> DB nucA qe dtok dmg-              | n == nucsC -> DB nucC qe dtok dmg-              | n == nucsG -> DB nucG qe dtok dmg-              | n == nucsT -> DB nucT qe dtok dmg-              | otherwise  -> DB nucA (Q 0) dtok dmg-      where-        !q   = case b_qual `V.unsafeIndex` i of Q 0xff -> Q 30 ; x -> x         -- quality; invalid (0xff) becomes 30-        !q'  | i >= B.length baq = q                                            -- no BAQ available-             | otherwise = Q (unQ q + (B.index baq i - 64))                     -- else correct for BAQ-        !qe  = min q' b_mapq                                                    -- use MAPQ as upper limit-        !dmg = if i+i > max_seq then i-max_seq else i -    get_seq' :: Int -> DamagedBase-    get_seq' i = case b_seq `V.unsafeIndex` i of                                -- nucleotide-            n | n == nucsA -> DB nucA qe dtok dmg nucsA-              | n == nucsC -> DB nucC qe dtok dmg nucsC-              | n == nucsG -> DB nucG qe dtok dmg nucsG-              | n == nucsT -> DB nucT qe dtok dmg nucsT-              | otherwise  -> DB nucA (Q 0) dtok dmg n+    get_seq :: Int -> (Nucleotides -> Nucleotides) -> DamagedBase+    get_seq i f = case b_seq `V.unsafeIndex` i of                                -- nucleotide+            n | n == nucsA -> DB nucA qe dtok dmg (f n)+              | n == nucsC -> DB nucC qe dtok dmg (f n)+              | n == nucsG -> DB nucG qe dtok dmg (f n)+              | n == nucsT -> DB nucT qe dtok dmg (f n)+              | otherwise  -> DB nucA (Q 0) dtok dmg (f n)       where         !q   = case b_qual `V.unsafeIndex` i of Q 0xff -> Q 30 ; x -> x         -- quality; invalid (0xff) becomes 30         !q'  | i >= B.length baq = q                                            -- no BAQ available@@ -176,7 +123,7 @@       where         -- We have to treat (MdNum 0), because samtools actually         -- generates(!) it all over the place and if not handled as a-        -- special case, it looks like an incinsistend MD field.+        -- special case, it looks like an inconsistent MD field.         drop_del n (MdDel ns : mds')             | n < length ns = MdDel (drop n ns) : mds'             | n > length ns = drop_del (n - length ns) mds'@@ -194,11 +141,11 @@     nextBase !wt !pos !is !ic !io mds = case mds of         MdNum   0 : mds' -> nextBase wt pos is ic io mds'         MdDel  [] : mds' -> nextBase wt pos is ic io mds'-        MdNum   1 : mds' -> nextBase' (get_seq' is      ) mds'-        MdNum   n : mds' -> nextBase' (get_seq' is      ) (MdNum (n-1) : mds')-        MdRep ref : mds' -> nextBase' (get_seq  is ref  ) mds'-        MdDel   _ : _    -> nextBase' (get_seq  is nucsN) mds-        [              ] -> nextBase' (get_seq  is nucsN) [ ]+        MdNum   1 : mds' -> nextBase' (get_seq is id   ) mds'+        MdNum   n : mds' -> nextBase' (get_seq is id   ) (MdNum (n-1) : mds')+        MdRep ref : mds' -> nextBase' (get_seq is $ const ref  ) mds'+        MdDel   _ : _    -> nextBase' (get_seq is $ const nucsN) mds+        [              ] -> nextBase' (get_seq is $ const nucsN) [ ]       where         nextBase' ref mds' = Base wt ref b_mapq $ nextIndel  [] [] (pos+1) (is+1) ic (io+1) mds' @@ -223,11 +170,9 @@             Mat :* cl | io == cl  -> nextIndel  ins     del   pos     is  (ic+1) 0 mds                       | otherwise -> indel del out $ nextBase (length del) pos is ic io mds -- ends up generating a 'Base'       where-        indel d o k = rlist o `seq` Indel d o k+        indel d o k = foldr seq (Indel d o k) o         out    = concat $ reverse ins-        isq cl = [ get_seq i gap | i <- [is..is+cl-1] ] : ins-        rlist [    ] = ()-        rlist (a:as) = a `seq` rlist as+        isq cl = [ get_seq i $ const gap | i <- [is..is+cl-1] ] : ins          -- We have to treat (MdNum 0), because samtools actually         -- generates(!) it all over the place and if not handled as a@@ -243,10 +188,11 @@ -- fitlering (so not very useful), but we keep them because it's easy to -- track them. -data CallStats = CallStats { read_depth       :: {-# UNPACK #-} !Int       -- number of contributing reads-                           , reads_mapq0      :: {-# UNPACK #-} !Int       -- number of (non-)contributing reads with MAPQ==0-                           , sum_mapq         :: {-# UNPACK #-} !Int       -- sum of map qualities of contributing reads-                           , sum_mapq_squared :: {-# UNPACK #-} !Int }     -- sum of squared map qualities of contributing reads+data CallStats = CallStats+    { read_depth       :: {-# UNPACK #-} !Int       -- number of contributing reads+    , reads_mapq0      :: {-# UNPACK #-} !Int       -- number of (non-)contributing reads with MAPQ==0+    , sum_mapq         :: {-# UNPACK #-} !Int       -- sum of map qualities of contributing reads+    , sum_mapq_squared :: {-# UNPACK #-} !Int }     -- sum of squared map qualities of contributing reads   deriving (Show, Eq, Generic)  instance Monoid CallStats where@@ -259,22 +205,6 @@                             , sum_mapq         = sum_mapq x + sum_mapq y                             , sum_mapq_squared = sum_mapq_squared x + sum_mapq_squared y } --- | Genotype likelihood values.  A variant call consists of a position,--- some measure of qualities, genotype likelihood values, and a--- representation of variants.  A note about the GL values:  @VCF@ would--- normalize them so that the smallest one becomes zero.  We do not do--- that here, since we might want to compare raw values for a model--- test.  We also store them in a 'Double' to make arithmetics easier.--- Normalization is appropriate when converting to @VCF@.------ If GL is given, we follow the same order used in VCF:--- \"the ordering of genotypes for the likelihoods is given by:--- F(j/k) = (k*(k+1)/2)+j.  In other words, for biallelic sites the--- ordering is: AA,AB,BB; for triallelic sites the ordering is:--- AA,AB,BB,AC,BC,CC, etc.\"--type GL = U.Vector Prob- newtype V_Nuc  = V_Nuc  (U.Vector Nucleotide)  deriving (Eq, Ord, Show) newtype V_Nucs = V_Nucs (U.Vector Nucleotides) deriving (Eq, Ord, Show) @@ -292,10 +222,10 @@ type IndelPile = [( Qual, ([Nucleotides], [DamagedBase]) )]   -- a list of indel variants  -- | Running pileup results in a series of piles.  A 'Pile' has the--- basic statistics of a 'VarCall', but no GL values and a pristine list--- of variants instead of a proper call.  We emit one pile with two--- 'BasePile's (one for each strand) and one 'IndelPile' (the one--- immediately following) at a time.+-- basic statistics of a 'VarCall', but no likelihood values and a+-- pristine list of variants instead of a proper call.  We emit one pile+-- with two 'BasePile's (one for each strand) and one 'IndelPile' (the+-- one immediately following) at a time.  data Pile' a b = Pile { p_refseq     :: {-# UNPACK #-} !Refseq                       , p_pos        :: {-# UNPACK #-} !Int@@ -308,27 +238,6 @@ -- | Raw pile.  Bases and indels are piled separately on forward and -- backward strands. type Pile = Pile' (BasePile, BasePile) (IndelPile, IndelPile)---- | Simple single population model.  'prob_div' is the fraction of--- homozygous divergent sites, 'prob_het' is the fraction of--- heterozygous variant sites among sites that are not homozygous--- divergent.-data SinglePop = SinglePop { prob_div :: !Double, prob_het :: !Double }---- | Computes posterior  genotype probabilities from likelihoods under--- the 'SinglePop' model.-{-# INLINE single_pop_posterior #-}-single_pop_posterior :: ( U.Unbox a, Ord a, Floating a )-                     => SinglePop -> Int -> U.Vector (Prob' a) -> U.Vector (Prob' a)-single_pop_posterior SinglePop{..} refix lks = U.map (/ U.sum v) v-  where-    priors = U.replicate (U.length lks)-                         ((1/3) * prob_het  * (1-prob_div))                                 -- hets-                U.// [ (refix, (1-prob_het) * (1-prob_div)) ]                               -- ref-                U.// [ (i,               (1/3) * prob_div ) | i <- hixes, i /= refix ]      -- homs--    v = U.zipWith (\l p -> l * toProb (realToFrac p)) lks priors-    hixes = takeWhile (< U.length lks) $ scanl (+) 0 [2..]  -- | The pileup enumeratee takes 'BamRaw's, decomposes them, interleaves -- the pieces appropriately, and generates 'Pile's.  The output will

@@ -1,3 +1,12 @@+-- | Parsers for BAM and SAM.+--+-- TONOTDO:+--+-- * Reader for gzipped\/bzipped\/bgzf'ed SAM.  Storing SAM is a bad idea,+--   so why would anyone ever want to compress, much less index it?+-- * Proper support for the "=" symbol.  It's completely alien to the+--   ususal representation of sequences.+ module Bio.Bam.Reader (     Block(..),     decompressBgzfBlocks,@@ -52,21 +61,13 @@ import qualified Data.Vector.Storable               as VS import qualified Data.Vector.Unboxed                as U --- ^ Parsers for BAM and SAM.  We employ an @Iteratee@ interface, and we--- strive to support everything possible in BAM.  The implementation of--- nucleotides is somewhat lacking:  the "=" symbol is not understood.------ TONOTDO:--- - Reader for gzipped/bzipped/bgzf'ed SAM.  Storing SAM is a bad idea,---   so why would anyone ever want to compress, much less index it?- type BamrawEnumeratee m b = Enumeratee' BamMeta Bytes [BamRaw] m b type BamEnumeratee m b = Enumeratee' BamMeta Bytes [BamRec] m b  isBamOrSam :: MonadIO m => Iteratee Bytes m (BamEnumeratee m a) isBamOrSam = maybe decodeSam wrap `liftM` isBam   where-    wrap enee it' = enee (\hdr -> mapStream unpackBam (it' hdr)) >>= lift . run+    wrap enee it' = enee (mapStream unpackBam . it') >>= lift . run   -- | Checks if a file contains BAM in any of the common forms,@@ -94,16 +95,17 @@     decodeSamLoop (meta_refs meta) (inner meta)  decodeSamLoop :: Monad m => Refs -> Enumeratee [Bytes] [BamRec] m a-decodeSamLoop refs inner = convStream (liftI parse_record) inner-  where !refs' = M.fromList $ zip [ nm | BamSQ { sq_name = nm } <- F.toList refs ] [toEnum 0..]-        ref x = M.lookupDefault invalidRefseq x refs'+decodeSamLoop refs = convStream (liftI parse_record)+  where+    !refs' = M.fromList $ zip [ nm | BamSQ { sq_name = nm } <- F.toList refs ] [toEnum 0..]+    ref x = M.lookupDefault invalidRefseq x refs' -        parse_record (EOF x) = icont parse_record x-        parse_record (Chunk []) = liftI parse_record-        parse_record (Chunk (l:ls)) | "@" `S.isPrefixOf` l = parse_record (Chunk ls)-        parse_record (Chunk (l:ls)) = case P.parseOnly (parseSamRec ref) l of-            Right  r -> idone [r] (Chunk ls)-            Left err -> icont parse_record (Just $ iterStrExc $ err ++ ", " ++ show l)+    parse_record (EOF x) = icont parse_record x+    parse_record (Chunk []) = liftI parse_record+    parse_record (Chunk (l:ls)) | "@" `S.isPrefixOf` l = parse_record (Chunk ls)+    parse_record (Chunk (l:ls)) = case P.parseOnly (parseSamRec ref) l of+        Right  r -> idone [r] (Chunk ls)+        Left err -> icont parse_record (Just $ iterStrExc $ err ++ ", " ++ show l)  -- | Parser for SAM that doesn't look for a header.  Has the advantage -- that it doesn't stall on a pipe that never delivers data.  Has the@@ -119,7 +121,7 @@                   <*> snum <*> sequ <*> quals <*> exts <*> pure 0   where     sep      = P.endOfInput <|> () <$ P.char '\t'-    word     = P.takeTill ((==) '\t') <* sep+    word     = P.takeTill ('\t' ==) <* sep     num      = P.decimal <* sep     num'     = P.decimal <* sep     snum     = P.signed P.decimal <* sep@@ -143,7 +145,7 @@      value    = P.char 'A' *> P.char ':' *> (Char <$>               anyWord8) <|>                P.char 'i' *> P.char ':' *> (Int  <$>     P.signed P.decimal) <|>-               P.char 'Z' *> P.char ':' *> (Text <$> P.takeTill ((==) '\t')) <|>+               P.char 'Z' *> P.char ':' *> (Text <$>   P.takeTill ('\t' ==)) <|>                P.char 'H' *> P.char ':' *> (Bin  <$>               hexarray) <|>                P.char 'f' *> P.char ':' *> (Float . realToFrac <$> P.double) <|>                P.char 'B' *> P.char ':' *> (@@ -186,7 +188,7 @@ -- @Iteratee@.) isBgzfBam  = do b <- isBgzf                 k <- if b then joinI $ enumInflate GZip defaultDecompressParams isPlainBam else return Nothing-                return $ (\_ -> (joinI . decompressBgzfBlocks . decodeBam)) `fmap` k+                return $ const (joinI . decompressBgzfBlocks . decodeBam) `fmap` k  isGzipBam  = do b <- isGzip                 k <- if b then joinI $ enumInflate GZip defaultDecompressParams isPlainBam else return Nothing@@ -212,14 +214,12 @@ concatDefaultInputs it0 = liftIO getArgs >>= \fs -> concatInputs fs it0  concatInputs :: (MonadIO m, MonadMask m) => [FilePath] -> Enumerator' BamMeta [BamRaw] m a-concatInputs [        ] = \k -> enumFd defaultBufSize stdInput (decodeAnyBam k) >>= run-concatInputs (fp0:fps0) = \k -> enum1 fp0 k >>= go fps0+concatInputs [       ] = enumFd defaultBufSize stdInput . decodeAnyBam >=> run+concatInputs (fp0:fps) = enum1 fp0 >=> go fps   where-    enum1 "-" k1 = enumFd   defaultBufSize stdInput (decodeAnyBam k1) >>= run-    enum1  fp k1 = enumFile defaultBufSize fp       (decodeAnyBam k1) >>= run--    go [       ] = return-    go (fp1:fps) = enum1 fp1 . const >=> go fps+    go = foldr (\fp1 -> (>=>) (enum1 fp1 . const)) return+    enum1 "-" = enumFd   defaultBufSize stdInput . decodeAnyBam >=> run+    enum1  fp = enumFile defaultBufSize fp       . decodeAnyBam >=> run  mergeDefaultInputs :: (MonadIO m, MonadMask m)     => (BamMeta -> Enumeratee [BamRaw] [BamRaw] (Iteratee [BamRaw] m) a)@@ -270,7 +270,7 @@                                    nm <- endianRead4 LSB >>= iGetString . fromIntegral                                    ln <- endianRead4 LSB                                    return $! acc |> BamSQ (S.init nm) (fromIntegral ln) []-                             ) Z.empty $ [1..nref]+                             ) Z.empty [1..nref]      -- Need to merge information from header into actual reference list.     -- The latter is the authoritative source for the *order* of the@@ -282,7 +282,7 @@         in meta { meta_refs = fmap (\s -> maybe s (mmerge' s) (M.lookup (sq_name s) tbl)) refs }      mmerge' l r | sq_length l == sq_length r = l { sq_other_shit = sq_other_shit l ++ sq_other_shit r }-               | otherwise                  = l -- contradiction in header, but we'll just ignore it+                | otherwise                  = l -- contradiction in header, but we'll just ignore it   {-# INLINE getBamRaw #-}

@@ -6,13 +6,6 @@ -- do not have support for ambiguity codes, and the "=" symbol is not -- understood. --- TODO:--- - Automatic creation of some kind of index.  If possible, this should---   be the standard index for sorted BAM and/or the newer CSI format.---   Optionally, a block index for slicing of large files, even unsorted---   ones.  Maybe an index by name and an index for group-sorted files.---   Sensible indices should be generated whenever a file is written.- module Bio.Bam.Rec (     BamRaw,     bamRaw,@@ -103,7 +96,7 @@         where             w = fromEnum op .|. shiftL num 4 --- | extracts the aligned length from a cigar line+-- | Extracts the aligned length from a cigar line. -- This gives the length of an alignment as measured on the reference, -- which is different from the length on the query or the length of the -- alignment.

@@ -1,10 +1,10 @@ module Bio.Bam.Regions where  import Bio.Bam.Header ( Refseq(..) )-import Data.List ( foldl' )-import qualified Data.IntMap as IM-import Prelude+import Bio.Prelude +import qualified Data.IntMap.Strict as IM+ data Region = Region { refseq :: !Refseq, start :: !Int, end :: !Int }   deriving (Eq, Ord, Show) @@ -32,7 +32,7 @@ addInt :: Int -> Int -> Subsequence -> Subsequence addInt b e (Subsequence m0) = Subsequence $ merge_into b e m0   where-    merge_into x y m = case lookupLT y m of+    merge_into x y m = case IM.lookupLT y m of         Just (u,v) | x < u && y <= v -> merge_into x v $ IM.delete u m    -- extend to the left                    | x < u           -> merge_into x y $ IM.delete u m    -- subsume                    | y <= v          -> m                                 -- subsumed@@ -40,11 +40,6 @@         _                            -> IM.insert  x y m                  -- no overlap  overlaps :: Int -> Int -> Subsequence -> Bool-overlaps b e (Subsequence m) = case lookupLT e m of+overlaps b e (Subsequence m) = case IM.lookupLT e m of         Just (_,v) -> b < v         Nothing    -> False--lookupLT :: IM.Key -> IM.IntMap a -> Maybe (IM.Key, a)-lookupLT k m | IM.null m1 = Nothing-             | otherwise  = Just $ IM.findMax m1-  where (m1,_) = IM.split k m

@@ -12,7 +12,7 @@  import qualified Data.ByteString        as B import qualified Data.ByteString.Char8  as T-import qualified Data.Map               as M+import qualified Data.Map.Strict        as M import qualified Data.Vector.Generic    as V import qualified Data.Vector.Storable   as VS import qualified Data.Vector.Unboxed    as U@@ -115,10 +115,10 @@     same_pos u v = b_cpos u == b_cpos v     b_cpos u = (b_rname u, b_pos u) -    nice_sort x = sortBy (comparing (V.length . b_seq . snd)) x+    nice_sort = sortBy $ comparing (V.length . b_seq . snd)      mapGroups f o = tryHead >>= maybe (return o) (\a -> eneeCheckIfDone (mg1 f a []) o)-    mg1 f a acc k = tryHead >>= \mb -> case mb of+    mg1 f a acc k = tryHead >>= \case                         Nothing -> return . k . Chunk . f $ a:acc                         Just b | same_pos a b -> mg1 f a (b:acc) k                                | otherwise -> eneeCheckIfDone (mg1 f b []) . k . Chunk . f $ a:acc@@ -209,7 +209,7 @@         (half_unaligned, half_aligned) = partition isUnmapped raw_half_pairs          mkMap :: Ord a => (BamRec -> a) -> [BamRec] -> (M.Map a (Int,Decision), [BamRec])-        mkMap f x = let m1 = M.map (\xs -> (length xs, collapse xs)) $ accumMap f id x+        mkMap f x = let m1 = M.map (length &&& collapse) $ accumMap f id x                     in (M.map (second fst) m1, concatMap (snd.snd) $ M.elems m1)          (pairs',r1)   = mkMap (\b -> (b_mate_pos b,   b_strand b, b_mate b)) pairs@@ -257,7 +257,7 @@          in case M.lookup str singles of             Nothing                    -> (              (m,v) : r,     unal )-            Just (n, Consensus      w) -> ( (n, add_xp_of w v) : r,     unal )      -- XXX do we need this w?!+            Just (n, Consensus      w) -> ( (n, add_xp_of w v) : r,     unal )             Just (n, Representative w) -> ( (n, add_xp_of w v) : r, w : unal )      -- No more pairs, delegate the problem@@ -412,11 +412,11 @@     xa' = nub' [ T.split ';' xas | Just (Text xas) <- map (lookup "XA" . b_exts) brs ]      modify_extensions es = foldr ($!) es $-        [ let vs = [ v | Just v <- map (lookup k . b_exts) brs ]+        [ let vs = mapMaybe (lookup k . b_exts) brs           in if null vs then id else updateE k $! maj vs | k <- do_maj ] ++         [ let vs = [ v | Just (Int v) <- map (lookup k . b_exts) brs ]           in if null vs then id else updateE k $! Int (rmsq vs) | k <- do_rmsq ] ++-        [ deleteE k | k <- useless ] +++        map deleteE useless ++         [ updateE "NM" $! Int nm'         , updateE "XP" $! Int (foldl' (\a b -> a `oplus` extAsInt 1 "XP" b) 0 brs)         , if null xa' then id else updateE "XA" $! (Text $ T.intercalate (T.singleton ';') xa')@@ -528,7 +528,7 @@   where     overhangs = not (isUnmapped b) && b_pos b < l && l < b_pos b + alignedLength (b_cigar b) -    do_wrap = case split_ecig (l - b_pos b) $ toECig (b_cigar b) (maybe [] id $ getMd b) of+    do_wrap = case split_ecig (l - b_pos b) $ toECig (b_cigar b) (fromMaybe [] $ getMd b) of                   (left,right) -> [ Right $ b { b_cigar = toCigar  left }            `setMD` left                                   , Left  $ b { b_cigar = toCigar right, b_pos = 0 } `setMD` right ] @@ -569,12 +569,10 @@ mask_all (Ins' n ec) = SMa' n $ mask_all ec mask_all (SMa' n ec) = SMa' n $ mask_all ec --- | Argh, this business with the CIGAR operations is a mess, it gets--- worse when combined with MD.  Okay, we will support CIGAR (no "=" and--- "X" operations) and MD.  If we have MD on input, we generate it on--- output, too.  And in between, we break everything into /very small/--- operations.  (Yes, the two terminating constructors are a weird--- hack.)+-- | Extended CIGAR.  This subsumes both the CIGAR string and the+-- optional MD field.  If we have MD on input, we generate it on output,+-- too.  And in between, we break everything into /very small/+-- operations.  data ECig = WithMD                      -- terminate, do generate MD field           | WithoutMD                   -- terminate, don't bother with MD@@ -589,7 +587,7 @@   toECig :: VS.Vector Cigar -> [MdOp] -> ECig-toECig cig md = go (VS.toList cig) md+toECig = go . VS.toList   where     go        cs  (MdNum  0:mds) = go cs mds     go        cs  (MdDel []:mds) = go cs mds

@@ -30,7 +30,7 @@ import qualified Data.Vector.Storable                   as W  -- | Trims from the 3' end of a sequence.--- @trim_3\' p b@ trims the 3' end of the sequence in @b@ at the+-- @trim_3' p b@ trims the 3' end of the sequence in @b@ at the -- earliest position such that @p@ evaluates to true on every suffix -- that was trimmed off.  Note that the 3' end may be the beginning of -- the sequence if it happens to be stored in reverse-complemented form.@@ -229,7 +229,7 @@ merged_seq :: (V.Vector v Nucleotides, V.Vector v Qual)            => Int -> v Nucleotides -> v Qual -> v Nucleotides -> v Qual -> v Nucleotides merged_seq l b_seq_r1 b_qual_r1 b_seq_r2 b_qual_r2 = V.concat-        [ V.take (l - len_r2) (b_seq_r1)+        [ V.take (l - len_r2) b_seq_r1         , V.zipWith4 zz          (V.take l $ V.drop (l - len_r2) b_seq_r1)                                  (V.take l $ V.drop (l - len_r2) b_qual_r1)                      (V.reverse $ V.take l $ V.drop (l - len_r1) b_seq_r2)@@ -246,7 +246,7 @@ merged_qual :: (V.Vector v Nucleotides, V.Vector v Qual)             => Word8 -> Int -> v Nucleotides -> v Qual -> v Nucleotides -> v Qual -> v Qual merged_qual qmax l b_seq_r1 b_qual_r1 b_seq_r2 b_qual_r2 = V.concat-        [ V.take (l - len_r2) (b_qual_r1)+        [ V.take (l - len_r2) b_qual_r1         , V.zipWith4 zz           (V.take l $ V.drop (l - len_r2) b_seq_r1)                                   (V.take l $ V.drop (l - len_r2) b_qual_r1)                       (V.reverse $ V.take l $ V.drop (l - len_r1) b_seq_r2)

@@ -1,3 +1,6 @@+-- | Printers for BAM and SAM.  BAM is properly supported, SAM can be+-- piped to standard output.+ module Bio.Bam.Writer (     IsBamRec(..),     encodeBamWith,@@ -29,12 +32,6 @@ import qualified Data.Vector.Unboxed                as U import qualified Data.Sequence                      as Z --- ^ Printers for BAM.  We employ an @Iteratee@ interface, and we strive--- to keep BAM records in their encoded form.  This is most compact and--- often faster, since it saves the time for repeated decoding and--- encoding, if that's not strictly needed.-- -- | write in SAM format to stdout -- This is useful for piping to other tools (say, AWK scripts) or for -- debugging.  No convenience function to send SAM to a file exists,@@ -118,11 +115,7 @@ pushBamRaw :: BamRaw -> BgzfTokens -> BgzfTokens pushBamRaw = TkLnString . raw_data --- | writes BAM encoded stuff to a file--- XXX This should(!) write indexes on the side---a simple block index--- for MapReduce style slicing, a standard BAM index or a name index--- would be possible.  When writing to a file, this makes even more--- sense than when writing to a @Handle@.+-- | Writes BAM encoded stuff to a file. writeBamFile :: IsBamRec r => FilePath -> BamMeta -> Iteratee [r] IO () writeBamFile fp meta =     C.bracket (liftIO $ openBinaryFile fp WriteMode)@@ -135,12 +128,7 @@ pipeBamOutput :: IsBamRec r => BamMeta -> Iteratee [r] IO () pipeBamOutput meta = encodeBamWith 0 meta =$ mapChunksM_ (liftIO . S.hPut stdout) --- | writes BAM encoded stuff to a @Handle@--- We generate BAM with dynamic blocks, then stream them out to the file.------ XXX This could write indexes on the side---a simple block index--- for MapReduce style slicing, a standard BAM index or a name index--- would be possible.+-- | Writes BAM encoded stuff to a 'Handle'. writeBamHandle :: (MonadIO m, IsBamRec r) => Handle -> BamMeta -> Iteratee [r] m () writeBamHandle hdl meta = encodeBamWith 6 meta =$ mapChunksM_ (liftIO . S.hPut hdl) 

@@ -52,11 +52,6 @@ import qualified Data.ByteString.Char8 as S import qualified Data.Vector.Unboxed   as U -#if __GLASGOW_HASKELL__ == 704-import Data.Vector.Generic          ( Vector(..) )-import Data.Vector.Generic.Mutable  ( MVector(..) )-#endif- -- | A nucleotide base.  We only represent A,C,G,T.  The contained -- 'Word8' ist guaranteed to be 0..3. newtype Nucleotide = N { unN :: Word8 } deriving ( Eq, Ord, Enum, Ix, Storable )@@ -299,7 +294,7 @@      readList ('-':cs) = readList cs     readList (c:cs) | isSpace c = readList cs-                    | otherwise = case readsPrec 0 (c:cs) of+                    | otherwise = case reads (c:cs) of                             [] -> [ ([],c:cs) ]                             xs -> [ (n:ns,r2) | (n,r1) <- xs, (ns,r2) <- readList r1 ]     readList [] = [([],[])]

@@ -1,6 +1,3 @@--- | Basically a reexport of "Data.Iteratee" less the names that clash--- with "Prelude" plus a handful of utilities.- module Bio.Iteratee (     iGetString,     iterLoop,@@ -102,19 +99,16 @@                         else it a >>= iterLoop it infixl 1 $== {-# INLINE ($==) #-}--- | Compose an 'Enumerator\'' with an 'Enumeratee', giving a new--- 'Enumerator\''.+-- | Compose an 'Enumerator'' with an 'Enumeratee', giving a new+-- 'Enumerator''. ($==) :: Monad m => Enumerator' hdr input m (Iteratee output m result)                  -> Enumeratee      input             output m result                  -> Enumerator' hdr                   output m result ($==) enum enee iter = run =<< enum (enee . iter) --- | Merge two 'Enumerator\''s into one.  The header provided by the--- inner 'Enumerator\'' is passed to the output iterator, the header--- provided by the outer 'Enumerator\'' is passed to the merging iteratee------ XXX  Something about those headers is unsatisfactory... there should---      be an unobtrusive way to combine headers.+-- | Merge two 'Enumerator''s into one.  The header provided by the+-- inner 'Enumerator'' is passed to the output iterator, the header+-- provided by the outer 'Enumerator'' is passed to the merging iteratee  {-# INLINE mergeEnums' #-} mergeEnums' :: (Nullable s2, Nullable s1, Monad m)@@ -259,9 +253,9 @@ pushQ a (QQ l  f b) = QQ (l+1) f (a:b)  popQ :: QQ a -> Maybe (a, QQ a)-popQ (QQ l (a:[]) b) = Just (a, QQ (l-1) (reverse b) [])-popQ (QQ l (a:fs) b) = Just (a, QQ (l-1) fs b) popQ (QQ _ [    ] _) = Nothing+popQ (QQ l [ a  ] b) = Just (a, QQ (l-1) (reverse b) [])+popQ (QQ l (a:fs) b) = Just (a, QQ (l-1) fs b)  cancelAll :: MonadIO m => QQ (Async a) -> m () cancelAll (QQ _ ff bb) = liftIO $ mapM_ cancel (ff ++ bb)@@ -302,7 +296,7 @@     go mv i         | i == n    = return n         | otherwise =-            tryHead >>= \x -> case x of+            tryHead >>= \case                 Nothing -> return i                 Just  a -> liftIO (VM.write mv i a) >> go mv (i+1) @@ -310,7 +304,7 @@ stream2vector :: (MonadIO m, VG.Vector v a) => Iteratee [a] m (v a) stream2vector = liftIO (VM.new 1024) >>= go 0   where-    go !i !mv = tryHead >>= \x -> case x of+    go !i !mv = tryHead >>= \case                   Nothing -> liftIO $ VG.unsafeFreeze $ VM.take i mv                   Just  a -> do mv' <- if VM.length mv == i then liftIO (VM.grow mv (VM.length mv)) else return mv                                 when (i `rem` 0x10000 == 0) $ liftIO performGC@@ -318,7 +312,7 @@                                 go (i+1) mv'  withFileFd :: (MonadIO m, MonadMask m) => FilePath -> (Fd -> m a) -> m a-withFileFd filepath iter = CMC.bracket+withFileFd filepath = CMC.bracket     (liftIO $ openFd filepath ReadOnly Nothing defaultFileFlags)-    (liftIO . closeFd) iter+    (liftIO . closeFd) 

@@ -15,7 +15,6 @@   -- ** Control functions   ,run   ,tryRun-  ,mapIteratee   ,ilift   ,ifold   -- ** Creating Iteratees@@ -168,9 +167,6 @@ instance NullPoint s => MonadTrans (Iteratee s) where   lift m = Iteratee $ \onDone _ -> m >>= flip onDone (Chunk emptyP) --- instance (MonadBase b m, Nullable s, NullPoint s) => MonadBase b (Iteratee s m) where-  -- liftBase = lift . liftBase- instance (MonadIO m, Nullable s, NullPoint s) => MonadIO (Iteratee s m) where   liftIO = lift . liftIO @@ -216,14 +212,6 @@     onCont' _ Nothing  = return $ maybeExc (toException EofException)     onCont' _ (Just e) = return $ maybeExc e     maybeExc e = maybe (Left (E.throw e)) Left (fromException e)---- |Transform a computation inside an @Iteratee@.-mapIteratee :: (NullPoint s, Monad n, Monad m) =>-  (m a -> n b)-  -> Iteratee s m a-  -> Iteratee s n b-mapIteratee f = lift . f . run-{-# DEPRECATED mapIteratee "This function will be removed, compare to 'ilift'" #-}  -- | Lift a computation in the inner monad of an iteratee. --

@@ -17,8 +17,6 @@ import Foreign.C.String                     ( withCAString ) import Foreign.C.Types                      ( CInt(..), CChar(..), CUInt(..), CULong(..) ) import Foreign.Marshal.Alloc                ( mallocBytes, free, allocaBytes )-import Foreign.Ptr                          ( nullPtr, castPtr, Ptr, plusPtr, minusPtr )-import Foreign.Storable                     ( peekByteOff, pokeByteOff )  import qualified Data.ByteString            as S import qualified Data.ByteString.Unsafe     as S

@@ -1,12 +1,8 @@ {-# LANGUAGE ForeignFunctionInterface #-}--- | Buffer builder to assemble Bgzf blocks.  The plan is to serialize+-- | Buffer builder to assemble Bgzf blocks.  The idea is to serialize -- stuff (BAM and BCF) into a buffer, then Bgzf chunks from the buffer. -- We use a large buffer, and we always make sure there is plenty of--- space in it (to avoid redundant checks).  Whenever a block is ready--- to be compressed, we stick it into a MVar.  When we run out of space,--- we simply use a new buffer.  Multiple threads grab pieces from the--- MVar, compress them, pass them downstream through another MVar.  A--- final thread restores the order and writes the blocks.+-- space in it (to avoid redundant checks).  module Bio.Iteratee.Builder (     BB(..),@@ -94,12 +90,12 @@     withForeignPtr arr1 $ \d ->         withForeignPtr (buffer b) $ \s ->              copyBytes d (plusPtr s (off b)) (used b - off b)-    return $ BB { buffer = arr1-                , size   = sz'-                , off    = 0-                , used   = used b - off b-                , mark   = if mark  b == maxBound then maxBound else mark  b - off b-                , mark2  = if mark2 b == maxBound then maxBound else mark2 b - off b }+    return BB{ buffer = arr1+             , size   = sz'+             , off    = 0+             , used   = used b - off b+             , mark   = if mark  b == maxBound then maxBound else mark  b - off b+             , mark2  = if mark2 b == maxBound then maxBound else mark2 b - off b }  compressChunk' :: Int -> ForeignPtr Word8 -> Int -> Int -> IO B.ByteString compressChunk' lv fptr off len =@@ -134,7 +130,7 @@             case tk of TkEnd -> liftI $ go bb k                        _     -> go1 bb k tk -        | otherwise = do+        | otherwise =             eneeCheckIfDone (flush_blocks tk bb { off = off bb + maxBlockSize }) $                 k $ Chunk [compressChunk' lv (buffer bb) (off bb) maxBlockSize] @@ -149,7 +145,7 @@ fillBuffer :: BB -> BgzfTokens -> IO (BB, BgzfTokens) fillBuffer bb0 tk = withForeignPtr (buffer bb0) (\p -> go_slowish p bb0 tk)   where-    go_slowish p bb tk1 = go_fast p bb (used bb) tk1+    go_slowish p bb = go_fast p bb (used bb)      go_fast p bb use tk1 = case tk1 of         -- no space?  not our job.@@ -178,8 +174,7 @@          TkString   s tk'             -- Too big, can't handle.  We will get progressively bigger-            -- buffers and eventually handle it; for very large strings,-            -- it works, but isn't ideal.  XXX+            -- buffers and eventually handle it.             | B.length s > size bb - use -> return (bb { used = use },tk')              | otherwise  -> do let ln = B.length s@@ -192,8 +187,7 @@          TkLnString s tk'             -- Too big, can't handle.  We will get progressively bigger-            -- buffers and eventually handle it; for very large strings,-            -- it works, but isn't ideal.  XXX+            -- buffers and eventually handle it.             | B.length s > size bb - use - 4 -> return (bb { used = use },tk')              | otherwise  -> do let ln = B.length s

@@ -226,7 +226,7 @@ -- This is provided as a higher-performance alternative to enumWords, and -- is equivalent to treating the stream as a Data.ByteString.Char8.ByteString. enumWordsBS :: Monad m => Enumeratee Bytes [Bytes] m a-enumWordsBS iter = convStream getter iter+enumWordsBS = convStream getter   where     getter = liftI step     lChar = isSpace . C.last

@@ -118,7 +118,7 @@   fromException = enumExceptionFromException  -- |Create an enumerator exception from a @String@.-data EnumStringException = EnumStringException String+newtype EnumStringException = EnumStringException String   deriving (Show, Typeable)  instance Exception EnumStringException where@@ -130,7 +130,7 @@ enStrExc = EnumException . EnumStringException  -- |The enumerator received an 'IterException' it could not handle.-data EnumUnhandledIterException = EnumUnhandledIterException IterException+newtype EnumUnhandledIterException = EnumUnhandledIterException IterException   deriving (Show, Typeable)  instance Exception EnumUnhandledIterException where@@ -177,7 +177,7 @@   fromIterException = Just  -- |A seek request within an @Iteratee@.-data SeekException = SeekException FileOffset+newtype SeekException = SeekException FileOffset   deriving (Typeable, Show)  instance Exception SeekException where@@ -197,7 +197,7 @@ instance IException EofException where  -- |An @Iteratee exception@ specified by a @String@.-data IterStringException = IterStringException String deriving (Typeable, Show)+newtype IterStringException = IterStringException String deriving (Typeable, Show)  instance Exception IterStringException where   toException   = iterExceptionToException

@@ -56,10 +56,10 @@ -- |The enumerator of a POSIX File Descriptor: a variation of @enumFd@ that -- supports RandomIO (seek requests). enumFdRandom :: MonadIO m => Int -> Fd -> Enumerator Bytes m a-enumFdRandom bs fd iter = enumFdCatch bs fd handler iter+enumFdRandom bs fd = enumFdCatch bs fd handler   where     handler (SeekException off) =-        Nothing <$ (liftIO . fdSeek fd AbsoluteSeek $ fromIntegral off)+        liftIO $ fdSeek fd AbsoluteSeek (fromIntegral off) >> return Nothing  enumFile' :: (MonadIO m, MonadMask m) =>   (Int -> Fd -> Enumerator s m a)

@@ -2,7 +2,6 @@             ,RankNTypes             ,FlexibleContexts             ,ScopedTypeVariables-            ,BangPatterns             ,DeriveDataTypeable #-}  -- |Monadic and General Iteratees:@@ -516,7 +515,7 @@ -- | Enumerate chunks from a list -- enumList :: (Monad m) => [s] -> Enumerator s m a-enumList chunks = go chunks+enumList = go  where   go [] i = return i   go xs' i = runIter i idoneM (onCont xs')@@ -544,8 +543,8 @@   (st -> m (Either SomeException ((Bool, st), s)))   -> st   -> Enumerator s m a-enumFromCallback c st =-  enumFromCallbackCatch c (\NotAnException -> return Nothing) st+enumFromCallback c =+  enumFromCallbackCatch c (\NotAnException -> return Nothing)  -- Dummy exception to catch in enumFromCallback -- This never gets thrown, but it lets us

@@ -770,7 +770,7 @@             onDone a str'@(Chunk c') =                 od (a, newLen - fromIntegral (length c')) str'             onDone a str'@(EOF _) = od (a, n) str'-            onCont f' mExc = oc (go newLen f') mExc+            onCont f' = oc (go newLen f') {-# INLINE countConsumed #-}  -- ------------------------------------------------------------------------@@ -808,8 +808,5 @@  -- | Folds a monadic function over an 'Iteratee'. foldStreamM :: Monad m => (b -> a -> m b) -> b -> Iteratee [a] m b-foldStreamM k = foldChunksM go-  where-    go b [   ] = return b-    go b (h:t) = k b h >>= \b' -> go b' t+foldStreamM = foldChunksM . foldM {-# INLINE foldStreamM #-}

@@ -6,12 +6,6 @@     module System.Posix.Files,     module System.Posix.IO,     module System.Posix.Types,-#if !MIN_VERSION_base_prelude(1,1,0)-    module Foreign.Storable,-    module Foreign.Ptr,-    module Foreign.ForeignPtr,-    module Foreign.StablePtr,-#endif      Bytes, LazyBytes,     HashMap,@@ -26,11 +20,6 @@ #endif #endif -#if !MIN_VERSION_base(4,7,0)-    isLeft,-    isRight,-#endif- #if !MIN_VERSION_base(4,8,0)     first,     second,@@ -49,10 +38,8 @@ import BasePrelude #if MIN_VERSION_base(4,9,0)                     hiding ( EOF, log1p, log1pexp, log1mexp, expm1 )-#elif MIN_VERSION_base(4,7,0)-                    hiding ( EOF ) #else-                    hiding ( EOF, block )+                    hiding ( EOF ) #endif  #if !MIN_VERSION_base(4,8,0)@@ -75,13 +62,6 @@ import System.Posix.IO import System.Posix.Types -#if !MIN_VERSION_base_prelude(1,1,0)-import Foreign.Storable-import Foreign.Ptr-import Foreign.ForeignPtr-import Foreign.StablePtr-#endif- import qualified Data.ByteString.Unsafe as B import qualified Data.ByteString.Lazy   as BL import qualified Data.ByteString.Char8  as S@@ -112,12 +92,6 @@ instance Unpack Text       where unpack = T.unpack instance Unpack String     where unpack = id -#if !MIN_VERSION_base(4,7,0)-isLeft, isRight :: Either a b -> Bool-isLeft = either (const False) (const True)-isRight = either (const True) (const False)-#endif- fdPut :: Fd -> Bytes -> IO () fdPut fd s = B.unsafeUseAsCStringLen s $ \(p,l) ->              throwErrnoIf_ (/= fromIntegral l) "fdPut" $@@ -128,7 +102,7 @@  withFd :: FilePath -> OpenMode -> Maybe FileMode -> OpenFileFlags        -> (Fd -> IO a) -> IO a-withFd fp om fm ff k = bracket (openFd fp om fm ff) closeFd k+withFd fp om fm ff = bracket (openFd fp om fm ff) closeFd  -- | Converts 'Bytes' into 'Text'.  This uses UTF8, but if there is an -- error, it pretends it was Latin1.  Evil as this is, it tends to Just

@@ -1,3 +1,14 @@+-- | Would you believe it?  The 2bit format stores blocks of Ns in a table at+-- the beginning of a sequence, then packs four bases into a byte.  So it+-- is neither possible nor necessary to store Ns in the main sequence, and+-- you would think they aren't stored there, right?  And they aren't.+-- Instead Ts are stored which the reader has to replace with Ns.+--+-- The sensible way to treat these is probably to just say there are two+-- kinds of implied annotation (repeats and large gaps for a typical+-- genome), which can be interpreted in whatever way fits.  And that's why+-- we have 'Mask' and 'getSubseqWith'.+ module Bio.TwoBit (         TwoBitFile(..),         TwoBitSequence(..),@@ -27,24 +38,11 @@ import           Data.Binary.Get import qualified Data.ByteString                as B import qualified Data.ByteString.Lazy           as L-import qualified Data.IntMap                    as I+import qualified Data.IntMap.Strict             as I import qualified Data.HashMap.Lazy              as M import qualified Data.Vector.Unboxed            as U import           System.Random --- ^ Would you believe it?  The 2bit format stores blocks of Ns in a table at--- the beginning of a sequence, then packs four bases into a byte.  So it--- is neither possible nor necessary to store Ns in the main sequence, and--- you would think they aren't stored there, right?  And they aren't.--- Instead Ts are stored which the reader has to replace with Ns.------ The sensible way to treat these is probably to just say there are two--- kinds of implied annotation (repeats and large gaps for a typical--- genome), which can be interpreted in whatever way fits.  And that's why--- we have 'Mask' and 'getSubseqWith'.------ TODO:  use binary search for the Int->Int mappings on the raw data?- data TwoBitFile = TBF {     tbf_raw :: B.ByteString,     -- This map is intentionally lazy.  May or may not be important.@@ -91,7 +89,7 @@      readBlockList = getWord32 >>= \n -> liftM2 zip (repM n getWord32) (repM n getWord32) --- | Repeat monadic action 'n' times.  Returns result in reverse(!)+-- | Repeat monadic action @n@ times.  Returns result in reverse(!) -- order, but doesn't build a huge list of thunks in memory. repM :: Monad m => Int -> m a -> m [a] repM n0 m = go [] n0@@ -150,8 +148,8 @@ -- lowercasing.  Here, we take a user supplied function to apply -- masking. getSubseqWith :: (Nucleotide -> Mask -> a) -> TwoBitFile -> Range -> [a]-getSubseqWith maskf tbf (Range { r_pos = Pos { p_seq = chr, p_start = start }, r_length = len }) = do-    let sq1 = maybe (error $ unpack chr ++ " doesn't exist") id $ M.lookup chr (tbf_seqs tbf)+getSubseqWith maskf tbf Range{ r_pos = Pos { p_seq = chr, p_start = start }, r_length = len } = do+    let sq1 = fromMaybe (error $ unpack chr ++ " doesn't exist") $ M.lookup chr (tbf_seqs tbf)     let go = getFwdSubseqWith tbf sq1     if start < 0         then reverse $ take len $ go (maskf . cmp_nt) (-start-len)@@ -162,8 +160,8 @@  -- | Works only in forward direction. getLazySubseq :: TwoBitFile -> Position -> [Nucleotide]-getLazySubseq tbf (Pos { p_seq = chr, p_start = start }) = do-    let sq1 = maybe (error $ unpack chr ++ " doesn't exist") id $ M.lookup chr (tbf_seqs tbf)+getLazySubseq tbf Pos{ p_seq = chr, p_start = start } = do+    let sq1 = fromMaybe (error $ unpack chr ++ " doesn't exist") $ M.lookup chr (tbf_seqs tbf)     let go  = getFwdSubseqWith tbf sq1     if start < 0         then error "sorry, can't go backwards"

@@ -1,7 +1,7 @@ {-# LANGUAGE ForeignFunctionInterface #-} module Bio.Util.MMap where -import Bio.Prelude hiding ( left, right, chr )+import Bio.Prelude import Data.ByteString.Internal ( fromForeignPtr ) import Foreign.C.Types 

@@ -1,19 +1,18 @@+-- | Random useful stuff I didn't know where to put.+ module Bio.Util.Numeric (     wilson, invnormcdf, choose,     estimateComplexity, showNum, showOOM,     log1p, expm1, (<#>),     log1mexp, log1pexp,-    lsum, llerp,-    sigmoid2, isigmoid2+    lsum, llerp                 ) where  import Data.List ( foldl1' ) import Data.Char ( intToDigit ) import Prelude --- ^ Random useful stuff I didn't know where to put.---- | calculates the Wilson Score interval.+-- | Calculates the Wilson Score interval. -- If @(l,m,h) = wilson c x n@, then @m@ is the binary proportion and -- @(l,h)@ it's @c@-confidence interval for @x@ positive examples out of -- @n@ observations.  @c@ is typically something like 0.05.@@ -33,9 +32,9 @@ showNum = triplets [] . reverse . show   where     triplets acc [] = acc-    triplets acc (a:[]) = a:acc-    triplets acc (a:b:[]) = b:a:acc-    triplets acc (a:b:c:[]) = c:b:a:acc+    triplets acc [a] = a:acc+    triplets acc [a,b] = b:a:acc+    triplets acc [a,b,c] = c:b:a:acc     triplets acc (a:b:c:s) = triplets (',':c:b:a:acc) s  showOOM :: Double -> String@@ -51,8 +50,8 @@                                             '0' : if s == '.' then [] else [s]                         | otherwise = findSuffix (y*0.001) ss --- Stolen from Lennart Augustsson's erf package, who in turn took it rom--- http://home.online.no/~pjacklam/notes/invnorm/ Accurate to about 1e-9.+-- Stolen from Lennart Augustsson's erf package, who in turn took it from+-- <http://home.online.no/~pjacklam/notes/invnorm/> Accurate to about 1e-9. invnormcdf :: (Ord a, Floating a) => a -> a invnormcdf p =     let a1 = -3.969683028665376e+01@@ -108,25 +107,31 @@ -- How many different things are there? -- -- Let the total number be @m@.  The copy number follows a Poisson--- distribution with paramter @\lambda@.  Let @z := e^{\lambda}@, then+-- distribution with paramter @\lambda@.  Let \( z := e^{\lambda} \), then -- we have: -----   P( 0 ) = e^{-\lambda} = 1/z---   P( 1 ) = \lambda e^{-\lambda} = ln z / z---   P(>=1) = 1 - e^{-\lambda} = 1 - 1/z------   singles = m ln z / z---   total   = m (1 - 1/z)------   D := total/singles = (1 - 1/z) * z / ln z---   f := z - 1 - D ln z = 0+-- \[+--   P( 0 ) = e^{-\lambda} = \frac{1}{z}                    \\+--   P( 1 ) = \lambda e^{-\lambda} = \frac{\ln z}{z}                   \\+--   P(\ge 1) = 1 - e^{-\lambda} = 1 - \frac{1}{z}                    \\+-- \]+-- \[+--   \mbox{singles} = m \frac{\ln z}{z}                   \\+--   \mbox{total}   = m \left( 1 - \frac{1}{z} \right)                  \\+-- \]+-- \[+--   D := \frac{\mbox{total}}{\mbox{singles}} = (1 - \frac{1}{z}) * \frac{z}{\ln z}                  \\+--   f := z - 1 - D \ln z = 0+-- \] -- -- To get @z@, we solve using Newton iteration and then substitute to -- get @m@: -----   df/dz = 1 - D/z---   z' := z - z (z - 1 - D ln z) / (z - D)---   m = singles * z /log z+-- \[+--   df/dz = 1 - D/z                                    \\+--   z' = z - \frac{ z (z - 1 - D \ln z) }{ z - D }     \\+--   m = \mbox{singles} * \frac{z}{\ln z}+-- \] -- -- It converges as long as the initial @z@ is large enough, and @10D@ -- (in the line for @zz@ below) appears to work well.@@ -144,7 +149,7 @@     m = fromIntegral singles * zz / log zz  --- | Computes @log (exp x + exp y)@ without leaving the log domain and+-- | Computes \( \ln \left( e^x + e^y \right) \) without leaving the log domain and -- hence without losing precision. infixl 5 <#> {-# INLINE (<#>) #-}@@ -152,7 +157,7 @@ x <#> y = if x >= y then x + log1pexp (y-x) else y + log1pexp (x-y)  -- | Computes @log (1+x)@ to a relative precision of @10^-8@ even for--- very small @x@.  Stolen from http://www.johndcook.com/cpp_log_one_plus_x.html+-- very small @x@.  Stolen from <http://www.johndcook.com/cpp_log_one_plus_x.html> {-# INLINE log1p #-} log1p :: (Floating a, Ord a) => a -> a log1p x | x < -1 = error "log1p: argument must be greater than -1"@@ -163,20 +168,20 @@         | otherwise = (1 - 0.5*x) * x  --- | Computes @exp x - 1@ to a relative precision of @10^-10@ even for--- very small @x@.  Stolen from http://www.johndcook.com/cpp_expm1.html+-- | Computes \( e^x - 1 \) to a relative precision of @10^-10@ even for+-- very small @x@.  Stolen from <http://www.johndcook.com/cpp_expm1.html> {-# INLINE expm1 #-} expm1 :: (Floating a, Ord a) => a -> a expm1 x | x > -0.00001 && x < 0.00001 = (1 + 0.5 * x) * x       -- Taylor approx         | otherwise                   = exp x - 1               -- direct eval --- | Computes @log (1 - exp x)@, following Martin Mächler.+-- | Computes \( \ln (1 - e^x) \), following Martin Mächler. {-# INLINE log1mexp #-} log1mexp :: (Floating a, Ord a) => a -> a log1mexp x | x > - log 2 = log (- expm1 x)            | otherwise   = log1p (- exp x) --- | Computes @log (1 + exp x)@, following Martin Mächler.+-- | Computes \( \ln (1 + e^x) \), following Martin Mächler. {-# INLINE log1pexp #-} log1pexp :: (Floating a, Ord a) => a -> a log1pexp x | x <=  -37 = exp x@@ -185,14 +190,14 @@            | otherwise = x  --- | Computes \( \log ( \sum_i e^{x_i} ) \) sensibly.  The list must be+-- | Computes \( \ln ( \sum_i e^{x_i} ) \) sensibly.  The list must be -- sorted in descending(!) order. {-# INLINE lsum #-} lsum :: (Floating a, Ord a) => [a] -> a-lsum xs = foldl1' (\x y -> if x >= y then x + log1pexp (y-x) else err) xs-    where err = error $ "lsum: argument list must be in descending order"+lsum = foldl1' (\x y -> if x >= y then x + log1pexp (y-x) else err)+    where err = error "lsum: argument list must be in descending order" --- | Computes \( \log \left( c e^x + (1-c) e^y \right) \).+-- | Computes \( \ln \left( c e^x + (1-c) e^y \right) \). {-# INLINE llerp #-} llerp :: (Floating a, Ord a) => a -> a -> a -> a llerp c x y | c <= 0.0  = y@@ -200,20 +205,9 @@             | x >= y    = log     c  + x + log1p ( (1-c)/c * exp (y-x) )        -- Hmm.             | otherwise = log1p (-c) + y + log1p ( c/(1-c) * exp (x-y) )        -- Hmm. --- | Binomial coefficient:  @n `choose` k == n! / ((n-k)! k!)@+-- | Binomial coefficient: \( \mbox{choose n k} = \frac{n!}{(n-k)! k!} \) {-# INLINE choose #-} choose :: Integral a => a -> a -> a-n `choose` k = product [n-k+1 .. n] `div` product [2..k]----- | Kind-of sigmoid function that maps the reals to the interval--- @[0,1)@.  Good to compute a probability without introducing boundary--- conditions.-sigmoid2 :: (Fractional a, Floating a) => a -> a-sigmoid2 x = y*y where y = (exp x - 1) / (exp x + 1)---- | Inverse of 'sigmoid2'.-isigmoid2 :: (Fractional a, Floating a) => a -> a-isigmoid2 y = log $ (1 + sqrt y) / (1 - sqrt y)+choose n k = product [n-k+1 .. n] `div` product [2..k]  

@@ -10,10 +10,9 @@  -- | Decompresses Gzip or Bgzf and passes everything else on.  In -- reality, it simply decompresses Gzip, and when done, looks for--- another Gzip stream.  Trailing garbage is returned as is, therefore,--- uncompressed files are passed through.  Since there is a small chance--- to attempt compression of an uncompressed stream, the original data--- is returned in case of an error.+-- another Gzip stream.  Since there is a small chance to attempt+-- decompression of an uncompressed stream, the original data is+-- returned in case of an error. decompressGzip :: L.ByteString -> L.ByteString decompressGzip s = case L.uncons s of     Just (31, s') -> case L.uncons s' of

@@ -1,10 +1,10 @@ #ifndef INCLUDED_MYERS_ALIGN #define INCLUDED_MYERS_ALIGN -enum myers_align_mode { -	myers_align_globally = 0,-	myers_align_is_prefix = 1,-	myers_align_has_prefix = 2 } ;+enum myers_align_mode {+    myers_align_globally = 0,+    myers_align_is_prefix = 1,+    myers_align_has_prefix = 2 } ;  //! \brief aligns two sequences in O(nd) time //! This alignment algorithm following Eugene W. Myers: "An O(ND)