diff --git a/BioInf/CMCompare.hs b/BioInf/CMCompare.hs new file mode 100644 --- /dev/null +++ b/BioInf/CMCompare.hs @@ -0,0 +1,354 @@ + +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE GeneralizedNewtypeDeriving #-} +{-# LANGUAGE StandaloneDeriving #-} +{-# LANGUAGE RecordWildCards #-} +{-# LANGUAGE DeriveDataTypeable #-} + +-- | This program compares two Infernal covariance models with each other. +-- Based on the Infernal CM scoring mechanism, a Link sequence and Link score +-- are calculated. The Link sequence is defined as the sequence scoring highest +-- in both models simultanuously. +-- +-- The complete algorithm is described in: +-- +-- "Christian Höner zu Siederdissen, and Ivo L. Hofacker. 2010. Discriminatory +-- power of RNA family models. Bioinformatics 26, no. 18: 453–59." +-- +-- +-- +-- +-- +-- NOTE always use coverage analysis to find out, if we really used all code +-- paths (in long models, if a path is not taken, there is a bug) + +-- NOTE when comparing hits with cmsearch, use the following commandline: +-- +-- cmsearch --no-null3 --cyk --fil-no-hmm --fil-no-qdb +-- +-- --no-null3 : important, the test sequence is so short that null3 can easily +-- generate scores that are way off! remember, we are interested in a sequence +-- that is typically embedded in something large +-- +-- --fil-no-hmm, --fil-no-qdb: do not use heuristics for speedup, they +-- sometimes hide results (in at least one case) +-- +-- (--toponly): if the comparison was done onesided +-- +-- (-g): if you want to compare globally + +module BioInf.CMCompare where + +import Control.Arrow (first,second,(***)) +import Control.Lens +import Control.Monad +import Data.Array.IArray +import Data.List (maximumBy,nub,sort) +import qualified Data.Map as M +import System.Console.CmdArgs +import System.Environment (getArgs) +import Text.Printf + +import Biobase.Primary +import Biobase.SElab.CM +import Biobase.SElab.CM.Import +import Biobase.SElab.Types + + + +-- * optimization functions + +-- | Type of the optimization functions. + +type Opt a = + ( CM -> StateID -> a -- E + , CM -> StateID -> BitScore -> a -> a -- lbegin + , CM -> StateID -> BitScore -> a -> a -- S + , CM -> StateID -> BitScore -> a -> a -- D + , CM -> StateID -> BitScore -> (Char,Char,BitScore) -> a -> a -- MP + , CM -> StateID -> BitScore -> (Char,BitScore) -> a -> a -- ML + , CM -> StateID -> BitScore -> (Char,BitScore) -> a -> a -- IL + , CM -> StateID -> BitScore -> (Char,BitScore) -> a -> a -- MR + , CM -> StateID -> BitScore -> (Char,BitScore) -> a -> a -- IR + , CM -> StateID -> a -> a -> a -- B + , [(a,a)] -> [(a,a)] -- optimization + , a -> String -- finalize, make pretty for output + ) + +-- | Calculates the cyk optimal score over both models. + +cykMaxiMin :: Opt BitScore +cykMaxiMin = (end,lbegin,start,delete,matchP,matchL,insertL,matchR,insertR,branch,opt,finalize) where + end _ _ = 0 + lbegin _ _ t s = t + s + start _ _ t s = t + s + delete _ _ t s = t + s + matchP _ _ t (_,_,e) s = t + e + s + matchL _ _ t (_,e) s = t + e + s + insertL _ _ t (_,e) s = t + e + s + matchR _ _ t (_,e) s = t + e + s + insertR _ _ t (_,e) s = t + e + s + branch _ _ s t = s + t + opt [] = [] + opt xs = [maximumBy (\(a,b) (c,d) -> (min a b) `compare` (min c d)) xs] -- (xs `using` parList rdeepseq)] + finalize s = show s + +-- | Return the nucleotide sequence leading to the score. uses an optional +-- endmarker to denote end states. the string is the same for both models. this +-- is the only Opt function, currently, for which this is true. + +rnaString :: Bool -> Opt [Char] +rnaString endmarker = (end,lbegin,start,delete,matchP,matchL,insertL,matchR,insertR,branch,opt,finalize) where + end _ _ = ['N' | endmarker] + lbegin _ _ _ s = s + start _ _ _ s = s + delete _ _ _ s = s + matchP _ _ _ (k1,k2,_) s = [k1] ++ s ++ [k2] + matchL _ _ _ (k,_) s = k : s + insertL _ _ _ (k,_) s = k : s + matchR _ _ _ (k,_) s = s ++ [k] + insertR _ _ _ (k,_) s = s ++ [k] + branch _ _ s t = s ++ t + opt = id + finalize s = if endmarker + then concatMap f s + else concatMap show s + f x + | x=='N' = "_" + | otherwise = show x + +-- | Dotbracket notation, again with an endmarker, to see the secondary +-- structure corresponding to the rnastring. + +dotBracket :: Bool -> Opt String +dotBracket endmarker = (end,lbegin,start,delete,matchP,matchL,insertL,matchR,insertR,branch,opt,finalize) where + end _ _ = ['_' | endmarker] + lbegin _ _ _ s = s + start _ _ _ s = s + delete _ _ _ s = s + matchP _ _ _ _ s = "(" ++ s ++ ")" + matchL _ _ _ _ s = '.' : s + insertL _ _ _ _ s = ',' : s + matchR _ _ _ _ s = s ++ "." + insertR _ _ _ _ s = s ++ "," + branch _ _ s t = s ++ t + opt = id + finalize s = s + +-- | Show the nodes which were visited to get the score. the last node can +-- occur multiple times. if it does, local end transitions were used. + +visitedNodes :: Opt [NodeID] +visitedNodes = (end,lbegin,start,delete,matchP,matchL,insertL,matchR,insertR,branch,opt,finalize) where + end cm k = [((cm^.states) M.! k) ^. nodeID] + lbegin cm k _ s = s + start cm k _ s = ((cm^.states) M.! k) ^. nodeID : s + delete cm k _ s = ((cm^.states) M.! k) ^. nodeID : s + matchP cm k _ _ s = ((cm^.states) M.! k) ^. nodeID : s + matchL cm k _ _ s = ((cm^.states) M.! k) ^. nodeID : s + insertL cm k _ _ s = ((cm^.states) M.! k) ^. nodeID : s + matchR cm k _ _ s = ((cm^.states) M.! k) ^. nodeID : s + insertR cm k _ _ s = ((cm^.states) M.! k) ^. nodeID : s + branch cm k s t = ((cm^.states) M.! k) ^. nodeID : (s ++ t) + opt = id -- NOTE do not sort, do not nub ! + finalize xs = (show $ map unNodeID xs) -- NOTE do not sort, do not nub ! + +-- | Detailed output of the different states, that were visited. + +extendedOutput :: Opt String +extendedOutput = (end,lbegin,start,delete,matchP,matchL,insertL,matchR,insertR,branch,opt,finalize) where + end cm sid = printf "E %5d %5d" (unStateID sid) (unNodeID $ ((cm^.states) M.! sid)^.nodeID) + lbegin cm sid t s = printf "lbegin %5d %5d %7.3f \n%s" (unStateID sid) (unNodeID $ ((cm^.states) M.! sid)^.nodeID) (unBitScore t) s + start cm sid t s = printf "S %5d %5d %7.3f \n%s" (unStateID sid) (unNodeID $ ((cm^.states) M.! sid)^.nodeID) (unBitScore t) s + delete cm sid t s = printf "D %5d %5d %7.3f \n%s" (unStateID sid) (unNodeID $ ((cm^.states) M.! sid)^.nodeID) (unBitScore t) s + matchP cm sid t (k1,k2,e) s = printf "MP %5d %5d %7.3f %7.3f %1s %1s\n%s" (unStateID sid) (unNodeID $ ((cm^.states) M.! sid)^.nodeID) (unBitScore t) (unBitScore e) (show k1) (show k2) s + matchL cm sid t (k,e) s = printf "ML %5d %5d %7.3f %7.3f %1s\n%s" (unStateID sid) (unNodeID $ ((cm^.states) M.! sid)^.nodeID) (unBitScore t) (unBitScore e) (show k) s + insertL cm sid t (k,e) s = printf "IL %5d %5d %7.3f %7.3f %1s\n%s" (unStateID sid) (unNodeID $ ((cm^.states) M.! sid)^.nodeID) (unBitScore t) (unBitScore e) (show k) s + matchR cm sid t (k,e) s = printf "MR %5d %5d %7.3f %7.3f %1s\n%s" (unStateID sid) (unNodeID $ ((cm^.states) M.! sid)^.nodeID) (unBitScore t) (unBitScore e) (show k) s + insertR cm sid t (k,e) s = printf "IR %5d %5d %7.3f %7.3f %1s\n%s" (unStateID sid) (unNodeID $ ((cm^.states) M.! sid)^.nodeID) (unBitScore t) (unBitScore e) (show k) s + branch cm sid s t = printf "B %5d %5d\n%s\n%s" (unStateID sid) (unNodeID $ ((cm^.states) M.! sid) ^. nodeID) s t + opt = id + finalize s = "\nLabel State Node Trans Emis\n\n" ++ s + +-- | Algebra product operation. + +(<*>) :: Eq a => Opt a -> Opt b -> Opt (a,b) +algA <*> algB = (end,lbegin,start,delete,matchP,matchL,insertL,matchR,insertR,branch,opt,finalize) where + (endA,lbeginA,startA,deleteA,matchPA,matchLA,insertLA,matchRA,insertRA,branchA,optA,finalizeA) = algA + (endB,lbeginB,startB,deleteB,matchPB,matchLB,insertLB,matchRB,insertRB,branchB,optB,finalizeB) = algB + end cm k = (endA cm k, endB cm k) + lbegin cm k t (sA,sB) = (lbeginA cm k t sA, lbeginB cm k t sB) + start cm k t (sA,sB) = (startA cm k t sA, startB cm k t sB) + delete cm k t (sA,sB) = (deleteA cm k t sA, deleteB cm k t sB) + matchP cm k t e (sA,sB) = (matchPA cm k t e sA, matchPB cm k t e sB) + matchL cm k t e (sA,sB) = (matchLA cm k t e sA, matchLB cm k t e sB) + insertL cm k t e (sA,sB) = (insertLA cm k t e sA, insertLB cm k t e sB) + matchR cm k t e (sA,sB) = (matchRA cm k t e sA, matchRB cm k t e sB) + insertR cm k t e (sA,sB) = (insertRA cm k t e sA, insertRB cm k t e sB) + branch cm k (sA,sB) (tA,tB) = (branchA cm k sA tA, branchB cm k sB tB) + opt xs = [((xl1,xl2),(xr1,xr2)) | (xl1,xr1) <- nub $ optA [(yl1,yr1) | ((yl1,yl2),(yr1,yr2)) <- xs] + , (xl2,xr2) <- optB [(yl2,yr2) | ((yl1,yl2),(yr1,yr2)) <- xs, (yl1,yr1) == (xl1,xr1)] + ] + finalize (sA,sB) = finalizeA sA ++ "\n" ++ finalizeB sB + + + +-- * The grammar for CM comparison. + +-- | Recursion in two CMs simultanously. + +recurse :: Bool -> Opt a -> CM -> CM -> Array (StateID,StateID) [(a,a)] +recurse fastIns (end,lbegin,start,delete,matchP,matchL,insertL,matchR,insertR,branch,opt,finalize) m1 m2 = locarr where + + loc k1 k2 + | otherwise = opt $ do + r <- arr ! (k1, k2) + return $ (lbegin m1 k1 lb1 *** lbegin m2 k2 lb2) r + where + lb1 = M.findWithDefault (BitScore (-10000)) k1 (m1^.localBegin) + lb2 = M.findWithDefault (BitScore (-10000)) k2 (m2^.localBegin) + + rec k1 k2 = let xyz = rec' k1 k2 + in xyz -- traceShow ("rec",k1,((m1^.states) M.! k1) ^. stateType,k2,((m2^.states) M.! k2) ^. stateType) xyz + rec' k1 k2 + -- + | t1 == E && t2 == E = [(end m1 k1, end m2 k2)] + -- + | t1 == S && t2 == S = opt $ do + (c1,tr1) <- s1 ^. transitions ++ [(ls1,le1)] + (c2,tr2) <- s2 ^. transitions ++ [(ls2,le2)] + r <- arr ! (c1, c2) + return $ (start m1 k1 tr1 *** start m2 k2 tr2) r + | t1 == D && t2 == D = opt $ do + (c1,tr1) <- s1 ^. transitions ++ [(ls1,le1)] + (c2,tr2) <- s2 ^. transitions ++ [(ls2,le2)] + r <- arr ! (c1, c2) + return $ (delete m1 k1 tr1 *** delete m2 k2 tr2) r + -- match pair emitting states + | t1 == MP && t2 == MP + = opt $ do + (c1,tr1) <- s1 ^. transitions ++ [(ls1,le1)] + (c2,tr2) <- s2 ^. transitions ++ [(ls2,le2)] + (e1,e2) <- zip (s1 ^. emits ^. pair) (s2 ^. emits ^. pair) + r <- arr ! (c1, c2) + return $ (matchP m1 k1 tr1 e1 *** matchP m2 k2 tr2 e2) r + -- match left emitting states + | t1 `elem` lstates && t2 `elem` lstates + = opt $ do + (c1,tr1) <- s1 ^. transitions ++ [(ls1,le1)] + (c2,tr2) <- s2 ^. transitions ++ [(ls2,le2)] + guard $ (not fastIns && (c1 /= k1 || c2 /= k2)) || (fastIns && c1/=k1 && c2/=k2) + (e1,e2) <- zip (s1 ^. emits ^. single) (s2 ^. emits ^. single) + r <- arr ! (c1, c2) + let f = if t1 == ML then matchL else insertL + let g = if t2 == ML then matchL else insertL + return $ (f m1 k1 tr1 e1 *** g m2 k2 tr2 e2) r + -- match right emitting states + | t1 `elem` rstates && t2 `elem` rstates + = opt $ do + (c1,tr1) <- s1 ^. transitions ++ [(ls1,le1)] + (c2,tr2) <- s2 ^. transitions ++ [(ls2,le2)] + guard $ (not fastIns && (c1 /= k1 || c2 /= k2)) || (fastIns && c1/=k1 && c2/=k2) + (e1,e2) <- zip (s1 ^. emits ^. single) (s2 ^. emits ^. single) + r <- arr ! (c1, c2) + let f = if t1 == MR then matchR else insertR + let g = if t2 == MR then matchR else insertR + return $ (f m1 k1 tr1 e1 *** g m2 k2 tr2 e2) r + -- if one state is E, we can only delete states, except for another S state, which will go into local end + -- it is not possible to use an emitting state on the right as those would require emitting on the left, too! + | t1 == E && t2 `elem` [D,S] = opt $ do + (c2,tr2) <- s2 ^. transitions ++ [(ls2,le2)] + r <- arr ! (k1,c2) + return $ if t2 == D then second (delete m2 k2 tr2) r else second (start m2 k2 tr2) r + -- the other way around with D,E + | t1 `elem` [D,S] && t2 == E = opt $ do + (c1,tr1) <- s1 ^. transitions ++ [(ls1,le2)] + r <- arr ! (c1,k2) + return $ if t1 == D then first (delete m1 k1 tr1) r else first (start m1 k1 tr1) r + -- two branching states + | t1 == B && t2 == B = opt $ + let + [(l1,_),(r1,_)] = s1 ^. transitions + [(l2,_),(r2,_)] = s2 ^. transitions + in + -- both branches are matched + do + (s1,s2) <- arr ! (l1,l2) -- left branch (m1,m2) + (t1,t2) <- arr ! (r1,r2) -- right branch (m1,m2) + return (branch m1 k1 s1 t1, branch m2 k2 s2 t2) -- (m1,m2) + ++ + do + (t1,s2) <- arr ! (r1,l2) -- match right branch of m1 with left branch of m2 + -- local ends for other branches + x <- arr ! (ls1,ls2) + let (s1,t2) = (delete m1 l1 le1 *** delete m2 l2 le2) x + return (branch m1 k1 s1 t1, branch m2 k2 s2 t2) + ++ + do + (s1,t2) <- arr ! (l1,r2) + x <- arr ! (ls1,ls2) + let (t1,s2) = (delete m1 l1 le1 *** delete m2 l2 le2) x + return (branch m1 k1 s1 t1, branch m2 k2 s2 t2) + -- branch - non-branch + | t1 == B && t2 /= B = opt $ + let + [(l,_), (r,_)] = s1 ^. transitions + in + do + (s1,s2) <- arr ! (l,k2) -- left branch and m2 + x <- arr ! (ls1,ls2) + -- dont do anything for ls2, since we do not have to + -- delete a branch in model 2. + let (t1,t2) = first (delete m1 r le1) x + return (branch m1 k1 s1 t1, branch m2 k2 s2 t2) + ++ + do + (t1,t2) <- arr ! (r,k2) -- right branch and m2 + x <- arr ! (ls1,ls2) + let (s1,s2) = first (delete m1 l le1) x -- delete left branch in m1 + return (branch m1 k1 s1 t1, branch m2 k2 s2 t2) + -- branch - non-branch + | t1 /= B && t2 == B = opt $ + let + [(l,_), (r,_)] = s2 ^. transitions + in + do + (s1,s2) <- arr ! (k1,l) + x <- arr ! (ls1,ls2) + let (t1,t2) = second (delete m2 r le2) x + return (branch m1 k1 s1 t1, branch m2 k2 s2 t2) + ++ + do + (t1,t2) <- arr ! (k1,r) + x <- arr ! (ls1,ls2) + let (s1,s2) = second (delete m2 l le2) x + return (branch m1 k1 s1 t1, branch m2 k2 s2 t2) + -- S state versus any + | t1 == S = opt $ do + (c1,tr1) <- s1 ^. transitions ++ [(ls1,le1)] + r <- arr ! (c1, k2) + return $ first (start m1 k1 tr1) r + -- S state versus any + | t2 == S = opt $ do + (c2,tr2) <- s2 ^. transitions ++ [(ls2,le2)] + r <- arr ! (k1, c2) + return $ second (start m2 k2 tr2) r + -- + | otherwise = [] + where + s1 = (m1 ^. states) M.! k1 + s2 = (m2 ^. states) M.! k2 + t1 = s1 ^. stateType + t2 = s2 ^. stateType + le1 = M.findWithDefault (BitScore (-10000)) k1 (m1^.localEnd) + le2 = M.findWithDefault (BitScore (-10000)) k2 (m2^.localEnd) + ls1 = sn1 + ls2 = sn2 + lstates = [ML,IL] + rstates = [MR,IR] + + locarr = (array ((0,0),(sn1,sn2)) [((k1,k2),loc k1 k2) | k1 <- [0 .. sn1], k2 <- [0 .. sn2]]) + arr = (array ((0,0),(sn1,sn2)) [((k1,k2),rec k1 k2) | k1 <- [0 .. sn1], k2 <- [0 .. sn2]]) `asTypeOf` locarr + sn1 = fst . M.findMax $ m1 ^. states + sn2 = fst . M.findMax $ m2 ^. states + diff --git a/CMCompare.cabal b/CMCompare.cabal --- a/CMCompare.cabal +++ b/CMCompare.cabal @@ -1,5 +1,5 @@ name: CMCompare -version: 0.0.1.2 +version: 0.0.1.4 author: Christian Hoener zu Siederdissen, Ivo L. Hofacker maintainer: choener@tbi.univie.ac.at copyright: Christian Hoener zu Siederdissen, Ivo L. Hofacker, 2010 @@ -10,29 +10,48 @@ license-file: LICENSE build-type: Simple stability: experimental -cabal-version: >= 1.4.0 +cabal-version: >= 1.6.0 description: Compares two Infernal covariance models. Returns the common MaxiMin score and the offending RNA sequence. High scores point toward low discriminative power of the two models. Based on: - "Discriminatory Power or RNA Family Models, Hoener zu - Siederdissen and Hofacker, 2010, accepted for eccb10" + . + "Christian Höner zu Siederdissen, and Ivo L. Hofacker. 2010. + Discriminatory power of RNA family models. Bioinformatics 26, + no. 18: 453–59" + . + extra-source-files: scripts/HighScoreEdges.sh scripts/NeighborGraph.sh -executable hsCMCompare +library + build-depends: - base >= 4 && < 5, - cmdargs == 0.6.4, - array, + base >= 4 && < 5 , + array == 0.4.* , + containers == 0.5.* , + lens == 3.* , - Biobase >= 0.1.0 && < 0.2.0, - HsTools >= 0.0.1.1 && < 0.0.2 + BiobaseXNA == 0.6.3.* , + BiobaseInfernal == 0.7.* + exposed-modules: + BioInf.CMCompare + + ghc-options: + -O2 -rtsopts -fllvm -optlo-O3 -optlo-inline -optlo-std-compile-opts + + + +executable CMCompare + + build-depends: + cmdargs >= 0.10 + main-is: CMCompare.hs ghc-options: - -O2 + -O2 -rtsopts -fllvm -optlo-O3 -optlo-inline -optlo-std-compile-opts diff --git a/CMCompare.hs b/CMCompare.hs --- a/CMCompare.hs +++ b/CMCompare.hs @@ -1,16 +1,23 @@ +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE GeneralizedNewtypeDeriving #-} {-# LANGUAGE StandaloneDeriving #-} {-# LANGUAGE RecordWildCards #-} {-# LANGUAGE DeriveDataTypeable #-} --- Based on: Discriminatory Power of RNA Family Models, Christian Hoener zu --- Siederdissen and Ivo Hofacker, 2010, accepted for eccb10: +-- | This program compares two Infernal covariance models with each other. +-- Based on the Infernal CM scoring mechanism, a Link sequence and Link score +-- are calculated. The Link sequence is defined as the sequence scoring highest +-- in both models simultanuously. -- --- Preprint: +-- The complete algorithm is described in: -- --- http://www.tbi.univie.ac.at/newpapers/abstracts/abstractTBI-p-2010-5.html - - - +-- "Christian Höner zu Siederdissen, and Ivo L. Hofacker. 2010. Discriminatory +-- power of RNA family models. Bioinformatics 26, no. 18: 453–59." +-- +-- +-- +-- +-- -- NOTE always use coverage analysis to find out, if we really used all code -- paths (in long models, if a path is not taken, there is a bug) @@ -29,335 +36,26 @@ -- -- (-g): if you want to compare globally --- {{{ module descriptor + module Main where -import Data.Array.IArray -import Text.Printf +import Control.Arrow (first,second,(***)) import Control.Monad -import Debug.Trace -import System.Environment (getArgs) +import Data.Array.IArray import System.Console.CmdArgs -import Data.List (maximumBy,nub,sort) -import Control.Arrow (first,second,(***)) +import Text.Printf -import Biobase.Infernal.CM -import Biobase.Infernal.CM.Import -import Biobase.RNA hiding (nucE) -import qualified Biobase.RNA as RNA -import Debug.Trace.Tools +import Biobase.SElab.CM +import Biobase.SElab.CM.Import +import Biobase.SElab.Types --- }}} +import BioInf.CMCompare --- * optimization functions --- {{{ type of the optimization functions -type StateID = Int -- TODO should this go into BiobaseCM? -type CM' = CM () () - -type Opt a = - ( CM' -> StateID -> a -- E - , CM' -> StateID -> Double -> a -> a -- lbegin - , CM' -> StateID -> Double -> a -> a -- S - , CM' -> StateID -> Double -> a -> a -- D - , CM' -> StateID -> Double -> Emission -> a -> a -- MP - , CM' -> StateID -> Double -> Emission -> a -> a -- ML - , CM' -> StateID -> Double -> Emission -> a -> a -- IL - , CM' -> StateID -> Double -> Emission -> a -> a -- MR - , CM' -> StateID -> Double -> Emission -> a -> a -- IR - , CM' -> StateID -> a -> a -> a -- B - , [(a,a)] -> [(a,a)] -- optimization - , a -> String -- finalize, make pretty for output - ) - --- }}} - --- {{{ optimization functions - --- | calculates the cyk optimal score over both models. - -cykMaxiMin :: Opt Double -cykMaxiMin = (end,lbegin,start,delete,matchP,matchL,insertL,matchR,insertR,branch,opt,finalize) where - end _ _ = 0 - lbegin _ _ t s = t + s - start _ _ t s = t + s - delete _ _ t s = t + s - matchP _ _ t (EmitP _ _ e) s = t + e + s - matchL _ _ t (EmitS _ e) s = t + e + s - insertL _ _ t (EmitS _ e) s = t + e + s - matchR _ _ t (EmitS _ e) s = t + e + s - insertR _ _ t (EmitS _ e) s = t + e + s - branch _ _ s t = s + t - opt [] = [] - opt xs = [maximumBy (\(a,b) (c,d) -> (min a b) `compare` (min c d)) xs] - finalize s = show s - --- | return the nucleotide sequence leading to the score. uses an optional --- endmarker to denote end states. the string is the same for both models. this --- is the only Opt function, currently, for which this is true. - -rnaString :: Bool -> Opt [Nucleotide] -rnaString endmarker = (end,lbegin,start,delete,matchP,matchL,insertL,matchR,insertR,branch,opt,finalize) where - end _ _ = [RNA.nucE | endmarker] - lbegin _ _ _ s = s - start _ _ _ s = s - delete _ _ _ s = s - matchP _ _ _ (EmitP k1 k2 _) s = [k1] ++ s ++ [k2] - matchL _ _ _ (EmitS k _) s = k : s - insertL _ _ _ (EmitS k _) s = k : s - matchR _ _ _ (EmitS k _) s = s ++ [k] - insertR _ _ _ (EmitS k _) s = s ++ [k] - branch _ _ s t = s ++ t - opt = id - finalize s = if endmarker - then concatMap f s - else concatMap show s - f x - | x==RNA.nucE = "_" - | otherwise = show x - --- | dotbracket notation, again with an endmarker, to see the secondary --- structure corresponding to the rnastring. - -dotBracket :: Bool -> Opt String -dotBracket endmarker = (end,lbegin,start,delete,matchP,matchL,insertL,matchR,insertR,branch,opt,finalize) where - end _ _ = ['_' | endmarker] - lbegin _ _ _ s = s - start _ _ _ s = s - delete _ _ _ s = s - matchP _ _ _ _ s = "(" ++ s ++ ")" - matchL _ _ _ _ s = '.' : s - insertL _ _ _ _ s = ',' : s - matchR _ _ _ _ s = s ++ "." - insertR _ _ _ _ s = s ++ "," - branch _ _ s t = s ++ t - opt = id - finalize s = s - --- | show the nodes which were visited to get the score. the last node can --- occur multiple times. if it does, local end transitions were used. - -visitedNodes :: Opt [Int] -visitedNodes = (end,lbegin,start,delete,matchP,matchL,insertL,matchR,insertR,branch,opt,finalize) where - end cm k = [snode (states cm ! k)] - lbegin cm k _ s = s - start cm k _ s = snode (states cm ! k) : s - delete cm k _ s = snode (states cm ! k) : s - matchP cm k _ _ s = snode (states cm ! k) : s - matchL cm k _ _ s = snode (states cm ! k) : s - insertL cm k _ _ s = snode (states cm ! k) : s - matchR cm k _ _ s = snode (states cm ! k) : s - insertR cm k _ _ s = snode (states cm ! k) : s - branch cm k s t = snode (states cm ! k) : (s ++ t) - opt = id -- NOTE do not sort, do not nub ! - finalize xs = "Nodes: " ++ show xs -- NOTE do not sort, do not nub ! - --- | detailed output of the different states, that were visited. - -extendedOutput :: Opt String -extendedOutput = (end,lbegin,start,delete,matchP,matchL,insertL,matchR,insertR,branch,opt,finalize) where - end cm sid = printf "E %5d %5d" sid (snode (states cm ! sid)) - lbegin cm sid t s = printf "lbegin %5d %5d %7.3f \n%s" sid (snode (states cm ! sid)) t s - start cm sid t s = printf "S %5d %5d %7.3f \n%s" sid (snode (states cm ! sid)) t s - delete cm sid t s = printf "D %5d %5d %7.3f \n%s" sid (snode (states cm ! sid)) t s - matchP cm sid t (EmitP k1 k2 e) s = printf "MP %5d %5d %7.3f %7.3f %1s %1s\n%s" sid (snode (states cm ! sid)) t e (show k1) (show k2) s - matchL cm sid t (EmitS k e) s = printf "ML %5d %5d %7.3f %7.3f %1s\n%s" sid (snode (states cm ! sid)) t e (show k) s - insertL cm sid t (EmitS k e) s = printf "IL %5d %5d %7.3f %7.3f %1s\n%s" sid (snode (states cm ! sid)) t e (show k) s - matchR cm sid t (EmitS k e) s = printf "MR %5d %5d %7.3f %7.3f %1s\n%s" sid (snode (states cm ! sid)) t e (show k) s - insertR cm sid t (EmitS k e) s = printf "IR %5d %5d %7.3f %7.3f %1s\n%s" sid (snode (states cm ! sid)) t e (show k) s - branch cm sid s t = printf "B %5d %5d\n%s\n%s" sid (snode (states cm ! sid)) s t - opt = id - finalize s = "\nLabel State Node Trans Emis\n\n" ++ s - -(<*>) :: Eq a => Opt a -> Opt b -> Opt (a,b) -algA <*> algB = (end,lbegin,start,delete,matchP,matchL,insertL,matchR,insertR,branch,opt,finalize) where - (endA,lbeginA,startA,deleteA,matchPA,matchLA,insertLA,matchRA,insertRA,branchA,optA,finalizeA) = algA - (endB,lbeginB,startB,deleteB,matchPB,matchLB,insertLB,matchRB,insertRB,branchB,optB,finalizeB) = algB - end cm k = (endA cm k, endB cm k) - lbegin cm k t (sA,sB) = (lbeginA cm k t sA, lbeginB cm k t sB) - start cm k t (sA,sB) = (startA cm k t sA, startB cm k t sB) - delete cm k t (sA,sB) = (deleteA cm k t sA, deleteB cm k t sB) - matchP cm k t e (sA,sB) = (matchPA cm k t e sA, matchPB cm k t e sB) - matchL cm k t e (sA,sB) = (matchLA cm k t e sA, matchLB cm k t e sB) - insertL cm k t e (sA,sB) = (insertLA cm k t e sA, insertLB cm k t e sB) - matchR cm k t e (sA,sB) = (matchRA cm k t e sA, matchRB cm k t e sB) - insertR cm k t e (sA,sB) = (insertRA cm k t e sA, insertRB cm k t e sB) - branch cm k (sA,sB) (tA,tB) = (branchA cm k sA tA, branchB cm k sB tB) - opt xs = [((xl1,xl2),(xr1,xr2)) | (xl1,xr1) <- nub $ optA [(yl1,yr1) | ((yl1,yl2),(yr1,yr2)) <- xs] - , (xl2,xr2) <- optB [(yl2,yr2) | ((yl1,yl2),(yr1,yr2)) <- xs, (yl1,yr1) == (xl1,xr1)] - ] - finalize (sA,sB) = finalizeA sA ++ "\n" ++ finalizeB sB - --- }}} - --- * recursion in two CMs simultanously - --- {{{ main recursion - -recurse :: Opt a -> CM' -> CM' -> Array (Int,Int) [(a,a)] -recurse (end,lbegin,start,delete,matchP,matchL,insertL,matchR,insertR,branch,opt,finalize) m1 m2 = locarr where - - loc k1 k2 - | cmType m1 == CMProb || cmType m2 == CMProb = error "both models need to be score type models" - | otherwise = opt $ do - r <- arr ! (k1, k2) - return $ (lbegin m1 k1 lb1 *** lbegin m2 k2 lb2) r - where - lb1 = localBegin m1 ! k1 - lb2 = localBegin m2 ! k2 - - rec k1 k2 - -- - | t1 == E && t2 == E = [(end m1 k1, end m2 k2)] - -- - | t1 == S && t2 == S = opt $ do - Transition c1 tr1 <- schildren s1 ++ [Transition ls1 le1 | acceptLE le1] - Transition c2 tr2 <- schildren s2 ++ [Transition ls2 le2 | acceptLE le2] - r <- arr ! (c1, c2) - return $ (start m1 k1 tr1 *** start m2 k2 tr2) r - | t1 == D && t2 == D = opt $ do - Transition c1 tr1 <- schildren s1 ++ [Transition ls1 le1 | acceptLE le1] - Transition c2 tr2 <- schildren s2 ++ [Transition ls2 le2 | acceptLE le2] - r <- arr ! (c1, c2) - return $ (delete m1 k1 tr1 *** delete m2 k2 tr2) r - -- match pair emitting states - | t1 == MP && t2 == MP - = opt $ do - Transition c1 tr1 <- schildren s1 ++ [Transition ls1 le1 | acceptLE le1] - Transition c2 tr2 <- schildren s2 ++ [Transition ls2 le2 | acceptLE le2] - (e1,e2) <- zip (semission s1) (semission s2) - r <- arr ! (c1, c2) - return $ (matchP m1 k1 tr1 e1 *** matchP m2 k2 tr2 e2) r - -- match left emitting states - | t1 `elem` lstates && t2 `elem` lstates - = opt $ do - Transition c1 tr1 <- schildren s1 ++ [Transition ls1 le1 | acceptLE le1] - Transition c2 tr2 <- schildren s2 ++ [Transition ls2 le2 | acceptLE le2] - guard $ c1 /= k1 || c2 /= k2 - (e1,e2) <- zip (semission s1) (semission s2) - r <- arr ! (c1, c2) - let f = if t1 == ML then matchL else insertL - let g = if t2 == ML then matchL else insertL - return $ (f m1 k1 tr1 e1 *** g m2 k2 tr2 e2) r - -- match right emitting states - | t1 `elem` rstates && t2 `elem` rstates - = opt $ do - Transition c1 tr1 <- schildren s1 ++ [Transition ls1 le1 | acceptLE le1] - Transition c2 tr2 <- schildren s2 ++ [Transition ls2 le2 | acceptLE le2] - guard $ c1 /= k1 || c2 /= k2 - (e1,e2) <- zip (semission s1) (semission s2) - r <- arr ! (c1, c2) - let f = if t1 == MR then matchR else insertR - let g = if t2 == MR then matchR else insertR - return $ (f m1 k1 tr1 e1 *** g m2 k2 tr2 e2) r - -- if one state is E, we can only delete states, except for another S state, which will go into local end - -- it is not possible to use an emitting state on the right as those would require emitting on the left, too! - | t1 == E && t2 `elem` [D,S] = opt $ do - Transition c2 tr2 <- schildren s2 ++ [Transition ls2 le2 | acceptLE le2] - r <- arr ! (k1,c2) - return $ if t2 == D then second (delete m2 k2 tr2) r else second (start m2 k2 tr2) r - -- the other way around with D,E - | t1 `elem` [D,S] && t2 == E = opt $ do - Transition c1 tr1 <- schildren s1 ++ [Transition ls1 le1 | acceptLE le1] - r <- arr ! (c1,k2) - return $ if t1 == D then first (delete m1 k1 tr1) r else first (start m1 k1 tr1) r - -- two branching states - | t1 == B && t2 == B = opt $ - let - [Branch l1, Branch r1] = schildren s1 - [Branch l2, Branch r2] = schildren s2 - in - -- both branches are matched - do - (s1,s2) <- arr ! (l1,l2) -- left branch (m1,m2) - (t1,t2) <- arr ! (r1,r2) -- right branch (m1,m2) - return (branch m1 k1 s1 t1, branch m2 k2 s2 t2) -- (m1,m2) - ++ - do - (t1,s2) <- arr ! (r1,l2) -- match right branch of m1 with left branch of m2 - -- local ends for other branches - x <- arr ! (ls1,ls2) - let (s1,t2) = (delete m1 l1 le1 *** delete m2 l2 le2) x - return (branch m1 k1 s1 t1, branch m2 k2 s2 t2) - ++ - do - (s1,t2) <- arr ! (l1,r2) - x <- arr ! (ls1,ls2) - let (t1,s2) = (delete m1 l1 le1 *** delete m2 l2 le2) x - return (branch m1 k1 s1 t1, branch m2 k2 s2 t2) - -- branch - non-branch - | t1 == B && t2 /= B = opt $ - let - [Branch l, Branch r] = schildren s1 - in - do - (s1,s2) <- arr ! (l,k2) -- left branch and m2 - x <- arr ! (ls1,ls2) - -- dont do anything for ls2, since we do not have to - -- delete a branch in model 2. - let (t1,t2) = first (delete m1 r le1) x - return (branch m1 k1 s1 t1, branch m2 k2 s2 t2) - ++ - do - (t1,t2) <- arr ! (r,k2) -- right branch and m2 - x <- arr ! (ls1,ls2) - let (s1,s2) = first (delete m1 l le1) x -- delete left branch in m1 - return (branch m1 k1 s1 t1, branch m2 k2 s2 t2) - -- branch - non-branch - | t1 /= B && t2 == B = opt $ - let - [Branch l, Branch r] = schildren s2 - in - do - (s1,s2) <- arr ! (k1,l) - x <- arr ! (ls1,ls2) - let (t1,t2) = second (delete m2 r le2) x - return (branch m1 k1 s1 t1, branch m2 k2 s2 t2) - ++ - do - (t1,t2) <- arr ! (k1,r) - x <- arr ! (ls1,ls2) - let (s1,s2) = second (delete m2 l le2) x - return (branch m1 k1 s1 t1, branch m2 k2 s2 t2) - -- S state versus any - | t1 == S = opt $ do - Transition c1 tr1 <- schildren s1 ++ [Transition ls1 le1 | acceptLE le1] - r <- arr ! (c1, k2) - return $ first (start m1 k1 tr1) r - -- S state versus any - | t2 == S = opt $ do - Transition c2 tr2 <- schildren s2 ++ [Transition ls2 le2 | acceptLE le2] - r <- arr ! (k1, c2) - return $ second (start m2 k2 tr2) r - -- - | otherwise = [] - where - s1 = states m1 ! k1 - s2 = states m2 ! k2 - t1 = stype s1 - t2 = stype s2 - le1 = localEnd m1 ! k1 - le2 = localEnd m2 ! k2 - ls1 = snd . bounds $ states m1 -- last state (E) - ls2 = snd . bounds $ states m2 - lstates = [ML,IL] - rstates = [MR,IR] - acceptLE x - | cmType m1 == CMScore && x > (-1/0) = True - | cmType m1 == CMProb && x /= 0 = True - | otherwise = False - - locarr = (array ((0,0),(sn1,sn2)) [((k1,k2),loc k1 k2) | k1 <- [0 .. sn1], k2 <- [0 .. sn2]]) `asTypeOf` arr - arr = (array ((0,0),(sn1,sn2)) [((k1,k2),rec k1 k2) | k1 <- [0 .. sn1], k2 <- [0 .. sn2]]) `asTypeOf` locarr - (_,sn1) = bounds $ states m1 - (_,sn2) = bounds $ states m2 - --- }}} - --- {{{ main - +-- * Handling user I/O. +-- -- TODO add an option to filter by minimal score (default: -1000000) if the -- filter is on, we get a result only, if the score is above the threshold @@ -368,6 +66,7 @@ , pend :: Double , endmarker :: Bool , nobeginilir :: Bool + , fastIns :: Bool , models :: [String] } deriving (Show,Data,Typeable) @@ -378,6 +77,7 @@ , pend = 0.05 &= help "aggregate local end probability" , endmarker = False &= help "add an endmarker into the rnastring to denote local ends" , nobeginilir = False &= help "trailing left or right nucleotides change the score" + , fastIns = False &= help "fast insertion heuristic" , models = def &= args -- &= help "path to exactly two covariance models" } &= summary "CMCompare: Discriminatory Power of RNA Family Models" &= help "(c) 2010, Christian Hoener zu Siederdissen and Ivo Hofacker\nchoener@tbi.univie.ac.at\nLicensed under the GPLv3\n" &= verbosity -- TODO put fixTransition in BiobaseInfernal, rename and whatnot @@ -388,6 +88,7 @@ -- second can only, if IL eats a nucleotide. for single CM search, this is -- taken care of by the DP algorithm which doesn't work here. +{- fixTransition :: CM' -> CM' fixTransition x = x{states = ss // [(0,rtnew)]} where ss = states x @@ -395,13 +96,14 @@ tr = schildren rt trnew = [Transition 1 0, Transition 2 0] ++ drop 2 tr rtnew = rt{schildren = trnew} +-} applyIf c f = if c then f else id -answers optf m1 m2 = map (finalize *** finalize) . opt . concat . elems $ recurse optf m1 m2 where +answers fastIns optf m1 m2 = map (finalize *** finalize) . opt . concat . elems $ recurse fastIns optf m1 m2 where (end,lbegin,start,delete,matchP,matchL,insertL,matchR,insertR,branch,opt,finalize) = optf -results optf m1 m2 = opt . concat . elems $ recurse optf m1 m2 where +results fastIns optf m1 m2 = opt . concat . elems $ recurse fastIns optf m1 m2 where (end,lbegin,start,delete,matchP,matchL,insertL,matchR,insertR,branch,opt,finalize) = optf main = do @@ -412,36 +114,20 @@ unless (length models == 2) $ do fail "give exactly two CMs" let [a,b] = models - Right [m1'] <- fromFile a - Right [m2'] <- fromFile b - let m1 = applyIf (not nobeginilir) fixTransition $ applyIf (not global) (cmMakeLocal pbegin pend) m1' - let m2 = applyIf (not nobeginilir) fixTransition $ applyIf (not global) (cmMakeLocal pbegin pend) m2' + [theA] <- fromFile a + [theB] <- fromFile b + let m1 = if nobeginilir then theA else makeLocal pbegin pend theA -- TODO !!! applyIf (not nobeginilir) fixTransition $ applyIf (not global) (cmMakeLocal pbegin pend) m1' + let m2 = if nobeginilir then theA else makeLocal pbegin pend theB -- TODO !!! applyIf (not nobeginilir) fixTransition $ applyIf (not global) (cmMakeLocal pbegin pend) m2' let pr = (\(x,y) -> putStrLn x >> putStrLn "" >> putStrLn y >> putStrLn "++++++++++++") when (quiet && not normal) $ do - let (a1,a2) = head $ results cykMaxiMin m1 m2 - printf "%s %s %10.3f %10.3f\n" a b a1 a2 + let (a1,a2) = head $ results fastIns cykMaxiMin m1 m2 + printf "%s %s %10.3f %10.3f\n" a b (unBitScore a1) (unBitScore a2) when (normal && not loud) $ do - let ((((cyk1,vn1),rna1),db1),(((cyk2,vn2),_),db2)) = head $ results (cykMaxiMin <*> visitedNodes <*> rnaString endmarker <*> dotBracket endmarker) m1 m2 - let rnaBoth = concatMap f rna1 where + let ((((cyk1,vn1),rna1),db1),(((cyk2,vn2),_),db2)) = head $ results fastIns (cykMaxiMin <*> visitedNodes <*> rnaString endmarker <*> dotBracket endmarker) m1 m2 + let rnaBoth = map f rna1 where f x - | x==RNA.nucE = "_" - | otherwise = show x - printf "%s %s %10.3f %10.3f %s %s %s %s %s\n" a b cyk1 cyk2 rnaBoth db1 db2 (show vn1) (show vn2) - when loud . mapM_ pr $ answers (cykMaxiMin <*> visitedNodes <*> rnaString endmarker <*> dotBracket endmarker <*> extendedOutput) m1 m2 - --- }}} - --- * Helper functions - --- {{{ - --- | summation in logspace - --- TODO time for a new library ;) - -logSum x y = h + log (1 + exp (l-h)) where - h = max x y - l = min x y - --- }}} + | x=='N' = '_' + | otherwise = x + printf "%s %s %10.3f %10.3f %s %s %s %s %s\n" a b (unBitScore cyk1) (unBitScore cyk2) rnaBoth db1 db2 (show $ map unNodeID vn1) (show $ map unNodeID vn2) + when loud . mapM_ pr $ answers fastIns (cykMaxiMin <*> visitedNodes <*> rnaString endmarker <*> dotBracket endmarker <*> extendedOutput) m1 m2