a50-0.2: src/A50.hs

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