BiobaseInfernal 0.6.2.0 → 0.7.0.0

raw patch · 30 files changed

+1067/−1299 lines, 30 filesdep +BiobaseXNAdep +attoparsec-conduitdep +bytestring-lexingdep −attoparsec-iterateedep −iterateedep −iteratee-compressdep ~PrimitiveArraydep ~attoparsecdep ~biocorePVP ok

version bump matches the API change (PVP)

Dependencies added: BiobaseXNA, attoparsec-conduit, bytestring-lexing, conduit, lens, primitive

Dependencies removed: attoparsec-iteratee, iteratee, iteratee-compress

Dependency ranges changed: PrimitiveArray, attoparsec, biocore, either-unwrap, transformers, tuple, vector

API changes (from Hackage documentation)

Files

• Biobase/Infernal.hs +0/−42
• Biobase/Infernal/Align.hs +0/−36
• Biobase/Infernal/Align/Import.hs +0/−69
• Biobase/Infernal/CM.hs +0/−67
• Biobase/Infernal/CM/Export.hs +0/−12
• Biobase/Infernal/CM/Import.hs +0/−131
• Biobase/Infernal/Clan.hs +0/−24
• Biobase/Infernal/Clan/Import.hs +0/−44
• Biobase/Infernal/Hit.hs +0/−36
• Biobase/Infernal/RfamFasta.hs +0/−68
• Biobase/Infernal/RfamFasta/Import.hs +0/−99
• Biobase/Infernal/TabularHit.hs +0/−39
• Biobase/Infernal/TabularHit/Import.hs +0/−57
• Biobase/Infernal/Taxonomy.hs +0/−32
• Biobase/Infernal/Taxonomy/Import.hs +0/−80
• Biobase/Infernal/Types.hs +0/−93
• Biobase/Infernal/VerboseHit.hs +0/−54
• Biobase/Infernal/VerboseHit/Export.hs +0/−106
• Biobase/Infernal/VerboseHit/Import.hs +0/−142
• Biobase/Infernal/VerboseHit/Internal.hs +0/−18
• Biobase/SElab/CM.hs +251/−0
• Biobase/SElab/CM/Import.hs +242/−0
• Biobase/SElab/HMM.hs +74/−0
• Biobase/SElab/HMM/Import.hs +178/−0
• Biobase/SElab/RfamNames.hs +22/−0
• Biobase/SElab/RfamNames/Import.hs +59/−0
• Biobase/SElab/Taxonomy.hs +36/−0
• Biobase/SElab/Taxonomy/Import.hs +63/−0
• Biobase/SElab/Types.hs +93/−0
• BiobaseInfernal.cabal +49/−50

@@ -1,42 +0,0 @@---- | Re-export the most import parts.--module Biobase.Infernal-  ( TabularHit(..)-  , thFromFile-  , eneeTabularHit-  , VerboseHit(..)-  , vhFromFile-  , eneeVerboseHit-  , vhEneeByteString-  , vhEneeByteStrings-  , SpeciesTaxonomy(..)-  , tFromFile-  , Clan(..)-  , cFromFile-  ) where--import Data.ByteString as BS-import Data.Iteratee as I--import Biobase.Infernal.TabularHit-import Biobase.Infernal.TabularHit.Import as TH-import Biobase.Infernal.VerboseHit-import Biobase.Infernal.VerboseHit.Import as VH-import Biobase.Infernal.VerboseHit.Export as VH-import Biobase.Infernal.Taxonomy-import Biobase.Infernal.Taxonomy.Import as T-import Biobase.Infernal.Clan-import Biobase.Infernal.Clan.Import as C--thFromFile = TH.fromFile--vhFromFile = VH.fromFile-vhEneeByteString :: Monad m => Enumeratee [VerboseHit] ByteString m a-vhEneeByteString = VH.eneeByteString-vhEneeByteStrings :: Monad m => Enumeratee [VerboseHit] [ByteString] m a-vhEneeByteStrings = VH.eneeByteStrings--tFromFile = T.fromFile--cFromFile = C.fromFile

@@ -1,36 +0,0 @@---- | "cmalign" provides two interesting results, bit scores of sequences--- aligned to the model and the alignments themselves.--module Biobase.Infernal.Align where--import Data.ByteString.Char8 (ByteString)--import Biobase.Infernal.Types------ | cmalign results, includes sequence scores if available.------ TODO stockholmAlignment, should be "biostockholm" (will be set after some--- fun iteratee tests). For now, the 'ByteString' holds everything needed to--- parse using biostockholm.--data Align = Align-  { modelIdentification :: ModelIdentification-  , sequenceScores      :: [SequenceScore]-  , stockholmAlignment  :: ByteString-  }---- | Individual sequence scores.------ TODO avgProbability should use Probability newtype--data SequenceScore = SequenceScore-  { sequenceName      :: !(ModelAccession,ModelIdentification,EmblAccession)  -- ^ sequence name, typically RFxxxxxx;RfamID;embl-accession-  , sLength           :: !Int       -- ^ aligned sequence length-  , totalBitScore     :: !BitScore  -- ^ total alignment bitscore-  , structureBitScore :: !BitScore  -- ^ structural score part-  , avgProbability    :: !Double    -- ^-  }-

