diff --git a/Biobase/SElab/CM.hs b/Biobase/SElab/CM.hs --- a/Biobase/SElab/CM.hs +++ b/Biobase/SElab/CM.hs @@ -1,9 +1,3 @@ -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE StandaloneDeriving #-} -{-# LANGUAGE TypeOperators #-} -{-# LANGUAGE GeneralizedNewtypeDeriving #-} -{-# LANGUAGE TemplateHaskell #-} -{-# LANGUAGE PackageImports #-} -- | Infernal CMs. -- @@ -12,240 +6,11 @@ -- TODO "fastCM :: CM -> FastCM" to make a data structure that is suitable for -- high-performance applications. -module Biobase.SElab.CM where - -import Control.Lens -import Data.ByteString.Char8 as BS -import Data.Ix (Ix) -import Data.Map as M -import Data.Primitive.Types -import Data.Vector as V -import Data.Vector.Unboxed as VU -import GHC.Base (quotInt,remInt) -import Prelude as P -import Data.List (genericLength) - -import Data.Array.Repa.Index - -import Data.Array.Repa.Index as R -import Data.Array.Repa.Shape as R -import Data.Array.Repa.ExtShape as R - -import Biobase.SElab.Types -import qualified Biobase.SElab.HMM as HMM - - - --- | Encode the CM versions we can parse - -data CMVersion - = Infernal10 BS.ByteString - | Infernal11 BS.ByteString - deriving (Eq,Ord,Show,Read) - --- | Encode CM node types. - -data NodeType - = BIF - | MATP - | MATL - | MATR - | BEGL - | BEGR - | ROOT - | END - deriving (Eq,Ord,Enum,Show,Read) - --- | Node IDs - -newtype NodeID = NodeID {unNodeID :: Int} - deriving (Eq,Ord,Show,Read) - --- | Encode CM state types. - -data StateType - = D - | MP - | ML - | MR - | IL - | IR - | S - | E - | B - | EL - deriving (Eq,Ord,Enum,Show,Read) - --- | State IDs - -newtype StateID = StateID {unStateID :: Int} - deriving (Eq,Ord,Show,Read,Prim,Ix,Enum,Num) - -illegalState = StateID $ -1 - --- | Certain states (IL,IR,ML,MR) emit a single nucleotide, one state emits a --- pair (MP), other states emit nothing. - -data Emits - = EmitsSingle { _single :: [(Char, BitScore)] } - | EmitsPair { _pair :: [(Char, Char, BitScore)] } - | EmitNothing - deriving (Eq,Ord,Show,Read) - -makeLenses ''Emits - --- | A single state. - -data State = State - { _stateID :: StateID -- ^ The ID of this state - , _nodeID :: NodeID -- ^ to which node does this state belong - , _nodeType :: NodeType -- ^ node type for this state - , _stateType :: StateType -- ^ type of the state - , _transitions :: [(StateID,BitScore)] -- ^ which transitions, id and bitscore - , _emits :: Emits -- ^ do we emit characters - } deriving (Eq,Ord,Show,Read) - -makeLenses ''State - --- | This is an Infernal covariance model. We have a number of blocks: --- --- - basic information like the name of the CM, accession number, etc. --- --- - advanced information: nodes and their states, and the states themselves. --- --- - unsorted information from the header / blasic block --- --- The 'CM' data structure is not suitable for high-performance applications. --- --- - score inequalities: trusted (lowest seed score) >= gathering (lowest full --- score) >= noise (random strings) --- --- --- --- Local entries into the CM. --- --- The "localBegin" lens returns a map of state id's. We either have just the --- root node (with the "S" state), or a set of states with type: MP,ML,MR,B. --- --- The "localEnd" lens on the other hand is the set of possible early exits --- from the model. - -data CM = CM - { _name :: Identification Rfam -- ^ name of model as in "tRNA" - , _accession :: Accession Rfam -- ^ RFxxxxx identification - , _version :: CMVersion -- ^ We can parse version 1.0 and 1.1 CMs - , _trustedCutoff :: BitScore -- ^ lowest score of any seed member - , _gathering :: BitScore -- ^ all scores at or above 'gathering' score are in the "full" alignment - , _noiseCutoff :: Maybe BitScore -- ^ highest score NOT included as member - , _nullModel :: VU.Vector BitScore -- ^ Null-model: categorical distribution on ACGU - - , _nodes :: M.Map NodeID (NodeType,[StateID]) -- ^ each node has a set of states - , _states :: M.Map StateID State -- ^ each state has a type, some emit characters, and some have children - - , _localBegin :: M.Map StateID BitScore -- ^ Entries into the CM. - , _localEnd :: M.Map StateID BitScore -- ^ Exits out of the CM. - - , _unsorted :: M.Map ByteString ByteString -- ^ all lines that are not handled. Multiline entries are key->multi-line entry - , _hmm :: Maybe HMM.HMM3 - } deriving (Show,Read) - -makeLenses ''CM - - - --- | Map of model names to individual CMs. - -type ID2CM = M.Map (Identification Rfam) CM - --- | Map of model accession numbers to individual CMs. - -type AC2CM = M.Map (Accession Rfam) CM - --- | Make a CM have local start/end behaviour, with "pbegin" and "pend" --- probabilities given. - -makeLocal :: Double -> Double -> CM -> CM -makeLocal pbegin pend cm = makeLocalEnd pend $ makeLocalBegin pbegin cm - --- | Insert all legal local beginnings, disable root node (and root states). --- The 'pbegin' probability the the total probability for local begins. The --- remaining "1-pbegin" is the probability to start with node 1. - -makeLocalBegin :: Double -> CM -> CM -makeLocalBegin pbegin cm = cm{_localBegin = lb} where - lb = M.fromList . P.map (\s -> (s^.stateID, if s^.nodeID==NodeID 1 then prob2Score 1 (1-pbegin) else prob2Score 1 (pbegin/l))) $ ss - l = genericLength ss - ss = P.filter (\s -> s^.stateType `P.elem` [MP,ML,MR,B]) . M.elems $ cm ^. states - --- | Insert all legal local ends. - -makeLocalEnd :: Double -> CM -> CM -makeLocalEnd pend cm = cm{_localEnd = le} where - le = M.fromList . P.map (\s -> (s^.stateID, prob2Score 1 (pend/l))) $ ss - l = genericLength ss - ss = P.filter (\s -> s^.stateType `P.elem` [MP,MP,MR,S] && s^.nodeType/=ROOT && notEnding s) . M.elems $ cm^.states - -- no local end, if the next node ends anyway - notEnding s = not . P.any (==E) . P.map ((^.stateType) . ((cm^.states) M.!) . fst) $ s^.transitions - - - --- Instances - -instance Shape sh => Shape (sh:.StateID) where - - rank (sh:._) = rank sh + 1 - {-# INLINE rank #-} - - zeroDim = zeroDim :. (StateID 0) - {-# INLINE zeroDim #-} - - unitDim = unitDim :. (StateID 1) - {-# INLINE unitDim #-} - - intersectDim (sh1 :. StateID n1) (sh2 :. StateID n2) = intersectDim sh1 sh2 :. StateID (min n1 n2) - {-# INLINE intersectDim #-} - - addDim (sh1 :. StateID n1) (sh2 :. StateID n2) = addDim sh1 sh2 :. StateID (n1+n2) - {-# INLINE addDim #-} - - size (sh :. StateID n) = R.size sh * n - {-# INLINE size #-} - - sizeIsValid (sh :. StateID n) - | R.size sh > 0 = n <= maxBound `div` R.size sh - | otherwise = False - {-# INLINE sizeIsValid #-} - - toIndex (sh1 :. StateID n1) (sh2 :. StateID n2) = toIndex sh1 sh2 * n1 + n2 - {-# INLINE toIndex #-} - - fromIndex (ds :. StateID d) n = fromIndex ds (n `quotInt` d) :. StateID r where - r | rank ds == 0 = n - | otherwise = n `remInt` d - {-# INLINE fromIndex #-} - - inShapeRange (sh1 :. StateID n1) (sh2 :. StateID n2) (zs :. StateID z) = (z >= n1) && (z < n2) && inShapeRange sh1 sh2 zs - {-# INLINE inShapeRange #-} - - listOfShape (sh :. StateID n) = n : listOfShape sh - {-# INLINE listOfShape #-} - - shapeOfList xx - = case xx of - [] -> error $ "shapeOfList empty in StateID" - (x:xs) -> shapeOfList xs :. StateID x - {-# INLINE shapeOfList #-} - - deepSeq (sh :. n) x = deepSeq sh (n `seq` x) - {-# INLINE deepSeq #-} - - - -instance ExtShape sh => ExtShape (sh:.StateID) where - - subDim (sh1 :. StateID n1) (sh2 :. StateID n2) = subDim sh1 sh2 :. StateID (n1-n2) - {-# INLINE subDim #-} +module Biobase.SElab.CM + ( module Biobase.SElab.CM.Types + , module Biobase.SElab.CM.Import + ) where - rangeList (sh1 :. StateID n1) (sh2 :. StateID n2) = [sh :. StateID n | sh <- rangeList sh1 sh2, n <- [n1 .. (n1+n2)] ] - {-# INLINE rangeList #-} +import Biobase.SElab.CM.Import (cmFromFile) +import Biobase.SElab.CM.Types diff --git a/Biobase/SElab/CM/Import.hs b/Biobase/SElab/CM/Import.hs --- a/Biobase/SElab/CM/Import.hs +++ b/Biobase/SElab/CM/Import.hs @@ -1,242 +1,246 @@ -{-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE DoAndIfThenElse #-} -{-# LANGUAGE ViewPatterns #-} -{-# LANGUAGE ViewPatterns #-} -{-# LANGUAGE TemplateHaskell #-} -{-# LANGUAGE PatternGuards #-} -{-# LANGUAGE BangPatterns #-} -{-# LANGUAGE NoMonomorphismRestriction #-} -{-# LANGUAGE OverloadedStrings #-} --- | Parses text-based covariance-model descriptions. +-- | Parses text-based covariance-model descriptions. This parser is +-- Utf8-aware. module Biobase.SElab.CM.Import where -import Control.Applicative -import Control.Arrow -import Control.Lens -import Control.Monad.IO.Class -import Control.Monad.IO.Class (liftIO, MonadIO) -import Control.Monad (unless) -import Data.Attoparsec.ByteString as AB -import Data.ByteString.Char8 as BS -import Data.ByteString.Lex.Double as BS -import Data.Char (isSpace,isAlpha,isDigit) -import Data.Conduit as C -import Data.Conduit.Attoparsec -import Data.Conduit.Binary as CB -import Data.Conduit.List as CL -import Data.Map as M -import Data.Maybe as M -import Data.Tuple.Select -import Data.Vector.Unboxed as VU (fromList) -import Prelude as P -import System.IO (stdout) +import Control.Applicative ( (<|>), pure, (<$>), (<$), (<*>), (<*) ) +import Control.DeepSeq (($!!)) +import Control.Lens +import Control.Monad (forM_) +import Control.Monad.IO.Class (MonadIO) +import Data.Attoparsec.ByteString (takeTill,count,many1,(),manyTill,option) +import Data.ByteString.Char8 (ByteString) +import Data.Default +import Data.Text.Encoding (decodeUtf8) +import Data.Text (Text) +import Data.Text (unpack) +import Data.Vector.Generic (fromList,empty,toList) +import Data.Vector.Generic.Lens +import Data.Vector (Vector) +import Debug.Trace +import qualified Data.Attoparsec.ByteString.Char8 as ABC +import qualified Data.ByteString.Char8 as BSC +import qualified Data.List as L +import qualified Data.Map as M +import qualified Data.Text as T +import qualified Data.Vector.Generic as VG +import System.FilePath (takeExtension) +import System.IO (stdin) -import Data.PrimitiveArray -import Data.PrimitiveArray.Zero +import Biobase.Primary.Letter +import Biobase.Primary.Nuc.RNA +import Biobase.Types.Accession (Accession(Accession),Rfam,retagAccession) +import Biobase.Types.Bitscore +import Data.PrimitiveArray hiding (fromList,map,toList) -import Biobase.SElab.CM -import Biobase.SElab.Types -import qualified Biobase.SElab.HMM as HMM -import qualified Biobase.SElab.HMM.Import as HMM +import Biobase.SElab.CM.ModelStructure +import Biobase.SElab.CM.Types +import Biobase.SElab.Common.Parser +import Biobase.SElab.HMM.Import (parseHMM) +import Biobase.SElab.HMM.Types (HMM) +import qualified Biobase.SElab.CM.Types as CM +import qualified Biobase.SElab.HMM.Types as HMM --- * Covariance model parsing. - --- ** Infernal 1.0 and 1.1 covariance model parser - -parseHeader = ($) <$ AB.string "INFERNAL" *> (Infernal10 <$ AB.string "-1" <|> Infernal11 <$ AB.string "1/a") <*> AB.takeByteString "INFERNAL line" +-- | Stream a 'ByteString' into 'CM's. +-- +-- NOTE Each CM is /always/ followed by the corresponding filter HMM. +-- (Infernal 1.1.1. at least) +-- +-- TODO this should yield @Either CM HMM@. Internally we check if part of +-- the stream @[... , CM, HMM, ...]@. We might want to provide a function +-- @mergeCmHmm@ that merges consecutive @CMs@ and @HMMs@ into the @CM@ but +-- still leaves the unmerged ones separately. Maybe we then want the +-- @these@ package, which has @a , b, (a,b)@ style data types. -lineParser p = CL.head >>= \x -> return . maybe (error "no more input") (either (\e -> error $ show (e,x)) id . AB.parseOnly p) $ x +--conduitCM :: (Monad m, MonadIO m, MonadThrow m) => Conduit ByteString m CM +--conduitCM = decodeUtf8 =$= conduitParserEither (parseCM "CM parser") =$= awaitForever (either (error . show) (yield . snd)) --- | Top-level parser for Infernal 1.0 and 1.1 human-readable covariance --- models. Reads all lines first, then builds up the CM. +-- | Simple convenience function for parsing HMM's without a lot of +-- fancyness. -parseCM1x :: (Monad m, MonadIO m) => Conduit ByteString m CM -parseCM1x = CB.lines =$= CL.sequence go where - go = do - hdr <- lineParser parseHeader - hs <- parseHeaders [] - ns <- parseNodes hdr [] - let nsMap = M.fromList . P.map (\n -> (sel2 n, (sel1 n, P.map (^. stateID) $ sel3 n))) $ ns - let ssMap = M.fromList . P.map ((^. stateID) &&& id) . P.concatMap (sel3) $ ns - lineParser $ (AB.string "//" "model end") - pk <- CL.peek - hmm <- case HMM.legalHMM pk of - True -> Just `fmap` HMM.parseHMM3 - False -> return Nothing - return CM - { _name = IDD $ hs M.! "NAME" - , _accession = ACC . readAccession . P.head . M.catMaybes $ P.map (`M.lookup` hs) ["ACC", "ACCESSION"] - , _version = hdr - , _trustedCutoff = BitScore . readBS $ hs M.! "TC" - , _gathering = BitScore . readBS $ hs M.! "GA" - , _noiseCutoff = (BitScore . readBS) `fmap` (M.lookup "NC" hs) - , _nullModel = VU.fromList . P.map readBitScore . BS.words $ hs M.! "NULL" +cmFromFile :: FilePath -> IO [CM] +cmFromFile fp = do + bs <- BSC.readFile fp + case ABC.parseOnly (many1 parseCM <* ABC.endOfInput) bs of + Left err -> error err + Right xs -> return xs - , _nodes = nsMap - , _states = ssMap +-- | Parser for covariance models (CMs). Will switch to specialized parsing +-- depending on the model version. - , _localBegin = flip M.singleton (BitScore 0) . (^.stateID) . P.head . P.filter (\s -> s^.stateType == S && s^.nodeID == NodeID 0 ) . M.elems $ ssMap - , _localEnd = M.empty +parseCM :: ABC.Parser CM +parseCM = do + pre <- parsePreCM + bdy <- parseCMBody pre + hm <- parseHMM + return $ set hmm (over HMM.accession retagAccession hm) bdy - , _unsorted = M.filter (not . flip P.elem ["NAME","ACCESSION","TC","GA","NC","NULL"]) hs - , _hmm = hmm - } +-- | -readBS = read . BS.unpack -readBitScore "*" = BitScore $ -1/0 -readBitScore x = BitScore . readBS $ x +parsePreCM :: ABC.Parser CM -- (CM, Text) +parsePreCM = do + v <- acceptedVersion + let cm' = version .~ v $ def + ls <- manyTill cmHeader ("CM" <|> "MODEL:") + let cm = L.foldl' (\a l -> l a) cm' ls + -- remainder <- ABC.takeText + return cm -- (cm, remainder) -readAccession xs - | BS.length xs /= 7 = error $ "can't read accession: " ++ BS.unpack xs - | "RF" == hdr && P.all isDigit tl = read tl - | otherwise = error $ "readAccession: " ++ BS.unpack xs - where (hdr,tl) = second BS.unpack . BS.splitAt 2 $ xs +parseCMBody :: CM -> ABC.Parser CM +parseCMBody cm = do + let v = cm^.version + nss <- case v of + -- parsing of 1.x versions + (vv,_) | "1." `T.isPrefixOf` vv -> manyTill node1x "//" + (vv,_) | "0.7" `T.isPrefixOf` vv -> manyTill node07 "//" + err -> error $ show err + ABC.endOfLine <|> pure () + buildCM nss cm def --- | Infernal 1.0 header parser. Greps all lines until the "MODEL:" line, then --- return lines to top-level parser. Parses three lines at once in case of --- "FT-" lines. +-- | We have all the parts, just need to fill up the optimized 'States' +-- data structure. +-- +-- TODO make sure that @ss@ is ordered by @sid@ and that there are no +-- missing states! -parseHeaders hs = do - p <- CL.head - case p of - (finishedHeader -> True) -> return . M.fromList - . P.map (second (BS.dropWhile isSpace) - . BS.break isSpace) - . P.reverse - $ hs - Nothing -> error $ "unexpected end of header, until here:" ++ (show $ P.reverse hs) - Just "" -> error "empty line" - Just l -> do ls <- if ("FT-" `isPrefixOf` l) then CL.take 2 else return [] - let lls = BS.concat $ l:ls - parseHeaders (lls:hs) +buildCM :: [((PInt () NodeIndex, Node),[(PInt () StateIndex, State)])] -> CM -> HMM Rfam -> ABC.Parser CM +buildCM nss cm cmhmm = do + let ns = M.fromList $ map fst nss + let ss = M.fromList $ concatMap snd nss + let cm' = set hmm cmhmm + $ set CM.cm (Left $ FlexibleModel { _fmStates = ss, _fmNodes = ns }) + $ cm + return $!! cm' -finishedHeader :: Maybe ByteString -> Bool -finishedHeader (Just x) = go x where - go "MODEL:" = True - go "CM" = True - go _ = False -finishedHeader _ = False +acceptedVersion :: ABC.Parser (T.Text,T.Text) +acceptedVersion = (new "new") <|> (old "old") "acceptedVersion" + where new = (,) <$ "INFERNAL1/a [" <*> (decodeUtf8 <$> ABC.takeTill (=='|') "x-|") <*> (eolT "|->") + old = (,"") <$ "INFERNAL-1 [" <*> (decodeUtf8 <$> ABC.takeTill (==']') "decode") <* (eolT "endOfLine") --- | Parses nodes. Will terminate on "//" which ends a CM. The state parser --- will just peek on "//", not remove it from the stream. +-- | Parse CM header information. -- --- A node is (node type, node id, set of states) +-- TODO not all header information is stored in the structure yet. -parseNodes hdr ns = do - p <- CL.peek - case (BS.dropWhile isAlpha `fmap` p) of - Nothing -> error "unexpected empty line" - Just "//" -> return . P.reverse $ ns - (isNode -> Just (ntype,nid)) -> do _ <- CL.head -- kill the line - ss <- parseStates hdr ntype nid [] - parseNodes hdr $ (ntype,nid,ss):ns +cmHeader :: ABC.Parser (CM -> CM) +cmHeader = ABC.choice + [ set name <$ "NAME" <*> eolT "name" + , set description <$ "DESC" <*> eolT "description" + , set statesInModel <$ "STATES" <*> eolN "states" + , set nodesInModel <$ "NODES" <*> eolN "nodes" + , set clen <$ "CLEN" <*> eolN "clen" + , set w <$ "W" <*> eolN "w" + , set alph <$ "ALPH" <*> eolT "alph" + , set referenceAnno <$ "RF" <*> eolB "rf" + , set consensusRes <$ "CONS" <*> eolB "cons" + , set alignColMap <$ "MAP" <*> eolB "map" + , set date <$ "DATE" <*> eolT "date" + , set pbegin <$ "PBEGIN" <*> eolD "pbegin" + , set pend <$ "PEND" <*> eolD "pend" + , set wbeta <$ "WBETA" <*> eolD "wbeta" + , set qdbBeta1 <$ "QDBBETA1" <*> eolD "qdbbeta1" + , set qdbBeta2 <$ "QDBBETA2" <*> eolD "qdbbeta2" + , set n2Omega <$ "N2OMEGA" <*> eolD "n2omega" + , set n3Omega <$ "N3OMEGA" <*> eolD "n3omega" + , set elseLF <$ "ELSELF" <*> eolD "elself" + , set nseq <$ "NSEQ" <*> eolN "nseq" + , set effn <$ "EFFN" <*> eolD "effn" + , set cksum <$ "CKSUM" <*> eolN "cksum" + , set ga <$ "GA" <*> eolD "ga" + , set tc <$ "TC" <*> eolD "tc" + , set accession . Accession <$ "ACC" <*> eolT "hmmAccession" + , ecm "ECMLC" + , ecm "ECMGC" + , ecm "ECMLI" + , ecm "ECMGI" + , set efp7gf <$> ((,) <$ "EFP7GF" <*> ssD <*> eolD "efp7gf") + , set nullModel . fromList . map Bitscore <$ "NULL" <*> ABC.count 4 ssD <* eolS "null" + , (\s -> over commandLineLog (|> decodeUtf8 s)) <$ "COM" <*> eolS "com" + , (\x -> over unknownLines (|> decodeUtf8 x)) <$> ABC.takeWhile1 (/='\n') <* ABC.take 1 + ] "cmHeader" + where + ecm s = (\a b c d e f -> set ecmlc (EValueParams a b c d e f)) <$ s <*> ssD <*> ssD <*> ssD <*> ssN <*> ssN <*> ssD <* eolS "ecm parser" --- | Parses all states for a node. We peek at the first line, then handle --- accordingly: if "//" the model will be done; is a node is coming up, return --- the state lines read until now. +-- | Parse a node together with the attached states. -parseStates hdr ntype nid xs = do - p <- CL.peek - case (BS.dropWhile isSpace `fmap` p) of - Nothing -> error "unexpected empty state" - Just "//" -> return . P.reverse $ xs - (isNode -> Just _) -> return . P.reverse $ xs - _ -> do Just x <- CL.head - let psx = parseState hdr ntype nid x - parseStates hdr ntype nid (psx:xs) +node1x :: ABC.Parser ((PInt () NodeIndex, Node), [(PInt () StateIndex, State)]) +node1x = + (\n ss -> (n & _2 . nodeStates .~ (VG.fromList $ map fst ss), ss)) <$> aNode <*> ABC.many1 (aState True) "node1x" + where + aNode = (\nty nid mapl mapr conl conr refl refr -> (nid, Node nty empty mapl mapr conl conr refl refr)) + <$ ABC.skipSpace <* ("[ " "[ ") + <*> anType <*> (ssN "node ID") + <* ABC.skipSpace <* "]" + <*> (ssN_ "mapL") <*> (ssN_ "mapR") + <*> (ssC "consL") <*> (ssC "consR") + <*> (ssC "rfL") <*> (eolC "rfR") "aNode" --- parseState :: ByteString -> State -parseState hdr ntype nid s - | P.null ws = error "parseState: no words" - | B == t = State { _stateID = StateID . readBS $ pn!!0 - , _stateType = t - , _nodeID = nid - , _nodeType = ntype - , _transitions = [ ( StateID . readBS $ pn!!3, 0) - , ( StateID . readBS $ pn!!4, 0) - ] - , _emits = EmitNothing - } - | otherwise = State { _stateID = StateID . readBS $ pn!!0 - , _stateType = t - , _nodeID = nid - , _nodeType = ntype - , _transitions = [ (StateID (i+k), readBitScore $ ts!!k) | k <- [0..n-1]] - , _emits = e - } +node07 :: ABC.Parser ((PInt () NodeIndex, Node), [(PInt () StateIndex, State)]) +node07 = + -- Update the node with the indices from all the states belonging with + -- this node. + (\n ss -> (n & _2 . nodeStates .~ (VG.fromList $ map fst ss), ss)) <$> aNode <*> ABC.many1 (aState False) "node07" where - ws = BS.words s - numPN = case hdr of - Infernal10 _ -> 5 - Infernal11 _ -> 9 -- last 4 values are QDB values ... - numTS = readBS $ pn!!4 - numES = case w of - "MP" -> 16 - (flip P.elem ["ML","MR","IL","IR"] -> True) -> 4 - _ -> 0 - ~([w],~(pn,~(ts,es))) = second (second (second (P.map readBitScore) . P.splitAt numTS) . P.splitAt numPN) . P.splitAt 1 $ ws - t = readBS w :: StateType - i = readBS $ pn!!3 - n = readBS $ pn!!4 - e = case t of - MP -> EmitsPair $ P.zipWith (\(c1,c2) k -> (c1,c2,k)) [ (c1,c2) | c1 <- "ACGU", c2 <- "ACGU" ] es - ((flip P.elem [ML,MR,IL,IR]) -> True) -> EmitsSingle $ P.zip "ACGU" es - _ -> EmitNothing + aNode = (\nty nix -> (nix, Node nty empty 0 0 '-' '-' '-' '-')) + <$ ABC.skipSpace <* ("[ " "[ ") + <*> anType <*> (ssN "node ID") + <* ABC.skipSpace <* "]" +-- | Parse the node type +anType :: ABC.Parser NodeType +anType = ABC.choice [ Bif <$ "BIF" , MatP <$ "MATP", MatL <$ "MATL", MatR <$ "MATR" + , BegL <$ "BEGL", BegR <$ "BEGR", Root <$ "ROOT", End <$ "END" ] "anType" -{- -parseState hdr ntype nid s - | B == t = State { _stateID = StateID . readBS $ ws!!1 - , _stateType = B - , _nodeID = nid - , _nodeType = ntype - , _transitions = [ ( StateID . readBS $ ws!!4, 0) - , ( StateID . readBS $ ws!!5, 0) - ] - , _emits = EmitNothing - } - | otherwise = State { _stateID = StateID . readBS $ ws!!1 - , _stateType = t -- stateTypeFromString . BS.unpack $ t - , _nodeID = nid - , _nodeType = ntype - , _transitions = [ (StateID (i+k), readBitScore $ ws!!(6+k)) - | k <- [0 .. n-1] ] - , _emits = e - } - where - last k = P.map readBitScore . P.reverse . P.take k . P.reverse $ ws - (t':_) = ws - n = readBS $ ws!!5 -- number of states - i = readBS $ ws!!4 -- first state - e = case t of - MP -> EmitsPair $ P.zipWith (\(c1,c2) k -> (c1,c2,k)) [ (c1,c2) | c1 <- "ACGU", c2 <- "ACGU" ] (last 16) - ((flip P.elem [ML,MR,IL,IR]) -> True) -> EmitsSingle . P.zip "ACGU" $ last 4 - _ -> EmitNothing --} +-- | Parsing of states. --- | Determine if a line is a node line ('Just'). If yes, we'll get the node --- type as string and the node identifier, too. +aState + :: Bool -- ^ Control if the four QDB parameters are present. Currently 'node1x' will parse those, 'node07' won't. + -> ABC.Parser (PInt () StateIndex, State) +aState parseQDB = do + ABC.skipSpace + _stateType <- asType "asType" + stateId <- ssN + _stateParents <- (\a b -> VG.fromList [a,b]) <$> ssZ <*> ssN + (c1,c2) <- (,) <$> ssZ <*> ssN + let chdn = if _stateType==B then [ PInt c1 , PInt c2 ] + else take c2 [ PInt c1 .. ] + _stateQDB <- if parseQDB + then QDB <$> ssN <*> ssN <*> ssN <*> ssN + else pure def + _stateTransitions <- if | _stateType == B -> pure $ fromList $ map (,0) chdn + | otherwise -> (fromList . zip chdn) <$> ABC.count c2 (Bitscore <$> ssD') + _stateEmissions <- if | emitsPair _stateType -> fromList <$> ABC.count 16 (Bitscore <$> ssD) + | emitsSingle _stateType -> fromList <$> ABC.count 4 (Bitscore <$> ssD) + | otherwise -> pure empty + eolS + return (stateId,State{..}) -isNode :: Maybe ByteString -> Maybe (NodeType, NodeID) -isNode (Just xs) - | BS.null xs = Nothing - | ("[":ntype:nid:"]":cm11) <- BS.words xs = Just (readBS ntype, NodeID . readBS $ nid) -isNode _ = Nothing +-- | Type of the state -fromFile :: FilePath -> IO [CM] -fromFile fp = do - runResourceT $ sourceFile fp $= parseCM1x $$ consume +asType :: ABC.Parser StateType +asType = ABC.choice [ D <$ "D" , MP <$ "MP", ML <$ "ML", MR <$ "MR" + , IL <$ "IL", IR <$ "IR", S <$ "S" , E <$ "E" , B <$ "B" ] -test :: IO () -test = do - xs10 <- runResourceT $ sourceFile "test10.cm" $= parseCM1x $$ consume -- sinkHandle stdout - xs11 <- runResourceT $ sourceFile "test11.cm" $= parseCM1x $$ consume -- sinkHandle stdout - print xs10 - print xs11 - return () +-- | Read a list of CMs from a given filename. The special filename @-@ +-- reads from @stdin@, while a suffix ending in @.gz@ will pipe through +-- @ungzip@ before parsing the contents. +--fromFile :: FilePath -> IO [CM] +--fromFile fp +-- | fp == "-" = runResourceT $ sourceHandle stdin $= conduitCM $$ consume +-- | takeExtension fp == ".gz" = runResourceT $ sourceFile fp $= ungzip $= conduitCM $$ consume +-- | otherwise = runResourceT $ sourceFile fp $= conduitCM $$ consume +-- +--test = do +-- cms11 <- fromFile "RF00563.cm" +-- --cms07 <- fromFile "rebecca-kirsch/split_split_chr3L_289_0.maf.gz.fa.cm.h1.3.h2.5" +-- forM_ cms11 $ \cm -> do +-- --print $ cm ^. unknownLines +-- --print $ makeLocal cm +-- --print $ (addLocalEnds cm) ^? nodes . vectorIx 1 . nstates . ix 0 . transitions +-- let q = (addLocalEnds cm) & nodes . vectorIx 1 . nodeMainState EntryState . transitions %~ id +-- mapM_ print $ VG.toList $ q ^. nodes . vectorIx 1 . nstates . ix 0 . transitions +-- diff --git a/Biobase/SElab/CM/ModelStructure.hs b/Biobase/SElab/CM/ModelStructure.hs new file mode 100644 --- /dev/null +++ b/Biobase/SElab/CM/ModelStructure.hs @@ -0,0 +1,668 @@ + +-- | Defines two model structures. One structure is designed to be easily +-- modifiable for working with a CM. The second is "static" but efficient +-- to use in applications. An isomorphism between the two structures is +-- provided. +-- +-- TODO Generalize to both, HMMs and CMs. This will require some thinking +-- on how to generalize everything from individual states to emission +-- systems. Emissions can probably be phantom-typed so that we know +-- emission orders, and other things. + +module Biobase.SElab.CM.ModelStructure where + +import Control.DeepSeq +import Control.Lens +import Data.Aeson (FromJSON,ToJSON) +import Data.Binary (Binary) +import Data.Default +import Data.Function (on) +import Data.Hashable (Hashable) +import Data.Ix (Ix) +import Data.Map (Map) +import Data.Serialize (Serialize) +import Data.Set (Set) +import Data.Vector.Unboxed.Deriving +import Debug.Trace +import GHC.Generics (Generic) +import qualified Data.List as L +import qualified Data.Map as M +import qualified Data.Vector as V +import qualified Data.Vector.Generic as VG +import qualified Data.Vector.Unboxed as VU +import Text.Read + +import Biobase.Primary.Letter +import Biobase.Primary.Nuc.RNA +import Biobase.Types.Bitscore +import Data.PrimitiveArray hiding (fromList,toList,map) + + + +-- * General things + +-- | Phantom-type a node index of @PInt@s + +data NodeIndex + +-- | Phantom-type a state index @PInt@s + +data StateIndex + +-- | The type of a node, efficiently encoded as an Int. + + + +newtype NodeType = NodeType Int + deriving (Eq,Ord,Generic,Ix) + +pattern Bif = NodeType 0 +pattern MatP = NodeType 1 +pattern MatL = NodeType 2 +pattern MatR = NodeType 3 +pattern BegL = NodeType 4 +pattern BegR = NodeType 5 +pattern Root = NodeType 6 +pattern End = NodeType 7 + +instance Binary NodeType +instance FromJSON NodeType +instance Hashable NodeType +instance Serialize NodeType +instance ToJSON NodeType +instance NFData NodeType + +instance Show NodeType where + show = \case + Bif -> "BIF" + MatP -> "MATP" + MatL -> "MATL" + MatR -> "MATR" + BegL -> "BEGL" + BegR -> "BEGR" + Root -> "ROOT" + End -> "END" + +instance Read NodeType where + readPrec = parens $ do + Ident s <- lexP + return $ case s of + "BIF" -> Bif + "MATP" -> MatP + "MATL" -> MatL + "MATR" -> MatR + "BEGL" -> BegL + "BEGR" -> BegR + "ROOT" -> Root + "END" -> End + _ -> error $ "read NodeType: " ++ s + +derivingUnbox "NodeType" + [t| NodeType -> Int |] [| \(NodeType n) -> n |] [| NodeType |] + + + +-- | Type of a state, a newtype wrapper for performance + +newtype StateType = StateType Int + deriving (Eq,Ord,Generic,Ix) + +pattern D = StateType 0 +pattern MP = StateType 1 +pattern ML = StateType 2 +pattern MR = StateType 3 +pattern IL = StateType 4 +pattern IR = StateType 5 +pattern S = StateType 6 +pattern E = StateType 7 +pattern B = StateType 8 +pattern EL = StateType 9 + +instance Binary StateType +instance FromJSON StateType +instance Hashable StateType +instance Serialize StateType +instance ToJSON StateType +instance NFData StateType + +instance Show StateType where + show = \case + D -> "D" + MP -> "MP" + ML -> "ML" + MR -> "MR" + IL -> "IL" + IR -> "IR" + S -> "S" + E -> "E" + B -> "B" + EL -> "EL" + (StateType e) -> "StateType " ++ show e + +instance Read StateType where + readPrec = parens $ do + Ident s <- lexP + return $ case s of + "D" -> D + "MP" -> MP + "ML" -> ML + "MR" -> MR + "IL" -> IL + "IR" -> IR + "S" -> S + "E" -> E + "B" -> B + "EL" -> EL + _ -> error $ "read StateType: " ++ s + +derivingUnbox "StateType" + [t| StateType -> Int |] [| \(StateType s) -> s |] [| StateType |] + +emitsSingle :: StateType -> Bool +emitsSingle s | s `elem` [ML,MR,IL,IR] = True + | otherwise = False +{-# Inline emitsSingle #-} + +emitsPair = (==) MP +{-# Inline emitsPair #-} + + + +-- * QDB parameters. + +-- | Query-dependent banding parameters. The four parameters are given in +-- increasing order. They are set to @-1@ if not given. + +data QDB = QDB + { _minExpSeqLenBeta2 :: ! Int + , _minExpSeqLenBeta1 :: ! Int + , _maxExpSeqLenBeta1 :: ! Int + , _maxExpSeqLenBeta2 :: ! Int + } + deriving (Eq,Show,Read,Generic) + +makeLenses ''QDB +makePrisms ''QDB + +instance Default QDB where + def = QDB + { _minExpSeqLenBeta2 = -1 + , _minExpSeqLenBeta1 = -1 + , _maxExpSeqLenBeta1 = -1 + , _maxExpSeqLenBeta2 = -1 + } + +instance Binary QDB +instance Serialize QDB +instance FromJSON QDB +instance ToJSON QDB +instance NFData QDB + +derivingUnbox "QDB" + [t| QDB -> (Int,Int,Int,Int) |] [| \(QDB a b c d) -> (a,b,c,d) |] [| \(a,b,c,d) -> QDB a b c d |] + + + +type Transitions b = VU.Vector (PInt () StateIndex, b) + + + +-- | A single state in a model. +-- +-- TODO Map (PInt () StateIndex) State + +data State = State + { _stateType :: ! StateType + -- ^ The type of the current state + , _stateParents :: ! (VU.Vector (PInt () StateIndex)) + -- ^ List of parents into this state + , _stateQDB :: ! QDB + -- ^ QDB information + , _stateTransitions :: ! (Transitions Bitscore) + -- ^ Into which children to we transition to + , _stateEmissions :: ! (VU.Vector Bitscore) + -- ^ Finally, emission scores, if given for this state. Different + -- stochastic models should interpret this differently! + -- For covariance models, the emission order is ACGU for single states + -- or AA,AC,AG,AU, CA,CC,CG,CU, GA,GC,GG,GU, UA,UC,UG,UU for pair + -- states. + -- TODO really only one entry? + } + deriving (Eq,Show,Read,Generic) + +makeLenses ''State +makePrisms ''State + +instance Default State where + def = State + { _stateType = StateType (-1) + , _stateParents = VG.empty + , _stateQDB = def + , _stateTransitions = VG.empty + , _stateEmissions = VG.empty + } + +instance Binary State +instance Serialize State +instance FromJSON State +instance ToJSON State +instance NFData State + + + +-- * High-performance structure for @State@s. Actual calculations are run +-- on these. + +-- | Encode all the information necessary to have *efficient* covariance +-- models. +-- +-- The index @PInt () StateIndex@ is the actual index type as given in +-- a model description. +-- +-- Transitions are encoded as a boxed vector of unboxed vectors. The outer +-- boxed vector is indexed by the current state. The inner unboxed vector +-- is indexed by the child number. For each child number we record the +-- target state and transition cost. +-- +-- TODO emissions pair/single +-- TODO local / global mode +-- TODO add QDB information here? +-- +-- TODO We need to modify how BiobaseXNA encodes RNA sequences (maybe ACGUN) +-- +-- TODO ugly but more efficient? Use just a single @Emit@ data structure? + +data States = States + { _statesType :: ! (Unboxed (PInt () StateIndex) StateType) + -- ^ Type of the state at the current index + , _statesParents :: ! (Boxed (PInt () StateIndex) (VU.Vector (PInt () StateIndex))) + -- ^ For each state, record which other states lead here + , _statesTransitions :: ! (Boxed (PInt () StateIndex) (Transitions Bitscore)) + -- ^ Transitions to a state, together with the transition score; + -- unpopulated transitions are set to @-1@. + -- TODO we have "forbidden" transitions. Consider how to handle these. + -- Easy solution is very low bitscores, maybe @-neginf@? + , _statesQDB :: ! (Unboxed (PInt () StateIndex) QDB) + , _statesEmitPair :: ! (Unboxed (Z:.PInt () StateIndex:.Letter RNA:.Letter RNA) Bitscore) + -- ^ Scores for the emission of a pair + , _statesEmitSingle :: ! (Unboxed (Z:.PInt () StateIndex:.Letter RNA) Bitscore) + -- ^ Scores for the emission of a single nucleotide + } + deriving (Eq,Show,Read,Generic) + +makeLenses ''States +makePrisms ''States + +instance Default States where + def = States + { _statesType = fromAssocs 0 0 (StateType $ -1) [] + , _statesParents = fromAssocs 0 0 VG.empty [] + , _statesTransitions = fromAssocs 0 0 VG.empty [] + , _statesQDB = fromAssocs 0 0 def [] + , _statesEmitPair = fromAssocs (Z:.0:.A:.A) (Z:.0:.A:.A) 0 [] + , _statesEmitSingle = fromAssocs (Z:.0:.A) (Z:.0:.A) 0 [] + } + +instance Binary States +instance Serialize States +instance FromJSON States +instance ToJSON States +instance NFData States + +-- | A pure getter to retrieve the last state + +sLastState :: Getter States (PInt () StateIndex) +sLastState = statesType . to bounds . to snd +{-# Inline sLastState #-} + + + +-- * Nodes for dynamically changeable models. + +-- | @Node@s are a high-level structure in covariance models, with each +-- node having one or more states as children. In addition, nodes carry +-- alignment-column based information. +-- +-- TODO @_nColL@ and @nColR@ should become @Index 1@ types. We'll do that +-- once we re-activate Stockholm file parsing. + +data Node = Node + { _nodeType :: ! NodeType + -- ^ Type of the node + , _nodeStates :: ! (V.Vector (PInt () StateIndex)) + -- ^ States associated with this node + , _nodeColL :: ! Int + -- ^ Column index in the corresponding Stockholm file + , _nodeColR :: ! Int + -- ^ Column index in the corresponding Stockholm file + , _nodeConL :: ! Char + -- ^ TODO + , _nodeConR :: ! Char + -- ^ TODO + , _nodeRefL :: ! Char + -- ^ TODO + , _nodeRefR :: ! Char + -- ^ TODO + } + deriving (Eq,Ord,Show,Read,Generic) + +makeLenses ''Node +makePrisms ''Node + +instance Binary Node +instance Serialize Node +instance FromJSON Node +instance ToJSON Node +instance NFData Node + +instance Default Node where + def = Node + { _nodeType = NodeType (-1) + , _nodeStates = VG.empty + , _nodeColL = -1 + , _nodeColR = -1 + , _nodeConL = '-' + , _nodeConR = '-' + , _nodeRefL = '-' + , _nodeRefR = '-' + } + + + +-- * High-performance structure for @Node@s. + +data Nodes = Nodes + { _nodesType :: ! (Unboxed (PInt () NodeIndex) NodeType) + , _nodesStates :: ! (Boxed (PInt () NodeIndex) (V.Vector (PInt () StateIndex))) + , _nodesColL :: ! (Unboxed (PInt () NodeIndex) Int) + , _nodesColR :: ! (Unboxed (PInt () NodeIndex) Int) + , _nodesConL :: ! (Unboxed (PInt () NodeIndex) Char) + , _nodesConR :: ! (Unboxed (PInt () NodeIndex) Char) + , _nodesRefL :: ! (Unboxed (PInt () NodeIndex) Char) + , _nodesRefR :: ! (Unboxed (PInt () NodeIndex) Char) + } + deriving (Eq,Show,Read,Generic) + +makeLenses ''Nodes +makePrisms ''Nodes + +instance Default Nodes where + def = Nodes + { _nodesType = fromAssocs 0 0 (NodeType $ -1) [] + , _nodesStates = fromAssocs 0 0 VG.empty [] + , _nodesColL = fromAssocs 0 0 (-1) [] + , _nodesColR = fromAssocs 0 0 (-1) [] + , _nodesConL = fromAssocs 0 0 '-' [] + , _nodesConR = fromAssocs 0 0 '-' [] + , _nodesRefL = fromAssocs 0 0 '-' [] + , _nodesRefR = fromAssocs 0 0 '-' [] + } + +instance Binary Nodes +instance Serialize Nodes +instance FromJSON Nodes +instance ToJSON Nodes +instance NFData Nodes + + + + +data StaticModel = StaticModel + { _smStates :: ! States + , _smNodes :: ! Nodes + } + deriving (Eq,Show,Read,Generic) + +instance Binary StaticModel +instance Serialize StaticModel +instance FromJSON StaticModel +instance ToJSON StaticModel +instance NFData StaticModel + +instance Default StaticModel where + def = StaticModel + { _smStates = def + , _smNodes = def + } + +-- | Model structure that is somewhat easy to modify. Before turning this +-- into a @StaticModel@, the model itself needs to be valid. + +data FlexibleModel = FlexibleModel + { _fmStates :: ! (Map (PInt () StateIndex) State) + , _fmNodes :: ! (Map (PInt () NodeIndex ) Node ) + } + deriving (Eq,Show,Read,Generic) + +instance Binary FlexibleModel +instance Serialize FlexibleModel +instance FromJSON FlexibleModel +instance ToJSON FlexibleModel +instance NFData FlexibleModel + +makeLenses ''FlexibleModel +makePrisms ''FlexibleModel + +instance Default FlexibleModel where + def = FlexibleModel + { _fmStates = def + , _fmNodes = def + } + + + +isValidModel :: FlexibleModel -> Bool +isValidModel = error "isvalidModel: write me!" + + + +-- * Isomorphisms between static and flexible models +-- +-- @flexibleToStatic . staticToFlexible == id@ +-- @staticToFlexible . flexibleToStatic == id@ + +-- | Make a flexible model static. +-- +-- TODO should *really* do some basic tests +-- +-- TODO this would be easier if we introduced hybrid arrays, not just +-- @Unboxed@ and @Boxed@. ... or if everything were unboxed. +-- +-- TODO needs to generalize over the actual model we are dealing with. This +-- includes how many characters to emit in pair and single. And the +-- underlying alphabet. +-- +-- TODO use @isValidModel@ for tests. + +flexibleToStatic :: FlexibleModel -> StaticModel +flexibleToStatic (FlexibleModel s n) + | True = StaticModel s' n' + where s' = States + { _statesType = fromAssocs 0 mix (StateType $ -1) $ zip ix $ s ^.. traverse . stateType + , _statesParents = fromAssocs 0 mix VG.empty $ zip ix $ s ^.. traverse . stateParents + , _statesTransitions = fromAssocs 0 mix VG.empty $ zip ix $ s ^.. traverse . stateTransitions + , _statesQDB = fromAssocs 0 mix def $ zip ix $ s ^.. traverse . stateQDB + -- + , _statesEmitPair = fromAssocs (Z:.0:.A:.A) (Z:.mix:.U:.U) def $ + [ (Z:.k:.n1:.n2,e) + | (k,es) <- zip ix $ s ^.. traverse . stateEmissions + , VG.length es == 16 + , ((n1,n2),e) <- zip ((,) <$> acgu <*> acgu) (VG.toList es) + ] + , _statesEmitSingle = fromAssocs (Z:.0:.A) (Z:.mix:.U) def $ + [ (Z:.k:.n1,e) + | (k,es) <- zip ix $ s ^.. traverse . stateEmissions + , VG.length es == 4 + , ((n1),e) <- zip acgu (VG.toList es) + ] + } where ix = M.keys s ; mix = maximum ix + n' = Nodes + { _nodesType = fromAssocs 0 mix (NodeType $ -1) $ zip ix $ n ^.. traverse . nodeType + , _nodesStates = fromAssocs 0 mix VG.empty $ zip ix $ n ^.. traverse . nodeStates + , _nodesColL = fromAssocs 0 mix (-1) $ zip ix $ n ^.. traverse . nodeColL + , _nodesColR = fromAssocs 0 mix (-1) $ zip ix $ n ^.. traverse . nodeColR + , _nodesConL = fromAssocs 0 mix '-' $ zip ix $ n ^.. traverse . nodeConL + , _nodesConR = fromAssocs 0 mix '-' $ zip ix $ n ^.. traverse . nodeConR + , _nodesRefL = fromAssocs 0 mix '-' $ zip ix $ n ^.. traverse . nodeRefL + , _nodesRefR = fromAssocs 0 mix '-' $ zip ix $ n ^.. traverse . nodeRefR + } where ix = M.keys n ; mix = maximum ix + +-- | Make static model flexible again. +-- +-- Static models are always (defined to be) valid models. +-- +-- TODO emission handling for generalized models + +staticToFlexible :: StaticModel -> FlexibleModel +staticToFlexible (StaticModel States{..} Nodes{..}) + = FlexibleModel s' n' + where s' = M.fromList $ map goS $ uncurry enumFromTo $ bounds _statesType + n' = M.fromList $ map goN $ uncurry enumFromTo $ bounds _nodesType + goS k = (k,) $ State + { _stateType = t + , _stateParents = _statesParents ! k + , _stateQDB = _statesQDB ! k + , _stateTransitions = _statesTransitions ! k + , _stateEmissions = if | emitsPair t -> VG.fromList [ _statesEmitPair ! (Z:.k:.i:.j) | (i,j) <- (,) <$> acgu <*> acgu ] + | emitsSingle t -> VG.fromList [ _statesEmitSingle ! (Z:.k:.i ) | i <- acgu ] + | otherwise -> VG.empty + } where t = _statesType ! k + goN k = (k,) $ Node + { _nodeType = _nodesType ! k + , _nodeStates = _nodesStates ! k + , _nodeColL = _nodesColL ! k + , _nodeColR = _nodesColR ! k + , _nodeConL = _nodesConL ! k + , _nodeConR = _nodesConR ! k + , _nodeRefL = _nodesRefL ! k + , _nodeRefR = _nodesRefR ! k + } + + + +-- * Local / Global mode conversion + +-- | The list of all nodes and states that can be the target of a local +-- begin. These are nodes with type @MatP@, @MatL@,@MatR@, or @Bif@. They +-- will not necessarily have been set this way. Targets of a local begin +-- are *never* @Root@ nodes and their states. + +internalEntries :: FlexibleModel -> [(PInt () NodeIndex, PInt () StateIndex)] +internalEntries FlexibleModel{..} = xs + where xs = concatMap givenN $ M.toList _fmNodes + givenN (n,Node{..}) + | _nodeType == MatP = [(n, getState MP _nodeStates)] + | _nodeType == MatL = [(n, getState ML _nodeStates)] + | _nodeType == MatR = [(n, getState MR _nodeStates)] + | _nodeType == Bif = [(n, getState B _nodeStates)] + | otherwise = [] + getState ty = head . filter ((==ty) . _stateType . (_fmStates M.!)) . VG.toList + +-- | The list of all nodes and states that can be the target of a local +-- end. +-- +-- Nodes that have and @End@ node following are excluded. + +internalExits :: FlexibleModel -> [(PInt () NodeIndex, PInt () StateIndex)] +internalExits FlexibleModel{..} = xs + where xs = concatMap givenN $ M.toList _fmNodes + givenN (n,Node{..}) + | _nodeType == MatP = [(n, getState MP _nodeStates) | noNextE _nodeStates ] + | _nodeType == MatL = [(n, getState ML _nodeStates) | noNextE _nodeStates ] + | _nodeType == MatR = [(n, getState MR _nodeStates) | noNextE _nodeStates ] + | _nodeType == BegL = [(n, getState S _nodeStates) | noNextE _nodeStates ] + | _nodeType == BegR = [(n, getState S _nodeStates) | noNextE _nodeStates ] + | otherwise = [] + getState ty = VG.head . VG.filter ((==ty) . _stateType . (_fmStates M.!)) + noNextE = VG.null . VG.filter ((==E) . _stateType . (_fmStates M.!)) + +-- | Create a new transition from a given state to another given state. +-- +-- Will die with an error if any of source or target state is not in the +-- model. + +insertTransition :: PInt () StateIndex -> PInt () StateIndex -> Bitscore -> FlexibleModel -> FlexibleModel +insertTransition frm to sc mdl + | fS <- M.lookup frm (mdl^.fmStates) + , tS <- M.lookup to (mdl^.fmStates) + = mdl + -- add the backlink + & fmStates . at to . _Just . stateParents %~ addParent frm + -- add the transition itself + & fmStates . at frm . _Just . stateTransitions %~ addTransition to sc + | otherwise = error $ "insertTransition: missing state(s)" + +-- | Given a state we come from (@frm@), insert into the vector of parents. + +addParent :: PInt () StateIndex -> VU.Vector (PInt () StateIndex) -> VU.Vector (PInt () StateIndex) +addParent frm = VG.fromList . L.nub . (frm:) . VG.toList + +-- | Adds a transition at the right position in the @Transitions@ vector. +-- +-- This operation takes @O(n^2)@ time for each insert! (Though @n@ is +-- typically @<=6@. + +addTransition :: VU.Unbox s => PInt () StateIndex -> s -> Transitions s -> Transitions s +addTransition to sc ts = VG.fromList . L.nubBy ((==) `on` fst) $ VG.toList xs ++ [(to,sc)] ++ VG.toList ys + where (xs,ys) = VG.partition (( FlexibleModel -> FlexibleModel +addLocalBegins b mdl = foldl go mdl $ (,) <$> ss <*> (map snd $ internalEntries mdl) + where + -- list of states to modify. Assumed to be @S 0@ to @IR 2@. + ss = mdl ^.. fmNodes . at 0 . traverse . nodeStates . traverse + go m (f,t) = insertTransition f t b m + +-- | Given a @CM@, add the necessary transitions to create local ends. + +addLocalEnds :: Bitscore -> FlexibleModel -> FlexibleModel +addLocalEnds b mdl' = foldl go mdl . map snd $ internalExits mdl + where + go m f = insertTransition f t b m + [(t,_)] = filter ((==EL) . _stateType . snd) . M.toList $ mdl ^. fmStates + mdl = createLocalEndState mdl' + +-- | Create the @EL@ state, together with its own node. + +createLocalEndState :: FlexibleModel -> FlexibleModel +createLocalEndState mdl + | null e = mdl & fmNodes . at (maxN+1) .~ Just n & fmStates . at (maxS+1) .~ Just s + | otherwise = mdl -- we already have an @EL@ state. + where + e = filter (==EL) $ mdl ^.. fmStates . traverse . stateType + (maxN, _) = M.findMax $ mdl ^. fmNodes + (maxS, _) = M.findMax $ mdl ^. fmStates + n = def & nodeType .~ End & nodeStates .~ VG.singleton maxS + s = def & stateType .~ EL + + +-- | Perform the necessary edge insertions to make a mode "local". If in +-- doubt, use @Just 0.05@ and @Just 0.05@ as parameters for the local +-- begins and local ends. +-- +-- TODO It holds that @makeLocal b e . makeLocal b e == makeLocal b e@. +-- (Provisionary; depending on how we shall go about modifying bitscores) +-- +-- TODO implement local ends part + +makeLocal + :: () + => Maybe Bitscore + -- ^ @Just@ the local begin bitscore, or @Nothing@ if local begins are + -- not desired. + -> Maybe Bitscore + -- ^ @Just@ the local end bitscore, or @Nothing@ if local ends are not + -- desired. + -> FlexibleModel + -> FlexibleModel +makeLocal mB mE = maybe id addLocalEnds mE . maybe id addLocalBegins mB + diff --git a/Biobase/SElab/CM/Types.hs b/Biobase/SElab/CM/Types.hs new file mode 100644 --- /dev/null +++ b/Biobase/SElab/CM/Types.hs @@ -0,0 +1,296 @@ + +-- | Efficient encoding of @Infernal 1.1@ covariance models. + +module Biobase.SElab.CM.Types where + +import Control.DeepSeq +import Control.Lens +import Data.Aeson (FromJSON,ToJSON) +import Data.Binary (Binary) +import Data.Default +import Data.Hashable (Hashable) +import Data.Ix (Ix) +import Data.List (genericLength) +import Data.Sequence (Seq) +import Data.Serialize (Serialize) +import Data.Text (Text) +import Data.Vector.Generic (empty,fromList,toList) +import Data.Vector.Generic.Lens +import Data.Vector.Unboxed.Deriving +import Data.Vector.Unboxed (Vector) +import Data.Word (Word32) +import GHC.Generics (Generic) +import qualified Data.Vector as V +import qualified Data.Vector.Generic as VG +import Text.Read +import Debug.Trace + +import Biobase.Primary.Letter +import Biobase.Primary.Nuc.RNA +import Biobase.Types.Accession +import Biobase.Types.Bitscore +import Data.PrimitiveArray hiding (fromList,toList) + +import Biobase.SElab.HMM.Types (HMM) +import Biobase.SElab.CM.ModelStructure + + + +-- | Extended CM information to calculate e-values + +data EValueParams = EValueParams + { _lambda :: !Double -- ^ λ>0 (lambda, slope) for exponential tails for local scores + , _tau :: !Double -- ^ τ (tau, location) for exponential tails for local scores + , _tau2 :: !Double -- ^ τ2 (tau, location again) for full histogram of all hits + , _dbSize :: !Int -- ^ database size in residues + , _numHits :: !Int -- ^ total number of non-overlapping hits + , _tailFit :: !Double -- ^ high-scoring tail fit + } + deriving (Eq,Show,Read,Generic) + +makeLenses ''EValueParams +makePrisms ''EValueParams + +instance Default EValueParams where + def = EValueParams + { _lambda = 0 + , _tau = 0 + , _tau2 = 0 + , _dbSize = 0 + , _numHits = 0 + , _tailFit = 0 + } + +instance Binary EValueParams +instance Serialize EValueParams +instance FromJSON EValueParams +instance ToJSON EValueParams +instance NFData EValueParams +instance Hashable EValueParams + + + + + + + +data EntryExit = EntryState | ExitState + deriving (Eq,Ord,Read,Show,Generic) + +{- + + + +-- | A pure getter to retrieve the last state + +sLastState :: Getter States (PInt () StateIndex) +sLastState = sStateType . to bounds . to snd +{-# Inline sLastState #-} + +-} + + + +-- | Covariance models for @Infernal@ non-coding RNA structures. + +data CM = CM + { _version :: (Text,Text) + , _name :: Text + , _accession :: Accession Rfam + , _description :: Text + , _clen :: Int + , _statesInModel :: Int + , _nodesInModel :: Int + , _w :: Int + , _referenceAnno :: Bool -- ^ have we picked up reference annotation from @GC RF@ lines in Stockholm? and integrated into match states? + , _consensusRes :: Bool -- ^ valid consensus residue annotation? + , _alignColMap :: Bool -- ^ if yes, we have map annotation in the main model annotating which multiple-alignment column a state came from + , _alph :: Text + , _date :: Text + , _commandLineLog :: Seq Text + , _pbegin :: Double + , _pend :: Double + , _wbeta :: Double + , _qdbBeta1 :: Double + , _qdbBeta2 :: Double + , _n2Omega :: Double + , _n3Omega :: Double + , _elseLF :: Double + , _nseq :: Int + , _effn :: Double + , _cksum :: Word32 + , _nullModel :: Vector Bitscore + , _ga :: Double + , _tc :: Double + , _efp7gf :: (Double,Double) + , _ecmlc :: EValueParams + , _ecmgc :: EValueParams + , _ecmli :: EValueParams + , _ecmgi :: EValueParams + , _cm :: Either FlexibleModel StaticModel +-- , _nodes :: V.Vector Node +-- , _states :: States + , _hmm :: HMM Rfam + , _unknownLines :: Seq Text + , _cmIsLocal :: Bool -- ^ @True@ if we are in local mode + } + deriving (Eq,Show,Read,Generic) + +makeLenses ''CM +makePrisms ''CM + +instance Default CM where + def = CM + { _version = ("","") + , _name = "" + , _accession = "" + , _description = "" + , _clen = 0 + , _statesInModel = 0 + , _nodesInModel = 0 + , _w = 0 + , _referenceAnno = False + , _consensusRes = False + , _alignColMap = False + , _alph = "" + , _date = "" + , _commandLineLog = def + , _pbegin = 0.05 + , _pend = 0.05 + , _wbeta = 0 + , _qdbBeta1 = 0 + , _qdbBeta2 = 0 + , _n2Omega = 0 + , _n3Omega = 0 + , _elseLF = 0 + , _nseq = 0 + , _effn = 0 + , _cksum = 0 + , _nullModel = empty + , _ga = 0 + , _tc = 0 + , _efp7gf = (0,0) + , _ecmlc = def + , _ecmgc = def + , _ecmli = def + , _ecmgi = def + , _cm = Left def +-- , _nodes = empty +-- , _states = def + , _hmm = def + , _unknownLines = def + , _cmIsLocal = False + } + +instance Binary CM +instance Serialize CM +instance FromJSON CM +instance ToJSON CM +instance NFData CM + + + +-- * Operations on CMs + +{- + +-- | Given an otherwise valid 'CM', build the efficient 'States' system. +-- +-- Normally, we will do @CM@ manipulations with the @CM@ itself, which is +-- much more highlevel, but slower due to the data structures used. +-- Building the actual 'States' structure provides a low-level +-- high-performance structure for computations. +-- +-- TODO check that @ss@ actually is sorted and dense. + +buildStatesFromCM :: CM -> States +buildStatesFromCM cm = States + { _sTransitions = fromList [ s^.transitions | s <- ss ] + , _sPairEmissions = fromAssocs (Z:.0:.A:.A) (Z:.maxState:.U:.U) def + . concatMap (\s -> [((Z:.s^.sid:.n1:.n2),e) | emitsPair (s^.sType), ((n1,n2),e) <- zip ((,) <$> acgu <*> acgu) (toList $ s^.emissions)]) $ ss + , _sSingleEmissions = fromAssocs (Z:.0:.A) (Z:.maxState:.U) def + . concatMap (\s -> [((Z:.s^.sid:.nt),e) | emitsSingle (s^.sType), (nt,e) <- zip acgu (toList $ s^.emissions)] ) $ ss + , _sStateType = fromAssocs 0 maxState (StateType $ -1) . Prelude.map ((,) <$> view sid <*> view sType) $ ss + } + where maxState = maximum $ ss^..folded.sid + ss = cm^..nodes.folded.nstates.folded + +-- | Determine if any of the children for a state is an @E@ state. + +hasEndNext :: CM -> State -> Bool +hasEndNext cm s = any (`elem` ss) kids + where kids = VG.toList . VG.map fst $ s^.transitions + ss = cm^..nodes.folded.nstates.folded.filtered ((==E) . view sType).sid + +-- | Create a local @CM@. +-- +-- A local @CM@ has transitions from the @Root 0@ states to all other @MATP +-- / MP@, @MATL / ML@, @MATR / MR@, @BIF / B@ states, except for the state +-- in @Node 1@, which already is connected. The probabilities are set to +-- @PBEGIN@ or @0.05@ in case @PBEGIN@ is not set, and divided by the +-- number of such states to move to. + +makeLocal :: CM -> CM +makeLocal cm = (addLocalEnds $ addLocalBegins cm) & cmIsLocal .~ True + +-- | Given a @CM@, add the necessary transitions to create local +-- beginnings. +-- +-- This is done by simply adding additional transitions from the @S 0@ +-- state and its companions @IL 1@ and @IR 2@. + +addLocalBegins :: CM -> CM +addLocalBegins cm = cm & nodes . vectorIx 0 . nstates . traverse . transitions %~ addbegs + where lbegs = drop 1 $ cm^..nodes.folded.(nodeMainState EntryState).sid -- the first node already has a main transition from @S 0@. + lbp = localBeginBitscore cm + addbegs :: Transitions Bitscore -> Transitions Bitscore + addbegs ts = ts VG.++ (VG.map (,lbp) $ fromList lbegs) + +-- | Calculate the actual local beginnings score. + +localBeginBitscore :: CM -> Bitscore +localBeginBitscore cm = prob2Score 1 $ cm^.pbegin / genericLength lbegs + where lbegs = drop 1 $ cm^..nodes.folded.(nodeMainState EntryState).sid + +-- | Add a single @local end@ state and transitions to this state. + +addLocalEnds :: CM -> CM +addLocalEnds cm = lendcm & states .~ (buildStatesFromCM lendcm) + where elsid = maximum (cm^..nodes.folded.nstates.folded.sid) + 1 + el = State -- the new local end state to be inserted. + { _sType = EL + , _sid = elsid + , _sParents = (-1,-1) -- TODO ??? + , _sqdb = (0,0,0,0) -- TODO ??? + , _transitions = empty + , _emissions = empty + } + ell = Node + { _nstates = VG.singleton el + , _ntype = End + , _nid = maximum (cm^..nodes.folded.nid) + 1 + , _nColL = -1 -- TODO all below + , _nColR = -1 + , _nConL = '-' + , _nConR = '-' + , _nRefL = '-' + , _nRefR = '-' + } + lep = localEndBitscore cm + addEnds :: Transitions Bitscore -> Transitions Bitscore + addEnds ts = ts `VG.snoc` (elsid,lep) + lendcm = (cm & nodes %~ (`VG.snoc` ell)) + & nodes.traverse.(nodeMainState ExitState).filtered (not . hasEndNext cm).transitions %~ addEnds + +-- | + +localEndBitscore :: CM -> Bitscore +localEndBitscore cm = prob2Score 1 $ cm^.pend / (genericLength $ cm^..nodes.folded.(nodeMainState ExitState).filtered (not . hasEndNext cm)) + +-- | + +makeGlobal :: CM -> CM +makeGlobal = undefined + +-} + diff --git a/Biobase/SElab/Common/Parser.hs b/Biobase/SElab/Common/Parser.hs new file mode 100644 --- /dev/null +++ b/Biobase/SElab/Common/Parser.hs @@ -0,0 +1,38 @@ +{-# LANGUAGE NoMonomorphismRestriction #-} +{-# LANGUAGE OverloadedStrings #-} + +-- | Common parser helpers. The @-S@ versions convert from @ByteString@ to +-- @Text@ here. + +module Biobase.SElab.Common.Parser where + +import Control.Applicative +import Data.Attoparsec.ByteString.Char8 +import Data.Char (isAlpha,isDigit) +import Data.Char.Util +import Data.Text.Encoding (decodeUtf8) +import Prelude hiding (takeWhile) + + + +-- * Helper functions + +ssN = skipSpace *> decimal +ssN_ = skipSpace *> ((-1) <$ "-" <|> decimal) +ssZ = skipSpace *> signed decimal +ssQ = skipSpace *> rational +ssD = skipSpace *> double "Double" +ssD' = skipSpace *> ((-999999) <$ "*" <|> double) "Double ('*' aware)" +ssS = skipSpace *> takeTill (\c' -> let c = c2w8 c' in isEndOfLine c || isHorizontalSpace c) +ssT = decodeUtf8 <$ skipSpace <*> takeTill (\c' -> let c = c2w8 c' in isEndOfLine c || isHorizontalSpace c) +ssC = skipSpace *> anyChar + +eolD = skipSpace *> double <* endOfLine "eolD" +eolN = skipSpace *> decimal <* endOfLine "eolN" +eolR = skipSpace *> rational <* endOfLine "eolR" +eolZ = skipSpace *> signed decimal <* endOfLine "eolZ" +eolC = skipSpace *> satisfy (not . isSpace) <* endOfLine "eolC" +eolS = takeWhile (isHorizontalSpace . c2w8) *> takeTill (isEndOfLine . c2w8) <* endOfLine "eolS" +eolT = decodeUtf8 <$ takeWhile (isHorizontalSpace . c2w8) <*> takeTill (isEndOfLine . c2w8) <* endOfLine "eolS" +eolB = skipSpace *> (True <$ "yes" <|> False <$ "no") <* endOfLine "eolB" + diff --git a/Biobase/SElab/HMM.hs b/Biobase/SElab/HMM.hs --- a/Biobase/SElab/HMM.hs +++ b/Biobase/SElab/HMM.hs @@ -1,74 +1,11 @@ -{-# LANGUAGE TemplateHaskell #-} -{-# LANGUAGE EmptyDataDecls #-} --- | HMMER3 HMMs. Since we do not understand HMMER3 HMMs yet, this is actually --- just a small ``throw-away'' parser to successfully parse Infernal 1.1 CMs. --- The next version should have a real working parser. --- --- TODO in the future, we should split parsing into just grabbing lines between --- HMMER and "//" and handling in-between. We need extraction of individual --- models and similar fun. - -module Biobase.SElab.HMM where - -import Data.ByteString.Char8 as BS -import Control.Lens - -import Biobase.SElab.Types - - - -data HMM - -data Alphabet - = Amino - | DNA - | RNA - | Coins - | Dice - | Custom - deriving (Eq,Show,Read) - --- | Negated natural logarithm of probability. --- --- TODO put into types stuff - -newtype NegLogProb = NLP Double - deriving (Show,Read) - --- | The nodes in an HMM. Starting with Node "0" for BEGIN. - -data Node = Node - { _nid :: Int - , _matchE :: [NegLogProb] -- [] for BEGIN - , _insertE :: [NegLogProb] -- insertions - , _trans :: [NegLogProb] -- transitions: B->M1 B->I0 B->D1 I0->M1 I0->I0 0.0 * ||| Mk->Mk+1 Mk->Ik Mk->Dk+1 Ik->Mk+1 Ik->Ik Dk->Mk+1 Dk->Dk+1 - } - deriving (Show,Read) - -makeLenses ''Node +-- | --- | The HMM3 data structure in ``slow mode''. --- --- TODO shouldn't this be "Identification Pfam" ? --- --- TODO maybe redo the whole "idd" idea and just keep the string? +module Biobase.SElab.HMM + ( module Biobase.SElab.HMM.Types + , module Biobase.SElab.HMM.Import + ) where -data HMM3 = HMM3 - { _version :: (ByteString,ByteString) - , _idd :: Identification HMM - , _acc :: Maybe (Accession HMM) - , _description :: Maybe ByteString - , _leng :: Int -- mandatory >0 count of match states - , _alph :: Alphabet - , _rf :: Bool - , _cs :: Bool - , _alignMap :: Bool - , _date :: ByteString - , _symAlph :: [ByteString] - , _transHeaders :: [ByteString] - , _compo :: [NegLogProb] - , _nodes :: [Node] - } deriving (Show,Read) +import Biobase.SElab.HMM.Import (hmmFromFile, parseHMM) +import Biobase.SElab.HMM.Types -makeLenses ''HMM3 diff --git a/Biobase/SElab/HMM/Import.hs b/Biobase/SElab/HMM/Import.hs --- a/Biobase/SElab/HMM/Import.hs +++ b/Biobase/SElab/HMM/Import.hs @@ -1,178 +1,157 @@ -{-# LANGUAGE DoAndIfThenElse #-} -{-# LANGUAGE OverloadedStrings #-} -{-# LANGUAGE NoMonomorphismRestriction #-} -- | Import HMMER3 HMM models. module Biobase.SElab.HMM.Import where -import Control.Monad.IO.Class (liftIO, MonadIO) -import Data.ByteString.Char8 as BS -import Data.ByteString.Lex.Double as BS -import Data.Conduit as C -import Data.Conduit.Binary as CB -import Data.Conduit.List as CL -import Control.Monad (unless) -import Prelude as P -import Control.Arrow +import Control.Applicative ( (<|>), pure, (<$>), (<$), (<*>), (*>), (<*) ) +import Control.Lens hiding ((|>)) +import Control.Monad +import Control.Monad.IO.Class (MonadIO) +import Data.Attoparsec.ByteString (count,many1,(),manyTill,option) +import Data.ByteString.Char8 (ByteString,unpack) +import Data.Char (isSpace,isAlpha,isDigit) +import Data.Char.Util +import Data.Default +import Data.Sequence ((|>)) +import Data.Text.Encoding (decodeUtf8) +import Data.Text (Text) +import Data.Vector.Unboxed (fromList) +import Debug.Trace +import qualified Data.Attoparsec.ByteString as AB +import qualified Data.Attoparsec.ByteString.Char8 as ABC +import qualified Data.ByteString.Char8 as BSC +import qualified Data.List as L import qualified Data.Map as M -import Data.Char (toLower) +import qualified Data.Text as T +import qualified Data.Vector.Unboxed as VU +import System.FilePath (takeExtension) +import Control.DeepSeq (($!!)) -import Biobase.SElab.HMM -import Biobase.SElab.Types +import Biobase.Primary +import Biobase.Types.Accession (Accession(..)) +import Biobase.Types.Bitscore +import Data.PrimitiveArray as PA hiding (map) +import Biobase.SElab.Common.Parser +import Biobase.SElab.HMM.Types --- * Different HMMer parsers --- ** HMMER3 / b --- | --- --- TODO not everything is currently being parsed. Notably the rf,cs,alignmap --- annotations. - --- parseHMM3 :: (Monad m, MonadIO m) => Conduit ByteString m HMM3 -parseHMM3 = go where - go = do - hdr' <- CL.head - unless (legalHMM hdr') . error $ "no legal HMM at header: " ++ show hdr' - let Just hdr = hdr' - hs <- headerMap `fmap` headerLines - (sas,ths) <- sathLines - let asize = P.length sas - c <- compoLine - n0 <- parseBegin asize - ns <- parseNodes asize - Just "//" <- CL.head - return $ HMM3 - { _version = second (BS.dropWhile (==' ')) . BS.span (/=' ') $ hdr - , _idd = IDD $ hs M.! "NAME" - , _acc = fmap (ACC . readBS) $ "AC" `M.lookup` hs - , _description = "DESC" `M.lookup` hs - , _leng = readBS $ hs M.! "LENG" - , _alph = readAlph $ hs M.! "ALPH" - , _rf = readBoolean $ M.findWithDefault "no" "RF" hs - , _cs = readBoolean $ M.findWithDefault "no" "CS" hs - , _alignMap = readBoolean $ M.findWithDefault "no" "MAP" hs - , _date = M.findWithDefault "" "DATE" hs - , _symAlph = sas - , _transHeaders = ths - , _compo = c - , _nodes = n0:ns - } - --- | Check, if we have a legal HMMER3 model. - -legalHMM :: Maybe ByteString -> Bool -legalHMM (Just s) - | w == "HMMER3/f" = True - | w == "HMMER3/b" = True - where (w:_) = BS.words s -legalHMM _ = False - - - --- * Helper functions - --- | Read boolean flags. - -readBoolean = f . BS.map toLower where - f "no" = False - f "yes" = True - f x = error $ "unknown boolean: " ++ show x - --- | Determine which alphabet is in use by the HMM. - -readAlph = f . BS.map toLower where - f "dna" = DNA - f "rna" = RNA - f "coins" = Coins - f "dice" = Dice - f "amino" = Amino - f "custom" = Custom - f a = error $ "unknown alph: " ++ show a - --- | Read from a bytestring into a structure. - -readBS = read . BS.unpack +-- | Simple convenience function for parsing HMM's without a lot of +-- fancyness. --- | create associative map of the key/value data. +hmmFromFile :: FilePath -> IO [HMM ()] +hmmFromFile fp = do + bs <- BSC.readFile fp + case AB.parseOnly (many1 parseHMM <* AB.endOfInput) bs of + Left err -> error err + Right xs -> return xs -headerMap xs = M.fromList . P.map f $ xs where - f = second (BS.dropWhile (==' ')) . BS.span (/=' ') +-- | +-- +-- NOTE the idea of filling with @999999@ is that if we run the HMM, then any +-- score bugs will yield weird results that show up immediately. --- | Parse the two beginning lines. +parseHMM :: ABC.Parser (HMM xfam) +parseHMM = do + pre <- parsePreHMM + bdy <- parseHMMBody pre + return bdy -parseBegin asize = do - Just i' <- CL.head - Just t' <- CL.head - return $ Node - 0 - [] - (P.map (readNLP . BS.unpack) $ BS.words i') - (P.map (readNLP . BS.unpack) $ BS.words t') +-- | Parse the header of an HMM, and return the partially filled HMM and +-- a ByteString with the non-parsed remainder. --- | Parse all individual nodes, except the first one, which uses 'parseBegin'. +parsePreHMM :: ABC.Parser (HMM xfam) -- (HMM xfam, Text) +parsePreHMM = do + v <- acceptedVersion + let hmm' = version .~ v $ def + ls <- hmmHeader `manyTill` "HMM" + let hmm = L.foldl' (\a l -> l a) hmm' ls + eolS + eolS + -- remainder <- ABC.takeText + return hmm -- (hmm, remainder) -parseNodes asize = go [] where - go xs = do - p <- CL.peek - case p of - (Just "//") -> return $ P.reverse xs - _ -> do Just m' <- CL.head - Just i' <- CL.head - Just t' <- CL.head - let (nid:m) = BS.words m' - let n = Node - (read . BS.unpack $ nid) - (P.map (readNLP . BS.unpack) $ P.take asize m) - (P.map (readNLP . BS.unpack) $ BS.words i') - (P.map (readNLP . BS.unpack) $ BS.words t') - go (n:xs) +parseHMMBody :: HMM xfam -> ABC.Parser (HMM xfam) +parseHMMBody hmm = do + l <- component0 + ls <- (component (length $ l^._2)) `manyTill` "//" + ABC.skipSpace + return + $!! set matchScores (PA.fromAssocs (Z:.0:.Letter 0) (Z:.(PInt $ length ls):.(Letter . subtract 1 . length $ l^._2)) 999999 + [((Z:.s:.k),Bitscore v) | (s,vs) <- zip [0..] (l^._2:map (view (_2._1)) ls), (k,v) <- zip [Letter 0 ..] vs ]) + $ set insertScores (PA.fromAssocs (Z:.0:.Letter 0) (Z:.(PInt $ length ls):.(Letter . subtract 1 . length $ l^._3)) 999999 + [((Z:.s:.k),Bitscore v) | (s,vs) <- zip [0..] (l^._3:map (view _3 ) ls), (k,v) <- zip [Letter 0 ..] vs ]) + $ set transitionScores (PA.fromAssocs (Z:.0:.Letter 0) (Z:.(PInt $ length ls):.(Letter . subtract 1 . length $ l^._4)) 999999 + [((Z:.s:.k),Bitscore v) | (s,vs) <- zip [0..] (l^._4:map (view _4 ) ls), (k,v) <- zip [Letter 0 ..] vs ]) + $ hmm --- | Read a HMMER negated log-probability. +acceptedVersion :: ABC.Parser (Text,Text) +acceptedVersion = (,) <$> (decodeUtf8 <$> vOk) <* ABC.skipSpace <*> eolT "accepted Version" where + vOk = "HMMER3/b" <|> "HMMER3/f" -readNLP :: String -> NegLogProb -readNLP = go where - go "*" = NLP $ 1/0 - go xs = NLP . read $ xs +hmmHeader :: ABC.Parser (HMM xfam -> HMM xfam) +hmmHeader = ABC.choice + [ set name <$ "NAME" <*> eolT "name" + , set accession . Accession <$ "ACC" <*> eolT "hmmAccession" + , set description <$ "DESC" <*> eolT "description" + , set modelLength <$ "LENG" <*> eolN "leng" + , set maxInstanceLen . Just <$ "MAXL" <*> eolN "maxl" + , set alphabet <$ "ALPH" <*> eolT "alph" + , set referenceAnno <$ "RF" <*> eolB "rf" + , set consensusStruc <$ "CS" <*> eolB "cs" + , set consensusRes <$ "CONS" <*> eolB "cons" + , set alignColMap <$ "MAP" <*> eolB "map" + , set modelMask <$ "MM" <*> eolB "mm" + , set date <$ "DATE" <*> eolT "date" + , set nseq . Just <$ "NSEQ" <*> eolN "nseq" + , set effnseq . Just <$ "EFFN" <*> eolD "effn" + , set chksum . Just <$ "CKSUM" <*> eolN "cksum" + , (\l r -> set gatheringTh (Just (l,r))) <$ "GA" <*> ssD <*> ssD <* eolS + , (\l r -> set trustedCutoff (Just (l,r))) <$ "TC" <*> ssD <*> ssD <* eolS + , (\l r -> set noiseCutoff (Just (l,r))) <$ "NC" <*> ssD <*> ssD <* eolS + , (\l r -> set msv (Just (l,r))) <$ "STATS LOCAL MSV" <*> ssD <*> ssD <* eolS + , (\l r -> set viterbi (Just (l,r))) <$ "STATS LOCAL VITERBI" <*> ssD <*> ssD <* eolS + , (\l r -> set forward (Just (l,r))) <$ "STATS LOCAL FORWARD" <*> ssD <*> ssD <* eolS + , (\s -> over commandLineLog (|> decodeUtf8 s)) <$ "COM" <*> eolS "com" + , (\x -> over unknownLines (|> decodeUtf8 x)) <$> ABC.takeWhile1 (/='\n') <* ABC.take 1 + ] "hmmHeader" --- | Read the optional COMPO line. +-- | TODO -compoLine = do - Just p <- CL.peek - case (BS.words p) of - ("COMPO":xs) -> CL.head >>= \_ -> return $ P.map (NLP . read . BS.unpack) xs - _ -> return [] +component0 :: ABC.Parser Component0 +component0 = (,,,) <$> ident <*> matches <*> inserts <*> moves "COMPO/0" where + ident = ABC.skipSpace *> (0 <$ "COMPO") "ident" -- optional + matches = manyTill ssD ABC.endOfLine "matches" -- optional + inserts = manyTill ssD ABC.endOfLine "inserts" + moves = count 7 ssD' <* ABC.endOfLine "moves" --- | Read the alphabet and transition lines. +-- | Parse components. Matches come with annotations. These depend on the specific model. -sathLines = do - Just sa' <- CL.head - Just th' <- CL.head - let (sa:sas) = BS.words sa' - let ths = BS.words th' - if sa == "HMM" - then return (sas,ths) - else error $ "NOT THE HMM symalph lines: " ++ show (sa:sas,ths) +component :: Int -> ABC.Parser Component +component k = (,,,) <$> ident <*> (matches "matches") <*> inserts <*> moves "component" where + ident = ABC.skipSpace *> (error "COMPO parsed in component" <$ "COMPO" <|> ABC.decimal) "ident" + matches = matchHMM <|> matchSubHMM "matchHMM-CM" + matchHMM = (,,,,,) <$> count k ssD <*> melMAP <*> pure ' ' <*> melRF <*> melCS <*> pure ' ' <* (ABC.endOfLine "eol") "matchHMM" + matchSubHMM = (,,,,,) <$> count k ssD <*> melMAP <*> melCONS <*> melRF <*> melCS <*> melStruc <* (ABC.endOfLine "eol") "matchSubHMM" -- SubHMM of a CM, not a CM! + inserts = count k ssD <* ABC.endOfLine "inserts" + moves = count 7 ssD' <* ABC.endOfLine "moves" + melMAP = skipHorizSpace *> (ABC.decimal <|> (0 <$ "-")) "melMAP" + melCONS = skipHorizSpace *> ABC.anyChar "melCONS" + melRF = skipHorizSpace *> ABC.anyChar "melRF" + melCS = skipHorizSpace *> ABC.anyChar "melCS" + melStruc = skipHorizSpace *> ABC.anyChar "melSTRUC" -- not defined in the Userguide! + skipHorizSpace = ABC.skipWhile (ABC.isHorizontalSpace . c2w8) --- | All the header lines until we see "HMM". +type Component0 = (Int, [Double] , [Double], [Double]) -headerLines = go [] where - go xs = do - p <- CL.peek - case p of - (Just x) | "HMM" `BS.isPrefixOf` x -> return $ P.reverse xs - | otherwise -> CL.drop 1 >> go (x:xs) - Nothing -> error $ "no more lines after: " ++ show (P.reverse xs) +-- | A Component line. Index (starting with 1, zero is COMPO). Then comes the +-- match line with the scores, a MAP annotation, consensus residue, reference +-- annotation, and consensus structure. HMMer HMMs don't have a consensus +-- structure. +type Component = (Int, ([Double],Int,Char,Char,Char,Char), [Double], [Double]) --- | Simple test for the HMMer parser. - -test :: IO () -test = do - xs <- runResourceT $ sourceFile "test.hmm" =$= CB.lines $= CL.sequence parseHMM3 $$ consume -- sinkHandle stdout - print xs diff --git a/Biobase/SElab/HMM/Types.hs b/Biobase/SElab/HMM/Types.hs new file mode 100644 --- /dev/null +++ b/Biobase/SElab/HMM/Types.hs @@ -0,0 +1,111 @@ + +-- | Efficient encoding of Hidden-Markov models as used by @HMMER 3@ and +-- @Infernal 1.1@. + +module Biobase.SElab.HMM.Types where + +import Control.DeepSeq +import Control.Lens +import Data.Aeson (FromJSON,ToJSON) +import Data.Binary (Binary) +import Data.Default +import Data.PrimitiveArray +import Data.Sequence (Seq) +import Data.Serialize (Serialize) +import Data.Text (Text) +import Data.Vector.Unboxed (Vector,empty) +import Data.Word (Word32) +import GHC.Generics (Generic) + +import Biobase.Primary +import Biobase.Types.Accession +import Biobase.Types.Bitscore +import Data.PrimitiveArray + + + +-- | Which node in the HMM are we in? + +data NodeIndex + +-- | An efficient encoding of Infernal HMM models. With @xfam@ phantom type +-- that will pin the type to @Pfam@ or @Rfam@ (or maybe others later). +-- +-- TODO Parsing has only been tested for HMMER3 and Infernal 1.1 + +data HMM xfam = HMM + { _version :: (Text,Text) -- ^ magic string aka @HMMER3/f@ (HMMer) or @HMMER3/i@ (Infernal), followed by the bracketed version info + , _name :: Text -- ^ the name of this HMM; for Rfam HMMs, same as CM + , _accession :: Accession xfam -- ^ a la @PF01234@ or @RF01234@. + , _description :: Text -- ^ one-line free text description + , _maxInstanceLen :: Maybe Int -- ^ upper on length at which an instance of the model is expected to be found + , _alphabet :: Text -- ^ alphabet type, case insensitive. We are interested in @amino@, @DNA@, @RNA@, but there are others as well. + , _modelLength :: Int -- ^ number of match states in the model + , _referenceAnno :: Bool -- ^ have we picked up reference annotation from @GC RF@ lines in Stockholm? and integrated into match states? + , _consensusRes :: Bool -- ^ valid consensus residue annotation? + , _consensusStruc :: Bool -- ^ picked up consensus annotation @@SS_cons@? + , _alignColMap :: Bool -- ^ if yes, we have map annotation in the main model annotating which multiple-alignment column a state came from + , _modelMask :: Bool -- ^ if yes, a model mask is active. Annotates columns which are set to background frequency instead of observed model frequency + , _gatheringTh :: Maybe (Double,Double) -- ^ gathering thresholds @ga1@ and @ga2@ + , _trustedCutoff :: Maybe (Double,Double) -- ^ trusted cutoffs @tc1@ and @tc2@ + , _noiseCutoff :: Maybe (Double,Double) -- ^ noise cutoffs @nc1@ and @nc2@ + , _date :: Text -- ^ model creation date + , _commandLineLog :: Seq Text -- ^ commands to build this model + , _nseq :: Maybe Int -- ^ number of sequences in multiple alignment + , _effnseq :: Maybe Double -- ^ effective number of sequences (after weighting) + , _chksum :: Maybe Word32 -- ^ checksum (TODO: replace Word32 with actual checksum newtype) + , _msv :: Maybe (Double,Double) -- ^ μ (mu) and λ (lambda) for gumbel distribution + , _viterbi :: Maybe (Double,Double) -- ^ μ (mu) and λ (lambda) for gumbel distribution + , _forward :: Maybe (Double,Double) -- ^ t (tau) and l (lambda) for exponential tails + , _matchMap :: Vector Int -- ^ match node alignment index + , _matchRef :: Vector Char -- ^ match node reference annotation + , _matchCons :: Vector Char -- ^ match node consensus annotation + , _matchScores :: Unboxed (Z:.PInt I NodeIndex:.Letter Unknown) Bitscore -- ^ + , _insertScores :: Unboxed (Z:.PInt I NodeIndex:.Letter Unknown) Bitscore -- ^ + , _transitionScores :: Unboxed (Z:.PInt I NodeIndex:.Letter Unknown) Bitscore -- ^ + , _unknownLines :: Seq Text -- ^ filled with header lines that can not be parsed + } deriving (Eq,Show,Read,Generic) + +makeLenses ''HMM +makePrisms ''HMM + +instance Default (HMM xfam) where + def = HMM + { _version = ("","") + , _name = "" + , _accession = "" + , _description = "" + , _maxInstanceLen = Nothing + , _alphabet = "" + , _modelLength = -1 + , _referenceAnno = False + , _consensusRes = False + , _consensusStruc = False + , _alignColMap = False + , _modelMask = False + , _gatheringTh = Nothing + , _trustedCutoff = Nothing + , _noiseCutoff = Nothing + , _date = "" + , _commandLineLog = def + , _nseq = Nothing + , _effnseq = Nothing + , _chksum = Nothing + , _msv = Nothing + , _viterbi = Nothing + , _forward = Nothing + , _matchMap = empty + , _matchRef = empty + , _matchCons = empty + , _matchScores = fromAssocs (Z:.0:.Letter 0) (Z:.0:.Letter 0) def [] + , _insertScores = fromAssocs (Z:.0:.Letter 0) (Z:.0:.Letter 0) def [] + , _transitionScores = fromAssocs (Z:.0:.Letter 0) (Z:.0:.Letter 0) def [] + , _unknownLines = def + } + +instance Binary (HMM xfam) +instance Serialize (HMM xfam) +instance FromJSON (HMM xfam) +instance ToJSON (HMM xfam) +instance NFData (HMM xfam) + diff --git a/Biobase/SElab/Model.hs b/Biobase/SElab/Model.hs new file mode 100644 --- /dev/null +++ b/Biobase/SElab/Model.hs @@ -0,0 +1,8 @@ + +module Biobase.SElab.Model + ( fromFile + , keepAllModels + ) where + +import Biobase.SElab.Model.Import + diff --git a/Biobase/SElab/Model/Import.hs b/Biobase/SElab/Model/Import.hs new file mode 100644 --- /dev/null +++ b/Biobase/SElab/Model/Import.hs @@ -0,0 +1,216 @@ + +-- | Import SElab HMM and CM models, and combine these together into full +-- CM models. +-- +-- TODO A pair @Right CM +++ Left HMM@ should become a @Right CM@, a lonely +-- @Left HMM@ should become a @Left HMM@, a lonely @Right CM@ should become +-- a @Right CM@. + +module Biobase.SElab.Model.Import where + +import Control.Applicative +import Control.DeepSeq +import Control.Lens (set,over,(^.),zoom) +import Control.Monad (liftM2) +import Control.Monad.Trans.Class (lift) +import Control.Monad.Trans.Writer.Strict +import Control.Monad (void,unless) +import Control.Monad (when,replicateM) +import Control.Parallel.Strategies (using,parList,rdeepseq,parMap) +import Data.ByteString (ByteString) +import Data.Maybe (catMaybes) +import Data.Monoid +import Data.String (IsString) +import Data.Text (Text) +import Data.Text.Encoding (decodeUtf8) +import Debug.Trace +import qualified Data.Attoparsec.ByteString.Char8 as ABC +import qualified Data.Attoparsec.Text as AT +import qualified Data.ByteString.Char8 as BS +import qualified Data.List as L +import qualified Data.Map.Strict as M +import qualified Data.Text as T +import qualified Data.Text.IO as T +import qualified Pipes as P +import qualified Pipes.Attoparsec as PA +import qualified Pipes.ByteString as PB +import qualified Pipes.GZip as PG +import qualified Pipes.Parse as PP +import qualified Pipes.Prelude as P +import qualified Pipes.Safe.Prelude as PSP +import System.FilePath (takeExtension) +import System.IO (stdin,withFile,IOMode(..)) + +import Biobase.Types.Accession +import Pipes.Split.ByteString + +import Biobase.SElab.CM.Import as CM +import Biobase.SElab.CM.Types as CM +import Biobase.SElab.HMM.Import as HMM +import Biobase.SElab.HMM.Types as HMM +import Biobase.SElab.Model.Types + + + +-- | Filter a model after the header, not the body, has been parsed. + +type PreFilterFun = Text -> Accession () -> Either (HMM ()) CM -> Bool + +-- | Filter a model after the full model has been parsed. This is the same +-- type as @PreFilterFun@. + +type PostFilterFun = Text -> Accession () -> Either (HMM ()) CM -> Bool + +-- TODO use a builder? + +newtype Log = Log { getLog :: Text } + deriving (Monoid,IsString) + +-- | The type of logger we use + +type Logger m = WriterT Log m + + +-- | Combine CMs with their HMMs. Assumes that each CM is followed by its +-- HMM. + +attachHMMs :: (Monad m) => PP.Producer Model (Logger m) r -> PP.Producer CM (Logger m) ((), PP.Producer Model (Logger m) r) +attachHMMs = PP.parsed go where + go :: (Monad m) => PP.Parser Model (Logger m) (Either () CM) + go = do + mcm <- PP.draw + case mcm of + Nothing -> return $ Left () + Just (Left hmm) -> do + lift . tell . Log $ "HMM: " <> (hmm^.HMM.name) <> " is an orphan\n" + go + Just (Right cm) -> do + mhm <- PP.draw + case mhm of + Nothing -> do + lift . tell . Log $ "CM: " <> (cm^.CM.name) <> " has no attached HMM and the stream is finished\n" + return $ Right cm + -- TODO actually check that these belong together + Just (Left hm) | (cm^.CM.name) == (hm^.HMM.name) -> do + -- cm and hmm belong together + return . Right $ set hmm (over HMM.accession retagAccession hm) cm + -- The HMM doesn't belong to our CM + Just (Left hm) -> do + lift . tell . Log $ "CM: " <> (cm^.CM.name) <> " and HMM: " <> (hm^.HMM.name) <> " do not belong together, dropping the HMM from the stream\n" + return . Right $ cm + Just (Right dup) -> do + lift . tell . Log $ "CM: " <> (cm^.CM.name) <> " has no attached HMM\n" + PP.unDraw $ Right dup + go + + + +-- | Parses @HMM@ and @CM@ models from Rfam. The filtering function takes +-- the model name and accession and allows for premature termination of the +-- parsing of the current model. + +parseSelectively :: (Monad m) + => PreFilterFun + -- ^ filter function for premodels + -> PostFilterFun + -- ^ filter function for full models + -> PP.Producer ByteString (Logger m) r +-- -> PP.Producer (Maybe (Either (HMM xfam) CM)) (Logger m) (r, PP.Producer ByteString (Logger m) r) + -> PP.Producer Model (Logger m) ((), PP.Producer ByteString (Logger m) r) +parseSelectively preFltr postFltr p + = PP.parsed go p + P.>-> P.concat + P.>-> P.filter (\mdl -> postFltr (mdl^.modelName) (mdl^.modelAccession) mdl) + where + -- | Parse either a CM or a HMM ... + go = do + p <- zoom (splitKeepEnd "//\n") parseMdl + -- TODO can be simplified now + case p of + Left () -> return $ Left () + Right x -> return $ Right x + -- parse models. + -- TODO @Either ()@ should become @Either (Maybe Error)@ and only @Left + -- Nothing@ will be error-free stop. + handleError err = do + da <- PP.drawAll + lift . tell . Log $ "could not parse:\n" + lift . tell . Log $ T.pack $ show err + lift . tell . Log . decodeUtf8 $ BS.concat da + lift . tell . Log $ "\n" + return $ Left () + parseMdl :: Monad m => PP.StateT (PP.Producer ByteString (Logger m) x) (Logger m) (Either () (Maybe (Either (HMM ()) CM))) + parseMdl = do + -- if @pre@ is Nothing, the underlying producer is exhausted. + -- if @pre@ is @Just $ Left x@, then we have a parse error. + -- if @pre@ is @Just $ Right y@, then we have a successful parse. In + -- this case, @y@ is either a @Left hmm@ or a @Right cm@. + pre <- PA.parse $ (Left <$> parsePreHMM) <|> (Right <$> parsePreCM) + case pre of + -- we have nothing left to parse and indicate this now + Nothing -> return $ Left () + Just (Left err) -> handleError err + Just (Right mdl) -> if preFltr (mdl^.modelName) (mdl^.modelAccession) mdl + then do + case mdl of + Left hmm -> do + h <- PA.parse $ parseHMMBody hmm + case h of + Nothing -> handleError "premature end of parsing in hmm body" + Just (Left err) -> handleError err + Just (Right hh) -> do + da <- PP.drawAll + -- TODO check if @da@ is empty? + return $ Right $ Just $ Left hh + Right cm -> do + c <- PA.parse $ parseCMBody cm + case c of + Nothing -> handleError "premature end of parsing in cm body" + Just (Left err) -> handleError err + Just (Right d) -> do + da <- PP.drawAll + return . Right . Just $ Right d + else PP.skipAll >> (return . Right $ Nothing) + + + +-- | Keep all models + +keepAllModels _ _ _ = True + + + +-- | Load a number of models from file. Including pre- and full-model +-- filtering. + +fromFile + :: PreFilterFun + -- ^ filter premodels before they are fully parsed. Full parsing is + -- costly. Use @\name acc hmmOrcm -> True@ if all models should be + -- loaded. + -> PostFilterFun + -- ^ Filter full models before they are combined into the CM-HMM pair. + -> Bool + -- ^ If true, than any error during parsing means termination of the + -- program. + -> FilePath + -- ^ input file name. Can be @-@ for stdin. If a file and the file ends + -- with @.gz@, the file is uncompressed on the ly. + -> IO [CM] +fromFile preFltr postFltr stopOnError fp + | fp == "-" = parse (PB.fromHandle stdin) + | takeExtension fp == ".gz" = withFile fp ReadMode $ \hdl -> parse (PG.decompress $ PB.fromHandle hdl) + | otherwise = withFile fp ReadMode $ \hdl -> parse (PB.fromHandle hdl) + where + parse source = do + ((xs,((),rmdr)),log) <- runWriterT . P.toListM' $ attachHMMs $ parseSelectively preFltr postFltr source + -- TODO log should be empty + -- TODO rmdr should be empty + let Log l = log + unless (T.null l) $ do + T.putStrLn l + T.putStrLn "There have been errors parsing the models!" + when stopOnError $ do + error "stopping here!" + return xs + diff --git a/Biobase/SElab/Model/Types.hs b/Biobase/SElab/Model/Types.hs new file mode 100644 --- /dev/null +++ b/Biobase/SElab/Model/Types.hs @@ -0,0 +1,36 @@ + +-- | Generic model wrapper. + +module Biobase.SElab.Model.Types where + +import Control.Lens +import Data.Text (Text) +import GHC.Generics (Generic) + +import Biobase.Types.Accession + +import Biobase.SElab.CM.Types (CM) +import Biobase.SElab.HMM.Types (HMM) +import qualified Biobase.SElab.CM.Types as CM +import qualified Biobase.SElab.HMM.Types as HMM + + + +-- | Generic model wrapper. + +type Model = Either (HMM ()) CM + +-- | A getter for the name of the model. Be it CM or HMM. + +modelName :: Getter Model Text +modelName = to f + where f m = case m of + Left hmm -> HMM._name hmm + Right cm -> CM._name cm + +modelAccession :: Getter Model (Accession ()) +modelAccession = to f + where f m = case m of + Left hmm -> retagAccession $ HMM._accession hmm + Right cm -> retagAccession $ CM._accession cm + diff --git a/Biobase/SElab/RfamNames.hs b/Biobase/SElab/RfamNames.hs deleted file mode 100644 --- a/Biobase/SElab/RfamNames.hs +++ /dev/null @@ -1,22 +0,0 @@ -{-# LANGUAGE TemplateHaskell #-} - --- | The database of Rfam "names". For each model, we get to know which --- sequences it is built of, what the AC of the species is, and its name (or --- ID). - -module Biobase.SElab.RfamNames where - -import Control.Lens - -import Biobase.SElab.Types - - -data ModelNames = ModelNames - { _modelAC :: !(Accession Rfam) - , _modelID :: !(Identification Rfam) - -- TODO this would have been the sequence info - , _speciesAC :: Maybe (Accession Species) - , _speciesID :: Maybe (Identification Species) - } deriving (Show) - -makeLenses ''ModelNames diff --git a/Biobase/SElab/RfamNames/Import.hs b/Biobase/SElab/RfamNames/Import.hs deleted file mode 100644 --- a/Biobase/SElab/RfamNames/Import.hs +++ /dev/null @@ -1,59 +0,0 @@ -{-# LANGUAGE OverloadedStrings #-} -{-# LANGUAGE BangPatterns #-} -{-# LANGUAGE NoMonomorphismRestriction #-} - -module Biobase.SElab.RfamNames.Import where - -import Control.Applicative -import Control.Lens -import Data.Attoparsec as A hiding (parse) -import Data.Attoparsec.Char8 as A8 hiding (parse) -import Data.Conduit as C -import Data.Conduit.Attoparsec -import Data.Conduit.Binary as CB -import Data.Conduit.List as CL -import Data.Conduit.Util as C -import Data.Either.Unwrap as E -import Data.Map (Map) -import qualified Data.ByteString.Char8 as BS -import qualified Data.Map as M - -import Biobase.SElab.RfamNames -import Biobase.SElab.Types - - - -parse = CB.lines - =$ CL.map (parseOnly mkRfamName) - =$ CL.filter isRight - =$ CL.map fromRight - =$ C.zipSinks mapIdRfamNames mapAcRfamNames -{-# INLINE parse #-} - -mkRfamName = f <$> rfamAC <* char ';' <*> rfamID <* char ';' <*> seqident <* spaces <*> specAC <* char ':' <*> specID where - f rfac rfid sid spac spid = ModelNames rfac rfid spac spid - rfamAC = ACC <$ string "RF" <*> decimal - rfamID = IDD <$> A8.takeTill (==';') - seqident = A8.takeTill isSpace - specAC = (fmap (ACC . read . BS.unpack) . maybeBS) <$> A8.takeTill (==':') - specID = (fmap IDD . maybeBS) <$> takeByteString - spaces = many1 space - maybeBS s - | BS.null s = Nothing - | otherwise = Just s -{-# INLINE mkRfamName #-} - -mapIdRfamNames = CL.fold f M.empty where - f !mp x = M.insertWith' (++) (x ^. modelID) [x] mp -{-# INLINE mapIdRfamNames #-} - -mapAcRfamNames = CL.fold f M.empty where - f !mp x = M.insertWith' (++) (x ^. modelAC) [x] mp -{-# INLINE mapAcRfamNames #-} - -fromFile :: String -> IO ( Map (Identification Rfam) [ModelNames] - , Map (Accession Rfam) [ModelNames] - ) -fromFile fname = do - runResourceT $ CB.sourceFile fname $$ parse -{-# NOINLINE fromFile #-} diff --git a/Biobase/SElab/Taxonomy.hs b/Biobase/SElab/Taxonomy.hs --- a/Biobase/SElab/Taxonomy.hs +++ b/Biobase/SElab/Taxonomy.hs @@ -3,29 +3,21 @@ -- | Infernal contains a taxonomy database. This is a simple module reflecting -- said database. +-- +-- See Task # 1effbfec-7f3b-4530-a8a7-b82fef1a0c12 -module Biobase.SElab.Taxonomy where +module Biobase.SElab.Taxonomy + ( module Biobase.Types.Taxonomy + , module Biobase.SElab.Taxonomy.Import + ) where -import Control.Lens -import Data.Char (toLower) -import qualified Data.ByteString.Char8 as BS +import Biobase.Types.Taxonomy -import Biobase.SElab.Types +import Biobase.SElab.Taxonomy.Import --- | For each species, we store the name and a classification list from most --- general (head) to most specific (last). The database comes with the NCBI --- taxon identifier (taxid). - -data Taxonomy = Taxonomy - { _accession :: !(Accession Species) - , _name :: !(Identification Species) - , _classification :: [Classification] - } deriving (Show) - -makeLenses ''Taxonomy - +{- -- | Given a name such as "Drosophila Melanogaster", returns "d.melanogaster". shortenName :: Identification Species -> Identification Species @@ -34,3 +26,5 @@ | [w] <- ws = IDD w | otherwise = IDD . BS.map toLower $ BS.take 1 (ws!!0) `BS.append` (BS.cons '.' $ ws!!1) where ws = BS.words xs +-} + diff --git a/Biobase/SElab/Taxonomy/Import.hs b/Biobase/SElab/Taxonomy/Import.hs --- a/Biobase/SElab/Taxonomy/Import.hs +++ b/Biobase/SElab/Taxonomy/Import.hs @@ -1,63 +1,70 @@ -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE NoMonomorphismRestriction #-} -{-# LANGUAGE BangPatterns #-} --- | Iteratee-based importer. Provides a simple "fromFile" function that --- produces both maps in one pass. +-- | Import an Rfam @taxonomy.txt@ file. Provides a simple "fromFile" +-- function that produces both maps in one pass. @fromFile@ will check the +-- file suffix and on @.gz@ suffixes decompress the input on-the-fly. module Biobase.SElab.Taxonomy.Import where -import Control.Applicative -import Control.Lens -import Data.Attoparsec as A hiding (parse) -import Data.Attoparsec.Char8 (char,decimal) -import Data.ByteString.Char8 as BS -import Data.Conduit as C -import Data.Conduit.Attoparsec -import Data.Conduit.Binary as CB -import Data.Conduit.List as CL -import Data.Conduit.Util as C -import Data.Either.Unwrap as E -import Data.List as L -import Data.Map as M -import qualified Data.Attoparsec.ByteString as AB hiding (parse) -import qualified Data.Attoparsec.Char8 as A8 +import Codec.Compression.GZip (decompress) +import Control.Applicative ((<|>)) +import Control.Arrow (second) +import Control.Monad +import Data.Attoparsec.Text.Lazy as AT +import Data.Char (isDigit) +import Data.HashMap.Strict (HashMap) +import Data.List (foldl') +import Data.Text.Lazy.Encoding (decodeUtf8) +import Data.Text.Lazy.IO as TL +import Data.Text (Text) +import Data.Vector (fromList) +import qualified Data.ByteString.Lazy as BL +import qualified Data.HashMap.Strict as HM +import qualified Data.Text as T +import qualified Data.Text.Lazy as TL +import System.FilePath (takeExtension) +import Data.String.Conversions.Monomorphic -import Biobase.SElab.Taxonomy -import Biobase.SElab.Types +import Biobase.Types.Accession (Accession(..),Species) +import Biobase.Types.Names +import Biobase.Types.Taxonomy -parse = CB.lines - =$ CL.map (parseOnly mkTaxonomy) - =$ CL.filter isRight - =$ CL.map fromRight - =$ C.zipSinks mapIdTaxonomy mapAcTaxonomy -{-# INLINE parse #-} +-- | Parse a single Taxon line. +-- +-- TODO there are unknown words at the end of each line. make those known -mkTaxonomy :: Parser Taxonomy -mkTaxonomy = f <$> ptaxid <* tab <*> pname <* tab <*> takeByteString where - f k n xs = let - cs = L.map (Classification . copy . BS.dropWhile (==' ')) . BS.split ';' . BS.init $ xs - in Taxonomy (ACC k) (IDD $ copy n) cs - ptaxid = decimal - pname = A8.takeWhile (/='\t') - tab = char '\t' -{-# INLINE mkTaxonomy #-} +parseTaxon :: Parser Taxon +parseTaxon = do + accession <- Accession <$> takeWhile1 isDigit "accession" + skipSpace "1st space" + species <- speciesName <$> takeWhile1 (/='\t') "species" + skipSpace "2nd space" + classification <- (fromList . map (,Unknown) . map fromST) <$> takeWhile1 (\z -> z/=';' && z/='.') `sepBy` "; " "classification" + unknowns <- manyTill anyChar (endOfInput <|> endOfLine) + return $ Taxon {..} -mapIdTaxonomy :: Monad m => GSink Taxonomy m (M.Map (Identification Species) Taxonomy) -mapIdTaxonomy = CL.fold f M.empty where - f !mp x = M.insert (x ^. name) x mp -{-# INLINE mapIdTaxonomy #-} +-- | Taxonomy according to @Infernal@ stored in two hashmaps. The first +-- from @Accession@ to @Taxon@, the second from species name to @Taxon@. -mapAcTaxonomy :: Monad m => GSink Taxonomy m (M.Map (Accession Species) Taxonomy) -mapAcTaxonomy = CL.fold f M.empty where - f !mp x = M.insert (x ^. accession) x mp -{-# INLINE mapAcTaxonomy #-} +type Taxonomy = ( HashMap (Accession Species) Taxon -- ^ find @Taxon@ via accession number + , HashMap SpeciesName Taxon -- ^ find @Taxon@ via species name + ) -fromFile :: String -> IO ( Map (Identification Species) Taxonomy - , Map (Accession Species) Taxonomy - ) -fromFile fname = do - runResourceT $ CB.sourceFile fname $$ parse -{-# NOINLINE fromFile #-} +-- | Parses the taxonomy.txt file. + +parseTaxonomy :: Parser Taxonomy +parseTaxonomy = foldl' go (HM.empty , HM.empty) <$> manyTill parseTaxon endOfInput + where go (!a, !s) x = (HM.insert (accession x) x a, HM.insert (species x) x s) + +-- | Read @taxonomy.txt.gz@ / @taxonomy.txt@ file into structure. + +fromFile :: FilePath -> IO Taxonomy +fromFile f = case takeExtension f of + ".gz" -> (go . decodeUtf8 . decompress) <$> BL.readFile f + _ -> go <$> TL.readFile f + where go txt = case AT.parse parseTaxonomy txt of + Done "" r -> r + Done ncon r -> error $ "unconsumed input " ++ f ++ ": " ++ (TL.unpack $ TL.take 1000 ncon) + Fail ncon ctx err -> error $ "error parsing " ++ f ++ ": " ++ show (ctx,err) + diff --git a/Biobase/SElab/Types.hs b/Biobase/SElab/Types.hs deleted file mode 100644 --- a/Biobase/SElab/Types.hs +++ /dev/null @@ -1,93 +0,0 @@ -{-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE GeneralizedNewtypeDeriving #-} -{-# LANGUAGE StandaloneDeriving #-} -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE EmptyDataDecls #-} - --- | Infernal Stockholm files and covariance models, and other related files --- use a bunch of different identifiers. We provide newtypes for more type --- safety. --- --- TODO Use (Bio.Core.Sequence.Offset) instead of Int for sequence info --- --- TODO move 'BitScore's, null models, probabilities into its own library. - -module Biobase.SElab.Types where - -import Control.Arrow -import Data.ByteString.Char8 as BS -import Data.Vector.Unboxed.Base -import Data.Vector.Generic as VG -import Data.Vector.Generic.Mutable as VGM -import Data.Vector.Unboxed as VU -import Data.Primitive.Types - - - --- * 'Accession' and string 'Identifier' with phantom types. - --- | Accession number, in the format of RFxxxxx, PFxxxxx, or CLxxxxx. We keep --- only the Int-part. A phantom type specifies which kind of accession number --- this is. For Species, we just have an index, it seems. - -newtype Accession t = ACC {unACC :: Int} - deriving (Eq,Ord,Read,Show) - --- | One word name for the family or clan. Phantom-typed with the correct type --- of model. Can be a longer name for species. - -newtype Identification t = IDD {unIDD :: ByteString} - deriving (Eq,Ord,Read,Show) - --- | Tag as being a clan. - -data Clan - --- | Tag as being a Pfam model. - -data Pfam - --- | Tag as being an Rfam model. Used for Stockholm and CM files. - -data Rfam - --- | Species have an accession number, too. - -data Species - - --- | Infernal bit score. Behaves like a double (deriving Num). --- --- Infernal users guide, p.42: log-odds score in log_2 (aka bits). --- --- S = log_2 (P(seq|CM) / P(seq|null)) - -newtype BitScore = BitScore {unBitScore :: Double} - deriving (Eq,Ord,Read,Show,Num,Prim) - -deriving instance Unbox BitScore -deriving instance VGM.MVector VU.MVector BitScore -deriving instance VG.Vector VU.Vector BitScore - --- | Given a null model and a probability, calculate the corresponding --- 'BitScore'. - -prob2Score :: Double -> Double -> BitScore -prob2Score null x - | x==0 = BitScore $ -10000 - | otherwise = BitScore $ log (x/null) / log 2 -{-# INLINE prob2Score #-} - --- | Given a null model and a 'BitScore' return the corresponding probability. - -score2Prob :: Double -> BitScore -> Double -score2Prob null (BitScore x) - | x<=(-9999) = 0 - | otherwise = null * exp (x * log 2) -{-# INLINE score2Prob #-} - --- | Classification names (taxonomic classification) - -newtype Classification = Classification {unClassification :: ByteString} - deriving (Eq,Ord,Read,Show) - diff --git a/BiobaseInfernal.cabal b/BiobaseInfernal.cabal --- a/BiobaseInfernal.cabal +++ b/BiobaseInfernal.cabal @@ -1,16 +1,18 @@ name: BiobaseInfernal -version: 0.7.1.0 +version: 0.8.1.0 author: Christian Hoener zu Siederdissen -maintainer: choener@tbi.univie.ac.at -homepage: http://www.tbi.univie.ac.at/~choener/ -copyright: Christian Hoener zu Siederdissen, 2011 +maintainer: choener@bioinf.uni-leipzig.de +homepage: https://github.com/choener/BiobaseInfernal +bug-reports: https://github.com/choener/BiobaseInfernal/issues +copyright: Christian Hoener zu Siederdissen, 2011 - 2017 category: Bioinformatics -synopsis: Infernal data structures and tools license: GPL-3 license-file: LICENSE build-type: Simple stability: experimental -cabal-version: >= 1.6.0 +cabal-version: >= 1.10.0 +tested-with: GHC == 7.10.3, GHC == 8.0.1 +synopsis: Infernal data structures and tools description: Provides import and export facilities for Infernal/Rfam data formats. We include Stockholm, CM, verbose Infernal results, @@ -20,78 +22,167 @@ annotations. This extension should be backward-compatible with standard-compliant parsers. . - This package uses Int's to store sequence position information. - Don't compile for 32bit. (And yes, this is a TODO, to change to - Int64). - . - . - . - Changes in 0.7.1.0 - . - * bumped dependencies, compatible to newest BiobaseXNA, PrimitiveArray - . - Changes in 0.7.0.0 - . - * work-in-progress release (some features missing) - . - * working CM parsing - . - * type defns have changed. using phantom types to specify what kind of model we are working with - . - * using conduit instead of iteratee + The @cmsearchFilter@ program provides filtering and coloring + options. extra-source-files: + changelog.md + README.md + + library - build-depends: - base >3 && <5, - attoparsec == 0.10.* , - attoparsec-conduit == 0.5.* , - biocore == 0.2 , - bytestring , - bytestring-lexing == 0.4.* , - conduit == 0.5.* , - containers, - either-unwrap == 1.1 , - lens == 3.* , - primitive >= 0.5 , - repa == 3.2.* , - transformers == 0.3.* , - tuple == 0.2.* , - vector >= 0.10 , - BiobaseXNA == 0.7.* , - PrimitiveArray == 0.5.* + build-depends: base >= 4.7 && < 5.0 + , aeson >= 0.9 + , attoparsec >= 0.12 + , binary >= 0.7 + , bytestring + , cereal >= 0.4 + , cereal-text >= 0.1 + , cereal-vector >= 0.2 + , containers + , data-default >= 0.5 + , deepseq >= 1.4 + , filepath >= 1.3 + , hashable >= 1.2 + , lens >= 4.0 + , parallel >= 3.0 + , pipes >= 4.0 + , pipes-attoparsec >= 0.5 + , pipes-bytestring >= 2.0 + , pipes-parse >= 3.0 + , pipes-safe >= 2.2 + , pipes-zlib >= 0.4 + , primitive >= 0.5 + , strict >= 0.3 + , string-conversions >= 0.4 + , text >= 1.0 + , text-binary >= 0.1 + , transformers + , tuple >= 0.3 + , unordered-containers >= 0.2 + , utf8-string >= 1.0 + , vector >= 0.10 + , vector-th-unbox >= 0.2 + , zlib >= 0.6 + -- own libraries (keep tight control over versioning?) + , BiobaseTypes == 0.1.2.* + , BiobaseXNA == 0.9.3.* + , DPutils == 0.0.1.* + , PrimitiveArray == 0.8.0.* exposed-modules: --- Biobase.Infernal --- Biobase.Infernal.Align --- Biobase.Infernal.Align.Import --- Biobase.Infernal.Clan --- Biobase.Infernal.Clan.Import + Biobase.SElab.Model + Biobase.SElab.Model.Import + Biobase.SElab.Model.Types Biobase.SElab.CM Biobase.SElab.CM.Import + Biobase.SElab.CM.ModelStructure + Biobase.SElab.CM.Types + Biobase.SElab.Common.Parser Biobase.SElab.HMM Biobase.SElab.HMM.Import --- Biobase.Infernal.Hit --- Biobase.Infernal.RfamFasta --- Biobase.Infernal.RfamFasta.Import - Biobase.SElab.RfamNames - Biobase.SElab.RfamNames.Import --- Biobase.Infernal.TabularHit --- Biobase.Infernal.TabularHit.Import + Biobase.SElab.HMM.Types Biobase.SElab.Taxonomy Biobase.SElab.Taxonomy.Import - Biobase.SElab.Types --- Biobase.Infernal.VerboseHit --- Biobase.Infernal.VerboseHit.Export --- Biobase.Infernal.VerboseHit.Import --- Biobase.Infernal.VerboseHit.Internal - + default-language: + Haskell2010 + default-extensions: BangPatterns + , DataKinds + , DeriveGeneric + , FlexibleContexts + , FlexibleInstances + , GADTs + , GeneralizedNewtypeDeriving + , LambdaCase + , MultiParamTypeClasses + , MultiWayIf + , OverloadedStrings + , ParallelListComp + , PatternSynonyms + , RankNTypes + , RecordWildCards + , ScopedTypeVariables + , TemplateHaskell + , TupleSections + , TypeFamilies + , TypeOperators + , TypeSynonymInstances + , UndecidableInstances ghc-options: -O2 -funbox-strict-fields + + +-- provides advanced cmsearch hit filtering functionality. + +executable cmsearchFilter + build-depends: base + , cmdargs >= 0.10 + -- + , BiobaseInfernal + hs-source-dirs: + src + main-is: + cmsearchFilter.hs + default-language: + Haskell2010 + ghc-options: + -O2 + + + +benchmark parsing + type: + exitcode-stdio-1.0 + build-depends: base + , criterion >= 1.1.0.0 + , lens + , text + , transformers + -- + , BiobaseInfernal + hs-source-dirs: + tests + main-is: + parsing.hs + default-language: + Haskell2010 + ghc-options: + -O2 -threaded -rtsopts "-with-rtsopts=-N -I0" + + + +test-suite properties + type: + exitcode-stdio-1.0 + main-is: + properties.hs + ghc-options: + -threaded -rtsopts -with-rtsopts=-N + hs-source-dirs: + tests + default-language: + Haskell2010 + default-extensions: OverloadedStrings + , ScopedTypeVariables + , TemplateHaskell + build-depends: base + , HUnit >= 1.2 + , lens + , QuickCheck + , tasty >= 0.11 + , tasty-hunit >= 0.9 + , tasty-quickcheck >= 0.8 + , tasty-th >= 0.1 + -- + , BiobaseInfernal + + + source-repository head type: git location: git://github.com/choener/BiobaseInfernal + diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,19 @@ +[![Build Status](https://travis-ci.org/choener/BiobaseInfernal.svg?branch=master)](https://travis-ci.org/choener/BiobaseInfernal) + +# BiobaseInfernal + +Low-level handling of Infernal file formats. Appropriate efficient Haskell data +structures for high-performance computations with Infernal covariance models. + +In particular, the *_states* structure provides a numerics-friendly interface +to the actual model. + + + +#### Contact + +Christian Hoener zu Siederdissen +Leipzig University, Leipzig, Germany +choener@bioinf.uni-leipzig.de +http://www.bioinf.uni-leipzig.de/~choener/ + diff --git a/changelog.md b/changelog.md new file mode 100644 --- /dev/null +++ b/changelog.md @@ -0,0 +1,27 @@ +0.8.1.0 +------- + +- parse CM files with version 0.7 as well +- ADPfusion-0.5.2 + +0.8.0.0 +------- + +- updated travis file to use stack + +0.7.1.0 +------- + +- bumped dependencies, compatible to newest BiobaseXNA, PrimitiveArray + +0.7.0.0 +------- + +- work-in-progress release (some features missing) + +- working CM parsing + +- type defns have changed. using phantom types to specify what kind of model +we are working with + +- using conduit instead of iteratee diff --git a/src/cmsearchFilter.hs b/src/cmsearchFilter.hs new file mode 100644 --- /dev/null +++ b/src/cmsearchFilter.hs @@ -0,0 +1,9 @@ +{-# LANGUAGE DeriveDataTypeable #-} + +module Main where + +import System.Console.CmdArgs + + + +main = return () diff --git a/tests/parsing.hs b/tests/parsing.hs new file mode 100644 --- /dev/null +++ b/tests/parsing.hs @@ -0,0 +1,90 @@ + +-- | We benchmark different stages of parsing. @small@, @medium@, and +-- @large@ are three sets of CMs with different file sizes. +-- +-- TODO +-- +-- - only extract header information +-- - single-threaded full parsing +-- - multi-threaded full parsing + +module Main where + +import Control.Lens +import Control.Monad.Trans.Writer.Strict +import Control.Monad (when) +import Criterion.Main +import Data.List as L +import qualified Data.Text as T +import qualified Data.Text.IO as T +import Data.Maybe (maybeToList) + +import Biobase.SElab.CM.Types as CM +import Biobase.SElab.HMM.Types as HMM +import Biobase.SElab.Model.Import +import Biobase.SElab.Model.Types + + + +{- + +-- | Parse only the header + +parseOnlyHeader :: Int -> String -> IO [PreModel] +parseOnlyHeader n file = do + (xs,log) <- runResourceT $ runWriterT $ sourceFile file $= ungzip $= preModels $= ccm n $$ consume + when (not $ T.null log) $ T.putStr log +-- print $ length xs +-- print $ xs ^.. ix 0 . _Right . _1 . CM.name +-- print $ xs ^.. ix 1 . _Left . _1 . HMM.name + return xs + +parseFull :: Int -> String -> IO [CM] +parseFull n file = do + (xs,log) <- runResourceT $ runWriterT $ sourceFile file $= ungzip $= preModels $= ccm n $= finalizeModels 1 $= attachHMMs $$ consume + when (not $ T.null log) $ T.putStr log + return xs + +parseFullPar :: Int -> String -> IO [CM] +parseFullPar n file = do + (xs,log) <- runResourceT $ runWriterT $ sourceFile file $= ungzip $= preModels $= ccm n $= finalizeModels 64 $= attachHMMs $$ consume + when (not $ T.null log) $ T.putStr log + return xs + +main :: IO () +main = defaultMain + -- only parse the header + [ bgroup "header" + [ bench "x 1" $ nfIO $ parseOnlyHeader 1 "./rfam-models/CL00001.cm.gz" + ] + -- full parsing, single-threaded + , bgroup "singlethreaded" + [ bench "small" $ nfIO $ parseFull 1 "./rfam-models/CL00001.cm.gz" + , bench "medium" $ nfIO $ parseFull 10 "./rfam-models/CL00001.cm.gz" + , bench "large" $ nfIO $ parseFull 100 "./rfam-models/CL00001.cm.gz" + ] + -- full parsing, multi-threaded + , bgroup "multithreaded" + [ bench "small" $ nfIO $ parseFullPar 1 "./rfam-models/CL00001.cm.gz" + , bench "medium" $ nfIO $ parseFullPar 10 "./rfam-models/CL00001.cm.gz" + , bench "large" $ nfIO $ parseFullPar 100 "./rfam-models/CL00001.cm.gz" + ] + ] + + + +ccm n = go where + go = do + xx <- await + yy <- await + case xx of + Nothing -> return () + Just x -> do + Prelude.mapM_ yield $ L.concat $ L.replicate n $ x : maybeToList yy + go + +-} + +main :: IO () +main = return () + diff --git a/tests/properties.hs b/tests/properties.hs new file mode 100644 --- /dev/null +++ b/tests/properties.hs @@ -0,0 +1,48 @@ + +module Main where + +import Control.Lens +import Data.Either (isLeft) +import Data.Monoid (mempty) +import Test.HUnit +import Test.QuickCheck.Modifiers +import Test.QuickCheck.Property +import Test.Tasty +import Test.Tasty.HUnit (testCase) +import Test.Tasty.QuickCheck (testProperty) +import Test.Tasty.TH + +import Biobase.SElab.CM.ModelStructure (flexibleToStatic, staticToFlexible) +import qualified Biobase.SElab.CM as CM +import qualified Biobase.SElab.HMM as HMM +import qualified Biobase.SElab.Model as Mdl + + + +-- | Import a HMMER 3 HMM and do some basic consistency checks + +case_HMM_import = do + hmms <- HMM.hmmFromFile "tests/test.hmm" + let h = head hmms + assertEqual "tests/test.hmm has a single HMMER 3 HMM" 1 $ length hmms + assertEqual "unknown lines:" mempty $ h ^. HMM.unknownLines + assertEqual "amino alphabet" "amino" $ h ^. HMM.alphabet + +case_CM__import = do + -- (cms, log) <- Mdl.fromFile "tests/test11.cm" 1 (const True) + cms <- CM.cmFromFile "tests/test11.cm" + let c = head cms + --assertEqual "Mdl.fromFile has empty log" "" log + assertEqual "tests/test.CM has a single Infernal 1.1 CM" 1 $ length cms + assertEqual "unknown lines in CM:" mempty $ c ^. CM.unknownLines + assertEqual "unknown lines in sub-HMM:" mempty $ c ^. CM.hmm . HMM.unknownLines + assertEqual "RNA alphabet" "RNA" $ c ^. CM.alph + assertBool "is flexible model" $ isLeft $ c ^. CM.cm + let Left flx = c ^. CM.cm -- just asserted + assertEqual "flexible . static" flx (staticToFlexible . flexibleToStatic $ flx) + + + +main :: IO () +main = $(defaultMainGenerator) +