diff --git a/README b/README --- a/README +++ b/README @@ -15,11 +15,12 @@ USAGE - rselect n [m] input.seq + rselect [-r] n [m] input.seq Selects n sequences from the file input.seq. If the optional m is given, this limits the selection to happen only from the first - m sequences in the file, which may be more efficient. + m sequences in the file, which may be more efficient. If -r is given, + the sequences will be reoriented randomly. The selected sequences are written to standard output, so you probably want to redirect them to a file. diff --git a/estreps.cabal b/estreps.cabal --- a/estreps.cabal +++ b/estreps.cabal @@ -1,5 +1,5 @@ Name: estreps -Version: 0.1 +Version: 0.3 License: GPL License-file: LICENSE Author: Ketil Malde @@ -7,12 +7,13 @@ Stability: Beta Category: Bioinformatics Synopsis: Repeats from ESTs -Description: rselect - select a random set of sequences from a FASTA file. - reps - extract exact k-word repeats based that occur in - sequences grouped in different clusters. +Description: * rselect - select a random set of sequences from a FASTA file, optinally + with random orientation (forward/reverse complement). + * reps - extract exact k-word repeats based that occur in + sequences grouped in different clusters. . The Darcs repository is at: . -Homepage: http://malde.org/~ketil/ +Homepage: http://blog.malde.org/ Tested-With: GHC==6.8.2 Build-Type: Simple diff --git a/src/RSelect.lhs b/src/RSelect.lhs --- a/src/RSelect.lhs +++ b/src/RSelect.lhs @@ -8,26 +8,28 @@ module Main where -import Data.Set +import Data.Set hiding (partition) import System.Environment (getArgs) import System.Random import System.IO +import Data.List (isPrefixOf,partition) -import Bio.Sequence.Fasta +import Bio.Sequence.Fasta (countSeqs) +import Bio.Sequence import Unigene main :: IO () main = do - as <- getArgs - case as of [m,f] -> select (read m) Nothing f - [m,n,f] -> select (read m) (Just (read n)) f - _ -> error "Usage: rselect n [m] input.seq > output.seq" + (opts,as) <- partition (isPrefixOf "-") `fmap` getArgs + case as of [m,f] -> select opts (read m) Nothing f + [m,n,f] -> select opts (read m) (Just (read n)) f + _ -> error "Usage: rselect [-r] n [m] input.seq > output.seq" -- | Select m out of maybe n sequences from f -- Not specifying a maximum increases memory consumption -select :: Integer -> Maybe Integer -> FilePath -> IO () -select m mn f = do +select :: [String] -> Integer -> Maybe Integer -> FilePath -> IO () +select os m mn f = do seed <- newStdGen hPutStrLn stderr ("Random seed: "++show seed) ss <- ugRead f @@ -35,14 +37,21 @@ n' <- countSeqs f return (fromIntegral n') Just n' -> return n' - let ds = if m>n then error "rselect: m is larger than n, that won't work." + + rds <- randoms `fmap` newStdGen :: IO [Bool] + let randomize_dirs = if "-r" `elem` os then zipWith maybe_reverse rds else id + + ds = if m>n then error "rselect: m is larger than n, that won't work." else if m==n then take (fromIntegral n) $ repeat 1 else if m<=n `div` 2 then diffs $ pick m $ randomRs (1,n) seed else diffs $ invert n 1 $ pick (n-m) $ randomRs (1,n) seed - out = Prelude.filter (not.Prelude.null) $ choose ds ss + out = Prelude.filter (not.Prelude.null) $ Prelude.map randomize_dirs $ choose ds ss n `seq` hPutStrLn stderr ("N = "++show n++", M = "++show m) - mapM_ (\o -> putStrLn "# cluster delimiter" >> hWriteFasta stdout o) out + mapM_ (\o -> putStrLn "# cluster delimiter" >> hWriteFasta stdout o) out +maybe_reverse :: Bool -> Sequence -> Sequence +maybe_reverse b = if b then revcompl else id + -- | choose elements according to a list of deltas choose :: Show a => [Integer] -> [[a]] -> [[a]] choose = choose' [] @@ -62,6 +71,7 @@ diffs' (a:b:cs) = if b<=a then error ("diffs needs a monotonically increasing sequence"++show (take 10 (a:b:cs))) else (b-a) : diffs' (b:cs) diffs' [_] = [] + diffs' [] = [] -- | convert a list to the list of missing values invert :: Integer -> Integer -> [Integer] -> [Integer] diff --git a/src/Repeats.lhs b/src/Repeats.lhs --- a/src/Repeats.lhs +++ b/src/Repeats.lhs @@ -33,10 +33,11 @@ import Data.Maybe import System.Environment (getArgs) -import Util +-- import Util -- type Word = Int +main :: IO () main = do args <- getArgs when (length args < 2 || length args > 3) diff --git a/src/Unigene.hs b/src/Unigene.hs --- a/src/Unigene.hs +++ b/src/Unigene.hs @@ -16,10 +16,11 @@ ugRead f = do s <- B.readFile f let ls = dropWhile B.null $ B.lines $ s - when (B.head (head ls) /= '#') (hPutStrLn stderr ("Warning: '"++f++"' does not look like Unigene format")) + when (B.head (head ls) /= '#') (hPutStrLn stderr ("'"++f++"' does not look like Unigene format")) return (clusters $ ls) clusters :: [B.ByteString] -> [[Sequence]] -clusters (l:ls) = let ls' = if B.head l == '#' then (l:ls) else (B.pack "# dummy UG hdr" :l:ls) - ss = groupBy (const (('#' /=) . B.head)) $ filter (not . B.null) ls' - in map (mkSeqs . tail) ss +clusters (l:ls) = if B.head l == '#' + then map (mkSeqs . tail) $ groupBy (const (('#' /=) . B.head)) $ filter (not . B.null) (l:ls) + else [mkSeqs (l:ls)] +clusters [] = [] diff --git a/src/Util.hs b/src/Util.hs --- a/src/Util.hs +++ b/src/Util.hs @@ -3,10 +3,11 @@ import System.IO import System.IO.Unsafe +countIO :: Int -> String -> [a] -> IO [a] countIO step msg xs = sequence $ map unsafeInterleaveIO $ ct 0 xs - where ct s xs = let (a,b) = splitAt (step-1) xs - next = s+step - output c = hPutStr stderr ((take (length (show c)+2+length msg) $ repeat ' ')++'\r':msg++show c) >> hFlush stderr + where ct s xs' = let (a,b) = splitAt (step-1) xs' + next = s+step + output c = hPutStr stderr ((take (length (show c)+2+length msg) $ repeat ' ')++'\r':msg++show c) >> hFlush stderr in if null b then map return a ++ [] else map return a ++ [output next >> return (head b)] ++ ct next (tail b)