a50 (empty) → 0.2

raw patch · 7 files changed

+321/−0 lines, 7 filesdep +basedep +biodep +cmdargssetup-changed

Dependencies added: base, bio, cmdargs, containers, directory, process

Files

README +62/−0
Setup.hs +5/−0
a50.cabal +35/−0
src/A50.hs +60/−0
src/Blat.hs +63/−0
src/Gnuplot.hs +44/−0
src/Options.hs +52/−0

+ README view

@@ -0,0 +1,62 @@+a50 is a tool for comparing genome assemblies, providing a bit more+information than the usual numeric statistics, like N50.  For a quick+overview of the options, use 'a50 --help'.+++  General usage+  -------------++To compare assemblies, you need two or more fasta-formatted files+containing contigs.  Running 'a50' with the files as arguments+produces a plot with one curve per input.  On the x-axis are the+contigs, ordered by size, and on the y-axis is the corresponding+cumulative size.  Generally, a better assembly has a steeper curve+(big contigs) which ends early (few contigs) near the y-value+corresponding to the expected genome size.++The plot is generated using 'gnuplot', so the gnuplot executable must+be in your $PATH.  You can pass various information to gnuplot on the+commandline, specifically you can use -f to specify output format+('terminal' in gnuplot lingo), -o to specify an output file, and -e to+produce horizontal lines, e.g. at the expected genome size.++For example:++        a50 -t pdf -o asm.pdf asm1.fasta asm2.fasta -e 8e+8++will plot the assemblies given in asm1.fasta and asm2.fasta against an+expected genome size of 800Mb in PDF format to the file asm.pdf.  If+no output or format is specified, gnuplot will display the graph in a+window.  If an output file, but no format is specified, a50 will try+to determine the format from the file name extension.+++  Using EST or transcripts as a reference+  ---------------------------------------++A similar way to compare assemblies is to measure the gene coverage of+each contig.  To do this, you need to have the genes available,+typically in the form of raw or assembled ESTs.  By specifying a set+of transcripts using the -E option, e.g:++        a50 asm.fasta asm2.fasta -E ests.fasta++will plot curves with contigs of the assemblies ordered by size along+the x-axis as before, but the y-value will be the total amount of+transcript data mapped to the contigs.++The mapping process runs BLAT, and stores the resulting PSL files in a+tempdir, typically /tmp, but overriden with the $TMPDIR variable or+specified with -T.  Only the best hit for each EST is retained, so+exons in other contigs are not counted.+ ++  Comparison to other measures+  ----------------------------++You can read various other measures off the graph fairly easily.+Total assembly size is the y-value at which the curve ends, and the+number of contigs in the assembly is the corresponding x-value.  An+assembly with greater average contig size will end to the left and/or+above of an assembly with shorter average contigs.  N50 is the slope+of the graph at the y-value corresponding to half the genome size.

+ Setup.hs view

@@ -0,0 +1,5 @@+#!/usr/bin/env runhaskell+ +import Distribution.Simple+main = defaultMain+

+ a50.cabal view

@@ -0,0 +1,35 @@+Name:           a50+Version:        0.2+License:        GPL++Author:         Ketil Malde+Maintainer:     Ketil Malde <ketil@malde.org>++Category:       Bioinformatics+Synopsis:       Compare genome assemblies+Description:    a50 - a simple tool for graphing genome coverage and fragmentation.+                .+                Reads files of contigs, and compares them by plotting each as a line in a graph. +		The x-axis represents contig number, the y-axis represents total (cumulative) size.+		An ideal assembly contains a few, large contigs, thus this curve should rise steeply, and +		stop early (but at the expected genome size).  Conversely, a poor assembly consisting of+		many small fragments will have a less steep curve extending far to the right.+		.+		The graphs produced by a50 gives a simple and easy to grasp comparison between assemblies,+		and yet produces a more detailed and informative view than the usual metrics like largest+		contig size or N50.+		.+                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+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

+ src/A50.hs view

