diff --git a/bioinformatics-toolkit.cabal b/bioinformatics-toolkit.cabal --- a/bioinformatics-toolkit.cabal +++ b/bioinformatics-toolkit.cabal @@ -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 diff --git a/src/Bio/Data/Bed/Types.hs b/src/Bio/Data/Bed/Types.hs --- a/src/Bio/Data/Bed/Types.hs +++ b/src/Bio/Data/Bed/Types.hs @@ -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 -> "+" diff --git a/src/Bio/Data/Bed/Utils.hs b/src/Bio/Data/Bed/Utils.hs --- a/src/Bio/Data/Bed/Utils.hs +++ b/src/Bio/Data/Bed/Utils.hs @@ -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 () diff --git a/src/Bio/Motif.hs b/src/Bio/Motif.hs --- a/src/Bio/Motif.hs +++ b/src/Bio/Motif.hs @@ -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 diff --git a/src/Bio/Motif/Search.hs b/src/Bio/Motif/Search.hs --- a/src/Bio/Motif/Search.hs +++ b/src/Bio/Motif/Search.hs @@ -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 #-} + -} diff --git a/tests/Tests/Motif.hs b/tests/Tests/Motif.hs --- a/tests/Tests/Motif.hs +++ b/tests/Tests/Motif.hs @@ -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