korfu 0.1.1 → 0.2
raw patch · 4 files changed
+312/−5 lines, 4 filesdep ~basedep ~bio
Dependency ranges changed: base, bio
Files
- korfu.cabal +3/−2
- src/Blastx.hs +232/−0
- src/Korfu.hs +3/−3
- src/Options.hs +74/−0
korfu.cabal view
@@ -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
+ src/Blastx.hs view
@@ -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)
src/Korfu.hs view
@@ -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
+ src/Options.hs view
@@ -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+