diff --git a/biosff.cabal b/biosff.cabal --- a/biosff.cabal +++ b/biosff.cabal @@ -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 diff --git a/src/Bio/Sequence/SFF.hs b/src/Bio/Sequence/SFF.hs --- a/src/Bio/Sequence/SFF.hs +++ b/src/Bio/Sequence/SFF.hs @@ -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" diff --git a/src/Flower/Main.hs b/src/Flower/Main.hs --- a/src/Flower/Main.hs +++ b/src/Flower/Main.hs @@ -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 diff --git a/src/Flower/Options.hs b/src/Flower/Options.hs --- a/src/Flower/Options.hs +++ b/src/Flower/Options.hs @@ -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'