flower 0.2 → 0.3
raw patch · 8 files changed
+364/−27 lines, 8 filesdep +cmdargsdep +containersdep +randomdep ~basedep ~binarydep ~bionew-component:exe:flowselectnew-component:exe:flowtnew-component:exe:frecovernew-component:exe:frename
Dependencies added: cmdargs, containers, random
Dependency ranges changed: base, binary, bio
Files
- flower.cabal +35/−7
- src/FRecover.hs +9/−0
- src/FRename.hs +31/−0
- src/FlowSelect.hs +100/−0
- src/Flower.hs +33/−12
- src/Flowt.hs +139/−0
- src/Metrics.hs +10/−3
- src/Print.hs +7/−5
flower.cabal view
@@ -1,27 +1,55 @@ Name: flower-Version: 0.2+Version: 0.3 License: GPL-+Cabal-Version: >= 1.6 Author: Ketil Malde Maintainer: Ketil Malde <ketil@malde.org> Category: Bioinformatics Synopsis: Analyze 454 flowgrams (.SFF files)-Description: flower - The FLOWgram ExtracteR+Description: flower - FLOWgram ExtractoR tools .- Reads files in SFF-format and produces various output, including sequences- with quality, or flowgram data in tabular format.+ The flower executable reads files in SFF-format and produces various output, + including sequences with quality, or flowgram data in tabular format. .+ The fselect executable extracts reads from SFF-files, generating a new+ SFF-file with a subset of the reads based on various criteria.+ .+ 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.+ . The Darcs repository is at <http://malde.org/~ketil/biohaskell/flower>. HomePage: http://malde.org/~ketil/biohaskell/flower-Build-Depends: bio >= 0.4, base >=3 && <4, array >= 0.1, bytestring >= 0.9.1, binary+Build-Depends: bio >= 0.4.2, base >=3 && <5, array >= 0.1, bytestring >= 0.9.1, binary == 0.4.*, random, cmdargs, containers Build-Type: Simple Tested-with: GHC==6.8.3 - -- Data-files: README+-- Data-files: README+ Executable: flower Main-Is: Flower.hs Other-Modules: Print, Metrics Hs-Source-Dirs: src Ghc-Options: ++Executable: flowselect+Main-Is: FlowSelect.hs+Other-Modules: Metrics+Hs-Source-Dirs: src+Extensions: ExistentialQuantification++Executable: frecover+Main-Is: FRecover.hs+Hs-Source-Dirs: src++Executable: frename+Main-Is: FRename.hs+Hs-Source-Dirs: src++Executable: flowt+Main-Is: Flowt.hs+Hs-Source-Dirs: src+Extensions: DeriveDataTypeable
+ src/FRecover.hs view
@@ -0,0 +1,9 @@++module Main where+import Bio.Sequence.SFF+import System.Environment (getArgs)++main = mapM_ recoverFile =<< getArgs++recoverFile f = writeSFF (f++"_recovered") =<< recoverSFF f+
+ src/FRename.hs view
@@ -0,0 +1,31 @@++{-|+ Rename reads in .SFF files to avoid name clashes.+ Apparently, reads with the same name crashes Newbler, and is+ in any case a bad idea. This ensures uniqueness by appending a serial number to each read name in a set of files.+-}++module Main where++import Bio.Sequence.SFF+import System.Environment (getArgs)+import qualified Data.ByteString.Char8 as B++main = do+ fs <- getArgs+ if null fs then putStrLn "Usage: frename file1.sff [file2.sff ...]"+ else renameSFFs fs+ +renameSFFs :: [FilePath] -> IO ()+renameSFFs = go 0 + where go _ [] = return ()+ go current (f:fs) = do+ (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 }}
+ src/FlowSelect.hs view
@@ -0,0 +1,100 @@+{-# LANGUAGE ExistentialQuantification #-}++module Main where++import Prelude hiding (LT,GT)++import Bio.Sequence.SFF+import System.Environment (getArgs)+import System.Random++import Metrics++main :: IO ()+main = do+ args <- getArgs+ (input,output,myfilter) <- parseArgs args+ SFF h rs <- readSFF input+ c <- writeSFF' output (SFF h $ myfilter rs)+ putStrLn ("Wrote "++show c++" reads.")++type FilterSFF = [ReadBlock] -> [ReadBlock]++parseArgs :: [String] -> IO (FilePath,FilePath,FilterSFF)+parseArgs [e,i] = do+ f <- case words e of + ["Rand",x] -> do+ ps <- randomRs (0,1) `fmap` newStdGen+ let t = read x :: Double+ return (map snd . filter ((<t).fst) . zip ps)+ ["Pick",n] -> undefined+ _ -> return (filter (apply (read e) . getChars))+ return (i,"selected.sff", f)++parseArgs _ = error "Usage: fselect <expression> <input.sff>"++-- | This structure represents selection parameters for one read+data Characteristics = Ch { k2, ee :: Double -- k-square, expected errors+ , ns, len, tlen :: Int -- lenght, trimmed length+ }++getChars :: ReadBlock -> Characteristics+getChars rb = let rh = read_header rb + in Ch { k2 = (/100) $ fromIntegral $ quals $ flowgram rb+ , ee = 0 -- fixme!+ , ns = n_count rb+ , len = fromIntegral $ num_bases rh+ , tlen = fromIntegral $ clip_qual_right rh - clip_qual_left rh + 1}++data FilterFunction = forall a . Ord a => LT (Characteristics -> a) a+ | forall a . Ord a => GT (Characteristics -> a) a++instance Show FilterFunction where show _ = "<filterfunction>"++instance Read FilterFunction where+ readsPrec i 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)]+ (c:"tlen":rest) -> [((lookupO c) tlen x,r) | (x,r) <- (reads $ unwords rest)]+ (c:"ncount":rest) -> [((lookupO c) ns x,r) | (x,r) <- (reads $ unwords rest)]+ _ -> error ("Couldn't parse FilterFunction: "++take 100 str)++lookupO "LT" = LT+lookupO "GT" = GT+lookupO x = error ("FilterFunction must be either LT or GT, was "++take 100 x)+++data Filter = Func FilterFunction + | And Filter Filter+ | Or Filter Filter+ | Not Filter + deriving Show++-- Okay, so we should really return/expect all parses here.+instance Read Filter where+ readsPrec i str = readParen False p str+ where p s = case words s of + "And":rest -> let [(a,r)] = reads (unwords rest)+ [(b,c)] = reads r+ in [(And a b,c)]+ "Or":rest -> let ((a,r):_) = reads (unwords rest)+ ((b,c):_) = reads r+ in [(Or a b,c)]+ "Not":rest -> let ((a,r):_) = reads (unwords rest)+ in [(Not a,r)]+ "Func":rest -> let ((a,r):_) = reads (unwords rest)+ in [(Func a,r)]+ _ -> [] -- error ("Couldn't parse "++take 100 s)++-- myFilter = LT k2 1.5++eval :: FilterFunction -> Characteristics -> Bool+eval (LT f a) c = f c < a+eval (GT f a) c = f c > a++apply :: Filter -> Characteristics -> Bool+apply (Func f) = eval f +apply (And f1 f2) = \r -> apply f1 r && apply f2 r+apply (Or f1 f2) = \r -> apply f1 r || apply f2 r+apply (Not f) = \r -> not (apply f r)
src/Flower.hs view
@@ -3,8 +3,10 @@ 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 Print @@ -28,10 +30,17 @@ main = do args <- getArgs let (opts,files) = partition (\p -> case p of ('-':_) -> True; _ -> False) args++ 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 opts of - ["-r"] -> hWriteFasta stdout . sffToSequence- ["-R"] -> writeFastaQual ("flower.fasta") ("flower.qual") . sffToSequence+ 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@@ -50,8 +59,12 @@ ++" -i output header information\n" ++" -s output a summary of each read" )- writer `seq` mapM_ (\f -> writer =<< readSFF f) files+ writer `seq` mapM_ (\f -> writer =<< reader f) files +-- 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+ -- ------------------------------------------------------------ -- The -i option: Print header info -- ------------------------------------------------------------@@ -69,23 +82,30 @@ -- | 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\tqual\ttrimqual"+ 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 -- todo: date and time are usually constants! sum1 :: ReadBlock -> Builder sum1 r = let rh = read_header r nb = num_bases rh- h = read_name rh + h = read_name rh+ tr = trim r (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+ 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 $ take (fromIntegral (qright-qleft)) $ drop (fromIntegral qleft) $ flowgram r), nl])+ ++ 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'@@ -102,7 +122,8 @@ fi :: Flow -> ByteString fi f | f <= 9999 && f >= 0 = farray!f- | otherwise = error ("Can't show a flow value of "++show f)+ | otherwise = let (i,r) = f `divMod` 100 in B.pack (show i++"."++show r) + -- error ("Can't show flow values outside [0..99.99] (You had: "++show f++")") farray :: Array Flow ByteString farray = listArray (0,9999) [B.pack (showFFloat (Just 2) i "") | i <- [0,0.01..99.99::Double]]@@ -168,7 +189,7 @@ ins1 ('T',i) = bump t i ins1 (x,_) = error ("Illegal character "++show x++" in flow!") bump ar i = readArray ar i >>= \x -> writeArray ar i (x+1)- mapM_ ins1 (zip (cycle fl) (concat scores))+ mapM_ ins1 (zip (cycle fl) (map (\x->if x>9999 || x<0 then 9999 else x) $ concat scores)) a' <- unsafeFreeze a c' <- unsafeFreeze c g' <- unsafeFreeze g@@ -177,5 +198,5 @@ showHist :: (Hist,Hist,Hist,Hist) -> String showHist (as,cs,gs,ts) = "Score\tA\tC\tG\tT\tsum\n" ++ - unlines [concat $ intersperse "\t" $ (showFFloat (Just 2) (fromIntegral sc/100::Double) "") : map show [as!sc,cs!sc,gs!sc,ts!sc, as!sc+cs!sc+gs!sc+ts!sc]+ unlines [concat $ intersperse "\t" $ showFFloat (Just 2) (fromIntegral sc/100::Double) "" : map show [as!sc,cs!sc,gs!sc,ts!sc, as!sc+cs!sc+gs!sc+ts!sc] | sc <- [0..9999]]
+ src/Flowt.hs view
@@ -0,0 +1,139 @@+{-# LANGUAGE DeriveDataTypeable #-}++module Main where++import Bio.Sequence.SFF hiding (trim)+import Bio.Util (countIO)+import qualified Data.ByteString as B+import qualified Data.ByteString.Char8 as BC++import qualified Data.IntMap as M+import Data.IntMap (IntMap)++-- import Bloom+import qualified Data.IntSet as S+import Data.IntSet (IntSet)+import Text.Printf++import System.IO+import Control.Monad (when)++import System.Console.CmdArgs++version :: String+version = "flowt v0, copyright 2009-2010 Ketil Malde"++type FingerPrints = IntSet+type DupMap = IntMap [ReadBlock]+ +data ResultList = Then ReadBlock ResultList | EndWith DupMap++splitRes :: ResultList -> ([ReadBlock], DupMap)+splitRes (Then x rs) = let (ys,e) = splitRes rs in (x:ys,e)+splitRes (EndWith e) = ([],e)++trim = id -- trimFromTo 4 10000 <- this trims to last base called position!++data Options = O { thresh :: Double+ , fplen :: Int+ , summarize :: FilePath+ , clusters :: Bool+ , 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."++vlog :: Bool -> String -> IO ()+vlog v s = when v (hPutStr stderr s >> hFlush stderr)++main :: IO ()+main = do+ opts <- cmdArgs version [modes]+ -- putStrLn $ show opts+ verb <- isLoud+ vlog verb "Building the fingerprint index"+ let sff = case input opts of [x] -> x; _ -> error "You need to specify a (single) file name!\n(or use 'flowt -h' for help)."+ dups <- mkbf opts sff+ vlog verb (seq dups "...done!")+ SFF hs rs <- readSFF sff+ let (uqs,ds) = splitRes $ filter_unique opts dups rs++ writeSFF' "unique.sff" =<< SFF hs `fmap`+ (if verb then countIO "Output unique: " "..done!" 100 else return) uqs++ when (not . null . summarize $ opts) (gensum (summarize opts) ds)++ let cg = concatMap (clusterGroups opts) (M.elems ds)+ writeSFF' "duplic.sff" =<< SFF hs `fmap` + (if verb then countIO "Output cluster reps: " "..done!" 100 else return)+ (map head cg)+ when (clusters opts) $ BC.writeFile "clusters.txt" $ + BC.unlines $ map (BC.unwords . map (read_name . read_header)) $ cg+ return ()+ +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]+ +filter_unique :: Options -> FingerPrints -> [ReadBlock] -> ResultList+filter_unique opts dups = go M.empty + where go dm (r:rs) = let fp = fingerprint (fplen opts) r + in if fp `S.member` dups then let dm' = myinsert fp r dm in dm' `seq` go dm' rs+ else r `Then` go dm rs+ go dm [] = EndWith dm+ myinsert fp r dm = let v = M.findWithDefault [] fp dm+ v' = r:v+ in v `seq` v' `seq` dm `seq` M.insert fp v' dm++mkbf :: Options -> FilePath -> IO FingerPrints+mkbf opts sff = do + SFF _ rs <- readSFF sff+ let go seen dup (fp:rest) = if fp `S.member` seen then go seen (S.insert fp dup) rest+ else go (S.insert fp seen) dup rest+ go _ dup [] = dup+ return $ go S.empty S.empty $ map (fingerprint $ fplen opts) rs++{-+mkbf sff = do + SFF _ rs <- readSFF sff+ return $ snd $ mkFilters 1000000 $ map fingerprint rs+-}++-- | Calculating a fingerprint - basically just a hash of the first 20 elements of the flow_index. On 64 bits, it is +-- possible to use more, this is a sensitivity/specificity tradeoff.+fingerprint :: Int -> ReadBlock -> Int+fingerprint fpl = foldr (\x y -> y*3+x-1) 0 . map fromIntegral . B.unpack . B.take fpl . B.filter (/=0) . flow_index . trim+ -- trim == drop 4 for the TCAG key, then 3^20 ~ 2^32 for range (or should that be 4^16?)++-- | Align each cluster and merge reads that appear to be from the same clone.+cluster :: Options -> DupMap -> [ReadBlock]+cluster opts = concatMap (map mergeCG . clusterGroups opts) . M.elems+ +mergeCG :: [ReadBlock] -> ReadBlock+mergeCG = head -- todo: build a consensus flowgram and base call it.++clusterGroups :: Options -> [ReadBlock] -> [[ReadBlock]]+clusterGroups opts (c:cs) = let (this,rest) = span ((<= thresh opts) . matches c) cs + in (c:this) : clusterGroups opts rest+clusterGroups _ [] = []++-- crude flowgram match check, valid range?+matches :: ReadBlock -> ReadBlock -> Double+matches old new = let f = map fromIntegral . flowgram . trim in dist (f old) (f new)++dist :: [Double] -> [Double] -> Double+dist = (.) sum . zipWith (\x y -> (6*(x-y)/(x+y+4))^(2::Int))++averageflow = go . map (map fromIntegral . flowgram . trim)+ where go [] = []+ go xs = avg (map head xs) : go (filter (not . null) $ map tail xs)+ avg xs = sum xs / fromIntegral (length xs)
src/Metrics.hs view
@@ -3,6 +3,7 @@ module Metrics where import Bio.Sequence.SFF+import Bio.Sequence.SeqData -- import Test.QuickCheck @@ -10,6 +11,12 @@ quals :: [Flow] -> Flow quals q = floor $ (100 - 2*(sqrt $ (/fromIntegral (length q)) $ sum $ map (fromIntegral . (^2) . (flip (-) 50) . (`mod` 100) . (+50)) $ q)) -prop_quals :: [Flow] -> Bool-prop_quals fs = let q = quals fs in q <= 100 && q >= 0-+-- | Count number of n's in the sequence+-- The algorithm for generating Ns is a bit opaque, and appears to depend on the magnitude +-- of the noise flow values. We chicken out, and just count the called sequence.+n_count :: ReadBlock -> Int+n_count r = length . filter isN . clip . toStr . bases $ r+ where isN x = x=='N' || x == 'n'+ clip = take (right-left+1) . drop left+ right = fromIntegral $ clip_qual_right (read_header r)+ left = fromIntegral $ clip_qual_left (read_header r)
src/Print.hs view
@@ -1,7 +1,7 @@ module Print ( Builder, toLazyByteString, mconcat, fromByteString, char- , putInt, putInt2, putDate, putTime+ , putInt, putInt2, putDate, putTime, putFix ) where import Data.Binary.Builder@@ -11,8 +11,6 @@ import Data.Array.Unboxed import Data.Char (ord) --- import Test.QuickCheck- char = singleton . fromIntegral . ord putInt :: Int -> Builder@@ -43,5 +41,9 @@ putTime h m s = mconcat [putInt2 h, col, putInt2 m, col, putInt2 s] where col = char ':' --prop_int i = toLazyByteString (putInt i) == LB.pack (show i)+-- a bit hackish, maybe?+putFix :: Int -> Int -> Builder+putFix 1 n = let (i,f) = n `divMod` 10 in putInt i `mappend` char '.' `mappend` putInt f+putFix 2 n = let (i,f) = n `divMod` 100 in putInt i `mappend` char '.' `mappend` putInt2 f+putFix 3 n = let (i,f) = n `divMod` 1000 in putInt i `mappend` char '.' `mappend` putInt3 f+putFix _ _ = error "putFix only supports up to three fractional decimals"