diff --git a/bamstats.cabal b/bamstats.cabal --- a/bamstats.cabal +++ b/bamstats.cabal @@ -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 -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 diff --git a/src/Bam.hs b/src/Bam.hs --- a/src/Bam.hs +++ b/src/Bam.hs @@ -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 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 @@ -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 []