diff --git a/Biobase/Turner.hs b/Biobase/Turner.hs --- a/Biobase/Turner.hs +++ b/Biobase/Turner.hs @@ -1,3 +1,9 @@ +{-# LANGUAGE RecordWildCards #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE StandaloneDeriving #-} +{-# LANGUAGE GeneralizedNewtypeDeriving #-} +{-# LANGUAGE NoMonomorphismRestriction #-} +{-# LANGUAGE TemplateHaskell #-} {-# LANGUAGE TypeOperators #-} -- | The 'Turner2004' data structure reflects the RNA (and DNA) energy @@ -6,26 +12,85 @@ -- 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 Data.ByteString -import Data.Map as M + +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 -import Data.PrimitiveArray.Unboxed.Zero +import Data.PrimitiveArray as PA +import Data.PrimitiveArray.Zero --- | The parameters. Turner parameters are set by the Import module for --- nucleotides n,a,c,g,u. All values that are not read (or are ".") will end up --- with a value > 100K. --- --- TODO use 'Energy' instead of 'Double' --- --- TODO specialized shape types for pairs? + +-- | The actual Turner parameters return energies in Double format. + +newtype Energy = Energy Double + deriving (Eq,Ord,Num,Read,Show) + +deriving instance Prim Energy +deriving instance VGM.MVector VU.MVector Energy +deriving instance VG.Vector VU.Vector Energy +deriving instance VU.Unbox Energy + +-- | The Turner model with 'Energy's. + +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 ByteString 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) @@ -33,40 +98,45 @@ type PPNNN = PPNN:.Nuc type PPNNNN = PPNNN:.Nuc -data Turner2004 = Turner2004 - { stack :: Arr0 PP Double - , dangle3 :: Arr0 PN Double - , dangle5 :: Arr0 PN Double - , hairpinL :: Arr0 DIM1 Double - , hairpinMM :: Arr0 PNN Double - , hairpinLookup :: M.Map ByteString Double - , hairpinGGG :: Double - , hairpinCslope :: Double - , hairpinCintercept :: Double - , hairpinC3 :: Double - , bulgeL :: Arr0 DIM1 Double - , bulgeSingleC :: Double - , iloop1x1 :: Arr0 PPNN Double - , iloop2x1 :: Arr0 PPNNN Double - , iloop2x2 :: Arr0 PPNNNN Double - , iloopMM :: Arr0 PNN Double - , iloop2x3MM :: Arr0 PNN Double - , iloop1xnMM :: Arr0 PNN Double - , iloopL :: Arr0 DIM1 Double - , multiMM :: Arr0 PNN Double - , ninio :: Double - , maxNinio :: Double - , multiOffset :: Double - , multiNuc :: Double - , multiHelix :: Double - , multiAsym :: Double - , multiStrain :: Double - , extMM :: Arr0 PNN Double - , coaxial :: Arr0 PP Double -- no intervening unpaired nucleotides - , coaxStack :: Arr0 PNN Double - , tStackCoax :: Arr0 PNN Double - , largeLoop :: Double - , termAU :: Double - , intermolecularInit :: Double - } deriving () +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 @@ -41,14 +41,15 @@ import Data.ByteString.Char8 as BS import Data.ByteString.Lex.Double import Data.Char -import Data.Iteratee as I -import Data.Iteratee.Char as I -import Data.Iteratee.IO as I +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.Primary import Biobase.Secondary @@ -72,7 +73,7 @@ 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 $ splitEvery 4 loop' + 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 @@ -85,41 +86,41 @@ coaxial' <- blockFromFile $ fp prefix ++ "coaxial" <.> suffix cstack' <- blockFromFile $ fp prefix ++ "coaxstack" <.> suffix tstack' <- blockFromFile $ fp prefix ++ "tstackcoax" <.> suffix - return Turner2004 - { 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 = 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 = L.head $ imisc' !! 8 - , hairpinCslope = L.head $ imisc' !! 9 - , hairpinCintercept = L.head $ imisc' !! 10 - , hairpinC3 = L.head $ imisc' !! 11 - , bulgeL = fromAssocs (Z:.0) (Z:.30) infE $ L.zip d1_30 bulgeL' - , bulgeSingleC = 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 = fromAssocs (Z:.0) (Z:.30) infE $ L.zip d1_30 iloopL' - , multiMM = fromAssocs minPBB maxPBB infE $ L.zip keysPBB multiMM' - , ninio = L.head $ imisc' !! 2 - , maxNinio = L.head $ imisc' !! 1 - , multiOffset = (imisc' !! 3) !! 0 - , multiNuc = (imisc' !! 3) !! 1 - , multiHelix = (imisc' !! 3) !! 2 - , multiAsym = L.head $ imisc' !! 5 - , multiStrain = 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 = L.head $ imisc' !! 0 - , termAU = L.head $ imisc' !! 7 - , intermolecularInit = L.head $ imisc' !! 12 + 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) @@ -156,67 +157,69 @@ 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)] -infE = 999999 :: Double - --- * Iteratee stuff - --- | Transform input stream into list of list of doubles +infE = Energy 999999 -eneeBlocks :: (Functor m, Monad m) => Enumeratee ByteString [[Double]] m a -eneeBlocks = enumLinesBS ><> mapStream f ><> I.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 +-- * Conduit stuff -- | extract values. "." - values are extracted as > 100k -values :: ByteString -> [Double] +values :: ByteString -> [Energy] values xs | BS.null ys = [] | "." `isPrefixOf` ys = infE : values (BS.drop 1 ys) | Just (d,zs) <- readDouble ys - = d : values zs + = Energy d : values zs where ys = BS.dropWhile isSpace xs -- | Iteratee to parse tabulated loops (hairpins). -iTabulated :: (Functor m, Monad m) => Iteratee ByteString m [(ByteString,Double)] -iTabulated = joinI $ enumLinesBS ><> I.filter (BS.all isSpace) $ g where +parseTabulated :: Monad m => Sink ByteString m [(ByteString,Energy)] +parseTabulated = C.lines =$ CL.filter (not . BS.all isSpace) =$ g where g = do - I.drop 2 - joinI $ mapStream f stream2stream + CL.drop 2 + xs <- CL.map f =$ consume + return xs f x - | Just (d,_) <- readDouble v = (k,d) + | 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 -- | Convenience function -blockFromFile :: FilePath -> IO [Double] +blockFromFile :: FilePath -> IO [Energy] blockFromFile fp = do - i <- enumFile 8192 fp . joinI $ eneeBlocks stream2list - xs <- run i + 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" +-- | Transform input stream into list of list of doubles + +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 + + + -- | Parses the miscloop table -- --- NOTO extra brownie points for miscloop.dat for providing data in a form that +-- NOTE extra brownie points for miscloop.dat for providing data in a form that -- does not conform to normal number encoding. -iMiscLoop :: (Functor m, Monad m) => Iteratee ByteString m [[Double]] -iMiscLoop = joinI $ enumLinesBS ><> I.groupBy (\x y -> not $ BS.null y) $ f where +parseMiscLoop :: Monad m => Sink ByteString m [[Double]] +parseMiscLoop = C.lines =$ CL.groupBy (\x y -> not $ BS.null y) =$ f where f = do - I.drop 1 - xs <- fmap (L.map (L.map (readD) . BS.words . L.last)) $ stream2list - return xs + CL.drop 1 + xs <- consume + return . L.map (L.map readD . BS.words . L.last) $ xs -- | Parses stupidly encoded values like ".6" and "-.0". @@ -231,12 +234,12 @@ -- | miscFromFile :: FilePath -> IO [[Double]] -miscFromFile fp = run =<< enumFile 8192 fp iMiscLoop +miscFromFile fp = runResourceT $ sourceFile fp $$ parseMiscLoop -- | -tabFromFile :: FilePath -> IO [(ByteString,Double)] -tabFromFile fp = run =<< enumFile 8192 fp iTabulated +tabFromFile :: FilePath -> IO [(ByteString,Energy)] +tabFromFile fp = runResourceT $ sourceFile fp $$ parseTabulated allEq [] = True allEq (x:xs) = L.all (==x) xs diff --git a/BiobaseTurner.cabal b/BiobaseTurner.cabal --- a/BiobaseTurner.cabal +++ b/BiobaseTurner.cabal @@ -1,9 +1,9 @@ name: BiobaseTurner -version: 0.2.2.4 +version: 0.3.0.0 author: Christian Hoener zu Siederdissen maintainer: choener@tbi.univie.ac.at homepage: http://www.tbi.univie.ac.at/~choener/ -copyright: Christian Hoener zu Siederdissen, 2010-2012 +copyright: Christian Hoener zu Siederdissen, 2010-2013 category: Bioinformatics synopsis: Import Turner RNA parameters license: GPL-3 @@ -19,32 +19,31 @@ The file structure is geared towards humans, not machines. If you need to be able to export, send a mail. . - This is the 2011 post-library split version, hence not - deprecated anymore. - . NOTE This is rather fragile as some files use different index - enumerations, which we handle rather... simplistically. We - cannot fix one asymmetry case in dnastack.dh, as we do not - change sources. + enumerations, which we handle rather... simplistically. . - BIG FAT WARNING indexing now depends on repa shapes and index - representations. (with a good reason coming soon). For now, - just assume that this has performance benefits. + 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 >3 && <5, - bytestring >= 0.9, - bytestring-lexing >= 0.4, - containers >= 0.4, - filepath >= 1, - iteratee >= 0.8.8, - split >= 0.1.4, - vector == 0.9.* , - BiobaseXNA == 0.6.2.5 , - PrimitiveArray == 0.2.2.0 + 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