RNAwolf 0.3.0.2 → 0.3.1.0
raw patch · 10 files changed
+115/−128 lines, 10 filesdep ~BiobaseTrainingDataPVP ok
version bump matches the API change (PVP)
Dependency ranges changed: BiobaseTrainingData
API changes (from Hackage documentation)
+ BioInf.RNAwolf: rnaWolfOptimal :: Tables -> Double
Files
- BioInf/Keys.hs +4/−48
- BioInf/PassiveAggressive.hs +27/−30
- BioInf/RNAwolf.hs +28/−8
- BioInf/RNAwolf/Bulge.hs +16/−14
- BioInf/RNAwolf/Extern.hs +10/−10
- BioInf/RNAwolf/Hairpin.hs +0/−2
- BioInf/RNAwolf/Interior.hs +11/−7
- BioInf/RNAwolf/Stem.hs +0/−1
- RNAwolf.cabal +8/−4
- RNAwolfTrain.hs +11/−4
BioInf/Keys.hs view
@@ -140,26 +140,6 @@ = [ Stem (((nI,nJ),ijExt),((nL,nK),swap23 klExt)) -- , PairDistance (j-i-1) ] ++ concatMap f xs- {-- -- triplet stem (right nucleotide shared)- | [SSTree ((k,l),klExt) _ _] <- xs- , i+1==k && j==l- , let nK = pri `vuIndex` k- = [ StemTriplet ( ((nI,nJ),ijExt)- , ((nJ,nK),swap23 klExt)- )--- , PairDistance (j-i-1)- ] ++ concatMap f xs- -}- {-- -- triplet stem (left nucleotide shared)- | [SSTree ((k,l),klExt) _ _] <- xs- , i==k && j-1==l- , let nL = pri `vuIndex` l- = [ StemTriplet (((nI,nL),klExt),((nI,nJ),ijExt)) -- shared nuc first (nI), then 5' (nL) first--- , PairDistance (j-i-1)- ] ++ concatMap f xs- -} -- interior loops | [SSTree ((k,l),klExt) _ _] <- xs , let lenI = k-i-1; lenJ = j-l-1; len = lenI+lenJ@@ -167,9 +147,10 @@ , let nL = pri `vuIndex` l; nK = pri `vuIndex` k; lkExt = swap23 klExt , let nLm1 = pri `vuIndex` (l-1); nKp1 = pri `vuIndex` (k+1) , let nLp1 = pri `vuIndex` (l+1); nKm1 = pri `vuIndex` (k-1)- = [ InteriorLength len+-- = [ InteriorClose (((nN,nN),(cis,wc,wc)),nN,nN) -- the closing pair (don't block this)+ = [ InteriorClose (((nI,nJ),ijExt),nIp1,nJm1) -- the closing pair (don't block this)+ , InteriorLength len , InteriorAsym $ abs (lenI-lenJ)- , InteriorClose (((nI,nJ),ijExt),nIp1,nJm1) , InteriorClose (((nL,nK),lkExt),nLp1,nKm1) -- , PairDistance (j-i-1) ] ++ concatMap f xs@@ -181,33 +162,9 @@ , let nK = pri `vuIndex` k; nL = pri `vuIndex` l; lkExt = swap23 klExt = [ BulgeLength len , BulgeClose ((nI,nJ),ijExt)- , BulgeClose ((nL,nK),lkExt)--- , PairDistance (j-i-1)- ] ++ concatMap f xs- {-- -- bulge triplet (left)- | [SSTree ((k,l),klExt) _ _] <- xs- , let lenJ = j-l-1- , i==k && lenJ>0- , lenJ<=P.maxLength- , let nL = pri `vuIndex` l- = [ BulgeLength lenJ- , BulgeTriplet (((nI,nL),klExt),((nI,nJ),ijExt))--- , PairDistance (j-i-1)- ] ++ concatMap f xs- -}- {-- -- bulge triplet (right)- | [SSTree ((k,l),klExt) _ _] <- xs- , let lenI = k-i-1- , j==l && lenI>0- , lenI<=P.maxLength- , let nK = pri `vuIndex` k- = [ BulgeLength lenI- , BulgeTriplet (((nJ,nI),swap23 ijExt),((nJ,nK),swap23 klExt))+-- , BulgeClose ((nL,nK),lkExt) -- , PairDistance (j-i-1) ] ++ concatMap f xs- -} -- close a multibranched loop -- -- TODO what about shared multibranched loops? (see sequence GCGGCACCGUCCGCUCAAACAAACGG in fr3d DB)@@ -229,7 +186,6 @@ ] ++ concatMap f xs | otherwise = concatMap f xs- {- xxx| otherwise = error $ show ("unknown features:", inp , SSTree ((i,j),ijExt) () xs) -} where nI = pri `vuIndex` i nJ = pri `vuIndex` j
BioInf/PassiveAggressive.hs view
@@ -21,6 +21,7 @@ import Data.Set as S import Control.Arrow import Data.Map as M+import Text.Printf import Biobase.TrainingData import BioInf.Keys@@ -41,57 +42,53 @@ defaultPA :: Double -> P.Params -> TrainingData -> (P.Params,Double,Double,[(Int,Double)]) defaultPA aggressiveness params td@TrainingData{..}--- | kScore+0.02 < pScore = error $ show (pScore,kScore,pOnly,kOnly,tau,changes)--- | pScore > kScore = error "foo" | L.null $ pOnly++kOnly = (params,0,1,[]) | sty >= 0.999 = (params,0,1,[])--- | otherwise = error $ show (pOnly,kOnly,kScore,pScore,tau,changes)- | otherwise = ( heck+ | otherwise = ( new , tau , sty , changes ) where- new1 = P.fromList . VU.toList $ VU.accum (\v pm -> v+pm) cur changes- new2 = P.fromList . VU.toList $ VU.accum (\v pm -> v+pm) (VU.fromList $ P.toList new1) []- heck- | P.toList new1 == P.toList new2 = new1- | otherwise = error "fuck" -- ignore this line ;-) (impressive, that you actually read this code!)+ -- create new vector+ new = P.fromList . VU.toList $ VU.accum (\v pm -> v+pm) cur changes pFeatures = featureVector primary predicted kFeatures = featureVector primary secondary pOnly = pFeatures L.\\ kFeatures kOnly = kFeatures L.\\ pFeatures numChanges = genericLength $ pOnly ++ kOnly+ changes = zip kOnly (repeat $ negate tau) ++ zip pOnly (repeat tau) cur = VU.fromList . P.toList $ params pScore = sum . L.map (cur VU.!) $ pFeatures kScore = sum . L.map (cur VU.!) $ kFeatures- pScore2 = sum . L.map (cur VU.!) $ pFeatures- kScore2 = sum . L.map (cur VU.!) $ kFeatures+ -- weight calculation tau- | abs ((kScore2 - pScore2) - (kScore-pScore)) > 0.1- = error $ "abs: \n" ++ z- | val < 0 = error $ "val<0 \n" ++ z- | sty >= 0.999 = 0- | otherwise = val -- 100 * val+ | kScore + epsilon < pScore+ = error $ "S(known) < S(predicted)\n" ++ errorKnownTooGood td cur kFeatures pFeatures+ | sty > 0.999+ && kScore+epsilon < pScore+ = error $ "S(known) < S(predicted)\n" ++ errorKnownTooGood td cur kFeatures pFeatures+ | sty >= 0.999+ = 0+ | otherwise+ = val where val = min aggressiveness $ (kScore - pScore + sqrt (1-sty)) / (numChanges ^ 2)- z = show ( kScore,pScore,kScore - pScore- , kScore2,pScore2, kScore2 - pScore2- ) ++ "\n" ++ primary ++ "\n" ++ (concat $ intersperse "\n" comments) ++ "\n" ++- ( L.concatMap (\x -> show x ++ "\n")- $ L.map (fun &&& (cur VU.!)) kOnly ) ++ " <<<\n" ++- ( L.concatMap (\x -> show x ++ "\n") - $ L.map (fun &&& (cur VU.!)) pOnly ) ++ " ALL\n" ++- ( L.concatMap (\x -> show x ++ "\n")- $ L.map (fun &&& (cur VU.!)) pFeatures)- fun i = let lol = vks M.! i in (lol, fun2 lol)- fun2 hc@(HairpinClose k) = P.hairpinClose params PA.! k- fun2 hl@(HairpinLength l) = P.hairpinLength params PA.! l- fun2 _ = (-1) sty = case fmeasure (mkConfusionMatrix td) of -- currently optimizing using F_1 Left _ -> 1 Right v -> v- changes = zip kOnly (repeat $ negate tau) ++ zip pOnly (repeat tau)+ -- special constants+ epsilon = 0.1++-- | In case that the known structure has a score 'epsilon' better than the+-- predicted, we have an error condition, as this should never be the case.++errorKnownTooGood TrainingData{..} curPs kFeatures pFeatures = z where+ z = printf "S(known) = %7.4f, S(pred) = %7.4f, S(known) - S(pred) = %7.4f\n"+ kScore pScore (kScore - pScore)+ ++ printf "%s\n%s\n" primary (concat $ intersperse "\n" comments)+ kScore = sum . L.map (curPs VU.!) $ kFeatures+ pScore = sum . L.map (curPs VU.!) $ pFeatures -- | Pull in the statistical interface. From the confusion matrix, we -- automagically get everything we need.
BioInf/RNAwolf.hs view
@@ -26,6 +26,7 @@ module BioInf.RNAwolf ( rnaWolf , rnaWolfBacktrack+ , rnaWolfOptimal ) where import Control.Monad@@ -82,6 +83,15 @@ -- required to calculate this way, as otherwise the shared nucleotide -- variants will fail. forM_ (mkIJ n) $ \(i,j) -> do+ -- Fill helper tables that need to be accessed to calculate "weak": interior, bulged+ -- FIXME+ -- fill the interior LOOP table (includes everything except the closing pair)+ writeM nInteLoopM (i,j) . minimumVU $ Int.fInteriorLoop ps inp nInte i j+ -- fill bulge LOOP table+ writeM nBulgLoopM (i,j) . minimumVU $ Bul.fBulgeLoop ps inp nBulg i j+ -- multibranched close helper table (should improve speed for MLs by 2x3x3)+ writeM nMultLoopM (i,j) . minimumVU $ Mul.fMlLoop ps inp nMbr nMbr1 i j+ -- and now, fill the weak table forM_ citr $ \ct -> forM_ wsh $ \eI -> forM_ wsh $ \eJ -> do -- weak table (everything is weak) let vHairpin = minimumVU $ Hp.fHairpin imi ps inp i j ct eI eJ@@ -94,16 +104,10 @@ writeM nStemM (i,j) . minimumVU $ Stem.fNstem ps inp eStem i j -- fill the inner interior table writeM nInteM (i,j) . minimumVU $ Int.fInteriorInner ps inp eStem i j- -- fill the interior LOOP table (includes everything except the closing pair)- writeM nInteLoopM (i,j) . minimumVU $ Int.fInteriorLoop ps inp nInte i j -- fill multibranch helper table writeM nMultM (i,j) . minimumVU $ Mul.fMlHelix ps inp eStem i j- -- multibranched close helper table (should improve speed for MLs by 2x3x3)- writeM nMultLoopM (i,j) . minimumVU $ Mul.fMlLoop ps inp nMbr nMbr1 i j -- fill bulge close helper table writeM nBulgM (i,j) . minimumVU $ Bul.fBulgeInner ps inp eStem i j- -- fill bulge LOOP table- writeM nBulgLoopM (i,j) . minimumVU $ Bul.fBulgeLoop ps inp nBulg i j -- one or more multibranched stems let vUnpaired = minimumVU $ Mul.fUnpairedRight ps inp nMbr i j let vStem = minimumVU $ Mul.fMlStem ps inp nMult i j@@ -162,15 +166,17 @@ , 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)) $ btE 0 n delta+ in filter ((<=0).snd) . map (second (\z -> finalScore + delta -z)) $ res where+ res = btE 0 n delta btE i j d = -- trace (show ("btE",i,j,d)) $ Ext.btOne ps inp nextn i j d ++ -- [1] 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 =+ btNS i j d = -- trace (show ("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 = -- trace (show ("btES",i,j,ct,eI,eJ,d)) $@@ -205,7 +211,21 @@ epsilon = 0.001 imi = map fst . filter ((==nIMI).snd) $ zip [0..] (VU.toList inp) +-- | Return the optimal energy. +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 -- * Helper functions
BioInf/RNAwolf/Bulge.hs view
@@ -21,13 +21,11 @@ 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+ | otherwise = VU.singleton ( (i,j)+ , nBulgLoop ! (i,j)+ + bulgeClose ! ((nI,nJ),(ct,eI,eJ))+ ) 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@@ -47,7 +45,6 @@ | 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@@ -71,7 +68,11 @@ , nBulg ! (k,l) + bulgeLength ! (max (k-i-1) (j-l-1)) )- kls = ks VU.++ ls+ kls = fBulgeLoopIndices i j++-- | Index generator for bulged loops++fBulgeLoopIndices i j = ks VU.++ ls where ks = VU.fromList [ (k,l) | k <- takeWhile (\k -> k-i-1<=maxLength) [i+2 .. j-4], let l = j-1 ]@@ -111,14 +112,15 @@ | 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))+ f ext =+ ( ij+ , eStem ! ij+-- + bulgeClose ! ((nJ,nI),ext) ) where nI = inp VU.! i nJ = inp VU.! j- kls = VU.fromList [ ((i,j),(ctIJ,eI,eJ))- | eI<-wsh, eJ<-wsh, ctIJ<-citr+ ij = ((i,j),ext)+ kls = VU.fromList [ (ct,eI,eJ)+ | eI<-wsh, eJ<-wsh, ct<-citr ] {-# INLINE fBulgeInner #-}
BioInf/RNAwolf/Extern.hs view
@@ -17,7 +17,6 @@ - -- * Unpaired left nucleotide -- | An external loop with an unpaired nucleotide to the left@@ -25,16 +24,17 @@ 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+ | otherwise = VU.singleton ( (i+1,j)+ , nExtn ! (i+1,j)+ ) 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.+--+-- FIXME In btLeftUnpaired, allow only non-empty structures on the right. It+-- would be nice to make the recursion scheme take care of that. btLeftUnpaired :: Params@@ -50,6 +50,7 @@ , testD d' , (x,z) <- btE k l d' , testD z+ , not $ null x -- FIXME ? we only allow left-unpaired "structures" to the left of a real structure ] where ehere = unNExtn nExtn !(i,j) n = VU.length inp -1@@ -103,11 +104,10 @@ 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+ | otherwise = VU.singleton ( (i,j)+ , 0+ ) where- s = ( (i,j)- , 0- ) n = VU.length inp -1 {-# INLINE fOne #-}
BioInf/RNAwolf/Hairpin.hs view
@@ -15,8 +15,6 @@ -- | 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...
BioInf/RNAwolf/Interior.hs view
@@ -13,8 +13,10 @@ import BioInf.Params import BioInf.RNAwolf.Types +import Debug.Trace + -- * Outer part -- | The outer part of an interior loop. Given a certain basepair type, add the@@ -31,6 +33,7 @@ -- + if j-i-1<=maxDistance then pairDistance ! (j-i-1) else 0 ) ijSc = interiorClose ! (((nI,nJ),(ct,eI,eJ)),nIp1,nJm1)+-- ijSc = interiorClose ! (((nN,nN),(cis,wc,wc)),nN,nN) nI = inp VU.! i nJ = inp VU.! j nIp1 = inp VU.! (i+1)@@ -46,20 +49,21 @@ -> NInteLoop -> NBT -- recursive backtracking function for loops -> ExtBT-btInteriorOuter ps inp (EStem eStem) nInteLoop btILoop i j ct eI eJ d =+btInteriorOuter ps inp (EStem eStem) nInteLoop btILoop i j ct eI eJ d = -- iltrc ("ilOuter",i,j,lol) $ [ (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+-- , testD z ] where ij = ((i,j),(ct,eI,eJ)) ehere = eStem!ij n = VU.length inp -1+ lol = VU.toList $ fInteriorOuter ps inp nInteLoop i j ct eI eJ +iltrc k x = trace (show (k,x)) x -- * Loop part@@ -94,14 +98,14 @@ -> NInte -> NBT -> NBT-btInteriorLoop ps inp (NInteLoop nInteLoop) nInte btIL i j d =+btInteriorLoop ps inp (NInteLoop nInteLoop) nInte btIL i j d = -- iltrc ("ilLoop",i,j) $ [ (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+-- , testD z ] where ehere = nInteLoop!(i,j) n = VU.length inp -1@@ -141,14 +145,14 @@ -> EStem -> ExtBT -> NBT-btInteriorInner ps inp (NInte nInte) eStem btES i j d =+btInteriorInner ps inp (NInte nInte) eStem btES i j d = -- iltrc ("ilInner",i,j) $ [ (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+-- , testD z ] where n = VU.length inp -1 ehere = nInte!(i,j)
BioInf/RNAwolf/Stem.hs view
@@ -45,7 +45,6 @@ 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'
RNAwolf.cabal view
@@ -1,5 +1,5 @@ name: RNAwolf-version: 0.3.0.2+version: 0.3.1.0 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/@@ -42,11 +42,15 @@ the runtimes down (as has been done for the extended loops without triplets). .- /We have recently split the Biohaskell libraries into smaller+ We have recently split the Biohaskell libraries into smaller individual libraries. In addition, stacking, intermediate arrays, fusion and newtype-wrapping did require a number of changes. Please send a mail, if you encounter strange behaviour- or bugs./+ or bugs.+ .+ Last Changes:+ .+ * fixed bugs introduced by bulge/interior/multi-loops Flag llvm description: build using llvm backend@@ -64,7 +68,7 @@ vector, PrimitiveArray, BiobaseXNA,- BiobaseTrainingData == 0.1.*,+ BiobaseTrainingData >= 0.1.2.2, StatisticalMethods exposed-modules: BioInf.Keys
RNAwolfTrain.hs view
@@ -112,9 +112,11 @@ putStr "# INFO history:" zipWithM_ (printf " %4d %4.2f") [1::Int ..] $ rhos++[rho] putStrLn ""+ {- print $ sum $ map abs $ P.toList newp print $ minimum $ P.toList newp print $ maximum $ P.toList newp+ -} putStrLn "======================================\n" writeFile (printf "%04d.db" k) . show $ newp return (newp,rhos++[rho])@@ -123,17 +125,22 @@ foldTD :: Options -> Int -> (P.Params,Double,Double,Int) -> (TrainingData,Int) -> IO (P.Params,Double,Double,Int) foldTD o@Options{..} total (!p,accChange,rhosum,cooptimality) (td@TrainingData{},k) = do- print $ length $ primary td let pri = mkPrimary $ primary td let tables = rnaWolf p pri- let bs' = let f x = td{predicted = x} in map (first f) . take (maybe 1 id maxLoss) $ rnaWolfBacktrack p pri 0.00001 tables- let bs = pure $ minimumBy (comparing (sensitivity . mkConfusionMatrix . fst)) bs'+ let bs' = let f x = td{predicted = x} in+ map (first f) + . take (maybe 1 id maxLoss)+ $ rnaWolfBacktrack p pri 0.001 tables+ printf "co-opts: %d\n" $ length bs'+-- mapM_ print bs'+-- print $ rnaWolfOptimal tables+ let bs = pure $ minimumBy (comparing (fmeasure . mkConfusionMatrix . fst)) bs' case bs of [(x,ddd)] -> do let fV = featureVector (primary x) (predicted x) let pVU = VU.fromList . P.toList $ p let sss = map (pVU VU.!) fV- when (abs (ddd - sum sss) > 0.0001) $ do+ when (abs (ddd - sum sss) > 0.001) $ do printf "SCORE DIFFERENCE, backtracking score: %f, sum features: %f\n" ddd (sum sss) -- , " ", map (vks M.!) fV, " ", sss) mapM_ print $ zip (map (vks M.!) fV) sss print "You have found a bug, now write choener to have him fix it!"