RNAdesign 0.1.0.0 → 0.1.1.0
raw patch · 10 files changed
+329/−315 lines, 10 filesdep +bytestringdep +file-embeddep ~BiobaseXNAdep ~PrimitiveArray
Dependencies added: bytestring, file-embed
Dependency ranges changed: BiobaseXNA, PrimitiveArray
Files
- BioInf/RNAdesign.hs +30/−175
- BioInf/RNAdesign/Assignment.hs +9/−9
- BioInf/RNAdesign/CandidateChain.hs +96/−0
- BioInf/RNAdesign/Graph.hs +2/−2
- BioInf/RNAdesign/LogMultinomial.hs +1/−0
- BioInf/RNAdesign/OptParser.hs +4/−4
- README.md +89/−0
- RNAdesign.cabal +13/−31
- RNAdesign.hs +68/−90
- changelog +17/−4
BioInf/RNAdesign.hs view
@@ -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
BioInf/RNAdesign/Assignment.hs view
@@ -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
+ BioInf/RNAdesign/CandidateChain.hs view
@@ -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+
BioInf/RNAdesign/Graph.hs view
@@ -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
BioInf/RNAdesign/LogMultinomial.hs view
@@ -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"
BioInf/RNAdesign/OptParser.hs view
@@ -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] ]
README.md view
@@ -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
RNAdesign.cabal view
@@ -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:
RNAdesign.hs view
@@ -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) )
changelog view
@@ -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