packages feed

bamstats 0.3 → 0.4

raw patch · 3 files changed

+39/−20 lines, 3 files

Files

bamstats.cabal view
@@ -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 <ketil@malde.org> Synopsis:      A program to extract various information from BAM alignmnet files.+Homepage:      http://blog.malde.org/posts/bamstats.html  Library     Hs-Source-Dirs:  src
src/Bam.hs view
@@ -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)
src/Bio/SamTools/Classify.hs view
@@ -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))