src/Bio/SamTools/BamIndex.hs

-- | This module provides an interface to sorted, indexed BAM
-- alignment files, which allow rapid extraction of alignments that lie
-- within one specific region of one sequence.
module Bio.SamTools.BamIndex
 ( 
 IdxHandle, idxFilename, idxHeader
 , open, close, withIndex
 , Query, qyHandle
 , query, next
 ) 
 where 

import Control.Concurrent.MVar
import Control.Exception
import Control.Monad
import Data.Int (Int64)
import Foreign.ForeignPtr
import Foreign.Ptr

import Bio.SamTools.Bam
import Bio.SamTools.Internal
import Bio.SamTools.LowLevel 
 
-- | Handle for fetching alignments by region from a sorted, indexed
-- BAM file.
data IdxHandle = IdxHandle { idxFilename :: !FilePath -- ^ Filename of sorted, indexed BAM file
 , bamindex :: !(MVar (Ptr BamFileInt, Ptr BamIndexInt))
 , idxHeader :: !Header -- ^ Target sequences
 }
 
-- | Open a sorted, indexed BAM file. 
open :: FilePath -> IO IdxHandle
open filename = do
 f <- bamOpen filename "r"
 when (f == nullPtr) $ ioError . userError 
 $ "Error opening BAM file " ++ show filename
 i <- bamIndexLoad filename
 when (i == nullPtr) $ ioError . userError 
 $ "Error opening index for BAM file " ++ show filename
 mv <- newMVar (f, i)
 bhdr <- bamHeaderRead f
 when (bhdr == nullPtr) $ ioError . userError $
 "Error reading header from BAM file " ++ show filename
 bamInitHeaderHash bhdr
 addMVarFinalizer mv (finalizeBamIndex mv)
 hdr <- liftM Header . newForeignPtr bamHeaderDestroyPtr $ bhdr
 return $ IdxHandle { idxFilename = filename
 , bamindex = mv
 , idxHeader = hdr
 }
 
close :: IdxHandle -> IO ()
close = finalizeBamIndex . bamindex

finalizeBamIndex :: MVar (Ptr BamFileInt, Ptr BamIndexInt) -> IO ()
finalizeBamIndex mv = modifyMVar mv $ \(f, i) -> do
 unless (f == nullPtr) $ bamClose f >> return ()
 unless (i == nullPtr) $ bamIndexDestroy i
 return ((nullPtr, nullPtr), ())

data Query = Query { qyHandle :: !IdxHandle
 , iter :: !(MVar (Ptr BamIterInt))
 }
 
withIndex :: FilePath -> (IdxHandle -> IO a) -> IO a
withIndex filename = bracket (open filename) close
 
query :: IdxHandle -> Int -> (Int64, Int64) -> IO Query
query inh tid (start, end) = withMVar (bamindex inh) $ \(_f, idx) -> do
 it <- bamIterQuery idx (fromIntegral tid) (fromIntegral start) (fromIntegral end)
 when (it == nullPtr) $ ioError . userError
 $ "Error starting BAM query: " ++ show (idxFilename inh, (tid, (start, end)))
 mv <- newMVar it
 addMVarFinalizer mv (finalizeBamIter mv)
 return $ Query { qyHandle = inh, iter = mv }
 
finalizeBamIter :: MVar (Ptr BamIterInt) -> IO ()
finalizeBamIter mv = modifyMVar mv $ \it -> do
 unless (it == nullPtr) $ bamIterDestroy it
 return (nullPtr, ())
 
next :: Query -> IO (Maybe Bam1)
next rgn = withMVar (bamindex . qyHandle $ rgn) $ \(f, _idx) -> 
 withMVar (iter rgn) $ \it -> do
 b <- bamInit1
 res <- bamIterRead f it b
 if res > 0
 then do bptr <- newForeignPtr bamDestroy1Ptr b
 return . Just $ Bam1 { ptrBam1 = bptr, header = idxHeader . qyHandle $ rgn }
 else do bamDestroy1 b
 if res < -1
 then ioError . userError $
 "Error reading BAM query from " ++ show (idxFilename . qyHandle $ rgn)
 else return Nothing