packages feed

HsHTSLib 1.9.2 → 1.9.2.2

raw patch · 6 files changed

+108/−18 lines, 6 filesbinary-added

Files

HsHTSLib.cabal view
@@ -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
src/Bio/HTS/BAM.chs view
@@ -3,8 +3,11 @@ {-# LANGUAGE LambdaCase #-}  module Bio.HTS.BAM-    ( getBamHeader+    ( -- * BAM file header+      getBamHeader     , getSortOrder++      -- * BAM IO     , streamBam     , sinkBam
src/Bio/HTS/Types.hs view
@@ -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--}
src/Bio/HTS/Utils.hs view
@@ -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 #-}
+ tests/data/paired_end.bam view

binary file changed (absent → 1365169 bytes)

+ tests/data/paired_end_dedup.bam view

binary file changed (absent → 287553 bytes)