@@ -0,0 +1,60 @@+module Main where++import Bio.Sequence+import Data.List (intersperse, foldl')++import qualified Data.IntMap as M++import Gnuplot+import Blat+import Options++main :: IO ()+main = do +  opts <- getArgs+  ss <- mapM sizes $ inputs opts+  case estref opts of +    "" -> do+      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)++zipLists :: [String] -> [[Int]] -> String+zipLists fs ss = unlines ((concat $ map ("#\t"++) fs) : go (map (++repeat 0) ([1..]:ss)))+  where go :: [[Int]] -> [String]+        go xs = let hs = map head xs+                    ts = map tail xs+                in if all (==0) $ tail hs then []+                   else (concat $ intersperse "\t" $ map myshow hs) : go ts++myshow :: Integral 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+  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+        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

@@ -0,0 +1,63 @@+--  blat DB QRY OUT -t= -d= -minScore= -extendThroughN+--  blat P_2010_10_11_10_15_26_runAssembly/454AllContigs.fna ../EST/lakselus.fasta -t=dna -q=rna tmp.psl -minScore=50++module Blat (module Bio.Alignment.PSL, runBlat, gen_result, interleave) where++import Bio.Alignment.PSL+import Control.Arrow (second)+import Data.List (sortBy)+import System.Directory (doesFileExist, findExecutable)+import Control.Monad (when)+import System.Process (runCommand, waitForProcess) +import System.Exit (ExitCode(..))++blat_cmd :: String+blat_cmd = "blat -t=dna -q=rna -minScore=50 -extendThroughN "++-- run blat+runBlat :: FilePath -> FilePath -> FilePath -> IO [PSL]+runBlat tmpdir asm ests = do+  let pslfile = tmpdir++"/"++basename asm++"_vs_"++basename ests++".psl"+      basename = reverse . takeWhile (/='/') . reverse+  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]+    e <- waitForProcess blat+    case e of ExitSuccess -> return ()+              _ -> error $ show e+  readPSL pslfile++-- this will only get you the sequences that have a match!+gen_result :: [PSL] -> [(Int,Int)]+gen_result = coverage . order . pslbest++interleave :: Integral i => [i] -> [(i,i)] -> [i]+interleave sz [] = map (const 0) sz+interleave (sz:szs) covs@((s1,c1):cs) +  | sz > s1   = 0 : interleave szs covs+  | sz == s1  = c1 : interleave szs cs++-- if this fails, something is wrong with my assumptions!+  | otherwise = error ("interleave failed: "++show (take 10 (sz:szs)) ++" "++show (take 10 covs))+interleave [] rs = error ("interleave failed - leftover contigs: "++show (take 10 rs)++"...")++coverage :: [PSL] -> [(Int,Int)]+coverage = map (second (sum . map (sum . blocksizes))) . tgroup+  where tgroup [] = []+        tgroup (p1:ps) = let (this,rest) = span ((tname p1 ==) . tname) ps+                         in (tsize p1, p1:this) : tgroup rest++order :: [PSL] -> [PSL]+order = sortBy (comparing (\c -> (negate $ tsize c,tname c)))++pslbest :: [PSL] -> [PSL]+pslbest = map (last . sortBy (comparing match)) . qgroup+  where +    qgroup [] = []+    qgroup (p1:ps) = let (this,rest) = span ((qname p1 ==) . qname) ps+                     in (p1:this) : qgroup rest++comparing :: (Ord a) => (t -> a) -> t -> t -> Ordering+comparing f a b = compare (f a) (f b)

+ src/Gnuplot.hs view

@@ -0,0 +1,44 @@+module Gnuplot where    ++import System.Process+import System.Exit+import System.IO+import Data.List (intersperse)++import Control.Monad (when)+import System.Directory (findExecutable)+import Control.Concurrent (forkIO)++gnuplot :: Num i => [String] -> [(String,[i])] -> [Int] -> IO ()+gnuplot preamble cols hlines = do+  fe <- findExecutable "gnuplot"+  when (fe == Nothing) $ error "Couldn't find the 'gnuplot' executable - aborting"+  (i,o,e,p) <- runInteractiveCommand "gnuplot -persist"+  _ <- 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)+    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++  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 "++-- | generate the plot command+mkplots :: [(String,[a])] -> String+mkplots cs = "plot " ++ (concat . intersperse "," . map mkp $ cs)+  where mkp :: (String,a) -> String+        mkp (s,_) = "'-' with lines title '"++s++"'"+

+ src/Options.hs view

@@ -0,0 +1,52 @@+{-# LANGUAGE DeriveDataTypeable #-}++module Options where++import System.Console.CmdArgs+import Control.Monad (when)+import System.Directory (getTemporaryDirectory)+import Foreign (unsafePerformIO)++tmpdefault :: FilePath+tmpdefault = unsafePerformIO getTemporaryDirectory++data Opt = Opt { outfile  :: FilePath+               , format :: String+               , expect   :: [String]+               , inputs   :: [FilePath]+               , estref   :: FilePath+               , tmpdir   :: FilePath+               } deriving (Typeable, Data, Show, Eq)++myopt :: Opt+myopt = Opt +  { outfile = def &= help "Output file, if applicable" &= typFile+  , format  = def &= help "Gnuplot output format (a.k.a. 'terminal')"+  , expect  = []  &= help "Expected genome size"+  , inputs  = def &= args &= typFile+  , estref  = def &= help "Reference transcripts" &= typFile &= name "E"+  , tmpdir  = tmpdefault &= help "Set temporary directory" &= typDir &= name "T"+  } &= summary "a50 - compare genome assemblies" &= program "a50"++getArgs :: IO Opt+getArgs = do+  o <- setOutputFormat `fmap` cmdArgs myopt+  when (null $ inputs o) $ error "Please specify one or more input files!"+  return o+  +setOutputFormat :: Opt -> Opt+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."+  | otherwise       = o+    where +      determineFormat fp =+        let ext = reverse . takeWhile (/='.') . reverse+        in case ext fp of +          "png" -> "png"+          "jpg" -> "jpg"+          "ps"  -> "postscript eps"+          "svg" -> "svg"+          "pdf" -> "pdf"+          _ -> error "Couldn't determine format from output file name.\nPlease specify format with -f."