BiobaseTurner 0.0.2.2 → 0.3.1.1
raw patch · 4 files changed
Files
- Biobase/Turner.hs +129/−28
- Biobase/Turner/Import.hs +194/−121
- Biobase/Turner/Tables.hs +0/−133
- BiobaseTurner.cabal +40/−19
Biobase/Turner.hs view
@@ -1,41 +1,142 @@-{-# LANGUAGE TypeSynonymInstances #-}+{-# LANGUAGE RecordWildCards #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE StandaloneDeriving #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE NoMonomorphismRestriction #-}+{-# LANGUAGE TemplateHaskell #-}+{-# LANGUAGE TypeOperators #-} --- | 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.+-- | The 'Turner2004' data structure reflects the RNA (and DNA) energy+-- parameters known as the Turner 2004 data set.+--+-- In general, have a look here:+-- <http://rna.urmc.rochester.edu/NNDB/turner04/index.html> where parameters+-- are explained.+--+-- TODO need a "Functor" instance over elements "e". Or alternatively, generic+-- programming to capture stuff going on in 'e' module Biobase.Turner where -import Biobase.RNA (Nucleotide)-import qualified Biobase.Turner.Tables as T +import Control.Lens+import Data.Array.Repa.Index+import Data.ByteString (ByteString(..))+import qualified Data.ByteString+import qualified Data.Map as M+import qualified Data.Vector as V+import qualified Data.Vector.Unboxed as VU+import qualified Data.Vector.Generic as VG+import qualified Data.Vector.Generic.Mutable as VGM+import Data.Primitive.Types +import Biobase.Primary+import Biobase.Secondary+import Data.PrimitiveArray as PA+import Data.PrimitiveArray.Zero --- | 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. +-- | The actual Turner parameters return energies in Double format. --- 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+newtype Energy = Energy Double+ deriving (Eq,Ord,Num,Read,Show) -type TurnerTables = T.Turner2004 Basepair Nucleotide Double+deriving instance Prim Energy+deriving instance VGM.MVector VU.MVector Energy+deriving instance VG.Vector VU.Vector Energy+deriving instance VU.Unbox Energy -type TurnerEntropy = TurnerTables-type TurnerEnthalpy = TurnerTables-type Temperature = Double-type TurnerSet = (Temperature,TurnerEntropy,TurnerEnthalpy)+-- | The Turner model with 'Energy's. -type Basepair = (Nucleotide,Nucleotide)+type Turner2004 = Turner2004Model Energy++-- | The Turner energy tables. Parametrized over the storing vector type 'v'+-- and the actual element type 'e'.++data Turner2004Model e = Turner2004Model+ { _stack :: !(Unboxed PP e)+ , _dangle3 :: !(Unboxed PN e)+ , _dangle5 :: !(Unboxed PN e)+ , _hairpinL :: !(VU.Vector e)+ , _hairpinMM :: !(Unboxed PNN e)+ , _hairpinLookup :: !(M.Map Primary e)+ , _hairpinGGG :: !e+ , _hairpinCslope :: !e+ , _hairpinCintercept :: !e+ , _hairpinC3 :: !e+ , _bulgeL :: !(VU.Vector e)+ , _bulgeSingleC :: !e+ , _iloop1x1 :: !(Unboxed PPNN e)+ , _iloop2x1 :: !(Unboxed PPNNN e)+ , _iloop2x2 :: !(Unboxed PPNNNN e)+ , _iloopMM :: !(Unboxed PNN e)+ , _iloop2x3MM :: !(Unboxed PNN e)+ , _iloop1xnMM :: !(Unboxed PNN e)+ , _iloopL :: !(VU.Vector e)+ , _multiMM :: !(Unboxed PNN e)+ , _ninio :: !e+ , _maxNinio :: !e+ , _multiOffset :: !e+ , _multiNuc :: !e+ , _multiHelix :: !e+ , _multiAsym :: !e+ , _multiStrain :: !e+ , _extMM :: !(Unboxed PNN e)+ , _coaxial :: !(Unboxed PP e) -- no intervening unpaired nucleotides+ , _coaxStack :: !(Unboxed PNN e)+ , _tStackCoax :: !(Unboxed PNN e)+ , _largeLoop :: !e+ , _termAU :: !e+ , _intermolecularInit :: !e+ } deriving (Show)++type PP = (Z:.Nuc:.Nuc:.Nuc:.Nuc)+type PN = (Z:.Nuc:.Nuc:.Nuc)+type PNN = (Z:.Nuc:.Nuc:.Nuc:.Nuc)+type PPNN = PP:.Nuc:.Nuc+type PPNNN = PPNN:.Nuc+type PPNNNN = PPNNN:.Nuc++makeLenses ''Turner2004Model++-- | Map a function over all 'e' elements.++emap :: (VU.Unbox e, VU.Unbox e') => (e -> e') -> Turner2004Model e -> Turner2004Model e'+emap f Turner2004Model{..} = Turner2004Model+ { _stack = PA.map f _stack+ , _dangle3 = PA.map f _dangle3+ , _dangle5 = PA.map f _dangle5+ , _hairpinL = VU.map f _hairpinL+ , _hairpinMM = PA.map f _hairpinMM+ , _hairpinLookup = M.map f _hairpinLookup+ , _hairpinGGG = f _hairpinGGG+ , _hairpinCslope = f _hairpinCslope+ , _hairpinCintercept = f _hairpinCintercept+ , _hairpinC3 = f _hairpinC3+ , _bulgeL = VU.map f _bulgeL+ , _bulgeSingleC = f _bulgeSingleC+ , _iloop1x1 = PA.map f _iloop1x1+ , _iloop2x1 = PA.map f _iloop2x1+ , _iloop2x2 = PA.map f _iloop2x2+ , _iloopMM = PA.map f _iloopMM+ , _iloop2x3MM = PA.map f _iloop2x3MM+ , _iloop1xnMM = PA.map f _iloop1xnMM+ , _iloopL = VU.map f _iloopL+ , _multiMM = PA.map f _multiMM+ , _ninio = f _ninio+ , _maxNinio = f _maxNinio+ , _multiOffset = f _multiOffset+ , _multiNuc = f _multiNuc+ , _multiHelix = f _multiHelix+ , _multiAsym = f _multiAsym+ , _multiStrain = f _multiStrain+ , _extMM = PA.map f _extMM+ , _coaxial = PA.map f _coaxial+ , _coaxStack = PA.map f _coaxStack+ , _tStackCoax = PA.map f _tStackCoax+ , _largeLoop = f _largeLoop+ , _termAU = f _termAU+ , _intermolecularInit = f _intermolecularInit+ }+
Biobase/Turner/Import.hs view
@@ -1,8 +1,10 @@+{-# LANGUAGE PatternGuards #-}+{-# LANGUAGE TupleSections #-} {-# LANGUAGE NoMonomorphismRestriction #-}+{-# LANGUAGE OverloadedStrings #-} --- | 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.+-- | Turner file parser. Returns a Turner2004 data structure. Requires an+-- annoying amount of boilerplate. -- -- How is 'stack' data stored: --@@ -32,145 +34,216 @@ -- --TODO not sure if dangle3/dangle5 are correctly split or if they should switch -module Biobase.Turner.Import- ( parseTurner- ) where+module Biobase.Turner.Import 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 Control.Arrow+import Data.Array.Repa.Index+import Data.ByteString.Char8 as BS+import Data.ByteString.Lex.Double+import Data.Char+import Data.Conduit as C+import Data.Conduit.Binary as C+import Data.Conduit.List as CL+import Data.List.Split+import Data.Map as M+import Data.Maybe (fromJust)+import qualified Data.List as L+import System.FilePath.Posix+import qualified Data.Vector.Unboxed as VU -import Biobase.RNA-import Data.Ix.Tuple+import Biobase.Primary+import Biobase.Secondary import Data.PrimitiveArray-import Data.PrimitiveArray.Ix-import Text.Parsec.Numbers -import Biobase.Turner.Tables+import Biobase.Turner --- * 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 a directory, fill in the 'Turner2004' data structure +fromDir :: FilePath -> Prefix -> Suffix -> IO Turner2004+fromDir fp prefix suffix = do+ stack' <- blockFromFile $ fp </> prefix ++ "stack" <.> suffix+ dangle' <- blockFromFile $ fp </> prefix ++ "dangle" <.> suffix+ loop' <- blockFromFile $ fp </> prefix ++ "loop" <.> suffix+ hairpinMM' <- blockFromFile $ fp </> prefix ++ "tstackh" <.> suffix+ hairpinLk3 <- tabFromFile $ fp </> prefix ++ "triloop" <.> suffix+ hairpinLk4 <- tabFromFile $ fp </> prefix ++ "tloop" <.> suffix+ hairpinLk6 <- tabFromFile $ fp </> prefix ++ "hexaloop" <.> suffix+ let (dangle3',dangle5') = L.splitAt (L.length dangle' `div` 2) dangle'+ let (_:iloopL':bulgeL':hairpinL':[]) = L.transpose $ chunksOf 4 loop'+ iloop1x1' <- blockFromFile $ fp </> prefix ++ "int11" <.> suffix+ iloop2x1' <- blockFromFile $ fp </> prefix ++ "int21" <.> suffix+ iloop2x2' <- blockFromFile $ fp </> prefix ++ "int22" <.> suffix+ iloopMM' <- blockFromFile $ fp </> prefix ++ "tstacki" <.> suffix+ iloop2x3MM' <- blockFromFile $ fp </> prefix ++ "tstacki23" <.> suffix+ iloop1xnMM' <- blockFromFile $ fp </> prefix ++ "tstacki1n" <.> suffix+ multiMM' <- blockFromFile $ fp </> prefix ++ "tstackm" <.> suffix+ imisc' <- miscFromFile $ fp </> prefix ++ "miscloop" <.> suffix+ extMM' <- blockFromFile $ fp </> prefix ++ "tstack" <.> suffix+ coaxial' <- blockFromFile $ fp </> prefix ++ "coaxial" <.> suffix+ cstack' <- blockFromFile $ fp </> prefix ++ "coaxstack" <.> suffix+ tstack' <- blockFromFile $ fp </> prefix ++ "tstackcoax" <.> suffix+ return Turner2004Model+ { _stack = fromAssocs minPP maxPP infE $ L.zip keysPP stack'+ , _dangle3 = fromAssocs minPB maxPB infE $ L.zip keysPB dangle3'+ , _dangle5 = fromAssocs minPB maxPB infE $ L.zip keysPB dangle5'+ , _hairpinL = VU.fromList $ infE : hairpinL' -- fromAssocs (Z:.0) (Z:.30) infE $ L.zip d1_30 hairpinL'+ , _hairpinMM = fromAssocs minPBB maxPBB infE $ L.zip keysPBB hairpinMM'+ , _hairpinLookup = M.fromList $ hairpinLk3 ++ hairpinLk4 ++ hairpinLk6+ , _hairpinGGG = Energy . L.head $ imisc' !! 8+ , _hairpinCslope = Energy . L.head $ imisc' !! 9+ , _hairpinCintercept = Energy . L.head $ imisc' !! 10+ , _hairpinC3 = Energy . L.head $ imisc' !! 11+ , _bulgeL = VU.fromList $ infE : bulgeL' -- fromAssocs (Z:.0) (Z:.30) infE $ L.zip d1_30 bulgeL'+ , _bulgeSingleC = Energy . L.head $ imisc' !! 13+ , _iloop1x1 = fromAssocs minPPBB maxPPBB infE $ L.zip keysPPBB iloop1x1'+ , _iloop2x1 = fromAssocs minPPBBB maxPPBBB infE $ L.zip keysPPBBB iloop2x1'+ , _iloop2x2 = fromAssocs minPPBBBB maxPPBBBB infE $ L.zip (if (prefix == "" || suffix == "dh") then keysPPBBBBrna else keysPPBBBBdna) iloop2x2'+ , _iloopMM = fromAssocs minPBB maxPBB infE $ L.zip keysPBB iloopMM'+ , _iloop2x3MM = fromAssocs minPBB maxPBB infE $ L.zip keysPBB iloop2x3MM'+ , _iloop1xnMM = fromAssocs minPBB maxPBB infE $ L.zip keysPBB iloop1xnMM'+ , _iloopL = VU.fromList $ infE : iloopL' -- fromAssocs (Z:.0) (Z:.30) infE $ L.zip d1_30 iloopL'+ , _multiMM = fromAssocs minPBB maxPBB infE $ L.zip keysPBB multiMM'+ , _ninio = Energy . L.head $ imisc' !! 2+ , _maxNinio = Energy . L.head $ imisc' !! 1+ , _multiOffset = Energy $ (imisc' !! 3) !! 0+ , _multiNuc = Energy $ (imisc' !! 3) !! 1+ , _multiHelix = Energy $ (imisc' !! 3) !! 2+ , _multiAsym = Energy . L.head $ imisc' !! 5+ , _multiStrain = Energy . L.head $ imisc' !! 6+ , _extMM = fromAssocs minPBB maxPBB infE $ L.zip keysPBB extMM'+ , _coaxial = fromAssocs minPP maxPP infE $ L.zip keysPP coaxial'+ , _coaxStack = fromAssocs minPBB maxPBB infE $ L.zip keysPBB cstack'+ , _tStackCoax = fromAssocs minPBB maxPBB infE $ L.zip keysPBB tstack'+ , _largeLoop = Energy . L.head $ imisc' !! 0+ , _termAU = Energy . L.head $ imisc' !! 7+ , _intermolecularInit = Energy . L.head $ imisc' !! 12+ } +minPP = Z:.nN:.nN:.nN:.nN -- (minP,minP)+maxPP = Z:.nU:.nU:.nU:.nU -- (maxP,maxP)+minP = Z:.nN:.nN -- (nN,nN)+maxP = Z:.nU:.nU -- (nU,nU)+minPB = minP:.nN -- (minP,nN)+maxPB = maxP:.nU -- (maxP,nU)+minPBB = minPB:.nN -- (minP,nN,nN)+maxPBB = maxPB:.nU -- (maxP,nU,nU)+minPPBB = minPP:.nN:.nN -- (minP,minP,(nN,nN))+maxPPBB = maxPP:.nU:.nU -- (maxP,maxP,(nU,nU))+minPPBBB = minPPBB:.nN -- (minP,minP,(nN,nN,nN))+maxPPBBB = maxPPBB:.nU -- (maxP,maxP,(nU,nU,nU))+minPPBBBB = minPPBBB:.nN -- (minP,minP,(nN,nN,nN,nN))+maxPPBBBB = maxPPBBB:.nU -- (maxP,maxP,(nU,nU,nU,nU)) --- | Given the base dir and a suffix (.dat/.dh most likely), parse the relevant--- files. The prefix can be used for, eg., dna file+d1_30 = L.map (Z:.) [1..30] -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- }+keysPP = [{- ((k1,k2),(k4,k3)) -} Z:.k1:.k2:.k4:.k3 | k1 <- acgu, k3 <- acgu, k2 <- acgu, k4 <- acgu]+keysPB = [{- ((k1,k2),k3) -} Z:.k1:.k2:.k3 | k1 <- acgu, k2 <- acgu, k3 <- acgu]+keysPBB = [ {- ((k1,k2),k3,k4) -} Z:.k1:.k2:.k3:.k4+ | k1 <- acgu, k3 <- acgu, k2 <- acgu, k4 <- acgu]+keysPPBB = [ {- ((k1,k2),(k4,k3),(k5,k6)) -} Z:.k1:.k2:.k4:.k3:.k5:.k6+ | (k1,k2) <- plist11, k5 <- acgu, (k3,k4) <- plist11, k6 <- acgu]+keysPPBBB = [ {- ((k1,k2),(k4,k3),(k5,k6,k7)) -} Z:.k1:.k2:.k4:.k3:.k5:.k6:.k7+ | (k1,k2) <- plist11, k6 <- acgu, k5 <- acgu, (k3,k4) <- plist11, k7 <- acgu]+keysPPBBBBrna = [ {- ((k1,k2),(k4,k3),(k5,k6,k7,k8)) -} Z:.k1:.k2:.k4:.k3:.k5:.k6:.k7:.k8+ | (k1,k2) <- plist22rna, (k3,k4) <- plist22rna, k5 <- acgu, k8 <- acgu, k6 <- acgu, k7 <- acgu]+keysPPBBBBdna = [ {- ((k1,k2),(k4,k3),(k5,k6,k7,k8)) -} Z:.k1:.k2:.k4:.k3:.k5:.k6:.k7:.k8+ | (k1,k2) <- plist22dna, (k3,k4) <- plist22dna, k5 <- acgu, k8 <- acgu, k6 <- acgu, k7 <- acgu] -grabB basedir fsuffix fname = do- p <- parseFromFile pBlockFile (basedir ++ "/" ++ fname ++ fsuffix)- case p of- (Right ans) -> return ans- (Left err) -> error $ show err+plist11 = [(nA,nU),(nC,nG),(nG,nC),(nU,nA),(nG,nU),(nU,nG)]+plist22rna = [(nA,nU),(nC,nG),(nG,nC),(nG,nU),(nU,nA),(nU,nG)]+plist22dna = [(nA,nT),(nC,nG),(nG,nC),(nT,nA),(nG,nT),(nT,nG)] -grabA basedir fsuffix fname = do- p <- parseFromFile pAssocFile (basedir ++ "/" ++ fname ++ fsuffix)- case p of- (Right ans) -> return ans- (Left err) -> error $ show err+infE = Energy 999999 +-- * Conduit stuff +-- | extract values. "." - values are extracted as > 100k --- * File parser: We do not want to understand the format, just extract all--- data.+values :: ByteString -> [Energy]+values xs+ | BS.null ys+ = []+ | "." `isPrefixOf` ys+ = infE : values (BS.drop 1 ys)+ | Just (d,zs) <- readDouble ys+ = Energy d : values zs+ where ys = BS.dropWhile isSpace xs --- | Blocks of data.+-- | Iteratee to parse tabulated loops (hairpins). -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"+parseTabulated :: Monad m => Sink ByteString m [(ByteString,Energy)]+parseTabulated = C.lines =$ CL.filter (not . BS.all isSpace) =$ g where+ g = do+ CL.drop 2+ xs <- CL.map f =$ consume+ return xs+ f x+ | Just (d,_) <- readDouble v = (k,Energy d)+ | otherwise = error $ "tabulated: <" ++ BS.unpack x ++ ">"+ where (k,v) = second (BS.dropWhile isSpace) . BS.span (not . isSpace) . BS.dropWhile isSpace $ x --- | A File with associations.+-- | Convenience function -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 ' ')+blockFromFile :: FilePath -> IO [Energy]+blockFromFile fp = do+ xs <- runResourceT $ sourceFile fp $$ parseBlocks =$ consume+ if (allEq $ L.map L.length xs)+ then return $ L.concat xs+ else error $ "in file: " ++ fp ++ " we have unequal line lengths" --- a line with crap-infoRow = [] <$ anyChar `manyTill` newline+-- | Transform input stream into list of list of doubles --- a rather stupid float-dotFloat = (\s a b -> read $ s:'0':a:b) <$> (char '-' <|> pure ' ') <*> char '.' <*> many1 digit+parseBlocks :: Monad m => Conduit ByteString m [Energy]+parseBlocks = C.lines =$= CL.map f =$= CL.filter (not . L.null) where+ f x+ | "5'" `isPrefixOf` y = []+ | "3'" `isPrefixOf` y = []+ | "." `isPrefixOf` y = values y+ | Just (d,xs) <- readDouble y = values y+ | otherwise = [] -- error $ BS.unpack x+ where y = BS.dropWhile isSpace x -nuc = "ACGUT"+++-- | Parses the miscloop table+--+-- NOTE extra brownie points for miscloop.dat for providing data in a form that+-- does not conform to normal number encoding.++parseMiscLoop :: Monad m => Sink ByteString m [[Double]]+parseMiscLoop = C.lines =$ CL.groupBy (\x y -> not $ BS.null y) =$ f where+ f = do+ CL.drop 1+ xs <- consume+ return . L.map (L.map readD . BS.words . L.last) $ xs++-- | Parses stupidly encoded values like ".6" and "-.0".++readD :: ByteString -> Double+readD xs+ | BS.null xs = error "readD: null"+ | BS.head xs == '.' = readD $ BS.cons '0' xs+ | "-." `BS.isPrefixOf` xs = readD $ "-0." `BS.append` BS.drop 2 xs+ | Just (d,_) <- readDouble xs = d+ | otherwise = error $ BS.unpack xs++-- |++miscFromFile :: FilePath -> IO [[Double]]+miscFromFile fp = runResourceT $ sourceFile fp $$ parseMiscLoop++-- |++tabFromFile :: FilePath -> IO [(Primary,Energy)]+tabFromFile fp = fmap (L.map (first mkPrimary)) . runResourceT $ sourceFile fp $$ parseTabulated++allEq [] = True+allEq (x:xs) = L.all (==x) xs++type Prefix = FilePath+type Suffix = FilePath+
− 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)
BiobaseTurner.cabal view
@@ -1,39 +1,60 @@ name: BiobaseTurner-version: 0.0.2.2+version: 0.3.1.1 author: Christian Hoener zu Siederdissen maintainer: choener@tbi.univie.ac.at-copyright: Christian Hoener zu Siederdissen, 2010+homepage: http://www.tbi.univie.ac.at/~choener/+copyright: Christian Hoener zu Siederdissen, 2010-2013 category: Bioinformatics-synopsis: (deprecated) RNA folding data structures.+synopsis: Import Turner RNA parameters license: GPL-3 license-file: LICENSE build-type: Simple stability: experimental-cabal-version: >= 1.4.0+cabal-version: >= 1.6.0 description:- Data structures and parsers for the Turner 2004 set of RNA parameters.- - http://rna.urmc.rochester.edu/NNDB/index.html+ A data structure for Mathews / Turner RNA and DNA energy+ parameters. This library currently only provides an importer,+ not export functions. There are two reasons: (i) We currently+ have no use-case where we need more than import facilities (ii)+ The file structure is geared towards humans, not machines. If+ you need to be able to export, send a mail.+ .+ NOTE This is rather fragile as some files use different index+ enumerations, which we handle rather... simplistically.+ .+ In principle, all parameters should be symmetric regarding the+ stem direction. However, there is one asymmetry case in+ dnastack.dh. We do not fix this problem as we do not change the+ source files. library build-depends:- base >=4 && <5,- containers,- parsec >=3 && <4,- primitive >=0.3 && <0.4,- filemanip >=0.3.3 && <0.4,- split >=0.1.2 && <0.2,-- Biobase >=0.0.2 && <0.0.3,- HsTools >=0.0.1.1 && <0.0.2,- ParsecTools >=0.0.2 && <0.0.3,- PrimitiveArray >=0.0.2 && <0.0.3+ base >3 && <5,+ bytestring >= 0.9 ,+ bytestring-lexing >= 0.4 ,+ conduit >= 0.5 ,+ containers >= 0.4 ,+ filepath >= 1 ,+ lens >= 3.8 ,+ primitive >= 0.5 ,+ repa >= 3.2 ,+ split >= 0.2 ,+ vector >= 0.10 ,+ BiobaseXNA >= 0.7 ,+ PrimitiveArray >= 0.5 exposed-modules: Biobase.Turner Biobase.Turner.Import- Biobase.Turner.Tables ghc-options:- -O2+ -O2 -funbox-strict-fields++++source-repository head+ type: git+ location: git://github.com/choener/BiobaseTurner+