bbi (empty) → 0.1.0
raw patch · 7 files changed
+531/−0 lines, 7 filesdep +basedep +bytestringdep +cerealsetup-changed
Dependencies added: base, bytestring, cereal, conduit, containers, mtl, zlib
Files
- LICENSE +30/−0
- Setup.hs +2/−0
- bbi.cabal +37/−0
- src/Data/BBI.hs +248/−0
- src/Data/BBI/BigBed.hs +24/−0
- src/Data/BBI/BigWig.hs +117/−0
- src/Data/BBI/Utils.hs +73/−0
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright (c) 2014, Kai Zhang++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++ * Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.++ * Redistributions in binary form must reproduce the above+ copyright notice, this list of conditions and the following+ disclaimer in the documentation and/or other materials provided+ with the distribution.++ * Neither the name of Kai Zhang nor the names of other+ contributors may be used to endorse or promote products derived+ from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ bbi.cabal view
@@ -0,0 +1,37 @@+-- Initial bbi.cabal generated by cabal init. For further documentation, +-- see http://haskell.org/cabal/users-guide/++name: bbi+version: 0.1.0+synopsis: Tools for reading Big Binary Indexed files, e.g., bigBed, bigWig+-- description: +license: BSD3+license-file: LICENSE+author: Kai Zhang+maintainer: kai@kzhang.org+copyright: (c) 2014 Kai Zhang+category: Bio+build-type: Simple+-- extra-source-files: +cabal-version: >=1.10++library+ exposed-modules:+ Data.BBI+ , 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+ , mtl+ , containers+ , zlib++ hs-source-dirs: src++ default-language: Haskell2010
+ src/Data/BBI.hs view
@@ -0,0 +1,248 @@+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE Rank2Types #-}++module Data.BBI+ ( BbiFile(..)+ , BbiFileHeader(..)+ , FileType(..)+ , openBbiFile+ , closeBbiFile+ , getBbiFileHeader+ , getChromId+ , overlappingBlocks+ , readBlocks+ ) where++import Codec.Compression.Zlib (decompress)+import Control.Applicative ((<$>))+import Control.Monad.State+import qualified Data.ByteString as B+import qualified Data.ByteString.Lazy as BL+import Data.Conduit+import Data.Maybe+import qualified Data.Map.Strict as M+import System.IO+import Data.BBI.Utils++data FileType = BigBed+ | BigWig+ deriving (Show)++data BbiFile = BbiFile+ { _filehandle :: Handle+ , _header :: !BbiFileHeader+ , _chromTree :: !(M.Map B.ByteString (Int, Int))+ , _rtree :: !RTreeHeader+ } deriving (Show)++data BbiFileHeader = BbiFileHeader+ { _filetype :: !FileType+ , _endian :: !Endianness+ , _version :: !Int+ , _zoomlevels :: !Int+ , _chrTreeOffset :: !Int+ , _fullDataOffset :: !Int+ , _fullIndexOffset :: !Int+ , _filedCount :: !Int+ , _definedFiledCount :: !Int+ , _autoSqlOffset :: !Int+ , _totalSummaryOffset :: !Int+ , _uncompressBufSize :: !Int+ , _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 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++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)++data RTreeHeader = RTreeHeader+ { _blockSize :: Int+ , _itemCount :: Int+ , _startChromIx :: Int+ , _startBase :: Int+ , _endChromIx :: Int+ , _endBase :: Int+ , _endFileOffset :: Int+ , _itemsPerSlot :: Int+ , _rTreeReserved :: Int+ } deriving (Show)++getRTreeHeader :: Endianness -> Handle -> Integer -> IO (Either String RTreeHeader)+getRTreeHeader e h i = do+ hSeek h AbsoluteSeek i+ magic <- hReadInt32 e h+ if magic /= rtmagic+ then return $ Left "incorrect RTree magic"+ 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+ where+ rtmagic = 0x2468ACE0++overlap :: Int -> Int -> Int -> Int -> Int -> Int -> Int -> Bool+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 fl chr = fmap fst . M.lookup chr . _chromTree $ fl+{-# INLINE getChromId #-}++overlappingBlocks :: BbiFile -> (Int, Int, Int) -> IO [(Int, Int)]+overlappingBlocks fl (cid, start, end) =+ if overlap cid start end (_startChromIx $ _rtree fl)+ (_startBase $ _rtree fl)+ (_endChromIx $ _rtree fl)+ (_endBase $ _rtree fl)+ then catMaybes <$> go ((fromIntegral . _fullIndexOffset . _header) fl + 48)+ else return []+ where+ h = _filehandle fl+ e = _endian . _header $ fl+ go i = do+ hSeek h AbsoluteSeek i+ (isLeaf, n) <- readNode + if isLeaf+ then replicateM n readLeafItem+ else do + pos <- fromIntegral <$> hTell h+ results <- mapM (readNonLeafItem . (\x -> pos + x*24)) [0..n-1]+ return $ concat results++ readNode = do isLeaf <- hReadBool h+ _ <- hReadBool h+ count <- hReadInt16 e h+ return (isLeaf, count)++ readLeafItem = do+ stCIx <- hReadInt32 e h+ st <- hReadInt32 e h+ edCIx <- hReadInt32 e h+ ed <- hReadInt32 e h+ offset <- hReadInt64 e h+ size <- hReadInt64 e h+ return $ if overlap cid start end stCIx st edCIx ed+ then Just (offset, size)+ else Nothing++ readNonLeafItem p = do + hSeek h AbsoluteSeek $ fromIntegral p+ stCIx <- hReadInt32 e h+ st <- hReadInt32 e h+ edCIx <- hReadInt32 e h+ ed <- hReadInt32 e h+ next <- hReadInt64 e h+ if overlap cid start end stCIx st edCIx ed+ then go $ fromIntegral next+ else return []+{-# INLINE overlappingBlocks #-}++readBlocks :: BbiFile -> [(Int, Int)] -> Source IO B.ByteString+readBlocks fl blks = forM_ blks $ \(offset, size) -> do+ bs <- lift $ do hSeek handle AbsoluteSeek $ fromIntegral offset+ BL.hGet handle size+ yield . BL.toStrict . decompress $ bs+ where+ handle = _filehandle fl+{-# INLINE readBlocks #-}++{-+streamBbi :: BbiFile -> Source IO (B.ByteString, Int, Int, B.ByteString)+streamBbi fl = mapM_ (query fl) allChroms+ where+ allChroms = map (\(chr, (_, size)) -> (chr, 0, size-1)) . M.toList . _chromTree $ fl+ -}
+ src/Data/BBI/BigBed.hs view
@@ -0,0 +1,24 @@+module Data.BBI.BigBed where++import qualified Data.ByteString as B+import Data.Conduit+import qualified Data.Conduit.List as CL+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+ 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 #-}
+ src/Data/BBI/BigWig.hs view
@@ -0,0 +1,117 @@+module Data.BBI.BigWig+ ( BWFile+ , openBWFile+ , closeBWFile+ , queryBWFile+ ) where++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 Data.Maybe+import Data.BBI+import Data.BBI.Utils++newtype BWFile = BWFile BbiFile++openBWFile :: FilePath -> IO BWFile+openBWFile fl = do bbi <- openBbiFile fl+ case (_filetype . _header) bbi of+ BigWig -> return $ BWFile bbi+ _ -> error "not a bigWig file"++closeBWFile :: BWFile -> IO ()+closeBWFile (BWFile fl) = closeBbiFile fl++-- | wig section type+data SectionType = VarStep -- 0-indexed in binary+ | FixedStep -- 0-indexed in binary+ | BedGraph -- 0-indexed+ deriving (Show)++data WigHeader = WigHeader+ { _chromId :: !Int+ , _chromStart :: !Int+ , _chromEnd :: !Int+ , _itemStep :: !Int+ , _itemSpan :: !Int+ , _type :: !SectionType+ , _reserve :: !Int+ , _itemCount :: !Int+ } deriving (Show)++-- | the query should be 0-indexed+queryBWFile :: BWFile -> (B.ByteString, Int, Int) -> Source IO (B.ByteString, Int, Int, Double)+queryBWFile (BWFile fl) (chr, start, end) = unless (isNothing cid) $ do+ blks <- lift $ overlappingBlocks fl (fromJust cid, start, end)+ readBlocks fl blks $= toWigRecords endi $= filter'+ where+ filter' = CL.mapMaybe $ \(c, s, e, v) ->+ let s' = max start s+ e' = min end e+ in if c == fromJust cid && s' < e'+ then Just (chr, s', e', realToFrac v)+ else Nothing+ cid = getChromId fl chr+ endi = _endian . _header $ fl+{-# INLINE queryBWFile #-}++-- | 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 ->+ let header = readWigHeader endi bs+ in case _type header of+ VarStep -> readVarStep header $ B.drop 24 bs+ FixedStep -> readFixedStep header $ B.drop 24 bs+ BedGraph -> readBedGraph header $ B.drop 24 bs+ where+ readVarStep h = loop 0+ where+ loop i x | i >= n = []+ | otherwise = let s = (readInt32 endi . B.take 4) x+ e = s + sp+ v = readFloat32 endi . B.take 4 . B.drop 4 $ x+ in (cid, s, e, v) : loop (i+1) (B.drop 8 x)+ n = _itemCount h+ sp = _itemSpan h+ cid = _chromId h++ readFixedStep h = loop start 0+ 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) + n = _itemCount h+ sp = _itemSpan h+ step = _itemStep h+ start = _chromStart h+ cid = _chromId h++ readBedGraph h = loop 0+ where+ loop i x | i >= n = []+ | otherwise = let s = readInt32 endi . B.take 4 $ x+ e = readInt32 endi . B.take 4 . B.drop 4 $ x+ v = readFloat32 endi . B.take 4 . B.drop 8 $ x+ in (cid, s, e, v) : loop (i+1) (B.drop 12 x)+ n = _itemCount h+ cid = _chromId h+{-# INLINE toWigRecords #-}++readWigHeader :: Endianness -> B.ByteString -> WigHeader+readWigHeader e s = WigHeader (readInt32 e . B.take 4 $ s)+ (readInt32 e . B.take 4 . B.drop 4 $ s)+ (readInt32 e . B.take 4 . B.drop 8 $ s)+ (readInt32 e . B.take 4 . B.drop 12 $ s)+ (readInt32 e . B.take 4 . B.drop 16 $ s)+ (f . readInt8 . B.take 1 . B.drop 20 $ s)+ (readInt8 . B.take 1 . B.drop 21 $ s)+ (readInt16 e . B.take 2 . B.drop 22 $ s)+ where+ f x | x == 1 = BedGraph+ | x == 2 = VarStep+ | x == 3 = FixedStep+ | otherwise = error "Unknown Wig seciton type" +{-# INLINE readWigHeader #-}
+ src/Data/BBI/Utils.hs view
@@ -0,0 +1,73 @@+module Data.BBI.Utils+ ( Endianness(..)+ , hReadBool+ , readInt64+ , hReadInt64+ , readInt32+ , hReadInt32+ , readInt16+ , hReadInt16+ , readInt8+ , hReadInt8+ , readFloat32+ , fromRight+ ) where++import qualified Data.ByteString as B+import Data.Serialize+import System.IO++data Endianness = LE+ | BE+ deriving (Show)++hReadBool :: Handle -> IO Bool+hReadBool h = do+ bs <- B.hGet h 1+ return . fromRight . runGet (fmap (toEnum . fromIntegral) getWord8) $ bs+{-# INLINE hReadBool #-}++readInt64 :: Endianness -> B.ByteString -> Int+readInt64 LE = fromIntegral . fromRight . runGet getWord64le+readInt64 BE = fromIntegral . fromRight . runGet getWord64be+{-# INLINE readInt64 #-}++hReadInt64 :: Endianness -> Handle -> IO Int+hReadInt64 e h = fmap (readInt64 e) . B.hGet h $ 8+{-# INLINE hReadInt64 #-}++readInt32 :: Endianness -> B.ByteString -> Int+readInt32 LE = fromIntegral . fromRight . runGet getWord32le+readInt32 BE = fromIntegral . fromRight . runGet getWord32be+{-# INLINE readInt32 #-}++hReadInt32 :: Endianness -> Handle -> IO Int+hReadInt32 e h = fmap (readInt32 e) . B.hGet h $ 4+{-# INLINE hReadInt32 #-}++readInt16 :: Endianness -> B.ByteString -> Int+readInt16 LE = fromIntegral . fromRight . runGet getWord16le+readInt16 BE = fromIntegral . fromRight . runGet getWord16be+{-# INLINE readInt16 #-}++hReadInt16 :: Endianness -> Handle -> IO Int+hReadInt16 e h = fmap (readInt16 e) . B.hGet h $ 2+{-# INLINE hReadInt16 #-}++hReadInt8 :: Handle -> IO Int+hReadInt8 h = fmap readInt8 . B.hGet h $ 1+{-# INLINE hReadInt8 #-}++readInt8 :: B.ByteString -> Int+readInt8 = fromIntegral . fromRight . runGet getWord8+{-# INLINE readInt8 #-}++readFloat32 :: Endianness -> B.ByteString -> Float+readFloat32 LE = fromRight . runGet getFloat32le+readFloat32 BE = fromRight . runGet getFloat32be+{-# INLINE readFloat32 #-}++fromRight :: Either String b -> b+fromRight (Right x) = x+fromRight (Left s) = error s+{-# INLINE fromRight #-}