diff --git a/flower.cabal b/flower.cabal --- a/flower.cabal +++ b/flower.cabal @@ -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 . -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 diff --git a/src/FRecover.hs b/src/FRecover.hs --- a/src/FRecover.hs +++ b/src/FRecover.hs @@ -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 diff --git a/src/FRename.hs b/src/FRename.hs --- a/src/FRename.hs +++ b/src/FRename.hs @@ -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 }} diff --git a/src/FlowSelect.hs b/src/FlowSelect.hs --- a/src/FlowSelect.hs +++ b/src/FlowSelect.hs @@ -27,7 +27,7 @@ ps <- randomRs (0,1) `fmap` newStdGen let t = read x :: Double return (map snd . filter (( undefined + ["Pick",_] -> undefined _ -> return (filter (apply (read e) . getChars)) return (i,"selected.sff", f) @@ -52,7 +52,7 @@ instance Show FilterFunction where show _ = "" 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 diff --git a/src/Flower.hs b/src/Flower.hs --- a/src/Flower.hs +++ b/src/Flower.hs @@ -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] [ ..]\n" - ++" -r output reads in Fasta format\n" - ++" -R output reads in Fasta format with associated .qual\n" - ++" (generates files instead of writing to )\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 diff --git a/src/Flowt.hs b/src/Flowt.hs --- a/src/Flowt.hs +++ b/src/Flowt.hs @@ -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) diff --git a/src/Fork.hs b/src/Fork.hs new file mode 100644 --- /dev/null +++ b/src/Fork.hs @@ -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 diff --git a/src/Metrics.hs b/src/Metrics.hs --- a/src/Metrics.hs +++ b/src/Metrics.hs @@ -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 diff --git a/src/Options.hs b/src/Options.hs new file mode 100644 --- /dev/null +++ b/src/Options.hs @@ -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 diff --git a/src/Print.hs b/src/Print.hs --- a/src/Print.hs +++ b/src/Print.hs @@ -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