module Bio.SamTools.Bam (
HeaderSeq(..)
, Header, nTargets, targetSeqList, targetSeq, targetSeqName, targetSeqLen, lookupTarget
, Bam1
, targetID, targetName, targetLen, position
, isPaired, isProperPair, isUnmap, isMateUnmap, isReverse, isMateReverse
, isRead1, isRead2, isSecondary, isQCFail, isDup
, cigars, queryName, queryLength, querySeq, queryQual
, mateTargetID, mateTargetName, mateTargetLen, matePosition, insertSize
, nMismatch, nHits, matchDesc
, refSpLoc, refSeqLoc
, InHandle, inHeader
, openTamInFile, openTamInFileWithIndex, openBamInFile
, closeInHandle
, withTamInFile, withTamInFileWithIndex, withBamInFile
, get1
, OutHandle, outHeader
, openTamOutFile, openBamOutFile
, closeOutHandle
, withTamOutFile, withBamOutFile
, put1
)
where
import Control.Concurrent.MVar
import Control.Exception
import Control.Monad
import Data.Bits
import qualified Data.ByteString as BSW
import qualified Data.ByteString.Char8 as BS
import Foreign hiding (new)
import Foreign.C.Types
import Foreign.C.String
import qualified Data.Vector as V
import Bio.SeqLoc.OnSeq
import qualified Bio.SeqLoc.SpliceLocation as SpLoc
import Bio.SeqLoc.Strand
import Bio.SamTools.Cigar
import Bio.SamTools.Internal
import Bio.SamTools.LowLevel
targetID :: Bam1 -> Maybe Int
targetID b = unsafePerformIO $ withForeignPtr (ptrBam1 b) $ liftM fromTID . getTID
where fromTID ctid | ctid < 0 = Nothing
| otherwise = Just $! fromIntegral ctid
targetName :: Bam1 -> Maybe BS.ByteString
targetName b = liftM (targetSeqName (header b)) $! targetID b
targetLen :: Bam1 -> Maybe Int64
targetLen b = liftM (targetSeqLen (header b)) $! targetID b
position :: Bam1 -> Maybe Int64
position b = unsafePerformIO $ withForeignPtr (ptrBam1 b) $ liftM fromPos . getPos
where fromPos cpos | cpos < 0 = Nothing
| otherwise = Just $! fromIntegral cpos
isFlagSet :: BamFlag -> Bam1 -> Bool
isFlagSet f b = unsafePerformIO $ withForeignPtr (ptrBam1 b) $ liftM isfset . getFlag
where isfset = (== f) . (.&. f)
isPaired :: Bam1 -> Bool
isPaired = isFlagSet flagPaired
isProperPair :: Bam1 -> Bool
isProperPair = isFlagSet flagProperPair
isUnmap :: Bam1 -> Bool
isUnmap = isFlagSet flagUnmap
isMateUnmap :: Bam1 -> Bool
isMateUnmap = isFlagSet flagMUnmap
isReverse :: Bam1 -> Bool
isReverse = isFlagSet flagReverse
isMateReverse :: Bam1 -> Bool
isMateReverse = isFlagSet flagMReverse
isRead1 :: Bam1 -> Bool
isRead1 = isFlagSet flagRead1
isRead2 :: Bam1 -> Bool
isRead2 = isFlagSet flagRead2
isSecondary :: Bam1 -> Bool
isSecondary = isFlagSet flagSecondary
isQCFail :: Bam1 -> Bool
isQCFail = isFlagSet flagQCFail
isDup :: Bam1 -> Bool
isDup = isFlagSet flagDup
cigars :: Bam1 -> [Cigar]
cigars b = unsafePerformIO $ withForeignPtr (ptrBam1 b) $ \p -> do
nc <- getNCigar p
liftM (map toCigar) $! peekArray nc . bam1Cigar $ p
queryName :: Bam1 -> BS.ByteString
queryName b = unsafePerformIO $ withForeignPtr (ptrBam1 b) (return . bam1QName)
queryLength :: Bam1 -> Maybe Int64
queryLength b = unsafePerformIO $ withForeignPtr (ptrBam1 b) $ liftM fromc . getLQSeq
where fromc clq | clq < 1 = Nothing
| otherwise = Just $! fromIntegral clq
querySeq :: Bam1 -> Maybe BS.ByteString
querySeq b = unsafePerformIO $ withForeignPtr (ptrBam1 b) $ \p ->
let seqarr = bam1Seq p
getQSeq l | l < 1 = return Nothing
| otherwise = return $! Just $!
fst (BS.unfoldrN (fromIntegral l) (\i -> if i==l then Nothing else Just (seqiToChar $ bam1Seqi seqarr $ i,i+1)) 0)
in getLQSeq p >>= getQSeq
queryQual :: Bam1 -> Maybe BS.ByteString
queryQual b = unsafePerformIO $ withForeignPtr (ptrBam1 b) $ \p ->
let getQQual l | l < 1 = return Nothing
| otherwise = do q0 <- peek $! bam1Qual p
if q0 == 0xff
then return Nothing
else liftM (Just . BSW.map (+ 33)) $!
BS.packCStringLen (castPtr $! bam1Qual p, fromIntegral l)
in getLQSeq p >>= getQQual
seqiToChar :: CUChar -> Char
seqiToChar = (chars V.!) . fromIntegral
where chars = emptyChars V.// [(1, 'A'), (2, 'C'), (4, 'G'), (8, 'T'), (15, 'N')]
emptyChars = V.generate 16 (\idx -> error $ "Unknown char " ++ show idx)
mateTargetID :: Bam1 -> Maybe Int
mateTargetID b = unsafePerformIO $ withForeignPtr (ptrBam1 b) $ liftM fromTID . getMTID
where fromTID ctid | ctid < 0 = Nothing
| otherwise = Just $! fromIntegral ctid
mateTargetName :: Bam1 -> Maybe BS.ByteString
mateTargetName b = liftM (targetSeqName (header b)) $! mateTargetID b
mateTargetLen :: Bam1 -> Maybe Int64
mateTargetLen b = liftM (targetSeqLen (header b)) $! mateTargetID b
matePosition :: Bam1 -> Maybe Int64
matePosition b = unsafePerformIO $ withForeignPtr (ptrBam1 b) $ liftM fromPos . getMPos
where fromPos cpos | cpos < 0 = Nothing
| otherwise = Just $! fromIntegral cpos
insertSize :: Bam1 -> Maybe Int64
insertSize b = unsafePerformIO $ withForeignPtr (ptrBam1 b) $ liftM fromISize . getISize
where fromISize cis | cis < 1 = Nothing
| otherwise = Just $! fromIntegral cis
matchDesc :: Bam1 -> Maybe BS.ByteString
matchDesc b = unsafePerformIO $ withForeignPtr (ptrBam1 b) $ \p ->
withCAString "MD" $ \mdstr ->
do md <- bamAuxGet p mdstr
if md == nullPtr
then return Nothing
else do cstr <- bamAux2Z md
if cstr == nullPtr
then return Nothing
else liftM Just . BS.packCString $ cstr
nHits :: Bam1 -> Maybe Int
nHits b = unsafePerformIO $ withForeignPtr (ptrBam1 b) $ \p ->
withCAString "NH" $ \nhstr ->
do nh <- bamAuxGet p nhstr
if nh == nullPtr
then return Nothing
else liftM Just $! liftM fromIntegral $! bamAux2i nh
nMismatch :: Bam1 -> Maybe Int
nMismatch b = unsafePerformIO $ withForeignPtr (ptrBam1 b) $ \p ->
withCAString "NM" $ \nmstr ->
do nm <- bamAuxGet p nmstr
if nm == nullPtr
then return Nothing
else liftM Just $! liftM fromIntegral $! bamAux2i nm
refSpLoc :: Bam1 -> Maybe SpLoc.SpliceLoc
refSpLoc b | isUnmap b = Nothing
| otherwise = liftM (stranded strand) $! liftM2 (cigarToSpLoc) (position b) (Just . cigars $ b)
where strand = if isReverse b then Minus else Plus
refSeqLoc :: Bam1 -> Maybe SpliceSeqLoc
refSeqLoc b = liftM2 OnSeq (liftM toSeqLabel $! targetName b) (refSpLoc b)
data InHandle = InHandle { inFilename :: !FilePath
, samfile :: !(MVar (Ptr SamFileInt))
, inHeader :: !Header
}
newInHandle :: FilePath -> Ptr SamFileInt -> IO InHandle
newInHandle filename fsam = do
when (fsam == nullPtr) $ ioError . userError $ "Error opening BAM file " ++ show filename
mv <- newMVar fsam
addMVarFinalizer mv (finalizeSamFile mv)
bhdr <- getSbamHeader fsam
when (bhdr == nullPtr) $ ioError . userError $
"Error reading header from BAM file " ++ show filename
hdr <- newHeader bhdr
return $ InHandle { inFilename = filename, samfile = mv, inHeader = hdr }
openTamInFile :: FilePath -> IO InHandle
openTamInFile filename = sbamOpen filename "r" nullPtr >>= newInHandle filename
openTamInFileWithIndex :: FilePath -> FilePath -> IO InHandle
openTamInFileWithIndex filename indexname
= withCString indexname (sbamOpen filename "r" . castPtr) >>= newInHandle filename
openBamInFile :: FilePath -> IO InHandle
openBamInFile filename = sbamOpen filename "rb" nullPtr >>= newInHandle filename
finalizeSamFile :: MVar (Ptr SamFileInt) -> IO ()
finalizeSamFile mv = modifyMVar mv $ \fsam -> do
unless (fsam == nullPtr) $ sbamClose fsam
return (nullPtr, ())
closeInHandle :: InHandle -> IO ()
closeInHandle = finalizeSamFile . samfile
withTamInFile :: FilePath -> (InHandle -> IO a) -> IO a
withTamInFile filename = bracket (openTamInFile filename) closeInHandle
withTamInFileWithIndex :: FilePath -> FilePath -> (InHandle -> IO a) -> IO a
withTamInFileWithIndex filename indexname = bracket (openTamInFileWithIndex filename indexname) closeInHandle
withBamInFile :: FilePath -> (InHandle -> IO a) -> IO a
withBamInFile filename = bracket (openBamInFile filename) closeInHandle
get1 :: InHandle -> IO (Maybe Bam1)
get1 inh = withMVar (samfile inh) $ \fsam -> do
b <- bamInit1
res <- sbamRead fsam b
if res < 0
then do bamDestroy1 b
if res < 1
then ioError . userError $ "Error reading from BAM file " ++ show (inFilename inh)
else return Nothing
else do bptr <- newForeignPtr bamDestroy1Ptr b
return . Just $ Bam1 { ptrBam1 = bptr, header = inHeader inh }
data OutHandle = OutHandle { outFilename :: !FilePath
, outfile :: !(MVar (Ptr SamFileInt))
, outHeader :: !Header
}
newOutHandle :: String -> FilePath -> Header -> IO OutHandle
newOutHandle mode filename hdr = do
fsam <- withForeignPtr (unHeader hdr) $ sbamOpen filename mode . castPtr
when (fsam == nullPtr) $ ioError . userError $ "Error opening BAM file " ++ show filename
mv <- newMVar fsam
addMVarFinalizer mv (finalizeSamFile mv)
return $ OutHandle { outFilename = filename, outfile = mv, outHeader = hdr }
openTamOutFile :: FilePath -> Header -> IO OutHandle
openTamOutFile = newOutHandle "wh"
openBamOutFile :: FilePath -> Header -> IO OutHandle
openBamOutFile = newOutHandle "wb"
closeOutHandle :: OutHandle -> IO ()
closeOutHandle = finalizeSamFile . outfile
withTamOutFile :: FilePath -> Header -> (OutHandle -> IO a) -> IO a
withTamOutFile filename hdr = bracket (openTamOutFile filename hdr) closeOutHandle
withBamOutFile :: FilePath -> Header -> (OutHandle -> IO a) -> IO a
withBamOutFile filename hdr = bracket (openBamOutFile filename hdr) closeOutHandle
put1 :: OutHandle -> Bam1 -> IO ()
put1 outh b = withMVar (outfile outh) $ \fsam ->
withForeignPtr (ptrBam1 b) $ \p ->
sbamWrite fsam p >>= handleRes
where handleRes res | res > 0 = return ()
| otherwise = ioError . userError $ "Error writing to BAM file " ++ show (outFilename outh)