bioinformatics-toolkit 0.4.1 → 0.5.0

raw patch · 20 files changed

+523/−705 lines, 20 filesdep +lensPVP ok

version bump matches the API change (PVP)

Dependencies added: lens

API changes (from Hackage documentation)

Files

• bioinformatics-toolkit.cabal +3/−4
• src/Bio/ChIPSeq.hs +60/−67
• src/Bio/ChIPSeq/FragLen.hs +12/−11
• src/Bio/Data/Bam.hs +9/−7
• src/Bio/Data/Bed.hs +142/−361
• src/Bio/Data/Bed/Types.hs +228/−0
• src/Bio/Data/Fasta.hs +6/−5
• src/Bio/GO/GREAT.hs +0/−87
• src/Bio/GO/Parser.hs +2/−2
• src/Bio/Motif.hs +1/−1
• src/Bio/Motif/Search.hs +8/−24
• src/Bio/RealWorld/GENCODE.hs +31/−21
• src/Bio/RealWorld/UCSC.hs +2/−2
• src/Bio/Seq/IO.hs +1/−1
• src/Bio/Utils/Functions.hs +2/−2
• src/Bio/Utils/Overlap.hs +8/−45
• tests/Tests/Bam.hs +3/−3
• tests/Tests/Bed.hs +5/−4
• tests/Tests/GREAT.hs +0/−56
• tests/test.hs +0/−2

@@ -1,5 +1,5 @@ name:                bioinformatics-toolkit-version:             0.4.1+version:             0.5.0 synopsis:            A collection of bioinformatics tools description:         A collection of bioinformatics tools license:             MIT@@ -28,11 +28,11 @@     Bio.ChIPSeq     Bio.ChIPSeq.FragLen     Bio.Data.Bed+    Bio.Data.Bed.Types     Bio.Data.Bam     Bio.Data.Fasta     Bio.GO     Bio.GO.Parser-    Bio.GO.GREAT     Bio.Motif     Bio.Motif.Alignment     Bio.Motif.Merge@@ -67,6 +67,7 @@     , http-conduit >=2.1.8     , hexpat     , IntervalMap >=0.5.0.0+    , lens     , matrices >=0.4.3     , mtl >=2.1.3.1     , math-functions@@ -97,7 +98,6 @@     , data-default-class     , conduit     , mtl-  default-language:    Haskell2010  test-suite tests   type: exitcode-stdio-1.0@@ -109,7 +109,6 @@     , Tests.ChIPSeq     , Tests.Motif     , Tests.Seq-    , Tests.GREAT     , Tests.Tools    default-language:    Haskell2010

