diff --git a/bioinformatics-toolkit.cabal b/bioinformatics-toolkit.cabal
--- a/bioinformatics-toolkit.cabal
+++ b/bioinformatics-toolkit.cabal
@@ -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
diff --git a/src/Bio/ChIPSeq.hs b/src/Bio/ChIPSeq.hs
deleted file mode 100644
--- a/src/Bio/ChIPSeq.hs
+++ /dev/null
@@ -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 #-}
diff --git a/src/Bio/Data/Bam.hs b/src/Bio/Data/Bam.hs
--- a/src/Bio/Data/Bam.hs
+++ b/src/Bio/Data/Bam.hs
@@ -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 #-}
diff --git a/src/Bio/Data/Bed.hs b/src/Bio/Data/Bed.hs
--- a/src/Bio/Data/Bed.hs
+++ b/src/Bio/Data/Bed.hs
@@ -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 #-}
diff --git a/src/Bio/Data/Bed/Types.hs b/src/Bio/Data/Bed/Types.hs
--- a/src/Bio/Data/Bed/Types.hs
+++ b/src/Bio/Data/Bed/Types.hs
@@ -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)
diff --git a/src/Bio/Data/Bed/Utils.hs b/src/Bio/Data/Bed/Utils.hs
new file mode 100644
--- /dev/null
+++ b/src/Bio/Data/Bed/Utils.hs
@@ -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 #-}
diff --git a/src/Bio/Data/Fastq.hs b/src/Bio/Data/Fastq.hs
--- a/src/Bio/Data/Fastq.hs
+++ b/src/Bio/Data/Fastq.hs
@@ -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
diff --git a/src/Bio/Motif.hs b/src/Bio/Motif.hs
--- a/src/Bio/Motif.hs
+++ b/src/Bio/Motif.hs
@@ -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
diff --git a/src/Bio/Motif/Merge.hs b/src/Bio/Motif/Merge.hs
--- a/src/Bio/Motif/Merge.hs
+++ b/src/Bio/Motif/Merge.hs
@@ -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
diff --git a/src/Bio/RealWorld/GDC.hs b/src/Bio/RealWorld/GDC.hs
new file mode 100644
--- /dev/null
+++ b/src/Bio/RealWorld/GDC.hs
@@ -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 #-}
diff --git a/src/Bio/RealWorld/Reactome.hs b/src/Bio/RealWorld/Reactome.hs
new file mode 100644
--- /dev/null
+++ b/src/Bio/RealWorld/Reactome.hs
@@ -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 =
+-}
diff --git a/src/Bio/Utils/BitVector.hs b/src/Bio/Utils/BitVector.hs
new file mode 100644
--- /dev/null
+++ b/src/Bio/Utils/BitVector.hs
@@ -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
diff --git a/tests/Tests/Bam.hs b/tests/Tests/Bam.hs
--- a/tests/Tests/Bam.hs
+++ b/tests/Tests/Bam.hs
@@ -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
diff --git a/tests/Tests/Bed.hs b/tests/Tests/Bed.hs
--- a/tests/Tests/Bed.hs
+++ b/tests/Tests/Bed.hs
@@ -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
diff --git a/tests/Tests/ChIPSeq.hs b/tests/Tests/ChIPSeq.hs
deleted file mode 100644
--- a/tests/Tests/ChIPSeq.hs
+++ /dev/null
@@ -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'
-              -}
diff --git a/tests/Tests/Motif.hs b/tests/Tests/Motif.hs
--- a/tests/Tests/Motif.hs
+++ b/tests/Tests/Motif.hs
@@ -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"
diff --git a/tests/data/motifs.meme b/tests/data/motifs.meme
new file mode 100644
--- /dev/null
+++ b/tests/data/motifs.meme
@@ -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