diff --git a/Biobase/Turner.hs b/Biobase/Turner.hs --- a/Biobase/Turner.hs +++ b/Biobase/Turner.hs @@ -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: +-- 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 + } + diff --git a/Biobase/Turner/Import.hs b/Biobase/Turner/Import.hs --- a/Biobase/Turner/Import.hs +++ b/Biobase/Turner/Import.hs @@ -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 + diff --git a/Biobase/Turner/Tables.hs b/Biobase/Turner/Tables.hs deleted file mode 100644 --- a/Biobase/Turner/Tables.hs +++ /dev/null @@ -1,133 +0,0 @@ -{-# LANGUAGE TemplateHaskell #-} -{-# LANGUAGE DeriveFunctor #-} -{-# LANGUAGE RecordWildCards #-} -{-# LANGUAGE TypeSynonymInstances #-} - --- | Turner Nearest Neighbor Energy Tables} - --- TODO (Read) instance? (Vector has no Read instance...) - -module Biobase.Turner.Tables where - -import qualified Data.Map as M - -import Data.PrimitiveArray -import Data.PrimitiveArray.Ix - - - --- | A data structure describing all fields as used by the Turner 2004 --- parameter set. Some fields are commented out until they are being used. - -data Turner2004 a b c = Turner2004 - { stack :: PrimArray (a,a) c - , dangle3 :: PrimArray (a,b) c - , dangle5 :: PrimArray (a,b) c - , hairpinL :: PrimArray Int c - , hairpinMM :: PrimArray (a,b,b) c - , hairpinLookup :: M.Map [b] c --- , hairpinGGG :: c --- , hairpinCslope :: c --- , hairpinCintercept :: c --- , hairpinC3 :: c - , bulgeL :: PrimArray Int c --- , bulgeSingleC :: c - , iloop1x1 :: PrimArray (a,a,b,b) c - , iloop1x2 :: PrimArray (a,a,b,b,b) c - , iloop2x2 :: PrimArray (a,a,b,b,b,b) c - , iloopMM :: PrimArray (a,b,b) c - , iloop2x3MM :: PrimArray (a,b,b) c - , iloop1xnMM :: PrimArray (a,b,b) c - , iloopL :: PrimArray Int c - , multiMM :: PrimArray (a,b,b) c - , ninio :: c - , maxNinio :: c - , multiOffset :: c - , multiNuc :: c - , multiHelix :: c --- , multiAsym :: c --- , multiStrain :: c - , extMM :: PrimArray (a,b,b) c --- , coaxStack :: PrimArray (a,a) c --- , coaxStackOpen :: PrimArray (a,b,b) c --- , coaxStackCont :: PrimArray (a,b,b) c - , largeLoop :: Double - , termAU :: c - , intermolecularInit :: c - } - - - --- | Map functions over the payload. - -dmap f Turner2004{..} = Turner2004 - { stack = amap f stack - , dangle3 = amap f dangle3 - , dangle5 = amap f dangle5 - , hairpinL = amap f hairpinL - , hairpinMM = amap f hairpinMM - , hairpinLookup = fmap f hairpinLookup --- , hairpinGGG = f hairpinGGG --- , hairpinCslope = f hairpinCslope --- , hairpinCintercept = f hairpinCintercept --- , hairpinC3 = f hairpinC3 - , bulgeL = amap f bulgeL --- , bulgeSingleC = f bulgeSingleC - , iloop1x1 = amap f iloop1x1 - , iloop1x2 = amap f iloop1x2 - , iloop2x2 = amap f iloop2x2 - , iloopMM = amap f iloopMM - , iloop2x3MM = amap f iloop2x3MM - , iloop1xnMM = amap f iloop1xnMM - , iloopL = amap f iloopL - , multiMM = amap f multiMM - , ninio = f ninio - , maxNinio = f maxNinio - , multiOffset = f multiOffset - , multiNuc = f multiNuc - , multiHelix = f multiHelix --- , multiAsym = f multiAsym --- , multiStrain = f multiStrain - , extMM = amap f extMM --- , coaxStack = amap f coaxStack --- , coaxStackOpen = amap f coaxStackOpen --- , coaxStackCont = amap f coaxStackCont - , largeLoop = largeLoop - , termAU = f termAU - , intermolecularInit = f intermolecularInit - } - - - --- | Zip two payloads. --- --- TODO right now, we have undefined behaviour when some arrrays are of different length - -dZipWith f t1 t2 = Turner2004 - { stack = zipWithPA f (stack t1) (stack t2) - , dangle3 = zipWithPA f (dangle3 t1) (dangle3 t2) - , dangle5 = zipWithPA f (dangle5 t1) (dangle5 t2) - , hairpinL = zipWithPA f (hairpinL t1) (hairpinL t2) - , hairpinMM = zipWithPA f (hairpinMM t1) (hairpinMM t2) - , hairpinLookup = undefined -- zipWithPA f (hairpinLookup t1) (hairpinLookup t2) - , bulgeL = zipWithPA f (bulgeL t1) (bulgeL t2) - , iloop1x1 = zipWithPA f (iloop1x1 t1) (iloop1x1 t2) - , iloop1x2 = zipWithPA f (iloop1x2 t1) (iloop1x2 t2) - , iloop2x2 = zipWithPA f (iloop2x2 t1) (iloop2x2 t2) - , iloopMM = zipWithPA f (iloopMM t1) (iloopMM t2) - , iloop2x3MM = zipWithPA f (iloop2x3MM t1) (iloop2x3MM t2) - , iloop1xnMM = zipWithPA f (iloop1xnMM t1) (iloop1xnMM t2) - , iloopL = zipWithPA f (iloopL t1) (iloopL t2) - , multiMM = zipWithPA f (multiMM t1) (multiMM t2) - , ninio = f (ninio t1) (ninio t2) - , maxNinio = f (maxNinio t1) (maxNinio t2) - , multiOffset = f (multiOffset t1) (multiOffset t2) - , multiNuc = f (multiNuc t1) (multiNuc t2) - , multiHelix = f (multiHelix t1) (multiHelix t2) - , extMM = zipWithPA f (extMM t1) (extMM t2) - , largeLoop = largeLoop t1 -- no zipping here! TODO make this a constant somewhere! - , termAU = f (termAU t1) (termAU t2) - , intermolecularInit = f (intermolecularInit t1) (intermolecularInit t2) - } - -zipWithPA f a1 a2 = let (l,u) = bounds a1 in fromList l u $ zipWith f (toList a1) (toList a2) diff --git a/BiobaseTurner.cabal b/BiobaseTurner.cabal --- a/BiobaseTurner.cabal +++ b/BiobaseTurner.cabal @@ -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 +