diff --git a/HsHTSLib.cabal b/HsHTSLib.cabal
--- a/HsHTSLib.cabal
+++ b/HsHTSLib.cabal
@@ -1,5 +1,5 @@
 name:                HsHTSLib
-version:             1.9.2
+version:             1.9.2.2
 synopsis:            Bindings to htslib.
 description:         This package provides bindings to htslib, a library
                      for processing high throughput DNA sequencing data.
@@ -7,7 +7,7 @@
 license-file:        LICENSE
 author:              Kai Zhang
 maintainer:          kai@kzhang.org
-copyright:           (c) 2016-2019 Kai Zhang
+copyright:           (c) 2016-2021 Kai Zhang
 category:            Bio
 build-type:          Simple
 cabal-version:       >=1.10
@@ -26,6 +26,8 @@
   tests/data/example.sam
   tests/data/single_end.bam
   tests/data/single_end_dedup.bam
+  tests/data/paired_end.bam
+  tests/data/paired_end_dedup.bam
 
 library
   ghc-options:         -Wall
@@ -42,6 +44,7 @@
     , bytestring-lexing
     , conduit >= 1.3.0
     , containers
+    , vector
 
   hs-source-dirs:      src
   build-tools:         c2hs >= 0.25.0
diff --git a/src/Bio/HTS/BAM.chs b/src/Bio/HTS/BAM.chs
--- a/src/Bio/HTS/BAM.chs
+++ b/src/Bio/HTS/BAM.chs
@@ -3,8 +3,11 @@
 {-# LANGUAGE LambdaCase #-}
 
 module Bio.HTS.BAM
-    ( getBamHeader
+    ( -- * BAM file header
+      getBamHeader
     , getSortOrder
+
+      -- * BAM IO
     , streamBam
     , sinkBam
 
diff --git a/src/Bio/HTS/Types.hs b/src/Bio/HTS/Types.hs
--- a/src/Bio/HTS/Types.hs
+++ b/src/Bio/HTS/Types.hs
@@ -27,8 +27,10 @@
 -- | The BAM format.
 newtype BAM = BAM { unbam :: ForeignPtr Bam1 }
 
+-- | The BAM file header.
 newtype BAMHeader = BAMHeader {unbamHeader :: ForeignPtr BamHdr}
 
+-- | Convert bam file header to string.
 showBamHeader :: BAMHeader -> B.ByteString
 showBamHeader header = unsafePerformIO $
     withForeignPtr (unbamHeader header) $ \ptr -> do
@@ -45,10 +47,12 @@
 
 newtype CIGAR = CIGAR [(Int, Char)]
 
+-- | Convert CIGAR to string.
 cigar2String :: CIGAR -> B.ByteString
 cigar2String (CIGAR c) = B.concat $
     concatMap (\(i, x) -> [fromJust $ packDecimal i, B.singleton x]) c
 
+-- | Read CIGAR from string.
 string2Cigar :: B.ByteString -> CIGAR
 string2Cigar c = CIGAR $ go c
   where
@@ -74,6 +78,7 @@
     , _sam_aux   :: [((Char,Char), AuxiliaryData)]
     }
 
+-- | Convert SAM to string.
 showSam :: SAM -> B.ByteString
 showSam SAM{..} = B.intercalate "\t" $
     [ _sam_qname, pack' _sam_flag, fromMaybe "*" _sam_rname, pack' _sam_pos
@@ -95,6 +100,7 @@
                    | AuxFloatArray [Float]
                    deriving (Show)
 
+-- | Convert aux data to string.
 showAuxiliaryData :: ((Char, Char), AuxiliaryData) -> B.ByteString
 showAuxiliaryData ((x1,x2), aux) = B.pack [x1,x2] <> aux'
   where
@@ -104,8 +110,3 @@
         AuxFloat x -> B.pack $ ":f:" <> show x
         AuxString x -> ":Z:" <> x
         _ -> error "Not implemented"
-
-{-
--- | SAM record flag
-newtype Flag = Flag Word16
--}
diff --git a/src/Bio/HTS/Utils.hs b/src/Bio/HTS/Utils.hs
--- a/src/Bio/HTS/Utils.hs
+++ b/src/Bio/HTS/Utils.hs
@@ -3,10 +3,16 @@
 {-# LANGUAGE RecordWildCards #-}
 {-# LANGUAGE PartialTypeSignatures #-}
 module Bio.HTS.Utils
-    ( markDupBy
-    , makeKey
+    ( -- * Mark duplicates
+      markDupBy
+    , Orientation(..)
     , BAMKey
-    , Orientation
+    , makeKey
+    , makeKeySingle
+    , makeKeyPair
+    
+      -- * Other utilities
+    , fragmentSizeDistr
     ) where
 
 import           Conduit
@@ -15,6 +21,8 @@
 import Data.Maybe
 import qualified Data.Sequence as S
 import qualified Data.ByteString.Char8 as B
+import qualified Data.Vector.Unboxed as U
+import qualified Data.Vector.Unboxed.Mutable as UM
 
 import           Bio.HTS.BAM
 import           Bio.HTS.Types
@@ -25,16 +33,73 @@
                      , _loc1 :: Int
                      , _orientation :: Orientation
                      , _barcode :: Maybe B.ByteString }
-            | Pair { _ref_id1 :: Int
-                   , _ref_id2 :: Int
-                   , _loc1 :: Int
-                   , _loc2 :: Int
+            | Pair { _ref_id1 :: Int  -- ^ ref id of this tag
+                   , _ref_id2 :: Int  -- ^ ref id of the paired tag
+                   , _loc1 :: Int     -- ^ location of this tag
+                   , _loc2 :: Int     -- ^ location of the paired tag
                    , _orientation :: Orientation
-                   , _leftmost :: Bool
+                   , _leftmost :: Bool   -- ^ Is this tag leftmost
                    , _barcode :: Maybe B.ByteString }
             deriving (Eq, Ord, Show)
 
-makeKey :: (BAM -> Maybe B.ByteString)   -- ^ Get Barcode
+-- | Generate fingerprint for single-end reads.
+makeKeySingle :: (BAM -> Maybe B.ByteString)   -- ^ Barcode extraction function
+              -> BAM
+              -> BAMKey
+makeKeySingle fn bam = Single ref1 (if isFwd1 then lloc1 else rloc1)
+    (if isFwd1 then F else R) $ fn bam
+  where
+    ref1 = refId bam
+    isFwd1 = not $ isRC flg
+    lloc1 = startLoc bam - fst clipped1 + 1
+    rloc1 = endLoc bam + snd clipped1
+    flg = flag bam
+    clipped1 = getClipped $ fromJust $ cigar bam
+{-# INLINE makeKeySingle #-}
+
+-- | Create a pair of keys (single-end and paired-end).
+makeKeyPair :: (BAM -> Maybe B.ByteString)   -- ^ Barcode extraction function
+            -> (BAM, BAM)
+            -> BAMKey
+makeKeyPair fn (bam1, bam2) = Pair ref1 ref2 loc1 loc2 orientation isLeftMost $ fn bam1
+  where
+    ref1 = refId bam1
+    lloc1 = startLoc bam1 - fst clipped1 + 1
+    rloc1 = endLoc bam1 + snd clipped1
+    clipped1 = getClipped $ fromJust $ cigar bam1
+    isFwd1 = not $ isRC $ flag bam1
+
+    ref2 = mateRefId bam2
+    lloc2 = startLoc bam2 - fst clipped2 + 1
+    rloc2 = endLoc bam2 + snd clipped2
+    clipped2 = getClipped $ fromJust $ cigar bam2
+    isFwd2 = not $ isRC $ flag bam2
+
+    isLeftMost
+        | ref1 /= ref2 = ref1 < ref2
+        | otherwise = if isFwd1 == isFwd2
+            then if isFwd1 then lloc1 <= lloc2 else rloc1 <= rloc2
+            else if isFwd1 then lloc1 <= rloc2 else rloc1 <= lloc2
+    orientation
+        | isLeftMost = if isFwd1 == isFwd2
+            then if isFwd1
+                then if isFirstSegment flg then FF else RR
+                else if isFirstSegment flg then RR else FF
+            else if isFwd1 then FR else RF
+        | otherwise = if isFwd1 == isFwd2
+            then if isFwd1
+                then if isFirstSegment flg then RR else FF
+                else if isFirstSegment flg then FF else RR
+            else if isFwd1 then RF else FR
+    loc1 | isFwd1 == isFwd2 = if isLeftMost then lloc1 else rloc1
+         | otherwise = if isFwd1 then lloc1 else rloc1
+    loc2 | isFwd1 == isFwd2 = if isLeftMost then rloc2 else lloc2
+         | otherwise = if isFwd1 then rloc2 else lloc2
+    flg = flag bam1
+
+
+-- | Create a pair of keys (single-end and paired-end).
+makeKey :: (BAM -> Maybe B.ByteString)   -- ^ Barcode extraction function
         -> BAM
         -> (BAMKey, BAMKey)
 makeKey fn bam = (single, pair)
@@ -94,6 +159,7 @@
 
 -- | Remove duplicated reads. Duplicates are determined by
 -- checking for matching keys. The Key is comprised of:
+--
 -- 1. Chromosome
 -- 2. Orientation (forward/reverse)
 -- 3. Unclipped Start(forward)/End(reverse)
@@ -102,7 +168,7 @@
 -- Keep the read that has a higher base quality sum (sum of all
 -- base qualities in the record above 15).
 markDupBy :: MonadIO m
-          => (BAM -> Maybe B.ByteString)   -- ^ Get Barcode
+          => (BAM -> Maybe B.ByteString)   -- ^ Barcode extraction function, if any.
           -> ConduitT BAM BAM m ()
 markDupBy bcFn = go (-1,-1) M.empty S.empty
   where
@@ -144,3 +210,20 @@
     isSorted (chr1, loc1) (chr2, loc2) = chr1 < chr2 ||
         (chr1 == chr2 && loc1 <= loc2)
 {-# INLINE markDupBy #-}
+
+
+-- | Compute fragment size distribution from paired end BAM records.
+fragmentSizeDistr :: PrimMonad m
+                  => Int    -- ^ Largest fragment size
+                  -> ConduitT BAM o m (U.Vector Double)
+fragmentSizeDistr n = do
+    vec <- lift $ UM.replicate n 0
+    mapM_C $ f vec
+    vec' <- lift $ U.unsafeFreeze vec
+    return $ U.map (/ (U.sum vec')) vec'
+  where
+    f v x | s >= n = return ()
+          | otherwise = UM.modify v (+1) s
+      where
+        s = abs $ tLen x
+{-# INLINE fragmentSizeDistr #-}
diff --git a/tests/data/paired_end.bam b/tests/data/paired_end.bam
new file mode 100644
Binary files /dev/null and b/tests/data/paired_end.bam differ
diff --git a/tests/data/paired_end_dedup.bam b/tests/data/paired_end_dedup.bam
new file mode 100644
Binary files /dev/null and b/tests/data/paired_end_dedup.bam differ