diff --git a/Biobase.cabal b/Biobase.cabal --- a/Biobase.cabal +++ b/Biobase.cabal @@ -1,9 +1,9 @@ name: Biobase -version: 0.1.0.1 +version: 0.3.1.0 author: Christian Hoener zu Siederdissen maintainer: choener@tbi.univie.ac.at -homepage: http://www.tbi.univie.ac.at/~choener/Haskell/ -copyright: Christian Hoener zu Siederdissen, 2010 +homepage: http://www.tbi.univie.ac.at/~choener/ +copyright: Christian Hoener zu Siederdissen, 2010,2011 category: Bioinformatics synopsis: Base library for bioinformatics license: GPL-3 @@ -14,49 +14,53 @@ description: Base library for bioinformatics providing the following features: . - * RNA primary and secondary structure + RNA sequences and energy files: . - * Infernal covariance models + * efficient format for RNA sequences, based on the vector package . - * Turner and Vienna energy files + * internal conversion from strings and bytestrings . + * secondary structure manipulation functions, im- and export of Vienna-dotbracket notation . . - RNA sequences and energy files: . - * efficient format for RNA sequences, based on the vector package + Utility classes: . - * import from strings, bytestrings, fasta files (using the bio library) + * algebraic ring class . - * secondary structure manipulation functions, im- and export of Vienna-dotbracket notation + * instances for Gibbs free energy, partition function probabilities, and scores . - * import Turner energy files + * conversion between different entities . - * cf. http://rna.urmc.rochester.edu/NNDB/index.html . - * im- and export of Vienna 2.0 energy files . + Bioinformatics data sources: . + * mostly just imports, some export functions are provided . - Covariance models: + * Fasta (please use the "bio" library by Ketil Malde) . - * import covariance models + * Mafft . - * basic manipulation of covariance models + * RNA secondary structure: Dot-Bracket, CT . - * Stockholm file format manipulation + * ViennaRNA energy files (with default parameters) . + * Turner energy files: cf. . + * RNAstrand database: dot-parenthesis (one file, many entries) . - Utility classes: + * FR3D: basepairs . - * algebraic ring class + * Infernal: covariance models (text-based) . - * instances for Gibbs free energy, partition function probabilities, and scores + * MC-Fold: NCM database . - * conversion between different entities . - * ready for the vector library + . + Important notes: + . + * In general, this library is not unicode safe (and could somebody tell me the reason for encoding DNA/RNA sequences in unicode?) @@ -68,49 +72,75 @@ build-depends: array, base >=4 && <5, + ghc-prim, bytestring, containers, + deepseq, directory >=1, file-embed >=0.0.3, filemanip >=0.3.3 && <0.4, + filepath, + mtl, parsec >=3 && <4, primitive >=0.3 && <0.4, + split >=0.1.2.3, tuple >=0.2.0.1, + utility-ht, vector >=0.7 && <0.8, - split >=0.1.2.3, + zlib, + either-unwrap, HsTools >=0.0.1.1 && <0.0.2, ParsecTools >=0.0.2 && <0.0.3, - PrimitiveArray >=0.0.2.1 && <0.0.3 + PrimitiveArray >=0.0.3.1 && <0.0.4, + vector-read-instances exposed-modules: Biobase.Constants - Biobase.Infernal.CM - Biobase.Infernal.CM.Export - Biobase.Infernal.CM.Import - Biobase.Infernal.Stockholm - Biobase.Infernal.Stockholm.Import + Biobase.DataSource.CT + Biobase.DataSource.CT.Import + Biobase.DataSource.DotParen + Biobase.DataSource.DotParen.Import + Biobase.DataSource.FR3D + Biobase.DataSource.FR3D.Import + Biobase.DataSource.InfernalCM + Biobase.DataSource.InfernalCM.Export + Biobase.DataSource.InfernalCM.Import + Biobase.DataSource.Mafft + Biobase.DataSource.Mafft.Import + Biobase.DataSource.MCFold + Biobase.DataSource.MCFold.Conversions + Biobase.DataSource.MCFold.Import + Biobase.DataSource.Turner + Biobase.DataSource.Turner.Import + Biobase.DataSource.Turner.Tables + Biobase.DataSource.Vienna + Biobase.DataSource.Vienna.Default + Biobase.DataSource.Vienna.Export.ViennaC + Biobase.DataSource.Vienna.Export.ViennaPar + Biobase.DataSource.Vienna.Import.Turner + Biobase.DataSource.Vienna.Import.ViennaPar + Biobase.DataSource.Vienna.Modification.NonStandard + Biobase.DataSource.Vienna.Modification.Partition + Biobase.DataSource.Vienna.Modification.Temperature Biobase.RNA + Biobase.RNA.Complex + Biobase.RNA.Hashes + Biobase.RNA.NucBounds + Biobase.RNA.Pairs + Biobase.RNA.ViennaPair Biobase.Structure Biobase.Structure.DotBracket Biobase.Structure.Shapes - Biobase.Turner - Biobase.Turner.Import - Biobase.Turner.Tables Biobase.Types.Convert Biobase.Types.Energy Biobase.Types.Partition Biobase.Types.Ring Biobase.Types.Score - Biobase.Vienna - Biobase.Vienna.Default - Biobase.Vienna.Export.ViennaC - Biobase.Vienna.Export.ViennaPar - Biobase.Vienna.Import.Turner - Biobase.Vienna.Import.ViennaPar - Biobase.Vienna.Modification.NonStandard - Biobase.Vienna.Modification.Partition - Biobase.Vienna.Modification.Temperature ghc-options: - -O2 + -Odph + + if impl(ghc > 6.13) + ghc-options: + -fllvm diff --git a/Biobase/Constants.hs b/Biobase/Constants.hs --- a/Biobase/Constants.hs +++ b/Biobase/Constants.hs @@ -6,27 +6,29 @@ --- * energy "Infinities" for Int's +-- * energy "Infinities" -- | maxBound on Int is problematic, if we try something like maxBound + -- maxBound which happens very often. 'eInf' is used for settings like lookup -- tables, etc while comparisons are against 'eMax'. -- Both values need to be small enough not to overflow in typical situations. -eInf :: Int -eInf = 10000000 - -eMax :: Int -eMax = eInf `div` 10 - - +class Constants a where + eInf :: a + eMax :: a --- * floating infinity +instance Constants Int where + eInf = 10000000 + eMax = 1000000 + {-# INLINE eInf #-} + {-# INLINE eMax #-} --- here, we choose a real infinity and let the ieee handle it +instance Constants Double where + eInf = 10000000 + eMax = 1000000 + {-# INLINE eInf #-} + {-# INLINE eMax #-} -fInf :: Double -fInf = 1/0 -- * some constants, required for scaling and Boltzmann diff --git a/Biobase/DataSource/CT.hs b/Biobase/DataSource/CT.hs new file mode 100644 --- /dev/null +++ b/Biobase/DataSource/CT.hs @@ -0,0 +1,38 @@ + +-- | Import CT type structures + +module Biobase.DataSource.CT where + +import qualified Data.Array.IArray as A +import Data.List.Split + + + +-- | A CT file. Encodes one or more individual RNAs within one complex. The +-- 'nucs' array encodes (1,n) as '0' has a special meaning. It is used in +-- 'prv', 'nxt' do denote no next or previous nucleotide in the chain. + +data CT = CT + { comments :: [String] + , lengthLine :: String + , nucs :: A.Array Int Nuc + } + deriving (Show,Read) + +-- | A nucleotide. Encodes character, previous and next nucleotide on the chain +-- and pairing partner. + +data Nuc = Nuc + { nuc :: Char + , prv :: Int + , nxt :: Int + , partner :: Int + } + deriving (Show,Read) + +-- | transforms a 'CT' object into our internal format. + +fromCT :: CT -> () +fromCT ct = res where + res = () + chns = split (keepDelimsR $ whenElt ((==0) . nxt)) . A.elems $ nucs ct diff --git a/Biobase/DataSource/CT/Import.hs b/Biobase/DataSource/CT/Import.hs new file mode 100644 --- /dev/null +++ b/Biobase/DataSource/CT/Import.hs @@ -0,0 +1,64 @@ + +-- | Parse a CT file. + +module Biobase.DataSource.CT.Import where + +import Control.Applicative +import Text.Parsec.Char +import Text.Parsec.Combinator hiding (optional) +import Text.Parsec.Error +import Text.Parsec hiding (many) +import Text.Parsec hiding ((<|>), many, optional) +import Text.Parsec.String +import qualified Data.Array.IArray as A + +import Biobase.RNA +import Biobase.RNA.Pairs +import Biobase.DataSource.CT + + + +-- | Parse a CT file. + +fromFile :: FilePath -> IO (Either String CT) +fromFile fname = do + x <- parseFromFile ctP fname + return $ case x of + Left err -> Left $ show err + Right res -> res + +-- * Parser + +commentP :: GenParser Char st String +commentP = char '#' *> anyChar `manyTill` newline + +commentsP = many commentP + +lenlineP :: GenParser Char st (Int,String) +lenlineP = (,) <$> snumP <*> anyChar `manyTill` newline + +-- | We parse six fields: position, nucleotide, next nuc, previous nuc, pairing +-- partner and (position again). CT8 adds two additional fields (unknown) and +-- (unknown). Fields in brackets () are not returned, just parsed. + +nucLineP :: GenParser Char st (Int,Nuc) +nucLineP = f <$> snumP <*> scharP <*> snumP <*> snumP <*> snumP <* snumP <* choice [newline, snumP *> snumP *> newline] where + f pos chr prv nxt partner = (pos, Nuc chr prv nxt partner) + +nucLinesP :: GenParser Char st [(Int,Nuc)] +nucLinesP = many (try nucLineP) + +-- | Parse a CT file with some error handling. + +ctP :: GenParser Char st (Either String CT) +ctP = f <$> commentsP <*> lenlineP <*> nucLinesP <* spaces <* eof where + f cs (l,ll) ns + | l /= length ns = Left $ "CT inscribed length of: " ++ show l ++ " doesn't match read entries in: " ++ unlines cs + | otherwise = Right . CT cs ll $ A.array (1,length ns) ns + +scharP = spaces *> anyChar + +snumP :: GenParser Char st Int +snumP = read <$ try spaces <*> numberP + +numberP = many1 digit diff --git a/Biobase/DataSource/DotParen.hs b/Biobase/DataSource/DotParen.hs new file mode 100644 --- /dev/null +++ b/Biobase/DataSource/DotParen.hs @@ -0,0 +1,18 @@ + +-- | A dot-parenthesis file. These are mostly used by the RNAstrand DB. +-- Admittedly, the format is very similar to dot-bracket. + +module Biobase.DataSource.DotParen where + +import qualified Data.ByteString.Lazy.Char8 as BS + + + +data DotParen = DotParen + { comments :: [BS.ByteString] + , sequence :: BS.ByteString + , dotparen :: BS.ByteString + } + deriving (Show) + +type DotParens = [DotParen] diff --git a/Biobase/DataSource/DotParen/Import.hs b/Biobase/DataSource/DotParen/Import.hs new file mode 100644 --- /dev/null +++ b/Biobase/DataSource/DotParen/Import.hs @@ -0,0 +1,113 @@ +{-# LANGUAGE RecordWildCards #-} + +-- | Imports a dot-parenthesis file. +-- +-- TODO This importer is not unicode-safe. It silently filters out all +-- sequence/structure pairs with unequal lengths. This is work for later. + +module Biobase.DataSource.DotParen.Import + ( fromFile + , fromFileOrFail + ) where + +import Control.Applicative +import Control.Monad (when) +import Data.Either.Unwrap (isLeft,fromLeft,fromRight) +import qualified Data.ByteString.Lazy.Char8 as BS +import System.Exit (exitFailure) +import Text.Parsec.ByteString.Lazy +import Text.Parsec hiding (many,(<|>)) +import Data.List (findIndices) + +import Biobase.DataSource.DotParen + + + +-- | Import a set of sequence/structure pairs from a file. +-- +-- TODO improve error handling + +fromFile :: FilePath -> IO (Either String DotParens) +fromFile fname = do + res <- parseFromFile pFile fname + check res + where + check res + | Left err <- res = return $ Left (show err) + {- + | Right ls <- res + , xs@(x:_) <- findIndices unEqualLengths ls = return . Left $ "unequal sequence / structure length for: " ++ show (map (ls!!) xs) + -} + | Right ls <- res = return . Right . filter (not . unEqualLengths) $ ls + unEqualLengths DotParen{..} = BS.length sequence /= BS.length dotparen + +-- | Read RnaStrand data or fail + +fromFileOrFail fname = do + xs' <- fromFile fname + when (isLeft xs') $ do + print "error reading source:" + putStrLn $ fromLeft xs' + exitFailure + return $ fromRight xs' + + + +-- * The parsers + +-- | parse a comment line + +pComment :: GenParser Char st BS.ByteString +pComment = BS.pack <$ char '#' <*> manyTill anyChar newline + +-- | parse all comments + +pComments :: GenParser Char st [BS.ByteString] +pComments = many pComment + +data DS = D BS.ByteString | S BS.ByteString +isS (S _) = True +isS _ = False +fromDS (S x) = x +fromDS (D x) = x + + +-- | parse a sequence line + +pSequence :: GenParser Char st DS +pSequence = (\c cs -> S $ BS.pack $ c:cs) <$> noneOf "#" <*> manyTill anyChar newline where + +-- | parse one dot-paren line + +pDotParen :: GenParser Char st DS +pDotParen = D . BS.pack <$> manyTill dps newline where + dps = oneOf "()[]<>{}." + +-- | parse either a dot-paren line or a sequence line + +pDS :: GenParser Char st DS +pDS = choice [try pDotParen, pSequence] + +-- | Returns a pair of (sequence,dot-paren) + +pDSs :: GenParser Char st (BS.ByteString,BS.ByteString) +pDSs = f <$> many1 pDS where + f xs = let + (ss,ds) = span isS xs + in (BS.concat $ map fromDS ss,BS.concat $ map fromDS ds) + +-- | Parse one block + +pBlock :: GenParser Char st DotParen +pBlock = (\cs (ss,ds) -> DotParen cs ss ds) <$> pComments <* many newline <*> pDSs + +-- | Parse many blocks + +pBlocks :: GenParser Char st DotParens +pBlocks = many pBlock + +-- | Parse file + +pFile :: GenParser Char st DotParens +pFile = pBlocks <* eof + diff --git a/Biobase/DataSource/FR3D.hs b/Biobase/DataSource/FR3D.hs new file mode 100644 --- /dev/null +++ b/Biobase/DataSource/FR3D.hs @@ -0,0 +1,120 @@ +{-# LANGUAGE RecordWildCards #-} + +-- | Holds the 'FR3DBasepairs' structure extracted from a FR3D basepairs file. +-- The other files are still on the TODO list +-- +-- TODO write data structures and parsers for the other files. +-- +-- TODO rewrite by reducing all those data-structures to something more +-- intelligent. Use 'Nucleotide's for encoding. Have a converter to 'Complex', +-- too. + +module Biobase.DataSource.FR3D where + +import Data.Map (Map) +import qualified Data.Map as M +import Data.Char (isDigit) +import Control.Monad.Error -- TODO GHC < 7.0.0 only ??? +import Data.List (nubBy) + +import qualified Biobase.RNA.Complex as C +import Biobase.RNA.Pairs (ExtPair) + + + +-- * FR3D data structures + +-- | Wrap everything up, giving us our data source. + +data FR3D = FR3D + { complex :: C.Complex ExtPair -- ^ parsed PDB/FR3D as normal 'Complex' + , fr3dBasepairs :: FR3DBasepairs -- ^ the source information + } deriving (Show) + +-- | A complete 'FR3DBasepairs' structure. We have the 'pdbID' from which the FR3D +-- came from, the 'versions' of different programs, the 'chains' or different +-- nucleotide sequences used in the structure, the 'pairs' themselves and other +-- 'comments' which could not be identified correctly. + +data FR3DBasepairs = FR3DBasepairs + { pdbID :: String + , versions :: String + , chains :: Map String String + , pairs :: [Pair] + , comments :: [String] -- unknown comments, maybe really just comments + } deriving (Read,Show) + +-- | Each 'Pair' is a triple of the 'interaction' between the two 'Nucleotide's +-- 'nuc1' and nuc2'. + +data Pair = Pair + { interaction :: Interaction + , nuc1 :: Nucleotide + , nuc2 :: Nucleotide + } deriving (Read,Show,Eq,Ord) + +-- | Each 'Nucleotide' is described by the 'nuc' itself (A,C,G,U) (something +-- else?), the pdb 'num'ber (whatever that is), the 'chain' id and the +-- 'pos'ition in the chain. + +data Nucleotide = Nucleotide + { nuc :: Char + , num :: String -- should be "20" but could be "20A" (grr) + , chain :: String -- could be "chain 1" or "chain E" or so + , pos :: Int + } deriving (Read,Show,Eq,Ord) + +type Interaction = String -- ? + + + +-- | fix the 'pos' in all 'Nucleotide's. Not actually clear, if this /DOES/ fix +-- the problem! + +fixLengths :: FR3DBasepairs -> FR3DBasepairs +fixLengths b = b {pairs = map f $ pairs b} where + f p = p{nuc1 = g $ nuc1 p, nuc2 = g $ nuc2 p} + g n + | pos n > (length $ chains b M.! chain n) = n{pos = read $ filter isDigit $ num n} + | otherwise = n + +-- | Check if a given 'FR3DBasepairs' structure is correct, or if there are errors + +checkFR3DBasepairs :: FR3DBasepairs -> Either String FR3DBasepairs +checkFR3DBasepairs b' = checkChains b' >>= checkLength >>= return where + checkChains b + | null cnucs = Right b + | otherwise = Left $ "CHAIN: " ++ show cnucs + where cnucs = filter (\n -> not $ chain n `M.member` chains b) ps + checkLength b + | null cnucs = Right b + | otherwise = Left $ "LENGTH: " ++ show b ++ show cnucs + where cnucs = filter (\n -> length (chains b M.! chain n) < pos n) ps -- positions start countaing at "1" + checkNucs b + | null cnucs = Right b + | otherwise = Left $ "NUCS: " ++ show cnucs + where cnucs = filter (\n -> (chains b M.! chain n !! (pos n -1)) == nuc n) ps + ps = concatMap (\x -> [nuc1 x, nuc2 x]) $ pairs b' + + + +-- * Transform 'FR3DBasepairs' into 'Complex' + +-- | Aquire a 'Complex' from the data read. +-- +-- TODO write me + +makeComplex :: FR3DBasepairs -> C.Complex ExtPair +makeComplex b@FR3DBasepairs{..} + | length sqnc /= length chns = error $ "makeComplex: " ++ show b + | otherwise = C.mkComplex "TODO write me" {- sqnc strc chns -} where + sqnc = concatMap snd $ M.toAscList chains + strc = nubBy eqT12 $ map (\(Pair z n1 n2) -> (pos n1-1,pos n2-1,z)) pairs + chns = concatMap (\(k,v) -> concat $ replicate (length v) k) $ M.toAscList chains + +eqT12 (a,b,_) (x,y,_) = (a,b)==(x,y) || (a,b)==(y,x) + +-- | Simple wrapper fixing up an 'FR3D' + +mkFR3D :: FR3DBasepairs -> FR3D +mkFR3D b = FR3D {complex=makeComplex b, fr3dBasepairs=b} diff --git a/Biobase/DataSource/FR3D/Import.hs b/Biobase/DataSource/FR3D/Import.hs new file mode 100644 --- /dev/null +++ b/Biobase/DataSource/FR3D/Import.hs @@ -0,0 +1,92 @@ + +-- Import / parse FR3D basepairs files. + +module Biobase.DataSource.FR3D.Import where + +import Control.Applicative +import Control.Monad +import Data.Either +import Data.List +import System.Directory +import System.FilePath +import System.IO.Unsafe +import Text.Parsec.Char +import Text.Parsec.Combinator hiding (optional) +import Text.Parsec.Error +import Text.Parsec hiding ((<|>), many, optional) +import Text.Parsec.String +import qualified Data.Map as M +import Control.Arrow (first) + +import Biobase.DataSource.FR3D + +-- * New default importer + +-- | import a number of 'FR3D' structures. + +fromDir :: FilePath -> IO ([FR3D], [ParseError]) +fromDir dirname = do + (bps,errs) <- importFR3DBasepairsFiles dirname + return (map mkFR3D bps,errs) + +fromFile :: FilePath -> IO (Either ParseError FR3D) +fromFile fname = do + bp <- importFR3DBasepairsFile fname + return $ mkFR3D `fmap` bp + + + +-- * Older functions + +-- | Parse a complete FR3D basepairs file and return a 'FR3DBasepairs' structure. +-- Does not set the 'pdbID', that one is encoded in the filename. (And in each +-- individual 'Pair'...). + +parseFR3DBasepairs :: GenParser Char st FR3DBasepairs +parseFR3DBasepairs = mk <$> comments <* pdbdesc <*> many (parsePair <* newline) where + comments = many commentline + commentline = char '#' *> manyTill anyChar newline + pdbdesc = string "PDB_ID" *> manyTill anyChar newline + mk cs bps = + let -- TODO this is too fragile! + (v:_:chns) = cs + f (x:y:_) = (x,y) + in FR3DBasepairs "" v (M.fromList $ map (f . words) chns) {-(map ((\[n,c] -> (read n,c)) . drop 1 . words) chns)-} bps cs + +-- | One line of connected nucleotides forms a 'Pair' + +parsePair :: GenParser Char st Pair +parsePair = Pair <$ mtas <*> interact <*> parseNucleotide <* space <*> parseNucleotide where + interact = mtas + +-- | One of the two 'Nucleotide's in a 'Pair' + +parseNucleotide :: GenParser Char st Nucleotide +parseNucleotide = Nucleotide <$> anyChar <* space <*> mtas <*> mtas <*> number where + number = read <$> many1 digit + +-- | Given a directory, try to import all files matching "_basepairs_FR3D.txt". +-- In the result (as,bs), the as are the successful parses, while bs failed. +-- For now (NOV 2010), bs is empty ;-) + +importFR3DBasepairsFiles :: FilePath -> IO ([FR3DBasepairs], [Text.Parsec.Error.ParseError]) +importFR3DBasepairsFiles d = do + de <- doesDirectoryExist d + if de + then do + cnts <- getDirectoryContents d + let cands = filter ("basepairs_FR3D.txt" `isSuffixOf`) cnts + (psL,psR) <- partitionEithers `fmap` mapM (pFF d) cands + return (psR, psL) + else return ([], []) + +pFF d c = unsafeInterleaveIO . fmap (fmap (\x -> x{pdbID = take 4 c})) . parseFromFile parseFR3DBasepairs $ d c + +-- | Try to import just one basepair file. + +importFR3DBasepairsFile :: FilePath -> IO (Either Text.Parsec.Error.ParseError FR3DBasepairs) +importFR3DBasepairsFile f = fmap (fmap (\x -> x{pdbID = take 4 f})) $ parseFromFile parseFR3DBasepairs f + +-- * helper functions + +mtas = manyTill anyChar space diff --git a/Biobase/DataSource/InfernalCM.hs b/Biobase/DataSource/InfernalCM.hs new file mode 100644 --- /dev/null +++ b/Biobase/DataSource/InfernalCM.hs @@ -0,0 +1,326 @@ + +-- +-- Stores Infernal Covariance Model data. Built for ease of use, not speed but +-- should work reasonably well. +-- + +-- TODO Data.Vector? +-- +-- TODO use generics? +-- +-- TODO functions to change to probabilities! +-- +-- TODO cmCanonize function! +-- +-- TODO put functions into their own modules, a bit of cleanup +-- +-- TODO add functions to insert a new node between two already existing nodes; think about how to handle BIF +-- +-- TODO add ability to remove node; think about how to handle BIF +-- +-- NOTE maybe BIF should not be insertable/removable right now? + +module Biobase.DataSource.InfernalCM where + + +import Data.Array.IArray +import Data.List (genericLength) + +import Biobase.RNA hiding (nucE) -- "E type" nucleotides do not happen in CMs! + + + +-- * Data types for Covariance Models + +-- {{{ Data types + +-- | A complete covariance model. Each node and each state can be tagged with +-- additional data. Typically, say after parsing, the tag will be (). + +data CM n s = CM + { nodes :: Array Int (Node n) + , states :: Array Int (State s) + , header :: [(String,String)] -- keeps the list of header entries sorted! + , localBegin :: Array Int Double + , localEnd :: Array Int Double + , cmType :: CMType + , nullModel :: Array Nucleotide Double + } deriving (Show) + +-- | Describes one node + +data Node n = Node + { nid :: Int + , ntype :: NodeType + , nparents :: [Int] -- TODO can there be more than one? + , nchildren :: [Int] + , nstates :: [Int] + , ntag :: n + } deriving (Show) + +-- | One state + +data State s = State + { sid :: Int + , stype :: StateType + , snode :: Int + , sparents :: [Int] + , schildren :: [Transition] + , semission :: [Emission] + , stag :: s + } deriving (Show) + +-- | CMType is important if we want to set localBegin / localEnd! + +data CMType = CMProb | CMScore + deriving (Show,Eq) + +-- | can emit either one nucleotide or a pair + +data Emission + = EmitS {eNuc :: Nucleotide, escore :: Double} + | EmitP { eNucL :: Nucleotide, eNucR :: Nucleotide, escore :: Double} + deriving (Show) + +-- | branches are transition without attached probability becaue both branches are always taken + +data Transition + = Branch {tchild :: Int} + | Transition {tchild :: Int, tscore :: Double} + deriving (Show) + +-- | the different node types + +data NodeType = MATP | MATL | MATR | BIF | ROOT | BEGL | BEGR | END + deriving (Read,Show,Eq,Ord,Enum,Bounded) + +-- | the different state types + +data StateType = MP | IL | IR | D | ML | MR | B | S | E + deriving (Read,Show,Eq,Ord,Enum,Bounded) + +-- }}} + + +-- * make a local model out of a global one + +-- | generate a local model with local begin prob and local end prob + +cmMakeLocal :: Double -> Double -> CM n s -> CM n s +cmMakeLocal pbegin pend cm = cmMakeLocalBegin pbegin $ cmMakeLocalEnd pend cm + + + +cmMakeLocalBegin :: Double -> CM n s -> CM n s +cmMakeLocalBegin pbegin cm = cm{localBegin = localBegin cm // changes} where + changes = rootS : (start : intern) + rootS = (0, prob2Score 0 1.0) -- root disabled! + start = (head . nstates $ nodes cm ! 1, prob2Score (1-pbegin) 1.0) -- the first state after "root 0" + intern = map (\k -> (sid $ nodeMainState cm k,prob2Score (pbegin / l) 1.0)) nds + nds = filter (localBeginPossible cm) . elems $ nodes cm + l = genericLength nds + + + +-- TODO have to change the transition score, too! + +cmMakeLocalEnd :: Double -> CM n s -> CM n s +cmMakeLocalEnd pend cm = cm{localEnd = localEnd cm // changes} where + changes = map (\k -> (sid $ nodeMainState cm k,prob2Score (pend / l) 1.0)) nds + nds = filter (localBeginPossible cm) . elems $ nodes cm + l = genericLength nds + + + +-- * Transform between score and probability mode + +-- | given a CM in score mode, change it to probability mode + +cmScore2Prob :: CM n s -> CM n s +cmScore2Prob cm' = if cmType cm' == CMProb then cm' else CM + (nodes cm) + (statesScore2Prob cm $ states cm) + (header cm) + (localBeginScore2Prob $ localBegin cm) + (localEndScore2Prob $ localEnd cm) + CMProb + nm + where + nm = amap (flip score2Prob 0.25) $ nullModel cm' + cm = cm' {nullModel = nm} + + +-- | Given a CM in prob mode, change to score mode + +cmProb2Score :: CM n s -> CM n s +cmProb2Score cm' = if cmType cm' == CMScore then cm' else CM + (nodes cm) + (statesProb2Score cm $ states cm) + (header cm) + (localBeginProb2Score $ localBegin cm) + (localEndProb2Score $ localEnd cm) + CMScore + nm + where + nm = amap (flip prob2Score 0.25) $ nullModel cm' + cm = cm' {nullModel = nm} + + + +-- | normalize all PROBabilities in a CM + +cmNormalizeProbabilities :: CM n s -> CM n s +cmNormalizeProbabilities cm + | cmType cm == CMScore = error "cannot normalize score-type CM" + | otherwise = cm -- TODO have to map normalization over all scores! + + + +-- {{{ CM score/prob conversion helpers + +statesScore2Prob :: CM n s -> Array Int (State s) -> Array Int (State s) +statesScore2Prob cm sA = amap f sA where + f s = s {schildren = map fT $ schildren s, semission = map fE $ semission s} + fT b@(Branch _) = b + fT (Transition k v) = Transition k (score2Prob v 1.0) + fE (EmitS k v) = EmitS k (score2Prob v $ nullModel cm ! k) + fE (EmitP k1 k2 v) = EmitP k1 k2 (score2Prob v $ (nullModel cm ! k1) * (nullModel cm ! k2)) + + + +localBeginScore2Prob :: Array Int Double -> Array Int Double +localBeginScore2Prob sA = amap f sA where + f s = score2Prob s 1.0 + + + +localEndScore2Prob :: Array Int Double -> Array Int Double +localEndScore2Prob sA = amap f sA where + f s = score2Prob s 1.0 + + + +statesProb2Score :: CM n s -> Array Int (State s) -> Array Int (State s) +statesProb2Score cm sA = amap f sA where + f s = s {schildren = map fT $ schildren s, semission = map fE $ semission s} + fT b@(Branch _) = b + fT (Transition k v) = Transition k (prob2Score v 1.0) + fE (EmitS k v) = EmitS k (prob2Score v $ nullModel cm ! k) + fE (EmitP k1 k2 v) = EmitP k1 k2 (prob2Score v $ (nullModel cm ! k1) * (nullModel cm ! k2)) + + + +localBeginProb2Score :: Array Int Double -> Array Int Double +localBeginProb2Score sA = amap f sA where + f s = prob2Score s 1.0 + + + +localEndProb2Score :: Array Int Double -> Array Int Double +localEndProb2Score sA = amap f sA where + f s = score2Prob s 1.0 + +-- }}} + + + + +-- * Helper Functions + +-- {{{ helper functions + +-- | extract the main state for each node (eg MP state for MATP node) + +-- TODO shouldn't this just be "head $ nstates n"? + +nodeMainState :: CM n s -> Node n -> State s +nodeMainState cm n = head $ filter ((==st) . stype) ss where + (Just st) = (ntype n) `lookup` nodeMainStateAssocs + ss = map (states cm !) $ nstates n + + +-- | Checks for each node, if it can be target of a local begin. + +localBeginPossible :: CM n s -> Node n -> Bool +localBeginPossible cm n = + if ntype n `elem` okNodes + && (not . any (==0) $ nparents n) -- nodes reachable from "root" (that is: node 1) have handled specially + then True + else False + where + okNodes = [MATP,MATL,MATR,BIF] + + + +-- | Checks for each node if it can lead to a local end. + +localEndPossible :: CM n s -> Node n -> Bool +localEndPossible cm n = + if ntype n `elem` okNodes + && (END /= (ntype $ nodes cm ! (nid n +1))) + then True + else False + where + okNodes = [MATP,MATL,MATR,BEGL,BEGR] + + + +-- | transform scores into probabilities, given a nullmodel for x + +-- TODO quickcheck! + +score2Prob x null + | x == (-1/0) = 0 + | otherwise = exp (x * log 2) * null + + + +-- | back into scores + +prob2Score x null + | x == 0 = (-1/0) + | otherwise = log (x / null) / log 2 + + + +-- | Transform a state, setting probabilities instead of scores. Requires CM +-- knowledge for background model. + +-- TODO actually use the background model + +stateScore2Prob :: CM n s -> State s -> State s +stateScore2Prob cm s = error "implement me" + + + +-- | Transform a state, setting scores instead of probabilities. + +stateProb2Score :: CM n s -> State s -> State s +stateProb2Score cm s = error "implement me" +transitionTargets :: [Transition] -> [Int] +transitionTargets xs = map f xs where + f (Branch x) = x + f (Transition x _) = x + + +nodeMainStateAssocs :: [(NodeType,StateType)] +nodeMainStateAssocs = + [ (MATP, MP) + , (MATL, ML) + , (MATR, MR) + , (BIF, B) + , (ROOT, S) + , (BEGL, S) + , (BEGR, S) + , (END, E) + ] + +-- }}} + + + +-- TODO score -> prob : exp (x / log 2) -- CHECK THIS!!! +-- score -> prob : exp (x * log 2) +-- prob -> score : (log x) / log 2 +-- TODO mkGlobal (needed?) +-- TODO use logsum : log (exp x + exp y) = x + log (1 + exp (y-x)), where x>y diff --git a/Biobase/DataSource/InfernalCM/Export.hs b/Biobase/DataSource/InfernalCM/Export.hs new file mode 100644 --- /dev/null +++ b/Biobase/DataSource/InfernalCM/Export.hs @@ -0,0 +1,76 @@ +{-# LANGUAGE RecordWildCards #-} + +-- | write an Infernal model. We aim for version 1.0.2 of Infernal but right +-- now, there is no explicit version management! We just return header +-- information as we read it. So we are compatible to all versions that can be +-- parsed but have no safeguards in place. +-- +-- NOTE: exported model can apparently be calibrated. +-- +-- TODO make sure that (Export.toString . Import.fromString == id); the other +-- way around is not absolutely necessary (except if the Infernal parser should +-- require this at some point) +-- +-- TODO put functions like stringS into another module? + +module Biobase.DataSource.InfernalCM.Export where + +import qualified Data.Array.IArray as A +import Text.Printf + +import Biobase.DataSource.InfernalCM + +import Biobase.DataSource.InfernalCM.Import -- TODO for testing only! + + + +-- | Main export function for CMs. Creates a string that is accepted by +-- Infernal + +toString :: CM n s -> String +toString cm@CM{..} = unlines $ hdr ++ ["MODEL:"] ++ nds ++ ["//"] + where + hdr = map (\(k,v) -> k ++ " " ++ v) header + nds = concatMap (nodeToString cm) $ A.elems nodes + + + +-- | export a specific node, used by 'toString' + +nodeToString :: CM n s -> Node n -> [String] +nodeToString cm Node{..} = printf "\t\t\t\t[ %s %d ]" (show ntype) nid : map (stateToString cm . (states cm A.!)) nstates + + + +-- | export a specific states, used by 'nodeToString' + +stateToString :: CM n s -> State s -> String +stateToString cm State{..} = printf " %2s %6d %6d %2d %6d %6d %-55s %s" (show stype) sid maxP numP minC numC tscores escores + where + maxP + | null sparents = -1 + | otherwise = maximum sparents + numP = length sparents + minC + | null schildren = -1 + | otherwise = minimum $ map tchild schildren + numC + | stype == B = maximum $ map tchild schildren + | otherwise = length schildren + tscores :: String + tscores + | stype == B = "" + | otherwise = concatMap (printf "%s " . stringS . tscore) schildren + escores :: String + escores + | null semission = "" + | otherwise = concatMap (printf "%.3f " . escore) semission + + + +-- | export scores for transitions with the Infernal-specific "-infinity" value +-- of "*" + +stringS x + | x == (-1/0) = " *" + | otherwise = printf "%.3f" x diff --git a/Biobase/DataSource/InfernalCM/Import.hs b/Biobase/DataSource/InfernalCM/Import.hs new file mode 100644 --- /dev/null +++ b/Biobase/DataSource/InfernalCM/Import.hs @@ -0,0 +1,197 @@ + +-- +-- A parser for Infernal's Covariance Models (CM). Should work with version 1.0 +-- of Infernal. +-- + +-- TODO use generics? +-- +-- TODO importing the cmcalibrate calibration data is currently broken! (data is parsed as separate lines) + +module Biobase.DataSource.InfernalCM.Import + ( fromFile + , fromString + ) where + +import Data.Maybe (fromJust) +import qualified Data.Array.IArray as A +import qualified Data.List as L +import Text.ParserCombinators.Parsec hiding (State) + +import Biobase.DataSource.InfernalCM +import Biobase.RNA hiding (nucE) + + + +-- {{{ Parsec + +myChar = alphaNum <|> char '-' <|> char '_' <|> char '.' +headerChar = myChar <|> char '[' <|> char ']' +floating = many1 $ digit <|> char '.' <|> char 'e' <|> char '-' + +headerLine = do + key <- many1 headerChar + --spaces + many1 $ char ' ' + val <- many1 (headerChar <|> char ' ' <|> char '/' <|> char ':') + newline + return (key,val) + +unknownNumbers = do + many1 $ char ' ' + ns <- floating `sepEndBy` (many1 $ char ' ') + newline + return ("",concat ns) + +theHeader = do + hs <- many $ (try headerLine <|> try unknownNumbers) + string "MODEL:" + newline + return hs + +-- }}} + +-- {{{ Node + +node = do + spaces + string "[ " + name <- many letter + spaces + num <- many digit >>= (return . read) + spaces + string "]" + newline + states <- many $ try (state num) + return $ (Node + { ntype = read name + , nid = num + , nstates = map sid states + , nparents = [] + , nchildren = [] + , ntag = () + } + , states) + +-- }}} + +-- {{{ stuff + +-- TODO put some of this into HsTools/Parsec stuff + +aNum = many1 $ char '-' <|> digit +addneginf :: String -> Double +addneginf "*" = (-1)/0 --"-1000000000.0" +addneginf x = read x + +-- }}} + +-- {{{ State + +state nodeID = do + spaces + name <- many1 myChar + spaces + sid <- aNum >>= (return . read) + spaces + plast <- aNum >>= (return . read) + spaces + pnum <- aNum >>= (return . read) + spaces + cfirst <- aNum >>= (return . read) + spaces + cnum <- aNum >>= (return . read) + probs <- ( (many $ myChar <|> char '*') `sepBy` (many1 $ char ' ') >>= + return . map addneginf . filter ((/=) "") ) + newline + let s = case name of + "B" -> State + { stype = read name + , sid = sid + , snode = nodeID + , sparents = [plast-pnum+1 .. plast] + , schildren = [Branch cfirst, Branch cnum] + , semission = [] + , stag = () + } + _ -> State + { stype = read name + , sid = sid + , snode = nodeID + , sparents = [plast-pnum+1 .. plast] + , schildren = zipWith Transition [cfirst .. cfirst+cnum-1] probs + , semission = let keep = drop cnum probs in + case length keep of + 0 -> [] + 4 -> zipWith EmitS acgu keep + 16 -> zipWith (\(k1,k2) v -> EmitP k1 k2 v) acguPairs keep + _ -> error $ "strange number of probabilities" ++ show (keep) + , stag = () + } + return s + +-- }}} + +-- {{{ Model + +models = do + ms <- many model + eof + return ms + +model = do + h <- theHeader + ns <- many node + let states = concatMap snd ns + let nodes = map (addPCinfo states . fst) ns + -- just add all the node parent / child info + string "//" + newline + -- eof -- removed, we want to be able to read concatenated models! + return $ CM + { nodes = A.array (0, length nodes -1) $ zip (map nid nodes) nodes + , states = A.array (0, length states -1) $ zip (map sid states) states + , header = h + , localBegin = A.array (0, length states -1) $ zip [0 .. length states -1] (0.0 : repeat (-1/0)) + , localEnd = A.array (0, length states -1) $ zip [0 .. length states -1] (repeat (-1/0)) + , cmType = CMScore + , nullModel = A.array (nucA,nucU) $ zip acgu (map read . words . fromJust $ "NULL" `lookup` h) -- TODO circumvents the whole parsing stuff! + } + +-- }}} + +-- {{{ Stuff + +addPCinfo states n = + let + s = nstates n + sp = L.nub $ L.sort $ concatMap (sparents . (states !!)) s + sc = L.nub $ L.sort $ concatMap (transitionTargets . schildren . (states !!)) s + np = (L.nub $ L.sort $ map (snode . (states !!)) sp) L.\\ [nid n] + nc = (L.nub $ L.sort $ map (snode . (states !!)) sc) L.\\ [nid n] + in + n {nparents = np, nchildren = nc} + + + +-- | Two types of parsing, once using a file and once by parsing a string. + +fromFile f = parseFromFile models f +fromString s = parse models "(stdin)" s + +-- }}} + +-- | Helper function to remove impossible state transitions (those that have +-- -infty score). + +-- TODO move to InfernalCM.hs and have it for Prob and Score both + +{- +canonize cm = cm {states = A.amap f $ states cm} where + f s = s {schildren = filter g $ schildren s} + g (Branch _) = True + g (Transition _ v) = + h (_,s) + | s == (-1)/0 = False + | otherwise = True +-} diff --git a/Biobase/DataSource/MCFold.hs b/Biobase/DataSource/MCFold.hs new file mode 100644 --- /dev/null +++ b/Biobase/DataSource/MCFold.hs @@ -0,0 +1,190 @@ +{-# LANGUAGE GeneralizedNewtypeDeriving #-} + +-- | The next iteration of MC-Fold-DP which is, hopefully, a bit faster. It can +-- be seen that the dsConnect and ddConnect data structures do not care about +-- the type of connection, only which kind of NCM is involved and the pair. +-- This follows the way MC-Fold integrates away the specific type of pair +-- family involved. + +module Biobase.DataSource.MCFold where + +import Control.DeepSeq +import Data.Ix +import Data.Ix.Tuple +import Data.PrimitiveArray +import Data.PrimitiveArray.Ix -- only needed for the ghci part? +import qualified Data.Vector.Unboxed as VU + +import Biobase.RNA +import Biobase.RNA.Hashes + + + +-- | This version of the motif database accepts ACGU+E as nucleotides. +-- Occurances of E lead to a maximization over all strings where the E are +-- replaced by ACGU -- if the corresponding function was used, otherwise all E +-- map to infinity (actually only a large number). +-- +-- The dsConnect and ddConnect database stores the combined information for +-- connections between different cycles. + +data MotifDB = MotifDB + -- The most important data, NCMs and connections between two of those + { sCycles :: [(Int,CycleArray)] -- ^ store length of cycles, cycle array for this length + , dCycles :: [((Int,Int),CycleArray)] -- ^ store (left,right) lengths, cycle array + , dsConnect :: [(((Int,Int),Int),ConnectArray)] + , ddConnect :: [(((Int,Int),(Int,Int)),ConnectArray)] + -- Following is other data, mostly from the importer + , rawPairScores :: PairScoreArray + , rawJunctionScores :: JunctionArray + , rawSHingeCounts :: [(Int,[RawHinge])] -- ^ length of cycle, hinge data + , rawDHingeCounts :: [((Int,Int),[RawHinge])] -- ^ length of cycle, length of cycle, hinge data + } deriving (Show) + + + +-- * Newtypes for different annotations + +newtype TransCis = TransCis Int + deriving (Eq,Ord,Ix,NFData,Read,Show,Enum) +(trans:cis:_) = map TransCis [0..] + +instance Bounded TransCis where + minBound = trans + maxBound = cis + +newtype AntiPara = AntiPara Int + deriving (Eq,Ord,Ix,NFData,Read,Show,Enum) +(anti:para:_) = map AntiPara [0..] + +instance Bounded AntiPara where + minBound = anti + maxBound = para + +newtype Slash = Slash Int + deriving (Eq,Ord,Ix,NFData,Read,Show,Enum) +(wsp:hsp:csp:ssp:bsp:_) = map Slash [0..] + +instance Bounded Slash where + minBound = wsp + maxBound = bsp + +newtype Prime = Prime Int + deriving (Eq,Ord,Ix,NFData,Read,Show,Enum) +(p5:p3:_) = map Prime [0..] + +instance Bounded Prime where + minBound = p5 + maxBound = p3 + +-- | log-scores as derived from probabilities. +-- +-- TODO apply (min eInf)? + +logScore :: Double -> Double +logScore x = (-0.606) * log x + +type CycleArray = PrimArray HashedPrimary Double +type ConnectArray = PrimArray (Nucleotide,Nucleotide) Double -- ,Slash,Slash,AntiPara,TransCis) Double +type RawHinge = (((Prime,Nucleotide,Nucleotide),(Slash,Slash,AntiPara,TransCis)),Double) +type PairScoreArray = PrimArray (Nucleotide,Nucleotide) Double +type JunctionArray = PrimArray (Int,Int) Double + + + +-- * A lot of static information + +-- | Known single-pair types. (Length of nucleotide string, 0-based index). + +-- 4 -> (5,2) +-- 5 -> (5,3) +-- 6 -> (5,4) +-- +-- found to be wrong +-- 5 -> (0,5) + +knownSingleNCM :: VU.Vector (Int,Int) +knownSingleNCM = VU.fromList + [ (4,2) + , (5,3) + , (6,4) + ] + +-- | Known double-pair types. ((Length of left...,Length of right nucleotide +-- string), 0-based index + +-- found to be right through hacking! +-- all are (..) to next type +-- NCM type -> junctions IV.! (a,b) +-- (a+1,b+1) for the junction.db +-- (2,2) -> (5,5) +-- (2,3) -> (8,2) +-- (3,2) -> (9,2) + +knownDoubleNCM :: VU.Vector ((Int,Int),Int) +knownDoubleNCM = VU.fromList + [ ((2,2), 5) + , ((3,3), 6) + , ((2,3), 7) + , ((3,2), 8) + , ((2,4), 9) + , ((4,2),10) + ] + +-- | This is a list of known hinge types. If it has to be regenerated at some +-- point, run the following command inside the MCfold DB: +-- +-- cat *hinge | grep -v "^\[" | awk '{print $1 " " $2 " " $3}' | sort | uniq | awk '{printf " , \"%s\"\n", $0}' + +knownHinges = + [ "B/B anti trans" + , "B/H anti trans" + , "B/H para cis" + , "B/S anti cis" + , "B/S anti trans" + , "B/S para trans" + , "B/W anti cis" + , "B/W anti trans" + , "B/W para cis" + , "B/W para trans" + , "C/H para trans" + , "H/B anti trans" + , "H/C para trans" + , "H/H anti cis" + , "H/H para trans" + , "H/S anti trans" + , "H/S para trans" + , "H/W anti cis" + , "H/W anti trans" + , "H/W para cis" + , "H/W para trans" + , "S/B anti cis" + , "S/B para trans" + , "S/C para cis" + , "S/H anti cis" + , "S/H anti trans" + , "S/H para cis" + , "S/H para trans" + , "S/S anti cis" + , "S/S para trans" + , "S/W anti cis" + , "S/W anti trans" + , "S/W para cis" + , "S/W para trans" + , "W/B anti cis" + , "W/B anti trans" + , "W/B para cis" + , "W/C anti trans" + , "W/H anti cis" + , "W/H anti trans" + , "W/H para cis" + , "W/H para trans" + , "W/S anti cis" + , "W/S anti trans" + , "W/S para cis" + , "W/S para trans" + , "W/W anti cis" + , "W/W anti trans" + , "W/W para cis" + , "W/W para trans" + ] diff --git a/Biobase/DataSource/MCFold/Conversions.hs b/Biobase/DataSource/MCFold/Conversions.hs new file mode 100644 --- /dev/null +++ b/Biobase/DataSource/MCFold/Conversions.hs @@ -0,0 +1,58 @@ + +-- | Performs operations on cyclearrays, most importantly calculating energies +-- for "E" cases. +-- +-- TODO MotifDB -> MotifDB should be the type. +-- +-- TODO two different cases (i) activate "E" with reasonable defaults (ii) +-- generic sparse data correction. If (ii) is used, it should come before (i) +-- in order to produce reasonable results for motifs with "E" + +module Biobase.DataSource.MCFold.Conversions + ( maxArray + ) where + +import qualified Data.Vector.Unboxed as VU + +import Data.PrimitiveArray + +import Biobase.Constants +import Biobase.DataSource.MCFold +import Biobase.RNA +import Biobase.RNA.Hashes +import Biobase.RNA.NucBounds + + + +-- * activate + +-- | Given a CycleArray, create entries for all keys containing at least one E. +-- Set said energy to maximum over all non-E containing neighbors. + +maxArray :: CycleArray -> Int -> CycleArray +maxArray arr k = fromAssocs l u eInf xs where + (l,u) = bounds arr + xs = map (\k -> let p = mkPrimary k in (mkHashedPrimary (minExtended,maxExtended) p,maxOver arr p)) $ genKeys k + +-- | Perform maximum calculation for one element + +maxOver :: CycleArray -> Primary -> Double +maxOver arr p + | VU.all (/=nucE) p = arr ! mkci p + | otherwise = maximum $ map ((arr !) . mkci . mkPrimary) keys + where + mkci = mkHashedPrimary (minExtended,maxExtended) + keys = mkKeys $ VU.toList p + +-- | Create all neighboring keys for a given one + +mkKeys :: [Nucleotide] -> [[Nucleotide]] +mkKeys [x] = if x==nucE then [[y] |y<-acgu] else [[x]] +mkKeys (x:xs) = [ head y:ys | y<-mkKeys [x], ys<-mkKeys xs ] + +-- | Generate all keys of a given size + +genKeys :: Int -> [[Nucleotide]] +genKeys k + | k==1 = [[y] | y<-eacgu] + | k> 1 = [ y:ys | y<-eacgu, ys<-genKeys (k-1) ] diff --git a/Biobase/DataSource/MCFold/Import.hs b/Biobase/DataSource/MCFold/Import.hs new file mode 100644 --- /dev/null +++ b/Biobase/DataSource/MCFold/Import.hs @@ -0,0 +1,226 @@ +{-# LANGUAGE ScopedTypeVariables #-} + +-- | Importer for the NCM database. +-- +-- NOTE This version reads only known single/double-stranded NCMs. +-- +-- NOTE "infinity" is a large constant, given by RNA.Constants.eInf + +module Biobase.DataSource.MCFold.Import + ( parseDir + ) where + +import Control.Applicative +import Control.Arrow (first, second) +import Control.Monad (liftM) +import Data.Either.Unwrap (fromRight) +import Data.Maybe (fromJust) +import qualified Data.ByteString as B +import qualified Data.Map as M +import qualified Data.Vector.Unboxed as VU +import System.FilePath (()) +import Text.Parsec.ByteString +import Text.Parsec hiding (many, (<|>)) +import Text.Parsec.Prim (parse) +import Text.Printf + +import Text.Parsec.Numbers +import Data.PrimitiveArray + +import Biobase.Constants +import Biobase.DataSource.MCFold +import Biobase.RNA +import Biobase.RNA.Hashes +import Biobase.RNA.NucBounds + + + +parseDir :: FilePath -> IO MotifDB +parseDir dir = do + let sList = map fst $ VU.toList knownSingleNCM + let dList = map fst $ VU.toList knownDoubleNCM + -- basic filenames + let (spFiles :: [FilePath]) = map ((dir ) . printf "1_strand_%d.db") sList + let (dsFiles :: [FilePath]) = map ((dir ) . uncurry (printf "2_strand_%d_%d.db")) dList + -- parse cycle DB + let parseWith parser suffix lst = liftM (zip lst) . mapM (liftM fromRight . parseFromFile parser . (++suffix)) + sCycles <- parseWith parseCycle ".cycle" sList spFiles + dCycles <- parseWith parseCycle ".cycle" dList dsFiles + -- parse hinge DB: parses the original counted data + sHinges <- parseWith parseHinge ".hinge" sList spFiles + dHinges <- parseWith parseHinge ".hinge" dList dsFiles + -- parse the DB with individual pair contributions, dropping the 5 integer-based values first + pairScores <- liftM (fromAssocs minExtended maxExtended eInf . zip [(n1,n2)|n1<-eacgu,n2<-eacgu] . drop 5 . fromRight) . parseFromFile parsePairScores $ (dir "hinge.db") -- strangely, the data is in hinge.db + -- parse contributions by individual junctions + junctionScores <- liftM (fromList (0,0) (21,21) . fromRight) . parseFromFile parseJunction $ (dir "junction.db") + return $ MotifDB + { sCycles = map (second (toNCMArray pairScores)) sCycles + , dCycles = map (second (toNCMArray pairScores)) dCycles + , dsConnect = [ mkDSConnect junctionScores pairScores d s | d <- dHinges, s <- sHinges ] + , ddConnect = [ mkDDConnect junctionScores pairScores d d | d <- dHinges, d <- dHinges ] + , rawPairScores = pairScores + , rawJunctionScores = junctionScores + , rawSHingeCounts = sHinges + , rawDHingeCounts = dHinges + } + +-- | Create the connections array for a given doubleNCM known by (l,r) length +-- and singleNCM. Such a connection is then only defined by the connecting +-- pair. +-- +-- TODO check if connect calculates correctly! + +mkDSConnect :: JunctionArray -> PairScoreArray -> ((Int,Int),[RawHinge]) -> (Int,[RawHinge]) -> (((Int,Int),Int),ConnectArray) +mkDSConnect junction pair ((lenLeft,lenRight),ds) (len,ss) = (((lenLeft,lenRight),len), fromAssocs l u z xs) where + l = minExtended + u = maxExtended + z = eInf + xs = [ ((n1,n2), connect n1 n2) | n1<-acgu, n2<-acgu ] + dsMap = M.map normalize . M.fromListWith (++) $ map (\((pnn,others),v) -> (pnn,[(others,v)])) $ ds + ssMap = M.map normalize . M.fromListWith (++) $ map (\((pnn,others),v) -> (pnn,[(others,v)])) $ ss + connect n1 n2 = + let + common = [ va*vb | (ka,va) <- M.findWithDefault [] (p3,n1,n2) dsMap, (kb,vb) <- M.findWithDefault [] (p5,n1,n2) ssMap, ka==kb ] + c = if null common + then eInf + else junction ! (((lenLeft,lenRight) `lookup` knownDoubleNCM),(len `lookup` knownSingleNCM)) -- junction score between two NCMs + + (logScore $ sum common) + in c + lookup k = snd . fromJust . VU.find ((==k).fst) + normalize xs + | null xs = xs + | otherwise = map (second (/total)) xs + where total = sum $ map snd xs + +-- | Same as above. Some annoying code duplication. +-- +-- TODO check if connect calculates correctly + +mkDDConnect :: JunctionArray -> PairScoreArray -> ((Int,Int),[RawHinge]) -> ((Int,Int),[RawHinge]) -> (((Int,Int),(Int,Int)),ConnectArray) +mkDDConnect junction pair ((lenLA,lenRA),dsA) ((lenLB,lenRB),dsB) = (((lenLA,lenRA),(lenLB,lenRB)), fromAssocs l u z xs) where + l = minExtended + u = maxExtended + z = eInf + xs = [ ((n1,n2), connect n1 n2) | n1<-acgu, n2<-acgu ] + dsAMap = M.map normalize . M.fromListWith (++) $ map (\((pnn,others),v) -> (pnn,[(others,v)])) $ dsA + dsBMap = M.map normalize . M.fromListWith (++) $ map (\((pnn,others),v) -> (pnn,[(others,v)])) $ dsB + connect n1 n2 = + let + common = [ va*vb | (ka,va) <- M.findWithDefault [] (p3,n1,n2) dsAMap, (kb,vb) <- M.findWithDefault [] (p5,n1,n2) dsBMap, ka==kb ] + c = if null common + then eInf + else junction ! (((lenLA,lenRA) `lookup` knownDoubleNCM),((lenLB,lenRB) `lookup` knownDoubleNCM)) -- junction score between two NCMs + + (logScore $ sum common) + in c + lookup k = snd . fromJust . VU.find ((==k).fst) + normalize xs + | null xs = xs + | otherwise = map (second (/total)) xs + where total = sum $ map snd xs + +-- | Generates the array of all NCMs. We add the score for the closing pair +-- (hinge.db) directly here. + +toNCMArray :: PairScoreArray -> [Cycle] -> PrimArray HashedPrimary Double +toNCMArray pairScores xs + | null xs = error "toNCMArray expects non-empty list of Cycles" + | otherwise = fromAssocs l u z $ map f xs + where + lp = replicate (length . fst $ head xs) (minExtended :: Nucleotide) + up = replicate (length . fst $ head xs) (maxExtended :: Nucleotide) + l = mkci $ mkPrimary lp + u = mkci $ mkPrimary up + z = eInf + mkci s = mkHashedPrimary (minExtended,maxExtended) s + f (s,v) = + let + p = mkPrimary s + pS = pairScores ! (VU.head p, VU.last p) + in (mkci p,v + pS) + + + +-- * Parsers for different files + +-- | parser for junction.db + +parseJunction :: GenParser B.ByteString st [Double] +parseJunction = many1 (doubles <* spaces) <* eof where + doubles = try (const eInf <$> string "+inf") <|> floatX + +-- | parser for hinge.db Returns scores for individual pairs. + +parsePairScores :: GenParser B.ByteString st [Double] +parsePairScores = (\a b -> a ++ b) <$> fl <*> ls where + fl = many1 (read <$> many1 digit <* spaces) -- the 5 integers in the top line + ls = many1 (doubles <* spaces) <* eof -- the 5x5 matrix of doubles + doubles = try (const eInf <$> string "+inf") <|> floatX + +-- | parser for NCM database (*.cycle) + +parseCycle :: GenParser B.ByteString st [Cycle] +parseCycle = many1 (cycle <* newline) <* eof where + cycle :: GenParser B.ByteString st Cycle + cycle = (\a b -> (a,b)) <$> many1 letter <* ws <*> floatX + +-- many ws characters + +ws = many $ char ' ' + +-- | helper function for the strange double values +-- +-- TODO move floatX into HsTools! + +floatX :: GenParser B.ByteString st Double +floatX = toDouble <$> prefix <*> many digit <*> char '.' <*> many digit where + prefix = char '+' <|> char '-' + toDouble '+' ds1 dot ds2 = read $ ds1 ++ [dot] ++ ds2 + toDouble '-' ds1 dot ds2 = read $ "-" ++ ds1 ++ [dot] ++ ds2 + +-- | new hinge parser + +parseHinge :: GenParser B.ByteString st [RawHinge] +parseHinge = concat <$> many hinge <* eof where + mkap "anti" = anti + mkap "para" = para + mktc "trans" = trans + mktc "cis" = cis + mksp x = case x of + 'S' -> ssp + 'H' -> hsp + 'W' -> wsp + 'B' -> bsp + 'C' -> csp + _ -> error $ "parseHinge: " ++ [x] + hinge :: GenParser B.ByteString st [RawHinge] + -- use monad format to check opening / closing equality + hinge = do + string "[ " + n53 <- digit + let primeType = if n53 == '5' then p5 else p3 + string "' " + b1 <- oneOf "ACGU" + string " " + b2 <- oneOf "ACGU" + string " ]" + newline + hts <- many (try hingetype) + return $ map (\((s1,s2,ap,tr),count) -> (((primeType,mkNuc b1,mkNuc b2),(s1,s2,ap,tr)),count)) hts + hingetype :: GenParser B.ByteString st ((Slash,Slash,AntiPara,TransCis),Double) + hingetype = mkht + <$> hw + <* space + <*> (try (string "anti") <|> string "para") + <* space + <*> (try (string "trans") <|> string "cis") + <* space + <*> many1 digit + <* newline + mkht (a:_:b:[]) ap tc cnt = ((mksp a, mksp b, mkap ap, mktc tc), read cnt) + hw = many1 (char '/' <|> letter) + + + +-- * Types + +type Cycle = (String,Double) diff --git a/Biobase/DataSource/Mafft.hs b/Biobase/DataSource/Mafft.hs new file mode 100644 --- /dev/null +++ b/Biobase/DataSource/Mafft.hs @@ -0,0 +1,50 @@ + +-- | Mafft data. +-- +-- NOTE: For the "-" strand, the start position is given by ("genome size" - +-- "start position" + "length") +-- +-- TODO mafft blocks contain "i" parts which i haven't seen yet. These have to +-- be parsed, too. + +module Biobase.DataSource.Mafft where + +import qualified Data.ByteString.Lazy.Char8 as BS + + + +-- | A complete Mafft file. + +data Mafft = Mafft + { comments :: [BS.ByteString] + , blocks :: [Block] + } + deriving (Show) + +-- | One Mafft block. +-- +-- THERE is "a ", "i ", "s " + +data Block = Block + { alignment :: BS.ByteString + , sequences :: [Sequence] + } + deriving (Show) + +-- | Encodes one sequence entry. + +data Sequence = Sequence + { species :: BS.ByteString + , start :: Int + , length :: Int + , strand :: PlusMinus + , genomeSize :: Int + , alignedSequence :: CompressedBS + } + deriving (Show) + +newtype PlusMinus = PlusMinus Char + deriving (Eq,Show) + +newtype CompressedBS = CompressedBS BS.ByteString + deriving (Eq,Show) diff --git a/Biobase/DataSource/Mafft/Import.hs b/Biobase/DataSource/Mafft/Import.hs new file mode 100644 --- /dev/null +++ b/Biobase/DataSource/Mafft/Import.hs @@ -0,0 +1,74 @@ + +-- | Imports a mafft file. The mafft file may be gzip-compressed, as long as +-- the filename contains a ".gz" suffix. Files without ".gz" suffix are treated +-- as uncrompressed. + +module Biobase.DataSource.Mafft.Import + ( fromFile + ) where + +import qualified Codec.Compression.GZip as GZip +import System.FilePath (takeExtension) +import qualified Data.ByteString.Lazy.Char8 as BS +import Text.Parsec hiding (many,(<|>)) +import Text.Parsec.ByteString.Lazy +import Control.Applicative + +import Biobase.DataSource.Mafft + + + +test = fromFile "tba.maf.gz" + +-- | Import 'Mafft' from file. + +fromFile :: FilePath -> IO (Either String Mafft) +fromFile fname = do + content <- if (takeExtension fname == ".gz") + then GZip.decompress `fmap` BS.readFile fname + else BS.readFile fname + res <- return $ runP pFile () fname content + case res of + Left err -> return . Left $ show err + Right m -> return . Right $ m + +-- | Parses one comment line. + +pComment :: GenParser Char st BS.ByteString +pComment = BS.pack <$ char '#' <*> manyTill anyChar newline + +-- | Parses many comment lines. + +pComments :: GenParser Char st [BS.ByteString] +pComments = many pComment + +-- | Parses one complete 'Block'. + +pBlock :: GenParser Char st Block +pBlock = Block <$> pAli <*> many1 pSequence where + pAli = BS.pack <$ char 'a' <*> manyTill anyChar newline "expecting alignment (score)" + +-- | Parses many 'Block's. + +pBlocks :: GenParser Char st [Block] +pBlocks = (try $ pBlock <* newline) `manyTill` (string "##eof maf" <* spaces <* eof) + +-- | The Mafft file parser. + +pFile :: GenParser Char st Mafft +pFile = Mafft <$> pComments <*> pBlocks + +-- | Parses one 'Sequence' line. + +pSequence :: GenParser Char st Sequence +pSequence = Sequence <$ string "s " <*> sname <*> sstart <*> slength <*> sstrand <*> sgenomesize <*> salisequence where + sname = BS.pack <$> manyTill alphaNum (many1 space) + sstart = read <$> number <* spaces + slength = read <$> number <* spaces + sstrand = PlusMinus <$> (char '+' <|> char '-') <* spaces + sgenomesize = read <$> number <* spaces + salisequence = CompressedBS . GZip.compress . BS.pack <$> manyTill anyChar newline + +-- | Parses an integral number. + +number = many1 digit diff --git a/Biobase/DataSource/Turner.hs b/Biobase/DataSource/Turner.hs new file mode 100644 --- /dev/null +++ b/Biobase/DataSource/Turner.hs @@ -0,0 +1,41 @@ +{-# LANGUAGE TypeSynonymInstances #-} + +-- | Mathews / Turner nearest neighbor energy tables for interactions between +-- RNA-RNA and DNA-DNA. Using lazy tables for easy filling of missing data. +-- Transform tables into uvectors for efficient access. + +module Biobase.DataSource.Turner where + +import Biobase.RNA (Nucleotide) +import qualified Biobase.DataSource.Turner.Tables as T + + + +-- | Data structure containing all necessary tables for the Turner 2004 model. +-- Note the temperature, using an additional function, rescaling becomes +-- possible. A Show and a Read instance are automatically derived but one +-- should use the Turner and/or Vienna import/export facilities. All data uses +-- complete tables, where any base can pair with any other. Missing data ("." +-- in tables) is given by 'Nothing', otherwise we have 'Just Int'. +-- +-- In general, have a look here: +-- where +-- parameters are explained. +-- +-- Conversion from Turner to Vienna is only possible in one direction as a +-- number of parameters are not used in Vienna RNA Folding. + + + +-- We are missing tables for DNA-RNA data: helixdr.dat, stackdr.dat, +-- stackdr.dh. Several other files have unknown data, too. These are: +-- dnadynalign*, fam\_hmm\_pars.dat, helix.dat, int22-exp.dh + +type TurnerTables = T.Turner2004 Basepair Nucleotide Double + +type TurnerEntropy = TurnerTables +type TurnerEnthalpy = TurnerTables +type Temperature = Double +type TurnerSet = (Temperature,TurnerEntropy,TurnerEnthalpy) + +type Basepair = (Nucleotide,Nucleotide) diff --git a/Biobase/DataSource/Turner/Import.hs b/Biobase/DataSource/Turner/Import.hs new file mode 100644 --- /dev/null +++ b/Biobase/DataSource/Turner/Import.hs @@ -0,0 +1,176 @@ +{-# LANGUAGE NoMonomorphismRestriction #-} + +-- | Turner file parser. Returns a Turner2004 data structure. We store data in +-- the same way it is stored in the ViennaRNA package. Pairs are tuples +-- however. +-- +-- How is 'stack' data stored: +-- +-- AX +-- UY -> ((A,U),(Y,X)) +-- +-- How 'iloop1x1' is stored: +-- +-- X +-- A G +-- U C -> ((A,U),(C,G),X,Y) +-- Y +-- +-- Now 'iloop1x2' is stored: +-- +-- X +-- A G +-- U C -> ((A,U),(C,G),X,C,Y), single (X) first, then 5' to 3' +-- YC +-- +-- 'iloop2x2' is stored: +-- +-- XY +-- A G +-- U C -> ((A,U),(C,G),X,Y,y,x), X-->Y then x<--y +-- xy +-- +--TODO not sure if dangle3/dangle5 are correctly split or if they should switch + +module Biobase.DataSource.Turner.Import + ( parseTurner + ) where + +import Control.Applicative +import Control.Arrow (first) +import Control.Monad +import Data.Either +import Data.List.Split (splitEvery) +import Data.List (transpose) +import qualified Data.Map as M +import System.FilePath.Find +import Text.Parsec hiding ((<|>), many) +import Text.Parsec.String + +import Biobase.RNA +import Data.Ix.Tuple +import Data.PrimitiveArray +import Data.PrimitiveArray.Ix +import Text.Parsec.Numbers + +import Biobase.DataSource.Turner.Tables + + + +-- * Associate every file with a parse. + +minp = (nucA,nucA) +maxp = (nucU,nucU) +minpb = (minp,nucA) +maxpb = (maxp,nucU) +minpp = (minp,minp) +maxpp = (maxp,maxp) +minppbb = (minp,minp,nucA,nucA) +maxppbb = (maxp,maxp,nucU,nucU) +minppbbb = (minp,minp,nucA,nucA,nucA) +maxppbbb = (maxp,maxp,nucU,nucU,nucU) +minppbbbb = (minp,minp,nucA,nucA,nucA,nucA) +maxppbbbb = (maxp,maxp,nucU,nucU,nucU,nucU) +minpbb = (minp,nucA,nucA) +maxpbb = (maxp,nucU,nucU) +-- (4,3) switched for vienna rna compatibility +keyspp = [((k1,k2),(k4,k3)) | k1 <- acgu, k3 <- acgu, k2 <- acgu, k4 <- acgu] +keyspb = [((k1,k2),k3) | k1 <- acgu, k2 <- acgu, k3 <- acgu] +keyspbb = [((k1,k2),k3,k4) | k1 <- acgu, k3 <- acgu, k2 <- acgu, k4 <- acgu] +plist11 = [(nucA,nucU),(nucC,nucG),(nucG,nucC),(nucU,nucA),(nucG,nucU),(nucU,nucG)] +plist22 = [(nucA,nucU),(nucC,nucG),(nucG,nucC),(nucG,nucU),(nucU,nucA),(nucU,nucG)] +-- (4,3) switched for vienna rna compatibility +keysppbb = [((k1,k2),(k4,k3),k5,k6) | (k1,k2) <- plist11, k5 <- acgu, (k3,k4) <- plist11, k6 <- acgu] +-- (4,3) switched for vienna rna compatibility +keysppbbb = [((k1,k2),(k4,k3),k5,k6,k7) | (k1,k2) <- plist11, k6 <- acgu, k5 <- acgu, (k3,k4) <- plist11, k7 <- acgu] +-- (4,3) switched for vienna rna compatibility, 5786 is 5'3' order top, bottom! +keysppbbbb = [((k1,k2),(k4,k3),k5,k6,k7,k8) | (k1,k2) <- plist22, (k3,k4) <- plist22, k5 <- acgu, k8 <- acgu, k6 <- acgu, k7 <- acgu] +z = 999999 + + + +-- | Given the base dir and a suffix (.dat/.dh most likely), parse the relevant +-- files. The prefix can be used for, eg., dna file + +parseTurner prefix basedir fsuffix = do + vstack <- grabB basedir fsuffix $ prefix ++ "stack" + vdangles <- grabB basedir fsuffix $ prefix ++ "dangle" + vlengths <- grabB basedir fsuffix $ prefix ++ "loop" + vhairmm <- grabB basedir fsuffix $ prefix ++ "tstackh" + viloopmm <- grabB basedir fsuffix $ prefix ++ "tstacki" + viloop23mm <- grabB basedir fsuffix $ prefix ++ "tstacki23" + viloop1nmm <- grabB basedir fsuffix $ prefix ++ "tstacki1n" + vmultimm <- grabB basedir fsuffix $ prefix ++ "tstackm" + vextmm <- grabB basedir fsuffix $ prefix ++ "tstack" + viloop11 <- grabB basedir fsuffix $ prefix ++ "int11" + viloop12 <- grabB basedir fsuffix $ prefix ++ "int21" + viloop22 <- grabB basedir fsuffix $ prefix ++ "int22" + vlookups <- fmap (map (first (map mkNuc)) . concat) $ mapM (\f -> grabA basedir fsuffix $ prefix ++ f) ["triloop","tloop","hexaloop"] + vmisc <- grabB basedir fsuffix $ prefix ++ "miscloop" + let (vdangle3,vdangle5) = splitAt (length vdangles `div` 2) vdangles + let (_:viloopl:vbulgel:vhairpinl:[]) = transpose $ splitEvery 4 vlengths + return $ Turner2004 + { stack = fromAssocs minpp maxpp z $ zip keyspp vstack + , dangle3 = fromAssocs minpb maxpb z $ zip keyspb vdangle3 + , dangle5 = fromAssocs minpb maxpb z $ zip keyspb vdangle5 + , hairpinL = fromAssocs 0 30 z $ zip [1..] vhairpinl + , bulgeL = fromAssocs 0 30 z $ zip [1..] vbulgel + , iloopL = fromAssocs 0 30 z $ zip [1..] viloopl + , hairpinMM = fromAssocs minpbb maxpbb z $ zip keyspbb vhairmm + , iloopMM = fromAssocs minpbb maxpbb z $ zip keyspbb viloopmm + , iloop2x3MM = fromAssocs minpbb maxpbb z $ zip keyspbb viloop23mm + , iloop1xnMM = fromAssocs minpbb maxpbb z $ zip keyspbb viloop1nmm + , multiMM = fromAssocs minpbb maxpbb z $ zip keyspbb vmultimm + , extMM = fromAssocs minpbb maxpbb z $ zip keyspbb vextmm + , hairpinLookup = M.fromList vlookups + , iloop1x1 = fromAssocs minppbb maxppbb z $ zip keysppbb viloop11 + , iloop1x2 = fromAssocs minppbbb maxppbbb z $ zip keysppbbb viloop12 + , iloop2x2 = fromAssocs minppbbbb maxppbbbb z $ zip keysppbbbb viloop22 + , ninio = vmisc !! 2 + , maxNinio = vmisc !! 1 + , multiOffset = vmisc !! 6 + , multiNuc = vmisc !! 7 + , multiHelix = vmisc !! 8 + , largeLoop = vmisc !! 0 + , termAU = vmisc !! 14 + , intermolecularInit = vmisc !! 19 + } + +grabB basedir fsuffix fname = do + p <- parseFromFile pBlockFile (basedir ++ "/" ++ fname ++ fsuffix) + case p of + (Right ans) -> return ans + (Left err) -> error $ show err + +grabA basedir fsuffix fname = do + p <- parseFromFile pAssocFile (basedir ++ "/" ++ fname ++ fsuffix) + case p of + (Right ans) -> return ans + (Left err) -> error $ show err + + + +-- * File parser: We do not want to understand the format, just extract all +-- data. + +-- | Blocks of data. + +pBlockFile = concat <$> many (try numberRow <|> infoRow) <* eof where + -- a line with a least one number + numberRow = many (char ' ') *> (number <* notFollowedBy (char '\'')) `sepEndBy1` (many $ char ' ') <* newline + -- a number is a dot or some extended floating point number + number = (try dotFloat) <|> (999999 <$ char '.') <|> parseExtFloat "floating point number" + +-- | A File with associations. + +pAssocFile = concat <$> many1 (try assocRow <|> infoRow) where + -- one assocs of string, value + assocRow = (\a b -> [(a,b)]) <$ many (char ' ') <*> many1 (oneOf nuc) <* many (char ' ') <*> parseExtFloat <* many (char ' ') + +-- a line with crap +infoRow = [] <$ anyChar `manyTill` newline + +-- a rather stupid float +dotFloat = (\s a b -> read $ s:'0':a:b) <$> (char '-' <|> pure ' ') <*> char '.' <*> many1 digit + +nuc = "ACGUT" diff --git a/Biobase/DataSource/Turner/Tables.hs b/Biobase/DataSource/Turner/Tables.hs new file mode 100644 --- /dev/null +++ b/Biobase/DataSource/Turner/Tables.hs @@ -0,0 +1,133 @@ +{-# LANGUAGE TemplateHaskell #-} +{-# LANGUAGE DeriveFunctor #-} +{-# LANGUAGE RecordWildCards #-} +{-# LANGUAGE TypeSynonymInstances #-} + +-- | Turner Nearest Neighbor Energy Tables} + +-- TODO (Read) instance? (Vector has no Read instance...) + +module Biobase.DataSource.Turner.Tables where + +import qualified Data.Map as M + +import Data.PrimitiveArray +import Data.PrimitiveArray.Ix + + + +-- | A data structure describing all fields as used by the Turner 2004 +-- parameter set. Some fields are commented out until they are being used. + +data Turner2004 a b c = Turner2004 + { stack :: PrimArray (a,a) c + , dangle3 :: PrimArray (a,b) c + , dangle5 :: PrimArray (a,b) c + , hairpinL :: PrimArray Int c + , hairpinMM :: PrimArray (a,b,b) c + , hairpinLookup :: M.Map [b] c +-- , hairpinGGG :: c +-- , hairpinCslope :: c +-- , hairpinCintercept :: c +-- , hairpinC3 :: c + , bulgeL :: PrimArray Int c +-- , bulgeSingleC :: c + , iloop1x1 :: PrimArray (a,a,b,b) c + , iloop1x2 :: PrimArray (a,a,b,b,b) c + , iloop2x2 :: PrimArray (a,a,b,b,b,b) c + , iloopMM :: PrimArray (a,b,b) c + , iloop2x3MM :: PrimArray (a,b,b) c + , iloop1xnMM :: PrimArray (a,b,b) c + , iloopL :: PrimArray Int c + , multiMM :: PrimArray (a,b,b) c + , ninio :: c + , maxNinio :: c + , multiOffset :: c + , multiNuc :: c + , multiHelix :: c +-- , multiAsym :: c +-- , multiStrain :: c + , extMM :: PrimArray (a,b,b) c +-- , coaxStack :: PrimArray (a,a) c +-- , coaxStackOpen :: PrimArray (a,b,b) c +-- , coaxStackCont :: PrimArray (a,b,b) c + , largeLoop :: Double + , termAU :: c + , intermolecularInit :: c + } + + + +-- | Map functions over the payload. + +dmap f Turner2004{..} = Turner2004 + { stack = amap f stack + , dangle3 = amap f dangle3 + , dangle5 = amap f dangle5 + , hairpinL = amap f hairpinL + , hairpinMM = amap f hairpinMM + , hairpinLookup = fmap f hairpinLookup +-- , hairpinGGG = f hairpinGGG +-- , hairpinCslope = f hairpinCslope +-- , hairpinCintercept = f hairpinCintercept +-- , hairpinC3 = f hairpinC3 + , bulgeL = amap f bulgeL +-- , bulgeSingleC = f bulgeSingleC + , iloop1x1 = amap f iloop1x1 + , iloop1x2 = amap f iloop1x2 + , iloop2x2 = amap f iloop2x2 + , iloopMM = amap f iloopMM + , iloop2x3MM = amap f iloop2x3MM + , iloop1xnMM = amap f iloop1xnMM + , iloopL = amap f iloopL + , multiMM = amap f multiMM + , ninio = f ninio + , maxNinio = f maxNinio + , multiOffset = f multiOffset + , multiNuc = f multiNuc + , multiHelix = f multiHelix +-- , multiAsym = f multiAsym +-- , multiStrain = f multiStrain + , extMM = amap f extMM +-- , coaxStack = amap f coaxStack +-- , coaxStackOpen = amap f coaxStackOpen +-- , coaxStackCont = amap f coaxStackCont + , largeLoop = largeLoop + , termAU = f termAU + , intermolecularInit = f intermolecularInit + } + + + +-- | Zip two payloads. +-- +-- TODO right now, we have undefined behaviour when some arrrays are of different length + +dZipWith f t1 t2 = Turner2004 + { stack = zipWithPA f (stack t1) (stack t2) + , dangle3 = zipWithPA f (dangle3 t1) (dangle3 t2) + , dangle5 = zipWithPA f (dangle5 t1) (dangle5 t2) + , hairpinL = zipWithPA f (hairpinL t1) (hairpinL t2) + , hairpinMM = zipWithPA f (hairpinMM t1) (hairpinMM t2) + , hairpinLookup = undefined -- zipWithPA f (hairpinLookup t1) (hairpinLookup t2) + , bulgeL = zipWithPA f (bulgeL t1) (bulgeL t2) + , iloop1x1 = zipWithPA f (iloop1x1 t1) (iloop1x1 t2) + , iloop1x2 = zipWithPA f (iloop1x2 t1) (iloop1x2 t2) + , iloop2x2 = zipWithPA f (iloop2x2 t1) (iloop2x2 t2) + , iloopMM = zipWithPA f (iloopMM t1) (iloopMM t2) + , iloop2x3MM = zipWithPA f (iloop2x3MM t1) (iloop2x3MM t2) + , iloop1xnMM = zipWithPA f (iloop1xnMM t1) (iloop1xnMM t2) + , iloopL = zipWithPA f (iloopL t1) (iloopL t2) + , multiMM = zipWithPA f (multiMM t1) (multiMM t2) + , ninio = f (ninio t1) (ninio t2) + , maxNinio = f (maxNinio t1) (maxNinio t2) + , multiOffset = f (multiOffset t1) (multiOffset t2) + , multiNuc = f (multiNuc t1) (multiNuc t2) + , multiHelix = f (multiHelix t1) (multiHelix t2) + , extMM = zipWithPA f (extMM t1) (extMM t2) + , largeLoop = largeLoop t1 -- no zipping here! TODO make this a constant somewhere! + , termAU = f (termAU t1) (termAU t2) + , intermolecularInit = f (intermolecularInit t1) (intermolecularInit t2) + } + +zipWithPA f a1 a2 = let (l,u) = bounds a1 in fromList l u $ zipWith f (toList a1) (toList a2) diff --git a/Biobase/DataSource/Vienna.hs b/Biobase/DataSource/Vienna.hs new file mode 100644 --- /dev/null +++ b/Biobase/DataSource/Vienna.hs @@ -0,0 +1,96 @@ +{-# LANGUAGE NoMonomorphismRestriction #-} + +-- | Turner energy data as used in the ViennaRNA package. +-- +-- TODO should use energy type for tables, not int! + +module Biobase.DataSource.Vienna where + +import Data.Map (Map,empty) + +import Data.Ix.Tuple +import Data.PrimitiveArray +import Data.PrimitiveArray.Ix + +import Biobase.Constants +import Biobase.DataSource.Turner.Tables +import Biobase.RNA +import Biobase.RNA.NucBounds +import Biobase.RNA.ViennaPair +import Biobase.Types.Energy +import Biobase.Types.Partition + + + +-- Type for Vienna Tables + +-- type ViennaEnergySet = (Temperature, ViennaTables, ViennaTables) + +type ViennaIntTables = Turner2004 ViennaPair Nucleotide Int +type ViennaEnergyTables = Turner2004 ViennaPair Nucleotide Energy +type ViennaPartitionTables = Turner2004 ViennaPair Nucleotide Partition + +type Temperature = Double + + + +-- | An empty Turner2004 set, with Vienna-style keys for tables. Unused fields +-- are commented out. Activate them in BiobaseTurner before doing the same +-- here! + +emptyTables :: ViennaIntTables +emptyTables = Turner2004 + --stack + { stack = emptyIV + --dangle + , dangle3 = emptyPB + , dangle5 = emptyPB + --hairpin + , hairpinL = emptyV + , hairpinMM = emptyIV + , hairpinLookup = empty + --eInf + --eInf + --eInf + --eInf + --bulge + , bulgeL = emptyV + --eInf + --internal + , iloop1x1 = emptyPPBB + , iloop1x2 = emptyPPBBB + , iloop2x2 = emptyPPBBBB + , iloopMM = emptyPBB + , iloop2x3MM = emptyPBB + , iloop1xnMM = emptyPBB + , iloopL = emptyV + --multibranch + , multiMM = emptyPBB + , ninio = eInf + , maxNinio = eInf + , multiOffset = eInf + , multiNuc = eInf + , multiHelix = eInf + --eInf + --eInf + --exterior + , extMM = emptyPBB + --coaxial stacking + --emptyIV + --emptyIV + --emptyIV + --scalar values + , largeLoop = eInf -- large loop is always double + , termAU = eInf + , intermolecularInit = eInf + } + +emptyIV = fromAssocs minExtended maxExtended eInf [] +emptyPB = fromAssocs minExtended (maxBound,nucU) eInf [] +emptyPBB = fromAssocs minExtended (maxBound,nucU,nucU) eInf [] +emptyPPBB = fromAssocs minExtended (maxBound,maxBound,nucU,nucU) eInf [] +emptyPPBBB = fromAssocs minExtended (maxBound,maxBound,nucU,nucU,nucU) eInf [] +emptyPPBBBB = fromAssocs minExtended (maxBound,maxBound,nucU,nucU,nucU,nucU) eInf [] + +-- TODO should be an empty a +emptyV = fromAssocs 0 30 eInf [] diff --git a/Biobase/DataSource/Vienna/Default.hs b/Biobase/DataSource/Vienna/Default.hs new file mode 100644 --- /dev/null +++ b/Biobase/DataSource/Vienna/Default.hs @@ -0,0 +1,27 @@ + +{-# LANGUAGE TemplateHaskell #-} + +-- | Embedding of a default set of parameters. These are encoded in Vienna style and +-- use the Turner2004 set of parameters. + +module Biobase.DataSource.Vienna.Default + ( turner2004GH + ) where + +import Data.FileEmbed +import qualified Data.ByteString.Char8 as B + +import Biobase.DataSource.Vienna.Import.ViennaPar +import Biobase.DataSource.Turner.Tables +import Biobase.DataSource.Vienna +import Biobase.Types.Energy + + + +defaultPars :: B.ByteString +defaultPars = $(embedFile "parfiles/Turner2004GH.par") + +-- Turner entropy and enthalpy tables. + +turner2004GH = (dmap Energy ta,dmap Energy tb) where + (ta,tb) = importPar $ B.unpack defaultPars diff --git a/Biobase/DataSource/Vienna/Export/ViennaC.hs b/Biobase/DataSource/Vienna/Export/ViennaC.hs new file mode 100644 --- /dev/null +++ b/Biobase/DataSource/Vienna/Export/ViennaC.hs @@ -0,0 +1,264 @@ +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TemplateHaskell #-} + +-- | Exports the Vienna energy tables as a number of C files (one .c file, a +-- number .h files). For this reason, files are written into a target +-- directory! + +module Biobase.DataSource.Vienna.Export.ViennaC + ( export + ) where + +import Prelude hiding (pi) +import Text.Printf +import System.Directory +import Data.List.Split +import Data.List +import qualified Data.Map as M +import Data.FileEmbed +import qualified Data.ByteString.Char8 as B + +import Biobase.DataSource.Turner.Tables +import Biobase.DataSource.Vienna +import Data.PrimitiveArray +import Biobase.RNA + +import Biobase.DataSource.Vienna.Default + +(g,h) = turner2004GH + + + +type Directory = FilePath + +-- | Export the files required for default values in the ViennaRNA 2.0 C +-- package. + +export :: ViennaIntTables -> ViennaIntTables -> Directory -> IO () +export trnrG trnrH target = mapM_ (\f -> f trnrG trnrH target) [exportC, export11, export21, export22] + + + +-- | Export the main C file + +exportC :: ViennaIntTables -> ViennaIntTables -> Directory -> IO () +exportC trnrG trnrH target = writeFile (target ++ "/energy_par.c") out where + out = unlines [ B.unpack $(embedFile "templates/cheader") , xs, stck, lps, mms, dngls, tris, tets, hexs,includes] + xs = unlines + [ pd "lxc37" $ largeLoop trnrG + , pi "ML_intern37" $ multiHelix trnrG + , pi "ML_interndH" $ multiHelix trnrH + , pi "ML_closing37" $ multiOffset trnrG + , pi "ML_closingdH" $ multiOffset trnrH + , pi "ML_BASE37" $ multiNuc trnrG + , pi "ML_BASEdH" $ multiNuc trnrH + , pi "MAX_NINIO" $ maxNinio trnrG + , pi "ninio37" $ ninio trnrG + , pi "niniodH" $ ninio trnrH + , pi "TerminalAU37" $ termAU trnrG + , pi "TerminalAUdH" $ termAU trnrH + , pi "DuplexInit37" $ intermolecularInit trnrG + , pi "DuplexInitdH" $ intermolecularInit trnrH + , pi "TripleC37" 100 + , pi "TripleCdH" 1860 + , pi "MultipleCA37" 30 + , pi "MultipleCAdH" 340 + , pi "MultipleCB37" 160 + , pi "MultipleCBdH" 760 + ] + stck = unlines + [ "PUBLIC int stack37[NBPAIRS+1][NBPAIRS+1] =" + , (block2 8 . toList $ stack trnrG) + , "PUBLIC int stackdH[NBPAIRS+1][NBPAIRS+1] =" + , (block2 8 . toList $ stack trnrH) + ] + lps = unlines + [ "PUBLIC int hairpin37[31] = {" ++ (concat . intersperse "," . map printNum . toList $ hairpinL trnrG) ++ "};" + , "PUBLIC int hairpindH[31] = {" ++ (concat . intersperse "," . map printNum . toList $ hairpinL trnrH) ++ "};" + , "PUBLIC int bulge37[31] = {" ++ (concat . intersperse "," . map printNum . toList $ bulgeL trnrG) ++ "};" + , "PUBLIC int bulgedH[31] = {" ++ (concat . intersperse "," . map printNum . toList $ bulgeL trnrH) ++ "};" + , "PUBLIC int internal_loop37[31] = {" ++ (concat . intersperse "," . map printNum . toList $ iloopL trnrG) ++ "};" + , "PUBLIC int internal_loopdH[31] = {" ++ (concat . intersperse "," . map printNum . toList $ iloopL trnrH) ++ "};" + ] + mms = unlines + [ writeMM "I" iloopMM trnrG trnrH + , writeMM "H" hairpinMM trnrG trnrH + , writeMM "M" multiMM trnrG trnrH + , writeMM "1nI" iloop1xnMM trnrG trnrH + , writeMM "23I" iloop2x3MM trnrG trnrH + , writeMM "Ext" extMM trnrG trnrH + ] + dngls = unlines + [ "PUBLIC int dangle3_37[NBPAIRS+1][5] =" + , block2 5 . toList $ dangle3 trnrG + , "PUBLIC int dangle3_dH[NBPAIRS+1][5] =" + , block2 5 . toList $ dangle3 trnrH + , "PUBLIC int dangle5_37[NBPAIRS+1][5] =" + , block2 5 . toList $ dangle5 trnrG + , "PUBLIC int dangle5_dH[NBPAIRS+1][5] =" + , block2 5 . toList $ dangle5 trnrH + ] + t3 = getLookup 3 trnrG trnrH + t4 = getLookup 4 trnrG trnrH + t6 = getLookup 6 trnrG trnrH + tris = writeTabbed 3 "Triloop" t3 + tets = writeTabbed 4 "Tetraloop" t4 + hexs = writeTabbed 6 "Hexaloop" t6 + + includes = unlines $ map ((++ "\"") . ("#include \""++)) + [ "intl11.h" + , "intl11dH.h" + , "intl21.h" + , "intl21dH.h" + , "intl22.h" + , "intl22dH.h" + ] + + + +-- | Export iloops 1x1 + +export11 :: ViennaIntTables -> ViennaIntTables -> Directory -> IO () +export11 trnrG trnrH target = writeFile (target ++ "/intl11.h") outG >> writeFile (target ++ "/intl11dH.h") outH where + outG = unlines + [ "PUBLIC int int11_37[NBPAIRS+1][NBPAIRS+1][5][5] =" + , (block4 8 5 5 . toList $ iloop1x1 trnrG) + ] + outH = unlines + [ "PUBLIC int int11_dH[NBPAIRS+1][NBPAIRS+1][5][5] =" + , (block4 8 5 5 . toList $ iloop1x1 trnrH) + ] + + + +-- | Export iloops 2x1 + +export21 :: ViennaIntTables -> ViennaIntTables -> Directory -> IO () +export21 trnrG trnrH target = writeFile (target ++ "/intl21.h") outG >> writeFile (target ++ "/intl21dH.h") outH where + outG = unlines + [ "PUBLIC int int21_37[NBPAIRS+1][NBPAIRS+1][5][5][5] =" + , (block5 8 5 5 5 . toList $ iloop1x2 trnrG) + ] + outH = unlines + [ "PUBLIC int int21_dH[NBPAIRS+1][NBPAIRS+1][5][5][5] =" + , (block5 8 5 5 5 . toList $ iloop1x2 trnrH) + ] + + + +-- | Export iloops 2x2 + +export22 :: ViennaIntTables -> ViennaIntTables -> Directory -> IO () +export22 trnrG trnrH target = writeFile (target ++ "/intl22.h") outG >> writeFile (target ++ "/intl22dH.h") outH where + outG = unlines + [ "PUBLIC int int22_37[NBPAIRS+1][NBPAIRS+1][5][5][5][5] =" + , (block6 8 5 5 5 5 . toList $ iloop2x2 trnrG) + ] + outH = unlines + [ "PUBLIC int int22_dH[NBPAIRS+1][NBPAIRS+1][5][5][5][5] =" + , (block6 8 5 5 5 5 . toList $ iloop2x2 trnrH) + ] + + + +-- * Helper functions + +writeTabbed (ksize :: Int) key (ns,gs,hs) = unlines + [ printf "PUBLIC char %ss[%d] =" key ((ksize+3)*40+1) + , ns + , ";" + ] ++ printf "PUBLIC int %s37[40] = " key ++ gs ++ ";\n" + ++ printf "PUBLIC int %sdH[40] = " key ++ hs ++ ";\n" + +getLookup k trnrG trnrH = (nsF,gs,hs) where + nsF = concat $ intersperse "\n" $ map ((++ " \"") . (" \"" ++) . map fromNuc) ns + ns = filter ((==(k+2)) . length) . M.keys $ hairpinLookup trnrG + gs = lkup trnrG + hs = lkup trnrH + lkup tbl = ("{"++) . (++"}") . concat . intersperse "," $ map (printNum . (hairpinLookup tbl M.!)) ns + +writeMM typ fun trnrG trnrH = unlines + [ printf "PUBLIC int mismatch%s37[NBPAIRS+1][5][5] =" typ + , block3 5 5 . toList $ fun trnrG + , printf "PUBLIC int mismatch%sdH[NBPAIRS+1][5][5] =" typ + , block3 5 5 . toList $ fun trnrH + ] + +block2 s2 xs + = (++ "};") + . ("{" ++) + . concat + . intersperse "\n," + . map ((++ "}") . ("{" ++) . concat . intersperse ",") + . splitEvery s2 + $ map printNum xs + +block3 s2 s3 xs + = (++ "};") + . ("{" ++) + . concat + . intersperse "\n," + . map ((++ "}") . ("{" ++) . concat . intersperse ",") + . splitEvery s2 + . map ((++ "}\n ") . ("{" ++) . concat . intersperse ",") + . splitEvery s3 + $ map printNum xs + +block4 s2 s3 s4 xs + = (++ "};") + . ("{" ++) + . concat + . intersperse "\n," + . map ((++ "}") . ("{" ++) . concat . intersperse ",") + . splitEvery s2 + . map ((++ "}\n ") . ("{" ++) . concat . intersperse ",") + . splitEvery s3 + . map ((++ "}\n ") . ("{" ++) . concat . intersperse ",") + . splitEvery s4 + $ map printNum xs + +block5 s2 s3 s4 s5 xs + = (++ "};") + . ("{" ++) + . concat + . intersperse "\n," + . map ((++ "}") . ("{" ++) . concat . intersperse ",") + . splitEvery s2 + . map ((++ "}\n ") . ("{" ++) . concat . intersperse ",") + . splitEvery s3 + . map ((++ "}\n ") . ("{" ++) . concat . intersperse ",") + . splitEvery s4 + . map ((++ "}\n ") . ("{" ++) . concat . intersperse ",") + . splitEvery s5 + $ map printNum xs + +block6 s2 s3 s4 s5 s6 xs + = (++ "};") + . ("{" ++) + . concat + . intersperse "\n," + . map ((++ "}") . ("{" ++) . concat . intersperse ",") + . splitEvery s2 + . map ((++ "}\n ") . ("{" ++) . concat . intersperse ",") + . splitEvery s3 + . map ((++ "}\n ") . ("{" ++) . concat . intersperse ",") + . splitEvery s4 + . map ((++ "}\n ") . ("{" ++) . concat . intersperse ",") + . splitEvery s5 + . map ((++ "}\n ") . ("{" ++) . concat . intersperse ",") + . splitEvery s6 + $ map printNum xs + + +pd :: String -> Double -> String +pd k v = printf "PUBLIC double %s=%f;" k v + +pi :: String -> Int -> String +pi k v = printf "PUBLIC int %s=%d;" k v + +printNum :: Int -> String +printNum n + | n >= maxN = " INF" + | otherwise = printf "%6d" n + where + maxN = 10000 diff --git a/Biobase/DataSource/Vienna/Export/ViennaPar.hs b/Biobase/DataSource/Vienna/Export/ViennaPar.hs new file mode 100644 --- /dev/null +++ b/Biobase/DataSource/Vienna/Export/ViennaPar.hs @@ -0,0 +1,140 @@ + +{-# LANGUAGE NoMonomorphismRestriction #-} + +module Biobase.DataSource.Vienna.Export.ViennaPar + ( export + ) where + +import Data.Tuple.Select +import Data.List.Split +import Text.Printf +import Control.Arrow +import Data.List (intersperse) +import qualified Data.Map as M + +import Data.PrimitiveArray +import Biobase.DataSource.Turner.Tables +import Biobase.RNA +import Biobase.RNA.ViennaPair +import Biobase.Constants + +import Biobase.DataSource.Vienna + + + +export :: ViennaIntTables -> ViennaIntTables -> String +export trnr trnrH = hdr ++ blocks ++ mlps ++ ninios ++ misc ++ triloops ++ tetra ++ hexa ++ "\n#END" where + hdr = "## RNAfold parameter file v2.0\n\n" + mlps = printf "# ML_params\n%7d %7d %7d %7d %7d %7d\n\n" + (multiNuc trnr) (multiNuc trnrH) + (multiOffset trnr) (multiOffset trnrH) + (multiHelix trnr) (multiHelix trnrH) + ninios = printf "# NINIO\n%7d %7d %7d\n\n" + (ninio trnr) (ninio trnrH) (maxNinio trnr) + misc = printf "# Misc\n %7d %7d %7d %7d\n\n" + (intermolecularInit trnr) (intermolecularInit trnrH) (termAU trnr) (termAU trnrH) + triloops = "# Triloops\n" ++ printHairpinAssocs 5 trnr trnrH ++ "\n" + tetra = "# Tetraloops\n" ++ printHairpinAssocs 6 trnr trnrH ++ "\n" + hexa = "# Hexaloops\n" ++ printHairpinAssocs 8 trnr trnrH ++ "\n" + blocks = concat $ zipWith (++) + -- entropy terms + [ printBlock "stack" 7 pp2lkey $ stack trnr + , printBlock "mismatch_hairpin" 5 pbb2lkey $ hairpinMM trnr + , printBlock "mismatch_interior" 5 pbb2lkey $ iloopMM trnr + , printBlock "mismatch_interior_1n" 5 pbb2lkey $ iloop1xnMM trnr + , printBlock "mismatch_interior_23" 5 pbb2lkey $ iloop2x3MM trnr + , printBlock "mismatch_multi" 5 pbb2lkey $ multiMM trnr + , printBlock "mismatch_exterior" 5 pbb2lkey $ extMM trnr + , printBlock "dangle5" 5 pb2lkey $ dangle5 trnr + , printBlock "dangle3" 5 pb2lkey $ dangle3 trnr + , printBlock "int11" 5 ppbb2lkey $ iloop1x1 trnr + , printBlock "int21" 5 ppbbb2lkey $ iloop1x2 trnr + , printBlock22 "int22" 4 ppbbbb2lkey $ iloop2x2 trnr + , printLinear "hairpin" 10 $ hairpinL trnr + , printLinear "bulge" 10 $ bulgeL trnr + , printLinear "interior" 10 $ iloopL trnr + ] + -- enthalpy terms + [ printBlockH "stack" 7 pp2lkey $ stack trnrH + , printBlockH "mismatch_hairpin" 5 pbb2lkey $ hairpinMM trnrH + , printBlockH "mismatch_interior" 5 pbb2lkey $ iloopMM trnrH + , printBlockH "mismatch_interior_1n" 5 pbb2lkey $ iloop1xnMM trnrH + , printBlockH "mismatch_interior_23" 5 pbb2lkey $ iloop2x3MM trnrH + , printBlockH "mismatch_multi" 5 pbb2lkey $ multiMM trnrH + , printBlockH "mismatch_exterior" 5 pbb2lkey $ extMM trnrH + , printBlockH "dangle5" 5 pb2lkey $ dangle5 trnrH + , printBlockH "dangle3" 5 pb2lkey $ dangle3 trnrH + , printBlockH "int11" 5 ppbb2lkey $ iloop1x1 trnrH + , printBlockH "int21" 5 ppbbb2lkey $ iloop1x2 trnrH + , printBlock22H "int22" 4 ppbbbb2lkey $ iloop2x2 trnrH + , printLinearH "hairpin" 10 $ hairpinL trnrH + , printLinearH "bulge" 10 $ bulgeL trnrH + , printLinearH "interior" 10 $ iloopL trnrH + ] + + + +-- * Helper functions + +-- | Show the key of the line, minus the changing last key + +showKey :: [(LKey,Int)] -> String +showKey xs = + " /* " ++ + (concat $ intersperse "," $ init $ (map show ps) ++ (map show ns)) ++ + " */" + where + (ps,ns) = fst $ head xs + +-- | Transform from tuple-based keys to a pair of list-based keys. + +type LKey = ([ViennaPair],[Nucleotide]) + +pb2lkey (p1,b1) = ([p1],[b1]) +pbb2lkey (p1,b1,b2) = ([p1],[b1,b2]) +pp2lkey (k1,k2) = ([k1,k2],[]) +ppbb2lkey (p1,p2,b1,b2) = ([p1,p2],[b1,b2]) +ppbbb2lkey (p1,p2,b1,b2,b3) = ([p1,p2],[b1,b2,b3]) +ppbbbb2lkey (p1,p2,b1,b2,b3,b4) = ([p1,p2],[b1,b2,b3,b4]) + +-- | Print a block. + +printBlock = printBlockG noNP where + noNP ((ps,ns),v) = not $ any (==vpNP) ps || any (==nucChain) ns + +printBlockH s = printBlock (s ++ "_enthalpies") + +printBlock22 = printBlockG noNSNPE where + noNSNPE ((ps,ns),v) = not $ any (==vpNP) ps || any (==vpNS) ps || any (==nucE) ns || any (==nucChain) ns + +printBlock22H s = printBlock22 (s ++ "_enthalpies") + +printBlockG fltr s k tolkey xs' = + let + xs = filter fltr $ map (first tolkey) $ assocs xs' + in + printf "# %s\n" s ++ + (concatMap printLine $ splitEvery k xs) ++ + "\n" + +printLine xs = + concatMap printVal xs ++ " " ++ showKey xs ++ + printf "\n" + +printVal (k,v) + | v > eMax = printf " INF" + | otherwise = printf "%6d" v + +-- | A linear block is more boring + +printLinear s k xs' = let xs = assocs xs' in + printf "# %s\n" s ++ + (concatMap (\ys -> concatMap printVal ys ++ "\n") $ splitEvery k xs) ++ + "\n" + +printHairpinAssocs l trnr trnrH = res where + res = concat $ zipWith (\(k,v) vH -> printf "%s %7d %7d\n" (concatMap show k) v vH) xs ys + xs = filter ((==l).length.fst) $ M.assocs $ hairpinLookup trnr + ys = map snd $ filter ((==l).length.fst) $ M.assocs $ hairpinLookup trnrH + +printLinearH s = printLinear (s ++ "_enthalpies") diff --git a/Biobase/DataSource/Vienna/Import/Turner.hs b/Biobase/DataSource/Vienna/Import/Turner.hs new file mode 100644 --- /dev/null +++ b/Biobase/DataSource/Vienna/Import/Turner.hs @@ -0,0 +1,72 @@ +{-# LANGUAGE RecordWildCards #-} + +-- | Import data from a Turner table structure. Those can be imported using the +-- Turner.Import module. + +module Biobase.DataSource.Vienna.Import.Turner + ( fromTurner + ) where + +import Control.Arrow +import Data.Tuple.All +import qualified Data.Map as M + +import Biobase.DataSource.Turner.Tables +import Biobase.DataSource.Turner +import Data.PrimitiveArray +import Biobase.Constants +import Biobase.RNA +import Biobase.RNA.ViennaPair + +import Biobase.DataSource.Vienna +-- | Convert from Turner keys to Vienna keys and all values from Double to Int + +fromTurner :: TurnerTables -> ViennaIntTables +fromTurner Turner2004{..} = Turner2004 + { stack = convert pp stack + , dangle3 = convert pb dangle3 + , dangle5 = convert pb dangle5 + , hairpinL = amap f2i hairpinL + , hairpinMM = convert pbb hairpinMM + , hairpinLookup = M.map f2i hairpinLookup + , bulgeL = amap f2i bulgeL + , iloop1x1 = convert uu iloop1x1 + , iloop1x2 = convert uu iloop1x2 + , iloop2x2 = convert uu iloop2x2 + , iloopMM = convert pbb iloopMM + , iloop2x3MM = convert pbb iloop2x3MM + , iloop1xnMM = convert pbb iloop1xnMM + , iloopL = amap f2i iloopL + , multiMM = convert pbb multiMM + , ninio = f2i ninio + , maxNinio = f2i maxNinio + , multiOffset = f2i multiOffset + , multiNuc = f2i multiNuc + , multiHelix = f2i multiHelix + , extMM = convert pbb extMM + , largeLoop = largeLoop * 100 -- scale to decacals + , termAU = f2i termAU + , intermolecularInit = f2i intermolecularInit + } + + + +-- * Helper functions. + +-- | Vienna values are Int's that have 2-digit precision. + +f2i :: Double -> Int +f2i d = round $ d * 100 + + +-- | Deep breath ;-) + +pp = both mkViennaPair +pb = first mkViennaPair +pbb t = upd1 (mkViennaPair $ sel1 t) t +uu t = upd2 (mkViennaPair $ sel2 t) $ upd1 (mkViennaPair $ sel1 t) t +both f = first f . second f + +-- | Do the conversion, both on keys and values. + +convert fk = fromAssocs minBound maxBound eInf . map (first fk . second f2i) . assocs diff --git a/Biobase/DataSource/Vienna/Import/ViennaPar.hs b/Biobase/DataSource/Vienna/Import/ViennaPar.hs new file mode 100644 --- /dev/null +++ b/Biobase/DataSource/Vienna/Import/ViennaPar.hs @@ -0,0 +1,124 @@ +{-# LANGUAGE PatternGuards #-} +{-# LANGUAGE RecordWildCards #-} + +module Biobase.DataSource.Vienna.Import.ViennaPar + ( importPar + ) where + +import Text.Parsec hiding ((<|>), many, optional) +import Text.Parsec.String +import Text.Parsec.Char +import Text.Parsec.Combinator hiding (optional) +import Control.Applicative +import Data.List +import qualified Data.Map as M + +import Text.Parsec.Numbers +import Biobase.Constants +import Biobase.RNA +import Biobase.RNA.ViennaPair +import Biobase.DataSource.Turner.Tables +import Data.PrimitiveArray + +import Biobase.DataSource.Vienna + + + +-- | Imports Vienna tables from a given string. + +importPar :: String -> (ViennaIntTables,ViennaIntTables) +importPar s + | Right _ <- p = foldr updateTables (trnr,trnrH) bs + | Left err <- p = error $ show err + where + p = parse blocks "(stdin)" s + Right bs = p + (bgs,bhs) = partition (not . isHblock) bs + trnr = foldr updateTable emptyTables bgs + trnrH = foldr updateTable emptyTables bhs + +-- | Apply all parsed blocks to a Turner2004 structure. + +updateTable :: Block -> ViennaIntTables -> ViennaIntTables +updateTable Block{..} trnr -- @Turner2004{..} + | bhdr =? "stack" = trnr{stack = fromAssocs minBound maxBound z $ zip pp ints} + | bhdr =? "dangle3" = trnr{dangle3 = fromAssocs minBound (vpNS,nucU) z $ zip pb ints} + | bhdr =? "dangle5" = trnr{dangle5 = fromAssocs minBound (vpNS,nucU) z $ zip pb ints} + | bhdr =? "hairpin" = trnr{hairpinL = fromList 0 30 ints} + | bhdr =? "mismatch_hairpin" = trnr{hairpinMM = fromAssocs minBound (vpNS,nucU,nucU) z $ zip pbb ints} + | bhdr =? "bulge" = trnr{bulgeL = fromList 0 30 ints} + | bhdr =? "int11" = trnr{iloop1x1 = fromAssocs minBound (vpNS,vpNS,nucU,nucU) z $ zip ppbb ints} + | bhdr =? "int21" = trnr{iloop1x2 = fromAssocs minBound (vpNS,vpNS,nucU,nucU,nucU) z $ zip ppbbb ints} + | bhdr =? "int22" = trnr{iloop2x2 = fromAssocs minBound (vpNS,vpNS,nucU,nucU,nucU,nucU) z $ zip ppbbbb ints} + | bhdr =? "mismatch_interior" = trnr{iloopMM = fromAssocs minBound (vpNS,nucU,nucU) z $ zip pbb ints} + | bhdr =? "mismatch_interior_23" = trnr{iloop2x3MM = fromAssocs minBound (vpNS,nucU,nucU) z $ zip pbb ints} + | bhdr =? "mismatch_interior_1n" = trnr{iloop1xnMM = fromAssocs minBound (vpNS,nucU,nucU) z $ zip pbb ints} + | bhdr =? "interior" = trnr{iloopL = fromList 0 30 ints} + | bhdr =? "mismatch_multi" = trnr{multiMM = fromAssocs minBound (vpNS,nucU,nucU) z $ zip pbb ints} + | bhdr =? "mismatch_exterior" = trnr{extMM = fromAssocs minBound (vpNS,nucU,nucU) z $ zip pbb ints} + | otherwise = trnr -- unknown block! +updateTable ABlock{..} trnr + | otherwise = trnr -- CORRECT! + +-- | update with info where no separation between G/H was done. + +updateTables Block{..} (trnr,trnrH) + | bhdr =? "ML_params" + = ( trnr {multiNuc = ints!!0, multiOffset = ints!!2, multiHelix = ints!!4} + , trnrH{multiNuc = ints!!1, multiOffset = ints!!3, multiHelix = ints!!5} + ) + | bhdr =? "NINIO" + = ( trnr {ninio = ints!!0, maxNinio = ints!!2} + , trnrH {ninio = ints!!1} + ) + | bhdr =? "Misc" + = ( trnr {intermolecularInit = ints!!0, termAU = ints!!2} + , trnrH {intermolecularInit = ints!!1, termAU = ints!!3} + ) + | otherwise = (trnr,trnrH) +updateTables ABlock{..} (trnr,trnrH) + | bhdr =? "Triloops" || bhdr =? "Tetraloops" || bhdr =? "Hexaloops" + = ( trnr {hairpinLookup = hairpinLookup trnr `M.union` (M.fromList $ map (\(k,v,w) -> (s2ns k,v)) ascs)} + , trnrH {hairpinLookup = hairpinLookup trnrH `M.union` (M.fromList $ map (\(k,v,w) -> (s2ns k,w)) ascs)} + ) + | otherwise = (trnr,trnrH) + +pp = [(x,y) | x<-cgnsP,y<-cgnsP] +pb = [(x,y) | x<-cgnsP,y<-eacgu] +pbb = [(x,y,z) | x<-cgnsP,y<-eacgu,z<-eacgu] +ppbb = [(p1,p2,b1,b2) | p1<-cgnsP,p2<-cgnsP,b1<-eacgu,b2<-eacgu] +ppbbb = [(p1,p2,b1,b2,b3) | p1<-cgnsP,p2<-cgnsP,b1<-eacgu,b2<-eacgu,b3<-eacgu] +ppbbbb = [(p1,p2,b1,b2,b3,b4) | p1<-cguaP,p2<-cguaP,b1<-acgu,b2<-acgu,b3<-acgu,b4<-acgu] +z = eInf +-- xs is the block header name, ys the constant string +xs =? ys = xs == ys || xs == ys++"_enthalpies" -- and $ zipWith (==) xs ys +s2ns = map mkNuc + +-- * Simple parser for the vienna 2.0 format. + +isHblock :: Block -> Bool +isHblock b = "enthalpies" `elem` (tails $ bhdr b) + +data Block + = Block {bhdr :: String, ints :: [Int]} + | ABlock {bhdr :: String, ascs :: [(String,Int,Int)]} + deriving (Show) + +blocks :: GenParser Char st [Block] +blocks = vrna2 *> spaces *> (try assocBlock <|> try block) `sepEndBy1` spaces <* string "#END" <* spaces <* eof + +vrna2 = string "## RNAfold parameter file v2.0" + +block :: GenParser Char st Block +block = Block <$> header <* spaces <*> num `sepEndBy1` (try comment <|> spaces) where + num = ((eInf <$ string "INF") <|> parseIntegral) + +assocBlock :: GenParser Char st Block +assocBlock = ABlock <$> header <* spaces <*> asc `sepEndBy1` spaces where + asc = (,,) <$> (many1 $ oneOf "ACGU") <* spaces <*> parseIntegral <* spaces <*> parseIntegral + +header :: GenParser Char st String +header = char '#' *> space *> (alphaNum <|> char '_') `manyTill` newline + +comment :: GenParser Char st () +comment = const () <$ spaces <*> between (string "/*") (string "*/") (many $ noneOf "/*") <* spaces diff --git a/Biobase/DataSource/Vienna/Modification/NonStandard.hs b/Biobase/DataSource/Vienna/Modification/NonStandard.hs new file mode 100644 --- /dev/null +++ b/Biobase/DataSource/Vienna/Modification/NonStandard.hs @@ -0,0 +1,122 @@ + +-- | This module provides facilities to modify values in the tables and return +-- the modified tables. The main reason for the modules' existence are the 'E' +-- nucleotides and 'NS' base pairs used by the ViennaRNA package to denote +-- unknown or missing nucleotides. For each table type, we have a small +-- function. + +module Biobase.DataSource.Vienna.Modification.NonStandard + ( nonStandardMax + ) where + +import Biobase.RNA +import Biobase.RNA.ViennaPair +import Biobase.DataSource.Vienna +import Biobase.DataSource.Turner.Tables +import Data.Ix.Tuple +import Biobase.Constants +import Biobase.Types.Energy + +import Data.PrimitiveArray + + +-- | Each table-changing operation makes the promise that all individual +-- changes in a batch are independent of the changes in the same batch + +changeTable f tbl = fromAssocs l u z . map (f tbl) $ assocs tbl where + (l,u) = bounds tbl + z = Energy eInf + + + +-- |* Replace each occurrance of a value with a strange key according to this +-- rule. The rule is to take the maximum over all sane keys in the +-- neighborhood. Yes, this takes long... + +fPPBBBBmax tbl (k@(p1,p2,b1,b2,b3,b4),v) + | p1==vpNP || p2==vpNP = (k,v) -- don't do anything + | otherwise = (k,maximum' + [ tbl ! (p1n,p2n,b1n,b2n,b3n,b4n) + | p1n <- genKeyP p1 + , p2n <- genKeyP p2 + , b1n <- genKeyB b1 + , b2n <- genKeyB b2 + , b3n <- genKeyB b3 + , b4n <- genKeyB b4 + ]) + +fPPBBBmax tbl (k@(p1,p2,b1,b2,b3),v) + | p1==vpNP || p2==vpNP = (k,v) + | otherwise = (k,maximum' + [ tbl ! (p1n,p2n,b1n,b2n,b3n) + | p1n <- genKeyP p1 + , p2n <- genKeyP p2 + , b1n <- genKeyB b1 + , b2n <- genKeyB b2 + , b3n <- genKeyB b3 + ]) + +fPPBBmax tbl (k@(p1,p2,b1,b2),v) + | p1==vpNP || p2==vpNP = (k,v) + | otherwise = (k,maximum' + [ tbl ! (p1n,p2n,b1n,b2n) + | p1n <- genKeyP p1 + , p2n <- genKeyP p2 + , b1n <- genKeyB b1 + , b2n <- genKeyB b2 + ]) + +fPPmax tbl (k@(p1,p2),v) + | p1==vpNP || p2==vpNP = (k,v) + | otherwise = (k,maximum' + [ tbl ! (p1n,p2n) + | p1n <- genKeyP p1 + , p2n <- genKeyP p2 + ]) + +fPBBmax tbl (k@(p1,b1,b2),v) + | p1==vpNP = (k,v) + | otherwise = (k,maximum' + [ tbl ! (p1n,b1n,b2n) + | p1n <- genKeyP p1 + , b1n <- genKeyB b1 + , b2n <- genKeyB b2 + ]) + +fPBmax tbl (k@(p1,b1),v) + | p1==vpNP = (k,v) + | otherwise = (k,maximum' + [ tbl ! (p1n,b1n) + | p1n <- genKeyP p1 + , b1n <- genKeyB b1 + ]) + +maximum' = Energy . maximum . map unEnergy + + +-- * Generate keys. If the key is of the non-standard kind, then we go fishing, +-- otherwise we keep the key itself. + +genKeyP k = if k==vpNS then cguaP else [k] +genKeyB k = if k==nucE then acgu else [k] + + + +-- | The max operation applied to all relevant tables. (This serves as a +-- pointer, which tables to manipulate). + +nonStandardMax :: ViennaEnergyTables -> ViennaEnergyTables +nonStandardMax tbl = tbl + { iloop1x1 = changeTable fPPBBmax $ iloop1x1 tbl + , iloop1x2 = changeTable fPPBBBmax $ iloop1x2 tbl + , iloop2x2 = changeTable fPPBBBBmax $ iloop2x2 tbl + , hairpinMM = changeTable fPBBmax $ hairpinMM tbl + , iloopMM = changeTable fPBBmax $ iloopMM tbl + , iloop1xnMM = changeTable fPBBmax $ iloop1xnMM tbl + , iloop2x3MM = changeTable fPBBmax $ iloop2x3MM tbl + , multiMM = changeTable fPBBmax $ multiMM tbl + , extMM = changeTable fPBBmax $ extMM tbl + , dangle3 = changeTable fPBmax $ dangle3 tbl + , dangle5 = changeTable fPBmax $ dangle5 tbl + , stack = changeTable fPPmax $ stack tbl + } diff --git a/Biobase/DataSource/Vienna/Modification/Partition.hs b/Biobase/DataSource/Vienna/Modification/Partition.hs new file mode 100644 --- /dev/null +++ b/Biobase/DataSource/Vienna/Modification/Partition.hs @@ -0,0 +1,51 @@ + +-- | Small module for converting Vienna energy tables to partition function +-- value tables. +-- +-- Temperature scaling should happen on the level of energy tables, followed by +-- conversion here. +-- +-- TODO after the switch to the Energy type, change too, the partitionTables +-- generation + +module Biobase.DataSource.Vienna.Modification.Partition + ( partitionTables + ) where + +import Biobase.Types.Partition +import Biobase.Types.Energy +import Biobase.DataSource.Turner.Tables +import Biobase.DataSource.Vienna +import Biobase.Constants + + +-- Testing + +{- +import Biobase.DataSource.Vienna.Default +import Data.PrimitiveArray +import Debug.Trace.Tools +import Biobase.RNA +import Biobase.Types.Ring + +(g,h) = turner2004GH + +-- x = toList $ stack $ partitionTables 37 g +x = {- toList $ -} (stack $ partitionTables 37 g) ! (vpGC,vpGC) +-} + +-- | Create partition tables out of normal energy tables. Note that these +-- values are not normalized. An additional term 1/Z(1/kt) would be required +-- which is missing until we are done calculating the partition function. +-- +-- All energy contributions are rescaled from dekacals back to kcal/mol. +-- +-- We explicitly set all probabilities =0 where the energy is >= eMax. + +partitionTables :: Temperature -> ViennaEnergyTables -> ViennaPartitionTables +partitionTables tempe trnr = dmap f trnr where + kt = (tempe + kelvinC0) * gasconst + -- convert energy into probability; encapsulate as a (logFloat) partition + f (Energy e') = {- traceVal (show (e',e,p)) -} p where + e = fromIntegral e' / 100 + p = if e' >= eMax then Partition 0 else Partition . exp $ -e / kt -- can NOT use mkPartition as these values are NOT NORMALIZED diff --git a/Biobase/DataSource/Vienna/Modification/Temperature.hs b/Biobase/DataSource/Vienna/Modification/Temperature.hs new file mode 100644 --- /dev/null +++ b/Biobase/DataSource/Vienna/Modification/Temperature.hs @@ -0,0 +1,34 @@ + +-- | Given entropy and enthalpy sets, allow modification of the folding +-- temperature (from the default 37 Celsius) + +module Biobase.DataSource.Vienna.Modification.Temperature + ( adjustTemperature + ) where + +import Biobase.DataSource.Vienna +import Biobase.DataSource.Turner.Tables +import Biobase.Constants + + + +-- temperature scaling is a generalized zip: let g and h be a scalar, t the +-- temperature. Then: g_t = h - (h-g)*((t+k_0)/t_m). Where k_0 = 273.15 (Kelvin +-- of 37 Celsius) and t_m = 37+k_0 +-- +-- TODO abstract over table type? +-- +-- TODO this does not take into account that dangles are supposed to _always_ +-- be beneficial. We should therefor transform the energies and then apply 'min +-- 0' onto all dangle types. +-- +-- TODO use EnergyTables (Int is just for im-/export) + +adjustTemperature :: ViennaIntTables -> ViennaIntTables -> Temperature -> ViennaIntTables +adjustTemperature trnrG trnrH cels = dZipWith f trnrG trnrH where + f vG vH = + let g = fromIntegral vG + h = fromIntegral vH + k_0 = kelvinC0 + t_m = 37 + k_0 + in round $ h - (h-g) * ((cels+k_0)/t_m) diff --git a/Biobase/Infernal/CM.hs b/Biobase/Infernal/CM.hs deleted file mode 100644 --- a/Biobase/Infernal/CM.hs +++ /dev/null @@ -1,326 +0,0 @@ - --- --- Stores Infernal Covariance Model data. Built for ease of use, not speed but --- should work reasonably well. --- - --- TODO Data.Vector? --- --- TODO use generics? --- --- TODO functions to change to probabilities! --- --- TODO cmCanonize function! --- --- TODO put functions into their own modules, a bit of cleanup --- --- TODO add functions to insert a new node between two already existing nodes; think about how to handle BIF --- --- TODO add ability to remove node; think about how to handle BIF --- --- NOTE maybe BIF should not be insertable/removable right now? - -module Biobase.Infernal.CM where - - -import Data.Array.IArray -import Data.List (genericLength) - -import Biobase.RNA hiding (nucE) -- "E type" nucleotides do not happen in CMs! - - - --- * Data types for Covariance Models - --- {{{ Data types - --- | A complete covariance model. Each node and each state can be tagged with --- additional data. Typically, say after parsing, the tag will be (). - -data CM n s = CM - { nodes :: Array Int (Node n) - , states :: Array Int (State s) - , header :: [(String,String)] -- keeps the list of header entries sorted! - , localBegin :: Array Int Double - , localEnd :: Array Int Double - , cmType :: CMType - , nullModel :: Array Nucleotide Double - } deriving (Show) - --- | Describes one node - -data Node n = Node - { nid :: Int - , ntype :: NodeType - , nparents :: [Int] -- TODO can there be more than one? - , nchildren :: [Int] - , nstates :: [Int] - , ntag :: n - } deriving (Show) - --- | One state - -data State s = State - { sid :: Int - , stype :: StateType - , snode :: Int - , sparents :: [Int] - , schildren :: [Transition] - , semission :: [Emission] - , stag :: s - } deriving (Show) - --- | CMType is important if we want to set localBegin / localEnd! - -data CMType = CMProb | CMScore - deriving (Show,Eq) - --- | can emit either one nucleotide or a pair - -data Emission - = EmitS {eNuc :: Nucleotide, escore :: Double} - | EmitP { eNucL :: Nucleotide, eNucR :: Nucleotide, escore :: Double} - deriving (Show) - --- | branches are transition without attached probability becaue both branches are always taken - -data Transition - = Branch {tchild :: Int} - | Transition {tchild :: Int, tscore :: Double} - deriving (Show) - --- | the different node types - -data NodeType = MATP | MATL | MATR | BIF | ROOT | BEGL | BEGR | END - deriving (Read,Show,Eq,Ord,Enum,Bounded) - --- | the different state types - -data StateType = MP | IL | IR | D | ML | MR | B | S | E - deriving (Read,Show,Eq,Ord,Enum,Bounded) - --- }}} - - --- * make a local model out of a global one - --- | generate a local model with local begin prob and local end prob - -cmMakeLocal :: Double -> Double -> CM n s -> CM n s -cmMakeLocal pbegin pend cm = cmMakeLocalBegin pbegin $ cmMakeLocalEnd pend cm - - - -cmMakeLocalBegin :: Double -> CM n s -> CM n s -cmMakeLocalBegin pbegin cm = cm{localBegin = localBegin cm // changes} where - changes = rootS : (start : intern) - rootS = (0, prob2Score 0 1.0) -- root disabled! - start = (head . nstates $ nodes cm ! 1, prob2Score (1-pbegin) 1.0) -- the first state after "root 0" - intern = map (\k -> (sid $ nodeMainState cm k,prob2Score (pbegin / l) 1.0)) nds - nds = filter (localBeginPossible cm) . elems $ nodes cm - l = genericLength nds - - - --- TODO have to change the transition score, too! - -cmMakeLocalEnd :: Double -> CM n s -> CM n s -cmMakeLocalEnd pend cm = cm{localEnd = localEnd cm // changes} where - changes = map (\k -> (sid $ nodeMainState cm k,prob2Score (pend / l) 1.0)) nds - nds = filter (localBeginPossible cm) . elems $ nodes cm - l = genericLength nds - - - --- * Transform between score and probability mode - --- | given a CM in score mode, change it to probability mode - -cmScore2Prob :: CM n s -> CM n s -cmScore2Prob cm' = if cmType cm' == CMProb then cm' else CM - (nodes cm) - (statesScore2Prob cm $ states cm) - (header cm) - (localBeginScore2Prob $ localBegin cm) - (localEndScore2Prob $ localEnd cm) - CMProb - nm - where - nm = amap (flip score2Prob 0.25) $ nullModel cm' - cm = cm' {nullModel = nm} - - --- | Given a CM in prob mode, change to score mode - -cmProb2Score :: CM n s -> CM n s -cmProb2Score cm' = if cmType cm' == CMScore then cm' else CM - (nodes cm) - (statesProb2Score cm $ states cm) - (header cm) - (localBeginProb2Score $ localBegin cm) - (localEndProb2Score $ localEnd cm) - CMScore - nm - where - nm = amap (flip prob2Score 0.25) $ nullModel cm' - cm = cm' {nullModel = nm} - - - --- | normalize all PROBabilities in a CM - -cmNormalizeProbabilities :: CM n s -> CM n s -cmNormalizeProbabilities cm - | cmType cm == CMScore = error "cannot normalize score-type CM" - | otherwise = cm -- TODO have to map normalization over all scores! - - - --- {{{ CM score/prob conversion helpers - -statesScore2Prob :: CM n s -> Array Int (State s) -> Array Int (State s) -statesScore2Prob cm sA = amap f sA where - f s = s {schildren = map fT $ schildren s, semission = map fE $ semission s} - fT b@(Branch _) = b - fT (Transition k v) = Transition k (score2Prob v 1.0) - fE (EmitS k v) = EmitS k (score2Prob v $ nullModel cm ! k) - fE (EmitP k1 k2 v) = EmitP k1 k2 (score2Prob v $ (nullModel cm ! k1) * (nullModel cm ! k2)) - - - -localBeginScore2Prob :: Array Int Double -> Array Int Double -localBeginScore2Prob sA = amap f sA where - f s = score2Prob s 1.0 - - - -localEndScore2Prob :: Array Int Double -> Array Int Double -localEndScore2Prob sA = amap f sA where - f s = score2Prob s 1.0 - - - -statesProb2Score :: CM n s -> Array Int (State s) -> Array Int (State s) -statesProb2Score cm sA = amap f sA where - f s = s {schildren = map fT $ schildren s, semission = map fE $ semission s} - fT b@(Branch _) = b - fT (Transition k v) = Transition k (prob2Score v 1.0) - fE (EmitS k v) = EmitS k (prob2Score v $ nullModel cm ! k) - fE (EmitP k1 k2 v) = EmitP k1 k2 (prob2Score v $ (nullModel cm ! k1) * (nullModel cm ! k2)) - - - -localBeginProb2Score :: Array Int Double -> Array Int Double -localBeginProb2Score sA = amap f sA where - f s = prob2Score s 1.0 - - - -localEndProb2Score :: Array Int Double -> Array Int Double -localEndProb2Score sA = amap f sA where - f s = score2Prob s 1.0 - --- }}} - - - - --- * Helper Functions - --- {{{ helper functions - --- | extract the main state for each node (eg MP state for MATP node) - --- TODO shouldn't this just be "head $ nstates n"? - -nodeMainState :: CM n s -> Node n -> State s -nodeMainState cm n = head $ filter ((==st) . stype) ss where - (Just st) = (ntype n) `lookup` nodeMainStateAssocs - ss = map (states cm !) $ nstates n - - --- | Checks for each node, if it can be target of a local begin. - -localBeginPossible :: CM n s -> Node n -> Bool -localBeginPossible cm n = - if ntype n `elem` okNodes - && (not . any (==0) $ nparents n) -- nodes reachable from "root" (that is: node 1) have handled specially - then True - else False - where - okNodes = [MATP,MATL,MATR,BIF] - - - --- | Checks for each node if it can lead to a local end. - -localEndPossible :: CM n s -> Node n -> Bool -localEndPossible cm n = - if ntype n `elem` okNodes - && (END /= (ntype $ nodes cm ! (nid n +1))) - then True - else False - where - okNodes = [MATP,MATL,MATR,BEGL,BEGR] - - - --- | transform scores into probabilities, given a nullmodel for x - --- TODO quickcheck! - -score2Prob x null - | x == (-1/0) = 0 - | otherwise = exp (x * log 2) * null - - - --- | back into scores - -prob2Score x null - | x == 0 = (-1/0) - | otherwise = log (x / null) / log 2 - - - --- | Transform a state, setting probabilities instead of scores. Requires CM --- knowledge for background model. - --- TODO actually use the background model - -stateScore2Prob :: CM n s -> State s -> State s -stateScore2Prob cm s = error "implement me" - - - --- | Transform a state, setting scores instead of probabilities. - -stateProb2Score :: CM n s -> State s -> State s -stateProb2Score cm s = error "implement me" -transitionTargets :: [Transition] -> [Int] -transitionTargets xs = map f xs where - f (Branch x) = x - f (Transition x _) = x - - -nodeMainStateAssocs :: [(NodeType,StateType)] -nodeMainStateAssocs = - [ (MATP, MP) - , (MATL, ML) - , (MATR, MR) - , (BIF, B) - , (ROOT, S) - , (BEGL, S) - , (BEGR, S) - , (END, E) - ] - --- }}} - - - --- TODO score -> prob : exp (x / log 2) -- CHECK THIS!!! --- score -> prob : exp (x * log 2) --- prob -> score : (log x) / log 2 --- TODO mkGlobal (needed?) --- TODO use logsum : log (exp x + exp y) = x + log (1 + exp (y-x)), where x>y diff --git a/Biobase/Infernal/CM/Export.hs b/Biobase/Infernal/CM/Export.hs deleted file mode 100644 --- a/Biobase/Infernal/CM/Export.hs +++ /dev/null @@ -1,76 +0,0 @@ -{-# LANGUAGE RecordWildCards #-} - --- | write an Infernal model. We aim for version 1.0.2 of Infernal but right --- now, there is no explicit version management! We just return header --- information as we read it. So we are compatible to all versions that can be --- parsed but have no safeguards in place. --- --- NOTE: exported model can apparently be calibrated. --- --- TODO make sure that (Export.toString . Import.fromString == id); the other --- way around is not absolutely necessary (except if the Infernal parser should --- require this at some point) --- --- TODO put functions like stringS into another module? - -module Biobase.Infernal.CM.Export where - -import qualified Data.Array.IArray as A -import Text.Printf - -import Biobase.Infernal.CM - -import Biobase.Infernal.CM.Import -- TODO for testing only! - - - --- | Main export function for CMs. Creates a string that is accepted by --- Infernal - -toString :: CM n s -> String -toString cm@CM{..} = unlines $ hdr ++ ["MODEL:"] ++ nds ++ ["//"] - where - hdr = map (\(k,v) -> k ++ " " ++ v) header - nds = concatMap (nodeToString cm) $ A.elems nodes - - - --- | export a specific node, used by 'toString' - -nodeToString :: CM n s -> Node n -> [String] -nodeToString cm Node{..} = printf "\t\t\t\t[ %s %d ]" (show ntype) nid : map (stateToString cm . (states cm A.!)) nstates - - - --- | export a specific states, used by 'nodeToString' - -stateToString :: CM n s -> State s -> String -stateToString cm State{..} = printf " %2s %6d %6d %2d %6d %6d %-55s %s" (show stype) sid maxP numP minC numC tscores escores - where - maxP - | null sparents = -1 - | otherwise = maximum sparents - numP = length sparents - minC - | null schildren = -1 - | otherwise = minimum $ map tchild schildren - numC - | stype == B = maximum $ map tchild schildren - | otherwise = length schildren - tscores :: String - tscores - | stype == B = "" - | otherwise = concatMap (printf "%s " . stringS . tscore) schildren - escores :: String - escores - | null semission = "" - | otherwise = concatMap (printf "%.3f " . escore) semission - - - --- | export scores for transitions with the Infernal-specific "-infinity" value --- of "*" - -stringS x - | x == (-1/0) = " *" - | otherwise = printf "%.3f" x diff --git a/Biobase/Infernal/CM/Import.hs b/Biobase/Infernal/CM/Import.hs deleted file mode 100644 --- a/Biobase/Infernal/CM/Import.hs +++ /dev/null @@ -1,194 +0,0 @@ - --- --- A parser for Infernal's Covariance Models (CM). Should work with version 1.0 --- of Infernal. --- - --- TODO use generics? --- --- TODO importing the cmcalibrate calibration data is currently broken! (data is parsed as separate lines) - -module Biobase.Infernal.CM.Import (fromFile,fromString) where - -import Data.Maybe (fromJust) -import qualified Data.Array.IArray as A -import qualified Data.List as L -import Text.ParserCombinators.Parsec hiding (State) - -import Biobase.Infernal.CM -import Biobase.RNA hiding (nucE) - - - --- {{{ Parsec - -myChar = alphaNum <|> char '-' <|> char '_' <|> char '.' -headerChar = myChar <|> char '[' <|> char ']' -floating = many1 $ digit <|> char '.' <|> char 'e' <|> char '-' - -headerLine = do - key <- many1 headerChar - --spaces - many1 $ char ' ' - val <- many1 (headerChar <|> char ' ' <|> char '/' <|> char ':') - newline - return (key,val) - -unknownNumbers = do - many1 $ char ' ' - ns <- floating `sepEndBy` (many1 $ char ' ') - newline - return ("",concat ns) - -theHeader = do - hs <- many $ (try headerLine <|> try unknownNumbers) - string "MODEL:" - newline - return hs - --- }}} - --- {{{ Node - -node = do - spaces - string "[ " - name <- many letter - spaces - num <- many digit >>= (return . read) - spaces - string "]" - newline - states <- many $ try (state num) - return $ (Node - { ntype = read name - , nid = num - , nstates = map sid states - , nparents = [] - , nchildren = [] - , ntag = () - } - , states) - --- }}} - --- {{{ stuff - --- TODO put some of this into HsTools/Parsec stuff - -aNum = many1 $ char '-' <|> digit -addneginf :: String -> Double -addneginf "*" = (-1)/0 --"-1000000000.0" -addneginf x = read x - --- }}} - --- {{{ State - -state nodeID = do - spaces - name <- many1 myChar - spaces - sid <- aNum >>= (return . read) - spaces - plast <- aNum >>= (return . read) - spaces - pnum <- aNum >>= (return . read) - spaces - cfirst <- aNum >>= (return . read) - spaces - cnum <- aNum >>= (return . read) - probs <- ( (many $ myChar <|> char '*') `sepBy` (many1 $ char ' ') >>= - return . map addneginf . filter ((/=) "") ) - newline - let s = case name of - "B" -> State - { stype = read name - , sid = sid - , snode = nodeID - , sparents = [plast-pnum+1 .. plast] - , schildren = [Branch cfirst, Branch cnum] - , semission = [] - , stag = () - } - _ -> State - { stype = read name - , sid = sid - , snode = nodeID - , sparents = [plast-pnum+1 .. plast] - , schildren = zipWith Transition [cfirst .. cfirst+cnum-1] probs - , semission = let keep = drop cnum probs in - case length keep of - 0 -> [] - 4 -> zipWith EmitS acgu keep - 16 -> zipWith (\(k1,k2) v -> EmitP k1 k2 v) acguPairs keep - _ -> error $ "strange number of probabilities" ++ show (keep) - , stag = () - } - return s - --- }}} - --- {{{ Model - -models = do - ms <- many model - eof - return ms - -model = do - h <- theHeader - ns <- many node - let states = concatMap snd ns - let nodes = map (addPCinfo states . fst) ns - -- just add all the node parent / child info - string "//" - newline - -- eof -- removed, we want to be able to read concatenated models! - return $ CM - { nodes = A.array (0, length nodes -1) $ zip (map nid nodes) nodes - , states = A.array (0, length states -1) $ zip (map sid states) states - , header = h - , localBegin = A.array (0, length states -1) $ zip [0 .. length states -1] (0.0 : repeat (-1/0)) - , localEnd = A.array (0, length states -1) $ zip [0 .. length states -1] (repeat (-1/0)) - , cmType = CMScore - , nullModel = A.array (nucA,nucU) $ zip acgu (map read . words . fromJust $ "NULL" `lookup` h) -- TODO circumvents the whole parsing stuff! - } - --- }}} - --- {{{ Stuff - -addPCinfo states n = - let - s = nstates n - sp = L.nub $ L.sort $ concatMap (sparents . (states !!)) s - sc = L.nub $ L.sort $ concatMap (transitionTargets . schildren . (states !!)) s - np = (L.nub $ L.sort $ map (snode . (states !!)) sp) L.\\ [nid n] - nc = (L.nub $ L.sort $ map (snode . (states !!)) sc) L.\\ [nid n] - in - n {nparents = np, nchildren = nc} - - - --- | Two types of parsing, once using a file and once by parsing a string. - -fromFile f = parseFromFile models f -fromString s = parse models "(stdin)" s - --- }}} - --- | Helper function to remove impossible state transitions (those that have --- -infty score). - --- TODO move to InfernalCM.hs and have it for Prob and Score both - -{- -canonize cm = cm {states = A.amap f $ states cm} where - f s = s {schildren = filter g $ schildren s} - g (Branch _) = True - g (Transition _ v) = - h (_,s) - | s == (-1)/0 = False - | otherwise = True --} diff --git a/Biobase/Infernal/Stockholm.hs b/Biobase/Infernal/Stockholm.hs deleted file mode 100644 --- a/Biobase/Infernal/Stockholm.hs +++ /dev/null @@ -1,16 +0,0 @@ - --- | The beginnings of a Stockholm File parser. - -module Biobase.Infernal.Stockholm where - --- | Stockholm format data. We have a set of sequences with sequence data, a --- set of column annotations and unknown data (actually: known, but we do not --- care). - -data Stockholm = Stockholm - { sequences :: [(String,String)] - , colAnnotations :: [(String,String)] -- #=GC tag - , exAnnotations :: [(String,String)] -- NOTE this is not in the Stockholm format! - , unknown :: [String] - } - deriving (Show) diff --git a/Biobase/Infernal/Stockholm/Import.hs b/Biobase/Infernal/Stockholm/Import.hs deleted file mode 100644 --- a/Biobase/Infernal/Stockholm/Import.hs +++ /dev/null @@ -1,35 +0,0 @@ - --- | A simple parser for Stockholm data. As we do not interpret most stuff --- right now, Parsec is not required. - -module Biobase.Infernal.Stockholm.Import where - -import Biobase.Infernal.Stockholm -import Data.List - -fromFile :: String -> IO Stockholm -fromFile fname = do - ls <- readFile fname >>= return . filter (not . null) . lines - let (colanno, rest1) = partition (\x -> length x >= 4 && (and $ zipWith (==) "#=GC" x)) ls - let (seqdata, rest2) = partition (\(x:_) -> x/='#' && x/='/') rest1 - let (exdata , unk) = partition (\x -> length x >= 4 && (and $ zipWith (==) "#=EX" x)) rest2 - {- - print colanno - print rest1 - print seqdata - print rest2 - print exdata - print unk - -} - return Stockholm - { sequences = map (mkPair . words) seqdata - , colAnnotations = map (mkPair . drop 1 . words) colanno - , exAnnotations = map (mkPair . drop 1 . words) exdata - , unknown = unk - } - - - -mkPair [] = error "empty list" -mkPair [x] = error $ "just one element: " ++ show x -mkPair (x:xs) = (x,unwords xs) diff --git a/Biobase/RNA.hs b/Biobase/RNA.hs --- a/Biobase/RNA.hs +++ b/Biobase/RNA.hs @@ -1,247 +1,168 @@ +{-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE GeneralizedNewtypeDeriving #-} +{-# LANGUAGE MultiParamTypeClasses #-} {-# LANGUAGE StandaloneDeriving #-} {-# LANGUAGE TypeSynonymInstances #-} -{-# LANGUAGE PatternGuards #-} -{-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE FlexibleInstances #-} -- | Encodings for nucleotides and pairs. Currently very ViennaRNA-centric. -- This might change over time. -- -- TODO do not export Nucleotide ctor? -- --- TODO consider adding "nucAmpersand" as a chain seperator. This could make --- some things easier. +-- TODO With GHC 7.2 we should be able to activate Read,Show for the overlapped +-- VU.Vector Nucleotide and get better string representations for 'Primary'. +-- +-- TODO nucE -> nucN ? +-- +-- TODO add all possible characters? module Biobase.RNA where +import Control.DeepSeq import Data.Char (toUpper) import Data.Ix -import qualified Data.Vector.Unboxed as VU +import Data.Primitive.Types +import qualified Data.ByteString.Char8 as BSS +import qualified Data.ByteString.Lazy.Char8 as BSL import qualified Data.Vector.Generic as VG import qualified Data.Vector.Generic.Mutable as VGM -import Data.Primitive.Types +import qualified Data.Vector.Unboxed as VU -import Data.PrimitiveArray -import Data.PrimitiveArray.Ix +import Biobase.RNA.NucBounds -- * single nucleotide encoding +-- | Represent nucleotides as machine 'Int's + newtype Nucleotide = Nucleotide Int - deriving (Eq,Ord,Ix) + deriving (Eq,Ord,Ix,NFData) -(nucE:nucA:nucC:nucG:nucU:_) = map Nucleotide [0..] +-- | We have the canonical nucleotides A,C,G,U, unknown nucleotides as E and a +-- chain symbol used to separate different nucleotide chains. -nucCharList :: [(Int,Char)] -nucCharList = zip [0..] "EACGU" -charNucList :: [(Char,Int)] -charNucList = zip "EACGU" [0..] +(nucE: nucA:nucC:nucG:nucU: nucChain:_) = map Nucleotide [0..] -deriving instance Prim Nucleotide -deriving instance VGM.MVector VU.MVector Nucleotide -deriving instance VG.Vector VU.Vector Nucleotide -deriving instance VU.Unbox Nucleotide +-- | Conversion to and from nucleotides. --- TODO maybe use 'assert' and do not catch errors in optimized operation?! +class MkNucleotide a where + mkNuc :: a -> Nucleotide + fromNuc :: Nucleotide -> a +-- | instance for characters + +instance MkNucleotide Char where + mkNuc c' + | Just n <- c `lookup` c2n = n + | c=='&' = nucChain + | otherwise = nucE + where + c2n = zip "ACGU" [nucA..nucU] + c = toUpper c' + fromNuc n + | Just c <- n `lookup` n2c = c + | n == nucChain = '&' + | otherwise = 'E' + where n2c = zip [nucA..nucU] "ACGU" + +-- | The enum instance. + instance Enum Nucleotide where toEnum x - | x>=0 && x<=4 = Nucleotide x + | x>=0 && x<=5 = Nucleotide x | otherwise = error $ "can't make to enum" ++ show x fromEnum (Nucleotide x) = x {-# INLINE toEnum #-} {-# INLINE fromEnum #-} +deriving instance Prim Nucleotide +deriving instance VGM.MVector VU.MVector Nucleotide +deriving instance VG.Vector VU.Vector Nucleotide +deriving instance VU.Unbox Nucleotide + +-- | The Bounded instance uses bounds that are explicitly "too small" as other +-- symbols have special meaning and should not be used with normal bounds. +-- +-- TODO think about better ways to handle this! + instance Bounded Nucleotide where minBound = nucE - maxBound = nucU + maxBound = nucChain +-- | The typical lookup table should use 'NucBounds' instead of 'Bounded' and +-- here 'minExtended' and 'maxExtended' as this covers all characters one would +-- typically expect. + +instance NucBounds Nucleotide where + minNormal = nucA + maxNormal = nucU + minExtended = nucE + maxExtended = nucU + instance Show Nucleotide where - show (Nucleotide x) - | Just s <- x `lookup` nucCharList = [s] + show n = [fromNuc n] instance Read Nucleotide where readsPrec p [] = [] readsPrec p (x:xs) | x ==' ' = readsPrec p xs - | Just n <- x `lookup` charNucList = [(Nucleotide n,xs)] + | Just n <- x `lookup` c2n = [(n,xs)] | otherwise = [] - -charToNucleotide :: Char -> Nucleotide -charToNucleotide = f . toUpper where - f x - | Just n <- x `lookup` charNucList = Nucleotide n - | otherwise = nucE - --- | Since this function is used rather often. - -char2nuc = charToNucleotide - -nucleotideToChar :: Nucleotide -> Char -nucleotideToChar (Nucleotide n) = f n where - f x - | Just c <- x `lookup` nucCharList = c + where + c2n = zip "EACGU&" [nucE .. nucChain] --- * Represent nucleotide sequences using a PrimArray. --- --- TODO switch back to Unboxed.Vector? Test speed of both approaches! - -type Primary = PrimArray Int Nucleotide - --- | --- --- TODO is this fast enough? +-- * Representing nucleotide chains, aka primary sequence -instance Eq Primary where - xs == ys = bounds xs == bounds ys && toList xs == toList ys +-- | We use plain unboxed vectors for the primary sequence -instance Show Primary where - show p = "mkPrimary " ++ unPrimary p +type Primary = VU.Vector Nucleotide -instance Read Primary where - readsPrec p xs - | "mkPrimary " == chk = [(mkPrimary ps,ys)] - | otherwise = error $ show (chk,others,ps,ys) - where - (chk,others) = splitAt 10 xs - (ps,ys) = span (`elem` "ACGUE") others +-- | Conversion from/to primary sequence class MkPrimary a where mkPrimary :: a -> Primary unPrimary :: Primary -> a -instance MkPrimary String where - mkPrimary xs = fromList 0 (length xs -1) $ map charToNucleotide xs - unPrimary v = map nucleotideToChar $ toList v - {-# INLINE mkPrimary #-} - {-# INLINE unPrimary #-} - -instance MkPrimary (VU.Vector Char) where - mkPrimary = mkPrimary . VU.toList - unPrimary = VU.fromList . unPrimary - {-# INLINE mkPrimary #-} - {-# INLINE unPrimary #-} - - - --- * Encoding of Watson-Crick and Wobble Pairs in the Vienna RNA package style. - -newtype ViennaPair = ViennaPair Int - deriving (Eq,Ord,Ix) - -(vpNP:vpCG:vpGC:vpGU:vpUG:vpAU:vpUA:vpNS:_) = map ViennaPair [0..] - -vpStrList :: [(Int,String)] -vpStrList = zip [0..] ["NP","CG","GC","GU","UG","AU","UA","NS"] -strVpList :: [(String,Int)] -strVpList = zip ["NP","CG","GC","GU","UG","AU","UA","NS"] [0..] - -deriving instance VGM.MVector VU.MVector ViennaPair -deriving instance VG.Vector VU.Vector ViennaPair -deriving instance VU.Unbox ViennaPair -deriving instance Prim ViennaPair - -instance Enum ViennaPair where - toEnum x - | x>=0 && x<=7 = ViennaPair x - | otherwise = error $ "can't make to enum" ++ show x - fromEnum (ViennaPair x) = x - {-# INLINE toEnum #-} - {-# INLINE fromEnum #-} +instance MkPrimary [Nucleotide] where + mkPrimary = VU.fromList + unPrimary = VU.toList -instance Bounded ViennaPair where - minBound = vpNP - maxBound = vpNS +instance MkPrimary String where + mkPrimary = VU.fromList . map mkNuc + unPrimary = map fromNuc . VU.toList -instance Show ViennaPair where - show (ViennaPair x) - | Just s <- x `lookup` vpStrList = s +instance MkPrimary BSL.ByteString where + mkPrimary = VU.fromList . map mkNuc . BSL.unpack + unPrimary = BSL.pack . map fromNuc . VU.toList -instance Read ViennaPair where - readsPrec p [] = [] - readsPrec p [x] = [] - readsPrec p (x:y:xs) - | x ==' ' = readsPrec p (y:xs) - | Just n <- (x:y:[]) `lookup` strVpList = [(ViennaPair n,xs)] - | otherwise = [] +instance MkPrimary BSS.ByteString where + mkPrimary = VU.fromList . map mkNuc . BSS.unpack + unPrimary = BSS.pack . map fromNuc . VU.toList --- * convenience lists - -acgu = [nucA .. nucU] -acguStr = map show acgu -acguPairs = [(b1,b2) | b1<-acgu, b2<-acgu] -eacgu = [nucE .. nucU] -npnsP = [vpNP .. vpNS] -cguaP = [vpCG .. vpUA] -cgnsP = [vpCG .. vpNS] +-- * Some convenience structures +eacgu = [nucE..nucU] +acgu = [nucA..nucU] +acguPairs = [(x,y)|x<-acgu,y<-acgu] --- * Transformations between Nucleotide and ViennaPair +-- * Currently impossible instances (changing with ghc 7.2) -pairTable = fromAssocs (nucE,nucE) (nucU,nucU) vpNP - [ ((nucC,nucG),vpCG) - , ((nucG,nucC),vpGC) - , ((nucG,nucU),vpGU) - , ((nucU,nucG),vpUG) - , ((nucA,nucU),vpAU) - , ((nucU,nucA),vpUA) - ] -{- NOINLINE pairTable #-} +{- TODO reactivate these instances with a later version of GHC +instance Show Primary where + show p = "mkPrimary " ++ unPrimary p -toPair :: Nucleotide -> Nucleotide -> ViennaPair -toPair x y = pairTable ! (x,y) -{- - | x==nucC&&y==nucG = vpCG - | x==nucG&&y==nucC = vpGC - | x==nucG&&y==nucU = vpGU - | x==nucU&&y==nucG = vpUG - | x==nucA&&y==nucU = vpAU - | x==nucU&&y==nucA = vpUA - | otherwise = vpNP +instance Read Primary where + readsPrec p xs + | "mkPrimary " == chk = [(mkPrimary ps,ys)] + | otherwise = error $ show (chk,others,ps,ys) + where + (chk,others) = splitAt 10 xs + (ps,ys) = span (`elem` "ACGUE") others -} -{-# INLINE toPair #-} - -pairs :: Nucleotide -> Nucleotide -> Bool -pairs a b = (>vpNP) $ toPair a b -{-# INLINE pairs #-} - -vienna2Tuple p - | p==vpCG = (nucC,nucG) - | p==vpGC = (nucG,nucC) - | p==vpGU = (nucG,nucU) - | p==vpUG = (nucU,nucG) - | p==vpAU = (nucA,nucU) - | p==vpUA = (nucU,nucA) - | otherwise = error "non-standard pairs can't be backcasted" -{-# INLINE vienna2Tuple #-} - -vienna2RevTuple p - | p==vpCG = (nucG,nucC) - | p==vpGC = (nucC,nucG) - | p==vpGU = (nucU,nucG) - | p==vpUG = (nucG,nucU) - | p==vpAU = (nucU,nucA) - | p==vpUA = (nucA,nucU) - | otherwise = error "non-standard pairs can't be backcasted" -{-# INLINE vienna2RevTuple #-} - -tuple2Vienna (b1,b2) - | b1==nucC&&b2==nucG = vpCG - | b1==nucG&&b2==nucC = vpGC - | b1==nucG&&b2==nucU = vpGU - | b1==nucU&&b2==nucG = vpUG - | b1==nucA&&b2==nucU = vpAU - | b1==nucU&&b2==nucA = vpUA - | otherwise = vpNS -{-# INLINE tuple2Vienna #-} - -tuple2RevVienna (a,b) = tuple2Vienna (b,a) -{-# INLINE tuple2RevVienna #-} diff --git a/Biobase/RNA/Complex.hs b/Biobase/RNA/Complex.hs new file mode 100644 --- /dev/null +++ b/Biobase/RNA/Complex.hs @@ -0,0 +1,72 @@ +{-# LANGUAGE RecordWildCards #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE TypeSynonymInstances #-} +{-# LANGUAGE MultiParamTypeClasses #-} + +-- | Internal representation of RNA complexes. Ideally, each DataSource should +-- define an instance of 'MkComplex' for easy importing of foreign data. +-- +-- TODO unify with Biobase.Structure.Complex + +module Biobase.RNA.Complex where + +import Control.DeepSeq +import qualified Data.Vector.Unboxed as VU + +import Biobase.RNA +import Biobase.RNA.Pairs + + + +-- | An RNA complex. One or more RNA molecules linearly encoding as a +-- 'sequence'. Structure is stored as an unboxed vector. The (unspecified) type +-- is most likely either a 'CwwPair' or an 'ExtPair'. The 'chain' encodes (for +-- more than one participating nucleotide sequence) to which chain each +-- nucleotide in the sequence belongs. +-- +-- NOTE 'chain' may be empty and there may be "ampersand" markers in the +-- 'sequence' which then define split positions. + +data (VU.Unbox a, Show a, Read a) => Complex a = Complex + { sequence :: Primary + , structure :: Pairs a + , chain :: VU.Vector Char + } + deriving (Read,Show) + +-- | Create 'Complex'es from data sources. + +class MkComplex a b where + mkComplex :: a -> Complex b + + + +-- * 'CwwPair' instances + +-- | Create a complex where only sequence is given. + +instance MkComplex String CwwPair where + mkComplex s = error "write me" -- Complex (mkPrimary s) (mkSecondary "") (VU.empty) + +-- | Default interactive instance. Given sequence and dot-bracket string, both +-- as strings, create an instance. + +instance MkComplex (String,String) CwwPair where + mkComplex (s,p) = error "write me" -- Complex (mkPrimary s) (mkSecondary p) (VU.empty) + + + +-- * 'ExtPair' instances + +instance MkComplex String ExtPair where + mkComplex s = error "write me" -- Complex (mkPrimary s) (mkExtended "") (VU.empty) + + + +-- * Instances for Complex + +-- | The NFData instance. Since these are unboxed vectors, we just need to make +-- sure that evaluation begins. + +instance (NFData a, Read a, Show a, VU.Unbox a) => NFData (Complex a) where + rnf Complex{..} = rnf (VU.length sequence) `seq` rnf structure `seq` rnf (VU.length chain) `seq` () diff --git a/Biobase/RNA/Hashes.hs b/Biobase/RNA/Hashes.hs new file mode 100644 --- /dev/null +++ b/Biobase/RNA/Hashes.hs @@ -0,0 +1,32 @@ +{-# LANGUAGE GeneralizedNewtypeDeriving #-} + +-- | Fast hash functions for 'Primary' sequences. + +module Biobase.RNA.Hashes where + +import Control.DeepSeq +import Control.Exception.Base (assert) +import Data.Ix +import qualified Data.Vector.Unboxed as VU + +import Biobase.RNA + +newtype HashedPrimary = HashedPrimary Int + deriving (Eq,Ord,Ix,NFData,Read,Show,Enum,Bounded) + +-- | Given a piece of primary sequence information, reduce it to an index. +-- +-- Will throw an assertion in debug code if 'ps' are not within bounds. Note +-- that "mkPrimary [minBound]" and "mkPrimary [minBound,minBound]" map to the +-- same index. Meaning that indices are only unique within the same length +-- group. Furthermore, indices with different (l,u)-bounds are not compatible +-- with each other. All indices start at 0. +-- +-- The empty input produces an index of 0. + +mkHashedPrimary :: (Nucleotide,Nucleotide) -> Primary -> HashedPrimary +mkHashedPrimary (l,u) ps = assert (VU.all (\p -> l<=p && p<=u) ps) $ HashedPrimary idx where + idx = VU.sum $ VU.zipWith f ps (VU.enumFromStepN (VU.length ps -1) (-1) (VU.length ps)) + f p c = (fromEnum p - fromEnum l) * (cnst^c) + cnst = fromEnum u - fromEnum l + 1 +{-# INLINE mkHashedPrimary #-} diff --git a/Biobase/RNA/NucBounds.hs b/Biobase/RNA/NucBounds.hs new file mode 100644 --- /dev/null +++ b/Biobase/RNA/NucBounds.hs @@ -0,0 +1,54 @@ + +-- | A special class of bounds for RNA/pair encodings that are used to index +-- into tables. We typically encode more in the alphabets than we want to use +-- to index, so in order to keep things simple, we have specialized bounds. + +module Biobase.RNA.NucBounds where + +import GHC.Tuple + + + +-- | 'minNormal' and 'maxNormal' encode for, say, ACGU; while 'minExtended' and +-- 'maxExtended' would allow 'N' as well. See Biobase.RNA and +-- Biobase.RNA.ViennaPair for instances. + +class Bounded a => NucBounds a where + minNormal :: a + maxNormal :: a + minExtended :: a + maxExtended :: a + + + +-- * Instances for tuples of size 2-6 + +instance (NucBounds a, NucBounds b) => NucBounds (a,b) where + minNormal = (minNormal, minNormal) + maxNormal = (maxNormal, maxNormal) + minExtended = (minExtended, minExtended) + maxExtended = (maxExtended, maxExtended) + +instance (NucBounds a, NucBounds b, NucBounds c) => NucBounds (a,b,c) where + minNormal = (minNormal, minNormal, minNormal) + maxNormal = (maxNormal, maxNormal, maxNormal) + minExtended = (minExtended, minExtended, minExtended) + maxExtended = (maxExtended, maxExtended, maxExtended) + +instance (NucBounds a, NucBounds b, NucBounds c, NucBounds d) => NucBounds (a,b,c,d) where + minNormal = (minNormal, minNormal, minNormal, minNormal) + maxNormal = (maxNormal, maxNormal, maxNormal, maxNormal) + minExtended = (minExtended, minExtended, minExtended, minExtended) + maxExtended = (maxExtended, maxExtended, maxExtended, maxExtended) + +instance (NucBounds a, NucBounds b, NucBounds c, NucBounds d, NucBounds e) => NucBounds (a,b,c,d,e) where + minNormal = (minNormal, minNormal, minNormal, minNormal, minNormal) + maxNormal = (maxNormal, maxNormal, maxNormal, maxNormal, maxNormal) + minExtended = (minExtended, minExtended, minExtended, minExtended, minExtended) + maxExtended = (maxExtended, maxExtended, maxExtended, maxExtended, maxExtended) + +instance (NucBounds a, NucBounds b, NucBounds c, NucBounds d, NucBounds e, NucBounds f) => NucBounds (a,b,c,d,e,f) where + minNormal = (minNormal, minNormal, minNormal, minNormal, minNormal, minNormal) + maxNormal = (maxNormal, maxNormal, maxNormal, maxNormal, maxNormal, maxNormal) + minExtended = (minExtended, minExtended, minExtended, minExtended, minExtended, minExtended) + maxExtended = (maxExtended, maxExtended, maxExtended, maxExtended, maxExtended, maxExtended) diff --git a/Biobase/RNA/Pairs.hs b/Biobase/RNA/Pairs.hs new file mode 100644 --- /dev/null +++ b/Biobase/RNA/Pairs.hs @@ -0,0 +1,124 @@ +{-# LANGUAGE ParallelListComp #-} +{-# LANGUAGE RecordWildCards #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE TypeSynonymInstances #-} + +-- | We want to be able to generate secondary structure representations from a +-- wide variety of inputs. +-- +-- TODO include DotBracket, CT, lists of pairs,... some lossy conversions from +-- FamPair to CwwPair, etc... + +module Biobase.RNA.Pairs where + +import Control.DeepSeq +import Data.Vector.Unboxed.Read +import qualified Data.Vector.Unboxed as VU +import Data.List (nub) + +-- | Canonical secondary structures, with the Watson-Crick cis pair family. + +class MkSecondary a where + mkSecondary :: a -> Pairs CwwPair + +-- | Extended secondary structures, with any kind of pair family for each pair. + +class MkExtended a where + mkExtended :: a -> Pairs ExtPair + +-- | Pair container. We store the sequence length as well in order to be able +-- to reconstruct the unpaired nucleotides without having to store them +-- explicitly. + +data (VU.Unbox a, Show a, Read a) => Pairs a = Pairs + { seqLength :: Int + , pairs :: VU.Vector a + } + deriving (Eq,Show,Read) + + + +-- * Canonical secondary structures + +-- | Most specific instance, just wrapping length and vector + +instance MkSecondary (Int,VU.Vector CwwPair) where + mkSecondary (l,xs) = Pairs l xs + +-- | Only conversion to vector from list + +instance MkSecondary (Int,[CwwPair]) where + mkSecondary (l,xs) = mkSecondary (l,VU.fromList xs) + +-- | Instance for MkSecondary + +instance MkSecondary String where + mkSecondary xs = mkPairList xs + + + +-- * Extended secondary structures + +-- | The most specific instance possible + +instance MkExtended (Int,VU.Vector ExtPair) where + mkExtended (l,xs) = Pairs l xs + +-- | Just turn list into a vector + +instance MkExtended (Int,[ExtPair]) where + mkExtended (l,xs) = mkExtended (l,VU.fromList xs) + +-- | Create from a string. Here, we do not have pair family information and so +-- assume that all pairs are of the cWW type. + +instance MkExtended String where + mkExtended xs = mkExtended $ mkSecondary xs + +-- | Create extended structure from CwwPair structure. +-- +-- TODO still undefined until we have isostericity included + +instance MkExtended (Pairs CwwPair) where + mkExtended Pairs{..} = Pairs seqLength (VU.map (\(i,j) -> (i,j,error "mkExtended needs to know about isostericity")) pairs) + + + +-- * A basic dot-bracket parser. Not a datasource due to it being universally needed. + +-- | We only return 'CwwPair's as dot-bracket pairs, others can not be encoded. + +mkPairList :: String -> Pairs CwwPair +mkPairList s = mkSecondary (length s, concatMap (\b -> mkPairListWith b s) bs) where + z = nub s + bs = filter (\[l,r] -> l `elem` z && r `elem` z) defBrackets + +-- | Generic dot-bracket parser accepting more than one bracket type. + +mkPairListWith :: String -> String -> [CwwPair] +mkPairListWith [l,r] s = f [] 0 s where + f _ _ [] = [] + f s k (x:xs) + | x==l = f (k:s) (k+1) xs + | x==r = (head s, k) : f (tail s) (k+1) xs + | otherwise = f s (k+1) xs + +-- | These are the default bracket types to consider + +defBrackets = + [ "()" + , "{}" + , "[]" + , "<>" + ] ++ + [ [b,s] | b <- ['A'..'Z'] | s <- ['a'..'z'] ] + + + +-- * Types and instances + +type CwwPair = (Int,Int) +type ExtPair = (Int,Int,Int) + +instance (NFData a, Read a, Show a, VU.Unbox a) => NFData (Pairs a) where + rnf Pairs{..} = rnf seqLength `seq` rnf (VU.length pairs) `seq` () diff --git a/Biobase/RNA/ViennaPair.hs b/Biobase/RNA/ViennaPair.hs new file mode 100644 --- /dev/null +++ b/Biobase/RNA/ViennaPair.hs @@ -0,0 +1,108 @@ +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE GeneralizedNewtypeDeriving #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE StandaloneDeriving #-} + +-- | Encoding of Watson-Crick and Wobble Pairs in the Vienna RNA package style. + +module Biobase.RNA.ViennaPair where + +import Data.Ix +import qualified Data.Vector.Generic as VG +import Data.Primitive.Types +import qualified Data.Vector.Generic.Mutable as VGM +import qualified Data.Vector.Unboxed as VU +import Data.Tuple.HT (swap) + +import Biobase.RNA +import Biobase.RNA.NucBounds + + + +-- | Use machine Ints internally + +newtype ViennaPair = ViennaPair Int + deriving (Eq,Ord,Ix) + +(vpNP:vpCG:vpGC:vpGU:vpUG:vpAU:vpUA:vpNS:_) = map ViennaPair [0..] + +class MkViennaPair a where + mkViennaPair :: a -> ViennaPair + fromViennaPair :: ViennaPair -> a + +instance MkViennaPair (Nucleotide,Nucleotide) where + mkViennaPair (b1,b2) + | b1==nucC&&b2==nucG = vpCG + | b1==nucG&&b2==nucC = vpGC + | b1==nucG&&b2==nucU = vpGU + | b1==nucU&&b2==nucG = vpUG + | b1==nucA&&b2==nucU = vpAU + | b1==nucU&&b2==nucA = vpUA + | otherwise = vpNS + fromViennaPair p + | p==vpCG = (nucC,nucG) + | p==vpGC = (nucG,nucC) + | p==vpGU = (nucG,nucU) + | p==vpUG = (nucU,nucG) + | p==vpAU = (nucA,nucU) + | p==vpUA = (nucU,nucA) + | otherwise = error "non-standard pairs can't be backcasted" + +deriving instance VGM.MVector VU.MVector ViennaPair +deriving instance VG.Vector VU.Vector ViennaPair +deriving instance VU.Unbox ViennaPair +deriving instance Prim ViennaPair + +instance Enum ViennaPair where + toEnum x + | x>=0 && x<=7 = ViennaPair x + | otherwise = error $ "can't make to enum" ++ show x + fromEnum (ViennaPair x) = x + {-# INLINE toEnum #-} + {-# INLINE fromEnum #-} + +instance Bounded ViennaPair where + minBound = vpNP + maxBound = vpNS + +instance NucBounds ViennaPair where + minNormal = vpCG + maxNormal = vpUA + minExtended = vpNP + maxExtended = vpNS + +instance Show ViennaPair where + show x + | Just s <- x `lookup` pairToString = s + +instance Read ViennaPair where + readsPrec p [] = [] + readsPrec p [x] = [] + readsPrec p (x:y:xs) + | x ==' ' = readsPrec p (y:xs) + | Just n <- (x:y:[]) `lookup` s2p = [(n,xs)] + | otherwise = [] + where s2p = (map swap pairToString) + + + +-- | reverse a vienna pair + +revPair :: ViennaPair -> ViennaPair +revPair p + | p==vpCG = vpGC + | p==vpGC = vpCG + | p==vpGU = vpUG + | p==vpUG = vpGU + | p==vpAU = vpUA + | p==vpUA = vpAU + | p==vpNP = vpNP + | p==vpNS = vpNS + + + +-- * Convenience structures + +cguaP = [vpCG..vpUA] +cgnsP = [vpCG..vpNS] +pairToString = [(vpCG,"CG"),(vpGC,"GC"),(vpUA,"UA"),(vpAU,"AU"),(vpGU,"GU"),(vpUG,"UG"),(vpNS,"NS"),(vpNP,"NP")] diff --git a/Biobase/Turner.hs b/Biobase/Turner.hs deleted file mode 100644 --- a/Biobase/Turner.hs +++ /dev/null @@ -1,41 +0,0 @@ -{-# LANGUAGE TypeSynonymInstances #-} - --- | Mathews / Turner nearest neighbor energy tables for interactions between --- RNA-RNA and DNA-DNA. Using lazy tables for easy filling of missing data. --- Transform tables into uvectors for efficient access. - -module Biobase.Turner where - -import Biobase.RNA (Nucleotide) -import qualified Biobase.Turner.Tables as T - - - --- | Data structure containing all necessary tables for the Turner 2004 model. --- Note the temperature, using an additional function, rescaling becomes --- possible. A Show and a Read instance are automatically derived but one --- should use the Turner and/or Vienna import/export facilities. All data uses --- complete tables, where any base can pair with any other. Missing data ('.' --- in tables) is given by 'Nothing', otherwise we have 'Just Int'. --- --- In general, have a look here: --- \url{http://rna.urmc.rochester.edu/NNDB/turner04/index.html} where --- parameters are explained. --- --- Conversion from Turner to Vienna is only possible in one direction as a --- number of parameters are not used in Vienna RNA Folding. - - - --- We are missing tables for DNA-RNA data: helixdr.dat, stackdr.dat, --- stackdr.dh. Several other files have unknown data, too. These are: --- dnadynalign*, fam\_hmm\_pars.dat, helix.dat, int22-exp.dh - -type TurnerTables = T.Turner2004 Basepair Nucleotide Double - -type TurnerEntropy = TurnerTables -type TurnerEnthalpy = TurnerTables -type Temperature = Double -type TurnerSet = (Temperature,TurnerEntropy,TurnerEnthalpy) - -type Basepair = (Nucleotide,Nucleotide) diff --git a/Biobase/Turner/Import.hs b/Biobase/Turner/Import.hs deleted file mode 100644 --- a/Biobase/Turner/Import.hs +++ /dev/null @@ -1,176 +0,0 @@ -{-# LANGUAGE NoMonomorphismRestriction #-} - --- | Turner file parser. Returns a Turner2004 data structure. We store data in --- the same way it is stored in the ViennaRNA package. Pairs are tuples --- however. --- --- How is 'stack' data stored: --- --- AX --- UY -> ((A,U),(Y,X)) --- --- How 'iloop1x1' is stored: --- --- X --- A G --- U C -> ((A,U),(C,G),X,Y) --- Y --- --- Now 'iloop1x2' is stored: --- --- X --- A G --- U C -> ((A,U),(C,G),X,C,Y), single (X) first, then 5' to 3' --- YC --- --- 'iloop2x2' is stored: --- --- XY --- A G --- U C -> ((A,U),(C,G),X,Y,y,x), X-->Y then x<--y --- xy --- ---TODO not sure if dangle3/dangle5 are correctly split or if they should switch - -module Biobase.Turner.Import - ( parseTurner - ) where - -import Control.Applicative -import Control.Arrow (first) -import Control.Monad -import Data.Either -import Data.List.Split (splitEvery) -import Data.List (transpose) -import qualified Data.Map as M -import System.FilePath.Find -import Text.Parsec hiding ((<|>), many) -import Text.Parsec.String - -import Biobase.RNA -import Data.Ix.Tuple -import Data.PrimitiveArray -import Data.PrimitiveArray.Ix -import Text.Parsec.Numbers - -import Biobase.Turner.Tables - - - --- * Associate every file with a parse. - -minp = (nucA,nucA) -maxp = (nucU,nucU) -minpb = (minp,nucA) -maxpb = (maxp,nucU) -minpp = (minp,minp) -maxpp = (maxp,maxp) -minppbb = (minp,minp,nucA,nucA) -maxppbb = (maxp,maxp,nucU,nucU) -minppbbb = (minp,minp,nucA,nucA,nucA) -maxppbbb = (maxp,maxp,nucU,nucU,nucU) -minppbbbb = (minp,minp,nucA,nucA,nucA,nucA) -maxppbbbb = (maxp,maxp,nucU,nucU,nucU,nucU) -minpbb = (minp,nucA,nucA) -maxpbb = (maxp,nucU,nucU) --- (4,3) switched for vienna rna compatibility -keyspp = [((k1,k2),(k4,k3)) | k1 <- acgu, k3 <- acgu, k2 <- acgu, k4 <- acgu] -keyspb = [((k1,k2),k3) | k1 <- acgu, k2 <- acgu, k3 <- acgu] -keyspbb = [((k1,k2),k3,k4) | k1 <- acgu, k3 <- acgu, k2 <- acgu, k4 <- acgu] -plist11 = [(nucA,nucU),(nucC,nucG),(nucG,nucC),(nucU,nucA),(nucG,nucU),(nucU,nucG)] -plist22 = [(nucA,nucU),(nucC,nucG),(nucG,nucC),(nucG,nucU),(nucU,nucA),(nucU,nucG)] --- (4,3) switched for vienna rna compatibility -keysppbb = [((k1,k2),(k4,k3),k5,k6) | (k1,k2) <- plist11, k5 <- acgu, (k3,k4) <- plist11, k6 <- acgu] --- (4,3) switched for vienna rna compatibility -keysppbbb = [((k1,k2),(k4,k3),k5,k6,k7) | (k1,k2) <- plist11, k6 <- acgu, k5 <- acgu, (k3,k4) <- plist11, k7 <- acgu] --- (4,3) switched for vienna rna compatibility, 5786 is 5'3' order top, bottom! -keysppbbbb = [((k1,k2),(k4,k3),k5,k6,k7,k8) | (k1,k2) <- plist22, (k3,k4) <- plist22, k5 <- acgu, k8 <- acgu, k6 <- acgu, k7 <- acgu] -z = 999999 - - - --- | Given the base dir and a suffix (.dat/.dh most likely), parse the relevant --- files. The prefix can be used for, eg., dna file - -parseTurner prefix basedir fsuffix = do - vstack <- grabB basedir fsuffix $ prefix ++ "stack" - vdangles <- grabB basedir fsuffix $ prefix ++ "dangle" - vlengths <- grabB basedir fsuffix $ prefix ++ "loop" - vhairmm <- grabB basedir fsuffix $ prefix ++ "tstackh" - viloopmm <- grabB basedir fsuffix $ prefix ++ "tstacki" - viloop23mm <- grabB basedir fsuffix $ prefix ++ "tstacki23" - viloop1nmm <- grabB basedir fsuffix $ prefix ++ "tstacki1n" - vmultimm <- grabB basedir fsuffix $ prefix ++ "tstackm" - vextmm <- grabB basedir fsuffix $ prefix ++ "tstack" - viloop11 <- grabB basedir fsuffix $ prefix ++ "int11" - viloop12 <- grabB basedir fsuffix $ prefix ++ "int21" - viloop22 <- grabB basedir fsuffix $ prefix ++ "int22" - vlookups <- fmap (map (first (map charToNucleotide)) . concat) $ mapM (\f -> grabA basedir fsuffix $ prefix ++ f) ["triloop","tloop","hexaloop"] - vmisc <- grabB basedir fsuffix $ prefix ++ "miscloop" - let (vdangle3,vdangle5) = splitAt (length vdangles `div` 2) vdangles - let (_:viloopl:vbulgel:vhairpinl:[]) = transpose $ splitEvery 4 vlengths - return $ Turner2004 - { stack = fromAssocs minpp maxpp z $ zip keyspp vstack - , dangle3 = fromAssocs minpb maxpb z $ zip keyspb vdangle3 - , dangle5 = fromAssocs minpb maxpb z $ zip keyspb vdangle5 - , hairpinL = fromAssocs 0 30 z $ zip [1..] vhairpinl - , bulgeL = fromAssocs 0 30 z $ zip [1..] vbulgel - , iloopL = fromAssocs 0 30 z $ zip [1..] viloopl - , hairpinMM = fromAssocs minpbb maxpbb z $ zip keyspbb vhairmm - , iloopMM = fromAssocs minpbb maxpbb z $ zip keyspbb viloopmm - , iloop2x3MM = fromAssocs minpbb maxpbb z $ zip keyspbb viloop23mm - , iloop1xnMM = fromAssocs minpbb maxpbb z $ zip keyspbb viloop1nmm - , multiMM = fromAssocs minpbb maxpbb z $ zip keyspbb vmultimm - , extMM = fromAssocs minpbb maxpbb z $ zip keyspbb vextmm - , hairpinLookup = M.fromList vlookups - , iloop1x1 = fromAssocs minppbb maxppbb z $ zip keysppbb viloop11 - , iloop1x2 = fromAssocs minppbbb maxppbbb z $ zip keysppbbb viloop12 - , iloop2x2 = fromAssocs minppbbbb maxppbbbb z $ zip keysppbbbb viloop22 - , ninio = vmisc !! 2 - , maxNinio = vmisc !! 1 - , multiOffset = vmisc !! 6 - , multiNuc = vmisc !! 7 - , multiHelix = vmisc !! 8 - , largeLoop = vmisc !! 0 - , termAU = vmisc !! 14 - , intermolecularInit = vmisc !! 19 - } - -grabB basedir fsuffix fname = do - p <- parseFromFile pBlockFile (basedir ++ "/" ++ fname ++ fsuffix) - case p of - (Right ans) -> return ans - (Left err) -> error $ show err - -grabA basedir fsuffix fname = do - p <- parseFromFile pAssocFile (basedir ++ "/" ++ fname ++ fsuffix) - case p of - (Right ans) -> return ans - (Left err) -> error $ show err - - - --- * File parser: We do not want to understand the format, just extract all --- data. - --- | Blocks of data. - -pBlockFile = concat <$> many (try numberRow <|> infoRow) <* eof where - -- a line with a least one number - numberRow = many (char ' ') *> (number <* notFollowedBy (char '\'')) `sepEndBy1` (many $ char ' ') <* newline - -- a number is a dot or some extended floating point number - number = (try dotFloat) <|> (999999 <$ char '.') <|> parseExtFloat "floating point number" - --- | A File with associations. - -pAssocFile = concat <$> many1 (try assocRow <|> infoRow) where - -- one assocs of string, value - assocRow = (\a b -> [(a,b)]) <$ many (char ' ') <*> many1 (oneOf nuc) <* many (char ' ') <*> parseExtFloat <* many (char ' ') - --- a line with crap -infoRow = [] <$ anyChar `manyTill` newline - --- a rather stupid float -dotFloat = (\s a b -> read $ s:'0':a:b) <$> (char '-' <|> pure ' ') <*> char '.' <*> many1 digit - -nuc = "ACGUT" diff --git a/Biobase/Turner/Tables.hs b/Biobase/Turner/Tables.hs deleted file mode 100644 --- a/Biobase/Turner/Tables.hs +++ /dev/null @@ -1,133 +0,0 @@ -{-# LANGUAGE TemplateHaskell #-} -{-# LANGUAGE DeriveFunctor #-} -{-# LANGUAGE RecordWildCards #-} -{-# LANGUAGE TypeSynonymInstances #-} - --- | Turner Nearest Neighbor Energy Tables} - --- TODO (Read) instance? (Vector has no Read instance...) - -module Biobase.Turner.Tables where - -import qualified Data.Map as M - -import Data.PrimitiveArray -import Data.PrimitiveArray.Ix - - - --- | A data structure describing all fields as used by the Turner 2004 --- parameter set. Some fields are commented out until they are being used. - -data Turner2004 a b c = Turner2004 - { stack :: PrimArray (a,a) c - , dangle3 :: PrimArray (a,b) c - , dangle5 :: PrimArray (a,b) c - , hairpinL :: PrimArray Int c - , hairpinMM :: PrimArray (a,b,b) c - , hairpinLookup :: M.Map [b] c --- , hairpinGGG :: c --- , hairpinCslope :: c --- , hairpinCintercept :: c --- , hairpinC3 :: c - , bulgeL :: PrimArray Int c --- , bulgeSingleC :: c - , iloop1x1 :: PrimArray (a,a,b,b) c - , iloop1x2 :: PrimArray (a,a,b,b,b) c - , iloop2x2 :: PrimArray (a,a,b,b,b,b) c - , iloopMM :: PrimArray (a,b,b) c - , iloop2x3MM :: PrimArray (a,b,b) c - , iloop1xnMM :: PrimArray (a,b,b) c - , iloopL :: PrimArray Int c - , multiMM :: PrimArray (a,b,b) c - , ninio :: c - , maxNinio :: c - , multiOffset :: c - , multiNuc :: c - , multiHelix :: c --- , multiAsym :: c --- , multiStrain :: c - , extMM :: PrimArray (a,b,b) c --- , coaxStack :: PrimArray (a,a) c --- , coaxStackOpen :: PrimArray (a,b,b) c --- , coaxStackCont :: PrimArray (a,b,b) c - , largeLoop :: Double - , termAU :: c - , intermolecularInit :: c - } - - - --- | Map functions over the payload. - -dmap f Turner2004{..} = Turner2004 - { stack = amap f stack - , dangle3 = amap f dangle3 - , dangle5 = amap f dangle5 - , hairpinL = amap f hairpinL - , hairpinMM = amap f hairpinMM - , hairpinLookup = fmap f hairpinLookup --- , hairpinGGG = f hairpinGGG --- , hairpinCslope = f hairpinCslope --- , hairpinCintercept = f hairpinCintercept --- , hairpinC3 = f hairpinC3 - , bulgeL = amap f bulgeL --- , bulgeSingleC = f bulgeSingleC - , iloop1x1 = amap f iloop1x1 - , iloop1x2 = amap f iloop1x2 - , iloop2x2 = amap f iloop2x2 - , iloopMM = amap f iloopMM - , iloop2x3MM = amap f iloop2x3MM - , iloop1xnMM = amap f iloop1xnMM - , iloopL = amap f iloopL - , multiMM = amap f multiMM - , ninio = f ninio - , maxNinio = f maxNinio - , multiOffset = f multiOffset - , multiNuc = f multiNuc - , multiHelix = f multiHelix --- , multiAsym = f multiAsym --- , multiStrain = f multiStrain - , extMM = amap f extMM --- , coaxStack = amap f coaxStack --- , coaxStackOpen = amap f coaxStackOpen --- , coaxStackCont = amap f coaxStackCont - , largeLoop = largeLoop - , termAU = f termAU - , intermolecularInit = f intermolecularInit - } - - - --- | Zip two payloads. --- --- TODO right now, we have undefined behaviour when some arrrays are of different length - -dZipWith f t1 t2 = Turner2004 - { stack = zipWithPA f (stack t1) (stack t2) - , dangle3 = zipWithPA f (dangle3 t1) (dangle3 t2) - , dangle5 = zipWithPA f (dangle5 t1) (dangle5 t2) - , hairpinL = zipWithPA f (hairpinL t1) (hairpinL t2) - , hairpinMM = zipWithPA f (hairpinMM t1) (hairpinMM t2) - , hairpinLookup = undefined -- zipWithPA f (hairpinLookup t1) (hairpinLookup t2) - , bulgeL = zipWithPA f (bulgeL t1) (bulgeL t2) - , iloop1x1 = zipWithPA f (iloop1x1 t1) (iloop1x1 t2) - , iloop1x2 = zipWithPA f (iloop1x2 t1) (iloop1x2 t2) - , iloop2x2 = zipWithPA f (iloop2x2 t1) (iloop2x2 t2) - , iloopMM = zipWithPA f (iloopMM t1) (iloopMM t2) - , iloop2x3MM = zipWithPA f (iloop2x3MM t1) (iloop2x3MM t2) - , iloop1xnMM = zipWithPA f (iloop1xnMM t1) (iloop1xnMM t2) - , iloopL = zipWithPA f (iloopL t1) (iloopL t2) - , multiMM = zipWithPA f (multiMM t1) (multiMM t2) - , ninio = f (ninio t1) (ninio t2) - , maxNinio = f (maxNinio t1) (maxNinio t2) - , multiOffset = f (multiOffset t1) (multiOffset t2) - , multiNuc = f (multiNuc t1) (multiNuc t2) - , multiHelix = f (multiHelix t1) (multiHelix t2) - , extMM = zipWithPA f (extMM t1) (extMM t2) - , largeLoop = largeLoop t1 -- no zipping here! TODO make this a constant somewhere! - , termAU = f (termAU t1) (termAU t2) - , intermolecularInit = f (intermolecularInit t1) (intermolecularInit t2) - } - -zipWithPA f a1 a2 = let (l,u) = bounds a1 in fromList l u $ zipWith f (toList a1) (toList a2) diff --git a/Biobase/Vienna.hs b/Biobase/Vienna.hs deleted file mode 100644 --- a/Biobase/Vienna.hs +++ /dev/null @@ -1,89 +0,0 @@ -{-# LANGUAGE NoMonomorphismRestriction #-} - --- | Turner energy data as used in the ViennaRNA package. --- --- TODO should use energy type for tables, not int! - -module Biobase.Vienna where - -import Data.Map (Map,empty) - -import Data.Ix.Tuple -import Data.PrimitiveArray.Ix -import Data.PrimitiveArray -import Biobase.Types.Partition -import Biobase.Types.Energy - -import Biobase.RNA -import Biobase.Turner.Tables -import Biobase.Constants - - - --- Type for Vienna Tables - --- type ViennaEnergySet = (Temperature, ViennaTables, ViennaTables) - -type ViennaIntTables = Turner2004 ViennaPair Nucleotide Int -type ViennaEnergyTables = Turner2004 ViennaPair Nucleotide Energy -type ViennaPartitionTables = Turner2004 ViennaPair Nucleotide Partition - -type Temperature = Double - - - --- | An empty Turner2004 set, with Vienna-style keys for tables. Unused fields --- are commented out. Activate them in BiobaseTurner before doing the same --- here! - -emptyTables :: ViennaIntTables -emptyTables = Turner2004 - --stack - { stack = emptyIV - --dangle - , dangle3 = emptyIV - , dangle5 = emptyIV - --hairpin - , hairpinL = emptyV - , hairpinMM = emptyIV - , hairpinLookup = empty - --eInf - --eInf - --eInf - --eInf - --bulge - , bulgeL = emptyV - --eInf - --internal - , iloop1x1 = emptyIV - , iloop1x2 = emptyIV - , iloop2x2 = emptyIV - , iloopMM = emptyIV - , iloop2x3MM = emptyIV - , iloop1xnMM = emptyIV - , iloopL = emptyV - --multibranch - , multiMM = emptyIV - , ninio = eInf - , maxNinio = eInf - , multiOffset = eInf - , multiNuc = eInf - , multiHelix = eInf - --eInf - --eInf - --exterior - , extMM = emptyIV - --coaxial stacking - --emptyIV - --emptyIV - --emptyIV - --scalar values - , largeLoop = fInf -- large loop is always double - , termAU = eInf - , intermolecularInit = eInf - } - -emptyIV = fromAssocs minBound maxBound eInf [] - --- TODO should be an empty a -emptyV = fromAssocs 0 30 eInf [] diff --git a/Biobase/Vienna/Default.hs b/Biobase/Vienna/Default.hs deleted file mode 100644 --- a/Biobase/Vienna/Default.hs +++ /dev/null @@ -1,27 +0,0 @@ - -{-# LANGUAGE TemplateHaskell #-} - --- | Embedding of a default set of parameters. These are encoded in Vienna style and --- use the Turner2004 set of parameters. - -module Biobase.Vienna.Default - ( turner2004GH - ) where - -import Data.FileEmbed -import qualified Data.ByteString.Char8 as B - -import Biobase.Vienna.Import.ViennaPar -import Biobase.Turner.Tables -import Biobase.Vienna -import Biobase.Types.Energy - - - -defaultPars :: B.ByteString -defaultPars = $(embedFile "parfiles/Turner2004GH.par") - --- Turner entropy and enthalpy tables. - -turner2004GH = (dmap Energy ta,dmap Energy tb) where - (ta,tb) = importPar $ B.unpack defaultPars diff --git a/Biobase/Vienna/Export/ViennaC.hs b/Biobase/Vienna/Export/ViennaC.hs deleted file mode 100644 --- a/Biobase/Vienna/Export/ViennaC.hs +++ /dev/null @@ -1,264 +0,0 @@ -{-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE TemplateHaskell #-} - --- | Exports the Vienna energy tables as a number of C files (one .c file, a --- number .h files). For this reason, files are written into a target --- directory! - -module Biobase.Vienna.Export.ViennaC - ( export - ) where - -import Prelude hiding (pi) -import Text.Printf -import System.Directory -import Data.List.Split -import Data.List -import qualified Data.Map as M -import Data.FileEmbed -import qualified Data.ByteString.Char8 as B - -import Biobase.Turner.Tables -import Biobase.Vienna -import Data.PrimitiveArray -import Biobase.RNA - -import Biobase.Vienna.Default - -(g,h) = turner2004GH - - - -type Directory = FilePath - --- | Export the files required for default values in the ViennaRNA 2.0 C --- package. - -export :: ViennaIntTables -> ViennaIntTables -> Directory -> IO () -export trnrG trnrH target = mapM_ (\f -> f trnrG trnrH target) [exportC, export11, export21, export22] - - - --- | Export the main C file - -exportC :: ViennaIntTables -> ViennaIntTables -> Directory -> IO () -exportC trnrG trnrH target = writeFile (target ++ "/energy_par.c") out where - out = unlines [ B.unpack $(embedFile "templates/cheader") , xs, stck, lps, mms, dngls, tris, tets, hexs,includes] - xs = unlines - [ pd "lxc37" $ largeLoop trnrG - , pi "ML_intern37" $ multiHelix trnrG - , pi "ML_interndH" $ multiHelix trnrH - , pi "ML_closing37" $ multiOffset trnrG - , pi "ML_closingdH" $ multiOffset trnrH - , pi "ML_BASE37" $ multiNuc trnrG - , pi "ML_BASEdH" $ multiNuc trnrH - , pi "MAX_NINIO" $ maxNinio trnrG - , pi "ninio37" $ ninio trnrG - , pi "niniodH" $ ninio trnrH - , pi "TerminalAU37" $ termAU trnrG - , pi "TerminalAUdH" $ termAU trnrH - , pi "DuplexInit37" $ intermolecularInit trnrG - , pi "DuplexInitdH" $ intermolecularInit trnrH - , pi "TripleC37" 100 - , pi "TripleCdH" 1860 - , pi "MultipleCA37" 30 - , pi "MultipleCAdH" 340 - , pi "MultipleCB37" 160 - , pi "MultipleCBdH" 760 - ] - stck = unlines - [ "PUBLIC int stack37[NBPAIRS+1][NBPAIRS+1] =" - , (block2 8 . toList $ stack trnrG) - , "PUBLIC int stackdH[NBPAIRS+1][NBPAIRS+1] =" - , (block2 8 . toList $ stack trnrH) - ] - lps = unlines - [ "PUBLIC int hairpin37[31] = {" ++ (concat . intersperse "," . map printNum . toList $ hairpinL trnrG) ++ "};" - , "PUBLIC int hairpindH[31] = {" ++ (concat . intersperse "," . map printNum . toList $ hairpinL trnrH) ++ "};" - , "PUBLIC int bulge37[31] = {" ++ (concat . intersperse "," . map printNum . toList $ bulgeL trnrG) ++ "};" - , "PUBLIC int bulgedH[31] = {" ++ (concat . intersperse "," . map printNum . toList $ bulgeL trnrH) ++ "};" - , "PUBLIC int internal_loop37[31] = {" ++ (concat . intersperse "," . map printNum . toList $ iloopL trnrG) ++ "};" - , "PUBLIC int internal_loopdH[31] = {" ++ (concat . intersperse "," . map printNum . toList $ iloopL trnrH) ++ "};" - ] - mms = unlines - [ writeMM "I" iloopMM trnrG trnrH - , writeMM "H" hairpinMM trnrG trnrH - , writeMM "M" multiMM trnrG trnrH - , writeMM "1nI" iloop1xnMM trnrG trnrH - , writeMM "23I" iloop2x3MM trnrG trnrH - , writeMM "Ext" extMM trnrG trnrH - ] - dngls = unlines - [ "PUBLIC int dangle3_37[NBPAIRS+1][5] =" - , block2 5 . toList $ dangle3 trnrG - , "PUBLIC int dangle3_dH[NBPAIRS+1][5] =" - , block2 5 . toList $ dangle3 trnrH - , "PUBLIC int dangle5_37[NBPAIRS+1][5] =" - , block2 5 . toList $ dangle5 trnrG - , "PUBLIC int dangle5_dH[NBPAIRS+1][5] =" - , block2 5 . toList $ dangle5 trnrH - ] - t3 = getLookup 3 trnrG trnrH - t4 = getLookup 4 trnrG trnrH - t6 = getLookup 6 trnrG trnrH - tris = writeTabbed 3 "Triloop" t3 - tets = writeTabbed 4 "Tetraloop" t4 - hexs = writeTabbed 6 "Hexaloop" t6 - - includes = unlines $ map ((++ "\"") . ("#include \""++)) - [ "intl11.h" - , "intl11dH.h" - , "intl21.h" - , "intl21dH.h" - , "intl22.h" - , "intl22dH.h" - ] - - - --- | Export iloops 1x1 - -export11 :: ViennaIntTables -> ViennaIntTables -> Directory -> IO () -export11 trnrG trnrH target = writeFile (target ++ "/intl11.h") outG >> writeFile (target ++ "/intl11dH.h") outH where - outG = unlines - [ "PUBLIC int int11_37[NBPAIRS+1][NBPAIRS+1][5][5] =" - , (block4 8 5 5 . toList $ iloop1x1 trnrG) - ] - outH = unlines - [ "PUBLIC int int11_dH[NBPAIRS+1][NBPAIRS+1][5][5] =" - , (block4 8 5 5 . toList $ iloop1x1 trnrH) - ] - - - --- | Export iloops 2x1 - -export21 :: ViennaIntTables -> ViennaIntTables -> Directory -> IO () -export21 trnrG trnrH target = writeFile (target ++ "/intl21.h") outG >> writeFile (target ++ "/intl21dH.h") outH where - outG = unlines - [ "PUBLIC int int21_37[NBPAIRS+1][NBPAIRS+1][5][5][5] =" - , (block5 8 5 5 5 . toList $ iloop1x2 trnrG) - ] - outH = unlines - [ "PUBLIC int int21_dH[NBPAIRS+1][NBPAIRS+1][5][5][5] =" - , (block5 8 5 5 5 . toList $ iloop1x2 trnrH) - ] - - - --- | Export iloops 2x2 - -export22 :: ViennaIntTables -> ViennaIntTables -> Directory -> IO () -export22 trnrG trnrH target = writeFile (target ++ "/intl22.h") outG >> writeFile (target ++ "/intl22dH.h") outH where - outG = unlines - [ "PUBLIC int int22_37[NBPAIRS+1][NBPAIRS+1][5][5][5][5] =" - , (block6 8 5 5 5 5 . toList $ iloop2x2 trnrG) - ] - outH = unlines - [ "PUBLIC int int22_dH[NBPAIRS+1][NBPAIRS+1][5][5][5][5] =" - , (block6 8 5 5 5 5 . toList $ iloop2x2 trnrH) - ] - - - --- * Helper functions - -writeTabbed (ksize :: Int) key (ns,gs,hs) = unlines - [ printf "PUBLIC char %ss[%d] =" key ((ksize+3)*40+1) - , ns - , ";" - ] ++ printf "PUBLIC int %s37[40] = " key ++ gs ++ ";\n" - ++ printf "PUBLIC int %sdH[40] = " key ++ hs ++ ";\n" - -getLookup k trnrG trnrH = (nsF,gs,hs) where - nsF = concat $ intersperse "\n" $ map ((++ " \"") . (" \"" ++) . map nucleotideToChar) ns - ns = filter ((==(k+2)) . length) . M.keys $ hairpinLookup trnrG - gs = lkup trnrG - hs = lkup trnrH - lkup tbl = ("{"++) . (++"}") . concat . intersperse "," $ map (printNum . (hairpinLookup tbl M.!)) ns - -writeMM typ fun trnrG trnrH = unlines - [ printf "PUBLIC int mismatch%s37[NBPAIRS+1][5][5] =" typ - , block3 5 5 . toList $ fun trnrG - , printf "PUBLIC int mismatch%sdH[NBPAIRS+1][5][5] =" typ - , block3 5 5 . toList $ fun trnrH - ] - -block2 s2 xs - = (++ "};") - . ("{" ++) - . concat - . intersperse "\n," - . map ((++ "}") . ("{" ++) . concat . intersperse ",") - . splitEvery s2 - $ map printNum xs - -block3 s2 s3 xs - = (++ "};") - . ("{" ++) - . concat - . intersperse "\n," - . map ((++ "}") . ("{" ++) . concat . intersperse ",") - . splitEvery s2 - . map ((++ "}\n ") . ("{" ++) . concat . intersperse ",") - . splitEvery s3 - $ map printNum xs - -block4 s2 s3 s4 xs - = (++ "};") - . ("{" ++) - . concat - . intersperse "\n," - . map ((++ "}") . ("{" ++) . concat . intersperse ",") - . splitEvery s2 - . map ((++ "}\n ") . ("{" ++) . concat . intersperse ",") - . splitEvery s3 - . map ((++ "}\n ") . ("{" ++) . concat . intersperse ",") - . splitEvery s4 - $ map printNum xs - -block5 s2 s3 s4 s5 xs - = (++ "};") - . ("{" ++) - . concat - . intersperse "\n," - . map ((++ "}") . ("{" ++) . concat . intersperse ",") - . splitEvery s2 - . map ((++ "}\n ") . ("{" ++) . concat . intersperse ",") - . splitEvery s3 - . map ((++ "}\n ") . ("{" ++) . concat . intersperse ",") - . splitEvery s4 - . map ((++ "}\n ") . ("{" ++) . concat . intersperse ",") - . splitEvery s5 - $ map printNum xs - -block6 s2 s3 s4 s5 s6 xs - = (++ "};") - . ("{" ++) - . concat - . intersperse "\n," - . map ((++ "}") . ("{" ++) . concat . intersperse ",") - . splitEvery s2 - . map ((++ "}\n ") . ("{" ++) . concat . intersperse ",") - . splitEvery s3 - . map ((++ "}\n ") . ("{" ++) . concat . intersperse ",") - . splitEvery s4 - . map ((++ "}\n ") . ("{" ++) . concat . intersperse ",") - . splitEvery s5 - . map ((++ "}\n ") . ("{" ++) . concat . intersperse ",") - . splitEvery s6 - $ map printNum xs - - -pd :: String -> Double -> String -pd k v = printf "PUBLIC double %s=%f;" k v - -pi :: String -> Int -> String -pi k v = printf "PUBLIC int %s=%d;" k v - -printNum :: Int -> String -printNum n - | n >= maxN = " INF" - | otherwise = printf "%6d" n - where - maxN = 10000 diff --git a/Biobase/Vienna/Export/ViennaPar.hs b/Biobase/Vienna/Export/ViennaPar.hs deleted file mode 100644 --- a/Biobase/Vienna/Export/ViennaPar.hs +++ /dev/null @@ -1,139 +0,0 @@ - -{-# LANGUAGE NoMonomorphismRestriction #-} - -module Biobase.Vienna.Export.ViennaPar - ( export - ) where - -import Data.Tuple.Select -import Data.List.Split -import Text.Printf -import Control.Arrow -import Data.List (intersperse) -import qualified Data.Map as M - -import Data.PrimitiveArray -import Biobase.Turner.Tables -import Biobase.RNA -import Biobase.Constants - -import Biobase.Vienna - - - -export :: ViennaIntTables -> ViennaIntTables -> String -export trnr trnrH = hdr ++ blocks ++ mlps ++ ninios ++ misc ++ triloops ++ tetra ++ hexa ++ "\n#END" where - hdr = "## RNAfold parameter file v2.0\n\n" - mlps = printf "# ML_params\n%7d %7d %7d %7d %7d %7d\n\n" - (multiNuc trnr) (multiNuc trnrH) - (multiOffset trnr) (multiOffset trnrH) - (multiHelix trnr) (multiHelix trnrH) - ninios = printf "# NINIO\n%7d %7d %7d\n\n" - (ninio trnr) (ninio trnrH) (maxNinio trnr) - misc = printf "# Misc\n %7d %7d %7d %7d\n\n" - (intermolecularInit trnr) (intermolecularInit trnrH) (termAU trnr) (termAU trnrH) - triloops = "# Triloops\n" ++ printHairpinAssocs 5 trnr trnrH ++ "\n" - tetra = "# Tetraloops\n" ++ printHairpinAssocs 6 trnr trnrH ++ "\n" - hexa = "# Hexaloops\n" ++ printHairpinAssocs 8 trnr trnrH ++ "\n" - blocks = concat $ zipWith (++) - -- entropy terms - [ printBlock "stack" 7 pp2lkey $ stack trnr - , printBlock "mismatch_hairpin" 5 pbb2lkey $ hairpinMM trnr - , printBlock "mismatch_interior" 5 pbb2lkey $ iloopMM trnr - , printBlock "mismatch_interior_1n" 5 pbb2lkey $ iloop1xnMM trnr - , printBlock "mismatch_interior_23" 5 pbb2lkey $ iloop2x3MM trnr - , printBlock "mismatch_multi" 5 pbb2lkey $ multiMM trnr - , printBlock "mismatch_exterior" 5 pbb2lkey $ extMM trnr - , printBlock "dangle5" 5 pb2lkey $ dangle5 trnr - , printBlock "dangle3" 5 pb2lkey $ dangle3 trnr - , printBlock "int11" 5 ppbb2lkey $ iloop1x1 trnr - , printBlock "int21" 5 ppbbb2lkey $ iloop1x2 trnr - , printBlock22 "int22" 4 ppbbbb2lkey $ iloop2x2 trnr - , printLinear "hairpin" 10 $ hairpinL trnr - , printLinear "bulge" 10 $ bulgeL trnr - , printLinear "interior" 10 $ iloopL trnr - ] - -- enthalpy terms - [ printBlockH "stack" 7 pp2lkey $ stack trnrH - , printBlockH "mismatch_hairpin" 5 pbb2lkey $ hairpinMM trnrH - , printBlockH "mismatch_interior" 5 pbb2lkey $ iloopMM trnrH - , printBlockH "mismatch_interior_1n" 5 pbb2lkey $ iloop1xnMM trnrH - , printBlockH "mismatch_interior_23" 5 pbb2lkey $ iloop2x3MM trnrH - , printBlockH "mismatch_multi" 5 pbb2lkey $ multiMM trnrH - , printBlockH "mismatch_exterior" 5 pbb2lkey $ extMM trnrH - , printBlockH "dangle5" 5 pb2lkey $ dangle5 trnrH - , printBlockH "dangle3" 5 pb2lkey $ dangle3 trnrH - , printBlockH "int11" 5 ppbb2lkey $ iloop1x1 trnrH - , printBlockH "int21" 5 ppbbb2lkey $ iloop1x2 trnrH - , printBlock22H "int22" 4 ppbbbb2lkey $ iloop2x2 trnrH - , printLinearH "hairpin" 10 $ hairpinL trnrH - , printLinearH "bulge" 10 $ bulgeL trnrH - , printLinearH "interior" 10 $ iloopL trnrH - ] - - - --- * Helper functions - --- | Show the key of the line, minus the changing last key - -showKey :: [(LKey,Int)] -> String -showKey xs = - " /* " ++ - (concat $ intersperse "," $ init $ (map show ps) ++ (map show ns)) ++ - " */" - where - (ps,ns) = fst $ head xs - --- | Transform from tuple-based keys to a pair of list-based keys. - -type LKey = ([ViennaPair],[Nucleotide]) - -pb2lkey (p1,b1) = ([p1],[b1]) -pbb2lkey (p1,b1,b2) = ([p1],[b1,b2]) -pp2lkey (k1,k2) = ([k1,k2],[]) -ppbb2lkey (p1,p2,b1,b2) = ([p1,p2],[b1,b2]) -ppbbb2lkey (p1,p2,b1,b2,b3) = ([p1,p2],[b1,b2,b3]) -ppbbbb2lkey (p1,p2,b1,b2,b3,b4) = ([p1,p2],[b1,b2,b3,b4]) - --- | Print a block. - -printBlock = printBlockG noNP where - noNP ((ps,ns),v) = not $ any (==vpNP) ps - -printBlockH s = printBlock (s ++ "_enthalpies") - -printBlock22 = printBlockG noNSNPE where - noNSNPE ((ps,ns),v) = not $ any (==vpNP) ps || any (==vpNS) ps || any (==nucE) ns - -printBlock22H s = printBlock22 (s ++ "_enthalpies") - -printBlockG fltr s k tolkey xs' = - let - xs = filter fltr $ map (first tolkey) $ assocs xs' - in - printf "# %s\n" s ++ - (concatMap printLine $ splitEvery k xs) ++ - "\n" - -printLine xs = - concatMap printVal xs ++ " " ++ showKey xs ++ - printf "\n" - -printVal (k,v) - | v > eMax = printf " INF" - | otherwise = printf "%6d" v - --- | A linear block is more boring - -printLinear s k xs' = let xs = assocs xs' in - printf "# %s\n" s ++ - (concatMap (\ys -> concatMap printVal ys ++ "\n") $ splitEvery k xs) ++ - "\n" - -printHairpinAssocs l trnr trnrH = res where - res = concat $ zipWith (\(k,v) vH -> printf "%s %7d %7d\n" (concatMap show k) v vH) xs ys - xs = filter ((==l).length.fst) $ M.assocs $ hairpinLookup trnr - ys = map snd $ filter ((==l).length.fst) $ M.assocs $ hairpinLookup trnrH - -printLinearH s = printLinear (s ++ "_enthalpies") diff --git a/Biobase/Vienna/Import/Turner.hs b/Biobase/Vienna/Import/Turner.hs deleted file mode 100644 --- a/Biobase/Vienna/Import/Turner.hs +++ /dev/null @@ -1,71 +0,0 @@ -{-# LANGUAGE RecordWildCards #-} - --- | Import data from a Turner table structure. Those can be imported using the --- Turner.Import module. - -module Biobase.Vienna.Import.Turner - ( fromTurner - ) where - -import Control.Arrow -import Data.Tuple.All -import qualified Data.Map as M - -import Biobase.Turner.Tables -import Biobase.Turner -import Data.PrimitiveArray -import Biobase.Constants -import Biobase.RNA - -import Biobase.Vienna --- | Convert from Turner keys to Vienna keys and all values from Double to Int - -fromTurner :: TurnerTables -> ViennaIntTables -fromTurner Turner2004{..} = Turner2004 - { stack = convert pp stack - , dangle3 = convert pb dangle3 - , dangle5 = convert pb dangle5 - , hairpinL = amap f2i hairpinL - , hairpinMM = convert pbb hairpinMM - , hairpinLookup = M.map f2i hairpinLookup - , bulgeL = amap f2i bulgeL - , iloop1x1 = convert uu iloop1x1 - , iloop1x2 = convert uu iloop1x2 - , iloop2x2 = convert uu iloop2x2 - , iloopMM = convert pbb iloopMM - , iloop2x3MM = convert pbb iloop2x3MM - , iloop1xnMM = convert pbb iloop1xnMM - , iloopL = amap f2i iloopL - , multiMM = convert pbb multiMM - , ninio = f2i ninio - , maxNinio = f2i maxNinio - , multiOffset = f2i multiOffset - , multiNuc = f2i multiNuc - , multiHelix = f2i multiHelix - , extMM = convert pbb extMM - , largeLoop = largeLoop * 100 -- scale to decacals - , termAU = f2i termAU - , intermolecularInit = f2i intermolecularInit - } - - - --- * Helper functions. - --- | Vienna values are Int's that have 2-digit precision. - -f2i :: Double -> Int -f2i d = round $ d * 100 - - --- | Deep breath ;-) - -pp = both tuple2Vienna -pb = first tuple2Vienna -pbb t = upd1 (tuple2Vienna $ sel1 t) t -uu t = upd2 (tuple2Vienna $ sel2 t) $ upd1 (tuple2Vienna $ sel1 t) t -both f = first f . second f - --- | Do the conversion, both on keys and values. - -convert fk = fromAssocs minBound maxBound eInf . map (first fk . second f2i) . assocs diff --git a/Biobase/Vienna/Import/ViennaPar.hs b/Biobase/Vienna/Import/ViennaPar.hs deleted file mode 100644 --- a/Biobase/Vienna/Import/ViennaPar.hs +++ /dev/null @@ -1,123 +0,0 @@ -{-# LANGUAGE PatternGuards #-} -{-# LANGUAGE RecordWildCards #-} - -module Biobase.Vienna.Import.ViennaPar - ( importPar - ) where - -import Text.Parsec hiding ((<|>), many, optional) -import Text.Parsec.String -import Text.Parsec.Char -import Text.Parsec.Combinator hiding (optional) -import Control.Applicative -import Data.List -import qualified Data.Map as M - -import Text.Parsec.Numbers -import Biobase.Constants -import Biobase.RNA -import Biobase.Turner.Tables -import Data.PrimitiveArray - -import Biobase.Vienna - - - --- | Imports Vienna tables from a given string. - -importPar :: String -> (ViennaIntTables,ViennaIntTables) -importPar s - | Right _ <- p = foldr updateTables (trnr,trnrH) bs - | Left err <- p = error $ show err - where - p = parse blocks "(stdin)" s - Right bs = p - (bgs,bhs) = partition (not . isHblock) bs - trnr = foldr updateTable emptyTables bgs - trnrH = foldr updateTable emptyTables bhs - --- | Apply all parsed blocks to a Turner2004 structure. - -updateTable :: Block -> ViennaIntTables -> ViennaIntTables -updateTable Block{..} trnr -- @Turner2004{..} - | bhdr =? "stack" = trnr{stack = fromAssocs minBound maxBound z $ zip pp ints} - | bhdr =? "dangle3" = trnr{dangle3 = fromAssocs minBound maxBound z $ zip pb ints} - | bhdr =? "dangle5" = trnr{dangle5 = fromAssocs minBound maxBound z $ zip pb ints} - | bhdr =? "hairpin" = trnr{hairpinL = fromList 0 30 ints} - | bhdr =? "mismatch_hairpin" = trnr{hairpinMM = fromAssocs minBound maxBound z $ zip pbb ints} - | bhdr =? "bulge" = trnr{bulgeL = fromList 0 30 ints} - | bhdr =? "int11" = trnr{iloop1x1 = fromAssocs minBound maxBound z $ zip ppbb ints} - | bhdr =? "int21" = trnr{iloop1x2 = fromAssocs minBound maxBound z $ zip ppbbb ints} - | bhdr =? "int22" = trnr{iloop2x2 = fromAssocs minBound maxBound z $ zip ppbbbb ints} - | bhdr =? "mismatch_interior" = trnr{iloopMM = fromAssocs minBound maxBound z $ zip pbb ints} - | bhdr =? "mismatch_interior_23" = trnr{iloop2x3MM = fromAssocs minBound maxBound z $ zip pbb ints} - | bhdr =? "mismatch_interior_1n" = trnr{iloop1xnMM = fromAssocs minBound maxBound z $ zip pbb ints} - | bhdr =? "interior" = trnr{iloopL = fromList 0 30 ints} - | bhdr =? "mismatch_multi" = trnr{multiMM = fromAssocs minBound maxBound z $ zip pbb ints} - | bhdr =? "mismatch_exterior" = trnr{extMM = fromAssocs minBound maxBound z $ zip pbb ints} - | otherwise = trnr -- unknown block! -updateTable ABlock{..} trnr - | otherwise = trnr -- CORRECT! - --- | update with info where no separation between G/H was done. - -updateTables Block{..} (trnr,trnrH) - | bhdr =? "ML_params" - = ( trnr {multiNuc = ints!!0, multiOffset = ints!!2, multiHelix = ints!!4} - , trnrH{multiNuc = ints!!1, multiOffset = ints!!3, multiHelix = ints!!5} - ) - | bhdr =? "NINIO" - = ( trnr {ninio = ints!!0, maxNinio = ints!!2} - , trnrH {ninio = ints!!1} - ) - | bhdr =? "Misc" - = ( trnr {intermolecularInit = ints!!0, termAU = ints!!2} - , trnrH {intermolecularInit = ints!!1, termAU = ints!!3} - ) - | otherwise = (trnr,trnrH) -updateTables ABlock{..} (trnr,trnrH) - | bhdr =? "Triloops" || bhdr =? "Tetraloops" || bhdr =? "Hexaloops" - = ( trnr {hairpinLookup = hairpinLookup trnr `M.union` (M.fromList $ map (\(k,v,w) -> (s2ns k,v)) ascs)} - , trnrH {hairpinLookup = hairpinLookup trnrH `M.union` (M.fromList $ map (\(k,v,w) -> (s2ns k,w)) ascs)} - ) - | otherwise = (trnr,trnrH) - -pp = [(x,y) | x<-cgnsP,y<-cgnsP] -pb = [(x,y) | x<-cgnsP,y<-eacgu] -pbb = [(x,y,z) | x<-cgnsP,y<-eacgu,z<-eacgu] -ppbb = [(p1,p2,b1,b2) | p1<-cgnsP,p2<-cgnsP,b1<-eacgu,b2<-eacgu] -ppbbb = [(p1,p2,b1,b2,b3) | p1<-cgnsP,p2<-cgnsP,b1<-eacgu,b2<-eacgu,b3<-eacgu] -ppbbbb = [(p1,p2,b1,b2,b3,b4) | p1<-cguaP,p2<-cguaP,b1<-acgu,b2<-acgu,b3<-acgu,b4<-acgu] -z = eInf --- xs is the block header name, ys the constant string -xs =? ys = xs == ys || xs == ys++"_enthalpies" -- and $ zipWith (==) xs ys -s2ns = map charToNucleotide - --- * Simple parser for the vienna 2.0 format. - -isHblock :: Block -> Bool -isHblock b = "enthalpies" `elem` (tails $ bhdr b) - -data Block - = Block {bhdr :: String, ints :: [Int]} - | ABlock {bhdr :: String, ascs :: [(String,Int,Int)]} - deriving (Show) - -blocks :: GenParser Char st [Block] -blocks = vrna2 *> spaces *> (try assocBlock <|> try block) `sepEndBy1` spaces <* string "#END" <* spaces <* eof - -vrna2 = string "## RNAfold parameter file v2.0" - -block :: GenParser Char st Block -block = Block <$> header <* spaces <*> num `sepEndBy1` (try comment <|> spaces) where - num = ((eInf <$ string "INF") <|> parseIntegral) - -assocBlock :: GenParser Char st Block -assocBlock = ABlock <$> header <* spaces <*> asc `sepEndBy1` spaces where - asc = (,,) <$> (many1 $ oneOf "ACGU") <* spaces <*> parseIntegral <* spaces <*> parseIntegral - -header :: GenParser Char st String -header = char '#' *> space *> (alphaNum <|> char '_') `manyTill` newline - -comment :: GenParser Char st () -comment = const () <$ spaces <*> between (string "/*") (string "*/") (many $ noneOf "/*") <* spaces diff --git a/Biobase/Vienna/Modification/NonStandard.hs b/Biobase/Vienna/Modification/NonStandard.hs deleted file mode 100644 --- a/Biobase/Vienna/Modification/NonStandard.hs +++ /dev/null @@ -1,121 +0,0 @@ - --- | This module provides facilities to modify values in the tables and return --- the modified tables. The main reason for the modules' existence are the 'E' --- nucleotides and 'NS' base pairs used by the ViennaRNA package to denote --- unknown or missing nucleotides. For each table type, we have a small --- function. - -module Biobase.Vienna.Modification.NonStandard - ( nonStandardMax - ) where - -import Biobase.RNA -import Biobase.Vienna -import Biobase.Turner.Tables -import Data.Ix.Tuple -import Biobase.Constants -import Biobase.Types.Energy - -import Data.PrimitiveArray - - --- | Each table-changing operation makes the promise that all individual --- changes in a batch are independent of the changes in the same batch - -changeTable f tbl = fromAssocs l u z . map (f tbl) $ assocs tbl where - (l,u) = bounds tbl - z = Energy eInf - - - --- |* Replace each occurrance of a value with a strange key according to this --- rule. The rule is to take the maximum over all sane keys in the --- neighborhood. Yes, this takes long... - -fPPBBBBmax tbl (k@(p1,p2,b1,b2,b3,b4),v) - | p1==vpNP || p2==vpNP = (k,v) -- don't do anything - | otherwise = (k,maximum' - [ tbl ! (p1n,p2n,b1n,b2n,b3n,b4n) - | p1n <- genKeyP p1 - , p2n <- genKeyP p2 - , b1n <- genKeyB b1 - , b2n <- genKeyB b2 - , b3n <- genKeyB b3 - , b4n <- genKeyB b4 - ]) - -fPPBBBmax tbl (k@(p1,p2,b1,b2,b3),v) - | p1==vpNP || p2==vpNP = (k,v) - | otherwise = (k,maximum' - [ tbl ! (p1n,p2n,b1n,b2n,b3n) - | p1n <- genKeyP p1 - , p2n <- genKeyP p2 - , b1n <- genKeyB b1 - , b2n <- genKeyB b2 - , b3n <- genKeyB b3 - ]) - -fPPBBmax tbl (k@(p1,p2,b1,b2),v) - | p1==vpNP || p2==vpNP = (k,v) - | otherwise = (k,maximum' - [ tbl ! (p1n,p2n,b1n,b2n) - | p1n <- genKeyP p1 - , p2n <- genKeyP p2 - , b1n <- genKeyB b1 - , b2n <- genKeyB b2 - ]) - -fPPmax tbl (k@(p1,p2),v) - | p1==vpNP || p2==vpNP = (k,v) - | otherwise = (k,maximum' - [ tbl ! (p1n,p2n) - | p1n <- genKeyP p1 - , p2n <- genKeyP p2 - ]) - -fPBBmax tbl (k@(p1,b1,b2),v) - | p1==vpNP = (k,v) - | otherwise = (k,maximum' - [ tbl ! (p1n,b1n,b2n) - | p1n <- genKeyP p1 - , b1n <- genKeyB b1 - , b2n <- genKeyB b2 - ]) - -fPBmax tbl (k@(p1,b1),v) - | p1==vpNP = (k,v) - | otherwise = (k,maximum' - [ tbl ! (p1n,b1n) - | p1n <- genKeyP p1 - , b1n <- genKeyB b1 - ]) - -maximum' = Energy . maximum . map unEnergy - - --- * Generate keys. If the key is of the non-standard kind, then we go fishing, --- otherwise we keep the key itself. - -genKeyP k = if k==vpNS then cguaP else [k] -genKeyB k = if k==nucE then acgu else [k] - - - --- | The max operation applied to all relevant tables. (This serves as a --- pointer, which tables to manipulate). - -nonStandardMax :: ViennaEnergyTables -> ViennaEnergyTables -nonStandardMax tbl = tbl - { iloop1x1 = changeTable fPPBBmax $ iloop1x1 tbl - , iloop1x2 = changeTable fPPBBBmax $ iloop1x2 tbl - , iloop2x2 = changeTable fPPBBBBmax $ iloop2x2 tbl - , hairpinMM = changeTable fPBBmax $ hairpinMM tbl - , iloopMM = changeTable fPBBmax $ iloopMM tbl - , iloop1xnMM = changeTable fPBBmax $ iloop1xnMM tbl - , iloop2x3MM = changeTable fPBBmax $ iloop2x3MM tbl - , multiMM = changeTable fPBBmax $ multiMM tbl - , extMM = changeTable fPBBmax $ extMM tbl - , dangle3 = changeTable fPBmax $ dangle3 tbl - , dangle5 = changeTable fPBmax $ dangle5 tbl - , stack = changeTable fPPmax $ stack tbl - } diff --git a/Biobase/Vienna/Modification/Partition.hs b/Biobase/Vienna/Modification/Partition.hs deleted file mode 100644 --- a/Biobase/Vienna/Modification/Partition.hs +++ /dev/null @@ -1,51 +0,0 @@ - --- | Small module for converting Vienna energy tables to partition function --- value tables. --- --- Temperature scaling should happen on the level of energy tables, followed by --- conversion here. --- --- TODO after the switch to the Energy type, change too, the partitionTables --- generation - -module Biobase.Vienna.Modification.Partition - ( partitionTables - ) where - -import Biobase.Types.Partition -import Biobase.Types.Energy -import Biobase.Turner.Tables -import Biobase.Vienna -import Biobase.Constants - - --- Testing - -{- -import Biobase.Vienna.Default -import Data.PrimitiveArray -import Debug.Trace.Tools -import Biobase.RNA -import Biobase.Types.Ring - -(g,h) = turner2004GH - --- x = toList $ stack $ partitionTables 37 g -x = {- toList $ -} (stack $ partitionTables 37 g) ! (vpGC,vpGC) --} - --- | Create partition tables out of normal energy tables. Note that these --- values are not normalized. An additional term 1/Z(1/kt) would be required --- which is missing until we are done calculating the partition function. --- --- All energy contributions are rescaled from dekacals back to kcal/mol. --- --- We explicitly set all probabilities =0 where the energy is >= eMax. - -partitionTables :: Temperature -> ViennaEnergyTables -> ViennaPartitionTables -partitionTables tempe trnr = dmap f trnr where - kt = (tempe + kelvinC0) * gasconst - -- convert energy into probability; encapsulate as a (logFloat) partition - f (Energy e') = {- traceVal (show (e',e,p)) -} p where - e = fromIntegral e' / 100 - p = if e' >= eMax then Partition 0 else Partition . exp $ -e / kt -- can NOT use mkPartition as these values are NOT NORMALIZED diff --git a/Biobase/Vienna/Modification/Temperature.hs b/Biobase/Vienna/Modification/Temperature.hs deleted file mode 100644 --- a/Biobase/Vienna/Modification/Temperature.hs +++ /dev/null @@ -1,34 +0,0 @@ - --- | Given entropy and enthalpy sets, allow modification of the folding --- temperature (from the default 37 Celsius) - -module Biobase.Vienna.Modification.Temperature - ( adjustTemperature - ) where - -import Biobase.Vienna -import Biobase.Turner.Tables -import Biobase.Constants - - - --- temperature scaling is a generalized zip: let g and h be a scalar, t the --- temperature. Then: g_t = h - (h-g)*((t+k_0)/t_m). Where k_0 = 273.15 (Kelvin --- of 37 Celsius) and t_m = 37+k_0 --- --- TODO abstract over table type? --- --- TODO this does not take into account that dangles are supposed to _always_ --- be beneficial. We should therefor transform the energies and then apply 'min --- 0' onto all dangle types. --- --- TODO use EnergyTables (Int is just for im-/export) - -adjustTemperature :: ViennaIntTables -> ViennaIntTables -> Temperature -> ViennaIntTables -adjustTemperature trnrG trnrH cels = dZipWith f trnrG trnrH where - f vG vH = - let g = fromIntegral vG - h = fromIntegral vH - k_0 = kelvinC0 - t_m = 37 + k_0 - in round $ h - (h-g) * ((cels+k_0)/t_m) diff --git a/templates/cheader b/templates/cheader --- a/templates/cheader +++ b/templates/cheader @@ -2,15 +2,14 @@ /* Automatically generated using the TurnerParser - TurnerParser (c) 2008,2009,2010 + TurnerParser (c) 2008,2009,2010,2011 Christian Hoener zu Siederdissen, TBI Vienna choener (at) tbi.univie.ac.at - The library enabling this can be found at: - http://hackage.haskell.org/package/BiobaseVienna - the program can be found at: - (sorry, not yet) - install using cabal: cabal install (sorry, not yet) + Important links: + http://hackage.haskell.org/package/Biobase + http://www.tbi.univie.ac.at/~choener/ + http://www.biohaskell.org/ */ /*