diff --git a/LICENSE b/LICENSE --- a/LICENSE +++ b/LICENSE @@ -1,4 +1,4 @@ -Copyright (c)2011, Felipe Lessa +Copyright (c)2011-2012, Felipe Lessa All rights reserved. diff --git a/benchmarks/benchmark_v0.1.hs b/benchmarks/benchmark_v0.1.hs new file mode 100644 --- /dev/null +++ b/benchmarks/benchmark_v0.1.hs @@ -0,0 +1,8 @@ +import Bio.Sequence.Stockholm +import Control.Monad.Exception.Synchronous (Exceptional(..)) +import qualified Data.ByteString.Lazy as L + +main = L.getContents >>= mapM_ printName . parseStockholm + where + printName (Success (Stockholm file _ _)) = print (findAnn AC file) + printName (Exception ()) = putStr "---------------------" diff --git a/benchmarks/benchmark_v0.2.hs b/benchmarks/benchmark_v0.2.hs new file mode 100644 --- /dev/null +++ b/benchmarks/benchmark_v0.2.hs @@ -0,0 +1,19 @@ +{-# LANGUAGE OverloadedStrings #-} +import Bio.Sequence.Stockholm +import Bio.Sequence.Stockholm.Stream +import Control.Monad.Exception.Synchronous (Exceptional(..)) +import System.IO (stdin) +import qualified Data.Conduit as C +import qualified Data.Conduit.Binary as CB +import qualified Data.Conduit.List as CL + +main = main_event + +main_doc = C.runResourceT $ CB.sourceHandle stdin C.$$ parseStockholm C.=$ CL.mapM_ printName + where + printName (Stockholm file _ _) = print (findAnn AC file) + +main_event = C.runResourceT $ CB.sourceHandle stdin C.$$ parseEvents C.=$ CL.mapM_ printName + where + printName (EvGF "AC" name) = print name + printName _ = return () diff --git a/biostockholm.cabal b/biostockholm.cabal --- a/biostockholm.cabal +++ b/biostockholm.cabal @@ -1,6 +1,6 @@ Name: biostockholm -Version: 0.1.0.1 -Synopsis: Reading and writing Stockholm files (multiple sequence alignment, used by Rfam and Infernal). +Version: 0.2 +Synopsis: Parsing and rendering of Stockholm files (used by Pfam, Rfam and Infernal). License: BSD3 License-file: LICENSE Author: Felipe Lessa @@ -8,8 +8,15 @@ Category: Bioinformatics Build-type: Simple Cabal-version: >=1.8 +Extra-source-files: + benchmarks/benchmark_v0.1.hs + benchmarks/benchmark_v0.2.hs + tests/runtests.hs Description: - Parsing and pretty printing of files in Stockholm 1.0 format. See: + Parsing and rendering of files in Stockholm 1.0 format. Among + the users of the Stockholm format are Pfam, Rfam and Infernal. + These files hold information about families of proteins or + non-coding RNAs. For more information, please see: . * . @@ -25,32 +32,39 @@ Hs-Source-Dirs: src Exposed-modules: Bio.Sequence.Stockholm + Bio.Sequence.Stockholm.Document + Bio.Sequence.Stockholm.Stream Ghc-Options: -Wall Build-depends: base >= 3 && < 5 , containers >= 0.2 && < 0.5 , bytestring == 0.9.* - , deepseq == 1.1.* - , explicit-exception == 0.1.* + , deepseq >= 1.1 && < 1.3 + , conduit == 0.1.* + , attoparsec == 0.10.* + , attoparsec-conduit == 0.0.* + , blaze-builder == 0.3.* + , blaze-builder-conduit == 0.0.* , biocore >= 0.1 && < 0.3 Test-suite runtests Type: exitcode-stdio-1.0 - Hs-Source-Dirs: src + Hs-Source-Dirs: tests Main-is: runtests.hs - Other-modules: - Bio.Sequence.Stockholm Ghc-Options: -Wall - Cpp-Options: -DTEST Build-depends: base >= 3 && < 5 , containers >= 0.2 && < 0.5 , bytestring == 0.9.* - , deepseq == 1.1.* - , explicit-exception == 0.1.* + , conduit == 0.1.* + , zlib-conduit + , transformers >= 0.2 , biocore >= 0.1 && < 0.3 , hspec == 0.9.* - , HUnit == 1.2.* + , HUnit + , QuickCheck + + , biostockholm diff --git a/src/Bio/Sequence/Stockholm.hs b/src/Bio/Sequence/Stockholm.hs --- a/src/Bio/Sequence/Stockholm.hs +++ b/src/Bio/Sequence/Stockholm.hs @@ -10,565 +10,93 @@ -- - module Bio.Sequence.Stockholm - (-- * Data types - Stockholm(..) - ,StockholmSeq(..) - ,Ann(..) - ,FileAnnotation(..) - ,SequenceAnnotation(..) - ,ColumnAnnotation(..) - ,InFile - ,InSeq - ,findAnn + ( -- * Data types + Stockholm(..) + , StockholmSeq(..) + , Ann(..) + , FileAnnotation(..) + , SequenceAnnotation(..) + , ColumnAnnotation(..) + , InFile + , InSeq + , findAnn - -- * Parsing - ,parseStockholm - ,StockholmExc(..) + -- * Parsing + , parseStockholm - -- * Printing - ,prettyPrintStockholm + -- * Printing + , renderStockholm -#ifdef TEST - -- * Test cases - ,test_Stockholm -#endif + -- * Lazy I/O + , lazyParseStockholm + , lazyRenderStockholm ) where -- from base -import Control.Applicative ((<$>)) -import Control.Arrow (second) -import Control.DeepSeq (NFData(..)) -import Control.Monad (mplus) -import Data.Char (isSpace) import Data.List (find) -import Data.Maybe (fromMaybe) -import Data.Typeable (Typeable) - --- from containers -import qualified Data.Map as M +import System.IO.Unsafe (unsafePerformIO) -- from bytestring -import qualified Data.ByteString.Lazy.Char8 as B -import Data.ByteString.Lazy.Char8 (ByteString) - --- from biocore -import Bio.Core.Sequence - --- from explicit-exception -import Control.Monad.Exception.Synchronous (Exceptional, throw) - - -#ifdef TEST -import Test.Hspec.Monadic -import Test.Hspec.HUnit () -import Test.HUnit -#endif - - --- | An Stockholm 1.0 formatted file represented in memory. -data Stockholm = Stockholm [Ann FileAnnotation] - [Ann (ColumnAnnotation InFile)] - [StockholmSeq] - deriving (Show, Eq, Typeable) - -instance NFData Stockholm where - rnf (Stockholm file clmn seqs) = rnf file `seq` rnf clmn `seq` rnf seqs - --- | A sequence in Stockholm 1.0 format. -data StockholmSeq = StSeq !SeqLabel - !SeqData - [Ann SequenceAnnotation] - [Ann (ColumnAnnotation InSeq)] - deriving (Eq, Typeable) - --- We don't derive Show to be able support biocore-0.1, which --- doesn't have Show instances for SeqLabel and SeqData. -instance Show StockholmSeq where - showsPrec prec (StSeq (SeqLabel l) (SeqData d) sa ca) = - showParen (prec > 10) $ - showString "StSeq (SeqLabel " . - showsPrec 11 l . - showString ") (SeqData " . - showsPrec 11 d . - (')':) . (' ':) . - showsPrec 11 sa . - (' ':) . - showsPrec 11 ca - - -instance NFData StockholmSeq where - rnf (StSeq _ _ sa ca) = rnf sa `seq` rnf ca - -instance BioSeq StockholmSeq where - seqlabel (StSeq sl _ _ _) = sl - seqdata (StSeq _ sd _ _) = sd - seqlength (StSeq _ sd _ _) = Offset $ B.length (unSD sd) - --- | A generic annotation. -data Ann d = Ann { feature :: !d - , text :: !ByteString - } - deriving (Show, Eq, Ord, Typeable) - -instance NFData (Ann d) where - -- already strict, default instance - --- | Possible file annotations. -data FileAnnotation = - AC -- ^ Accession number: Accession number in form PFxxxxx.version or PBxxxxxx. - | ID -- ^ Identification: One word name for family. - | DE -- ^ Definition: Short description of family. - | AU -- ^ Author: Authors of the entry. - | SE -- ^ Source of seed: The source suggesting the seed members belong to one family. - | GA -- ^ Gathering method: Search threshold to build the full alignment. - | TC -- ^ Trusted Cutoff: Lowest sequence score and domain score of match in the full alignment. - | NC -- ^ Noise Cutoff: Highest sequence score and domain score of match not in full alignment. - | TP -- ^ Type: Type of family (presently Family, Domain, Motif or Repeat). - | SQ -- ^ Sequence: Number of sequences in alignment. - | AM -- ^ Alignment Method: The order ls and fs hits are aligned to the model to build the full align. - - | DC -- ^ Database Comment: Comment about database reference. - | DR -- ^ Database Reference: Reference to external database. - | RC -- ^ Reference Comment: Comment about literature reference. - | RN -- ^ Reference Number: Reference Number. - | RM -- ^ Reference Medline: Eight digit medline UI number. - | RT -- ^ Reference Title: Reference Title. - | RA -- ^ Reference Author: Reference Author - | RL -- ^ Reference Location: Journal location. - | PI -- ^ Previous identifier: Record of all previous ID lines. - | KW -- ^ Keywords: Keywords. - | CC -- ^ Comment: Comments. - | NE -- ^ Pfam accession: Indicates a nested domain. - | NL -- ^ Location: Location of nested domains - sequence ID, start and end of insert. - - | F_Other !ByteString -- ^ Other file annotation. - deriving (Show, Eq, Ord, Typeable) - --- | Possible column annotations. Phantom type can be 'InFile' --- or 'InSeq'. -data ColumnAnnotation a = - SS -- ^ Secondary structure. - | SA -- ^ Surface accessibility. - | TM -- ^ TransMembrane. - | PP -- ^ Posterior probability. - | LI -- ^ LIgand binding. - | AS -- ^ Active site. - | PAS -- ^ AS - Pfam predicted. - | SAS -- ^ AS - from SwissProt. - | IN -- ^ INtron (in or after). - - | C_Other !ByteString -- ^ Other column annotation. - deriving (Show, Eq, Ord, Typeable) - --- | Phantom type for 'ColumnAnnotation's of the whole file. -data InFile --- | Phantom type for 'ColumnAnnotation's of a single sequence. -data InSeq - --- | Possible sequence annotations. -data SequenceAnnotation = - S_AC -- ^ Accession number - | S_DE -- ^ Description - | S_DR -- ^ Database reference - | OS -- ^ Organism (species) - | OC -- ^ Organism classification (clade, etc.) - | LO -- ^ Look (Color, etc.) - - | S_Other !ByteString -- ^ Other sequence annotation. - deriving (Show, Eq, Ord, Typeable) - - --- | Class used internally below just to simplify the code. -class IsAnnotation a where - parseAnn :: ByteString -> a - showAnn :: a -> ByteString - -mkParseAnn :: (ByteString -> a -> ByteString) -> [(ByteString, a)] - -> (ByteString -> a) -> ByteString -> a -mkParseAnn modify anns mkOther = - let annots = M.fromList anns - in \feat -> let featMod = modify feat (error "mkParseAnn: never here") - in fromMaybe (mkOther feat) $ M.lookup featMod annots - -mkShowAnn :: Ord a => (ByteString -> a -> ByteString) -> [(ByteString, a)] - -> (a -> Maybe ByteString) -> a -> ByteString -mkShowAnn modify anns fromOther = - let annots = M.fromList [(a,b) | (b,a) <- anns] - in \ann -> fromMaybe (error "mkShowAnn: never here 2") $ - fromOther ann `mplus` (mod' <$> M.lookup ann annots) - where mod' = flip modify (error "mkShowAnn: never here 1") - -instance IsAnnotation SequenceAnnotation where - parseAnn = mkParseAnn const seqAnns S_Other - showAnn = mkShowAnn const seqAnns f - where f (S_Other o) = Just o - f _ = Nothing - -seqAnns :: [(ByteString, SequenceAnnotation)] -seqAnns = [("LO",LO), ("OC",OC), ("OS",OS), - ("AC",S_AC), ("DE",S_DE), ("DR",S_DR)] - - -instance IsAnnotation FileAnnotation where - parseAnn = mkParseAnn const fileAnns F_Other - showAnn = mkShowAnn const fileAnns f - where f (F_Other o) = Just o - f _ = Nothing - -fileAnns :: [(ByteString, FileAnnotation)] -fileAnns = [("AC",AC), ("AM",AM), ("AU",AU), ("CC",CC), - ("DC",DC), ("DE",DE), ("DR",DR), ("GA",GA), - ("ID",ID), ("KW",KW), ("NC",NC), ("NE",NE), - ("NL",NL), ("PI",PI), ("RA",RA), ("RC",RC), - ("RL",RL), ("RM",RM), ("RN",RN), ("RT",RT), - ("SE",SE), ("SQ",SQ), ("TC",TC), ("TP",TP)] - - -instance ClmnAnnLoc a => IsAnnotation (ColumnAnnotation a) where - parseAnn = mkParseAnn removeSuffix clmnAnns C_Other - where - removeSuffix feat phantom = - let suffix = clmnAnnSuffix phantom - (f, s) = B.splitAt (B.length feat - B.length suffix) feat - in if suffix == s then f else "" - showAnn = mkShowAnn addSuffix clmnAnns f - where - f (C_Other o) = Just o - f _ = Nothing - addSuffix feat phantom = feat `B.append` clmnAnnSuffix phantom - -clmnAnns :: [(ByteString, ColumnAnnotation a)] -clmnAnns = [("AS",AS), ("IN",IN), ("LI",LI), ("PAS",PAS), ("PP",PP), - ("SA",SA), ("SAS",SAS), ("SS",SS), ("TM",TM)] - -class ClmnAnnLoc a where - clmnAnnSuffix :: b a -> ByteString -instance ClmnAnnLoc InSeq where - clmnAnnSuffix _ = "" -instance ClmnAnnLoc InFile where - clmnAnnSuffix _ = "_cons" - -parseAnn' :: IsAnnotation a => ByteString -> ByteString -> Ann a -parseAnn' = Ann . parseAnn - - +import qualified Data.ByteString as B +import qualified Data.ByteString.Lazy as L +-- from conduit +import qualified Data.Conduit as C +import qualified Data.Conduit.Lazy as CZ +import qualified Data.Conduit.List as CL -#ifdef TEST -test_parseAnnots :: Specs -test_parseAnnots = - describe "parse*" $ do - it "1" $ parseAnn "AC" @?= AC - it "2" $ parseAnn "SQ" @?= SQ - it "3a" $ parseAnn "SS" @?= (SS :: ColumnAnnotation InSeq) - it "3b" $ parseAnn "SS" @?= (C_Other "SS" :: ColumnAnnotation InFile) - it "4a" $ parseAnn "SS_cons" @?= (SS :: ColumnAnnotation InFile) - it "4b" $ parseAnn "SS_cons" @?= (C_Other "SS_cons" :: ColumnAnnotation InSeq) - it "5a" $ parseAnn "SS_CONS" @?= (C_Other "SS_CONS" :: ColumnAnnotation InSeq) - it "5b" $ parseAnn "SS_CONS" @?= (C_Other "SS_CONS" :: ColumnAnnotation InFile) - it "6" $ parseAnn "LO" @?= LO -#endif +-- from this package +import Bio.Sequence.Stockholm.Stream +import Bio.Sequence.Stockholm.Document -- | Find an annotation. For example, you may use @'findAnn' 'SS'@ --- to find the secondary of an Stockholm file. -findAnn :: Eq d => d -> [Ann d] -> Maybe ByteString +-- to find the secondary on an Stockholm file. +findAnn :: Eq d => d -> [Ann d] -> Maybe L.ByteString findAnn x = fmap text . find ((== x) . feature) --- | Exceptions that may happen while parsing a Stockholm file. -class StockholmExc e where - -- | File is empty. - emptyFileExc :: e - - -- | Header is missing. - headerExc :: e - - -- | Malformed annotation. The line is passed as argument. - malformedAnnExc :: ByteString -> e - - -- | Unknown annotation type. - unknownAnnTypeExc :: Char -> e - - -- | Malformed sequence line. The line is passed as argument. - malformedSeqDataExc :: ByteString -> e - -instance StockholmExc () where - emptyFileExc = () - headerExc = () - malformedAnnExc _ = () - unknownAnnTypeExc _ = () - malformedSeqDataExc _ = () - -instance StockholmExc ByteString where - emptyFileExc = "parseStockholm: empty file." - headerExc = "parseStockholm: header is missing." - malformedAnnExc line = - B.concat ["parseStockholm: malformed annotation '", line, "'."] - unknownAnnTypeExc typ = - B.concat ["parseStockholm: unknown annotation type '", B.pack [typ], "'."] - malformedSeqDataExc line = - B.concat ["parseStockholm: malformed sequence data line '", line, "'."] - - - - - - - - - - - - - - - - --- | Some kind of annotation -data ParseAnnRet = - FileAnn !(Ann FileAnnotation) - | FileColAnn !(Ann (ColumnAnnotation InFile)) - | SeqAnn !ByteString !(Ann SequenceAnnotation) - | SeqColAnn !ByteString !(Ann (ColumnAnnotation InSeq)) - --- | @parseAnnotation line@ tries to parse a line as an annotation. -parseAnnotation :: (StockholmExc e) => ByteString -> Exceptional e ParseAnnRet -parseAnnotation line - | not (B.isPrefixOf "#=G" line) || B.length line < 5 = - throw (malformedAnnExc line) -parseAnnotation line = - let Just (typ, rest) = B.uncons $ B.drop 3 line - (word1, text1) = second dropSpace . B.break isSpace $ dropSpace rest - (word2, text2) = second dropSpace . B.break isSpace $ text1 - in case typ of - 'F' -> return $ FileAnn (parseAnn' word1 text1) - 'C' -> return $ FileColAnn (parseAnn' word1 text1) - 'S' -> return $ SeqAnn word1 (parseAnn' word2 text2) - 'R' -> return $ SeqColAnn word1 (parseAnn' word2 text2) - _ -> throw (unknownAnnTypeExc typ) - -dropSpace :: ByteString -> ByteString -dropSpace = B.dropWhile isSpace - - --- | @parseSeqData line@ tries to parse a line as some data of a sequence. -parseSeqData :: (StockholmExc e) => ByteString - -> Exceptional e (SeqLabel, SeqData) -parseSeqData str = case B.words str of - [ident, sq] -> return (SeqLabel ident, SeqData sq) - _ -> throw (malformedSeqDataExc str) - - --- | @parseStockholm@ parses a file in Stockholm 1.0 format. --- --- Each file must be completely read before it is used because --- the Stockholm format allows information to be given in any --- part of the file. However, there may be multiple --- \"Stockholm files\" concatenated in a single \"filesystem --- file\". These multiple files are read independently, which --- is why we return a list of 'Exceptional'@s@. +-- | @parseStockholm@ parses a stream of files in Stockholm 1.0 +-- format. -- --- If you prefer to read the whole file in one go, use --- @'sequence' (parseStockholm input)@, which will fail if any --- family fails. -parseStockholm :: (StockholmExc e) => ByteString - -> [Exceptional e Stockholm] -parseStockholm = map parseStockholm' . split . - filter (not . B.all isSpace) . B.lines - where - split [] = [] - split xs = let ~(y, ys) = break (B.isPrefixOf "//") xs - in y : split (tail ys) - - -parseStockholm' :: (StockholmExc e) => [ByteString] - -> Exceptional e Stockholm -parseStockholm' = header . filter (not . B.null) - where - -- Find - header (h:hs) - | h == stockholm = do (annots, seqs) <- go (emptyPA, M.empty) hs - return (makeStockholm annots seqs) - | otherwise = throw headerExc - where stockholm = "# STOCKHOLM 1.0" - header [] = throw emptyFileExc - - -- End of file - go acc [] = return acc - - -- Annotation - go (!annots, !seqs) (line:ls) | B.take 2 line == "#=" = do - annot <- parseAnnotation line - go (insertPA annot annots, seqs) ls - - -- Comment - go acc (l:ls) | B.head l == '#' = go acc ls - - -- Otherwise a sequence - go (!annots, !seqs) (line:ls) = do - seqData <- parseSeqData line - go (annots, insertDM seqData seqs) ls - - -type DiffMap a b = M.Map a [b] - -insertDM :: Ord a => (a, b) -> DiffMap a b -> DiffMap a b -insertDM (key, val) = M.insertWith' (\_ old -> val:old) key [val] - -finishDM :: (b -> ByteString) -> DiffMap a b -> M.Map a ByteString -finishDM f = fmap (B.concat . map f . reverse) - -type AnnMap d = DiffMap d ByteString - -insertAnn :: Ord d => Ann d -> AnnMap d -> AnnMap d -insertAnn (Ann key val) = insertDM (key, val) - -finishAnn :: AnnMap d -> [Ann d] -finishAnn m = [Ann a b | (a, b) <- M.toList (finishDM id m)] - -type SeqAnnMap d = M.Map ByteString (AnnMap d) - -insertSM :: Ord d => ByteString -> Ann d -> SeqAnnMap d -> SeqAnnMap d -insertSM sq ann = M.alter (just . insertAnn ann . fromMaybe M.empty) sq - where - just !x = Just x - -finishSM :: SeqAnnMap d -> M.Map ByteString [Ann d] -finishSM = fmap finishAnn - -data PartialAnns = - PartialAnns { paFileAnns :: !(AnnMap FileAnnotation) - , paFileColAnns :: !(AnnMap (ColumnAnnotation InFile)) - , paSeqAnns :: !(SeqAnnMap SequenceAnnotation) - , paSeqColAnns :: !(SeqAnnMap (ColumnAnnotation InSeq)) - } - -emptyPA :: PartialAnns -emptyPA = PartialAnns M.empty M.empty M.empty M.empty - -insertPA :: ParseAnnRet -> PartialAnns -> PartialAnns -insertPA (FileAnn ann) pa = pa { paFileAnns = insertAnn ann (paFileAnns pa) } -insertPA (FileColAnn ann) pa = pa { paFileColAnns = insertAnn ann (paFileColAnns pa) } -insertPA (SeqAnn sq ann) pa = pa { paSeqAnns = insertSM sq ann (paSeqAnns pa) } -insertPA (SeqColAnn sq ann) pa = pa { paSeqColAnns = insertSM sq ann (paSeqColAnns pa) } - - - --- | Glue everything in place, as the Stockholm format lets --- everything be everywhere and split in any number of parts. -makeStockholm :: PartialAnns -> DiffMap SeqLabel SeqData -> Stockholm -makeStockholm annots seqsDM = - let fileAnns_ = finishAnn (paFileAnns annots) - fileColAnns = finishAnn (paFileColAnns annots) - seqAnns_ = finishSM (paSeqAnns annots) - seqColAnns = finishSM (paSeqColAnns annots) - - stseqs = [StSeq label (SeqData dt) (f sq seqAnns_) (f sq seqColAnns) - | (label@(SeqLabel sq), dt) <- M.toList (finishDM unSD seqsDM)] - where - f = M.findWithDefault [] - in Stockholm fileAnns_ fileColAnns stseqs - +-- Each file must be completely read before it is used because +-- the Stockholm format allows information to be given in any +-- part of the file. However, there may be multiple \"Stockholm +-- files\" concatenated in a single \"filesystem file\". These +-- multiple files are read independently. If you need to process +-- large Stockholm files, consider using the streaming interface +-- on "Bio.Sequence.Stockholm.Stream". +parseStockholm :: C.ResourceThrow m => C.Conduit B.ByteString m Stockholm +parseStockholm = parseEvents C.=$= parseDoc --- | Pretty-prints an Stockholm file. We follow Rfam preferences --- and do not wrap lines. -prettyPrintStockholm :: Stockholm -> B.ByteString -prettyPrintStockholm (Stockholm file clmn seqs) = - let showAnnF :: IsAnnotation a => Char -> Ann a -> (ByteString, ByteString) - showAnnF t ann = (B.concat [B.pack ("#=G" ++ t : " "), - showAnn (feature ann)], text ann) - showAnnS :: IsAnnotation a => ByteString -> Char -> Ann a -> (ByteString, ByteString) - showAnnS s t ann = (B.unwords [B.pack ("#=G" ++ [t]), s, - showAnn (feature ann)], text ann) - - fileLines = map (showAnnF 'F') file - clmnLines = map (showAnnF 'C') clmn - sequences = do - StSeq (SeqLabel name) (SeqData seqd) sa ca <- seqs - (name, seqd) : map (showAnnS name 'R') ca - ++ map (showAnnS name 'S') sa - - allLines = fileLines ++ sequences ++ clmnLines - firstColLen = maximum $ map (B.length . fst) allLines - mkLine (col1, col2) = B.concat [col1, B.replicate n ' ', col2] - where n = 1 + firstColLen - B.length col1 - in B.unlines ("# STOCKHOLM 1.0" : map mkLine allLines ++ ["//"]) - +-- | Pretty prints an Stockholm file. +renderStockholm :: C.ResourceUnsafeIO m => C.Conduit Stockholm m B.ByteString +renderStockholm = renderDoc C.=$= renderEvents -#ifdef TEST -stockFile :: B.ByteString -stockFile = B.unlines [ - "# STOCKHOLM 1.0", - "#=GF AU Infernal 1.0", - "", - "#=GS Purine1 DE Number 1 :)", - "Purine1 AAAAUUGAAUAUCGUUUUACUUGUUUAUGUC-GUGAAU-UGGCAC-GACG", - "Purine2 AAAAUUUAAUAA-GAAGCACUCAUAUAAUCCCGAGAAUAUGGCUCGGGAG", - "Purine3 UGGCAGUAACUAGCGUCACUUCGUAUAACCCCAGUGAUAUGGAUUGGGGG", - "#=GC SS_cons :::::::::::::::::((((((((,,,<<<-<<<_______>>>->>>,", - "", - "# We may have comments =)", - "", - "Purine1 UUUCUACAAGGUG-CCGGAA--CACCUAACAAUAAGUAAGUCAGCAGUGA", - "Purine2 UCUCUACCGAACAACCGUAAAUUGUUCGACUAUGAGUGAAAGUGUACCUA", - "Purine3 UCUCUACCAGGAACCAAUAA--AUCCUGAUUACGAAGAGUUUAGUGCUUU", - "#=GC SS_cons ,,,,,,,<<<<<<_________>>>>>>,,))))))))::::::::::::", - "", - "Purine1 GAU", - "Purine2 GGG", - "Purine3 AGU", - "#=GC SS_cons :::", - "// "] - -purine1, purine2, purine3 :: SeqData -ss_cons :: ByteString -purine1 = SeqData "AAAAUUGAAUAUCGUUUUACUUGUUUAUGUC-GUGAAU-UGGCAC-GACGUUUCUACAAGGUG-CCGGAA--CACCUAACAAUAAGUAAGUCAGCAGUGAGAU" -purine2 = SeqData "AAAAUUUAAUAA-GAAGCACUCAUAUAAUCCCGAGAAUAUGGCUCGGGAGUCUCUACCGAACAACCGUAAAUUGUUCGACUAUGAGUGAAAGUGUACCUAGGG" -purine3 = SeqData "UGGCAGUAACUAGCGUCACUUCGUAUAACCCCAGUGAUAUGGAUUGGGGGUCUCUACCAGGAACCAAUAA--AUCCUGAUUACGAAGAGUUUAGUGCUUUAGU" -ss_cons = ":::::::::::::::::((((((((,,,<<<-<<<_______>>>->>>,,,,,,,,<<<<<<_________>>>>>>,,)))))))):::::::::::::::" - -result :: [Stockholm] -result = [Stockholm file clmn seqs] - where - file = [Ann AU "Infernal 1.0"] - clmn = [Ann SS ss_cons] - seqs = [mkStock "Purine1" purine1 [Ann S_DE "Number 1 :)"], - mkStock "Purine2" purine2 [], - mkStock "Purine3" purine3 []] - mkStock name data_ sa = StSeq name data_ sa [] - -stockFile2 :: B.ByteString -stockFile2 = B.unlines [stockFile, stockFile] - -result2 :: [Stockholm] -result2 = result ++ result - -returnExc :: [a] -> [Exceptional B.ByteString a] -returnExc = map return - -test_parseStockholm :: Specs -test_parseStockholm = - describe "parseStockholm" $ do - it "correctly parses test file 1" $ parseStockholm stockFile @?= returnExc result - it "correctly parses test file 2" $ parseStockholm stockFile2 @?= returnExc result2 - -test_prettyPrintStockholm :: Specs -test_prettyPrintStockholm = - describe "parseStockholm/prettyPrintStockholm" $ do - it "parses printed test file 1" $ parseStockholm (func result) @?= returnExc result - it "parses printed test file 2" $ parseStockholm (func result2) @?= returnExc result2 - where func = B.unlines . map prettyPrintStockholm -#endif +-- | Use lazy I/O to parse a stream of files in Stockholm 1.0 +-- format. We recommend using 'parseStockholm'. +lazyParseStockholm :: L.ByteString -> [Stockholm] +lazyParseStockholm lbs = + unsafePerformIO $ + C.runResourceT $ + CZ.lazyConsume $ + CL.sourceList (L.toChunks lbs) C.$= + parseStockholm +{-# NOINLINE lazyParseStockholm #-} -#ifdef TEST -test_Stockholm :: Specs -test_Stockholm = describe "Bio.Sequence.Stockholm" $ do - test_parseAnnots - test_parseStockholm - test_prettyPrintStockholm -#endif +-- | Use lazy I/O to render a list of 'Stockholm'@s@ into a +-- stream of files in Stockholm 1.0 format. We recommend using +-- 'renderStockholm'. +lazyRenderStockholm :: [Stockholm] -> L.ByteString +lazyRenderStockholm stos = + L.fromChunks $ + unsafePerformIO $ + C.runResourceT $ + CZ.lazyConsume $ + CL.sourceList stos C.$= + renderStockholm +{-# NOINLINE lazyRenderStockholm #-} diff --git a/src/Bio/Sequence/Stockholm/Document.hs b/src/Bio/Sequence/Stockholm/Document.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Sequence/Stockholm/Document.hs @@ -0,0 +1,427 @@ +{-# LANGUAGE BangPatterns, CPP, EmptyDataDecls, DeriveDataTypeable, OverloadedStrings #-} +-- | Take low-level 'Event'@s@ and turn them high-level data +-- structures. +module Bio.Sequence.Stockholm.Document + ( -- * Data types + Stockholm(..) + , StockholmSeq(..) + , Ann(..) + , FileAnnotation(..) + , SequenceAnnotation(..) + , ColumnAnnotation(..) + , InFile + , InSeq + + -- * Conduits + , parseDoc + , renderDoc + ) + where + +-- from base +import Control.Applicative ((<$>)) +import Control.DeepSeq (NFData(..)) +import Control.Monad (mplus) +import Data.Maybe (fromMaybe) +import Data.Typeable (Typeable) + +-- from containers +import qualified Data.Map as M + +-- from bytestring +import qualified Data.ByteString.Char8 as B +import qualified Data.ByteString.Lazy.Char8 as L + +-- from biocore +import Bio.Core.Sequence + +-- from conduit +import qualified Data.Conduit as C +import qualified Data.Conduit.List as CL + +-- from this package +import Bio.Sequence.Stockholm.Stream + + + +---------------------------------------------------------------------- +-- Types + +-- | An Stockholm 1.0 formatted file represented in memory. +data Stockholm = Stockholm [Ann FileAnnotation] + [Ann (ColumnAnnotation InFile)] + [StockholmSeq] + deriving (Show, Eq, Ord, Typeable) + +instance NFData Stockholm where + rnf (Stockholm file clmn seqs) = rnf file `seq` rnf clmn `seq` rnf seqs + + +-- | A sequence in Stockholm 1.0 format. +data StockholmSeq = StSeq !SeqLabel + !SeqData + [Ann SequenceAnnotation] + [Ann (ColumnAnnotation InSeq)] + deriving (Eq, Ord, Typeable) + +-- We don't derive Show to be able support biocore-0.1, which +-- doesn't have Show instances for SeqLabel and SeqData. +instance Show StockholmSeq where + showsPrec prec (StSeq (SeqLabel l) (SeqData d) sa ca) = + showParen (prec > 10) $ + showString "StSeq (SeqLabel " . + showsPrec 11 l . + showString ") (SeqData " . + showsPrec 11 d . + (')':) . (' ':) . + showsPrec 11 sa . + (' ':) . + showsPrec 11 ca + +instance NFData StockholmSeq where + rnf (StSeq _ _ sa ca) = rnf sa `seq` rnf ca + +instance BioSeq StockholmSeq where + seqlabel (StSeq sl _ _ _) = sl + seqdata (StSeq _ sd _ _) = sd + seqlength (StSeq _ sd _ _) = Offset $ L.length (unSD sd) + + +-- | A generic annotation. +data Ann d = Ann { feature :: !d + , text :: !L.ByteString + } + deriving (Show, Eq, Ord, Typeable) + +instance NFData (Ann d) where + -- already strict, default instance + + +-- | Possible file annotations. +data FileAnnotation = + AC -- ^ Accession number: Accession number in form PFxxxxx.version or PBxxxxxx. + | ID -- ^ Identification: One word name for family. + | DE -- ^ Definition: Short description of family. + | AU -- ^ Author: Authors of the entry. + | SE -- ^ Source of seed: The source suggesting the seed members belong to one family. + | GA -- ^ Gathering method: Search threshold to build the full alignment. + | TC -- ^ Trusted Cutoff: Lowest sequence score and domain score of match in the full alignment. + | NC -- ^ Noise Cutoff: Highest sequence score and domain score of match not in full alignment. + | TP -- ^ Type: Type of family (presently Family, Domain, Motif or Repeat). + | SQ -- ^ Sequence: Number of sequences in alignment. + | AM -- ^ Alignment Method: The order ls and fs hits are aligned to the model to build the full align. + + | DC -- ^ Database Comment: Comment about database reference. + | DR -- ^ Database Reference: Reference to external database. + | RC -- ^ Reference Comment: Comment about literature reference. + | RN -- ^ Reference Number: Reference Number. + | RM -- ^ Reference Medline: Eight digit medline UI number. + | RT -- ^ Reference Title: Reference Title. + | RA -- ^ Reference Author: Reference Author + | RL -- ^ Reference Location: Journal location. + | PI -- ^ Previous identifier: Record of all previous ID lines. + | KW -- ^ Keywords: Keywords. + | CC -- ^ Comment: Comments. + | NE -- ^ Pfam accession: Indicates a nested domain. + | NL -- ^ Location: Location of nested domains - sequence ID, start and end of insert. + + | F_Other !B.ByteString -- ^ Other file annotation. + deriving (Show, Eq, Ord, Typeable) + + +-- | Possible column annotations. Phantom type can be 'InFile' +-- or 'InSeq'. +data ColumnAnnotation a = + SS -- ^ Secondary structure. + | SA -- ^ Surface accessibility. + | TM -- ^ TransMembrane. + | PP -- ^ Posterior probability. + | LI -- ^ LIgand binding. + | AS -- ^ Active site. + | PAS -- ^ AS - Pfam predicted. + | SAS -- ^ AS - from SwissProt. + | IN -- ^ INtron (in or after). + + | C_Other !B.ByteString -- ^ Other column annotation. + deriving (Show, Eq, Ord, Typeable) + + +-- | Phantom type for 'ColumnAnnotation's of the whole file. +data InFile + + +-- | Phantom type for 'ColumnAnnotation's of a single sequence. +data InSeq + + +-- | Possible sequence annotations. +data SequenceAnnotation = + S_AC -- ^ Accession number + | S_DE -- ^ Description + | S_DR -- ^ Database reference + | OS -- ^ Organism (species) + | OC -- ^ Organism classification (clade, etc.) + | LO -- ^ Look (Color, etc.) + + | S_Other !B.ByteString -- ^ Other sequence annotation. + deriving (Show, Eq, Ord, Typeable) + + + + +---------------------------------------------------------------------- +-- Parsing and showing features + + +-- | Parse a feature. +type ParseFeature a = B.ByteString -> a + + +-- | Show a feature. +type ShowFeature a = a -> B.ByteString + + +-- | Helper to create 'ParseFeature'@s@. +mkParseFeature :: (B.ByteString -> a -> B.ByteString) + -> [(B.ByteString, a)] + -> (B.ByteString -> a) + -> ParseFeature a +mkParseFeature modify anns mkOther = + let annots = M.fromList anns + in \feat -> let featMod = modify feat (error "mkParseFeature: never here") + in fromMaybe (mkOther feat) $ M.lookup featMod annots + + +-- | Helper to create 'ShowFeature'@s@. +mkShowFeature :: Ord a => + (B.ByteString -> a -> B.ByteString) + -> [(B.ByteString, a)] + -> (a -> Maybe B.ByteString) + -> ShowFeature a +mkShowFeature modify anns fromOther = + let annots = M.fromList [(a,b) | (b,a) <- anns] + in \ann -> fromMaybe (error "mkShowFeature: never here 2") $ + fromOther ann `mplus` (mod' <$> M.lookup ann annots) + where mod' = flip modify (error "mkShowFeature: never here 1") + + +-- | Parse and show sequence annotations. +parseSeqFeature :: ParseFeature SequenceAnnotation +showSeqFeature :: ShowFeature SequenceAnnotation +(parseSeqFeature, showSeqFeature) = + ( mkParseFeature const seqFeatures S_Other + , mkShowFeature const seqFeatures f ) + where + f (S_Other o) = Just o + f _ = Nothing + + seqFeatures = [("LO",LO), ("OC",OC), ("OS",OS), + ("AC",S_AC), ("DE",S_DE), ("DR",S_DR)] + + +-- | Parse and show file annotations. +parseFileFeature :: ParseFeature FileAnnotation +showFileFeature :: ShowFeature FileAnnotation +(parseFileFeature, showFileFeature) = + ( mkParseFeature const fileFeatures F_Other + , mkShowFeature const fileFeatures f ) + where + f (F_Other o) = Just o + f _ = Nothing + + fileFeatures = [("AC",AC), ("AM",AM), ("AU",AU), ("CC",CC), + ("DC",DC), ("DE",DE), ("DR",DR), ("GA",GA), + ("ID",ID), ("KW",KW), ("NC",NC), ("NE",NE), + ("NL",NL), ("PI",PI), ("RA",RA), ("RC",RC), + ("RL",RL), ("RM",RM), ("RN",RN), ("RT",RT), + ("SE",SE), ("SQ",SQ), ("TC",TC), ("TP",TP)] + + +-- | Parse and show column annotations. +parseClmnFeature :: ClmnFeatureLoc a => ParseFeature (ColumnAnnotation a) +parseClmnFeature = mkParseFeature removeSuffix clmnFeatures C_Other + where + removeSuffix feat phantom = + let suffix = clmnFeatureSuffix phantom + (f, s) = B.splitAt (B.length feat - B.length suffix) feat + in if suffix == s then f else "" + +showClmnFeature :: ClmnFeatureLoc a => ShowFeature (ColumnAnnotation a) +showClmnFeature = mkShowFeature addSuffix clmnFeatures f + where + f (C_Other o) = Just o + f _ = Nothing + + addSuffix feat phantom = feat `B.append` clmnFeatureSuffix phantom + +clmnFeatures :: [(B.ByteString, ColumnAnnotation a)] +clmnFeatures = [("AS",AS), ("IN",IN), ("LI",LI), ("PAS",PAS), ("PP",PP), + ("SA",SA), ("SAS",SAS), ("SS",SS), ("TM",TM)] + +class ClmnFeatureLoc a where + clmnFeatureSuffix :: b a -> B.ByteString +instance ClmnFeatureLoc InSeq where + clmnFeatureSuffix _ = "" +instance ClmnFeatureLoc InFile where + clmnFeatureSuffix _ = "_cons" + + + +---------------------------------------------------------------------- +-- Specilized Maps. + +type DiffMap a b = M.Map a [b] + +insertDM :: Ord a => (a, b) -> DiffMap a b -> DiffMap a b +insertDM (key, val) = M.insertWith' (\_ old -> val:old) key [val] + +finishDM :: (b -> L.ByteString) -> DiffMap a b -> M.Map a L.ByteString +finishDM f = fmap (L.concat . map f . reverse) + +type AnnMap d = DiffMap d L.ByteString + +insertAnn :: Ord d => Ann d -> AnnMap d -> AnnMap d +insertAnn (Ann key val) = insertDM (key, val) + +finishAnn :: AnnMap d -> [Ann d] +finishAnn m = [Ann a b | (a, b) <- M.toList (finishDM id m)] + +type SeqAnnMap d = M.Map B.ByteString (AnnMap d) + +insertSM :: Ord d => B.ByteString -> Ann d -> SeqAnnMap d -> SeqAnnMap d +insertSM sq ann = M.alter (just . insertAnn ann . fromMaybe M.empty) sq + where + just !x = Just x + +finishSM :: SeqAnnMap d -> M.Map B.ByteString [Ann d] +finishSM = fmap finishAnn + +data PartialAnns = + PartialAnns { paFileAnns :: !(AnnMap FileAnnotation) + , paFileColAnns :: !(AnnMap (ColumnAnnotation InFile)) + , paSeqAnns :: !(SeqAnnMap SequenceAnnotation) + , paSeqColAnns :: !(SeqAnnMap (ColumnAnnotation InSeq)) + } + +emptyPA :: PartialAnns +emptyPA = PartialAnns M.empty M.empty M.empty M.empty + +insertPA_GF :: Ann (FileAnnotation ) -> PartialAnns -> PartialAnns +insertPA_GC :: Ann (ColumnAnnotation InFile) -> PartialAnns -> PartialAnns +insertPA_GS :: B.ByteString -> Ann (SequenceAnnotation ) -> PartialAnns -> PartialAnns +insertPA_GR :: B.ByteString -> Ann (ColumnAnnotation InSeq ) -> PartialAnns -> PartialAnns +insertPA_GF ann pa = pa { paFileAnns = insertAnn ann (paFileAnns pa) } +insertPA_GC ann pa = pa { paFileColAnns = insertAnn ann (paFileColAnns pa) } +insertPA_GS sq ann pa = pa { paSeqAnns = insertSM sq ann (paSeqAnns pa) } +insertPA_GR sq ann pa = pa { paSeqColAnns = insertSM sq ann (paSeqColAnns pa) } + + + + +---------------------------------------------------------------------- +-- [Event <-> Document] conversion functions + + +-- | Conduit that parses 'Event'@s@ into documents 'Stockholm'. +parseDoc :: C.Resource m => C.Conduit Event m Stockholm +parseDoc = C.conduitState LookingForHeader push close + where + + -- FIXME: Nice exceptions + + close LookingForHeader = return [] + close (InsideStockholm annots seqs) = return [makeStockholm annots seqs] + + push state (EvComment _) = + return (state, C.Producing []) + + push LookingForHeader EvHeader = + continue (emptyPA, M.empty) + push LookingForHeader x = + fail $ "parseDoc: unexpected " ++ show x ++ " before header" + + push (InsideStockholm _ _) EvHeader = + fail "parseDoc: unexpected header" + push (InsideStockholm annots seqs) EvEnd = + return (LookingForHeader, C.Producing [makeStockholm annots seqs]) + push (InsideStockholm annots seqs) (EvSeqData label data_) = + continue (annots, insertDM (label, data_) seqs) + push (InsideStockholm annots seqs) (EvGF feat data_) = + continue (insertPA_GF (Ann (parseFileFeature feat) data_) annots, seqs) + push (InsideStockholm annots seqs) (EvGC feat data_) = + continue (insertPA_GC (Ann (parseClmnFeature feat) data_) annots, seqs) + push (InsideStockholm annots seqs) (EvGS sq feat data_) = + continue (insertPA_GS sq (Ann (parseSeqFeature feat) data_) annots, seqs) + push (InsideStockholm annots seqs) (EvGR sq feat data_) = + continue (insertPA_GR sq (Ann (parseClmnFeature feat) data_) annots, seqs) + + continue (annots, seqs) = return (InsideStockholm annots seqs, C.Producing []) + {-# INLINE continue #-} + +data ParseDoc = LookingForHeader + | InsideStockholm + { pdAnnots :: {-# UNPACK #-} !PartialAnns + , pdSeqs :: !(DiffMap B.ByteString L.ByteString) + } + + +-- | Glue everything into place, as the Stockholm format lets +-- everything be everywhere and split in any number of parts. +makeStockholm :: PartialAnns -> DiffMap B.ByteString L.ByteString -> Stockholm +makeStockholm annots seqsDM = + let fileAnns_ = finishAnn (paFileAnns annots) + fileColAnns = finishAnn (paFileColAnns annots) + seqAnns_ = finishSM (paSeqAnns annots) + seqColAnns = finishSM (paSeqColAnns annots) + + stseqs = [StSeq (SeqLabel $ l sq) (SeqData dt) (f sq seqAnns_) (f sq seqColAnns) + | (sq, dt) <- M.toList (finishDM id seqsDM)] + where + f = M.findWithDefault [] + l = L.fromChunks . return + in Stockholm fileAnns_ fileColAnns stseqs + + +-- | Conduit that renders 'Stockholm'@s@ into 'Event'@s@. +renderDoc :: C.Resource m => C.Conduit Stockholm m Event +renderDoc = CL.concatMap toEvents + where + toEvents (Stockholm file clmn seqs) = + (EvHeader:) $ + toEventsFileAnns file $ + toEventsSeqs seqs $ + toEventsFileClmn clmn $ + [EvEnd] + + toEventsFileAnns [] = id + toEventsFileAnns (a:as) = + (EvGF (showFileFeature $ feature a) (text a) :) . + toEventsFileAnns as + + toEventsFileClmn [] = id + toEventsFileClmn (a:as) = + wrap (EvGC (showClmnFeature $ feature a)) (text a) . + toEventsFileClmn as + + toEventsSeqs (StSeq (SeqLabel name) (SeqData seqd) sa ca : xs) = + wrap (EvSeqData name') seqd . + toEventsSeqAnns name' sa . + toEventsSeqClmn name' ca . + toEventsSeqs xs + where name' = B.concat $ L.toChunks name + toEventsSeqs [] = id + + toEventsSeqAnns _ [] = id + toEventsSeqAnns n (a:as) = + (EvGS n (showSeqFeature $ feature a) (text a) :) . + toEventsSeqAnns n as + + toEventsSeqClmn _ [] = id + toEventsSeqClmn n (a:as) = + wrap (EvGR n (showClmnFeature $ feature a)) (text a) . + toEventsSeqClmn n as + + wrap :: (L.ByteString -> b) -> L.ByteString -> [b] -> [b] + wrap mk bs = case L.splitAt 70 bs of + (x, "") -> (mk x :) + (x, xs) -> (mk x :) . wrap mk xs diff --git a/src/Bio/Sequence/Stockholm/Stream.hs b/src/Bio/Sequence/Stockholm/Stream.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Sequence/Stockholm/Stream.hs @@ -0,0 +1,135 @@ +{-# LANGUAGE BangPatterns, CPP, EmptyDataDecls, DeriveDataTypeable, OverloadedStrings #-} +-- | Parsing of an Stockholm 1.0 file into a stream of events. +module Bio.Sequence.Stockholm.Stream + ( -- * Streams + Event(..) + , parseEvents + , renderEvents + ) + where + +-- from base +import Control.Applicative +import Data.Monoid (mappend) + +-- from bytestring +import qualified Data.ByteString as B +import qualified Data.ByteString.Lazy as L + +-- from conduit +import qualified Data.Conduit as C +import qualified Data.Conduit.Binary as CB +import qualified Data.Conduit.List as CL + +-- from attoparsec +import qualified Data.Attoparsec as A +import qualified Data.Attoparsec.Char8 as A8 + +-- from attoparsec-conduit +import Data.Conduit.Attoparsec (sinkParser) + +-- from blaze-builder +import qualified Blaze.ByteString.Builder as Blaze + +-- from blaze-builder-conduit +import Data.Conduit.Blaze (builderToByteString) + + +-- | An event (roughly a line in the file). +data Event = EvHeader + -- ^ @# STOCKHOLM 1.0@ + | EvEnd + -- ^ @\/\/@ + | EvComment L.ByteString + -- ^ @# ....@ + | EvSeqData B.ByteString L.ByteString + -- ^ @seqlabel seqdata@ + | EvGF B.ByteString L.ByteString + -- ^ @#GF feature data@ + | EvGC B.ByteString L.ByteString + -- ^ @#GC feature data@ + | EvGS B.ByteString B.ByteString L.ByteString + -- ^ @#GS seqlabel feature data@ + | EvGR B.ByteString B.ByteString L.ByteString + -- ^ @#GR seqlabel feature data@ + deriving (Eq, Ord, Show) + + +-- | Parse an 'Event' after 'EvHeader'. +eventParser :: A.Parser Event +eventParser = hash *> (ann <|> comment) + <|> end + <|> seqdata + where + word = A.takeTill A8.isHorizontalSpace <* spaces + tillNextLine = A.takeLazyByteString + + hash = A8.char '#' + comment = EvComment <$> tillNextLine + end = EvEnd <$ A8.string "//" <* spaces + seqdata = EvSeqData <$> word <*> tillNextLine + + ann = A8.string "=G" *> (gf <|> gc <|> gs <|> gr) + where + gf = EvGF <$ A8.char 'F' <* spaces <*> word <*> tillNextLine + gc = EvGC <$ A8.char 'C' <* spaces <*> word <*> tillNextLine + gs = EvGS <$ A8.char 'S' <* spaces <*> word <*> word <*> tillNextLine + gr = EvGR <$ A8.char 'R' <* spaces <*> word <*> word <*> tillNextLine + +-- | Parse 'EvHeader'. +headerParser :: A.Parser Event +headerParser = EvHeader <$ A8.char '#' <* spaces <* mystring "STOCKHOLM 1.0" <* spaces + where + mystring (x:xs) = A8.char x *> mystring xs + mystring [] = pure () + +spaces :: A.Parser () +spaces = A.skipWhile A8.isHorizontalSpace + + +-- | Conduit that parses a file into events. +parseEvents :: C.ResourceThrow m => C.Conduit B.ByteString m Event +parseEvents = C.sequenceSink LookingForHeader go + where + go LookingForHeader = do + dropSpaces + let emit = C.Emit InsideStockholm . (:[]) + insideLine C.=$ sinkParser $ C.Stop <$ A8.endOfInput + <|> emit <$> headerParser + + go InsideStockholm = do + dropSpaces + event <- insideLine C.=$ sinkParser eventParser + let newState = case event of + EvEnd -> LookingForHeader + _ -> InsideStockholm + return $ C.Emit newState [event] + + dropSpaces = CB.dropWhile A8.isSpace_w8 + insideLine = CB.takeWhile (/= 10) + +data ParseEvents = LookingForHeader | InsideStockholm + + +-- | Pretty print an event. +eventPrinter :: Event -> Blaze.Builder +eventPrinter ev = + case ev of + EvHeader -> bs "# STOCKHOLM 1.0\n" + EvEnd -> bs "//\n" + EvComment comment -> bs "#" <> lbs comment <> n + EvSeqData seqlabel seqdata -> bs seqlabel <> s <> lbs seqdata <> n + EvGF feature data_ -> bs "#=GF " <> bs feature <> s <> lbs data_ <> n + EvGC feature data_ -> bs "#=GC " <> bs feature <> s <> lbs data_ <> n + EvGS seqlabel feature data_ -> bs "#=GS " <> bs seqlabel <> s <> bs feature <> s <> lbs data_ <> n + EvGR seqlabel feature data_ -> bs "#=GR " <> bs seqlabel <> s <> bs feature <> s <> lbs data_ <> n + where bs = Blaze.fromByteString + lbs = Blaze.fromLazyByteString + (<>) = mappend + s = bs " " + n = bs "\n" + + +-- | Conduit that pretty prints an event stream into a file. +renderEvents :: C.ResourceUnsafeIO m => C.Conduit Event m B.ByteString +renderEvents = CL.map eventPrinter C.=$= builderToByteString diff --git a/tests/runtests.hs b/tests/runtests.hs new file mode 100644 --- /dev/null +++ b/tests/runtests.hs @@ -0,0 +1,317 @@ +{-# LANGUAGE OverloadedStrings, ScopedTypeVariables #-} + +-- from base +import Control.Applicative +import Control.Monad (zipWithM_) +import Data.List (sort) +import System.Environment (getEnv) +import System.IO.Unsafe (unsafePerformIO) + +-- from containers +import qualified Data.Map as M + +-- from bytestring +import qualified Data.ByteString.Char8 as B +import qualified Data.ByteString.Lazy.Char8 as L + +-- from biocore +import Bio.Core.Sequence + +-- from transformers +import Control.Monad.IO.Class (liftIO) + +-- from conduit +import qualified Data.Conduit as C +import qualified Data.Conduit.Binary as CB +import qualified Data.Conduit.Lazy as CZ +import qualified Data.Conduit.List as CL + +-- from zlib-conduit +import Data.Conduit.Zlib (ungzip) + +-- from QuickCheck +import Test.QuickCheck + +-- from hspec +import Test.Hspec.Monadic +import Test.Hspec.HUnit () +import Test.Hspec.QuickCheck (prop) +import Test.HUnit + +-- from this package +import Bio.Sequence.Stockholm +import Bio.Sequence.Stockholm.Document +import Bio.Sequence.Stockholm.Stream + + +main :: IO () +main = + hspecX $ do + describe "parseStockholm" $ do + it "correctly parses test file 1" $ do + ret <- strictParse stockFile + ret @?= result + it "correctly parses test file 2" $ do + ret <- strictParse stockFile2 + ret @?= result2 + + describe "renderStockholm/parseStockholm" $ do + it "parses rendered test file 1" $ do + rendered <- strictRender result + again <- strictParse rendered + again @?= result + it "parses rendered test file 2" $ do + rendered <- strictRender result2 + again <- strictParse rendered + again @?= result2 + prop "passes QuickCheck property" $ \(sto :: [Stockholm]) -> + unsafePerformIO $ do + rendered <- strictRender sto + again <- strictParse rendered + return (canonical again == canonical sto) + + describe "parseEvents" $ do + it "is able to parse gzipped RFAM_FULL" $ do + rfamFp <- getEnv "RFAM_FULL" + C.runResourceT $ + CB.sourceFile rfamFp C.$$ ungzip C.=$ parseEvents C.=$ CL.sinkNull + + describe "parseStockholm/renderStockholm/parseStockholm" $ do + it "roundtrips RFAM_SEED" $ do + rfamFp <- getEnv "RFAM_SEED" + C.runResourceT $ do + parsed1 <- CZ.lazyConsume $ CB.sourceFile rfamFp C.$= parseStockholm + parsed2 <- CZ.lazyConsume $ CL.sourceList parsed1 C.$= renderStockholm C.$= parseStockholm + liftIO (zipWithM_ (@?=) parsed2 parsed1) + + describe "renderEvents/parseEvents" $ do + prop "passes QuickCheck property" $ \(events :: [Event]) -> + unsafePerformIO $ do + rendered <- C.runResourceT $ CL.sourceList events C.$$ renderEvents C.=$ CL.consume + again <- C.runResourceT $ CL.sourceList rendered C.$$ parseEvents C.=$ CL.consume + return (canonical again == canonical events) + + -- -- Needs a better Arbitrary instance for [Event], since + -- -- eventList may generate, for example, annotations for + -- -- sequences that don't exist. + -- describe "parseDoc/renderDoc" $ do + -- prop "passes QuickCheck property" $ forAll eventList $ \(events :: [Event]) -> + -- unsafePerformIO $ do + -- rendered <- C.runResourceT $ CL.sourceList events C.$$ parseDoc C.=$ CL.consume + -- again <- C.runResourceT $ CL.sourceList rendered C.$$ renderDoc C.=$ CL.consume + -- return (canonical again == canonical events) + + describe "renderDoc/parseDoc" $ do + prop "passes QuickCheck property" $ \(docs :: [Stockholm]) -> + unsafePerformIO $ do + rendered <- C.runResourceT $ CL.sourceList docs C.$$ renderDoc C.=$ CL.consume + again <- C.runResourceT $ CL.sourceList rendered C.$$ parseDoc C.=$ CL.consume + return (canonical again == canonical docs) + + + + + +---------------------------------------------------------------------- + + +stockFile :: L.ByteString +stockFile = L.unlines [ + "# STOCKHOLM 1.0", + "#=GF AU Infernal 1.0", + "", + "#=GS Purine1 DE Number 1 :)", + "Purine1 AAAAUUGAAUAUCGUUUUACUUGUUUAUGUC-GUGAAU-UGGCAC-GACG", + "Purine2 AAAAUUUAAUAA-GAAGCACUCAUAUAAUCCCGAGAAUAUGGCUCGGGAG", + "Purine3 UGGCAGUAACUAGCGUCACUUCGUAUAACCCCAGUGAUAUGGAUUGGGGG", + "#=GC SS_cons :::::::::::::::::((((((((,,,<<<-<<<_______>>>->>>,", + "", + "# We may have comments =)", + "", + "Purine1 UUUCUACAAGGUG-CCGGAA--CACCUAACAAUAAGUAAGUCAGCAGUGA", + "Purine2 UCUCUACCGAACAACCGUAAAUUGUUCGACUAUGAGUGAAAGUGUACCUA", + "Purine3 UCUCUACCAGGAACCAAUAA--AUCCUGAUUACGAAGAGUUUAGUGCUUU", + "#=GC SS_cons ,,,,,,,<<<<<<_________>>>>>>,,))))))))::::::::::::", + "", + "Purine1 GAU", + "Purine2 GGG", + "Purine3 AGU", + "#=GC SS_cons :::", + "// "] + +purine1, purine2, purine3 :: SeqData +ss_cons :: L.ByteString +purine1 = SeqData "AAAAUUGAAUAUCGUUUUACUUGUUUAUGUC-GUGAAU-UGGCAC-GACGUUUCUACAAGGUG-CCGGAA--CACCUAACAAUAAGUAAGUCAGCAGUGAGAU" +purine2 = SeqData "AAAAUUUAAUAA-GAAGCACUCAUAUAAUCCCGAGAAUAUGGCUCGGGAGUCUCUACCGAACAACCGUAAAUUGUUCGACUAUGAGUGAAAGUGUACCUAGGG" +purine3 = SeqData "UGGCAGUAACUAGCGUCACUUCGUAUAACCCCAGUGAUAUGGAUUGGGGGUCUCUACCAGGAACCAAUAA--AUCCUGAUUACGAAGAGUUUAGUGCUUUAGU" +ss_cons = ":::::::::::::::::((((((((,,,<<<-<<<_______>>>->>>,,,,,,,,<<<<<<_________>>>>>>,,)))))))):::::::::::::::" + +result :: [Stockholm] +result = [Stockholm file clmn seqs] + where + file = [Ann AU "Infernal 1.0"] + clmn = [Ann SS ss_cons] + seqs = [mkStock "Purine1" purine1 [Ann S_DE "Number 1 :)"], + mkStock "Purine2" purine2 [], + mkStock "Purine3" purine3 []] + mkStock name data_ sa = StSeq name data_ sa [] + +stockFile2 :: L.ByteString +stockFile2 = L.unlines [stockFile, stockFile] + +result2 :: [Stockholm] +result2 = result ++ result + +sourceLBS :: C.Resource m => L.ByteString -> C.Source m B.ByteString +sourceLBS = CL.sourceList . L.toChunks + +strictParse :: L.ByteString -> IO [Stockholm] +strictParse lbs = C.runResourceT $ + sourceLBS lbs C.$= + parseStockholm C.$$ + CL.consume + +strictRender :: [Stockholm] -> IO L.ByteString +strictRender stos = fmap L.fromChunks $ + C.runResourceT $ + CL.sourceList stos C.$= + renderStockholm C.$$ + CL.consume + + +---------------------------------------------------------------------- + + +instance Arbitrary Event where + arbitrary = frequency + [ (3, EvComment <$> arbitrary) + , (10, EvSeqData <$> seqlabel <*> seqdata) + , (2, EvGF <$> feature <*> arbitrary) + , (2, EvGC <$> feature <*> arbitrary) + , (2, EvGS <$> seqlabel <*> feature <*> arbitrary) + , (2, EvGR <$> seqlabel <*> feature <*> arbitrary) + ] + where seqlabel = strict . unSL <$> arbitrary + seqdata = strict . unSD <$> arbitrary + feature = B.pack <$> listOf1 (elements alpha) + strict = B.concat . L.toChunks + +eventList :: Gen [Event] +eventList = sized $ \s -> frequency [ (100, single) + , (s, (++) <$> single <*> (resize (s `div` 2) eventList)) ] + where + single = (\xs -> EvHeader : xs ++ [EvEnd]) <$> listOf arbitrary + +instance Arbitrary Stockholm where + arbitrary = sized $ \s -> resize (min s 15) $ + Stockholm <$> arbitrary <*> arbitrary <*> arbitrary + shrink (Stockholm fileanns clmnanns stseqs) = + Stockholm <$> shrink fileanns <*> shrink clmnanns <*> shrink stseqs + + +instance Arbitrary StockholmSeq where + arbitrary = StSeq <$> arbitrary <*> arbitrary <*> arbitrary <*> arbitrary + shrink (StSeq label data_ seqanns clmnanns) = + StSeq label data_ <$> shrink seqanns <*> shrink clmnanns + +instance Arbitrary d => Arbitrary (Ann d) where + arbitrary = Ann <$> arbitrary <*> arbitrary + +instance Arbitrary FileAnnotation where + arbitrary = annArbitraryHelper list F_Other + where + list = [ AC, ID, DE, AU, SE, GA, TC, NC, TP, SQ, AM, DC + , DR, RC, RN, RM, RT, RA, RL, PI, KW, CC, NE, NL ] + +instance Arbitrary (ColumnAnnotation a) where + arbitrary = annArbitraryHelper list C_Other + where + list = [SS, SA, TM, PP, LI, AS, PAS, SAS, IN] + +instance Arbitrary SequenceAnnotation where + arbitrary = annArbitraryHelper list S_Other + where + list = [S_AC, S_DE, S_DR, OS, OC, LO] + +annArbitraryHelper :: Arbitrary b => [a] -> (b -> a) -> Gen a +annArbitraryHelper list other = + frequency $ (1, other <$> arbitrary) : + [(5, pure x) | x <- list] + +instance Arbitrary SeqLabel where + arbitrary = SeqLabel <$> (L.cons <$> elements alpha <*> (L.filter (/= ' ') <$> arbitrary)) + +instance Arbitrary SeqData where + arbitrary = SeqData . L.pack <$> listOf1 (elements ['A', 'T', 'C', 'G']) + +instance Arbitrary B.ByteString where + arbitrary = B.pack <$> (c3 <$> elements alpha + <*> listOf1 (elements $ alpha ++ " .?!|:[]{}") + <*> elements alpha) + where c3 a b c = a : b ++ [c] + +instance Arbitrary L.ByteString where + arbitrary = L.fromChunks <$> arbitrary + +alpha :: String +alpha = ['a'..'z'] ++ ['A'..'Z'] + + +---------------------------------------------------------------------- + +class Ord a => Canonical a where + canonical :: a -> a + canonical = id + + canonicalList :: [a] -> [a] + canonicalList = sort . map canonical + +instance Canonical FileAnnotation where +instance Canonical (ColumnAnnotation a) where +instance Canonical SequenceAnnotation where +instance Canonical SeqLabel where +instance Canonical SeqData where +instance Canonical B.ByteString where +instance Canonical L.ByteString where + +instance Canonical Event where + canonicalList = concat . map (glue . sort) . separate + where + separate (EvHeader : xs) = + case break (== EvEnd) xs of + (before, EvEnd : after) -> (EvHeader : before ++ [EvEnd]) : separate after + (before, rest) -> (EvHeader : before) : separate rest + separate (x:xs) = [x] : separate xs + separate [] = [] + (<>) = B.append + glue (EvSeqData sq1 data1 : EvSeqData sq2 data2 : xs) + | sq1 == sq2 = glue (EvSeqData sq1 (data1 <> data2) : xs) + glue (EvGF feat1 data1 : EvGF feat2 data2 : xs) + | feat1 == feat2 = glue (EvGF feat1 (data1 <> data2) : xs) + glue (EvGC feat1 data1 : EvGC feat2 data2 : xs) + | feat1 == feat2 = glue (EvGC feat1 (data1 <> data2) : xs) + glue (EvGS sq1 feat1 data1 : EvGS sq2 feat2 data2 : xs) + | sq1 == sq2 && feat1 == feat2 = glue (EvGS sq1 feat1 (data1 <> data2) : xs) + glue (EvGR sq1 feat1 data1 : EvGR sq2 feat2 data2 : xs) + | sq1 == sq2 && feat1 == feat2 = glue (EvGR sq1 feat1 (data1 <> data2) : xs) + glue (x1:x2:xs) = x1 : glue (x2:xs) + glue rest = rest + + +instance Canonical Stockholm where + canonical (Stockholm fileanns clmnanns stseqs) = + Stockholm (canonical fileanns) (canonical clmnanns) (canonical stseqs) + +instance Canonical StockholmSeq where + canonical (StSeq label data_ seqanns clmnanns) = + StSeq (canonical label) (canonical data_) (canonical seqanns) (canonical clmnanns) + +instance (Ord a, Canonical a) => Canonical [a] where + canonical = canonicalList + +instance Ord d => Canonical (Ann d) where + canonicalList = map mk . M.toList . toMap . map unMk + where + mk = uncurry Ann + unMk (Ann f d) = (f, d) + toMap = M.fromListWith (flip L.append)