diff --git a/samtools-iteratee.cabal b/samtools-iteratee.cabal --- a/samtools-iteratee.cabal +++ b/samtools-iteratee.cabal @@ -1,12 +1,12 @@ Name: samtools-iteratee -Version: 0.1.1 +Version: 0.2.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 +Category: Bioinformatics Build-type: Simple Cabal-version: >=1.4 @@ -17,7 +17,7 @@ Library Exposed-modules: Bio.SamTools.Iteratee - Build-depends: base >= 4.3 && < 5, bytestring >= 0.9 && < 0.10, haskell98, samtools >= 0.1.1, + Build-depends: base >= 4.2 && < 5, bytestring >= 0.9 && < 0.10, haskell98, samtools >= 0.1.1, transformers, iteratee >= 0.8 && < 0.9 Hs-Source-Dirs: src @@ -26,17 +26,7 @@ Buildable: False Main-Is: BamFilter.hs Other-Modules: Bio.SamTools.Iteratee - Build-depends: base >= 4.3 && < 5, bytestring >= 0.9 && < 0.10, haskell98, samtools >= 0.1.1, + Build-depends: base >= 4.2 && < 5, bytestring >= 0.9 && < 0.10, haskell98, samtools >= 0.1.1, 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 >= 0.1.1, - transformers, iteratee >= 0.8 && < 0.9, monads-tf, vector Hs-Source-Dirs: src Ghc-options: -Wall diff --git a/src/BamCount.hs b/src/BamCount.hs deleted file mode 100644 --- a/src/BamCount.hs +++ /dev/null @@ -1,75 +0,0 @@ -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 diff --git a/src/BamFilter.hs b/src/BamFilter.hs --- a/src/BamFilter.hs +++ b/src/BamFilter.hs @@ -3,6 +3,7 @@ import Control.Applicative import Control.Monad.Reader +import qualified Data.ByteString.Char8 as BS import Data.Maybe import System.Console.GetOpt @@ -12,27 +13,31 @@ import qualified Data.Iteratee as Iter import qualified Bio.SamTools.Bam as Bam -import qualified Bio.SamTools.Iteratee as Bam +import qualified Bio.SamTools.Iteratee as BamIter 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" + handleOpt (_, _, []) = usage "Specify exactly one BAM file" 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 - + Bam.withBamOutFile (confOutput conf) (Bam.inHeader hin) $ \hout -> + let handleBam bam | confIsWanted conf bam = Bam.put1 hout bam + | otherwise = return () + in BamIter.enumInHandle hin (Iter.mapM_ handleBam) >>= Iter.run + +confIsWanted :: Conf -> Bam.Bam1 -> Bool +confIsWanted conf | confPerfect conf = isPerfect + | otherwise = const True + +isPerfect :: Bam.Bam1 -> Bool +isPerfect = maybe False (== 0) . Bam.nMismatch + data Conf = Conf { confOutput :: !FilePath , confPerfect :: !Bool } deriving (Show) @@ -46,8 +51,8 @@ argOutput _ = Nothing optDescrs :: [OptDescr Arg] -optDescrs = [ Option ['o'] ["output"] (ReqArg ArgOutput "OUTFILE") "Output filename" - , Option [] ["perfect"] (NoArg ArgPerfect) "Filter perfect (NM:i:0) reads" +optDescrs = [ Option ['o'] ["output"] (ReqArg ArgOutput "OUTFILE") "Output filename" + , Option [] ["perfect"] (NoArg ArgPerfect) "Pass only perfect alignments" ] argsToConf :: [Arg] -> Either String Conf @@ -55,4 +60,4 @@ where conf = Conf <$> findOutput <*> (ReaderT $ return . elem ArgPerfect) - findOutput = ReaderT $ maybe (Left "No output file") return . listToMaybe . mapMaybe argOutput + findOutput = ReaderT $ maybe (Left "No output filename") return . listToMaybe . mapMaybe argOutput