diff --git a/korfu.cabal b/korfu.cabal
--- a/korfu.cabal
+++ b/korfu.cabal
@@ -1,5 +1,5 @@
 Name:           korfu
-Version:        0.1.1
+Version:        0.2
 License:        GPL
 Category:	Bioinformatics
 Author:         Ketil Malde
@@ -15,7 +15,8 @@
 Build-Type:	Simple
 
 Executable:     korfu
-Build-Depends:  base>=3&&<4.2, bio >= 0.3 && < 0.4, haskell98, simpleargs, bytestring
+Build-Depends:  base>=3 && <5, bio >= 0.4, haskell98, simpleargs, bytestring
 Hs-Source-Dirs: src
 Main-Is:        Korfu.hs
+Other-Modules:	Options, Blastx
 Ghc-Options:    -Wall -fglasgow-exts -O2 -funbox-strict-fields
diff --git a/src/Blastx.hs b/src/Blastx.hs
new file mode 100644
--- /dev/null
+++ b/src/Blastx.hs
@@ -0,0 +1,232 @@
+{- | 
+   parse blastX and try to determine correct reading frame
+   
+   Each blastx match (HSP) adds its bitscore? /length? to the 
+   first position of each codon.
+
+   Todo: smooth (-1,2,-1), then add -k for past stop codons, -k,+k for aug/kozaks
+
+   Use DP to determine optimal path through "frame"-space.
+   What is frame-shift cost?  Direction-shift? Termination?
+
+   Bonus points:
+     extrapolated start/stop
+     Poly-A, Kozak, STP, AUG, ...
+
+   Problem: protein switch: follow good/strong protein in one frame, then switch to a different,
+   weaker match to (artificially) extend the frame.
+
+-}
+module Blastx (Weights, Frame
+              , frames, select_dir -- , indet
+              ,genFasta, genCDS
+              ) where
+
+import Bio.Sequence
+import Bio.Alignment.BlastFlat
+
+import Data.List (groupBy,maximumBy)
+import Text.Printf
+import qualified Data.ByteString.Lazy.Char8 as B
+
+import qualified Debug.Trace (trace)
+trace :: String -> a -> a
+trace m = id -- Debug.Trace.trace m -- switch for debugging output
+
+-- todo: use ST/U/Arrays for efficiency
+type Weights = [Double]
+type Frame   = (Double,[Int])
+
+genFasta :: (Sequence Nuc, (Frame,Frame)) -> [Sequence Nuc]
+genFasta (s,frames@((sc_f,fs),(sc_r,rs))) = 
+    let d = select_dir s frames
+        score = printf "%s: %.2f/%.2f CDS: " (show d) sc_f sc_r
+        cds_f = if sc_f > 0 then printf "%d..%d" (head fs) (last fs) else "-"
+        cds_r = if sc_r > 0 then printf "%d..%d" (last rs) (head rs) else "-"
+    in -- trace (show sc_f++" : "++show fs++"\n"++show sc_r++" : "++show rs++"\n") $
+       return $ tag ((if d==Rev then revcompl else id) s) (score++cds_f++"/"++cds_r) 
+
+genCDS :: (Sequence Nuc, (Frame,Frame)) -> [Sequence Amino]
+genCDS (s,frames@((sc_f,fs),(sc_r,rs))) = 
+    trace (show sc_f++" : "++show fs++"\n"++show sc_r++" : "++show rs++"\n") $
+    case select_dir s frames of
+                  Indet -> []
+                  Fwd -> [cds_fwd frames s]
+                  Rev -> [cds_rev frames s]
+--                  Vague -> [cds_fwd fs s,cds_rev fs s]
+
+-- todo: appears to give correct results, but...
+cds_fwd, cds_rev :: (Frame,Frame) -> Sequence Nuc -> Sequence Amino
+cds_fwd ((score,is),_) s = tag (Seq (seqlabel s) sdata Nothing) 
+   (printf "FWD trans: %.2f CDS: %d..%d l=%d" score (last is) (head is) (seqlength s))
+   where sdata = toIUPAC $ concatMap (take 1 . translate s . fromIntegral) $ reverse is
+
+cds_rev (_,(score,is)) s = trace (show (sl-last is)++".."++show (sl-head is)++": "++show is) $
+   tag (Seq (seqlabel s) sdata Nothing)
+   (printf "REV trans: %.2f CDS: %d..%d l=%d" score (sl-head is) (sl-last is) sl)
+    where sl = fromIntegral $ seqlength s
+          sdata = toIUPAC $ concatMap (take 1 . translate (revcompl s) . fromIntegral) $ reverse is
+
+frames :: [Sequence Nuc] -> [BlastFlat] -> [(Sequence Nuc,(Frame,Frame))]
+frames ss bfs = joins ss (map merge $ partitions bfs)
+    where merge (p:ps) = let (wf,wr) = mergeW $ map genweights (p:ps)
+                         in (query p, (wf,wr))
+
+joins :: [Sequence Nuc] -> [(SeqData,(Weights,Weights))] -> [(Sequence Nuc,(Frame,Frame))]
+joins ss [] = map (\s -> (s,no_frame)) ss
+joins (s:ss) all@((fs,ws):rest) 
+    | fs /= seqheader s   = (s,no_frame) : joins ss all
+    | otherwise           = let wf' = zipWith (+) (fst ws) (patterns 10 s)
+                                -- this is WRONG:
+                                wr' = zipWith (+) (snd ws) (patterns 10 $ revcompl s)
+                            in  (s,(opt_score 10 wf',opt_score 10 wr')) : joins ss rest
+
+no_frame = ((0,[]),(0,[]))
+
+-- | Partition BlastFlats by query name
+partitions :: [BlastFlat] -> [[BlastFlat]]
+partitions = groupBy ((==) `on` query)
+    where f `on` g = \x y -> f (g x) (g y)
+
+{-- RECURRENCE
+
+opt_score (w1:w2:w3:w4:ws) = 
+    max | w1 + opt_score (w4:ws)
+        | w0 - p + opt_score (w3:w4:ws)
+        | w2 - p + opt_score ws
+
+Perhaps easier:
+  s1,s2,s3 = w1,w2,w3
+  si = wi+max (w(i-3),-p+w(i-2),-p+w(i-4))
+
+Todo: add codon bias
+      add amino comp(?)
+      add slight positive weight overall to prefer long ORFs?
+-}
+
+-- TODO: xxx-AUG-yyy- -- xxx is -p , AUG is +p-epsilon? (or just hit-bits?)
+-- starting at AUG gives bonus, passing (including) it is neutral
+-- but...we shouldn't encourage frame shifts at the AUG.
+-- xxx-STP-yyy  -- xxx is -p STP is +1.5p, yyy is -p (+eps - rather force two frame shifts)
+-- should be better to frame-shift before the STP than to include it...
+
+patterns :: Double -> Sequence Nuc -> Weights
+patterns p s = zipWith (+) (codon_adjust p s) (polyA_adjust s)
+
+-- | Add START+STOP bonus/penalties.
+--   Currently, a severe penalty for STP codons
+codon_adjust :: Double -> Sequence Nuc -> Weights
+codon_adjust p s = scores (f1,f2,f3)
+    where
+      [f1,f2,f3] = map (translate s) [0,1,2]
+      scores ([],[],[]) = []
+      scores (a:as,bs,cs) = (case a of STP -> -2*p
+                                       _   ->    0) : scores (bs,cs,as)
+
+-- increase scores before poly-A, decrease it afterwards
+polyA_adjust :: Sequence Nuc -> Weights
+polyA_adjust s = case pApos s of Just i' -> let i = fromIntegral i' 
+                                            in replicate i positive ++ replicate (sl-i) negative
+                                 Nothing -> replicate sl zero
+    where positive =  0.1
+          negative = -0.05
+          zero     =  0.01
+          sl       = fromIntegral $ seqlength s
+
+-- | Find a position whose suffix has > t As within w distance
+pApos :: Sequence Nuc -> Maybe Offset
+pApos s@(Seq _ d _) = go w0 w
+    where w0 = fromIntegral . length . filter isA . take (fromIntegral w) . toStr $ d
+          isA c = c == 'a' || c == 'A'
+          go :: Offset -> Offset -> Maybe Offset
+          go count pos 
+              | pos >= l  = Nothing
+              | count > t = Just pos
+              | otherwise = let next = if isA (s!pos) then 1 else 0
+                                prev = if isA (s!(pos-w)) then 1 else 0
+                            in go (count+next-prev) (pos+1)
+          l = seqlength s
+          t = 20 -- threshold
+          w = 25 -- window size
+
+-- | Calculate optimal score, and include the coordinates for the codons.
+--   Uses a kind of DP, retaining only four values (but all coords for best hit so far).
+opt_score :: Double -> [Double] -> Frame
+opt_score p (w1:w2:w3:w4:ws) = osc 4 (maximum' [i0,i1,i2,i3]) (i0,i1,i2,i3) ws
+    where (i0,i1,i2,i3) = ((w1,[0]),(w2,[1]),(w3,[2])
+                          , if w1 >= w2-p then (w4 + w1,[3,0])
+                            else (w4+w2-p,[3,1]))
+
+          -- osc params: curidx maxidx cache rest
+          osc i m (s0,s1,s2,s3) (w:ws) = 
+              (if (i-1) `rem` 3==0 then trace (printf "%d :\t%.1f\t%.1f\t%.1f\tmax %.1f,%d..%d" (i-1) (fst s0) (fst s1) (fst s2) (fst m) (head $ snd m) (last $ snd m)) else id) $
+              let (s,is) = maximum' [s1,(pm s0),(pm s2),(0,[])]
+              in osc (i+1) (if s>fst m then (s,is) else m) (s1,s2,s3,(w+s,i:is)) ws
+          osc _i m (_0,s1,s2,s3) [] = trace (printf "%d :\t%.1f\t%.1f\t%.1f\n" _i (fst s1) (fst s2) (fst s3)) $
+                                      maximum' [m,s1,s2,s3]
+
+          pm (s,is) = (s-p,is)
+          maximum' = maximumBy (compare `on` fst)
+          f `on` g = \x y -> f (g x) (g y)
+
+opt_score _ _ = (0,[]) -- less than four nucleotides is too short a sequence :-)
+
+{-
+-- | Subtract the opposing view from a set of weights (why does this matter?)
+convolve :: (Weights,Weights) -> (Weights,Weights)
+convolve = unzip . c . zip 
+    where c xs@(x:_) = c1 x : map c2 (takeWhile (>=2 . length) tails xs)
+          c1 ((a,b):(c,d):_) = (a-(b+c+d)/3, b-(a+c+d)/3)
+          c2 ((a,b):(c,d):(e,f):_) = let z = a+b+e+f 
+                                     in (c-(d+z)/5, d-(c+z)/5)
+          c2 z@[_,_] = c1 $ reverse z
+-}
+
+-- | Merge weights
+mergeW :: [(Weights,Weights)] -> (Weights,Weights)
+mergeW = foldr (sumW max) (repeat (-0.0001),repeat (-0.0001))
+
+sumW :: (Double -> Double -> Double) -> (Weights,Weights) -> (Weights,Weights) -> (Weights,Weights)
+sumW op (as,bs) (cs,ds) = (zipWith op as cs, zipWith op bs ds)
+
+-- | Generate appropriate weights for a single BLAST hits.
+--   BLAST coordinates are for the FWD strand only, starts on 1, and the hit is inclusive
+genweights :: BlastFlat -> (Weights,Weights)
+genweights bf = let Frame s _ = aux bf
+                    matchlen = q_to bf - q_from bf + 1
+                    b = 3*bits bf/fromIntegral matchlen
+                    res = replicate (q_from bf-1) 0 ++ 
+                          take matchlen (case s of Plus -> cycle [b,0,0]
+                                                   Minus -> cycle [0,0,b]) ++
+                          replicate (qlength bf - q_to bf) 0
+                in case s of Plus -> (res, replicate (qlength bf) 0)
+                             Minus -> (replicate (qlength bf) 0, reverse res)
+
+-- use appendHeader from new biolib (todo!)
+tag :: Sequence t -> String -> Sequence t
+tag (Seq h d _) t = (Seq (B.append h (B.pack (" "++t))) d Nothing)
+
+{-
+-- | What to do if no data is available
+indet :: Sequence -> Sequence
+indet s = tag s "INDET"
+
+fwd, rev, vague :: (Frame,Frame) -> Sequence -> Sequence
+fwd ((f,(i:is)),(r,_)) s = tag s 
+                       (printf "FWD: %.2f/%.2f CDS: %d..%d" f r (last (i:is)) i)
+rev ((f,_),(r,(i:is))) s = tag (revcompl s) 
+                       (printf "REV: %.2f/%.2f CDS: %d..%d" r f (sl-i) (sl-last (i:is)))
+    where sl = fromIntegral $ seqlength s
+vague ((f,(fi:fis)),(r,(ri:ris))) s = 
+    tag s (printf "WEAK: %.2f/%.2f CDSs: %d..%d/%d..%d" f r  (last (fi:fis)) fi (last (ri:ris)) ri)
+vague x s = error $ show x -- todo: happens if only score in one dir, but less than threshold
+-}
+
+-- | Orient based on BLAST data
+select_dir :: Sequence Nuc -> (Frame,Frame) -> Dir
+select_dir s scs@((sc_f,_if),(sc_r,_ir))
+    | sc_f-5 > sc_r = Fwd
+    | sc_r-5 > sc_f = Rev
+    | otherwise     = Indet
+
+data Dir = Indet | Fwd | Rev deriving (Show,Eq)
diff --git a/src/Korfu.hs b/src/Korfu.hs
--- a/src/Korfu.hs
+++ b/src/Korfu.hs
@@ -16,7 +16,7 @@
 main :: IO ()
 main = do
   opts <- getOptions
