diff --git a/app/viewSeq.hs b/app/viewSeq.hs --- a/app/viewSeq.hs +++ b/app/viewSeq.hs @@ -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 diff --git a/bioinformatics-toolkit.cabal b/bioinformatics-toolkit.cabal --- a/bioinformatics-toolkit.cabal +++ b/bioinformatics-toolkit.cabal @@ -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: diff --git a/data/1.fasta b/data/1.fasta deleted file mode 100644 --- a/data/1.fasta +++ /dev/null @@ -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 diff --git a/src/Bio/ChIPSeq.hs b/src/Bio/ChIPSeq.hs --- a/src/Bio/ChIPSeq.hs +++ b/src/Bio/ChIPSeq.hs @@ -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 diff --git a/src/Bio/Data/Bed.hs b/src/Bio/Data/Bed.hs --- a/src/Bio/Data/Bed.hs +++ b/src/Bio/Data/Bed.hs @@ -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 #-} diff --git a/src/Bio/GO/GREAT.hs b/src/Bio/GO/GREAT.hs --- a/src/Bio/GO/GREAT.hs +++ b/src/Bio/GO/GREAT.hs @@ -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 diff --git a/src/Bio/Motif.hs b/src/Bio/Motif.hs --- a/src/Bio/Motif.hs +++ b/src/Bio/Motif.hs @@ -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) diff --git a/src/Bio/Seq/IO.hs b/src/Bio/Seq/IO.hs --- a/src/Bio/Seq/IO.hs +++ b/src/Bio/Seq/IO.hs @@ -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: -- \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 = "" -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 diff --git a/src/Bio/Utils/Functions.hs b/src/Bio/Utils/Functions.hs --- a/src/Bio/Utils/Functions.hs +++ b/src/Bio/Utils/Functions.hs @@ -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 diff --git a/tests/Tests/Bed.hs b/tests/Tests/Bed.hs --- a/tests/Tests/Bed.hs +++ b/tests/Tests/Bed.hs @@ -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 diff --git a/tests/Tests/GREAT.hs b/tests/Tests/GREAT.hs new file mode 100644 --- /dev/null +++ b/tests/Tests/GREAT.hs @@ -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 diff --git a/tests/Tests/Motif.hs b/tests/Tests/Motif.hs --- a/tests/Tests/Motif.hs +++ b/tests/Tests/Motif.hs @@ -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 diff --git a/tests/test.hs b/tests/test.hs --- a/tests/test.hs +++ b/tests/test.hs @@ -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 ]