FastQ.hs

{- |
   Module: Bio.Sequence.FastQ

   Support the FastQ format that combines sequence and quality. See:

   * <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 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>

   * <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
    (
     -- * 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
import qualified Data.ByteString.Lazy.Char8 as B
import qualified Data.ByteString.Lazy as BB
import Data.List (unfoldr)

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

hReadFastQ :: Handle -> IO [Sequence Nuc]
hReadFastQ h = (go . B.lines) `fmap` B.hGetContents h

go :: [B.ByteString] -> [Sequence Nuc]
go = map (either error id) . unfoldr parse

-- | 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) =
    case (B.uncons h1,B.uncons h2) of
      -- The fast path: four-line format
      (Just ('@',h1name), Just ('+',h2name))
          | h1name == h2name || B.null h2name
            -> Just (Right $ Seq h1name sd (Just (BB.map (subtract 33) sq)), rest)
          | otherwise
            -> Just (Left $ "Bio.Sequence.FastQ: name mismatch:" ++ showStanza, rest)
      _     -> Just (Left $ "Bio.Sequence.FastQ: illegal FastQ format:" ++ showStanza, rest)
    where showStanza = unlines $ map B.unpack [ h1, sd, h2, sq ]
parse [] = Nothing
parse fs = let showStanza = unlines (map B.unpack fs)
               err = Left $ "Bio.Sequence.FastQ: illegal number of lines in FastQ format: " ++ showStanza
           in Just (err, [])

writeFastQ :: FilePath -> [Sequence Nuc] -> IO ()
writeFastQ f = B.writeFile f . B.concat . map unparse

hWriteFastQ :: Handle -> [Sequence Nuc] -> IO ()
hWriteFastQ h = B.hPut h . B.concat . map unparse

-- helper function for writing
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)
                                   ++" doesn not have quality data!")