diff --git a/BioInf/RNAfold.hs b/BioInf/RNAfold.hs deleted file mode 100644 --- a/BioInf/RNAfold.hs +++ /dev/null @@ -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 (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 deleted file mode 100644 --- a/BioInf/RNAfold/Combinators.hs +++ /dev/null @@ -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 (~~#) #-} - diff --git a/BioInf/RNAfold/Energy.hs b/BioInf/RNAfold/Energy.hs deleted file mode 100644 --- a/BioInf/RNAfold/Energy.hs +++ /dev/null @@ -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 #-} - diff --git a/BioInf/RNAfold/Library.hs b/BioInf/RNAfold/Library.hs deleted file mode 100644 --- a/BioInf/RNAfold/Library.hs +++ /dev/null @@ -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 Scalar Bool -empty (Z:.i:.j) = Scalar $ i==j - diff --git a/BioInf/RNAfold/QuickCheck.hs b/BioInf/RNAfold/QuickCheck.hs deleted file mode 100644 --- a/BioInf/RNAfold/QuickCheck.hs +++ /dev/null @@ -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 - diff --git a/BioInf/ViennaRNA.hs b/BioInf/ViennaRNA.hs new file mode 100644 --- /dev/null +++ b/BioInf/ViennaRNA.hs @@ -0,0 +1,2 @@ + +module BioInf.ViennaRNA where diff --git a/BioInf/ViennaRNA/Energy.hs b/BioInf/ViennaRNA/Energy.hs new file mode 100644 --- /dev/null +++ b/BioInf/ViennaRNA/Energy.hs @@ -0,0 +1,117 @@ +{-# LANGUAGE PatternGuards #-} +{-# LANGUAGE NoMonomorphismRestriction #-} +{-# LANGUAGE BangPatterns #-} + +module BioInf.ViennaRNA.Energy where + +import Data.Vector.Fusion.Stream.Monadic as SM +import qualified Data.Vector.Unboxed as VU +import Control.Lens +import Data.Array.Repa.Index +import Prelude as P +import qualified Data.Map as M + +import Data.PrimitiveArray as PA hiding ((!)) +import Data.PrimitiveArray.Zero as PA +import qualified Data.PrimitiveArray as PA +import Biobase.Turner +import Biobase.Vienna +import Biobase.Primary + +import BioInf.ViennaRNA.Signature + +import Debug.Trace + + + +mfe :: Monad m => Signature m Deka Deka +mfe = (hairpin,interior,multi,blockStem,blockUnpair,compsBR,compsBC,structW,structCS,structWS,structOpen,h) where + hairpin ener l lp xs rp r + | len <= 6 + , Just e <- (l `VU.cons` xs `VU.snoc` r) `M.lookup` _hairpinLookup ener = e + | len < 3 = huge + | len == 3 = (ener^.hairpinL) VU.! len + tAU + | len < 31 = (ener^.hairpinL) VU.! len + ener^.hairpinMM!(Z:.l:.r:.lp:.rp) + | otherwise = huge + where + !len = VU.length xs + !tAU = if (l,r) == (nC,nG) || (l,r) == (nG,nC) then Deka 0 else ener^.termAU + interior ener l ls li w ri rs r + | lls==0 && lrs==0 -- stack + = w + _stack ener ! (Z:.l:.r:.ri:.li) -- left, right, right inner, left inner + | lls==1 && lrs==0 || lls==0 && lrs==1 -- stack with slip + = w + _stack ener ! (Z:.l:.r:.ri:.li) + _bulgeL ener VU.! 1 + | lls==1 && lrs==1 + = w + _iloop1x1 ener ! (Z:.l:.r:.ri:.li:.lH:.rL) + | lls==1 && lrs==2 + = w + _iloop2x1 ener ! (Z:.l:.r:.ri:.li:.lH:.rH:.rL) + | lls==2 && lrs==1 + = w + _iloop2x1 ener ! (Z:.l:.r:.ri:.li:.rH:.lH:.lL) + | lls==2 && lrs==2 + = w + _iloop2x2 ener ! (Z:.l:.r:.ri:.li:.lH:.lL:.rH:.rL) + | min lls lrs == 2 && max lls lrs == 3 + = w + _iloop2x3MM ener ! (Z:.l:.r:.lH:.lL) + _iloop2x3MM ener ! (Z:.ri:.li:.rL:.rH) + _iloopL ener VU.! 5 + _ninio ener + | lls==0 && lrs > 1 && lrs <= 30 + = w + tAU + _bulgeL ener VU.! lrs + tUA + | lrs==0 && lls > 1 && lls <= 30 + = w + tAU + _bulgeL ener VU.! lls + tUA + | lrs==1 && lls > 2 && lls <= 30 + = w + _iloop1xnMM ener ! (Z:.li:.ri:.lL:.rH) + _iloop1xnMM ener ! (Z:.r:.l:.rL:.lH) + _iloopL ener VU.! lls + min (_ninio ener *. (lls-1)) (_maxNinio ener) + | lls==1 && lrs > 2 && lrs <= 30 + = w + _iloop1xnMM ener ! (Z:.li:.ri:.lL:.rH) + _iloop1xnMM ener ! (Z:.r:.l:.rL:.lH) + _iloopL ener VU.! lrs + min (_ninio ener *. (lrs-1)) (_maxNinio ener) + | lls>0 && lrs>0 && lls+lrs <= 30 -- TODO missing support for length constraints ? + = w + _iloopMM ener ! (Z:.l:.r:.lH:.rL) + _iloopMM ener ! (Z:.ri:.li:.rH:.lL) + _iloopL ener VU.! (lls+lrs) + min (_ninio ener *. (abs $ lls - lrs)) (_maxNinio ener) + | otherwise = huge -- NOTE later on, we should never get this score + where + !lls = VU.length ls + !lrs = VU.length rs + !tAU = if (l,r) `P.elem` [(nC,nG), (nG,nC)] then Deka 0 else ener^.termAU + !tUA = if (li,ri) `P.elem` [(nC,nG), (nG,nC)] then Deka 0 else ener^.termAU + lH = VU.unsafeHead ls + lL = VU.unsafeLast ls + rH = VU.unsafeHead rs + rL = VU.unsafeLast rs + multi ener l li b c ri r + = b + c + _multiMM ener ! (Z:.r:.l:.ri:.li) + _multiHelix ener + _multiOffset ener where + blockStem ener lo l s r ro + = s + _multiMM ener ! (Z:.l:.r:.lo:.ro) + _multiHelix ener + blockUnpair ener c b + = b + _multiNuc ener + compsBR ener b reg + = let Deka nuc = _multiNuc ener in b + (Deka $ nuc * (VU.length reg)) + compsBC ener b c + = b + c + structW ener w + = w + structCS ener c w + = w + structWS ener w s + = w + s + structOpen ener r + = 0 + h = foldl' min huge + {-# INLINE hairpin #-} + {-# INLINE interior #-} + {-# INLINE multi #-} + {-# INLINE blockStem #-} + {-# INLINE blockUnpair #-} + {-# INLINE compsBR #-} + {-# INLINE compsBC #-} + {-# INLINE structW #-} + {-# INLINE structCS #-} + {-# INLINE structWS #-} + {-# INLINE structOpen #-} + {-# INLINE h #-} +{-# INLINE mfe #-} + +huge = Deka 999999 +{-# INLINE huge #-} + +infixl 8 ! +(!) = (PA.!) +{-# INLINE (!) #-} + +(*.) :: Deka -> Int -> Deka +(Deka k) *. n = Deka $ k*n +{-# INLINE (*.) #-} + diff --git a/BioInf/ViennaRNA/Eval.hs b/BioInf/ViennaRNA/Eval.hs new file mode 100644 --- /dev/null +++ b/BioInf/ViennaRNA/Eval.hs @@ -0,0 +1,95 @@ +{-# LANGUAGE PatternGuards #-} + +-- Direct evaluation of the energy of a given structure. The RNAfold-based +-- variant finds the optimal subset of base pairs that conform to the given +-- structure, this algorithm gives the energy of exactly the given structure. + +module BioInf.ViennaRNA.Eval where + +import Data.Vector.Fusion.Util (Id(..)) +import qualified Data.Vector.Unboxed as VU +import Text.Printf + +import Biobase.Primary +import Biobase.Secondary +import Biobase.Secondary.Diagrams +import Biobase.Vienna + +import BioInf.ViennaRNA.Signature +import BioInf.ViennaRNA.Energy + +import Debug.Trace + + + +rnaEval ener s d1s = flatten $ eval mfe ener s d1s + +flatten :: SSTree PairIdx Structure -> (Deka, [String]) +flatten = f where + unDeka (Deka e) = e + f (SSExt l (External e) xs) = + let etot = e + sum (map fst ys) + ys = map f xs + here = printf "External loop: %d" (unDeka e) + in (etot, here : concatMap snd ys) + f (SSTree _ (Hairpin e l us r) [] ) = (e, [printf "Hairpin loop: %d" (unDeka e)]) + f (SSTree _ (Interior e l ls ll rr rs r) [y]) = + let etot = e + fst (f y) + in (etot, printf "Interior loop: %d" (unDeka e) : snd (f y)) + f (SSTree _ (Multi e ll l r rr) ys) = + let etot = e + sum (map (fst . f) ys) + in (etot, printf "Multi loop: %d %s" (unDeka e) (concatMap show [ll,l,r,rr]) : concatMap (snd . f) ys) + {- + f (SSTree p e xs) = + let etot = e + sum (map fst ys) + ys = map f xs + here + | null xs = printf "Hairpin loop: %d" (unDeka e) + | [_] <- xs = printf "Interior loop: %d" (unDeka e) + | otherwise = printf "Multibranched loop: %d" (unDeka e) + in (etot, here : concatMap snd ys) + -} + +data Structure + = External Deka + | Hairpin Deka Nuc Primary Nuc + | Interior Deka Nuc Primary Nuc Nuc Primary Nuc + | Multi Deka Nuc Nuc Nuc Nuc + +eval :: Signature Id Deka Deka -> Vienna2004 -> Primary -> D1Secondary -> SSTree PairIdx Structure +eval efun ener s d1s = annotateWithEnergy t where + t = d1sTree d1s + (hairpin,interior,multi,blockStem,blockUnpair,compsBR,compsBC,structW,structCS,structWS,structOpen,h) = efun + annotateWithEnergy :: SSTree PairIdx () -> SSTree PairIdx Structure + annotateWithEnergy (SSExt l () xs) = SSExt l e (map annotateWithEnergy xs) where + e = External 0 -- TODO sum of all external loop energies + annotateWithEnergy err@(SSTree (i,j) () xs) + -- hairpin + | null xs + = let pri = VU.slice (i+1) (j-i-1) s in SSTree (i,j) (Hairpin (hairpin ener si sii pri jjs sj) si pri sj) [] + -- interior loop + | [SSTree (k,l) () _] <- xs + = let kks = s VU.! k; sll = s VU.! l + e = interior ener si ls kks 0 sll rs sj + ls = VU.slice (i+1) (k-i-1) s + rs = VU.slice (l+1) (j-l-1) s + in SSTree (i,j) (Interior e si ls kks sll rs sj) (map annotateWithEnergy xs) + -- multibranched loop + | otherwise + = let e = multi ener si sii 0 0 jjs sj + + sum (map bStem xs) + + sum (map (\c -> blockUnpair ener c 0) cs) + cs = [] -- TODO all unpaired nucleotides + bStem (SSTree (k,l) () _) = + let kks = s VU.! (k-1) + sk = s VU.! k + sl = s VU.! l + sll = s VU.! (l+1) + in blockStem ener kks sk 0 sl sll + in SSTree (i,j) (Multi e si sii jjs sj) (map annotateWithEnergy xs) + where + si = s VU.! i + sj = s VU.! j + sii = s VU.! (i+1) + jjs = s VU.! (j-1) + diff --git a/BioInf/ViennaRNA/Fold.hs b/BioInf/ViennaRNA/Fold.hs new file mode 100644 --- /dev/null +++ b/BioInf/ViennaRNA/Fold.hs @@ -0,0 +1,192 @@ +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE BangPatterns #-} + +module BioInf.ViennaRNA.Fold where + +import Data.Vector.Fusion.Util (Id (..)) +import Data.Vector.Fusion.Stream.Monadic as SM +import qualified Data.Vector.Unboxed as VU +import Data.Array.Repa.Index +import Control.Monad +import Control.Monad.ST +import System.IO.Unsafe +import Prelude as P hiding (Maybe(..)) +import Data.Strict.Maybe +import Data.Strict.Tuple + +import Biobase.Secondary.Diagrams +import Data.Array.Repa.Index.Subword +import ADP.Fusion +import ADP.Fusion.Table +import Biobase.Vienna +import Biobase.Primary +import Data.PrimitiveArray as PA hiding ((!)) +import Data.PrimitiveArray.Zero as PA + +import BioInf.ViennaRNA.Signature +import BioInf.ViennaRNA.Energy + + + +basepairing :: Primary -> Subword -> Bool +basepairing inp (Subword(i:.j)) = i+1 Subword -> Bool +structureConstrains Nothing !_ = True +structureConstrains !(Just (D1S c)) (Subword (i:.j)) = (i Subword -> Bool +structC inp (Subword(i:.j)) = VU.length inp == j +{-# INLINE structC #-} + +-- TODO need to fix sized regions, then we are good to go -- performance-wise +-- +-- TODO backtracking +-- +-- TODO struct table +-- +-- TODO restrict structs to a linear band + +gRNAfold ener (hairpin,interior,multi,blockStem,blockUnpair,compsBR,compsBC,structW,structCS,structWS,structOpen,h) weak block comps struct cs inp = + ( weak , + hairpin ener <<< c % pr % hr % pl % c ||| + interior ener <<< c % ir % pr % weak % pl % ir % c ||| + multi ener <<< c % pl % block % comps % pl % c `check` (basepairing inp) `check` (structureConstrains cs) ... h + , block , + blockStem ener <<< pl % c % weak % c % pr ||| + blockUnpair ener <<< c % block ... h + , comps , + compsBR ener <<< block % r ||| + compsBC ener <<< block % comps ... h + , struct , +-- structW ener <<< weak ||| -- TODO peak left/right with default ; not needed anymore + structCS ener <<< c % struct ||| + structWS ener <<< weak % struct ||| -- peak here for weak, too + structOpen ener <<< r `check` (structC inp) ... h + ) where c = chr inp + r = region inp + pr = peekR inp + pl = peekL inp + hr = sregion 3 30 inp + ir = sregion 0 20 inp + {-# INLINE c #-} + {-# INLINE r #-} + {-# INLINE pr #-} + {-# INLINE pl #-} + {-# INLINE hr #-} + {-# INLINE ir #-} +{-# INLINE gRNAfold #-} + + + +pretty :: Monad m => Signature m String (SM.Stream m String) +pretty = (hairpin,interior,multi,blockStem,blockUnpair,compsBR,compsBC,structW,structCS,structWS,structOpen,h) where + hairpin _ _ _ r _ _ = "(" P.++ (P.replicate (VU.length r) '.') P.++ ")" + interior _ _ l _ w _ r _ = "(" P.++ (P.replicate (VU.length l) '.') P.++ w P.++ (P.replicate (VU.length r) '.') P.++ ")" + multi _ _ _ b c _ _ = "(" P.++ b P.++ c P.++ ")" + blockStem _ _ _ w _ _ = w + blockUnpair _ _ b = "." P.++ b + compsBR _ b r = b P.++ (P.replicate (VU.length r) '.') + compsBC _ b c = b P.++ c + structW _ w = w + structCS _ _ w = "." P.++ w + structWS _ w s = w P.++ s + structOpen _ r = P.replicate (VU.length r) '.' + h = return . id + +type CombSignature m e b = Signature m (e, m (SM.Stream m b)) (SM.Stream m b) + + +(<**) + :: (Monad m, Eq b, Eq e) -- , Show e, Show (m [b])) + => Signature m e e + -> Signature m b (SM.Stream m b) + -> CombSignature m e b +(<**) f s = (hairpin,interior,multi,blockStem,blockUnpair,compsBR,compsBC,structW,structCS,structWS,structOpen,h) where + (hairpinF,interiorF,multiF,blockStemF,blockUnpairF,compsBRF,compsBCF,structWF,structCSF,structWSF,structOpenF,hF) = f + (hairpinS,interiorS,multiS,blockStemS,blockUnpairS,compsBRS,compsBCS,structWs,structCSS,structWSS,structOpenS,hS) = s + + xs >>>= f = xs >>= return . SM.map f + ccm2 xs ys f = xs >>= \xx -> ys >>= \yy -> return $ SM.concatMap (\x -> SM.map (\y -> f x y) yy) xx + + hairpin ener l lp xs rp r = (hairpinF ener l lp xs rp r, return $ SM.singleton $ hairpinS ener l lp xs rp r) + interior ener l ls li (wF,wS) ri rs r = (interiorF ener l ls li wF ri rs r, wS >>>= \w -> interiorS ener l ls li w ri rs r) + multi ener l li (bF,bS) (cF,cS) ri r = (multiF ener l li bF cF ri r, ccm2 bS cS $ \b c -> multiS ener l li b c ri r) + blockStem ener lo l (sF,sS) r ro = (blockStemF ener lo l sF r ro, sS >>>= \s -> blockStemS ener lo l s r ro) + blockUnpair ener c (bF,bS) = (blockUnpairF ener c bF, bS >>>= \s -> blockUnpairS ener c s) + compsBR ener (bF,bS) reg = (compsBRF ener bF reg, bS >>>= \s -> compsBRS ener s reg) + compsBC ener (bF,bS) (cF,cS) = (compsBCF ener bF cF, ccm2 bS cS $ \b c -> compsBCS ener b c) + structW ener (wF,wS) = (structWF ener wF, wS >>>= \w -> structWs ener w) + structCS ener c (wF,wS) = (structCSF ener c wF, wS >>>= \w -> structCSS ener c w) + structWS ener (wF,wS) (sF,sS) = (structWSF ener wF sF, ccm2 wS sS $ \w s -> structWSS ener w s) + structOpen ener r = (structOpenF ener r, return . SM.singleton $ structOpenS ener r) + h xs = do + hfs <- hF $ SM.map P.fst xs + let phfs = SM.concatMapM P.snd . SM.filter ((hfs==) . P.fst) $ xs + hS phfs + +rnaFoldConstrained :: Vienna2004 -> Primary -> D1Secondary -> (Deka,[String]) +rnaFoldConstrained ener inp s = (struct ! (Z:.subword 0 n), bt) where + (_,Z:.Subword (_:.n)) = bounds weak + len = VU.length inp + (weak,block,comps,struct) = unsafePerformIO (rnaFoldFill ener (Just s) inp) + bt = backtrack ener (Just s) inp (weak,block,comps,struct) +{-# NOINLINE rnaFoldConstrained #-} + +rnaFold :: Vienna2004 -> Primary -> (Deka,[String]) +rnaFold ener inp = (struct ! (Z:.subword 0 n), bt) where + (_,Z:.Subword (_:.n)) = bounds weak + len = VU.length inp + (weak,block,comps,struct) = unsafePerformIO (rnaFoldFill ener Nothing inp) + bt = backtrack ener Nothing inp (weak,block,comps,struct) +{-# NOINLINE rnaFold #-} + +rnaFoldFill :: Vienna2004 -> Maybe (D1Secondary) -> Primary -> IO (PA.Unboxed (Z:.Subword) Deka, PA.Unboxed (Z:.Subword) Deka, PA.Unboxed (Z:.Subword) Deka, PA.Unboxed (Z:.Subword) Deka) +rnaFoldFill !ener !cs !inp = do + let n = VU.length inp + !weak' <- newWithM (Z:.subword 0 0) (Z:.subword 0 n) huge + !block' <- newWithM (Z:.subword 0 0) (Z:.subword 0 n) huge + !comps' <- newWithM (Z:.subword 0 0) (Z:.subword 0 n) huge + !struc' <- newWithM (Z:.subword 0 0) (Z:.subword 0 n) 0 + fillTables $ gRNAfold ener mfe (mTblSw NonEmptyT weak') (mTblSw NonEmptyT block') (mTblSw NonEmptyT comps') (mTblSw NonEmptyT struc') cs inp + weakF <- freeze weak' + blockF <- freeze block' + compsF <- freeze comps' + strucF <- freeze struc' + return (weakF,blockF,compsF,strucF) +{-# NOINLINE rnaFoldFill #-} + +fillTables (MTbl _ weak, weakF, MTbl _ block, blockF, MTbl _ comps, compsF, MTbl _ struc, strucF) = do + let (_,Z:.Subword (0:.n)) = boundsM weak + forM_ [n,n-1..0] $ \i -> forM_ [i..n] $ \j -> do + weakF (subword i j) >>= writeM weak (Z:.subword i j) + blockF (subword i j) >>= writeM block (Z:.subword i j) + compsF (subword i j) >>= writeM comps (Z:.subword i j) + strucF (subword i j) >>= writeM struc (Z:.subword i j) +{-# INLINE fillTables #-} + +-- * backtracking + +backtrack ener cs (inp :: Primary) (weak :: PA.Unboxed (Z:.Subword) Deka, block :: PA.Unboxed (Z:.Subword) Deka, comps :: PA.Unboxed (Z:.Subword) Deka, struct :: PA.Unboxed (Z:.Subword) Deka) = unId . SM.toList . unId $ sF $ subword 0 n where + n = VU.length inp + w :: SwBtTbl Id Deka String + w = btTbl NonEmptyT weak (wF :: Subword -> Id (SM.Stream Id String)) + b :: SwBtTbl Id Deka String + b = btTbl NonEmptyT block (bF :: Subword -> Id (SM.Stream Id String)) + c :: SwBtTbl Id Deka String + c = btTbl NonEmptyT comps (cF :: Subword -> Id (SM.Stream Id String)) + s :: SwBtTbl Id Deka String + s = btTbl NonEmptyT struct (sF :: Subword -> Id (SM.Stream Id String)) + (_,wF,_,bF,_,cF,_,sF) = gRNAfold ener (mfe <** pretty) w b c s cs inp +{-# INLINE backtrack #-} + diff --git a/BioInf/ViennaRNA/Signature.hs b/BioInf/ViennaRNA/Signature.hs new file mode 100644 --- /dev/null +++ b/BioInf/ViennaRNA/Signature.hs @@ -0,0 +1,38 @@ + +module BioInf.ViennaRNA.Signature where + +import Data.Vector.Fusion.Stream.Monadic as SM + +import Biobase.Primary +import Biobase.Vienna + + + +type Signature m a r = + -- weak / hairpin + ( Vienna2004 -> Nuc -> Nuc -> Primary -> Nuc -> Nuc -> a + -- weak / interior + , Vienna2004 -> Nuc -> Primary -> Nuc -> a -> Nuc -> Primary -> Nuc -> a + -- weak / multibranch + , Vienna2004 -> Nuc -> Nuc -> a -> a -> Nuc -> Nuc -> a + -- block / multistem + , Vienna2004 -> Nuc -> Nuc -> a -> Nuc -> Nuc -> a + -- block / unpaired + , Vienna2004 -> Nuc -> a -> a + -- comps / block region + , Vienna2004 -> a -> Primary -> a + -- comps / block comps + , Vienna2004 -> a -> a -> a + -- struct / weak + , Vienna2004 -> a -> a + -- struct / char-struct + , Vienna2004 -> Nuc -> a -> a + -- struct / weak-struct + , Vienna2004 -> a -> a -> a + -- struct / open + , Vienna2004 -> Primary -> a + -- all / objective + , Stream m a -> m r + ) + + diff --git a/RNAFold.cabal b/RNAFold.cabal --- a/RNAFold.cabal +++ b/RNAFold.cabal @@ -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 diff --git a/RNAFold.hs b/RNAFold.hs deleted file mode 100644 --- a/RNAFold.hs +++ /dev/null @@ -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 diff --git a/src/RNAEval.hs b/src/RNAEval.hs new file mode 100644 --- /dev/null +++ b/src/RNAEval.hs @@ -0,0 +1,91 @@ +{-# LANGUAGE RecordWildCards #-} +{-# LANGUAGE DeriveDataTypeable #-} + +-- | RNAEval tool. + +module Main where + + + +import System.Console.CmdArgs + +import Biobase.Primary +import Biobase.Secondary.Diagrams +import Biobase.Vienna +import qualified Biobase.Turner.Import as TI + +import BioInf.ViennaRNA.Fold +import BioInf.ViennaRNA.Eval + + + +data Options + = Eval + { params :: String + } + | ConstrainedFold + { params :: String + } + deriving (Show,Data,Typeable) + +oEval = Eval + { params = "./params" &= help "Turner 2004 RNA parameters (defaults to ./params)" + } + +oConstrainedFold = ConstrainedFold + { params = "../params" + } + +main = do + o <- cmdArgs $ modes [oEval &= auto, oConstrainedFold] + xs <- fmap lines getContents + tm <- fmap turnerToVienna $ TI.fromDir (params o) "" ".dat" + case o of + Eval{..} -> mapM_ (doEval tm) $ toPairs xs + ConstrainedFold{..} -> mapM_ (doCF tm) $ toPairs xs + +toPairs (x1:x2:xs) = (x1,x2) : toPairs xs +toPairs [x] = error "single last line remaining" +toPairs [] = [] + +doEval tm (inp,str) = do + print $ length inp + print $ rnaEval tm (mkPrimary inp) (mkD1S str) + +doCF tm (inp,str) = do + print $ length inp + print $ rnaFoldConstrained tm (mkPrimary inp) (mkD1S str) + + + +test inp str = do + tm <- fmap turnerToVienna $ TI.fromDir "./params" "" ".dat" + doEval tm (inp,str) + +tests = mapM_ (uncurry test) + [ ( "CCUGACUGGCGUUGACAUAUGGUU" + , ".......(((((......)).)))" + ) + , ( "CUGGGGGUGACAUCCCCCC" + , "..(((((......)).)))" + ) + , ( "GGCGUUGACAUAUGGUU" + , "(((((......)).)))" + ) + , ( "GGGGUUGACAUACCCCC" + , "(((((......)).)))" + ) + , ( "GGCGUUGACAUAUGUU" + , "(((((......)))))" + ) + , ( "GGGGGUGACAUCCCCC" + , "(((((......)))))" + ) + , ( "GGGGGUGACCCCC" + , "(((((...)))))" + ) + , ( "CCUGACUGGCGUUGACAUAUGGUUGCUUGAGCGUAGCCAGGUGUUGGUGGUCCAGUGCAUCAAGGUGCCGUCGGAUCGGAUACUUGGCUUUGCUUAGAUU" + , ".......(((((......)).)))(.(((((((.(((((((.(((((((.....((((......)))))))).)))......))))))).))))))).)." + ) + ] + diff --git a/src/RNAFold.hs b/src/RNAFold.hs new file mode 100644 --- /dev/null +++ b/src/RNAFold.hs @@ -0,0 +1,37 @@ +{-# LANGUAGE DeriveDataTypeable #-} +{-# LANGUAGE RecordWildCards #-} + +-- | Simple wrapper around the rnafold library. + +module Main where + + + +import System.Console.CmdArgs + +import Biobase.Primary +import Biobase.Vienna +import qualified Biobase.Turner.Import as TI + +import BioInf.ViennaRNA.Fold + + + +data Options = Options + { params :: String + } deriving (Show,Data,Typeable) + +options = Options + { params = "./params" &= help "Turner 2004 RNA parameters (defaults to ./params)" + } + +main = do + Options{..} <- cmdArgs options + xs <- fmap lines getContents + tm <- fmap turnerToVienna $ TI.fromDir params "" ".dat" + mapM_ (run' tm) xs + +run' tm inp = do + print $ length inp + print $ rnaFold tm (mkPrimary inp) +