diff --git a/BioInf/RNAdesign.hs b/BioInf/RNAdesign.hs --- a/BioInf/RNAdesign.hs +++ b/BioInf/RNAdesign.hs @@ -1,91 +1,49 @@ -{-# LANGUAGE RankNTypes #-} {-# LANGUAGE NoMonomorphismRestriction #-} {-# LANGUAGE RecordWildCards #-} module BioInf.RNAdesign where +import Control.Arrow (first,second) +import Control.Monad.Primitive +import Control.Monad.Primitive.Class +import Data.List (nub,group,sort,(\\),genericLength) +import Data.Tuple.Select (sel1) import qualified Data.Array.IArray as A -import System.IO.Unsafe -import Control.Monad.IO.Class -import Control.Monad.Primitive.Class -import Control.Monad.Primitive -import System.Random.MWC.Monad -import Control.Monad -import qualified Data.Vector.Unboxed as VU -import qualified Data.Vector as V -import Data.List (sort,group) import qualified Data.Map as M -import qualified Data.Vector.Fusion.Stream.Monadic as SM -import Control.Arrow -import System.IO.Unsafe -- TODO remove -import Data.List -import Data.Tuple.Select - -import Biobase.Primary -import Biobase.Secondary.Diagrams -import Biobase.Secondary -import Biobase.Vienna -import qualified BioInf.ViennaRNA.Bindings as RNA -- NOTE removes the ability to call into ghci! -import BioInf.ViennaRNA.Eval +import qualified Data.Vector as V +import qualified Data.Vector.Unboxed as VU +import System.IO.Unsafe (unsafePerformIO) +import System.Random.MWC.Monad -import BioInf.RNAdesign.Graph -import BioInf.RNAdesign.OptParser -import BioInf.RNAdesign.Assignment -import BioInf.RNAdesign.LogMultinomial +import Biobase.Primary +import Biobase.Primary.IUPAC +import Biobase.Secondary.Diagrams +import Biobase.Secondary (PairIdx(..)) +import Biobase.Vienna +import qualified BioInf.ViennaRNA.Bindings as RNA -import Debug.Trace +import BioInf.RNAdesign.Assignment +import BioInf.RNAdesign.CandidateChain +import BioInf.RNAdesign.Graph +import BioInf.RNAdesign.LogMultinomial +import BioInf.RNAdesign.OptParser --- A single candidate, with its sequence and the score, this sequence receives. --- Candidates are ordered by their scores. - -data Candidate = Candidate - { candidate :: Primary - , score :: Score - } deriving (Eq,Show) - -instance Ord Candidate where - (Candidate _ a) <= (Candidate _ b) = ropt a <= ropt b - --- | Create an initial, legal, candidate. Give it a really bad score. - -mkInitial :: (MonadPrim m, PrimMonad m) => Int -> DesignProblem -> Rand m Candidate -mkInitial l dp = do - let z = VU.replicate l nA - c <- foldM mutateOneAssignment z $ assignments dp - return $ Candidate c (Score [] 999999) - -{- --- | Sum probabilities over base pairs in the structural constraints - -sumProbStructures :: Primary -> [D1Secondary] -> Double -sumProbStructures inp ss = s where - s = sum $ map ((bp A.!) . first (+1) . second (+1)) ps - ps = concatMap snd (map fromD1S ss :: [(Int,[PairIdx])]) - bp = let (_,_,bp') = unsafePerformIO (RNA.part $ concatMap show $ VU.toList inp) in bp' - -sumProbNotStructures :: Primary -> [D1Secondary] -> Double -sumProbNotStructures inp ss = undefined - -probabilityDefect inp str = s where - s = sum (map (bp A.!) ps) - sum (map (bp A.!) ups) - ups = [ (i,j) | i<-[1..l], j<-[i..l] ] \\ ps - (l,ps) = second (map (first (+1) . second (+1))) $ fromD1S str :: (Int,[PairIdx]) - bp = let (_,_,bp') = unsafePerformIO (RNA.part $ concatMap show $ VU.toList inp) in bp' --} +-- | probabilityDefectAll inp ss = s where ca :: A.Array (Int,Int) Double ca = A.amap (\c -> c / n) . A.accumArray (+) 0 ((1,1),(l,l)) $ zip ps (repeat 1) n = genericLength ss --- s = sum (map (abs . (n-) . (bp A.!)) ps) + sum (map (bp A.!) ups) s = sum (map (\ix -> abs $ ca A.! ix - bp A.! ix) ps) + sum (map (bp A.!) ups) l = VU.length inp ups = [ (i,j) | i<-[1..l], j<-[i..l] ] \\ ps ps = map (first (+1) . second (+1)) $ concatMap snd (map fromD1S ss :: [(Int,[PairIdx])]) bp = let (_,_,bp') = unsafePerformIO (RNA.part $ concatMap show $ VU.toList inp) in bp' +-- | + ensembleDefect inp str = s where s = n - 2 * sps - sus n = fromIntegral $ VU.length inp @@ -105,9 +63,6 @@ sops = [ ("eos" , \k -> unsafePerformIO $ RNA.eos (concatMap show (VU.toList inp)) (fromD1S $ secs !! (k-1))) , ("ed" , \k -> ensembleDefect inp (secs !! (k-1))) -- ensemble defect --- , ("pdef" , \k -> probabilityDefect inp (secs !! (k-1))) --- [ ("EOS",\k -> let (Deka e) = fst $ rnaEval ener inp $ secs !! (k-1) in fromIntegral e / 100) --- , ("PF" ... ] mops = [ ("sum",sum) @@ -118,44 +73,29 @@ [ ("Ged" , probabilityDefectAll inp secs) -- global ensemble defect a la ``me'' , ("gibbs" , sel1 . unsafePerformIO $ RNA.part (concatMap show (VU.toList inp))) , ("mfe" , fst . unsafePerformIO $ RNA.mfe (concatMap show (VU.toList inp))) --- , ("Pin" , sumProbStructures inp secs) ] props = [ ("logMN", \ps -> lmn ps inp) ] +-- | + lmn ps inp = logMultinomial l p c where l = VU.length inp p = VU.fromList ps cM = M.fromList . map (\z -> (head z, length z)) . group . sort $ VU.toList inp c = VU.fromList $ map (\z -> M.findWithDefault 0 z cM) acgu -data Score = Score - { eoss :: [Deka] - , ropt :: Double - } deriving (Eq,Show,Read) - -instance Ord Score where - (Score _ a) <= (Score _ b) = a<=b +-- | scoreSequence :: String -> Vienna2004 -> DesignProblem -> Primary -> Score scoreSequence optfun ener DesignProblem{..} s = score where - score = Score - { eoss = error "don't call this" -- map (fst . rnaEval ener s) structures - , ropt = resolveOpt optfun ener s structures - } - --- | This structure defines a "design problem" - -data DesignProblem = DesignProblem - { structures :: [D1Secondary] - , assignments :: [Assignment] - } deriving (Eq,Read,Show) + score = Score $ resolveOpt optfun ener s structures -- | Given a set of structures, create the set of independent graphs and -- assignment possibilities. -mkDesignProblem :: Int -> [String] -> [String] -> DesignProblem +mkDesignProblem :: Int -> [String] -> String -> DesignProblem mkDesignProblem asnLimit xs scs = dp where dp = DesignProblem { structures = map mkD1S xs @@ -163,93 +103,8 @@ } gs = independentGraphs xs as = map (allCandidates asnLimit sv) gs - ss = M.map fixup . M.unionsWith ((nub .) . (++)) $ map (M.fromList . zip [0..] . (map ((:[]). mkNuc))) scs + --ss = M.map fixup . M.unionsWith ((nub .) . (++)) $ map (M.fromList . zip [0..] . (map ((:[]). mkNuc))) scs + ss = M.map fixup . M.fromList . zip [0..] . map (map mkNuc . fromSymbol) $ scs sv = V.fromList $ map (\k -> M.findWithDefault acgu k ss) [0 .. length (head xs) - 1] fixup zs = filter (/=nN) $ if (all (==nN) zs) then acgu else zs - -unfoldStreamNew - :: forall m . (MonadPrim m, PrimMonad m) - => Int -> Int -> Int -> (Primary -> Score) -> (Candidate -> Candidate -> Rand m Bool) -> DesignProblem -> Candidate -> SM.Stream (Rand m) Candidate -unfoldStreamNew burnin number thin score f dp = go where - go s = SM.map snd -- remove remaining indices from stream - . SM.take number -- take the number of sequences we want - . SM.filter ((==0) . flip mod thin . fst) -- keep only every thin'th sequence - . SM.indexed -- add index - . SM.drop burnin -- drop the burnin sequences - . SM.drop 1 -- drop original input - . SM.scanlM' (mutateOneAssignmentCandidateWith score f) s -- starting with 's', mutate s further and further using cycled assignments - $ SM.unfoldr (Just . first head . splitAt 1) (cycle $ assignments dp) -- create inifinite cycled assignments - -unfoldStream - :: forall m . (MonadPrim m, PrimMonad m) - => Int -> Int -> Int -> (Primary -> Primary -> Rand m Bool) -> DesignProblem -> Primary -> SM.Stream (Rand m) Primary -unfoldStream burnin number thin f dp = go where - go s = SM.map snd -- remove remaining indices from stream - . SM.take number -- take the number of sequences we want - . SM.filter ((==0) . flip mod thin . fst) -- keep only every thin'th sequence - . SM.indexed -- add index - . SM.drop burnin -- drop the burnin sequences - . SM.drop 1 -- drop original input - . SM.scanlM' (mutateOneAssignmentWith f) s -- starting with 's', mutate s further and further using cycled assignments - $ SM.unfoldr (Just . first head . splitAt 1) (cycle $ assignments dp) -- create inifinite cycled assignments - --- | Mutate the sequence in a candidate - -mutateOneAssignmentCandidateWith - :: (MonadPrim m, PrimMonad m) - => (Primary -> Score) -> (Candidate -> Candidate -> Rand m Bool) -> Candidate -> Assignment -> Rand m Candidate -mutateOneAssignmentCandidateWith score f old Assignment{..} = do - i <- uniformR (0,V.length assignment -1) -- inclusive range for Int - let cs = VU.zip columns (assignment V.! i) - let nw = VU.update (candidate old) cs - let new = Candidate nw (score nw) - b <- f old new - return $ if b then new else old - --- | Mutate the sequence using one assignment with evaluation function. - -mutateOneAssignmentWith - :: (MonadPrim m, PrimMonad m) - => (Primary -> Primary -> Rand m Bool) -> Primary -> Assignment -> Rand m Primary -mutateOneAssignmentWith f old Assignment{..} = do - i <- uniformR (0,V.length assignment -1) -- inclusive range for Int - let cs = VU.zip columns (assignment V.! i) - let new = VU.update old cs - b <- f old new - return $ if b then new else old - --- | Create a number of sequences, thinning the list of candidates to yield --- more independent candidates. The optimization function is used to make the --- choice between emitting the current candidate again and selecting a new one. - -generateSequences - :: (MonadPrim m, PrimMonad m) - => Int -> Int -> (Primary -> Primary -> Rand m Bool) -> DesignProblem -> Primary -> Rand m [Primary] -generateSequences number thin f dp s = go number thin s where - go n t s - | n < 1 = return [] - | t == 0 = do s' <- mutateSequence f dp s - ss <- go (n-1) thin s' - return $ s' : ss - | otherwise = mutateSequence f dp s >>= go n (t-1) - --- | Mutate a sequence using the possible assignments. - -mutateSequence - :: (MonadPrim m, PrimMonad m) - => (Primary -> Primary -> Rand m Bool) -> DesignProblem -> Primary -> Rand m Primary -mutateSequence f dp old = do - new <- foldM mutateOneAssignment old $ assignments dp - b <- f old new - return $ if b then new else old - --- | Mutate the sequence using one assignment. - -mutateOneAssignment - :: (MonadPrim m, PrimMonad m) - => Primary -> Assignment -> Rand m Primary -mutateOneAssignment s Assignment{..} = do - i <- uniformR (0,V.length assignment -1) -- inclusive range for Int - let cs = VU.zip columns (assignment V.! i) - return $ VU.update s cs diff --git a/BioInf/RNAdesign/Assignment.hs b/BioInf/RNAdesign/Assignment.hs --- a/BioInf/RNAdesign/Assignment.hs +++ b/BioInf/RNAdesign/Assignment.hs @@ -3,15 +3,16 @@ module BioInf.RNAdesign.Assignment where -import Control.Arrow -import Data.Graph.Inductive.Graph -import Data.List (nub,sortBy,sort,genericLength) -import Data.Ord +import Control.Arrow +import Control.Lens +import Control.Lens.Tuple +import Data.Function +import Data.Graph.Inductive.Graph +import Data.Graph.Inductive.Query +import Data.List (nub,sortBy,sort,genericLength) +import Data.Ord import qualified Data.Vector as V import qualified Data.Vector.Unboxed as VU -import Control.Lens.Tuple -import Control.Lens -import Data.Function import Biobase.Primary import Biobase.Secondary.Vienna @@ -19,10 +20,9 @@ import BioInf.RNAdesign.Graph -import Data.Graph.Inductive.Query -import Debug.Trace +-- | data Assignment = Assignment { columns :: VU.Vector Int diff --git a/BioInf/RNAdesign/CandidateChain.hs b/BioInf/RNAdesign/CandidateChain.hs new file mode 100644 --- /dev/null +++ b/BioInf/RNAdesign/CandidateChain.hs @@ -0,0 +1,96 @@ +{-# LANGUAGE NoMonomorphismRestriction #-} +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE RecordWildCards #-} + +module BioInf.RNAdesign.CandidateChain where + +import Control.Arrow (first) +import Control.Monad (foldM) +import Control.Monad.Primitive +import Control.Monad.Primitive.Class +import Data.Function (on) +import qualified Data.Vector as V +import qualified Data.Vector.Fusion.Stream.Monadic as SM +import qualified Data.Vector.Unboxed as VU +import System.Random.MWC.Monad + +import Biobase.Primary +import Biobase.Secondary.Diagrams +import Biobase.Vienna + +import BioInf.RNAdesign.Assignment (Assignment(..)) + + + +-- | A single candidate, with its sequence and the score, this sequence +-- receives. Candidates are ordered by their scores. + +data Candidate = Candidate + { candidate :: Primary + , score :: Score + } deriving (Eq,Show) + +instance Ord Candidate where + (<=) = (<=) `on` score + +-- | The likelihood score we get. +-- +-- TODO replace Score Likelihood / LogLikelihood (once we switch to the more +-- generic MCMC library) + +newtype Score = Score { unScore :: Double } + deriving (Eq,Ord,Show,Read) + +-- | This structure defines a "design problem" + +data DesignProblem = DesignProblem + { structures :: [D1Secondary] + , assignments :: [Assignment] + } deriving (Eq,Read,Show) + +-- | Create an initial, legal, candidate. + +mkInitial :: (MonadPrim m, PrimMonad m) => (Primary -> Score) -> Int -> DesignProblem -> Rand m Candidate +mkInitial scoring l dp = do + let z = VU.replicate l nA + foldM (mutateOneAssignmentWith scoring (\_ _ -> return True)) (Candidate z (scoring z)) $ assignments dp + +-- | Create a stream of 'Candidate's from an initial candidate. + +unfoldStream + :: forall m . (MonadPrim m, PrimMonad m) + => Int -> Int -> Int -> (Primary -> Score) -> (Candidate -> Candidate -> Rand m Bool) -> DesignProblem -> Candidate + -> SM.Stream (Rand m) Candidate +unfoldStream burnin number thin score f dp = go where + go s = SM.map snd -- remove remaining indices from stream + . SM.take number -- take the number of sequences we want + . SM.filter ((==0) . flip mod thin . fst) -- keep only every thin'th sequence + . SM.indexed -- add index + . SM.drop burnin -- drop the burnin sequences + . SM.drop 1 -- drop original input + . SM.scanlM' (mutateOneAssignmentWith score f) s -- starting with 's', mutate s further and further using cycled assignments + $ SM.unfoldr (Just . first head . splitAt 1) + (cycle $ assignments dp) -- create inifinite cycled assignments + +-- | Mutate the current (or "old") sequence under the possible 'assignment's as +-- prescribed by 'Assignment'. The modifying assignment is selected uniformly. +-- The monadic @old -> new -> Rand m Bool@ function chooses between the old and +-- the new candidate. It can be used to, e.g., allow always choosing "new" if +-- it is better, but choosing "new" as well if some stochastic value (hence +-- dependence on @Rand m@) indicates so. + +mutateOneAssignmentWith + :: (MonadPrim m, PrimMonad m) + => (Primary -> Score) -- ^ the likelihood function, gives a sequence a score + -> (Candidate -> Candidate -> Rand m Bool) -- ^ choose between old and new sequence (monadic, stochastic) + -> Candidate -- ^ "old" / current sequence + -> Assignment -- ^ possible assignments for the sequence + -> Rand m Candidate -- ^ the "new" sequence +mutateOneAssignmentWith score f old Assignment{..} = do + i <- uniformR (0,V.length assignment -1) -- inclusive range for Int + let cs = VU.zip columns (assignment V.! i) + let nw = VU.update (candidate old) cs + let new = Candidate nw (score nw) + b <- f old new + return $ if b then new else old + diff --git a/BioInf/RNAdesign/Graph.hs b/BioInf/RNAdesign/Graph.hs --- a/BioInf/RNAdesign/Graph.hs +++ b/BioInf/RNAdesign/Graph.hs @@ -1,12 +1,12 @@ module BioInf.RNAdesign.Graph where -import Data.Graph.Inductive.Graph +import Control.Arrow (first,second) import Data.Graph.Inductive.Basic +import Data.Graph.Inductive.Graph import Data.Graph.Inductive.PatriciaTree import Data.Graph.Inductive.Query import Data.List (nub,partition) -import Control.Arrow (first,second) import Biobase.Secondary.Diagrams diff --git a/BioInf/RNAdesign/LogMultinomial.hs b/BioInf/RNAdesign/LogMultinomial.hs --- a/BioInf/RNAdesign/LogMultinomial.hs +++ b/BioInf/RNAdesign/LogMultinomial.hs @@ -4,6 +4,7 @@ import qualified Data.Vector.Unboxed as VU + logMultinomial :: Int -> VU.Vector Double -> VU.Vector Int -> Double logMultinomial n' ps xs' | VU.length ps /= VU.length xs' = error "logMultinomial: P-vector and count-vector of unequal length" diff --git a/BioInf/RNAdesign/OptParser.hs b/BioInf/RNAdesign/OptParser.hs --- a/BioInf/RNAdesign/OptParser.hs +++ b/BioInf/RNAdesign/OptParser.hs @@ -11,12 +11,12 @@ ( parseOptString ) where +import Control.Applicative import Text.Parsec.Expr import Text.Parsec hiding ((<|>)) import Text.Parsec.Language -import Text.Parsec.Token -import Control.Applicative import Text.Parsec.String +import Text.Parsec.Token import Text.Parsec.Numbers @@ -58,7 +58,7 @@ g x = f (read x) mkMultiOp :: NumSecStructs -> (SingleOp,MultiOp) -> GenParser Char st Double -mkMultiOp nss ((s,sf),(m,mf)) = {- (\xs -> error $ show (xs, map sf xs, mf $ map sf xs)) <$ -} (\xs -> mf $ map sf xs) <$ +mkMultiOp nss ((s,sf),(m,mf)) = (\xs -> mf $ map sf xs) <$ string m <* string "(" <* string s <* string "," <*> secs <* string ")" where secs = try ([1..nss] <$ string "all") <|> map read <$> many1 digit `sepBy1` string "," @@ -73,7 +73,7 @@ parseGlobalOp gops = choice $ map (try . mkGlobalOp) gops optable = [ [prefix "-" negate, prefix "+" id] - , [binary "^" (**) AssocLeft] --, binary "**" (**) AssocLeft] + , [binary "^" (**) AssocLeft] , [binary "*" (*) AssocLeft, binary "/" (/) AssocLeft] , [binary "+" (+) AssocLeft, binary "-" (-) AssocLeft] ] diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -1,1 +1,90 @@ +RNAdesign +========= + +The RNAdesign program solves the multi-target RNA sequence design problem. You +can give one or more structural targets for which a single compatible sequence +is designed. + +PAPER +===== + +Christian Hoener zu Siederdissen, Stefan Hammer, Ingrid Abfalter, Ivo L. Hofacker, Christoph Flamm, Peter F. Stadler. +Computational Design of RNAs with Complex Energy Landscapes. +2013. Biopolymers. 99, no. 12. 99. 1124–36. http://dx.doi.org/10.1002/bip.22337. + +Contact +======= + +choener@tbi.univie.ac.at + + + +HOW TO USE RNAdesign +==================== + +RNAdesign designs RNA sequences given one or more structural targets. The +program offers a variety of optimization functions that each can be used to +optimize candidate sequence towards a certain goal, say, minimal ensemble +defect or small energetic distance to another target structure. + + +RNAdesign input +--------------- + +Structural targets are given via stdin, preferably via an input file. Below is +a the small tri-stable from our paper, which you should then pipe to RNAdesign: +"echo tri-stable.dat | RNAdesign" + +"cat tri-stable.dat:" + + # a tri-stable example target. (optional comment) + ((((....))))....((((....))))........ + ........((((....((((....))))....)))) + ((((((((....))))((((....))))....)))) + # below follows a simple (and optional) sequence constraint. + CKNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNB + +The input may contain many comments lines, starting with a hash "#" and at most +one sequence constraint line. All of these lines are optional, except of course +for the structural constraints. + + +Optimization functions +---------------------- + +Depending on the actual design you are looking for, you'll want to modify the +optimization function. Below, the different options available are detailed. By +giving a complex "--optfun", many different design goals can be tried. + +A good optimization goal is (as an example for three targets): + +--optfun "eos(1)+eos(2)+eos(3) - 3 * gibbs + 1 * ((eos(1)-eos(2))^2 + (eos(1)-eos(3))^2 + (eos(2)-eos(3))^2)" + +This way, the sequence produces close-to-mfe foldings with the targets (left) +and the targets are close together in terms of energy. (1 * ) scales the two +terms according to user choice. + +### binary, combining: + ++ - * / :: the four basic operations +^ :: (^) generalized power function + +### binary, apply function to many targets: + +sum max min :: run function over set of targets: sum(eos,1,2) or sum(eos,all) + +### unary, apply to single target: + +eos :: energy of a structure: eos(1) +ed :: ensemble defect of a structure: ed(3) + +### nullary, constant for the current sequence: + +Ged :: global, weighted ensemble defect: Ged +gibbs :: gibbs free energy of sequence +mfe :: minimum free energy of sequence + +### special: + +logMN :: requires four parameters logMN(0.2,0.3,0.3,0.2) penalizes according to given mono-nucleotide distribution in order of ACGU diff --git a/RNAdesign.cabal b/RNAdesign.cabal --- a/RNAdesign.cabal +++ b/RNAdesign.cabal @@ -1,6 +1,7 @@ name: RNAdesign -version: 0.1.0.0 -author: Christian Hoener zu Siederdissen +version: 0.1.1.0 +author: Christian Hoener zu Siederdissen, 2013-2014 +copyright: Christian Hoener zu Siederdissen, Stefan Hammer, Ingrid Abfalter, Ivo L. Hofacker, Christoph Flamm, Peter F. Stadler, 2013-2014 maintainer: choener@tbi.univie.ac.at category: Bioinformatics synopsis: Multi-target RNA sequence design @@ -27,11 +28,11 @@ . If you find this program useful, please cite: . + @ Christian Hoener zu Siederdissen, Stefan Hammer, Ingrid Abfalter, Ivo L. Hofacker, Christoph Flamm, Peter F. Stadler - . - @Computational design of RNAs with complex energy landscapes@ - . - Biopolymers, 99, 12, 1124-1136, 2013, Wiley + Computational design of RNAs with complex energy landscapes + 2013. Biopolymers. 99, no. 12. 99. 1124–36. http://dx.doi.org/10.1002/bip.22337 + @ @@ -58,14 +59,15 @@ fgl-extras-decompositions >= 0.1.0.0 , BiobaseTurner >= 0.3.1.1 , BiobaseVienna >= 0.3 , - BiobaseXNA >= 0.7.0.1 , + BiobaseXNA >= 0.8.1 , ParsecTools >= 0.0.2 && < 0.0.3 , - PrimitiveArray >= 0.5 , + PrimitiveArray >= 0.5.3 , RNAFold >= 1.99.3.3 , ViennaRNA-bindings >= 0.1.0.0 exposed-modules: BioInf.RNAdesign BioInf.RNAdesign.Assignment + BioInf.RNAdesign.CandidateChain BioInf.RNAdesign.Graph BioInf.RNAdesign.LogMultinomial BioInf.RNAdesign.OptParser @@ -74,29 +76,9 @@ executable RNAdesign build-depends: - base >= 4 && < 5 , - array >= 0.4 , - cmdargs == 0.10.* , - containers , - fgl >= 5.4 , - lens >= 3.9 , - monad-primitive >= 0.1 , - mwc-random-monad >= 0.6 , - parallel >= 3.2 , - parsec >= 3 , - primitive >= 0.5 , - random >= 1.0 , - transformers >= 0.3 , - tuple >= 0.2 , - vector >= 0.10 , - fgl-extras-decompositions >= 0.1.0.0 , - BiobaseTurner >= 0.3.1.1 , - BiobaseVienna >= 0.3 , - BiobaseXNA >= 0.7.0.1 , - ParsecTools >= 0.0.2 && < 0.0.3 , - PrimitiveArray >= 0.5 , - RNAFold >= 1.99.3.3 , - ViennaRNA-bindings >= 0.1.0.0 + bytestring >= 0.10 , + cmdargs == 0.10.* , + file-embed >= 0.0.6 main-is: RNAdesign.hs ghc-options: diff --git a/RNAdesign.hs b/RNAdesign.hs --- a/RNAdesign.hs +++ b/RNAdesign.hs @@ -1,3 +1,5 @@ +{-# LANGUAGE DoAndIfThenElse #-} +{-# LANGUAGE TemplateHaskell #-} {-# LANGUAGE QuasiQuotes #-} {-# LANGUAGE BangPatterns #-} {-# LANGUAGE ScopedTypeVariables #-} @@ -5,139 +7,115 @@ {-# LANGUAGE NoMonomorphismRestriction #-} {-# LANGUAGE RecordWildCards #-} --- | +-- | Given one or more structural constraints, and possibly sequence +-- constraints for certain columns, design a sequence which is optimal +-- according to a user-defined optimization function. Optimization works via a +-- Markov Chain. module Main where -import System.Console.CmdArgs -import Data.List -import Data.Char (isAlpha) -import Control.Monad -import System.Random.MWC.Monad -import System.Random.MWC.Distributions.Monad -import qualified Data.Vector.Unboxed as VU -import Text.Printf -import Data.Ord +import Control.Arrow +import Control.Monad +import Data.Char (isAlpha) +import Data.FileEmbed +import Data.Function +import Data.List +import Data.Ord +import Data.Version (showVersion) +import qualified Data.ByteString.Char8 as BS import qualified Data.Vector.Fusion.Stream.Monadic as SM -import Control.Arrow -import Data.Function -import System.IO +import qualified Data.Vector.Unboxed as VU +import System.Console.CmdArgs +import System.IO +import System.Random.MWC.Distributions.Monad +import System.Random.MWC.Monad +import Text.Printf -import Biobase.Primary -import Biobase.Vienna +import Biobase.Primary +import Biobase.Vienna import qualified Biobase.Turner.Import as TI import BioInf.RNAdesign +import BioInf.RNAdesign.CandidateChain import BioInf.RNAdesign.Assignment -import Debug.Trace +import Paths_RNAdesign (version) --- * configuration +-- * Configuration data Config = Config - { number :: Int - , thin :: Int - , burnin :: Int - , scale :: Double - , optfun :: String - , veclen :: Int - , turner :: String --- , exhaustive :: Bool -- TODO want to think about this for number of structures > 3, IF the total sequence space size is less than say 100.000 - , initial :: String + { number :: Int + , thin :: Int + , burnin :: Int + , scale :: Double + , optfun :: String + , veclen :: Int + , turner :: String + , initial :: String , sequenceConstraints :: Bool - , explore :: Bool + , explore :: Bool + , showManual :: Bool } deriving (Show,Data,Typeable) config = Config - { number = 50 &= help "Number of candidate sequences to generate (50)" - , thin = 50 &= help "keep only every n'th sequence (50)" - , burnin = 100 &= help "remove the first burnin sequences (100)" - , scale = 1 &= help "acceptance scale parameter (1); exponentially distributed with mean 'scale^(-1)' (smaller scale means longer jumps)" - , optfun = "sum(eos,all)" &= help "optimization function, \"sum(eos,all)\" tries to minimize the sum of the energies" - , veclen = 1000000 &= help "multiple structure constraints lead to large connected components, veclen restricts the number of component solutions to store." - , turner = "./params" &= help "directory containing the Turner 2004 RNA energy tables (with a default of \"./params/\"" --- , exhaustive = False &= help "exhaustively search the nucleotide space" - , initial = "" &= help "start from this initial sequence" - , explore = def &= help "explore sequences, do not sort of nub list" - , sequenceConstraints = def &= help "activate sequence constraints" - } &= help shortHelp &= details longHelp &= summary "RNAdesign v0.0.2, (C) Christian Hoener zu Siederdissen 2013, choener@tb.univie.ac.at" &= program "RNAdesign" + { number = 50 &= help "Number of candidate sequences to generate (50)" + , thin = 50 &= help "keep only every n'th sequence (50)" + , burnin = 100 &= help "remove the first burnin sequences (100)" + , scale = 1 &= help "acceptance scale parameter (1); exponentially distributed with mean 'scale^(-1)' (smaller scale means longer jumps)" + , optfun = "sum(eos,all)" &= help "optimization function, \"sum(eos,all)\" tries to minimize the sum of the energies" + , veclen = 1000000 &= help "multiple structure constraints lead to large connected components, veclen restricts the number of component solutions to store." + , turner = "./params" &= help "directory containing the Turner 2004 RNA energy tables (with a default of \"./params/\"" + , initial = "" &= help "start from this initial sequence" + , explore = def &= help "explore sequences, do not sort of nub list" + , sequenceConstraints = def &= help "activate sequence constraints" + , showManual = def &= help "" + } &= help shortHelp + &= details [] + &= summary ("RNAdesign " ++ showVersion version ++ " (C) Christian Hoener zu Siederdissen 2013--2014, choener@tbi.univie.ac.at") + &= program "RNAdesign" shortHelp = "The defaults work acceptably well and produce a results extremely fast. " -longHelp = - [ "RNAdesign designs RNA sequences given one or more structural targets. The" - , "program offers a variety of optimization functions that each can be used to" - , "optimize candidate sequence towards a certain goal, say, minimal ensemble" - , "defect or small energetic distance to another target structure. By giving a" - , "complex \"--optun\", many different design goals can be tried. The following" - , "functions are available:" - , "binary, combining:" - , "+ - * / :: the four basic operations" - , "^ :: (^) generalized power function" - , "" - , "binary, apply function to many targets:" - , "sum max min :: run function over set of targets: sum(eos,1,2) or sum(eos,all)" - , "" - , "unary, apply to single target:" - , "eos :: energy of a structure: eos(1)" - , "ed :: ensemble defect of a structure: ed(3)" - , "nullary, constant for the current sequence:" - , "Ged :: global, weighted ensemble defect: Ged" - , "gibbs :: gibbs free energy of sequence" - , "mfe :: minimum free energy of sequence" - , "" - , "special:" - , "logMN :: requires four parameters logMN(0.2,0.3,0.3,0.2) penalizes according to given mono-nucleotide distribution" - , "" - , "A good optimization goal is (as an example for three targets):" - , "--optfun \"eos(1)+eos(2)+eos(3) - 3*gibbs +" - , " 1 * ((eos(1)-eos(2))^2 + (eos(1)-eos(3))^2 + (eos(2)-eos(3))^2)\"" - , "This way, the sequence produces close-to-mfe foldings with the targets (left) and the targets are close together in terms of energy. (1 *) scales the two terms according to user choice." - , "\n\n\n" - , "If you find this tool useful, please cite:" - , "" - , "Christian Hoener zu Siederdissen, Stefan Hammer, Ingrid Abfalter, Ivo L. Hofacker, Christoph Flamm, Peter F. Stadler." - , "A Graph Coloring Approach to Designing Multi-Stable Nucleic Acid Sequences." - , "submitted, 2013." - , "" - , "Contact: choener@tbi.univie.ac.at" - , "Given one or more structures in dot-bracket format of the same length, returns a compatible assignment of nucleotides." - , "Compatible nucleotides are those that allow folding of the sequence into all given structures." - ] +embeddedManual = $(embedFile "README.md") main = do hSetBuffering stdout NoBuffering hSetBuffering stderr NoBuffering cmds@Config{..} <- cmdArgs config + if showManual + then BS.putStrLn embeddedManual + else do turner <- fmap turnerToVienna $ TI.fromDir turner "" ".dat" - strs' <- fmap lines $ getContents + strs' <- fmap (filter ((/="#") . take 1) . lines) $ getContents let (scs,strs) = partition (any isAlpha) . filter ((">"/=) . take 1) $ strs' unless (length strs > 0) $ error "no structures given!" let l = length $ head strs unless (all ((l==) . length) strs) $ error "structures of different size detected" - let dp = mkDesignProblem veclen strs (if sequenceConstraints then scs else []) + unless (not sequenceConstraints || sequenceConstraints && length scs<=1) $ error "sequence constraint error" + let dp = mkDesignProblem veclen strs (if sequenceConstraints && length scs==1 then head scs else "") let defOpt old new = let oldS = scoreSequence optfun turner dp old newS = scoreSequence optfun turner dp new in do t <- exponential scale - return $ ropt newS <= ropt oldS || t >= ropt newS - ropt oldS + return $ unScore newS <= unScore oldS || t >= unScore newS - unScore oldS let calcScore = scoreSequence optfun turner dp let walk old new = do t <- exponential scale - let sn = ropt $ score new - let so = ropt $ score old + let sn = unScore $ score new + let so = unScore $ score old return $ sn <= so || t >= sn - so let ini = if null initial -- start from initial sequence or generate one from the ensemble - then mkInitial l dp - else return $ Candidate (mkPrimary initial) (Score [] 999999) - xs <- runWithSystemRandom . asRandIO $ (ini >>= SM.toList . unfoldStreamNew burnin number thin calcScore walk dp) + then mkInitial calcScore l dp + else let pri = mkPrimary initial + in return $ Candidate pri (calcScore pri) + xs <- runWithSystemRandom . asRandIO $ (ini >>= SM.toList . unfoldStream burnin number thin calcScore walk dp) let pna = product . map numAssignments $ assignments dp - printf "# Size of sequence space: %d %s\n\n" pna {-(product . map numAssignments $ assignments dp)-} (show . map numAssignments $ assignments dp) + printf "# Size of sequence space: %d %s\n\n" pna (show . map numAssignments $ assignments dp) unless (pna>0) $ error "empty sequence space, aborting!" - mapM_ (\ys -> printf "%s %4d %8.2f\n" (concatMap show . VU.toList . candidate . head $ ys) (length ys) (ropt . score $ head ys)) + mapM_ (\ys -> printf "%s %4d %8.2f\n" (concatMap show . VU.toList . candidate . head $ ys) (length ys) (unScore . score $ head ys)) . ( if explore then map (:[]) - else ( sortBy (comparing (ropt . score . head)) + else ( sortBy (comparing (unScore . score . head)) . groupBy ((==) `on` candidate) . sortBy (comparing candidate) ) diff --git a/changelog b/changelog --- a/changelog +++ b/changelog @@ -1,7 +1,20 @@ +0.1.1.0 + +- IUPAC nomenclature for sequence constraints +- --showmanual will now show README.md, while --help shows shorter help + 0.1.0.0 - * uses new ViennaRNA bindings - * added correct name +- major cleanup of source +- preparation for MCMC library transition +- small typos fixed + +0.1.0.0 + +- uses new ViennaRNA bindings +- added correct name + 0.0.2.1 - * post-publication version - * allows continuous Markovian walk for special applications + +- post-publication version +- allows continuous Markovian walk for special applications