diff --git a/Bio/Alignment/BlastXML.hs b/Bio/Alignment/BlastXML.hs --- a/Bio/Alignment/BlastXML.hs +++ b/Bio/Alignment/BlastXML.hs @@ -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" diff --git a/Bio/Sequence.hs b/Bio/Sequence.hs --- a/Bio/Sequence.hs +++ b/Bio/Sequence.hs @@ -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] diff --git a/Bio/Sequence/FastQ.hs b/Bio/Sequence/FastQ.hs --- a/Bio/Sequence/FastQ.hs +++ b/Bio/Sequence/FastQ.hs @@ -5,26 +5,37 @@ * - 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: * * - Currently, we only support the non-Solexa FastQ, adding\/subtracting 33 for - the quality values. + * + 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!") diff --git a/Bio/Sequence/Fasta.hs b/Bio/Sequence/Fasta.hs --- a/Bio/Sequence/Fasta.hs +++ b/Bio/Sequence/Fasta.hs @@ -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) diff --git a/Bio/Sequence/SFF.hs b/Bio/Sequence/SFF.hs --- a/Bio/Sequence/SFF.hs +++ b/Bio/Sequence/SFF.hs @@ -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) = diff --git a/Bio/Sequence/SeqData.hs b/Bio/Sequence/SeqData.hs --- a/Bio/Sequence/SeqData.hs +++ b/Bio/Sequence/SeqData.hs @@ -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. diff --git a/Bio/Sequence/TwoBit.hs b/Bio/Sequence/TwoBit.hs --- a/Bio/Sequence/TwoBit.hs +++ b/Bio/Sequence/TwoBit.hs @@ -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 diff --git a/bio.cabal b/bio.cabal --- a/bio.cabal +++ b/bio.cabal @@ -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