diff --git a/flower.cabal b/flower.cabal --- a/flower.cabal +++ b/flower.cabal @@ -1,5 +1,5 @@ Name: flower -Version: 0.1.2 +Version: 0.2 License: GPL Author: Ketil Malde @@ -15,13 +15,13 @@ The Darcs repository is at . HomePage: http://malde.org/~ketil/biohaskell/flower -Build-Depends: bio >= 0.3.5 && <0.4, base >= 3 && < 4, array >= 0.1, bytestring >= 0.9.1, binary +Build-Depends: bio >= 0.4, base >=3 && <4, array >= 0.1, bytestring >= 0.9.1, binary Build-Type: Simple Tested-with: GHC==6.8.3 -- Data-files: README Executable: flower Main-Is: Flower.hs -Other-Modules: Print, Statistics +Other-Modules: Print, Metrics Hs-Source-Dirs: src Ghc-Options: diff --git a/src/Flower.hs b/src/Flower.hs --- a/src/Flower.hs +++ b/src/Flower.hs @@ -11,6 +11,7 @@ import System.IO (stdout) import System.Environment (getArgs) import Numeric (showFFloat) +import Data.Char (toLower) import Data.List (intersperse, partition) import Data.ByteString.Char8 (pack,unpack,ByteString) import qualified Data.ByteString.Char8 as B @@ -21,33 +22,54 @@ import Data.Array.ST import Control.Monad.ST +import Metrics + main :: IO () main = do args <- getArgs let (opts,files) = partition (\p -> case p of ('-':_) -> True; _ -> False) args - case opts of - ["-r"] -> mapM_ (\f -> hWriteFasta stdout . sffToSequence =<< readSFF f) files - ["-R"] -> mapM_ (\f -> writeFastaQual (f++".fasta") (f++".qual") . sffToSequence =<< readSFF f) files - ["-q"] -> mapM_ (\f -> hWriteFastQ stdout . sffToSequence =<< readSFF f) files - ["-f"] -> L1.putStrLn . L1.fromChunks . intersperse (B.pack "\n") . concat =<< mapM showflow files - ["-h"] -> mapM_ (\f -> (putStr . sffToHistogram) =<< readSFF f) files - ["-s"] -> mapM_ (\f -> summarize =<< readSFF f) files + writer :: SFF -> IO () + writer = case opts of + ["-r"] -> hWriteFasta stdout . sffToSequence + ["-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 - _ -> error ("Usage: flower -[f|q|r|R] [ ..]\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" - ++" -s output a summary of each read" - ) + _ -> 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 =<< readSFF f) files +-- ------------------------------------------------------------ +-- 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) + ] + -- ---------------------------------------------------------- +-- 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\tx_loc\ty_loc\tlen\tqual" + putStrLn "# name........\tdate......\ttime....\treg\ttrim_l\ttrim_r\tx_loc\ty_loc\tlen\tqual\ttrimqual" L1.putStrLn . toLazyByteString . mconcat . map sum1 $ rs -- todo: date and time are usually constants! @@ -55,44 +77,52 @@ sum1 r = let rh = read_header r nb = num_bases rh h = read_name rh - rn = decodeReadName h - (y,m,d) = date rn - reg = region rn - (hh,mm,ss) = time rn - in mconcat ([fromByteString h, tb, putDate y m d, tb, putTime hh mm ss, tb, putInt2 reg - ,tb, putInt (fromIntegral $ x_loc rn), tb, putInt (fromIntegral $ y_loc rn), tb, putInt (fromIntegral nb) - ,tb, fromByteString (fi $ quals $ flowgram r), nl]) - --- | Take the fractional parts of the flows, and sum their squares -quals :: [Flow] -> Flow -quals q = floor $ (*(100/fromIntegral (length q))) $ sqrt $ sum $ map (fromIntegral . (^2) . (flip (-) 50) . (`mod` 100) . (+50)) $ q + (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) + 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]) -tb, nl :: Builder +tb, nl, q :: Builder tb = char '\t' nl = char '\n' +q = char '?' --- these are clumsy, since we just might need the file name -showflow :: FilePath -> IO [ByteString] -showflow f = return . {- map (\s -> B.concat [B.pack f,t,s]) . -} showrun =<< readSFF f +-- ---------------------------------------------------------- +-- 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 = (!) farray +fi f | f <= 9999 && f >= 0 = farray!f + | otherwise = error ("Can't show a flow value of "++show f) farray :: Array Flow ByteString -farray = listArray (0,10000) [B.pack (showFFloat (Just 2) i "") | i <- [0,0.01..99.99::Double]] - -showrun :: SFF -> [ByteString] -showrun (SFF h rs) = concatMap (showread h) rs +farray = listArray (0,9999) [B.pack (showFFloat (Just 2) i "") | i <- [0,0.01..99.99::Double]] tab :: ByteString tab = B.pack "\t" showread :: CommonHeader -> ReadBlock -> [ByteString] -showread h rd = let rn = read_name $ read_header rd +showread h rd = let rh = read_header rd + rn = read_name rh + maskFlows = mask rh 1 qgroups . unpack qgroups = qgroup (B1.unpack $ flow_index rd) (L1.unpack $ quality rd) - format p c v q = B.concat [rn,tab,B.pack (show p),tab,B.pack [c],tab,fi v,tab,B.pack (show q)] - in zipWith4 format [(1::Int)..] (unpack $ flow h) (flowgram rd) qgroups + format p c v q = B.concat [rn,tab,B.pack (show p),tab,B.pack [c],tab,fi v,tab,B.pack (init $ drop 1 $ show q)] + 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 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 + 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 zipWith4 _ _ _ _ _ = [] @@ -107,8 +137,21 @@ in q1 : qgroup irest qrest qgroup (i:is) qs = [] : qgroup (i-1:is) qs -sffToHistogram :: SFF -> String + +-- ---------------------------------------------------------- +-- 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/Metrics.hs b/src/Metrics.hs new file mode 100644 --- /dev/null +++ b/src/Metrics.hs @@ -0,0 +1,15 @@ +-- Calculate various characteristics on sequence quality + +module Metrics where + +import Bio.Sequence.SFF + +-- import Test.QuickCheck + +-- | 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)) + +prop_quals :: [Flow] -> Bool +prop_quals fs = let q = quals fs in q <= 100 && q >= 0 + diff --git a/src/Statistics.hs b/src/Statistics.hs deleted file mode 100644 --- a/src/Statistics.hs +++ /dev/null @@ -1,44 +0,0 @@ --- | Yet another simple module for implementing statistics stuff. - -module Statistics (normal, normals, stdnormal, stdnormals, module System.Random) where -import System.Random - -stdnormal :: StdGen -> (Double,StdGen) -stdnormal = normal 0 1 - -stdnormals :: StdGen -> [Double] -stdnormals = normals 0 1 - -normal :: Double -> Double -> StdGen -> (Double,StdGen) -normal mu sigma = \g -> let (x,g') = randomR (0,1) g in (invcumnorm mu sigma x,g') - -normals :: Double -> Double -> StdGen -> [Double] -normals mu sigma = \g -> let (x,g') = normal mu sigma g in x : normals mu sigma g' - -lognormal :: Double -> Double -> StdGen -> Double -lognormal mu sigma = undefined - --- support - -invcumnorm mu sigma z = mu + search (-limit*sigma) (limit*sigma) - where search a b = let c = (a+b)/2 - cn = cumnorm 0 sigma c - in if abs (z - cn) < 10*epsilon || abs (a-b) < epsilon then c - else if cn > z then search a c - else search c b - -cumstdnorm :: Double -> Double -cumstdnorm x = 0.5*(1+erf (x/sqrt 2)) - -cumnorm :: Double -> Double -> Double -> Double -cumnorm mu sigma x = 0.5*(1+erf((x-mu)/(sigma*sqrt 2))) - --- taylor expansion, see wikipedia "error function". Tested within the range (-limit..limit) -erf :: Double -> Double -erf x | x>limit = 1 - | x< negate limit = 0 - | otherwise = (2/sqrt pi)*sum (reverse $ takeWhile ((>=epsilon).abs) [ ((-1)**n*x**(2*n+1)) / (fac n*(2*n+1)) | n <- [0..]]) - -epsilon = 0.0000000001 -limit = 4.4 :: Double -fac x = product [2..x]