diff --git a/bioinformatics-toolkit.cabal b/bioinformatics-toolkit.cabal
--- a/bioinformatics-toolkit.cabal
+++ b/bioinformatics-toolkit.cabal
@@ -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
diff --git a/src/Bio/ChIPSeq.hs b/src/Bio/ChIPSeq.hs
--- a/src/Bio/ChIPSeq.hs
+++ b/src/Bio/ChIPSeq.hs
@@ -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 #-}
diff --git a/src/Bio/ChIPSeq/FragLen.hs b/src/Bio/ChIPSeq/FragLen.hs
--- a/src/Bio/ChIPSeq/FragLen.hs
+++ b/src/Bio/ChIPSeq/FragLen.hs
@@ -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 #-}
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
@@ -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 #-}
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
@@ -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 #-}
diff --git a/src/Bio/Data/Bed/Types.hs b/src/Bio/Data/Bed/Types.hs
new file mode 100644
--- /dev/null
+++ b/src/Bio/Data/Bed/Types.hs
@@ -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)
diff --git a/src/Bio/Data/Fasta.hs b/src/Bio/Data/Fasta.hs
--- a/src/Bio/Data/Fasta.hs
+++ b/src/Bio/Data/Fasta.hs
@@ -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
diff --git a/src/Bio/GO/GREAT.hs b/src/Bio/GO/GREAT.hs
deleted file mode 100644
--- a/src/Bio/GO/GREAT.hs
+++ /dev/null
@@ -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 #-}
diff --git a/src/Bio/GO/Parser.hs b/src/Bio/GO/Parser.hs
--- a/src/Bio/GO/Parser.hs
+++ b/src/Bio/GO/Parser.hs
@@ -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
diff --git a/src/Bio/Motif.hs b/src/Bio/Motif.hs
--- a/src/Bio/Motif.hs
+++ b/src/Bio/Motif.hs
@@ -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
diff --git a/src/Bio/Motif/Search.hs b/src/Bio/Motif/Search.hs
--- a/src/Bio/Motif/Search.hs
+++ b/src/Bio/Motif/Search.hs
@@ -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
diff --git a/src/Bio/RealWorld/GENCODE.hs b/src/Bio/RealWorld/GENCODE.hs
--- a/src/Bio/RealWorld/GENCODE.hs
+++ b/src/Bio/RealWorld/GENCODE.hs
@@ -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 #-}
diff --git a/src/Bio/RealWorld/UCSC.hs b/src/Bio/RealWorld/UCSC.hs
--- a/src/Bio/RealWorld/UCSC.hs
+++ b/src/Bio/RealWorld/UCSC.hs
@@ -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
diff --git a/src/Bio/Seq/IO.hs b/src/Bio/Seq/IO.hs
--- a/src/Bio/Seq/IO.hs
+++ b/src/Bio/Seq/IO.hs
@@ -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
diff --git a/src/Bio/Utils/Functions.hs b/src/Bio/Utils/Functions.hs
--- a/src/Bio/Utils/Functions.hs
+++ b/src/Bio/Utils/Functions.hs
@@ -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
 
 
diff --git a/src/Bio/Utils/Overlap.hs b/src/Bio/Utils/Overlap.hs
--- a/src/Bio/Utils/Overlap.hs
+++ b/src/Bio/Utils/Overlap.hs
@@ -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
diff --git a/tests/Tests/Bam.hs b/tests/Tests/Bam.hs
--- a/tests/Tests/Bam.hs
+++ b/tests/Tests/Bam.hs
@@ -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) )
diff --git a/tests/Tests/Bed.hs b/tests/Tests/Bed.hs
--- a/tests/Tests/Bed.hs
+++ b/tests/Tests/Bed.hs
@@ -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)
diff --git a/tests/Tests/GREAT.hs b/tests/Tests/GREAT.hs
deleted file mode 100644
--- a/tests/Tests/GREAT.hs
+++ /dev/null
@@ -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)
diff --git a/tests/test.hs b/tests/test.hs
--- a/tests/test.hs
+++ b/tests/test.hs
@@ -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
     ]