packages feed

RNAwolf 0.3.0.0 → 0.3.0.1

raw patch · 10 files changed

+1178/−1 lines, 10 files

Files

+ BioInf/RNAwolf/Bulge.hs view
@@ -0,0 +1,143 @@+{-# LANGUAGE RecordWildCards #-}++module BioInf.RNAwolf.Bulge where++import qualified Data.Vector.Unboxed as VU++import Biobase.Primary+import Biobase.Secondary+import Data.PrimitiveArray+import Data.PrimitiveArray.Ix++import BioInf.Params+import BioInf.RNAwolf.Types++++-- * Outer part of a bulge++-- | The outer closing pair of a bulge loop (one unpaired region).++fBulgeOuter :: BaseF (NBulgLoop -> ExtFeatures (VU.Vector (PairIdx,Double)))+fBulgeOuter Params{..} inp (NBulgLoop nBulgLoop) i j ct eI eJ+  | i<0 || j>n = error $ "fBulgeOuter: " ++ show (i,j)+  | otherwise = VU.singleton s+  where+    s =+      ( (i,j)+      , nBulgLoop ! (i,j)+      + bulgeClose ! ((nI,nJ),(ct,eI,eJ))+      )+    nI = inp VU.! i+    nJ = inp VU.! j+    n = VU.length inp -1+{-# INLINE fBulgeOuter #-}++-- | Outer part of a normal bulge++btBulgeOuter+  :: Params+  -> Primary+  -> EStem+  -> NBulgLoop+  -> NBT+  -> ExtBT+btBulgeOuter ps inp (EStem eStem) nBulgLoop btBULoop i j ct eI eJ d =+  [ (ij:x,z)+  | i>=0,i<j,j<=n+  , ((k,l),enext) <- VU.toList $ fBulgeOuter ps inp nBulgLoop i j ct eI eJ+  , let d' = newD d ehere enext+--  , trc' ("btESbulge",ij,d') $ testD d'+  , testD d'+  , (x,z) <- btBULoop k l d'+  ] where+      ij = ((i,j),(ct,eI,eJ))+      ehere = eStem!ij+      n = VU.length inp -1++++-- * Loop part of a bulge++-- | The loop-part of bulges. Increases speed by 2x++fBulgeLoop :: BaseF (NBulg -> Features (VU.Vector (PairIdx,Double)))+fBulgeLoop Params{..} inp (NBulg nBulg) i j+  | j-i<=6 = VU.empty+  | otherwise = VU.map f kls+  where+    f (k,l) =+      ( (k,l)+      , nBulg ! (k,l)+      + bulgeLength ! (max (k-i-1) (j-l-1))+      )+    kls = ks VU.++ ls+    ks = VU.fromList [ (k,l)+                     | k <- takeWhile (\k -> k-i-1<=maxLength) [i+2 .. j-4], let l = j-1+                     ]+    ls = VU.fromList [ (k,l)+                     | let k = i+1, l <- takeWhile (\l -> j-l-1<=maxLength) [j-2,j-3 .. i+4]+                     ]+{-# INLINE fBulgeLoop #-}++-- | Backtrack the bulge loop part.++btBulgeLoop+  :: Params+  -> Primary+  -> NBulgLoop+  -> NBulg+  -> NBT+  -> NBT+btBulgeLoop ps inp (NBulgLoop nBulgLoop) nBulg btBU i j d =+  [ (x,z)+  | i>=0,i<j,j<=n+  , ((k,l),enext) <- VU.toList $ fBulgeLoop ps inp nBulg i j+  , let d' = newD d ehere enext+  , testD d'+  , (x,z) <- btBU k l d'+  , testD z+  ] where+      ehere = nBulgLoop!(i,j)+      n = VU.length inp -1+++-- * Inner part of a bulge++-- | Inner part of a bulge to speed up bulge calculations++fBulgeInner :: BaseF (EStem -> Features (VU.Vector (ExtPairIdx,Double)))+fBulgeInner Params{..} inp (EStem eStem) i j+  | j-i<2 = VU.empty+  | otherwise = VU.map f kls+  where+    f ijExt@((i,j),(ctIJ,eI,eJ)) =+      ( ijExt+      , eStem ! ijExt+      + bulgeClose ! ((nJ,nI),(ctIJ,eJ,eI))+      ) where nI = inp VU.! i+              nJ = inp VU.! j+    kls = VU.fromList [ ((i,j),(ctIJ,eI,eJ))+                      | eI<-wsh, eJ<-wsh, ctIJ<-citr+                      ]+{-# INLINE fBulgeInner #-}++-- | Backtrack the inner part of a bulge++btBulgeInner+  :: Params+  -> Primary+  -> NBulg+  -> EStem+  -> ExtBT+  -> NBT+btBulgeInner ps inp (NBulg nBulg) eStem btES i j d =+  [ (x,z)+  | i>=0,i<j,j<=n+  , ((_,(eI,eJ,ct)),enext) <- VU.toList $ fBulgeInner ps inp eStem i j+  , let d' = newD d ehere enext+  , testD d'+  , (x,z) <- btES i j eI eJ ct d'+  ] where+      n = VU.length inp -1+      ehere = nBulg!(i,j)
+ BioInf/RNAwolf/Extern.hs view
@@ -0,0 +1,172 @@+{-# LANGUAGE RecordWildCards #-}++-- | External loops are complete substructures, of which zero to many sit on+-- the chain of nucleotides.++module BioInf.RNAwolf.Extern where++import qualified Data.Vector.Unboxed as VU++import Biobase.Primary+import Biobase.Secondary+import Data.PrimitiveArray+import Data.PrimitiveArray.Ix++import BioInf.Params+import BioInf.RNAwolf.Types+++++-- * Unpaired left nucleotide++-- | An external loop with an unpaired nucleotide to the left++fLeftUnpaired :: BaseF (NExtn -> Features (VU.Vector (PairIdx,Double)))+fLeftUnpaired Params{..} inp (NExtn nExtn) i j+  | i<0 || j>n || i>=j = error $ "Extern.fLeftUnpaired: " ++ show (i,j)+  | otherwise = VU.singleton s+  where+    s =+      ( (i+1,j)+      , nExtn ! (i+1,j)+      )+    n = VU.length inp -1+{-# INLINE fLeftUnpaired #-}++-- | Backtracking a structure with an unpaired nucleotide to the left.++btLeftUnpaired+  :: Params+  -> Primary+  -> NExtn+  -> NBT+  -> NBT+btLeftUnpaired ps inp nExtn btE i j d =+  [ (x,z)+  | i>=0,i<j,j<=n+  , ((k,l),enext) <- VU.toList $ fLeftUnpaired ps inp nExtn i j+  , let d' = newD d ehere enext+  , testD d'+  , (x,z) <- btE k l d'+  , testD z+  ] where+      ehere = unNExtn nExtn !(i,j)+      n = VU.length inp -1++++-- * exactly one stem with indices (i,k), i<k<=j++-- | Energy for exactly one stem at (i,k)++fStem :: BaseF (NStem -> Features (VU.Vector (PairIdx,Double)))+fStem Params{..} inp (NStem nStem) i j+  | i<0 || j>n || i>=j = error $ "Extern.fStem: " ++ show (i,j)+  | otherwise = VU.map f $ VU.enumFromN (i+1) (j-i)+  where+    f k = ( (i,k)+          , nStem !(i,k)+          )+    n = VU.length inp -1+{-# INLINE fStem #-}++-- | Backtrack one stem with right index k.++btStem+  :: Params+  -> Primary+  -> NExtn+  -> NStem+  -> NBT+  -> NBT+btStem ps inp nExtn nStem btNStem i j d =+  [ (x,z)+  | i>=0,i<j,j<=n+  , ((_,k),enext) <- VU.toList $ fStem ps inp nStem i j+  , let d' = newD d ehere enext+  , testD d'+  , (x,z) <- btNStem i k d'+  , testD z+  ] where+      ehere = unNExtn nExtn !(i,j)+      n = VU.length inp -1++++-- * The neutral element for energy. This is an unpaired stretch between (i,j)++-- | This one is important as otherwise, some stretches of nucleotides will+-- always have to be paired. (Obviously, I forgot to add this one for a+-- time...)++fOne :: BaseF (Features (VU.Vector (PairIdx,Double)))+fOne Params{..} inp i j+  | i<0 || j>n || i>j = error $ "Extern.fOne: " ++ show (i,j)+  | otherwise = VU.singleton s+  where+    s = ( (i,j)+        , 0+        )+    n = VU.length inp -1+{-# INLINE fOne #-}++btOne+  :: Params+  -> Primary+  -> NExtn+  -> NBT+btOne ps inp nExtn i j d =+  [ (x,z)+  | i>=0,i<=j,j<=n+  , (_,enext) <- VU.toList $ fOne ps inp i j+  , let d' = newD d ehere enext+  , testD d'+  , let x = []+  , let z = d'+  , testD z+  ] where+      ehere = unNExtn nExtn !(i,j)+      n = VU.length inp -1++++-- * Two or more stems in the external structure++-- | External structures with more than one stem have a NStem on the left and+-- an external NExtn structure on the right.++fStems :: BaseF (NStem -> NExtn -> Features (VU.Vector (Int,Double)))+fStems Params{..} inp (NStem nStem) (NExtn nExtn) i j+  | i<0 || j>n || i>j = error $ "Extern.fStems: " ++ show (i,j)+  | otherwise = VU.map f $ VU.enumFromN (i+1) (j-i-1)+  where+    f k = ( k+          , nStem !(i,k) + nExtn !(k+1,j)+          )+    n = VU.length inp -1+{-# INLINE fStems #-}++-- | Backtracking of an external structure with more than one stem++btStems+  :: Params+  -> Primary+  -> NStem+  -> NExtn+  -> NBT+  -> NBT+  -> NBT+btStems ps inp nStem nExtn btNS btE i j d =+  [ (x++y,z)+  | i>=0,i<j,j<=n+  , (k,enext) <- VU.toList $ fStems ps inp nStem nExtn i j+  , let d' = newD d ehere enext+  , testD d'+  , (x,z') <- btNS i k d'+  , testD z'+  , (y,z) <- btE (k+1) j z'+  , testD z+  ] where+      ehere = unNExtn nExtn !(i,j)+      n = VU.length inp -1
+ BioInf/RNAwolf/Hairpin.hs view
@@ -0,0 +1,59 @@+{-# LANGUAGE RecordWildCards #-}++module BioInf.RNAwolf.Hairpin where++import qualified Data.Vector.Unboxed as VU++import Biobase.Primary+import Biobase.Secondary+import Data.PrimitiveArray+import Data.PrimitiveArray.Ix++import BioInf.Params+import BioInf.RNAwolf.Types++-- | A hairpin is a number of 0 or more unpaired nucleotides, enclosed by the+-- nucleotides (i,j) which pair.+--+-- TODO should we allow loops with more than 30 unpaired nucleotides?+--+-- TODO should we allow hairpins with no unpaired nucleotides in the pin? They+-- do occur, but only under special circumstances which we should model+-- differently...+--+-- TODO re-allow IMI++fHairpin :: [Int] -> BaseF (ExtFeatures (VU.Vector (ExtPairIdx,Double)))+fHairpin imi Params{..} inp i j ct eI eJ+--  | checkIMI imi    = VU.singleton (k,interMolInit)+  | j-i<3           = VU.empty+  | j-i-1>maxLength = VU.empty+  | otherwise       = VU.singleton (k,v)+  where+    k = ((i,j),(ct,eI,eJ))+    v = 0+      + hairpinLength ! (j-i-1)+      + hairpinClose  ! (((inp VU.! i,inp VU.! j),(ct,eI,eJ)),inp VU.! (i+1), inp VU.! (j-1))+    checkIMI [] = False+    checkIMI (x:xs) = i<x && j>x || checkIMI xs+--      + pairDistance  ! (j-i-1)+{-# INLINE fHairpin #-}++-- | Backtracking hairpins.++btHairpin+  :: Params+  -> Primary+  -> EStem+  -> ExtBT+btHairpin ps inp (EStem eStem) i j ct eI eJ d =+  [ ([ij],d')+  | i>=0,j-i>=3,j<=n+  , (_,enext) <- VU.toList $ fHairpin imi ps inp i j ct eI eJ+  , let d' = newD d ehere enext+  , testD d'+  ] where+      ij = ((i,j),(ct,eI,eJ))+      ehere = eStem!ij+      n = VU.length inp -1+      imi = map fst . filter ((==nIMI).snd) $ zip [0..] (VU.toList inp)
+ BioInf/RNAwolf/Interior.hs view
@@ -0,0 +1,196 @@++{-# LANGUAGE RecordWildCards #-}++module BioInf.RNAwolf.Interior where++import qualified Data.Vector.Unboxed as VU++import Biobase.Primary+import Biobase.Secondary+import Data.PrimitiveArray+import Data.PrimitiveArray.Ix++import BioInf.Params+import BioInf.RNAwolf.Types++++-- * Outer part++-- | The outer part of an interior loop. Given a certain basepair type, add the+-- cost from the unpaired part.++fInteriorOuter :: BaseF (NInteLoop -> ExtFeatures (VU.Vector (PairIdx,Double)))+fInteriorOuter Params{..} inp (NInteLoop nInteLoop) i j ct eI eJ+  | j-i<4     = VU.empty+  | otherwise = VU.map f $ VU.singleton (i,j)+  where+    f (k,l) = ( (i,j)+              , nInteLoop ! (i,j)+              + ijSc+--              + if j-i-1<=maxDistance then pairDistance ! (j-i-1) else 0+              )+    ijSc = interiorClose ! (((nI,nJ),(ct,eI,eJ)),nIp1,nJm1)+    nI   = inp VU.! i+    nJ   = inp VU.! j+    nIp1 = inp VU.! (i+1)+    nJm1 = inp VU.! (j-1)+{-# INLINE fInteriorOuter #-}++-- |++btInteriorOuter+  :: Params+  -> Primary+  -> EStem+  -> NInteLoop+  -> NBT -- recursive backtracking function for loops+  -> ExtBT+btInteriorOuter ps inp (EStem eStem) nInteLoop btILoop i j ct eI eJ d =+  [ (ij:x,z) -- interior loop+  | i>=0,i<j,j<=n+  , ((k,l),enext) <- VU.toList $ fInteriorOuter ps inp nInteLoop i j ct eI eJ+  , i<k && l<j+  , let d' = newD d ehere enext+  , testD d'+  , (x,z) <- btILoop k l d'+  , testD z+  ] where+      ij = ((i,j),(ct,eI,eJ))+      ehere = eStem!ij+      n = VU.length inp -1++++-- * Loop part++-- | Performs the interior loop calculations between (i,j) "outer" and (k,l)+-- "inner" part. The score based on the unpaired nucleotides is independent of+-- both, the outer and the inner basepair type.+--+-- NOTE / TODO -- fusion enabled for this function (due to it taking 50% of the+-- time), full fusion is still dependent on other factors and needs to be+-- checked (in particular, we still have allocation events)++fInteriorLoop :: BaseF (NInte -> Features (VU.Vector (PairIdx,Double)))+fInteriorLoop Params{..} inp (NInte nInte) i j+  | j-i<4 = VU.empty+  | otherwise = VU.map f kls+  where+    f (k,l) = ( (k,l)+              , nInte ! (k,l)+              + interiorLength ! (lenI+lenJ)+              + interiorAsym ! (abs $ lenI - lenJ)+              ) where lenI = k-i-1; lenJ = j-l-1+    kls  = VU.map (\(dI,dJ) -> (i+dI,j-dJ)) $ fInteriorKLs i j+{-# INLINE fInteriorLoop #-}++-- | Backtrack the unpaired loop region++btInteriorLoop+  :: Params+  -> Primary+  -> NInteLoop+  -> NInte+  -> NBT+  -> NBT+btInteriorLoop ps inp (NInteLoop nInteLoop) nInte btIL i j d =+  [ (x,z)+  | i>=0,i<j,j<=n+  , ((k,l),enext) <- VU.toList $ fInteriorLoop ps inp nInte i j+  , let d' = newD d ehere enext+  , testD d'+  , (x,z) <- btIL k l d'+  , testD z+  ] where+      ehere = nInteLoop!(i,j)+      n = VU.length inp -1++++-- * Inner part++-- | This opens up an interior loop. For each index (i,j) we minimize over all+-- possible basepair types.++fInteriorInner :: BaseF (EStem -> Features (VU.Vector (ExtPairIdx,Double)))+fInteriorInner Params{..} inp (EStem eStem) i j+  | j-i<2 = VU.empty+  | i==0 || j+1==VU.length inp = VU.empty+  | otherwise = VU.map f kls+  where+    f ijExt@((i,j),(ctIJ,eI,eJ)) =+      ( ijExt+      , eStem ! ijExt+      + interiorClose ! (((nJ,nI),(ctIJ,eJ,eI)),nJp1,nIm1)+      ) where nI = inp VU.! i+              nJ = inp VU.! j+              nIm1 = inp VU.! (i-1)+              nJp1 = inp VU.! (j+1)+    kls = VU.fromList [ ((i,j),(ctIJ,eI,eJ))+                      | eI<-wsh, eJ<-wsh, ctIJ<-citr+                      ]+{-# INLINE fInteriorInner #-}++-- | Backtrack from an NInte result to the corresponding EStem parts++btInteriorInner+  :: Params+  -> Primary+  -> NInte+  -> EStem+  -> ExtBT+  -> NBT+btInteriorInner ps inp (NInte nInte) eStem btES i j d =+  [ (x,z)+  | i>=0,i<j,j<=n+  , ((_,(eI,eJ,ct)),enext) <- VU.toList $ fInteriorInner ps inp eStem i j+  , let d' = newD d ehere enext+  , testD d'+  , (x,z) <- btES i j eI eJ ct d'+  , testD z+  ] where+      n = VU.length inp -1+      ehere = nInte!(i,j)++++-- * Helper functions++-- | Since backtracking interior loops is mostly selfcontained, we encapsulate+-- the above three functions -- which we can't do easily with the forward+-- calculations as they actually have to save on runtime.++{-+btInteriorComplete+  :: Params+  -> Primary+  -> EStem+  -> NInteLoop+  -> NInte+  -> ExtBT+  -> ExtBT+btInteriorComplete ps pri eStem nInteLoop nInte btExtStem i j ct eI eJ d =+  btInteriorOuter ps pri eStem btiloop i j ct eI eJ d+  where btiloop i j d = btInteriorLoop  ps pri nInteLoop nInte btinner i j d+        btinner = btInteriorInner ps pri nInte     eStem btExtStem+-}++-- | Given the outer indices (i,j), produces delta_i and delta_j so that+-- i+delta_i and j-delta_j are the inner indices. 'fInteriorKLs' should fuse+-- and should make sure that l-k>=4 is always true (maxd). Furthermore the+-- maximal unpaired length of both sides combined is determined by 'maxLength'.+--+-- TODO better name than 'maxLength'++fInteriorKLs :: Int -> Int -> VU.Vector (Int,Int)+fInteriorKLs i j = didjs where+  didjs = VU.unfoldr mkDs (4,2)+  mkDs (d,s)+    | d>maxd = Nothing+    | s>=d-2 = Just ((d-s,s),(d+1,2))+    | otherwise = Just ((d-s,s),(d,s+1))+  {-# INLINE mkDs #-}+  maxd = min maxLength (j-i-4)+{-# INLINE fInteriorKLs #-}+
+ BioInf/RNAwolf/Multibranched.hs view
@@ -0,0 +1,350 @@+{-# LANGUAGE RecordWildCards #-}++-- | Functions for handling non-triplet multibranched loops.+--+-- TODO We can do the loop-splitting thing again to speed up multibranched+-- closing by x10.++module BioInf.RNAwolf.Multibranched where++import qualified Data.Vector.Unboxed as VU++import Biobase.Primary+import Biobase.Secondary+import Data.PrimitiveArray+import Data.PrimitiveArray.Ix++import BioInf.Params+import BioInf.RNAwolf.Types++++-- * An unpaired nucleotide to the right of an NMbr structure.++-- | Energy for having the rightmost nucleotide (at j) unpaired in NMBr.++fUnpairedRight :: BaseF (NMbr -> Features (VU.Vector (PairIdx,Double)))+fUnpairedRight Params{..} inp (NMbr nMbr) i j+  | i<0 || j>n = error $ "Multibranched.fUnpairedRight: " ++ show (i,j)+  | i==j       = VU.empty+  | otherwise  = VU.singleton s+  where+    s = ( (i,j-1)+        , nMbr !(i,j-1)+        )+    n = VU.length inp -1++-- | Backtrack in NMbr if the nucleotide at j is unpaired.++btUnpairedRight+  :: Params+  -> Primary+  -> NMbr+  -> NBT+  -> NBT+btUnpairedRight ps inp nMbr btM i j d =+  [ (x,z)+  | i>=0,i<j,j<=n+  , ((_,k),enext) <- VU.toList $ fUnpairedRight ps inp nMbr i j+  , let d' = newD d ehere enext+  , testD d'+  , (x,z) <- btM i k d'+  , testD z+  ] where+      ehere = unNMbr nMbr !(i,j)+      n = VU.length inp -1+++-- * An unpaired nucleotide to the right of an NMbr1 structure.++-- | Energy for having the rightmost nucleotide (at j) unpaired in NMBr1.++fUnpairedRight1 :: BaseF (NMbr1 -> Features (VU.Vector (PairIdx,Double)))+fUnpairedRight1 Params{..} inp (NMbr1 nMbr1) i j+  | i<0 || j>n = error $ "Multibranched.fUnpairedRight: " ++ show (i,j)+  | i==j       = VU.empty+  | otherwise  = VU.singleton s+  where+    s = ( (i,j-1)+        , nMbr1 !(i,j-1)+        )+    n = VU.length inp -1++-- | Backtrack NMbr1 if the nucleotide at j is unpaired.++btUnpairedRight1+  :: Params+  -> Primary+  -> NMbr1+  -> NBT+  -> NBT+btUnpairedRight1 ps inp nMbr1 btM1 i j d =+  [ (x,z)+  | i>=0,i<j,j<=n+  , ((_,k),enext) <- VU.toList $ fUnpairedRight1 ps inp nMbr1 i j+  , let d' = newD d ehere enext+  , testD d'+  , (x,z) <- btM1 i k d'+  , testD z+  ] where+      ehere = unNMbr1 nMbr1 !(i,j)+      n = VU.length inp -1++++-- * EStem to Helix start++-- | A multibranched helix (except the closing one). (i,j) are closed by a+-- basepair. Backtracking into the EStem reveals the type of pairing.++fMlHelix :: BaseF (EStem -> Features (VU.Vector (ExtPairIdx,Double)))+fMlHelix Params{..} inp (EStem eStem) i j+  | i==0 || j==VU.length inp -1 = VU.empty -- TODO not required ?!+  | otherwise = VU.map f exts+  where+    f ext@(ct,eI,eJ) =+      ( ijExt+      , eStem ! ijExt+      + mbClose ! (((nJ,nI),(ct,eJ,eI)),inp VU.! (j+1) ,inp VU.! (i-1))+      + multiHelix+      ) where+          ijExt = ((i,j),(ct,eI,eJ))+          nI = inp VU.! i+          nJ = inp VU.! j+    exts = VU.fromList [ (ct,eI,eJ)+                       | j-i>=2, i>0, j+1<VU.length inp+                       , eI<-wsh, eJ<-wsh, ct<-citr+                       ]+{-# INLINE fMlHelix #-}++-- | Backtracks from (i,j) in NMult into the extended-pairing EStem.++btMlHelix+  :: Params+  -> Primary+  -> NMult+  -> EStem+  -> ExtBT+  -> NBT+btMlHelix ps inp (NMult nMult) eStem btES i j d =+  [ (x,z)+  | i>0,i<j,j<n -- correct boundaries since we access elements at (i-1) and (j+1)+  , ((_,(eI,eJ,ct)),enext) <- VU.toList $ fMlHelix ps inp eStem i j+  , let d' = newD d ehere enext+  , testD d'+  , (x,z) <- btES i j eI eJ ct d'+  , testD z+  ] where+      ehere = nMult !(i,j)+      n = VU.length inp -1++++-- * Closes a multibranched loop and produces an extended structure at+-- ((i,j),ext)++-- | Closes a multibranch loop.+--+-- TODO make completely triplet compliant++fMlClose :: BaseF (NMultLoop -> ExtFeatures (VU.Vector (PairIdx,Double)))+fMlClose Params{..} inp (NMultLoop nMultLoop) i j ct eI eJ+  | i>=j      = VU.empty+  | otherwise = VU.singleton s+  where+    s = ( (i,j)+        , nMultLoop !(i,j)+        + mlc+        )+    mlc = 0+        + multiBranched+        + multiHelix+        + mbClose ! ( ((inp VU.! i, inp VU.! j),(ct,eI,eJ))+                    , inp VU.! (i+1)+                    , inp VU.! (j-1)+                    )+  --      + if j-i-1<=maxDistance then pairDistance ! (j-i-1) else 0+{-# INLINE fMlClose #-}++-- | Backtrack from and extended annotation (ij,ext) into the helper table+-- NMultLoop.++btMlClose+  :: Params+  -> Primary+  -> EStem+  -> NMultLoop+  -> NBT+  -> ExtBT+btMlClose ps inp (EStem eStem) nMultLoop btMultLoop i j ct eI eJ d =+  [ (ij:x,z)+  | i>=0,i<j,j<=n+  , (_,enext) <- VU.toList $ fMlClose ps inp nMultLoop i j ct eI eJ+  , let d' = newD d ehere enext+  , testD d'+  , (x,z) <- btMultLoop i j d'+  , testD z+  ] where+      ij = ((i,j),(ct,eI,eJ))+      ehere = eStem!ij+      n = VU.length inp -1++++-- * Multibranched loop helper table++-- | Multibranched loop helper function that combines "at least one stem" with+-- "exactly one stem" but does not add the closing energy from (i,j).++fMlLoop :: BaseF (NMbr -> NMbr1 -> Features (VU.Vector (Int,Double)))+fMlLoop Params{..} inp (NMbr nMbr) (NMbr1 nMbr1) i j = VU.map f ks where+  f k = ( k+        , nMbr!(i+1,k)+        + nMbr1!(k+1,j-1)+        )+  ks  = VU.enumFromN (i+2) (j-i-3) -- == [i+2 .. j-2]+{-# INLINE fMlLoop #-}++-- | Backtracking the multibranched loop.++btMlLoop+  :: Params+  -> Primary+  -> NMultLoop+  -> NMbr+  -> NMbr1+  -> NBT+  -> NBT+  -> NBT+btMlLoop ps inp (NMultLoop nMultLoop) nMbr nMbr1 btM btM1 i j d =+  [ (x++y,z)+  | i>=0,i<j,j<=n+  , (k,enext) <- VU.toList $ fMlLoop ps inp nMbr nMbr1 i j+  , let d' = newD d ehere enext+  , testD d'+  , i+1<k+  , (x,z') <- btM (i+1) k d'+  , testD z'+  , k+1<j-1+  , (y,z)  <- btM1 (k+1) (j-1) z'+  , testD z+  ] where+      ehere = nMultLoop !(i,j)+      n = VU.length inp -1++++-- * Backtracking of a multibranched stem with unpaired nucleotides to the+-- left.++-- | Backtrack a single stem in NMbr, where the stem has zero or more unpaired+-- nucleotides to the left.++fMlStem :: BaseF (NMult -> Features (VU.Vector (Int,Double)))+fMlStem Params{..} inp (NMult nMult) i j = VU.map f ks where+  f k = ( k+        , nMult!(k,j)+        )+  ks  = VU.enumFromN i (j-i-1) -- == [i..j-2]+{-# INLINE fMlStem #-}++-- | Backtrack by trying to find a multilooped helix.++btMlStem+  :: Params+  -> Primary+  -> NMbr+  -> NMult+  -> NBT+  -> NBT+btMlStem ps inp (NMbr nMbr) nMult btMH i j d =+  [ (x,z) -- stem at (k,j)+  | i>=0,i<j,j<=n+  , (k,enext) <- VU.toList $ fMlStem ps inp nMult i j+  , let d' = newD d ehere enext+  , testD d'+  , (x,z) <- btMH k j d'+  , testD z+  ] where+      ehere = nMbr!(i,j)+      n = VU.length inp -1++++-- * Backtracking of at least two stems by finding one or more stems to the+-- left and exactly one stem to the right.++-- | Add a stem to a multibranch table containing already at least one stem.++fMlStems :: BaseF (NMbr -> NMult -> Features (VU.Vector (Int,Double)))+fMlStems Params{..} inp (NMbr nMbr) (NMult nMult) i j = VU.map f ks where+  f k = ( k+        , nMbr!(i,k)+        + nMult!(k+1,j)+        )+  ks  = VU.enumFromN (i+2) (j-i-4) -- == [i+2..j-3]+{-# INLINE fMlStems #-}++-- | Backtrack by finding the splitting index between an NMbr composite+-- structure and a single multibranched stem NMult (which can contain unpaired+-- nucleotides to the left).++btMlStems+  :: Params+  -> Primary+  -> NMbr+  -> NMult+  -> NBT+  -> NBT+  -> NBT+btMlStems ps inp nMbr nMult btM btMH i j d =+  [ (x++y,z) -- nMbr ++ nStem+  | i>=0,i<j,j<=n+  , (k,enext) <- VU.toList $ fMlStems ps inp nMbr nMult i j+  , let d' = newD d ehere enext+  , testD d'+  , (x,z') <- btM i k d'+  , testD z'+  , (y,z) <- btMH (k+1) j z'+  , testD z+  ] where+      ehere = unNMbr nMbr !(i,j)+      n = VU.length inp -1++++-- * Backtrack a single stem in NMbr1. This stem is closed at (i,j).++-- | Add a single stem to a multibranch table containing zero stems already.+--+-- TODO this would be equal to mlHelix, unify!++fMl1Stem :: BaseF (NMult -> Features (VU.Vector ((Int,Int),Double)))+fMl1Stem Params{..} inp (NMult nMult) i j = VU.singleton s where+  s = ( (i,j)+      , nMult!(i,j)+      )+{-# INLINE fMl1Stem #-}++-- | Backtrack a single stem closed at (i,j) for NMbr1. Takes the route through+-- NMult which solves for the exact pairtype.++btMl1Stem+  :: Params+  -> Primary+  -> NMbr1+  -> NMult+  -> NBT+  -> NBT+btMl1Stem ps inp (NMbr1 nMbr1) nMult btMH i j d =+  [ (x,z)+  | i>=0,i<j,j<=n+  , ((k,l),enext) <- VU.toList $ fMl1Stem ps inp nMult i j+  , i==k && j==l+  , let d' = newD d ehere enext+  , testD d'+  , (x,z) <- btMH i j d'+  , testD z+  ] where+      ehere = nMbr1!(i,j)+      n = VU.length inp -1
+ BioInf/RNAwolf/Stem.hs view
@@ -0,0 +1,101 @@++{-# LANGUAGE RecordWildCards #-}++module BioInf.RNAwolf.Stem where++import qualified Data.Vector.Unboxed as VU++import Biobase.Primary+import Biobase.Secondary+import Data.PrimitiveArray+import Data.PrimitiveArray.Ix++import BioInf.Params+import BioInf.RNAwolf.Types++++-- * Collapses an extended stem into a normal stem, or allows going back from a+-- normal stem to an extended stem.++-- | A normal stem is created by taking the minimum over all possible basepairs+-- of the extended stem.++fNstem :: BaseF (EStem -> Features (VU.Vector (ExtPairIdx,Double)))+fNstem Params{..} inp (EStem eStem) i j+  | i<0 || j>n = error $ "Stem.fNstem: " ++ show (i,j)+  | j-i<=2      = VU.empty+  | otherwise   = VU.map f $ VU.fromList [ (ct,eI,eJ) | ct<-citr, eI<-wsh, eJ<-wsh ]+  where+    f ext = ( idx+            , eStem !idx+            ) where idx = ((i,j),ext)+    n = VU.length inp -1+{-# INLINE fNstem #-}++-- | Backtrack from a normal stem back into the extended stem.++btNstem+  :: Params+  -> Primary+  -> NStem+  -> EStem+  -> ExtBT+  -> NBT+btNstem ps inp nStem eStem btES i j d =+  [ (x,z)+  | i>=0,i<j,j<=n+  , ct <- citr, eI <- wsh, eJ <- wsh+  , ((_,(ct,eI,eJ)),enext) <- VU.toList $ fNstem ps inp eStem i j+  , let d' = newD d ehere enext+  , testD d'+  , (x,z) <- btES i j ct eI eJ d'+  , testD z+  ] where+     ehere = unNStem nStem !(i,j)+     n = VU.length inp -1++++-- * Stacking with extended basepair information++-- | A stem is extended by another pair. The score contribution is dependent on+-- the previous pair. Note that for score lookup purposes, the inner pair is+-- switched.++fStem :: BaseF (EStem -> ExtFeatures (VU.Vector (ExtPairIdx,Double)))+fStem Params{..} inp (EStem eStem) i j ct eI eJ+  | j-i<3     = VU.empty+  | otherwise = VU.map f $ eKLs+  where+    f (ctKL,eK,eL) = ( ((k,l),(ctKL,eK,eL))+                     , eStem ! ((k,l),(ctKL,eK,eL))+                     + stem ! (ijExt,((inp VU.! l,inp VU.! k),(ctKL,eL,eK)))+--                     + if j-i-1<=maxDistance then pairDistance ! (j-i-1) else 0+                     )+    ijExt = ((inp VU.! i, inp VU.! j),(ct,eI,eJ))+    k = i+1+    l = j-1+    eKLs = VU.fromList [ (ctKL,eK,eL) | eK<-wsh, eL<-wsh, ctKL<-citr ]+    {-# NOINLINE eKLs #-} -- speed-up for small input sizes+{-# INLINE fStem #-}++-- | Stem backtracking.++btStem+  :: Params+  -> Primary+  -> EStem+  -> ExtBT -- recursive backtracking function+  -> ExtBT+btStem ps inp eStem btES i j ct eI eJ d =+  [ (ij:x,z)+  | i>=0,i<j,j<=n+  , (((k,l),(ctKL,eK,eL)),enext) <- VU.toList $ fStem ps inp eStem i j ct eI eJ+  , let d' = newD d ehere enext+  , testD d'+  , (x,z) <- btES k l ctKL eK eL d'+  ] where+      ij = ((i,j),(ct,eI,eJ))+      ehere = unEStem eStem !ij+      n = VU.length inp -1
+ BioInf/RNAwolf/TripletBulge.hs view
@@ -0,0 +1,17 @@++{-# LANGUAGE RecordWildCards #-}++module BioInf.RNAwolf.TripletBulge where++import qualified Data.Vector.Unboxed as VU++import Biobase.Primary+import Biobase.Secondary+import Data.PrimitiveArray+import Data.PrimitiveArray.Ix++import BioInf.Params+import BioInf.RNAwolf.Types+++
+ BioInf/RNAwolf/TripletStem.hs view
@@ -0,0 +1,64 @@++{-# LANGUAGE RecordWildCards #-}++module BioInf.RNAwolf.TripletStem where++import qualified Data.Vector.Unboxed as VU++import Biobase.Primary+import Biobase.Secondary+import Data.PrimitiveArray+import Data.PrimitiveArray.Ix++import BioInf.Params+import BioInf.RNAwolf.Types+++{-+      let vTStem    = 999999 -- minimumVU $ fTripletStem ps inp eStem      i j ct eI eJ+-}++{-+-- | Triplet stems have a shared nucleotide.+--   _C+-- A-_B where A is paired with both B and C.+--+-- If A,B is paired with cWW, then A,C can only use the S or H edge.++fTripletStem :: BaseF (ExtTable -> ExtFeatures (VU.Vector (ExtPairIdx,Double)))+fTripletStem Params{..} inp eStem i j ct eI eJ+  | j-i<2     = VU.empty+  | otherwise = {- VU.map fI iShared VU.++ -} VU.map fJ jShared+  where+    fI klExt@((k,l),(ctKL,eK,eL))+      = ( klExt+        , eStem ! klExt+--        + stemTriplet ! ( ((nI,nL),(ctKL,eK,eL)), ((nI,nJ),(ct,eI,eJ)) )+--        + if j-i-1<=maxDistance then pairDistance ! (j-i-1) else 0+        ) where nL = inp VU.! l+    fJ klExt@((k,l),(ctKL,eK,eL))+      = ( klExt+        , eStem ! klExt+        + stemTriplet ! ( ((nI,nJ),((ct,eI,eJ))), ((nJ,nK),(ctKL,eJ,eK)) ) --l==j+--        + if j-i-1<=maxDistance then pairDistance ! (j-i-1) else 0+        ) where nK = inp VU.! k+    nI = inp VU.! i+    nJ = inp VU.! j+    iShared = VU.fromList [ ((k,l),(du,eK,eL))+                          | let k=i+                          , let l=j-1+                          , eK<-[wc,sugar,hoogsteen]+                          , eK/=eI+                          , eL<-[wc,sugar,hoogsteen]+                          , du<-[cis,trans]+                          ]+    jShared = VU.fromList [ ((k,l),(du,eK,eL))+                          | let k=i+1+                          , let l=j+                          , eK<-[wc,sugar,hoogsteen]+                          , eL<-[wc,sugar,hoogsteen]+                          , du<-[cis,trans]+                          , eL/=eJ+                          ]+{- INLINE fTripletStem #-}+-}
+ BioInf/RNAwolf/Types.hs view
@@ -0,0 +1,65 @@++module BioInf.RNAwolf.Types where++import Biobase.Primary+import Biobase.Secondary+import Data.PrimitiveArray++import BioInf.Params++++newD d here next = d - (next - here)+testD d = d>=0++-- | Should really go into BiobaseXNA++wsh = [wc,sugar,hoogsteen]++-- | Should really go into BiobaseXNA++citr = [cis,trans]+++++type ExtBT = Int -> Int -> CTisomerism -> Edge -> Edge -> Double -> BTAnswer+type NBT = Int -> Int -> Double -> BTAnswer++type BTAnswer = [([ExtPairIdx],Double)]++type Table    = PrimArray PairIdx    Double+type ExtTable = PrimArray ExtPairIdx Double++type BaseF a = Params -> Primary -> a+type ExtFeatures a = Int -> Int -> CTisomerism -> Edge -> Edge -> a+type Features a = Int -> Int -> a++type Tables = ( EStem+              , NStem+              , NInte+              , NInteLoop+              , NBulg+              , NBulgLoop+              , NMult+              , NMbr+              , NMbr1+              , NMultLoop+              , NExtn+              )++-- ** Newtype wrappers for all tables.+--+-- NOTE Don't ever not newtype-wrap or you will hurt your brain.++newtype EStem = EStem {unEStem :: ExtTable}+newtype NStem = NStem {unNStem :: Table}+newtype NInte = NInte {unNInte :: Table}+newtype NInteLoop = NInteLoop {unNInteLoop :: Table}+newtype NMult = NMult {unNMult :: Table}+newtype NBulg = NBulg {unNBulg :: Table}+newtype NBulgLoop = NBulgLoop {unBulgLoop :: Table}+newtype NMbr  = NMbr  {unNMbr :: Table}+newtype NMbr1 = NMbr1 {unNMbr1 :: Table}+newtype NExtn = NExtn {unNExtn :: Table}+newtype NMultLoop = NMultLoop {unMultLoop :: Table}
RNAwolf.cabal view
@@ -1,5 +1,5 @@ name:           RNAwolf-version:        0.3.0.0+version:        0.3.0.1 author:         Christian Hoener zu Siederdissen, Stephan H Bernhart, Peter F Stadler, Ivo L Hofacker copyright:      Christian Hoener zu Siederdissen, 2010-2011 homepage:       http://www.tbi.univie.ac.at/software/rnawolf/@@ -73,6 +73,16 @@     BioInf.Params.Import     BioInf.PassiveAggressive     BioInf.RNAwolf+    BioInf.RNAwolf.Bulge+    BioInf.RNAwolf.Extern+    BioInf.RNAwolf.Hairpin+    BioInf.RNAwolf.Interior+    BioInf.RNAwolf.Multibranched+    BioInf.RNAwolf.Stem+    BioInf.RNAwolf.TripletBulge+    BioInf.RNAwolf.TripletStem+    BioInf.RNAwolf.Types+   ghc-options:     -O2   if flag(llvm)