flower-0.6: src/FlowSelect.hs

{-# LANGUAGE ExistentialQuantification #-}

module Main where

import Prelude hiding (LT,GT)

import Bio.Sequence.SFF
import System.Environment (getArgs)
import System.Random

import Metrics

main :: IO ()
main = do
  args <- getArgs
  (input,output,myfilter) <- parseArgs args
  SFF h rs <- readSFF input
  c <- writeSFF' output (SFF h $ myfilter rs)
  putStrLn ("Wrote "++show c++" reads.")

type FilterSFF = [ReadBlock] -> [ReadBlock]

parseArgs :: [String] -> IO (FilePath,FilePath,FilterSFF)
parseArgs [e,i] = do
          f <- case words e of 
                 ["Rand",x] -> do
                              ps <- randomRs (0,1) `fmap` newStdGen
                              let t = read x :: Double
                              return (map snd . filter ((<t).fst) . zip ps)
                 ["Pick",_] -> undefined
                 _  -> return (filter (apply (read e) . getChars))
          return (i,"selected.sff", f)

parseArgs _ = error "Usage: fselect <expression> <input.sff>"

-- | This structure represents selection parameters for one read
data Characteristics = Ch { k2, ee :: Double -- k-square, expected errors
                          , ns, len, tlen :: Int -- lenght, trimmed length
                          }

getChars :: ReadBlock -> Characteristics
getChars rb = let rh = read_header rb 
              in Ch { k2 = (/100) $ fromIntegral $ quals $ flowgram rb
                    , ee = 0 -- fixme!
                    , ns = n_count rb
                    , len  = fromIntegral $ num_bases rh
                    , tlen = fromIntegral $ clip_qual_right rh - clip_qual_left rh + 1}

data FilterFunction = forall a . Ord a => LT (Characteristics -> a) a
                    | forall a . Ord a => GT (Characteristics -> a) a

instance Show FilterFunction where show _ = "<filterfunction>"

instance Read FilterFunction where
    readsPrec _ str = case words str of 
                        (c:"k2":rest) -> [((lookupO c) k2 x,r) | (x,r) <- (reads $ unwords rest)]
                        (c:"ee":rest) -> [((lookupO c) ee x,r) | (x,r) <- (reads $ unwords rest)]
                        (c:"len":rest) -> [((lookupO c) len x,r) | (x,r) <- (reads $ unwords rest)]
                        (c:"tlen":rest) -> [((lookupO c) tlen x,r) | (x,r) <- (reads $ unwords rest)]
                        (c:"ncount":rest) -> [((lookupO c) ns x,r) | (x,r) <- (reads $ unwords rest)]
                        _ -> error ("Couldn't parse FilterFunction: "++take 100 str)

lookupO :: Ord a => String -> (Characteristics -> a) -> a -> FilterFunction
lookupO "LT" = LT
lookupO "GT" = GT
lookupO x = error ("FilterFunction must be either LT or GT, was "++take 100 x)


data Filter = Func FilterFunction 
            | And Filter Filter
            | Or Filter Filter
            | Not Filter 
              deriving Show

-- Okay, so we should really return/expect all parses here.
instance Read Filter where
    readsPrec _ str = readParen False p str
        where p s = case words s of 
                      "And":rest -> let [(a,r)] = reads (unwords rest)
                                        [(b,c)] = reads r
                                    in [(And a b,c)]
                      "Or":rest -> let ((a,r):_) = reads (unwords rest)
                                       ((b,c):_) = reads r
                                    in [(Or a b,c)]
                      "Not":rest -> let ((a,r):_) = reads (unwords rest)
                                    in [(Not a,r)]
                      "Func":rest -> let ((a,r):_) = reads (unwords rest)
                                     in [(Func a,r)]
                      _ -> [] -- error ("Couldn't parse "++take 100 s)

-- myFilter = LT k2 1.5

eval :: FilterFunction -> Characteristics -> Bool
eval (LT f a) c = f c < a
eval (GT f a) c = f c > a

apply :: Filter -> Characteristics -> Bool
apply (Func f)    = eval f 
apply (And f1 f2) = \r -> apply f1 r && apply f2 r
apply (Or  f1 f2) = \r -> apply f1 r || apply f2 r
apply (Not f)     = \r -> not (apply f r)