biosff (empty) → 0.1
raw patch · 10 files changed
+1109/−0 lines, 10 filesdep +arraydep +basedep +binarysetup-changed
Dependencies added: array, base, binary, biocore, bytestring, cmdargs, mtl
Files
- Setup.hs +3/−0
- biosff.cabal +32/−0
- src/Bio/Sequence/SFF.hs +459/−0
- src/Bio/Sequence/SFF_filters.hs +144/−0
- src/Bio/Sequence/SFF_name.hs +86/−0
- src/Flower/Fork.hs +14/−0
- src/Flower/Main.hs +245/−0
- src/Flower/Metrics.hs +23/−0
- src/Flower/Options.hs +52/−0
- src/Flower/Print.hs +51/−0
+ Setup.hs view
@@ -0,0 +1,3 @@+#!/usr/bin/env runhaskell+import Distribution.Simple+main = defaultMain
+ biosff.cabal view
@@ -0,0 +1,32 @@+Name: biosff+Version: 0.1+Synopsis: Library and executables for working with SFF files+Description: The library contains the functionality for reading and writing+ SFF files (sequencing data from 454 and Ion Torrent). It duplicates+ code from (and is incompatible with) the "bio" library.+Homepage: http://biohaskell.org/+License: LGPL+Author: Ketil Malde+Maintainer: ketil@malde.org+Stability: Experimental+Category: Bioinformatics+Build-type: Simple+Cabal-version: >=1.2++Library+ Exposed-modules: Bio.Sequence.SFF+ Other-modules: Bio.Sequence.SFF_name, Bio.Sequence.SFF_filters+ Build-depends: base >= 3 && < 5, biocore >= 0.1, binary < 0.5, bytestring, array+ Hs-Source-Dirs: src+ Ghc-Options: -Wall++Executable flower+ Main-Is: Main.hs+ Other-Modules: Fork, Options, Metrics, Print+ Build-Depends: base >= 3 && < 5, cmdargs, mtl >= 2+ Hs-Source-Dirs: src, src/Flower+ Ghc-Options: -Wall++-- Executable flowt++-- Executable flowselect
+ src/Bio/Sequence/SFF.hs view
@@ -0,0 +1,459 @@+{- | Read and write the SFF file format used by+ Roche\/454 sequencing to store flowgram data.++ A flowgram is a series of values (intensities) representing homopolymer runs of+ A,G,C, and T in a fixed cycle, and usually displayed as a histogram.++ This file is based on information in the Roche FLX manual. Among other sources for information about+ the format, are The Staden Package, which contains an io_lib with a C routine for parsing this format.+ According to comments in the sources, the io_lib implementation is based on a file+ called getsff.c, which I've been unable to track down. Other software parsing SFFs + are QIIME, sff_extract, and Celera's sffToCa.++ It is believed that all values are stored big endian.+-}++module Bio.Sequence.SFF ( SFF(..), CommonHeader(..)+ , ReadHeader(..), ReadBlock(..)+ , readSFF, writeSFF, writeSFF', recoverSFF+ -- , sffToSequence, rbToSequence+ , trim, trimFromTo -- , trimKey+ , baseToFlowPos, flowToBasePos+ , trimFlows+ , test, convert, flowgram+ , masked_bases, cumulative_index+ , packFlows, unpackFlows+ , Flow, Qual, Index, SeqData, QualData+ , ReadName (..), decodeReadName, encodeReadName+ ) where++import Bio.Core.Sequence+import Bio.Sequence.SFF_name++import Data.Int+import qualified Data.ByteString.Lazy as LB+import qualified Data.ByteString.Lazy.Char8 as LBC+import qualified Data.ByteString as B+import qualified Data.ByteString.Char8 as BC+import Data.ByteString (ByteString)+import Control.Monad (when,replicateM,replicateM_)++import Data.List (intersperse)+import Data.Binary+import Data.Binary.Get (getByteString,getLazyByteString)+import qualified Data.Binary.Get as G+import Data.Binary.Put (putByteString,putLazyByteString)+import Data.Char (toUpper, toLower)+import Text.Printf (printf)+import System.IO++-- | The type of flowgram value+type Flow = Int16+type Index = Word8++-- Global variables holding static information+-- | An SFF file always start with this magic number.+magic :: Int32+magic = 0x2e736666++-- | Version is always 1.+versions :: [Int32]+versions = [1]++-- | Read an SFF file.+readSFF :: FilePath -> IO SFF+readSFF f = return . decode =<< LB.readFile f++{-+-- | Extract the read without the initial (TCAG) key.+trimKey :: CommonHeader -> Sequence Nuc -> Maybe (Sequence Nuc)+trimKey ch (Seq n s q) = let (k,s2) = LB.splitAt (fromIntegral $ key_length ch) s+ in if LBC.map toLower k==LBC.map toLower (LB.fromChunks [key ch]) + then Just $ Seq n s2 (liftM (LB.drop (fromIntegral $ key_length ch)) q)+ else Nothing -- error ("Couldn't match key in sequence "++LBC.unpack n++" ("++LBC.unpack k++" vs. "++BC.unpack (key ch)++")!")+-}++instance BioSeq ReadBlock where+ seqlabel rb = SeqLabel $ LB.fromChunks [read_name $ read_header rb]+ seqdata rb = let h = read_header rb+ (left,right) = (clip_qual_left h, clip_qual_right h)+ (a,b) = LB.splitAt (fromIntegral right) $ unSD $ bases rb+ (c,d) = LB.splitAt (fromIntegral left-1) a+ in SeqData $ LBC.concat [LBC.map toLower c, LBC.map toUpper d,LBC.map toLower b]+ seqlength rb = fromIntegral $ num_bases $ read_header rb++instance BioSeqQual ReadBlock where+ seqqual = quality++{- -- | Extract the sequences from an 'SFF' data structure.+sffToSequence :: SFF -> [Sequence Nuc]+sffToSequence (SFF _ rs) = map rbToSequence rs++-- | Extract the sequence information from a 'ReadBlock'.+rbToSequence :: ReadBlock -> Sequence Nuc+rbToSequence r = Seq (LB.fromChunks [read_name h ,BC.pack (" qclip: "++show left ++".."++show right)])+ (seqdata r)+ (Just $ seqqual r)+ where h = read_header r+ (left,right) = (clip_qual_left h, clip_qual_right h)+-}++-- | Trim a 'ReadBlock' limiting the number of flows. If writing to+-- an SFF file, make sure you update the 'CommonHeader' accordingly.+-- See @examples/Flx.hs@ for how to use this. +trimFlows :: Integral i => i -> ReadBlock -> ReadBlock+trimFlows l rb = rb { read_header = rh { num_bases = fromIntegral n+ , clip_qual_right = min cqr $ fromIntegral n+ }+ , flow_data = B.take (2*fromIntegral l) (flow_data rb)+ , flow_index = B.take n (flow_index rb)+ , bases = SeqData $ LB.take (fromIntegral n) (unSD $ bases rb)+ , quality = QualData $ LB.take (fromIntegral n) (unQD $ quality rb)+ }+ where n = (flowToBasePos rb l)-1+ rh = read_header rb+ cqr = clip_qual_right rh++-- trimming the flowgram is necessary, but how to deal with the shift in flow+-- sequence - i.e. what to do when trimming "splits" a flow into trimmed/untrimmed bases?++-- | Trim a read to specific sequence position, inclusive bounds.+trimFromTo :: (Integral i) => i -> i -> ReadBlock -> ReadBlock+trimFromTo x r rd = let+ l = x-1+ trim_seq = LB.drop (fromIntegral l) . LB.take (fromIntegral r)+ trim_seq' = B.drop (fromIntegral l) . B.take (fromIntegral r)+ trim_flw = B.drop ((2*) $ fromIntegral $ baseToFlowPos rd l) . B.take ((2*) $ fromIntegral $ baseToFlowPos rd r)+ new_flw = trim_flw (flow_data rd)+ padding = B.replicate (B.length (flow_data rd) - B.length new_flw) 0+ rh = read_header rd+ [r',l'] = map fromIntegral [r,l]+ rh' = rh { num_bases = fromIntegral (r'-l')+ , clip_qual_left = max 0 $ clip_qual_left rh-l'+ , clip_qual_right = min (clip_qual_right rh-l') (r'-l'+1)+ }+ in rd { read_header = rh'+ , flow_data = B.concat [new_flw, padding]+ , flow_index = trim_seq' (flow_index rd)+ , bases = SeqData $ trim_seq $ unSD $ bases rd+ , quality = QualData $ trim_seq $ unQD $ quality rd+ }++-- | Trim a read according to clipping information+trim :: ReadBlock -> ReadBlock+trim rb = let rh = read_header rb in trimFromTo (clip_qual_left rh) (clip_qual_right rh) rb++-- | Convert a flow position to the corresponding sequence position+flowToBasePos :: Integral i => ReadBlock -> i -> Int+flowToBasePos rd fp = length $ takeWhile (<=fp) $ scanl (+) 0 $ map fromIntegral $ B.unpack $ flow_index rd++-- | Convert a sequence position to the corresponding flow position+baseToFlowPos :: Integral i => ReadBlock -> i -> Int+baseToFlowPos rd sp = sum $ map fromIntegral $ B.unpack $ B.take (fromIntegral sp) $ flow_index rd++-- | Read an SFF file, but be resilient against errors.+recoverSFF :: FilePath -> IO SFF+recoverSFF f = return . unRecovered . decode =<< LB.readFile f++-- | Write an 'SFF' to the specified file name+writeSFF :: FilePath -> SFF -> IO ()+writeSFF = encodeFile++-- | Write an 'SFF' to the specified file name, but go back and+-- update the read count. Useful if you want to output a lazy+-- stream of 'ReadBlock's. Returns the number of reads written.+writeSFF' :: FilePath -> SFF -> IO Int+writeSFF' f (SFF hs rs) = do+ h <- openFile f WriteMode+ LBC.hPut h $ encode hs+ c <- writeReads h (fromIntegral $ flow_length hs) rs+ hSeek h AbsoluteSeek 20+ LBC.hPut h $ encode c+ hClose h+ return $ fromIntegral c++-- | Write 'ReadBlock's to a file handle.+writeReads :: Handle -> Int -> [ReadBlock] -> IO Int32+writeReads _ _ [] = return 0+writeReads h i (r:rs) = do+ LBC.hPut h $ encode (RBI i r)+ c <- writeReads h i rs+ return $! (c+1)++data RBI = RBI Int ReadBlock++-- | Wrapper for ReadBlocks since they need additional information+instance Binary RBI where + put (RBI c r) = do+ putRB c r+ get = undefined+ +-- --------------------------------------------------+-- | test serialization by output'ing the header and first two reads +-- in an SFF, and the same after a decode + encode cycle.+test :: FilePath -> IO ()+test file = do + (SFF h rs) <- readSFF file + let sff = (SFF h (take 2 rs))+ putStrLn $ show $ sff+ putStrLn ""+ putStrLn $ show $ (decode $ encode sff :: SFF)++-- --------------------------------------------------+-- | Convert a file by decoding it and re-encoding it+-- This will lose the index (which isn't really necessary)+convert :: FilePath -> IO ()+convert file = writeSFF (file++".out") =<< readSFF file++-- | Generalized function for padding+pad :: Integral a => a -> Put+pad x = replicateM_ (fromIntegral x) (put zero) where zero = 0 :: Word8 ++-- | Generalized function to skip padding+skip :: Integral a => a -> Get ()+skip = G.skip . fromIntegral++-- | The data structure storing the contents of an SFF file (modulo the index)+data SFF = SFF !CommonHeader [ReadBlock]++instance Show SFF where + show (SFF h rs) = (show h ++ "Read Blocks:\n\n" ++ concatMap show rs)++instance Binary SFF where+ get = do+ -- Parse CommonHeader+ chead <- get+ -- Get the ReadBlocks+ rds <- replicateM (fromIntegral (num_reads chead))+ (do + rh <- get :: Get ReadHeader+ getRB chead rh+ )+ return (SFF chead rds)++ put (SFF hd rds) = do+ put hd+ mapM_ (put . RBI (fromIntegral $ flow_length hd)) rds++-- | Helper function for decoding a 'ReadBlock'.+{-# INLINE getRB #-}+getRB :: CommonHeader -> ReadHeader -> Get ReadBlock+getRB chead rh = do+ let nb = fromIntegral $ num_bases rh+ nb' = fromIntegral $ num_bases rh+ fl = fromIntegral $ flow_length chead+ fg <- getByteString (2*fl)+ fi <- getByteString nb+ bs <- getLazyByteString nb'+ qty <- getLazyByteString nb'+ let l = (fl*2+nb*3) `mod` 8+ when (l > 0) (skip (8-l))+ return (ReadBlock rh fg fi (SeqData bs) (QualData qty))++-- | A ReadBlock can't be an instance of Binary directly, since it depends on+-- information from the CommonHeader.+putRB :: Int -> ReadBlock -> Put+putRB fl rb = do+ put (read_header rb)+ putByteString (flow_data rb)+ -- ensure that flowgram has correct lenght+ replicateM_ (2*fl-B.length (flow_data rb)) (put (0::Word8))+ putByteString (flow_index rb)+ putLazyByteString (unSD $ bases rb)+ putLazyByteString (unQD $ quality rb)+ let nb = fromIntegral $ num_bases $ read_header rb+ l = (fl*2+nb*3) `mod` 8+ when (l > 0) (pad (8-l))++-- | Unpack the flow_data field into a list of flow values+unpackFlows :: ByteString -> [Flow]+unpackFlows = dec . map fromIntegral . B.unpack + where dec (d1:d2:rest) = d1*256+d2 : dec rest+ dec [] = []+ dec _ = error "odd flowgram length?!"++-- | Pack a list of flows into the corresponding binary structure (the flow_data field)+packFlows :: [Flow] -> ByteString+packFlows = B.pack . map fromIntegral . merge + where merge (x:xs) = let (a,b) = x `divMod` 256 in a:b:merge xs+ merge [] = []++-- ----------------------------------------------------------+-- | SFF has a 31-byte common header+--+-- The format is open to having the index anywhere between reads,+-- we should really keep count and check for each read. In practice, it+-- seems to be places after the reads.+-- +-- The following two fields are considered part of the header, but as+-- they are static, they are not part of the data structure+--+-- @ +-- magic :: Word32 -- 0x2e736666, i.e. the string \".sff\"+-- version :: Word32 -- 0x00000001+-- @+data CommonHeader = CommonHeader {+ index_offset :: Int64 -- ^ Points to a text(?) section+ , index_length, num_reads :: Int32+ , key_length, flow_length :: Int16+ , flowgram_fmt :: Word8+ , flow, key :: ByteString + }++instance Show CommonHeader where+ show (CommonHeader io il nr kl fl fmt f k) =+ "Common Header:\n\n" ++ (unlines $ map (" "++) + ["index_off:\t"++show io ++"\tindex_len:\t"++show il+ ,"num_reads:\t"++show nr+ ,"key_len:\t" ++show kl ++ "\tflow_len:\t"++show fl+ ,"format\t:" ++show fmt+ ,"flow\t:" ++BC.unpack f+ ,"key\t:" ++BC.unpack k+ , ""+ ])++instance Binary CommonHeader where+ get = do { m <- get ; when (m /= magic) $ error (printf "Incorrect magic number - got %8x, expected %8x" m magic)+ ; v <- get ; when (not (v `elem` versions)) $ error (printf "Unexpected version - got %d, supported are: %s" v (unwords $ map show versions))+ ; io <- get ; ixl <- get ; nrd <- get+ ; chl <- get ; kl <- get ; fl <- get ; fmt <- get+ ; fw <- getByteString (fromIntegral fl)+ ; k <- getByteString (fromIntegral kl)+ ; skip (chl-(31+fl+kl)) -- skip to boundary+ ; return (CommonHeader io ixl nrd kl fl fmt fw k)+ }++ put ch = let CommonHeader io il nr kl fl fmt f k = ch { index_offset = 0 } in+ do { let cl = 31+fl+kl+ l = cl `mod` 8+ padding = if l > 0 then 8-l else 0+ ; put magic; put (last versions); put io; put il; put nr; put (cl+padding); put kl; put fl; put fmt+ ; putByteString f; putByteString k+ ; pad padding -- skip to boundary+ }++-- ---------------------------------------------------------- +-- | Each Read has a fixed read header, containing various information.+data ReadHeader = ReadHeader {+ name_length :: Int16+ , num_bases :: Int32+ , clip_qual_left, clip_qual_right+ , clip_adapter_left, clip_adapter_right :: Int16+ , read_name :: ByteString+}++instance Show ReadHeader where+ show (ReadHeader nl nb cql cqr cal car rn) =+ (" Read Header:\n" ++) $ unlines $ map (" "++) + [ "name_len:\t"++show nl, "num_bases:\t"++show nb+ , "clip_qual:\t"++show cql++"..."++show cqr+ , "clip_adap:\t"++show cal++"..."++show car+ , "read name:\t"++BC.unpack rn+ , "" + ]++instance Binary ReadHeader where+ get = do+ { rhl <- get; nl <- get; nb <- get+ ; cql <- get; cqr <- get ; cal <- get ; car <- get+ ; n <- getByteString (fromIntegral nl)+ ; skip (rhl - (16+ nl))+ ; return (ReadHeader nl nb cql cqr cal car n)+ }+ put (ReadHeader nl nb cql cqr cal car rn) = + do { let rl = 16+nl+ l = rl `mod` 8+ padding = if l > 0 then 8-l else 0+ ; put (rl+padding); put nl; put nb; put cql; put cqr; put cal; put car+ ; putByteString rn + ; pad padding+ }++-- ----------------------------------------------------------+-- | This contains the actual flowgram for a single read.+data ReadBlock = ReadBlock {+ read_header :: ! ReadHeader+ -- The data block+ , flow_data :: ! ByteString -- nb! use unpackFlows for this+ , flow_index :: ! ByteString+ , bases :: ! SeqData+ , quality :: ! QualData+ }++-- | Helper function to access the flowgram+flowgram :: ReadBlock -> [Flow]+flowgram = unpackFlows . flow_data++-- | Extract the sequence with masked bases in lower case+masked_bases :: ReadBlock -> SeqData+masked_bases rb = let+ l = fromIntegral $ clip_qual_left $ read_header rb+ r = fromIntegral $ clip_qual_right $ read_header rb+ SeqData s = bases rb+ in SeqData $ LBC.concat [ LBC.map toLower $ LBC.take (l-1) s+ , LBC.take r (LBC.drop (l-1) s)+ , LBC.map toLower $ LBC.drop r s]++-- | Extract the index as absolute coordinates, not relative.+cumulative_index :: ReadBlock -> [Int]+cumulative_index = scanl1 (+) . map fromIntegral . B.unpack . flow_index++instance Show ReadBlock where+ show (ReadBlock h f i (SeqData b) (QualData q)) =+ show h ++ unlines (map (" "++) + ["flowgram:\t"++show (unpackFlows f)+ , "index:\t"++(concat . intersperse " " . map show . B.unpack) i+ , "bases:\t"++LBC.unpack b+ , "quality:\t"++(concat . intersperse " " . map show . LB.unpack) q+ , ""+ ])++-- ------------------------------------------------------------+-- | RSFF wraps an SFF to provide an instance of Binary with some more error checking.+data RSFF = RSFF { unRecovered :: SFF }++instance Binary RSFF where + get = do+ -- Parse CommonHeader+ chead <- get+ -- Get the first read block+ r1 <- do rh <- get + getRB chead rh+ -- Get subsequent read blocks+ rds <- replicateM (fromIntegral (num_reads chead))+ (do rh <- getSaneHeader (take 4 $ BC.unpack $ read_name $ read_header r1)+ getRB chead rh)+ return (RSFF $ SFF chead (r1:rds))+ put = error "You should not serialize an RSFF"++-- | This allows us to decode the constant parts of the read header for verifying its correcness.+data PartialReadHeader = PartialReadHeader {+ _pread_header_lenght :: Int16 + , _pname_length :: Int16+ , _pnum_bases :: Int32+ , _pclip_qual_left, _pclip_qual_right+ , _clip_adapter_left, _pclip_adapter_right :: Int16+ , _pread_name :: ByteString -- length four+}++instance Binary PartialReadHeader where+ get = do { rhl <- get; nl <- get; nb <- get; ql <- get; qr <- get; al <- get; ar <- get; rn <- getByteString 4 + ; return (PartialReadHeader rhl nl nb ql qr al ar rn) }+ put = error "You should not serialize a PartialReadHeader"++-- | Ensure that the header we're decoding matches our expectations.+getSaneHeader :: String -> Get ReadHeader+getSaneHeader prefix = do+ buf <- getLazyByteString 20+ decodeSaneH prefix buf ++-- | Decode a 'ReadHeader', verifying that the data make sense.+decodeSaneH :: String -> LBC.ByteString -> Get ReadHeader+decodeSaneH prefix buf = do+ let PartialReadHeader rhl nl _nb _ql _qr _al _ar rn = decode buf+ if rhl >= 20 && nl > 0 && all id (zipWith (==) prefix (BC.unpack rn))+ then do buf2 <- getLazyByteString (fromIntegral rhl-20)+ return (decode $ LB.concat [buf,buf2])+ else do x <- getLazyByteString 1 -- error "skip one byte, try again"+ decodeSaneH prefix (LBC.concat [buf,x])+
+ src/Bio/Sequence/SFF_filters.hs view
@@ -0,0 +1,144 @@+-- | This implements a number of filters used in the Titanium pipeline, +-- based on published documentation.+module Bio.Sequence.SFF_filters where++import Bio.Sequence.SFF (ReadBlock(..), ReadHeader(..)+ , flowToBasePos, flowgram, cumulative_index)++import Bio.Core.Sequence+import qualified Data.ByteString.Lazy as B+import qualified Data.ByteString as SB+import qualified Data.ByteString.Lazy.Char8 as BL+import Data.List (tails)+import Data.Char (toUpper)++-- Ti uses a set of filters, described in the (something) manual.+-- (GS Run Processor Application, section 3.2.2++)++-- ** Discarding filters++-- | DiscardFilters determine whether a read is to be retained or discarded+type DiscardFilter = ReadBlock -> Bool -- True to retain, False to discard++-- | This filter discards empty sequences.+discard_empty :: DiscardFilter+discard_empty rb = num_bases (read_header rb) >= 5++-- | Discard sequences that don't have the given key tag (typically TCAG) at the start+-- of the read.+discard_key :: String -> DiscardFilter+discard_key key rb = (map toUpper key==) $ take (length key) $ BL.unpack $ unSD $ bases rb++-- | 3.2.2.1.2 The "dots" filter discards sequences where the last positive flow is +-- before flow 84, and flows with >5% dots (i.e. three successive noise values) +-- before the last postitive flow. The percentage can be given as a parameter.+discard_dots :: Double -> DiscardFilter+discard_dots p rb = let dotcount = SB.length $ SB.filter (>3) $ flow_index rb+ in fromIntegral dotcount / fromIntegral (BL.length $ unSD $ bases rb) < p+ && last (cumulative_index rb) >= 84++-- | 3.2.2.1.3 The "mixed" filter discards sequences with more than 70% positive flows. +-- Also, discard with <30% noise, >20% middle (0.45..0.75) or <30% positive.+discard_mixed :: DiscardFilter+discard_mixed rb = let fs = dropWhile (<50) . reverse . flowgram $ rb+ fl = dlength fs+ in and+ [ (dlength (filter (>50) fs) / fl) < 0.7 -- 70% positive+ , (dlength (filter (<45) fs) / fl) > 0.3 -- 30% noise+ , (dlength (filter (>75) fs) / fl) > 0.3 -- 30% postivie+ , (dlength (filter (\f -> f<=75 && f>=45) fs) / fl) < 0.2+ ]++-- | Discard a read if the number of untrimmed flows is less than n (n=186 for Titanium)+discard_length :: Int -> DiscardFilter+discard_length n rb = length (flowgram rb) >= n++-- ** Trimming filters++-- | TrimFilters modify the read, typically trimming it for quality+type TrimFilter = ReadBlock -> ReadBlock++-- | 3.2.2.1.4 Signal intensity trim - trim back until <3% borderline flows (0.5..0.7).+-- Then trim borderline values or dots from the end (use a window).+trim_sigint :: TrimFilter+trim_sigint rb = clipSeq rb (sigint rb)++-- n counts the "bad" flow values, m counts flow position+sigint :: ReadBlock -> Int+sigint rb = let bs = drop 1 $ scanl (\(n,m,_) f -> if f >= 50 && f <= 70 then (n+1,m+1,f) else (n,m+1,f)) (0,0,0) $ flowgram rb + xs = dropWhile (\(_,_,f) -> f<=70) + $ dropWhile (\(n,m,_)->(1000*n) `div` m > (30::Int)) + $ reverse bs+ in case xs of [] -> error "no sequence left?"+ ((_,m,_):_) -> flowToBasePos rb m++-- | 3.2.2.1.5 Primer filter +-- This looks for the B-adaptor at the end of the read. The 454 implementation isn't very+-- effective at finding mutated adaptors.+trim_primer :: String -> TrimFilter+trim_primer s rb = clipSeq rb (find_primer s rb)++find_primer :: String -> ReadBlock -> Int+find_primer s rb = go (num_bases (read_header rb) - 10)+ where go i | i <= 5 = fromIntegral (num_bases $ read_header rb)+ | match i = fromIntegral i+ | otherwise = go (i-1)+ match j = s' `B.isPrefixOf` B.drop (fromIntegral j) (unSD $ bases rb)+ s' = BL.pack $ map toUpper $ take 14 s++-- 3.2.2.1.6 Trimback valley filter is ignored, we don't understand the description.++-- | 3.2.2.1.7 Quality score trimming trims using a 10-base window until a Q20 average is found.+trim_qual20 :: Int -> TrimFilter+trim_qual20 w rs = clipSeq rs $ qual20 w rs++qual20 :: Int -> ReadBlock -> Int+qual20 w rs = (fromIntegral $ num_bases $ read_header rs)+ - (length . takeWhile (<20) . map (avg . take w) . tails . reverse . B.unpack $ unQD $ quality rs)++-- ** Utility functions++-- | List length as a double (eliminates many instances of fromIntegral)+dlength :: [a] -> Double+dlength = fromIntegral . length++-- | Calculate average of a list+avg :: Integral a => [a] -> Double+avg xs = sum (map fromIntegral xs) / dlength xs++-- | Translate a number of flows to position in sequence, and update clipping data accordingly+clipFlows :: ReadBlock -> Int -> ReadBlock+clipFlows rb n = clipSeq rb (flowToBasePos rb n)++-- | Update clip_qual_right if more severe than previous value+clipSeq :: ReadBlock -> Int -> ReadBlock+clipSeq rb n' = let n = fromIntegral n' + rh = read_header rb + in if clip_qual_right rh <= n then rb else rb { read_header = rh {clip_qual_right = n }}++-- ** Data++-- Celera docs, at http://sourceforge.net/apps/mediawiki/wgs-assembler/index.php?title=SffToCA++-- These are used for mate-pair libraries, should be located around the middle of the read:++flx_linker = "GTTGGAACCGAAAGGGTTTGAATTCAAACCCTTTCGGTTCCAAC" -- Celera+ti_linker = "TCGTATAACTTCGTATAATGTATGCTATACGAAGTTATTACG" -- 20K cod jump++-- ti_linker and this: "AGCATATTGAAGCATATTACATACGATATGCTTCAATAATGC"+-- from "GS FLX Titanium 3 kb Span Paired End Library Preparation Method Manual April 2009"+-- ftp://ftp.genome.ou.edu/pub/for_broe/titanium/++-- These are used at the end of RNA (cDNA) sequences, after the poly-A tail:++rna_adapter = "ggcgggcgatgtctcgtctgagcgggctggcaaggc" -- cod transcripts?+rna_adapter2 = "ttcgcagtgagtgacaggctagtagctgagcgggctggcaaggc" -- Cod_c.sff+rna_adapter3 = "gacggggcggatgtctcgtctgagcgggcgtggcaaggc" -- COD1.sff++-- These are used at the end of DNA sequencing reads:++rapid_adapter = "agtcgtggaggcaaggcacacagggatagg" -- sea louse reads, key GACT+ti_adapter_b = "ctgagactgccaaggcacacagggggatagg" -- sea bass and l.s.Ca+ ++
+ src/Bio/Sequence/SFF_name.hs view
@@ -0,0 +1,86 @@+module Bio.Sequence.SFF_name where++import qualified Data.ByteString.Char8 as B+import Data.ByteString.Char8 (ByteString, pack)+import Data.Array.Unboxed+import Data.Char (ord)++-- | Read names encode various information, as per this struct.+data ReadName = ReadName { date :: (Int,Int,Int)+ , time :: (Int,Int,Int)+ , region :: Int+ , x_loc, y_loc :: Int } deriving Show++-- ----------------------------------------------------------+-- Decoding++decodeReadName :: ByteString -> Maybe ReadName+decodeReadName b = do t <- decodeDate $ B.take 6 b+ r <- fst `fmap` (B.readInt $ B.take 2 $ B.drop 7 b)+ l <- decodeLocation $ B.drop 9 b+ return $ ReadName { date = (\[y,m,d] -> (y,m,d)) (take 3 t)+ , time = (\[hh,mm,ss] -> (hh,mm,ss)) (drop 3 t)+ , region = r+ , x_loc = fst l, y_loc = snd l }++decodeLocation :: ByteString -> Maybe (Int,Int)+decodeLocation l = (`divMod` 4096) `fmap` decode36 l++decodeDate :: ByteString -> Maybe [Int]+decodeDate d = (fixyear . reverse . (`divMods` [60,60,24,32,13])) =<< decode36 d+ where fixyear (i:is) = Just (2000+i:is)+ fixyear [] = Nothing++-- ----------------------------------------------------------+-- Encoding++encodeReadName :: ReadName -> ByteString+encodeReadName r = B.concat [ encodeDate (date r) (time r) + , encodeRegion (region r)+ , encodeLocation (x_loc r) (y_loc r)]++encodeLocation :: Int -> Int -> ByteString+encodeLocation = undefined++encodeRegion :: Int -> ByteString+encodeRegion = undefined++encodeDate :: (Int,Int,Int) -> (Int,Int,Int) -> ByteString+encodeDate = undefined++-- ----------------------------------------------------------++divMods :: Int -> [Int] -> [Int]+divMods x (i:is) = let (a,b) = x `divMod` i+ in b : divMods a is+divMods x [] = [x]++-- ----------------------------------------------------------+-- Decoding base36 strings++decode36 :: ByteString -> Maybe Int+decode36 s = (foldr1 (\a b -> b*36+a) . reverse) `fmap` (mapM decCh . B.unpack $ s)++{-+decode36' = dec 0+ where dec i b = case uncons b of Just (c,rest) -> dec (i*36+fromJust (decCh c)) rest+ Nothing -> i+ fromJust (Just z) = z+-}++decCh :: Char -> Maybe Int+decCh x | x >= 'A' && x <= 'Z' = Just (ord x - ord 'A')+ | x >= '0' && x <= '9' = Just (26 + ord x - ord '0')+ | otherwise = Nothing -- error ("decode36: can't decode "++show x)++encode36 :: Int -> ByteString+encode36 = pack . map (b36!) . reverse . enc+ where+ enc 0 = []+ enc i = let (a,b) = i `divMod` 36+ in b : enc a++b36 :: UArray Int Char+b36 = listArray (0,35) (['A'..'Z']++['0'..'9'])++
+ src/Flower/Fork.hs view
@@ -0,0 +1,14 @@+module Fork where++import Control.Concurrent+import Control.Exception++-- | Spawn a set of threads and wait for them to complete.+forkAndWait :: [IO ()] -> IO ()+forkAndWait actions = mapM myForkIO actions >>= mapM_ takeMVar+ where+ myForkIO :: IO () -> IO (MVar ())+ myForkIO io = do+ mvar <- newEmptyMVar+ _ <- forkIO (io `finally` putMVar mvar ())+ return mvar
+ src/Flower/Main.hs view
@@ -0,0 +1,245 @@+-- FlowEr - FLOWgram ExtractoR+module Main (main) where++import Bio.Sequence.SFF+import Bio.Sequence.SFF_filters+import Bio.Core++import Print+import Text.Printf++import System.IO (stdout, Handle, openFile, IOMode(..), hClose, hPutStrLn)++import Numeric (showFFloat)+import Data.Char (toLower)+import Data.List (intersperse)+import Data.ByteString.Char8 (unpack,ByteString)+import qualified Data.ByteString.Char8 as B+import qualified Data.ByteString as B1+import qualified Data.ByteString.Lazy.Char8 as L+import qualified Data.ByteString.Lazy as L1++import Data.Array.Unboxed+import Data.Array.ST+import Control.Monad.ST+import Control.Monad.State++import Metrics+import qualified Options as O+import Options (Opts)+import Fork++main :: IO ()+main = do+ opts <- O.getArgs+ when (null $ O.inputs opts) $ error "Please provide an input file - or use --help for more information."+ forkAndWait $ buildActions opts++type Action = IO ()+type Trimmer = ReadBlock -> ReadBlock++buildActions :: Opts -> [Action]+buildActions o = let+ inp = mapM readSFF (O.inputs o)+ tr = map (mkTrimmer o)+ ch (SFF h _) = h+ rs (SFF _ r) = r+ in snd $ flip runState [] $ do+ on (O.info o) (\h -> mapM_ (hPutStrLn h . getHeader . ch) =<< inp)+ on (O.fasta o) (\h -> mapM_ (L1.hPut h . L1.concat . map toFasta . tr . rs) =<< inp)+ on (O.fqual o) (\h -> mapM_ (L1.hPut h . L1.concat . map toFastaQual . tr . rs) =<< inp)+ on (O.text o) (\h -> mapM_ (hPutStrLn h . dumpText . tr . rs) =<< inp)+ on (O.fastq o) (\h -> mapM_ (L1.hPut h . L1.concat . map toFastQ . tr . rs) =<< inp)+ on (O.summarize o) (\h -> mapM_ (L1.hPut h . summarize . tr . rs) =<< inp) -- should we trim?+ on (O.filters o) (\h -> mapM_ (L1.hPut h . sum_filters . rs) =<< inp)+ on (O.histogram o) (\h -> mapM_ (\(SFF c r) -> hPutStrLn h . showHist . histogram (B.unpack $ flow c) . map flowgram . tr $ r) =<< inp)+ on (O.flowgram o) (\h -> mapM_ (\(SFF c r) -> L1.hPut h . L1.fromChunks . intersperse (B.pack "\n") . concatMap (showread c) $ r) =<< inp)++on :: Maybe FilePath -> (Handle -> Action) -> State [Action] ()+on Nothing _ = return ()+on (Just f) act = modify $ (:) $ case f of + "-" -> act stdout+ _ -> do h <- openFile f WriteMode+ act h+ hClose h++mkTrimmer :: Opts -> Trimmer+mkTrimmer o = case (O.trimKey o, O.trim o) of+ (True,True) -> error "Please specify only one of --trim and --trimkey"+ (True,False) -> \r -> trimFromTo 5 (num_bases $ read_header r) r+ (False,True) -> trim+ (False,False) -> id++-- ------------------------------------------------------------+-- No option - dump as text format+-- ------------------------------------------------------------+dumpText :: [ReadBlock] -> String+dumpText rs = concat . map toText $ rs+ where toText :: ReadBlock -> String+ toText r = concat [ gt, B.unpack (read_name rh), nl+ , maybe "" ((\s->info++s++nl) . formatRN) $ decodeReadName (read_name rh)+ , let (lf,rt) = (clip_adapter_right rh, clip_adapter_left rh) + in if lf /= 0 || rt /= 0 then adapter ++ show lf ++ sp++ show rt else ""+ , clip, show (clip_qual_left rh), sp, show (clip_qual_right rh), nl+ , flows, B.unpack $ B.unwords $ map fi $ flowgram r, nl+ , idx, unwords $ map show $ cumulative_index' r, nl+ , base, L.unpack (unSD $ seqdata r), nl+ , qual, unwords $ map show $ L1.unpack (unQD $ quality r), nl+ ]+ where rh = read_header r+ gt = ">"+ nl = "\n"+ sp = " "+ info = " Info: \t"+ clip = " Clip: \t"+ adapter = " Adap: \t"+ flows = " Flows:\t"+ idx = " Index:\t"+ base = " Bases:\t"+ qual = " Quals:\t"+ formatRN (ReadName (yr,mo,dy) (h,m,s) r' x y) = + printf "%4d-%02d-%02d %02d:%02d:%02d R%d (%d,%d)" yr mo dy h m s r' x y++cumulative_index' :: ReadBlock -> [Int]+cumulative_index' = scanl1 (+) . map fromIntegral . B1.unpack . flow_index++-- ------------------------------------------------------------+-- The -i option: Print header info+-- ------------------------------------------------------------+getHeader :: CommonHeader -> String+getHeader 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 :: [ReadBlock] -> L.ByteString+summarize rs = do+ L.concat [ L.pack "# name........\tdate......\ttime....\treg\ttrim_l\ttrim_r\tx_loc\ty_loc\tlen\tK2\ttrimK2\tncount\tavgQ\ttravgQ\n"+ , 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+ tr = trim r+ tb, nl, q :: Builder+ tb = char '\t'+ nl = char '\n'+ q = char '?'+ + (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 qs = let l = fromIntegral (L1.length qs)+ in if l>0 then putFix 2 $ sum (map fromIntegral $ L1.unpack qs) * 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 $ flowgram tr)+ , putInt (n_count r)+ , avg_qual $ unQD $ quality r, avg_qual $ unQD $ quality tr]) ++ [nl]++-- ----------------------------------------------------------+-- The --filters option, summarize filters+-- ----------------------------------------------------------+sum_filters :: [ReadBlock] -> L1.ByteString+sum_filters rs = toLazyByteString $ mconcat (header:map sumf1 rs)+ where + header = fromByteString $ B.pack "# name..... \tlength \tl_trim \tr_trim \tE K D M L\tSig Q20 Adp\n"+ sumf1 rb = let+ rh = read_header rb+ rn = read_name rh+ nb = fromIntegral $ num_bases rh+ (cl,cr) = (fromIntegral $ clip_qual_left rh, fromIntegral $ clip_qual_right rh)+ dfs = mconcat $ intersperse (char ' ') $+ map (\f -> if f rb then char '+' else char ' ') + [discard_empty, discard_key "tcag", discard_dots 0.05, discard_mixed, discard_length 186]+ tfs = mconcat $ intersperse (char ' ') $ map (\f -> putInt3 (f rb))+ [sigint, qual20 10, find_primer rapid_adapter]+ in mconcat (intersperse (char '\t') [fromByteString rn, putInt nb, putInt cl, putInt cr, dfs, tfs]++[char '\n'])++-- ----------------------------------------------------------+-- The -F option: Output the sequence of flows, one flow per line+-- ----------------------------------------------------------++fi :: Flow -> ByteString+fi f | f <= 9999 && f >= 0 = farray!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]]++tab :: ByteString+tab = B.pack "\t"++showread :: CommonHeader -> ReadBlock -> [ByteString]+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) (map Qual $ L1.unpack $ unQD $ 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 (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 :: ReadHeader -> Int -> [[a]] -> [Char] -> [Char]+mask _ _ _ [] = [] -- 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+mask _ _ _ _ = error "internal error in 'mask'"++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 _ _ _ _ _ = []++-- | Take the unpacked index_offsets and quality values, and return +-- a list of groups of quality values, each group corresponding to a flow value. +-- Flow values < 0.5 result in empty groups.+qgroup :: [Index] -> [Qual] -> [[Qual]]+qgroup [] [] = let rest = []:rest in rest+qgroup is@(1:_) qs = let (iz,irest) = span (==0) (tail is)+ (q1,qrest) = splitAt (length iz+1) qs+ in q1 : qgroup irest qrest+qgroup (i:is) qs = [] : qgroup (i-1:is) qs+qgroup _ _ = error "internal error in 'qgroup'"++-- ----------------------------------------------------------+-- The -h option: Output a histogram of flow values+-- ----------------------------------------------------------++type Hist = UArray Flow Int++histogram :: String -> [[Flow]] -> (Hist,Hist,Hist,Hist)+histogram fl scores = runST $ do + let zero = newArray (0,9999) 0 :: ST s (STUArray s Flow Int)+ a <- zero+ c <- zero+ g <- zero+ t <- zero+ let ins1 ('A',i) = bump a i+ ins1 ('C',i) = bump c i+ ins1 ('G',i) = bump g i+ 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) (map (\x->if x>9999 || x<0 then 9999 else x) $ concat scores))+ a' <- unsafeFreeze a+ c' <- unsafeFreeze c+ g' <- unsafeFreeze g+ t' <- unsafeFreeze t+ return (a',c',g',t')++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]+ | sc <- [0..9999]]
+ src/Flower/Metrics.hs view
@@ -0,0 +1,23 @@+-- Calculate various characteristics on sequence quality++module Metrics where++import Bio.Core+import Bio.Sequence.SFF+import qualified Data.ByteString.Lazy.Char8 as B++-- import Test.QuickCheck++-- | Take the fractional parts of the flows, and sum their squares (the "K²" metric)+quals :: [Flow] -> Flow+quals q = floor $ ((100::Double) - 2*(sqrt $ (/fromIntegral (length q)) $ sum $ map (fromIntegral . (^(2::Integer)) . (flip (-) 50) . (`mod` 100) . (+50)) $ q))++-- | 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 . B.unpack . unSD . 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/Flower/Options.hs view
@@ -0,0 +1,52 @@+{-# LANGUAGE DeriveDataTypeable #-}++module Options where++import System.Console.CmdArgs+import Control.Monad (when)+import Data.Maybe (isJust)++data Opts = Opts + { trimKey :: Bool+ , trim :: Bool+ , summarize :: Maybe FilePath+ , filters :: Maybe FilePath+ , info :: Maybe FilePath+ , fasta :: Maybe FilePath+ , fqual :: Maybe FilePath+ , fastq :: Maybe FilePath+ , flowgram :: Maybe FilePath+ , histogram :: Maybe FilePath+ , inputs :: [FilePath]+ , text :: Maybe FilePath+ } deriving (Data,Typeable, Show, Eq)++optdef :: Ann+optdef = opt ("-"::String)++opts :: Opts+opts = Opts+ { trimKey = False &= help "Trim only the TCAG key sequence"+ , trim = False &= help "Trim quality using clipping information" &= name "t"+ , summarize = def &= help "Output per sequence summary information" &= typFile &= optdef+ , filters = def &= help "Output filtering information" &= typFile &= optdef+ , info = def &= help "Output brief overview of the contents" &= typFile &= optdef+ , fasta = def &= help "Output FASTA-formatted sequences" &= typFile &= name "f" &= optdef+ , fqual = def &= help "Output phred qualities" &= typFile &= name "q" &= optdef+ , fastq = def &= help "Output FastQ-formatted sequence and Sanger quality" &= typFile &= name "Q" &= optdef+ , flowgram = def &= help "Output flowgram information in tabular form" &= typFile &= name "F" &= optdef+ , histogram = def &= help "Output histogram of flow values" &= typFile &= name "h" &= optdef+ , text = def &= help "Output SFF information as text (default)" &= typFile &= name "T" &= optdef+ , inputs = def &= args &= typFile+ } + &= summary "flower v0.7 - Extract information from SFF files" + &= program "flower"++getArgs :: IO Opts+getArgs = do+ o <- cmdArgs opts + -- print o+ let outs = filter isJust $ map ($o) [summarize,filters,info,fasta,fqual,fastq,flowgram,histogram,text]+ when ((length $ filter (==Just "-") $ outs) > 1) $ error "If you specify more than one output format, you need to specify output files"+ let o' = if null outs then o { text = Just "-" } else o+ return o'
+ src/Flower/Print.hs view
@@ -0,0 +1,51 @@+module Print + (+ Builder, toLazyByteString, mconcat, fromByteString, char+ , putInt, putInt2, putInt3, putDate, putTime, putFix+ ) where ++import Data.Binary.Builder+import Data.Monoid+import Data.ByteString.Char8 (ByteString, pack)+import Data.Array.Unboxed+import Data.Char (ord)++char :: Char -> Builder+char = singleton . fromIntegral . ord++putInt :: Int -> Builder+putInt y = if y < 0 then char '-' `append` putInt' (negate y)+ else putInt' y+ where putInt' x = case x `divMod` 1000 of+ (0,r) -> fromByteString (ints!r)+ (a,r) -> putInt a `append` putInt3 r++-- zero padded ints++putInt2 :: Int -> Builder+putInt2 x | x<100 = fromByteString (int2s!x)+ | otherwise = fromByteString (pack "xx")++putInt3 :: Int -> Builder+putInt3 x | x<1000 = fromByteString (int3s!x)+ | otherwise = fromByteString (pack "xxx")++ints, int2s, int3s :: 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])++putDate :: Int -> Int -> Int -> Builder+putDate y m d = mconcat [putInt y, dash, putInt2 m, dash, putInt2 d]+ where dash = char '-'++putTime :: Int -> Int -> Int -> Builder+putTime h m s = mconcat [putInt2 h, col, putInt2 m, col, putInt2 s]+ where col = char ':'++-- 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"