@@ -1,69 +0,0 @@-{-# LANGUAGE OverloadedStrings #-}---- | Parses "cmalign" results.------ NOTE have not tested if this works with multiple results in a file, but--- could ;-)--module Biobase.Infernal.Align.Import where--import Data.Iteratee as I-import Data.Iteratee.Char as I-import Data.Iteratee.IO as I-import Data.Iteratee.ZLib as IZ-import Data.ByteString.Char8 as BS-import Prelude as P--import Biobase.Infernal.Align-import Biobase.Infernal.Types------ | Transforms bytestring to list of 'Align' data.--eneeAlign :: (Monad m) => Enumeratee ByteString [Align] m a-eneeAlign = enumLinesBS ><> convStream go where-  go = do-    -- lets start with some comment lines-    cs <- I.takeWhile (("#" ==) . BS.take 1)-    -- there should be score lines now-    ss <- I.takeWhile (\s -> "# STOCKHOLM 1.0" /= s && (not $ BS.null s))-    -- Stockholm lines-    xs <- I.takeWhile (/="//")-    x <- I.head-    return [Align-      { modelIdentification = ModelIdentification ""-      , sequenceScores = P.map mkScore ss-      , stockholmAlignment = BS.unlines $ xs++[x]-      }]---- | Creates the required sequence score.--mkScore s = SequenceScore-  { sequenceName = undefined $ ws!!0-  , sLength = read . BS.unpack $ ws!!1-  , totalBitScore = BitScore . read . BS.unpack $ ws!!2-  , structureBitScore = BitScore . read . BS.unpack $ ws!!3-  , avgProbability = read . BS.unpack $ ws!!4-  } where ws = BS.words s---- | Convenience function creating all maps.--fromFileZip :: FilePath -> IO [Align]-fromFileZip fp = run =<< ( enumFile 8192 fp-                         . joinI-                         . enumInflate GZipOrZlib defaultDecompressParams-                         . joinI-                         . eneeAlign-                         $ stream2stream-                         )---- | Convenience function creating all maps.--fromFile :: FilePath -> IO [Align]-fromFile fp = run =<< ( enumFile 8192 fp-                      . joinI-                      . eneeAlign-                      $ stream2stream-                      )-

@@ -1,67 +0,0 @@---- | Infernal CMs.--module Biobase.Infernal.CM where--import Data.ByteString as BS-import Data.Map as M-import Data.Vector as V-import Data.Vector.Unboxed as VU--import Data.PrimitiveArray-import Data.PrimitiveArray.Ix--import Biobase.Infernal.Types------ | A datatype representing Infernal covariance models. This is a new--- representation that is incompatible with the one once found in "Biobase".--- The most important difference is that lookups are mapped onto efficient data--- structures, currently "PrimitiveArray".------ [1] Each "State" of a covariance model has up to 6 transition scores, hence--- we need s*6 cells for transitions.------ [2] Each "State" of a covariance has up to 16 emission scores, so we have--- s*16 cells for emissions, with unused cells set to a really high score.------ On top of these basic structures, we then place additional high-level--- constructs.------ [3] 'paths' are allowed transitions. This can safe a check, if the--- transition is encoded with a forbidden score.------ [4] 'localBegin' and 'localEnd' are local entry and exit strategies. A--- 'localBegin' is a transition score to certain states, all such transitions--- are in 'begins'. A 'localEnd' is a transition score to a local end state.------ NOTE that trustedCutoff > gathering > noiseCutoff------ TODO as with other projects, we should not use Double's but "Score" and--- "Probability" newtypes.--data CM = CM-  { name          :: ModelIdentification  -- ^ name of model as in "tRNA"-  , accession     :: ModelAccession       -- ^ RFxxxxx identification-  , trustedCutoff :: BitScore -- ^ lowest score of true member-  , gathering     :: BitScore -- ^ all scores at or above 'gathering' score are in the "full" alignment-  , noiseCutoff   :: Maybe BitScore -- ^ highest score NOT included as member-  , transition :: PrimArray (Int,Int) Double-  , emission :: PrimArray (Int,Int) Double-  , paths :: V.Vector (VU.Vector Double)-  , localBegin :: VU.Vector Double-  , begins :: VU.Vector Int-  , localEnd :: VU.Vector (Double)-  , nodes :: V.Vector (VU.Vector Int)-  }-  deriving (Show)---- | Map of model names to individual CMs.--type ID2CM = M.Map ModelIdentification CM---- | Map of model accession numbers to individual CMs.--type AC2CM = M.Map ModelAccession CM-

@@ -1,12 +0,0 @@---- | Transforms the internal representation of a CM back into a version that--- can be used by Infernal.------ Note that models are transformed into 'ByteString' as-is, the exporter does--- not make sure that probabilities add to one, that we write out probabilities--- instead of scores, and so on.------ TODO some of the notes above will become less problematic once we use--- newtypes, as a probability-CM will not be accepted by the exporter by then.--module Biobase.Infernal.CM.Export where

@@ -1,131 +0,0 @@-{-# LANGUAGE BangPatterns #-}-{-# LANGUAGE NoMonomorphismRestriction #-}-{-# LANGUAGE OverloadedStrings #-}---- | Iteratee-based parsing of Infernal covariance models.------ TODO does not create working CMs yet. Only partial key/value parsing is--- implemented.--module Biobase.Infernal.CM.Import where--import Control.Arrow-import Control.Monad (unless)-import Data.ByteString.Char8 as BS-import Data.Iteratee as I-import Data.Iteratee.Char as I-import Data.Iteratee.IO as I-import Data.Iteratee.Iteratee as I-import Data.Iteratee.ListLike as I-import Data.Iteratee.ZLib as IZ-import Data.Map as M-import Prelude as P-import Control.Monad.IO.Class (liftIO, MonadIO)--import Data.PrimitiveArray-import Data.PrimitiveArray.Ix--import Biobase.Infernal.CM-import Biobase.Infernal.Types------ * iteratee stuff---- | iteratee-based parsing of human-readable CMs.--eneeCM :: (Monad m) => Enumeratee ByteString [CM] m a-eneeCM = enumLinesBS ><> convStream f where-  f = do-    -- initial (mostly key/value) data-    hs' <- I.takeWhile (/="MODEL:")-    let hs = M.fromList . P.map (second (BS.dropWhile (==' ')) . BS.break (==' ')) $ hs'-    -- model begins-    mb <- I.tryHead-    unless (mb == Just "MODEL:") . error $ "model error: " ++ show (hs,mb,"head")-    -- nodes-    ns <- iterNodes-    -- model ends-    me <- I.tryHead-    unless (me == Just "//") . error $ "model error: " ++ show (hs,me,"tail")-    return . (:[]) $ CM-      { name = ModelIdentification $ hs M.! "NAME"-      , accession = ModelAccession . bsRead . BS.drop 2 $ hs M.! "ACCESSION"-      , gathering = BitScore . bsRead $ hs M.! "GA"-      , trustedCutoff = BitScore . bsRead $ hs M.! "TC"-      , noiseCutoff = let x = hs M.! "NC" in if x == "undefined" then Nothing else Just . BitScore . bsRead $ x-      , transition = error "not implemented yet"-      , emission = error "not implemented yet"-      , paths = error "not implemented yet"-      , localBegin = error "not implemented yet"-      , begins = error "not implemented yet"-      , localEnd = error "not implemented yet"-      , nodes = error "not implemented yet"-      } where bsRead = read . BS.unpack--iterNodes :: (Monad m) => Iteratee [ByteString] m [Node]-iterNodes = do-  hdr' <- I.head-  let (ishdr,(hdr,nidx)) = isNodeHeader hdr'-  unless ishdr $ error $ show hdr'-  xs <- I.takeWhile (fst . isState)-  pk <- I.peek-  let n = Node-            { nodeHeader = hdr-            , nodeIndex = nidx-            }-  case pk of-    Just "//" -> return []-    Just x-      | (True,_) <- isNodeHeader x -> do-          ns <- iterNodes-          return $ n:ns-    e -> error $ show e--data Node = Node-  { nodeHeader :: ByteString-  , nodeIndex :: Int-  }--isNodeHeader :: ByteString -> (Bool,(ByteString,Int))-isNodeHeader xs = (isnh,(hdr,nidx)) where-  isnh = BS.elem '[' xs && BS.elem ']' xs-  [hdr,nidx'] = BS.words . BS.init . BS.takeWhile (/=']') . BS.drop 1 . BS.dropWhile (/='[') $ xs-  nidx = read . BS.unpack $ nidx'--isState :: ByteString -> (Bool,ByteString)-isState xs'-  | P.null xs = (False,"")-  | P.head xs `P.elem` [ "[", "//" ] = (False,"")-  | P.head xs `P.elem` [ "S", "IL", "IR", "MATR", "MR", "D", "MP", "ML", "B", "E" ] = (True,"")-  | otherwise = error $ show xs-  where-    xs = BS.words xs'---- * convenience functions---- | Read covariance models from file. This parser reads one or more CMs from--- file.--fromFile :: FilePath -> IO (ID2CM, AC2CM)-fromFile fp = run =<< ( enumFile 8192 fp-                      . joinI-                      . eneeCM-                      $ I.zip (mkMap name) (mkMap accession)-                      )---- | Read covariance models from a compressed file.--fromFileZip :: FilePath -> IO (ID2CM, AC2CM)-fromFileZip fp = run =<< ( enumFile 8192 fp-                         . joinI-                         . enumInflate GZipOrZlib defaultDecompressParams-                         . joinI-                         . eneeCM-                         $ I.zip (mkMap name) (mkMap accession)-                         )---- | map creation helper--mkMap f = I.foldl' (\ !m x -> M.insert (f x) x m) M.empty-

@@ -1,24 +0,0 @@-{-# LANGUAGE OverloadedStrings #-}---- | Rfam clans are a set of biologically related Rfam families. This module--- provides simple abstraction methods and loaders from file and ByteString.------ TODO load and parse with enumerator--module Biobase.Infernal.Clan where--import Data.ByteString.Char8 (ByteString)--import Biobase.Infernal.Types------ | Simple Rfam clan data.--data Clan = Clan-  { cAccession  :: !ClanAccession       -- ^ result of the "AC    CL00001" line, keeping "1" in this case.-  , cIdentifier :: !ClanIdentification  -- ^ the "ID    tRNA" line, keeping "tRNA".-  , cMembers    :: ![ModelAccession]    -- ^ all the "MB    RF00005;", "MB    RF00023;" lines, keeping "[5,23]".-  , cStrings    :: ![ByteString]        -- ^ all lines of each clan, without any processing (except being in lines).-  } deriving (Read,Show,Eq)-

@@ -1,44 +0,0 @@-{-# LANGUAGE PatternGuards #-}-{-# LANGUAGE OverloadedStrings #-}---- | Importing clan data is probably never time-critical as the total file size--- is extremely small. Should this ever change, swap to iteratee. The file is--- read /strictly/.--module Biobase.Infernal.Clan.Import where--import qualified Data.ByteString.Char8 as BS-import Data.List--import Biobase.Infernal.Clan-import Biobase.Infernal.Types------ | Import the complete data from an uncompressed source file.--fromFile :: FilePath -> IO [Clan]-fromFile fp = fromByteString `fmap` BS.readFile fp---- | Transform a bytestring into a list of 'Clan's.--fromByteString :: BS.ByteString -> [Clan]-fromByteString s = map mkClan-                 . groupBy (\x y -> "AC"/=(head . BS.words $y))-                 . BS.lines-                 $ s---- | Given a list of bytestrings, create one Clan.------ TODO return Maybe, make crash-safe (not really high on the list...)--mkClan :: [BS.ByteString] -> Clan-mkClan xs = Clan-  { cAccession  = ClanAccession . f . BS.drop 2 . (!!1) . BS.words . head . filter ((=="AC") . BS.take 2) $ xs-  , cIdentifier = ClanIdentification . (!!1) . BS.words . head . filter ((=="ID") . BS.take 2) $ xs-  , cMembers    = map (ModelAccession . f . BS.drop 2 . BS.init . (!!1)) . filter ((=="MB") . (!!0)) . map BS.words $ xs-  , cStrings    = xs-  } where-      f s-        | Just (k, _) <- BS.readInt s = k-        | otherwise = error $ "mkClan: " ++ BS.unpack s

@@ -1,36 +0,0 @@---- | Accessors for Infernal hits.------ TODO modelStartStop pair? same for target?------ TODO newtypes for these returns?--module Biobase.Infernal.Hit where--import Data.ByteString.Char8 (ByteString)--import Biobase.Infernal.Types------ | Generalized accessors for VerboseHit's and TabularHit's.--class Hit a where-  -- | Model name (like 5S_rRNA).-  model       :: a -> ModelIdentification-  -- | Target name, typically the scaffold or chromosome where the hit occurs.-  target      :: a -> Scaffold-  -- | Start of submodel.-  modelStart  :: a -> Int-  -- | Stop of submodel.-  modelStop   :: a -> Int-  -- | Start of substring in target.-  targetStart :: a -> Int-  -- | Stop of substring in target.-  targetStop  :: a -> Int-  -- | Bit score of the hit of model in target.-  bitScore    :: a -> BitScore-  -- | Evalue, expectation of bit score of higher in target sequence of length.-  evalue      :: a -> Double-  -- | G/C content in target.-  gcPercent   :: a -> Int

@@ -1,68 +0,0 @@-{-# LANGUAGE OverloadedStrings #-}-{-# LANGUAGE RecordWildCards #-}---- | The Rfam.fasta.gz file provides useful information: (1) conversion between--- Rfam accession and Rfam identifier, (2) species accession, (3) name of said--- species, and (4) the sequence fasta file.--module Biobase.Infernal.RfamFasta where--import Bio.Core.Sequence-import Data.ByteString.Char8 as BS-import Data.Map as M-import qualified Data.ByteString.Lazy.Char8 as BSL-import Text.Printf--import Biobase.Infernal.Types------ | Rfam FASTA entry.--data RfamFasta = RfamFasta-  { modelAccession    :: !ModelAccession      -- ^ Rfam accession number RFxxxxx (the xxxxx part).-  , modelIdentifier   :: !ModelIdentification -- ^ Rfam identifier (like 5S_rRNA).-  , sequenceAccession :: !EmblAccession       -- ^ EMBL sequence accession identifier and position.-  , speciesAccession  :: !SpeciesAccession    -- ^ Rfam species accession.-  , speciesName       :: !SpeciesName         -- ^ Species name.-  , fastaData         :: !StrictSeqData       -- ^ FASTA data-  } deriving (Show)---- | Since RfamFasta entries are just fasta entries...--instance BioSeq RfamFasta where-  seqlabel RfamFasta{..}  = SeqLabel . BSL.fromChunks $ [BS.concat-    [ BS.pack . printf "RF%05d" . unModelAccession $ modelAccession-    , ";"-    , unModelIdentification modelIdentifier-    , ";"-    , let (a,b,c) = unEmblAccession sequenceAccession in BS.concat [a, "/", BS.pack $ show b, "-", BS.pack $ show c]-    , "   "-    , BS.pack . show . unSpeciesAccession $ speciesAccession-    , ":"-    , unSpeciesName speciesName-    ] ]-  seqdata RfamFasta{..}   = SeqData . BSL.fromChunks $ [unStrictSeqData fastaData]-  seqlength RfamFasta{..} = Offset . fromInteger . toInteger . BS.length . unStrictSeqData $ fastaData------ * Some in-memory lookup systems.---- | Model accession to model identifier--type ModelAC2ID = Map ModelAccession ModelIdentification---- | Model identifier to model accession--type ModelID2AC = Map ModelIdentification ModelAccession---- | Model accession and sequence accession to 'RfamFasta' entry (and model--- accession to all entries for this accession).--type ACAC2RfamFasta = Map ModelAccession (Map EmblAccession RfamFasta)---- | Model identifier and sequence accession to 'RfamFasta' entry.--type IDAC2RfamFasta = Map ModelIdentification (Map EmblAccession RfamFasta)-

@@ -1,99 +0,0 @@-{-# LANGUAGE OverloadedStrings #-}-{-# LANGUAGE BangPatterns #-}---- | Imports an Rfam Fasta file and provides simultaneous export to four--- different data structures for lookups.--module Biobase.Infernal.RfamFasta.Import where--import Control.Arrow ((***))-import Data.ByteString.Char8 as BS-import Data.Iteratee as I-import Data.Iteratee.Char as I-import Data.Iteratee.IO as I-import Data.Iteratee.ZLib as IZ-import Data.Map as M-import Prelude as P--import Biobase.Infernal.RfamFasta-import Biobase.Infernal.Types------ | Enumeratee for RfamFasta entries from a ByteString.--eneeRfamFasta :: (Monad m) => Enumeratee ByteString [RfamFasta] m a-eneeRfamFasta = enumLinesBS ><> convStream f where-  f = do-    th <- I.tryHead-    case th of-      Nothing -> error "huh?"-      Just h  -> do-                   let (ana,sps) = (BS.split ';' *** BS.split ':' . BS.dropWhile (==' ')) . BS.break (==' ') $ h-                   fs <- I.takeWhile (\s -> ">" /= BS.take 1 s)-                   return . (:[]) $ RfamFasta-                     { modelAccession    = ModelAccession . read . P.drop 2 . unpack $ ana!!0-                     , modelIdentifier   = ModelIdentification $ ana!!1-                     , sequenceAccession = mkEmblAccession $ ana!!2-                     -- , speciesAC = maybe (error $ "ERROR: " ++ show (unpack $ sps!!0,unpack s)) fst . readInt $ sps!!0-                     , speciesAccession  = SpeciesAccession . maybe (-1) fst . readInt $ sps!!0-                     , speciesName = SpeciesName $ sps!!1-                     , fastaData = StrictSeqData . BS.copy . BS.concat $ fs-                     }------ * In-memory lookup---- | Create a mapping between rfam family accession numbers and rfam family--- names.--iModelAC2ID :: (Monad m) => Iteratee [RfamFasta] m ModelAC2ID-iModelAC2ID = I.foldl' f M.empty where-  f !m x = insertWith' const (modelAccession x) (modelIdentifier x) m---- | Create a mapping between rfam family names and rfam family accession--- numbers.--iModelID2AC :: (Monad m) => Iteratee [RfamFasta] m ModelID2AC-iModelID2AC = I.foldl' f M.empty where-  f !m x = insertWith' const (modelIdentifier x) (modelAccession x) m---- | Provides a mapping between (Rfam accession, sequence accession) and the--- complete 'RfamFasta'.--iACAC2RfamFasta :: (Monad m) => Iteratee [RfamFasta] m ACAC2RfamFasta-iACAC2RfamFasta = I.foldl' f M.empty where-  f !m x = insertWith' union (modelAccession x) (M.singleton (sequenceAccession x) x) m---- | Provides a mapping between (Rfam name, sequence accession) and the complete--- 'RfamFasta'.--iIDAC2RfamFasta :: (Monad m) => Iteratee [RfamFasta] m IDAC2RfamFasta-iIDAC2RfamFasta = I.foldl' f M.empty where-  f !m x = insertWith' union (modelIdentifier x) (M.singleton (sequenceAccession x) x) m------ * File reading.---- | Convenience function creating all maps.--fromFileZip :: FilePath -> IO (ModelAC2ID, ModelID2AC, ACAC2RfamFasta, IDAC2RfamFasta)-fromFileZip fp = run =<< ( enumFile 8192 fp-                         . joinI-                         . enumInflate GZipOrZlib defaultDecompressParams-                         . joinI-                         . eneeRfamFasta-                         $ I.zip4 iModelAC2ID iModelID2AC iACAC2RfamFasta iIDAC2RfamFasta-                         )---- | Convenience function creating all maps.--fromFile :: FilePath -> IO (ModelAC2ID, ModelID2AC, ACAC2RfamFasta, IDAC2RfamFasta)-fromFile fp = run =<< ( enumFile 8192 fp-                      . joinI-                      . eneeRfamFasta-                      $ I.zip4 iModelAC2ID iModelID2AC iACAC2RfamFasta iIDAC2RfamFasta-                      )-

@@ -1,39 +0,0 @@---- | Simple tabular hits as returned by Infernal.--module Biobase.Infernal.TabularHit where--import Data.ByteString.Char8 as BS--import Biobase.Infernal.Hit-import Biobase.Infernal.Types------ | Tabular Infernal hits. See Biobase.Infernal.Hit for description of the--- individual fields.--data TabularHit = TabularHit-  { thModel       :: !ModelIdentification-  , thTarget      :: !Scaffold-  , thTargetStart :: !Int-  , thTargetStop  :: !Int-  , thModelStart  :: !Int-  , thModelStop   :: !Int-  , thBitScore    :: !BitScore-  , thEvalue      :: !Double-  , thGCpercent   :: !Int-  } deriving (Read,Show)---- | Generalized accessors.--instance Hit TabularHit where-  model       = thModel-  target      = thTarget-  modelStart  = thModelStart-  modelStop   = thModelStop-  targetStart = thTargetStart-  targetStop  = thTargetStop-  bitScore    = thBitScore-  evalue      = thEvalue-  gcPercent   = thGCpercent

@@ -1,57 +0,0 @@-{-# LANGUAGE OverloadedStrings #-}---- Importing tabular hits is rather easy, as they are one entry per line.--module Biobase.Infernal.TabularHit.Import where--import Data.ByteString.Char8 as BS-import Data.Iteratee as I-import Data.Iteratee.Iteratee as I-import Data.Iteratee.ListLike as I-import Data.Iteratee.Char as I-import Data.Either.Unwrap-import Data.Attoparsec as A hiding (takeTill)-import Data.Attoparsec.Char8 as A-import Control.Applicative-import Data.Iteratee.IO as I--import Biobase.Infernal.TabularHit-import Biobase.Infernal.Types------ | Transform a stream into tabular hits.--eneeTabularHit :: (Functor m, Monad m) => Enumeratee ByteString [TabularHit] m a-eneeTabularHit = enumLinesBS ><> I.filter (\x -> not $ BS.null x || isPrefixOf "#" x) ><> mapStream f where-  f = fromRight . parseOnly p-  mkTH mName tName tStart tStop qStart qStop bScore eValue gc = TabularHit-    (ModelIdentification tName)-    (Scaffold tName)-    tStart-    tStop-    qStart-    qStop-    (BitScore bScore)-    eValue-    gc-  p = mkTH <$> pString  -- model name-           <*> pString  -- target name-           <*> pDecimal -- target start-           <*> pDecimal -- target stop-           <*> pDecimal -- query start-           <*> pDecimal -- query stop-           <*> pDouble  -- bit score-           <*> pDouble  -- evalue-           <*> pDecimal -- gc content-  pString  = A.skipSpace *> A.takeTill A.isSpace-  pDecimal = A.skipSpace *> A.decimal-  pDouble  = A.skipSpace *> A.double---- | Convenience function to load from file and return a big list of tabular--- hits.--fromFile :: FilePath -> IO [TabularHit]-fromFile fp = do-  i <- enumFile 8192 fp . joinI $ eneeTabularHit stream2stream-  run i

@@ -1,32 +0,0 @@-{-# LANGUAGE PatternGuards #-}---- | Infernal contains a taxonomy database. This is a simple module reflecting--- said database.--module Biobase.Infernal.Taxonomy where--import qualified Data.ByteString.Char8 as BS-import Data.Char (toLower)--import Biobase.Infernal.Types------ | 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 SpeciesTaxonomy = SpeciesTaxonomy-  { stAccession      :: !SpeciesAccession-  , stName           :: !SpeciesName-  , stClassification :: ![Classification]-  } deriving (Show)---- | Given a name such as "Drosophila Melanogaster", returns "d.melanogaster".--shortenName :: SpeciesName -> SpeciesName-shortenName (SpeciesName xs)-  | null ws   = SpeciesName xs-  | [w] <- ws = SpeciesName w-  | otherwise = SpeciesName . BS.map toLower $ BS.take 1 (ws!!0) `BS.append` (BS.cons '.' $ ws!!1)-  where ws = BS.words xs

@@ -1,80 +0,0 @@-{-# LANGUAGE BangPatterns #-}---- | Iteratee-based importer. Provides a simple "fromFile" function that--- produces both maps in one pass.--module Biobase.Infernal.Taxonomy.Import where--import Control.Applicative-import Data.Attoparsec as A-import Data.Attoparsec.Char8 as A8-import Data.Attoparsec.Iteratee-import Data.ByteString.Char8 as BS-import Data.Either.Unwrap as E-import Data.Iteratee as I-import Data.Iteratee.Char as I-import Data.Iteratee.IO as I-import Data.Iteratee.ListLike as I-import Data.List as L-import Data.Map as M--import Biobase.Infernal.Taxonomy-import Biobase.Infernal.Types------ | Provide name-based lookup as the most-common usage scenario.------ TODO there are 9 duplicates in the names, let's find them and see what is--- going on--iSpeciesMap :: Monad m => Iteratee [SpeciesTaxonomy] m (M.Map SpeciesName SpeciesTaxonomy)-iSpeciesMap = I.foldl' f M.empty where-  f !m x = M.insert (stName x) x m---- | And a map based on taxon id--iTaxIdMap :: Monad m => Iteratee [SpeciesTaxonomy] m (M.Map SpeciesAccession SpeciesTaxonomy)-iTaxIdMap = I.foldl' f M.empty where-  f !m x = M.insert (stAccession x) x m---- | Imports taxonomy data.--eneeSpecies :: Monad m => Enumeratee ByteString [Either String SpeciesTaxonomy] m a-eneeSpecies = enumLinesBS ><> mapStream (parseOnly mkSpecies)---- | Given a 'ByteString', create a species entry.------ NOTE The taxonomy format is, for each species, a line consisting of: taxid ---- tab - species name - tab - semicolon separated list of classification names--- - dot - end of line.--mkSpecies :: Parser SpeciesTaxonomy-mkSpecies = 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 SpeciesTaxonomy (SpeciesAccession k) (SpeciesName $ copy n) cs-  ptaxid   = decimal-  pname    = A8.takeWhile (/='\t')-  tab      = char '\t'---- | Convenience function: given a taxonomy file, produce both maps simultanously.--fromFile :: FilePath -> IO (M.Map SpeciesName SpeciesTaxonomy, M.Map SpeciesAccession SpeciesTaxonomy)-fromFile fp = do-  i <- enumFile 8192 fp-    . joinI-    . (eneeSpecies ><> I.filter isRight ><> mapStream fromRight)-    $ I.zip iSpeciesMap iTaxIdMap-  run i---- * Testing--{--test :: IO ()-test = do-  (s,t) <- fromFile "/home/choener/tmp/taxonomy"-  print $ M.size s-  print $ M.size t-  return ()--}

@@ -1,93 +0,0 @@-{-# LANGUAGE GeneralizedNewtypeDeriving #-}---- | All these different accession numbers and identifiers are confusing,--- newtype's to the rescue.------ TODO some of these names might have to change in the future...------ TODO Use INT64 instead of Int...--module Biobase.Infernal.Types where--import Control.Arrow-import Data.ByteString.Char8 as BS------ * Rfam Clans---- | Clan accession identifier--newtype ClanAccession = ClanAccession {unClanAccession :: Int}-  deriving (Eq,Ord,Read,Show)---- | Clan model name--newtype ClanIdentification = ClanIdentification {unClanIdentification :: ByteString}-  deriving (Eq,Ord,Read,Show)------ * Covariance models or Stockholm multiple alignments.---- | The numeric identifier of a covarience model or Stockholm multiple--- alignment as in RFxxxxx.--newtype ModelAccession = ModelAccession {unModelAccession :: Int}-  deriving (Eq,Ord,Read,Show)---- | String identifier of a covariance model or Stockholm multiple alignment as--- in "5S_rRNA".--newtype ModelIdentification = ModelIdentification {unModelIdentification :: ByteString}-  deriving (Eq,Ord,Read,Show)------ * Individual sequence information---- | EMBL sequence accession based on sequence accession and sequence start to--- stop. (Should this then be RfamSequenceAccession?)--newtype EmblAccession = EmblAccession {unEmblAccession :: (ByteString,Int,Int)}-  deriving (Eq,Ord,Read,Show)---- | Simple function to create 'EmblAccession' from a 'ByteString'.--mkEmblAccession :: ByteString -> EmblAccession-mkEmblAccession s = EmblAccession (sid,start,stop) where-  (sid,(Just (start,_),Just (stop,_))) = second ((BS.readInt *** (BS.readInt . BS.drop 1)) . BS.span (/='-') . BS.drop 1) . BS.span (/='/') $ s---- | Numeric species accession number.--newtype SpeciesAccession = SpeciesAccession {unSpeciesAccession :: Int}-  deriving (Eq,Ord,Read,Show)---- | String name for species.--newtype SpeciesName = SpeciesName {unSpeciesName :: ByteString}-  deriving (Eq,Ord,Read,Show)---- | Strict FASTA data.--newtype StrictSeqData = StrictSeqData {unStrictSeqData :: ByteString}-  deriving (Eq,Ord,Read,Show)---- | Classification names (taxonomic classification)--newtype Classification = Classification {unClassification :: ByteString}-  deriving (Eq,Ord,Read,Show)------ * More generic newtypes, sequence identification, etc---- | Identifies a certain scaffold or chromosome where a hit occurs--newtype Scaffold = Scaffold {unScaffold :: ByteString}-  deriving (Eq,Ord,Read,Show)---- | Infernal bit score. Behaves like a double (deriving Num).--newtype BitScore = BitScore {unBitScore :: Double}-  deriving (Eq,Ord,Read,Show,Num)

@@ -1,54 +0,0 @@-{-# LANGUAGE RecordWildCards #-}-{-# OPTIONS_GHC -funbox-strict-fields #-}---- | Provides a datatype for cmsearch verbose output. The Import/Export system--- now allows for primitive annotations using "##" as the first two characters.--- Annotations are only accepted for individual hits.------ TODE biocore / Strand for strand information?--module Biobase.Infernal.VerboseHit where--import Data.ByteString.Char8 as BS-import Text.Printf--import Biobase.Infernal.Hit-import Biobase.Infernal.Types------ | Captures a complete alignment--data VerboseHit = VerboseHit-  { vhTargetStart :: !Int          -- ^ part of target sequence (start counting at 1)-  , vhTargetStop  :: !Int-  , vhModelStart  :: !Int          -- ^ which part of the CM/stk do we align to-  , vhModelStop   :: !Int          -- ^ which part of the CM/stk do we align to-  , vhModel       :: !ModelIdentification   -- ^ the CM for this alignment-  , vhStrand      :: !Strand       -- ^ should be either '+' or '-'-  , vhBitScore    :: !BitScore     -- ^ bit score-  , vhEvalue      :: !Double       -- ^ number of hits we expect to find with 'score' or higher for 'targetSequence' length-  , vhPvalue      :: !Double       -- ^ ?-  , vhGCpercent   :: !Int          -- ^ ?-  , vhTarget      :: !Scaffold     -- ^ scaffold, chromosome, ... (the name of the sequence, not the sequence data!)-  , vhWuss        :: !ByteString   -- ^ fancy secondary structure annotation using wuss notation-  , vhConsensus   :: !ByteString   -- ^ query consensus (upper: highly, lower: weak/no)-  , vhScoring     :: !ByteString   -- ^ represents where positive and negative scores come from-  , vhSequence    :: !ByteString   -- ^ the target sequence which aligns to the model-  , vhAnnotation  :: ![ByteString] -- ^ any annotations that could be associated (# lines)-  } deriving (Show,Read)--type Strand = Char---- | Generalized accessors.--instance Hit VerboseHit where-  model       = vhModel-  target      = vhTarget-  modelStart  = vhModelStart-  modelStop   = vhModelStop-  targetStart = vhTargetStart-  targetStop  = vhTargetStop-  bitScore    = vhBitScore-  evalue      = vhEvalue-  gcPercent   = vhGCpercent

@@ -1,106 +0,0 @@-{-# LANGUAGE NoMonomorphismRestriction #-}-{-# LANGUAGE RankNTypes #-}-{-# LANGUAGE RecordWildCards #-}-{-# LANGUAGE OverloadedStrings #-}---- | Exports VerboseHit results back into text. As a likely scenario is a--- pipeline where hits are to be filtered out, this provides enumeratee's that--- handle additional annotations as required by the file format for CMs,--- scaffolds, and strand information. If you just need a way to show the data,--- use printVerboseHit.--module Biobase.Infernal.VerboseHit.Export where--import Control.Monad.Trans.Class (lift)-import Data.ByteString.Char8 as BS-import Data.Iteratee as I-import Data.Maybe-import Prelude as P-import Text.Printf--import Biobase.Infernal.Types-import Biobase.Infernal.VerboseHit-import Biobase.Infernal.VerboseHit.Internal------ | Transforms a list of verbose hits into a bytestring.------ TODO How to append the last line "//" to the finished stream, if at least--- one element was printed?--eneeByteString :: Monad m => Enumeratee [VerboseHit] ByteString m a-eneeByteString = eneeByteStrings ><> mapChunks BS.concat---- | This transformer keeps a 1-1 relationship between each 'VerboseHit' and--- bytestring representation. Useful for merging different streams, if--- individual 'VerboseHit's are to be annotated.--eneeByteStrings :: Monad m => Enumeratee [VerboseHit] [ByteString] m a-eneeByteStrings = unfoldConvStream f (AliGo BS.empty BS.empty '?' []) where-  f acc = do-    h <- I.head-    let na = newAcc acc h-    p <- I.peek-    return ( fst na-           , return . BS.unlines $ snd na ++ P.map (append "##") (vhAnnotation h)  ++ [showVerboseHit h] ++ maybe ["//"] (const []) p-           )---- | Given the current state "a" and verbose hit "h", determine if any state--- switches have to be emitted.--newAcc a@(AliGo{..}) h@VerboseHit{..}-  | otherwise = ( AliGo (unModelIdentification vhModel) (unScaffold vhTarget) vhStrand [], ls )-  where ls = [ "//" | aliCM /= BS.empty && bCM ] ++-             [ "CM: " `BS.append` unModelIdentification vhModel | bCM ] ++-             [ ">" `BS.append` unScaffold vhTarget `BS.append` "\n" | bCM || bSc] ++-             [ str `BS.append` " strand results:\n" | bCM || bSc || bSt ]-        bCM = aliCM /= unModelIdentification vhModel-        bSc = aliScaffold /= unScaffold vhTarget-        bSt = aliStrand /= vhStrand-        str-          | vhStrand == '+' = "Plus"-          | vhStrand == '-' = "Minus"-          | otherwise       = "Unknown"------ | Convert a 'VerboseHit' to a string, ready for printing as in the input--- file.--showVerboseHit :: VerboseHit -> BS.ByteString-showVerboseHit VerboseHit{..} = BS.unlines-  [ BS.pack $ printf " Query = %d - %d, Target = %d - %d"-                vhModelStart vhModelStop vhTargetStart vhTargetStop-  , BS.pack $ printf " Score = %.2f, E = %f, P = %.4e, GC = %d"-                (unBitScore vhBitScore) vhEvalue vhPvalue vhGCpercent-  , ""-  , ws11 `BS.append` vhWuss-  , (BS.pack $ printf "%10d " vhModelStart)-    `BS.append` vhConsensus-    `BS.append` (BS.pack $ printf " %d" vhModelStop)-  , ws11 `BS.append` vhScoring-  , (BS.pack $ printf "%10d " vhTargetStart)-    `BS.append` vhSequence-    `BS.append` (BS.pack $ printf " %d" vhTargetStop)-  ] where-    ws11 = BS.pack $ P.replicate 11 ' '----{----import Biobase.Infernal.VerboseHit.Import--test = do-  xs <- fromFile "/home/choener/tmp/infernal-1.0.2/tutorial/tmp.res"-  i <- enumList [xs] $ joinI $ eneeByteString stream2stream-  ys <- run i-  {--  BS.putStrLn ys-  print $ BS.length ys-  print $ P.length $ BS.lines ys-  -}-  BS.putStrLn $ BS.take 1000 ys-  return ()--}-

@@ -1,142 +0,0 @@-{-# LANGUAGE DoAndIfThenElse #-}-{-# LANGUAGE ScopedTypeVariables #-}-{-# LANGUAGE OverloadedStrings #-}-{-# LANGUAGE RecordWildCards #-}---- | Enumeratee that transforms a stream of 'ByteString's into a stream of--- 'VerboseHit's.--module Biobase.Infernal.VerboseHit.Import-  ( eneeVerboseHit-  , fromFile-  ) where--import Control.Applicative-import Control.Monad as M-import Data.Attoparsec as A-import Data.Attoparsec.Char8 as A8-import Data.Attoparsec.Iteratee as EAP-import Data.ByteString.Char8 as BS-import Data.Either.Unwrap-import Data.Iteratee as I-import Data.Iteratee.Char as I-import Data.Iteratee.IO as I-import Data.Iteratee.Iteratee as I-import Data.Iteratee.ListLike as I-import Data.Tuple.Select-import Prelude as P--import Biobase.Infernal.Types-import Biobase.Infernal.VerboseHit-import Biobase.Infernal.VerboseHit.Internal------ | Transforms a stream into verbose hits. We need to keep a state in the--- accumulator to keep track of the current CM, scaffold and strand.--eneeVerboseHit :: (Functor m, Monad m) => Enumeratee BS.ByteString [VerboseHit] m a-eneeVerboseHit = enumLinesBS ><> I.filter (not . BS.null) ><> unfoldConvStream f (AliGo BS.empty BS.empty '?' []) where-  f acc = do-    h' <- tryHead-    case h' of-      Nothing -> return (acc, [])-      (Just h)-        | "##"   `isPrefixOf` h -> return (acc{aliAnnotation = aliAnnotation acc ++ [BS.drop 2 h]},[])-        | "CM: " `isPrefixOf` h -> return (acc{aliCM = BS.copy $ BS.drop 4 h, aliAnnotation = []}, [])-        | ">"    `isPrefixOf` h -> return (acc{aliScaffold = BS.copy $ BS.drop 1 h, aliAnnotation = []}, [])-        | "Plus strand results"  `isInfixOf` h -> return (acc{aliStrand = '+', aliAnnotation = []}, [])-        | "Minus strand results" `isInfixOf` h -> return (acc{aliStrand = '-', aliAnnotation = []}, [])-        | " Query" `isInfixOf` h -> do-            x <- qs h (aliCM acc) (aliScaffold acc) (aliStrand acc) (aliAnnotation acc)-            return (acc{aliAnnotation = []},x)-        | otherwise -> return (acc,[])---- | Parses one CM query result.--qs :: Monad m => ByteString -> ByteString -> ByteString -> Char -> [ByteString] -> Iteratee [ByteString] m [VerboseHit]-qs query cm scaf pm anno = do-  let q = fromRight . parseOnly qt $ query-  s <- I.head >>= return . fromRight . parseOnly sepg-  l <- fourLines $ sel4 q-  return . pure $ VerboseHit-    { vhTarget = Scaffold scaf-    , vhModel = ModelIdentification cm-    , vhStrand = pm-    , vhModelStart = sel1 q-    , vhModelStop = sel2 q-    , vhTargetStart = sel3 q-    , vhTargetStop = sel4 q-    , vhBitScore = BitScore $ sel1 s-    , vhEvalue = sel2 s-    , vhPvalue = sel3 s-    , vhGCpercent = sel4 s-    , vhWuss = cpy $ l!!0-    , vhConsensus = cpy $ l!!1-    , vhScoring = cpy $ l!!2-    , vhSequence = cpy $ l!!3-    , vhAnnotation = anno-    }-  where-    cpy = BS.copy . BS.concat-    qt = (,,,) <$ A.string " Query = "   <*> A8.decimal <* A.string " - " <*> A8.decimal-               <* A.string ", Target = " <*> A8.decimal <* A.string " - " <*> A8.decimal-               <?> "qt"-    sepg = (,,,) <$ A.string " Score = " <*> A8.double-                 <* A.string ", E = "    <*> A8.double-                 <* A.string ", P = "    <*> A8.double-                 <* A.string ", GC = " <* A8.skipSpace <*> A8.decimal-                 <?> "sepg"---- | Parses multiple four-line elements.--fourLinesOld to = do-  I.dropWhile BS.null-  ls <- joinI $ I.take 4 stream2stream-  let ws = BS.length . BS.takeWhile isSpace . P.head $ ls-  let cs = BS.length . BS.dropWhile isSpace . P.head $ ls-  let xs = P.map (BS.take cs . BS.drop ws) ls-  if (P.length ls /= 4) ||-     ("-" == (P.last . BS.words . P.last $ ls)) ||-     (to == (read . BS.unpack . P.last . BS.words . P.last $ ls))-  then return . P.map (:[]) $ xs-  else fourLines to >>= return . (P.zipWith (:) xs)--fourLines to = do-  I.dropWhile BS.null-  mp <- I.peek-  case mp of-    Nothing -> return $ P.replicate 4 []-    Just p-      | "//" `isInfixOf` p-      || "CM" `isInfixOf` p-      || "Query" `isInfixOf` p-      || ">" `isPrefixOf` p-      || "strand" `isInfixOf` p-      || "#" `isPrefixOf` p-      -> return $ P.replicate 4 []-      | otherwise-      -> do ls <- joinI $ I.take 4 stream2stream-            let ws = BS.length . BS.takeWhile isSpace . P.head $ ls-            let cs = BS.length . BS.dropWhile isSpace . P.head $ ls-            let xs = P.map (BS.take cs . BS.drop ws) ls-            fourLines to >>= return . (P.zipWith (:) xs)---- | Convenience function: read all results into a single list.--fromFile :: FilePath -> IO [VerboseHit]-fromFile fp = do-  i <- enumFile 8192 fp . joinI $ eneeVerboseHit stream2list-  run i------ How to use this enumeratee.--{--test = do-  i <- enumFile 8192 "test.vh" $ joinI $ eneeVerboseHit stream2list-  xs <- run i-  P.mapM_ (\x -> print x >> P.putStrLn "\n\n\n") xs-  print $ P.length xs--}

@@ -1,18 +0,0 @@---- | Shared, internal stuff.--module Biobase.Infernal.VerboseHit.Internal where--import Data.ByteString.Char8 as BS------ | State for import and export functions--data AliGo = AliGo-  { aliCM :: ByteString-  , aliScaffold :: ByteString-  , aliStrand :: Char-  , aliAnnotation :: [ByteString]-  } deriving (Show)-

@@ -0,0 +1,251 @@+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE StandaloneDeriving #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE TemplateHaskell #-}+{-# LANGUAGE PackageImports #-}++-- | Infernal CMs.+--+-- TODO order of nucleotides? ACGU?+--+-- 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 "PrimitiveArray" Data.Array.Repa.Index++import Data.Array.Repa.Index as R+import Data.Array.Repa.Shape as R+import Data.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 #-}++  rangeList (sh1 :. StateID n1) (sh2 :. StateID n2) = [sh :. StateID n | sh <- rangeList sh1 sh2, n <- [n1 .. (n1+n2)] ]+  {-# INLINE rangeList #-}+

@@ -0,0 +1,242 @@+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE DoAndIfThenElse #-}+{-# LANGUAGE ViewPatterns #-}+{-# LANGUAGE ViewPatterns #-}+{-# LANGUAGE TemplateHaskell #-}+{-# LANGUAGE PatternGuards #-}+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE NoMonomorphismRestriction #-}+{-# LANGUAGE OverloadedStrings #-}++-- | Parses text-based covariance-model descriptions.++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 Data.PrimitiveArray+import Data.PrimitiveArray.Zero++import Biobase.SElab.CM+import Biobase.SElab.Types+import qualified Biobase.SElab.HMM as HMM+import qualified Biobase.SElab.HMM.Import 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"++lineParser p = CL.head >>= \x -> return . maybe (error "no more input") (either (\e -> error $ show (e,x)) id . AB.parseOnly p) $ x++-- | Top-level parser for Infernal 1.0 and 1.1 human-readable covariance+-- models. Reads all lines first, then builds up the CM.++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"++      , _nodes  = nsMap+      , _states = ssMap++      , _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++      , _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++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++-- | 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.++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)++finishedHeader :: Maybe ByteString -> Bool+finishedHeader (Just x) = go x where+  go "MODEL:" = True+  go "CM" = True+  go _ = False+finishedHeader _ = False++-- | Parses nodes. Will terminate on "//" which ends a CM. The state parser+-- will just peek on "//", not remove it from the stream.+--+-- A node is (node type, node id, set of states)++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++-- | 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.++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)++-- 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+                      }+  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++++{-+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+-}++-- | Determine if a line is a node line ('Just'). If yes, we'll get the node+-- type as string and the node identifier, too.++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++fromFile :: FilePath -> IO [CM]+fromFile fp = do+  runResourceT $ sourceFile fp $= parseCM1x $$ consume++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 ()+

@@ -0,0 +1,74 @@+{-# 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?++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)++makeLenses ''HMM3

@@ -0,0 +1,178 @@+{-# 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 qualified Data.Map as M+import Data.Char (toLower)++import Biobase.SElab.HMM+import Biobase.SElab.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++-- | create associative map of the key/value data.++headerMap xs = M.fromList . P.map f $ xs where+  f = second (BS.dropWhile (==' ')) . BS.span (/=' ')++-- | Parse the two beginning lines.++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 all individual nodes, except the first one, which uses 'parseBegin'.++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)++-- | Read a HMMER negated log-probability.++readNLP :: String -> NegLogProb+readNLP = go where+  go "*" = NLP $ 1/0+  go xs  = NLP . read $ xs++-- | Read the optional COMPO line.++compoLine = do+  Just p <- CL.peek+  case (BS.words p) of+    ("COMPO":xs) -> CL.head >>= \_ -> return $ P.map (NLP . read . BS.unpack) xs+    _ -> return []++-- | Read the alphabet and transition lines.++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)++-- | All the header lines until we see "HMM".++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)++++-- | 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+

@@ -0,0 +1,22 @@+{-# 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

@@ -0,0 +1,59 @@+{-# 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 #-}

@@ -0,0 +1,36 @@+{-# LANGUAGE TemplateHaskell #-}+{-# LANGUAGE PatternGuards #-}++-- | Infernal contains a taxonomy database. This is a simple module reflecting+-- said database.++module Biobase.SElab.Taxonomy where++import Control.Lens+import Data.Char (toLower)+import qualified Data.ByteString.Char8 as BS++import Biobase.SElab.Types++++-- | 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+shortenName (IDD xs)+  | null ws   = IDD xs+  | [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

@@ -0,0 +1,63 @@+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE NoMonomorphismRestriction #-}+{-# LANGUAGE BangPatterns #-}++-- | Iteratee-based importer. Provides a simple "fromFile" function that+-- produces both maps in one pass.++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 Biobase.SElab.Taxonomy+import Biobase.SElab.Types++++parse =  CB.lines+      =$ CL.map (parseOnly mkTaxonomy)+      =$ CL.filter isRight+      =$ CL.map fromRight+      =$ C.zipSinks mapIdTaxonomy mapAcTaxonomy+{-# INLINE parse #-}++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 #-}++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 #-}++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 #-}++fromFile :: String -> IO ( Map (Identification Species) Taxonomy+                         , Map (Accession Species) Taxonomy+                         )+fromFile fname = do+  runResourceT $ CB.sourceFile fname $$ parse+{-# NOINLINE fromFile #-}

@@ -0,0 +1,93 @@+{-# 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)+

@@ -1,5 +1,5 @@ name:           BiobaseInfernal-version:        0.6.2.0+version:        0.7.0.0 author:         Christian Hoener zu Siederdissen maintainer:     choener@tbi.univie.ac.at homepage:       http://www.tbi.univie.ac.at/~choener/@@ -10,7 +10,7 @@ license-file:   LICENSE build-type:     Simple stability:      experimental-cabal-version:  >= 1.4.0+cabal-version:  >= 1.6.0 description:                 Provides import and export facilities for Infernal/Rfam data                 formats. We include Stockholm, CM, verbose Infernal results,@@ -26,68 +26,67 @@                 .                 .                 .-                Changes in 0.6.2.0-                .-                * added CM parsing (implementation and interface subject to-                  change)-                .-                Changes in 0.6.1.0-                .-                * added cmalign results parser-                .-                Changes in 0.6.0.1-                .-                * haddock should finish now+                Changes in 0.7.0.0                 .-                Changes in 0.6.0.0+                * work-in-progress release (some features missing)                 .-                * multiple changes to data representation (mostly newtypes) and-                  documentation+                * working CM parsing                 .-                * load the Rfam.fasta.gz file (and prepare lookup structures)+                * type defns have changed. using phantom types to specify what kind of model we are working with                 .-                * partial biocore integration+                * using conduit instead of iteratee ++ extra-source-files:  library   build-depends:     base >3 && <5,-    attoparsec,-    attoparsec-iteratee,-    biocore,-    bytestring,+    attoparsec          == 0.10.*   ,+    attoparsec-conduit  == 0.5.*    ,+    biocore             == 0.2      ,+    bytestring                      ,+    bytestring-lexing   == 0.4.*    ,+    conduit             == 0.5.*    ,     containers,-    either-unwrap,-    iteratee,-    iteratee-compress,-    transformers,-    tuple,-    vector,-    PrimitiveArray+    either-unwrap       == 1.1      ,+    lens                == 3.*      ,+    primitive           >= 0.5      ,+    transformers        == 0.3.*    ,+    tuple               == 0.2.*    ,+    vector              >= 0.10     ,+    BiobaseXNA          == 0.6.3.*  ,+    PrimitiveArray      == 0.4.0.0    exposed-modules:-    Biobase.Infernal-    Biobase.Infernal.Align-    Biobase.Infernal.Align.Import-    Biobase.Infernal.Clan-    Biobase.Infernal.Clan.Import-    Biobase.Infernal.CM-    Biobase.Infernal.CM.Export-    Biobase.Infernal.CM.Import-    Biobase.Infernal.Hit-    Biobase.Infernal.RfamFasta-    Biobase.Infernal.RfamFasta.Import-    Biobase.Infernal.TabularHit-    Biobase.Infernal.TabularHit.Import-    Biobase.Infernal.Taxonomy-    Biobase.Infernal.Taxonomy.Import-    Biobase.Infernal.Types-    Biobase.Infernal.VerboseHit-    Biobase.Infernal.VerboseHit.Export-    Biobase.Infernal.VerboseHit.Import-    Biobase.Infernal.VerboseHit.Internal+--    Biobase.Infernal+--    Biobase.Infernal.Align+--    Biobase.Infernal.Align.Import+--    Biobase.Infernal.Clan+--    Biobase.Infernal.Clan.Import+    Biobase.SElab.CM+    Biobase.SElab.CM.Import+    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.Taxonomy+    Biobase.SElab.Taxonomy.Import+    Biobase.SElab.Types+--    Biobase.Infernal.VerboseHit+--    Biobase.Infernal.VerboseHit.Export+--    Biobase.Infernal.VerboseHit.Import+--    Biobase.Infernal.VerboseHit.Internal    ghc-options:     -O2 -funbox-strict-fields +source-repository head+  type: git+  location: git://github.com/choener/BiobaseInfernal