flower 0.5 → 0.6
raw patch · 10 files changed
+173/−121 lines, 10 filesdep +mtldep ~biodep ~cmdargs
Dependencies added: mtl
Dependency ranges changed: bio, cmdargs
Files
- flower.cabal +19/−13
- src/FRecover.hs +2/−2
- src/FRename.hs +5/−6
- src/FlowSelect.hs +4/−3
- src/Flower.hs +77/−82
- src/Flowt.hs +13/−13
- src/Fork.hs +14/−0
- src/Metrics.hs +1/−1
- src/Options.hs +37/−0
- src/Print.hs +1/−1
flower.cabal view
@@ -1,5 +1,5 @@ Name: flower-Version: 0.5+Version: 0.6 License: GPL Cabal-Version: >= 1.6 Author: Ketil Malde@@ -10,7 +10,9 @@ Description: flower - FLOWgram ExtractoR tools . The flower executable reads files in SFF-format and produces various output, - including sequences with quality, or flowgram data in tabular format.+ including sequences with quality, or flowgram data in tabular format. By default,+ it outputs a textual representation of the data in the SFF-file, much like+ 'sffinfo' from Roche does. . The fselect executable extracts reads from SFF-files, generating a new SFF-file with a subset of the reads based on various criteria.@@ -18,38 +20,42 @@ Sometimes SFF files will appear to be corrupted, with all-zero blocks in the file. The frecover program ignores these and tries to resync with the file after an invalid region. This was likely a one-time bug in the 454 software, so this program- is probably not so useful any more.- .+ is probably not so useful any more.+ . The Darcs repository is at <http://malde.org/~ketil/biohaskell/flower>. -HomePage: http://malde.org/~ketil/biohaskell/flower-Build-Depends: bio >= 0.4.2, base >=3 && <5, array >= 0.1, bytestring >= 0.9.1, binary == 0.4.*, random, cmdargs <= 0.1.1, containers+HomePage: http://blog.malde.org/index.php/flower/+Build-Depends: bio >= 0.4.7, base >=3 && <5, array >= 0.1, bytestring >= 0.9.1, binary == 0.4.*, random, cmdargs >= 0.5, containers, mtl Build-Type: Simple-Tested-with: GHC==6.8.3+Tested-with: GHC==6.8.3, GHC==6.12.1 -- Data-files: README Executable: flower Main-Is: Flower.hs-Other-Modules: Print, Metrics+Other-Modules: Print, Metrics, Options, Fork Hs-Source-Dirs: src-Ghc-Options: +Ghc-Options: -Wall Executable: flowselect Main-Is: FlowSelect.hs Other-Modules: Metrics Hs-Source-Dirs: src+Ghc-Options: -Wall Extensions: ExistentialQuantification Executable: frecover Main-Is: FRecover.hs Hs-Source-Dirs: src+Ghc-Options: -Wall Executable: frename Main-Is: FRename.hs Hs-Source-Dirs: src+Ghc-Options: -Wall -Executable: flowt-Main-Is: Flowt.hs-Hs-Source-Dirs: src-Extensions: DeriveDataTypeable+Executable: flowt+Main-Is: Flowt.hs+Hs-Source-Dirs: src+Ghc-Options: -Wall -fno-warn-unused-do-bind+Extensions: DeriveDataTypeable
src/FRecover.hs view
@@ -3,7 +3,7 @@ import Bio.Sequence.SFF import System.Environment (getArgs) +main :: IO () main = mapM_ recoverFile =<< getArgs--recoverFile f = writeSFF (f++"_recovered") =<< recoverSFF f+ where recoverFile f = writeSFF (f++"_recovered") =<< recoverSFF f
src/FRename.hs view
@@ -11,6 +11,7 @@ import System.Environment (getArgs) import qualified Data.ByteString.Char8 as B +main :: IO () main = do fs <- getArgs if null fs then putStrLn "Usage: frename file1.sff [file2.sff ...]"@@ -23,9 +24,7 @@ (SFF h rs) <- readSFF f writeSFF ("r_"++f) (SFF h $ renameFrom current rs) go (current+num_reads h) fs- -renameFrom i rs = zipWith update [i..] rs- where update j r = let h = read_header r- rn = B.concat [read_name h, B.pack "_", B.pack (show j)]- in r { read_header = h { name_length = fromIntegral $ B.length rn,- read_name = rn }}+ renameFrom i rs = zipWith update [i..] rs+ where update j r = let h = read_header r+ rn = B.concat [read_name h, B.pack "_", B.pack (show j)]+ in r { read_header = h { name_length = fromIntegral $ B.length rn, read_name = rn }}
src/FlowSelect.hs view
@@ -27,7 +27,7 @@ ps <- randomRs (0,1) `fmap` newStdGen let t = read x :: Double return (map snd . filter ((<t).fst) . zip ps)- ["Pick",n] -> undefined+ ["Pick",_] -> undefined _ -> return (filter (apply (read e) . getChars)) return (i,"selected.sff", f) @@ -52,7 +52,7 @@ instance Show FilterFunction where show _ = "<filterfunction>" instance Read FilterFunction where- readsPrec i str = case words str of + 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)]@@ -60,6 +60,7 @@ (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)@@ -73,7 +74,7 @@ -- Okay, so we should really return/expect all parses here. instance Read Filter where- readsPrec i str = readParen False p str+ 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
src/Flower.hs view
@@ -3,20 +3,20 @@ module Main (main) where import Bio.Sequence.SFF-import Bio.Sequence.SFF_filters import Bio.Sequence.Fasta-import Bio.Sequence.FastQ-import Bio.Util (countIO)+import Bio.Sequence.FastQ (hWriteSangerQ)+-- import Bio.Util (countIO) import Print import Text.Printf -import System.IO (stdout)-import System.Environment (getArgs)+import System.IO (stdout, Handle, openFile, IOMode(..), hClose, hPutStrLn)+-- import System.Environment (getArgs)+ import Numeric (showFFloat) import Data.Char (toLower)-import Data.List (intersperse, partition)-import Data.ByteString.Char8 (pack,unpack,ByteString)+import Data.List (intersperse)+import Data.ByteString.Char8 (unpack,ByteString) import qualified Data.ByteString.Char8 as B import qualified Data.ByteString as B1 import qualified Data.ByteString.Lazy.Char8 as L@@ -25,53 +25,62 @@ import Data.Array.Unboxed import Data.Array.ST import Control.Monad.ST+import Control.Monad.State import Metrics+import qualified Options as O+import Options (Opts)+import Fork main :: IO () main = do- args <- getArgs- let (opts,files) = partition (\p -> case p of ('-':_) -> True; _ -> False) args+ opts <- O.getArgs+ when (null $ O.inputs opts) $ error "Please provide an input file!"+ SFF h rs <- readSFF (O.inputs opts)+ let acts = case buildActions opts h of + [] -> [dumpText] + as -> as+ writeInfo (O.info opts) h+ forkAndWait $ map ($ map (mkTrimmer opts) rs) acts - reader :: String -> IO SFF- reader f = if "-v" `elem` opts then do SFF h rs <- readSFF f - rs' <- countIO "reads: " "done" 20 rs- return (SFF h rs')- else readSFF f - writer :: SFF -> IO ()- writer = case filter (/="-v") opts of - ["-r"] -> \s@(SFF ch _) -> hWriteFasta stdout . trim_keys ch . sffToSequence $ s- ["-R"] -> writeFastaQual "flower.fasta" "flower.qual" . sffToSequence- ["-q"] -> hWriteFastQ stdout . sffToSequence- ["-f"] -> L1.putStrLn . L1.fromChunks . intersperse (B.pack "\n") . showflow- ["-h"] -> putStr . sffToHistogram- ["-H"] -> putStr . sffToHistogramClip- ["-i"] -> putStr . getHeader- ["-s"] -> summarize- [] -> dumpText- _ -> error ("Usage: flower -[f|h|H|i|q|r|R|s] <file.sff> [<file2.sff> ..]\n"- ++" -r output reads in Fasta format\n"- ++" -R output reads in Fasta format with associated .qual\n"- ++" (generates files instead of writing to <stdout>)\n"- ++" -q output in FastQ format\n"- ++" -f output the flowgram in tabular format\n"- ++" -h output a histogram table of flow values\n"- ++" -H output a histogram of flows after clipping\n"- ++" -i output header information\n"- ++" -s output a summary of each read"- )- writer `seq` mapM_ (\f -> writer =<< reader f) files+type Action = [ReadBlock] -> IO ()+type Trimmer = ReadBlock -> ReadBlock+type Info = Maybe FilePath --- trim keys if they match, eliminate reads that don't.-trim_keys _ [] = []-trim_keys ch (s:ss) = maybe id ((:) . id) (trimKey ch s) $ trim_keys ch ss+buildActions :: Opts -> CommonHeader -> [Action]+buildActions o ch = snd $ flip runState [] $ do+ on (O.fasta o) (\h -> hWriteFasta h . map rbToSequence) -- todo: add qual!+ on (O.fastq o) (\h -> hWriteSangerQ h . map rbToSequence)+ on (O.summarize o) (\h -> L1.hPut h . summarize)+ on (O.flowgram o) (\h -> L1.hPut h . L1.fromChunks . intersperse (B.pack "\n") . concatMap (showread ch))+ on (O.histogram o) (\h -> hPutStrLn h . showHist . histogram (B.unpack $ flow ch) . map flowgram) +on :: Maybe FilePath -> (Handle -> Action) -> State [Action] ()+on Nothing _ = return ()+on (Just f) act = modify $ (:) $ case f of + "-" -> act stdout+ _ -> \rs -> do h <- openFile f WriteMode+ act h rs+ hClose h++mkTrimmer :: Opts -> Trimmer+mkTrimmer o = case (O.trimKey o, O.trim o) of+ (True,True) -> error "Please specify only one of --trim and --trimkey"+ (True,False) -> \r -> trimFromTo 4 (num_bases $ read_header r) r+ (False,True) -> trim+ (False,False) -> id++writeInfo :: Info -> CommonHeader -> IO ()+writeInfo Nothing = const (return ())+writeInfo (Just "-") = putStr . getHeader+writeInfo (Just f) = writeFile f . getHeader+ -- ------------------------------------------------------------ -- No option - dump as text format -- -------------------------------------------------------------dumpText :: SFF -> IO ()-dumpText (SFF _ rs) = putStr . concat . map toText $ rs+dumpText :: [ReadBlock] -> IO ()+dumpText rs = putStr . concat . map toText $ rs where toText :: ReadBlock -> String toText r = concat [ gt, B.unpack (read_name rh), nl , maybe "" ((\s->info++s++nl) . formatRN) $ decodeReadName (read_name rh)@@ -79,7 +88,7 @@ in if lf /= 0 || rt /= 0 then adapter ++ show lf ++ sp++ show rt else "" , clip, show (clip_qual_left rh), sp, show (clip_qual_right rh), nl , flows, unwords $ map show $ flowgram r, nl- , index, unwords $ map show $ cumulative_index' r, nl+ , idx, unwords $ map show $ cumulative_index' r, nl , base, L.unpack (masked_bases' r), nl , qual, unwords $ map show $ L1.unpack (quality r), nl ]@@ -91,11 +100,11 @@ clip = " Clip: \t" adapter = " Adap: \t" flows = " Flows:\t"- index = " Index:\t"+ idx = " Index:\t" base = " Bases:\t" qual = " Quals:\t"- formatRN (ReadName (yr,mo,dy) (h,m,s) r x y) = - printf "%4d-%02d-%02d %02d:%02d:%02d R%d (%d,%d)" yr mo dy h m s r x y+ formatRN (ReadName (yr,mo,dy) (h,m,s) r' x y) = + printf "%4d-%02d-%02d %02d:%02d:%02d R%d (%d,%d)" yr mo dy h m s r' x y -- Get the bases, but lower case masked bases. Also in biolib's SFF.hs (sans prime),@@ -115,22 +124,22 @@ -- ------------------------------------------------------------ -- The -i option: Print header info -- -------------------------------------------------------------getHeader :: SFF -> String-getHeader (SFF h _) = unlines ["Index: \t" ++ show (index_offset h,index_length h)- ,"Num_reads:\t" ++ show (num_reads h)- ,"Num_flows:\t" ++ show (flow_length h)- ,"Key: \t" ++ unpack (key h)- ]+getHeader :: CommonHeader -> String+getHeader h = unlines ["Index: \t" ++ show (index_offset h,index_length h)+ ,"Num_reads:\t" ++ show (num_reads h)+ ,"Num_flows:\t" ++ show (flow_length h)+ ,"Key: \t" ++ unpack (key h)+ ] -- ---------------------------------------------------------- -- The -s option: Summarize each read on one line -- ---------------------------------------------------------- -- | Summarize each read on one line of output-summarize :: SFF -> IO ()-summarize (SFF _rh rs) = do- putStrLn "# name........\tdate......\ttime....\treg\ttrim_l\ttrim_r\tx_loc\ty_loc\tlen\tK2\ttrimK2\tncount\tavgQ\ttravgQ\tf:sgnt\tf:q20"- L1.putStrLn . toLazyByteString . mconcat . map sum1 $ rs+summarize :: [ReadBlock] -> L1.ByteString+summarize rs = do+ L.concat [ L.pack "# name........\tdate......\ttime....\treg\ttrim_l\ttrim_r\tx_loc\ty_loc\tlen\tK2\ttrimK2\tncount\tavgQ\ttravgQ\n"+ , toLazyByteString . mconcat . map sum1 $ rs] -- todo: date and time are usually constants! sum1 :: ReadBlock -> Builder@@ -138,35 +147,30 @@ nb = num_bases rh h = read_name rh tr = trim r+ tb, nl, q :: Builder+ tb = char '\t'+ nl = char '\n'+ q = char '?'+ (rndec1,rndec2) = case decodeReadName h of Just rn -> let ((y,m,d),reg,(hh,mm,ss)) = (date rn,region rn,time rn) in ([putDate y m d, putTime hh mm ss, putInt2 reg] ,[putInt (fromIntegral $ x_loc rn), putInt (fromIntegral $ y_loc rn)]) Nothing -> ([q,q,q],[q,q]) (qleft,qright) = (clip_qual_left rh, clip_qual_right rh)- avg_qual q = let l = fromIntegral (L1.length q)- in if l>0 then putFix 2 $ sum (map fromIntegral $ L1.unpack q) * 100 `div` l+ avg_qual qs = let l = fromIntegral (L1.length qs)+ in if l>0 then putFix 2 $ sum (map fromIntegral $ L1.unpack qs) * 100 `div` l else putFix 2 0 in mconcat $ intersperse tb ([fromByteString h] ++ rndec1 ++ [putInt (fromIntegral qleft), putInt (fromIntegral qright)] ++ rndec2 ++ [putInt (fromIntegral nb) , fromByteString (fi $ quals $ flowgram r), fromByteString (fi $ quals $ flowgram tr) , putInt (n_count r)- , avg_qual $ quality r, avg_qual $ quality tr- , putInt (flowToBasePos r $ sigint r), putInt (qual20 r)]) ++ [nl]--tb, nl, q :: Builder-tb = char '\t'-nl = char '\n'-q = char '?'+ , avg_qual $ quality r, avg_qual $ quality tr]) ++ [nl] -- ---------------------------------------------------------- -- The -f option: Output the sequence of flows, one flow per line -- ---------------------------------------------------------- --- | output a list of flows-showflow :: SFF -> [ByteString]-showflow (SFF h rs) = concatMap (showread h) rs- fi :: Flow -> ByteString fi f | f <= 9999 && f >= 0 = farray!f | otherwise = let (i,r) = f `divMod` 100 in B.pack (show i++"."++show r) @@ -187,9 +191,11 @@ in zipWith4 format [(1::Int)..] (maskFlows $ flow h) (flowgram rd) qgroups -- lower case based on the clip_qual values-mask rh p _ [] = [] -- qgroups are infinite+mask :: ReadHeader -> Int -> [[a]] -> [Char] -> [Char]+mask _ _ _ [] = [] -- qgroups are infinite mask rh p (q1:qs) (c:cs) = c' : mask rh (p+length q1) qs cs where c' = if fromIntegral p < clip_qual_left rh || fromIntegral p > clip_qual_right rh then toLower c else c+mask _ _ _ _ = error "internal error in 'mask'" zipWith4 :: (a -> b -> c -> d -> e) -> [a] -> [b] -> [c] -> [d] -> [e] zipWith4 f (a:as) (b:bs) (c:cs) (d:ds) = f a b c d : zipWith4 f as bs cs ds@@ -204,22 +210,11 @@ (q1,qrest) = splitAt (length iz+1) qs in q1 : qgroup irest qrest qgroup (i:is) qs = [] : qgroup (i-1:is) qs-+qgroup _ _ = error "internal error in 'qgroup'" -- ---------------------------------------------------------- -- The -h and -H options: Output a histogram of flow values -- -------------------------------------------------------------- | Generate a histogram of flow values from an SFF file-sffToHistogram, sffToHistogramClip :: SFF -> String-sffToHistogram (SFF h rs) = showHist . histogram (B.unpack $ flow h) . map flowgram $ rs-sffToHistogramClip (SFF h rs) = showHist . histogram (B.unpack $ flow h) . map clip_flowgram $ rs--clip_flowgram rd = let (l,r) = (fromIntegral $ clip_qual_left (read_header rd)-1, fromIntegral $ clip_qual_right (read_header rd))- ps = take r $ B1.unpack $ flow_index rd- p1 = fromIntegral $ sum $ take l ps- p2 = fromIntegral $ sum $ drop l ps- in take p2 $ drop p1 $ flowgram rd type Hist = UArray Flow Int
src/Flowt.hs view
@@ -20,9 +20,6 @@ import System.Console.CmdArgs -version :: String-version = "flowt v0.4, copyright 2009-2010 Ketil Malde"- type FingerPrints = IntSet type DupMap = IntMap [ReadBlock] @@ -32,6 +29,7 @@ splitRes (Then x rs) = let (ys,e) = splitRes rs in (x:ys,e) splitRes (EndWith e) = ([],e) +trim :: ReadBlock -> ReadBlock trim = id -- trimFromTo 4 10000 <- this trims to last base called position! data Options = O { thresh :: Double@@ -41,21 +39,21 @@ , input :: [FilePath] } deriving (Data,Typeable,Show,Eq) -modes :: Mode Options-modes = mode $ O { thresh = 50 &= text "similarity threshold"- , fplen = 20 &= text "fingerprint size"- , summarize = def &= empty "-" & text "output cluster summary"- , clusters = True &= text "output complete clusters"- , input = def &= args & typFile }- &= prog "flowt"- & text "Filter out reads from duplicate clones in 454 sequencing."+mymode :: Options+mymode = O { thresh = 50 &= help "similarity threshold"+ , fplen = 20 &= help "fingerprint size"+ , summarize = def &= {- empty "-" &= -} help "output cluster summary"+ , clusters = True &= help "output complete clusters"+ , input = def &= args &= typFile }+ &= program "flowt"+ &= summary "flowt v0.6 - filter out reads from duplicate clones in 454 sequencing." vlog :: Bool -> String -> IO () vlog v s = when v (hPutStr stderr s >> hFlush stderr) main :: IO () main = do- opts <- cmdArgs version [modes]+ opts <- cmdArgs mymode -- putStrLn $ show opts verb <- isLoud vlog verb "Building the fingerprint index"@@ -78,11 +76,12 @@ BC.unlines $ map (BC.unwords . map (read_name . read_header)) $ cg return () +gensum :: String -> IntMap [ReadBlock] -> IO () gensum f ds = write $ unlines $ map showcluster $ M.assocs ds where write = if f == "-" then putStrLn else writeFile f showcluster (k,v) = let a = averageflow v in printf "%16x" k++":\t" ++ show (length v) ++ unwords (map (BC.unpack . read_name . read_header) v) ++ "\n" ++ - unlines [let f = flowgram (trim x) in printf "%6.1f " (dist a $ map fromIntegral f) ++ concatMap (printf "%3d ") f | x <- v]+ unlines [let fs = flowgram (trim x) in printf "%6.1f " (dist a $ map fromIntegral fs) ++ concatMap (printf "%3d ") fs | x <- v] filter_unique :: Options -> FingerPrints -> [ReadBlock] -> ResultList filter_unique opts dups = go M.empty @@ -133,6 +132,7 @@ dist :: [Double] -> [Double] -> Double dist = (.) sum . zipWith (\x y -> (6*(x-y)/(x+y+4))^(2::Int)) +averageflow :: [ReadBlock] -> [Double] averageflow = go . map (map fromIntegral . flowgram . trim) where go [] = [] go xs = avg (map head xs) : go (filter (not . null) $ map tail xs)
+ src/Fork.hs view
@@ -0,0 +1,14 @@+module Fork where++import Control.Concurrent+import Control.Exception++-- | Spawn a set of threads and wait for them to complete.+forkAndWait :: [IO ()] -> IO ()+forkAndWait actions = mapM myForkIO actions >>= mapM_ takeMVar+ where+ myForkIO :: IO () -> IO (MVar ())+ myForkIO io = do+ mvar <- newEmptyMVar+ _ <- forkIO (io `finally` putMVar mvar ())+ return mvar
src/Metrics.hs view
@@ -9,7 +9,7 @@ -- | Take the fractional parts of the flows, and sum their squares (the "K²" metric) quals :: [Flow] -> Flow-quals q = floor $ (100 - 2*(sqrt $ (/fromIntegral (length q)) $ sum $ map (fromIntegral . (^2) . (flip (-) 50) . (`mod` 100) . (+50)) $ q))+quals q = floor $ ((100::Double) - 2*(sqrt $ (/fromIntegral (length q)) $ sum $ map (fromIntegral . (^(2::Integer)) . (flip (-) 50) . (`mod` 100) . (+50)) $ q)) -- | Count number of n's in the sequence -- The algorithm for generating Ns is a bit opaque, and appears to depend on the magnitude
+ src/Options.hs view
@@ -0,0 +1,37 @@+{-# LANGUAGE DeriveDataTypeable #-}++module Options where++import System.Console.CmdArgs++data Opts = Opts + { trimKey :: Bool+ , trim :: Bool+ , summarize :: Maybe FilePath+ , info :: Maybe FilePath+ , fasta :: Maybe FilePath+ , fastq :: Maybe FilePath+ , flowgram :: Maybe FilePath+ , histogram :: Maybe FilePath+ , inputs :: FilePath+ , text :: Maybe FilePath+ } deriving (Data,Typeable, Show, Eq)++opts :: Opts+opts = Opts+ { trimKey = False &= help "Trim only the TCAG key sequence"+ , trim = False &= help "Trim quality using clipping information" &= name "t"+ , summarize = def &= help "Output per sequence summary information" &= typFile+ , info = def &= help "Output brief overview of the contents" &= typFile + , fasta = def &= help "Output FASTA-formatted sequences" &= typFile &= name "f"+ , fastq = def &= help "Output FastQ-formatted sequence and quality" &= typFile &= name "q"+ , flowgram = def &= help "Output flowgram information in tabular form" &= typFile &= name "F"+ , histogram = def &= help "Output histogram of flow values" &= typFile &= name "h"+ , text = def &= help "Output SFF information as text (default)" &= typFile &= name "T"+ , inputs = def &= args &= typFile+ } + &= summary "flower v0.6 - Extract information from SFF files" + &= program "flower"++getArgs :: IO Opts+getArgs = cmdArgs opts
src/Print.hs view
@@ -7,10 +7,10 @@ import Data.Binary.Builder import Data.Monoid import Data.ByteString.Char8 (ByteString, pack)-import qualified Data.ByteString.Lazy.Char8 as LB import Data.Array.Unboxed import Data.Char (ord) +char :: Char -> Builder char = singleton . fromIntegral . ord putInt :: Int -> Builder