bio 0.4.7 → 0.4.8

raw patch · 8 files changed

+130/−46 lines, 8 filesdep +directorydep ~tagsoupPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependencies added: directory

Dependency ranges changed: tagsoup

API changes (from Hackage documentation)

Files

• Bio/Alignment/BlastXML.hs +8/−6
• Bio/Sequence.hs +26/−3
• Bio/Sequence/FastQ.hs +45/−13
• Bio/Sequence/Fasta.hs +14/−6
• Bio/Sequence/SFF.hs +28/−6
• Bio/Sequence/SeqData.hs +2/−2
• Bio/Sequence/TwoBit.hs +3/−6
• bio.cabal +4/−4

@@ -20,6 +20,8 @@ import Control.Parallel import Data.List (isPrefixOf) +type STag = Tag String+ -- | Parse BLAST results in XML format readXML :: FilePath -> IO [BlastResult] readXML fp = do @@ -38,7 +40,7 @@                   in  rest `par` first : if null rest then [] else breaks p rest breaks _ []     = [] -getFrom :: [Tag] -> String -> String+getFrom :: [STag] -> String -> String getFrom list tag = let xs = sections (isTagOpenName tag) list                     in if null xs || null (head xs) || (null . drop 1 . head) xs                        then error ("Couldn't find tag '"++show tag++"' in\n"++showSome list)@@ -47,11 +49,11 @@                              x -> error ("Unexpeced tag: "++ show x)  -- Use pattern match since 'length' is strict, defeating the purpose.-showSome :: [Tag] -> String+showSome :: [STag] -> String showSome a@(_:_:_:_:_:_:_) = (init . show . take 5 $ a)++" ... ]" showSome a                 = show a -xml2br :: [Tag] -> [[Tag]] -> BlastResult+xml2br :: [STag] -> [[STag]] -> BlastResult xml2br h is = BlastResult { blastprogram = get "BlastOutput_program"                           , blastversion = bv                           , blastdate = bd @@ -64,7 +66,7 @@     where (bv,bd) = B.break (=='[') $ get "BlastOutput_version"           get = B.pack . getFrom h -iter2rec :: [[Tag]] -> BlastRecord+iter2rec :: [[STag]] -> BlastRecord iter2rec (i:hs) = BlastRecord                { query = B.pack $ get "Iteration_query-def"               , qlength = read $ get "Iteration_query-len"@@ -73,7 +75,7 @@     where get = getFrom i iter2rec [] = error "iter2rec: got empty list of sections!" -hit2hit :: [Tag] -> BlastHit+hit2hit :: [STag] -> BlastHit hit2hit hs = BlastHit               { subject = B.pack $ get "Hit_def"              , slength = read $ get "Hit_len"@@ -82,7 +84,7 @@     where get = getFrom hs  -hsp2match :: [Tag] -> BlastMatch+hsp2match :: [STag] -> BlastMatch hsp2match ms = BlastMatch                { bits   = read $ get "Hsp_bit-score"                , e_val  = read $ get "Hsp_evalue"

@@ -38,6 +38,8 @@     -- Combines sequence data and quality in one file.     -- Warning: Solexa uses a different formula for the quality values!     , readFastQ, writeFastQ, hReadFastQ, hWriteFastQ+    , readSangerQ, writeSangerQ, hReadSangerQ, hWriteSangerQ+    , readIllumina, writeIllumina, hReadIllumina, hWriteIllumina      -- ** The phd file format ("Bio.Sequence.Phd")     -- | These contain base (nucleotide) calling information,@@ -66,16 +68,37 @@ import Bio.Sequence.FastQ import Bio.Sequence.Phd import Bio.Sequence.TwoBit+import Bio.Sequence.SFF  -- sequence-oriented stuff import Bio.Sequence.Entropy import Bio.Sequence.HashWord +import Control.Monad (filterM)+import System.Directory (doesFileExist)+ -- | Read nucleotide sequences in any format - Fasta, SFF, FastQ, 2bit, PHD...+--   Todo: read quality automatically if available readNuc :: FilePath -> IO [Sequence Nuc]-readNuc = undefined-          -- check file contents-          -- magic number+readNuc fp  +  | ext `elem` ["fasta", "fna", "fa", "fst"] = do +       ps <- findQual fp +       ss <- (case ps of [q] -> readFastaQual fp q+                         []  -> readFasta fp+                         qs  -> error ("Ambigous quality file for "++show fp++": "++show qs))+       return (map castSeq ss)+  | ext == "2bit"                     = read2Bit                               $ fp+  | ext == "sff"                      = fmap sffToSequence . readSFF           $ fp+  | ext `elem` ["fq","fastq"]         = readFastQ                              $ fp+  | ext2 == "phd"                     = fmap return . readPhd                  $ fp -- only a single sequence+  -- "ace" ?+  | otherwise                         = error "readNuc: unknown file suffix!"+  where+    ext = reverse . takeWhile (/='.') . reverse $ fp+    ext2 = reverse . takeWhile (/='.') . dropWhile (=='.') . dropWhile (/='.') . reverse $ fp+    basename = reverse . dropWhile (=='.') . dropWhile (/= '.') . reverse+    findQual = filterM doesFileExist . qualnames+    qualnames f = [f++".qual",basename f++".qual"]  -- | Read protein sequences in any supported format (i.e. Fasta) readProt :: FilePath -> IO [Sequence Amino]

@@ -5,26 +5,37 @@     * <http://www.bioperl.org/wiki/FASTQ_sequence_format> -   Of course, this is yet another vaguely defined pseudo-standard with -   conflicting definitions.  Of course Solexa had to go and invent a different, -   but indistinguishably so, way to do it:+   Of course, this is yet another vaguely defined pseudo-standard with+   conflicting definitions.  Of course Solexa had to go and invent not one, but two +   different, and indistinguishably so, ways to do it:     * <http://www.bcgsc.ca/pipermail/ssrformat/2007-March/000137.html>     * <http://maq.sourceforge.net/fastq.shtml> -   Currently, we only support the non-Solexa FastQ, adding\/subtracting 33 for -   the quality values.+   * <http://en.wikipedia.org/wiki/FASTQ_format> +   Sanger-style FastQ-format is supported with the (h)read/writeSangerQ functions,+   and the new Illumina/Solexa-style with (h)read/writeIllumina.+    As far as I know, FastQ is only used for nucleotide sequences, never amino acid. -} -module Bio.Sequence.FastQ ++module Bio.Sequence.FastQ     (      -- * Reading FastQ     readFastQ, hReadFastQ, parse      -- * Writing FastQ     , writeFastQ, hWriteFastQ, unparse++    -- * use Sanger-style quality information+    , readSangerQ, hReadSangerQ+    , writeSangerQ, hWriteSangerQ+    +    -- * use Illumina (>v1.3)-style quality information+    , readIllumina, hReadIllumina+    , writeIllumina, hWriteIllumina     ) where  import System.IO@@ -34,8 +45,29 @@  import Bio.Sequence.SeqData +{-# DEPRECATED readFastQ, hReadFastQ, writeFastQ, hWriteFastQ  "FastQ assumes Sanger-style quality info use {read,write}SangerQ or -Illumina instead" #-}++readSangerQ, readIllumina :: FilePath -> IO [Sequence Nuc]+readSangerQ = readFastQ+readIllumina f = addQual (negate 31) `fmap` readFastQ f++hReadSangerQ, hReadIllumina :: Handle -> IO [Sequence Nuc]+hReadSangerQ = hReadFastQ+hReadIllumina h = addQual (negate 31) `fmap` hReadFastQ h++writeSangerQ, writeIllumina :: FilePath -> [Sequence Nuc] -> IO ()+writeSangerQ = writeFastQ+writeIllumina f = writeFastQ f . addQual 31 ++hWriteSangerQ, hWriteIllumina :: Handle -> [Sequence Nuc] -> IO ()+hWriteSangerQ = hWriteFastQ+hWriteIllumina h = hWriteFastQ h . addQual 31++addQual :: Qual -> [Sequence Nuc] -> [Sequence Nuc]+addQual q = map (\(Seq h d mq) -> (Seq h d $ BB.map (+q) `fmap` mq))+ readFastQ :: FilePath -> IO [Sequence Nuc]-readFastQ f = (go . B.lines) `fmap` B.readFile f +readFastQ f = (go . B.lines) `fmap` B.readFile f  hReadFastQ :: Handle -> IO [Sequence Nuc] hReadFastQ h = (go . B.lines) `fmap` B.hGetContents h@@ -46,7 +78,7 @@ -- | Parse one FastQ entry, suitable for using in 'unfoldr' over --   'B.lines' from a file parse :: [B.ByteString] -> Maybe (Either String (Sequence Nuc), [B.ByteString])-parse (h1:sd:h2:sq:rest) = +parse (h1:sd:h2:sq:rest) =     case (B.uncons h1,B.uncons h2) of       (Just ('@',h1name), Just ('+',h2name))           | h1name == h2name -> Just (Right $ Seq h1name sd (Just (BB.map (subtract 33) sq)), rest)@@ -58,15 +90,15 @@                err = Left $ "Bio.Sequence.FastQ: illegal number of lines in FastQ format: " ++ showStanza            in Just (err, []) -writeFastQ :: FilePath -> [Sequence a] -> IO ()+writeFastQ :: FilePath -> [Sequence Nuc] -> IO () writeFastQ f = B.writeFile f . B.concat . map unparse -hWriteFastQ :: Handle -> [Sequence a] -> IO ()+hWriteFastQ :: Handle -> [Sequence Nuc] -> IO () hWriteFastQ h = B.hPut h . B.concat . map unparse  -- helper function for writing-unparse :: Sequence a -> B.ByteString-unparse (Seq h sd (Just sq)) = +unparse :: Sequence Nuc -> B.ByteString+unparse (Seq h sd (Just sq)) =     B.unlines [B.cons '@' h, sd, B.cons '+' h, BB.map (+33) sq]-unparse (Seq h _ Nothing) = error ("Bio.Sequence.FastQ: sequence " ++ show (B.unpack h) +unparse (Seq h _ Nothing) = error ("Bio.Sequence.FastQ: sequence " ++ show (B.unpack h)                                    ++" doesn not have quality data!")

@@ -28,7 +28,7 @@ ) where  --- import Data.Char (isSpace)+import Data.Char (chr) -- isSpace import Data.List (groupBy,intersperse) import System.IO @@ -59,8 +59,13 @@  -- | Lazily read sequences from a FASTA-formatted file readFasta :: FilePath -> IO [Sequence Unknown]-readFasta f = (mkSeqs . B.lines) `fmap` B.readFile f+readFasta f = (mkSeqs . blines) `fmap` B.readFile f +-- | Replacement for 'lines' that gobbles control-M's+-- Some tools, like CLC, likes to add these to the end of each line.+blines :: B.ByteString -> [B.ByteString]+blines = B.lines . B.filter (/=Data.Char.chr 13)+ -- | Write sequences to a FASTA-formatted file. --   Line length is 60. writeFasta :: FilePath -> [Sequence a] -> IO ()@@ -71,7 +76,7 @@  -- | Read quality data for sequences to a file. readQual :: FilePath -> IO [Sequence Unknown]-readQual f = (mkQual . B.lines) `fmap` B.readFile f+readQual f = (mkQual . blines) `fmap` B.readFile f  -- | Write quality data for sequences to a file. writeQual :: FilePath -> [Sequence a] -> IO ()@@ -115,7 +120,7 @@  -- | Lazily read sequence from handle hReadFasta :: Handle -> IO [Sequence Unknown]-hReadFasta h = (mkSeqs . B.lines) `fmap` B.hGetContents h+hReadFasta h = (mkSeqs . blines) `fmap` B.hGetContents h  -- | Write sequences in FASTA format to a handle. hWriteFasta :: Handle -> [Sequence a] -> IO ()@@ -153,7 +158,10 @@ mkSeqs = map mkSeq . blocks  mkSeq :: [ByteString] -> Sequence Unknown-mkSeq (l:ls) = Seq (B.drop 1 l) (B.concat $ takeWhile isSeq ls) Nothing+mkSeq (l:ls) +  -- maybe check this?  | B.length l < 2 || isSpace (B.head $ B.tail l) +  --  = error "Trying to read sequence without a name...and failing."+  | otherwise = Seq (B.drop 1 l) (B.concat $ takeWhile isSeq ls) Nothing     where isSeq s = (not . B.null) s && ((flip elem) (['A'..'Z']++['a'..'z']) . B.head) s mkSeq [] = error "empty input to mkSeq" @@ -177,6 +185,6 @@ countSeqs :: FilePath -> IO Int countSeqs f = do   ss <- B.readFile f-  let hdrs = filter (('>'==).B.head) $ filter (not . B.null) $ B.lines ss+  let hdrs = filter (('>'==).B.head) $ filter (not . B.null) $ blines ss   return (length hdrs) 

@@ -14,9 +14,11 @@ module Bio.Sequence.SFF ( SFF(..), CommonHeader(..)                         , ReadHeader(..), ReadBlock(..)                         , readSFF, writeSFF, writeSFF', recoverSFF-                        , sffToSequence, trim, trimFromTo, trimKey+                        , sffToSequence, rbToSequence+                        , trim, trimFromTo, trimKey                         , baseToFlowPos, flowToBasePos                         , test, convert, flowgram+                        , masked_bases, cumulative_index                         , packFlows, unpackFlows                         , Flow, Qual, Index, SeqData, QualData                         , ReadName (..), decodeReadName, encodeReadName@@ -63,9 +65,11 @@                              else Nothing -- error ("Couldn't match key in sequence "++LBC.unpack n++" ("++LBC.unpack k++" vs. "++BC.unpack (key ch)++")!")  sffToSequence :: SFF -> [Sequence Nuc]-sffToSequence (SFF ch rs) = map r2s rs-    where r2s r = clip (read_header r, bases r, quality r)-          clip (h, s, q) = let (left,right) = (clip_qual_left h,clip_qual_right h)+sffToSequence (SFF _ rs) = map rbToSequence rs++rbToSequence :: ReadBlock -> Sequence Nuc+rbToSequence r = clip (read_header r, bases r, quality r)+    where clip (h, s, q) = let (left,right) = (clip_qual_left h,clip_qual_right h)                                split x = let (a,b) = LB.splitAt (fromIntegral right) x                                               (c,d) = LB.splitAt (fromIntegral left-1) a                                          in [c,d,b]@@ -81,7 +85,13 @@ trimFromTo l r rd = let 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)-                    in rd { read_header = read_header rd -- FIXME: Update num_bases?+                        rh = read_header rd+                        [r',l'] = map fromIntegral [r,l]+                        rh' = rh { num_bases = fromIntegral (r'-l'+1)+                                 , 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 = trim_flw (flow_data rd)                           , flow_index = trim_seq' (flow_index rd)                           , bases = trim_seq (bases rd)@@ -203,7 +213,7 @@   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))+  replicateM_ (2*fl-B.length (flow_data rb)) (put (0::Word8))   putByteString (flow_index rb)   putLazyByteString (bases rb)   putLazyByteString (quality rb)@@ -327,6 +337,18 @@  flowgram :: ReadBlock -> [Flow] flowgram = unpackFlows . flow_data++masked_bases :: ReadBlock -> SeqData+masked_bases rb = let+  l = fromIntegral $ clip_qual_left $ read_header rb+  r = fromIntegral $ clip_qual_right $ read_header rb+  s = bases rb+  in 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]++cumulative_index :: ReadBlock -> [Int]+cumulative_index = scanl1 (+) . map fromIntegral . B.unpack . flow_index  instance Show ReadBlock where     show (ReadBlock h f i b q) =

@@ -282,7 +282,7 @@  -- | Return sequence label (first word of header) seqlabel :: Sequence a -> SeqData-seqlabel (Seq l _ _) = head (B.words l)+seqlabel (Seq l _ _) = case B.words l of (x:_) -> x; [] -> B.empty  -- | Return full header. seqheader :: Sequence a -> SeqData@@ -306,7 +306,7 @@ --   all but the sequence label (i.e. first word). appendHeader, setHeader :: Sequence a -> String -> Sequence a appendHeader (Seq h d q) t = (Seq (B.append h (B.pack (" "++t))) d q)-setHeader (Seq h d q) t = (Seq (B.append (head $ B.words h) (B.pack (" "++t))) d q)+setHeader s@(Seq _ d q) t = (Seq (B.append (seqlabel s) (B.pack (" "++t))) d q)  -- | Returns a sequence with all internal storage freshly copied and -- with sequence and quality data present as a single chunk.  

@@ -283,7 +283,7 @@   -- | Parse a (lazy) ByteString as sequences in the 2bit format.-decode2Bit :: B.ByteString -> [Sequence Unknown]+decode2Bit :: B.ByteString -> [Sequence Nuc] decode2Bit cs = let                      -- decode to (Header, [Entry]) from ByteString                     (Entries h es) = decode cs :: Entries@@ -422,17 +422,14 @@ unSRLE :: SRLE -> SR unSRLE (SRLE x) = x --- -- | Read sequences from a file in 2bit format and  -- | unmarshall/deserialize into Sequence format.-read2Bit  :: FilePath -> IO [Sequence Unknown]+read2Bit  :: FilePath -> IO [Sequence Nuc] read2Bit f = B.readFile f >>= return . decode2Bit  -- | Read sequences from a file handle in the 2bit format and -- | unmarshall/deserialze into Sequence format.-hRead2Bit :: Handle   -> IO [Sequence Unknown]+hRead2Bit :: Handle   -> IO [Sequence Nuc] hRead2Bit h = B.hGetContents h >>= return . decode2Bit  

@@ -1,9 +1,9 @@ Name:                bio-Version:             0.4.7+Version:             0.4.8 License:             LGPL License-file:        LICENSE Author:              Ketil Malde-Maintainer:          ketil@ii.uib.no+Maintainer:          ketil@malde.org  Category:            Bioinformatics Synopsis:            A bioinformatics library@@ -26,8 +26,8 @@  Tested-With:         GHC==6.12.1 Build-Type:          Simple-Build-Depends:       base>=4 && <5, QuickCheck>=2, binary >=0.4 && <0.5, tagsoup>=0.4 && <0.8, bytestring >= 0.9.1,-                     containers, array, parallel, parsec, random, old-time, mtl+Build-Depends:       base>=4 && <5, QuickCheck>=2, binary >=0.4 && <0.5, tagsoup>=0.8, bytestring >= 0.9.1,+                     containers, array, parallel, parsec, random, old-time, mtl, directory  Data-Files:          README