RNAFold 0.0.2.1 → 1.99.3.4
raw patch · 14 files changed
Files
- BioInf/RNAEval.hs +0/−122
- BioInf/RNAFold.hs +0/−115
- BioInf/RNAFold/Energy.hs +0/−92
- BioInf/RNAFold/EnergyInt.hs +0/−235
- BioInf/RNAFold/Functions.hs +0/−575
- BioInf/ViennaRNA.hs +2/−0
- BioInf/ViennaRNA/Energy.hs +117/−0
- BioInf/ViennaRNA/Eval.hs +95/−0
- BioInf/ViennaRNA/Fold.hs +192/−0
- BioInf/ViennaRNA/Signature.hs +38/−0
- README.md +39/−0
- RNAFold.cabal +85/−48
- src/RNAEval.hs +91/−0
- src/RNAFold.hs +37/−0
− BioInf/RNAEval.hs
@@ -1,122 +0,0 @@---- | Given a sequence and a structure, evaluate the energy.------ TODO switch to 'SSTree [Energy]' type!--module BioInf.RNAEval where--import Text.Printf-import qualified Data.Vector.Unboxed as VU--import Biobase.RNA-import Biobase.Structure-import Biobase.Structure.DotBracket-import Biobase.Types.Ring-import Biobase.Vienna-import Data.PrimitiveArray--import BioInf.RNAFold.EnergyInt-import BioInf.RNAFold.Functions------ | Sum up a complete (sub-) tree.--rnaEval :: ViennaTables -> String -> String -> Int-rnaEval vna pri' sec' = treeSum $ annotateWithEnergy vna pri sst where- sec = dotbracketToPairlist sec'- sst = toSSTree sec- pri = mkPrimary pri'------ | Evaluate the energy of a secondary structure tree with sequence. We abuse--- the normal folding functions with a dummy table full of (one :: Energy).--- This is probably slower than another method but quickly written.------ TODO this is basically crapfuck ;-) Should use the FoldFunctions directly--- instead of that strange table. Should not xxxOpt either.--annotateWithEnergy :: ViennaTables -> Primary -> SSTree () -> SSTree [Int]-annotateWithEnergy trnr pri t = f t where- -- extern part- f (SSExt l _ []) = SSExt l [one] []- f (SSExt l _ xs) =- let- es = map ((\(i,j) -> externalLoopOpt trnr pri dummy i j) . treeIJ) xs- in- SSExt l es (map f xs)- -- hairpin- f (SSTree i j _ []) = SSTree i j [hairpinOpt trnr pri i j] []- -- 1-loop- f (SSTree i j _ [x])- -- stack- | i+1==k&&j-1==l- = SSTree i j [stackOpt trnr pri dummy i j] [f x]- | (di,dj) `VU.elem` tabbedInteriorLoopDistances i j- = SSTree i j [tabbedInteriorLoopOpt trnr pri dummy i j] [f x]- | di==0&&dj>1- = SSTree i j [bulgeLOpt trnr pri dummy i j] [f x]- | di>1&&dj==0- = SSTree i j [bulgeROpt trnr pri dummy i j] [f x]- | di==1&&dj>1- = SSTree i j [interior1xnLOpt trnr pri dummy i j] [f x]- | di>1&&dj==1- = SSTree i j [interior1xnROpt trnr pri dummy i j] [f x]- | ds+dl>30- = error "not handling loops with size >30"- -- large interiorr loop- | otherwise = SSTree i j [largeInteriorLoopOpt trnr pri dummy i j] [f x]- where- (k,l) = treeIJ x- di = k-i-1- dj = j-l-1- ds = min di dj- dl = max di dj- -- multiple loop- f (SSTree i j _ xs) =- let- cl = multibranchCloseOpt trnr pri dummy dummy i j- ls = map ((\(k,l) -> multibranchIJLoopOpt trnr pri dummy k l) . treeIJ) xs- in- SSTree i j (cl:ls) (map f xs)-- dummy = fromAssocs (0,0) (n,n) one []- n = snd $ bounds pri------ | convert an annotated tree into strings that explain each score.------ TODO this is a stupid name--explainTree :: Primary -> SSTree [Int] -> [String]-explainTree inp sst = f sst where- -- | exterior loop- f (SSExt _ _ []) = [""]- f (SSExt _ _ xs) = this : concatMap f xs where- eners = map (\(SSTree _ _ [e] _) -> e) xs- ener = sum eners- this = printf "%-20s %6d %s" "External loop" ener (show eners)- -- | hairpin- f (SSTree i j [e] []) = [this] where- this = printf "%-20s (%4d,%4d) %4s %6d" "Hairpin loop" (i+1) (j+1) (show $ pair inp i j) e- f (SSTree i j [e] [x@(SSTree k l _ _)]) = this : f x where- this = printf "%-20s (%4d,%4d) %4s (%4d,%4d) %4s %6d" "Interior" (i+1) (j+1) (show $ pair inp i j) (k+1) (l+1) (show $ pair inp k l) e- f (SSTree i j es xs) = concatMap f xs ++ [this] where- this = printf "%-20s (%4d,%4d) %4s %6d %6d, %s" "Multi loop" (i+1) (j+1) (show $ pair inp i j) (sum es) (head es) (show $ tail es)------ | input sequence indices--treeIJ (SSTree i j _ _) = (i,j)----- | Run a sum over the tree. (foldl (+), the tree annotations)------ TODO that should go into Biobase.Structure as a generic walk over the tree--treeSum t = f t where- f (SSExt _ es xs) = ringProductL $ es ++ map f xs- f (SSTree _ _ es xs) = ringProductL $ es ++ map f xs
− BioInf/RNAFold.hs
@@ -1,115 +0,0 @@---- | ViennaRNA folding based on an algebraic ring structure. This should--- combine the goals of few lines of codes, multiple different folding--- functions and extensibility.------ NOTE Assume that you want '-d 3' for folding with dangles. Then you can just--- instanciate the folding functions, replacing only those functions where the--- folding changes based on the new dangle options.------ NOTE compile with: -fno-method-sharing--module BioInf.RNAFold where--import Control.Monad-import Control.Monad.ST--import Biobase.RNA-import Biobase.Types.Ring-import Data.PrimitiveArray-import Biobase.Structure--import BioInf.RNAFold.Functions--import Debug.Trace.Tools-import Debug.Trace----- | Folding works on unboxed values of a Ring-type for which a FoldFunctions--- instance does exist. By default, we have this for Energy values. Again, we--- use a class as we could be interested in probabilistic backtracking or--- something like that.--type ResultTables a =- ( Table a -- weak structures- , Table a -- strong structures- , Table a -- exactly one component- , Table a -- one or more components- , Table a -- complete external structures- )--type Pairlist = [(Int,Int)]--class (FoldFunctions a) => Fold a where-- fold :: TurnerTables a -> Primary -> (ResultTables a)- foldST :: TurnerTables a -> Primary -> ST s (ResultTables a)- backtrack :: TurnerTables a -> Primary -> (ResultTables a) -> a -> [(Secondary,a)]-- -- | We have a default instance for folding based on Rings-- fold trnr inp = runST $ foldST trnr inp- {-# INLINE fold #-}-- foldST trnr inp = do- let n = snd $ bounds inp- (weakM,weak) <- mkTable n- (strongM,strong) <- mkTable n- (externM,extern) <- mkTableWith one n- (mbr1M,mbr1) <- mkTable n- (mbrM,mbr) <- mkTable n- forM_ [n,n-1..0] $ \i -> forM_ [i,i+1..n] $ \j -> do- let pIJ = pair inp i j- when (pIJ/=vpNP&&i+3<j) $ do- -- weak table- let hpVal = {-# SCC "hpVal" #-} hairpinOpt trnr inp i j- let ilVal = {-# SCC "ilVal" #-} ringSumL- [ largeInteriorLoopOpt trnr inp strong i j- , tabbedInteriorLoopOpt trnr inp strong i j- , bulgeLOpt trnr inp strong i j- , bulgeROpt trnr inp strong i j- , interior1xnLOpt trnr inp strong i j- , interior1xnROpt trnr inp strong i j- ]- let mbVal = {-# SCC "mbVal" #-} multibranchCloseOpt trnr inp mbr mbr1 i j- writeM weakM (i,j) $ ringSumL [hpVal,ilVal,mbVal]- -- strong table- when (i+5<j) $ do- let stValW = {-# SCC "stValW" #-} stackOpt trnr inp weak i j- let stValS = {-# SCC "stValS" #-} stackOpt trnr inp strong i j- writeM strongM (i,j) $ ringSumL [stValW,stValS]- -- multibranch loops- when (i>0&&j<n) $ do- -- M1- let mbr1ValS = {-# SCC "mbr1ValS" #-} multibranchIJLoopOpt trnr inp strong i j- let mbr1ValU = {-# SCC "mbr1ValU" #-} multibranchUnpairedJOpt trnr inp mbr1 i j- writeM mbr1M (i,j) $ ringSumL [mbr1ValS,mbr1ValU]- -- M- let mbrValU = {-# SCC "mbrValU" #-} multibranchUnpairedJOpt trnr inp mbr i j- let mbrValS = {-# SCC "mbrValS" #-} multibranchKJHelixOpt trnr inp strong i j- let mbrValMS = {-# SCC "mbrValMS" #-} multibranchAddKJHelixOpt trnr inp mbr strong i j- writeM mbrM (i,j) $ ringSumL [mbrValU,mbrValS,mbrValMS]- -- fill only part of the F array- let j=n- forM_ [n-6,n-7..0] $ \i -> do- let extUP = {-# SCC "extUP" #-} if i<n then extern ! (i+1,j) else zero- let extStr = {-# SCC "extStr" #-} externalLoopOpt trnr inp strong i j- let extAddL = {-# SCC "extAddL" #-} externalAddLoopOpt trnr inp strong extern i j- writeM externM (i,j) $ ringSumL [extUP,extStr,extAddL,one] -- always add 'one' as the open chain should always be contained- return (weak,strong,mbr1,mbr,extern)- {-# INLINE foldST #-}----- * Helper functions---- Create one 2dim - IxTable with default value.--mkTable n = mkTableWith zero n---- Create one IxTable with user-supplied value.--mkTableWith v n = do- tM <- fromAssocsM (0,0) (n,n) v []- t <- unsafeFreezeM tM- return (tM,t)-
− BioInf/RNAFold/Energy.hs
@@ -1,92 +0,0 @@-{-# LANGUAGE StandaloneDeriving #-}--module BioInf.RNAFold.Energy- ( FoldFunctions (..)- , Fold (..)- ) where--import BioInf.RNAFold-import Biobase.Types.Energy-import Biobase.Types.Ring-import BioInf.RNAFold.Functions----instance FoldFunctions Energy--instance Fold Energy where- backtrack trnr inp tbls = error "write me"-{-- backtrack trnr inp (weak,strong,mbr1,mbr,extern) = ext 0 n delta where- n = VU.length inp -1- delta = one :: Energy- overallBest = extern `unsafeIndex` (0,n)- ext i j d = let bestE = extern `unsafeIndex` (i,j) in- [ []- | i==j- , overallBest == one- ] ++ -- gives us the unfolded sequence- [ r- | i<j- , r <- ext (i+1) j d- ] ++- [ r- | i<j- , ((k,l),ce) <- externalLoopIdx trnr inp strong i j- , let dNew = ce `rmult` d `rmult` neg bestE- , dNew <= one- , r <- str k l d- ] ++- [ r++s- | i<j- , (k,ce) <- externalAddLoopIdx trnr inp strong extern i j- , let dNew = ce `rmult` d `rmult` neg bestE- , dNew <= one- , r <- str i k d- , s <- ext (k+1) j d ++ [[]] -- not 'd' but what 'str' leaves us with!, the [[]] only if the str-part alone works out- ]- str i j d = let bestE = strong `unsafeIndex` (i,j) in- [ (i,j):r- | i+2<j- , ((k,l),ce) <- stackIdx trnr inp strong i j- , let dNew = ce `rmult` d `rmult` neg bestE- , dNew <= one- , r <- str k l d- ] ++- [ (i,j):r- | i+2<j- , ((k,l),ce) <- stackIdx trnr inp weak i j- , let dNew = ce `rmult` d `rmult` neg bestE- , dNew <= one- , r <- wea k l d- ]- wea :: Int -> Int -> a -> [Pairlist]- wea i j d = let bestE = weak `unsafeIndex` (i,j) in- [ [(i,j)]- | i+3<j- , (ce :: Energy) <- hairpinIdx trnr inp i j -- needs to be qualified as we give no other table- ]--}------ | (DEBUG) Print an array.--{--printArr arr = do- --let (IT.IxTable (0,0) (_,n) _) = arr- let n = snd $ snd $ bounds arr- forM_ [0..n] $ \i -> do- putStrLn ""- forM [0..n] $ \j -> do- if j>=i- then do- let v = arr `unsafeIndex` (i,j)- if isZero v- then printf "%5s" "_i_"- else printf "%5i" (unEnergy $ arr `unsafeIndex` (i,j))- else do- putStr " "- putStrLn ""- return ()--}
− BioInf/RNAFold/EnergyInt.hs
@@ -1,235 +0,0 @@-{-# LANGUAGE StandaloneDeriving #-}---- | Temporary hackery until all base libraries understand newtype Energy. And--- yes, for testing too.--module BioInf.RNAFold.EnergyInt- ( FoldFunctions (..)- , Fold (..)- ) where--import Biobase.Types.Ring-import Biobase.RNA-import Biobase.Constants-import Biobase.Structure.DotBracket-import Biobase.Structure--import Data.PrimitiveArray--import BioInf.RNAFold-import BioInf.RNAFold.Functions---- for testing only-{--import Biobase.Vienna.Default-import Debug.Trace.Tools-import Data.List.Split-import Text.Printf-import Control.Monad-import Biobase.Turner.Tables--}---instance Ring Int where- (.+.) = min- (.*.) = (+)- neg = negate- one = 0- zero = eInf- isZero x = x>eMax- n .^. k = n * k- (.^^.) = error "write me"- {-# INLINE (.+.) #-}- {-# INLINE (.*.) #-}- {-# INLINE neg #-}- {-# INLINE one #-}- {-# INLINE zero #-}- {-# INLINE isZero #-}- {-# INLINE (.^.) #-}--instance FoldFunctions Int where- calcTermAU tAU p = if p/=vpCG && p/=vpGC then tAU else 0- calcNinio maxNno nno l = min maxNno (nno * l)- calcLargeLoop l = floor $ 108.856 * log (fromIntegral l / 30)- {-# INLINE calcTermAU #-}- {-# INLINE calcNinio #-}- {-# INLINE calcLargeLoop #-}---- TODO detach backtrack so that all instances that use this kind of--- backtracking can immediately use it!--instance Fold Int where- backtrack trnr inp (weak,strong,mbr1,mbr,extern) delta = filter ((<=0) . snd) $ map f $ ext 0 n delta where- f (xs,z) = (Secondary (n+1) xs,optE+delta-z)- n = snd $ bounds inp- optE = extern!(0,n)- newD d here next = d - (next - here)- --- -- All exterior loop decompositions- --- ext i j d = let ehere = extern!(i,j) in- [ (x,z) -- shorter ext- | i<j- , let d' = newD d ehere (extern!(i+1,j))- , d'>=0- , (x,z) <- ext (i+1) j d'- ] ++- [ (x,z) -- loop at (i,j)- | i<j- , (_,enext) <- externalLoopIdx trnr inp strong i j -- _ = (i,j)- , let d' = newD d ehere enext- , d'>=0- , (x,z) <- str i j d'- ] ++- [ (x++y,z)- | i<j- , (k,enext) <- externalAddLoopIdx trnr inp strong extern i j- , let d' = newD d ehere enext- , d'>=0- , (x,z') <- str i k d'- , (y,z) <- ext (k+1) j z' ++ [([],z') | z'>=0 && extern!(k+1,j)==0]- ]- --- -- A strong stem on top of strong of weak- --- str i j d = let ehere = strong!(i,j) in- [ ((i,j):x,z)- | i<j- , (_,enext) <- stackIdx trnr inp strong i j- , let d' = newD d ehere enext- , d'>=0- , (x,z) <- str (i+1) (j-1) d'- ] ++- [ ((i,j):x,z)- | i<j- , (_,enext) <- stackIdx trnr inp weak i j- , let d' = newD d ehere enext- , d'>=0- , (x,z) <- wea (i+1) (j-1) d'- ]- wea i j d = let ehere = weak!(i,j) in- --- -- A hairpin closes at (i,j)- --- [ ([(i,j)],d')- | i<j- , enext <- hairpinIdx trnr inp i j- , let d' = newD d ehere enext- , d'>=0- ] ++- --- -- interior loops- --- [ ((i,j):x,z)- | i<j- , ((k,l),enext) <- (- largeInteriorLoopIdx trnr inp strong i j ++- tabbedInteriorLoopIdx trnr inp strong i j ++- bulgeLIdx trnr inp strong i j ++- bulgeRIdx trnr inp strong i j ++- interior1xnLIdx trnr inp strong i j ++- interior1xnRIdx trnr inp strong i j- )- , let d' = newD d ehere enext- , d'>=0- , (x,z) <- str k l d'- ] ++- --- -- multibranch loops closed at (i,j)- --- [ ((i,j):x++y,z)- | i<j- , (k,enext) <- multibranchCloseIdx trnr inp mbr mbr1 i j- , let d' = newD d ehere enext- , d'>=0- , (x,z') <- mul (i+1) k d'- , (y,z) <- mul1 (k+1) (j-1) z' -- z' is what 'mul' left us- ]- mul i j d = let ehere = mbr!(i,j) in- --- -- unpaired nucleotide on the J side- --- [ (x,z)- | i<j- , ((k,l),enext) <- multibranchUnpairedJIdx trnr inp mbr i j- , let d' = newD d ehere enext- , d'>=0- , (x,z) <- mul k l d'- ] ++- --- -- a helix (k,j)- --- [ (x,z)- | i<j- , (k,enext) <- multibranchKJHelixIdx trnr inp strong i j- , let d' = newD d ehere enext- , d'>=0- , (x,z) <- str k j d'- ] ++- --- -- add a helix to an existing structure- --- [ (x++y,z)- | i<j- , (k,enext) <- multibranchAddKJHelixIdx trnr inp mbr strong i j- , let d' = newD d ehere enext- , d'>=0- , (x,z') <- mul i k d'- , (y,z) <- str (k+1) j z'- ]- mul1 i j d = let ehere = mbr1!(i,j) in- --- -- Simply the strong part in a multibranch loop- --- [ (x,z)- | i<j- , (_,enext) <- multibranchIJLoopIdx trnr inp strong i j- , let d' = newD d ehere enext- , d'>=0- , (x,z) <- str i j d'- ] ++- --- -- unpaired nucleotide on the J side- --- [ (x,z)- | i<j- , ((k,l),enext) <- multibranchUnpairedJIdx trnr inp mbr1 i j- , let d' = newD d ehere enext- , d'>=0- , (x,z) <- mul1 k l d'- ]--{--test inp' k =- let- (_,n) = bounds inp- bts = backtrack trnr inp tbls k- inp = mkPrimary inp'- trnr = fst turner2004GH- tbls@(weak,strong,mbr1,mbr,extern) = fold trnr inp- optE = extern ! (0,n)- f x- | x > 999 = " x"- | x < -999 = "-999"- | otherwise = printf "%4d" x- g t = do- let ls = splitEvery (n+1) $ map f $ toList t- putStr " "- mapM_ (printf "%4d") [0 :: Int .. (length $ head ls) -1]- putStrLn ""- zipWithM_ (\k xs -> printf "%4d" k >> mapM_ putStr xs >> putStrLn"") [0 :: Int ..] ls- in do- print "weak"- g weak- print "strong"- g strong- print "mbr L"- g mbr- print "mbr1 R"- g mbr1- print "extern"- g extern- print optE- putStrLn inp'- mapM_ (\(s,e) -> (putStr $ dotbracket s) >> (putStrLn $ " " ++ show e)) bts--}
− BioInf/RNAFold/Functions.hs
@@ -1,575 +0,0 @@-{-# LANGUAGE PatternGuards #-}-{-# LANGUAGE BangPatterns #-}-{-# LANGUAGE RecordWildCards #-}---- | These functions are implementations of RNA secondary structure folding as--- described in "Bompfuenewere et al., 2006, Variations on RNA folding and--- alignment"------ We have all the facilities needed for folding with the RNA parameter of--- Turner 2004 http://rna.urmc.rochester.edu/NNDB/turner04/ but consider only--- "double dangles" which correspond to the ViennaRNA package option "-d2".--- They are a bit easier to implement and are what is used for partition--- function calculations. In addition, it seems unlikely to see a statiscally--- relevant improvement in predication with "-d1" or "-d3".------ All functions work on an algebraic ring structure. This should make it--- easier to implement certain functionality without having to rewrite all the--- functions given here. Try deriving a new 'Ring' instance first and see if it--- /just works/.------ These functions do quite well, performancewise. GHC-HEAD with "-Odph" and--- "-fllvm" takes 14.4s, while the highly optimized viennaRNA package (the yet--- unpublished 2.0 version) takes about 1-2s on a sequence of length 1000.------ NOTE For GHC <= 6.12.3 you should copy the default instances into your--- instance BLA, otherwise the resulting code will be slow. Or you could just--- wait for the new GHC to arrive! The new one produces good code without such--- stuff.------ TODO single nucleotide bulges:--- http://rna.urmc.rochester.edu/NNDB/turner04/bulge.html , check what Vienna--- 2.0 does!--module BioInf.RNAFold.Functions- ( FoldFunctions (..)- , Table- , TurnerTables- , ringSumL- , pair- , riap- , ringProductL- , tabbedInteriorLoopDistances- ) where--import qualified Data.Vector.Unboxed as VU-import Control.Exception (assert)-import Data.List (foldl')-import qualified Data.Map as M--import Biobase.RNA-import Biobase.Turner.Tables-import Biobase.Types.Ring-import Biobase.Vienna-import Data.PrimitiveArray-import Data.Primitive.Types--import Debug.Trace.Tools--type Cell = (Int,Int)-type Table a = PrimArray Cell a-type TurnerTables a = Turner2004 ViennaPair Nucleotide a------ | The folding functions. It could happen that we need different folding--- functions with the same type, hence the class-based approach. The default--- instance uses the usual ring methods.--class (Show a, Ring a, VU.Unbox a, Prim a) => FoldFunctions a where-- stackOpt :: TurnerTables a -> Primary -> Table a -> Int -> Int -> a- stackIdx :: TurnerTables a -> Primary -> Table a -> Int -> Int -> [(Cell,a)]-- hairpinOpt :: TurnerTables a -> Primary -> Int -> Int -> a- hairpinIdx :: TurnerTables a -> Primary -> Int -> Int -> [a]-- largeInteriorLoopOpt :: TurnerTables a -> Primary -> Table a -> Int -> Int -> a- largeInteriorLoopIdx :: TurnerTables a -> Primary -> Table a -> Int -> Int -> [(Cell,a)]-- tabbedInteriorLoopOpt :: TurnerTables a -> Primary -> Table a -> Int -> Int -> a- tabbedInteriorLoopIdx :: TurnerTables a -> Primary -> Table a -> Int -> Int -> [(Cell,a)]-- bulgeLOpt :: TurnerTables a -> Primary -> Table a -> Int -> Int -> a- bulgeLIdx :: TurnerTables a -> Primary -> Table a -> Int -> Int -> [(Cell,a)]-- bulgeROpt :: TurnerTables a -> Primary -> Table a -> Int -> Int -> a- bulgeRIdx :: TurnerTables a -> Primary -> Table a -> Int -> Int -> [(Cell,a)]-- interior1xnLOpt :: TurnerTables a -> Primary -> Table a -> Int -> Int -> a- interior1xnLIdx :: TurnerTables a -> Primary -> Table a -> Int -> Int -> [(Cell,a)]-- interior1xnROpt :: TurnerTables a -> Primary -> Table a -> Int -> Int -> a- interior1xnRIdx :: TurnerTables a -> Primary -> Table a -> Int -> Int -> [(Cell,a)]-- multibranchIJLoopOpt :: TurnerTables a -> Primary -> Table a -> Int -> Int -> a- multibranchIJLoopIdx :: TurnerTables a -> Primary -> Table a -> Int -> Int -> [(Cell,a)]-- multibranchUnpairedJOpt :: TurnerTables a -> Primary -> Table a -> Int -> Int -> a- multibranchUnpairedJIdx :: TurnerTables a -> Primary -> Table a -> Int -> Int -> [(Cell,a)]-- multibranchKJHelixOpt :: TurnerTables a -> Primary -> Table a -> Int -> Int -> a- multibranchKJHelixIdx :: TurnerTables a -> Primary -> Table a -> Int -> Int -> [(Int,a)]-- multibranchAddKJHelixOpt :: TurnerTables a -> Primary -> Table a -> Table a -> Int -> Int -> a- multibranchAddKJHelixIdx :: TurnerTables a -> Primary -> Table a -> Table a -> Int -> Int -> [(Int,a)]-- multibranchCloseOpt :: TurnerTables a -> Primary -> Table a -> Table a -> Int -> Int -> a- multibranchCloseIdx :: TurnerTables a -> Primary -> Table a -> Table a -> Int -> Int -> [(Int,a)]-- externalLoopOpt :: TurnerTables a -> Primary -> Table a -> Int -> Int -> a- externalLoopIdx :: TurnerTables a -> Primary -> Table a -> Int -> Int -> [(Cell,a)]-- externalAddLoopOpt :: TurnerTables a -> Primary -> Table a -> Table a -> Int -> Int -> a- externalAddLoopIdx :: TurnerTables a -> Primary -> Table a -> Table a -> Int -> Int -> [(Int,a)] -- return only 'k' index-- -- | Calculate the ninio asymmetric malus. Can not be written using ring- -- functions alone as a 'min' or 'max' functions is required.-- calcNinio :: a -> a -> Int -> a-- -- | Applies a terminal AU/GU penalty, where required.- --- -- TODO shouldn't this be just: if CG||GC then one else termAU?-- calcTermAU :: a -> ViennaPair -> a -- ^ Apply terminal AU penalty-- -- | large hairpin loops >30 require special calculations that involve- -- 'floor', rounding and other stuff that can not be handled by the Ring- -- class alone-- calcLargeLoop :: Int -> a-- -- NOTE Copy these functions into your own instance. In most cases, your are- -- now done and get optimized loops. Each function that you do not copy uses- -- the version here. This can lead to less efficient code.- --- -- NOTE Fixed in GHC 6.13. The default instances should yield superb code!-- stackOpt trnr inp tbl i j =- VU.foldl' (.+.) zero $ stackBase trnr inp tbl i j- hairpinOpt trnr inp i j =- VU.foldl' (.+.) zero $ hairpinBase trnr inp i j- largeInteriorLoopOpt trnr inp strong i j =- VU.foldl' (.+.) zero $ largeInteriorLoopBase trnr inp strong i j- tabbedInteriorLoopOpt trnr inp strong i j =- VU.foldl' (.+.) zero $ tabbedInteriorLoopBase trnr inp strong i j- bulgeLOpt trnr inp strong i j =- VU.foldl' (.+.) zero $ bulgeLBase trnr inp strong i j- bulgeROpt trnr inp strong i j =- VU.foldl' (.+.) zero $ bulgeRBase trnr inp strong i j- interior1xnLOpt trnr inp strong i j =- VU.foldl' (.+.) zero $ interior1xnLBase trnr inp strong i j- interior1xnROpt trnr inp strong i j =- VU.foldl' (.+.) zero $ interior1xnRBase trnr inp strong i j- multibranchIJLoopOpt trnr inp strong i j =- VU.foldl' (.+.) zero $ multibranchIJLoopBase trnr inp strong i j- multibranchUnpairedJOpt trnr inp mtable i j =- VU.foldl' (.+.) zero $ multibranchUnpairedJBase trnr inp mtable i j- multibranchKJHelixOpt trnr inp strong i j =- VU.foldl' (.+.) zero $ multibranchKJHelixBase trnr inp strong i j- multibranchAddKJHelixOpt trnr inp table strong i j =- VU.foldl' (.+.) zero $ multibranchAddKJHelixBase trnr inp table strong i j- multibranchCloseOpt trnr inp m m1 i j =- VU.foldl' (.+.) zero $ multibranchCloseBase trnr inp m m1 i j- externalLoopOpt trnr inp strong i j =- VU.foldl' (.+.) zero $ externalLoopBase trnr inp strong i j- externalAddLoopOpt trnr inp strong extern i j =- VU.foldl' (.+.) zero $ externalAddLoopBase trnr inp strong extern i j-- -- NOTE copy and instanciate these functions if you have to work with many- -- candidate sequences. Otherwise you probably do not need the speed-up from- -- these functions. This stuff is for backtracking mostly as you get lists- -- of indices with attached scores.- --- -- NOTE With GHC 6.13 we should get optimized code anyways!-- stackIdx trnr inp tbl i j =- [((i+1,j-1),stackOpt trnr inp tbl i j)]- hairpinIdx trnr inp i j =- VU.toList $ hairpinBase trnr inp i j- largeInteriorLoopIdx trnr inp strong i j =- VU.toList $ VU.zip (interiorLoopIndices i j) (largeInteriorLoopBase trnr inp strong i j)- tabbedInteriorLoopIdx trnr inp strong i j =- VU.toList $ VU.zip (tabbedInteriorLoopIndices i j) $ tabbedInteriorLoopBase trnr inp strong i j- bulgeLIdx trnr inp strong i j =- VU.toList $ VU.zip (VU.map (\k -> (i+1+k,j-1)) . uncurry VU.enumFromN $ bulgeLimit i j) $ bulgeLBase trnr inp strong i j- bulgeRIdx trnr inp strong i j =- VU.toList $ VU.zip (VU.map (\k -> (i+1,j-1-k)) . uncurry VU.enumFromN $ bulgeLimit i j) $ bulgeRBase trnr inp strong i j- interior1xnLIdx trnr inp strong i j =- VU.toList $ VU.zip (VU.map (\k -> (i+1+k,j-2)) $ uncurry VU.enumFromN $ iloop1xnLimit i j) $ interior1xnLBase trnr inp strong i j- interior1xnRIdx trnr inp strong i j =- VU.toList $ VU.zip (VU.map (\k -> (i+2,j-1-k)) $ uncurry VU.enumFromN $ iloop1xnLimit i j) $ interior1xnRBase trnr inp strong i j- multibranchIJLoopIdx trnr inp strong i j =- VU.toList $ VU.zip (VU.singleton (i,j)) $ multibranchIJLoopBase trnr inp strong i j- multibranchUnpairedJIdx trnr inp mtable i j =- VU.toList $ VU.zip (VU.singleton (i,j-1)) $ multibranchUnpairedJBase trnr inp mtable i j- multibranchKJHelixIdx trnr inp strong i j =- VU.toList $ VU.zip (uncurry VU.enumFromN $ multibranchKJHelixLimit i j) $ multibranchKJHelixBase trnr inp strong i j- multibranchCloseIdx trnr inp m m1 i j = -- TODO this is wrong!- VU.toList $ VU.zip (uncurry VU.enumFromN $ multibranchCloseLimit i j) $ multibranchCloseBase trnr inp m m1 i j- multibranchAddKJHelixIdx trnr inp table strong i j =- VU.toList $ VU.zip (uncurry VU.enumFromN $ multibranchAddKJHelixLimit i j) $ multibranchAddKJHelixBase trnr inp table strong i j- externalLoopIdx trnr inp strong i j =- VU.toList $ VU.singleton ((i,j), externalLoopOpt trnr inp strong i j)- externalAddLoopIdx trnr inp strong extern i j =- VU.toList $ VU.zip (uncurry VU.enumFromN $ externalAddLoopLimit i j) $ externalAddLoopBase trnr inp strong extern i j------ * We provide a default instance working on rings.---- | Simple hairpin loops. If (i,j) pairs then, first, try to do a lookup of--- the exact sequence of the hairpin in a tabulated list. If this fails then,--- second, try length 3 hairpins and so on.--hairpinBase :: (FoldFunctions a) => TurnerTables a -> Primary -> Int -> Int -> VU.Vector a-hairpinBase Turner2004{..} inp i j = VU.singleton go where- go- | pIJ==vpNP || l<3 = zero- | l<=6, Just v <- s `M.lookup` hairpinLookup = v- | l==3 = (hairpinL ! l) -- .*. tAU -- apparantly not anymore according to NNDB- | l>=31 = (hairpinL ! 30) .*. (hairpinMM ! (pIJ,bI,bJ)) .*. (calcLargeLoop l)- | otherwise = {-traceVal (show (i,j,pIJ,bI,bJ,hairpinL!l,tAU,hairpinMM!(pIJ,bI,bJ))) $ -} (hairpinL ! l) .*. (hairpinMM ! (pIJ,bI,bJ))- l = j-i-1- s = assert (i>=0 && checkBounds inp j) $ [inp ! k | k <- [i..j-1]]- bI = inp ! (i+1)- bJ = inp ! (j-1)- pIJ = pair inp i j- tAU = calcTermAU termAU pIJ-{-# INLINE hairpinBase #-}---- | Extend a stem by another pair.--stackBase :: (FoldFunctions a) => TurnerTables a -> Primary -> Table a -> Int -> Int -> VU.Vector a-stackBase Turner2004{..} inp tbl i j = VU.singleton $ go (i+1,j-1) where- pIJ = pair inp i j- go (k,l)- | pIJ==vpNP || pKL==vpNP || isZero tE- = zero- | otherwise- = tE .*. ( stack ! (pIJ,pKL))- where- pKL = riap inp k l- tE = tbl ! (k,l)-{-# INLINE stackBase #-}---- | These are special cases of interior loops. Short iloops, such as 1x2 loops--- are completely tabulated. Look each one up here.------ TODO needs 1-bulges as a special case--tabbedInteriorLoopBase :: (Show a, FoldFunctions a) => TurnerTables a -> Primary -> Table a -> Int -> Int -> VU.Vector a-tabbedInteriorLoopBase Turner2004{..} inp strong i j = res where- res = VU.map (\(k,l) -> (strong ! (k,l)) .*. (ftabbed (k-i-1,j-l-1) (k,l))) ili- ili = tabbedInteriorLoopIndices i j- ftabbed (di,dj) (k,l)- | ds==0 && dl==1 = bulgeL!1 .*. stack!(pIJ,pKL) -- no termAU, since the helical stack continues (nndb)- | ds==1 && dl==1 = iloop1x1 ! (pIJ,pKL,bI,bJ)- | di==1 && dj==2 = iloop1x2 ! (pIJ,pKL,bI,bL,bJ)- | di==2 && dj==1 = iloop1x2 ! (pIJ,pKL,bJ,bI,bK)- | ds==2 && dl==2 = iloop2x2 ! (pIJ,pKL,bI,bK,bL,bJ)- | ds==2 && dl==3 = iloop2x3MM!(pIJ,bI,bJ) .*. iloop2x3MM ! (pKL,bL,bK) .*. iloopL ! 5 .*. ninio- where- pKL = riap inp k l- bK = inp ! (k-1)- bL = inp ! (l+1)- ds = min di dj- dl = max di dj- pIJ = pair inp i j- bI = inp ! (i+1)- bJ = inp ! (j-1)-{-# INLINE tabbedInteriorLoopBase #-}---- | A large number of iloops up to a total maximum loop size of 30 have to be--- checked. There are about 300 cases for j-i>>30. In addition, this loop does--- not like fusion very much and does many different lookups.------ TODO Any performance increase here should yield substantial improvements for--- the overall performance of folding algorithms.------ TODO performance improvement: Split mismatch calculation into its own table--- in the caller. Callee gets an additional argument.------ TODO didjs needs to go back into a *Limit function.--largeInteriorLoopBase :: (FoldFunctions a) => TurnerTables a -> Primary -> Table a -> Int -> Int -> VU.Vector a-largeInteriorLoopBase Turner2004{..} inp strong i j = res where- res = VU.map (\(di,dj) ->- ijmm .*.- (strong ! (i+di,j-dj)) .*.- (iloopL ! (di+dj-2)) .*. -- -2 because di,dj is total IL length +2- (iloopMM ! (riap inp (i+di) (j-dj), inp ! (i+di-1), inp ! (j-dj+1))) .*.- calcNinio maxNinio ninio (abs $ di-dj)- ) didjs- -- NOTE NO FUSION! :-(- -- TODO check this!- didjs = interiorLoopDistances i j- {-- didjs = traceVal (show (i,j)) $ VU.unfoldr (\(di,dj) ->- if di>maxc- then Nothing- else if di+dj>=nexc- then Just ((di,dj),(di+1,3 ))- else Just ((di,dj),(di ,dj+1))- ) (3,5) -}- -- constant- {-- maxc = min 27 (j-i-16)- nexc = min 30 (j-i-13)- -}- ijmm = iloopMM ! (pIJ,bI,bJ)- pIJ = pair inp i j- bI = inp ! (i+1)- bJ = inp ! (j-1)-{-# INLINE largeInteriorLoopBase #-}---- | A bulge on the left side of the inner closing pair. 1-bulges are treated--- as tabulated interior loops.------ (..((...)))--- 01234567890--bulgeLBase :: (FoldFunctions a) => TurnerTables a -> Primary -> Table a -> Int -> Int -> VU.Vector a-bulgeLBase Turner2004{..} inp strong i j = res where- res = VU.map (\k ->- strong ! (i+1+k,j-1) .*.- bulgeL ! k .*.- calcTermAU termAU (riap inp (i+1+k) (j-1)) .*.- tAUij- ) . uncurry VU.enumFromN $ bulgeLimit i j- tAUij = calcTermAU termAU $ pair inp i j-{-# INLINE bulgeLBase #-}---- | A bulge on the right side of the inner closing pair. Mirroring bulgeLBase------ ((~)...)--bulgeRBase :: (FoldFunctions a) => TurnerTables a -> Primary -> Table a -> Int -> Int -> VU.Vector a-bulgeRBase Turner2004{..} inp strong i j = res where- res = VU.map (\k ->- strong ! (i+1,j-1-k) .*.- bulgeL ! k .*.- calcTermAU termAU (riap inp (i+1) (j-1-k)) .*.- tAUij- ) . uncurry VU.enumFromN $ bulgeLimit i j- tAUij = calcTermAU termAU $ pair inp i j-{-# INLINE bulgeRBase #-}---- | 1xn iloops work a bit like bulges. They are more complicated because we--- have to consider 'ninio' values and terminal mismatches.------ TODO how big an improvement would 'iloopL1xn' be with integrated 'ninio'--- values?------ (...(~).)--interior1xnLBase Turner2004{..} inp strong i j = res where- res = VU.map (\k ->- strong ! (i+1+k,j-2) .*.- iloopL ! (k+1) .*.- iloop1xnMM ! (riap inp (i+1+k) (j-2), inp ! (i+k), inp ! (j-1)) .*.- calcNinio maxNinio ninio (k-1) .*.- pIJmm- ) . uncurry VU.enumFromN $ iloop1xnLimit i j- pIJmm = iloop1xnMM ! (pair inp i j, inp ! (i+1), inp ! (j-1))-{-# INLINE interior1xnLBase #-}---- | A mirror to the above.------ (.(~)...)--interior1xnRBase Turner2004{..} inp strong i j = res where- res = VU.map (\k ->- strong ! (i+2,j-1-k) .*.- iloopL ! (k+1) .*.- iloop1xnMM ! (riap inp (i+2) (j-1-k), inp ! (j-k), inp ! (i+1)) .*.- calcNinio maxNinio ninio (k-1) .*.- pIJmm- ) . uncurry VU.enumFromN $ iloop1xnLimit i j- pIJmm = iloop1xnMM ! (pair inp i j, inp ! (i+1), inp ! (j-1))-{-# INLINE interior1xnRBase #-}---- | Close a multibranch loop by trying to combine one element of 'm' with one--- of 'm1'.------ <[[...]][[...]]>--- 0123456789012345--- 1--multibranchCloseBase :: (FoldFunctions a) => TurnerTables a -> Primary -> Table a -> Table a -> Int -> Int -> VU.Vector a-multibranchCloseBase Turner2004{..} inp m m1 i j = res where- res = VU.map (\k ->- m ! (i+1,k) .*.- m1 ! (k+1,j-1) .*.- ijmm .*.- mbcl- ) . uncurry VU.enumFromN $ multibranchCloseLimit i j- ijmm = multiMM ! (pIJ,bJ,bI)- mbcl = multiOffset .*. multiHelix- pIJ = riap inp i j- bI = inp ! (i+1)- bJ = inp ! (j-1)-{-# INLINE multibranchCloseBase #-}---- | Adds a multibranch loop at exactly (i,j).--multibranchIJLoopBase Turner2004{..} inp strong i j = res where- res = VU.singleton $ strong ! (i,j) .*. mbrhlx .*. mbrmm- mbrhlx = multiHelix- mbrmm = multiMM ! (pIJ,bI,bJ) -- TODO correct orientation?- pIJ = pair inp i j- bI = inp ! (i-1)- bJ = inp ! (j+1)-{-# INLINE multibranchIJLoopBase #-}---- | Trivial function that adds a single unpaired nucleotide within a--- multibranch loop.--multibranchUnpairedJBase Turner2004{..} inp mtable i j = res where- res = VU.singleton $ mtable ! (i,j-1) .*. mbrup- mbrup = multiNuc-{-# INLINE multibranchUnpairedJBase #-}---- | Adds a multibranch loop at (k,j), with k-i unpaired nucleotides to the--- left of 'k'.--multibranchKJHelixBase :: (FoldFunctions a) => TurnerTables a -> Primary -> Table a -> Int -> Int -> VU.Vector a-multibranchKJHelixBase Turner2004{..} inp strong i j = res where- res = VU.map (\k ->- multiNuc .^. (k-i) .*.- strong ! (k,j) .*.- multiHelix .*.- multiMM ! (pair inp k j, inp ! (k-1), inp ! (j+1))- ) . uncurry VU.enumFromN $ multibranchKJHelixLimit i j-{-# INLINE multibranchKJHelixBase #-}---- |--multibranchAddKJHelixBase Turner2004{..} inp table strong i j = res where- res = VU.map (\k ->- table ! (i,k) .*.- strong ! (k+1,j) .*.- multiMM ! (pair inp (k+1) j, inp ! k, inp ! (j+1))- ) . uncurry VU.enumFromN $ multibranchAddKJHelixLimit i j-{-# INLINE multibranchAddKJHelixBase #-}---- | [[...]]--- 0123456--externalLoopBase Turner2004{..} inp strong i j = res where- res = VU.singleton $ strong ! (i,j) .*. mm .*. tAU- n = snd $ bounds inp- pIJ = pair inp i j- bI = inp ! (i-1)- bJ = inp ! (j+1)- tAU = calcTermAU termAU pIJ- mm- | i>0&&j<n = extMM ! (pIJ,bI,bJ)- | i>0 = dangle5 ! (pIJ,bI)- | j<n = dangle3 ! (pIJ,bJ)- | otherwise = one-{-# INLINE externalLoopBase #-}---- |--externalAddLoopBase :: (FoldFunctions a) => TurnerTables a -> Primary -> Table a -> Table a -> Int -> Int -> VU.Vector a-externalAddLoopBase trnr@Turner2004{..} inp strong extern i j = res where- res = VU.map (\k ->- externalLoopOpt trnr inp strong i k .*.- extern ! (k+1,j)- ) . uncurry VU.enumFromN $ externalAddLoopLimit i j-{-# INLINE externalAddLoopBase #-}------ * Helper Functions.---- TODO We definitely need to check that this works!--pair :: Primary -> Int -> Int -> ViennaPair-pair inp i j- = assert (checkBounds inp i && checkBounds inp j)- $ toPair (inp `unsafeIndex` i) (inp `unsafeIndex` j)-{-# INLINE pair #-}--riap inp i j- = assert (i>=0 && j>=0 && checkBounds inp i && checkBounds inp j)- $ toPair (inp ! j) (inp ! i)-{-# INLINE riap #-}--ringSum :: (Ring a, VU.Unbox a) => VU.Vector a -> a-ringSum v = VU.foldl' (.+.) zero v-{- INLINE ringSum #-}--ringSumL :: (Ring a, VU.Unbox a) => [a] -> a-ringSumL v = foldl' (.+.) zero v-{-# INLINE ringSumL #-}--ringProduct :: (Ring a, VU.Unbox a) => VU.Vector a -> a-ringProduct v = VU.foldl' (.*.) one v-{-# INLINE ringProduct #-}--ringProductL :: (Ring a, VU.Unbox a) => [a] -> a-ringProductL v = foldl' (.*.) one v---- | Explicit index generator for interior loops. Does not create indices for:--- - bulges 0 k--- - 1xn loops 1 k--- - 2x3 loops 2 3, 3 2--- - tabulated 1 1, 1 2, 2 1, 2 2--interiorLoopDistances i j =- VU.concatMap (- \d -> VU.map (\d' -> (d',d-d')) -- written as a tuple- $ VU.enumFromN 3 (d-5)) -- for each distance, all possible left/right combinations- $ VU.enumFromN 8 (min 23 (j-i-13)) -- diagonal distance or number of unpaired nucleotides -2.-{-# INLINE interiorLoopDistances #-}---- WTF ?! the function below is 10.000x slower than the function above!--interiorLoopIndices :: Int -> Int -> VU.Vector (Int,Int)-interiorLoopIndices !i !j = VU.map (\(k,l) -> (i+k,j-l)) $ interiorLoopDistances i j-{-# INLINE interiorLoopIndices #-}------- TODO filter 'ili' to accept only indices that are not too close. Does--tabbedInteriorLoopDistances :: Int -> Int -> VU.Vector (Int,Int)-tabbedInteriorLoopDistances i j- | j-i>=8 = VU.fromList [(0,1),(1,0),(1,1),(1,2),(2,1),(2,2),(2,3),(3,2)]- | otherwise = VU.empty-{-# INLINE tabbedInteriorLoopDistances #-}---- | Yes, therer is the (+1,+1) and yes, this is inconsistent with the other function--tabbedInteriorLoopIndices i j = VU.map (\(di,dj) -> (i+di+1,j-dj-1)) $ tabbedInteriorLoopDistances i j-{-# INLINE tabbedInteriorLoopIndices #-}---- * All limits depending on (i,j) should be here.---- | Bulge limitations--bulgeLimit :: Int -> Int -> (Int,Int)-bulgeLimit i j = (2,min 29 $ j-i-9)-{-# INLINE bulgeLimit #-}---- | 1xn iloop limits--iloop1xnLimit :: Int -> Int -> (Int,Int)-iloop1xnLimit i j = (3,min 26 $ j-i-10)-{-# INLINE iloop1xnLimit #-}---- | closing of a multibranch loop--multibranchCloseLimit :: Int -> Int -> (Int,Int)-multibranchCloseLimit i j = (i+1,j-i-2) -- (i+7, j-i-14)-{-# INLINE multibranchCloseLimit #-}---- | add mb helix--multibranchAddKJHelixLimit :: Int -> Int -> (Int,Int)-multibranchAddKJHelixLimit i j = (i+1,j-i-1)-{-# INLINE multibranchAddKJHelixLimit #-}---- | mb helix--multibranchKJHelixLimit :: Int -> Int -> (Int,Int)-multibranchKJHelixLimit i j = (i,j-i)-{-# INLINE multibranchKJHelixLimit #-}---- | add external loop--externalAddLoopLimit :: Int -> Int -> (Int,Int)-externalAddLoopLimit i j = (i+5,j-i-5)-{-# INLINE externalAddLoopLimit #-}
+ BioInf/ViennaRNA.hs view
@@ -0,0 +1,2 @@++module BioInf.ViennaRNA where
+ BioInf/ViennaRNA/Energy.hs view
@@ -0,0 +1,117 @@+{-# LANGUAGE PatternGuards #-}+{-# LANGUAGE NoMonomorphismRestriction #-}+{-# LANGUAGE BangPatterns #-}++module BioInf.ViennaRNA.Energy where++import Data.Vector.Fusion.Stream.Monadic as SM+import qualified Data.Vector.Unboxed as VU+import Control.Lens+import Data.Array.Repa.Index+import Prelude as P+import qualified Data.Map as M++import Data.PrimitiveArray as PA hiding ((!))+import Data.PrimitiveArray.Zero as PA+import qualified Data.PrimitiveArray as PA+import Biobase.Turner+import Biobase.Vienna+import Biobase.Primary++import BioInf.ViennaRNA.Signature++import Debug.Trace++++mfe :: Monad m => Signature m Deka Deka+mfe = (hairpin,interior,multi,blockStem,blockUnpair,compsBR,compsBC,structW,structCS,structWS,structOpen,h) where+ hairpin ener l lp xs rp r+ | len <= 6+ , Just e <- (l `VU.cons` xs `VU.snoc` r) `M.lookup` _hairpinLookup ener = e+ | len < 3 = huge+ | len == 3 = (ener^.hairpinL) VU.! len + tAU+ | len < 31 = (ener^.hairpinL) VU.! len + ener^.hairpinMM!(Z:.l:.r:.lp:.rp)+ | otherwise = huge+ where+ !len = VU.length xs+ !tAU = if (l,r) == (nC,nG) || (l,r) == (nG,nC) then Deka 0 else ener^.termAU+ interior ener l ls li w ri rs r+ | lls==0 && lrs==0 -- stack+ = w + _stack ener ! (Z:.l:.r:.ri:.li) -- left, right, right inner, left inner+ | lls==1 && lrs==0 || lls==0 && lrs==1 -- stack with slip+ = w + _stack ener ! (Z:.l:.r:.ri:.li) + _bulgeL ener VU.! 1+ | lls==1 && lrs==1+ = w + _iloop1x1 ener ! (Z:.l:.r:.ri:.li:.lH:.rL)+ | lls==1 && lrs==2+ = w + _iloop2x1 ener ! (Z:.l:.r:.ri:.li:.lH:.rH:.rL)+ | lls==2 && lrs==1+ = w + _iloop2x1 ener ! (Z:.l:.r:.ri:.li:.rH:.lH:.lL)+ | lls==2 && lrs==2+ = w + _iloop2x2 ener ! (Z:.l:.r:.ri:.li:.lH:.lL:.rH:.rL)+ | min lls lrs == 2 && max lls lrs == 3+ = w + _iloop2x3MM ener ! (Z:.l:.r:.lH:.lL) + _iloop2x3MM ener ! (Z:.ri:.li:.rL:.rH) + _iloopL ener VU.! 5 + _ninio ener+ | lls==0 && lrs > 1 && lrs <= 30+ = w + tAU + _bulgeL ener VU.! lrs + tUA+ | lrs==0 && lls > 1 && lls <= 30+ = w + tAU + _bulgeL ener VU.! lls + tUA+ | lrs==1 && lls > 2 && lls <= 30+ = w + _iloop1xnMM ener ! (Z:.li:.ri:.lL:.rH) + _iloop1xnMM ener ! (Z:.r:.l:.rL:.lH) + _iloopL ener VU.! lls + min (_ninio ener *. (lls-1)) (_maxNinio ener)+ | lls==1 && lrs > 2 && lrs <= 30+ = w + _iloop1xnMM ener ! (Z:.li:.ri:.lL:.rH) + _iloop1xnMM ener ! (Z:.r:.l:.rL:.lH) + _iloopL ener VU.! lrs + min (_ninio ener *. (lrs-1)) (_maxNinio ener)+ | lls>0 && lrs>0 && lls+lrs <= 30 -- TODO missing support for length constraints ?+ = w + _iloopMM ener ! (Z:.l:.r:.lH:.rL) + _iloopMM ener ! (Z:.ri:.li:.rH:.lL) + _iloopL ener VU.! (lls+lrs) + min (_ninio ener *. (abs $ lls - lrs)) (_maxNinio ener)+ | otherwise = huge -- NOTE later on, we should never get this score+ where+ !lls = VU.length ls+ !lrs = VU.length rs+ !tAU = if (l,r) `P.elem` [(nC,nG), (nG,nC)] then Deka 0 else ener^.termAU+ !tUA = if (li,ri) `P.elem` [(nC,nG), (nG,nC)] then Deka 0 else ener^.termAU+ lH = VU.unsafeHead ls+ lL = VU.unsafeLast ls+ rH = VU.unsafeHead rs+ rL = VU.unsafeLast rs+ multi ener l li b c ri r+ = b + c + _multiMM ener ! (Z:.r:.l:.ri:.li) + _multiHelix ener + _multiOffset ener where+ blockStem ener lo l s r ro+ = s + _multiMM ener ! (Z:.l:.r:.lo:.ro) + _multiHelix ener+ blockUnpair ener c b+ = b + _multiNuc ener+ compsBR ener b reg+ = let Deka nuc = _multiNuc ener in b + (Deka $ nuc * (VU.length reg))+ compsBC ener b c+ = b + c+ structW ener w+ = w+ structCS ener c w+ = w+ structWS ener w s+ = w + s+ structOpen ener r+ = 0+ h = foldl' min huge+ {-# INLINE hairpin #-}+ {-# INLINE interior #-}+ {-# INLINE multi #-}+ {-# INLINE blockStem #-}+ {-# INLINE blockUnpair #-}+ {-# INLINE compsBR #-}+ {-# INLINE compsBC #-}+ {-# INLINE structW #-}+ {-# INLINE structCS #-}+ {-# INLINE structWS #-}+ {-# INLINE structOpen #-}+ {-# INLINE h #-}+{-# INLINE mfe #-}++huge = Deka 999999+{-# INLINE huge #-}++infixl 8 !+(!) = (PA.!)+{-# INLINE (!) #-}++(*.) :: Deka -> Int -> Deka+(Deka k) *. n = Deka $ k*n+{-# INLINE (*.) #-}+
+ BioInf/ViennaRNA/Eval.hs view
@@ -0,0 +1,95 @@+{-# LANGUAGE PatternGuards #-}++-- Direct evaluation of the energy of a given structure. The RNAfold-based+-- variant finds the optimal subset of base pairs that conform to the given+-- structure, this algorithm gives the energy of exactly the given structure.++module BioInf.ViennaRNA.Eval where++import Data.Vector.Fusion.Util (Id(..))+import qualified Data.Vector.Unboxed as VU+import Text.Printf++import Biobase.Primary+import Biobase.Secondary+import Biobase.Secondary.Diagrams+import Biobase.Vienna++import BioInf.ViennaRNA.Signature+import BioInf.ViennaRNA.Energy++import Debug.Trace++++rnaEval ener s d1s = flatten $ eval mfe ener s d1s++flatten :: SSTree PairIdx Structure -> (Deka, [String])+flatten = f where+ unDeka (Deka e) = e+ f (SSExt l (External e) xs) =+ let etot = e + sum (map fst ys)+ ys = map f xs+ here = printf "External loop: %d" (unDeka e)+ in (etot, here : concatMap snd ys)+ f (SSTree _ (Hairpin e l us r) [] ) = (e, [printf "Hairpin loop: %d" (unDeka e)])+ f (SSTree _ (Interior e l ls ll rr rs r) [y]) =+ let etot = e + fst (f y)+ in (etot, printf "Interior loop: %d" (unDeka e) : snd (f y))+ f (SSTree _ (Multi e ll l r rr) ys) =+ let etot = e + sum (map (fst . f) ys)+ in (etot, printf "Multi loop: %d %s" (unDeka e) (concatMap show [ll,l,r,rr]) : concatMap (snd . f) ys)+ {-+ f (SSTree p e xs) =+ let etot = e + sum (map fst ys)+ ys = map f xs+ here+ | null xs = printf "Hairpin loop: %d" (unDeka e)+ | [_] <- xs = printf "Interior loop: %d" (unDeka e)+ | otherwise = printf "Multibranched loop: %d" (unDeka e)+ in (etot, here : concatMap snd ys)+ -}++data Structure+ = External Deka+ | Hairpin Deka Nuc Primary Nuc+ | Interior Deka Nuc Primary Nuc Nuc Primary Nuc+ | Multi Deka Nuc Nuc Nuc Nuc++eval :: Signature Id Deka Deka -> Vienna2004 -> Primary -> D1Secondary -> SSTree PairIdx Structure+eval efun ener s d1s = annotateWithEnergy t where+ t = d1sTree d1s+ (hairpin,interior,multi,blockStem,blockUnpair,compsBR,compsBC,structW,structCS,structWS,structOpen,h) = efun+ annotateWithEnergy :: SSTree PairIdx () -> SSTree PairIdx Structure+ annotateWithEnergy (SSExt l () xs) = SSExt l e (map annotateWithEnergy xs) where+ e = External 0 -- TODO sum of all external loop energies+ annotateWithEnergy err@(SSTree (i,j) () xs)+ -- hairpin+ | null xs+ = let pri = VU.slice (i+1) (j-i-1) s in SSTree (i,j) (Hairpin (hairpin ener si sii pri jjs sj) si pri sj) []+ -- interior loop+ | [SSTree (k,l) () _] <- xs+ = let kks = s VU.! k; sll = s VU.! l+ e = interior ener si ls kks 0 sll rs sj+ ls = VU.slice (i+1) (k-i-1) s+ rs = VU.slice (l+1) (j-l-1) s+ in SSTree (i,j) (Interior e si ls kks sll rs sj) (map annotateWithEnergy xs)+ -- multibranched loop+ | otherwise+ = let e = multi ener si sii 0 0 jjs sj+ + sum (map bStem xs)+ + sum (map (\c -> blockUnpair ener c 0) cs)+ cs = [] -- TODO all unpaired nucleotides+ bStem (SSTree (k,l) () _) =+ let kks = s VU.! (k-1)+ sk = s VU.! k+ sl = s VU.! l+ sll = s VU.! (l+1)+ in blockStem ener kks sk 0 sl sll+ in SSTree (i,j) (Multi e si sii jjs sj) (map annotateWithEnergy xs)+ where+ si = s VU.! i+ sj = s VU.! j+ sii = s VU.! (i+1)+ jjs = s VU.! (j-1)+
+ BioInf/ViennaRNA/Fold.hs view
@@ -0,0 +1,192 @@+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE BangPatterns #-}++module BioInf.ViennaRNA.Fold where++import Data.Vector.Fusion.Util (Id (..))+import Data.Vector.Fusion.Stream.Monadic as SM+import qualified Data.Vector.Unboxed as VU+import Data.Array.Repa.Index+import Control.Monad+import Control.Monad.ST+import System.IO.Unsafe+import Prelude as P hiding (Maybe(..))+import Data.Strict.Maybe+import Data.Strict.Tuple++import Biobase.Secondary.Diagrams+import Data.Array.Repa.Index.Subword+import ADP.Fusion+import ADP.Fusion.Table+import Biobase.Vienna+import Biobase.Primary+import Data.PrimitiveArray as PA hiding ((!))+import Data.PrimitiveArray.Zero as PA++import BioInf.ViennaRNA.Signature+import BioInf.ViennaRNA.Energy++++basepairing :: Primary -> Subword -> Bool+basepairing inp (Subword(i:.j)) = i+1<j && f (inp VU.! i) (inp VU.! (j-1)) where+ f l r = l==nC && r==nG+ || l==nG && r==nC+ || l==nA && r==nU+ || l==nU && r==nA+ || l==nG && r==nU+ || l==nU && r==nG+ {-# INLINE f #-}+{-# INLINE basepairing #-}++structureConstrains :: Maybe D1Secondary -> Subword -> Bool+structureConstrains Nothing !_ = True+structureConstrains !(Just (D1S c)) (Subword (i:.j)) = (i<j) && (VU.unsafeIndex c i == j-1)+{-# INLINE structureConstrains #-}++structC :: Primary -> Subword -> Bool+structC inp (Subword(i:.j)) = VU.length inp == j+{-# INLINE structC #-}++-- TODO need to fix sized regions, then we are good to go -- performance-wise+--+-- TODO backtracking+--+-- TODO struct table+--+-- TODO restrict structs to a linear band++gRNAfold ener (hairpin,interior,multi,blockStem,blockUnpair,compsBR,compsBC,structW,structCS,structWS,structOpen,h) weak block comps struct cs inp =+ ( weak ,+ hairpin ener <<< c % pr % hr % pl % c |||+ interior ener <<< c % ir % pr % weak % pl % ir % c |||+ multi ener <<< c % pl % block % comps % pl % c `check` (basepairing inp) `check` (structureConstrains cs) ... h+ , block ,+ blockStem ener <<< pl % c % weak % c % pr |||+ blockUnpair ener <<< c % block ... h+ , comps ,+ compsBR ener <<< block % r |||+ compsBC ener <<< block % comps ... h+ , struct ,+-- structW ener <<< weak ||| -- TODO peak left/right with default ; not needed anymore+ structCS ener <<< c % struct |||+ structWS ener <<< weak % struct ||| -- peak here for weak, too+ structOpen ener <<< r `check` (structC inp) ... h+ ) where c = chr inp+ r = region inp+ pr = peekR inp+ pl = peekL inp+ hr = sregion 3 30 inp+ ir = sregion 0 20 inp+ {-# INLINE c #-}+ {-# INLINE r #-}+ {-# INLINE pr #-}+ {-# INLINE pl #-}+ {-# INLINE hr #-}+ {-# INLINE ir #-}+{-# INLINE gRNAfold #-}++++pretty :: Monad m => Signature m String (SM.Stream m String)+pretty = (hairpin,interior,multi,blockStem,blockUnpair,compsBR,compsBC,structW,structCS,structWS,structOpen,h) where+ hairpin _ _ _ r _ _ = "(" P.++ (P.replicate (VU.length r) '.') P.++ ")"+ interior _ _ l _ w _ r _ = "(" P.++ (P.replicate (VU.length l) '.') P.++ w P.++ (P.replicate (VU.length r) '.') P.++ ")"+ multi _ _ _ b c _ _ = "(" P.++ b P.++ c P.++ ")"+ blockStem _ _ _ w _ _ = w+ blockUnpair _ _ b = "." P.++ b+ compsBR _ b r = b P.++ (P.replicate (VU.length r) '.')+ compsBC _ b c = b P.++ c+ structW _ w = w+ structCS _ _ w = "." P.++ w+ structWS _ w s = w P.++ s+ structOpen _ r = P.replicate (VU.length r) '.'+ h = return . id++type CombSignature m e b = Signature m (e, m (SM.Stream m b)) (SM.Stream m b)+++(<**)+ :: (Monad m, Eq b, Eq e) -- , Show e, Show (m [b]))+ => Signature m e e+ -> Signature m b (SM.Stream m b)+ -> CombSignature m e b+(<**) f s = (hairpin,interior,multi,blockStem,blockUnpair,compsBR,compsBC,structW,structCS,structWS,structOpen,h) where+ (hairpinF,interiorF,multiF,blockStemF,blockUnpairF,compsBRF,compsBCF,structWF,structCSF,structWSF,structOpenF,hF) = f+ (hairpinS,interiorS,multiS,blockStemS,blockUnpairS,compsBRS,compsBCS,structWs,structCSS,structWSS,structOpenS,hS) = s+ + xs >>>= f = xs >>= return . SM.map f+ ccm2 xs ys f = xs >>= \xx -> ys >>= \yy -> return $ SM.concatMap (\x -> SM.map (\y -> f x y) yy) xx++ hairpin ener l lp xs rp r = (hairpinF ener l lp xs rp r, return $ SM.singleton $ hairpinS ener l lp xs rp r)+ interior ener l ls li (wF,wS) ri rs r = (interiorF ener l ls li wF ri rs r, wS >>>= \w -> interiorS ener l ls li w ri rs r)+ multi ener l li (bF,bS) (cF,cS) ri r = (multiF ener l li bF cF ri r, ccm2 bS cS $ \b c -> multiS ener l li b c ri r)+ blockStem ener lo l (sF,sS) r ro = (blockStemF ener lo l sF r ro, sS >>>= \s -> blockStemS ener lo l s r ro)+ blockUnpair ener c (bF,bS) = (blockUnpairF ener c bF, bS >>>= \s -> blockUnpairS ener c s)+ compsBR ener (bF,bS) reg = (compsBRF ener bF reg, bS >>>= \s -> compsBRS ener s reg)+ compsBC ener (bF,bS) (cF,cS) = (compsBCF ener bF cF, ccm2 bS cS $ \b c -> compsBCS ener b c)+ structW ener (wF,wS) = (structWF ener wF, wS >>>= \w -> structWs ener w)+ structCS ener c (wF,wS) = (structCSF ener c wF, wS >>>= \w -> structCSS ener c w)+ structWS ener (wF,wS) (sF,sS) = (structWSF ener wF sF, ccm2 wS sS $ \w s -> structWSS ener w s)+ structOpen ener r = (structOpenF ener r, return . SM.singleton $ structOpenS ener r)+ h xs = do+ hfs <- hF $ SM.map P.fst xs+ let phfs = SM.concatMapM P.snd . SM.filter ((hfs==) . P.fst) $ xs+ hS phfs++rnaFoldConstrained :: Vienna2004 -> Primary -> D1Secondary -> (Deka,[String])+rnaFoldConstrained ener inp s = (struct ! (Z:.subword 0 n), bt) where+ (_,Z:.Subword (_:.n)) = bounds weak+ len = VU.length inp+ (weak,block,comps,struct) = unsafePerformIO (rnaFoldFill ener (Just s) inp)+ bt = backtrack ener (Just s) inp (weak,block,comps,struct)+{-# NOINLINE rnaFoldConstrained #-}++rnaFold :: Vienna2004 -> Primary -> (Deka,[String])+rnaFold ener inp = (struct ! (Z:.subword 0 n), bt) where+ (_,Z:.Subword (_:.n)) = bounds weak+ len = VU.length inp+ (weak,block,comps,struct) = unsafePerformIO (rnaFoldFill ener Nothing inp)+ bt = backtrack ener Nothing inp (weak,block,comps,struct)+{-# NOINLINE rnaFold #-}++rnaFoldFill :: Vienna2004 -> Maybe (D1Secondary) -> Primary -> IO (PA.Unboxed (Z:.Subword) Deka, PA.Unboxed (Z:.Subword) Deka, PA.Unboxed (Z:.Subword) Deka, PA.Unboxed (Z:.Subword) Deka)+rnaFoldFill !ener !cs !inp = do+ let n = VU.length inp+ !weak' <- newWithM (Z:.subword 0 0) (Z:.subword 0 n) huge+ !block' <- newWithM (Z:.subword 0 0) (Z:.subword 0 n) huge+ !comps' <- newWithM (Z:.subword 0 0) (Z:.subword 0 n) huge+ !struc' <- newWithM (Z:.subword 0 0) (Z:.subword 0 n) 0+ fillTables $ gRNAfold ener mfe (mTblSw NonEmptyT weak') (mTblSw NonEmptyT block') (mTblSw NonEmptyT comps') (mTblSw NonEmptyT struc') cs inp+ weakF <- freeze weak'+ blockF <- freeze block'+ compsF <- freeze comps'+ strucF <- freeze struc'+ return (weakF,blockF,compsF,strucF)+{-# NOINLINE rnaFoldFill #-}++fillTables (MTbl _ weak, weakF, MTbl _ block, blockF, MTbl _ comps, compsF, MTbl _ struc, strucF) = do+ let (_,Z:.Subword (0:.n)) = boundsM weak+ forM_ [n,n-1..0] $ \i -> forM_ [i..n] $ \j -> do+ weakF (subword i j) >>= writeM weak (Z:.subword i j)+ blockF (subword i j) >>= writeM block (Z:.subword i j)+ compsF (subword i j) >>= writeM comps (Z:.subword i j)+ strucF (subword i j) >>= writeM struc (Z:.subword i j)+{-# INLINE fillTables #-}++-- * backtracking++backtrack ener cs (inp :: Primary) (weak :: PA.Unboxed (Z:.Subword) Deka, block :: PA.Unboxed (Z:.Subword) Deka, comps :: PA.Unboxed (Z:.Subword) Deka, struct :: PA.Unboxed (Z:.Subword) Deka) = unId . SM.toList . unId $ sF $ subword 0 n where+ n = VU.length inp+ w :: SwBtTbl Id Deka String+ w = btTbl NonEmptyT weak (wF :: Subword -> Id (SM.Stream Id String))+ b :: SwBtTbl Id Deka String+ b = btTbl NonEmptyT block (bF :: Subword -> Id (SM.Stream Id String))+ c :: SwBtTbl Id Deka String+ c = btTbl NonEmptyT comps (cF :: Subword -> Id (SM.Stream Id String))+ s :: SwBtTbl Id Deka String+ s = btTbl NonEmptyT struct (sF :: Subword -> Id (SM.Stream Id String))+ (_,wF,_,bF,_,cF,_,sF) = gRNAfold ener (mfe <** pretty) w b c s cs inp+{-# INLINE backtrack #-}+
+ BioInf/ViennaRNA/Signature.hs view
@@ -0,0 +1,38 @@++module BioInf.ViennaRNA.Signature where++import Data.Vector.Fusion.Stream.Monadic as SM++import Biobase.Primary+import Biobase.Vienna++++type Signature m a r =+ -- weak / hairpin+ ( Vienna2004 -> Nuc -> Nuc -> Primary -> Nuc -> Nuc -> a+ -- weak / interior+ , Vienna2004 -> Nuc -> Primary -> Nuc -> a -> Nuc -> Primary -> Nuc -> a+ -- weak / multibranch+ , Vienna2004 -> Nuc -> Nuc -> a -> a -> Nuc -> Nuc -> a+ -- block / multistem+ , Vienna2004 -> Nuc -> Nuc -> a -> Nuc -> Nuc -> a+ -- block / unpaired+ , Vienna2004 -> Nuc -> a -> a+ -- comps / block region+ , Vienna2004 -> a -> Primary -> a+ -- comps / block comps+ , Vienna2004 -> a -> a -> a+ -- struct / weak+ , Vienna2004 -> a -> a+ -- struct / char-struct+ , Vienna2004 -> Nuc -> a -> a+ -- struct / weak-struct+ , Vienna2004 -> a -> a -> a+ -- struct / open+ , Vienna2004 -> Primary -> a+ -- all / objective+ , Stream m a -> m r+ )++
+ README.md view
@@ -0,0 +1,39 @@++ViennaRNA RNAfold v2, MFE variant+using the ADPfusion library++++Introduction+============++This algorithm is the second, and much larger, test case for ADPfusion. We+implement "RNAfold v2" in the MFE variant using "-d2" dangles. Both a library+version and an executable are created. The "RNAFold" binary expects single+sequences, one per line. Backtracking tracks all co-optimal structures.++++Installation+============++A simple "cabal update && cabal-dev install RNAFold" should be enough.++++Runtime notes+=============++Using Haskell and ADPfusion, we come to within x3-x4 for this package. Between+the initial test case / submission (in 0.0.0.3) I have traded in some+performance improvements for much better readability in BioInf.RNAfold.Energy.+The C version of RNAfold employs some other methods to improve performance.+Consider:++base -~+ inner-1 +~- base+base -~+ inner-2 +~- base++where it is advantageous to calculate the outer basepair only once, not twice+as we are doing. It is probably better to try to improve the handling of+fusioned code and/or final assembler generation than finding calculations+common to different parts of CFG's.
RNAFold.cabal view
@@ -1,64 +1,101 @@ name: RNAFold-version: 0.0.2.1-author: Christian Hoener zu Siederdissen (Haskell), Ivo L. Hofacker et al (ViennaRNA)+version: 1.99.3.4+author: Christian Hoener zu Siederdissen (Haskell), Ivo L. Hofacker et al (ViennaRNA), 2010-2013+copyright: Christian Hoener zu Siederdissen, 2010-2013+homepage: http://www.tbi.univie.ac.at/~choener/adpfusion maintainer: choener@tbi.univie.ac.at-copyright: Christian Hoener zu Siederdissen, 2010 category: Bioinformatics-synopsis: RNA secondary structure prediction license: GPL-3 license-file: LICENSE build-type: Simple stability: experimental-cabal-version: >= 1.4.0+cabal-version: >= 1.8.0+synopsis: RNA secondary structure prediction description:- Provides the folding functions as used in the ViennaRNA- package. Here, they are in Haskell form to be used by Haskell- programs.+ RNAfold v2 using the ADPfusion library. The RNAfold algorithm+ is used to determine how fast we can be compared to a highly+ optimized C program. .- - This is a release aimed at testing Data.Vector- - Expect major performance issues with GHC < 6.13!+ Please use GHC 7.6 or newer.+ .+ NOTE I'd like to rename this package to RNAfold, like the C+ implementation. Do not install "globally", especially if you+ normally use RNAfold from the ViennaRNA package, for obvious+ reasons. +Extra-Source-Files:+ README.md+ library build-depends:- base >=4 && <5,- containers,- vector >=0.7,- primitive >=0.3,+ base >=4&&<5 ,+ cmdargs >= 0.10 ,+ containers ,+ deepseq >= 1.3 ,+ lens >= 3.8 ,+ primitive >= 0.5 ,+ repa >= 3.2 ,+ strict >= 0.3.2 ,+ vector >= 0.10 ,+ ADPfusion >= 0.2.0.0 ,+ BiobaseTurner >= 0.3.1.1 ,+ BiobaseVienna >= 0.3 ,+ BiobaseXNA >= 0.7 ,+ PrimitiveArray >= 0.5+ exposed-modules:+ BioInf.ViennaRNA+ BioInf.ViennaRNA.Energy+ BioInf.ViennaRNA.Eval+ BioInf.ViennaRNA.Fold+ BioInf.ViennaRNA.Signature+ ghc-options:+ -Odph+ -funbox-strict-fields+ -funfolding-use-threshold100+ -funfolding-keeness-factor100+ -fllvm -optlo-O3 -optlo-inline -optlo-std-compile-opts - Biobase >=0.0.2,- BiobaseTurner >=0.0.2,- BiobaseVienna >=0.0.2,- BiobaseTypes >=0.0.2,- HsTools >=0.0.1.1,- PrimitiveArray >=0.0.2 && <0.0.3+executable RNAFold+ build-depends:+ base >= 4 && < 5 ,+ cmdargs >= 0.10 ,+ BiobaseTurner >= 0.3 ,+ BiobaseVienna >= 0.3 ,+ BiobaseXNA >= 0.7 ,+ RNAFold+ main-is:+ RNAFold.hs+ hs-source-dirs:+ src+ ghc-options:+ -rtsopts+ -Odph+ -funbox-strict-fields+ -funfolding-use-threshold100+ -funfolding-keeness-factor100+ -fllvm -optlo-O3 -optlo-inline -optlo-std-compile-opts - exposed-modules:- BioInf.RNAFold,- BioInf.RNAFold.Energy,- BioInf.RNAFold.EnergyInt,- BioInf.RNAEval,- BioInf.RNAFold.Functions+executable RNAEval+ build-depends:+ base >= 4 && < 5 ,+ cmdargs >= 0.10 ,+ BiobaseTurner >= 0.3 ,+ BiobaseVienna >= 0.3 ,+ BiobaseXNA >= 0.7 ,+ RNAFold+ main-is:+ RNAEval.hs+ hs-source-dirs:+ src+ ghc-options:+ -rtsopts+ -Odph+ -funbox-strict-fields+ -funfolding-use-threshold100+ -funfolding-keeness-factor100+ -fllvm -optlo-O3 -optlo-inline -optlo-std-compile-opts - if impl(ghc > 6.13)- ghc-options:- -Odph- -fllvm- -fforce-recomp- -funbox-strict-fields- -fllvm- -optlo-O3- -optlc-O3- -fdicts-cheap- -fspec-constr- -funbox-strict-fields- -funfolding-use-threshold=100- -funfolding-creation-threshold=100- else- ghc-options:-Odph+source-repository head+ type: git+ location: git://github.com/choener/RNAfold --- -fno-method-sharing -- performance does not improve!--- -fdicts-cheap--- -fspec-constr--- -funbox-strict-fields--- -funfolding-use-threshold=100--- -funfolding-creation-threshold=100
+ src/RNAEval.hs view
@@ -0,0 +1,91 @@+{-# LANGUAGE RecordWildCards #-}+{-# LANGUAGE DeriveDataTypeable #-}++-- | RNAEval tool.++module Main where++++import System.Console.CmdArgs++import Biobase.Primary+import Biobase.Secondary.Diagrams+import Biobase.Vienna+import qualified Biobase.Turner.Import as TI++import BioInf.ViennaRNA.Fold+import BioInf.ViennaRNA.Eval++++data Options+ = Eval+ { params :: String+ }+ | ConstrainedFold+ { params :: String+ }+ deriving (Show,Data,Typeable)++oEval = Eval+ { params = "./params" &= help "Turner 2004 RNA parameters (defaults to ./params)"+ }++oConstrainedFold = ConstrainedFold+ { params = "../params"+ }++main = do+ o <- cmdArgs $ modes [oEval &= auto, oConstrainedFold]+ xs <- fmap lines getContents+ tm <- fmap turnerToVienna $ TI.fromDir (params o) "" ".dat"+ case o of+ Eval{..} -> mapM_ (doEval tm) $ toPairs xs+ ConstrainedFold{..} -> mapM_ (doCF tm) $ toPairs xs++toPairs (x1:x2:xs) = (x1,x2) : toPairs xs+toPairs [x] = error "single last line remaining"+toPairs [] = []++doEval tm (inp,str) = do+ print $ length inp+ print $ rnaEval tm (mkPrimary inp) (mkD1S str)++doCF tm (inp,str) = do+ print $ length inp+ print $ rnaFoldConstrained tm (mkPrimary inp) (mkD1S str)++++test inp str = do+ tm <- fmap turnerToVienna $ TI.fromDir "./params" "" ".dat"+ doEval tm (inp,str)++tests = mapM_ (uncurry test)+ [ ( "CCUGACUGGCGUUGACAUAUGGUU"+ , ".......(((((......)).)))"+ )+ , ( "CUGGGGGUGACAUCCCCCC"+ , "..(((((......)).)))"+ )+ , ( "GGCGUUGACAUAUGGUU"+ , "(((((......)).)))"+ )+ , ( "GGGGUUGACAUACCCCC"+ , "(((((......)).)))"+ )+ , ( "GGCGUUGACAUAUGUU"+ , "(((((......)))))"+ )+ , ( "GGGGGUGACAUCCCCC"+ , "(((((......)))))"+ )+ , ( "GGGGGUGACCCCC"+ , "(((((...)))))"+ )+ , ( "CCUGACUGGCGUUGACAUAUGGUUGCUUGAGCGUAGCCAGGUGUUGGUGGUCCAGUGCAUCAAGGUGCCGUCGGAUCGGAUACUUGGCUUUGCUUAGAUU"+ , ".......(((((......)).)))(.(((((((.(((((((.(((((((.....((((......)))))))).)))......))))))).))))))).)."+ )+ ]+
+ src/RNAFold.hs view
@@ -0,0 +1,37 @@+{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE RecordWildCards #-}++-- | Simple wrapper around the rnafold library.++module Main where++++import System.Console.CmdArgs++import Biobase.Primary+import Biobase.Vienna+import qualified Biobase.Turner.Import as TI++import BioInf.ViennaRNA.Fold++++data Options = Options+ { params :: String+ } deriving (Show,Data,Typeable)++options = Options+ { params = "./params" &= help "Turner 2004 RNA parameters (defaults to ./params)"+ }++main = do+ Options{..} <- cmdArgs options+ xs <- fmap lines getContents+ tm <- fmap turnerToVienna $ TI.fromDir params "" ".dat"+ mapM_ (run' tm) xs++run' tm inp = do+ print $ length inp+ print $ rnaFold tm (mkPrimary inp)+