diff --git a/MutationOrder.cabal b/MutationOrder.cabal --- a/MutationOrder.cabal +++ b/MutationOrder.cabal @@ -1,5 +1,5 @@ name: MutationOrder -version: 0.0.0.1 +version: 0.0.0.2 author: Maria Beatriz Walter Costa, Christian Hoener zu Siederdissen, 2017 copyright: Christian Hoener zu Siederdissen, 2017 homepage: https://github.com/choener/MutationOrder @@ -111,6 +111,7 @@ build-depends: base , bytestring , cmdargs >= 0.10 + , file-embed >= 0.0.8 , filepath -- , MutationOrder @@ -119,6 +120,7 @@ default-extensions: BangPatterns , DeriveDataTypeable , RecordWildCards + , TemplateHaskell main-is: MutationOrder.hs default-language: diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -2,32 +2,94 @@ Determine the most likely order of mutations from one RNA sequence to another. -1. Walter Costa, Maria Beatriz and Hoener zu Siederdissen, Christian and Tulpan, Dan and Stadler, Peter F. and Nowick, Katja + Walter Costa, Maria Beatriz and Hoener zu Siederdissen, Christian and Tulpan, Dan and Stadler, Peter F. and Nowick, Katja *Uncovering the Structural Evolution of the Human Accelerated Region 1* 2017, submitted [preprint](http://www.bioinf.uni-leipzig.de/~choener/pdfs/wal-hoe-2017.pdf) +# General information + +Given two RNA sequences, one ancestral, and one extant, we want to determine +the most likely path of evolution under different measures of fitness. + +This program produces the (i) maximum-likelihood path, (ii) all end +probabilities, (iii) all start-end probabilities, (iv) all edge probabilities, +and (v) the maximum expected accuracy path for these two RNA sequences. + +In detail: +(i) gives the optimal path(s) for the fitness function +(ii) gives for each nucleotide polymorphism, how likely it is, that this mutation was introduced *last* +(iii) looks at all pairs of (first mutation, last mutation) and gives the probability that these two mutations are the begin and end of the chain of mutations +(iv) yields for all pairs of nodes (i -> j) the probability that this path occurs, over the whole ensemble of all possible paths +(v) produces the path of maximal weight using the probabilities produced in (iv) + # Usage instructions -We assume that you have two Fasta files, *from.fa* and *to.fa* but they can be -named however is convenient (say *chimp.fa* and *human.fa*). Each file has to -contain exactly one sequence and both sequences have to be of the same length. +## example usage +We assume that you have two Fasta files, *chimp_118.fa* and *human_118.fa* but +they can be named however is convenient. Each file has to contain exactly one +sequence and both sequences have to be of the same length. + +For testing with chimp and human, the provided chimp-human.json.gz database +should be used, otherwise the initial foldings will be recalculated. All +required files are available under 'Binaries' at the bottom of the page. + +In case, you don't want or can't use the provided work database, run +./MutationOrder with --verbose + We then run -```./MutationOrder --workdb from-to.json.gz --scoretype pairdistcen --onlypositive --outputprefix test``` + ./MutationOrder --workdb chimp-human.json.gz --scoretype pairdistcen --onlypositive --outputprefix test chimp_118.fa human_118.fa This will generate ```test.run```, ```test-edge.eps```, and -```test-meaorder.eps```. The ```test.run``` file provides extensive output of -the optimal path, the first-last probabilities, the edge probabilities, and the -mea output. The two ```eps``` files give a graphical representation of the edge +```test-meaorder.eps```. + +The ```test.run``` file provides extensive output of the optimal path, the +first-last probabilities, the edge probabilities, and the mea output. This +conforms to (i) -- (v) mentioned above. + +The two ```eps``` files give a graphical representation of the edge probabilities, for the ```meaorder``` in order of the path of maximum expected accuracy. -The ```--scoretype``` allows for ```mfe```, ```centroid```, ```pairdistcen```, -and ```pairdistmfe```, which analyse possible evoluationary paths according to -mfe energy, centroid energy, smallest base pair distances for each step in the -```cen```troid or ```mfe``` case. +The work database collects intermediate structures and their folding and is +only created once. The initial run will, however, take some time. I.e. for +'HAR1' this requires 1-4 hours depending on the machine. Further runs complete +*much* faster. In minutes for HAR1. + +## Command-line options + + --help provides short help + --verbose will show folding steps during the initial run + + -w + --workdb=ITEM the database where to store intermediate foldings + -t + --temperature=NUM annealing temperature. Values close to 0 favor optimal paths. The default is 1.0 + --fillweight=FILLWEIGHT provides logarithmic and linear fill styles for probability plots. The full style always fills the box + --fillstyle=FILLSTYLE normally, boxes are sized, but all in the same color. This changes the opacity of the color as well. Does not work well for eps files + --cooptcount=INT how many co-optimals to count for (the count in the .run file is produced differently) + --cooptprint=INT how many co-optimals to actually print + --outprefix=ITEM how to prefix all output files + --scoretype=SCORETYPE choose 'mfe', 'centroid', 'pairdistmfe', or 'pairdistcen' for the evaluation of each mutational step + --positivesquared square positive energies to penalize bad moves + --onlypositive minimize only over penalties, not energy gains + --equalstart each possible mutation is selected with equal probability as the initial one + --posscaled=NUM,NUM in =x,y will exponentiate all numbers >=x by the constant y. For value k>=x, we have k^y + --lkupfile=ITEM developer option to feed the initial work database with known foldings (usable but very raw and undocumented. needs 5-line rnafold output) + --showmanual will show this manual + +The allowed score types are: + + mfe +which optimizes based on the minimum free energy of each intermediate sequence + centroid +which instead looks at the energy of the centroid structure + pairdistmfe +which minimizes the base pair distance between following mutations using mfe structures + pairdistcen +which minimizes the base pair distance between following mutations using centroid structures diff --git a/changelog.md b/changelog.md --- a/changelog.md +++ b/changelog.md @@ -1,3 +1,9 @@ +0.0.0.2 +------- + +- updated readme +- included readme as manual + 0.0.0.1 ------- diff --git a/src/MutationOrder.hs b/src/MutationOrder.hs --- a/src/MutationOrder.hs +++ b/src/MutationOrder.hs @@ -3,10 +3,12 @@ module Main where -import System.Console.CmdArgs -import System.FilePath +import Control.Monad +import Data.FileEmbed import qualified Data.ByteString.Char8 as BS -import Control.Monad +import System.Console.CmdArgs +import System.Exit (exitSuccess, exitFailure) +import System.FilePath import BioInf.MutationOrder import BioInf.MutationOrder.RNA (createRNAlandscape) @@ -19,7 +21,7 @@ deriving (Show,Data,Typeable) data Options - = Options + = MutationOrder { infiles :: [FilePath] , workdb :: FilePath , temperature :: Double @@ -34,27 +36,29 @@ , equalStart :: Bool , posscaled :: Maybe (Double,Double) , lkupfile :: Maybe FilePath + , showmanual :: Bool } | GenSequences { infiles :: [FilePath] } deriving (Show,Data,Typeable) -oOptions = Options +oMutationOrder = MutationOrder { infiles = def &= args , workdb = "work.db" &= help "name of the database to store intermediates in" , temperature = 1.0 &= help "lower temperatures favor the more optimal paths, defaults to 1.0" - , fillweight = FWlog - , fillstyle = FSfull - , cooptcount = 100000 - , cooptprint = 2 - , outprefix = "tmp" + , fillweight = FWlog &= help "scale options for probability plots: fwlog, fwlinear, fwfill" + , fillstyle = FSfull &= help "fill options for probability plots: fsopacitylog, fsopacitylinear, fsfull" + , cooptcount = 1000 &= help "how many co-optimals to count" + , cooptprint = 2 &= help "how many co-optimals to actually print out" + , outprefix = "tmp" &= help "prefix for output files" , scoretype = Centroid &= help "choose 'mfe', 'centroid', 'pairdistmfe', or 'pairdistcen' for the evaluation of each mutational step" , positivesquared = False &= help "square positive energies to penalize worse structures" , onlypositive = False &= help "minimize only over penalties, not energy gains" - , equalStart = False - , posscaled = Nothing - , lkupfile = Nothing + , equalStart = False &= help "run mea with equal start probabilities" + , posscaled = Nothing &= help "--posscaled=x,y scale all values >= x by using y as exponent" + , lkupfile = Nothing &= help "developer option: if an RNAfold file with foldings exists, then use it" + , showmanual = False &= help "shows the manual" } oGenSequences = GenSequences @@ -63,9 +67,9 @@ main :: IO () main = do - o <- cmdArgs $ modes [oOptions, oGenSequences] &= verbosity + o <- cmdArgs $ modes [oMutationOrder &= auto, oGenSequences] &= verbosity case o of - Options{} -> mainProgram o + MutationOrder{} -> mainProgram o GenSequences{} -> genSequences o genSequences o = do @@ -76,8 +80,17 @@ forM_ ls $ \(k,sq) -> BS.putStrLn sq return () +embeddedManual = $(embedFile "README.md") + mainProgram oOptions = do - let Options{..} = oOptions + let MutationOrder{..} = oOptions + when showmanual $ do + BS.putStrLn embeddedManual + exitSuccess + unless (length infiles == 2) $ do + BS.putStrLn embeddedManual + putStrLn "\n\n\nThis program expects exactly two equal-length fasta files as input" + exitFailure isL <- isLoud let fwdScaleFunction = (if positivesquared then squaredPositive else id)