bioinformatics-toolkit 0.7.0 → 0.8.0

raw patch · 6 files changed

+115/−102 lines, 6 filesPVP ok

version bump matches the API change (PVP)

API changes (from Hackage documentation)

Files

• bioinformatics-toolkit.cabal +1/−1
• src/Bio/Data/Bed/Types.hs +15/−12
• src/Bio/Data/Bed/Utils.hs +50/−65
• src/Bio/Motif.hs +12/−9
• src/Bio/Motif/Search.hs +27/−15
• tests/Tests/Motif.hs +10/−0

@@ -1,5 +1,5 @@ name:                bioinformatics-toolkit-version:             0.7.0+version:             0.8.0 synopsis:            A collection of bioinformatics tools description:         A collection of bioinformatics tools license:             MIT

@@ -56,7 +56,7 @@     chromStart :: Lens' b Int     chromEnd :: Lens' b Int     name :: Lens' b (Maybe B.ByteString)-    score :: Lens' b (Maybe Double)+    score :: Lens' b (Maybe Int)     strand :: Lens' b (Maybe Bool)      -- | Return the size of a bed region.@@ -90,7 +90,7 @@     , _bed_chromStart :: !Int     , _bed_chromEnd   :: !Int     , _bed_name       :: !(Maybe B.ByteString)-    , _bed_score      :: !(Maybe Double)+    , _bed_score      :: !(Maybe Int)     , _bed_strand     :: !(Maybe Bool)  -- ^ True: "+", False: "-"     } deriving (Eq, Show, Read) @@ -122,21 +122,24 @@         getName x | x == "." = Nothing                   | otherwise = Just x         getScore x | x == "." = Nothing-                   | otherwise = Just . readDouble $ x+                   | s >= 0 && s <= 1000 = Just s+                   | otherwise = error "score must be in [0, 1000]."+          where+            s = readInt x         getStrand str | str == "-" = Just False                       | str == "+" = Just True                       | otherwise = Nothing     {-# INLINE fromLine #-}      toLine (BED f1 f2 f3 f4 f5 f6) = B.intercalate "\t"-        [ f1, (B.pack.show) f2, (B.pack.show) f3, fromMaybe "." f4, score'-        , strand' ]+        [ f1, fromJust $ packDecimal f2, fromJust $ packDecimal f3+        , fromMaybe "." f4, score', strand' ]       where         strand' | f6 == Just True = "+"                 | f6 == Just False = "-"                 | otherwise = "."         score' = case f5 of-                     Just x -> toShortest x+                     Just x -> fromJust $ packDecimal x                      _      -> "."     {-# INLINE toLine #-} @@ -180,7 +183,7 @@     , _npStart  :: !Int     , _npEnd    :: !Int     , _npName   :: !(Maybe B.ByteString)-    , _npScore  :: !Double+    , _npScore  :: !Int     , _npStrand :: !(Maybe Bool)     , _npSignal  :: !Double     , _npPvalue :: !(Maybe Double)@@ -207,7 +210,7 @@      fromLine l = NarrowPeak a (readInt b) (readInt c)         (if d == "." then Nothing else Just d)-        (readDouble e)+        (readInt e)         (if f == "." then Nothing else if f == "+" then Just True else Just False)         (readDouble g)         (readDoubleNonnegative h)@@ -219,7 +222,7 @@      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+        , fromJust $ packDecimal e         , case f of             Nothing   -> "."             Just True -> "+"@@ -239,7 +242,7 @@     , _bpStart  :: !Int     , _bpEnd    :: !Int     , _bpName   :: !(Maybe B.ByteString)-    , _bpScore  :: !Double+    , _bpScore  :: !Int     , _bpStrand :: !(Maybe Bool)     , _bpSignal  :: !Double     , _bpPvalue :: !(Maybe Double)@@ -264,7 +267,7 @@      fromLine l = BroadPeak a (readInt b) (readInt c)         (if d == "." then Nothing else Just d)-        (readDouble e)+        (readInt e)         (if f == "." then Nothing else if f == "+" then Just True else Just False)         (readDouble g)         (readDoubleNonnegative h)@@ -275,7 +278,7 @@      toLine (BroadPeak a b c d e f g h i) = B.intercalate "\t"         [ a, fromJust $ packDecimal b, fromJust $ packDecimal c, fromMaybe "." d-        , toShortest e+        , fromJust $ packDecimal e         , case f of             Nothing   -> "."             Just True -> "+"

@@ -1,14 +1,15 @@ {-# LANGUAGE FlexibleContexts      #-} {-# LANGUAGE OverloadedStrings     #-} {-# LANGUAGE LambdaCase #-}+{-# LANGUAGE RecordWildCards #-} {-# LANGUAGE BangPatterns #-}  module Bio.Data.Bed.Utils     ( fetchSeq     , fetchSeq'-    , motifScan-    , getMotifScore-    , getMotifPValue+    , CutoffMotif(..)+    , mkCutoffMotif+    , scanMotif      , monoColonalize     , BaseMap(..)     , baseMap@@ -39,7 +40,7 @@  import           Bio.Data.Bed import           Bio.Data.Bed.Types-import           Bio.Motif                    (Bkgd (..), Motif (..))+import           Bio.Motif                    (Bkgd (..), Motif (..), CDF, PWM) import qualified Bio.Motif                    as Motif import qualified Bio.Motif.Search             as Motif import           Bio.Seq hiding (length)@@ -64,71 +65,55 @@ fetchSeq' g beds = runConduit $ yieldMany beds .| fetchSeq g .| sinkList {-# INLINE fetchSeq' #-} --- | Identify motif binding sites-motifScan :: (BEDLike b, MonadIO m)-          => Genome -> [Motif] -> Bkgd -> Double -> ConduitT b BED m ()-motifScan g motifs bg p = awaitForever $ \bed -> do-    r <- liftIO $ getSeq g (bed^.chrom, bed^.chromStart, bed^.chromEnd)-    case r of-        Left _    -> liftIO $ hPutStrLn stderr $-            "Warning: no sequence for region: " ++ show -                (bed^.chrom, bed^.chromStart, bed^.chromEnd)-        Right dna -> mapM_ (getTFBS dna (bed^.chrom, bed^.chromStart)) motifs'-  where-    getTFBS dna (chr, s) (nm, (pwm, cutoff), (pwm', cutoff')) = toProducer-        ( (Motif.findTFBS bg pwm (dna :: DNA IUPAC) cutoff True .|-            mapC (\i -> bed & chromStart +~ i & chromEnd +~ i & strand .~ Just True)) >>-          (Motif.findTFBS bg pwm' dna cutoff' True .|-            mapC (\i -> bed & chromStart +~ i & chromEnd +~ i & strand .~ Just False)) )-      where-        n = Motif.size pwm-        bed = asBed chr s (s+n) & name .~ Just nm-    motifs' = flip map motifs $ \(Motif nm pwm) ->-        let cutoff = Motif.pValueToScore p bg pwm-            cutoff' = Motif.pValueToScore p bg pwm'-            pwm' = Motif.rcPWM pwm-        in (nm, (pwm, cutoff), (pwm', cutoff'))-{-# INLINE motifScan #-}+-- | Motif with predefined cutoff score. All necessary intermediate data+-- structure for motif scanning are stored.+data CutoffMotif = CutoffMotif+    { _motif_name :: B.ByteString+    , _motif_pwm :: PWM+    , _motif_sigma :: U.Vector Double+    , _motif_pwm_rc :: PWM+    , _motif_sigma_rc :: U.Vector Double+    , _background :: Bkgd+    , _cutoff :: Double+    , _cdf :: CDF } --- | Retrieve motif matching scores-getMotifScore :: MonadIO m-              => Genome -> [Motif] -> Bkgd -> ConduitT BED BED m ()-getMotifScore g motifs bg = awaitForever $ \bed -> do-    r <- liftIO $ getSeq g (bed^.chrom, bed^.chromStart, bed^.chromEnd)-    let r' = case bed^.strand 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")-                        (fromJust $ bed^.name) motifMap-                sc = Motif.score bg pwm (dna :: DNA IUPAC)-            yield $ score .~ Just sc $ bed+mkCutoffMotif :: Bkgd+              -> Double  -- ^ p-value+              -> Motif -> CutoffMotif+mkCutoffMotif bg p motif = CutoffMotif (_name motif) (_pwm motif) sigma pwm'+    sigma' bg sc $ Motif.truncateCDF (1 - p * 10) cdf   where-    motifMap = M.fromListWith (error "found motif with same name") $-        map (\(Motif nm pwm) -> (nm, pwm)) motifs-{-# INLINE getMotifScore #-}+    cdf = Motif.scoreCDF bg $ _pwm motif+    sc = Motif.cdf' cdf $ 1 - p+    sigma = Motif.optimalScoresSuffix bg $ _pwm motif+    pwm' = Motif.rcPWM $ _pwm motif+    sigma' = Motif.optimalScoresSuffix bg pwm' -getMotifPValue :: Monad m-               => Maybe Double   -- ^ whether to truncate the motif score CDF.-                                 -- Doing this will significantly reduce memory-                                 -- usage without sacrifice accuracy.-               -> [Motif] -> Bkgd -> ConduitT BED BED m ()-getMotifPValue truncation motifs bg = mapC $ \bed ->-    let nm = fromJust $ bed^.name-        sc = fromJust $ bed^.score-        d = M.lookupDefault (error "can't find motif with given name")-                nm motifMap-        p = 1 - Motif.cdf d sc-     in score .~ Just p $ bed+-- | Motif score is in [0, 1000]: ( 1 / (1 + exp (-(-logP - 5))) ) * 1000.+scanMotif :: (BEDLike b, MonadIO m)+          => Genome -> [CutoffMotif] -> ConduitT b BED m ()+scanMotif g motifs = awaitForever $ \bed -> do+    let (chr, start, end) = (bed^.chrom, bed^.chromStart, bed^.chromEnd)+    liftIO (getSeq g (chr, start, end)) >>= \case+        Left _    -> liftIO $ hPutStrLn stderr $+            "Warning: no sequence for region: " ++ show (chr, start, end)+        Right dna -> forM_ motifs $ \CutoffMotif{..} -> do+            let mkBed str (i, sc) = BED chr (start + i) (start + i + n)+                    (Just $ _motif_name) (Just $ toAffinity $ 1 - Motif.cdf _cdf sc)+                    (Just str)+                n = Motif.size _motif_pwm+            -- Scan forward strand+            Motif.findTFBSWith _motif_sigma _background _motif_pwm+                (dna :: DNA IUPAC) _cutoff True .| mapC (mkBed True)+            -- Scan reverse strand+            Motif.findTFBSWith _motif_sigma_rc _background _motif_pwm_rc+                dna _cutoff True .| mapC (mkBed False)   where-    motifMap = M.fromListWith (error "getMotifPValue: found motif with same name") $-        map (\(Motif nm pwm) -> (nm, compressCDF $ Motif.scoreCDF bg pwm)) motifs-    compressCDF = case truncation of-        Nothing -> id-        Just x  -> Motif.truncateCDF x-{-# INLINE getMotifPValue #-}+    toAffinity x' = round $ sc * 1000+      where+        sc = 1 / (1 + exp (-(x - 5)))+        x = negate $ logBase 10 $ max 1e-20 x'+{-# INLINE scanMotif #-}  -- | process a sorted BED stream, keep only mono-colonal tags monoColonalize :: Monad m => ConduitT BED BED m ()

@@ -252,15 +252,18 @@ cdf' :: CDF -> Double -> Double cdf' (CDF v) p     | p > 1 || p < 0 = error "p must be in [0,1]"-    | otherwise = let i = binarySearchBy cmp v p-                  in case () of-                      _ | i >= n -> 1/0-                        | i == 0 -> if p == snd (U.head v) then fst (U.head v) else undefined-                        | otherwise -> let (a, a') = v U.! (i-1)-                                           (b, b') = v U.! i-                                           α = (b' - p) / (b' - a')-                                           β = (p - a') / (b' - a')-                                       in α * a + β * b+    | otherwise =+        let i = binarySearchBy cmp v p+        in case () of+          _ | i >= n -> 1/0+            | i == 0 -> if p == snd (U.head v)+                then fst (U.head v)+                else error "cdf': impossible!"+            | otherwise -> let (a, a') = v U.! (i-1)+                               (b, b') = v U.! i+                               α = (b' - p) / (b' - a')+                               β = (p - a') / (b' - a')+                           in α * a + β * b   where     cmp (_,a) = compare a     n = U.length v

@@ -3,12 +3,10 @@  module Bio.Motif.Search     ( findTFBS-    , findTFBS'+    , findTFBSWith     , findTFBSSlow     , maxMatchSc-    , SpaceDistribution(..)-    , spaceConstraint-    , spaceConstraintHelper+    , optimalScoresSuffix     ) where  import           Conduit@@ -40,21 +38,35 @@          -> Double          -> Bool    -- ^ whether to skip ambiguous sequences. Recommend: True                     -- in most cases-         -> ConduitT i Int m ()-findTFBS bg pwm dna thres skip = loop 0+         -> ConduitT i (Int, Double) m ()+findTFBS bg pwm dna thres skip = findTFBSWith sigma bg pwm dna thres skip   where+    sigma = optimalScoresSuffix bg pwm+{-# INLINE findTFBS #-}++findTFBSWith :: Monad m+             => U.Vector Double   -- ^ best possible match score of suffixes+             -> Bkgd+             -> PWM+             -> DNA a+             -> Double+             -> Bool    -- ^ whether to skip ambiguous sequences. Recommend: True+                        -- in most cases+             -> ConduitT i (Int, Double) m ()+findTFBSWith sigma bg pwm dna thres skip = loop 0+  where     loop !i | i >= l - n + 1 = return ()-            | otherwise = do let (d, _) = searchFn bg pwm sigma dna i thres+            | otherwise = do let (d, sc) = searchFn bg pwm sigma dna i thres                              if d == n - 1-                                then yield i >> loop (i+1)+                                then yield (i, sc) >> loop (i+1)                                 else loop (i+1)-    sigma = optimalScoresSuffix bg pwm     l = Seq.length dna     n = size pwm     searchFn | skip = lookAheadSearch'              | otherwise = lookAheadSearch-{-# INLINE findTFBS #-}+{-# INLINE findTFBSWith #-} +{- -- | a parallel version of findTFBS findTFBS' :: Bkgd           -> PWM@@ -78,14 +90,15 @@     searchFn | skip = lookAheadSearch'              | otherwise = lookAheadSearch {-# INLINE findTFBS' #-}+-}  -- | use naive search-findTFBSSlow :: Monad m => Bkgd -> PWM -> DNA a -> Double -> ConduitT i Int m ()+findTFBSSlow :: Monad m => Bkgd -> PWM -> DNA a -> Double -> ConduitT i (Int, Double) m () findTFBSSlow bg pwm dna thres = scores' bg pwm dna .| loop 0   where     loop i = do v <- await                 case v of-                    Just v' ->  if v' >= thres then yield i >> loop (i+1)+                    Just v' ->  if v' >= thres then yield (i, v') >> loop (i+1)                                                else loop (i+1)                     _ -> return () {-# INLINE findTFBSSlow #-}@@ -193,6 +206,7 @@     n = size pwm {-# INLINE lookAheadSearch' #-} +{- data SpaceDistribution = SpaceDistribution     { _motif1   :: Motif     , _nSites1  :: (Int, Int)@@ -254,15 +268,13 @@         nRF = map (nOverlap rv1 fw2' w) $ reverse rs {-# INLINE spaceConstraintHelper #-} -{- computePValue :: Double -> [Int] -> [(Int, Double)] computePValue p xs = zip xs $ map (pValue n p) xs   where     n = foldl' (+) 0 xs-{-# INLINE computePValue #-}  pValue :: Int -> Double -> Int -> Double pValue n p x | n > 2000 = complCumulative (poisson (fromIntegral n* p)) $ fromIntegral x              | otherwise = complCumulative (binomial n p) $ fromIntegral x-{-# INLINE pValue #-}+ -}

@@ -43,6 +43,7 @@     , testCase "Max matching score" maxScTest     , testCase "pValue calculation" pValueTest     , testCase "CDF truncate test" cdfTruncateTest+    , testCase "Reverse Complentary" reverseComplentary     ]  @@ -86,3 +87,12 @@     assertEqual "CDF truncate" expect actual   where     pValueToScore' p bg pwm = cdf' (truncateCDF 0.999 $ scoreCDF bg pwm) $ 1 - p++reverseComplentary :: Assertion+reverseComplentary = do+    ms <- motifs+    let expect = map (\x -> map approx $ scores def (_pwm x) dna) ms+        actual = map (\x -> map approx $ reverse $ scores def (rcPWM $ _pwm x) $ rc dna) ms+    assertEqual "reverse complentary" expect actual+  where+    approx x = round $ 100000 * x