bioinformatics-toolkit 0.2.3 → 0.2.4
raw patch · 12 files changed
+48/−141 lines, 12 filesdep −samtools
Dependencies removed: samtools
Files
- bioinformatics-toolkit.cabal +4/−7
- src/Bio/ChIPSeq.hs +3/−4
- src/Bio/Data/Bam.hs +0/−62
- src/Bio/Data/Bed.hs +10/−3
- src/Bio/Motif.hs +8/−0
- src/Bio/RealWorld/GENCODE.hs +1/−1
- tests/Tests/Bam.hs +0/−21
- tests/Tests/ChIPSeq.hs +0/−29
- tests/Tests/Motif.hs +22/−10
- tests/data/example.bam binary
- tests/data/example.bam.bai binary
- tests/test.hs +0/−4
bioinformatics-toolkit.cabal view
@@ -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
src/Bio/ChIPSeq.hs view
@@ -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
− src/Bio/Data/Bam.hs
@@ -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 #-}
src/Bio/Data/Bed.hs view
@@ -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 #-}
src/Bio/Motif.hs view
@@ -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
src/Bio/RealWorld/GENCODE.hs view
@@ -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
− tests/Tests/Bam.hs
@@ -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'
− tests/Tests/ChIPSeq.hs
@@ -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'
tests/Tests/Motif.hs view
@@ -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
− tests/data/example.bam
binary file changed (689964 → absent bytes)
− tests/data/example.bam.bai
binary file changed (1973832 → absent bytes)
tests/test.hs view
@@ -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