BiobaseInfernal 0.7.1.0 → 0.8.1.0
raw patch · 22 files changed
+2138/−969 lines, 22 filesdep +BiobaseInfernaldep +BiobaseTypesdep +DPutilsdep −attoparsec-conduitdep −biocoredep −bytestring-lexingdep ~BiobaseXNAdep ~PrimitiveArraydep ~attoparsecnew-component:exe:cmsearchFilter
Dependencies added: BiobaseInfernal, BiobaseTypes, DPutils, HUnit, QuickCheck, aeson, binary, cereal, cereal-text, cereal-vector, cmdargs, criterion, data-default, deepseq, filepath, hashable, parallel, pipes, pipes-attoparsec, pipes-bytestring, pipes-parse, pipes-safe, pipes-zlib, strict, string-conversions, tasty, tasty-hunit, tasty-quickcheck, tasty-th, text, text-binary, unordered-containers, utf8-string, vector-th-unbox, zlib
Dependencies removed: attoparsec-conduit, biocore, bytestring-lexing, conduit, either-unwrap, repa
Dependency ranges changed: BiobaseXNA, PrimitiveArray, attoparsec, base, lens, transformers, tuple
Files
- Biobase/SElab/CM.hs +6/−241
- Biobase/SElab/CM/Import.hs +211/−207
- Biobase/SElab/CM/ModelStructure.hs +668/−0
- Biobase/SElab/CM/Types.hs +296/−0
- Biobase/SElab/Common/Parser.hs +38/−0
- Biobase/SElab/HMM.hs +7/−70
- Biobase/SElab/HMM/Import.hs +127/−148
- Biobase/SElab/HMM/Types.hs +111/−0
- Biobase/SElab/Model.hs +8/−0
- Biobase/SElab/Model/Import.hs +216/−0
- Biobase/SElab/Model/Types.hs +36/−0
- Biobase/SElab/RfamNames.hs +0/−22
- Biobase/SElab/RfamNames/Import.hs +0/−59
- Biobase/SElab/Taxonomy.hs +11/−17
- Biobase/SElab/Taxonomy/Import.hs +58/−51
- Biobase/SElab/Types.hs +0/−93
- BiobaseInfernal.cabal +152/−61
- README.md +19/−0
- changelog.md +27/−0
- src/cmsearchFilter.hs +9/−0
- tests/parsing.hs +90/−0
- tests/properties.hs +48/−0
Biobase/SElab/CM.hs view
@@ -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
Biobase/SElab/CM/Import.hs view
@@ -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+--
+ Biobase/SElab/CM/ModelStructure.hs view
@@ -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 ((<to) . fst) ts++-- | Given a @CM@, add the necessary transitions to create local+-- beginnings.+--+-- Local beginnings are created by adding transitions from the @S 0@ state+-- to the main states of each node. Activating local ends does not modify+-- any existing transition or emission probabilities.+--+-- This will add @S 0@, @IL 1@ and @IR 2@ as parent state to all nodes.++addLocalBegins :: Bitscore -> 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+
+ Biobase/SElab/CM/Types.hs view
@@ -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++-}+
+ Biobase/SElab/Common/Parser.hs view
@@ -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"+
Biobase/SElab/HMM.hs view
@@ -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
Biobase/SElab/HMM/Import.hs view
@@ -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
+ Biobase/SElab/HMM/Types.hs view
@@ -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)+
+ Biobase/SElab/Model.hs view
@@ -0,0 +1,8 @@++module Biobase.SElab.Model+ ( fromFile+ , keepAllModels+ ) where++import Biobase.SElab.Model.Import+
+ Biobase/SElab/Model/Import.hs view
@@ -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+
+ Biobase/SElab/Model/Types.hs view
@@ -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+
− Biobase/SElab/RfamNames.hs
@@ -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
− Biobase/SElab/RfamNames/Import.hs
@@ -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 #-}
Biobase/SElab/Taxonomy.hs view
@@ -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+-}+
Biobase/SElab/Taxonomy/Import.hs view
@@ -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)+
− Biobase/SElab/Types.hs
@@ -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)-
BiobaseInfernal.cabal view
@@ -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+
+ README.md view
@@ -0,0 +1,19 @@+[](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/ +
+ changelog.md view
@@ -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
+ src/cmsearchFilter.hs view
@@ -0,0 +1,9 @@+{-# LANGUAGE DeriveDataTypeable #-}++module Main where++import System.Console.CmdArgs++++main = return ()
+ tests/parsing.hs view
@@ -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 ()+
+ tests/properties.hs view
@@ -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)+