diff --git a/ChangeLog.md b/ChangeLog.md --- a/ChangeLog.md +++ b/ChangeLog.md @@ -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. diff --git a/HsHTSLib.cabal b/HsHTSLib.cabal --- a/HsHTSLib.cabal +++ b/HsHTSLib.cabal @@ -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 diff --git a/cbits/hs_htslib.c b/cbits/hs_htslib.c --- a/cbits/hs_htslib.c +++ b/cbits/hs_htslib.c @@ -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; } + + + diff --git a/src/Bio/HTS.chs b/src/Bio/HTS.chs deleted file mode 100644 --- a/src/Bio/HTS.chs +++ /dev/null @@ -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 #-} diff --git a/src/Bio/HTS.hs b/src/Bio/HTS.hs new file mode 100644 --- /dev/null +++ b/src/Bio/HTS.hs @@ -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 diff --git a/src/Bio/HTS/BAM.chs b/src/Bio/HTS/BAM.chs new file mode 100644 --- /dev/null +++ b/src/Bio/HTS/BAM.chs @@ -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 #-} diff --git a/src/Bio/HTS/Internal.chs b/src/Bio/HTS/Internal.chs --- a/src/Bio/HTS/Internal.chs +++ b/src/Bio/HTS/Internal.chs @@ -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 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 @@ -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 diff --git a/src/Bio/HTS/Utils.hs b/src/Bio/HTS/Utils.hs new file mode 100644 --- /dev/null +++ b/src/Bio/HTS/Utils.hs @@ -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 #-} diff --git a/tests/Main.hs b/tests/Main.hs --- a/tests/Main.hs +++ b/tests/Main.hs @@ -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 diff --git a/tests/data/single_end.bam b/tests/data/single_end.bam new file mode 100644 Binary files /dev/null and b/tests/data/single_end.bam differ diff --git a/tests/data/single_end_dedup.bam b/tests/data/single_end_dedup.bam new file mode 100644 Binary files /dev/null and b/tests/data/single_end_dedup.bam differ