biohazard 1.0.0 → 1.0.1

raw patch · 17 files changed

+328/−204 lines, 17 filesdep +semigroupsdep −binarydep −directorydep −filepathdep ~asyncdep ~basePVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependencies added: semigroups

Dependencies removed: binary, directory, filepath, random, scientific

Dependency ranges changed: async, base

API changes (from Hackage documentation)

Files

• biohazard.cabal +9/−16
• src/Bio/Adna.hs +19/−13
• src/Bio/Bam/Header.hs +13/−7
• src/Bio/Bam/Index.hs +13/−9
• src/Bio/Bam/Pileup.hs +7/−4
• src/Bio/Bam/Rec.hs +33/−45
• src/Bio/Bam/Trim.hs +5/−5
• src/Bio/Base.hs +7/−7
• src/Bio/Iteratee/Base.hs +18/−14
• src/Bio/Iteratee/Bgzf.hsc +6/−1
• src/Bio/Iteratee/Exception.hs +1/−1
• src/Bio/Iteratee/IO.hs +27/−24
• src/Bio/Iteratee/Iteratee.hs +10/−17
• src/Bio/Iteratee/List.hs +6/−8
• src/Bio/Prelude.hs +16/−15
• src/Bio/TwoBit.hs +37/−18
• src/Bio/Util/Storable.hs +101/−0

@@ -1,5 +1,5 @@ Name:                biohazard-Version:             1.0.0+Version:             1.0.1 Synopsis:            bioinformatics support library Description:         This is a collection of modules I separated from                      various bioinformatics tools.@@ -15,7 +15,7 @@  Cabal-version:       >= 1.10 Build-type:          Simple-Tested-With:         GHC == 7.8.4, GHC == 7.10.3, GHC == 8.0.2, GHC == 8.2.1+Tested-with:         GHC == 7.8.4, GHC == 7.10.3, GHC == 8.0.2, GHC == 8.2.1, GHC == 8.4.1  source-repository head   type:     git@@ -52,22 +52,18 @@                        Bio.TwoBit,                        Bio.Util.MMap,                        Bio.Util.Numeric,+                       Bio.Util.Storable,                        Bio.Util.Zlib -  Build-depends:       async                    >= 2.0 && < 2.2,+  Build-depends:       async                    >= 2.0 && < 2.3,                        attoparsec               >= 0.10 && < 0.14,-                       base                     >= 4.7 && < 4.11,+                       base                     >= 4.7 && < 4.12,                        base-prelude             == 1.2.*,-                       binary                   >= 0.7 && < 0.9,                        bytestring               >= 0.10.2 && < 0.11,                        containers               == 0.5.*,-                       directory                >= 1.2 && < 1.4,                        exceptions               >= 0.6 && < 0.9,-                       filepath                 >= 1.3 && < 1.5,                        hashable                 >= 1.0 && < 1.3,                        primitive                >= 0.5 && < 0.7,-                       random                   >= 1.0 && < 1.2,-                       scientific               == 0.3.*,                        stm                      == 2.4.*,                        text                     >= 1.0 && < 1.3,                        transformers             >= 0.4.1 && < 0.6,@@ -78,6 +74,9 @@                        vector-th-unbox          == 0.2.*,                        zlib                     == 0.6.* +  if !impl(ghc >= 8.0)+    build-depends: semigroups == 0.18.*+   Ghc-options:         -Wall    Default-Language:    Haskell2010@@ -97,16 +96,10 @@                        DeriveGeneric,                        ExistentialQuantification,                        ForeignFunctionInterface,-                       FunctionalDependencies,                        GeneralizedNewtypeDeriving,-                       PatternGuards,                        Rank2Types,-                       ScopedTypeVariables,                        TemplateHaskell,-                       TupleSections, -                       TypeFamilies,-                       TypeOperators,-                       UndecidableInstances+                       TypeFamilies    Hs-source-dirs:      src   Include-dirs:        src/cbits

