packages feed

bioinformatics-toolkit 0.2.0 → 0.2.1

raw patch · 13 files changed

+492/−308 lines, 13 filesdep ~http-conduit

Dependency ranges changed: http-conduit

Files

app/viewSeq.hs view
@@ -8,8 +8,8 @@ main :: IO () main = do     [fl, chr, start, end] <- getArgs-    g <- pack fl-    s <- getSeqs g [(B.pack chr, read start, read end)] :: IO [Either String (DNA IUPAC)]-    case head s of-        Left err -> error err-        Right x -> B.putStrLn . toBS $ x+    withGenome fl $ \g -> do+        s <- getSeq g (B.pack chr, read start, read end) :: IO (Either String (DNA IUPAC))+        case s of+            Left err -> error err+            Right x -> B.putStrLn . toBS $ x
bioinformatics-toolkit.cabal view
@@ -2,7 +2,7 @@ -- further documentation, see http://haskell.org/cabal/users-guide/  name:                bioinformatics-toolkit-version:             0.2.0+version:             0.2.1 synopsis:            A collection of bioinformatics tools description:         A collection of bioinformatics tools license:             MIT@@ -63,7 +63,7 @@     , containers >=0.5     , data-default-class     , double-conversion-    , http-conduit+    , http-conduit >=2.1.8     , hexpat     , matrices >=0.4.3     , mtl >=2.1.3.1@@ -144,6 +144,7 @@     , Tests.ChIPSeq     , Tests.Motif     , Tests.Seq+    , Tests.GREAT    default-language:    Haskell2010   build-depends:
− data/1.fasta 
@@ -1,30 +0,0 @@-MEME version 4--ALPHABET= ACGT--strands: + ---Background letter frequencies-A 0.25 C 0.25 G 0.25 T 0.25--MOTIF tree0.2_0_GGAKATWATMTC-letter-probability matrix: alength= 4 w= 12 nsites= 1 E= 0- 0.001 0.001 0.997 0.001- 0.001 0.003 0.995 0.001- 0.997 0.001 0.001 0.001- 0.001 0.001 0.474 0.524- 0.997 0.001 0.001 0.001- 0.997 0.001 0.001 0.001- 0.25 0.25 0.25 0.25- 0.25 0.25 0.25 0.25--MOTIF tree0.2_1_CTCCCTAAAATG-letter-probability matrix: alength= 4 w= 12 nsites= 1 E= 0- 0.25 0.25 0.25 0.25- 0.25 0.25 0.25 0.25- 0.001 0.001 0.997 0.001- 0.001 0.003 0.995 0.001- 0.997 0.001 0.001 0.001- 0.001 0.001 0.474 0.524- 0.997 0.001 0.001 0.001- 0.997 0.001 0.001 0.001
src/Bio/ChIPSeq.hs view
@@ -79,7 +79,7 @@     vec <- lift $ GM.replicate l 0     n <- foldMC (count vec) (0 :: Int)     let factor = fromIntegral n / 1e9-    lift $ liftM (G.imap (\i x -> x / factor / (fromIntegral . size) (regions V.! i)))+    lift $ liftM (G.imap (\i x -> x / factor / (fromIntegral . bedSize) (regions V.! i)))          $ G.unsafeFreeze vec   where     count v nTags tag = do@@ -246,7 +246,7 @@     peaks' = map f peaks     f b = let chr = chrom b               c = (chromStart b + chromEnd b) `div` 2-          in asBed chr (c-r) (c+r)+          in BED3 chr (c-r) (c+r)     mergeFn xs = BED (chrom $ head xs) lo hi Nothing (Just $ fromIntegral $ length xs) Nothing       where         lo = minimum . map chromStart $ xs
src/Bio/Data/Bed.hs view
@@ -15,9 +15,19 @@  module Bio.Data.Bed     ( BEDLike(..)++    -- * BED6 format+    , BED(..)+    -- * BED3 format+    , BED3(..)+    -- * NarrowPeak format+    , NarrowPeak(..)+     , BEDTree     , bedToTree     , sortedBedToTree+    , intersecting+    , isIntersected     , splitBed     , splitBedBySize     , splitBedBySizeLeft@@ -33,14 +43,23 @@     , mergeBedWith     , mergeSortedBed     , mergeSortedBedWith-    -- * BED6 format-    , BED(..)-    -- * BED3 format-    , BED3(..)+    , splitOverlapped+    , hReadBed+    , hReadBed'+    , readBed+    , readBed'+    , hWriteBed+    , hWriteBed'+    , writeBed+    , writeBed'++    -- * Utilities     , fetchSeq     , fetchSeq'+    , motifScan+    , getMotifScore+    , getMotifPValue     , compareBed-    , convert     ) where  import Control.Arrow ((***))@@ -52,89 +71,203 @@ import qualified Data.Foldable as F import qualified Data.HashMap.Strict as M import qualified Data.IntervalMap.Strict as IM-import Data.List (groupBy)-import Data.Maybe (fromMaybe)+import Data.List (groupBy, sortBy)+import Data.Maybe (fromMaybe, fromJust) import qualified Data.Vector as V import qualified Data.Vector.Algorithms.Intro as I import System.IO+import Data.ByteString.Lex.Integral (packDecimal)+import Data.Double.Conversion.ByteString (toShortest)+import Data.Ord (comparing) +import Bio.Motif hiding (_name)+import Bio.Motif.Search import Bio.Seq import Bio.Seq.IO import Bio.Utils.Misc ( binBySizeLeft, binBySize, binBySizeOverlap, bins                       , readInt, readDouble ) --- | a class representing BED-like data, e.g., BED3, BED6 and BED12. BED format+-- | A class representing BED-like data, e.g., BED3, BED6 and BED12. BED format -- uses 0-based index (see documentation). class BEDLike b where-    -- | construct bed record from chromsomoe, start location and end location+    -- | Construct bed record from chromsomoe, start location and end location     asBed :: B.ByteString -> Int -> Int -> b -    -- | convert bytestring to bed format+    -- | Convert bytestring to bed format     fromLine :: B.ByteString -> b -    -- | convert bed to bytestring+    -- | Convert bed to bytestring     toLine :: b -> B.ByteString -    -- | field accessor+    -- | Field accessor     chrom :: b -> B.ByteString     chromStart :: b -> Int     chromEnd :: b -> Int+    bedName :: b -> Maybe B.ByteString+    bedScore :: b -> Maybe Double+    bedStrand :: b -> Maybe Bool -    toBed3 :: b -> BED3-    toBed3 bed = BED3 (chrom bed) (chromStart bed) (chromEnd bed)+    convert :: BEDLike b' => b' -> b+    convert bed = asBed (chrom bed) (chromStart bed) (chromEnd bed)+    {-# INLINE convert #-} -    -- | return the size of a bed region-    size :: b -> Int-    size bed = chromEnd bed - chromStart bed+    -- | Return the size of a bed region.+    bedSize :: b -> Int+    bedSize bed = chromEnd bed - chromStart bed -    hReadBed :: MonadIO m => Handle -> Source m b-    hReadBed h = do eof <- liftIO $ hIsEOF h-                    unless eof $ do-                        line <- liftIO $ B.hGetLine h-                        yield $ fromLine line-                        hReadBed h-    {-# INLINE hReadBed #-}+    {-# MINIMAL asBed, fromLine, toLine, chrom, chromStart, chromEnd,+                bedName, bedScore, bedStrand #-} -    -- | non-streaming version-    hReadBed' :: MonadIO m => Handle -> m [b]-    hReadBed' h = hReadBed h $$ sinkList-    {-# INLINE hReadBed' #-}+-- * BED6 format -    readBed :: MonadIO m => FilePath -> Source m b-    readBed fl = do handle <- liftIO $ openFile fl ReadMode-                    hReadBed handle-                    liftIO $ hClose handle-    {-# INLINE readBed #-}+-- | BED6 format, as described in http://genome.ucsc.edu/FAQ/FAQformat.html#format1.7+data BED = BED+    { _chrom :: !B.ByteString+    , _chromStart :: {-# UNPACK #-} !Int+    , _chromEnd :: {-# UNPACK #-} !Int+    , _name :: !(Maybe B.ByteString)+    , _score :: !(Maybe Double)+    , _strand :: !(Maybe Bool)  -- ^ True: "+", False: "-"+    } deriving (Eq, Show, Read) -    -- | non-streaming version-    readBed' :: MonadIO m => FilePath -> m [b]-    readBed' fl = readBed fl $$ sinkList-    {-# INLINE readBed' #-}+instance Default BED where+    def = BED+        { _chrom = ""+        , _chromStart = 0+        , _chromEnd = 0+        , _name = Nothing+        , _score = Nothing+        , _strand = Nothing+        } -    hWriteBed :: MonadIO m => Handle -> Sink b m ()-    hWriteBed handle = do-        x <- await-        case x of-            Nothing -> return ()-            Just bed -> (liftIO . B.hPutStrLn handle . toLine) bed >> hWriteBed handle-    {-# INLINE hWriteBed #-}+instance BEDLike BED where+    asBed chr s e = BED chr s e Nothing Nothing Nothing -    hWriteBed' :: MonadIO m => Handle -> [b] -> m ()-    hWriteBed' handle beds = yieldMany beds $$ hWriteBed handle-    {-# INLINE hWriteBed' #-}+    fromLine l = evalState (f (B.split '\t' l)) 1+      where+        f :: [B.ByteString] -> State Int BED+        f [] = do i <- get+                  if i <= 3 then error "Read BED fail: Incorrect number of fields"+                            else return def+        f (x:xs) = do+            i <- get+            put (i+1)+            bed <- f xs+            case i of+                1 -> return $ bed {_chrom = x}+                2 -> return $ bed {_chromStart = readInt x}+                3 -> return $ bed {_chromEnd = readInt x}+                4 -> return $ bed {_name = guard' x}+                5 -> return $ bed {_score = getScore x}+                6 -> return $ bed {_strand = getStrand x}+                _ -> return def -    writeBed :: MonadIO m => FilePath -> Sink b m ()-    writeBed fl = do handle <- liftIO $ openFile fl WriteMode-                     hWriteBed handle-                     liftIO $ hClose handle-    {-# INLINE writeBed #-}+        guard' x | x == "." = Nothing+                 | otherwise = Just x+        getScore x | x == "." = Nothing+                   | otherwise = Just . readDouble $ x+        getStrand str | str == "-" = Just False+                      | str == "+" = Just True+                      | otherwise = Nothing+    {-# INLINE fromLine #-} -    writeBed' :: MonadIO m => FilePath -> [b] -> m ()-    writeBed' fl beds = yieldMany beds $$ writeBed fl-    {-# INLINE writeBed' #-}+    toLine (BED f1 f2 f3 f4 f5 f6) = B.intercalate "\t"+        [ f1, (B.pack.show) f2, (B.pack.show) f3, fromMaybe "." f4, score'+        , strand' ]+      where+        strand' | f6 == Just True = "+"+                | f6 == Just False = "-"+                | otherwise = "."+        score' = case f5 of+                     Just x -> (B.pack.show) x+                     _ -> "."+    {-# INLINE toLine #-} -    {-# MINIMAL asBed, fromLine, toLine, chrom, chromStart, chromEnd #-}+    chrom = _chrom+    chromStart = _chromStart+    chromEnd = _chromEnd+    bedName = _name+    bedScore = _score+    bedStrand = _strand +    convert bed = BED (chrom bed) (chromStart bed) (chromEnd bed) (bedName bed)+                      (bedScore bed) (bedStrand bed)++-- * BED3 format++data BED3 = BED3 !B.ByteString !Int !Int deriving (Eq, Show, Read)++instance BEDLike BED3 where+    asBed = BED3++    fromLine l = case B.split '\t' l of+                    (a:b:c:_) -> BED3 a (readInt b) $ readInt c+                    _ -> error "Read BED fail: Incorrect number of fields"+    {-# INLINE fromLine #-}++    toLine (BED3 a b c) = B.intercalate "\t"+        [a, fromJust $ packDecimal b, fromJust $ packDecimal c]+    {-# INLINE toLine #-}++    chrom (BED3 f1 _ _) = f1+    chromStart (BED3 _ f2 _) = f2+    chromEnd (BED3 _ _ f3) = f3+    bedName = const Nothing+    bedScore = const Nothing+    bedStrand = const Nothing++-- | ENCODE narrowPeak format: https://genome.ucsc.edu/FAQ/FAQformat.html#format12+data NarrowPeak = NarrowPeak+    { _npChrom :: !B.ByteString+    , _npStart :: !Int+    , _npEnd :: !Int+    , _npName :: !(Maybe B.ByteString)+    , _npScore :: !Double+    , _npStrand :: !(Maybe Bool)+    , _npSigal :: !Double+    , _npPvalue :: !(Maybe Double)+    , _npQvalue :: !(Maybe Double)+    , _npPeak :: !(Maybe Int)+    } deriving (Eq, Show, Read)++instance BEDLike NarrowPeak where+    asBed chr s e = NarrowPeak chr s e Nothing 0 Nothing 0 Nothing Nothing Nothing++    fromLine l = NarrowPeak a (readInt b) (readInt c)+        (if d == "." then Nothing else Just d)+        (readDouble e)+        (if f == "." then Nothing else if f == "+" then Just True else Just False)+        (readDouble g)+        (if readDouble h < 0 then Nothing else Just $ readDouble h)+        (if readDouble i < 0 then Nothing else Just $ readDouble i)+        (if readInt j < 0 then Nothing else Just $ readInt j)+      where+        (a:b:c:d:e:f:g:h:i:j:_) = B.split '\t' l+    {-# INLINE fromLine #-}++    toLine (NarrowPeak a b c d e f g h i j) = B.intercalate "\t"+        [ a, fromJust $ packDecimal b, fromJust $ packDecimal c, fromMaybe "." d+        , toShortest e+        , case f of+            Nothing -> "."+            Just True -> "+"+            _ -> "-"+        , toShortest g, fromMaybe "-1" $ fmap toShortest h+        , fromMaybe "-1" $ fmap toShortest i+        , fromMaybe "-1" $ fmap (fromJust . packDecimal) j+        ]+    {-# INLINE toLine #-}++    chrom = _npChrom+    chromStart = _npStart+    chromEnd = _npEnd+    bedName = _npName+    bedScore = Just . _npScore+    bedStrand = _npStrand++    convert bed = NarrowPeak (chrom bed) (chromStart bed) (chromEnd bed) (bedName bed)+        (fromMaybe 0 $ bedScore bed) (bedStrand bed) 0 Nothing Nothing Nothing+ type BEDTree a = M.HashMap B.ByteString (IM.IntervalMap Int a)  -- | convert a set of bed records to interval tree, with combining function for@@ -170,6 +303,17 @@         return v {-# INLINE bedToTree #-} +intersecting :: BEDLike b => BEDTree a -> b -> IM.IntervalMap Int a+intersecting tree x = IM.intersecting (M.lookupDefault IM.empty chr tree) interval+  where+    chr = chrom x+    interval = IM.IntervalCO (chromStart x) $ chromEnd x+{-# INLINE intersecting #-}++isIntersected :: BEDLike b => BEDTree a -> b -> Bool+isIntersected tree = not . IM.null . intersecting tree+{-# INLINE isIntersected #-}+ -- | split a bed region into k consecutive subregions, discarding leftovers splitBed :: BEDLike b => Int -> b -> [b] splitBed k bed = map (uncurry (asBed chr)) . bins k $ (s, e)@@ -281,7 +425,7 @@ {-# INLINE mergeBed #-}  mergeBedWith :: (BEDLike b, Monad m)-              => ([b] -> b) -> [b] -> Source m b+             => ([b] -> a) -> [b] -> Source m a mergeBedWith f = mergeSortedBedWith f . sortBed {-# INLINE mergeBedWith #-} @@ -295,7 +439,7 @@ {-# INLINE mergeSortedBed #-}  mergeSortedBedWith :: (BEDLike b, Monad m)-                   => ([b] -> b) -> Sorted (V.Vector b) -> Source m b+                   => ([b] -> a) -> Sorted (V.Vector b) -> Source m a mergeSortedBedWith mergeFn (Sorted beds) = do     (_, r) <- V.foldM' f acc0 . V.tail $ beds     yield $ mergeFn r@@ -313,124 +457,147 @@         e' = chromEnd bed {-# INLINE mergeSortedBedWith #-} --- * BED6 format+-- | Split overlapped regions into non-overlapped regions. The input must be overlapped.+-- This function is usually used with `mergeBedWith`.+splitOverlapped :: BEDLike b => ([b] -> a) -> [b] -> [(BED3, a)]+splitOverlapped fun xs = filter ((>0) . bedSize . fst) $+    evalState (F.foldrM f [] $ init xs') x0+  where+    x0 = (\(a,b) -> (fromEither a, M.singleton (chromStart b, chromEnd b) b)) $ last xs'+    xs' = sortBy (comparing (fromEither . fst)) $ concatMap+        ( \x -> [(Left $ chromStart x, x), (Right $ chromEnd x, x)] ) xs+    f (i, x) acc = do+        (j, set) <- get+        let bed = (BED3 chr (fromEither i) j, fun $ M.elems set)+            set' = case i of+                Left _ -> M.delete (chromStart x, chromEnd x) set+                Right _ -> M.insert (chromStart x, chromEnd x) x set+        put (fromEither i, set')+        return (bed:acc)+    fromEither (Left x) = x+    fromEither (Right x) = x+    chr = chrom $ head xs+{-# INLINE splitOverlapped #-} --- | BED6 format, as described in http://genome.ucsc.edu/FAQ/FAQformat.html#format1.7-data BED = BED-    { _chrom :: !B.ByteString-    , _chromStart :: {-# UNPACK #-} !Int-    , _chromEnd :: {-# UNPACK #-} !Int-    , _name :: !(Maybe B.ByteString)-    , _score :: !(Maybe Double)-    , _strand :: !(Maybe Bool)  -- ^ True: "+", False: "-"-    } deriving (Eq, Show, Read)+-- | Read records from a bed file handler in a streaming fashion.+hReadBed :: (BEDLike b, MonadIO m) => Handle -> Source m b+hReadBed h = do eof <- liftIO $ hIsEOF h+                unless eof $ do+                    line <- liftIO $ B.hGetLine h+                    yield $ fromLine line+                    hReadBed h+{-# INLINE hReadBed #-} -instance Default BED where-    def = BED-        { _chrom = ""-        , _chromStart = 0-        , _chromEnd = 0-        , _name = Nothing-        , _score = Nothing-        , _strand = Nothing-        }+-- | Non-streaming version.+hReadBed' :: (BEDLike b, MonadIO m) => Handle -> m [b]+hReadBed' h = hReadBed h $$ sinkList+{-# INLINE hReadBed' #-} -instance BEDLike BED where-    asBed chr s e = BED chr s e Nothing Nothing Nothing+-- | Read records from a bed file in a streaming fashion.+readBed :: (BEDLike b, MonadIO m) => FilePath -> Source m b+readBed fl = do handle <- liftIO $ openFile fl ReadMode+                hReadBed handle+                liftIO $ hClose handle+{-# INLINE readBed #-} -    fromLine l = evalState (f (B.split '\t' l)) 1-      where-        f :: [B.ByteString] -> State Int BED-        f [] = do i <- get-                  if i <= 3 then error "Read BED fail: Incorrect number of fields"-                            else return def-        f (x:xs) = do-            i <- get-            put (i+1)-            bed <- f xs-            case i of-                1 -> return $ bed {_chrom = x}-                2 -> return $ bed {_chromStart = readInt x}-                3 -> return $ bed {_chromEnd = readInt x}-                4 -> return $ bed {_name = guard' x}-                5 -> return $ bed {_score = getScore x}-                6 -> return $ bed {_strand = getStrand x}-                _ -> return def+-- | Non-streaming version.+readBed' :: (BEDLike b, MonadIO m) => FilePath -> m [b]+readBed' fl = readBed fl $$ sinkList+{-# INLINE readBed' #-} -        guard' x | x == "." = Nothing-                 | otherwise = Just x-        getScore x | x == "." = Nothing-                   | otherwise = Just . readDouble $ x-        getStrand str | str == "-" = Just False-                      | str == "+" = Just True-                      | otherwise = Nothing-    {-# INLINE fromLine #-}+hWriteBed :: (BEDLike b, MonadIO m) => Handle -> Sink b m ()+hWriteBed handle = do+    x <- await+    case x of+        Nothing -> return ()+        Just bed -> (liftIO . B.hPutStrLn handle . toLine) bed >> hWriteBed handle+{-# INLINE hWriteBed #-} -    toLine (BED f1 f2 f3 f4 f5 f6) = B.intercalate "\t" [ f1-                                                        , (B.pack.show) f2-                                                        , (B.pack.show) f3-                                                        , fromMaybe "." f4-                                                        , score'-                                                        , strand'-                                                        ]-      where-        strand' | f6 == Just True = "+"-                | f6 == Just False = "-"-                | otherwise = "."-        score' = case f5 of-                     Just x -> (B.pack.show) x-                     _ -> "."-    {-# INLINE toLine #-}+hWriteBed' :: (BEDLike b, MonadIO m) => Handle -> [b] -> m ()+hWriteBed' handle beds = yieldMany beds $$ hWriteBed handle+{-# INLINE hWriteBed' #-} -    chrom = _chrom-    chromStart = _chromStart-    chromEnd = _chromEnd+writeBed :: (BEDLike b, MonadIO m) => FilePath -> Sink b m ()+writeBed fl = do handle <- liftIO $ openFile fl WriteMode+                 hWriteBed handle+                 liftIO $ hClose handle+{-# INLINE writeBed #-} +writeBed' :: (BEDLike b, MonadIO m) => FilePath -> [b] -> m ()+writeBed' fl beds = yieldMany beds $$ writeBed fl+{-# INLINE writeBed' #-}+ -- | retreive sequences fetchSeq :: (BioSeq DNA a, MonadIO m) => Genome -> Conduit BED m (Either String (DNA a))-fetchSeq g = do gH <- liftIO $ gHOpen g-                table <- liftIO $ readIndex gH-                conduitWith gH table-                liftIO $ gHClose gH+fetchSeq g = mapMC f   where-    conduitWith h index' = do-        bed <- await-        case bed of-            Just (BED chr start end _ _ isForward) -> do-                dna <- liftIO $ getSeq h index' (chr, start, end)-                case isForward of-                    Just False -> yield $ fmap rc dna-                    _ -> yield dna-                conduitWith h index'-            _ -> return ()+    f (BED chr start end _ _ isForward) = do+        dna <- liftIO $ getSeq g (chr, start, end)+        return $ case isForward of+            Just False -> rc <$> dna+            _ -> dna {-# INLINE fetchSeq #-}  fetchSeq' :: (BioSeq DNA a, MonadIO m) => Genome -> [BED] -> m [Either String (DNA a)] fetchSeq' g beds = yieldMany beds $= fetchSeq g $$ sinkList {-# INLINE fetchSeq' #-} --- * BED3 format--data BED3 = BED3 !B.ByteString !Int !Int deriving (Eq, Show, Read)--instance Default BED3 where-    def = BED3 "" 0 0--instance BEDLike BED3 where-    asBed = BED3--    fromLine l = case B.split '\t' l of-                    (a:b:c:_) -> BED3 a (readInt b) $ readInt c-                    _ -> error "Read BED fail: Incorrect number of fields"-    {-# INLINE fromLine #-}--    toLine (BED3 f1 f2 f3) = B.intercalate "\t" [f1, (B.pack.show) f2, (B.pack.show) f3]-    {-# INLINE toLine #-}+-- | Identify motif binding sites+motifScan :: (BEDLike b, MonadIO m)+          => Genome -> [Motif] -> Bkgd -> Double -> Conduit b m BED+motifScan g motifs bg p = awaitForever $ \bed -> do+    let chr = chrom bed+        s = chromStart bed+        e = chromEnd bed+    r <- liftIO $ getSeq g (chr, s, e)+    case r of+        Left _ -> return ()+        Right dna -> mapM_ (getTFBS dna (chr, s)) motifs'+  where+    getTFBS dna (chr, s) (nm, (pwm, cutoff), (pwm', cutoff')) = toProducer+        ( (findTFBS bg pwm (dna :: DNA IUPAC) cutoff True =$=+            mapC (\i -> BED chr (s+i) (s+i+n) (Just nm) Nothing $ Just True)) >>+          (findTFBS bg pwm' dna cutoff' True =$=+            mapC (\i -> BED chr (s+i) (s+i+n) (Just nm) Nothing $ Just False)) )+      where+        n = Bio.Motif.size pwm+    motifs' = flip map motifs $ \(Motif nm pwm) ->+        let cutoff = pValueToScore p bg pwm+            cutoff' = pValueToScore p bg pwm'+            pwm' = rcPWM pwm+        in (nm, (pwm, cutoff), (pwm', cutoff'))+{-# INLINE motifScan #-} -    chrom (BED3 f1 _ _) = f1-    chromStart (BED3 _ f2 _) = f2-    chromEnd (BED3 _ _ f3) = f3+-- | Retrieve motif matching scores+getMotifScore :: MonadIO m+              => Genome -> [Motif] -> Bkgd -> Conduit BED m BED+getMotifScore g motifs bg = awaitForever $ \(BED chr s e (Just nm) _ isForward) -> do+    r <- liftIO $ getSeq g (chr, s, e)+    let r' = case isForward of+            Just False -> rc <$> r+            _ -> r+    case r' of+        Left _ -> return ()+        Right dna -> do+            let pwm = M.lookupDefault (error "can't find motif with given name")+                      nm motifMap+                sc = score bg pwm (dna :: DNA IUPAC)+            yield $ BED chr s e (Just nm) (Just sc) isForward+  where+    motifMap = M.fromListWith (error "found motif with same name") $+        map (\(Motif nm pwm) -> (nm, pwm)) motifs+{-# INLINE getMotifScore #-} -convert :: (BEDLike b1, BEDLike b2) => b1 -> b2-convert b = asBed (chrom b) (chromStart b) (chromEnd b)-{-# INLINE convert #-}+getMotifPValue :: MonadIO m+               => [Motif] -> Bkgd -> Conduit BED m BED+getMotifPValue motifs bg = awaitForever $ \bed -> do+    let sc = fromJust $ _score bed+        nm = fromJust $ _name bed+        d = M.lookupDefault (error "can't find motif with given name")+                nm motifMap+        p = 1 - cdf d sc+    yield bed{_score = Just p}+  where+    motifMap = M.fromListWith (error "found motif with same name") $+        map (\(Motif nm pwm) -> (nm, scoreCDF bg pwm)) motifs+{-# INLINE getMotifPValue #-}
src/Bio/GO/GREAT.hs view
@@ -5,22 +5,22 @@     , enrichedTerms     ) where -import Control.Monad.Primitive-import qualified Data.ByteString.Char8 as B-import Conduit-import Data.Default.Class-import qualified Data.HashMap.Strict as M-import qualified Data.IntervalMap as IM-import qualified Data.Vector as V-import Data.Vector.Algorithms.Intro (sortBy)-import Data.List (sort, group)-import Data.Ord (comparing)-import Data.IntervalMap.Interval (lowerBound, upperBound)-import Data.Maybe+import           Conduit+import           Control.Monad.Primitive+import qualified Data.ByteString.Char8        as B+import           Data.Default.Class+import           Data.Function                (on)+import qualified Data.HashMap.Strict          as M+import qualified Data.IntervalMap             as IM+import           Data.List                    (foldl', group, sort, sortBy)+import           Data.Maybe+import           Data.Ord                     (comparing)+import qualified Data.Vector                  as V+import qualified Data.Vector.Algorithms.Intro as I -import Bio.Data.Bed-import Bio.GO-import Bio.Utils.Functions+import           Bio.Data.Bed+import           Bio.GO+import           Bio.Utils.Functions  -- | how to associate genomic regions with genes data AssocRule = BasalPlusExtension Int Int Int@@ -28,33 +28,31 @@                | OneNearest  instance Default AssocRule where-    def = BasalPlusExtension 5000 1000 50000+    def = BasalPlusExtension 5000 1000 1000000 -type Gene = ( B.ByteString  -- ^ chromosome-            , Int           -- ^ tss-            , Bool          -- ^ is forward stranded-            , [GOId]-            )+type Gene a = ( ( B.ByteString  -- ^ chromosome+              , Int           -- ^ tss+              , Bool          -- ^ is forward stranded+              ), a)  -- | given a gene list and the rule, compute the rulatory domain for each gene-getRegulatoryDomains :: AssocRule -> [Gene] -> BEDTree [GOId]-getRegulatoryDomains (BasalPlusExtension up dw ext) gs =-    bedToTree undefined $ flip map basal $ \(BED3 chr s e, go) ->-        let intervals = fromJust . M.lookup chr $ basalTree-            s' = let im = IM.intersecting intervals $ IM.OpenInterval (s-ext) s-                 in if IM.null im-                     then s - ext-                     else min s $ maximum $ map upperBound $ IM.keys im-            e' = let im = IM.intersecting intervals $ IM.OpenInterval e (e+ext)-                 in if IM.null im-                     then e + ext-                     else max e $ minimum $ map lowerBound $ IM.keys im-        in (BED3 chr s' e', go)+getRegulatoryDomains :: AssocRule -> [Gene a] -> [(BED3, a)]+getRegulatoryDomains (BasalPlusExtension up dw ext) genes = (extendTail r ext, a) : rs   where-    basalTree = bedToTree (error "encounter redundant regions") basal-    basal = flip map gs $ \(chr,tss,str,go) ->-        if str then (BED3 chr (tss - up) (tss + dw), go)-               else (BED3 chr (tss - dw) (tss + up), go)+    (rs, Just (r,a)) = foldl' f ([], Nothing) $ sortBy (compareBed `on` fst) basal+    f (acc, Nothing) (b,x) = (acc, Just (extendHead b ext, x))+    f (acc, Just (b', x')) (b,x)+        | chrom b' /= chrom b = ( (extendTail b' ext, x') : acc+                                , Just (extendHead b ext, x) )+        | chromEnd b' >= chromStart b = ((b',x') : acc, Just (b,x))+        | otherwise = let ext' = min ext $ (chromStart b - chromEnd b') `div` 2+                      in ((extendTail b' ext', x') : acc, Just (extendHead b ext', x))+    extendHead (BED3 chr s e) l | s - l >= 0 = BED3 chr (s-l) e+                                | otherwise = BED3 chr 0 e+    extendTail (BED3 chr s e) l = BED3 chr s (e+l)+    basal = flip map genes $ \((chr, tss, str), x) ->+        if str then (BED3 chr (tss - up) (tss + dw), x)+               else (BED3 chr (tss - dw) (tss + up), x) getRegulatoryDomains _ _ = undefined {-# INLINE getRegulatoryDomains #-} @@ -130,5 +128,5 @@     v <- V.unsafeThaw $ V.fromList $ M.toList $ flip M.mapWithKey table1 $ \t c ->         let k = fromMaybe (error "x") $ M.lookup t table2         in (1 - hyperquick c k n1 n2, fromIntegral (c*n2) / fromIntegral (n1*k))-    sortBy (comparing snd) v+    I.sortBy (comparing snd) v     V.unsafeFreeze v
src/Bio/Motif.hs view
@@ -267,7 +267,7 @@     loop (prev,scFn) i         | i < n =             let (lo,hi) = minMax (1/0,-1/0) 0-                nBin' = min 100000 $ ceiling $ (hi - lo) / precision+                nBin' = min 200000 $ ceiling $ (hi - lo) / precision                 step = (hi - lo) / fromIntegral nBin'                 idx x = let j = truncate $ (x - lo) / step                         in if j >= nBin' then nBin' - 1 else j@@ -299,7 +299,7 @@                  in minMax (foldr min l [s1,s2,s3,s4],foldr max h [s1,s2,s3,s4]) (x+1)             | otherwise = minMax (l,h) (x+1)     toCDF (v, scFn) = CDF $ V.imap (\i x -> (scFn i, x)) $ V.scanl1 (+) v-    precision = 0.002+    precision = 1e-4     n = size pwm     log' x | x == 0 = log 0.001            | otherwise = log x@@ -350,7 +350,7 @@               1 -> do when (startOfPwm x) $ put (2, str ++ [B.words x !! 7])                       go xs               2 -> let x' = B.dropWhile (== ' ') x-                   in if B.null x'+                   in if B.null x' || "URL" `B.isPrefixOf` x'                          then do put (0, [])                                  r <- go xs                                  return (toMotif str : r)
src/Bio/Seq/IO.hs view
@@ -2,13 +2,10 @@  module Bio.Seq.IO     ( Genome-    , GenomeH-    , gHOpen-    , gHClose-    , pack-    , getSeqs+    , openGenome+    , closeGenome+    , withGenome     , getSeq-    , readIndex     , getChrom     , getChrSizes     , mkIndex@@ -16,6 +13,7 @@  import Bio.Seq import Bio.Utils.Misc (readInt)+import Control.Exception (bracket) import qualified Data.ByteString.Char8 as B import qualified Data.HashMap.Lazy as M import Data.List.Split@@ -24,9 +22,9 @@ -- | The first 2048 bytes are header. Header consists of a magic string, -- followed by chromosome information. Example: -- <HASKELLBIOINFORMATICS>\nCHR1 START SIZE-newtype Genome = G FilePath+data Genome = Genome !Handle !IndexTable -newtype GenomeH = GH Handle+type IndexTable = M.HashMap B.ByteString (Int, Int)  headerSize :: Int headerSize = 2048@@ -34,84 +32,70 @@ magic :: String magic = "<HASKELLBIOINFORMATICS>" -pack :: FilePath -> IO Genome-pack fl = withFile fl ReadMode f-  where f h = do l <- hGetLine h-                 if l == magic-                    then return $ G fl-                    else error "Bio.Seq.Query.pack: Incorrect format"--gHOpen :: Genome -> IO GenomeH-gHOpen (G fl) = do h <- openFile fl ReadMode-                   return $ GH h+openGenome :: FilePath -> IO Genome+openGenome fl = do+    h <- openFile fl ReadMode+    sig <- hGetLine h+    if sig == magic+        then Genome h <$> readIndex h+        else error "Bio.Seq.Query.openGenome: Incorrect format"+  where+    readIndex h = do+        header <- B.hGetLine h+        return $ M.fromList . map f . chunksOf 3 . B.words $ header+      where+        f [k, v, l] = (k, (readInt v, readInt l))+        f _ = error "error"+{-# INLINE openGenome #-} -gHClose :: GenomeH -> IO ()-gHClose (GH h) = hClose h+closeGenome :: Genome -> IO ()+closeGenome (Genome h _) = hClose h+{-# INLINE closeGenome #-} -type IndexTable = M.HashMap B.ByteString (Int, Int)+withGenome :: FilePath -> (Genome -> IO a) -> IO a+withGenome fl fn = bracket (openGenome fl) closeGenome fn+{-# INLINE withGenome #-}  type Query = (B.ByteString, Int, Int) -- (chr, start, end), zero-based index, half-close-half-open -getSeqs :: BioSeq s a => Genome -> [Query] -> IO [Either String (s a)]-getSeqs g querys = do gH <- gHOpen g-                      index <- readIndex gH-                      r <- mapM (getSeq gH index) querys-                      gHClose gH-                      return r-{-# INLINE getSeqs #-}--getSeq :: BioSeq s a => GenomeH -> IndexTable -> Query -> IO (Either String (s a))-getSeq (GH h) index (chr, start, end) = case M.lookup chr index of-    Just (chrStart, chrSize) -> if end > chrSize-        then return $ Left $ "Bio.Seq.getSeq: out of index: " ++-                 show end ++ ">" ++ show chrSize-        else do-            hSeek h AbsoluteSeek $ fromIntegral $ headerSize + chrStart + start-            (Right . fromBS) <$> B.hGet h (end - start)+getSeq :: BioSeq s a => Genome -> Query -> IO (Either String (s a))+getSeq (Genome h index) (chr, start, end) = case M.lookup chr index of+    Just (chrStart, chrSize) ->+        if end > chrSize+            then return $ Left $ "Bio.Seq.getSeq: out of index: " ++ show end +++                ">" ++ show chrSize+            else do+                hSeek h AbsoluteSeek $ fromIntegral $ headerSize + chrStart + start+                (Right . fromBS) <$> B.hGet h (end - start)     _ -> return $ Left $ "Bio.Seq.getSeq: Cannot find " ++ show chr {-# INLINE getSeq #-} -getChrom :: Genome -> B.ByteString -> IO (Maybe (DNA IUPAC))-getChrom g chr = do-    chrSize <- getChrSizes g-    case lookup chr chrSize of-        Just s -> do [Right dna] <- getSeqs g [(chr, 0, s)]-                     return $ Just dna-        _ -> return Nothing+getChrom :: Genome -> B.ByteString -> IO (Either String (DNA IUPAC))+getChrom g chr = case lookup chr chrSize of+    Just s -> getSeq g (chr, 0, s)+    _ -> return $ Left "Unknown chromosome"+  where+    chrSize = getChrSizes g {-# INLINE getChrom #-} -getChrSizes :: Genome -> IO [(B.ByteString, Int)]-getChrSizes g = do gh <- gHOpen g-                   table <- readIndex gh-                   gHClose gh-                   return . map (\(k, (_, l)) -> (k, l)) . M.toList $ table+getChrSizes :: Genome -> [(B.ByteString, Int)]+getChrSizes (Genome _ table) = map (\(k, (_, l)) -> (k, l)) . M.toList $ table {-# INLINE getChrSizes #-} -readIndex :: GenomeH -> IO IndexTable-readIndex (GH h) = do header <- B.hGetLine h >> B.hGetLine h-                      return $ M.fromList . map f . chunksOf 3 . B.words $ header-  where-    f [k, v, l] = (k, (readInt v, readInt l))-    f _ = error "error"-{-# INLINE readIndex #-}- -- | indexing a genome. mkIndex :: [FilePath]    -- ^ fasta files representing individual chromosomes         -> FilePath      -- ^ output file         -> IO ()-mkIndex fls outFl = do-    outH <- openFile outFl WriteMode+mkIndex fls outFl = withFile outFl WriteMode $ \outH -> do     hPutStr outH $ magic ++ "\n" ++ replicate 2024 '#'  -- header size: 1024     chrs <- mapM (write outH) fls     hSeek outH AbsoluteSeek 24     B.hPutStrLn outH $ mkHeader chrs-    hClose outH   where-    write handle fl = do h <- openFile fl ReadMode-                         fastaHeader <- B.hGetLine h-                         n <- loop 0 h-                         hClose h-                         return (B.tail fastaHeader, n)+    write handle fl = withFile fl ReadMode $ \h -> do+        fastaHeader <- B.hGetLine h+        n <- loop 0 h+        return (B.tail fastaHeader, n)       where         loop !n h' = do eof <- hIsEOF h'                         if eof then return n
src/Bio/Utils/Functions.hs view
@@ -33,7 +33,7 @@ import qualified Data.Vector as V import qualified Data.Vector.Generic as G import qualified Data.Vector.Unboxed as U-import Statistics.Sample (meanVarianceUnb, mean)+import Statistics.Sample (meanVarianceUnb) import Statistics.Function (sortBy)  -- | inverse hyperbolic sine transformation
tests/Tests/Bed.hs view
@@ -3,6 +3,8 @@ module Tests.Bed (tests) where  import Bio.Data.Bed+import Data.List (sortBy)+import Data.Function (on) import Test.Tasty import Test.Tasty.HUnit import qualified Data.Vector as V@@ -11,6 +13,7 @@ tests = testGroup "Test: Bio.Data.Bed"     [ testCase "sortBed" sortBedTest     , testCase "split" splitBedTest+    , testCase "splitOverlapped" splitOverlappedTest     ]  sortBedTest :: Assertion@@ -31,3 +34,28 @@     s3' = map f [ ( 0, 20), (10, 30), (20, 40), (30, 50), (40, 60)                 , (50, 70), (60, 80), (70, 90), (80, 100) ]     f (a,b) = asBed "chr1" a b++splitOverlappedTest :: Assertion+splitOverlappedTest = expect @=? result+  where+    input = map (\(a,b) -> BED3 "chr1" a b)+        [ (0, 100)+        , (10, 20)+        , (50, 150)+        , (120, 160)+        , (155, 200)+        , (155, 220)+        ]+    expect = map (\((a,b), x) -> (BED3 "chr1" a b, x))+        [ ((0, 10), 1)+        , ((10, 20), 2)+        , ((20, 50), 1)+        , ((50, 100), 2)+        , ((100, 120), 1)+        , ((120, 150), 2)+        , ((150, 155), 1)+        , ((155, 160), 3)+        , ((160, 200), 2)+        , ((200, 220), 1)+        ]+    result = sortBy (compareBed `on` fst) $ splitOverlapped length input
+ tests/Tests/GREAT.hs view
@@ -0,0 +1,34 @@+{-# LANGUAGE OverloadedStrings #-}+module Tests.GREAT (tests) where++import           Bio.Data.Bed+import           Bio.GO.GREAT+import           Data.List+import           Data.Ord++import           Test.Tasty+import           Test.Tasty.HUnit++testData = (genes, results)+  where+    genes = [ (("chr1", 1000, True), 0)+            , (("chr1", 6000, True), 1)+            , (("chr2", 6000, True), 2)+            , (("chr2", 2000000, False), 3)+            ]+    results = [ (BED3 "chr1" 0 2000, 0)+              , (BED3 "chr1" 1000 1007000, 1)+              , (BED3 "chr2" 0 1003000, 2)+              , (BED3 "chr2" 1003000 3005000, 3)+              ]++tests :: TestTree+tests = testGroup "Test: Bio.GO.GREAT"+    [ testCase "getRegulatoryDomains" testAssoc+    ]++testAssoc :: Assertion+testAssoc = do+    let result = sortBy (comparing snd) $ getRegulatoryDomains+             (BasalPlusExtension 5000 1000 1000000) $ fst testData+    snd testData @=? result
tests/Tests/Motif.hs view
@@ -50,7 +50,7 @@ findTFBSTest = do     ms <- motifs     let (Motif _ pwm) = head ms-    expect <- findTFBS def pwm dna (0.6 * optimalScore def pwm) $$ CL.consume+    expect <- findTFBS def pwm dna (0.6 * optimalScore def pwm) True $$ CL.consume     actual <- findTFBSSlow def pwm dna (0.6 * optimalScore def pwm) $$ CL.consume     assertEqual "findTFBS" expect actual
tests/test.hs view
@@ -3,6 +3,7 @@ import qualified Tests.Motif as Motif import qualified Tests.ChIPSeq as ChIPSeq import qualified Tests.Seq as Seq+import qualified Tests.GREAT as GREAT import Test.Tasty  main :: IO ()@@ -12,4 +13,5 @@     , Seq.tests     , ChIPSeq.tests     , Motif.tests+    , GREAT.tests     ]