CMCompare 0.0.1.2 → 0.0.1.4
raw patch · 3 files changed
+426/−367 lines, 3 filesdep +BiobaseInfernaldep +BiobaseXNAdep +containersdep −Biobasedep −HsToolsdep ~arraydep ~cmdargsnew-component:exe:CMCompare
Dependencies added: BiobaseInfernal, BiobaseXNA, containers, lens
Dependencies removed: Biobase, HsTools
Dependency ranges changed: array, cmdargs
Files
- BioInf/CMCompare.hs +354/−0
- CMCompare.cabal +30/−11
- CMCompare.hs +42/−356
+ BioInf/CMCompare.hs view
@@ -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."+--+-- <http://bioinformatics.oxfordjournals.org/content/26/18/i453.long>+--+--+--+-- 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+
CMCompare.cabal view
@@ -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"+ .+ <http://bioinformatics.oxfordjournals.org/content/26/18/i453.long> 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
CMCompare.hs view
@@ -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."+--+-- <http://bioinformatics.oxfordjournals.org/content/26/18/i453.long>+--+--+-- -- 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