diff --git a/README b/README
new file mode 100644
--- /dev/null
+++ b/README
@@ -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.
diff --git a/Setup.hs b/Setup.hs
new file mode 100644
--- /dev/null
+++ b/Setup.hs
@@ -0,0 +1,5 @@
+#!/usr/bin/env runhaskell
+ 
+import Distribution.Simple
+main = defaultMain
+
diff --git a/a50.cabal b/a50.cabal
new file mode 100644
--- /dev/null
+++ b/a50.cabal
@@ -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
diff --git a/src/A50.hs b/src/A50.hs
new file mode 100644
--- /dev/null
+++ b/src/A50.hs
@@ -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
diff --git a/src/Blat.hs b/src/Blat.hs
new file mode 100644
--- /dev/null
+++ b/src/Blat.hs
@@ -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)
diff --git a/src/Gnuplot.hs b/src/Gnuplot.hs
new file mode 100644
--- /dev/null
+++ b/src/Gnuplot.hs
@@ -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++"'"
+
diff --git a/src/Options.hs b/src/Options.hs
new file mode 100644
--- /dev/null
+++ b/src/Options.hs
@@ -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."