@@ -506,6 +506,7 @@   instance Monoid a => Monoid (DmgStats a) where+    mappend = merge_dmg_stats mappend     mempty = DmgStats { basecompo5 = empty_compo                       , basecompo3 = empty_compo                       , substs5    = empty_subst@@ -523,19 +524,24 @@         empty_compo = [ (nuc, U.empty) | nuc <- [Just nucA, Just nucC, Just nucG, Just nucT, Nothing] ]         empty_subst = [ (n1 :-> n2, U.empty) | n1 <- [nucA..nucT], n2 <- [nucA..nucT] ] -    a `mappend` b = DmgStats { basecompo5 = zipWith s1 (basecompo5 a) (basecompo5 b)-                             , basecompo3 = zipWith s1 (basecompo3 a) (basecompo3 b)-                             , substs5    = zipWith s2 (substs5    a) (substs5    b)-                             , substs3    = zipWith s2 (substs3    a) (substs3    b)-                             , substs5d5  = zipWith s2 (substs5d5  a) (substs5d5  b)-                             , substs3d5  = zipWith s2 (substs3d5  a) (substs3d5  b)-                             , substs5d3  = zipWith s2 (substs5d3  a) (substs5d3  b)-                             , substs3d3  = zipWith s2 (substs3d3  a) (substs3d3  b)-                             , substs5dd  = zipWith s2 (substs5dd  a) (substs5dd  b)-                             , substs3dd  = zipWith s2 (substs3dd  a) (substs3dd  b)-                             , substs5cpg = zipWith s2 (substs5cpg a) (substs5cpg b)-                             , substs3cpg = zipWith s2 (substs3cpg a) (substs3cpg b)-                             , stats_more = mappend    (stats_more a) (stats_more b) }+instance Semigroup a => Semigroup (DmgStats a) where+    (<>) = merge_dmg_stats (<>)++merge_dmg_stats :: (a -> a -> a) -> DmgStats a -> DmgStats a -> DmgStats a+merge_dmg_stats plus a b =+    DmgStats { basecompo5 = zipWith s1 (basecompo5 a) (basecompo5 b)+             , basecompo3 = zipWith s1 (basecompo3 a) (basecompo3 b)+             , substs5    = zipWith s2 (substs5    a) (substs5    b)+             , substs3    = zipWith s2 (substs3    a) (substs3    b)+             , substs5d5  = zipWith s2 (substs5d5  a) (substs5d5  b)+             , substs3d5  = zipWith s2 (substs3d5  a) (substs3d5  b)+             , substs5d3  = zipWith s2 (substs5d3  a) (substs5d3  b)+             , substs3d3  = zipWith s2 (substs3d3  a) (substs3d3  b)+             , substs5dd  = zipWith s2 (substs5dd  a) (substs5dd  b)+             , substs3dd  = zipWith s2 (substs3dd  a) (substs3dd  b)+             , substs5cpg = zipWith s2 (substs5cpg a) (substs5cpg b)+             , substs3cpg = zipWith s2 (substs3cpg a) (substs3cpg b)+             , stats_more = plus       (stats_more a) (stats_more b) }       where         s1 (x, u) (z, v) | x /= z    = error "Mismatch in zip.  This is a bug."                          | U.null u  = (x, v)

@@ -106,11 +106,14 @@  instance Monoid BamMeta where     mempty = BamMeta mempty noRefs [] []-    a `mappend` b = BamMeta { meta_hdr = meta_hdr a `mappend` meta_hdr b-                            , meta_refs = meta_refs a >< meta_refs b-                            , meta_other_shit = nub $ meta_other_shit a ++ meta_other_shit b-                            , meta_comment = nub $ meta_comment a ++ meta_comment b }+    mappend = (<>) +instance Semigroup BamMeta where+    a <> b = BamMeta { meta_hdr = meta_hdr a `mappend` meta_hdr b+                     , meta_refs = meta_refs a >< meta_refs b+                     , meta_other_shit = nub $ meta_other_shit a ++ meta_other_shit b+                     , meta_comment = nub $ meta_comment a ++ meta_comment b }+ data BamHeader = BamHeader {         hdr_version :: (Int, Int),         hdr_sorting :: !BamSorting,@@ -119,9 +122,12 @@  instance Monoid BamHeader where     mempty = BamHeader (1,0) Unknown []-    a `mappend` b = BamHeader { hdr_version = hdr_version a `min` hdr_version b-                              , hdr_sorting = let u = hdr_sorting a ; v = hdr_sorting b in if u == v then u else Unknown-                              , hdr_other_shit = hdr_other_shit a ++ hdr_other_shit b }+    mappend = (<>)++instance Semigroup BamHeader where+    a <> b = BamHeader { hdr_version = hdr_version a `min` hdr_version b+                       , hdr_sorting = let u = hdr_sorting a ; v = hdr_sorting b in if u == v then u else Unknown+                       , hdr_other_shit = hdr_other_shit a ++ hdr_other_shit b }  data BamSQ = BamSQ {         sq_name :: Seqid,

@@ -18,8 +18,7 @@ import Bio.Bam.Regions              ( Region(..), Subsequence(..) ) import Bio.Iteratee import Bio.Prelude-import System.Directory             ( doesFileExist )-import System.FilePath              ( dropExtension, takeExtension, (<.>) )+import System.Posix.Files           ( fileExist )  import qualified Bio.Bam.Regions                as R import qualified Data.IntMap.Strict             as M@@ -145,15 +144,20 @@ -- two, and finally in the file itself.  The first file that exists and -- can actually be parsed, is used. readBamIndex :: FilePath -> IO (BamIndex ())-readBamIndex fp | takeExtension fp == ".bai" = enumFile defaultBufSize fp readBaiIndex >>= run-                | takeExtension fp == ".csi" = enumFile defaultBufSize fp readBaiIndex >>= run-                | otherwise = tryIx               (fp <.> "bai") $-                              tryIx (dropExtension fp <.> "bai") $-                              tryIx               (fp <.> "csi") $-                              tryIx (dropExtension fp <.> "csi") $+readBamIndex fp | ".bai" `isSuffixOf` fp = enumFile defaultBufSize fp readBaiIndex >>= run+                | ".csi" `isSuffixOf` fp = enumFile defaultBufSize fp readBaiIndex >>= run+                | otherwise = tryIx               (fp ++ ".bai") $+                              tryIx (dropExtension fp ++  "bai") $+                              tryIx               (fp ++ ".csi") $+                              tryIx (dropExtension fp ++  "csi") $                               enumFile defaultBufSize fp readBaiIndex >>= run   where-    tryIx f k = do e <- doesFileExist f+    dropExtension p = reverse $ (if null b then "." ++ f else b) ++ d+      where+        (f,d) = break (=='/') $ reverse p+        b     = dropWhile (/='.') f++    tryIx f k = do e <- fileExist f                    if e then do r <- enumFile defaultBufSize f readBaiIndex >>= tryRun                                 case r of Right                     ix -> return ix                                           Left (IterStringException _) -> k

@@ -200,10 +200,13 @@                             , reads_mapq0      = 0                             , sum_mapq         = 0                             , sum_mapq_squared = 0 }-    mappend x y = CallStats { read_depth       = read_depth x + read_depth y-                            , reads_mapq0      = reads_mapq0 x + reads_mapq0 y-                            , sum_mapq         = sum_mapq x + sum_mapq y-                            , sum_mapq_squared = sum_mapq_squared x + sum_mapq_squared y }+    mappend     = (<>)++instance Semigroup CallStats where+    x <> y = CallStats { read_depth       = read_depth x + read_depth y+                       , reads_mapq0      = reads_mapq0 x + reads_mapq0 y+                       , sum_mapq         = sum_mapq x + sum_mapq y+                       , sum_mapq_squared = sum_mapq_squared x + sum_mapq_squared y }  newtype V_Nuc  = V_Nuc  (U.Vector Nucleotide)  deriving (Eq, Ord, Show) newtype V_Nucs = V_Nucs (U.Vector Nucleotides) deriving (Eq, Ord, Show)

@@ -50,6 +50,7 @@ import Bio.Bam.Header import Bio.Iteratee import Bio.Prelude+import Bio.Util.Storable  import Control.Monad.Primitive      ( unsafePrimToPrim, unsafeInlineIO ) import Foreign.C.Types              ( CInt(..), CSize(..) )@@ -81,21 +82,12 @@     sizeOf    _ = 4     alignment _ = 1 -    peek p = do w0 <- peekByteOff p 0 :: IO Word8-                w1 <- peekByteOff p 1 :: IO Word8-                w2 <- peekByteOff p 2 :: IO Word8-                w3 <- peekByteOff p 3 :: IO Word8-                let w = fromIntegral w0 `shiftL`  0 .|.  fromIntegral w1 `shiftL`  8 .|.-                        fromIntegral w2 `shiftL` 16 .|.  fromIntegral w3 `shiftL` 24+    peek p = do w <- fromIntegral <$> peekUnalnWord32LE p                 return $ toEnum (w .&. 0xf) :* shiftR w 4 -    poke p (op :* num) = do pokeByteOff p 0 (fromIntegral $ shiftR w  0 :: Word8)-                            pokeByteOff p 1 (fromIntegral $ shiftR w  8 :: Word8)-                            pokeByteOff p 2 (fromIntegral $ shiftR w 16 :: Word8)-                            pokeByteOff p 3 (fromIntegral $ shiftR w 24 :: Word8)-        where-            w = fromEnum op .|. shiftL num 4+    poke p (op :* num) = pokeUnalnWord32LE p . fromIntegral $ fromEnum op .|. shiftL num 4 + -- | 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@@ -235,13 +227,13 @@ {-# INLINE[1] unpackBam #-} unpackBam :: BamRaw -> BamRec unpackBam br = BamRec {-        b_rname =      Refseq $ getInt32  0,-        b_pos   =               getInt32  4,-        b_mapq  =           Q $ getInt8   9,-        b_flag  =               getInt16 14,-        b_mrnm  =      Refseq $ getInt32 20,-        b_mpos  =               getInt32 24,-        b_isize = fromIntegral (getInt32 28 :: Int32),+        b_rname =      Refseq $ getWord32  0,+        b_pos   =               getInt32   4,+        b_mapq  =           Q $ getInt8    9,+        b_flag  =               getInt16  14,+        b_mrnm  =      Refseq $ getWord32 20,+        b_mpos  =               getInt32  24,+        b_isize =               getInt32  28,          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),@@ -257,29 +249,27 @@         off_q = off_s + (l_seq + 1) `div` 2         off_e = off_q +  l_seq -        l_read_name = getInt8   8 - 1-        l_seq       = getInt32 16-        l_cigar     = getInt16 12+        l_read_name = getInt8    8 - 1+        l_seq       = getWord32 16+        l_cigar     = getInt16  12          getInt8 :: Num a => Int -> a         getInt8  o = fromIntegral (B.unsafeIndex (raw_data br) o) -        getInt16 :: (Num a, Bits a) => Int -> a-        getInt16 o = fromIntegral (B.unsafeIndex (raw_data br) o) .|.-                     fromIntegral (B.unsafeIndex (raw_data br) $ o+1) `shiftL`  8+        getInt16 :: Num a => Int -> a+        getInt16 o = unsafeDupablePerformIO $ B.unsafeUseAsCString (raw_data br) $+                     fmap fromIntegral . peekUnalnWord16LE . flip plusPtr o -        getInt32 :: (Num a, Bits a) => Int -> a-        getInt32 o = fromIntegral (        B.unsafeIndex (raw_data br) $ o+0 )             .|.-                     fromIntegral (        B.unsafeIndex (raw_data br) $ o+1 ) `shiftL`  8 .|.-                     fromIntegral (        B.unsafeIndex (raw_data br) $ o+2 ) `shiftL` 16 .|.-                     fromIntegral (signed (B.unsafeIndex (raw_data br) $ o+3)) `shiftL` 24+        getWord32 :: Num a => Int -> a+        getWord32 o = unsafeDupablePerformIO $ B.unsafeUseAsCString (raw_data br) $+                      fmap fromIntegral . peekUnalnWord32LE . flip plusPtr o -        signed :: Word8 -> Int8-        signed = fromIntegral+        -- ensures proper sign extension+        getInt32 :: Num a => Int -> a+        getInt32 o = fromIntegral (getWord32 o :: Int32) --- | A collection of extension fields.  The key is actually only two @Char@s, but that proved impractical.--- (Hmm... we could introduce a Key type that is a 16 bit int, then give--- it an @instance IsString@... practical?)+-- | A collection of extension fields.  A 'BamKey' is actually two ASCII+-- characters. type Extensions = [( BamKey, Ext )]  -- | Deletes all occurences of some extension field.@@ -337,16 +327,14 @@      getInt 'C' s | S.length s >= 1 = fromIntegral (fromIntegral (B.index s 0) :: Word8)     getInt 'c' s | S.length s >= 1 = fromIntegral (fromIntegral (B.index s 0) ::  Int8)-    getInt 'S' s | S.length s >= 2 = fromIntegral                               (i :: Word16)-        where i = fromIntegral (B.index s 0) .|. fromIntegral (B.index s 1) `shiftL` 8-    getInt 's' s | S.length s >= 2 = fromIntegral                               (i ::  Int16)-        where i = fromIntegral (B.index s 0) .|. fromIntegral (B.index s 1) `shiftL` 8-    getInt 'I' s | S.length s >= 4 = fromIntegral                               (i :: Word32)-        where i = fromIntegral (B.index s 0)             .|. fromIntegral (B.index s 1) `shiftL`  8 .|.-                  fromIntegral (B.index s 2) `shiftL` 16 .|. fromIntegral (B.index s 3) `shiftL` 24-    getInt 'i' s | S.length s >= 4 = fromIntegral                               (i ::  Int32)-        where i = fromIntegral (B.index s 0)             .|. fromIntegral (B.index s 1) `shiftL`  8 .|.-                  fromIntegral (B.index s 2) `shiftL` 16 .|. fromIntegral (B.index s 3) `shiftL` 24+    getInt 'S' s | S.length s >= 2 = fromIntegral                         (i :: Word16)+        where i = unsafeDupablePerformIO $ B.unsafeUseAsCString s $ peekUnalnWord16LE+    getInt 's' s | S.length s >= 2 = fromIntegral            (fromIntegral i ::  Int16)+        where i = unsafeDupablePerformIO $ B.unsafeUseAsCString s $ peekUnalnWord16LE+    getInt 'I' s | S.length s >= 4 = fromIntegral                         (i :: Word32)+        where i = unsafeDupablePerformIO $ B.unsafeUseAsCString s $ peekUnalnWord32LE+    getInt 'i' s | S.length s >= 4 = fromIntegral            (fromIntegral i ::  Int32)+        where i = unsafeDupablePerformIO $ B.unsafeUseAsCString s $ peekUnalnWord32LE     getInt _ _ = 0      getFloat s = unsafeDupablePerformIO $ alloca $ \buf ->

@@ -192,8 +192,8 @@  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) = [ r1', r2'     ]-    | qual1 < lowq                           = [ 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 ]@@ -279,9 +279,9 @@ -- and two qualities). trimBam :: Int -> Int -> [W.Vector Nucleotides] -> BamRec -> [BamRec] trimBam lowq highq ads1 r1-    | V.null (b_seq r1)              = [ r1'      ]+    | V.null (b_seq r1)              = [          ]     | mlen == 0 && qual1 >= highq    = [          ]-    | qual1 < lowq                   = [ r1'      ]+    | qual1 < lowq || mlen < 0       = [ r1'      ]     | qual1 >= highq                 = [      r1t ]     | otherwise = map flag_alternative [ r1', r1t ]   where@@ -417,7 +417,7 @@  {-# INLINE twoMins #-} twoMins :: (Bounded a, Ord a) => a -> Int -> (Int -> a) -> (a,Int,a)-twoMins a0 imax f = go a0 0 maxBound 0 0+twoMins a0 imax f = go a0 (-1) maxBound 0 0   where     go !m1 !i1 !m2 !i2 !i         | i == imax = (m1,i1,m2)

@@ -46,8 +46,7 @@ import Bio.Util.Numeric             ( log1pexp, log1mexp ) import Data.ByteString.Internal     ( c2w, w2c ) import Data.Vector.Unboxed.Deriving ( derivingUnbox )-import System.Directory             ( doesFileExist )-import System.FilePath              ( (</>), isAbsolute, splitSearchPath )+import System.Posix.Files           ( fileExist )  import qualified Data.ByteString.Char8 as S import qualified Data.Vector.Unboxed   as U@@ -360,10 +359,11 @@ -- | Finds a file by searching the environment variable BIOHAZARD like a -- PATH. findAuxFile :: FilePath -> IO FilePath-findAuxFile fn | isAbsolute fn = return fn-               | otherwise = go . maybe ["."] splitSearchPath . lookup "BIOHAZARD" =<< getEnvironment+findAuxFile fn | "/" `isPrefixOf` fn = return fn+               | otherwise = go . fromMaybe "." . lookup "BIOHAZARD" =<< getEnvironment   where-    go [    ] = return fn-    go (p:ps) = doesFileExist (p </> fn) >>=-                bool (return $ p </> fn) (go ps)+    go "" = return fn+    go pp = let (p,ps) = break (==':') pp+            in fileExist (p ++ "/" ++ fn) >>=+               bool (return $ p ++ "/" ++ fn) (go $ drop 1 ps) 

@@ -1,4 +1,4 @@-{-# LANGUAGE TypeFamilies,UndecidableInstances,Rank2Types,ExistentialQuantification #-}+{-# LANGUAGE Rank2Types #-}  -- |Monadic Iteratees: -- incremental input parsers, processors and transformers@@ -76,11 +76,16 @@   (EOF (Just e1)) == (EOF (Just e2)) = typeOf e1 == typeOf e2   _ == _                             = False +instance Semigroup c => Semigroup (Stream c) where+  EOF mErr <>        _ = EOF mErr+  _        <> EOF mErr = EOF mErr+  Chunk s1 <> Chunk s2 = Chunk (s1 <> s2)+ instance Monoid c => Monoid (Stream c) where-  mempty                        = Chunk mempty-  mappend (EOF mErr)         _  = EOF mErr-  mappend         _  (EOF mErr) = EOF mErr-  mappend (Chunk s1) (Chunk s2) = Chunk (s1 `mappend` s2)+  mempty                      = Chunk mempty+  EOF mErr `mappend`        _ = EOF mErr+  _        `mappend` EOF mErr = EOF mErr+  Chunk s1 `mappend` Chunk s2 = Chunk (mappend s1 s2)  -- |Map a function over a stream. instance Functor Stream where@@ -154,15 +159,14 @@     => Iteratee s m a     -> (a -> Iteratee s m b)     -> Iteratee s m b-bindIteratee = self-    where-        self m f = Iteratee $ \onDone onCont ->-             let m_done a (Chunk s)-                   | nullC s     = runIter (f a) onDone onCont-                 m_done a stream = runIter (f a) (const . flip onDone stream) f_cont-                   where f_cont k Nothing = runIter (k stream) onDone onCont-                         f_cont k e       = onCont k e-             in runIter m m_done (onCont . (flip self f .))+bindIteratee m f = Iteratee $ \onDone onCont ->+    let m_done a (Chunk s)+          | nullC s     = runIter (f a) onDone onCont+        m_done a stream = runIter (f a) (const . flip onDone stream) f_cont+          where f_cont k Nothing = runIter (k stream) onDone onCont+                f_cont k e       = onCont k e+    in runIter m m_done (onCont . (flip bindIteratee f .))+  instance NullPoint s => MonadTrans (Iteratee s) where   lift m = Iteratee $ \onDone _ -> m >>= flip onDone (Chunk emptyP)

@@ -34,9 +34,14 @@  instance Monoid Block where     mempty = Block 0 S.empty-    mappend (Block x s) (Block _ t) = Block x (s `S.append` t)+    mappend = (<>)     mconcat [] = Block 0 S.empty     mconcat bs@(Block x _:_) = Block x $ S.concat [s|Block _ s <- bs]++instance Semigroup Block where+    Block x s <> Block _ t = Block x (s `S.append` t)+    sconcat (Block x b :| bs) = Block x . sconcat $ b :| map block_contents bs+  -- | "Decompresses" a plain file.  What's actually happening is that the -- offset in the input stream is tracked and added to the @Bytes@s

@@ -1,4 +1,4 @@-{-# LANGUAGE DeriveDataTypeable, ExistentialQuantification #-}+{-# LANGUAGE ExistentialQuantification #-}  -- |Monadic and General Iteratees: -- Messaging and exception handling.

@@ -10,18 +10,18 @@   ,enumFileRandom   -- * FileDescriptor based enumerators for monadic iteratees   ,enumFd-  ,enumFdCatch   ,enumFdRandom ) where  import Bio.Iteratee.Iteratee-import Bio.Prelude hiding ( bracket )+import Bio.Prelude hiding ( bracket, loop ) import Control.Monad.Catch import Control.Monad.IO.Class-import Data.ByteString.Internal (createAndTrim) import System.IO (SeekMode(..)) +import qualified Data.ByteString as B+ -- | Default buffer size in elements.  This was 1024 in "Data.Iteratee", -- which is obviously too small.  Since we often want to merge many -- files, a read should take more time than a seek.  This sets the@@ -29,37 +29,40 @@ defaultBufSize :: Int defaultBufSize = 2*1024*1024 --- --------------------------------------------------------------------------- Binary Random IO enumerators--makefdCallback :: MonadIO m => Int -> Fd -> st -> m (Either SomeException ((Bool, st), Bytes))-makefdCallback bufsize fd st = do-  s <- liftIO . createAndTrim bufsize $ \p ->-       fromIntegral <$> fdReadBuf fd (castPtr p) (fromIntegral bufsize)-  return $ Right ((True, st), s)- -- |The enumerator of a POSIX File Descriptor.  This version enumerates -- over the entire contents of a file, in order, unless stopped by -- the iteratee.  In particular, seeking is not supported. enumFd :: MonadIO m => Int -> Fd -> Enumerator Bytes m a-enumFd bufsize fd = enumFromCallback (makefdCallback bufsize fd) ()+enumFd bufsize fd = loop+  where+    loop iter = runIter iter idoneM onCont --- |A variant of enumFd that catches exceptions raised by the @Iteratee@.-enumFdCatch- :: (IException e, MonadIO m)-    => Int-    -> Fd-    -> (e -> m (Maybe EnumException))-    -> Enumerator Bytes m a-enumFdCatch bufsize fd handler = enumFromCallbackCatch (makefdCallback bufsize fd) handler ()+    onCont k j@(Just _) = return (icont k j)+    onCont k   Nothing  = do+        s <- liftIO $ fdGet bufsize fd+        if B.null s then return $ liftI k+                    else loop . k $ Chunk s + -- |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 = enumFdCatch bs fd handler+enumFdRandom bs fd = loop   where-    handler (SeekException off) =-        liftIO $ fdSeek fd AbsoluteSeek (fromIntegral off) >> return Nothing+    loop iter = runIter iter idoneM onCont++    onCont k Nothing  = do+        s <- liftIO $ fdGet bs fd+        if B.null s then return $ liftI k+                    else loop . k $ Chunk s++    onCont k j@(Just e) = case fromException e of+      Just (SeekException off) -> do+                   liftIO . void $ fdSeek fd AbsoluteSeek (fromIntegral off)+                   loop $ liftI k+      Nothing -> return (icont k j)++  enumFile' :: (MonadIO m, MonadMask m) =>   (Int -> Fd -> Enumerator s m a)

@@ -1,8 +1,4 @@-{-# LANGUAGE KindSignatures-            ,RankNTypes-            ,FlexibleContexts-            ,ScopedTypeVariables-            ,DeriveDataTypeable #-}+{-# LANGUAGE Rank2Types, FlexibleContexts #-}  -- |Monadic and General Iteratees: -- incremental input parsers, processors and transformers@@ -187,7 +183,7 @@ -- The converters show a different way of composing two iteratees: -- `vertical' rather than `horizontal' -type Enumeratee sFrom sTo (m :: * -> *) a =+type Enumeratee sFrom sTo m a =   Iteratee sTo m a   -> Iteratee sFrom m (Iteratee sTo m a) @@ -543,17 +539,14 @@   (st -> m (Either SomeException ((Bool, st), s)))   -> st   -> Enumerator s m a-enumFromCallback c =-  enumFromCallbackCatch c (\NotAnException -> return Nothing)---- Dummy exception to catch in enumFromCallback--- This never gets thrown, but it lets us--- share plumbing-data NotAnException = NotAnException- deriving (Show, Typeable)--instance Exception NotAnException where-instance IException NotAnException where+enumFromCallback c = loop+  where+    loop st iter = runIter iter idoneM (onCont st)+    check k (True,  st') = loop st' . k . Chunk+    check k (False,_st') = return . k . Chunk+    onCont st k Nothing  = c st >>=+        either (return . k . EOF . Just) (uncurry (check k))+    onCont _st k j = return (icont k j)  -- |Create an enumerator from a callback function with an exception handler. -- The exception handler is called if an iteratee reports an exception.

@@ -1,5 +1,3 @@-{-# LANGUAGE TupleSections, ScopedTypeVariables #-}- -- |Monadic Iteratees: -- incremental input parsers, processors and transformers --@@ -602,8 +600,8 @@     (b', y') <- lift $ runIter y0 od oc     case checkDone a' b' of       Just (Right (a,b,s))  -> idone (a,b) s  -- 's' may be EOF, needs to stay-      Just (Left (Left a))  -> liftM (a,) y'-      Just (Left (Right b)) -> liftM (,b) x'+      Just (Left (Left a))  -> liftM ((,) a) y'+      Just (Left (Right b)) -> liftM (flip (,) b) x'       Nothing               -> liftI (step x' y')   where     step x y (Chunk xs) | nullC xs = liftI (step x y)@@ -612,8 +610,8 @@       (b', y') <- lift $ (\i -> runIter i od oc) =<< enumPure1Chunk xs y       case checkDone a' b' of         Just (Right (a,b,s))  -> idone (a,b) s-        Just (Left (Left a))  -> liftM (a,) y'-        Just (Left (Right b)) -> liftM (,b) x'+        Just (Left (Left a))  -> liftM ((,) a) y'+        Just (Left (Right b)) -> liftM (flip (,) b) x'         Nothing               -> liftI (step x' y')     step x y (EOF err) = joinIM $ case err of       Nothing -> (liftM2.liftM2) (,) (enumEof   x) (enumEof   y)@@ -683,7 +681,7 @@     (a', x') <- lift $ runIter i1 od oc     (_,  y') <- lift $ runIter i2 od oc     case a' of-      Just (a, s) -> flip idone s =<< lift (liftM (a,) $ run i2)+      Just (a, s) -> flip idone s =<< lift (liftM ((,) a) $ run i2)       Nothing     -> go x' y'   where     od a s = return (Just (a, s), idone a s)@@ -761,7 +759,7 @@ countConsumed :: (Monad m, Integral n) => Iteratee [el] m a -> Iteratee [el] m (a, n) countConsumed i = go 0 (const i) (Chunk emptyP)   where-    go !n f str@(EOF _) = (, n) `liftM` f str+    go !n f str@(EOF _) = flip (,) n `liftM` f str     go !n f str@(Chunk c) = Iteratee rI       where         newLen = n + fromIntegral (length c)

@@ -1,7 +1,9 @@-{-# LANGUAGE CPP, TypeOperators #-}+{-# LANGUAGE CPP #-} module Bio.Prelude (     module Bio.Base,     module BasePrelude,+    module Data.List.NonEmpty,+    module Data.Semigroup,     module System.IO,     module System.Posix.Files,     module System.Posix.IO,@@ -14,11 +16,6 @@     IntSet,     Text, LazyText,     Pair(..),-#ifndef __HADDOCK__-#ifdef __GLASGOW_HASKELL__-    (:!:),-#endif-#endif  #if !MIN_VERSION_base(4,8,0)     first,@@ -30,6 +27,7 @@      Hashable(..),     Unpack(..),+    fdGet,     fdPut,     fdPutLazy,     withFd@@ -37,9 +35,9 @@  import BasePrelude #if MIN_VERSION_base(4,9,0)-                    hiding ( EOF, log1p, log1pexp, log1mexp, expm1 )+                    hiding ( (<>), EOF, log1p, log1pexp, log1mexp, expm1 ) #else-                    hiding ( EOF )+                    hiding ( (<>), EOF ) #endif  #if !MIN_VERSION_base(4,8,0)@@ -49,6 +47,9 @@  import Bio.Base import Data.ByteString     ( ByteString )+import Data.ByteString.Internal ( createAndTrim )+import Data.List.NonEmpty  ( NonEmpty(..) )+import Data.Semigroup      ( Semigroup(..) ) import Data.Text           ( Text ) import Data.Hashable       ( Hashable(..) ) import Data.HashMap.Strict ( HashMap )@@ -77,13 +78,6 @@ -- | A strict pair. data Pair a b = !a :!: !b deriving(Eq, Ord, Show, Read, Bounded, Ix) -#ifndef __HADDOCK__-#ifdef __GLASGOW_HASKELL__--- This gives a nicer syntax for the type but only works in GHC for now.-type (:!:) = Pair-#endif-#endif- -- | Class of things that can be unpacked into 'String's.  Kind of the -- opposite of 'IsString'. class Unpack s where unpack :: s -> String@@ -91,6 +85,13 @@ instance Unpack ByteString where unpack = S.unpack instance Unpack Text       where unpack = T.unpack instance Unpack String     where unpack = id++-- | @fdGet bs fd@ reads up to @bs@ 'Bytes' from file descriptor @Fd@.+-- Returns an empty 'Bytes' at end of file.+fdGet :: Int -> Fd -> IO Bytes+fdGet bs fd =+    createAndTrim bs $ \p ->+        fromIntegral <$> fdReadBuf fd (castPtr p) (fromIntegral bs)  fdPut :: Fd -> Bytes -> IO () fdPut fd s = B.unsafeUseAsCStringLen s $ \(p,l) ->

@@ -35,13 +35,14 @@  import           Bio.Prelude hiding ( left, right, chr ) import           Bio.Util.MMap-import           Data.Binary.Get+import           Bio.Util.Storable+import           Control.Monad.Trans.State import qualified Data.ByteString                as B-import qualified Data.ByteString.Lazy           as L+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           System.Random+import           Foreign.C.Types ( CChar )  data TwoBitFile = TBF {     tbf_raw :: B.ByteString,@@ -60,11 +61,13 @@ openTwoBit :: FilePath -> IO TwoBitFile openTwoBit fp = do         raw <- unsafeMMapFile fp-        return $ flip runGet (L.fromChunks [raw]) $ do+        B.unsafeUseAsCString raw $ \praw ->+        -- return $ flip runGet (L.fromChunks [raw]) $ do+            flip evalStateT praw $ do                     sig <- getWord32be-                    getWord32 <- case sig of-                            0x1A412743 -> return $ fromIntegral `fmap` getWord32be-                            0x4327411A -> return $ fromIntegral `fmap` getWord32le+                    getWord32 <- case sig :: Word32 of+                            0x1A412743 -> return getWord32be+                            0x4327411A -> return getWord32le                             _          -> fail $ "invalid .2bit signature " ++ showHex sig []                      version <- getWord32@@ -73,19 +76,34 @@                     nseqs <- getWord32                     _reserved <- getWord32 -                    TBF raw <$> foldM (\ix _ -> do !key <- getWord8 >>= getByteString . fromIntegral+                    TBF raw <$> foldM (\ix _ -> do !key <- getWord8 >>= getByteString                                                    !off <- getWord32                                                    return $! M.insert key (mkBlockIndex raw getWord32 off) ix                                       ) M.empty [1..nseqs] +type Get = StateT (Ptr CChar) IO++getWord8, getWord32be, getWord32le :: Num a => Get a+getWord8    = StateT $ \p -> peekUnalnWord32BE  p >>= \w -> return (fromIntegral w, plusPtr p 1)+getWord32be = StateT $ \p -> peekUnalnWord32BE  p >>= \w -> return (fromIntegral w, plusPtr p 4)+getWord32le = StateT $ \p -> peekUnalnWord32LE  p >>= \w -> return (fromIntegral w, plusPtr p 4)++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 raw getWord32 ofs = runGet getBlock $ L.fromChunks [B.drop ofs raw]+mkBlockIndex raw getWord32 ofs =+    unsafePerformIO $+    B.unsafeUseAsCString raw $ \praw ->+    evalStateT getBlock (plusPtr praw ofs)   where-    getBlock = do ds <- getWord32+    getBlock = do p0 <- get+                  ds <- getWord32                   nb <- readBlockList                   mb <- readBlockList-                  len <- getWord32 >> bytesRead-                  return $! TBS (I.fromList nb) (I.fromList mb) (ofs + fromIntegral len) ds+                  _  <- getWord32+                  p1 <- get+                  return $! TBS (I.fromList nb) (I.fromList mb) (ofs + minusPtr p1 p0) ds      readBlockList = getWord32 >>= \n -> liftM2 zip (repM n getWord32) (repM n getWord32) @@ -222,11 +240,12 @@ -- On a 32bit platform, this will fail for genomes larger than 1G bases. -- However, if you're running this code on a 32bit platform, you have -- bigger problems to worry about.-getRandomSeq :: RandomGen g => TwoBitFile                   -- ^ 2bit file-                            -> Int                          -- ^ desired length-                            -> g                            -- ^ RNG-                            -> ((Range, [Nucleotide]), g)   -- ^ position, sequence, new RNG-getRandomSeq tbf len = draw+getRandomSeq :: TwoBitFile                   -- ^ 2bit file+             -> Int                          -- ^ desired length+             -> (Int -> g -> (Int, g))       -- ^ draw random int below limit+             -> g                            -- ^ RNG+             -> ((Range, [Nucleotide]), g)   -- ^ position, sequence, new RNG+getRandomSeq tbf len rndInt = draw   where     names = getSeqnames tbf     lengths = map (getSeqLength tbf) names@@ -236,7 +255,7 @@     draw g0 | good      = ((r', sq), gn)             | otherwise = draw gn       where-        (p0, gn) = randomR (0, 2*total-1) g0+        (p0, gn) = rndInt (2*total) g0         p = p0 `shiftR` 1         Just ((o,s),_) = I.maxViewWithKey $ fst $ I.split (p+1) frags         r' = (if odd p0 then id else reverseRange) $ clampPosition tbf $ Range (Pos s (p-o)) len

@@ -0,0 +1,101 @@+{-# LANGUAGE CPP #-}+-- | Utilities to read multibyte quantities from arbitrary positions.+module Bio.Util.Storable+    ( peekWord8+    , peekUnalnWord16LE+    , peekUnalnWord16BE+    , peekUnalnWord32LE+    , peekUnalnWord32BE+    , pokeUnalnWord32LE+    ) where++#if __GLASGOW_HASKELL__ >= 710+#define HAVE_BYTESWAP_PRIMOPS+#endif++#if i386_HOST_ARCH || x86_64_HOST_ARCH+#define MEM_UNALIGNED_OPS+#endif++import Bio.Prelude++#if defined(HAVE_BYTESWAP_PRIMOPS)+import GHC.Word ( byteSwap16, byteSwap32 )+#endif++peekWord8 :: Ptr a -> IO Word8+peekWord8 = peek . castPtr++#if defined(MEM_UNALIGNED_OPS) && defined(WORDS_BIGENDIAN) && defined(HAVE_BYTESWAP_PRIMOPS)+peekUnalnWord16LE :: Ptr a -> IO Word16+peekUnalnWord16LE = fmap byteSwap16 . peek . castPtr++peekUnalnWord32LE :: Ptr a -> IO Word32+peekUnalnWord32LE = fmap byteSwap32 . peek . castPtr++pokeUnalnWord32LE :: Ptr a -> Word32 -> IO ()+pokeUnalnWord32LE p w = poke (castPtr p) (byteSwap32 w)++#elif defined(MEM_UNALIGNED_OPS) && !defined(WORDS_BIGENDIAN)+peekUnalnWord16LE :: Ptr a -> IO Word16+peekUnalnWord16LE = peek . castPtr++peekUnalnWord32LE :: Ptr a -> IO Word32+peekUnalnWord32LE = peek . castPtr++pokeUnalnWord32LE :: Ptr a -> Word32 -> IO ()+pokeUnalnWord32LE p w = poke (castPtr p) w++#else+peekUnalnWord16LE :: Ptr a -> IO Word16+peekUnalnWord16LE p = do+    x <- fromIntegral <$> peekWord8 (plusPtr p 0)+    y <- fromIntegral <$> peekWord8 (plusPtr p 1)+    return $! x .|. unsafeShiftL y 8++peekUnalnWord32LE :: Ptr a -> IO Word32+peekUnalnWord32LE p = do+    x <- fromIntegral <$> peekWord8 (plusPtr p 0)+    y <- fromIntegral <$> peekWord8 (plusPtr p 1)+    z <- fromIntegral <$> peekWord8 (plusPtr p 2)+    w <- fromIntegral <$> peekWord8 (plusPtr p 3)+    return $! x .|. unsafeShiftL y 8 .|. unsafeShiftL z 16 .|. unsafeShiftL w 24++pokeUnalnWord32LE :: Ptr a -> Word32 -> IO ()+pokeUnalnWord32LE p w = do pokeByteOff p 0 (fromIntegral $ shiftR w  0 :: Word8)+                           pokeByteOff p 1 (fromIntegral $ shiftR w  8 :: Word8)+                           pokeByteOff p 2 (fromIntegral $ shiftR w 16 :: Word8)+                           pokeByteOff p 3 (fromIntegral $ shiftR w 24 :: Word8)+#endif+++#if defined(MEM_UNALIGNED_OPS) && !defined(WORDS_BIGENDIAN) && defined(HAVE_BYTESWAP_PRIMOPS)+peekUnalnWord16BE :: Ptr a -> IO Word16+peekUnalnWord16BE = fmap byteSwap16 . peek . castPtr++peekUnalnWord32BE :: Ptr a -> IO Word32+peekUnalnWord32BE = fmap byteSwap32 . peek . castPtr++#elif defined(MEM_UNALIGNED_OPS) && defined(WORDS_BIGENDIAN)+peekUnalnWord16BE :: Ptr a -> IO Word16+peekUnalnWord16BE = peek . castPtr++peekUnalnWord32BE :: Ptr a -> IO Word32+peekUnalnWord32BE = peek . castPtr++#else+peekUnalnWord16BE :: Ptr a -> IO Word16+peekUnalnWord16BE p = do+    x <- fromIntegral <$> peekWord8 (plusPtr p 0)+    y <- fromIntegral <$> peekWord8 (plusPtr p 1)+    return $! y .|. unsafeShiftL x 8++peekUnalnWord32BE :: Ptr a -> IO Word32+peekUnalnWord32BE p = do+    x <- fromIntegral <$> peekWord8 (plusPtr p 0)+    y <- fromIntegral <$> peekWord8 (plusPtr p 1)+    z <- fromIntegral <$> peekWord8 (plusPtr p 2)+    w <- fromIntegral <$> peekWord8 (plusPtr p 3)+    return $! w .|. unsafeShiftL z 8 .|. unsafeShiftL y 16 .|. unsafeShiftL x 24+#endif+