diff --git a/bamstats.cabal b/bamstats.cabal --- a/bamstats.cabal +++ b/bamstats.cabal @@ -1,5 +1,5 @@ Name: bamstats -Version: 0.3 +Version: 0.4 License: GPL Cabal-Version: >= 1.6 Build-Type: Simple @@ -7,6 +7,7 @@ Author: Ketil Malde Maintainer: Ketil Malde Synopsis: A program to extract various information from BAM alignmnet files. +Homepage: http://blog.malde.org/posts/bamstats.html Library Hs-Source-Dirs: src diff --git a/src/Bam.hs b/src/Bam.hs --- a/src/Bam.hs +++ b/src/Bam.hs @@ -6,6 +6,7 @@ import Bio.SamTools.Classify import Bio.SamTools.Bam import System.Console.CmdArgs +import Control.Monad (when) data Options = Stats { numrd :: Maybe Int, inputs :: [FilePath] } | Hist { numrd :: Maybe Int, bins, maxdist :: Int, inputs :: [FilePath], plot :: Maybe String } @@ -27,12 +28,12 @@ numrd = Nothing &= help "max number of reads to include" -- , bins = 200 &= help "number of bins" &= typ "INT" , mindist = 100 &= help "minimum insert size to count" &= typ "INT" - , delta = 1.01 &= help "relative size of next bin" &= typ "FLOAT" - , inputs = [] &= args &= typ "BAM file(s)" + , delta = 1.01 &= help "relative size of next bin" &= typ "FLOAT" + , inputs = [] &= args &= typ "BAM file(s)" } &= help "output approximate quantiles instead of the full histogram" -dump = Dump { numrd = Just 100 &= help "max number of reads (default: 100)" - , inputs = [] &= args &= typ "BAM file(s)" - } &= help "Dump alignments in the different classes." +dump = Dump { numrd = Just 100 &= help "max number of reads (default: 100)" + , inputs = [] &= args &= typ "BAM file(s)" + } &= help "Dump alignments in the different classes." histgen :: Options -> Hist histgen o = H 0 [(b,0) | b <- enumFromThenTo size (2*size) (maxdist o)] @@ -43,10 +44,11 @@ main :: IO () main = do - o <- cmdArgs $ modes [clfy,hst,quants,dump] - &= help "Extract information from BAM files" + o <- cmdArgs $ modes [clfy,hst,quants,dump] + &= help "Extract information from BAM files" &= program "bam" &= summary "bam v0.3, ©2012 Ketil Malde" - let geninp f = (case numrd o of Just x -> take x; Nothing -> id) `fmap` readBams f + when (null $ inputs o) $ error "No input files specified" + let geninp f = do {putStrLn ("## Input file: "++f) ; (case numrd o of Just x -> take x; Nothing -> id) `fmap` readBams f} genout = putStrLn . case o of Stats {} -> display . classify (CS 0 0 0 0 0) Hist {} -> (case plot o of Nothing -> display Just str -> genplot str) . classify (histgen o) diff --git a/src/Bio/SamTools/Classify.hs b/src/Bio/SamTools/Classify.hs --- a/src/Bio/SamTools/Classify.hs +++ b/src/Bio/SamTools/Classify.hs @@ -10,7 +10,7 @@ import Data.Maybe (isNothing, fromJust) import Text.Printf (printf) -data Stats a = Class { total :: !Int, innies, outies, lefties, righties :: !a } deriving Show +data Stats a = Class { total, unmapped, orphans, splits :: !Int, innies, outies, lefties, righties :: !a } deriving Show -- -------------------------------------------------- -- The generic collection framework @@ -19,22 +19,31 @@ -- | Pretty-print a 'Stats' value. display :: Insertable x => Stats x -> String display c = unlines $ map (intercalate "\t") - [ "Alignment" : dispheader (innies c) + [ "#Alignment" : dispheader (innies c) , "innies " : disp1 t (innies c) , "outies " : disp1 t (outies c) , "lefties " : disp1 t (lefties c) , "righties " : disp1 t (righties c) - , ["Total reads: ", show t]] + , summarize c t] where t = total c +summarize :: Stats a -> Int -> [String] +summarize c t = [printf "\nTotal reads: %7d\n" t + ++" unmapped: "++show1 unmapped++"\n" + ++" orphans: "++show1 orphans++"\n" + ++" split pairs: "++show2 splits] + where show1 f = printf "%7d" (f c)++percent (f c) + show2 f = printf "%7d" (f c`div`2)++percent (f c) + percent x = printf " (%.1f%%)" (100.0*fromIntegral x/fromIntegral t::Double) + showQuants :: Stats Hist -> String showQuants cs = unlines $ map (intercalate "\t") - [ "Alignment" : header (innies cs) + [ "#Alignment" : header (innies cs) , "innies " : quants t (innies cs) , "outies " : quants t (outies cs) , "lefties " : quants t (lefties cs) , "righties " : quants t (righties cs) - , ["Total reads: ", show t]] + , summarize cs t] where t = total cs percentiles = [0.05,0.25,0.5,0.75,0.95] :: [Double] quants tot h = printf "%14d" (hcount h) @@ -76,7 +85,7 @@ -- | Extract info from alignments classify :: Insertable x => x -> [Bam1] -> Stats x -classify def = foldl' (class1 . bump) (Class 0 def def def def) +classify def = foldl' (class1 . bump) (Class 0 0 0 0 def def def def) -- | Update count bump :: Stats a -> Stats a @@ -85,20 +94,27 @@ -- | Update data structure with a single alignment class1 :: Insertable x => Stats x -> Bam1 -> Stats x class1 c0 b - | isUnmapped b = c0 + | isNothing (insertSize b) = add_unmapped b c0 | isOpposite b = (if firstUpstream b then add_innie else add_outie) b c0 | otherwise = (if firstUpstream b then add_rightie else add_leftie) b c0 -add_innie, add_outie, add_rightie, add_leftie :: Insertable x => Bam1 -> Stats x -> Stats x +add_innie, add_outie, add_rightie, add_leftie, add_unmapped + :: Insertable x => Bam1 -> Stats x -> Stats x add_innie b c = c { innies = insert b (innies c) } add_outie b c = c { outies = insert b (outies c) } add_rightie b c = c { righties = insert b (righties c) } add_leftie b c = c { lefties = insert b (lefties c) } -isUnmapped, isOpposite, firstUpstream, isBefore :: Bam1 -> Bool -isUnmapped = isNothing . insertSize +add_unmapped b c + | isUnmap b = c { unmapped = unmapped c + 1 } + | isMateUnmap b = c { orphans = orphans c + 1 } + -- both mateTgtID and tgtId should be Just something + | mateTargetID b /= targetID b = c { splits = splits c + 1 } + | otherwise = c + +isOpposite, firstUpstream, isBefore :: Bam1 -> Bool -- isUnmap b || isMateUnmap b || mateTargetID b /= targetID b -- not the same -- insertSize is Nothing if negative (downstream read) isOpposite b = isReverse b /= isMateReverse b @@ -130,7 +146,7 @@ kt = (x4 - 4*m*x3 + 6*m2*x2 - 4*m3*x + n*m*m3)/(s*s*s*s*n) - 3 statistics :: Stats ClassStats -> Stats Statistics -statistics (Class t i o l r) = (Class t (mkstats i) (mkstats o) (mkstats l) (mkstats r)) -- todo: Functor instance? +statistics (Class t ump orp sp i o l r) = (Class t ump orp sp (mkstats i) (mkstats o) (mkstats l) (mkstats r)) -- todo: Functor instance? instance Insertable ClassStats where insert b (CS c s s2 s3 s4) = CS (c+1) (s+d) (s2+d^(2::Int)) (s3+d^(3::Int)) (s4+d^(4::Int))