diff --git a/flower.cabal b/flower.cabal --- a/flower.cabal +++ b/flower.cabal @@ -1,5 +1,5 @@ Name: flower -Version: 0.7 +Version: 0.7.2 License: GPL Cabal-Version: >= 1.6 Author: Ketil Malde @@ -8,20 +8,21 @@ Category: Bioinformatics Synopsis: Analyze 454 flowgrams (.SFF files) Description: flower - FLOWgram ExtractoR tools + . + This is now obsolete, please install the biosff library + (which includes the 'flower' exectuable) instead. . The flower executable reads files in SFF-format and produces various output, 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. - . - 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 flowselect executable extracts reads from SFF-files, generating a new + SFF-file with a subset of the reads based on various quality criteria. + . + The flowt program removes (artificial) duplicates from SFF files. It's currently a work + in progress, but included if you'd like to play with it. It's faster than other approaches + (e.g. CD-HIT), and ought to be more sensitive and specific, but this needs to be proven. HomePage: http://biohaskell.org/Applications/Flower -- Source-Repository: http//malde.org/~ketil/biohaskell/flower Build-Type: Simple @@ -44,18 +45,6 @@ Build-Depends: bio >= 0.4.9, random, base >= 3 && < 5 Extensions: ExistentialQuantification -Executable frecover - Main-Is: FRecover.hs - Hs-Source-Dirs: src - Build-Depends: bio >= 0.4.9, base >= 3 && < 5 - Ghc-Options: -Wall - -Executable frename - Main-Is: FRename.hs - Hs-Source-Dirs: src - Build-Depends: bio >= 0.4.9, base >= 3 && < 5, bytestring >= 0.9.1 - Ghc-Options: -Wall - Executable flowt Main-Is: Flowt.hs Hs-Source-Dirs: src @@ -63,5 +52,5 @@ Ghc-Options: -Wall -fno-warn-unused-do-bind else Ghc-Options: -Wall - Build-Depends: bio >= 0.4.9, base >= 3 && < 5, cmdargs >= 0.5, containers, bytestring >= 0.9.1 + Build-Depends: bio > 0.4.9, base >= 3 && < 5, cmdargs >= 0.5, containers, bytestring >= 0.9.1 Extensions: DeriveDataTypeable diff --git a/src/FRecover.hs b/src/FRecover.hs deleted file mode 100644 --- a/src/FRecover.hs +++ /dev/null @@ -1,9 +0,0 @@ - -module Main where -import Bio.Sequence.SFF -import System.Environment (getArgs) - -main :: IO () -main = mapM_ recoverFile =<< getArgs - where recoverFile f = writeSFF (f++"_recovered") =<< recoverSFF f - diff --git a/src/FRename.hs b/src/FRename.hs deleted file mode 100644 --- a/src/FRename.hs +++ /dev/null @@ -1,30 +0,0 @@ - -{-| - 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 :: IO () -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 }} diff --git a/src/Flowt.hs b/src/Flowt.hs --- a/src/Flowt.hs +++ b/src/Flowt.hs @@ -29,8 +29,10 @@ 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! (use trimKey?) trim :: ReadBlock -> ReadBlock -trim = id -- trimFromTo 4 10000 <- this trims to last base called position! (use trimKey?) +trim rb = let rh = read_header rb in trimFromTo (clip_qual_left rh) (clip_qual_right rh) rb data Options = O { thresh :: Double , fplen :: Int @@ -79,9 +81,9 @@ 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 + showcluster (k,v) = let a = map fromIntegral $ flowgram $ trim $ head v -- averageflow v in printf "%16x" k++":\t" ++ show (length v) ++ unwords (map (BC.unpack . read_name . read_header) v) ++ "\n" ++ - unlines [let fs = flowgram (trim x) in printf "%6.1f " (dist a $ map fromIntegral fs) ++ concatMap (printf "%3d ") fs | x <- v] + unlines [let fs = flowgram (trim x) in printf "%8.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 @@ -101,12 +103,6 @@ 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 @@ -129,9 +125,14 @@ matches :: ReadBlock -> ReadBlock -> Double matches old new = let f = map fromIntegral . flowgram . trim in dist (f old) (f new) +-- Flowgram distance, not including zeroes at end of flowgram dist :: [Double] -> [Double] -> Double -dist = (.) sum . zipWith (\x y -> (6*(x-y)/(x+y+4))^(2::Int)) - +dist as bs = go (reverse as) (reverse bs) + where go (x:xs) (y:ys) = + if x==0 || y==0 then go xs ys + else sum $ zipWith (\i j -> (6*(i-j)/(i+j+4))^(2::Int)) (x:xs) (y:ys) + go _ _ = 0 + averageflow :: [ReadBlock] -> [Double] averageflow = go . map (map fromIntegral . flowgram . trim) where go [] = [] diff --git a/src/Print.hs b/src/Print.hs --- a/src/Print.hs +++ b/src/Print.hs @@ -18,7 +18,7 @@ else putInt' y where putInt' x = case x `divMod` 1000 of (0,r) -> fromByteString (ints!r) - (a,r) -> putInt a `append` putInt3 r + (a,r) -> putInt a `append` putInt3_0 r -- zero padded ints @@ -30,10 +30,15 @@ putInt3 x | x<1000 = fromByteString (int3s!x) | otherwise = fromByteString (pack "xxx") -ints, int2s, int3s :: Array Int ByteString +putInt3_0 :: Int -> Builder +putInt3_0 x | x<1000 = fromByteString (int3s0!x) + | otherwise = fromByteString (pack "xxx") + +ints, int2s, int3s, int3s0 :: Array Int ByteString ints = listArray (0,999) [pack (show i) | i <- [0..999::Int]] int2s = listArray (0,99) (map (pack . ('0':) . show) [0..9::Int]++[pack (show i) | i <- [10..99::Int]]) int3s = listArray (0,999) (map (pack . (' ':) . (' ':) . show) [0..9::Int] ++ map (pack . (' ':) . show ) [10..99::Int] ++ map (pack . show) [100..999::Int]) +int3s0 = listArray (0,999) (map (pack . ('0':) . ('0':) . show) [0..9::Int] ++ map (pack . ('0':) . show ) [10..99::Int] ++ map (pack . show) [100..999::Int]) putDate :: Int -> Int -> Int -> Builder putDate y m d = mconcat [putInt y, dash, putInt2 m, dash, putInt2 d]