bbi 0.1.0 → 0.1.1

raw patch · 6 files changed

+293/−152 lines, 6 filesdep +bbidep +bioinformatics-toolkitdep +randomdep ~basedep ~bytestringdep ~conduitPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependencies added: bbi, bioinformatics-toolkit, random, tasty, tasty-golden, tasty-hunit, vector

Dependency ranges changed: base, bytestring, conduit

API changes (from Hackage documentation)

Files

• bbi.cabal +29/−9
• src/Data/BBI.hs +165/−119
• src/Data/BBI/BigBed.hs +54/−18
• src/Data/BBI/BigWig.hs +5/−6
• tests/Tests/BigBed.hs +33/−0
• tests/tests.hs +7/−0

@@ -1,18 +1,14 @@--- Initial bbi.cabal generated by cabal init.  For further documentation, --- see http://haskell.org/cabal/users-guide/- name:                bbi-version:             0.1.0+version:             0.1.1 synopsis:            Tools for reading Big Binary Indexed files, e.g., bigBed, bigWig--- description:         +-- description: license:             BSD3 license-file:        LICENSE author:              Kai Zhang maintainer:          kai@kzhang.org-copyright:           (c) 2014 Kai Zhang+copyright:           (c) 2014-2019 Kai Zhang category:            Bio build-type:          Simple--- extra-source-files:   cabal-version:       >=1.10  library@@ -21,13 +17,12 @@     , Data.BBI.BigBed     , Data.BBI.BigWig     , Data.BBI.Utils-  -- other-modules:           build-depends:       base >=4.7 && <5.0     , bytestring >=0.10.0.0     , cereal >=0.4.0-    , conduit >=1.2.0.0+    , conduit     , mtl     , containers     , zlib@@ -35,3 +30,28 @@   hs-source-dirs:      src    default-language:    Haskell2010++test-suite tests+  type: exitcode-stdio-1.0+  hs-source-dirs: tests+  main-is: tests.hs+  other-modules:+      Tests.BigBed++  default-language:    Haskell2010+  build-depends:+      base+    , bytestring+    , conduit+    , random+    , vector+    , tasty+    , tasty-golden+    , tasty-hunit+    , mtl+    , bioinformatics-toolkit+    , bbi++source-repository  head+  type: git+  location: https://github.com/kaizhang/bbi.git

