bioinformatics-toolkit 0.5.1 → 0.6.0
raw patch · 17 files changed
+773/−523 lines, 17 filesdep +attoparsecdep +conduit-extradep ~HsHTSLibdep ~matrices
Dependencies added: attoparsec, conduit-extra
Dependency ranges changed: HsHTSLib, matrices
Files
- bioinformatics-toolkit.cabal +24/−18
- src/Bio/ChIPSeq.hs +0/−248
- src/Bio/Data/Bam.hs +43/−43
- src/Bio/Data/Bed.hs +0/−97
- src/Bio/Data/Bed/Types.hs +5/−1
- src/Bio/Data/Bed/Utils.hs +373/−0
- src/Bio/Data/Fastq.hs +48/−43
- src/Bio/Motif.hs +32/−31
- src/Bio/Motif/Merge.hs +1/−1
- src/Bio/RealWorld/GDC.hs +32/−0
- src/Bio/RealWorld/Reactome.hs +51/−0
- src/Bio/Utils/BitVector.hs +63/−0
- tests/Tests/Bam.hs +16/−7
- tests/Tests/Bed.hs +27/−2
- tests/Tests/ChIPSeq.hs +0/−31
- tests/Tests/Motif.hs +4/−1
- tests/data/motifs.meme +54/−0
bioinformatics-toolkit.cabal view
@@ -1,12 +1,12 @@ name: bioinformatics-toolkit-version: 0.5.1+version: 0.6.0 synopsis: A collection of bioinformatics tools description: A collection of bioinformatics tools license: MIT license-file: LICENSE author: Kai Zhang maintainer: kai@kzhang.org-copyright: (c) 2014-2018 Kai Zhang+copyright: (c) 2014-2019 Kai Zhang category: Bio build-type: Simple extra-source-files: README.md@@ -20,15 +20,16 @@ tests/data/peaks.sorted.bed tests/data/example_intersect_peaks.bed tests/data/motifs.fasta+ tests/data/motifs.meme library hs-source-dirs: src ghc-options: -Wall exposed-modules:- Bio.ChIPSeq Bio.ChIPSeq.FragLen Bio.Data.Bed Bio.Data.Bed.Types+ Bio.Data.Bed.Utils Bio.Data.Bam Bio.Data.Fasta Bio.Data.Fastq@@ -42,11 +43,14 @@ Bio.RealWorld.ENCODE Bio.RealWorld.Ensembl Bio.RealWorld.GENCODE+ Bio.RealWorld.GDC Bio.RealWorld.ID+ Bio.RealWorld.Reactome Bio.RealWorld.UCSC Bio.RealWorld.Uniprot Bio.Seq Bio.Seq.IO+ Bio.Utils.BitVector Bio.Utils.Functions Bio.Utils.Misc Bio.Utils.Overlap@@ -56,30 +60,32 @@ base >=4.11 && <5.0 , aeson , aeson-pretty- , bytestring >=0.10- , bytestring-lexing >=0.5+ , attoparsec+ , bytestring >= 0.10+ , bytestring-lexing >= 0.5 , case-insensitive , clustering- , conduit >=1.3.0- , containers >=0.5+ , conduit >= 1.3.0+ , conduit-extra+ , containers >= 0.5 , data-ordlist , data-default-class , double-conversion- , HsHTSLib >=1.3.2.3- , http-conduit >=2.1.8+ , HsHTSLib >= 1.9.2+ , http-conduit >= 2.1.8 , hexpat- , IntervalMap >=0.5.0.0+ , IntervalMap >= 0.5.0.0 , lens- , matrices >=0.4.3- , mtl >=2.1.3.1+ , matrices >= 0.5.0+ , mtl >= 2.1.3.1 , math-functions- , parallel >=3.2+ , parallel >= 3.2 , primitive , split- , statistics >=0.13.2.1- , text >=0.11- , transformers >=0.3.0.0- , unordered-containers >=0.2+ , statistics >= 0.13.2.1+ , text >= 0.11+ , transformers >= 0.3.0.0+ , unordered-containers >= 0.2 , word8 , vector , vector-algorithms@@ -108,7 +114,6 @@ other-modules: Tests.Bed , Tests.Bam- , Tests.ChIPSeq , Tests.Motif , Tests.Seq , Tests.Tools@@ -120,6 +125,7 @@ , random , vector , data-default-class+ , lens , tasty , tasty-golden , tasty-hunit
− src/Bio/ChIPSeq.hs
@@ -1,248 +0,0 @@-{-# LANGUAGE BangPatterns #-}-{-# LANGUAGE FlexibleContexts #-}-{-# LANGUAGE OverloadedStrings #-}--module Bio.ChIPSeq- ( monoColonalize- , rpkmBed- , rpkmSortedBed- , countTagsBinBed- , countTagsBinBed'- , tagCountDistr- , peakCluster- ) where--import Conduit-import Control.Monad (forM, forM_, liftM)-import Control.Monad.Primitive (PrimMonad)-import qualified Data.Foldable as F-import Data.Function (on)-import qualified Data.HashMap.Strict as M-import qualified Data.IntervalMap as IM-import Control.Lens ((^.), (&), (.~))-import Data.Maybe (fromJust, fromMaybe)-import qualified Data.Vector as V-import qualified Data.Vector.Algorithms.Intro as I-import qualified Data.Vector.Generic as G-import qualified Data.Vector.Generic.Mutable as GM-import qualified Data.Vector.Unboxed as U--import Bio.Data.Bed---- | process a sorted BED stream, keep only mono-colonal tags-monoColonalize :: Monad m => ConduitT BED BED m ()-monoColonalize = do- x <- headC- case x of- Just b -> yield b >> concatMapAccumC f b- Nothing -> return ()- where- f cur prev = case compareBed prev cur of- GT -> error $- "Input is not sorted: " ++ show prev ++ " > " ++ show cur- LT -> (cur, [cur])- _ -> if prev^.strand == cur^.strand then (cur, []) else (cur, [cur])-{-# INLINE monoColonalize #-}---- | calculate RPKM on a set of unique regions. Regions (in bed format) would be kept in--- memory but not tag file.--- RPKM: Readcounts per kilobase per million reads. Only counts the starts of tags-rpkmBed :: (PrimMonad m, BEDLike b, G.Vector v Double)- => [b] -> ConduitT BED o m (v Double)-rpkmBed regions = do- v <- lift $ do v' <- V.unsafeThaw . V.fromList . zip [0..] $ regions- I.sortBy (compareBed `on` snd) v'- V.unsafeFreeze v'- let (idx, sortedRegions) = V.unzip v- n = G.length idx- rc <- rpkmSortedBed $ Sorted sortedRegions-- lift $ do- result <- GM.new n- G.sequence_ . G.imap (\x i -> GM.unsafeWrite result i (rc U.! x)) $ idx- G.unsafeFreeze result-{-# INLINE rpkmBed #-}---- | calculate RPKM on a set of regions. Regions must be sorted. The Sorted data--- type is used to remind users to sort their data.-rpkmSortedBed :: (PrimMonad m, BEDLike b, G.Vector v Double)- => Sorted (V.Vector b) -> ConduitT BED o m (v Double)-rpkmSortedBed (Sorted regions) = do- 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)))- $ G.unsafeFreeze vec- where- count v nTags tag = do- let p | tag^.strand == Just True = tag^.chromStart- | tag^.strand == Just False = tag^.chromEnd - 1- | otherwise = error "Unkown strand"- xs = concat $ IM.elems $- IM.containing (M.lookupDefault IM.empty (tag^.chrom) intervalMap) p- addOne v xs- return $ succ nTags-- intervalMap = sortedBedToTree (++) . Sorted . G.toList . G.zip regions .- G.map return . G.enumFromN 0 $ l- addOne v' = mapM_ $ \x -> GM.unsafeRead v' x >>= GM.unsafeWrite v' x . (+1)- l = G.length regions-{-# INLINE rpkmSortedBed #-}---- | divide each region into consecutive bins, and count tags for each bin and--- return the number of all tags. Note: a tag is considered to be overlapped--- with a region only if the starting position of the tag is in the region. For--- the common sense overlapping, use countTagsBinBed'.-countTagsBinBed :: (Integral a, PrimMonad m, G.Vector v a, BEDLike b)- => Int -- ^ bin size- -> [b] -- ^ regions- -> ConduitT BED o m ([v a], Int)-countTagsBinBed k beds = do- initRC <- lift $ forM beds $ \bed -> do- let start = bed^.chromStart- end = bed^.chromEnd- num = (end - start) `div` k- index i = (i - start) `div` k- v <- GM.replicate num 0- return (v, index)-- sink 0 $ V.fromList initRC- where- sink !nTags vs = do- tag <- await- case tag of- Just bed -> do- let p | bed^.strand == Just True = bed^.chromStart- | bed^.strand == Just False = bed^.chromEnd - 1- | otherwise = error "profiling: unkown strand"- overlaps = concat $ IM.elems $- IM.containing (M.lookupDefault IM.empty (bed^.chrom) intervalMap) p- lift $ forM_ overlaps $ \x -> do- let (v, idxFn) = vs `G.unsafeIndex` x- i = let i' = idxFn p- l = GM.length v- in if i' >= l then l - 1 else i'- GM.unsafeRead v i >>= GM.unsafeWrite v i . (+1)- sink (nTags+1) vs-- _ -> do rc <- lift $ mapM (G.unsafeFreeze . fst) $ G.toList vs- return (rc, nTags)-- intervalMap = bedToTree (++) $ zip beds $ map return [0..]-{-# INLINE countTagsBinBed #-}---- | Same as countTagsBinBed, except that tags are treated as complete intervals--- instead of single points.-countTagsBinBed' :: (Integral a, PrimMonad m, G.Vector v a, BEDLike b1, BEDLike b2)- => Int -- ^ bin size- -> [b1] -- ^ regions- -> ConduitT b2 o m ([v a], Int)-countTagsBinBed' k beds = do- initRC <- lift $ forM beds $ \bed -> do- let start = bed^.chromStart- end = bed^.chromEnd- num = (end - start) `div` k- index i = (i - start) `div` k- v <- GM.replicate num 0- return (v, index)-- sink 0 $ V.fromList initRC- where- sink !nTags vs = do- tag <- await- case tag of- Just bed -> do- let chr = bed^.chrom- start = bed^.chromStart- end = bed^.chromEnd- overlaps = concat $ IM.elems $ IM.intersecting- (M.lookupDefault IM.empty chr intervalMap) $ IM.IntervalCO start end- lift $ forM_ overlaps $ \x -> do- let (v, idxFn) = vs `G.unsafeIndex` x- lo = let i = idxFn start- in if i < 0 then 0 else i- hi = let i = idxFn end- l = GM.length v- in if i >= l then l - 1 else i- forM_ [lo..hi] $ \i ->- GM.unsafeRead v i >>= GM.unsafeWrite v i . (+1)- sink (nTags+1) vs-- _ -> do rc <- lift $ mapM (G.unsafeFreeze . fst) $ G.toList vs- return (rc, nTags)-- intervalMap = bedToTree (++) $ zip beds $ map return [0..]-{-# 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-rpkmBam fl = do- nTags <- lift $ readBam fl $$ foldMC (\acc bam -> return $- if isUnmap bam then acc else acc + 1) 0.0- handle <- lift $ BI.open fl- conduit nTags handle- where- conduit n h = do- x <- await- case x of- Nothing -> lift $ BI.close h- Just bed -> do let chr = chrom bed- s = chromStart bed- e = chromEnd bed- rc <- lift $ viewBam h (chr, s, e) $$ readCount s e- yield $ rc * 1e9 / n / fromIntegral (e-s)- conduit n h- readCount l u = foldMC f 0.0- where- f acc bam = do let p1 = fromIntegral . fromJust . position $ bam- rl = fromIntegral . fromJust . queryLength $ bam- p2 = p1 + rl - 1- return $ if isReverse bam- then if l <= p2 && p2 < u then acc + 1- else acc- else if l <= p1 && p1 < u then acc + 1- else acc-{-# INLINE rpkmBam #-}--}--tagCountDistr :: PrimMonad m => G.Vector v Int => ConduitT BED o m (v Int)-tagCountDistr = loop M.empty- where- loop m = do- x <- await- case x of- Just bed -> do- let p | fromMaybe True (bed^.strand) = bed^.chromStart- | otherwise = 1 - bed^.chromEnd- case M.lookup (bed^.chrom) m of- Just table -> loop $ M.insert (bed^.chrom) (M.insertWith (+) p 1 table) m- _ -> loop $ M.insert (bed^.chrom) (M.fromList [(p,1)]) m- _ -> lift $ do- vec <- GM.replicate 100 0- F.forM_ m $ \table ->- F.forM_ table $ \v -> do- let i = min 99 v- GM.unsafeRead vec i >>= GM.unsafeWrite vec i . (+1)- G.unsafeFreeze vec-{-# INLINE tagCountDistr #-}---- | cluster peaks-peakCluster :: (BEDLike b, Monad m)- => [b] -- ^ peaks- -> Int -- ^ radius- -> Int -- ^ cutoff- -> ConduitT o BED m ()-peakCluster peaks r th = mergeBedWith mergeFn peaks' .| filterC g- where- peaks' = map f peaks- f b = let c = (b^.chromStart + b^.chromEnd) `div` 2- in asBed (b^.chrom) (c-r) (c+r) :: BED3- mergeFn xs = asBed (head xs ^. chrom) lo hi & score .~ Just (fromIntegral $ length xs)- where- lo = minimum $ map (^.chromStart) xs- hi = maximum $ map (^.chromEnd) xs- g b = fromJust (b^.score) >= fromIntegral th-{-# INLINE peakCluster #-}
src/Bio/Data/Bam.hs view
@@ -1,67 +1,67 @@ {-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE LambdaCase #-} module Bio.Data.Bam- ( Bam- , HeaderState- , withBamFile- , readBam- , writeBam+ ( BAM+ , getBamHeader+ , streamBam+ , sinkBam+ , bamToBedC , bamToBed+ , bamToFastqC+ , bamToFastq , sortedBamToBedPE ) where import Bio.Data.Bed+import Bio.Data.Fastq import Bio.HTS-import Bio.HTS.Types (Bam, FileHeader (..))+import Bio.HTS.Types (BAM, BAMHeader, SortOrder(..)) import Conduit import Control.Lens ((&), (.~))-import Control.Monad.Reader (ask, lift) -- | Convert bam record to bed record. Unmapped reads will be discarded.-bamToBed :: ConduitT Bam BED HeaderState ()-bamToBed = mapMC bamToBed1 .| concatC-{-# INLINE bamToBed #-}+bamToBedC :: MonadIO m => BAMHeader -> ConduitT BAM BED m ()+bamToBedC header = mapC (bamToBed header) .| concatC+{-# INLINE bamToBedC #-} +-- | Convert bam record to fastq record.+bamToFastqC :: Monad m => ConduitT BAM Fastq m ()+bamToFastqC = mapC bamToFastq .| concatC+{-# INLINE bamToFastqC #-}+ -- | Convert pairedend bam file to bed. the bam file must be sorted by names, -- e.g., using "samtools sort -n". This condition is checked from Bam header.-sortedBamToBedPE :: ConduitT Bam (BED, BED) HeaderState ()-sortedBamToBedPE = do- maybeBam <- await- case maybeBam of- Nothing -> return ()- Just b' -> do- leftover b'- sortOrd <- getSortOrder <$> lift ask- case sortOrd of- Queryname -> loopBedPE .| concatC- _ -> error "Bam file must be sorted by NAME."+sortedBamToBedPE :: Monad m => BAMHeader -> ConduitT BAM (BED, BED) m ()+sortedBamToBedPE header = case getSortOrder header of+ Queryname -> loopBedPE .| concatC+ _ -> error "Bam file must be sorted by NAME." where- loopBedPE = do- pair <- (,) <$$> await <***> await- case pair of- Nothing -> return ()- Just (bam1, bam2) -> if qName bam1 /= qName bam2- then error "Adjacent records have different query names. Aborted."- else do- bed1 <- lift $ bamToBed1 bam1- bed2 <- lift $ bamToBed1 bam2- yield $ (,) <$> bed1 <*> bed2- loopBedPE+ loopBedPE = (,) <$$> await <***> await >>= \case+ Nothing -> return ()+ Just (bam1, bam2) -> if queryName bam1 /= queryName bam2+ then error "Adjacent records have different query names. Aborted."+ else do+ yield $ (,) <$> bamToBed header bam1 <*> bamToBed header bam2+ loopBedPE where (<$$>) = fmap . fmap (<***>) = (<*>) . fmap (<*>) {-# INLINE sortedBamToBedPE #-} --bamToBed1 :: Bam -> HeaderState (Maybe BED)-bamToBed1 bam = do- BamHeader hdr <- lift ask- return $- (\chr -> asBed chr start end & name .~ nm & score .~ sc & strand .~ str)- <$> getChr hdr bam+-- | Convert BAM to BED.+bamToBed :: BAMHeader -> BAM -> Maybe BED+bamToBed header bam = mkBed <$> refName header bam where- start = fromIntegral $ position bam- end = fromIntegral $ endPos bam- nm = Just $ qName bam+ mkBed chr = asBed chr start end &+ name .~ nm & score .~ sc & strand .~ str+ start = startLoc bam+ end = endLoc bam+ nm = Just $ queryName bam str = Just $ not $ isRev bam sc = Just $ fromIntegral $ mapq bam-{-# INLINE bamToBed1 #-}+{-# INLINE bamToBed #-}++-- | Convert BAM to Fastq.+bamToFastq :: BAM -> Maybe Fastq+bamToFastq bam = Fastq (queryName bam) <$> getSeq bam <*> qualityS bam+{-# INLINE bamToFastq #-}
src/Bio/Data/Bed.hs view
@@ -25,7 +25,6 @@ , splitBedBySize , splitBedBySizeLeft , splitBedBySizeOverlap- , Sorted(..) , sortBed , intersectBed , intersectBedWith@@ -46,12 +45,6 @@ , writeBed , writeBed' - -- * Utilities- , fetchSeq- , fetchSeq'- , motifScan- , getMotifScore- , getMotifPValue , compareBed ) where @@ -65,18 +58,12 @@ import qualified Data.HashMap.Strict as M import qualified Data.IntervalMap.Strict as IM import Data.List (groupBy, sortBy)-import Data.Maybe (fromJust) import Data.Ord (comparing) import qualified Data.Vector as V import qualified Data.Vector.Algorithms.Intro as I import System.IO import Bio.Data.Bed.Types-import Bio.Motif (Bkgd (..), Motif (..))-import qualified Bio.Motif as Motif-import qualified Bio.Motif.Search as Motif-import Bio.Seq-import Bio.Seq.IO import Bio.Utils.Misc (binBySize, binBySizeLeft, binBySizeOverlap, bins) @@ -150,9 +137,6 @@ binBySizeOverlap k o (bed^.chromStart, bed^.chromEnd) {-# INLINE splitBedBySizeOverlap #-} --- | a type to imply that underlying data structure is sorted-newtype Sorted b = Sorted {fromSorted :: b} deriving (Show, Read, Eq)- -- | Compare bed records using only the chromosome, start and end positions. -- Unlike the ``compare'' from the Ord type class, this function can compare -- different types of BED data types.@@ -315,84 +299,3 @@ writeBed' :: (BEDConvert b, MonadIO m) => FilePath -> [b] -> m () writeBed' fl beds = runConduit $ yieldMany beds .| writeBed fl {-# INLINE writeBed' #-}---- | retreive sequences-fetchSeq :: (BioSeq DNA a, MonadIO m)- => Genome- -> ConduitT BED (Either String (DNA a)) m ()-fetchSeq g = mapMC f- where- f bed = do- dna <- liftIO $ getSeq g (bed^.chrom, bed^.chromStart, bed^.chromEnd)- return $ case bed^.strand of- Just False -> rc <$> dna- _ -> dna-{-# INLINE fetchSeq #-}--fetchSeq' :: (BioSeq DNA a, MonadIO m) => Genome -> [BED] -> m [Either String (DNA a)]-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 _ -> return ()- 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 #-}---- | 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- where- motifMap = M.fromListWith (error "found motif with same name") $- map (\(Motif nm pwm) -> (nm, pwm)) motifs-{-# INLINE getMotifScore #-}--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- 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 #-}
src/Bio/Data/Bed/Types.hs view
@@ -19,6 +19,7 @@ , _bed , _data , BEDTree+ , Sorted(..) ) where import Control.Lens@@ -135,7 +136,7 @@ | f6 == Just False = "-" | otherwise = "." score' = case f5 of- Just x -> (B.pack.show) x+ Just x -> toShortest x _ -> "." {-# INLINE toLine #-} @@ -314,3 +315,6 @@ {-# INLINE toLine #-} type BEDTree a = M.HashMap B.ByteString (IM.IntervalMap Int a)++-- | a type to imply that underlying data structure is sorted+newtype Sorted b = Sorted {fromSorted :: b} deriving (Show, Read, Eq)
+ src/Bio/Data/Bed/Utils.hs view
@@ -0,0 +1,373 @@+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE BangPatterns #-}++module Bio.Data.Bed.Utils+ ( fetchSeq+ , fetchSeq'+ , motifScan+ , getMotifScore+ , getMotifPValue+ , monoColonalize+ , BaseMap(..)+ , baseMap+ , queryBaseMap+ , rpkmBed+ , rpkmSortedBed+ , countTagsBinBed+ , countTagsBinBed'+ , tagCountDistr+ , peakCluster+ ) where++import Conduit+import Control.Lens+import Control.Monad.State.Strict+import qualified Data.ByteString.Char8 as B+import qualified Data.Foldable as F+import Data.Function (on)+import qualified Data.HashMap.Strict as M+import qualified Data.IntervalMap.Strict as IM+import Data.Maybe (fromJust, fromMaybe)+import qualified Data.Vector as V+import qualified Data.Vector.Algorithms.Intro as I+import qualified Data.Vector.Generic as G+import qualified Data.Vector.Generic.Mutable as GM+import qualified Data.Vector.Unboxed as U++import Bio.Data.Bed+import Bio.Data.Bed.Types+import Bio.Motif (Bkgd (..), Motif (..))+import qualified Bio.Motif as Motif+import qualified Bio.Motif.Search as Motif+import Bio.Seq hiding (length)+import Bio.Seq.IO+import qualified Bio.Utils.BitVector as BV+++-- | retreive sequences+fetchSeq :: (BioSeq DNA a, MonadIO m)+ => Genome+ -> ConduitT BED (Either String (DNA a)) m ()+fetchSeq g = mapMC f+ where+ f bed = do+ dna <- liftIO $ getSeq g (bed^.chrom, bed^.chromStart, bed^.chromEnd)+ return $ case bed^.strand of+ Just False -> rc <$> dna+ _ -> dna+{-# INLINE fetchSeq #-}++fetchSeq' :: (BioSeq DNA a, MonadIO m) => Genome -> [BED] -> m [Either String (DNA a)]+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 _ -> return ()+ 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 #-}++-- | 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+ where+ motifMap = M.fromListWith (error "found motif with same name") $+ map (\(Motif nm pwm) -> (nm, pwm)) motifs+{-# INLINE getMotifScore #-}++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+ 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 #-}++-- | process a sorted BED stream, keep only mono-colonal tags+monoColonalize :: Monad m => ConduitT BED BED m ()+monoColonalize = do+ x <- headC+ case x of+ Just b -> yield b >> concatMapAccumC f b+ Nothing -> return ()+ where+ f cur prev = case compareBed prev cur of+ GT -> error $+ "Input is not sorted: " ++ show prev ++ " > " ++ show cur+ LT -> (cur, [cur])+ _ -> if prev^.strand == cur^.strand then (cur, []) else (cur, [cur])+{-# INLINE monoColonalize #-}++newtype BaseMap = BaseMap (M.HashMap B.ByteString BV.BitVector)++-- | Count the tags (starting positions) at each position in the genome.+baseMap :: PrimMonad m+ => [(B.ByteString, Int)] -- ^ chromosomes and their sizes+ -> ConduitT BED o m BaseMap+baseMap chrs = do+ bvs <- lift $ fmap M.fromList $ forM chrs $ \(chr, n) -> do+ bv <- BV.zeros n+ return (chr, bv)++ mapM_C $ \bed -> case M.lookup (bed^.chrom) bvs of+ Nothing -> return ()+ Just bv -> if fromMaybe True $ bed^.strand+ then BV.set bv $ bed^.chromStart+ else BV.set bv $ bed^.chromEnd++ lift $ fmap BaseMap $ sequence $ fmap BV.unsafeFreeze bvs ++queryBaseMap :: BEDLike b => b -> BaseMap -> Maybe [Bool]+queryBaseMap bed (BaseMap bm) = case M.lookup (bed^.chrom) bm of+ Nothing -> Nothing+ Just bv ->+ let res = map (bv BV.!) [bed^.chromStart .. bed^.chromEnd - 1]+ in case bed^.strand of+ Just False -> Just $ reverse res+ _ -> Just res++-- | calculate RPKM on a set of unique regions. Regions (in bed format) would be kept in+-- memory but not tag file.+-- RPKM: Readcounts per kilobase per million reads. Only counts the starts of tags+rpkmBed :: (PrimMonad m, BEDLike b, G.Vector v Double)+ => [b] -> ConduitT BED o m (v Double)+rpkmBed regions = do+ v <- lift $ do v' <- V.unsafeThaw . V.fromList . zip [0..] $ regions+ I.sortBy (compareBed `on` snd) v'+ V.unsafeFreeze v'+ let (idx, sortedRegions) = V.unzip v+ n = G.length idx+ rc <- rpkmSortedBed $ Sorted sortedRegions++ lift $ do+ result <- GM.new n+ G.sequence_ . G.imap (\x i -> GM.unsafeWrite result i (rc U.! x)) $ idx+ G.unsafeFreeze result+{-# INLINE rpkmBed #-}++-- | calculate RPKM on a set of regions. Regions must be sorted. The Sorted data+-- type is used to remind users to sort their data.+rpkmSortedBed :: (PrimMonad m, BEDLike b, G.Vector v Double)+ => Sorted (V.Vector b) -> ConduitT BED o m (v Double)+rpkmSortedBed (Sorted regions) = do+ 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)))+ $ G.unsafeFreeze vec+ where+ count v nTags tag = do+ let p | tag^.strand == Just True = tag^.chromStart+ | tag^.strand == Just False = tag^.chromEnd - 1+ | otherwise = error "Unkown strand"+ xs = concat $ IM.elems $+ IM.containing (M.lookupDefault IM.empty (tag^.chrom) intervalMap) p+ addOne v xs+ return $ succ nTags++ intervalMap = sortedBedToTree (++) . Sorted . G.toList . G.zip regions .+ G.map return . G.enumFromN 0 $ l+ addOne v' = mapM_ $ \x -> GM.unsafeRead v' x >>= GM.unsafeWrite v' x . (+1)+ l = G.length regions+{-# INLINE rpkmSortedBed #-}++-- | divide each region into consecutive bins, and count tags for each bin and+-- return the number of all tags. Note: a tag is considered to be overlapped+-- with a region only if the starting position of the tag is in the region. For+-- the common sense overlapping, use countTagsBinBed'.+countTagsBinBed :: (Integral a, PrimMonad m, G.Vector v a, BEDLike b)+ => Int -- ^ bin size+ -> [b] -- ^ regions+ -> ConduitT BED o m ([v a], Int)+countTagsBinBed k beds = do+ initRC <- lift $ forM beds $ \bed -> do+ let start = bed^.chromStart+ end = bed^.chromEnd+ num = (end - start) `div` k+ index i = (i - start) `div` k+ v <- GM.replicate num 0+ return (v, index)++ sink 0 $ V.fromList initRC+ where+ sink !nTags vs = do+ tag <- await+ case tag of+ Just bed -> do+ let p | bed^.strand == Just True = bed^.chromStart+ | bed^.strand == Just False = bed^.chromEnd - 1+ | otherwise = error "profiling: unkown strand"+ overlaps = concat $ IM.elems $+ IM.containing (M.lookupDefault IM.empty (bed^.chrom) intervalMap) p+ lift $ forM_ overlaps $ \x -> do+ let (v, idxFn) = vs `G.unsafeIndex` x+ i = let i' = idxFn p+ l = GM.length v+ in if i' >= l then l - 1 else i'+ GM.unsafeRead v i >>= GM.unsafeWrite v i . (+1)+ sink (nTags+1) vs++ _ -> do rc <- lift $ mapM (G.unsafeFreeze . fst) $ G.toList vs+ return (rc, nTags)++ intervalMap = bedToTree (++) $ zip beds $ map return [0..]+{-# INLINE countTagsBinBed #-}++-- | Same as countTagsBinBed, except that tags are treated as complete intervals+-- instead of single points.+countTagsBinBed' :: (Integral a, PrimMonad m, G.Vector v a, BEDLike b1, BEDLike b2)+ => Int -- ^ bin size+ -> [b1] -- ^ regions+ -> ConduitT b2 o m ([v a], Int)+countTagsBinBed' k beds = do+ initRC <- lift $ forM beds $ \bed -> do+ let start = bed^.chromStart+ end = bed^.chromEnd+ num = (end - start) `div` k+ index i = (i - start) `div` k+ v <- GM.replicate num 0+ return (v, index)++ sink 0 $ V.fromList initRC+ where+ sink !nTags vs = do+ tag <- await+ case tag of+ Just bed -> do+ let chr = bed^.chrom+ start = bed^.chromStart+ end = bed^.chromEnd+ overlaps = concat $ IM.elems $ IM.intersecting+ (M.lookupDefault IM.empty chr intervalMap) $ IM.IntervalCO start end+ lift $ forM_ overlaps $ \x -> do+ let (v, idxFn) = vs `G.unsafeIndex` x+ lo = let i = idxFn start+ in if i < 0 then 0 else i+ hi = let i = idxFn end+ l = GM.length v+ in if i >= l then l - 1 else i+ forM_ [lo..hi] $ \i ->+ GM.unsafeRead v i >>= GM.unsafeWrite v i . (+1)+ sink (nTags+1) vs++ _ -> do rc <- lift $ mapM (G.unsafeFreeze . fst) $ G.toList vs+ return (rc, nTags)++ intervalMap = bedToTree (++) $ zip beds $ map return [0..]+{-# 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+rpkmBam fl = do+ nTags <- lift $ readBam fl $$ foldMC (\acc bam -> return $+ if isUnmap bam then acc else acc + 1) 0.0+ handle <- lift $ BI.open fl+ conduit nTags handle+ where+ conduit n h = do+ x <- await+ case x of+ Nothing -> lift $ BI.close h+ Just bed -> do let chr = chrom bed+ s = chromStart bed+ e = chromEnd bed+ rc <- lift $ viewBam h (chr, s, e) $$ readCount s e+ yield $ rc * 1e9 / n / fromIntegral (e-s)+ conduit n h+ readCount l u = foldMC f 0.0+ where+ f acc bam = do let p1 = fromIntegral . fromJust . position $ bam+ rl = fromIntegral . fromJust . queryLength $ bam+ p2 = p1 + rl - 1+ return $ if isReverse bam+ then if l <= p2 && p2 < u then acc + 1+ else acc+ else if l <= p1 && p1 < u then acc + 1+ else acc+{-# INLINE rpkmBam #-}+-}++tagCountDistr :: PrimMonad m => G.Vector v Int => ConduitT BED o m (v Int)+tagCountDistr = loop M.empty+ where+ loop m = do+ x <- await+ case x of+ Just bed -> do+ let p | fromMaybe True (bed^.strand) = bed^.chromStart+ | otherwise = 1 - bed^.chromEnd+ case M.lookup (bed^.chrom) m of+ Just table -> loop $ M.insert (bed^.chrom) (M.insertWith (+) p 1 table) m+ _ -> loop $ M.insert (bed^.chrom) (M.fromList [(p,1)]) m+ _ -> lift $ do+ vec <- GM.replicate 100 0+ F.forM_ m $ \table ->+ F.forM_ table $ \v -> do+ let i = min 99 v+ GM.unsafeRead vec i >>= GM.unsafeWrite vec i . (+1)+ G.unsafeFreeze vec+{-# INLINE tagCountDistr #-}++-- | cluster peaks+peakCluster :: (BEDLike b, Monad m)+ => [b] -- ^ peaks+ -> Int -- ^ radius+ -> Int -- ^ cutoff+ -> ConduitT o BED m ()+peakCluster peaks r th = mergeBedWith mergeFn peaks' .| filterC g+ where+ peaks' = map f peaks+ f b = let c = (b^.chromStart + b^.chromEnd) `div` 2+ in asBed (b^.chrom) (c-r) (c+r) :: BED3+ mergeFn xs = asBed (head xs ^. chrom) lo hi & score .~ Just (fromIntegral $ length xs)+ where+ lo = minimum $ map (^.chromStart) xs+ hi = maximum $ map (^.chromEnd) xs+ g b = fromJust (b^.score) >= fromIntegral th+{-# INLINE peakCluster #-}
src/Bio/Data/Fastq.hs view
@@ -1,16 +1,21 @@+{-# LANGUAGE LambdaCase #-}+{-# LANGUAGE OverloadedStrings #-} module Bio.Data.Fastq ( Fastq(..) , parseFastqC- , parseFastqUnsafeC , fastqToByteString+ , qualitySummary , trimPolyA ) where import Conduit import Control.Monad (when) import qualified Data.ByteString.Char8 as B+import qualified Data.ByteString.Lazy as BL import qualified Data.ByteString as BS import Data.Maybe (isJust)+import Data.Attoparsec.ByteString+import Data.Conduit.Attoparsec -- | A FASTQ file normally uses four lines per sequence. --@@ -27,54 +32,54 @@ data Fastq = Fastq { fastqSeqId :: B.ByteString , fastqSeq :: B.ByteString- , fastqSeqInfo :: B.ByteString , fastqSeqQual :: B.ByteString } deriving (Show, Eq) -parseFastqC :: Monad m => ConduitT B.ByteString Fastq m ()-parseFastqC = linesUnboundedAsciiC .| conduit- where- conduit = do- l1 <- await- l2 <- await- l3 <- await- l4 <- await- case mkFastqRecord <$> l1 <*> l2 <*> l3 <*> l4 of- Nothing -> when (isJust l1) $ error "file ends prematurely"- Just x -> yield x >> conduit+parseFastqC :: MonadThrow m => ConduitT B.ByteString Fastq m ()+parseFastqC = conduitParser fastqParser .| mapC snd {-# INLINE parseFastqC #-} -parseFastqUnsafeC :: Monad m => ConduitT B.ByteString Fastq m ()-parseFastqUnsafeC = linesUnboundedAsciiC .| conduit- where- conduit = do- l1 <- await- l2 <- await- l3 <- await- l4 <- await- case mkFastqRecordUnsafe <$> l1 <*> l2 <*> l3 <*> l4 of- Nothing -> when (isJust l1) $ error "file ends prematurely"- Just x -> yield x >> conduit-{-# INLINE parseFastqUnsafeC #-}+fastqParser :: Parser Fastq+fastqParser = do+ ident <- word8 64 *> takeTill (==10)+ sequence <- BS.filter (/=10) <$> takeTill (==43)+ skip (/=10)+ score <- BS.filter (/=10) <$> takeTill (==64)+ return $ Fastq ident sequence score+{-# INLINE fastqParser #-} -fastqToByteString :: Fastq -> [B.ByteString]-fastqToByteString (Fastq a b c d) = ['@' `B.cons` a, b, '+' `B.cons` c, d]+fastqToByteString :: Fastq -> B.ByteString+fastqToByteString (Fastq a b c) = "@" <> a <> "\n" <> b <> "\n+\n" <> c {-# INLINE fastqToByteString #-} --- | Make Fastq record from Bytestrings, without sanity check.-mkFastqRecordUnsafe :: B.ByteString -- ^ First line- -> B.ByteString -- ^ Second line- -> B.ByteString -- ^ Third line- -> B.ByteString -- ^ Fourth line- -> Fastq-mkFastqRecordUnsafe l1 l2 l3 l4 = Fastq (B.tail l1) l2 (B.tail l3) l4-{-# INLINE mkFastqRecordUnsafe #-}+-- | Get the mean and variance of quality scores at every position.+qualitySummary :: Monad m => ConduitT Fastq o m [(Double, Double)]+qualitySummary = mapC (map fromIntegral . decodeQualSc) .| meanVarianceC -mkFastqRecord :: B.ByteString -- ^ First line- -> B.ByteString -- ^ Second line- -> B.ByteString -- ^ Third line- -> B.ByteString -- ^ Fourth line- -> Fastq+meanVarianceC :: Monad m => ConduitT [Double] o m [(Double, Double)]+meanVarianceC = peekC >>= \case+ Nothing -> error "Empty input"+ Just x -> fst <$> foldlC f (replicate (length x) (0,0), 0)+ where+ f (acc, n) xs = let acc' = zipWith g acc xs in (acc', n')+ where+ n' = n + 1+ g (m, s) x = (m', s')+ where+ m' = m + d / fromIntegral n'+ s' = s + d * (x - m')+ d = x - m+{-# INLINE meanVarianceC #-}++decodeQualSc :: Fastq -> [Int]+decodeQualSc = map (fromIntegral . (\x -> x - 33)) . BS.unpack .fastqSeqQual+{-# INLINE decodeQualSc #-}++pError :: Int -> Double+pError x = 10 ** (negate (fromIntegral x) / 10)+{-# INLINE pError #-}++{- mkFastqRecord l1 l2 l3 l4 = Fastq (parseLine1 l1) (parseLine2 l2) (parseLine3 l3) (parseLine4 l4) where@@ -96,12 +101,12 @@ where f b = let b' = fromIntegral b :: Int in b' >= 33 && b' <= 126-{-# INLINE mkFastqRecord #-}+-} -- | Remove trailing 'A' trimPolyA :: Int -> Fastq -> Fastq-trimPolyA n f@(Fastq a b c d)- | B.length trailing >= n = Fastq a b' c $ B.take (B.length b') d+trimPolyA n f@(Fastq a b c)+ | B.length trailing >= n = Fastq a b' $ B.take (B.length b') c | otherwise = f where (b', trailing) = B.spanEnd (=='A') b
src/Bio/Motif.hs view
@@ -49,6 +49,7 @@ import Numeric.MathFunctions.Constants (m_epsilon) import Prelude hiding (sum) import Text.Printf (printf)+import Statistics.Function (minMax) import Bio.Seq import Bio.Utils.Functions (binarySearchBy)@@ -270,44 +271,44 @@ truncateCDF x (CDF v) = CDF $ U.filter ((>=x) . snd) v {-# INLINE truncateCDF #-} --- | Approximate the cdf of motif matching scores+-- | Approximate the cdf of motif matching scores using dynamic programming.+-- Algorithm:+-- Scan the PWM from left to right. For each position $i$, compute a score+-- density function $s_i$ such that $s_i(x)$ is the total number of sequences+-- with score $x$. scoreCDF :: Bkgd -> PWM -> CDF scoreCDF (BG (a,c,g,t)) pwm = toCDF $ loop (U.singleton 1, const 0) 0 where loop (prev,scFn) i- | i < n =- let (lo,hi) = minMax (1/0,-1/0) 0- 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- v = U.create $ do+ | i >= n = (prev, scFn)+ | lo < hi =+ let v = U.create $ do new <- UM.replicate nBin' 0- flip U.imapM_ prev $ \x p ->- when (p /= 0) $ do- let idx_a = idx $ sc + log' (M.unsafeIndex (_mat pwm) (i,0)) - log a- idx_c = idx $ sc + log' (M.unsafeIndex (_mat pwm) (i,1)) - log c- idx_g = idx $ sc + log' (M.unsafeIndex (_mat pwm) (i,2)) - log g- idx_t = idx $ sc + log' (M.unsafeIndex (_mat pwm) (i,3)) - log t- sc = scFn x- new `UM.read` idx_a >>= UM.write new idx_a . (a * p + )- new `UM.read` idx_c >>= UM.write new idx_c . (c * p + )- new `UM.read` idx_g >>= UM.write new idx_g . (g * p + )- new `UM.read` idx_t >>= UM.write new idx_t . (t * p + )+ flip U.imapM_ prev $ \x prob ->+ when (prob /= 0) $ do+ let sc = scFn x+ UM.modify new (a * prob +) $ getIdx $ sc + a_at i+ UM.modify new (c * prob +) $ getIdx $ sc + c_at i+ UM.modify new (g * prob +) $ getIdx $ sc + g_at i+ UM.modify new (t * prob +) $ getIdx $ sc + t_at i return new in loop (v, \x -> (fromIntegral x + 0.5) * step + lo) (i+1)- | otherwise = (prev, scFn)+ | otherwise = loop (prev, scFn) (i+1) where- minMax (l,h) x- | x >= U.length prev = (l,h)- | prev U.! x /= 0 =- let sc = scFn x- s1 = sc + log' (M.unsafeIndex (_mat pwm) (i,0)) - log a- s2 = sc + log' (M.unsafeIndex (_mat pwm) (i,1)) - log c- s3 = sc + log' (M.unsafeIndex (_mat pwm) (i,2)) - log g- s4 = sc + log' (M.unsafeIndex (_mat pwm) (i,3)) - log t- in minMax (foldr min l [s1,s2,s3,s4],foldr max h [s1,s2,s3,s4]) (x+1)- | otherwise = minMax (l,h) (x+1)+ getIdx x = let j = truncate $ (x - lo) / step+ in if j >= nBin' then nBin' - 1 else j+ lo = lo' + min'+ hi = hi' + max'+ nBin' = min 200000 $ ceiling $ (hi - lo) / precision+ step = (hi - lo) / fromIntegral nBin'+ lo' = scFn $ fst . U.head $ U.dropWhile ((==0) . snd) $ U.indexed prev+ hi' = scFn $ fst . U.head $ U.dropWhile ((==0) . snd) $ U.reverse $+ U.indexed prev+ (min', max') = minMax $ U.fromList [a_at i, c_at i, g_at i, t_at i]+ a_at i = log' (M.unsafeIndex (_mat pwm) (i,0)) - log a+ c_at i = log' (M.unsafeIndex (_mat pwm) (i,1)) - log c+ g_at i = log' (M.unsafeIndex (_mat pwm) (i,2)) - log g+ t_at i = log' (M.unsafeIndex (_mat pwm) (i,3)) - log t toCDF (v, scFn) = CDF $ compressCDF $ U.imap (\i x -> (scFn i, x)) $ U.scanl1 (+) v compressCDF v = U.ifilter f v where@@ -315,7 +316,7 @@ f i (_, x) | i == 0 || i == len = True | otherwise = x - snd (v `U.unsafeIndex` (i-1)) > m_epsilon || snd (v `U.unsafeIndex` (i+1)) - x > m_epsilon- precision = 1e-4+ precision = 1e-5 n = size pwm log' x | x == 0 = log 0.001 | otherwise = log x
src/Bio/Motif/Merge.hs view
@@ -18,7 +18,7 @@ import Control.Monad.ST (runST, ST) import qualified Data.ByteString.Char8 as B import Data.List (dropWhileEnd)-import qualified Data.Matrix.Symmetric.Mutable as MSU+import qualified Data.Matrix.Symmetric.Generic.Mutable as MSU import qualified Data.Matrix.Unboxed as MU import Data.Maybe import qualified Data.Vector as V
+ src/Bio/RealWorld/GDC.hs view
@@ -0,0 +1,32 @@+-- NIH Genomic Data Commons+{-# LANGUAGE OverloadedStrings #-}++module Bio.RealWorld.GDC+ (downloadData) where++import Conduit+import qualified Data.Text as T+import qualified Data.ByteString.Char8 as B+import Data.Maybe (fromJust)+import Network.HTTP.Conduit++baseurl :: String+baseurl = "https://api.gdc.cancer.gov/"++-- | Download data+downloadData :: String -- ^ UUID+ -> FilePath -- ^ Output dir+ -> IO FilePath+downloadData uuid dir = do+ request <- parseRequest url+ manager <- newManager tlsManagerSettings+ runResourceT $ do+ response <- http request manager+ let filename = T.unpack $ snd $ T.breakOnEnd "filename=" $ T.pack $+ B.unpack $ fromJust $ lookup "Content-Disposition" $+ responseHeaders response+ runConduit $ responseBody response .| sinkFileBS (dir ++ "/" ++ filename)+ return filename+ where+ url = baseurl ++ "data/" ++ uuid+{-# INLINE downloadData #-}
+ src/Bio/RealWorld/Reactome.hs view
@@ -0,0 +1,51 @@+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE DeriveGeneric #-}++module Bio.RealWorld.Reactome+ ( getPathways+ ) where++import Data.Aeson+import GHC.Generics (Generic)+import qualified Data.Text as T+import Network.HTTP.Simple++base :: String+base = "https://reactome.org/ContentService"++data Obj = Obj+ { className :: Maybe T.Text+ , dbId :: Int + , displayName :: T.Text+ , schemaClass :: Maybe T.Text+ , stId :: Maybe T.Text+ , stIdVersion :: Maybe T.Text+ } deriving (Show, Generic)++instance ToJSON Obj+instance FromJSON Obj+ +-- | All Reactome top level pathways+getPathways :: String -> IO [Obj]+getPathways species = do+ req <- parseRequest url+ response <- httpJSON req+ return $ getResponseBody response+ where+ url = base ++ "/data/pathways/top/" ++ species++{-+pathwayAnalysis :: [B.ByteString] -- ^ A list of identifiers+ -> IO B.ByteString+pathwayAnalysis ids = do+ initReq <- parseRequest base+ let request = urlEncodedBody [] initReq+ manager <- newManager tlsManagerSettings+ r <- fmap M.fromList $ runResourceT $ do+ response <- http request manager+ runConduit $ responseBody response .| linesUnboundedAsciiC .|+ (dropC 1 >> mapC ((\[a,b] -> (a,b)) . B.split '\t')) .| sinkList+ return $ map (flip M.lookup r) ids+ where+ params =+-}
+ src/Bio/Utils/BitVector.hs view
@@ -0,0 +1,63 @@+module Bio.Utils.BitVector+ ( BitVector+ , BitMVector+ , size+ , (!)+ , set+ , clear+ , unsafeFreeze+ , zeros+ , toList+ ) where++import qualified Data.Vector.Unboxed as U+import Control.Monad.Primitive+import qualified Data.Vector.Unboxed.Mutable as UM+import Data.Word+import Data.Bits+import Text.Printf (printf)++data BitVector = BitVector Int (U.Vector Word8)++data BitMVector s = BitMVector Int (UM.MVector s Word8)++size :: BitVector -> Int+size (BitVector n _) = n++(!) :: BitVector -> Int -> Bool+(!) = index++index :: BitVector -> Int -> Bool+index (BitVector n v) idx+ | idx >= n = error $ printf "index out of bounds (%d,%d)" idx n+ | otherwise = testBit (v `U.unsafeIndex` i) j+ where+ i = idx `div` 8+ j = idx `mod` 8++set :: PrimMonad m => BitMVector (PrimState m) -> Int -> m ()+set (BitMVector _ mv) idx = UM.modify mv ((flip setBit) j) i+ where+ i = idx `div` 8+ j = idx `mod` 8++clear :: PrimMonad m => BitMVector (PrimState m) -> Int -> m ()+clear (BitMVector _ mv) idx = UM.modify mv ((flip clearBit) j) i+ where+ i = idx `div` 8+ j = idx `mod` 8++unsafeFreeze :: PrimMonad m => BitMVector (PrimState m) -> m BitVector+unsafeFreeze (BitMVector n mv) = U.unsafeFreeze mv >>= return . BitVector n++zeros :: PrimMonad m => Int -> m (BitMVector (PrimState m))+zeros n = UM.replicate n' 0 >>= return . BitMVector n+ where+ n' = if j == 0 then i else i + 1+ i = n `div` 8+ j = n `mod` 8++toList :: BitVector -> [Bool]+toList bv = flip map [0..n-1] $ \i -> bv ! i+ where+ n = size bv
tests/Tests/Bam.hs view
@@ -24,23 +24,32 @@ bamIOTest = do goldenVsFile "BAM Read/Write Test" input output io where- io = withBamFile input $ \h -> runConduit $ readBam h .| writeBam output+ io = do+ header <- getBamHeader input+ runResourceT $ runConduit $ streamBam input .| sinkBam output header input = "tests/data/example.bam" output = "tests/data/example_copy.bam" bamToBedTest :: Assertion bamToBedTest = do bed <- readBed' "tests/data/example.bed"- bed' <- withBamFile "tests/data/example.bam" $ \h ->- runConduit $ readBam h .| bamToBed .| sinkList- (bed == bed') @? "bamToBedTest"+ bed' <- do+ let input = "tests/data/example.bam" + header <- getBamHeader input+ runResourceT $ runConduit $ streamBam input .| bamToBedC header .| sinkList+ forM_ (zip bed bed') $ \(a,b) -> + if a == b then return () else error $ show (a,b)+ (bed == bed') @? show (head bed, head bed') sortedBamToBedPETest :: Assertion sortedBamToBedPETest = do bedpe <- readBedPE "tests/data/pairedend.bedpe" :: IO [(BED3, BED3)]- bedpe' <- withBamFile "tests/data/pairedend.bam" $ \h -> runConduit $- readBam h .| sortedBamToBedPE .|- mapC (\(x,y) -> (convert x, convert y)) .| sinkList+ bedpe' <- do+ let input = "tests/data/pairedend.bam"+ header <- getBamHeader input+ runResourceT $ runConduit $ streamBam input .|+ sortedBamToBedPE header .|+ mapC (\(x,y) -> (convert x, convert y)) .| sinkList forM_ (zip bedpe bedpe') $ \(b1, b2) -> (b1 == b2 || b1 == swap b2) @? show (b1,b2) where readBedPE fl = do
tests/Tests/Bed.hs view
@@ -3,12 +3,19 @@ module Tests.Bed (tests) where import Bio.Data.Bed+import Bio.Data.Bed.Types+import Bio.Data.Bed.Utils+import Bio.Utils.BitVector+import Control.Lens import Conduit import Data.Function (on)-import Data.List (sortBy)+import Data.List (sortBy, sort) import qualified Data.Vector as V import Test.Tasty import Test.Tasty.HUnit+import qualified Data.HashMap.Strict as M+import Data.Ord+import Data.Maybe tests :: TestTree tests = testGroup "Test: Bio.Data.Bed"@@ -16,6 +23,7 @@ , testCase "split" splitBedTest , testCase "splitOverlapped" splitOverlappedTest , testCase "intersectBed" intersectBedTest+ , testCase "baseMap" baseMapTest ] sortBedTest :: Assertion@@ -67,5 +75,22 @@ intersectBedTest = do expect <- readBed' "tests/data/example_intersect_peaks.bed" :: IO [BED3] peaks <- readBed' "tests/data/peaks.bed" :: IO [BED3]- result <- readBed "tests/data/example.bed" =$= intersectBed peaks $$ sinkList+ result <- runConduit $ readBed "tests/data/example.bed" .| intersectBed peaks .| sinkList expect @=? result++baseMapTest :: Assertion+baseMapTest = do+ BaseMap bv <- runConduit $ readBed "tests/data/example.bed" .|+ baseMap [("chr1", 300000000)]+ let res = M.lookupDefault undefined "chr1" $+ fmap (map fst . filter snd . zip [0..] . toList) bv+ expect <- runConduit $ readBed "tests/data/example.bed" .|+ concatMapC f .| sinkList+ sort expect @=? sort res+ where+ f :: BED -> Maybe Int+ f bed = if bed^.chrom == "chr1"+ then Just $ if fromJust (bed^.strand)+ then bed^.chromStart+ else bed^.chromEnd+ else Nothing
− tests/Tests/ChIPSeq.hs
@@ -1,31 +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"- [ - ]--{--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
@@ -31,7 +31,10 @@ _ -> undefined motifs :: IO [Motif]-motifs = readFasta' "tests/data/motifs.fasta"+motifs = do+ m1 <- readFasta' "tests/data/motifs.fasta"+ -- m2 <- readMEME "tests/data/motifs.meme"+ return m1 tests :: TestTree tests = testGroup "Test: Bio.Motif"
+ tests/data/motifs.meme view
@@ -0,0 +1,54 @@+MEME version 4++ALPHABET= ACGT++strands: + -++Background letter frequencies+A 0.303 C 0.183 G 0.209 T 0.306 ++MOTIF lexA+letter-probability matrix: alength= 4 w= 12 nsites= 14 E= 3.2e-035 + 0.000000 0.000000 1.000000 0.000000 + 0.000000 0.000000 0.000000 1.000000 + 0.857143 0.000000 0.071429 0.071429 + 0.000000 0.071429 0.000000 0.928571 + 0.857143 0.000000 0.071429 0.071429 + 0.142857 0.000000 0.000000 0.857143 + 0.571429 0.071429 0.214286 0.142857 + 0.285714 0.285714 0.000000 0.428571 + 1.000000 0.000000 0.000000 0.000000 + 0.285714 0.214286 0.000000 0.500000 + 0.428571 0.500000 0.000000 0.071429 + 0.000000 1.000000 0.000000 0.000000 ++MOTIF YPL133C+757.pfm+letter-probability matrix: alength= 4 w= 12 nsites= 1 E= 0+ 0.25 0.25 0.25 0.25+ 0 0.5 0 0.5+ 0 1 0 0+ 0 1 0 0+ 0.333333 0 0.333333 0.333333+ 0 0 0 1+ 0 0 0 1+ 1 0 0 0+ 0 1 0 0+ 0 1 0 0+ 0 0 1 0+ 0.25 0.25 0.25 0.25++MOTIF 7 V_ELK1_01++letter-probability matrix: alength= 4 w= 12 nsites= 4 E= 0+ 0.250000 0.250000 0.250000 0.250000 + 0.500000 0.000000 0.250000 0.250000 + 0.500000 0.250000 0.250000 0.000000 + 0.750000 0.000000 0.250000 0.000000 + 0.000000 0.750000 0.250000 0.000000 + 0.000000 0.000000 1.000000 0.000000 + 1.000000 0.000000 0.000000 0.000000 + 0.000000 0.250000 0.000000 0.750000 + 0.250000 0.250000 0.250000 0.250000 + 0.250000 0.750000 0.000000 0.000000 + 0.250000 0.250000 0.500000 0.000000 + 0.250000 0.250000 0.000000 0.500000