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/RNAfold.hs b/BioInf/RNAfold.hs new file mode 100644 --- /dev/null +++ b/BioInf/RNAfold.hs @@ -0,0 +1,335 @@ +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE TupleSections #-} +{-# LANGUAGE NoMonomorphismRestriction #-} +{-# LANGUAGE RecordWildCards #-} + +module BioInf.RNAfold where + +import Control.Monad +import Control.Monad.Primitive +import Control.Monad.ST +import Data.Array.Repa.Index +import qualified Data.Vector.Fusion.Stream.Monadic as S +import qualified Data.Vector.Fusion.Stream as P +import qualified Data.Vector.Unboxed as VU +import Control.Monad.State.Lazy +import Control.Arrow (first,second,(***)) +import qualified Data.List as L +import Control.Exception (assert) +import Data.Strict.Tuple hiding (fst,snd) + +import Biobase.Primary +import Biobase.Secondary.Vienna + +import Data.PrimitiveArray +import Data.PrimitiveArray.Unboxed.Zero + +import ADP.Fusion.Monadic +import ADP.Fusion.Monadic.Internal + +import Debug.Trace +import Text.Printf +import GHC.Exts + +import Biobase.Primary +import Biobase.Secondary.Vienna +import Biobase.Vienna +import Biobase.Vienna.Default + +import BioInf.RNAfold.Combinators +import BioInf.RNAfold.Energy +import BioInf.RNAfold.Library + + + +-- | + +testRNAfold :: String -> (Int,[String]) +testRNAfold inp' = struct `seq` bt `seq` (struct!(Z:.0:.n),bt) where + (_,Z:._:.n) = bounds struct + ener = fst turnerRNA2004 + inp = mkPrimary inp' + tbls@(weak,block,comps,struct) = runST (rnafold ener inp) + bt = btRNAfold ener inp tbls +{-# NOINLINE testRNAfold #-} + +-- | +testInput = "cccacccaaagggaaaaggg" +test = testRNAfold testInput + +rnafold :: Vienna2004 -> Primary -> ST s + ( Arr0 DIM2 Int + , Arr0 DIM2 Int + , Arr0 DIM2 Int + , Arr0 DIM2 Int + ) +rnafold ener inp = do + + let !n = let (_,Z:.l) = bounds inp in l+1 + let base = base' inp + {-# INLINE base #-} + let baseLr = baseLr' inp + {-# INLINE baseLr #-} + let baselR = baselR' inp + {-# INLINE baselR #-} + let basepairing = basepairing' inp + {-# INLINE basepairing #-} + let stackpairing = stackpairing' inp + {-# INLINE stackpairing #-} + let reglen = reglen' inp + {-# INLINE reglen #-} + let primary = primary' inp + {-# INLINE primary #-} + let primaryPR = primaryPR' inp + {-# INLINE primaryPR #-} + let primaryPL = primaryPL' inp + {-# INLINE primaryPL #-} + -- this is a bit unfortunate, but otherwise we get type inference problems + let hS :: S.Stream (ST s) Int -> ScalarM (ST s Int) + hS = ScalarM . S.foldl' min (999999::Int) + {-# INLINE hS #-} + + weak <- fromAssocsM (Z:.0:.0) (Z:.n:.n) 999999 [] + block <- fromAssocsM (Z:.0:.0) (Z:.n:.n) 999999 [] + comps <- fromAssocsM (Z:.0:.0) (Z:.n:.n) 999999 [] + struct <- fromAssocsM (Z:.0:.0) (Z:.n:.n) 0 [] + + fillTables + weak ( + multiOF ener <<< baseLr -~+ (multiIF ener <<< block +~+ comps ... hS) +~- baselR ||| + iloopOF ener <<< baseLr -~+ (iloopIF ener <<< primary #~~ weak ~~# primary ... hS) +~- baselR ||| + iloop1NF ener <<< primary ---~+ weak +~@ primary ||| + iloopN1F ener <<< primary @~+ weak +~--- primary ||| + bulgeRF ener <<< baseLr -~+ weak +~* primary ||| + bulgeLF ener <<< primary *~+ weak +~- baselR ||| + tinyloopF ener <<< primaryPR &~+ weak +~& primaryPL ||| + hairpinF ener <<< baseLr -~+ primary +~- baselR ... h `with` basepairing ) + block ( + adjustStream n (justStemF ener <<< baseLr -~+ weak +~- baselR) ||| + regionStemF ener <<< base -~+ block ... h ) + comps ( + bsF <<< block +~+ reglen ||| + bcF <<< block +~+ comps ||| + iD <<< block ... h ) + struct ( + iD <<< weak ||| + rSF <<< base -~~ struct ||| + cmF <<< weak +~+ struct ||| + nilF <<< empty ... h `with` constrained (\(Z:.i:.j) -> j==n) ) + + weak' <- freeze weak + block' <- freeze block + comps' <- freeze comps + struct' <- freeze struct + + return + ( weak' + , block' + , comps' + , struct' + ) +{-# INLINE rnafold #-} + +infixl 9 *~+, +~*, &~+, +~&, ---~+, +~@, +~---, @~+ + +(*~+) = makeLeft_MinRight (3,31) 1 +{-# INLINE (*~+) #-} + +(+~*) = makeMinLeft_Right 1 (3,31) +{-# INLINE (+~*) #-} + +(&~+) = makeLeft_MinRight (1,4) 1 +{-# INLINE (&~+) #-} + +(+~&) = makeMinLeft_Right 1 (1,4) +{-# INLINE (+~&) #-} + +(---~+) = makeLeft_MinRight (3,3) 1 +{-# INLINE (---~+) #-} + +(+~@) = makeMinLeft_Right 1 (5,31) +{-# INLINE (+~@) #-} + +(+~---) = makeMinLeft_Right 1 (3,3) +{-# INLINE (+~---) #-} + +(@~+) = makeLeft_MinRight (5,31) 1 +{-# INLINE (@~+) #-} + +-- * backtracking +-- +-- For now, we replicate the grammar as the optimizer is rather fragile + +btRNAfold + :: Vienna2004 + -> Primary + -> (Arr0 DIM2 Int, Arr0 DIM2 Int, Arr0 DIM2 Int, Arr0 DIM2 Int) + -> [String] +btRNAfold ener inp (weak,block,comps,struct) = structG (Z:.0:.n) where + + !n = let (_,Z:.l) = bounds inp in l+1 + base = base' inp + baseLr = baseLr' inp + baselR = baselR' inp + reglen = reglen' inp + basepairing = basepairing' inp + primary = primary' inp + primaryPR = primaryPR' inp + primaryPL = primaryPL' inp + + -- + + weak' :: DIM2 -> Scalar (Int, [String]) + weak' ij = Scalar (weak!ij, weakG ij) + + block' :: DIM2 -> Scalar (Int, [String]) + block' ij = Scalar (block!ij, blockG ij) + + comps' :: DIM2 -> Scalar (Int, [String]) + comps' ij = Scalar (comps!ij, compsG ij) + + struct' :: DIM2 -> Scalar (Int, [String]) + struct' ij = Scalar (struct!ij, structG ij) + + -- + + weakG :: DIM2 -> [String] + weakG = ( + multiBT ener <<< baseLr -~+ block' +~+ comps' +~- baselR ||| + iloopBT ener <<< baseLr -~+ primary #~~ weak' ~~# primary +~- baselR ||| + iloop1NBT ener <<< primary ---~+ weak' +~@ primary ||| + iloopN1BT ener <<< primary @~+ weak' +~--- primary ||| + bulgeRBT ener <<< baseLr -~+ weak' +~* primary ||| + bulgeLBT ener <<< primary *~+ weak' +~- baselR ||| + tinyloopBT ener <<< primaryPR &~+ weak' +~& primaryPL ||| + hairpinBT ener <<< baseLr -~+ primary +~- baselR ..@ (hBT weak) `withBT` basepairing + ) + + blockG :: DIM2 -> [String] + blockG = ( + adjustStreamBT n (justStemBT ener <<< baseLr -~+ weak' +~- baselR) ||| + regionStemBT ener <<< base -~+ block' ..@ (hBT block) + ) + + compsG :: DIM2 -> [String] + compsG = ( + bsBT <<< block' +~+ reglen ||| + bcBT <<< block' +~+ comps' ||| + iDBT <<< block' ..@ (hBT comps) + ) + + structG :: DIM2 -> [String] + structG = ( + iDBT <<< weak' ||| + rSBT <<< base -~~ struct' ||| + cmBT <<< weak' +~+ struct' ||| + nilBT <<< empty ..@ (hBT struct `withBT` constrained (\(Z:.i:.j) -> j==n) ) + ) + + multiBT ener l (b,bS) (c,cS) r = + let e = multiOF ener l (multiIF ener b c) r + in (e, ["("++x++y++")" | x<-bS, y<-cS]) + iloopBT ener lo ls@(_:!:li:!:lj) (w,wS) rs@(_:!:ri:!:rj) ro = + let e = iloopOF ener lo (iloopIF ener ls w rs) ro + in (e, L.map (\s -> "("++replicate (lj-li) '.'++"("++s++")"++replicate (rj-ri) '.'++")") wS) + iloop1NBT ener ls (w,wS) rs@(_:!:ri:!:rj) = + let e = iloop1NF ener ls w rs + in (e, L.map (\s -> "(.("++s++")"++replicate (rj-ri-1) '.'++")") wS) + iloopN1BT ener ls@(_:!:li:!:lj) (w,wS) rs = + let e = iloopN1F ener ls w rs + in (e, L.map (\s -> "("++replicate (lj-li-1) '.'++"("++s++").)") wS) + bulgeRBT ener ls (w,wS) rs@(_:!:ri:!:rj) = + let e = bulgeRF ener ls w rs + in (e, L.map (\s -> "("++s++"."++replicate (rj-ri-1) '.'++")") wS) + bulgeLBT ener ls@(_:!:li:!:lj) (w,wS) rs = + let e = bulgeLF ener ls w rs + in (e, L.map (\s -> "("++replicate (lj-li-1) '.'++"."++s++")") wS) + hairpinBT ener llp reg@(xs:!:i:!:j) rpr = + let e = hairpinF ener llp reg rpr + in (e, ["(" ++ replicate (j-i+1) '.' ++ ")"]) + tinyloopBT ener ls@(_:!:li:!:lj) (w,sW) rs@(_:!:ri:!:rj) = + let e = tinyloopF ener ls w rs + in (e, L.map (\s -> "("++replicate (lj-li-1) '.'++s++replicate (rj-ri-1) '.'++")") sW) + + regionStemBT ener nc (w,sW) = + let e = regionStemF ener nc w + in (e, L.map (\s -> "." ++ s) sW) + justStemBT ener llp (w,sW) rpr = + let e = justStemF ener llp w rpr + in (e, sW) + + bcBT (b,bW) (c,cW) = + let e = b+c + in (e,[ x++y | x<-bW, y<-cW ]) + bsBT (b,bW) reg = (b, L.map (++ replicate reg '.') bW) + iDBT = id + + ssBT len = (ssF len, [replicate len '.']) + rSBT n (w,wS) = + let e = rSF n w + in (e, map ("."++) wS) + cmBT (w,wS) (s,sS) = + let e = cmF w s + in (e, [x++y | x<-wS, y<-sS]) + nilBT b = if b then (nilF b, [""]) else (nilF b, []) + + hBT tbl ij = L.concatMap snd . L.filter ((tbl!ij==).fst) . P.toList + + +-- * different energy functions (very simplified signature) + + + +-- * Functions that will be part of a bioinformatics DP library + + + +-- ** + +adjustStream :: Int -> (DIM2 -> S.Stream (ST s) Int) -> DIM2 -> S.Stream (ST s) Int +adjustStream !n sgen (Z:.i:.j) + | i>0 && j (DIM2 -> P.Stream elm) -> DIM2 -> P.Stream elm +adjustStreamBT !n sgen (Z:.i:.j) + | i>0 && j0 && j MArr0 (PrimState m) DIM2 Int -> (DIM2 -> m Int) + -> MArr0 (PrimState m) DIM2 Int -> (DIM2 -> m Int) + -> MArr0 (PrimState m) DIM2 Int -> (DIM2 -> m Int) + -> MArr0 (PrimState m) DIM2 Int -> (DIM2 -> m Int) + -> m () +fillTables aT aF bT bF cT cF dT dF = do + let (_,Z:.n:._) = boundsM aT + forM_ [n,n-1 .. 0] $ \i -> forM_ [i..n] $ \j -> do + let ij = Z:.i:.j + aF ij >>= writeM aT ij + bF ij >>= writeM bT ij + cF ij >>= writeM cT ij + dF ij >>= writeM dT ij +{-# INLINE fillTables #-} + diff --git a/BioInf/RNAfold/Combinators.hs b/BioInf/RNAfold/Combinators.hs new file mode 100644 --- /dev/null +++ b/BioInf/RNAfold/Combinators.hs @@ -0,0 +1,52 @@ +{-# LANGUAGE TemplateHaskell #-} +{-# LANGUAGE NoMonomorphismRestriction #-} + +-- | RNAfold combinators, extracted for quickcheck + +module BioInf.RNAfold.Combinators + ( (~~#) + , (#~~) + ) where + +import Data.Array.Repa.Index +import qualified Data.Vector.Fusion.Stream as S +import Data.List + +import qualified ADP.Fusion as F +import qualified ADP.Fusion.Monadic as M +import qualified ADP.Fusion.Monadic.Internal as F +import ADP.Fusion.Monadic (makeLeft_MinRight) +import ADP.Fusion.Monadic.Internal (Box(..)) + + + +infixl 9 #~~, ~~# + +-- | The structure on the left is a subword with size 2-28. The maximal size +-- could be 30 but since the two combinators are linked, 29,30 would fail +-- anyways. + +(#~~) = makeLeft_MinRight (3,30) 1 +{-# INLINE (#~~) #-} + +-- | The structure on the right is a subword with size 2-30, however we inspect +-- the stack an reduce the maximal size. + +(~~#) xs ys = Box mk step xs ys where + minT = 8 -- minimal total size of region + minC = 3 -- minimal number of nuc's on the right + maxC = 32 + {-# INLINE mk #-} + mk (z:.k:.j,a,b) = let (_:.i) = z + cnsmd = k-i -- consumed part + l = max k (j-maxC+cnsmd) + in return (z:.k:.l:.j,a,b) + {-# INLINE step #-} + step (z:.k:.l:.j,a,b) + | l<=j-(max 0 $ minT - cnsmd) && l+minC<=j + = return $ S.Yield (z:.k:.l:.j,a,b) (z:.k:.l+1:.j,a,b) + | otherwise = return $ S.Done + where cnsmd = k-i + (_:.i) = z +{-# INLINE (~~#) #-} + diff --git a/BioInf/RNAfold/Energy.hs b/BioInf/RNAfold/Energy.hs new file mode 100644 --- /dev/null +++ b/BioInf/RNAfold/Energy.hs @@ -0,0 +1,252 @@ +{-# LANGUAGE PackageImports #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE RecordWildCards #-} + +-- | A set of energy functions that are modelled after the ViennaRNA package, +-- Version 2 (with d=2). +-- +-- As part of the design, we could have (i) either continued giving many +-- parameters or (ii) have fewer parameters that require input of the +-- (Primary,Index,Index) type. Since the compilation speed of the grammars +-- using these functions depends on the number of arguments (in case (i) +-- compilation takes minutes!), this approach has benefits for testing the +-- fusion library. +-- +-- This means compilation is a lot faster, but runtime is not @2.8x@ slower but +-- @3.5x@ slower. + +module BioInf.RNAfold.Energy where + +import Control.Exception (assert) +import Data.Strict.Tuple hiding (fst,snd) +import qualified Data.Vector.Fusion.Stream.Monadic as S + +import Biobase.Primary +import Biobase.Secondary.Vienna +import Biobase.Vienna +import Data.PrimitiveArray +import "PrimitiveArray" Data.Array.Repa.Index + +import Debug.Trace + + + +-- | Hairpin structures. Hairpins with less than 3 unpaired nucleotides are +-- forbidden. +-- +-- NOTE @(xs,i,j)@ is indeed *only* the unpaired stretch. Hence, the length is +-- @j-i+1@, as given. +-- +-- TODO Activate tabulated hairpin structures. + +hairpinF :: Vienna2004 -> (Nuc :!: Nuc) -> (Primary :!: Int :!: Int) -> (Nuc :!: Nuc) -> Int +hairpinF Vienna2004{..} (l :!: lp) (xs :!: i :!: j) (rp :!: r) = hairpinType where + {-# INLINE hairpinType #-} + hairpinType + -- TODO use sliceEq for tabulated hairpins + {- + | len <= 6, Just v <- find tabulated hairpin + -} + | len < 3 = 999999 + | len == 3 = hairpinL!(Z:.len) + tAU + | len > 31 = hairpinL!(Z:.30) + hairpinMM!(Z:.p:.lp:.rp) + llp + | otherwise = hairpinL!(Z:.len) + hairpinMM!(Z:.p:.lp:.rp) + p = mkViennaPair (l,r) + len = j-i+1 + tAU = if p/=vpCG && p/=vpGC then termAU else 0 + llp = floor $ 108.856 * log (fromIntegral len / 30) +{-# INLINE hairpinF #-} + +-- | Tiny loops are small interior loops. This includes canonical stacks +-- without any unpaired nucleotides, and small, tabulated interior loops. + +tinyloopF :: Vienna2004 -> (Primary :!: Int :!: Int) -> Int -> (Primary :!: Int :!: Int) -> Int +tinyloopF Vienna2004{..} (ls :!: li :!: lj) w (rs :!: ri :!: rj) = assert (assertL && assertR) $ loopType where + assertL = let (_,Z:.n) = bounds ls in li>=0 && lj<=n + assertR = let (_,Z:.n) = bounds rs in ri>=0 && rj<=n + {-# INLINE loopType #-} + loopType + | dl==0 || dr==0 = error $ "bug in tinyloop: " ++ show (li,lj,[ls!(Z:.k) | k<-[li..lj]],ri,rj,[rs!(Z:.l) | l<-[ri..rj]]) + -- normal stack + | dl==1 && dr==1 = w+stack!(Z:.op:.ip) + -- one intervening unpaired nucleotide doesn't break the stack + | dl==1 && dr==2 = w+stack!(Z:.op:.ip)+bulgeL!(Z:.1) + | dl==2 && dr==1 = w+stack!(Z:.op:.ip)+bulgeL!(Z:.1) + -- 1x1 symmetric interior loop + | dl==2 && dr==2 = w+iloop1x1!(Z:.op:.ip:.bI:.bJ) + -- 1x2 interior loop + | dl==2 && dr==3 = w+iloop2x1!(Z:.op:.ip:.bI:.bL:.bJ) + -- 2x1 interior loop + | dl==3 && dr==2 = w+iloop2x1!(Z:.op:.ip:.bJ:.bI:.bK) + -- 2x2 interior loop + | dl==3 && dr==3 = w+iloop2x2!(Z:.op:.ip:.bI:.bK:.bL:.bJ) + -- 2x3 interior loops with specialized mismatches + | dl==3 && dr==4 + || dl==4 && dr==3 = w+iloop2x3MM!(Z:.op:.bI:.bJ) + iloop2x3MM!(Z:.ip:.bL:.bK) + iloopL!(Z:.5) + ninio + | otherwise = 999999 -- overlaps with big interior loops + op = mkViennaPair (ls!(Z:.li),rs!(Z:.rj)) + ip = mkViennaPair (rs!(Z:.ri),ls!(Z:.lj)) + bI = ls!(Z:.li+1) + bJ = rs!(Z:.rj-1) + bK = ls!(Z:.lj-1) + bL = rs!(Z:.ri+1) + dl = lj-li + dr = rj-ri +{-# INLINE tinyloopF #-} + +-- | A left bulge @(....[[...]])@ which four unpaired nucleotides in the bulge. +-- the left bulge @ls@ will be given six nucleotides (note, @ls@ is the +-- complete input, use @li@ and @lj@ as the first and last included nucleotide +-- index), the two outer ones being for the outer and inner loop. On the right, +-- we have @rp@ and @r@ which are nucleotides. @ls!(Z:.li)@ and @r@ form the +-- outer Vienna pair. @rp@ and @ls!(Z:.lj)@ form the inner pair. + +bulgeLF :: Vienna2004 -> (Primary :!: Int :!: Int) -> Int -> (Nuc :!: Nuc) -> Int +bulgeLF Vienna2004{..} (ls :!: li :!: lj) w (rp :!: r) = assert (lj-li<=30) $ w + tAUlr + lenE + tAUrpcp where + tAUlr = terminalAU termAU lr + tAUrpcp = terminalAU termAU rpcp + lr = mkViennaPair (ls!(Z:.li),r) + rpcp = mkViennaPair (rp,ls!(Z:.lj)) + lenE = bulgeL!(Z:.lj-li-1) +{-# INLINE bulgeLF #-} + +-- | A right bulge @([[...]]....)@. See 'bulgeLF' for how this works. + +bulgeRF :: Vienna2004 -> (Nuc :!: Nuc) -> Int -> (Primary :!: Int :!: Int) -> Int +bulgeRF Vienna2004{..} (l :!: lp) w (rs :!: ri :!: rj) = assert (rj-ri<=30) $ w + tAUlr + lenE + tAUcplp where + tAUlr = terminalAU termAU lr + tAUcplp = terminalAU termAU cplp + lr = mkViennaPair (l,rs!(Z:.rj)) + cplp = mkViennaPair (rs!(Z:.ri),lp) + lenE = bulgeL!(Z:.rj-ri-1) +{-# iNLINE bulgeRF #-} + +-- | An interior loop with @N@ unpaired nucleotides to the left and @1@ +-- unpaired nucleotide to the right. The regions @ls@ and @rs@ each have 2 +-- nucleotides more than are unpaired. These first and last nucleotides form +-- the last paired or first pairs in the stacks around the loop. + +iloopN1F :: Vienna2004 -> (Primary :!: Int :!: Int) -> Int -> (Primary :!: Int :!: Int) -> Int +iloopN1F Vienna2004{..} (ls:!:li:!:lj) w (rs:!:ri:!:rj) + = assert (lj-li-1 <=29 && rj-ri == 2) + $ w + outerMM + lenE + innerMM + owninio where + poLR = mkViennaPair (ls!(Z:.li), rs!(Z:.rj)) + piRL = mkViennaPair (rs!(Z:.ri), ls!(Z:.lj)) + outerMM = iloop1xnMM!(Z:.poLR:.ls!(Z:.li+1):.rs!(Z:.ri+1)) + innerMM = iloop1xnMM!(Z:.piRL:.rs!(Z:.ri+1):.ls!(Z:.lj-1)) + lenE = iloopL!(Z:.(lL+1)) + owninio = min maxNinio (ninio * (lL-1)) + lL = lj-li-1 + lR = rj-ri-1 +{-# INLINE iloopN1F #-} + +-- | 1xN interior loops. + +iloop1NF :: Vienna2004 -> (Primary :!: Int :!: Int) -> Int -> (Primary :!: Int :!: Int) -> Int +iloop1NF Vienna2004{..} (ls :!: li :!: lj) w (rs :!: ri :!: rj) + = assert (lj-li == 2 && rj-ri-1 <=29) + $ w + outerMM + lenE + innerMM + owninio where + poLR = mkViennaPair (ls!(Z:.li), rs!(Z:.rj)) + piRL = mkViennaPair (rs!(Z:.ri), ls!(Z:.lj)) + outerMM = iloop1xnMM!(Z:.poLR:.ls!(Z:.li+1):.rs!(Z:.rj-1)) + innerMM = iloop1xnMM!(Z:.piRL:.rs!(Z:.ri+1):.ls!(Z:.li+1)) + lenE = iloopL!(Z:.(lR+1)) + owninio = min maxNinio (ninio * (lR-1)) + lL = lj-li-1 + lR = rj-ri-1 +{-# INLINE iloop1NF #-} + +iloopIF :: Vienna2004 -> (Primary :!: Int :!: Int) -> Int -> (Primary :!: Int :!: Int) -> Int +iloopIF Vienna2004{..} (ls :!: li :!: lj) w (rs :!: ri :!: rj) = w + iloopMM!(Z:.p:.bR:.bL) + iloopL!(Z:.len) + owninio where + p = mkViennaPair (rs!(Z:.ri), ls!(Z:.lj)) + bL = ls!(Z:.lj-1) + bR = rs!(Z:.ri+1) + len = (lj-li)+(rj-ri) + owninio = min maxNinio (ninio * (abs $ (lj-li) - (rj-ri))) +{-# INLINE iloopIF #-} + +-- | + +iloopOF :: Vienna2004 -> (Nuc :!: Nuc) -> Int -> (Nuc :!: Nuc) -> Int +iloopOF Vienna2004{..} (l :!: lp) iloopif (rp :!: r) = {- CORE "iloopOF" -} iloopif + iloopMM!(Z:.op:.lp:.rp) + where op = mkViennaPair (l,r) +{-# INLINE iloopOF #-} + +-- | + +multiIF :: Vienna2004 -> Int -> Int -> Int +multiIF Vienna2004{..} b c = b+c +{-# INLINE multiIF #-} + +-- | + +multiOF :: Vienna2004 -> (Nuc :!: Nuc) -> Int -> (Nuc :!: Nuc) -> Int +multiOF Vienna2004{..} (l :!: lp) multiif (rp :!: r) = multiif + multiMM!(Z:.p:.lp:.rp) + multiHelix + multiOffset + terminalAU termAU p + where p = mkViennaPair (l,r) +{-# INLINE multiOF #-} + +-- | + +regionStemF :: Vienna2004 -> Nuc -> Int -> Int +regionStemF Vienna2004{..} _ w = w + multiNuc +{-# INLINE regionStemF #-} + +-- | + +justStemF :: Vienna2004 -> (Nuc :!: Nuc) -> Int -> (Nuc :!: Nuc) -> Int +justStemF Vienna2004{..} (l :!: lp) w (rp :!: r) = w + multiMM!(Z:.p:.l:.r) + where p = mkViennaPair (lp,rp) +{-# INLINE justStemF #-} + +-- | + +bsF :: Int -> Int -> Int +bsF b reg = b +{-# INLINE bsF #-} + +-- | + +rSF :: Nuc -> Int -> Int +rSF nuc w = w +{-# INLINE rSF #-} + +-- | + +bcF :: Int -> Int -> Int +bcF b c = b + c +{-# INLINE bcF #-} + +-- | + +ssF :: Int -> Int +ssF reg = 0 +{-# INLINE ssF #-} + +-- | + +cmF :: Int -> Int -> Int +cmF w s = w+s +{-# INLINE cmF #-} + +-- | + +nilF :: Bool -> Int +nilF e = if e then 0 else 999999 +{-# INLINE nilF #-} + +-- | + +iD = id +{-# INLINE iD #-} + +-- | + +h :: Monad m => S.Stream m Int -> m Int +h = S.foldl' min (999999::Int) +{-# INLINE h #-} + +-- | + +terminalAU termAU p = if p/=vpCG && p/=vpGC then termAU else 0 +{-# INLINE terminalAU #-} + diff --git a/BioInf/RNAfold/Library.hs b/BioInf/RNAfold/Library.hs new file mode 100644 --- /dev/null +++ b/BioInf/RNAfold/Library.hs @@ -0,0 +1,140 @@ +{-# LANGUAGE TypeOperators #-} + +-- | A library of helpers for ADPfusion algorithms. + +module BioInf.RNAfold.Library where + +import Control.Exception (assert) +import Data.Array.Repa.Index +import Debug.Trace +import qualified Data.Vector.Unboxed as VU +import Data.Strict.Tuple hiding (fst,snd) + +import Biobase.Primary +import Biobase.Secondary.Vienna +import Data.PrimitiveArray + +import ADP.Fusion.Monadic +import ADP.Fusion.Monadic.Internal + + + +-- | + +base' :: Primary -> DIM2 -> (Scalar Nuc) +base' inp (Z:.i:.j) = Scalar $ index inp (Z:.i) +{-# INLINE base' #-} + +-- | Nucleotide, second one to the right. The assertion allows a size of one or +-- two to capture special cases of looking outside of the bounds (used by +-- @justStemF@ in RNAfold). + +baseLr' :: Primary -> DIM2 -> (Scalar (Nuc :!: Nuc)) +baseLr' inp (Z:.i:.j) + = assert (let (_,Z:.n) = bounds inp in (i+1==j || i+2==j) && i>=0 && i+1<=n) + . Scalar $ (index inp (Z:.i) :!: index inp (Z:.i+1)) +{-# INLINE baseLr' #-} + +-- | + +baselR' :: Primary -> DIM2 -> (Scalar (Nuc :!: Nuc)) +baselR' inp (Z:.i:.j) + = assert (let (_,Z:.n) = bounds inp in i+1==j && i>0 && i<=n) + . Scalar $ (index inp (Z:.i-1) :!: index inp (Z:.i)) +{-# INLINE baselR' #-} + +-- | + +region' :: VU.Vector Nuc -> DIM2 -> (Scalar (VU.Vector Nuc)) +region' inp (Z:.i:.j) + = assert (let n = VU.length inp -1 in i<=j && i>=0 && j<=n+1) + . Scalar $ VU.unsafeSlice i (j-i) inp +{-# INLINE region' #-} + +-- | A 'Primary' together with the lowest included nucleotide and the highest +-- included nucleotide. + +primary' :: Primary -> DIM2 -> (Scalar (Primary :!: Int :!: Int)) +primary' inp (Z:.i:.j) + = assert (let (_,Z:.u) = bounds inp in i>=0 && j-1<=u && i<=j) + $ Scalar (inp :!: i :!: j-1) +{-# INLINE primary' #-} + +-- | A 'Primary' together with the lowest included nucleotide and the highest +-- included nucleotide. + +primaryPR' :: Primary -> DIM2 -> (Scalar (Primary :!: Int :!: Int)) +primaryPR' inp (Z:.i:.j) + = assert (let (_,Z:.u) = bounds inp in i>=0 && j-2<=u && i<=j) + $ Scalar (inp :!: i :!: j) +{-# INLINE primaryPR' #-} + +-- | A 'Primary' together with the lowest included nucleotide and the highest +-- included nucleotide. + +primaryPL' :: Primary -> DIM2 -> (Scalar (Primary :!: Int :!: Int)) +primaryPL' inp (Z:.i:.j) + = assert (let (_,Z:.u) = bounds inp in i>0 && j-1<=u && i<=j) + $ Scalar (inp :!: i-1 :!: j-1) +{-# INLINE primaryPL' #-} + +-- | Vector of nucleotides peeking one nucleotide to the left. + +regionpl' :: VU.Vector Nuc -> DIM2 -> (Scalar (VU.Vector Nuc)) +regionpl' inp (Z:.i:.j) + = assert (let n = VU.length inp -1 in i-1<=j && i>0 && j<=n+1) + . Scalar $ VU.unsafeSlice (i-1) (j-i+1) inp +{-# INLINE regionpl' #-} + +-- | Vector of nucleotides peeaking one nucleotide to the right. + +regionpr' :: VU.Vector Nuc -> DIM2 -> (Scalar (VU.Vector Nuc)) +regionpr' inp (Z:.i:.j) + = assert (let n = VU.length inp -1 in i<=j && i>=0 && j+1<=n+1) + . Scalar $ VU.unsafeSlice i (j-i+1) inp +{-# INLINE regionpr' #-} + +-- | Tests if (i,j) is a valid base pair. + +basepairing' inp (Z:.i:.j) = tf + where p = mkViennaPair (inp!(Z:.i), inp!(Z:.j-1)) + Z:.n = snd . bounds $ inp + tf = i>=0 && j>0 && i<=j && j-1<=n && i<=n && p /= vpNS +{-# INLINE basepairing' #-} + +stackpairing' inp k (Z:.i:.j) = tf + where ps = [ mkViennaPair (inp!(Z:.i+l), inp!(Z:.j-1-l)) | l <-[0..k-1] ] + Z:.n = snd . bounds $ inp + tf = i>=k-1 && j>k-1 && i+k-1<=j-k+1 && j-k<=n && i+k-1<=n && all (/=vpNS) ps +{-# INLINE stackpairing' #-} + +constrained cns ij = cns ij +{-# INLINE constrained #-} + +-- | + +reglen' :: Primary -> DIM2 -> Scalar Int +reglen' inp (Z:.i:.j) + = assert (let (Z:.o,Z:.n) = bounds inp in i>=0 && (j<=n+1 || n== -1) && i<=j) + . Scalar $ j-i +{-# INLINE reglen' #-} + +-- | + +reglenpl' :: Primary -> DIM2 -> Scalar (Nuc,Nuc,Int) +reglenpl' inp (Z:.i:.j) + = assert (let (_,Z:.n) = bounds inp in i>0 && j<=n && i<=j) + . Scalar $ (index inp (Z:.i-1),index inp (Z:.i),j-i) +{-# INLINE reglenpl' #-} + +reglenpr' :: Primary -> DIM2 -> Scalar (Int,Nuc,Nuc) +reglenpr' inp (Z:.i:.j) + = assert (let (_,Z:.n) = bounds inp in i>=0 && j Scalar Bool +empty (Z:.i:.j) = Scalar $ i==j + diff --git a/BioInf/RNAfold/QuickCheck.hs b/BioInf/RNAfold/QuickCheck.hs new file mode 100644 --- /dev/null +++ b/BioInf/RNAfold/QuickCheck.hs @@ -0,0 +1,31 @@ + +module BioInf.RNAfold.QuickCheck where + +import Test.QuickCheck +import Test.QuickCheck.All + +import ADP.Fusion.QuickCheck.Arbitrary -- hiding (options,customCheck,allProps) + + + +-- * property checking + +fCombined (i,j) = S.toList $ (,,) F.<<< fRegion #~~ fRegion ~~# fRegion F.... id $ Z:.i:.j + +bCombined (i,j) = [ ( (i,k),(k,l),(l,j) ) + | k <- [i..j] + , l <- [k..j] + , k-i >= 3 + , j-l >= 3 + , (k-i) + (j-l) >= 8 + , (k-i) + (j-l) <= 32 + ] + +prop_Combined (Small i, Small j) = fCombined (i,j) == bCombined (i,j) + +options = stdArgs {maxSuccess = 1000} + +customCheck = quickCheckWithResult options + +allProps = $forAllProperties customCheck + 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,102 @@ name: RNAFold -version: 0.0.2.1 -author: Christian Hoener zu Siederdissen (Haskell), Ivo L. Hofacker et al (ViennaRNA) +version: 1.99.1.0 +author: Christian Hoener zu Siederdissen (Haskell), Ivo L. Hofacker et al (ViennaRNA), 2010-2012 +copyright: Christian Hoener zu Siederdissen, 2010-2012 +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.6.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! + If possible, build using the GHC llvm backend, and GHC-7.2.2. + GHC-7.4.x produces very bad code on my system, please benchmark + using 7.2.2. + . + 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. + . + NOTE I am reluctant to call this v2 for now. +Extra-Source-Files: + BioInf/RNAfold/QuickCheck.hs + README.md + +Flag llvm + description: build using llvm backend + default: True + + + library build-depends: - base >=4 && <5, - containers, - vector >=0.7, - primitive >=0.3, + base >= 4 && < 5, + mtl, + strict, + primitive == 0.4.* , + vector == 0.9.* , + PrimitiveArray == 0.2.1.1 , + BiobaseVienna == 0.2.2.3 , + BiobaseXNA == 0.6.2.2 , + ADPfusion == 0.0.1.0 + exposed-modules: + BioInf.RNAfold + BioInf.RNAfold.Combinators + BioInf.RNAfold.Energy + BioInf.RNAfold.Library + ghc-options: + -Odph + -funbox-strict-fields + -fspec-constr + -funfolding-use-threshold100 + -funfolding-keeness-factor100 + if flag (llvm) + ghc-options: + -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 - exposed-modules: - BioInf.RNAFold, - BioInf.RNAFold.Energy, - BioInf.RNAFold.EnergyInt, - BioInf.RNAEval, - BioInf.RNAFold.Functions - if impl(ghc > 6.13) +executable RNAFold + build-depends: + base >= 4 && < 5, + mtl, + strict, + primitive == 0.4.* , + vector == 0.9.* , + PrimitiveArray == 0.2.1.1 , + BiobaseVienna == 0.2.2.3 , + BiobaseXNA == 0.6.2.2 , + ADPfusion == 0.0.1.0 + main-is: + RNAFold.hs + other-modules: + BioInf.RNAfold + BioInf.RNAfold.Combinators + BioInf.RNAfold.Energy + BioInf.RNAfold.Library + ghc-options: + -rtsopts + -Odph + -funbox-strict-fields + -fspec-constr + -funfolding-use-threshold100 + -funfolding-keeness-factor100 + if flag (llvm) 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 + -fllvm -optlo-O3 -optlo-inline -optlo-std-compile-opts --- -fno-method-sharing -- performance does not improve! --- -fdicts-cheap --- -fspec-constr --- -funbox-strict-fields --- -funfolding-use-threshold=100 --- -funfolding-creation-threshold=100 + + +source-repository head + type: git + location: git://github.com/choener/RNAfold + diff --git a/RNAFold.hs b/RNAFold.hs new file mode 100644 --- /dev/null +++ b/RNAFold.hs @@ -0,0 +1,18 @@ + +-- | Simple wrapper around STrnafold + +module Main where + +import BioInf.RNAfold + + + +main = do + xs <- fmap lines getContents + mapM_ doRNAfold xs + +doRNAfold inp = do + let (e,bs) = testRNAfold inp + putStrLn inp + print e + mapM_ putStrLn bs