diff --git a/BioInf/RNAEval.hs b/BioInf/RNAEval.hs deleted file mode 100644 --- a/BioInf/RNAEval.hs +++ /dev/null @@ -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 diff --git a/BioInf/RNAFold.hs b/BioInf/RNAFold.hs deleted file mode 100644 --- a/BioInf/RNAFold.hs +++ /dev/null @@ -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+30&&j do - let extUP = {-# SCC "extUP" #-} if i Int -> a -> [Pairlist] - wea i j d = let bestE = weak `unsafeIndex` (i,j) in - [ [(i,j)] - | i+3 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 () --} diff --git a/BioInf/RNAFold/EnergyInt.hs b/BioInf/RNAFold/EnergyInt.hs deleted file mode 100644 --- a/BioInf/RNAFold/EnergyInt.hs +++ /dev/null @@ -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=0 - , (x,z) <- ext (i+1) j d' - ] ++ - [ (x,z) -- loop at (i,j) - | i=0 - , (x,z) <- str i j d' - ] ++ - [ (x++y,z) - | i=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=0 - , (x,z) <- str (i+1) (j-1) d' - ] ++ - [ ((i,j):x,z) - | i=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=0 - ] ++ - -- - -- interior loops - -- - [ ((i,j):x,z) - | i=0 - , (x,z) <- str k l d' - ] ++ - -- - -- multibranch loops closed at (i,j) - -- - [ ((i,j):x++y,z) - | i=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=0 - , (x,z) <- mul k l d' - ] ++ - -- - -- a helix (k,j) - -- - [ (x,z) - | i=0 - , (x,z) <- str k j d' - ] ++ - -- - -- add a helix to an existing structure - -- - [ (x++y,z) - | i=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=0 - , (x,z) <- str i j d' - ] ++ - -- - -- unpaired nucleotide on the J side - -- - [ (x,z) - | i=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 --} diff --git a/BioInf/RNAFold/Functions.hs b/BioInf/RNAFold/Functions.hs deleted file mode 100644 --- a/BioInf/RNAFold/Functions.hs +++ /dev/null @@ -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&&j0 = dangle5 ! (pIJ,bI) - | j 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 #-} diff --git a/BioInf/ViennaRNA.hs b/BioInf/ViennaRNA.hs new file mode 100644 --- /dev/null +++ b/BioInf/ViennaRNA.hs @@ -0,0 +1,2 @@ + +module BioInf.ViennaRNA where diff --git a/BioInf/ViennaRNA/Energy.hs b/BioInf/ViennaRNA/Energy.hs new file mode 100644 --- /dev/null +++ b/BioInf/ViennaRNA/Energy.hs @@ -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 (*.) #-} + diff --git a/BioInf/ViennaRNA/Eval.hs b/BioInf/ViennaRNA/Eval.hs new file mode 100644 --- /dev/null +++ b/BioInf/ViennaRNA/Eval.hs @@ -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) + diff --git a/BioInf/ViennaRNA/Fold.hs b/BioInf/ViennaRNA/Fold.hs new file mode 100644 --- /dev/null +++ b/BioInf/ViennaRNA/Fold.hs @@ -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 Subword -> Bool +structureConstrains Nothing !_ = True +structureConstrains !(Just (D1S c)) (Subword (i:.j)) = (i 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 #-} + diff --git a/BioInf/ViennaRNA/Signature.hs b/BioInf/ViennaRNA/Signature.hs new file mode 100644 --- /dev/null +++ b/BioInf/ViennaRNA/Signature.hs @@ -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 + ) + + diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -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. diff --git a/RNAFold.cabal b/RNAFold.cabal --- a/RNAFold.cabal +++ b/RNAFold.cabal @@ -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 diff --git a/src/RNAEval.hs b/src/RNAEval.hs new file mode 100644 --- /dev/null +++ b/src/RNAEval.hs @@ -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" + , ".......(((((......)).)))(.(((((((.(((((((.(((((((.....((((......)))))))).)))......))))))).))))))).)." + ) + ] + diff --git a/src/RNAFold.hs b/src/RNAFold.hs new file mode 100644 --- /dev/null +++ b/src/RNAFold.hs @@ -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) +