samtools-iteratee (empty) → 0.1
raw patch · 6 files changed
+254/−0 lines, 6 filesdep +basedep +bytestringdep +haskell98setup-changed
Dependencies added: base, bytestring, haskell98, iteratee, monads-tf, samtools, transformers, vector
Files
- LICENSE +21/−0
- Setup.hs +2/−0
- samtools-iteratee.cabal +42/−0
- src/BamCount.hs +75/−0
- src/BamFilter.hs +58/−0
- src/Bio/SamTools/Iteratee.hs +56/−0
+ LICENSE view
@@ -0,0 +1,21 @@+The MIT License++Copyright (c) 2011 Nicholas Ingolia++Permission is hereby granted, free of charge, to any person obtaining a copy+of this software and associated documentation files (the "Software"), to deal+in the Software without restriction, including without limitation the rights+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell+copies of the Software, and to permit persons to whom the Software is+furnished to do so, subject to the following conditions:++The above copyright notice and this permission notice shall be included in+all copies or substantial portions of the Software.++THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN+THE SOFTWARE.
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ samtools-iteratee.cabal view
@@ -0,0 +1,42 @@+Name: samtools-iteratee+Version: 0.1+Synopsis: Iteratee interface to SamTools library+Description: Iteratee interface to SamTools library+License: MIT+License-file: LICENSE+Author: Nicholas Ingolia+Maintainer: nick@ingolia.org+Category: Bio++Build-type: Simple+Cabal-version: >=1.4++Flag Utilities+ Description: Build utility programs+ Default: False++Library+ Exposed-modules: Bio.SamTools.Iteratee+ Build-depends: base >= 4.3 && < 5, bytestring >= 0.9 && < 0.10, haskell98, samtools, + transformers, iteratee >= 0.8 && < 0.9+ Hs-Source-Dirs: src++Executable bam-filter+ if !flag(Utilities)+ Buildable: False+ Main-Is: BamFilter.hs+ Other-Modules: Bio.SamTools.Iteratee+ Build-depends: base >= 4.3 && < 5, bytestring >= 0.9 && < 0.10, haskell98, samtools, + transformers, iteratee >= 0.8 && < 0.9, monads-tf+ Hs-Source-Dirs: src+ Ghc-options: -Wall++Executable bam-count+ if !flag(Utilities)+ Buildable: False+ Main-Is: BamCount.hs+ Other-Modules: Bio.SamTools.Iteratee+ Build-depends: base >= 4.3 && < 5, bytestring >= 0.9 && < 0.10, haskell98, samtools, + transformers, iteratee >= 0.8 && < 0.9, monads-tf, vector+ Hs-Source-Dirs: src+ Ghc-options: -Wall
+ src/BamCount.hs view
@@ -0,0 +1,75 @@+module Main+ where++import Control.Applicative+import Control.Monad.Reader+import qualified Data.ByteString.Char8 as BS+import Data.Int (Int64)+import Data.List+import Data.Maybe++import System.Console.GetOpt+import System.Environment+import System.IO++import qualified Data.Iteratee as Iter+import qualified Data.Vector.Unboxed as U+import qualified Data.Vector.Unboxed.Mutable as UM++import qualified Bio.SamTools.Bam as Bam+import qualified Bio.SamTools.Iteratee as Bam++main :: IO ()+main = getArgs >>= handleOpt . getOpt RequireOrder optDescrs+ where handleOpt (_, _, errs@(_:_)) = usage (unlines errs)+ handleOpt (args, [bam], []) = either usage (doBamCount bam) $ argsToConf args+ handleOpt (_, _, []) = usage "Specify exactly one BAM input"+ usage errs = do prog <- getProgName+ hPutStr stderr $ usageInfo prog optDescrs+ hPutStrLn stderr errs++doBamCount :: FilePath -> Conf -> IO ()+doBamCount bam conf = Bam.withBamInFile bam $ \hin -> do+ tctio <- UM.replicate (Bam.nTargets . Bam.inHeader $ hin) 0+ let countBam b = flip (maybe (return ())) (Bam.targetID b) $ \tid -> do + n0 <- UM.read tctio tid+ UM.write tctio tid $! succ n0+ bamiter = Iter.joinI $ Iter.filter (wanted conf) $ Iter.mapM_ countBam+ Bam.enumInHandle hin bamiter >>= Iter.run+ U.freeze tctio >>= writeFile (confOutput conf) . targetTable (Bam.inHeader hin)+ +targetTable :: Bam.Header -> U.Vector Int64 -> String+targetTable h = unlines . U.ifoldr addTargetLine []+ where addTargetLine idx n = (targetLine :)+ where targetLine = intercalate "\t" [ BS.unpack $ Bam.targetSeqName h idx+ , show $ Bam.targetSeqLen h idx+ , show n+ ]++wanted :: Conf -> Bam.Bam1 -> Bool+wanted conf | confPerfect conf = maybe False (== 0) . Bam.nMismatch+ | otherwise = const True+ +data Conf = Conf { confOutput :: !FilePath+ , confPerfect :: !Bool+ } deriving (Show)++data Arg = ArgOutput { unArgOutput :: !String }+ | ArgPerfect+ deriving (Show, Read, Eq, Ord)++argOutput :: Arg -> Maybe String+argOutput (ArgOutput del) = Just del+argOutput _ = Nothing++optDescrs :: [OptDescr Arg]+optDescrs = [ Option ['o'] ["output"] (ReqArg ArgOutput "OUTFILE") "Output filename"+ , Option [] ["perfect"] (NoArg ArgPerfect) "Filter perfect (NM:i:0) reads"+ ]++argsToConf :: [Arg] -> Either String Conf+argsToConf = runReaderT conf+ where conf = Conf <$> + findOutput <*>+ (ReaderT $ return . elem ArgPerfect)+ findOutput = ReaderT $ maybe (Left "No output file") return . listToMaybe . mapMaybe argOutput
+ src/BamFilter.hs view
@@ -0,0 +1,58 @@+module Main+ where++import Control.Applicative+import Control.Monad.Reader+import Data.Maybe++import System.Console.GetOpt+import System.Environment+import System.IO++import qualified Data.Iteratee as Iter++import qualified Bio.SamTools.Bam as Bam+import qualified Bio.SamTools.Iteratee as Bam++main :: IO ()+main = getArgs >>= handleOpt . getOpt RequireOrder optDescrs+ where handleOpt (_, _, errs@(_:_)) = usage (unlines errs)+ handleOpt (args, [bam], []) = either usage (doBamFilter bam) $ argsToConf args+ handleOpt (_, _, []) = usage "Specify exactly one BAM input"+ usage errs = do prog <- getProgName+ hPutStr stderr $ usageInfo prog optDescrs+ hPutStrLn stderr errs++doBamFilter :: FilePath -> Conf -> IO ()+doBamFilter bam conf = Bam.withBamInFile bam $ \hin ->+ Bam.withBamOutFile (confOutput conf) (Bam.inHeader hin) $ \hout ->+ let bamiter = Iter.joinI $ Iter.filter (wanted conf) $ Iter.mapM_ (Bam.put1 hout)+ in Bam.enumInHandle hin bamiter >>= Iter.run+ +wanted :: Conf -> Bam.Bam1 -> Bool+wanted conf | confPerfect conf = maybe False (== 0) . Bam.nMismatch+ | otherwise = const True+ +data Conf = Conf { confOutput :: !FilePath+ , confPerfect :: !Bool+ } deriving (Show)++data Arg = ArgOutput { unArgOutput :: !String }+ | ArgPerfect+ deriving (Show, Read, Eq, Ord)++argOutput :: Arg -> Maybe String+argOutput (ArgOutput del) = Just del+argOutput _ = Nothing++optDescrs :: [OptDescr Arg]+optDescrs = [ Option ['o'] ["output"] (ReqArg ArgOutput "OUTFILE") "Output filename"+ , Option [] ["perfect"] (NoArg ArgPerfect) "Filter perfect (NM:i:0) reads"+ ]++argsToConf :: [Arg] -> Either String Conf+argsToConf = runReaderT conf+ where conf = Conf <$> + findOutput <*>+ (ReaderT $ return . elem ArgPerfect)+ findOutput = ReaderT $ maybe (Left "No output file") return . listToMaybe . mapMaybe argOutput
+ src/Bio/SamTools/Iteratee.hs view
@@ -0,0 +1,56 @@+module Bio.SamTools.Iteratee+ (+ enumInHandle+ , enumTam, enumTamWithIndex, enumBam+ , enumQuery, enumIndexRegion, enumBamRegion+ ) +where+ +import Control.Exception+import Control.Monad.IO.Class+import qualified Data.ByteString.Char8 as BS++import qualified Bio.SamTools.Bam as Bam+import qualified Bio.SamTools.BamIndex as BamIndex+import qualified Data.Iteratee as Iter+ +enumInHandle :: Bam.InHandle -> Iter.Enumerator [Bam.Bam1] IO a+enumInHandle inh i0 = step i0+ where step iter = Bam.get1 inh >>= maybe eof next+ where eof = return iter+ next b = Iter.runIter iter Iter.idoneM onCont+ where onCont k Nothing = step . k $ Iter.Chunk [b]+ onCont k e = return $ Iter.icont k e+ +enumTam :: FilePath -> Iter.Enumerator [Bam.Bam1] IO a+enumTam inname i0 = bracket (Bam.openTamInFile inname) Bam.closeInHandle $ \h ->+ enumInHandle h i0+ +enumTamWithIndex :: FilePath -> FilePath -> Iter.Enumerator [Bam.Bam1] IO a+enumTamWithIndex inname idxname i0 + = bracket (Bam.openTamInFileWithIndex inname idxname) Bam.closeInHandle $ \h ->+ enumInHandle h i0+ +enumBam :: FilePath -> Iter.Enumerator [Bam.Bam1] IO a+enumBam inname i0 = bracket (Bam.openBamInFile inname) Bam.closeInHandle $ \h ->+ enumInHandle h i0+ +enumQuery :: BamIndex.Query -> Iter.Enumerator [Bam.Bam1] IO a+enumQuery q i0 = step i0+ where step iter = BamIndex.next q >>= maybe eof next+ where eof = return iter+ next b = Iter.runIter iter Iter.idoneM onCont+ where onCont k Nothing = step . k $ Iter.Chunk [b]+ onCont k e = return $ Iter.icont k e++enumIndexRegion :: BamIndex.IdxHandle -> Int -> (Int, Int) -> Iter.Enumerator [Bam.Bam1] IO a+enumIndexRegion h tid bnds i0 = do+ q <- liftIO $ BamIndex.query h tid bnds+ enumQuery q i0+ +enumBamRegion :: FilePath -> BS.ByteString -> (Int, Int) -> Iter.Enumerator [Bam.Bam1] IO a+enumBamRegion inname seqname bnds i0 = bracket (BamIndex.open inname) BamIndex.close $ \h -> do+ tid <- lookupTid . BamIndex.idxHeader $ h+ enumIndexRegion h tid bnds i0+ where lookupTid h = maybe noTid return $! Bam.lookupTarget h seqname+ noTid = ioError . userError $ "Target " ++ show seqname ++ " not found in " ++ show inname