diff --git a/flower.cabal b/flower.cabal --- a/flower.cabal +++ b/flower.cabal @@ -1,27 +1,55 @@ Name: flower -Version: 0.2 +Version: 0.3 License: GPL - +Cabal-Version: >= 1.6 Author: Ketil Malde Maintainer: Ketil Malde 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 . 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 diff --git a/src/FRecover.hs b/src/FRecover.hs new file mode 100644 --- /dev/null +++ b/src/FRecover.hs @@ -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 + diff --git a/src/FRename.hs b/src/FRename.hs new file mode 100644 --- /dev/null +++ b/src/FRename.hs @@ -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 }} diff --git a/src/FlowSelect.hs b/src/FlowSelect.hs new file mode 100644 --- /dev/null +++ b/src/FlowSelect.hs @@ -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 (( undefined + _ -> return (filter (apply (read e) . getChars)) + return (i,"selected.sff", f) + +parseArgs _ = error "Usage: fselect " + +-- | 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 _ = "" + +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) diff --git a/src/Flower.hs b/src/Flower.hs --- a/src/Flower.hs +++ b/src/Flower.hs @@ -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]] diff --git a/src/Flowt.hs b/src/Flowt.hs new file mode 100644 --- /dev/null +++ b/src/Flowt.hs @@ -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) diff --git a/src/Metrics.hs b/src/Metrics.hs --- a/src/Metrics.hs +++ b/src/Metrics.hs @@ -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) diff --git a/src/Print.hs b/src/Print.hs --- a/src/Print.hs +++ b/src/Print.hs @@ -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"