bamstats 0.2 → 0.3

raw patch · 3 files changed

+33/−13 lines, 3 filesPVP ok

version bump matches the API change (PVP)

API changes (from Hackage documentation)

Files

• bamstats.cabal +1/−1
• src/Bam.hs +5/−3
• src/Bio/SamTools/Classify.hs +27/−9

@@ -1,5 +1,5 @@ Name:          bamstats-Version:       0.2+Version:       0.3 License:       GPL Cabal-Version: >= 1.6 Build-Type:    Simple

@@ -8,7 +8,7 @@ import System.Console.CmdArgs  data Options = Stats { numrd :: Maybe Int, inputs :: [FilePath] }-             | Hist { numrd :: Maybe Int, bins, maxdist :: Int, inputs :: [FilePath] }+             | Hist { numrd :: Maybe Int, bins, maxdist :: Int, inputs :: [FilePath], plot :: Maybe String }              | Quants { numrd :: Maybe Int, mindist :: Int, delta :: Double, inputs :: [FilePath] }              | Dump { numrd :: Maybe Int, inputs :: [FilePath] }              deriving (Data,Typeable,Read,Show)@@ -21,6 +21,7 @@             , bins  = 20      &= help "number of bins" &= typ "INT"             , maxdist = 1000  &= help "max insert size to count" &= typ "INT"             , inputs = []     &= args &= typ "BAM file(s)"+            , plot = Nothing  &= help "gnuplot output, with optional gnuplot commands" &= opt ""             } &= help "Collect a histogram of insert sizes" quants = Quants {               numrd = Nothing &= help "max number of reads to include"@@ -44,10 +45,11 @@ main = do   o <- cmdArgs $ modes [clfy,hst,quants,dump]         &= help "Extract information from BAM files" -       &= program "bam" &= summary "bam v0.0, ©2012 Ketil Malde"+       &= 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       genout = putStrLn . case o of Stats {} -> display . classify (CS 0 0 0 0 0)-                                    Hist {}  -> display . classify (histgen o)+                                    Hist {}  -> (case plot o of Nothing -> display+                                                                Just str -> genplot str) . classify (histgen o)                                     Quants {} -> showQuants . classify (quantgen o)                                     Dump {}  -> display . classify (Bams [])   mapM_ (\f -> genout =<< geninp f) $ inputs o

@@ -4,6 +4,7 @@  module Bio.SamTools.Classify where +import Prelude hiding (sum) import Bio.SamTools.Bam import Data.List (foldl', intercalate) import Data.Maybe (isNothing, fromJust)@@ -27,14 +28,14 @@   where t = total c  showQuants :: Stats Hist -> String-showQuants c = unlines $ map (intercalate "\t")-  [  "Alignment" : header (innies c)-  ,  "innies   " : quants t (innies c)-  ,  "outies   " : quants t (outies c)-  ,  "lefties  " : quants t (lefties c)-  ,  "righties " : quants t (righties c)+showQuants cs = unlines $ map (intercalate "\t")+  [  "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]]-  where t = total c+  where t = total cs         percentiles = [0.05,0.25,0.5,0.75,0.95] :: [Double]         quants tot h = printf "%14d" (hcount h)               : printf "%5.2f%%" (200*fromIntegral (hcount h)/fromIntegral tot::Double)@@ -46,12 +47,28 @@         go (f:fs) sum prev ((b,c):rest) =            let next = sum+c in            if next >= f-          then let b' = b - round (fromIntegral (next-f)/fromIntegral c*fromIntegral (b-fst prev))+          then let b' = b - round (fromIntegral (next-f)/fromIntegral c*fromIntegral (b-fst prev)::Double)                in printf "%7d" b' : (go fs sum prev ((b,c):rest))           else     go (f:fs) next (b,c) rest         go [] _ _ _ = []         go _ _ _ [] = error "ran out of reads!?" +genplot :: String -> Stats Hist -> String+genplot cmds h = preamble ++ cmds ++ "\n" ++ plot+  where preamble = unlines ["set style data boxes"+                           ,"set style fill solid border 0"+                           ,"set xlabel 'insert length'"+                           ,"set ylabel 'read count'"]+        plot = "plot '-' ti 'innies', '-' ti 'outies', '-' ti 'lefties', '-' ti 'righties'\n"+                  ++concat [ plot1 $ buckets $ innies h+                           , plot1 $ map (\(x,y)->(-x,-y)) $ buckets $ outies h+                           , plot1 $ map (\(x,y)->(-x, y)) $ buckets $ lefties h+                           , plot1 $ map (\(x,y)->( x,-y)) $ buckets $ righties h]+        plot1 = unlines . go 0+        go prev ((b,c):rest@(_:_)) = (show ((prev+b) `div` 2) ++ " " ++ show c) : go (b+1) rest+        go _ [_] = ["e"] -- ignore last bucket, which is a catch-all+        go _ [] = error "no data"+         class Insertable x where   insert :: Bam1 -> x -> x   disp1  :: Int  -> x -> [String]@@ -135,7 +152,8 @@   insert b h = H (hcount h + 1) (go (fromIntegral $ fromJust $ insertSize b) (buckets h))     where go x ((b1,v):bs@(_:_)) =             if x > b1 then (b1,v):go x bs else inc (b1,v):bs-          go x [(b1,v)] = [inc (b1,v)]+          go _ [(b1,v)] = [inc (b1,v)]+          go _ [] = error "this never happens"           inc (x,y) = let y' = y+1 in y' `seq` (x,y')   disp1 tot h = printf "%14d" (hcount h)               : printf "%5.2f%%" (200*fromIntegral (hcount h)/fromIntegral tot::Double)