packages feed

flower 0.1.2 → 0.2

raw patch · 4 files changed

+102/−88 lines, 4 filesdep ~bio

Dependency ranges changed: bio

Files

flower.cabal view
@@ -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 <http://malde.org/~ketil/biohaskell/flower>.  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:
src/Flower.hs view
@@ -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] <file.sff> [<file2.sff> ..]\n"-                ++"  -r  output reads in Fasta format\n"-                ++"  -R  output reads in Fasta format with associated .qual\n"-                ++"      (generates files instead of writing to <stdout>)\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] <file.sff> [<file2.sff> ..]\n"+                             ++"  -r  output reads in Fasta format\n"+                             ++"  -R  output reads in Fasta format with associated .qual\n"+                             ++"      (generates files instead of writing to <stdout>)\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
+ src/Metrics.hs view
@@ -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+
− src/Statistics.hs 
@@ -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]