a50 0.2 → 0.4
raw patch · 6 files changed
+84/−52 lines, 6 files
Files
- a50.cabal +10/−9
- src/A50.hs +29/−26
- src/Blat.hs +3/−2
- src/Gnuplot.hs +18/−13
- src/Options.hs +3/−2
- src/Sort.hs +21/−0
a50.cabal view
@@ -1,5 +1,5 @@ Name: a50-Version: 0.2+Version: 0.4 License: GPL Author: Ketil Malde@@ -21,15 +21,16 @@ . The Darcs repository is at <http://malde.org/~ketil/biohaskell/a50>. --- HomePage: http://blog.malde.org/area50-Build-Depends: base>=3 && <5, bio > 0.4, process, containers, cmdargs >= 0.5, directory+HomePage: http://blog.malde.org/index.php/a50-a-graphical-comparison-of-genome-assemblies Build-Type: Simple Cabal-Version: >= 1.2.3 Tested-with: GHC==6.12.1- Data-files: README-Executable: a50-Main-Is: A50.hs-Other-Modules: Gnuplot, Options, Blat-Hs-Source-Dirs: src-Ghc-Options: -Wall++Executable a50+ Main-Is: A50.hs+ Other-Modules: Gnuplot, Options, Blat, Sort+ Hs-Source-Dirs: src+ Ghc-Options: -Wall+ Build-Depends: base>=3 && <5, bio > 0.4, process, containers, cmdargs >= 0.5, directory+
src/A50.hs view
@@ -1,13 +1,12 @@ module Main where -import Bio.Sequence-import Data.List (intersperse, foldl')--import qualified Data.IntMap as M+import Data.List (intersperse) import Gnuplot import Blat-import Options+import Sort (sizes, sort, SizeTable)+import Options (Opt(..), getArgs)+import Control.Monad (when) main :: IO () main = do @@ -15,16 +14,32 @@ ss <- mapM sizes $ inputs opts case estref opts of "" -> do+ when (format opts /= "plot") $ n50 (map read $ expect opts) $ zip (inputs opts) ss mkplot opts $ map (each 10) $ map ((scanl1 (+)) . sort) ss est -> do ps <- mapM (\asm -> fmap gen_result $ runBlat (tmpdir opts) asm est) (inputs opts) mkplot opts $ map (each 10) $ map (scanl1 (+)) $ zipWith interleave (map sort ss) ps -- putStr $ zipLists fs (map ((scanl1 (+)) . sort) ss) -each :: Int -> [a] -> [a]-each _ [] = []-each n (x:xs) = x : each n (drop (n-1) xs)+n50 :: [Int] -> [(FilePath,SizeTable)] -> IO ()+n50 [e] iss = do+ let genOut (f,st) = putStrLn (f ++ "\t"++go [e`div`4,e`div`2,3*e`div`4] (scanl1 (+) $ sort st))+ go (e0:es) (z1:z2:zs) = + if z2 > e0 then show (z2-z1)++"\t"++go es (z1:z2:zs)+ else go (e0:es) (z2:zs)+ go [] (z1:z2:zs) = z1 `seq` go [] (z2:zs)+ go es [z1] = concatMap (const "-\t") es ++ show z1++" total size"+ go _ [] = error "shouldn't happen"+ putStrLn ("Assembly \tn25\tn50\tn75\t\testimated size="++show e)+ mapM_ genOut iss+n50 _ _ = return () +each :: Int -> [a] -> [(Int,a)]+each k = go 0+ where+ go _ [] = []+ go i (x:xs) = (i,x) : go (i+k) (drop (k-1) xs)+ zipLists :: [String] -> [[Int]] -> String zipLists fs ss = unlines ((concat $ map ("#\t"++) fs) : go (map (++repeat 0) ([1..]:ss))) where go :: [[Int]] -> [String]@@ -33,28 +48,16 @@ in if all (==0) $ tail hs then [] else (concat $ intersperse "\t" $ map myshow hs) : go ts -myshow :: Integral i => i -> String+myshow :: (Integral i, Show i) => i -> String myshow 0 = "" myshow n = show n -mkplot :: Opt -> [[Int]] -> IO ()-mkplot o ns = gnuplot [conf,outp,labels,tics] (zip (inputs o) ns) es+mkplot :: Opt -> [[(Int,Int)]] -> IO ()+mkplot o ns = let my_out = if format o == "plot" then gnudat else gnuplot+ in my_out [conf,outp,labels,tics] (zip (inputs o) ns) es where labels = "set ylabel 'cumulative size'; set xlabel 'contig number'"- tics = "set format y '%.0s%c'; set format x '%.0f0'"- conf = if null $ format o then "" else "set terminal "++format o+ tics = "set format y '%.1s%c'; set format x '%.0f'"+ conf = if (null $ format o) || format o == "plot" then "" else "set terminal "++format o outp = if null $ outfile o then "" else "set out '"++outfile o++"'" es = map read $ expect o -sizes :: FilePath -> IO [Int]-sizes f = map (fromIntegral . seqlength) `fmap` readFasta f--sort :: [Int] -> [Int]-sort = concatMap (\(x,c) -> replicate c (fromIntegral $ negate x)) . M.toAscList . freqs---- equivalent to 'M.fromList . map (\x->(fromIntegral $ negate x,1))', --- except for not blowing the stack-freqs :: [Int] -> M.IntMap Int-freqs = foldl' ins M.empty . map fromIntegral . map negate- where ins m x = case M.lookup x m of - Just v -> v `seq` M.insert x (v+1) m- Nothing -> M.insert x 1 m
src/Blat.hs view
@@ -19,11 +19,12 @@ runBlat tmpdir asm ests = do let pslfile = tmpdir++"/"++basename asm++"_vs_"++basename ests++".psl" basename = reverse . takeWhile (/='/') . reverse+ quote str = "\""++str++"\"" dfe <- doesFileExist pslfile when (not dfe) $ do fe <- findExecutable "blat" when (fe == Nothing) $ error "Couldn't find the 'blat' executable - aborting"- blat <- runCommand $ unwords [blat_cmd,asm,ests,pslfile]+ blat <- runCommand $ unwords (blat_cmd : map quote [asm,ests,pslfile]) e <- waitForProcess blat case e of ExitSuccess -> return () _ -> error $ show e@@ -33,7 +34,7 @@ gen_result :: [PSL] -> [(Int,Int)] gen_result = coverage . order . pslbest -interleave :: Integral i => [i] -> [(i,i)] -> [i]+interleave :: (Integral i, Show i) => [i] -> [(i,i)] -> [i] interleave sz [] = map (const 0) sz interleave (sz:szs) covs@((s1,c1):cs) | sz > s1 = 0 : interleave szs covs
src/Gnuplot.hs view
@@ -9,7 +9,7 @@ import System.Directory (findExecutable) import Control.Concurrent (forkIO) -gnuplot :: Num i => [String] -> [(String,[i])] -> [Int] -> IO ()+gnuplot :: (Show i, Num i) => [String] -> [(String,[(Int,i)])] -> [Double] -> IO () gnuplot preamble cols hlines = do fe <- findExecutable "gnuplot" when (fe == Nothing) $ error "Couldn't find the 'gnuplot' executable - aborting"@@ -17,24 +17,29 @@ _ <- forkIO (hGetContents o >>= hPutStr stdout) -- ugly hack to limit error output _ <- forkIO $ do - let go ~(l1:l2:ls) =do putStrLn l1- if ('^' `elem` l1) - then do - putStrLn l2- else go (l2:ls)+ let go ~(l1:l2:ls) = do putStrLn l1+ if ('^' `elem` l1)+ then do+ putStrLn l2+ else go (l2:ls) go =<< fmap lines (hGetContents e)-- - hPutStr i $ unlines $ preamble- hPutStrLn i (mkplots cols ++ concatMap ((',':) . show) hlines)- mapM_ (\col -> do {hPutStr i . unlines . map show . snd $ col; hPutStrLn i "e"}) cols-+ genplot i preamble cols hlines hClose i- x <- waitForProcess p case x of ExitSuccess -> return () ExitFailure j -> hPutStrLn stderr (errmsg++show j) >> return () where errmsg = "'gnuplot' failed with exit code "++gnudat :: (Show i, Num i) => [String] -> [(String,[(Int,i)])] -> [Double] -> IO ()+gnudat = genplot stdout++-- | Write gnuplot instructions to a handle (todo: make pure)+genplot :: (Show i, Num i) => Handle -> [String] -> [(String,[(Int,i)])] -> [Double] -> IO ()+genplot i preamble cols hlines = do+ hPutStr i $ unlines $ preamble+ let show' (x,y) = show x ++ "\t" ++ show y+ hPutStrLn i (mkplots cols ++ concatMap ((',':) . show) hlines)+ mapM_ (\col -> do {hPutStr i . unlines . map show' . snd $ col; hPutStrLn i "e"}) cols -- | generate the plot command mkplots :: [(String,[a])] -> String
src/Options.hs view
@@ -5,7 +5,7 @@ import System.Console.CmdArgs import Control.Monad (when) import System.Directory (getTemporaryDirectory)-import Foreign (unsafePerformIO)+import System.IO.Unsafe (unsafePerformIO) tmpdefault :: FilePath tmpdefault = unsafePerformIO getTemporaryDirectory@@ -38,7 +38,8 @@ setOutputFormat o | null (format o) && null (outfile o) = o | null (format o) = o { format = determineFormat $ outfile o }- | null (outfile o) = error "Please specify an output file (-o) when you specify format."+ | format o /= "plot" && null (outfile o) = error "Please specify an output file (-o) when you specify format."+ | format o == "plot" && not (null $ outfile o) = error "'plot' format isn't supported with '-o',\nand only writes to standard output." | otherwise = o where determineFormat fp =
+ src/Sort.hs view
@@ -0,0 +1,21 @@+module Sort where++import qualified Data.IntMap as M+import Bio.Sequence+import Data.List (foldl')++type SizeTable = M.IntMap Int++sizes :: FilePath -> IO SizeTable+sizes f = (freqs . map (fromIntegral . seqlength)) `fmap` readFasta f++sort :: SizeTable -> [Int]+sort = concatMap (\(x,c) -> replicate c (fromIntegral $ negate x)) . M.toAscList++-- equivalent to 'M.fromList . map (\x->(fromIntegral $ negate x,1))', +-- except for not blowing the stack+freqs :: [Int] -> SizeTable+freqs = foldl' ins M.empty . map fromIntegral . map negate+ where ins m x = case M.lookup x m of + Just v -> v `seq` M.insert x (v+1) m+ Nothing -> M.insert x 1 m