HsHTSLib 1.9.1 → 1.9.2
raw patch · 12 files changed
+738/−504 lines, 12 filesdep +containersdep −mtlbinary-added
Dependencies added: containers
Dependencies removed: mtl
Files
- ChangeLog.md +8/−2
- HsHTSLib.cabal +9/−2
- cbits/hs_htslib.c +4/−0
- src/Bio/HTS.chs +0/−487
- src/Bio/HTS.hs +9/−0
- src/Bio/HTS/BAM.chs +513/−0
- src/Bio/HTS/Internal.chs +5/−0
- src/Bio/HTS/Types.hs +25/−7
- src/Bio/HTS/Utils.hs +146/−0
- tests/Main.hs +19/−6
- tests/data/single_end.bam binary
- tests/data/single_end_dedup.bam binary
ChangeLog.md view
@@ -1,5 +1,11 @@ # Revision history for HsHTSLib -## 1.9.0 -- YYYY-mm-dd+## 1.9.2 -- 2019-03-18 -* htslib-1.9+* Name changes: seqName -> queryName, getChrId -> refId, getChr -> refName, mateChrId -> mateRefId, mateChr -> mateRefName.+* Use 1-based coordinate system in SAM format.+* Add MarkdupBy++## 1.9.1 -- 2019-03-18++* Upgrade the C library to v1.9.
HsHTSLib.cabal view
@@ -1,5 +1,5 @@ name: HsHTSLib-version: 1.9.1+version: 1.9.2 synopsis: Bindings to htslib. description: This package provides bindings to htslib, a library for processing high throughput DNA sequencing data.@@ -14,13 +14,18 @@ extra-source-files: README.md ChangeLog.md++ -- C header files htslib-1.9/*.h htslib-1.9/htslib/*.h htslib-1.9/cram/*.h htslib-1.9/os/*.h + -- TEST DATA tests/data/example.bam tests/data/example.sam+ tests/data/single_end.bam+ tests/data/single_end_dedup.bam library ghc-options: -Wall@@ -28,12 +33,15 @@ Bio.HTS Bio.HTS.Internal Bio.HTS.Types+ Bio.HTS.BAM+ Bio.HTS.Utils build-depends: base >= 4.10 && < 5.0 , bytestring , bytestring-lexing , conduit >= 1.3.0+ , containers hs-source-dirs: src build-tools: c2hs >= 0.25.0@@ -110,7 +118,6 @@ , tasty-golden , tasty-hunit , conduit- , mtl source-repository head type: git
cbits/hs_htslib.c view
@@ -46,3 +46,7 @@ int bam_get_l_aux_(bam1_t *b) { return bam_get_l_aux(b); } +void bam_mark_dup(bam1_t *b) { b->core.flag |= BAM_FDUP; }+++
− src/Bio/HTS.chs
@@ -1,487 +0,0 @@-{-# LANGUAGE ForeignFunctionInterface #-}-{-# LANGUAGE OverloadedStrings #-}-{-# LANGUAGE LambdaCase #-}--module Bio.HTS- ( getBamHeader- , getSortOrder- , streamBam- , sinkBam-- -- * Field accessors- , getChrId- , getChr- , startLoc- , endLoc- , readLen- , isRev- , flag- , mapq- , getSeq- , seqName- , qualityS- , quality- , cigar- , mateChr- , mateChrId- , mateStartLoc- , tLen- , auxData- , queryAuxData- , bamToSam-- -- * Modify BAM- , appendAux-- -- * Flag interpretation- , hasMultiSegments- , isProperAligned- , isUnmapped- , isNextUnmapped- , isRC- , isNextRC- , isFirstSegment- , isLastSegment- , isSecondary- , isBadQual- , isDup- , isSupplementary- ) where--import Conduit-import Data.Bits (testBit)-import qualified Data.ByteString.Char8 as B-import qualified Data.ByteString as BS-import Data.Int-import Data.Word-import Foreign.C.String-import Foreign.ForeignPtr-import Foreign.Marshal.Array-import Foreign.Marshal.Alloc-import Foreign.Storable (poke, peekByteOff)-import Foreign.Ptr-import System.IO-import System.IO.Unsafe (unsafePerformIO)--import Bio.HTS.Internal-import Bio.HTS.Types--#include "htslib/sam.h"--getBamHeader :: FilePath -> IO BAMHeader-getBamHeader input = withHTSFile input ReadMode $ \hts ->- getBgzf hts >>= bamHdrRead >>= newForeignPtr bamHdrDestroy >>=- return . BAMHeader-{-# INLINE getBamHeader #-}-foreign import ccall unsafe "&bam_hdr_destroy"- bamHdrDestroy :: FunPtr (Ptr BamHdr -> IO ())---- | Get the sort information.-getSortOrder :: BAMHeader -> SortOrder-getSortOrder header = case lookup "SO" fields of- Just "unknown" -> Unknown- Just "unsorted" -> Unsorted- Just "queryname" -> Queryname- Just "coordinate" -> Coordinate- _ -> Unknown- where- fields = map (f . B.split ':') $ tail $ B.split '\t' $ head $ B.lines $- showBamHeader header- f [a,b] = (a,b)- f _ = error "Auxiliary field parsing failed!"---- | Create a Bam stream from the file.-streamBam :: MonadResource m => FilePath -> ConduitT i BAM m ()-streamBam input = bracketP (htsOpen input "r") htsClose $ \hts -> do- bgzf <- liftIO $ getBgzf hts- _ <- liftIO $ bamHdrRead bgzf >>= newForeignPtr bamHdrDestroy- loop bgzf- where- loop bgzf = do- ptr <- liftIO $ callocBytes {# sizeof bam1_t #}- code <- liftIO $ bamRead1 bgzf ptr- case () of- _ | code > 0 -> liftIO (newForeignPtr bamDestory1 ptr) >>=- yield . BAM >> loop bgzf- | code == -1 -> return ()- | code == -2 -> error "truncated file"- | otherwise -> error $ "read bam failed with code: " ++ show code-{-# INLINE streamBam #-}--foreign import ccall unsafe "&bam_destroy1"- bamDestory1 :: FunPtr (Ptr Bam1 -> IO ())---- | Write Bam records to a file.-sinkBam :: MonadResource m => FilePath -> BAMHeader -> ConduitT BAM o m ()-sinkBam output header = bracketP (htsOpen output "wb") htsClose $ \hts -> do- bgzf <- liftIO $ getBgzf hts- liftIO (withForeignPtr (unbamHeader header) $ bamHdrWrite bgzf) >>= \case- 0 -> mapM_C $ \bam -> liftIO $ do- code <- withForeignPtr (unbam bam) (bamWrite1 bgzf)- if code < 0- then error "bam_write1 failed"- else return ()- _ -> error "'bam_hdr_write' failed."----- | Return the chromosome id.-getChrId :: BAM -> Int-getChrId = unsafePerformIO . flip withForeignPtr fun . unbam- where- fun = fmap fromIntegral . {#get bam1_t->core.tid #}-{-# INLINE getChrId #-}---- | Return the chromosome name given the bam file header.-getChr :: BAMHeader -> BAM -> Maybe B.ByteString-getChr header bam- | chr < 0 = Nothing- | otherwise = Just $ unsafePerformIO $- withForeignPtr (unbamHeader header) $ \h ->- bamChr h chr >>= B.packCString- where- chr = getChrId bam-{-# INLINE getChr #-}---- | Return the 0-based starting location.-startLoc :: BAM -> Int-startLoc = unsafePerformIO . flip withForeignPtr fun . unbam- where- fun = fmap fromIntegral . {#get bam1_t->core.pos #}-{-# INLINE startLoc #-}---- | For a mapped read, this is just position + cigar2rlen.--- For an unmapped read (either according to its flags or if it has no cigar--- string), we return position + 1 by convention.-endLoc :: BAM -> Int-endLoc = unsafePerformIO . flip withForeignPtr bamEndpos .unbam-{-# INLINE endLoc #-}---- | Return the query length (read length).-readLen :: BAM -> Int-readLen = unsafePerformIO . flip withForeignPtr fun . unbam- where- fun = fmap fromIntegral . {#get bam1_t->core.l_qseq #}-{-# INLINE readLen #-}---- | Whether the query is on the reverse strand.-isRev :: BAM -> Bool-isRev = unsafePerformIO . flip withForeignPtr bamIsRev . unbam-{-# INLINE isRev #-}---- | Return the flag.-flag :: BAM -> Word16-flag = unsafePerformIO . flip withForeignPtr fn . unbam- where- fn = fmap fromIntegral . {#get bam1_t->core.flag #}-{-# INLINE flag #-}---- | MAPping Quality. It equals −10 log10 Pr{mapping position is wrong},--- rounded to the nearest integer. A value 255 indicates that the--- mapping quality is not available.-mapq :: BAM -> Word8-mapq = unsafePerformIO . flip withForeignPtr fn . unbam- where- fn = fmap fromIntegral . {#get bam1_t->core.qual #}-{-# INLINE mapq #-}---- | Return the DNA sequence.-getSeq :: BAM -> Maybe B.ByteString-getSeq = unsafePerformIO . flip withForeignPtr fn . unbam- where- fn b = fromIntegral <$> {#get bam1_t->core.l_qseq #} b >>= \case- 0 -> return Nothing- n -> allocaArray n $ \str -> do- bamGetSeq b str n- Just <$> B.packCStringLen (str, n)-{-# INLINE getSeq #-}---- | Get the name of the query.-seqName :: BAM -> B.ByteString-seqName = unsafePerformIO . flip withForeignPtr fn . unbam- where- fn b = {#get bam1_t->data #} b >>= B.packCString . castPtr-{-# INLINE seqName #-}---- | Human readable quality score which is: Phred base quality + 33.-qualityS :: BAM -> Maybe B.ByteString-qualityS = fmap (BS.map (+33)) . quality---- | Phred base quality (a sequence of 0xFF if absent).-quality :: BAM -> Maybe B.ByteString-quality = unsafePerformIO . flip withForeignPtr fn . unbam- where- fn b = fromIntegral <$> {#get bam1_t->core.l_qseq #} b >>= \case- 0 -> return Nothing- n -> allocaArray n $ \str -> bamGetQual b str n >>= \case- 0 -> Just <$> B.packCStringLen (str, n)- _ -> return Nothing-{-# INLINE quality #-}--cigar :: BAM -> Maybe [(Int, Char)]-cigar = unsafePerformIO . flip withForeignPtr fn . unbam- where- fn b = fromIntegral <$> {#get bam1_t->core.n_cigar #} b >>= \case- 0 -> return Nothing- n -> allocaArray n $ \num -> allocaArray n $ \str -> do- bamGetCigar b num str n- num' <- peekArray (fromIntegral n) num- str' <- peekArray (fromIntegral n) str- return $ Just $ zip (map fromIntegral num') $- map castCCharToChar str'-{-# INLINE cigar #-}--mateChrId :: BAM -> Int-mateChrId = unsafePerformIO . flip withForeignPtr fn . unbam- where- fn = fmap fromIntegral . {#get bam1_t->core.mtid #}-{-# INLINE mateChrId #-}--mateChr :: BAMHeader -> BAM -> Maybe B.ByteString-mateChr header bam- | chr < 0 = Nothing- | otherwise = Just $ unsafePerformIO $- withForeignPtr (unbamHeader header) $ \h ->- bamChr h chr >>= B.packCString- where- chr = mateChrId bam-{-# INLINE mateChr #-}---- | 0-based-mateStartLoc :: BAM -> Int-mateStartLoc = unsafePerformIO . flip withForeignPtr fn . unbam- where- fn = fmap fromIntegral . {#get bam1_t->core.mpos #}-{-# INLINE mateStartLoc #-}--tLen :: BAM -> Int-tLen = unsafePerformIO . flip withForeignPtr fn . unbam- where- fn = fmap fromIntegral . {#get bam1_t->core.isize #}-{-# INLINE tLen #-}--auxData :: BAM -> [((Char, Char), AuxiliaryData)]-auxData bam = unsafePerformIO $ withForeignPtr (unbam bam) $ \b -> do- l <- bamGetLAux b- aux <- bamGetAux b - go aux l- where- go ptr i- | i <= 0 = return []- | otherwise = do- name <- (,) <$> peekByteOff ptr 0 <*> peekByteOff ptr 1- castCCharToChar <$> peekByteOff ptr 2 >>= \case- 'A' -> do- r <- AuxChar . castCCharToChar <$> peekByteOff ptr 3- rs <- go (plusPtr ptr 4) $ i - 4- return $ (name, r) : rs- -- int8_t- 'c' -> do- r <- AuxInt . fromIntegral <$> (peekByteOff ptr 3 :: IO Int8)- rs <- go (plusPtr ptr 4) $ i - 4- return $ (name, r) : rs- -- uint8_t- 'C' -> do- r <- AuxInt . fromIntegral <$> (peekByteOff ptr 3 :: IO Word8)- rs <- go (plusPtr ptr 4) $ i - 4- return $ (name, r) : rs- -- int16_t- 's' -> do- r <- AuxInt . fromIntegral <$> (peekByteOff ptr 3 :: IO Int16)- rs <- go (plusPtr ptr 5) $ i - 5- return $ (name, r) : rs- -- uint16_t- 'S' -> do- r <- AuxInt . fromIntegral <$> (peekByteOff ptr 3 :: IO Word16)- rs <- go (plusPtr ptr 5) $ i - 5- return $ (name, r) : rs- -- int32_t- 'i' -> do- r <- AuxInt . fromIntegral <$> (peekByteOff ptr 3 :: IO Int32)- rs <- go (plusPtr ptr 7) $ i - 7- return $ (name, r) : rs- -- uint32_t- 'I' -> do- r <- AuxInt . fromIntegral <$> (peekByteOff ptr 3 :: IO Word32)- rs <- go (plusPtr ptr 7) $ i - 7- return $ (name, r) : rs- 'f' -> do- r <- AuxFloat <$> peekByteOff ptr 3- rs <- go (plusPtr ptr 7) $ i - 7- return $ (name, r) : rs- 'Z' -> do- str <- B.packCString (plusPtr ptr 3)- let l = B.length str + 1 + 3- rs <- go (plusPtr ptr l) $ i - l- return $ (name, AuxString str) : rs- 'H' -> do - str <- B.packCString (plusPtr ptr 3)- let l = B.length str + 1 + 3- rs <- go (plusPtr ptr l) $ i - l- return $ (name, AuxByteArray str) : rs- 'B' -> do- n <- fromIntegral <$> (peekByteOff ptr 4 :: IO Int32)- castCCharToChar <$> peekByteOff ptr 3 >>= \case- 'c' -> do- r <- AuxIntArray . map fromIntegral <$>- (peekArray n $ plusPtr ptr 8 :: IO [Int8])- let l = 8 + n- rs <- go (plusPtr ptr l) $ i - l- return $ (name, r) : rs- 'C' -> do- r <- AuxIntArray . map fromIntegral <$>- (peekArray n $ plusPtr ptr 8 :: IO [Word8])- let l = 8 + n- rs <- go (plusPtr ptr l) $ i - l- return $ (name, r) : rs- 's' -> do- r <- AuxIntArray . map fromIntegral <$>- (peekArray n $ plusPtr ptr 8 :: IO [Int16])- let l = 8 + 2 * n- rs <- go (plusPtr ptr l) $ i - l- return $ (name, r) : rs- 'S' -> do- r <- AuxIntArray . map fromIntegral <$>- (peekArray n $ plusPtr ptr 8 :: IO [Word16])- let l = 8 + 2 * n- rs <- go (plusPtr ptr l) $ i - l- return $ (name, r) : rs- 'i' -> do- r <- AuxIntArray . map fromIntegral <$>- (peekArray n $ plusPtr ptr 8 :: IO [Int32])- let l = 8 + 4 * n- rs <- go (plusPtr ptr l) $ i - l- return $ (name, r) : rs- 'I' -> do- r <- AuxIntArray . map fromIntegral <$>- (peekArray n $ plusPtr ptr 8 :: IO [Word32])- let l = 8 + 4 * n- rs <- go (plusPtr ptr l) $ i - l- return $ (name, r) : rs- 'f' -> do- r <- AuxFloatArray <$>- (peekArray n $ plusPtr ptr 8 :: IO [Float])- let l = 8 + 4 * n- rs <- go (plusPtr ptr l) $ i - l- return $ (name, r) : rs- x -> error $ "Unknown auxiliary record type: " ++ [x]- x -> error $ "Unknown auxiliary record type: " ++ [x]-{-# INLINE auxData #-}--queryAuxData :: (Char, Char) -> BAM -> Maybe AuxiliaryData-queryAuxData (x1,x2) bam = unsafePerformIO $- withForeignPtr (unbam bam) $ \b -> do- ptr <- bamAuxGet b [x1,x2]- if ptr == nullPtr- then return Nothing- else Just <$> getAuxData1 ptr-{-# INLINE queryAuxData #-}--getAuxData1 :: Ptr () -> IO AuxiliaryData-getAuxData1 ptr = castCCharToChar <$> peekByteOff ptr 0 >>= \case- 'A' -> AuxChar . castCCharToChar <$> peekByteOff ptr 1- 'c' -> AuxInt . fromIntegral <$> (peekByteOff ptr 1 :: IO Int8)- 'C' -> AuxInt . fromIntegral <$> (peekByteOff ptr 1 :: IO Word8)- 's' -> AuxInt . fromIntegral <$> (peekByteOff ptr 1 :: IO Int16)- 'S' -> AuxInt . fromIntegral <$> (peekByteOff ptr 1 :: IO Word16)- 'i' -> AuxInt . fromIntegral <$> (peekByteOff ptr 1 :: IO Int32)- 'I' -> AuxInt . fromIntegral <$> (peekByteOff ptr 1 :: IO Word32)- 'f' -> AuxFloat <$> peekByteOff ptr 1- 'Z' -> AuxString <$> B.packCString (plusPtr ptr 1)- 'H' -> AuxByteArray <$> B.packCString (plusPtr ptr 1)- 'B' -> do- n <- fromIntegral <$> (peekByteOff ptr 2 :: IO Int32)- peekByteOff ptr 1 >>= \case- 'c' -> AuxIntArray . map fromIntegral <$>- (peekArray n $ plusPtr ptr 6 :: IO [Int8])- 'C' -> AuxIntArray . map fromIntegral <$>- (peekArray n $ plusPtr ptr 6 :: IO [Word8])- 's' -> AuxIntArray . map fromIntegral <$>- (peekArray n $ plusPtr ptr 6 :: IO [Int16])- 'S' -> AuxIntArray . map fromIntegral <$>- (peekArray n $ plusPtr ptr 6 :: IO [Word16])- 'i' -> AuxIntArray . map fromIntegral <$>- (peekArray n $ plusPtr ptr 6 :: IO [Int32])- 'I' -> AuxIntArray . map fromIntegral <$>- (peekArray n $ plusPtr ptr 6 :: IO [Word32])- 'f' -> AuxFloatArray <$>- (peekArray n $ plusPtr ptr 6 :: IO [Float])- x -> error $ "Unknown auxiliary record type: " ++ [x]- x -> error $ "Unknown auxiliary record type: " ++ [x]-{-# INLINE getAuxData1 #-}- --- | Convert Bam record to Sam record.-bamToSam :: BAMHeader -> BAM -> SAM-bamToSam h b = SAM (seqName b) (flag b) (getChr h b) (startLoc b) (mapq b)- (cigar b) (mateChr h b) (mateStartLoc b) (tLen b) (getSeq b) (quality b)- (auxData b)---- | Template having multiple segments in sequencing-hasMultiSegments :: Word16 -> Bool-hasMultiSegments f = testBit f 0---- | Each segment properly aligned according to the aligner-isProperAligned :: Word16 -> Bool-isProperAligned f = testBit f 1---- | Segment unmapped-isUnmapped :: Word16 -> Bool-isUnmapped f = testBit f 2---- | Next segment in the template unmapped-isNextUnmapped :: Word16 -> Bool-isNextUnmapped f = testBit f 3---- | SEQ being reverse complemented-isRC :: Word16 -> Bool-isRC f = testBit f 4---- | SEQ of the next segment in the template being reverse complemented-isNextRC :: Word16 -> Bool-isNextRC f = testBit f 5---- | The first segment in the template-isFirstSegment :: Word16 -> Bool-isFirstSegment f = testBit f 6---- | The last segment in the template-isLastSegment :: Word16 -> Bool-isLastSegment f = testBit f 7---- | Secondary alignment-isSecondary :: Word16 -> Bool-isSecondary f = testBit f 8---- | Not passing filters, such as platform/vendor quality controls-isBadQual :: Word16 -> Bool-isBadQual f = testBit f 9---- | PCR or optical duplicate-isDup :: Word16 -> Bool-isDup f = testBit f 10---- | Supplementary alignment-isSupplementary :: Word16 -> Bool-isSupplementary f = testBit f 11-------------------------------------------------------------------------------------- Modify BAM ------------------------------------------------------------------------------------- | Append tag data to a bam record.-appendAux :: (Char, Char) -- ^ Tag- -> AuxiliaryData -- ^ Data- -> BAM- -> IO ()-appendAux (x1,x2) aux bam = withForeignPtr (unbam bam) $ \b -> do- code <- case aux of- AuxChar d -> with d $ bamAuxAppend b [x1,x2] 'A' 1- AuxInt d -> with d $ bamAuxAppend b [x1,x2] 'i' 4- AuxFloat d -> with d $ bamAuxAppend b [x1,x2] 'f' 4- AuxString d -> B.useAsCString d $- bamAuxAppend b [x1,x2] 'Z' (B.length d + 1) . castPtr- _ -> error "Not implemented"- if code == 0 then return () else error "Append aux failed"- where- with x fun = alloca $ \ptr -> poke ptr x >> fun (castPtr ptr)-{-# INLINE appendAux #-}
+ src/Bio/HTS.hs view
@@ -0,0 +1,9 @@+module Bio.HTS+ ( module Bio.HTS.Types+ , module Bio.HTS.BAM+ , module Bio.HTS.Utils+ ) where++import Bio.HTS.Types+import Bio.HTS.BAM+import Bio.HTS.Utils
+ src/Bio/HTS/BAM.chs view
@@ -0,0 +1,513 @@+{-# LANGUAGE ForeignFunctionInterface #-}+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE LambdaCase #-}++module Bio.HTS.BAM+ ( getBamHeader+ , getSortOrder+ , streamBam+ , sinkBam++ -- * Field accessors+ , refId+ , refName+ , startLoc+ , endLoc+ , readLen+ , isRev+ , flag+ , mapq+ , getSeq+ , queryName+ , qualityS+ , quality+ , sumQual+ , cigar+ , mateRefId+ , mateRefName+ , mateStartLoc+ , tLen+ , auxData+ , queryAuxData+ , bamToSam++ -- * Flag interpretation+ , hasMultiSegments+ , isProperAligned+ , isUnmapped+ , isNextUnmapped+ , isRC+ , isNextRC+ , isFirstSegment+ , isLastSegment+ , isSecondary+ , isBadQual+ , isDup+ , isSupplementary++ -- * Modify BAM+ , appendAux+ , setDup+ ) where++import Conduit+import Data.Bits (testBit)+import qualified Data.ByteString.Char8 as B+import qualified Data.ByteString as BS+import Data.Int+import Data.Word+import Foreign.C.String+import Foreign.ForeignPtr+import Foreign.Marshal.Array+import Foreign.Marshal.Alloc+import Foreign.Storable (poke, peekByteOff)+import Foreign.Ptr+import System.IO+import System.IO.Unsafe (unsafePerformIO)++import Bio.HTS.Internal+import Bio.HTS.Types++#include "htslib/sam.h"++getBamHeader :: FilePath -> IO BAMHeader+getBamHeader input = withHTSFile input ReadMode $ \hts ->+ getBgzf hts >>= bamHdrRead >>= newForeignPtr bamHdrDestroy >>=+ return . BAMHeader+{-# INLINE getBamHeader #-}+foreign import ccall unsafe "&bam_hdr_destroy"+ bamHdrDestroy :: FunPtr (Ptr BamHdr -> IO ())++-- | Get the sort information.+getSortOrder :: BAMHeader -> SortOrder+getSortOrder header = case lookup "SO" fields of+ Just "unknown" -> Unknown+ Just "unsorted" -> Unsorted+ Just "queryname" -> Queryname+ Just "coordinate" -> Coordinate+ _ -> Unknown+ where+ fields = map (f . B.split ':') $ tail $ B.split '\t' $ head $ B.lines $+ showBamHeader header+ f [a,b] = (a,b)+ f _ = error "Auxiliary field parsing failed!"++-- | Create a Bam stream from the file.+streamBam :: MonadResource m => FilePath -> ConduitT i BAM m ()+streamBam input = bracketP (htsOpen input "r") htsClose $ \hts -> do+ bgzf <- liftIO $ getBgzf hts+ _ <- liftIO $ bamHdrRead bgzf >>= newForeignPtr bamHdrDestroy+ loop bgzf+ where+ loop bgzf = do+ ptr <- liftIO $ callocBytes {# sizeof bam1_t #}+ code <- liftIO $ bamRead1 bgzf ptr+ case () of+ _ | code > 0 -> liftIO (newForeignPtr bamDestory1 ptr) >>=+ yield . BAM >> loop bgzf+ | code == -1 -> return ()+ | code == -2 -> error "truncated file"+ | otherwise -> error $ "read bam failed with code: " ++ show code+{-# INLINE streamBam #-}++foreign import ccall unsafe "&bam_destroy1"+ bamDestory1 :: FunPtr (Ptr Bam1 -> IO ())++-- | Write Bam records to a file.+sinkBam :: MonadResource m => FilePath -> BAMHeader -> ConduitT BAM o m ()+sinkBam output header = bracketP (htsOpen output "wb") htsClose $ \hts -> do+ bgzf <- liftIO $ getBgzf hts+ liftIO (withForeignPtr (unbamHeader header) $ bamHdrWrite bgzf) >>= \case+ 0 -> mapM_C $ \bam -> liftIO $ do+ code <- withForeignPtr (unbam bam) (bamWrite1 bgzf)+ if code < 0+ then error "bam_write1 failed"+ else return ()+ _ -> error "'bam_hdr_write' failed."+++-- | Reference sequence ID, −1 <= refId < n_ref; -1 for a read+-- without a mapping position.+refId :: BAM -> Int+refId = unsafePerformIO . flip withForeignPtr fun . unbam+ where+ fun = fmap fromIntegral . {#get bam1_t->core.tid #}+{-# INLINE refId #-}++-- | Return the reference sequence name (chromosome name) given the+-- bam file header.+refName :: BAMHeader -> BAM -> Maybe B.ByteString+refName header bam+ | chr < 0 = Nothing+ | otherwise = Just $ unsafePerformIO $+ withForeignPtr (unbamHeader header) $ \h ->+ bamChr h chr >>= B.packCString+ where+ chr = refId bam+{-# INLINE refName #-}++-- | Return the 0-based leftmost coordinate.+startLoc :: BAM -> Int+startLoc = unsafePerformIO . flip withForeignPtr fun . unbam+ where+ fun = fmap fromIntegral . {#get bam1_t->core.pos #}+{-# INLINE startLoc #-}++-- | Return the end location according to a 0-based coordinate system.+-- For a mapped read, this is just position + cigar2rlen.+-- For an unmapped read (either according to its flags or if it has no cigar+-- string), we return position + 1 by convention.+endLoc :: BAM -> Int+endLoc = unsafePerformIO . flip withForeignPtr bamEndpos .unbam+{-# INLINE endLoc #-}++-- | Return the query length (read length).+readLen :: BAM -> Int+readLen = unsafePerformIO . flip withForeignPtr fun . unbam+ where+ fun = fmap fromIntegral . {#get bam1_t->core.l_qseq #}+{-# INLINE readLen #-}++-- | Whether the query is on the reverse strand.+isRev :: BAM -> Bool+isRev = unsafePerformIO . flip withForeignPtr bamIsRev . unbam+{-# INLINE isRev #-}++-- | Return the flag.+flag :: BAM -> Word16+flag = unsafePerformIO . flip withForeignPtr fn . unbam+ where+ fn = fmap fromIntegral . {#get bam1_t->core.flag #}+{-# INLINE flag #-}++-- | MAPping Quality. It equals −10 log10 Pr{mapping position is wrong},+-- rounded to the nearest integer. A value 255 indicates that the+-- mapping quality is not available.+mapq :: BAM -> Word8+mapq = unsafePerformIO . flip withForeignPtr fn . unbam+ where+ fn = fmap fromIntegral . {#get bam1_t->core.qual #}+{-# INLINE mapq #-}++-- | Return the DNA sequence.+getSeq :: BAM -> Maybe B.ByteString+getSeq = unsafePerformIO . flip withForeignPtr fn . unbam+ where+ fn b = fromIntegral <$> {#get bam1_t->core.l_qseq #} b >>= \case+ 0 -> return Nothing+ n -> allocaArray n $ \str -> do+ bamGetSeq b str n+ Just <$> B.packCStringLen (str, n)+{-# INLINE getSeq #-}++-- | Get the name of the query.+queryName :: BAM -> B.ByteString+queryName = unsafePerformIO . flip withForeignPtr fn . unbam+ where+ fn b = {#get bam1_t->data #} b >>= B.packCString . castPtr+{-# INLINE queryName #-}++-- | Human readable quality score which is: Phred base quality + 33.+qualityS :: BAM -> Maybe B.ByteString+qualityS = fmap (BS.map (+33)) . quality++-- | Phred base quality (a sequence of 0xFF if absent).+quality :: BAM -> Maybe B.ByteString+quality = unsafePerformIO . flip withForeignPtr fn . unbam+ where+ fn b = fromIntegral <$> {#get bam1_t->core.l_qseq #} b >>= \case+ 0 -> return Nothing+ n -> allocaArray n $ \str -> bamGetQual b str n >>= \case+ 0 -> Just <$> B.packCStringLen (str, n)+ _ -> return Nothing+{-# INLINE quality #-}++-- | Sum of quality scores that above certain threshold (e.g., 15).+sumQual :: Int -- ^ Threshold+ -> BAM+ -> Maybe Int+sumQual th bam = BS.foldl' f 0 <$> quality bam + where+ f acc x | fromIntegral x < th = acc+ | otherwise = fromIntegral x + acc+{-# INLINE sumQual #-}++-- | Get the CIGAR string.+cigar :: BAM -> Maybe CIGAR+cigar = unsafePerformIO . flip withForeignPtr fn . unbam+ where+ fn b = fromIntegral <$> {#get bam1_t->core.n_cigar #} b >>= \case+ 0 -> return Nothing+ n -> allocaArray n $ \num -> allocaArray n $ \str -> do+ bamGetCigar b num str n+ num' <- peekArray (fromIntegral n) num+ str' <- peekArray (fromIntegral n) str+ return $ Just $ CIGAR $ zip (map fromIntegral num') $+ map castCCharToChar str'+{-# INLINE cigar #-}++-- | Ref-ID of the next segment (the paired read).+mateRefId :: BAM -> Int+mateRefId = unsafePerformIO . flip withForeignPtr fn . unbam+ where+ fn = fmap fromIntegral . {#get bam1_t->core.mtid #}+{-# INLINE mateRefId #-}++mateRefName :: BAMHeader -> BAM -> Maybe B.ByteString+mateRefName header bam+ | chr < 0 = Nothing+ | otherwise = Just $ unsafePerformIO $+ withForeignPtr (unbamHeader header) $ \h ->+ bamChr h chr >>= B.packCString+ where+ chr = mateRefId bam+{-# INLINE mateRefName #-}++-- | 0-based+mateStartLoc :: BAM -> Int+mateStartLoc = unsafePerformIO . flip withForeignPtr fn . unbam+ where+ fn = fmap fromIntegral . {#get bam1_t->core.mpos #}+{-# INLINE mateStartLoc #-}++-- | Observed Template length.+-- If all segments are mapped to the same reference, the unsigned+-- observed template length equals the number of bases from the+-- leftmost mapped base to the rightmost mapped base. +tLen :: BAM -> Int+tLen = unsafePerformIO . flip withForeignPtr fn . unbam+ where+ fn = fmap fromIntegral . {#get bam1_t->core.isize #}+{-# INLINE tLen #-}++auxData :: BAM -> [((Char, Char), AuxiliaryData)]+auxData bam = unsafePerformIO $ withForeignPtr (unbam bam) $ \b -> do+ l <- bamGetLAux b+ aux <- bamGetAux b + go aux l+ where+ go ptr i+ | i <= 0 = return []+ | otherwise = do+ name <- (,) <$> peekByteOff ptr 0 <*> peekByteOff ptr 1+ castCCharToChar <$> peekByteOff ptr 2 >>= \case+ 'A' -> do+ r <- AuxChar . castCCharToChar <$> peekByteOff ptr 3+ rs <- go (plusPtr ptr 4) $ i - 4+ return $ (name, r) : rs+ -- int8_t+ 'c' -> do+ r <- AuxInt . fromIntegral <$> (peekByteOff ptr 3 :: IO Int8)+ rs <- go (plusPtr ptr 4) $ i - 4+ return $ (name, r) : rs+ -- uint8_t+ 'C' -> do+ r <- AuxInt . fromIntegral <$> (peekByteOff ptr 3 :: IO Word8)+ rs <- go (plusPtr ptr 4) $ i - 4+ return $ (name, r) : rs+ -- int16_t+ 's' -> do+ r <- AuxInt . fromIntegral <$> (peekByteOff ptr 3 :: IO Int16)+ rs <- go (plusPtr ptr 5) $ i - 5+ return $ (name, r) : rs+ -- uint16_t+ 'S' -> do+ r <- AuxInt . fromIntegral <$> (peekByteOff ptr 3 :: IO Word16)+ rs <- go (plusPtr ptr 5) $ i - 5+ return $ (name, r) : rs+ -- int32_t+ 'i' -> do+ r <- AuxInt . fromIntegral <$> (peekByteOff ptr 3 :: IO Int32)+ rs <- go (plusPtr ptr 7) $ i - 7+ return $ (name, r) : rs+ -- uint32_t+ 'I' -> do+ r <- AuxInt . fromIntegral <$> (peekByteOff ptr 3 :: IO Word32)+ rs <- go (plusPtr ptr 7) $ i - 7+ return $ (name, r) : rs+ 'f' -> do+ r <- AuxFloat <$> peekByteOff ptr 3+ rs <- go (plusPtr ptr 7) $ i - 7+ return $ (name, r) : rs+ 'Z' -> do+ str <- B.packCString (plusPtr ptr 3)+ let l = B.length str + 1 + 3+ rs <- go (plusPtr ptr l) $ i - l+ return $ (name, AuxString str) : rs+ 'H' -> do + str <- B.packCString (plusPtr ptr 3)+ let l = B.length str + 1 + 3+ rs <- go (plusPtr ptr l) $ i - l+ return $ (name, AuxByteArray str) : rs+ 'B' -> do+ n <- fromIntegral <$> (peekByteOff ptr 4 :: IO Int32)+ castCCharToChar <$> peekByteOff ptr 3 >>= \case+ 'c' -> do+ r <- AuxIntArray . map fromIntegral <$>+ (peekArray n $ plusPtr ptr 8 :: IO [Int8])+ let l = 8 + n+ rs <- go (plusPtr ptr l) $ i - l+ return $ (name, r) : rs+ 'C' -> do+ r <- AuxIntArray . map fromIntegral <$>+ (peekArray n $ plusPtr ptr 8 :: IO [Word8])+ let l = 8 + n+ rs <- go (plusPtr ptr l) $ i - l+ return $ (name, r) : rs+ 's' -> do+ r <- AuxIntArray . map fromIntegral <$>+ (peekArray n $ plusPtr ptr 8 :: IO [Int16])+ let l = 8 + 2 * n+ rs <- go (plusPtr ptr l) $ i - l+ return $ (name, r) : rs+ 'S' -> do+ r <- AuxIntArray . map fromIntegral <$>+ (peekArray n $ plusPtr ptr 8 :: IO [Word16])+ let l = 8 + 2 * n+ rs <- go (plusPtr ptr l) $ i - l+ return $ (name, r) : rs+ 'i' -> do+ r <- AuxIntArray . map fromIntegral <$>+ (peekArray n $ plusPtr ptr 8 :: IO [Int32])+ let l = 8 + 4 * n+ rs <- go (plusPtr ptr l) $ i - l+ return $ (name, r) : rs+ 'I' -> do+ r <- AuxIntArray . map fromIntegral <$>+ (peekArray n $ plusPtr ptr 8 :: IO [Word32])+ let l = 8 + 4 * n+ rs <- go (plusPtr ptr l) $ i - l+ return $ (name, r) : rs+ 'f' -> do+ r <- AuxFloatArray <$>+ (peekArray n $ plusPtr ptr 8 :: IO [Float])+ let l = 8 + 4 * n+ rs <- go (plusPtr ptr l) $ i - l+ return $ (name, r) : rs+ x -> error $ "Unknown auxiliary record type: " ++ [x]+ x -> error $ "Unknown auxiliary record type: " ++ [x]+{-# INLINE auxData #-}++queryAuxData :: (Char, Char) -> BAM -> Maybe AuxiliaryData+queryAuxData (x1,x2) bam = unsafePerformIO $+ withForeignPtr (unbam bam) $ \b -> do+ ptr <- bamAuxGet b [x1,x2]+ if ptr == nullPtr+ then return Nothing+ else Just <$> getAuxData1 ptr+{-# INLINE queryAuxData #-}++getAuxData1 :: Ptr () -> IO AuxiliaryData+getAuxData1 ptr = castCCharToChar <$> peekByteOff ptr 0 >>= \case+ 'A' -> AuxChar . castCCharToChar <$> peekByteOff ptr 1+ 'c' -> AuxInt . fromIntegral <$> (peekByteOff ptr 1 :: IO Int8)+ 'C' -> AuxInt . fromIntegral <$> (peekByteOff ptr 1 :: IO Word8)+ 's' -> AuxInt . fromIntegral <$> (peekByteOff ptr 1 :: IO Int16)+ 'S' -> AuxInt . fromIntegral <$> (peekByteOff ptr 1 :: IO Word16)+ 'i' -> AuxInt . fromIntegral <$> (peekByteOff ptr 1 :: IO Int32)+ 'I' -> AuxInt . fromIntegral <$> (peekByteOff ptr 1 :: IO Word32)+ 'f' -> AuxFloat <$> peekByteOff ptr 1+ 'Z' -> AuxString <$> B.packCString (plusPtr ptr 1)+ 'H' -> AuxByteArray <$> B.packCString (plusPtr ptr 1)+ 'B' -> do+ n <- fromIntegral <$> (peekByteOff ptr 2 :: IO Int32)+ peekByteOff ptr 1 >>= \case+ 'c' -> AuxIntArray . map fromIntegral <$>+ (peekArray n $ plusPtr ptr 6 :: IO [Int8])+ 'C' -> AuxIntArray . map fromIntegral <$>+ (peekArray n $ plusPtr ptr 6 :: IO [Word8])+ 's' -> AuxIntArray . map fromIntegral <$>+ (peekArray n $ plusPtr ptr 6 :: IO [Int16])+ 'S' -> AuxIntArray . map fromIntegral <$>+ (peekArray n $ plusPtr ptr 6 :: IO [Word16])+ 'i' -> AuxIntArray . map fromIntegral <$>+ (peekArray n $ plusPtr ptr 6 :: IO [Int32])+ 'I' -> AuxIntArray . map fromIntegral <$>+ (peekArray n $ plusPtr ptr 6 :: IO [Word32])+ 'f' -> AuxFloatArray <$>+ (peekArray n $ plusPtr ptr 6 :: IO [Float])+ x -> error $ "Unknown auxiliary record type: " ++ [x]+ x -> error $ "Unknown auxiliary record type: " ++ [x]+{-# INLINE getAuxData1 #-}+ +-- | Convert Bam record to Sam record.+bamToSam :: BAMHeader -> BAM -> SAM+bamToSam h b = SAM (queryName b) (flag b) (refName h b) (startLoc b + 1) (mapq b)+ (cigar b) (mateRefName h b) (mateStartLoc b + 1) (tLen b) (getSeq b) (quality b)+ (auxData b)++-- | Template having multiple segments in sequencing+hasMultiSegments :: Word16 -> Bool+hasMultiSegments f = testBit f 0++-- | Each segment properly aligned according to the aligner+isProperAligned :: Word16 -> Bool+isProperAligned f = testBit f 1++-- | Segment unmapped+isUnmapped :: Word16 -> Bool+isUnmapped f = testBit f 2++-- | Next segment in the template unmapped+isNextUnmapped :: Word16 -> Bool+isNextUnmapped f = testBit f 3++-- | SEQ being reverse complemented+isRC :: Word16 -> Bool+isRC f = testBit f 4++-- | SEQ of the next segment in the template being reverse complemented+isNextRC :: Word16 -> Bool+isNextRC f = testBit f 5++-- | The first segment in the template+isFirstSegment :: Word16 -> Bool+isFirstSegment f = testBit f 6++-- | The last segment in the template+isLastSegment :: Word16 -> Bool+isLastSegment f = testBit f 7++-- | Secondary alignment+isSecondary :: Word16 -> Bool+isSecondary f = testBit f 8++-- | Not passing filters, such as platform/vendor quality controls+isBadQual :: Word16 -> Bool+isBadQual f = testBit f 9++-- | PCR or optical duplicate+isDup :: Word16 -> Bool+isDup f = testBit f 10++-- | Supplementary alignment+isSupplementary :: Word16 -> Bool+isSupplementary f = testBit f 11+++--------------------------------------------------------------------------------+-- Modify BAM +--------------------------------------------------------------------------------++-- | Append tag data to a bam record.+appendAux :: (Char, Char) -- ^ Tag+ -> AuxiliaryData -- ^ Data+ -> BAM+ -> IO ()+appendAux (x1,x2) aux bam = withForeignPtr (unbam bam) $ \b -> do+ code <- case aux of+ AuxChar d -> with d $ bamAuxAppend b [x1,x2] 'A' 1+ AuxInt d -> with d $ bamAuxAppend b [x1,x2] 'i' 4+ AuxFloat d -> with d $ bamAuxAppend b [x1,x2] 'f' 4+ AuxString d -> B.useAsCString d $+ bamAuxAppend b [x1,x2] 'Z' (B.length d + 1) . castPtr+ _ -> error "Not implemented"+ if code == 0 then return () else error "Append aux failed"+ where+ with x fun = alloca $ \ptr -> poke ptr x >> fun (castPtr ptr)+{-# INLINE appendAux #-}++-- | Turn on duplicate flag.+setDup :: BAM -> IO ()+setDup bam = withForeignPtr (unbam bam) bamMarkDup+{-# INLINE setDup #-}
src/Bio/HTS/Internal.chs view
@@ -30,6 +30,7 @@ , bamGetLAux , bamAuxGet , bamAuxAppend+ , bamMarkDup ) where import Control.Exception (bracket)@@ -51,6 +52,7 @@ void bam_get_cigar_(bam1_t *b, int *num, char *str, uint32_t l); uint8_t* bam_get_aux_(bam1_t *b); int bam_get_l_aux_(bam1_t *b);+void bam_mark_dup(bam1_t *b); #endc -- | Opaque data representing the hts file.@@ -124,6 +126,9 @@ {#fun bam_aux_append as ^ { castPtr `Ptr Bam1', `String', `Char', `Int', castPtr `Ptr ()' } -> `CInt' #}++{#fun bam_mark_dup as ^ { castPtr `Ptr Bam1' } -> `()' #}+ -------------------------------------------------------------------------------- -- BAM Header
src/Bio/HTS/Types.hs view
@@ -7,6 +7,9 @@ , showBamHeader , SortOrder(..) , SAM(..)+ , CIGAR(..)+ , cigar2String+ , string2Cigar , AuxiliaryData(..) , showSam ) where@@ -40,14 +43,29 @@ | Coordinate deriving (Show, Eq) --- | The SAM format.+newtype CIGAR = CIGAR [(Int, Char)]++cigar2String :: CIGAR -> B.ByteString+cigar2String (CIGAR c) = B.concat $+ concatMap (\(i, x) -> [fromJust $ packDecimal i, B.singleton x]) c++string2Cigar :: B.ByteString -> CIGAR+string2Cigar c = CIGAR $ go c+ where+ go x = case readDecimal x of+ Nothing -> error $ "parse cigar fail: " ++ show c+ Just (n, remain) -> if B.length remain == 1+ then [(n, B.head remain)]+ else (n, B.head remain) : go (B.tail remain)+ +-- | The SAM format. The SAM format uses 1-based coordinate system. data SAM = SAM { _sam_qname :: !B.ByteString , _sam_flag :: !Word16 , _sam_rname :: !(Maybe B.ByteString) , _sam_pos :: !Int , _sam_mapq :: !Word8- , _sam_cigar :: !(Maybe [(Int, Char)])+ , _sam_cigar :: !(Maybe CIGAR) , _sam_rnext :: !(Maybe B.ByteString) , _sam_pnext :: !Int , _sam_tlen :: !Int@@ -58,12 +76,12 @@ showSam :: SAM -> B.ByteString showSam SAM{..} = B.intercalate "\t" $- [ _sam_qname, pack' _sam_flag, fromMaybe "*" _sam_rname, pack' $ _sam_pos + 1- , pack' _sam_mapq, fromMaybe "*" $ f <$> _sam_cigar, fromMaybe "*" _sam_rnext- , pack' $ _sam_pnext + 1, pack' _sam_tlen, fromMaybe "*" _sam_seq- , fromMaybe "*" $ BS.map (+33) <$> _sam_qual ] ++ map showAuxiliaryData _sam_aux+ [ _sam_qname, pack' _sam_flag, fromMaybe "*" _sam_rname, pack' _sam_pos+ , pack' _sam_mapq, fromMaybe "*" $ cigar2String <$> _sam_cigar+ , fromMaybe "*" _sam_rnext, pack' _sam_pnext, pack' _sam_tlen+ , fromMaybe "*" _sam_seq, fromMaybe "*" $ BS.map (+33) <$> _sam_qual+ ] ++ map showAuxiliaryData _sam_aux where- f = B.concat . concatMap (\(i, x) -> [pack' i, B.singleton x]) pack' :: Integral a => a -> B.ByteString pack' = fromJust . packDecimal
+ src/Bio/HTS/Utils.hs view
@@ -0,0 +1,146 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE LambdaCase #-}+{-# LANGUAGE RecordWildCards #-}+{-# LANGUAGE PartialTypeSignatures #-}+module Bio.HTS.Utils+ ( markDupBy+ , makeKey+ , BAMKey+ , Orientation+ ) where++import Conduit+import qualified Data.Map.Strict as M+import Data.List+import Data.Maybe+import qualified Data.Sequence as S+import qualified Data.ByteString.Char8 as B++import Bio.HTS.BAM+import Bio.HTS.Types++data Orientation = F | R | FF | FR | RR | RF deriving (Eq, Ord, Show)++data BAMKey = Single { _ref_id1 :: Int+ , _loc1 :: Int+ , _orientation :: Orientation+ , _barcode :: Maybe B.ByteString }+ | Pair { _ref_id1 :: Int+ , _ref_id2 :: Int+ , _loc1 :: Int+ , _loc2 :: Int+ , _orientation :: Orientation+ , _leftmost :: Bool+ , _barcode :: Maybe B.ByteString }+ deriving (Eq, Ord, Show)++makeKey :: (BAM -> Maybe B.ByteString) -- ^ Get Barcode+ -> BAM+ -> (BAMKey, BAMKey)+makeKey fn bam = (single, pair)+ where+ single = Single ref1 (if isFwd1 then lloc1 else rloc1)+ (if isFwd1 then F else R) bc+ pair = Pair ref1 ref2 loc1 loc2 orientation isLeftMost bc+ bc = fn bam+ ref1 = refId bam+ lloc1 = startLoc bam - fst clipped1 + 1+ rloc1 = endLoc bam + snd clipped1+ clipped1 = getClipped $ fromJust $ cigar bam+ isFwd1 = not $ isRC flg+ ref2 = mateRefId bam+ lloc2 = mateStartLoc bam - fst clipped2 + 1+ rloc2 = mateStartLoc bam + ciglen cig + snd clipped2+ clipped2 = getClipped cig+ isFwd2 = not $ isNextRC flg+ cig = case queryAuxData ('M', 'C') bam of+ Just (AuxString x) -> string2Cigar x+ _ -> error "No MC tag. Please run samtools fixmate on file first."+ 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 bam+ ciglen (CIGAR c) = foldl' f 0 c+ where f acc (n,x) = if x `elem` "MDN=X" then n + acc else acc++-- | Get the number of clips at both ends.+getClipped :: CIGAR -> (Int, Int)+getClipped (CIGAR c) | length c <= 1 = (0,0)+ | otherwise = (getC $ head c, getC $ last c)+ where+ getC x = case x of+ (n, 'S') -> n+ (n, 'H') -> n+ _ -> 0+{-# INLINE getClipped #-}+++-- | 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)+-- 4. Barcode+--+-- 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+ -> ConduitT BAM BAM m ()+markDupBy bcFn = go (-1,-1) M.empty S.empty+ where+ go prev keyMap readBuf = await >>= \case+ Nothing -> mapM_ (\(_,_,x) -> yield x) readBuf+ Just bam -> markDup prev keyMap readBuf bam >>=+ (\(a,b,c) -> go a b c)+ markDup prev keyMap readBuf bam+ | isUnmapped flg = yield bam >> return (prev, keyMap, readBuf)+ | not (isSorted prev cur) = error $+ "bad coordinate order: " ++ show (prev, cur)+ | isDup flg = return (cur, keyMap, readBuf')+ | otherwise = do+ keyMap' <- addAndMark key bam keyMap+ update keyMap' readBuf' cur+ where+ readBuf' = readBuf S.|> (key, _loc1 singleKey, bam)+ cur = (refId bam, startLoc bam)+ key | not (hasMultiSegments flg) || isNextUnmapped flg = singleKey+ | otherwise = pairKey+ (singleKey, pairKey) = makeKey bcFn bam+ flg = flag bam+ update keyMap readBuf (chr, loc) = mapM_ (\(_,_,b) -> yield b) exclude >>+ return ((chr, loc), keyMap', kept)+ where+ keyMap' = foldl' (\m (k,_,_) -> M.delete k m) keyMap exclude+ (exclude, kept) = S.breakl toKeep readBuf+ toKeep (key, pos,_) = _ref_id1 key == chr && pos + max_len > loc+ addAndMark key bam keyMap = liftIO $ case M.lookup key keyMap of+ Nothing -> return $ M.insert key (bam, sc) keyMap+ Just (old, old_sc) -> if sc > old_sc -- || (sc == old_sc && queryName bam < queryName old)+ then setDup old >> return (M.insert key (bam, sc) keyMap)+ else setDup bam >> return keyMap+ where+ sc = case queryAuxData ('m', 's') bam of+ Just (AuxInt x) -> x + fromJust (sumQual 15 bam)+ _ -> fromJust $ sumQual 15 bam+ max_len = 500+ isSorted (chr1, loc1) (chr2, loc2) = chr1 < chr2 ||+ (chr1 == chr2 && loc1 <= loc2)+{-# INLINE markDupBy #-}
tests/Main.hs view
@@ -1,19 +1,18 @@ {-# LANGUAGE LambdaCase #-} {-# LANGUAGE FlexibleContexts #-} import Test.Tasty-import Test.Tasty.Golden+import Test.Tasty.Golden import Conduit-import Control.Monad.Reader import Bio.HTS-import Bio.HTS.Types main :: IO () main = defaultMain $ testGroup "Main"- [ tests ]+ [ bamToSamTest+ , markDupTest ] -tests :: TestTree-tests = goldenVsFile "BAM Read/Write Test" expect output $+bamToSamTest :: TestTree+bamToSamTest = goldenVsFile "BAM Read/Write Test" expect output $ bamFileToSamFile input output where input = "tests/data/example.bam"@@ -27,3 +26,17 @@ header <- getBamHeader input runResourceT $ runConduit $ streamBam input .| mapC (showSam . bamToSam header) .| unlinesAsciiC .| sinkFile output++markDupTest :: TestTree+markDupTest = testGroup "Mark duplicates"+ [ markDup "Single end" "tests/data/single_end_dedup.bam" + "tests/data/single_end.bam" "dedup.bam" + , markDup "Paired end" "tests/data/paired_end_dedup.bam" + "tests/data/paired_end.bam" "dedup.bam" + ]+ where+ markDup name expect input output = goldenVsFile name expect output $ do+ header <- getBamHeader input+ runResourceT $ runConduit $ streamBam input .|+ markDupBy (const Nothing) .|+ filterC (not . isDup . flag) .| sinkBam output header
+ tests/data/single_end.bam view
binary file changed (absent → 1096842 bytes)
+ tests/data/single_end_dedup.bam view
binary file changed (absent → 849691 bytes)