@@ -1,5 +1,6 @@ {-# LANGUAGE OverloadedStrings #-} {-# LANGUAGE Rank2Types #-}+{-# LANGUAGE LambdaCase #-}  module Data.BBI     ( BbiFile(..)@@ -18,7 +19,7 @@ import Control.Monad.State import qualified Data.ByteString as B import qualified Data.ByteString.Lazy as BL-import Data.Conduit+import Conduit import Data.Maybe import qualified Data.Map.Strict as M import System.IO@@ -28,13 +29,45 @@               | BigWig     deriving (Show) -data BbiFile = BbiFile-    { _filehandle :: Handle-    , _header :: !BbiFileHeader-    , _chromTree :: !(M.Map B.ByteString (Int, Int))-    , _rtree :: !RTreeHeader-    } deriving (Show)+-- | Get the file type.+getFileType :: Handle -> IO (FileType, Endianness)+getFileType h = do+    hSeek h AbsoluteSeek 0+    magic <- B.hGet h 4+    let magicBE = readInt32 BE magic+        magicLE = readInt32 LE magic+    return $ case () of+        _ | magicBE == bigBedMagic -> (BigBed, BE)+          | magicBE == bigWigMagic -> (BigWig, BE)+          | magicLE == bigBedMagic -> (BigBed, LE)+          | magicLE == bigWigMagic -> (BigWig, LE)+          | otherwise -> error "Not a bigBed or bigWig file!"+  where+    bigWigMagic = 0x888FFC26+    bigBedMagic = 0x8789F2EB+{-# INLINE getFileType #-} ++-- | The header for big index file. Here are the Byte-by-byte details:+-- --------------------  ----  ----  -------------------------------------------+-- Name                  Size  Type	 Description+-- --------------------  ----  ----  -------------------------------------------+-- magic                 4     uint  0x888FFC26 for BigWig, 0x8789F2EB for+--                                   BigBed.  If byte-swapped, all+--                                   numbers in file are byte-swapped.+-- version               2     uint	 Currently 3. +-- zoomLevels            2     uint	 Number of different zoom summary resolutions.+-- chromosomeTreeOffset	 8     uint	 Offset in file to chromosome B+ tree index.+-- fullDataOffset        8     uint	 Offset to main (unzoomed) data. Points+--                                   specifically to the dataCount.+-- fullIndexOffset       8     uint	 Offset to R tree index of items.+-- fieldCount            2     uint	 Number of fields in BED file. (0 for BigWig)+-- definedFieldCount     2     uint	 Number of fields that are predefined BED fields.+-- autoSqlOffset         8     uint	 Offset to zero-terminated string with .as spec.+-- totalSummaryOffset    8     uint	 Offset to overall file summary data block.+-- uncompressBufSize     4     uint	 Maximum size of decompression buffer+--                                   needed (nonzero on compressed files).+-- reserved              8     uint	 Reserved for future expansion. Currently 0. data BbiFileHeader = BbiFileHeader     { _filetype :: !FileType     , _endian :: !Endianness@@ -51,96 +84,30 @@     , _reserved :: !Int     } deriving (Show) -openBbiFile :: FilePath -> IO BbiFile-openBbiFile fl = do-    h <- openFile fl ReadMode-    header <- fromRight <$> getBbiFileHeader h-    t <- fromRight <$> getChromTreeAsList (_endian header) h (fromIntegral $ _chrTreeOffset header)-    rtree <- fromRight <$> getRTreeHeader (_endian header) h (fromIntegral $ _fullIndexOffset header)-    return $ BbiFile h header (M.fromList t) rtree--closeBbiFile :: BbiFile -> IO ()-closeBbiFile = hClose . _filehandle--getBbiFileHeader :: Handle -> IO (Either String BbiFileHeader)+getBbiFileHeader :: Handle -> IO BbiFileHeader getBbiFileHeader h = do-    ft <- getFileType h-    case ft of-        Right (t, e) -> do-            hSeek h AbsoluteSeek 4-            v <- hReadInt16 e h-            zl <- hReadInt16 e h-            ct <- hReadInt64 e h-            fd <- hReadInt64 e h-            fi <- hReadInt64 e h-            fc <- hReadInt16 e h-            df <- hReadInt16 e h-            as <- hReadInt64 e h-            ts <- hReadInt64 e h-            ub <- hReadInt32 e h-            r <- hReadInt64 e h-            return . Right $ BbiFileHeader t e v zl ct fd fi fc df as ts ub r-        Left e -> return $ Left e+    (ft, endi) <- getFileType h +    hSeek h AbsoluteSeek 4+    BbiFileHeader <$> return ft <*> return endi <*> hReadInt16 endi h <*>+        hReadInt16 endi h <*> hReadInt64 endi h <*> hReadInt64 endi h <*>+        hReadInt64 endi h <*> hReadInt16 endi h <*> hReadInt16 endi h <*>+        hReadInt64 endi h <*> hReadInt64 endi h <*> hReadInt32 endi h <*>+        hReadInt64 endi h+{-# INLINE getBbiFileHeader #-} -getFileType :: Handle -> IO (Either String (FileType, Endianness))-getFileType h = do-    hSeek h AbsoluteSeek 0-    magic <- B.hGet h 4-    let magicBE = readInt32 BE magic-        magicLE = readInt32 LE magic-    return $ case () of-        _ | magicBE == bigBedMagic -> Right (BigBed, BE)-          | magicBE == bigWigMagic -> Right (BigWig, BE)-          | magicLE == bigBedMagic -> Right (BigBed, LE)-          | magicLE == bigWigMagic -> Right (BigWig, LE)-          | otherwise -> Left "not a bigBed/bigWig file"-  where-    bigWigMagic = 0x888FFC26-    bigBedMagic = 0x8789F2EB-{-# INLINE getFileType #-} -getChromTreeAsList :: Endianness -> Handle -> Integer -> IO (Either String [(B.ByteString, (Int, Int))])-getChromTreeAsList e h offset = do-    hSeek h AbsoluteSeek offset-    magic <- hReadInt32 e h-    if magic /= chromTreeMagic-       then return $ Left "wrong chromosome tree header"-       else do-           _ <- hReadInt32 e h -- blockSize, not used-           keySize <- hReadInt32 e h-           _ <- hReadInt32 e h -- valSize, not used-           _ <- hReadInt64 e h -- itemCount, not used-           _ <- hReadInt64 e h -- reserved, not used-           Right <$> traverseTree keySize-  where-    traverseTree ks = go-      where-        go = do (isLeaf, n) <- readNode-                if isLeaf-                   then replicateM n $ readLeafItem ks-                   else do-                       xs <- replicateM n $ readNonLeafItem ks-                       rs <- forM xs $ \(_, off) -> do-                           hSeek h AbsoluteSeek $ fromIntegral off-                           go-                       return $ concat rs--    chromTreeMagic = 0x78CA8C91--    readNode = do isLeaf <- hReadBool h-                  _ <- hReadBool h-                  count <- hReadInt16 e h-                  return (isLeaf, count)--    readLeafItem n = do key <- B.filter (/=0) <$> B.hGet h n-                        chromId <- hReadInt32 e h-                        chromSize <- hReadInt32 e h-                        return (key, (chromId, chromSize))--    readNonLeafItem n = do key <- B.filter (/=0) <$> B.hGet h n-                           childOffset <- hReadInt64 e h-                           return (key, childOffset)-+-- | R tree index header.+-- Name Size	Type	Description+-- magic	4	uint	0x2468ACE0.  If byte-swapped all numbers in index are byte-swapped.+-- blockSize	4	uint	Number of children per block (not byte size of block). +-- itemCount	8	uint	The number of chromosomes/contigs.+-- startChromIx	4	uint	ID of first chromosome in index.+-- startBase	4	uint	Position of first base in index.+-- endChromIx	4	uint	ID of last chromosome in index.+-- endBase	4	uint	Position of last base in index.+-- endFileOffset	8	uint	Position in file where data being indexed ends.+-- itemsPerSlot	4	uint	Number of items pointed to by leaves of index.+-- Reserved	4	uint	Reserved for future expansion. Currently 0. data RTreeHeader = RTreeHeader     { _blockSize :: Int     , _itemCount :: Int@@ -153,36 +120,110 @@     , _rTreeReserved :: Int     } deriving (Show) -getRTreeHeader :: Endianness -> Handle -> Integer -> IO (Either String RTreeHeader)+getRTreeHeader :: Endianness+               -> Handle+               -> Integer   -- ^ Location of R Tree header in the file.+               -> IO RTreeHeader getRTreeHeader e h i = do-    hSeek h AbsoluteSeek i-    magic <- hReadInt32 e h-    if magic /= rtmagic-       then return $ Left "incorrect RTree magic"+    magic <- hSeek h AbsoluteSeek i >> hReadInt32 e h+    if magic /= 0x2468ACE0+       then error "incorrect RTree magic"+       else RTreeHeader <$> hReadInt32 e h <*> hReadInt64 e h <*>+           hReadInt32 e h <*> hReadInt32 e h <*> hReadInt32 e h <*>+           hReadInt32 e h <*> hReadInt64 e h <*> hReadInt32 e h <*>+           hReadInt32 e h+++-- | Big index file format. Here are the Byte-by-byte details:+-- --------------  ------  -----------------------------------------------------+-- Name            Size    Description+-- --------------  ------  -----------------------------------------------------+-- bbiHeader       64      Contains high-level information about file and+--                         offsets to various parts of file.+-- zoomHeaders     N*24    One for each level of zoom built into file.+-- autoSql         Varies  Zero-terminated string in autoSql format describing+--                         formats.  Optional, not used in BigWig.+-- totalSummary    40      Statistical summary of entire file.+-- chromosomeTree  Varies  B+ tree-formatted index of chromosomes, their sizes,+--                         and a unique ID for each.+-- dataCount       4       Number of records in data. For BigWig this+--                         corresponds to the number of sections, not the+--                         number of floating point values.+-- data            Varies  Possibly compressed data in format specific for+--                         file type.+-- index           Varies  R tree index of data.+-- zoomInfo        Varies  One for each zoom level.+data BbiFile = BbiFile+    { _filehandle :: Handle+    , _header :: !BbiFileHeader+    , _chromTree :: !(M.Map Chromosome (ChromID, ChromSize))+    , _rtree :: !RTreeHeader+    } deriving (Show)++type Chromosome = B.ByteString+type ChromSize = Int+type ChromID = Int++openBbiFile :: FilePath -> IO BbiFile+openBbiFile fl = do+    h <- openFile fl ReadMode+    header <- getBbiFileHeader h+    t <- getChromTreeAsList (_endian header) h (fromIntegral $ _chrTreeOffset header)+    rtree <- getRTreeHeader (_endian header) h (fromIntegral $ _fullIndexOffset header)+    return $ BbiFile h header (M.fromList t) rtree++closeBbiFile :: BbiFile -> IO ()+closeBbiFile = hClose . _filehandle++getChromTreeAsList :: Endianness+                   -> Handle+                   -> Integer   -- ^ Location+                   -> IO [(B.ByteString, (Int, Int))]+getChromTreeAsList e h offset = do+    magic <- hSeek h AbsoluteSeek offset >> hReadInt32 e h+    if magic /= 0x78CA8C91+       then error "wrong chromosome tree header"        else do-           bs <- hReadInt32 e h-           ic <- hReadInt64 e h-           sc <- hReadInt32 e h-           sb <- hReadInt32 e h-           ec <- hReadInt32 e h-           eb <- hReadInt32 e h-           eo <- hReadInt64 e h-           ips <- hReadInt32 e h-           r <- hReadInt32 e h-           return $ Right $ RTreeHeader bs ic sc sb ec eb eo ips r+           _ <- hReadInt32 e h -- blockSize, not used+           keySize <- hReadInt32 e h+           _ <- hReadInt32 e h -- valSize, not used+           _ <- hReadInt64 e h -- itemCount, not used+           _ <- hReadInt64 e h -- reserved, not used+           traverseTree keySize   where-    rtmagic = 0x2468ACE0+    traverseTree ks = go+      where+        go = readNode >>= \case+            (True, n) -> replicateM n $ readLeafItem ks+            (False, n) -> fmap concat $+                replicateM n (readNonLeafItem ks) >>= mapM ( \(_, off) -> do+                    hSeek h AbsoluteSeek $ fromIntegral off+                    go )+    readNode = do+        isLeaf <- hReadBool h+        _ <- hReadBool h+        count <- hReadInt16 e h+        return (isLeaf, count)+    readLeafItem n = do+        key <- B.filter (/=0) <$> B.hGet h n+        chromId <- hReadInt32 e h+        chromSize <- hReadInt32 e h+        return (key, (chromId, chromSize))+    readNonLeafItem n = do+        key <- B.filter (/=0) <$> B.hGet h n+        childOffset <- hReadInt64 e h+        return (key, childOffset)  overlap :: Int -> Int -> Int -> Int -> Int -> Int -> Int -> Bool-overlap qChr qStart qEnd rStartChr rStart rEndChr rEnd = +overlap qChr qStart qEnd rStartChr rStart rEndChr rEnd =     (qChr, qStart) < (rEndChr, rEnd) && (qChr, qEnd) > (rStartChr, rStart) {-# INLINE overlap #-} -getChromId :: BbiFile -> B.ByteString -> Maybe Int+getChromId :: BbiFile -> B.ByteString -> Maybe ChromID getChromId fl chr = fmap fst . M.lookup chr . _chromTree $ fl {-# INLINE getChromId #-} -overlappingBlocks :: BbiFile -> (Int, Int, Int) -> IO [(Int, Int)]+overlappingBlocks :: BbiFile -> (ChromID, Int, Int) -> IO [Block] overlappingBlocks fl (cid, start, end) =     if overlap cid start end (_startChromIx $ _rtree fl)                              (_startBase $ _rtree fl)@@ -195,10 +236,10 @@     e = _endian . _header $ fl     go i = do         hSeek h AbsoluteSeek i-        (isLeaf, n) <- readNode +        (isLeaf, n) <- readNode         if isLeaf            then replicateM n readLeafItem-           else do +           else do                pos <- fromIntegral <$> hTell h                results <- mapM (readNonLeafItem . (\x -> pos + x*24)) [0..n-1]                return $ concat results@@ -219,7 +260,7 @@                     then Just (offset, size)                     else Nothing -    readNonLeafItem p = do +    readNonLeafItem p = do         hSeek h AbsoluteSeek $ fromIntegral p         stCIx <- hReadInt32 e h         st <- hReadInt32 e h@@ -231,10 +272,15 @@            else return [] {-# INLINE overlappingBlocks #-} -readBlocks :: BbiFile -> [(Int, Int)] -> Source IO B.ByteString+type Block = (Offset, BlockSize)+type Offset = Int+type BlockSize = Int++-- | Read Blocks from the file. +readBlocks :: BbiFile -> [Block] -> Source IO B.ByteString readBlocks fl blks = forM_ blks $ \(offset, size) -> do-    bs <- lift $ do hSeek handle AbsoluteSeek $ fromIntegral offset-                    BL.hGet handle size+    bs <- lift $ hSeek handle AbsoluteSeek (fromIntegral offset) >>+        BL.hGet handle size     yield . BL.toStrict . decompress $ bs   where     handle = _filehandle fl

@@ -1,24 +1,60 @@-module Data.BBI.BigBed where+module Data.BBI.BigBed+    ( BBedFile+    , openBBedFile+    , closeBBedFile+    , query+    ) where  import qualified Data.ByteString as B-import Data.Conduit-import qualified Data.Conduit.List as CL+import qualified Data.ByteString.Char8 as BC+import qualified Data.Map as M+import Conduit+ import Data.BBI import Data.BBI.Utils -toBedRecords :: Monad m => Endianness -> Conduit B.ByteString m (Int, Int, Int, B.ByteString)-toBedRecords e = do-    d <- await-    case d of-        Nothing -> return ()-        Just bs -> go bs+newtype BBedFile = BBedFile BbiFile++openBBedFile :: FilePath -> IO BBedFile+openBBedFile fl = do+    bbi <- openBbiFile fl+    case (_filetype . _header) bbi of+        BigBed -> return $ BBedFile bbi+        _ -> error "not a bigbed file"++closeBBedFile :: BBedFile -> IO ()+closeBBedFile (BBedFile fl) = closeBbiFile fl++-- | Find all regions that are overlapped with the query in the bigbed file.+query :: (B.ByteString, Int, Int)   -- ^ (Chr, start, end)+      -> BBedFile+      -> ConduitT () (B.ByteString, Int, Int, B.ByteString) IO ()+query (chr, s, e) (BBedFile fl) = case getChromId fl chr of+    Just i -> do+        blks <- lift $ overlappingBlocks fl (i, s, e)+        readBlocks fl blks .| concatMapC (map f . decodeBlock endi) .| filterC+          (\(_,s',e',_) -> isOverlapped (s',e') (s,e))+    Nothing -> return ()   where-    go s | B.null s = return ()-         | otherwise = do-             let chr = readInt32 e . B.take 4 $ s-                 start = readInt32 e . B.take 4 . B.drop 4 $ s-                 end = readInt32 e . B.take 4 . B.drop 8 $ s-                 (rest, remain) = B.break (==0) . B.drop 12 $ s-             yield (chr, start, end, rest)-             go $ B.tail remain-{-# INLINE toBedRecords #-}+    endi = _endian $ _header fl+    id2Chr = M.fromList $ map (\(a,(b,_)) -> (b,a) ) $ M.toList $ _chromTree fl+    f (a,b,c,d) = (M.findWithDefault undefined a id2Chr, b, c, d)++-- | Whether two intervals are overlapped.+isOverlapped :: (Int, Int) -> (Int, Int) -> Bool+isOverlapped (lo1, hi1) (lo2, hi2) = not (lo2 > hi1 || lo1 > hi2)+{-# INLINE isOverlapped #-}++-- | Decode blocks in BigBed.+decodeBlock :: Endianness+            -> B.ByteString +            -> [(Int, Int, Int, B.ByteString)]+decodeBlock e blk+    | B.null blk = []+    | otherwise = (chr, start, end, rest) : decodeBlock e (B.tail remain)+  where+    chr = readInt32 e $ B.take 4 blk+    start = readInt32 e $ B.take 4 $ B.drop 4 blk+    end = readInt32 e $ B.take 4 $ B.drop 8 blk+    (rest, remain) = BC.break (=='\0') $ B.drop 12 blk+{-# INLINE decodeBlock #-}

@@ -8,8 +8,7 @@ import Control.Monad (unless) import Control.Monad.Trans (lift) import qualified Data.ByteString as B-import Data.Conduit-import qualified Data.Conduit.List as CL+import Conduit import Data.Maybe import Data.BBI import Data.BBI.Utils@@ -48,7 +47,7 @@     blks <- lift $ overlappingBlocks fl (fromJust cid, start, end)     readBlocks fl blks $= toWigRecords endi $= filter'   where-    filter' = CL.mapMaybe $ \(c, s, e, v) ->+    filter' = concatMapC $ \(c, s, e, v) ->                  let s' = max start s                      e' = min end e                  in if c == fromJust cid && s' < e'@@ -60,7 +59,7 @@  -- | convert bytestring to wig style record, but the output is 0-indexed toWigRecords :: Monad m => Endianness -> Conduit B.ByteString m (Int, Int, Int, Float)-toWigRecords endi = CL.concatMap $ \bs ->+toWigRecords endi = concatMapC $ \bs ->     let header = readWigHeader endi bs     in case _type header of         VarStep -> readVarStep header $ B.drop 24 bs@@ -82,7 +81,7 @@       where         loop st i x | i >= n = []                     | otherwise = let v = readFloat32 endi . B.take 4 $ x-                                  in (cid, st, st+sp, v) : loop (st+step) (i+1) (B.drop 4 x) +                                  in (cid, st, st+sp, v) : loop (st+step) (i+1) (B.drop 4 x)         n = _itemCount h         sp = _itemSpan h         step = _itemStep h@@ -113,5 +112,5 @@     f x | x == 1 = BedGraph         | x == 2 = VarStep         | x == 3 = FixedStep-        | otherwise = error "Unknown Wig seciton type" +        | otherwise = error "Unknown Wig seciton type" {-# INLINE readWigHeader #-}

@@ -0,0 +1,33 @@+{-# LANGUAGE OverloadedStrings #-}++module Tests.BigBed+    ( tests ) where++import Bio.Data.Bed+import Bio.Data.Bed.Types+import Conduit++import Test.Tasty+import Test.Tasty.HUnit++import Data.BBI.BigBed++tests :: TestTree+tests = testGroup "Test: BigBed"+    [ testCase "intersect" intersectTest+    ]++intersectTest :: Assertion+intersectTest = do+    bed <- readBed "tests/data/sample.bed"+    bbed <- openBBedFile "tests/data/sample.bigbed"+    compareIntersection (BED3 "chr1" 1 10000000) bed bbed++compareIntersection :: BED3 -> [BED3] -> BBedFile -> Assertion+compareIntersection x@(BED3 chr s e) bed bbed = do+    observed <- runConduit $ query (chr, s, e) bbed .| mapC f .| sinkList+    fromSorted (sortBed observed) @?= fromSorted (sortBed expected)+  where+    expected = runIdentity $ runConduit $+        yieldMany bed .| intersectBed [x] .| sinkList+    f (a,b,c,_) = BED3 a b c

@@ -0,0 +1,7 @@+import qualified Tests.BigBed as BB+import Test.Tasty++main :: IO ()+main = defaultMain $ testGroup "Main"+    [ BB.tests+    ]