diff --git a/LICENSE b/LICENSE
--- a/LICENSE
+++ b/LICENSE
@@ -1,6 +1,6 @@
 The MIT License (MIT)
 
-Copyright (c) 2014-2016 Kai Zhang
+Copyright (c) 2014-2020 Kai Zhang
 
 Permission is hereby granted, free of charge, to any person obtaining a copy
 of this software and associated documentation files (the "Software"), to deal
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.9.0
+version:             0.9.1
 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-2019 Kai Zhang
+copyright:           (c) 2014-2020 Kai Zhang
 category:            Bio
 build-type:          Simple
 extra-source-files:  README.md
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
@@ -35,7 +35,8 @@
     , mergeBedWith
     , mergeSortedBed
     , mergeSortedBedWith
-    , splitOverlapped
+--    , splitOverlapped
+    , countOverlapped
 
     -- * IO
     , streamBed
@@ -58,7 +59,7 @@
 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.List                    (groupBy, sortBy, group)
 import           Data.Ord                     (comparing)
 import           Data.Conduit.Zlib           (gzip, ungzip, multiple)
 import qualified Data.Vector                  as V
@@ -236,14 +237,16 @@
             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) . size . fst) $
-    evalState (F.foldrM f [] $ init xs') x0
+    evalState (F.foldrM f [] $ init anchors) $
+        (\(a,b) -> (fromEither a, M.singleton (b^.chromStart, b^.chromEnd) b)) $
+        last anchors
   where
-    x0 = (\(a,b) -> (fromEither a, M.singleton (b^.chromStart, b^.chromEnd) b)) $ last xs'
-    xs' = sortBy (comparing (fromEither . fst)) $ concatMap
+    anchors = sortBy (comparing (fromEither . fst)) $ concatMap
         ( \x -> [(Left $ x^.chromStart, x), (Right $ x^.chromEnd, x)] ) xs
     f (i, x) acc = do
         (j, s) <- get
@@ -257,7 +260,32 @@
     fromEither (Right x) = x
     chr = head xs ^. chrom
 {-# INLINE splitOverlapped #-}
+-}
 
+-- | Split overlapped regions into non-overlapped regions. The input must be overlapped.
+-- This function is usually used with `mergeBedWith`.
+countOverlapped :: BEDLike b => [b] -> [(BED3, Int)]
+countOverlapped xs = reverse $ (\(_,_,x) -> x) $
+    F.foldl' go (p0, c0, []) $ tail anchors
+  where
+    (Left p0, c0) = head anchors
+    go (x, accum, result) (Left p, count)
+        | x == p = (p, accum + count, result)
+        | otherwise = (p, accum + count, (asBed chr x $ p, accum):result)
+    go (x, accum, result) (Right p, count) =
+      (p+1, accum - count, (asBed chr x $ p + 1, accum):result)
+    anchors = map (\x -> (head x, length x)) $ group $ sortBy cmp $ concatMap
+        (\x -> [Left $ x^.chromStart, Right $ x^.chromEnd - 1]) xs
+    cmp (Left x) (Left y) = compare x y
+    cmp (Right x) (Right y) = compare x y
+    cmp (Left x) (Right y) = case compare x y of
+        EQ -> LT
+        o -> o
+    cmp (Right x) (Left y) = case compare x y of
+        EQ -> GT
+        o -> o
+    chr = head xs ^. chrom
+{-# INLINE countOverlapped #-}
 
 streamBed :: (MonadResource m, BEDConvert b, MonadIO m)
           => FilePath -> ConduitT i b m () 
diff --git a/tests/Tests/Bed.hs b/tests/Tests/Bed.hs
--- a/tests/Tests/Bed.hs
+++ b/tests/Tests/Bed.hs
@@ -57,9 +57,11 @@
         , (120, 160)
         , (155, 200)
         , (155, 220)
+        , (0, 10)
+        , (0, 10)
         ]
     expect = map (\((a,b), x) -> (asBed "chr1" a b, x))
-        [ ((0, 10), 1)
+        [ ((0, 10), 3)
         , ((10, 20), 2)
         , ((20, 50), 1)
         , ((50, 100), 2)
@@ -70,7 +72,8 @@
         , ((160, 200), 2)
         , ((200, 220), 1)
         ]
-    result = sortBy (compareBed `on` fst) $ splitOverlapped length input
+    --result = sortBy (compareBed `on` fst) $ splitOverlapped length input
+    result = sortBy (compareBed `on` fst) $ countOverlapped input
 
 mergeBedTest :: Assertion
 mergeBedTest = expect @=? result