-  ss <- readFasta $ fasta opts 
+  ss <- map castToNuc `fmap` (readFasta $ fasta opts)
   if gen_all opts 
      then do 
        case cds_out opts of
@@ -42,10 +42,10 @@
        when (isJust g) $ hClose (fromJust g)
        when (isJust h) $ hClose (fromJust h)
 
-all_frames :: Sequence -> [Sequence]
+all_frames :: Sequence Nuc -> [Sequence Amino]
 all_frames s = 
     map (\i -> (Seq (fromStr (toStr (seqlabel s)++"_F+"++show i))  (toIUPAC $ translate s i) Nothing)) [0,1,2]
     ++ map (\i -> (Seq (fromStr (toStr (seqlabel s)++"_F-"++show i)) (toIUPAC $ translate (revcompl s) i) Nothing)) [0,1,2]
 
-process :: (a -> [Sequence]) -> Handle -> a -> IO ()
+process :: (a -> [Sequence t]) -> Handle -> a -> IO ()
 process fmt h = hWriteFasta h . fmt
diff --git a/src/Options.hs b/src/Options.hs
new file mode 100644
--- /dev/null
+++ b/src/Options.hs
@@ -0,0 +1,74 @@
+module Options where
+
+import System.Environment (getArgs)
+import System.Console.GetOpt
+import System (exitWith,ExitCode(..))
+import Control.Monad
+
+-- import Bio.Util (countIO)
+
+getOptions :: IO Options
+getOptions  = do as <- getArgs
+                 let (optf,non,err) = getOpt Permute options as
+                 when (not (null err && null non)) $ error (usage err)
+                 if null as
+                     then normalUsage
+                     else parseargs optf
+
+normalUsage :: IO a
+normalUsage = do
+  putStrLn $ usageInfo (usagemsg++synopsis) options
+  exitWith ExitSuccess
+
+usage :: [String] -> String
+usage errs = usageInfo (concat errs ++ usagemsg) options
+
+usagemsg, synopsis :: String
+usagemsg = "Usage: korfu -i <file> -x <file.xml> [-R <file>] [-C <file>]\n"
+synopsis = "Characterize ESTs using BLASTX hits and sequence heuristics.\n"
+
+data Options = Opts { redir_out, cds_out :: FilePath
+                    , default_fwd :: Bool
+                    , gen_all :: Bool
+                    , fasta,blastx :: FilePath }
+
+options :: [OptDescr (Options -> IO Options)]
+options = 
+    [ -- Output options
+      Option ['R'] ["redir-out"]     (ReqArg (\arg opt -> return (if gen_all opt then incompat
+                                                                  else opt { redir_out = arg })) "File") 
+                 "File for sequence output,\nreoriented to predicted 5' to 3' direction."
+    , Option ['C'] ["cds-out"] (ReqArg (\arg opt -> return opt { cds_out = arg }) "File") 
+                 "File for predicted CDS output."
+-- todo:
+--    , Option ['r'] ["reverse-default"]  (NoArg (\opt -> return opt { default_fwd = False }))
+--                 "Reverse sequences without good orientation predictions"
+
+    , Option [] ["all"] (NoArg (\opt -> return (if null (blastx opt) then opt { gen_all = True } 
+                                                else incompat))) "Generate all frames"
+
+    -- Input options
+    , Option ['i'] ["sequence-file"] (ReqArg (\arg opt -> return opt { fasta = arg }) "File")
+                 "File for sequence input (Fasta-format)."
+    , Option ['x'] ["blast-xml-file"] (ReqArg (\arg opt -> return (if gen_all opt then incompat
+                                                                   else opt { blastx = arg })) "File")
+                 "File for Blastx result input\n(XML-format, use 'blastall -m 7')."
+
+    -- Other options
+    , Option ['h'] ["help"] (NoArg (\_ -> normalUsage)) "Print usage information."
+    ]
+        where incompat = error "'--all' is incompatible with '-x' and '-R'."
+
+defaultopts :: Options
+defaultopts = Opts { redir_out = "" 
+                   , cds_out   = ""
+                   , default_fwd = True
+                   , gen_all = False
+                   , fasta = error "No sequence input specified"
+                   , blastx = ""
+                   }
+
+parseargs :: [Options -> IO Options] -> IO Options
+parseargs args = foldl (>>=) (return defaultopts) args >>= verify
+    where verify = return . id
+