RNAFold 1.99.1.3 → 1.99.3.4
raw patch · 14 files changed
+627/−891 lines, 14 filesdep +BiobaseTurnerdep +RNAFolddep +cmdargsdep −mtldep ~ADPfusiondep ~BiobaseViennadep ~BiobaseXNAnew-component:exe:RNAEval
Dependencies added: BiobaseTurner, RNAFold, cmdargs, containers, deepseq, lens, repa
Dependencies removed: mtl
Dependency ranges changed: ADPfusion, BiobaseVienna, BiobaseXNA, PrimitiveArray, primitive, vector
Files
- BioInf/RNAfold.hs +0/−344
- BioInf/RNAfold/Combinators.hs +0/−52
- BioInf/RNAfold/Energy.hs +0/−252
- BioInf/RNAfold/Library.hs +0/−140
- BioInf/RNAfold/QuickCheck.hs +0/−31
- BioInf/ViennaRNA.hs +2/−0
- BioInf/ViennaRNA/Energy.hs +117/−0
- BioInf/ViennaRNA/Eval.hs +95/−0
- BioInf/ViennaRNA/Fold.hs +192/−0
- BioInf/ViennaRNA/Signature.hs +38/−0
- RNAFold.cabal +55/−54
- RNAFold.hs +0/−18
- src/RNAEval.hs +91/−0
- src/RNAFold.hs +37/−0
− BioInf/RNAfold.hs
@@ -1,344 +0,0 @@-{-# LANGUAGE ScopedTypeVariables #-}-{-# 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 :: MArr0 s DIM2 Int <- fromAssocsM (Z:.0:.0) (Z:.n:.n) 999999 []- block :: MArr0 s DIM2 Int <- fromAssocsM (Z:.0:.0) (Z:.n:.n) 999999 []- comps :: MArr0 s DIM2 Int <- fromAssocsM (Z:.0:.0) (Z:.n:.n) 999999 []- struct :: MArr0 s DIM2 Int <- fromAssocsM (Z:.0:.0) (Z:.n:.n) 0 []-- let iif = iloopIF ener <<< primary #~~ weak ~~# primary ... hS- {-# INLINE [0] iif #-}-- let mif = multiIF ener <<< block +~+ comps ... hS- {-# INLINE [0] mif #-}-- fillTables- weak (- -- multiOF ener <<< baseLr -~+ (multiIF ener <<< block +~+ comps ... hS) +~- baselR |||- multiOF ener <<< baseLr -~+ mif +~- baselR |||- -- iloopOF ener <<< baseLr -~+ (iloopIF ener <<< primary #~~ weak ~~# primary ... hS) +~- baselR |||- iloopOF ener <<< baseLr -~+ iif +~- 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<n = sgen (Z:.i-1:.j+1)- | otherwise = S.empty-{-# INLINE adjustStream #-}--adjustStreamBT :: Int -> (DIM2 -> P.Stream elm) -> DIM2 -> P.Stream elm-adjustStreamBT !n sgen (Z:.i:.j)- | i>0 && j<n = sgen (Z:.i-1:.j+1)- | otherwise = P.empty-{-# INLINE adjustStreamBT #-}--infixl 6 .!.-(.!.) stream (h,n) (Z:.i:.j)- | i>0 && j<n = h $ stream (Z:.i-1:.j+1)- | otherwise = h $ S.empty-{-# INLINE (.!.) #-}---- |--infixl 5 `with`-with xs cond ij = if cond ij then xs ij else return 999999-{-# INLINE with #-}--infixl 5 `withBT`-withBT xs cond ij = if cond ij then xs ij else return []---- |--fillTables- :: PrimMonad m- => 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 #-}-
− BioInf/RNAfold/Combinators.hs
@@ -1,52 +0,0 @@-{-# 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 (~~#) #-}-
− BioInf/RNAfold/Energy.hs
@@ -1,252 +0,0 @@-{-# 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 #-}-
− BioInf/RNAfold/Library.hs
@@ -1,140 +0,0 @@-{-# 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<n && i<=j)- . Scalar $ (j-i,index inp (Z:.j-1), index inp (Z:.j))-{-# INLINE reglenpr' #-}---- | True, if the subword at ij is empty.--empty :: DIM2 -> Scalar Bool-empty (Z:.i:.j) = Scalar $ i==j-
− BioInf/RNAfold/QuickCheck.hs
@@ -1,31 +0,0 @@--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-
+ BioInf/ViennaRNA.hs view
@@ -0,0 +1,2 @@++module BioInf.ViennaRNA where
+ BioInf/ViennaRNA/Energy.hs view
@@ -0,0 +1,117 @@+{-# LANGUAGE PatternGuards #-}+{-# LANGUAGE NoMonomorphismRestriction #-}+{-# LANGUAGE BangPatterns #-}++module BioInf.ViennaRNA.Energy where++import Data.Vector.Fusion.Stream.Monadic as SM+import qualified Data.Vector.Unboxed as VU+import Control.Lens+import Data.Array.Repa.Index+import Prelude as P+import qualified Data.Map as M++import Data.PrimitiveArray as PA hiding ((!))+import Data.PrimitiveArray.Zero as PA+import qualified Data.PrimitiveArray as PA+import Biobase.Turner+import Biobase.Vienna+import Biobase.Primary++import BioInf.ViennaRNA.Signature++import Debug.Trace++++mfe :: Monad m => Signature m Deka Deka+mfe = (hairpin,interior,multi,blockStem,blockUnpair,compsBR,compsBC,structW,structCS,structWS,structOpen,h) where+ hairpin ener l lp xs rp r+ | len <= 6+ , Just e <- (l `VU.cons` xs `VU.snoc` r) `M.lookup` _hairpinLookup ener = e+ | len < 3 = huge+ | len == 3 = (ener^.hairpinL) VU.! len + tAU+ | len < 31 = (ener^.hairpinL) VU.! len + ener^.hairpinMM!(Z:.l:.r:.lp:.rp)+ | otherwise = huge+ where+ !len = VU.length xs+ !tAU = if (l,r) == (nC,nG) || (l,r) == (nG,nC) then Deka 0 else ener^.termAU+ interior ener l ls li w ri rs r+ | lls==0 && lrs==0 -- stack+ = w + _stack ener ! (Z:.l:.r:.ri:.li) -- left, right, right inner, left inner+ | lls==1 && lrs==0 || lls==0 && lrs==1 -- stack with slip+ = w + _stack ener ! (Z:.l:.r:.ri:.li) + _bulgeL ener VU.! 1+ | lls==1 && lrs==1+ = w + _iloop1x1 ener ! (Z:.l:.r:.ri:.li:.lH:.rL)+ | lls==1 && lrs==2+ = w + _iloop2x1 ener ! (Z:.l:.r:.ri:.li:.lH:.rH:.rL)+ | lls==2 && lrs==1+ = w + _iloop2x1 ener ! (Z:.l:.r:.ri:.li:.rH:.lH:.lL)+ | lls==2 && lrs==2+ = w + _iloop2x2 ener ! (Z:.l:.r:.ri:.li:.lH:.lL:.rH:.rL)+ | min lls lrs == 2 && max lls lrs == 3+ = w + _iloop2x3MM ener ! (Z:.l:.r:.lH:.lL) + _iloop2x3MM ener ! (Z:.ri:.li:.rL:.rH) + _iloopL ener VU.! 5 + _ninio ener+ | lls==0 && lrs > 1 && lrs <= 30+ = w + tAU + _bulgeL ener VU.! lrs + tUA+ | lrs==0 && lls > 1 && lls <= 30+ = w + tAU + _bulgeL ener VU.! lls + tUA+ | lrs==1 && lls > 2 && lls <= 30+ = w + _iloop1xnMM ener ! (Z:.li:.ri:.lL:.rH) + _iloop1xnMM ener ! (Z:.r:.l:.rL:.lH) + _iloopL ener VU.! lls + min (_ninio ener *. (lls-1)) (_maxNinio ener)+ | lls==1 && lrs > 2 && lrs <= 30+ = w + _iloop1xnMM ener ! (Z:.li:.ri:.lL:.rH) + _iloop1xnMM ener ! (Z:.r:.l:.rL:.lH) + _iloopL ener VU.! lrs + min (_ninio ener *. (lrs-1)) (_maxNinio ener)+ | lls>0 && lrs>0 && lls+lrs <= 30 -- TODO missing support for length constraints ?+ = w + _iloopMM ener ! (Z:.l:.r:.lH:.rL) + _iloopMM ener ! (Z:.ri:.li:.rH:.lL) + _iloopL ener VU.! (lls+lrs) + min (_ninio ener *. (abs $ lls - lrs)) (_maxNinio ener)+ | otherwise = huge -- NOTE later on, we should never get this score+ where+ !lls = VU.length ls+ !lrs = VU.length rs+ !tAU = if (l,r) `P.elem` [(nC,nG), (nG,nC)] then Deka 0 else ener^.termAU+ !tUA = if (li,ri) `P.elem` [(nC,nG), (nG,nC)] then Deka 0 else ener^.termAU+ lH = VU.unsafeHead ls+ lL = VU.unsafeLast ls+ rH = VU.unsafeHead rs+ rL = VU.unsafeLast rs+ multi ener l li b c ri r+ = b + c + _multiMM ener ! (Z:.r:.l:.ri:.li) + _multiHelix ener + _multiOffset ener where+ blockStem ener lo l s r ro+ = s + _multiMM ener ! (Z:.l:.r:.lo:.ro) + _multiHelix ener+ blockUnpair ener c b+ = b + _multiNuc ener+ compsBR ener b reg+ = let Deka nuc = _multiNuc ener in b + (Deka $ nuc * (VU.length reg))+ compsBC ener b c+ = b + c+ structW ener w+ = w+ structCS ener c w+ = w+ structWS ener w s+ = w + s+ structOpen ener r+ = 0+ h = foldl' min huge+ {-# INLINE hairpin #-}+ {-# INLINE interior #-}+ {-# INLINE multi #-}+ {-# INLINE blockStem #-}+ {-# INLINE blockUnpair #-}+ {-# INLINE compsBR #-}+ {-# INLINE compsBC #-}+ {-# INLINE structW #-}+ {-# INLINE structCS #-}+ {-# INLINE structWS #-}+ {-# INLINE structOpen #-}+ {-# INLINE h #-}+{-# INLINE mfe #-}++huge = Deka 999999+{-# INLINE huge #-}++infixl 8 !+(!) = (PA.!)+{-# INLINE (!) #-}++(*.) :: Deka -> Int -> Deka+(Deka k) *. n = Deka $ k*n+{-# INLINE (*.) #-}+
+ BioInf/ViennaRNA/Eval.hs view
@@ -0,0 +1,95 @@+{-# LANGUAGE PatternGuards #-}++-- Direct evaluation of the energy of a given structure. The RNAfold-based+-- variant finds the optimal subset of base pairs that conform to the given+-- structure, this algorithm gives the energy of exactly the given structure.++module BioInf.ViennaRNA.Eval where++import Data.Vector.Fusion.Util (Id(..))+import qualified Data.Vector.Unboxed as VU+import Text.Printf++import Biobase.Primary+import Biobase.Secondary+import Biobase.Secondary.Diagrams+import Biobase.Vienna++import BioInf.ViennaRNA.Signature+import BioInf.ViennaRNA.Energy++import Debug.Trace++++rnaEval ener s d1s = flatten $ eval mfe ener s d1s++flatten :: SSTree PairIdx Structure -> (Deka, [String])+flatten = f where+ unDeka (Deka e) = e+ f (SSExt l (External e) xs) =+ let etot = e + sum (map fst ys)+ ys = map f xs+ here = printf "External loop: %d" (unDeka e)+ in (etot, here : concatMap snd ys)+ f (SSTree _ (Hairpin e l us r) [] ) = (e, [printf "Hairpin loop: %d" (unDeka e)])+ f (SSTree _ (Interior e l ls ll rr rs r) [y]) =+ let etot = e + fst (f y)+ in (etot, printf "Interior loop: %d" (unDeka e) : snd (f y))+ f (SSTree _ (Multi e ll l r rr) ys) =+ let etot = e + sum (map (fst . f) ys)+ in (etot, printf "Multi loop: %d %s" (unDeka e) (concatMap show [ll,l,r,rr]) : concatMap (snd . f) ys)+ {-+ f (SSTree p e xs) =+ let etot = e + sum (map fst ys)+ ys = map f xs+ here+ | null xs = printf "Hairpin loop: %d" (unDeka e)+ | [_] <- xs = printf "Interior loop: %d" (unDeka e)+ | otherwise = printf "Multibranched loop: %d" (unDeka e)+ in (etot, here : concatMap snd ys)+ -}++data Structure+ = External Deka+ | Hairpin Deka Nuc Primary Nuc+ | Interior Deka Nuc Primary Nuc Nuc Primary Nuc+ | Multi Deka Nuc Nuc Nuc Nuc++eval :: Signature Id Deka Deka -> Vienna2004 -> Primary -> D1Secondary -> SSTree PairIdx Structure+eval efun ener s d1s = annotateWithEnergy t where+ t = d1sTree d1s+ (hairpin,interior,multi,blockStem,blockUnpair,compsBR,compsBC,structW,structCS,structWS,structOpen,h) = efun+ annotateWithEnergy :: SSTree PairIdx () -> SSTree PairIdx Structure+ annotateWithEnergy (SSExt l () xs) = SSExt l e (map annotateWithEnergy xs) where+ e = External 0 -- TODO sum of all external loop energies+ annotateWithEnergy err@(SSTree (i,j) () xs)+ -- hairpin+ | null xs+ = let pri = VU.slice (i+1) (j-i-1) s in SSTree (i,j) (Hairpin (hairpin ener si sii pri jjs sj) si pri sj) []+ -- interior loop+ | [SSTree (k,l) () _] <- xs+ = let kks = s VU.! k; sll = s VU.! l+ e = interior ener si ls kks 0 sll rs sj+ ls = VU.slice (i+1) (k-i-1) s+ rs = VU.slice (l+1) (j-l-1) s+ in SSTree (i,j) (Interior e si ls kks sll rs sj) (map annotateWithEnergy xs)+ -- multibranched loop+ | otherwise+ = let e = multi ener si sii 0 0 jjs sj+ + sum (map bStem xs)+ + sum (map (\c -> blockUnpair ener c 0) cs)+ cs = [] -- TODO all unpaired nucleotides+ bStem (SSTree (k,l) () _) =+ let kks = s VU.! (k-1)+ sk = s VU.! k+ sl = s VU.! l+ sll = s VU.! (l+1)+ in blockStem ener kks sk 0 sl sll+ in SSTree (i,j) (Multi e si sii jjs sj) (map annotateWithEnergy xs)+ where+ si = s VU.! i+ sj = s VU.! j+ sii = s VU.! (i+1)+ jjs = s VU.! (j-1)+
+ BioInf/ViennaRNA/Fold.hs view
@@ -0,0 +1,192 @@+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE BangPatterns #-}++module BioInf.ViennaRNA.Fold where++import Data.Vector.Fusion.Util (Id (..))+import Data.Vector.Fusion.Stream.Monadic as SM+import qualified Data.Vector.Unboxed as VU+import Data.Array.Repa.Index+import Control.Monad+import Control.Monad.ST+import System.IO.Unsafe+import Prelude as P hiding (Maybe(..))+import Data.Strict.Maybe+import Data.Strict.Tuple++import Biobase.Secondary.Diagrams+import Data.Array.Repa.Index.Subword+import ADP.Fusion+import ADP.Fusion.Table+import Biobase.Vienna+import Biobase.Primary+import Data.PrimitiveArray as PA hiding ((!))+import Data.PrimitiveArray.Zero as PA++import BioInf.ViennaRNA.Signature+import BioInf.ViennaRNA.Energy++++basepairing :: Primary -> Subword -> Bool+basepairing inp (Subword(i:.j)) = i+1<j && f (inp VU.! i) (inp VU.! (j-1)) where+ f l r = l==nC && r==nG+ || l==nG && r==nC+ || l==nA && r==nU+ || l==nU && r==nA+ || l==nG && r==nU+ || l==nU && r==nG+ {-# INLINE f #-}+{-# INLINE basepairing #-}++structureConstrains :: Maybe D1Secondary -> Subword -> Bool+structureConstrains Nothing !_ = True+structureConstrains !(Just (D1S c)) (Subword (i:.j)) = (i<j) && (VU.unsafeIndex c i == j-1)+{-# INLINE structureConstrains #-}++structC :: Primary -> Subword -> Bool+structC inp (Subword(i:.j)) = VU.length inp == j+{-# INLINE structC #-}++-- TODO need to fix sized regions, then we are good to go -- performance-wise+--+-- TODO backtracking+--+-- TODO struct table+--+-- TODO restrict structs to a linear band++gRNAfold ener (hairpin,interior,multi,blockStem,blockUnpair,compsBR,compsBC,structW,structCS,structWS,structOpen,h) weak block comps struct cs inp =+ ( weak ,+ hairpin ener <<< c % pr % hr % pl % c |||+ interior ener <<< c % ir % pr % weak % pl % ir % c |||+ multi ener <<< c % pl % block % comps % pl % c `check` (basepairing inp) `check` (structureConstrains cs) ... h+ , block ,+ blockStem ener <<< pl % c % weak % c % pr |||+ blockUnpair ener <<< c % block ... h+ , comps ,+ compsBR ener <<< block % r |||+ compsBC ener <<< block % comps ... h+ , struct ,+-- structW ener <<< weak ||| -- TODO peak left/right with default ; not needed anymore+ structCS ener <<< c % struct |||+ structWS ener <<< weak % struct ||| -- peak here for weak, too+ structOpen ener <<< r `check` (structC inp) ... h+ ) where c = chr inp+ r = region inp+ pr = peekR inp+ pl = peekL inp+ hr = sregion 3 30 inp+ ir = sregion 0 20 inp+ {-# INLINE c #-}+ {-# INLINE r #-}+ {-# INLINE pr #-}+ {-# INLINE pl #-}+ {-# INLINE hr #-}+ {-# INLINE ir #-}+{-# INLINE gRNAfold #-}++++pretty :: Monad m => Signature m String (SM.Stream m String)+pretty = (hairpin,interior,multi,blockStem,blockUnpair,compsBR,compsBC,structW,structCS,structWS,structOpen,h) where+ hairpin _ _ _ r _ _ = "(" P.++ (P.replicate (VU.length r) '.') P.++ ")"+ interior _ _ l _ w _ r _ = "(" P.++ (P.replicate (VU.length l) '.') P.++ w P.++ (P.replicate (VU.length r) '.') P.++ ")"+ multi _ _ _ b c _ _ = "(" P.++ b P.++ c P.++ ")"+ blockStem _ _ _ w _ _ = w+ blockUnpair _ _ b = "." P.++ b+ compsBR _ b r = b P.++ (P.replicate (VU.length r) '.')+ compsBC _ b c = b P.++ c+ structW _ w = w+ structCS _ _ w = "." P.++ w+ structWS _ w s = w P.++ s+ structOpen _ r = P.replicate (VU.length r) '.'+ h = return . id++type CombSignature m e b = Signature m (e, m (SM.Stream m b)) (SM.Stream m b)+++(<**)+ :: (Monad m, Eq b, Eq e) -- , Show e, Show (m [b]))+ => Signature m e e+ -> Signature m b (SM.Stream m b)+ -> CombSignature m e b+(<**) f s = (hairpin,interior,multi,blockStem,blockUnpair,compsBR,compsBC,structW,structCS,structWS,structOpen,h) where+ (hairpinF,interiorF,multiF,blockStemF,blockUnpairF,compsBRF,compsBCF,structWF,structCSF,structWSF,structOpenF,hF) = f+ (hairpinS,interiorS,multiS,blockStemS,blockUnpairS,compsBRS,compsBCS,structWs,structCSS,structWSS,structOpenS,hS) = s+ + xs >>>= f = xs >>= return . SM.map f+ ccm2 xs ys f = xs >>= \xx -> ys >>= \yy -> return $ SM.concatMap (\x -> SM.map (\y -> f x y) yy) xx++ hairpin ener l lp xs rp r = (hairpinF ener l lp xs rp r, return $ SM.singleton $ hairpinS ener l lp xs rp r)+ interior ener l ls li (wF,wS) ri rs r = (interiorF ener l ls li wF ri rs r, wS >>>= \w -> interiorS ener l ls li w ri rs r)+ multi ener l li (bF,bS) (cF,cS) ri r = (multiF ener l li bF cF ri r, ccm2 bS cS $ \b c -> multiS ener l li b c ri r)+ blockStem ener lo l (sF,sS) r ro = (blockStemF ener lo l sF r ro, sS >>>= \s -> blockStemS ener lo l s r ro)+ blockUnpair ener c (bF,bS) = (blockUnpairF ener c bF, bS >>>= \s -> blockUnpairS ener c s)+ compsBR ener (bF,bS) reg = (compsBRF ener bF reg, bS >>>= \s -> compsBRS ener s reg)+ compsBC ener (bF,bS) (cF,cS) = (compsBCF ener bF cF, ccm2 bS cS $ \b c -> compsBCS ener b c)+ structW ener (wF,wS) = (structWF ener wF, wS >>>= \w -> structWs ener w)+ structCS ener c (wF,wS) = (structCSF ener c wF, wS >>>= \w -> structCSS ener c w)+ structWS ener (wF,wS) (sF,sS) = (structWSF ener wF sF, ccm2 wS sS $ \w s -> structWSS ener w s)+ structOpen ener r = (structOpenF ener r, return . SM.singleton $ structOpenS ener r)+ h xs = do+ hfs <- hF $ SM.map P.fst xs+ let phfs = SM.concatMapM P.snd . SM.filter ((hfs==) . P.fst) $ xs+ hS phfs++rnaFoldConstrained :: Vienna2004 -> Primary -> D1Secondary -> (Deka,[String])+rnaFoldConstrained ener inp s = (struct ! (Z:.subword 0 n), bt) where+ (_,Z:.Subword (_:.n)) = bounds weak+ len = VU.length inp+ (weak,block,comps,struct) = unsafePerformIO (rnaFoldFill ener (Just s) inp)+ bt = backtrack ener (Just s) inp (weak,block,comps,struct)+{-# NOINLINE rnaFoldConstrained #-}++rnaFold :: Vienna2004 -> Primary -> (Deka,[String])+rnaFold ener inp = (struct ! (Z:.subword 0 n), bt) where+ (_,Z:.Subword (_:.n)) = bounds weak+ len = VU.length inp+ (weak,block,comps,struct) = unsafePerformIO (rnaFoldFill ener Nothing inp)+ bt = backtrack ener Nothing inp (weak,block,comps,struct)+{-# NOINLINE rnaFold #-}++rnaFoldFill :: Vienna2004 -> Maybe (D1Secondary) -> Primary -> IO (PA.Unboxed (Z:.Subword) Deka, PA.Unboxed (Z:.Subword) Deka, PA.Unboxed (Z:.Subword) Deka, PA.Unboxed (Z:.Subword) Deka)+rnaFoldFill !ener !cs !inp = do+ let n = VU.length inp+ !weak' <- newWithM (Z:.subword 0 0) (Z:.subword 0 n) huge+ !block' <- newWithM (Z:.subword 0 0) (Z:.subword 0 n) huge+ !comps' <- newWithM (Z:.subword 0 0) (Z:.subword 0 n) huge+ !struc' <- newWithM (Z:.subword 0 0) (Z:.subword 0 n) 0+ fillTables $ gRNAfold ener mfe (mTblSw NonEmptyT weak') (mTblSw NonEmptyT block') (mTblSw NonEmptyT comps') (mTblSw NonEmptyT struc') cs inp+ weakF <- freeze weak'+ blockF <- freeze block'+ compsF <- freeze comps'+ strucF <- freeze struc'+ return (weakF,blockF,compsF,strucF)+{-# NOINLINE rnaFoldFill #-}++fillTables (MTbl _ weak, weakF, MTbl _ block, blockF, MTbl _ comps, compsF, MTbl _ struc, strucF) = do+ let (_,Z:.Subword (0:.n)) = boundsM weak+ forM_ [n,n-1..0] $ \i -> forM_ [i..n] $ \j -> do+ weakF (subword i j) >>= writeM weak (Z:.subword i j)+ blockF (subword i j) >>= writeM block (Z:.subword i j)+ compsF (subword i j) >>= writeM comps (Z:.subword i j)+ strucF (subword i j) >>= writeM struc (Z:.subword i j)+{-# INLINE fillTables #-}++-- * backtracking++backtrack ener cs (inp :: Primary) (weak :: PA.Unboxed (Z:.Subword) Deka, block :: PA.Unboxed (Z:.Subword) Deka, comps :: PA.Unboxed (Z:.Subword) Deka, struct :: PA.Unboxed (Z:.Subword) Deka) = unId . SM.toList . unId $ sF $ subword 0 n where+ n = VU.length inp+ w :: SwBtTbl Id Deka String+ w = btTbl NonEmptyT weak (wF :: Subword -> Id (SM.Stream Id String))+ b :: SwBtTbl Id Deka String+ b = btTbl NonEmptyT block (bF :: Subword -> Id (SM.Stream Id String))+ c :: SwBtTbl Id Deka String+ c = btTbl NonEmptyT comps (cF :: Subword -> Id (SM.Stream Id String))+ s :: SwBtTbl Id Deka String+ s = btTbl NonEmptyT struct (sF :: Subword -> Id (SM.Stream Id String))+ (_,wF,_,bF,_,cF,_,sF) = gRNAfold ener (mfe <** pretty) w b c s cs inp+{-# INLINE backtrack #-}+
+ BioInf/ViennaRNA/Signature.hs view
@@ -0,0 +1,38 @@++module BioInf.ViennaRNA.Signature where++import Data.Vector.Fusion.Stream.Monadic as SM++import Biobase.Primary+import Biobase.Vienna++++type Signature m a r =+ -- weak / hairpin+ ( Vienna2004 -> Nuc -> Nuc -> Primary -> Nuc -> Nuc -> a+ -- weak / interior+ , Vienna2004 -> Nuc -> Primary -> Nuc -> a -> Nuc -> Primary -> Nuc -> a+ -- weak / multibranch+ , Vienna2004 -> Nuc -> Nuc -> a -> a -> Nuc -> Nuc -> a+ -- block / multistem+ , Vienna2004 -> Nuc -> Nuc -> a -> Nuc -> Nuc -> a+ -- block / unpaired+ , Vienna2004 -> Nuc -> a -> a+ -- comps / block region+ , Vienna2004 -> a -> Primary -> a+ -- comps / block comps+ , Vienna2004 -> a -> a -> a+ -- struct / weak+ , Vienna2004 -> a -> a+ -- struct / char-struct+ , Vienna2004 -> Nuc -> a -> a+ -- struct / weak-struct+ , Vienna2004 -> a -> a -> a+ -- struct / open+ , Vienna2004 -> Primary -> a+ -- all / objective+ , Stream m a -> m r+ )++
RNAFold.cabal view
@@ -1,7 +1,7 @@ name: RNAFold-version: 1.99.1.3-author: Christian Hoener zu Siederdissen (Haskell), Ivo L. Hofacker et al (ViennaRNA), 2010-2012-copyright: Christian Hoener zu Siederdissen, 2010-2012+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 category: Bioinformatics@@ -9,90 +9,91 @@ license-file: LICENSE build-type: Simple stability: experimental-cabal-version: >= 1.6.0+cabal-version: >= 1.8.0 synopsis: RNA secondary structure prediction description: 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. .- 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.+ 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.- .- 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,- mtl >= 2,- strict >= 0.3.2,- primitive == 0.4.* ,- vector == 0.9.* ,- PrimitiveArray == 0.2.2.0 ,- BiobaseVienna == 0.2.2.5 ,- BiobaseXNA == 0.6.2.5 ,- ADPfusion == 0.0.1.2+ 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.RNAfold- BioInf.RNAfold.Combinators- BioInf.RNAfold.Energy- BioInf.RNAfold.Library+ BioInf.ViennaRNA+ BioInf.ViennaRNA.Energy+ BioInf.ViennaRNA.Eval+ BioInf.ViennaRNA.Fold+ BioInf.ViennaRNA.Signature ghc-options:- -O2+ -Odph -funbox-strict-fields- -fspec-constr- -fspec-constr-threshold=500- -fspec-constr-count=10 -funfolding-use-threshold100 -funfolding-keeness-factor100- -fsimplifier-phases10- if flag (llvm)- ghc-options:- -fllvm -optlo-O3 -optlo-inline -optlo-std-compile-opts--+ -fllvm -optlo-O3 -optlo-inline -optlo-std-compile-opts executable RNAFold build-depends:--- we get all build depends from the library dependencies+ base >= 4 && < 5 ,+ cmdargs >= 0.10 ,+ BiobaseTurner >= 0.3 ,+ BiobaseVienna >= 0.3 ,+ BiobaseXNA >= 0.7 ,+ RNAFold main-is: RNAFold.hs- other-modules:- BioInf.RNAfold- BioInf.RNAfold.Combinators- BioInf.RNAfold.Energy- BioInf.RNAfold.Library+ hs-source-dirs:+ src ghc-options: -rtsopts- -O2+ -Odph -funbox-strict-fields- -fspec-constr- -fspec-constr-threshold=500- -fspec-constr-count=10 -funfolding-use-threshold100 -funfolding-keeness-factor100- -fsimplifier-phases10- if flag (llvm)- ghc-options:- -fllvm -optlo-O3 -optlo-inline -optlo-std-compile-opts-+ -fllvm -optlo-O3 -optlo-inline -optlo-std-compile-opts +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 source-repository head type: git
− RNAFold.hs
@@ -1,18 +0,0 @@---- | 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
+ src/RNAEval.hs view
@@ -0,0 +1,91 @@+{-# LANGUAGE RecordWildCards #-}+{-# LANGUAGE DeriveDataTypeable #-}++-- | RNAEval tool.++module Main where++++import System.Console.CmdArgs++import Biobase.Primary+import Biobase.Secondary.Diagrams+import Biobase.Vienna+import qualified Biobase.Turner.Import as TI++import BioInf.ViennaRNA.Fold+import BioInf.ViennaRNA.Eval++++data Options+ = Eval+ { params :: String+ }+ | ConstrainedFold+ { params :: String+ }+ deriving (Show,Data,Typeable)++oEval = Eval+ { params = "./params" &= help "Turner 2004 RNA parameters (defaults to ./params)"+ }++oConstrainedFold = ConstrainedFold+ { params = "../params"+ }++main = do+ o <- cmdArgs $ modes [oEval &= auto, oConstrainedFold]+ xs <- fmap lines getContents+ tm <- fmap turnerToVienna $ TI.fromDir (params o) "" ".dat"+ case o of+ Eval{..} -> mapM_ (doEval tm) $ toPairs xs+ ConstrainedFold{..} -> mapM_ (doCF tm) $ toPairs xs++toPairs (x1:x2:xs) = (x1,x2) : toPairs xs+toPairs [x] = error "single last line remaining"+toPairs [] = []++doEval tm (inp,str) = do+ print $ length inp+ print $ rnaEval tm (mkPrimary inp) (mkD1S str)++doCF tm (inp,str) = do+ print $ length inp+ print $ rnaFoldConstrained tm (mkPrimary inp) (mkD1S str)++++test inp str = do+ tm <- fmap turnerToVienna $ TI.fromDir "./params" "" ".dat"+ doEval tm (inp,str)++tests = mapM_ (uncurry test)+ [ ( "CCUGACUGGCGUUGACAUAUGGUU"+ , ".......(((((......)).)))"+ )+ , ( "CUGGGGGUGACAUCCCCCC"+ , "..(((((......)).)))"+ )+ , ( "GGCGUUGACAUAUGGUU"+ , "(((((......)).)))"+ )+ , ( "GGGGUUGACAUACCCCC"+ , "(((((......)).)))"+ )+ , ( "GGCGUUGACAUAUGUU"+ , "(((((......)))))"+ )+ , ( "GGGGGUGACAUCCCCC"+ , "(((((......)))))"+ )+ , ( "GGGGGUGACCCCC"+ , "(((((...)))))"+ )+ , ( "CCUGACUGGCGUUGACAUAUGGUUGCUUGAGCGUAGCCAGGUGUUGGUGGUCCAGUGCAUCAAGGUGCCGUCGGAUCGGAUACUUGGCUUUGCUUAGAUU"+ , ".......(((((......)).)))(.(((((((.(((((((.(((((((.....((((......)))))))).)))......))))))).))))))).)."+ )+ ]+
+ src/RNAFold.hs view
@@ -0,0 +1,37 @@+{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE RecordWildCards #-}++-- | Simple wrapper around the rnafold library.++module Main where++++import System.Console.CmdArgs++import Biobase.Primary+import Biobase.Vienna+import qualified Biobase.Turner.Import as TI++import BioInf.ViennaRNA.Fold++++data Options = Options+ { params :: String+ } deriving (Show,Data,Typeable)++options = Options+ { params = "./params" &= help "Turner 2004 RNA parameters (defaults to ./params)"+ }++main = do+ Options{..} <- cmdArgs options+ xs <- fmap lines getContents+ tm <- fmap turnerToVienna $ TI.fromDir params "" ".dat"+ mapM_ (run' tm) xs++run' tm inp = do+ print $ length inp+ print $ rnaFold tm (mkPrimary inp)+