module BioInf.RNAwolf
( rnaWolf
, rnaWolfBacktrack
, rnaWolfOptimal
) where
import Control.Monad
import Control.Monad.ST
import qualified Data.Vector.Unboxed as VU
import Control.Arrow
import Data.PrimitiveArray
import Data.PrimitiveArray.Ix
import Biobase.Primary
import Biobase.Secondary
import BioInf.Params
import BioInf.RNAwolf.Types
import qualified BioInf.RNAwolf.Bulge as Bul
import qualified BioInf.RNAwolf.Extern as Ext
import qualified BioInf.RNAwolf.Hairpin as Hp
import qualified BioInf.RNAwolf.Interior as Int
import qualified BioInf.RNAwolf.Multibranched as Mul
import qualified BioInf.RNAwolf.Stem as Stem
import qualified BioInf.RNAwolf.TripletBulge as TrB
import qualified BioInf.RNAwolf.TripletStem as TrS
import Debug.Trace
rnaWolf :: Params -> Primary -> Tables
rnaWolf ps inp = runST $ foldST ps inp
foldST :: Params -> Primary -> ST s Tables
foldST ps inp = do
let n = VU.length inp 1
let imi = map fst . filter ((==nIMI).snd) $ zip [0..] (VU.toList inp)
(eStemM,eStem) <- second EStem `fmap` mkExtTable n
(nStemM,nStem) <- second NStem `fmap` mkTable n
(nInteM,nInte) <- second NInte `fmap` mkTable n
(nMultM,nMult) <- second NMult `fmap` mkTable n
(nBulgM,nBulg) <- second NBulg `fmap` mkTable n
(nMbrM ,nMbr ) <- second NMbr `fmap` mkTable n
(nMbr1M,nMbr1) <- second NMbr1 `fmap` mkTable n
(nExtnM,nExtn) <- second NExtn `fmap` mkTableWith 0 n
(nBulgLoopM,nBulgLoop) <- second NBulgLoop `fmap` mkTable n
(nInteLoopM,nInteLoop) <- second NInteLoop `fmap` mkTable n
(nMultLoopM,nMultLoop) <- second NMultLoop `fmap` mkTable n
forM_ (mkIJ n) $ \(i,j) -> do
writeM nInteLoopM (i,j) . minimumVU $ Int.fInteriorLoop ps inp nInte i j
writeM nBulgLoopM (i,j) . minimumVU $ Bul.fBulgeLoop ps inp nBulg i j
writeM nMultLoopM (i,j) . minimumVU $ Mul.fMlLoop ps inp nMbr nMbr1 i j
forM_ citr $ \ct -> forM_ wsh $ \eI -> forM_ wsh $ \eJ -> do
let vHairpin = minimumVU $ Hp.fHairpin imi ps inp i j ct eI eJ
let vStem = minimumVU $ Stem.fStem ps inp eStem i j ct eI eJ
let vInterior = minimumVU $ Int.fInteriorOuter ps inp nInteLoop i j ct eI eJ
let vMlClose = minimumVU $ Mul.fMlClose ps inp nMultLoop i j ct eI eJ
let vBulge = minimumVU $ Bul.fBulgeOuter ps inp nBulgLoop i j ct eI eJ
writeM eStemM ((i,j),(ct,eI,eJ)) $ minimum [vHairpin,vStem,vInterior,vMlClose,vBulge]
writeM nStemM (i,j) . minimumVU $ Stem.fNstem ps inp eStem i j
writeM nInteM (i,j) . minimumVU $ Int.fInteriorInner ps inp eStem i j
writeM nMultM (i,j) . minimumVU $ Mul.fMlHelix ps inp eStem i j
writeM nBulgM (i,j) . minimumVU $ Bul.fBulgeInner ps inp eStem i j
let vUnpaired = minimumVU $ Mul.fUnpairedRight ps inp nMbr i j
let vStem = minimumVU $ Mul.fMlStem ps inp nMult i j
let vStems = minimumVU $ Mul.fMlStems ps inp nMbr nMult i j
writeM nMbrM (i,j) $ minimum [vUnpaired, vStem, vStems]
let vUnpaired = minimumVU $ Mul.fUnpairedRight1 ps inp nMbr1 i j
let vStem = minimumVU $ Mul.fMl1Stem ps inp nMult i j
writeM nMbr1M (i,j) $ minimum [vUnpaired,vStem]
let j=n
forM_ [n2,n3..0] $ \i -> do
let unp = minimumVU $ Ext.fLeftUnpaired ps inp nExtn i j
let es = minimumVU $ Ext.fStem ps inp nStem i j
let esl = minimumVU $ Ext.fStems ps inp nStem nExtn i j
let one = minimumVU $ Ext.fOne ps inp i j
writeM nExtnM (i,j) $ minimum [unp,esl,es,one]
return ( eStem
, nStem
, nInte
, nInteLoop
, nBulg
, nBulgLoop
, nMult
, nMbr
, nMbr1
, nMultLoop
, nExtn)
rnaWolfBacktrack :: Params -> Primary -> Double -> Tables -> [([ExtPairIdx],Double)]
rnaWolfBacktrack ps inp delta ( estem@(EStem eStem)
, nstem@(NStem nStem)
, ninte@(NInte nInte)
, ninteloop@(NInteLoop nInteLoop)
, nbulg@(NBulg nBulg)
, nbulgloop@(NBulgLoop nBulgLoop)
, nmult@(NMult nMult)
, nmbr@(NMbr nMbr)
, nmbr1@(NMbr1 nMbr1)
, nmultloop@(NMultLoop nMultLoop)
, nextn@(NExtn nExtn)
)
| null res = [([],0)]
| otherwise = let finalScore = nExtn ! (0,n)
in filter ((<=0).snd) . map (second (\z -> finalScore + delta z)) $ res
where
res = btE 0 n delta
btE i j d =
Ext.btOne ps inp nextn i j d ++
Ext.btLeftUnpaired ps inp nextn btE i j d ++
Ext.btStem ps inp nextn nstem btNS i j d ++
Ext.btStems ps inp nstem nextn btNS btE i j d
btNS i j d =
Stem.btNstem ps inp nstem estem btES i j d
btES :: Int -> Int -> CTisomerism -> Edge -> Edge -> Double -> [([ExtPairIdx],Double)]
btES i j ct eI eJ d =
Hp.btHairpin ps inp estem i j ct eI eJ d ++
Int.btInteriorOuter ps inp estem ninteloop btILoop i j ct eI eJ d ++
Bul.btBulgeOuter ps inp estem nbulgloop btBULoop i j ct eI eJ d ++
Mul.btMlClose ps inp estem nmultloop btMultLoop i j ct eI eJ d ++
Stem.btStem ps inp estem btES i j ct eI eJ d
btILoop i j d =
Int.btInteriorLoop ps inp ninteloop ninte btIL i j d
btIL i j d =
Int.btInteriorInner ps inp ninte estem btES i j d
btBULoop i j d =
Bul.btBulgeLoop ps inp nbulgloop nbulg btBU i j d
btBU i j d =
Bul.btBulgeInner ps inp nbulg estem btES i j d
btMH i j d =
Mul.btMlHelix ps inp nmult estem btES i j d
btMultLoop i j d =
Mul.btMlLoop ps inp nmultloop nmbr nmbr1 btM btM1 i j d
btM i j d =
Mul.btUnpairedRight ps inp nmbr btM i j d ++
Mul.btMlStem ps inp nmbr nmult btMH i j d ++
Mul.btMlStems ps inp nmbr nmult btM btMH i j d
btM1 i j d = let ehere = nMbr1!(i,j) in
Mul.btUnpairedRight1 ps inp nmbr1 btM1 i j d ++
Mul.btMl1Stem ps inp nmbr1 nmult btMH i j d
newD d here next = d (next here)
testD d = d>=0
n = VU.length inp 1
epsilon = 0.001
imi = map fst . filter ((==nIMI).snd) $ zip [0..] (VU.toList inp)
rnaWolfOptimal :: Tables -> Double
rnaWolfOptimal ( estem@(EStem eStem)
, nstem@(NStem nStem)
, ninte@(NInte nInte)
, ninteloop@(NInteLoop nInteLoop)
, nbulg@(NBulg nBulg)
, nbulgloop@(NBulgLoop nBulgLoop)
, nmult@(NMult nMult)
, nmbr@(NMbr nMbr)
, nmbr1@(NMbr1 nMbr1)
, nmultloop@(NMultLoop nMultLoop)
, nextn@(NExtn nExtn)
) = nExtn ! (0,n) where n = snd . snd $ bounds nExtn
minimumVU xs = VU.foldl' (\(!acc) (!k,!v) -> min acc v) 999999 xs
mkTable n = mkTableWith 9999999 n
mkTableWith z n = do
tM <- fromAssocsM (0,0) (n,n) z []
t <- unsafeFreezeM tM
return (tM,t)
mkExtTable n = mkExtTableWith 9999999 n
mkExtTableWith z n = do
tM <- fromAssocsM ((0,0),(cis,wc,wc)) ((n,n),(trans,hoogsteen,hoogsteen)) z []
t <- unsafeFreezeM tM
return (tM,t)
mkIJ n = [ (i,j) | d <- [0..n], j<-[n,n1..0], let i=jd, j>=0, i>=0 ]
trc k x = trace (show (k,x)) x
trc' k x = trace (show k) x
trci' c k x = if c then trace (show k) x else x