bamstats 0.1 → 0.2

raw patch · 3 files changed

+49/−14 lines, 3 filesPVP ok

version bump matches the API change (PVP)

API changes (from Hackage documentation)

Files

• bamstats.cabal +2/−2
• src/Bam.hs +21/−8
• src/Bio/SamTools/Classify.hs +26/−4

@@ -1,12 +1,12 @@ Name:          bamstats-Version:       0.1+Version:       0.2 License:       GPL Cabal-Version: >= 1.6 Build-Type:    Simple Category:      Bioinformatics Author:        Ketil Malde Maintainer:    Ketil Malde <ketil@malde.org>-Synopsis:      A program to extract various information from BAM alignmnet files+Synopsis:      A program to extract various information from BAM alignmnet files.  Library     Hs-Source-Dirs:  src

@@ -9,10 +9,11 @@  data Options = Stats { numrd :: Maybe Int, inputs :: [FilePath] }              | Hist { numrd :: Maybe Int, bins, maxdist :: Int, inputs :: [FilePath] }+             | Quants { numrd :: Maybe Int, mindist :: Int, delta :: Double, inputs :: [FilePath] }              | Dump { numrd :: Maybe Int, inputs :: [FilePath] }              deriving (Data,Typeable,Read,Show) -clfy, hst, dump :: Options+clfy, hst, quants, dump :: Options clfy = Stats { numrd = Nothing &= help "max number of reads (default: all)"                 , inputs = []     &= args &= typ "BAM file(s)"                 } &= help "Calculate statistics on insert sizes"@@ -21,20 +22,32 @@             , maxdist = 1000  &= help "max insert size to count" &= typ "INT"             , inputs = []     &= args &= typ "BAM file(s)"             } &= help "Collect a histogram of insert sizes"+quants = Quants {+              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)"                +            } &= 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."-       +            } &= help "Dump alignments in the different classes."       ++histgen :: Options -> Hist histgen o = H 0 [(b,0) | b <- enumFromThenTo size (2*size) (maxdist o)]   where size = (maxdist o `div` bins o) +quantgen :: Options -> Hist+quantgen o = H 0 [(round b,0) | b <- scanl1 (+) $ iterate (* delta o) (fromIntegral $ mindist o)]+ main :: IO () main = do-  o <- cmdArgs $ modes [clfy,hst,dump] +  o <- cmdArgs $ modes [clfy,hst,quants,dump]         &= help "Extract information from BAM files" -       &= program "bam" &= summary "bam v0.0, ©2011 Ketil Malde"+       &= program "bam" &= summary "bam v0.0, ©2012 Ketil Malde"   let geninp f = (case numrd o of Just x -> take x; Nothing -> id) `fmap` readBams f-      genout = case o of Stats {} -> putStrLn . display . classify (cdef :: ClassStats)-                         Hist {}  -> putStrLn . display . classify (histgen o)-                         Dump {}  -> putStrLn . display . classify (cdef :: Collect)+      genout = putStrLn . case o of Stats {} -> display . classify (CS 0 0 0 0 0)+                                    Hist {}  -> display . classify (histgen o)+                                    Quants {} -> showQuants . classify (quantgen o)+                                    Dump {}  -> display . classify (Bams [])   mapM_ (\f -> genout =<< geninp f) $ inputs o

@@ -26,11 +26,36 @@   , ["Total reads: ", show t]]   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)+  , ["Total reads: ", show t]]+  where t = total c+        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)+              : if hcount h == 0 then map (const "    N/A") percentiles+                else go (map (round . (*fromIntegral (hcount h))) percentiles) 0 (0,0) (buckets h)+        header _h = "         count":"  prop":map (printf "%6.0f%%" . (*100)) percentiles+        go :: [Int] -> Int -> (Int,Int) -> [(Int,Int)] -> [String]+        go (0:fs) sum prev inp = "      0" : go fs sum prev inp+        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))+               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!?"+ class Insertable x where   insert :: Bam1 -> x -> x   disp1  :: Int  -> x -> [String]   dispheader :: x -> [String]-  cdef :: x    -- | Extract info from alignments classify :: Insertable x => x -> [Bam1] -> Stats x@@ -94,7 +119,6 @@   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))      where d = maybe (error ("no insert size?\n"++show b)) fromIntegral $ insertSize b   dispheader _ = ["         count","   prop","   mean","  stdev","   skew","   kurt"]-  cdef = CS 0 0 0 0 0   disp1 tot cs = printf "%14d" (ccount cs)                 : printf "%5.2f%%" (200*fromIntegral (ccount cs)/fromIntegral tot::Double)                : map (printf "%7.1f") [mean s, stdev s, skew s, kurt s]@@ -117,7 +141,6 @@               : printf "%5.2f%%" (200*fromIntegral (hcount h)/fromIntegral tot::Double)               : map (printf "%7d" . snd) (buckets h)   dispheader h = "         count":"  prop":(map (printf "%7d" . fst) $ init $ buckets h)++["    >"]-  cdef = undefined  -- -------------------------------------------------- -- Just "Collect" all the Bams in different classes@@ -129,4 +152,3 @@   insert b (Bams bs) = Bams (b:bs)   disp1 _ (Bams bs) = [unlines $ ("":map show bs)]   dispheader = const ["foo"]-  cdef = Bams []