biosff 0.1 → 0.2
raw patch · 4 files changed
+54/−21 lines, 4 filesdep ~binary
Dependency ranges changed: binary
Files
- biosff.cabal +12/−3
- src/Bio/Sequence/SFF.hs +21/−8
- src/Flower/Main.hs +16/−7
- src/Flower/Options.hs +5/−3
biosff.cabal view
@@ -1,5 +1,5 @@ Name: biosff-Version: 0.1+Version: 0.2 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@@ -11,12 +11,21 @@ Stability: Experimental Category: Bioinformatics Build-type: Simple-Cabal-version: >=1.2+Cabal-version: >=1.6 +Source-repository head+ Type: darcs+ Location: http://malde.org/~ketil/biohaskell/biosff++source-repository this+ Type: darcs+ Location: http://malde.org/~ketil/biohaskell/biosff+ Tag: 0.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+ Build-depends: base >= 3 && < 5, biocore >= 0.1, binary, bytestring, array Hs-Source-Dirs: src Ghc-Options: -Wall
src/Bio/Sequence/SFF.hs view
@@ -25,6 +25,7 @@ , packFlows, unpackFlows , Flow, Qual, Index, SeqData, QualData , ReadName (..), decodeReadName, encodeReadName+ , putRB, getRB ) where import Bio.Core.Sequence@@ -40,7 +41,7 @@ import Data.List (intersperse) import Data.Binary-import Data.Binary.Get (getByteString,getLazyByteString)+import Data.Binary.Get (getByteString,getLazyByteString,runGetState) import qualified Data.Binary.Get as G import Data.Binary.Put (putByteString,putLazyByteString) import Data.Char (toUpper, toLower)@@ -62,8 +63,21 @@ -- | Read an SFF file. readSFF :: FilePath -> IO SFF-readSFF f = return . decode =<< LB.readFile f+readSFF f = do+ file <- LB.readFile f+ let (header, remaining, _) = runGetState (get::Get CommonHeader) file 0+ blocks = getBlocks header (fromIntegral $ num_reads header) remaining+ return (SFF header blocks) +getBlocks :: CommonHeader -> Int -> LB.ByteString -> [ReadBlock]+getBlocks header n str+ | n == 0 = []+ | otherwise = case runGetState (getBlock $ fromIntegral $ flow_length $ header) str 0 of+ (block, remaining, _) -> block : getBlocks header (n-1) remaining++getBlock :: Int -> Get ReadBlock+getBlock flows = get >>= getRB flows+ {- -- | Extract the read without the initial (TCAG) key. trimKey :: CommonHeader -> Sequence Nuc -> Maybe (Sequence Nuc)@@ -227,7 +241,7 @@ rds <- replicateM (fromIntegral (num_reads chead)) (do rh <- get :: Get ReadHeader- getRB chead rh+ getRB (fromIntegral $ flow_length chead) rh ) return (SFF chead rds) @@ -237,11 +251,10 @@ -- | Helper function for decoding a 'ReadBlock'. {-# INLINE getRB #-}-getRB :: CommonHeader -> ReadHeader -> Get ReadBlock-getRB chead rh = do+getRB :: Int -> ReadHeader -> Get ReadBlock+getRB fl 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'@@ -418,11 +431,11 @@ chead <- get -- Get the first read block r1 <- do rh <- get - getRB chead rh+ getRB (fromIntegral $ flow_length 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)+ getRB (fromIntegral $ flow_length chead) rh) return (RSFF $ SFF chead (r1:rds)) put = error "You should not serialize an RSFF"
src/Flower/Main.hs view
@@ -52,7 +52,8 @@ 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.histogram o) (\h -> mapM_ (\(SFF c r) -> hPutStrLn h . showHist 9999 ["A","C","G","T"] . histogram (B.unpack $ flow c) . map flowgram . tr $ r) =<< inp)+ on (O.histpos o) (\h -> mapM_ (\(SFF _ r) -> hPutStrLn h . showHist 549 [] . histpos 549 . 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] ()@@ -219,7 +220,7 @@ type Hist = UArray Flow Int -histogram :: String -> [[Flow]] -> (Hist,Hist,Hist,Hist)+histogram :: String -> [[Flow]] -> [Hist] histogram fl scores = runST $ do let zero = newArray (0,9999) 0 :: ST s (STUArray s Flow Int) a <- zero@@ -237,9 +238,17 @@ c' <- unsafeFreeze c g' <- unsafeFreeze g t' <- unsafeFreeze t- return (a',c',g',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]] +showHist :: Int -> [String] -> [Hist] -> String+showHist mx hd hs = (if not (null hd) then concat (intersperse "\t" ("Score":hd)) ++ "\n" else "") +++ unlines [concat $ intersperse "\t" $ showFFloat (Just 2) (fromIntegral sc/100::Double) "" : [show (h!sc) | h <- hs] | sc <- [0..fromIntegral mx]]++histpos :: Int -> [ReadBlock] -> [Hist]+histpos mx scores = runST $ do + let mx' = fromIntegral mx+ zero = newArray (0,mx') 0 :: ST s (STUArray s Flow Int)+ hs <- sequence $ replicate (length $ flowgram $ head scores) zero+ let bump ar i = when (i>=0 && i<= mx') (readArray ar i >>= \x -> writeArray ar i (x+1))+ sequence_ $ concatMap (zipWith bump hs . flowgram) scores+ mapM unsafeFreeze hs
src/Flower/Options.hs view
@@ -17,6 +17,7 @@ , fastq :: Maybe FilePath , flowgram :: Maybe FilePath , histogram :: Maybe FilePath+ , histpos :: Maybe FilePath , inputs :: [FilePath] , text :: Maybe FilePath } deriving (Data,Typeable, Show, Eq)@@ -35,18 +36,19 @@ , 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+ , histogram = def &= help "Output histogram of flow values by nucleotide" &= typFile &= name "h" &= optdef+ , histpos = def &= help "Output histogram of flow values by flow cycle" &= 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" + &= summary "flower (biosff v0.2) - 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]+ let outs = filter isJust $ map ($o) [summarize,filters,info,fasta,fqual,fastq,flowgram,histogram,histpos,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'