packages feed

Biobase 0.1.0.1 → 0.3.1.0

raw patch · 52 files changed

+3502/−2154 lines, 52 filesdep +deepseqdep +either-unwrapdep +filepathdep ~PrimitiveArray

Dependencies added: deepseq, either-unwrap, filepath, ghc-prim, mtl, utility-ht, vector-read-instances, zlib

Dependency ranges changed: PrimitiveArray

Files

Biobase.cabal view
@@ -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. <http://rna.urmc.rochester.edu/NNDB/index.html>                 .+                * 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
Biobase/Constants.hs view
@@ -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
+ Biobase/DataSource/CT.hs view
@@ -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
+ Biobase/DataSource/CT/Import.hs view
@@ -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
+ Biobase/DataSource/DotParen.hs view
@@ -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]
+ Biobase/DataSource/DotParen/Import.hs view
@@ -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+
+ Biobase/DataSource/FR3D.hs view
@@ -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}
+ Biobase/DataSource/FR3D/Import.hs view
@@ -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
+ Biobase/DataSource/InfernalCM.hs view
@@ -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
+ Biobase/DataSource/InfernalCM/Export.hs view
@@ -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
+ Biobase/DataSource/InfernalCM/Import.hs view
@@ -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+-}
+ Biobase/DataSource/MCFold.hs view
@@ -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"+  ]
+ Biobase/DataSource/MCFold/Conversions.hs view
@@ -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) ]
+ Biobase/DataSource/MCFold/Import.hs view
@@ -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)
+ Biobase/DataSource/Mafft.hs view
@@ -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)
+ Biobase/DataSource/Mafft/Import.hs view
@@ -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
+ Biobase/DataSource/Turner.hs view
@@ -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:+-- <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)
+ Biobase/DataSource/Turner/Import.hs view
@@ -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"
+ Biobase/DataSource/Turner/Tables.hs view
@@ -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)
+ Biobase/DataSource/Vienna.hs view
@@ -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 []
+ Biobase/DataSource/Vienna/Default.hs view
@@ -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
+ Biobase/DataSource/Vienna/Export/ViennaC.hs view
@@ -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
+ Biobase/DataSource/Vienna/Export/ViennaPar.hs view
@@ -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")
+ Biobase/DataSource/Vienna/Import/Turner.hs view
@@ -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
+ Biobase/DataSource/Vienna/Import/ViennaPar.hs view
@@ -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
+ Biobase/DataSource/Vienna/Modification/NonStandard.hs view
@@ -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+  }
+ Biobase/DataSource/Vienna/Modification/Partition.hs view
@@ -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
+ Biobase/DataSource/Vienna/Modification/Temperature.hs view
@@ -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)
− Biobase/Infernal/CM.hs 
@@ -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
− Biobase/Infernal/CM/Export.hs 
@@ -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
− Biobase/Infernal/CM/Import.hs 
@@ -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--}
− Biobase/Infernal/Stockholm.hs 
@@ -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)
− Biobase/Infernal/Stockholm/Import.hs 
@@ -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)
Biobase/RNA.hs view
@@ -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 #-}
+ Biobase/RNA/Complex.hs view
@@ -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` ()
+ Biobase/RNA/Hashes.hs view
@@ -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 #-}
+ Biobase/RNA/NucBounds.hs view
@@ -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)
+ Biobase/RNA/Pairs.hs view
@@ -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` ()
+ Biobase/RNA/ViennaPair.hs view
@@ -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")]
− Biobase/Turner.hs 
@@ -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)
− Biobase/Turner/Import.hs 
@@ -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"
− Biobase/Turner/Tables.hs 
@@ -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)
− Biobase/Vienna.hs 
@@ -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 []
− Biobase/Vienna/Default.hs 
@@ -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
− Biobase/Vienna/Export/ViennaC.hs 
@@ -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
− Biobase/Vienna/Export/ViennaPar.hs 
@@ -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")
− Biobase/Vienna/Import/Turner.hs 
@@ -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
− Biobase/Vienna/Import/ViennaPar.hs 
@@ -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
− Biobase/Vienna/Modification/NonStandard.hs 
@@ -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-  }
− Biobase/Vienna/Modification/Partition.hs 
@@ -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
− Biobase/Vienna/Modification/Temperature.hs 
@@ -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)
templates/cheader view
@@ -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/ */  /*