diff --git a/Biobase/Infernal/Align.hs b/Biobase/Infernal/Align.hs new file mode 100644 --- /dev/null +++ b/Biobase/Infernal/Align.hs @@ -0,0 +1,36 @@ + +-- | "cmalign" provides two interesting results, bit scores of sequences +-- aligned to the model and the alignments themselves. + +module Biobase.Infernal.Align where + +import Data.ByteString.Char8 (ByteString) + +import Biobase.Infernal.Types + + + +-- | cmalign results, includes sequence scores if available. +-- +-- TODO stockholmAlignment, should be "biostockholm" (will be set after some +-- fun iteratee tests). For now, the 'ByteString' holds everything needed to +-- parse using biostockholm. + +data Align = Align + { modelIdentification :: ModelIdentification + , sequenceScores :: [SequenceScore] + , stockholmAlignment :: ByteString + } + +-- | Individual sequence scores. +-- +-- TODO avgProbability should use Probability newtype + +data SequenceScore = SequenceScore + { sequenceName :: !(ModelAccession,ModelIdentification,EmblAccession) -- ^ sequence name, typically RFxxxxxx;RfamID;embl-accession + , sLength :: !Int -- ^ aligned sequence length + , totalBitScore :: !BitScore -- ^ total alignment bitscore + , structureBitScore :: !BitScore -- ^ structural score part + , avgProbability :: !Double -- ^ + } + diff --git a/Biobase/Infernal/Align/Import.hs b/Biobase/Infernal/Align/Import.hs new file mode 100644 --- /dev/null +++ b/Biobase/Infernal/Align/Import.hs @@ -0,0 +1,69 @@ +{-# LANGUAGE OverloadedStrings #-} + +-- | Parses "cmalign" results. +-- +-- NOTE have not tested if this works with multiple results in a file, but +-- could ;-) + +module Biobase.Infernal.Align.Import where + +import Data.Iteratee as I +import Data.Iteratee.Char as I +import Data.Iteratee.IO as I +import Data.Iteratee.ZLib as IZ +import Data.ByteString.Char8 as BS +import Prelude as P + +import Biobase.Infernal.Align +import Biobase.Infernal.Types + + + +-- | Transforms bytestring to list of 'Align' data. + +eneeAlign :: (Monad m) => Enumeratee ByteString [Align] m a +eneeAlign = enumLinesBS ><> convStream go where + go = do + -- lets start with some comment lines + cs <- I.takeWhile (("#" ==) . BS.take 1) + -- there should be score lines now + ss <- I.takeWhile (\s -> "# STOCKHOLM 1.0" /= s && (not $ BS.null s)) + -- Stockholm lines + xs <- I.takeWhile (/="//") + x <- I.head + return [Align + { modelIdentification = ModelIdentification "" + , sequenceScores = P.map mkScore ss + , stockholmAlignment = BS.unlines $ xs++[x] + }] + +-- | Creates the required sequence score. + +mkScore s = SequenceScore + { sequenceName = undefined $ ws!!0 + , sLength = read . BS.unpack $ ws!!1 + , totalBitScore = BitScore . read . BS.unpack $ ws!!2 + , structureBitScore = BitScore . read . BS.unpack $ ws!!3 + , avgProbability = read . BS.unpack $ ws!!4 + } where ws = BS.words s + +-- | Convenience function creating all maps. + +fromFileZip :: FilePath -> IO [Align] +fromFileZip fp = run =<< ( enumFile 8192 fp + . joinI + . enumInflate GZipOrZlib defaultDecompressParams + . joinI + . eneeAlign + $ stream2stream + ) + +-- | Convenience function creating all maps. + +fromFile :: FilePath -> IO [Align] +fromFile fp = run =<< ( enumFile 8192 fp + . joinI + . eneeAlign + $ stream2stream + ) + diff --git a/Biobase/Infernal/CM.hs b/Biobase/Infernal/CM.hs --- a/Biobase/Infernal/CM.hs +++ b/Biobase/Infernal/CM.hs @@ -3,15 +3,18 @@ module Biobase.Infernal.CM where +import Data.ByteString as BS +import Data.Map as M import Data.Vector as V import Data.Vector.Unboxed as VU -import Data.ByteString as BS import Data.PrimitiveArray import Data.PrimitiveArray.Ix +import Biobase.Infernal.Types + -- | A datatype representing Infernal covariance models. This is a new -- representation that is incompatible with the one once found in "Biobase". -- The most important difference is that lookups are mapped onto efficient data @@ -33,15 +36,17 @@ -- 'localBegin' is a transition score to certain states, all such transitions -- are in 'begins'. A 'localEnd' is a transition score to a local end state. -- +-- NOTE that trustedCutoff > gathering > noiseCutoff +-- -- TODO as with other projects, we should not use Double's but "Score" and -- "Probability" newtypes. data CM = CM - { name :: ByteString - , accession :: Int - , ga :: Double - , tc :: Double - , nc :: Double + { name :: ModelIdentification -- ^ name of model as in "tRNA" + , accession :: ModelAccession -- ^ RFxxxxx identification + , trustedCutoff :: BitScore -- ^ lowest score of true member + , gathering :: BitScore -- ^ all scores at or above 'gathering' score are in the "full" alignment + , noiseCutoff :: Maybe BitScore -- ^ highest score NOT included as member , transition :: PrimArray (Int,Int) Double , emission :: PrimArray (Int,Int) Double , paths :: V.Vector (VU.Vector Double) @@ -51,3 +56,12 @@ , nodes :: V.Vector (VU.Vector Int) } deriving (Show) + +-- | Map of model names to individual CMs. + +type ID2CM = M.Map ModelIdentification CM + +-- | Map of model accession numbers to individual CMs. + +type AC2CM = M.Map ModelAccession CM + diff --git a/Biobase/Infernal/CM/Import.hs b/Biobase/Infernal/CM/Import.hs --- a/Biobase/Infernal/CM/Import.hs +++ b/Biobase/Infernal/CM/Import.hs @@ -1,34 +1,131 @@ +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE NoMonomorphismRestriction #-} +{-# LANGUAGE OverloadedStrings #-} -- | Iteratee-based parsing of Infernal covariance models. +-- +-- TODO does not create working CMs yet. Only partial key/value parsing is +-- implemented. module Biobase.Infernal.CM.Import where +import Control.Arrow +import Control.Monad (unless) +import Data.ByteString.Char8 as BS import Data.Iteratee as I import Data.Iteratee.Char as I import Data.Iteratee.IO as I import Data.Iteratee.Iteratee as I import Data.Iteratee.ListLike as I -import Data.ByteString as BS +import Data.Iteratee.ZLib as IZ +import Data.Map as M +import Prelude as P +import Control.Monad.IO.Class (liftIO, MonadIO) +import Data.PrimitiveArray +import Data.PrimitiveArray.Ix + import Biobase.Infernal.CM +import Biobase.Infernal.Types -- * iteratee stuff --- | +-- | iteratee-based parsing of human-readable CMs. eneeCM :: (Monad m) => Enumeratee ByteString [CM] m a eneeCM = enumLinesBS ><> convStream f where f = do - th <- tryHead - return undefined + -- initial (mostly key/value) data + hs' <- I.takeWhile (/="MODEL:") + let hs = M.fromList . P.map (second (BS.dropWhile (==' ')) . BS.break (==' ')) $ hs' + -- model begins + mb <- I.tryHead + unless (mb == Just "MODEL:") . error $ "model error: " ++ show (hs,mb,"head") + -- nodes + ns <- iterNodes + -- model ends + me <- I.tryHead + unless (me == Just "//") . error $ "model error: " ++ show (hs,me,"tail") + return . (:[]) $ CM + { name = ModelIdentification $ hs M.! "NAME" + , accession = ModelAccession . bsRead . BS.drop 2 $ hs M.! "ACCESSION" + , gathering = BitScore . bsRead $ hs M.! "GA" + , trustedCutoff = BitScore . bsRead $ hs M.! "TC" + , noiseCutoff = let x = hs M.! "NC" in if x == "undefined" then Nothing else Just . BitScore . bsRead $ x + , transition = error "not implemented yet" + , emission = error "not implemented yet" + , paths = error "not implemented yet" + , localBegin = error "not implemented yet" + , begins = error "not implemented yet" + , localEnd = error "not implemented yet" + , nodes = error "not implemented yet" + } where bsRead = read . BS.unpack +iterNodes :: (Monad m) => Iteratee [ByteString] m [Node] +iterNodes = do + hdr' <- I.head + let (ishdr,(hdr,nidx)) = isNodeHeader hdr' + unless ishdr $ error $ show hdr' + xs <- I.takeWhile (fst . isState) + pk <- I.peek + let n = Node + { nodeHeader = hdr + , nodeIndex = nidx + } + case pk of + Just "//" -> return [] + Just x + | (True,_) <- isNodeHeader x -> do + ns <- iterNodes + return $ n:ns + e -> error $ show e +data Node = Node + { nodeHeader :: ByteString + , nodeIndex :: Int + } + +isNodeHeader :: ByteString -> (Bool,(ByteString,Int)) +isNodeHeader xs = (isnh,(hdr,nidx)) where + isnh = BS.elem '[' xs && BS.elem ']' xs + [hdr,nidx'] = BS.words . BS.init . BS.takeWhile (/=']') . BS.drop 1 . BS.dropWhile (/='[') $ xs + nidx = read . BS.unpack $ nidx' + +isState :: ByteString -> (Bool,ByteString) +isState xs' + | P.null xs = (False,"") + | P.head xs `P.elem` [ "[", "//" ] = (False,"") + | P.head xs `P.elem` [ "S", "IL", "IR", "MATR", "MR", "D", "MP", "ML", "B", "E" ] = (True,"") + | otherwise = error $ show xs + where + xs = BS.words xs' + -- * convenience functions -- | Read covariance models from file. This parser reads one or more CMs from -- file. -fromFile :: FilePath -> [CM] -fromFile = undefined +fromFile :: FilePath -> IO (ID2CM, AC2CM) +fromFile fp = run =<< ( enumFile 8192 fp + . joinI + . eneeCM + $ I.zip (mkMap name) (mkMap accession) + ) + +-- | Read covariance models from a compressed file. + +fromFileZip :: FilePath -> IO (ID2CM, AC2CM) +fromFileZip fp = run =<< ( enumFile 8192 fp + . joinI + . enumInflate GZipOrZlib defaultDecompressParams + . joinI + . eneeCM + $ I.zip (mkMap name) (mkMap accession) + ) + +-- | map creation helper + +mkMap f = I.foldl' (\ !m x -> M.insert (f x) x m) M.empty + diff --git a/BiobaseInfernal.cabal b/BiobaseInfernal.cabal --- a/BiobaseInfernal.cabal +++ b/BiobaseInfernal.cabal @@ -1,5 +1,5 @@ name: BiobaseInfernal -version: 0.6.0.1 +version: 0.6.2.0 author: Christian Hoener zu Siederdissen maintainer: choener@tbi.univie.ac.at homepage: http://www.tbi.univie.ac.at/~choener/ @@ -26,6 +26,15 @@ . . . + Changes in 0.6.2.0 + . + * added CM parsing (implementation and interface subject to + change) + . + Changes in 0.6.1.0 + . + * added cmalign results parser + . Changes in 0.6.0.1 . * haddock should finish now @@ -44,9 +53,9 @@ library build-depends: base >3 && <5, - biocore, attoparsec, attoparsec-iteratee, + biocore, bytestring, containers, either-unwrap, @@ -59,6 +68,8 @@ exposed-modules: Biobase.Infernal + Biobase.Infernal.Align + Biobase.Infernal.Align.Import Biobase.Infernal.Clan Biobase.Infernal.Clan.Import Biobase.Infernal.CM