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.3 +version: 0.2.4 synopsis: A collection of bioinformatics tools description: A collection of bioinformatics tools license: MIT @@ -17,8 +17,8 @@ data-files: tests/data/*.bed tests/data/*.sorted.bed - tests/data/*.bam - tests/data/*.bam.bai +-- tests/data/*.bam +-- tests/data/*.bam.bai tests/data/*.meme tests/data/*.fasta @@ -29,9 +29,9 @@ exposed-modules: Bio.ChIPSeq Bio.ChIPSeq.FragLen - Bio.Data.Bam Bio.Data.Bed -- Bio.Data.BigWig +-- Bio.Data.Bam Bio.Data.Fasta Bio.GO Bio.GO.Parser @@ -78,7 +78,6 @@ , parallel >=3.2 , primitive , palette - , samtools , split , statistics >=0.13.2.1 , text >=0.11 @@ -149,8 +148,6 @@ main-is: test.hs other-modules: Tests.Bed - , Tests.Bam - , Tests.ChIPSeq , Tests.Motif , Tests.Seq , Tests.GREAT diff --git a/src/Bio/ChIPSeq.hs b/src/Bio/ChIPSeq.hs --- a/src/Bio/ChIPSeq.hs +++ b/src/Bio/ChIPSeq.hs @@ -1,19 +1,17 @@ {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE OverloadedStrings #-} {-# LANGUAGE BangPatterns #-} + module Bio.ChIPSeq ( monoColonalize , rpkmBed , rpkmSortedBed - , rpkmBam , countTagsBinBed , countTagsBinBed' , tagCountDistr , peakCluster ) where -import Bio.SamTools.Bam -import qualified Bio.SamTools.BamIndex as BI import Control.Monad (liftM, forM_, forM) import Control.Monad.Primitive (PrimMonad) import Conduit @@ -28,7 +26,6 @@ import qualified Data.Vector.Generic as G import qualified Data.Vector.Generic.Mutable as GM -import Bio.Data.Bam import Bio.Data.Bed -- | process a sorted BED stream, keep only mono-colonal tags @@ -184,6 +181,7 @@ {-# INLINE countTagsBinBed' #-} +{- -- | calculate RPKM using BAM file (*.bam) and its index file (*.bam.bai), using -- constant space rpkmBam :: BEDLike b => FilePath -> Conduit b IO Double @@ -214,6 +212,7 @@ else if l <= p1 && p1 < u then acc + 1 else acc {-# INLINE rpkmBam #-} +-} tagCountDistr :: PrimMonad m => G.Vector v Int => Sink BED m (v Int) tagCountDistr = loop M.empty diff --git a/src/Bio/Data/Bam.hs b/src/Bio/Data/Bam.hs deleted file mode 100644 --- a/src/Bio/Data/Bam.hs +++ /dev/null @@ -1,62 +0,0 @@ --------------------------------------------------------------------------------- --- | --- Module : $Header$ --- Copyright : (c) 2014 Kai Zhang --- License : MIT - --- Maintainer : kai@kzhang.org --- Stability : experimental --- Portability : portable - --- functions for processing BED files --------------------------------------------------------------------------------- - -module Bio.Data.Bam - ( readBam - , bamToBed - , viewBam - ) where - -import Bio.SamTools.Bam -import Bio.SamTools.BamIndex -import Control.Monad.Trans.Class (lift) -import qualified Data.ByteString.Char8 as B -import Conduit (Source, Conduit, yield, concatMapC) -import Data.Maybe (fromJust) -import Bio.Data.Bed - -readBam :: FilePath -> Source IO Bam1 -readBam fl = do handle <- lift $ openBamInFile fl - go handle - where - go h = do x <- lift $ get1 h - case x of - Nothing -> lift $ closeInHandle h - Just bam -> yield bam >> go h -{-# INLINE readBam #-} - -bamToBed :: Monad m => Conduit Bam1 m BED -bamToBed = concatMapC f - where - f bam =case targetName bam of - Just chr -> - let start = fromIntegral . fromJust . position $ bam - end = start + (fromIntegral . fromJust . queryLength) bam - nm = Just . queryName $ bam - strand = Just . not . isReverse $ bam - in Just $ BED chr start end nm Nothing strand - _ -> Nothing -{-# INLINE bamToBed #-} - -viewBam :: IdxHandle -> (B.ByteString, Int, Int) -> Source IO Bam1 -viewBam handle (chr, s, e) = case lookupTarget (idxHeader handle) chr of - Nothing -> return () - Just chrId -> do - q <- lift $ query handle chrId (fromIntegral s,fromIntegral e) - go q - where - go q' = do r <- lift $ next q' - case r of - Nothing -> return () - Just bam -> yield bam >> go q' -{-# INLINE viewBam #-} 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 @@ -598,8 +598,12 @@ map (\(Motif nm pwm) -> (nm, pwm)) motifs {-# INLINE getMotifScore #-} -getMotifPValue :: Monad m => [Motif] -> Bkgd -> Conduit BED m BED -getMotifPValue motifs bg = mapC $ \bed -> +getMotifPValue :: Monad m + => Maybe Double -- ^ whether to truncate the motif score CDF. + -- Doing this will significantly reduce memory + -- usage without sacrifice accuracy. + -> [Motif] -> Bkgd -> Conduit BED m BED +getMotifPValue truncation motifs bg = mapC $ \bed -> let nm = fromJust $ bedName bed sc = fromJust $ bedScore bed d = M.lookupDefault (error "can't find motif with given name") @@ -608,5 +612,8 @@ in bed{_score = Just p} where motifMap = M.fromListWith (error "getMotifPValue: found motif with same name") $ - map (\(Motif nm pwm) -> (nm, scoreCDF bg pwm)) motifs + map (\(Motif nm pwm) -> (nm, compressCDF $ scoreCDF bg pwm)) motifs + compressCDF = case truncation of + Nothing -> id + Just x -> truncateCDF x {-# INLINE getMotifPValue #-} diff --git a/src/Bio/Motif.hs b/src/Bio/Motif.hs --- a/src/Bio/Motif.hs +++ b/src/Bio/Motif.hs @@ -1,5 +1,6 @@ {-# LANGUAGE BangPatterns #-} {-# LANGUAGE OverloadedStrings #-} + module Bio.Motif ( PWM(..) , size @@ -17,6 +18,7 @@ , CDF(..) , cdf , cdf' + , truncateCDF , scoreCDF , pValueToScore , pValueToScoreExact @@ -224,6 +226,7 @@ pValueToScore :: Double -> Bkgd -> PWM -> Double pValueToScore p bg pwm = cdf' (scoreCDF bg pwm) $ 1 - p +-- | The cumulative distribution function in the form of (x, P(X <= x)) newtype CDF = CDF (U.Vector (Double, Double)) deriving (Read, Show) -- P(X <= x) @@ -259,6 +262,11 @@ cmp (_,a) = compare a n = U.length v {-# INLINE cdf' #-} + +-- | Truncate the CDF by a value, in order to reduce the memory usage. +truncateCDF :: Double -> CDF -> CDF +truncateCDF x (CDF v) = CDF $ U.filter ((>=x) . snd) v +{-# INLINE truncateCDF #-} -- approximate the cdf of motif matching scores scoreCDF :: Bkgd -> PWM -> CDF diff --git a/src/Bio/RealWorld/GENCODE.hs b/src/Bio/RealWorld/GENCODE.hs --- a/src/Bio/RealWorld/GENCODE.hs +++ b/src/Bio/RealWorld/GENCODE.hs @@ -35,7 +35,7 @@ where f l | B.head l == '#' || f3 /= "gene" = Nothing | otherwise = Just $ Gene (mk $ getField "gene_name") (getField "gene_id") f1 - (readInt f4) (readInt f5) (f7=="+") + (readInt f4 - 1) (readInt f5) (f7=="+") where [f1,_,f3,f4,f5,_,f7,_,f9] = B.split '\t' l fields = map (B.break (==' ') . strip) $ B.split ';' f9 diff --git a/tests/Tests/Bam.hs b/tests/Tests/Bam.hs deleted file mode 100644 --- a/tests/Tests/Bam.hs +++ /dev/null @@ -1,21 +0,0 @@ -{-# LANGUAGE OverloadedStrings #-} - -module Tests.Bam (tests) where - -import Bio.Data.Bed -import Bio.Data.Bam -import Conduit -import Test.Tasty -import Test.Tasty.HUnit -import qualified Data.Vector as V - -tests :: TestTree -tests = testGroup "Test: Bio.Data.Bam" - [ testCase "bamToBed" bamToBedTest - ] - -bamToBedTest :: Assertion -bamToBedTest = do - bed <- readBed' "tests/data/example.bed" - bed' <- readBam "tests/data/example.bam" =$= bamToBed $$ sinkList - bed @=? bed' diff --git a/tests/Tests/ChIPSeq.hs b/tests/Tests/ChIPSeq.hs deleted file mode 100644 --- a/tests/Tests/ChIPSeq.hs +++ /dev/null @@ -1,29 +0,0 @@ -module Tests.ChIPSeq (tests) where - -import Bio.Data.Bed -import Bio.ChIPSeq -import Data.Conduit -import qualified Data.Conduit.List as CL -import Test.Tasty -import Test.Tasty.HUnit -import qualified Data.Vector as V -import Text.Printf - -peaks :: IO [BED3] -peaks = readBed' "tests/data/peaks.bed" - -tags :: Source IO BED -tags = readBed "tests/data/example.bed" - -tests :: TestTree -tests = testGroup "Test: Bio.ChIPSeq" - [ testCase "rpkm" testRPKM - ] - -testRPKM :: Assertion -testRPKM = do regions <- peaks - r1 <- tags $$ rpkmBed regions - r2 <- CL.sourceList regions $= rpkmBam "tests/data/example.bam" $$ CL.consume - let r1' = map (printf "%0.6f") . V.toList $ r1 :: [String] - r2' = map (printf "%0.6f") r2 - r1' @=? r2' diff --git a/tests/Tests/Motif.hs b/tests/Tests/Motif.hs --- a/tests/Tests/Motif.hs +++ b/tests/Tests/Motif.hs @@ -2,19 +2,21 @@ module Tests.Motif (tests) where -import Test.Tasty -import Test.Tasty.HUnit -import System.Random -import Data.Default.Class -import qualified Data.Conduit.List as CL import qualified Data.ByteString.Char8 as B -import Data.Conduit +import Data.Conduit +import qualified Data.Conduit.List as CL +import Data.Default.Class +import System.Random +import Test.Tasty +import Test.Tasty.HUnit -import Bio.Seq -import Bio.Motif -import Bio.Motif.Search -import Bio.Data.Fasta +import Bio.Data.Fasta +import Bio.Motif +import Bio.Motif.Search +import Bio.Seq +import qualified Data.Vector.Unboxed as U + dna :: DNA Basic dna = fromBS $ B.pack $ map f $ take 5000 $ randomRs (0, 3) (mkStdGen 2) where @@ -35,6 +37,7 @@ testCase "TFBS scanning" findTFBSTest , testCase "Max matching score" maxScTest , testCase "pValue calculation" pValueTest + , testCase "CDF truncate test" cdfTruncateTest ] @@ -69,3 +72,12 @@ assertEqual "pValueToScore" expect actual where approx x = round $ 10 * x + +cdfTruncateTest :: Assertion +cdfTruncateTest = do + ms <- motifs + let expect = map (pValueToScore 1e-4 def . _pwm) ms + actual = map (pValueToScore' 1e-4 def . _pwm) ms + assertEqual "CDF truncate" expect actual + where + pValueToScore' p bg pwm = cdf' (truncateCDF 0.999 $ scoreCDF bg pwm) $ 1 - p diff --git a/tests/data/example.bam b/tests/data/example.bam deleted file mode 100644 Binary files a/tests/data/example.bam and /dev/null differ diff --git a/tests/data/example.bam.bai b/tests/data/example.bam.bai deleted file mode 100644 Binary files a/tests/data/example.bam.bai and /dev/null differ diff --git a/tests/test.hs b/tests/test.hs --- a/tests/test.hs +++ b/tests/test.hs @@ -1,7 +1,5 @@ import qualified Tests.Bed as Bed -import qualified Tests.Bam as Bam 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 qualified Tests.Tools as Tools @@ -10,9 +8,7 @@ main :: IO () main = defaultMain $ testGroup "Main" [ Bed.tests - , Bam.tests , Seq.tests - , ChIPSeq.tests , Motif.tests , GREAT.tests , Tools.tests