bio-0.3.3: Bio/Alignment/ACE.hs
{- |
Read ACE format assembly files
These are typically output by sequence assembly tools,
like CAP3 or Phrap.
Documented in the section labelled \"ACE FILE FORMAT\" at
<http://bozeman.mbt.washington.edu/consed/distributions/README.14.0.txt>
Briefly: each field is a line starting with a two letter code,
in some cases followed by data lines termintated by a blank line.
-- AS contigs reads
-- CO contig_name bases reads segments compl (CAP3: segments=0)
-- sequence
-- BQ
-- base_qualities
-- AF read1 compl padded_start_consensus (negatives meaning?)
-- AF read2 ..
-- BS segments
-- RD read1 bases info_items info_tags (latter two set to 0 by CAP3)
-- sequence
-- QA read1 qual_start qual_end align_start align_end
-- DS (phred header? left empty by CAP3)
-- RD read2 ...
-}
{-# LANGUAGE CPP #-}
module Bio.Alignment.ACE (readACE, writeACE, Assembly(..)
,Gaps,extractGaps, insertGaps, ptest
)
where
#include <interlude.h>
import Interlude hiding (lines,words,readFile,unwords -- ByteString clashes
,reads -- Assembly clash
)
import Bio.Sequence.SeqData
import qualified Data.ByteString.Lazy.Char8 as B
import Data.ByteString.Lazy.Char8 (ByteString,words,pack,unpack,readFile,unwords)
import Text.ParserCombinators.Parsec
import Text.ParserCombinators.Parsec.Pos (newPos)
import Control.Monad (liftM) -- ,when?
import Data.Char (chr)
data Dir = Fwd | Rev deriving (Eq,Show)
data Assembly = Asm { contig :: (Sequence,Gaps), reads :: [(Offset,Dir,Sequence,Gaps)] }
deriving Show
type Gaps = [Offset]
type Str = ByteString
-- | ACE header lines with parameters
-- The tokenizer (scanner) should convert input into a list of these,
-- which in turn can be parsed by Parsec
data ACE = AS Str Str | CO Str Str Str Str Str
| BQ | AF Str Str Str | BS [Str]
| RD Str Str Str Str
| QA Str Str Str Str Str
| DS [Str] | Other [Str] | Empty
deriving (Eq)
instance Show ACE where
show (AS x y) = "AS "++uw [x,y]
show (CO a b c d e) = "CO "++uw [a,b,c,d,e]
show BQ = "BQ"
show (AF a b c) = "AF "++uw [a,b,c]
show (RD a b c d) = "RD "++uw [a,b,c,d]
show (QA a b c d e) = "QA "++uw [a,b,c,d,e]
show (DS ss) = "DS "++uw ss
show (Other ss) = uw ss
show Empty = "(blank)"
show _ = "unknown ACE string"
uw = unpack . unwords
-- | The Parsec parser type
type AceParser a = GenParser (SourcePos,ACE) () a
-- | Parse a single token, primitive parser
parse1 :: (ACE -> Maybe a) -> AceParser a
parse1 p = token sho pos pred
where sho (_,t) = show t
pos (n,_) = n
pred (_,t) = p t
-- | Test parser p on a list of ACE elements
ptest :: Show a => String -> AceParser a -> [ACE] -> IO ()
ptest m p = parseTest p . source m
-- | Add SourcePoses to a stream of ACEs.
source :: String -> [ACE] -> [(SourcePos,ACE)]
source m = zip (iterate (\sp -> incSourceLine sp 1) (newPos m 1 0))
-- | Parse a complete ACE file as a set of assemblies.
ace :: AceParser [[Assembly]]
ace = many1 ace1
ace1 :: AceParser [Assembly]
ace1 = do
as
many blank
many (do ctg >>= asm) -- apparently, CAP3 outputs empty assemblies (AS 0 0)
-- | parse the initial header
as :: AceParser (Int,Int)
as = parse1 (\t -> case t of AS cs rs -> do c <- B.readInt cs
r <- B.readInt rs
return (fst c,fst r)
_ -> Nothing) <?> "AS <int> <int>"
blank = parse1 (\t -> case t of Empty -> Just ()
_ -> Nothing) <?> "empty line"
-- | Gaps are coded as '*'s, this function removes them, and returns
-- the sequence along with the list of gap positions.
extractGaps :: Str -> (Str,Gaps)
extractGaps str = (B.filter (/='*') str,B.elemIndices '*' str)
-- todo: faster to lift concat out of the inner loop?
insertGaps :: Char -> (Str,Gaps) -> Str
insertGaps c (str,gaps) = go str B.empty 0 gaps
where go str acc p (next:rest) = let (a,b) = B.splitAt (next-p) str
in go b (B.concat [acc,a,pack [c]]) (next+1) rest
go str acc _ [] = B.append acc str
-- | parse the contig and quality information (CO, BQ)
ctg :: AceParser (Sequence,Gaps)
ctg = do
name <- co
sd <- sdata
let (sd',gaps) = extractGaps sd
many blank
bq
sq <- qdata
many blank
-- todo: gaps?
return (Seq name sd' (Just sq),gaps)
co, sdata, qdata :: AceParser Str
co = parse1 (\t -> case t of CO name a b c _compl -> do
_bs <- B.readInt a
_rds <- B.readInt b
_seg <- B.readInt c
return name
_ -> Nothing) <?> "CO name <int> <int> <int> bool"
sdata = do return . B.concat =<< many1 sdata1
where sdata1 = parse1 (\t -> case t of Other sd -> Just (unwords sd)
_ -> Nothing) <?> "sequence data"
qdata = do return . B.concat =<< many1 qdata1
where qdata1 = parse1 (\t -> case t of Other sd -> liftM (pack . map chr) (readInts sd)
_ -> Nothing) <?> "sequence data"
-- | Read a list of Ints in the Maybe monad
readInts :: [ByteString] -> Maybe [Int]
readInts [] = Just []
readInts (x:xs) = do (a,_) <- B.readInt x
as <- readInts xs
return $ (a:as)
bq :: AceParser ()
bq = parse1 (\t -> case t of BQ -> Just (); _ -> Nothing) <?> "BQ"
-- | Given the CO info, get the AFS'es
asm :: (Sequence,Gaps) -> AceParser Assembly
asm cg = do
many blank
afs <- many1 af
_bss <- many bs
many blank
rds cg afs
-- | Parse a list of AFS, followed by actual read, and merge them
-- afs :: Sequence -> AceParser [Sequence] -- plus some auxiliary info?
-- todo: better error handling!
af :: AceParser (Str,Dir,Offset)
af = parse1 (\t -> case t of AF a b c -> Just (a,f b,readInt' c)
_ -> Nothing) <?> "AF name compl pad_start"
where f b = case unpack b of "U" -> Fwd; "C" -> Rev
readInt' x = case B.readInt x of Just (a,_) -> fromIntegral a
bs :: AceParser (Int,Int,Str)
bs = parse1 (\t -> case t of BS [x,y,n] -> Just (readInt' x, readInt' y,n)
_ -> Nothing) <?> "BS x y name"
where readInt' i = case B.readInt i of Just (a,_) -> fromIntegral a
rds :: (Sequence,Gaps) -> [(Str,Dir,Offset)] -> AceParser Assembly
rds cg xs = do
r <- many1 rseq
-- todo: check the number and merge with the afs
let f (_name,d,off) (s,gs) = (off,d,s,gs)
return $ Asm { contig = cg, reads = zipWith f xs r }
rseq :: AceParser (Sequence,Gaps)
rseq = do
(rn,_len,_,_) <- rd
(s,gaps) <- return . extractGaps =<< sdata
-- when (B.length s == fromIntegral len) $ (fail "Incorrect sequence length!")
-- todo: fix gaps!
many1 blank
qa
ds
many blank
return (Seq rn s Nothing,gaps) -- huh?
-- | parse each read (RD, QA, DS)
rd :: AceParser (Str,Int,Int,Int)
rd = parse1 (\t -> case t of RD a b c d -> do [x,y,z] <- readInts [b,c,d]
return (a,x,y,z)
_ -> Nothing) <?> "RD <string> <int> <int> <int>"
qa :: AceParser ()
qa = parse1 (\_ -> Just ())
ds :: AceParser ()
ds = parse1 (\t -> case t of DS _ -> Just (); _ -> Nothing) <?> "DS"
-- ----------------------------------------------------------
-- Convert lines into tokens
tokenize :: ByteString -> [ACE]
tokenize = map tokenize1 . B.lines
-- Tokenise a single line
-- todo: error on incorrect (partial) format, error reports with line number
tokenize1 :: ByteString -> ACE
tokenize1 l | B.null l = Empty
| otherwise = let (h:ws) = words l
in case (unpack h,ws) of
("AS",[cs,rs]) -> AS cs rs
("CO",[nm,bs,rs,segs,compl]) -> CO nm bs rs segs compl
("BQ",[]) -> BQ
("AF",[a,b,c]) -> AF a b c
("BS",_) -> BS ws
("RD",[a,b,c,d]) -> RD a b c d
("QA",[a,b,c,d,e]) -> QA a b c d e
("DS",_) -> DS ws
_ -> Other (h:ws)
-- | Reading an ACE file.
readACE :: FilePath -> IO [[Assembly]]
readACE f = parseit =<< B.readFile f
where parseit = \s -> case (parse ace f . source f . tokenize) s of
Left e -> fail (show e)
Right a -> return a
-- formatError msg = error ("readACE: incorrect format in "++msg)
writeACE :: FilePath -> [Assembly] -> IO ()
writeACE = undefined
-- todo: hWrite etc