@@ -1,6 +1,6 @@-{-# LANGUAGE FlexibleContexts #-}-{-# LANGUAGE OverloadedStrings #-}-{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE BangPatterns          #-}+{-# LANGUAGE FlexibleContexts      #-}+{-# LANGUAGE OverloadedStrings     #-}  module Bio.ChIPSeq     ( monoColonalize@@ -12,48 +12,43 @@     , peakCluster     ) where -import Control.Monad (liftM, forM_, forM)-import Control.Monad.Primitive (PrimMonad)-import Conduit-import Data.Function (on)-import qualified Data.Foldable as F-import qualified Data.HashMap.Strict as M-import qualified Data.IntervalMap as IM-import Data.Maybe (fromJust)-import qualified Data.Vector as V-import qualified Data.Vector.Unboxed as U+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.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  -- | process a sorted BED stream, keep only mono-colonal tags-monoColonalize :: Monad m => Conduit BED m BED+monoColonalize :: Monad m => ConduitT BED BED m () monoColonalize = do-    x <- await-    F.forM_ x loop+    x <- headC+    case x of+        Just b -> yield b >> concatMapAccumC f b+        Nothing -> return ()   where-    loop prev = do-        x <- await-        case x of-            Nothing -> yield prev-            Just current ->-                case () of-                   _ | compareBed prev current == GT ->-                         error $ "Bio.ChIPSeq.monoColonalize: Input is not sorted: " ++ show prev ++ " > " ++ show current-                     | chromStart prev == chromStart current &&-                       chromEnd prev == chromEnd current &&-                       chrom prev == chrom current &&-                       _strand prev == _strand current -> loop prev-                     | otherwise -> yield prev >> loop current+    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] -> Sink BED m (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'@@ -71,21 +66,20 @@ -- | 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) -> Sink BED m (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 . bedSize) (regions V.! i)))+    lift $ liftM (G.imap (\i x -> x / factor / (fromIntegral . size) (regions V.! i)))          $ G.unsafeFreeze vec   where     count v nTags tag = do-        let chr = chrom tag-            p | _strand tag == Just True = chromStart tag-              | _strand tag == Just False = chromEnd tag - 1+        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 chr intervalMap) p+                IM.containing (M.lookupDefault IM.empty (tag^.chrom) intervalMap) p         addOne v xs         return $ succ nTags @@ -102,11 +96,11 @@ countTagsBinBed :: (Integral a, PrimMonad m, G.Vector v a, BEDLike b)            => Int   -- ^ bin size            -> [b]   -- ^ regions-           -> Sink BED m ([v a], Int)+           -> ConduitT BED o m ([v a], Int) countTagsBinBed k beds = do     initRC <- lift $ forM beds $ \bed -> do-        let start = chromStart bed-            end = chromEnd bed+        let start = bed^.chromStart+            end = bed^.chromEnd             num = (end - start) `div` k             index i = (i - start) `div` k         v <- GM.replicate num 0@@ -117,12 +111,12 @@     sink !nTags vs = do         tag <- await         case tag of-            Just (BED chr start end _ _ strand) -> do-                let p | strand == Just True = start-                      | strand == Just False = end - 1+            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 chr intervalMap) p+                        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@@ -142,11 +136,11 @@ countTagsBinBed' :: (Integral a, PrimMonad m, G.Vector v a, BEDLike b1, BEDLike b2)                  => Int   -- ^ bin size                  -> [b1]   -- ^ regions-                 -> Sink b2 m ([v a], Int)+                 -> ConduitT b2 o m ([v a], Int) countTagsBinBed' k beds = do     initRC <- lift $ forM beds $ \bed -> do-        let start = chromStart bed-            end = chromEnd bed+        let start = bed^.chromStart+            end = bed^.chromEnd             num = (end - start) `div` k             index i = (i - start) `div` k         v <- GM.replicate num 0@@ -158,9 +152,9 @@         tag <- await         case tag of             Just bed -> do-                let chr = chrom bed-                    start = chromStart bed-                    end = chromEnd bed+                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@@ -214,18 +208,18 @@ {-# INLINE rpkmBam #-} -} -tagCountDistr :: PrimMonad m => G.Vector v Int => Sink BED m (v Int)+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 chr s e _ _ (Just str)) -> do-                let p | str = s-                      | otherwise = 1 - e-                case M.lookup chr m of-                    Just table -> loop $ M.insert chr (M.insertWith (+) p 1 table) m-                    _ -> loop $ M.insert chr (M.fromList [(p,1)]) m+            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 ->@@ -240,16 +234,15 @@             => [b]   -- ^ peaks             -> Int   -- ^ radius             -> Int   -- ^ cutoff-            -> Source m BED+            -> ConduitT o BED m () peakCluster peaks r th = mergeBedWith mergeFn peaks' .| filterC g   where     peaks' = map f peaks-    f b = let chr = chrom b-              c = (chromStart b + chromEnd b) `div` 2-          in BED3 chr (c-r) (c+r)-    mergeFn xs = BED (chrom $ head xs) lo hi Nothing (Just $ fromIntegral $ length xs) Nothing+    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 (_score b) >= fromIntegral th+        lo = minimum $ map (^.chromStart) xs+        hi = maximum $ map (^.chromEnd) xs+    g b = fromJust (b^.score) >= fromIntegral th {-# INLINE peakCluster #-}

@@ -5,13 +5,14 @@     , naiveCCWithSmooth     ) where -import qualified Data.ByteString.Char8 as B-import qualified Data.HashSet as S-import qualified Data.HashMap.Strict as M-import Data.List (foldl', maximumBy)-import Data.Ord (comparing)-import Bio.Data.Bed-import Control.Parallel.Strategies (parMap, rpar)+import           Bio.Data.Bed+import           Control.Lens                ((^.))+import           Control.Parallel.Strategies (parMap, rpar)+import qualified Data.ByteString.Char8       as B+import qualified Data.HashMap.Strict         as M+import qualified Data.HashSet                as S+import           Data.List                   (foldl', maximumBy)+import           Data.Ord                    (comparing)  -- | estimate fragment length for a ChIP-seq experiment fragLength :: (Int, Int) -> [BED] -> Int@@ -25,10 +26,10 @@ fromBED = map toSet . M.toList . M.fromListWith f . map parseLine     where         parseLine :: BED -> (B.ByteString, ([Int], [Int]))-        parseLine x = case _strand x of-            Just True -> (_chrom x, ([_chromStart x], []))-            Just False -> (_chrom x, ([], [_chromEnd x]))-            _ -> error "Unknown Strand!"+        parseLine x = case x^.strand of+            Just True  -> (x^.chrom, ([x^.chromStart], []))+            Just False -> (x^.chrom, ([], [x^.chromEnd]))+            _          -> error "Unknown Strand!"         f (a,b) (a',b') = (a++a', b++b')         toSet (chr, (forwd, rev)) = (chr, (S.fromList forwd, S.fromList rev)) {-# INLINE fromBED #-}

@@ -11,18 +11,19 @@  import           Bio.Data.Bed import           Bio.HTS-import           Bio.HTS.Types             (Bam, FileHeader (..))+import           Bio.HTS.Types        (Bam, FileHeader (..)) import           Conduit+import           Control.Lens         ((&), (.~)) import           Control.Monad.Reader (ask, lift)  -- | Convert bam record to bed record. Unmapped reads will be discarded.-bamToBed :: Conduit Bam HeaderState BED+bamToBed :: ConduitT Bam BED HeaderState () bamToBed = mapMC bamToBed1 .| concatC {-# INLINE bamToBed #-}  -- | 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 :: Conduit Bam HeaderState (BED, BED)+sortedBamToBedPE :: ConduitT Bam (BED, BED) HeaderState () sortedBamToBedPE = do     maybeBam <- await     case maybeBam of@@ -32,9 +33,8 @@             sortOrd <- getSortOrder <$> lift ask             case sortOrd of                 Queryname -> loopBedPE .| concatC-                _ -> error "Bam file must be sorted by NAME."+                _         -> error "Bam file must be sorted by NAME."   where-    loopBedPE :: Conduit Bam HeaderState (Maybe (BED, BED))     loopBedPE = do         pair <- (,) <$$> await <***> await         case pair of@@ -55,11 +55,13 @@ bamToBed1 :: Bam -> HeaderState (Maybe BED) bamToBed1 bam = do     BamHeader hdr <- lift ask-    return $ (\chr -> BED chr start end nm sc strand) <$> getChr hdr bam+    return $+        (\chr -> asBed chr start end & name .~ nm & score .~ sc & strand .~ str)+        <$> getChr hdr bam   where     start = fromIntegral $ position bam     end = fromIntegral $ endPos bam     nm = Just $ qName bam-    strand = Just $ not $ isRev bam+    str = Just $ not $ isRev bam     sc = Just $ fromIntegral $ mapq bam {-# INLINE bamToBed1 #-}

@@ -1,15 +1,19 @@-{-# LANGUAGE OverloadedStrings #-}-{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleContexts      #-}+{-# LANGUAGE OverloadedStrings     #-}  module Bio.Data.Bed     ( BEDLike(..)--    -- * BED6 format-    , BED(..)-    -- * BED3 format-    , BED3(..)-    -- * NarrowPeak format-    , NarrowPeak(..)+    , BEDConvert(..)+    , BED+    , BED3+    , NarrowPeak+    , npSignal+    , npPvalue+    , npQvalue+    , npPeak+    , BEDExt(..)+    , _bed+    , _data      , BEDTree     , bedToTree@@ -51,220 +55,30 @@     , compareBed     ) where -import Control.Arrow ((***))-import Control.Monad.State.Strict-import qualified Data.ByteString.Char8 as B-import Conduit-import Data.Default.Class (Default(..))-import Data.Function (on)-import qualified Data.Foldable as F-import qualified Data.HashMap.Strict as M-import qualified Data.IntervalMap.Strict as IM-import Data.List (groupBy, sortBy)-import Data.Maybe (fromMaybe, fromJust)-import qualified Data.Vector as V+import           Conduit+import           Control.Arrow                ((***))+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.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 Data.ByteString.Lex.Integral (packDecimal)-import Data.Double.Conversion.ByteString (toShortest)-import Data.Ord (comparing)--import Bio.Motif hiding (_name)-import Bio.Motif.Search-import Bio.Seq-import Bio.Seq.IO-import Bio.Utils.Misc ( binBySizeLeft, binBySize, binBySizeOverlap, bins-                      , readInt, readDouble )---- | A class representing BED-like data, e.g., BED3, BED6 and BED12. BED format--- uses 0-based index (see documentation).-class BEDLike b where-    -- | Construct bed record from chromsomoe, start location and end location-    asBed :: B.ByteString -> Int -> Int -> b--    -- | Convert bytestring to bed format-    fromLine :: B.ByteString -> b--    -- | Convert bed to bytestring-    toLine :: b -> B.ByteString--    -- | Field accessor-    chrom :: b -> B.ByteString-    chromStart :: b -> Int-    chromEnd :: b -> Int-    bedName :: b -> Maybe B.ByteString-    bedScore :: b -> Maybe Double-    bedStrand :: b -> Maybe Bool--    convert :: BEDLike b' => b' -> b-    convert bed = asBed (chrom bed) (chromStart bed) (chromEnd bed)-    {-# INLINE convert #-}--    -- | Return the size of a bed region.-    bedSize :: b -> Int-    bedSize bed = chromEnd bed - chromStart bed--    {-# MINIMAL asBed, fromLine, toLine, chrom, chromStart, chromEnd,-                bedName, bedScore, bedStrand #-}---- * BED6 format---- | BED6 format, as described in http://genome.ucsc.edu/FAQ/FAQformat.html#format1.7-data BED = BED-    { _chrom :: !B.ByteString-    , _chromStart :: {-# UNPACK #-} !Int-    , _chromEnd :: {-# UNPACK #-} !Int-    , _name :: !(Maybe B.ByteString)-    , _score :: !(Maybe Double)-    , _strand :: !(Maybe Bool)  -- ^ True: "+", False: "-"-    } deriving (Eq, Show, Read)--instance Ord BED where-    compare (BED x1 x2 x3 x4 x5 x6) (BED y1 y2 y3 y4 y5 y6) =-        compare (x1,x2,x3,x4,x5,x6) (y1,y2,y3,y4,y5,y6)--instance Default BED where-    def = BED-        { _chrom = ""-        , _chromStart = 0-        , _chromEnd = 0-        , _name = Nothing-        , _score = Nothing-        , _strand = Nothing-        }--instance BEDLike BED where-    asBed chr s e = BED chr s e Nothing Nothing Nothing--    fromLine l = evalState (f (B.split '\t' l)) 1-      where-        f :: [B.ByteString] -> State Int BED-        f [] = do i <- get-                  if i <= 3 then error "Read BED fail: Incorrect number of fields"-                            else return def-        f (x:xs) = do-            i <- get-            put (i+1)-            bed <- f xs-            case i of-                1 -> return $ bed {_chrom = x}-                2 -> return $ bed {_chromStart = readInt x}-                3 -> return $ bed {_chromEnd = readInt x}-                4 -> return $ bed {_name = guard' x}-                5 -> return $ bed {_score = getScore x}-                6 -> return $ bed {_strand = getStrand x}-                _ -> return def--        guard' x | x == "." = Nothing-                 | otherwise = Just x-        getScore x | x == "." = Nothing-                   | otherwise = Just . readDouble $ 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' ]-      where-        strand' | f6 == Just True = "+"-                | f6 == Just False = "-"-                | otherwise = "."-        score' = case f5 of-                     Just x -> (B.pack.show) x-                     _ -> "."-    {-# INLINE toLine #-}--    chrom = _chrom-    chromStart = _chromStart-    chromEnd = _chromEnd-    bedName = _name-    bedScore = _score-    bedStrand = _strand--    convert bed = BED (chrom bed) (chromStart bed) (chromEnd bed) (bedName bed)-                      (bedScore bed) (bedStrand bed)---- * BED3 format--data BED3 = BED3 !B.ByteString !Int !Int deriving (Eq, Show, Read)--instance Ord BED3 where-    compare (BED3 x1 x2 x3) (BED3 y1 y2 y3) = compare (x1,x2,x3) (y1,y2,y3)--instance BEDLike BED3 where-    asBed = BED3--    fromLine l = case B.split '\t' l of-                    (a:b:c:_) -> BED3 a (readInt b) $ readInt c-                    _ -> error "Read BED fail: Incorrect number of fields"-    {-# INLINE fromLine #-}--    toLine (BED3 a b c) = B.intercalate "\t"-        [a, fromJust $ packDecimal b, fromJust $ packDecimal c]-    {-# INLINE toLine #-}--    chrom (BED3 f1 _ _) = f1-    chromStart (BED3 _ f2 _) = f2-    chromEnd (BED3 _ _ f3) = f3-    bedName = const Nothing-    bedScore = const Nothing-    bedStrand = const Nothing---- | ENCODE narrowPeak format: https://genome.ucsc.edu/FAQ/FAQformat.html#format12-data NarrowPeak = NarrowPeak-    { _npChrom :: !B.ByteString-    , _npStart :: !Int-    , _npEnd :: !Int-    , _npName :: !(Maybe B.ByteString)-    , _npScore :: !Double-    , _npStrand :: !(Maybe Bool)-    , _npSigal :: !Double-    , _npPvalue :: !(Maybe Double)-    , _npQvalue :: !(Maybe Double)-    , _npPeak :: !(Maybe Int)-    } deriving (Eq, Show, Read)--instance BEDLike NarrowPeak where-    asBed chr s e = NarrowPeak chr s e Nothing 0 Nothing 0 Nothing Nothing Nothing--    fromLine l = NarrowPeak a (readInt b) (readInt c)-        (if d == "." then Nothing else Just d)-        (readDouble e)-        (if f == "." then Nothing else if f == "+" then Just True else Just False)-        (readDouble g)-        (if readDouble h < 0 then Nothing else Just $ readDouble h)-        (if readDouble i < 0 then Nothing else Just $ readDouble i)-        (if readInt j < 0 then Nothing else Just $ readInt j)-      where-        (a:b:c:d:e:f:g:h:i:j:_) = B.split '\t' l-    {-# INLINE fromLine #-}--    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-        , case f of-            Nothing -> "."-            Just True -> "+"-            _ -> "-"-        , toShortest g, fromMaybe "-1" $ fmap toShortest h-        , fromMaybe "-1" $ fmap toShortest i-        , fromMaybe "-1" $ fmap (fromJust . packDecimal) j-        ]-    {-# INLINE toLine #-}--    chrom = _npChrom-    chromStart = _npStart-    chromEnd = _npEnd-    bedName = _npName-    bedScore = Just . _npScore-    bedStrand = _npStrand--    convert bed = NarrowPeak (chrom bed) (chromStart bed) (chromEnd bed) (bedName bed)-        (fromMaybe 0 $ bedScore bed) (bedStrand bed) 0 Nothing Nothing Nothing+import           System.IO -type BEDTree a = M.HashMap B.ByteString (IM.IntervalMap Int a)+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)  -- | Convert a set of sorted bed records to interval tree, with combining -- function for equal keys.@@ -272,38 +86,28 @@                 => (a -> a -> a)                 -> Sorted (f (b, a))                 -> BEDTree a-sortedBedToTree f (Sorted xs) =-      M.fromList-    . map ((head *** IM.fromAscListWith f) . unzip)-    . groupBy ((==) `on` fst)-    . map (\(bed, x) -> (chrom bed, (IM.IntervalCO (chromStart bed) (chromEnd bed), x)))-    . F.toList-    $ xs+sortedBedToTree f (Sorted xs) = M.fromList $+    map ((head *** IM.fromAscListWith f) . unzip) $ groupBy ((==) `on` fst) $+    map (\(b, x) -> (b^.chrom, (IM.IntervalCO (b^.chromStart) (b^.chromEnd), x))) $+    F.toList xs {-# INLINE sortedBedToTree #-}  bedToTree :: BEDLike b           => (a -> a -> a)           -> [(b, a)]           -> BEDTree a-bedToTree f xs =-      M.fromList-    . map ((head *** IM.fromAscListWith f) . unzip)-    . groupBy ((==) `on` fst)-    . map (\(bed, x) -> (chrom bed, (IM.IntervalCO (chromStart bed) (chromEnd bed), x)))-    . V.toList-    $ xs'-  where-    xs' = V.create $ do-        v <- V.unsafeThaw . V.fromList $ xs+bedToTree f xs = M.fromList $ map ((head *** IM.fromAscListWith f) . unzip) $+    groupBy ((==) `on` fst) $+    map (\(b, x) -> (b^.chrom, (IM.IntervalCO (b^.chromStart) (b^.chromEnd), x))) $+    V.toList $ V.create $ do+        v <- V.unsafeThaw $ V.fromList xs         I.sortBy (compareBed `on` fst) v         return v {-# INLINE bedToTree #-}  intersecting :: BEDLike b => BEDTree a -> b -> IM.IntervalMap Int a-intersecting tree x = IM.intersecting (M.lookupDefault IM.empty chr tree) interval-  where-    chr = chrom x-    interval = IM.IntervalCO (chromStart x) $ chromEnd x+intersecting tree x = IM.intersecting (M.lookupDefault IM.empty (x^.chrom) tree) $+    IM.IntervalCO (x^.chromStart) $ x^.chromEnd {-# INLINE intersecting #-}  isIntersected :: BEDLike b => BEDTree a -> b -> Bool@@ -311,56 +115,39 @@ {-# INLINE isIntersected #-}  sizeOverlapped :: (BEDLike b1, BEDLike b2) => b1 -> b2 -> Int-sizeOverlapped x y | chr1 /= chr2 = 0-                   | overlap < 0 = 0-                   | otherwise = overlap+sizeOverlapped b1 b2 | b1^.chrom /= b2^.chrom = 0+                     | overlap < 0 = 0+                     | otherwise = overlap   where-    overlap = minimum [e1 - s2, e2 - s1, e1 - s1, e2 - s2]-    chr1 = chrom x-    s1 = chromStart x-    e1 = chromEnd x-    chr2 = chrom y-    s2 = chromStart y-    e2 = chromEnd y-+    overlap = minimum [ b1^.chromEnd - b2^.chromStart+                      , b2^.chromEnd - b1^.chromStart+                      , b1^.chromEnd - b1^.chromStart+                      , b2^.chromEnd - b2^.chromStart ]  -- | split a bed region into k consecutive subregions, discarding leftovers-splitBed :: BEDLike b => Int -> b -> [b]-splitBed k bed = map (uncurry (asBed chr)) . bins k $ (s, e)-  where-    chr = chrom bed-    s = chromStart bed-    e = chromEnd bed+splitBed :: BEDConvert b => Int -> b -> [b]+splitBed k bed = map (uncurry (asBed (bed^.chrom))) $+    bins k (bed^.chromStart, bed^.chromEnd) {-# INLINE splitBed #-}  -- | split a bed region into consecutive fixed size subregions, discarding leftovers-splitBedBySize :: BEDLike b => Int -> b -> [b]-splitBedBySize k bed = map (uncurry (asBed chr)) . binBySize k $ (s, e)-  where-    chr = chrom bed-    s = chromStart bed-    e = chromEnd bed+splitBedBySize :: BEDConvert b => Int -> b -> [b]+splitBedBySize k bed = map (uncurry (asBed (bed^.chrom))) $+    binBySize k (bed^.chromStart, bed^.chromEnd) {-# INLINE splitBedBySize #-}  -- | split a bed region into consecutive fixed size subregions, including leftovers-splitBedBySizeLeft :: BEDLike b => Int -> b -> [b]-splitBedBySizeLeft k bed = map (uncurry (asBed chr)) . binBySizeLeft k $ (s, e)-  where-    chr = chrom bed-    s = chromStart bed-    e = chromEnd bed+splitBedBySizeLeft :: BEDConvert b => Int -> b -> [b]+splitBedBySizeLeft k bed = map (uncurry (asBed (bed^.chrom))) $+    binBySizeLeft k (bed^.chromStart, bed^.chromEnd) {-# INLINE splitBedBySizeLeft #-} -splitBedBySizeOverlap :: BEDLike b+splitBedBySizeOverlap :: BEDConvert b                       => Int     -- ^ bin size                       -> Int     -- ^ overlap size                       -> b -> [b]-splitBedBySizeOverlap k o bed = map (uncurry (asBed chr)) .-    binBySizeOverlap k o $ (s, e)-  where-    chr = chrom bed-    s = chromStart bed-    e = chromEnd bed+splitBedBySizeOverlap k o bed = map (uncurry (asBed (bed^.chrom))) $+    binBySizeOverlap k o (bed^.chromStart, bed^.chromEnd) {-# INLINE splitBedBySizeOverlap #-}  -- | a type to imply that underlying data structure is sorted@@ -370,10 +157,8 @@ -- Unlike the ``compare'' from the Ord type class, this function can compare -- different types of BED data types. compareBed :: (BEDLike b1, BEDLike b2) => b1 -> b2 -> Ordering-compareBed x y = compare x' y'-  where-    x' = (chrom x, chromStart x, chromEnd x)-    y' = (chrom y, chromStart y, chromEnd y)+compareBed b1 b2 = compare (b1^.chrom, b1^.chromStart, b1^.chromEnd)+                           (b2^.chrom, b2^.chromStart, b2^.chromEnd) {-# INLINE compareBed #-}  -- | sort BED, first by chromosome (alphabetical order), then by chromStart, last by chromEnd@@ -385,7 +170,7 @@ {-# INLINE sortBed #-}  -- | return records in A that are overlapped with records in B-intersectBed :: (BEDLike b1, BEDLike b2, Monad m) => [b2] -> Conduit b1 m b1+intersectBed :: (BEDLike b1, BEDLike b2, Monad m) => [b2] -> ConduitT b1 b1 m () intersectBed b = intersectSortedBed b'   where     b' = sortBed b@@ -393,59 +178,55 @@  -- | return records in A that are overlapped with records in B intersectSortedBed :: (BEDLike b1, BEDLike b2, Monad m)-                   => Sorted (V.Vector b2) -> Conduit b1 m b1+                   => Sorted (V.Vector b2) -> ConduitT b1 b1 m () intersectSortedBed (Sorted b) = filterC (not . IM.null . intersecting tree)   where     tree = sortedBedToTree (\_ _ -> ()) . Sorted $ V.map (\x -> (x,())) b {-# INLINE intersectSortedBed #-}  intersectBedWith :: (BEDLike b1, BEDLike b2, Monad m)-                 => ([b2] -> a)+                 => (b1 -> [b2] -> a)                  -> [b2]-                 -> Conduit b1 m (b1, a)+                 -> ConduitT b1 a m () intersectBedWith fn = intersectSortedBedWith fn . sortBed {-# INLINE intersectBedWith #-}  intersectSortedBedWith :: (BEDLike b1, BEDLike b2, Monad m)-                       => ([b2] -> a)+                       => (b1 -> [b2] -> a)                        -> Sorted (V.Vector b2)-                       -> Conduit b1 m (b1, a)-intersectSortedBedWith fn (Sorted b) = mapC f+                       -> ConduitT b1 a m ()+intersectSortedBedWith fn (Sorted b) = mapC $ \input -> fn input+    $ concat $ IM.elems $ intersecting tree input   where-    f bed = (bed, fn $ concat $ IM.elems $ intersecting tree bed)     tree = sortedBedToTree (++) $ Sorted $ V.map (\x -> (x, [x])) b {-# INLINE intersectSortedBedWith #-}  isOverlapped :: (BEDLike b1, BEDLike b2) => b1 -> b2 -> Bool-isOverlapped bed1 bed2 = chr == chr' && not (e <= s' || e' <= s)-  where-    chr = chrom bed1-    s = chromStart bed1-    e = chromEnd bed1-    chr' = chrom bed2-    s' = chromStart bed2-    e' = chromEnd bed2+isOverlapped b1 b2 = b1^.chrom == b2^.chrom &&+    not (b1^.chromEnd <= b2^.chromStart || b2^.chromEnd <= b1^.chromStart) -mergeBed :: (BEDLike b, Monad m) => [b] -> Source m b+mergeBed :: (BEDConvert b, Monad m) => [b] -> ConduitT i b m () mergeBed = mergeSortedBed . sortBed {-# INLINE mergeBed #-}  mergeBedWith :: (BEDLike b, Monad m)-             => ([b] -> a) -> [b] -> Source m a+             => ([b] -> a) -> [b] -> ConduitT i a m () mergeBedWith f = mergeSortedBedWith f . sortBed {-# INLINE mergeBedWith #-} -mergeSortedBed :: (BEDLike b, Monad m) => Sorted (V.Vector b) -> Source m b+mergeSortedBed :: (BEDConvert b, Monad m)+               => Sorted (V.Vector b)+               -> ConduitT i b m () mergeSortedBed = mergeSortedBedWith f   where-    f xs = asBed (chrom $ head xs) lo hi+    f xs = asBed (head xs ^. chrom) lo hi       where-        lo = minimum . map chromStart $ xs-        hi = maximum . map chromEnd $ xs+        lo = minimum $ map (^.chromStart) xs+        hi = maximum $ map (^.chromEnd) xs {-# INLINE mergeSortedBed #-}  mergeSortedBedWith :: (BEDLike b, Monad m)-                   => ([b] -> a) -> Sorted (V.Vector b) -> Source m a+                   => ([b] -> a) -> Sorted (V.Vector b) -> ConduitT i a m () mergeSortedBedWith mergeFn (Sorted beds)     | V.null beds = return ()     | otherwise = do@@ -453,42 +234,42 @@         yield $ mergeFn r   where     x0 = V.head beds-    acc0 = ((chrom x0, chromStart x0, chromEnd x0), [x0])+    acc0 = ((x0^.chrom, x0^.chromStart, x0^.chromEnd), [x0])     f ((chr,lo,hi), acc) bed         | chr /= chr' || s' > hi = yield (mergeFn acc) >>                                    return ((chr',s',e'), [bed])         | e' > hi = return ((chr',lo,e'), bed:acc)         | otherwise = return ((chr,lo,hi), bed:acc)       where-        chr' = chrom bed-        s' = chromStart bed-        e' = chromEnd bed+        chr' = bed^.chrom+        s' = bed^.chromStart+        e' = bed^.chromEnd {-# INLINE mergeSortedBedWith #-}  -- | Split overlapped regions into non-overlapped regions. The input must be overlapped. -- This function is usually used with `mergeBedWith`. splitOverlapped :: BEDLike b => ([b] -> a) -> [b] -> [(BED3, a)]-splitOverlapped fun xs = filter ((>0) . bedSize . fst) $+splitOverlapped fun xs = filter ((>0) . size . fst) $     evalState (F.foldrM f [] $ init xs') x0   where-    x0 = (\(a,b) -> (fromEither a, M.singleton (chromStart b, chromEnd b) b)) $ last xs'+    x0 = (\(a,b) -> (fromEither a, M.singleton (b^.chromStart, b^.chromEnd) b)) $ last xs'     xs' = sortBy (comparing (fromEither . fst)) $ concatMap-        ( \x -> [(Left $ chromStart x, x), (Right $ chromEnd x, x)] ) xs+        ( \x -> [(Left $ x^.chromStart, x), (Right $ x^.chromEnd, x)] ) xs     f (i, x) acc = do         (j, set) <- get-        let bed = (BED3 chr (fromEither i) j, fun $ M.elems set)+        let bed = (asBed chr (fromEither i) j, fun $ M.elems set)             set' = case i of-                Left _ -> M.delete (chromStart x, chromEnd x) set-                Right _ -> M.insert (chromStart x, chromEnd x) x set+                Left _  -> M.delete (x^.chromStart, x^.chromEnd) set+                Right _ -> M.insert (x^.chromStart, x^.chromEnd) x set         put (fromEither i, set')         return (bed:acc)-    fromEither (Left x) = x+    fromEither (Left x)  = x     fromEither (Right x) = x-    chr = chrom $ head xs+    chr = head xs ^. chrom {-# INLINE splitOverlapped #-}  -- | Read records from a bed file handler in a streaming fashion.-hReadBed :: (BEDLike b, MonadIO m) => Handle -> Source m b+hReadBed :: (BEDConvert b, MonadIO m) => Handle -> ConduitT i b m () hReadBed h = do eof <- liftIO $ hIsEOF h                 unless eof $ do                     line <- liftIO $ B.hGetLine h@@ -497,23 +278,23 @@ {-# INLINE hReadBed #-}  -- | Non-streaming version.-hReadBed' :: (BEDLike b, MonadIO m) => Handle -> m [b]+hReadBed' :: (BEDConvert b, MonadIO m) => Handle -> m [b] hReadBed' h = runConduit $ hReadBed h .| sinkList {-# INLINE hReadBed' #-}  -- | Read records from a bed file in a streaming fashion.-readBed :: (BEDLike b, MonadIO m) => FilePath -> Source m b+readBed :: (BEDConvert b, MonadIO m) => FilePath -> ConduitT i b m () readBed fl = do handle <- liftIO $ openFile fl ReadMode                 hReadBed handle                 liftIO $ hClose handle {-# INLINE readBed #-}  -- | Non-streaming version.-readBed' :: (BEDLike b, MonadIO m) => FilePath -> m [b]+readBed' :: (BEDConvert b, MonadIO m) => FilePath -> m [b] readBed' fl = runConduit $ readBed fl .| sinkList {-# INLINE readBed' #-} -hWriteBed :: (BEDLike b, MonadIO m) => Handle -> Sink b m ()+hWriteBed :: (BEDConvert b, MonadIO m) => Handle -> ConduitT b o m () hWriteBed handle = do     x <- await     case x of@@ -521,29 +302,31 @@         Just bed -> (liftIO . B.hPutStrLn handle . toLine) bed >> hWriteBed handle {-# INLINE hWriteBed #-} -hWriteBed' :: (BEDLike b, MonadIO m) => Handle -> [b] -> m ()+hWriteBed' :: (BEDConvert b, MonadIO m) => Handle -> [b] -> m () hWriteBed' handle beds = runConduit $ yieldMany beds .| hWriteBed handle {-# INLINE hWriteBed' #-} -writeBed :: (BEDLike b, MonadIO m) => FilePath -> Sink b m ()+writeBed :: (BEDConvert b, MonadIO m) => FilePath -> ConduitT b o m () writeBed fl = do handle <- liftIO $ openFile fl WriteMode                  hWriteBed handle                  liftIO $ hClose handle {-# INLINE writeBed #-} -writeBed' :: (BEDLike b, MonadIO m) => FilePath -> [b] -> m ()-writeBed' fl beds = yieldMany beds $$ writeBed fl+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 -> Conduit BED m (Either String (DNA a))+fetchSeq :: (BioSeq DNA a, MonadIO m)+         => Genome+         -> ConduitT BED (Either String (DNA a)) m () fetchSeq g = mapMC f   where-    f (BED chr start end _ _ isForward) = do-        dna <- liftIO $ getSeq g (chr, start, end)-        return $ case isForward of+    f bed = do+        dna <- liftIO $ getSeq g (bed^.chrom, bed^.chromStart, bed^.chromEnd)+        return $ case bed^.strand of             Just False -> rc <$> dna-            _ -> dna+            _          -> dna {-# INLINE fetchSeq #-}  fetchSeq' :: (BioSeq DNA a, MonadIO m) => Genome -> [BED] -> m [Either String (DNA a)]@@ -552,45 +335,43 @@  -- | Identify motif binding sites motifScan :: (BEDLike b, MonadIO m)-          => Genome -> [Motif] -> Bkgd -> Double -> Conduit b m BED+          => Genome -> [Motif] -> Bkgd -> Double -> ConduitT b BED m () motifScan g motifs bg p = awaitForever $ \bed -> do-    let chr = chrom bed-        s = chromStart bed-        e = chromEnd bed-    r <- liftIO $ getSeq g (chr, s, e)+    r <- liftIO $ getSeq g (bed^.chrom, bed^.chromStart, bed^.chromEnd)     case r of-        Left _ -> return ()-        Right dna -> mapM_ (getTFBS dna (chr, s)) motifs'+        Left _    -> return ()+        Right dna -> mapM_ (getTFBS dna (bed^.chrom, bed^.chromStart)) motifs'   where     getTFBS dna (chr, s) (nm, (pwm, cutoff), (pwm', cutoff')) = toProducer-        ( (findTFBS bg pwm (dna :: DNA IUPAC) cutoff True =$=-            mapC (\i -> BED chr (s+i) (s+i+n) (Just nm) Nothing $ Just True)) >>-          (findTFBS bg pwm' dna cutoff' True =$=-            mapC (\i -> BED chr (s+i) (s+i+n) (Just nm) Nothing $ Just False)) )+        ( (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 = Bio.Motif.size pwm+        n = Motif.size pwm+        bed = asBed chr s (s+n) & name .~ Just nm     motifs' = flip map motifs $ \(Motif nm pwm) ->-        let cutoff = pValueToScore p bg pwm-            cutoff' = pValueToScore p bg pwm'-            pwm' = rcPWM 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 -> Conduit BED m BED-getMotifScore g motifs bg = awaitForever $ \(BED chr s e (Just nm) _ isForward) -> do-    r <- liftIO $ getSeq g (chr, s, e)-    let r' = case isForward of+              => 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+            _          -> r     case r' of         Left _ -> return ()         Right dna -> do             let pwm = M.lookupDefault (error "can't find motif with given name")-                      nm motifMap-                sc = score bg pwm (dna :: DNA IUPAC)-            yield $ BED chr s e (Just nm) (Just sc) isForward+                        (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@@ -600,18 +381,18 @@                => 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+               -> [Motif] -> Bkgd -> ConduitT BED BED m () getMotifPValue truncation motifs bg = mapC $ \bed ->-    let nm = fromJust $ bedName bed-        sc = fromJust $ bedScore 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 - cdf d sc-     in bed{_score = Just p}+        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 $ scoreCDF bg pwm)) motifs+        map (\(Motif nm pwm) -> (nm, compressCDF $ Motif.scoreCDF bg pwm)) motifs     compressCDF = case truncation of         Nothing -> id-        Just x -> truncateCDF x+        Just x  -> Motif.truncateCDF x {-# INLINE getMotifPValue #-}

@@ -0,0 +1,228 @@+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE TemplateHaskell #-}++module Bio.Data.Bed.Types where++import           Control.Lens+import qualified Data.ByteString.Char8             as B+import           Data.ByteString.Lex.Integral      (packDecimal)+import           Data.Default.Class                (Default (..))+import           Data.Double.Conversion.ByteString (toShortest)+import qualified Data.HashMap.Strict               as M+import qualified Data.IntervalMap.Strict           as IM+import           Data.Maybe                        (fromJust, fromMaybe)+import           Data.Monoid                       ((<>))++import           Bio.Utils.Misc                    (readDouble, readInt)++-- | A class representing BED-like data, e.g., BED3, BED6 and BED12. BED format+-- uses 0-based index (see documentation).+class BEDLike b where+    -- | Field lens+    chrom :: Lens' b B.ByteString+    chromStart :: Lens' b Int+    chromEnd :: Lens' b Int+    name :: Lens' b (Maybe B.ByteString)+    score :: Lens' b (Maybe Double)+    strand :: Lens' b (Maybe Bool)++    -- | Return the size of a bed region.+    size :: b -> Int+    size bed = bed^.chromEnd - bed^.chromStart+    {-# INLINE size #-}++    {-# MINIMAL chrom, chromStart, chromEnd, name, score, strand #-}++class BEDLike b => BEDConvert b where+    -- | Construct bed record from chromsomoe, start location and end location+    asBed :: B.ByteString -> Int -> Int -> b++    -- | Convert bytestring to bed format+    fromLine :: B.ByteString -> b++    -- | Convert bed to bytestring+    toLine :: b -> B.ByteString++    convert :: BEDLike b' => b' -> b+    convert bed = asBed (bed^.chrom) (bed^.chromStart) (bed^.chromEnd)+    {-# INLINE convert #-}++    {-# MINIMAL asBed, fromLine, toLine #-}++-- * BED6 format++-- | BED6 format, as described in http://genome.ucsc.edu/FAQ/FAQformat.html#format1.7+data BED = BED+    { _bed_chrom      :: !B.ByteString+    , _bed_chromStart :: !Int+    , _bed_chromEnd   :: !Int+    , _bed_name       :: !(Maybe B.ByteString)+    , _bed_score      :: !(Maybe Double)+    , _bed_strand     :: !(Maybe Bool)  -- ^ True: "+", False: "-"+    } deriving (Eq, Show, Read)++instance Ord BED where+    compare (BED x1 x2 x3 x4 x5 x6) (BED y1 y2 y3 y4 y5 y6) =+        compare (x1,x2,x3,x4,x5,x6) (y1,y2,y3,y4,y5,y6)++instance BEDLike BED where+    chrom = lens _bed_chrom (\bed x -> bed { _bed_chrom = x })+    chromStart = lens _bed_chromStart (\bed x -> bed { _bed_chromStart = x })+    chromEnd = lens _bed_chromEnd (\bed x -> bed { _bed_chromEnd = x })+    name = lens _bed_name (\bed x -> bed { _bed_name = x })+    score = lens _bed_score (\bed x -> bed { _bed_score = x })+    strand = lens _bed_strand (\bed x -> bed { _bed_strand = x })++instance BEDConvert BED where+    asBed chr s e = BED chr s e Nothing Nothing Nothing++    fromLine l = f $ take 6 $ B.split '\t' l+      where+        f [f1,f2,f3,f4,f5,f6] = BED f1 (readInt f2) (readInt f3) (getName f4)+            (getScore f5) (getStrand f6)+        f [f1,f2,f3,f4,f5] = BED f1 (readInt f2) (readInt f3) (getName f4)+            (getScore f5) Nothing+        f [f1,f2,f3,f4] = BED f1 (readInt f2) (readInt f3) (getName f4)+            Nothing Nothing+        f [f1,f2,f3] = asBed f1 (readInt f2) (readInt f3)+        f _ = error "Read BED fail: Not enough fields!"+        getName x | x == "." = Nothing+                  | otherwise = Just x+        getScore x | x == "." = Nothing+                   | otherwise = Just . readDouble $ 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' ]+      where+        strand' | f6 == Just True = "+"+                | f6 == Just False = "-"+                | otherwise = "."+        score' = case f5 of+                     Just x -> (B.pack.show) x+                     _      -> "."+    {-# INLINE toLine #-}++    convert bed = BED (bed^.chrom) (bed^.chromStart) (bed^.chromEnd) (bed^.name)+                      (bed^.score) (bed^.strand)++-- * BED3 format++data BED3 = BED3+    { _bed3_chrom       :: !B.ByteString+    , _bed3_chrom_start :: !Int+    , _bed3_chrom_end   :: !Int+    } deriving (Eq, Show, Read)++instance Ord BED3 where+    compare (BED3 x1 x2 x3) (BED3 y1 y2 y3) = compare (x1,x2,x3) (y1,y2,y3)++instance BEDLike BED3 where+    chrom = lens _bed3_chrom (\bed x -> bed { _bed3_chrom = x })+    chromStart = lens _bed3_chrom_start (\bed x -> bed { _bed3_chrom_start = x })+    chromEnd = lens _bed3_chrom_end (\bed x -> bed { _bed3_chrom_end = x })+    name = lens (const Nothing) (\bed _ -> bed)+    score = lens (const Nothing) (\bed _ -> bed)+    strand = lens (const Nothing) (\bed _ -> bed)++instance BEDConvert BED3 where+    asBed = BED3++    fromLine l = case B.split '\t' l of+                    (a:b:c:_) -> BED3 a (readInt b) $ readInt c+                    _ -> error "Read BED fail: Incorrect number of fields"+    {-# INLINE fromLine #-}++    toLine (BED3 a b c) = B.intercalate "\t"+        [a, fromJust $ packDecimal b, fromJust $ packDecimal c]+    {-# INLINE toLine #-}++-- | ENCODE narrowPeak format: https://genome.ucsc.edu/FAQ/FAQformat.html#format12+data NarrowPeak = NarrowPeak+    { _npChrom  :: !B.ByteString+    , _npStart  :: !Int+    , _npEnd    :: !Int+    , _npName   :: !(Maybe B.ByteString)+    , _npScore  :: !Double+    , _npStrand :: !(Maybe Bool)+    , _npSignal  :: !Double+    , _npPvalue :: !(Maybe Double)+    , _npQvalue :: !(Maybe Double)+    , _npPeak   :: !(Maybe Int)+    } deriving (Eq, Show, Read)++makeLensesFor [ ("_npSignal", "npSignal")+              , ("_npPvalue", "npPvalue")+              , ("_npQvalue", "npQvalue")+              , ("_npPeak", "npPeak")+              ] ''NarrowPeak++instance BEDLike NarrowPeak where+    chrom = lens _npChrom (\bed x -> bed { _npChrom = x })+    chromStart = lens _npStart (\bed x -> bed { _npStart = x })+    chromEnd = lens _npEnd (\bed x -> bed { _npEnd = x })+    name = lens _npName (\bed x -> bed { _npName = x })+    score = lens (Just . _npScore) (\bed x -> bed { _npScore = fromJust x })+    strand = lens _npStrand (\bed x -> bed { _npStrand = x })++instance BEDConvert NarrowPeak where+    asBed chr s e = NarrowPeak chr s e Nothing 0 Nothing 0 Nothing Nothing Nothing++    fromLine l = NarrowPeak a (readInt b) (readInt c)+        (if d == "." then Nothing else Just d)+        (readDouble e)+        (if f == "." then Nothing else if f == "+" then Just True else Just False)+        (readDouble g)+        (if readDouble h < 0 then Nothing else Just $ readDouble h)+        (if readDouble i < 0 then Nothing else Just $ readDouble i)+        (if readInt j < 0 then Nothing else Just $ readInt j)+      where+        (a:b:c:d:e:f:g:h:i:j:_) = B.split '\t' l+    {-# INLINE fromLine #-}++    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+        , case f of+            Nothing   -> "."+            Just True -> "+"+            _         -> "-"+        , toShortest g, fromMaybe "-1" $ fmap toShortest h+        , fromMaybe "-1" $ fmap toShortest i+        , fromMaybe "-1" $ fmap (fromJust . packDecimal) j+        ]+    {-# INLINE toLine #-}++    convert bed = NarrowPeak (bed^.chrom) (bed^.chromStart) (bed^.chromEnd) (bed^.name)+        (fromMaybe 0 $ bed^.score) (bed^.strand) 0 Nothing Nothing Nothing++data BEDExt bed a = BEDExt+    { _ext_bed :: bed+    , _ext_data :: a+    } deriving (Eq, Show, Read)++makeLensesFor [("_ext_bed", "_bed"), ("_ext_data", "_data")] ''BEDExt++instance BEDLike bed => BEDLike (BEDExt bed a) where+    chrom = _bed . chrom+    chromStart = _bed . chromStart+    chromEnd = _bed . chromEnd+    name = _bed . name+    score = _bed . score+    strand = _bed . strand++instance (Default a, Read a, Show a, BEDConvert bed) => BEDConvert (BEDExt bed a) where+    asBed chr s e = BEDExt (asBed chr s e) def++    fromLine l = let (a, b) = B.breakEnd (=='\t') l+                 in BEDExt (fromLine $ B.init a) $ read $ B.unpack b+    {-# INLINE fromLine #-}++    toLine (BEDExt bed a) = toLine bed <> "\t" <> B.pack (show a)+    {-# INLINE toLine #-}++type BEDTree a = M.HashMap B.ByteString (IM.IntervalMap Int a)

@@ -16,12 +16,12 @@     -- to a specific data type     fromFastaRecord :: (B.ByteString, [B.ByteString]) -> f -    readFasta :: FilePath -> Source (ResourceT IO) f-    readFasta fl = fastaReader fl =$= mapC fromFastaRecord+    readFasta :: FilePath -> ConduitT i f (ResourceT IO) ()+    readFasta fl = fastaReader fl .| mapC fromFastaRecord      -- | non-stream version, read whole file in memory     readFasta' :: FilePath -> IO [f]-    readFasta' fl = runResourceT $ readFasta fl $$ sinkList+    readFasta' fl = runResourceT $ runConduit $ readFasta fl .| sinkList     {-# MINIMAL fromFastaRecord #-}  instance BioSeq s a => FastaLike (s a) where@@ -34,8 +34,9 @@     fromFastaRecord (name, mat) = Motif name (toPWM mat)     {-# INLINE fromFastaRecord #-} -fastaReader :: FilePath -> Source (ResourceT IO) (B.ByteString, [B.ByteString])-fastaReader fl = sourceFile fl =$= linesUnboundedAsciiC =$= loop []+fastaReader :: FilePath+            -> ConduitT i (B.ByteString, [B.ByteString]) (ResourceT IO) ()+fastaReader fl = sourceFile fl .| linesUnboundedAsciiC .| loop []   where     loop acc = do         x <- await

@@ -1,87 +0,0 @@-{-# LANGUAGE BangPatterns #-}-module Bio.GO.GREAT-    ( AssocRule(..)-    , defaultAssocRule-    , getRegulatoryDomains-    , read3DContact-    , get3DRegulatoryDomains-    ) where--import           Conduit-import           Control.Monad         (forM_)-import qualified Data.ByteString.Char8 as B-import           Data.Function         (on)-import qualified Data.IntervalMap      as IM-import           Data.List             (foldl', sortBy)-import           Data.List.Ordered     (nubSort)-import           Data.Maybe            (fromJust, isNothing)--import           Bio.Data.Bed-import           Bio.Utils.Misc        (readInt)---- | How to associate genomic regions with genes-data AssocRule =-    BasalPlusExtension Int Int Int  -- ^ Upstream, downstream and extension-  | TwoNearest                      -- ^ Not implemented-  | OneNearest                      -- ^ Not implemented---- | The default rule is defined as -5k ~ +1k around TSSs,--- plus up to 1000k extension.-defaultAssocRule :: AssocRule-defaultAssocRule = BasalPlusExtension 5000 1000 1000000---- | A Gene consists of the chromosome name, TSS, strandness and an associated value.-type Gene a = ((B.ByteString, Int, Bool), a)---- | Given a gene list and the rule, compute the rulatory domain for each gene-getRegulatoryDomains :: AssocRule -> [Gene a] -> [(BED3, a)]-getRegulatoryDomains _ [] = error "No gene available for domain assignment!"-getRegulatoryDomains (BasalPlusExtension up dw ext) genes = zip-    (loop $ [Nothing] ++ map Just basal ++ [Nothing]) names-  where-    loop (a:b:c:rest) = fn a b c : loop (b:c:rest)-    loop _            = []-    fn left (Just (BED3 chr s e)) right = BED3 chr leftPos rightPos-      where-        leftPos-            | isNothing left || chr /= chrom (fromJust left) = max (s - ext) 0-            | otherwise = min s $ max (s - ext) $ chromEnd $ fromJust left-        rightPos-            | isNothing right || chr /= chrom (fromJust right) = e + ext   -- TODO: bound check-            | otherwise = max e $ min (e + ext) $ chromStart $ fromJust right-    (basal, names) = unzip $ sortBy (compareBed `on` fst) $ flip map genes $-        \((chr, tss, str), x) -> if str then (BED3 chr (tss - up) (tss + dw), x)-                                        else (BED3 chr (tss - dw) (tss + up), x)-getRegulatoryDomains _ _ = undefined-{-# INLINE getRegulatoryDomains #-}---- | Read 3D contacts from a file, where each line contains 6 fields separated--- by Tabs corresponding to 2 interacting loci. Example:--- chr1 [TAB] 11 [TAB] 100 [TAB] chr2 [TAB] 23 [TAB] 200-read3DContact :: FilePath -> Source (ResourceT IO) (BED3, BED3)-read3DContact input = sourceFileBS input =$= linesUnboundedAsciiC =$= mapC f-  where-    f x = let (chr1:s1:e1:chr2:s2:e2:_) = B.split '\t' x-          in (BED3 chr1 (readInt s1) (readInt e1), BED3 chr2 (readInt s2) (readInt e2))---- | 3D regulatory domains of a gene include the gene's promoter and regions--- that contact with the gene's promoter.-get3DRegulatoryDomains :: (Monad m, Ord a)-                       => [Gene a]   -- ^ Genes-                       -> Int        -- ^ Upstream-                       -> Int        -- ^ Downstream-                       -> Conduit (BED3, BED3) m (BED3, a)-get3DRegulatoryDomains genes up dw = concatRegions >> promoters-  where-    promoters = forM_ genes $ \((chr, tss, str), x) ->-        if str then yield (BED3 chr (gateZero $ tss - up) (tss + dw), x)-               else yield (BED3 chr (gateZero $ tss - dw) (tss + up), x)-    concatRegions = concatMapC $ \(locA, locB) ->-        map ((,) locB) (intersect locA) ++ map ((,) locA) (intersect locB)-    basal = bedToTree (++) $ flip map genes $ \((chr, tss, str), x) ->-        if str then (BED3 chr (tss - up) (tss + dw), [x])-               else (BED3 chr (tss - dw) (tss + up), [x])-    intersect = nubSort . concat . IM.elems . intersecting basal-    gateZero x | x < 0 = 0-               | otherwise = x-{-# INLINE get3DRegulatoryDomains #-}

@@ -75,8 +75,8 @@  -- | GO Annotation File (GAF) Format 2.1 Parser. For details read: -- http://geneontology.org/page/go-annotation-file-gaf-format-21.-readGAF :: FilePath -> Source (ResourceT IO) GAF-readGAF input = sourceFileBS input =$= linesUnboundedAsciiC =$=+readGAF :: FilePath -> ConduitT i GAF (ResourceT IO) ()+readGAF input = sourceFileBS input .| linesUnboundedAsciiC .|     (dropWhileC isCom >> mapC parseLine)   where     isCom l = B.head l == '!' || B.null l

@@ -132,7 +132,7 @@ {-# INLINE scores #-}  -- | A streaming version of scores.-scores' :: Monad m => Bkgd -> PWM -> DNA a -> Source m Double+scores' :: Monad m => Bkgd -> PWM -> DNA a -> ConduitT i Double m () scores' bg p@(PWM _ pwm) dna = go 0   where     go i | i < n - len + 1 = do yield $ scoreHelp bg p $ B.take len $ B.drop i s

@@ -11,8 +11,7 @@     , spaceConstraintHelper     ) where -import           Conduit                     (Source, await, mapC, sinkVector,-                                              yield, ($$), (=$=))+import           Conduit import           Control.Arrow               ((***)) import           Control.Monad.ST            (runST) import           Control.Parallel.Strategies (parMap, rdeepseq)@@ -41,7 +40,7 @@          -> Double          -> Bool    -- ^ whether to skip ambiguous sequences. Recommend: True                     -- in most cases-         -> Source m Int+         -> ConduitT i Int m () findTFBS bg pwm dna thres skip = loop 0   where     loop !i | i >= l - n + 1 = return ()@@ -81,9 +80,8 @@ {-# INLINE findTFBS' #-}  -- | use naive search-findTFBSSlow :: Monad m => Bkgd -> PWM -> DNA a -> Double -> Source m Int-findTFBSSlow bg pwm dna thres = scores' bg pwm dna =$=-                                       loop 0+findTFBSSlow :: Monad m => Bkgd -> PWM -> DNA a -> Double -> ConduitT i Int m ()+findTFBSSlow bg pwm dna thres = scores' bg pwm dna .| loop 0   where     loop i = do v <- await                 case v of@@ -184,7 +182,7 @@               'C' -> log $! matchC / c               'G' -> log $! matchG / g               'T' -> log $! matchT / t-              _ -> -1 / 0+              _   -> -1 / 0         matchA = addSome $ M.unsafeIndex (_mat pwm) (d,0)         matchC = addSome $ M.unsafeIndex (_mat pwm) (d,1)         matchG = addSome $ M.unsafeIndex (_mat pwm) (d,2)@@ -195,20 +193,6 @@     n = size pwm {-# INLINE lookAheadSearch' #-} -{--data MotifCompo = MotifCompo-    { _motif1        :: Motif-    , _motif2        :: Motif-    , _sameDirection :: Bool-    , _spacing       :: Int-    , _score         :: Double-    }--instance Show MotifCompo where-    show (MotifCompo m1 m2 isSame sp sc) = intercalate "\t"-        [B.unpack $ _name m1, B.unpack $ _name m2, show isSame, show sp, show sc]-        -}- data SpaceDistribution = SpaceDistribution     { _motif1   :: Motif     , _nSites1  :: (Int, Int)@@ -235,9 +219,9 @@     motifs = nubBy ((==) `on` _name) $ concatMap (\(a,b) -> [a,b]) pairs     findSites (Motif _ pwm) = (fwd, rev)       where-        fwd = runST $ findTFBS bg pwm dna cutoff True $$ sinkVector-        rev = runST $ findTFBS bg pwm' dna cutoff' True =$=-              mapC (+ (size pwm - 1)) $$ sinkVector+        fwd = runST $ runConduit $ findTFBS bg pwm dna cutoff True .| sinkVector+        rev = runST $ runConduit $ findTFBS bg pwm' dna cutoff' True .|+              mapC (+ (size pwm - 1)) .| sinkVector         cutoff = pValueToScore th bg pwm         cutoff' = pValueToScore th bg pwm'         pwm' = rcPWM pwm

@@ -3,43 +3,53 @@ module Bio.RealWorld.GENCODE     ( Gene(..)     , readGenes-    , readGenes'-    , parseGenes     ) where  import           Conduit import qualified Data.ByteString.Char8 as B-import           Data.Maybe            (fromJust) import           Data.CaseInsensitive  (CI, mk)+import qualified Data.HashMap.Strict   as M+import           Data.Maybe            (fromJust)  import           Bio.Utils.Misc        (readInt) +-- | GTF's position is 1-based, but here we convert it to 0-based indexing. data Gene = Gene-    { geneName   :: !(CI B.ByteString)-    , geneId     :: !B.ByteString-    , geneChrom  :: !B.ByteString-    , geneStart  :: !Int-    , geneEnd    :: !Int-    , geneStrand :: !Bool+    { geneName        :: !(CI B.ByteString)+    , geneId          :: !B.ByteString+    , geneChrom       :: !B.ByteString+    , geneLeft        :: !Int+    , geneRight       :: !Int+    , geneStrand      :: !Bool+    , geneTranscripts :: [(Int, Int)]     } deriving (Show)  -- | Read gene information from Gencode GTF file-readGenes :: MonadResource m => FilePath -> Source m Gene-readGenes input = sourceFile input =$= parseGenes--readGenes' :: FilePath -> IO [Gene]-readGenes' input = runResourceT $ readGenes input $$ sinkList--parseGenes :: Monad m => Conduit B.ByteString m Gene-parseGenes = linesUnboundedAsciiC =$= concatMapC f+readGenes :: FilePath -> IO [Gene]+readGenes input = do+    (genes, transcripts) <- runResourceT $ runConduit $ sourceFile input .|+        linesUnboundedAsciiC .| foldlC f (M.empty, M.empty)+    return $ M.elems $ M.foldlWithKey'+        (\m k v -> M.adjust (\g -> g{geneTranscripts=v}) k m)+        genes transcripts   where-    f l | B.head l == '#' || f3 /= "gene" = Nothing-        | otherwise = Just $ Gene (mk $ getField "gene_name") (getField "gene_id") f1-            (readInt f4 - 1) (readInt f5) (f7=="+")+    f (genes, transcripts) l+        | B.head l == '#' = (genes, transcripts)+        | f3 == "gene" =+            let g = Gene (mk $ getField "gene_name") i f1 (readInt f4 - 1)+                    (readInt f5 - 1) (f7=="+") []+                i = getField "gene_id"+            in (M.insert i g genes, transcripts)+        | f3 == "transcript" =+            let t = (readInt f4 - 1, readInt f5 - 1)+                i = getField "gene_id"+                updateFn Nothing = Just [t]+                updateFn (Just x) = Just $ t : x+            in (genes, M.alter updateFn i transcripts)+        | otherwise = (genes, transcripts)       where         [f1,_,f3,f4,f5,_,f7,_,f9] = B.split '\t' l         fields = map (B.break (==' ') . strip) $ B.split ';' f9         getField x = B.init $ B.drop 2 $ fromJust $ lookup x fields     strip = fst . B.spanEnd isSpace . B.dropWhile isSpace     isSpace = (== ' ')-{-# INLINE parseGenes #-}

@@ -37,7 +37,7 @@ getJunction g = (_chrom g, U.map fst $ _introns g)  -- | read genes from UCSC "knownGenes.tsv"-readUCSCGenes :: FilePath -> Source IO UCSCGene+readUCSCGenes :: FilePath -> ConduitT i UCSCGene IO () readUCSCGenes fl = do     handle <- liftIO $ openFile fl ReadMode     _ <- liftIO $ B.hGetLine handle   -- header@@ -54,7 +54,7 @@ {-# INLINE readUCSCGenes #-}  readUCSCGenes' :: FilePath -> IO [UCSCGene]-readUCSCGenes' fl = readUCSCGenes fl $$ sinkList+readUCSCGenes' fl = runConduit $ readUCSCGenes fl .| sinkList {-# INLINE readUCSCGenes' #-}  readGeneFromLine :: B.ByteString -> UCSCGene

@@ -98,7 +98,7 @@     B.hPutStr outH $ B.unlines [magic, mkHeader chrs]     mapM_ (B.hPutStr outH) dnas   where-    readSeq fl = runResourceT $ fastaReader fl =$= conduit $$ sinkList+    readSeq fl = runResourceT $ runConduit $ fastaReader fl .| conduit .| sinkList       where         conduit = awaitForever $ \(chrName, seqs) -> do             let dna = B.concat seqs

@@ -85,8 +85,8 @@             in loop k' ak' bk' espk'         | otherwise = 1 - (ak / bk - epsk / 2)     s = foldl' (\s' k -> 1 + s' * invJm _n x _N k) 1.0 [x..m-2]-    invJm _n x _N m = ( 1 - fromIntegral x / fromIntegral (m+1) ) /-                          ( 1 - fromIntegral (_n-1-x) / fromIntegral (_N-1-m) )+    invJm _n _x _N _m = ( 1 - fromIntegral _x / fromIntegral (_m+1) ) /+                          ( 1 - fromIntegral (_n-1-_x) / fromIntegral (_N-1-_m) )     e = 1e-20  

@@ -6,6 +6,7 @@  import           Bio.Data.Bed import           Conduit+import           Control.Lens                ((^.)) import           Control.Monad import qualified Data.ByteString.Char8       as B import           Data.Function@@ -24,55 +25,17 @@         create xs = (fst.fst.head $ xs, IM.fromDistinctAscList.map f $ xs) {-# INLINE toMap #-} -{---- | coverages of bins--- FIXME: Too ugly coverage :: [BED]  -- ^ genomic locus in BED format-         -> [BED]  -- ^ reads in BED format-         -> (V.Vector Double, Int)-coverage bin tags = getResult (V.create (VM.replicate (n+1) 0 >>= go tags))-  where-    getResult v = (V.zipWith normalize (V.slice 0 n v) featWidth, v V.! n)-    go ts v = do-        forM_ ts (\t -> do-            let set = M.lookup (t^.chrom) featMap-                s = t^.chromStart-                e = t^.chromEnd-                b = (s, e)-                l = s - e + 1-                intervals = case set of-                    Just iMap -> IM.intersecting iMap.toInterval $ b-                    _ -> []-            forM_ intervals (\interval -> do-                let i = snd interval-                    nucl = overlap b . fst $ interval-                VM.write v i . (+nucl) =<< VM.read v i-                )-            VM.write v n . (+l) =<< VM.read v n-            )-        return v-    featMap = toMap.map (\x -> (x^.chrom, (x^.chromStart, x^.chromEnd))) $ bin-    featWidth = V.fromList.map (\x -> x^.chromEnd - x^.chromStart) $ bin-    n = length bin-    overlap (l, u) (IM.ClosedInterval l' u')-        | l' >= l = if u' <= u then u'-l'+1 else u-l'+1-        | otherwise = if u' <= u then u'-l+1 else u-l+1-    overlap _ _ = 0-    normalize a b = fromIntegral a / fromIntegral b-    -}--coverage :: [BED]  -- ^ genomic locus in BED format-         -> Source IO BED  -- ^ reads in BED format+         -> ConduitT () BED IO ()  -- ^ reads in BED format          -> IO (V.Vector Double, Int)-coverage bin tags = liftM getResult $ tags $$ sink+coverage bin tags = liftM getResult $ runConduit $ tags .| sink   where-    sink :: Sink BED IO (V.Vector Int)     sink = do         v <- lift $ VM.replicate (n+1) 0         mapM_C $ \t -> do-                let set = M.lookup (_chrom t) featMap-                    s = _chromStart t-                    e = _chromEnd t+                let set = M.lookup (t^.chrom) featMap+                    s = t^.chromStart+                    e = t^.chromEnd                     b = (s, e)                     l = e - s + 1                     intervals = case set of@@ -86,8 +49,8 @@                 VM.write v n . (+l) =<< VM.read v n         lift $ V.freeze v     getResult v = (V.zipWith normalize (V.slice 0 n v) featWidth, v V.! n)-    featMap = toMap.map (\x -> (_chrom x, (_chromStart x, _chromEnd x))) $ bin-    featWidth = V.fromList.map (\x -> _chromEnd x - _chromStart x) $ bin+    featMap = toMap.map (\x -> (x^.chrom, (x^.chromStart, x^.chromEnd))) $ bin+    featWidth = V.fromList $ map size bin     n = length bin     overlap (l, u) (IM.ClosedInterval l' u')         | l' >= l = if u' <= u then u'-l'+1 else u-l'+1

@@ -37,7 +37,7 @@  sortedBamToBedPETest :: Assertion sortedBamToBedPETest = do-    bedpe <- readBedPE "tests/data/pairedend.bedpe"+    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@@ -46,5 +46,5 @@     readBedPE fl = do         c <- B.readFile fl         return $ map (f . B.split '\t') $ B.lines c-    f (f1:f2:f3:f4:f5:f6:_) = ( BED3 f1 (readInt f2) (readInt f3)-                              , BED3 f4 (readInt f5) (readInt f6) )+    f (f1:f2:f3:f4:f5:f6:_) = ( asBed f1 (readInt f2) (readInt f3)+                              , asBed f4 (readInt f5) (readInt f6) )

@@ -27,7 +27,7 @@ splitBedTest :: Assertion splitBedTest = (s1', s2', s3') @=? (s1, s2, s3)   where-    bed = BED3 "chr1" 0 99+    bed = asBed "chr1" 0 99     s1 = splitBedBySize 20 bed     s1' = map f [(0, 20), (20, 40), (40, 60), (60, 80)]     s2 = splitBedBySizeLeft 20 bed@@ -35,12 +35,13 @@     s3 = splitBedBySizeOverlap 20 10 bed     s3' = map f [ ( 0, 20), (10, 30), (20, 40), (30, 50), (40, 60)                 , (50, 70), (60, 80), (70, 90), (80, 100) ]-    f (a,b) = asBed "chr1" a b+    f (a,b) = asBed "chr1" a b :: BED3  splitOverlappedTest :: Assertion splitOverlappedTest = expect @=? result   where-    input = map (\(a,b) -> BED3 "chr1" a b)+    input :: [BED3]+    input = map (\(a,b) -> asBed "chr1" a b)         [ (0, 100)         , (10, 20)         , (50, 150)@@ -48,7 +49,7 @@         , (155, 200)         , (155, 220)         ]-    expect = map (\((a,b), x) -> (BED3 "chr1" a b, x))+    expect = map (\((a,b), x) -> (asBed "chr1" a b, x))         [ ((0, 10), 1)         , ((10, 20), 2)         , ((20, 50), 1)

@@ -1,56 +0,0 @@-{-# LANGUAGE OverloadedStrings #-}-module Tests.GREAT (tests) where--import           Bio.Data.Bed-import           Bio.GO.GREAT-import           Conduit-import           Control.Monad.Identity (runIdentity)-import           Data.Function          (on)-import           Data.List-import           Data.Ord--import           Test.Tasty-import           Test.Tasty.HUnit--genes = [ (("chr1", 1000, True), 0)-        , (("chr1", 6000, True), 1)-        , (("chr2", 6000, True), 2)-        , (("chr2", 2000000, False), 3)-        , (("chr2", 2020000, True), 4)-        ]--results = [ (BED3 "chr1" 0 2000, 0)-          , (BED3 "chr1" 1000 1007000, 1)-          , (BED3 "chr2" 0 1007000, 2)-          , (BED3 "chr2" 999000 2015000, 3)-          , (BED3 "chr2" 2005000 3021000, 4)-          ]--loops = [ (BED3 "chr1" 900 1200, BED3 "chr1" 10000 11000)-        , (BED3 "chr2" 0 100, BED3 "chr2" 5500 6500)-        ]--results3D = [ (BED3 "chr1" 0 2000, 0)-            , (BED3 "chr1" 1000 7000, 1)-            , (BED3 "chr2" 1000 7000, 2)-            , (BED3 "chr2" 1999000 2005000, 3)-            , (BED3 "chr1" 10000 11000, 0)-            , (BED3 "chr1" 10000 11000, 1)-            , (BED3 "chr2" 0 100, 2)-            , (BED3 "chr2" 2015000 2021000,4)-            ]--tests :: TestTree-tests = testGroup "Test: Bio.GO.GREAT"-    [ testCase "getRegulatoryDomains" testAssoc-    , testCase "get3DRegulatoryDomains" testAssoc3D-    ]--testAssoc :: Assertion-testAssoc = results @=? sortBy (comparing snd)-    (getRegulatoryDomains (BasalPlusExtension 5000 1000 1000000) genes)--testAssoc3D :: Assertion-testAssoc3D = sortBy (compareBed `on` fst) results3D @=?-    sortBy (compareBed `on` fst) (runIdentity $ yieldMany loops =$=-        get3DRegulatoryDomains genes 5000 1000 $$ sinkList)

@@ -2,7 +2,6 @@ import qualified Tests.Bam as Bam import qualified Tests.Motif as Motif import qualified Tests.Seq as Seq-import qualified Tests.GREAT as GREAT import qualified Tests.Tools as Tools import Test.Tasty @@ -12,6 +11,5 @@     , Bam.tests     , Seq.tests     , Motif.tests-    , GREAT.tests     , Tools.tests     ]