/PhyloSim.R
R | 6741 lines | 2944 code | 500 blank | 3297 comment | 491 complexity | 599a6b1f6c5fd46221d5d716d177dc82 MD5 | raw file
Possible License(s): GPL-3.0
- ##
- ## Copyright 2009 Botond Sipos
- ## See the package description for licensing information.
- ##
- ##########################################################################/**
- #
- # @RdocClass PhyloSim
- # \alias{phylosim}
- #
- # @title "The PhyloSim class"
- #
- # \description{
- #
- # PhyloSim is an extensible object-oriented framework for the Monte Carlo simulation
- # of sequence evolution written in 100 percent \code{R}.
- # It is built on the top of the \code{\link{R.oo}} and \code{\link{ape}} packages and uses
- # Gillespie's direct method to simulate substitutions, insertions and deletions.
- #
- # Key features offered by the framework:
- # \itemize{
- # \item Simulation of the evolution of a set of discrete characters with arbitrary states evolving
- # by a continuous-time Markov process with an arbitrary rate matrix.
- # \item Explicit implementations of the most popular substitution models (for nucleotides, amino acids and codons).
- # \item Simulation under the popular models of among-sites rate variation, like the gamma (+G) and invariants plus gamma (+I+G) models.
- # \item The possibility to simulate with arbitrarily complex patterns of among-sites rate variation by setting the site specific rates according to any \code{R} expression.
- # \item Simulation with one or more separate insertion and/or deletion processes acting on the sequences, which sample the insertion/deletion length from an arbitrary discrete distribution or an \code{R} expression (so all the probability distributions implemented in \code{R} are readily available for this purpose).
- # \item Simulation of the effects of variable functional constraints over the sites by site-process-specific insertion and deletion tolerance parameters, which determine the rejection probability of a proposed insertion/deletion.
- # \item The possibility of having a different set of processes and site-process-specific parameters for every site, which allow for an arbitrary number of partitions in the simulated data.
- # \item Simulation of heterotachy and other cases of non-homogeneous evolution by allowing the user to set "node hook" functions altering the site properties at internal nodes of the phylogeny.
- # \item The possibility to export the counts of various events ("branch statistics") as phylo objects (see \code{\link{exportStatTree.PhyloSim}}).
- # }
- #
- # General notes:
- # \itemize{
- # \item The \code{Sequence} objects have no "immortal links". The simulation
- # is aborted if the sequence length shrinks to zero. It is up to the user
- # to choose sensible indel rates and sequence lengths to prevent that.
- # \item The sites near the beginning and end of the sequences have less sites proposing
- # insertion and deletion events around the so the insertion and deletion processes
- # have an "edge effect". The user can simulate
- # realistic flanking sequences to alleviate the edge effect in the simulation settings where
- # it may be an issue.
- # }
- #
- # Notes on performance:
- # \itemize{
- # \item The pure \code{R} implementation offers felxibility, but also comes
- # with a slower simulation speed. If the \code{PSIM_FAST} object is present in the environment, a "fast & careless"
- # mode is enabled. In this mode most of the error checking is skipped, increasing the speed.
- # It is recomended that simulations are only run in fast mode if you are sure that the simulation
- # settings are free from errors. It is probably a good practice to set up the simulations in normal mode
- # with short sequences and enable fast mode when running the actual simulation with long sequences.
- # \item Please note, that no "branch statistics" are saved in fast mode.
- # \item Logging also has a negative impact on performance, so it's not a good idea to run
- # large simulations with the logging enabled.
- # \item The time needed to run a simulation depends not only on the number of the sites,
- # but also on the length of the tree.
- # \item Constructing \code{Sequence} objects with large number of sites is expensive. Avoid doing
- # that inside a cycle.
- # \item In the case of \code{Sequence} objects with a large number of sites (more than 10 000) the
- # amount of available memory can be limiting as well.
- # }
- #
- # The examples below demonstrate only some more common simulation settings,
- # the framework offers much more flexibility. See the package
- # vignette (\code{vignette("PhyloSim",package="phylosim")}) and the
- # examples directory (\url{http://github.com/sbotond/phylosim/tree/master/examples/})
- # for additional examples.
- #
- # @classhierarchy
- # }
- #
- #
- # \references{
- # Gillespie, DT (1977) Exact stochastic simulation of coupled chemical reactions -
- # J. Phys. Chem. 81 (25):2340-2361 \url{http://dx.doi.org/10.1021/j100540a008}
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{phylo}{A rooted phylo object, constructed by the APE package.}
- # \item{root.seq}{A valid Sequence object with Process objects attached. Used as the starting sequence during simulation.}
- # \item{name}{The name of the object (a character vector of length one).}
- # \item{log.file}{Name of the file used for logging.}
- # \item{log.level}{An integer specifying the verbosity of logging (see \code{\link{setLogLevel.PhyloSim}}).}
- # \item{...}{Not used.}
- # }
- #
- # \section{Fields and Methods}{
- # @allmethods
- # }
- #
- # \examples{
- # set.seed(1)
- # ## The following examples demonstrate
- # ## the typical use of the framework.
- # ## See the package vignette and
- # ## \url{http://github.com/sbotond/phylosim/tree/master/examples/}
- #
- # ## The ll() method gives information about the methods defined
- # ## in the immediate class of an object.
- # ## Useful when exploring the framework.
- #
- # s<-Sequence()
- # ll(s)
- # ll(PhyloSim())
- # ll(GTR())
- #
- # ## Example 1 - A short simulation:
- # ## simulate nucleotide seqeunces and display
- # ## the resulting alignment matrix.
- #
- # Simulate(
- # PhyloSim(phy=rcoal(3),
- # root=NucleotideSequence(string="ATGC", proc=list(list(JC69())) ) )
- # )$alignment
- #
- # # Construct a phylo object for the following
- # # simulations, scale total tree length to 1:
- #
- # tmp<-PhyloSim(phylo=rcoal(3))
- # scaleTree(tmp,1/tmp$treeLength)
- # tmp$treeLength
- # t<-tmp$phylo
- #
- # ## Example 3 - simulating rate variation,
- # ## insertions and deletions.
- # ## See the examples/example_3_clean.R file
- # ## in the phylosim GitHub repository.
- #
- # # construct a GTR process object
- # gtr<-GTR(
- # name="MyGTR",
- # rate.params=list(
- # "a"=1, "b"=2, "c"=3,
- # "d"=1, "e"=2, "f"=3
- # ),
- # base.freqs=c(2,2,1,1)/6
- # )
- # # get object summary
- # summary(gtr)
- # # display a bubble plot
- # plot(gtr)
- #
- # # construct root sequence object
- # s<-NucleotideSequence(length=20)
- # # attach process via virtual field
- # s$processes<-list(list(gtr))
- # # sample states from the equilibrium
- # # distribution of the attached processes
- # sampleStates(s)
- # # create among-site rate variation by sampling
- # # the "rate.multiplier" site-process-specific parameter
- # # from a discrete gamma distribution (GTR+G).
- # plusGamma(s,gtr,shape=0.1)
- # # make the range 11:12 invariable
- # setRateMultipliers(s,gtr,0,11:12)
- # # get the rate multipliers for s and gtr
- # getRateMultipliers(s,gtr)
- #
- # construct a deletion process object
- # # proposing lengths in the range 1:3
- # d<-DiscreteDeletor(
- # rate=0.1,
- # name="MyDel",
- # sizes=c(1:3),
- # probs=c(3/6,2/6,1/6)
- # )
- # # get object
- # summary(d)
- # # plot deletion length distribution
- # plot(d)
- # # attach deletion process d to sequence s
- # attachProcess(s,d)
- # # create a region rejecting all deletions
- # setDeletionTolerance(s,d,0,11:12)
- #
- # # construct an insertion process object
- # # proposing lengths in the range 1:3
- # i<-DiscreteInsertor(
- # rate=0.1,
- # name="MyDel",
- # sizes=c(1:2),
- # probs=c(1/2,1/2),
- # template.seq=NucleotideSequence(length=1,processes=list(list(JC69())))
- # )
- # # states will be sampled from the JC69 equilibrium distribution
- # # get object
- # summary(i)
- # # plot insertion length distribution
- # plot(i)
- # # attach insertion process i to sequence s
- # attachProcess(s,i)
- # # create a region rejecting all insertions
- # setInsertionTolerance(s,i,0,11:12)
- #
- # # plot total site rates
- # plot(s)
- # # construct simulation object
- # sim<-PhyloSim(root.seq=s, phylo=t)
- # # get object summary
- # summary(sim)
- # # plot tree
- # plot(sim)
- # # run simulation
- # Simulate(sim)
- # # get the list of recorded per-branch event counts
- # getBranchEvents(sim)
- # # export the number of substitutions as a phylo object
- # subst<-exportStatTree(sim,"substitution")
- # # plot the exported phylo object
- # plot(subst)
- # # plot tree and alignment
- # plot(sim)
- # # save and display alingment
- # file<-paste("PhyloSim_dummy_fasta_",Sys.getpid(),".fas",sep="");
- # saveAlignment(sim,file=file);
- # # print out the Fasta file
- # cat(paste(scan(file=file,what=character(),sep="\n"),collapse="\n"));cat("\n");
- # # delete Fasta file
- # unlink(file);
- #
- # # See \url{http://github.com/sbotond/phylosim/tree/master/examples/examples_class.R}
- # # for the full list of PhyloSim constructor examples.
- #
- # ## See the package vignette and
- # ## the GitHub repository for even more examples.
- # }
- #
- # @author
- #
- # \seealso{
- # \code{\link{PSRoot} \link{Alphabet} \link{AminoAcidAlphabet}
- # \link{AminoAcidSequence} \link{AminoAcidSubst}
- # \link{AnyAlphabet} \link{BinaryAlphabet} \link{BinarySequence}
- # \link{BinarySubst} \link{BrownianInsertor} \link{CodonAlphabet}
- # \link{CodonSequence} \link{CodonUNREST} \link{ContinuousDeletor}
- # \link{ContinuousInsertor} \link{cpREV} \link{DiscreteDeletor}
- # \link{DiscreteInsertor} \link{Event} \link{F81} \link{F84}
- # \link{FastFieldDeletor} \link{GeneralDeletor}
- # \link{GeneralInDel} \link{GeneralInsertor} \link{GeneralSubstitution}
- # \link{GTR} \link{GY94} \link{HKY} \link{JC69} \link{JTT} \link{JTT.dcmut}
- # \link{K80} \link{K81} \link{LG} \link{mtArt} \link{mtMam} \link{mtREV24}
- # \link{MtZoa} \link{NucleotideAlphabet} \link{NucleotideSequence} \link{PAM}
- # \link{PAM.dcmut} \link{PhyloSim} \link{Process} \link{QMatrix} \link{Sequence}
- # \link{Site} \link{T92} \link{TN93} \link{UNREST} \link{WAG}}
- # }
- #
- #*/###########################################################################
- setConstructorS3(
- "PhyloSim",
- function(
- phylo=NA,
- root.seq=NA,
- name=NA,
- log.file=NA,
- log.level=-1, # no loggin is performed by default
- ...
- ) {
- this<-PSRoot();
- this<-extend(this,
- "PhyloSim",
- .name="Anonymous",
- .phylo=NA,
- .root.sequence=NA,
- .sequences=list(), # references to the sequence objects
- .node.hooks=list(), # references to the node hook functions.
- .branch.stats=list(), # branch statistics.
- .alignment=NA, # the resulting alignment in fasat format.
- .log.file=NA, # the name of the log file.
- .log.connection=NA, # connection for the log file.
- .log.level=NA # log level
- );
- if(!all(is.na(phylo))){
- this$phylo<-phylo;
- }
- if(!all(is.na(root.seq))){
- this$rootSeq<-root.seq;
- }
-
- if(!missing(name)){
- this$name<-name;
- }
- if(!missing(log.file)){
- this$logFile<-log.file;
- } else {
- # Setting default log file:
- tmp<-this$id;
- tmp<-gsub(":","_",tmp);
- this$logFile<-paste(tmp,".log",sep="");
- }
- # Setting log level:
- this$logLevel<-log.level;
- return(this);
- },
- enforceRCC=TRUE
- );
- ##
- ## Method: checkConsistency
- ##
- ###########################################################################/**
- #
- # @RdocMethod checkConsistency
- #
- # @title "Check object consistency"
- #
- # \description{
- # @get "title".
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{An object.}
- # \item{...}{Not used.}
- # }
- #
- #
- # \value{
- # Returns an invisible TRUE if no inconsistencies found in the object, throws
- # an error otherwise.
- # }
- #
- # @author
- #
- # \seealso{
- # @seeclass
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "checkConsistency",
- class="PhyloSim",
- function(
- this,
- ...
- ){
- may.fail<-function(this) {
-
- # Checking the name:
- this$name<-this$name;
- # Checking the phylo object:
- if (!any(is.na(this$.phylo)) & !is.phylo(this$.phylo) ){
- throw("The phylo object is invalid!\n");
- }
- # Checking the log level:
- if(!is.numeric(this$.log.level) | (length(this$.log.level) != 1) ){
- throw("The log level must be numeric vector of length 1!\n");
- }
- # Checking lof file:
- if(!is.character(this$.log.file) | (length(this$.log.level) != 1) ){
- throw("The log file must be charcter vector of length 1!\n");
- }
- # Checking the sequences:
- for (seq in this$.sequences){
- if(is.Sequence(seq)){
- checkConsistency(seq);
- }
- }
- # Checking node hooks:
- for (hook in this$.node.hooks){
- if(!is.null(hook) & !is.function(hook)){
- throw("Invalid node hook found!\n");
- }
- }
- # Checking the alignment:
- if(!any(is.na(this$.alignment))){
- .checkAlignmentConsistency(this, this$.alignment);
- }
-
- }
- tryCatch(may.fail(this));
- return(invisible(TRUE));
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- )
- ##
- ## Method: is.phylo.default
- ##
- ###########################################################################/**
- #
- # @RdocDefault is.phylo
- #
- # @title "Check if an object is an instance of the phylo class"
- #
- # \description{
- # @get "title".
- # Phylo objects are created by the \pkg{APE} package. This method just return the value of \code{inherits(this,"phylo")}.
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{An object.}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # TRUE or FALSE.
- # }
- #
- # \examples{
- # # load APE
- # library(ape);
- # # create some objects
- # o1<-Object();
- # o2<-rcoal(3);
- # # check if they are phylo objects
- # is.phylo(o1);
- # is.phylo(o2);
- #
- # }
- #
- # @author
- #
- # \seealso{
- # The \pkg{ape} package.
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "is.phylo",
- class="default",
- function(
- this,
- ...
- ){
- inherits(this,"phylo");
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- )
- ##
- ## Method: setPhylo
- ##
- ###########################################################################/**
- #
- # @RdocMethod setPhylo
- #
- # @title "Set the phylo object for a PhyloSim object"
- #
- # \description{
- # @get "title".
- #
- # The internal structure of the provided phylo object is reordered in a cladeweise fashion.
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{A PhyloSim object.}
- # \item{value}{A phylo object created by the \pkg{ape} package.}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # A phylo object or FALSE.
- # }
- #
- # \examples{
- # #create a PhyloSim object
- # sim<-PhyloSim();
- # # creat a phylo object
- # tree<-rcoal(3);
- # # get/set phylo object
- # setPhylo(sim,tree);
- # getPhylo(sim,tree);
- # # get/set phylo object via virtual field
- # sim$tree<-rcoal(5);
- # sim$tree;
- # }
- #
- # @author
- #
- # \seealso{
- # The PhyloSim class, the \pkg{ape} package.
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "setPhylo",
- class="PhyloSim",
- function(
- this,
- value,
- ...
- ){
- if(missing(value)){
- throw("No object provided!\n");
- }
- else if(!is.phylo(value)){
- throw("The new value must be a \"phylo\" object!\n");
- }
- else if(!is.rooted(value)){
- throw("The new value must be a rooted \"phylo\" object!\n");
- }
- else {
- .checkTipLabels(value);
- this$.phylo<-value;
- this$.phylo<-reorder(this$.phylo, order="cladewise");
- for (i in this$nodes){
- this$.sequences[[i]]<-NA;
- }
- return(this$.phylo);
- }
- return(FALSE);
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: .checkTipLabels
- ##
- setMethodS3(
- ".checkTipLabels",
- class="phylo",
- function(
- this,
- ...
- ){
- for(label in this$tip.label){
- if(length(grep("^Node \\d+$",label,perl=TRUE,value=FALSE)) > 0){
- throw("Sorry, but the node labels matching \"Node \\d+\" are reserved for internal nodes! Blaming label: ",label,".\n");
- }
- else if(length(grep("^Root node \\d+$",label,perl=TRUE,value=FALSE)) > 0){
- throw("Sorry, but the node labels matching \"Root node \\d+\" are reserved for the root node! Blaming label: ",label,".\n");
- }
-
- }
- },
- private=TRUE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- )
- ##
- ## Method: getPhylo
- ##
- ###########################################################################/**
- #
- # @RdocMethod getPhylo
- #
- # @title "Get the phylo object aggregated in a PhyloSim object"
- #
- # \description{
- # @get "title".
- #
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{A PhyloSim object.}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # A phylo object or NA.
- # }
- #
- # \examples{
- # #create a PhyloSim object
- # sim<-PhyloSim();
- # # creat a phylo object
- # tree<-rcoal(3);
- # # get/set phylo object
- # setPhylo(sim,tree);
- # getPhylo(sim,tree);
- # # get/set phylo object via virtual field
- # sim$tree<-rcoal(5);
- # sim$tree;
- # }
- #
- # @author
- #
- # \seealso{
- # The PhyloSim class, the \pkg{ape} package.
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "getPhylo",
- class="PhyloSim",
- function(
- this,
- ...
- ){
- this$.phylo;
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- )
- ##
- ## Method: setRootSeq
- ##
- ###########################################################################/**
- #
- # @RdocMethod setRootSeq
- #
- # @title "Set the root sequence for a PhyloSim object"
- #
- # \description{
- # @get "title".
- #
- # The root sequence will be used as a starting point for the simulation. The phylo object must be set before
- # trying to set the root sequence object.
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{A PhyloSim object.}
- # \item{value}{A valid Sequence object.}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # The root Sequence object if succesfull, FALSE otherwise.
- # }
- #
- # \examples{
- # # create some objects
- # sim<-PhyloSim(phylo=rcoal(3));
- # seq<-NucleotideSequence(string="ATGCC");
- # # set/get root sequence
- # setRootSeq(sim, seq);
- # getRootSeq(sim, seq);
- # # set/get root sequence via virtual field
- # sim$rootSeq<-BinarySequence(string="111000111000");
- # sim$rootSeq;
- #
- # }
- #
- # @author
- #
- # \seealso{
- # @seeclass Sequence Process
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "setRootSeq",
- class="PhyloSim",
- function(
- this,
- value,
- ...
- ){
- if(missing(value)){
- throw("No object provided!\n");
- }
- else if(!is.Sequence(value)){
- throw("The new value must be a sequence object!\n");
- }
- else {
- this$.root.sequence<-clone(value);
- this$.root.sequence$name<-paste("Root node",this$rootNode);
-
- # Call garbage collection:
- gc();
- gc();
- return(this$.root.sequence);
- }
- return(FALSE);
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- )
- ##
- ## Method: getRootSeq
- ##
- ###########################################################################/**
- #
- # @RdocMethod getRootSeq
- #
- # @title "Get the root sequence aggregated by a PhyloSim object"
- #
- # \description{
- # @get "title".
- #
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{A PhyloSim object.}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # The root Sequence object or NA.
- # }
- #
- # \examples{
- # # create some objects
- # sim<-PhyloSim(phylo=rcoal(3));
- # seq<-NucleotideSequence(string="ATGCC");
- # # set/get root sequence
- # setRootSeq(sim, seq);
- # getRootSeq(sim, seq);
- # # set/get root sequence via virtual field
- # sim$rootSeq<-BinarySequence(string="111000111000");
- # sim$rootSeq;
- #
- # }
- #
- # @author
- #
- # \seealso{
- # @seeclass Sequence Process
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "getRootSeq",
- class="PhyloSim",
- function(
- this,
- ...
- ){
- this$.root.sequence;
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- )
- ##
- ## Method: as.character.PhyloSim
- ##
- ###########################################################################/**
- #
- # @RdocMethod as.character
- #
- # @title "Return the character representation of a PhyloSim object"
- #
- # \description{
- # @get "title".
- #
- # The character representation is the identifier of the PhyloSim object as returned by the \code{getId} method.
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{x}{A PhyloSim object.}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # A character vector of length one.
- # }
- #
- # \examples{
- # # create a PhyloSim object
- # o<-PhyloSim(name="MySim");
- # # get character representation
- # as.character(o)
- # }
- #
- # @author
- #
- # \seealso{
- # @seeclass
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "as.character",
- class="PhyloSim",
- function(
- x,
- ...
- ){
- return(getId(x));
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- )
- ##
- ## Method: getName
- ##
- ###########################################################################/**
- #
- # @RdocMethod getName
- #
- # @title "Get the name of a PhyloSim object"
- #
- # \description{
- # @get "title".
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{A PhyloSim object.}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # A character vector of length one.
- # }
- #
- # \examples{
- # # create a PhyloSim object
- # o<-PhyloSim();
- # # set/get name
- # setName(o,"MySim");
- # getName(o,"MySim");
- # # set/get name via virtual field
- # o$name<-"George";
- # o$name
- # }
- #
- # @author
- #
- # \seealso{
- # @seeclass
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "getName",
- class="PhyloSim",
- function(
- this,
- ...
- ){
- this$.name;
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: setName
- ##
- ###########################################################################/**
- #
- # @RdocMethod setName
- #
- # @title "Set the name of a PhyloSim object"
- #
- # \description{
- # @get "title".
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{A PhyloSim object.}
- # \item{new.name}{A character vector of length one.}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # The new name.
- # }
- #
- # \examples{
- # # create a PhyloSim object
- # o<-PhyloSim();
- # # set/get name
- # setName(o,"MySim");
- # getName(o,"MySim");
- # # set/get name via virtual field
- # o$name<-"George";
- # o$name
- # }
- #
- # @author
- #
- # \seealso{
- # @seeclass
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "setName",
- class="PhyloSim",
- function(
- this,
- new.name,
- ...
- ){
- this$.name<-as.character(new.name);
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: getId
- ##
- ###########################################################################/**
- #
- # @RdocMethod getId
- #
- # @title "Get the unique identifier of a PhyloSim object"
- #
- # \description{
- # @get "title".
- # The unique identifier is the concatenation of the class, the object name as returned by getName() and the object hash
- # as returned by hashCode().
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{A PhyloSim object.}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # A character vector of length one.
- # }
- #
- # \examples{
- # # create a PhyloSim object
- # o<-PhyloSim(name="MySim");
- # # get id
- # getId(o);
- # # get id via virtual field
- # o$id;
- # }
- #
- # @author
- #
- # \seealso{
- # @seeclass
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "getId",
- class="PhyloSim",
- function(
- this,
- ...
- ){
- this.class<-class(this)[1];
- id<-paste(this.class,this$.name,hashCode(this),sep=":");
- return(id);
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: setId
- ##
- ###########################################################################/**
- #
- # @RdocMethod setId
- #
- # @title "Forbidden action: setting the unique identifier of a PhyloSim object"
- #
- # \description{
- # @get "title".
- #
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{An object.}
- # \item{value}{Not used.}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # Throws an error.
- # }
- #
- # @author
- #
- # \seealso{
- # @seeclass
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "setId",
- class="PhyloSim",
- function(
- this,
- value,
- ...
- ){
- throw("Id is generated automatically and it cannot be set!\n");
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: Simulate
- ##
- ###########################################################################/**
- #
- # @RdocMethod Simulate
- #
- # @title "Run a simulation according to a PhyloSim object"
- #
- # \description{
- # @get "title".
- #
- # The phylo object and the root sequence must be set before attempting to run a simulation.
- # Also the bigRate of the root sequence must not be NA or zero, so at least one sane
- # Process object must be attached to the root sequence object.
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{A PhyloSim object.}
- # \item{quiet}{TRUE or FALSE (default).}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # The PhyloSim object (invisible).
- # }
- #
- # \examples{
- # # Create a PhyloSim object.
- # # Provide the phylo object
- # # and the root sequence.
- # sim<-PhyloSim(
- # name="TinySim",
- # phylo=rcoal(3),
- # root.seq=NucleotideSequence(string="ATGC",processes=list(list(JC69())))
- # );
- # # Run the simulation
- # Simulate(sim);
- # # Print the resulting sequences
- # sim$sequences
- # # Print the resulting alignment
- # sim$alignment
- # }
- #
- # @author
- #
- # \seealso{
- # @seeclass
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "Simulate",
- class="PhyloSim",
- function(
- this,
- quiet=FALSE,
- ...
- ){
- if(!is.phylo(this$.phylo)){
- throw("Cannot simulate because the phylo object is not set or it is invalid!\n");
- }
- # Check for the root sequence:
- else if(!is.Sequence(this$.root.sequence)){
- throw("Cannot simulate because the root sequence is not set or it is invalid!\n");
- }
- # Check bigRate validity:
- else if(is.na(this$.root.sequence$bigRate)){
- throw("Cannot simulate because the bigRate of the root sequence is NA!\n");
- }
- else{
- # Warn for zero bigRate:
- if(this$.root.sequence$bigRate == 0){
- warning("The bigRate of the root sequence is zero! You are running a pointless simulation!\n");
- }
- if(exists(x="PSIM_FAST")){
- if(!quiet){
- cat("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
- cat("!! WARNING: fast & careless mode is on, most of the error checking is omitted! !!\n");
- cat("!! Please note that this also disables the saving of branch statistics. !!\n");
- cat("!! You can go back to normal mode by deleting the PSIM_FAST object. !!\n");
- cat("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
- }
- Log(this,"WARNING: fast & careless mode is on, most of the error checking is omitted!");
- }
-
- # Attach root sequence to root node:
- Log(this,paste("Attaching root sequence ",this$.root.sequence$id,sep=""));
- attachSeqToNode(this, node=getRootNode(this),seq=this$.root.sequence);
- # Write protecting the root sequence:
- Log(this,paste("Write protecting root sequence ",this$.root.sequence$id,sep=""));
- this$.root.sequence$writeProtected<-TRUE;
- # Traverse the tree and simulate:
- Log(this,paste("Starting simulation on the object",this$id));
- edge.counter<-1;
- n.edges<-this$nedges;
- for(edge in 1:n.edges){
- if(!quiet){ cat("Simulating edge",edge,"of", n.edges,"\n");}
- Log(this,paste("Starting to simulate edge",edge,"of",n.edges));
- .simulateEdge(this,number=edge);
- edge.counter<-edge.counter+1;
- }
- }
- Log(this, "Simulation finished, building alignment!\n");
- this$.alignment<-.recoverAlignment(this);
- # Flush the log connection:
- if(!is.na(this$.log.connection)){
- close(this$.log.connection);
- }
- # Call the garbage collector:
- gc();
- gc();
- return(invisible(this));
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: .simulateEdge
- ##
- setMethodS3(
- ".simulateEdge",
- class="PhyloSim",
- function(
- this,
- number=NA,
- ...
- ){
- # Get edge:
- edge<-getEdge(this, number);
- # Get parent node:
- start.seq<-getSeqFromNode(this, edge[[1,"from"]]);
- # Evolve sequence:
- new.seq<-.evolveBranch(this, start.seq=start.seq, branch.length=edge[1,"length"], old.node=edge[[1,"from"]],new.node=edge[[1,"to"]], branch.number=number);
- # Write protect the sequence:
- new.seq$writeProtected<-TRUE;
- # Attach sequence to children node:
- attachSeqToNode(this, node=edge[1,"to"], seq=new.seq);
- },
- private=TRUE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: .evolveBranch
- ##
- setMethodS3(
- ".evolveBranch",
- class="PhyloSim",
- function(
- this,
- start.seq=NA,
- branch.length=NA,
- old.node=NA,
- new.node=NA,
- branch.number=NA,
- ...
- ){
- if(!exists(x="PSIM_FAST")){
- if(missing(start.seq)){
- throw("No starting sequence provided!\n");
- }
- else if(missing(branch.length)){
- throw("No branch length provided!\n");
- }
- else if(!is.numeric(branch.length)){
- throw("The branch length must be numeric!\n");
- }
- }
- if(.checkSeq(this, start.seq) ){
-
- # Cloning the starting sequence:
- seq<-clone(start.seq);
-
- # Set the name of the sequence object:
- if(is.tip(this, new.node)){
- seq$name<-this$tipLabels[[new.node]];
- }
- else {
- seq$name<-paste("Node",new.node);
- }
- .GillespieDirect(this, seq=seq, branch.length=branch.length, branch.number=branch.number);
-
- # Call the node hook if exists:
- hook<-this$.node.hooks[[as.character(new.node)]];
- if(!is.null(hook) & is.function(hook)){
- Log(this,paste("Calling node hook for node",new.node));
- seq<-hook(seq=seq);
- if(!is.Sequence(seq)){
- throw("Node hook returned an invalid sequence object!\n");
- }
- else if(is.na(seq$bigRate)){
- throw("Node hook returned sequence with NA bigRate!\n");
- }
- else if(seq$bigRate == 0.0){
- throw("Node hook returned sequence with zero bigRate!\n");
- }
- else{
- checkConsistency(seq, omit.sites=TRUE);
- }
- }
- # Return the resulting sequence object:
- return(seq);
-
- }
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: .GillespieDirect
- ##
- setMethodS3(
- ".GillespieDirect",
- class="PhyloSim",
- function(
- this,
- seq=NA,
- branch.length=NA,
- branch.number=NA,
- ...
- ){
- Debug(this, paste("Branch length is",branch.length));
-
- # Initialize time:
- time<-0.0;
-
- # Sample the next waiting time until
- # the branch length is consumed:
- while( (time<-time + rexp(n=1, rate=(big.rate<-getBigRate(seq)))) <= branch.length){
- # Generate a random number between zero and the bigRate:
- E<-runif(n=1,min=0,max=big.rate);
- # Identify the target site:
- site.number<-which(.getCumulativeRatesFast(seq) >= E)[[1]];
- # Get the events from the target site:
- site<-seq$.sites[[site.number]];
- site$.position<-site.number;
- events<-getEvents(site);
- site$.position<-NULL;
-
- # Get the rates:
- rates<-double();
- for(e in events){
- rates<-c(rates,e$.rate);
- }
-
- # Calculate the corresponding cumulative rates:
- if(site.number > 1){
- rates<-cumsum(c(seq$.cumulative.rates[[site.number - 1]], rates));
- }
- else {
- rates<-cumsum(c(0.0, rates));
- }
- # Pick the event:
- event.number<-which(rates >= E)[[1]] - 1;
- event<-events[[event.number]];
- # Log the event:
- Log(this,paste("Performing event [",event$.name,"] at position",event$.position,"generated by the process",event$.process$.id));
- # Perform the event:
- event.details<-Perform(event);
- Debug(this,paste("Remaining branch length is",(branch.length-time) ));
-
- # Log event details:
-
- # Log deletion event details:
- if(event$.name == "Deletion"){
- Log(this,paste("The process",event$.process,"proposed to delete range",paste(event.details$range,collapse="--"),". Accepted:",event.details$accepted));
- }
- # Log insertion event details:
- else if(event$.name == "Insertion"){
- message<-paste("The process ",event$.process," proposed insertion at position ",event.details$position,". Accepted: ",event.details$accepted,sep="");
- if(event.details$accepted == TRUE){
- message<-paste(message,"."," Insert length was ",event.details$length,sep="");
- }
- Log(this, message);
- }
- # Update branch statistics if not in fast mode:
- if(!exists(x="PSIM_FAST")){
- .UpdateBranchStats(this,event,event.details, branch.number);
- }
- # Abort if sequence length shrunk to zero:
- if(seq$.length == 0){
- message<-paste("Terminating the simulation because the length of the sequence ",seq$name," shrunk to zero! Please be more careful when tuning the indel rates!\n");
- Log(this, message);
- throw(message);
- }
- } #/while
- # Calling garbage collection:
- gc();
- gc();
- return(seq);
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: attachSeqToNode
- ##
- ###########################################################################/**
- #
- # @RdocMethod attachSeqToNode
- #
- # @title "Assotiate a Sequence object with a given node of a phylo object aggregated by a PhyloSim object"
- #
- # \description{
- # @get "title".
- #
- # This method is mainly used internally.
- #
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{A PhyloSim object.}
- # \item{node}{Node identifier.}
- # \item{seq}{A Sequence object.}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # The PhyloSim object (invisible).
- # }
- #
- #
- # @author
- #
- # \seealso{
- # @seeclass
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "attachSeqToNode",
- class="PhyloSim",
- function(
- this,
- node=NA,
- seq=NA,
- ...
- ){
- if(!is.phylo(this$.phylo)){
- throw("The phylo object is not set, sequence to node is not possible!\n");
- }
- if(missing(node)){
- throw("No node specified!\n");
- }
- else if(missing(seq)){
- throw("No sequence object given");
- }
- else if(.checkNode(this,node) & .checkSeq(this, seq)){
-
- if(is.Sequence(this$.sequences[[node]])){
- throw("The node has already an attached sequence. Detach that before trying to attach a new one!\n");
- }
- else {
- this$.sequences[[as.numeric(node)]]<-seq;
- return(invisible(this));
- }
- }
-
- return(invisible(this));
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: attachHookToNode
- ##
- ###########################################################################/**
- #
- # @RdocMethod attachHookToNode
- #
- # @title "Attach a callback function to a given node of a phylo object aggregated by a PhyloSim object"
- #
- # \description{
- # @get "title".
- #
- # A "node hook" is a function which accepts a Sequence object through the named argument "seq" and returns a
- # Sequence object. The node hook function must accept any object which inherits from the \code{Sequence} class!
- #
- # After simulating the branch leading to the node, the resulting Sequence object is passed
- # to the node hook and the returned object is used to simulate the downstream branches.
- #
- # By using node hooks the attached processes can be replaced during simulation, hence enabling the simulation of
- # non-homogeneous sequence evolution.
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{A PhyloSim object.}
- # \item{node}{Node identifier.}
- # \item{fun}{A function (see above).}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # The PhyloSim object (invisible).
- # }
- #
- # \examples{
- # # Create a PhyloSim object.
- # # Provide the phylo object
- # # and the root sequence.
- # sim<-PhyloSim(
- # name="TinySim",
- # phylo=rcoal(3),
- # root.seq=NucleotideSequence(string="ATGC",processes=list(list(JC69())))
- # );
- # # create a node hook function
- # hook<-function(seq=NA){
- # # replace the substitution process with F84
- # if(inherits(seq,"NucleotideSequence")){
- # cat("Replacing JC69 with F84.\n");
- # seq$processes<-list(list(F84(rate.params=list("Kappa" = 2))));
- # }
- # return(seq);
- # }
- # # attach hook function to node 5
- # attachHookToNode(sim,5,hook);
- # # Run the simulation
- # Simulate(sim);
- # # Check if the processes have been truly replaced
- # lapply(sim$sequences, getUniqueProcesses.Sequence)
- # # Print the resulting alignment
- # sim$alignment
- # }
- #
- # @author
- #
- # \seealso{
- # @seeclass
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "attachHookToNode",
- class="PhyloSim",
- function(
- this,
- node=NA,
- fun=NA,
- ...
- ){
- if(!is.phylo(this$.phylo)){
- throw("The phylo object is not set, attaching node hook is not possible!\n");
- }
- if(missing(node)){
- throw("No node specified!\n");
- }
- else if(missing(fun)){
- throw("No function given!");
- }
- else if(!is.function(fun)){
- throw("The argument \"fun\" must be a function!\n");
- }
- else if( length(intersect(names(formals(fun)), "seq")) == 0 ){
- throw("The function argument must have a an argument named \"seq\"");
- }
- else if(!is.Sequence(fun(Sequence(length=1)))){
- throw("The insert hook function must return a Sequence object!\n");
- }
- else if( .checkNode(this,node) ){
- if(is.function(this$.node.hooks[[as.character(node)]])){
- throw("The node has already an attached node hook. Detach that before trying to attach a new one!\n");
- }
- else {
- this$.node.hooks[[as.character(node)]]<-fun;
- return(invisible(this));
- }
- }
-
- return(invisible(this));
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: .checkNode
- ##
- setMethodS3(
- ".checkNode",
- class="PhyloSim",
- function(
- this,
- node=NA,
- ...
- ){
- if(missing(node)){
- throw("No node specified!\n");
- } else if( length(intersect(node, getNodes(this))) != 1){
- throw("The specified node is invalid!\n");
- }
- else {
- return(TRUE);
- }
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: .checkSeq
- ##
- setMethodS3(
- ".checkSeq",
- class="PhyloSim",
- function(
- this,
- seq=NA,
- ...
- ){
- if(missing(seq)){
- throw("No sequence specified!\n");
- } else if(!is.Sequence(seq)){
- throw("The sequence object is invalid!\n");
- }
- else {
- return(TRUE);
- }
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: detachSeqFromNode
- ##
- ###########################################################################/**
- #
- # @RdocMethod detachSeqFromNode
- #
- # @title "Detach a Sequence object from a given node of a phylo object aggregated by a PhyloSim object"
- #
- # \description{
- # @get "title".
- #
- # This method is mainly used internally.
- #
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{A PhyloSim object.}
- # \item{node}{Node identifier.}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # The PhyloSim object (invisible).
- # }
- #
- # @author
- #
- # \seealso{
- # @seeclass
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "detachSeqFromNode",
- class="PhyloSim",
- function(
- this,
- node=NA,
- ...
- ){
- if(missing(node)){
- throw("No node specified!\n");
- }
- else if( .checkNode(this,node) ){
-
- this$.sequences[[as.numeric(node)]]<-NA;
- }
-
- return(invisible(this));
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: detachHookFromNode
- ##
- ###########################################################################/**
- #
- # @RdocMethod detachHookFromNode
- #
- # @title "Detach a node hook function from a given node of a phylo object aggregated by a PhyloSim object"
- #
- # \description{
- # @get "title".
- #
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{A PhyloSim object.}
- # \item{node}{Node identifier.}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # The PhyloSim object (invisible).
- # }
- #
- # \examples{
- # # Create a PhyloSim object.
- # # Provide the phylo object
- # # and the root sequence.
- # sim<-PhyloSim(
- # name="TinySim",
- # phylo=rcoal(3),
- # root.seq=NucleotideSequence(string="ATGC",processes=list(list(JC69())))
- # );
- # # create a node hook function
- # hook<-function(seq=NA){
- # # replace the substitution process with F84
- # if(inherits(seq,"NucleotideSequence")){
- # cat("Replacing JC69 with F84.\n");
- # seq$processes<-list(list(F84(rate.params=list("Kappa" = 2))));
- # }
- # return(seq);
- # }
- # # attach hook function to node 5
- # attachHookToNode(sim,5,hook);
- # # detach hook from node 5
- # detachHookFromNode(sim,5);
- # # Run the simulation again
- # Simulate(sim); # You should not see the message printed out by the "hook" function.
- #
- # }
- #
- # @author
- #
- # \seealso{
- # attachHookToNode PhyloSim Simulate.PhyloSim
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "detachHookFromNode",
- class="PhyloSim",
- function(
- this,
- node=NA,
- ...
- ){
- if(missing(node)){
- throw("No node specified!\n");
- }
- else if( .checkNode(this,node) ){
-
- this$.node.hooks[[as.character(node)]]<-NA;
- }
-
- return(invisible(this));
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: getSeqFromNode
- ##
- ###########################################################################/**
- #
- # @RdocMethod getSeqFromNode
- #
- # @title "Get the Sequence object associated with a given node of a phylo object aggregated by a PhyloSim object"
- #
- # \description{
- # @get "title".
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{A PhyloSim object.}
- # \item{node}{Node identifier.}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # A Sequence object.
- # }
- #
- # \examples{
- # # Create a PhyloSim object.
- # # Provide the phylo object
- # # and the root sequence.
- # sim<-PhyloSim(
- # name="TinySim",
- # phylo=rcoal(3),
- # root.seq=NucleotideSequence(string="ATG",processes=list(list(JC69())))
- # );
- # # get the sequence associated with node 5
- # getSeqFromNode(sim,5) # Should be NA
- # # Run the simulation
- # Simulate(sim)
- # # try again
- # getSeqFromNode(sim,5)
- # }
- #
- # @author
- #
- # \seealso{
- # @seeclass
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "getSeqFromNode",
- class="PhyloSim",
- function(
- this,
- node=NA,
- ...
- ){
- if(missing(node)){
- throw("No node specified!\n");
- }
- else if( .checkNode(this,node) ){
-
- return(this$.sequences[[as.numeric(node)]]);
- }
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: getSequences
- ##
- ###########################################################################/**
- #
- # @RdocMethod getSequences
- #
- # @title "Gets all the Sequence objects associated with the nodes of a phylo object aggregated by a PhyloSim object"
- #
- # \description{
- # @get "title".
- #
- # The order of the Sequence objects in the returned list reflects the identifiers of the associated nodes.
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{A PhyloSim object.}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # A list of sequence objects.
- # }
- #
- # \examples{
- # # Create a PhyloSim object.
- # # Provide the phylo object
- # # and the root sequence.
- # sim<-PhyloSim(
- # name="TinySim",
- # phylo=rcoal(3),
- # root.seq=NucleotideSequence(string="ATG",processes=list(list(JC69())))
- # );
- # # run the simulation
- # Simulate(sim)
- # # get all the associated sequence objects
- # getSequences(sim)
- # # get the sequence associated with node 3
- # # via virtual field
- # sim$sequences[[3]]
- # }
- #
- # @author
- #
- # \seealso{
- # @seeclass
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "getSequences",
- class="PhyloSim",
- function(
- this,
- ...
- ){
- slist<-list();
- for (node in getNodes(this)){
- slist[[node]]<-getSeqFromNode(this, node=node);
- }
- return(slist);
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: setSequences
- ##
- ###########################################################################/**
- #
- # @RdocMethod setSequences
- #
- # @title "Forbidden action: setting the Sequence objects associated with the nodes of a phylo object aggregated by a PhyloSim object"
- #
- # \description{
- # @get "title".
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{An object.}
- # \item{value}{Not used.}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # Throws an error.
- # }
- #
- # @author
- #
- # \seealso{
- # @seeclass
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "setSequences",
- class="PhyloSim",
- function(
- this,
- value,
- ...
- ){
- virtualAssignmentForbidden(this);
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: getAlignment
- ##
- ###########################################################################/**
- #
- # @RdocMethod getAlignment
- #
- # @title "Get the alignment stored in a PhyloSim object"
- #
- # \description{
- # @get "title".
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{A PhyloSim object.}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # The alignment as a matrix. Gap are represented by strings composed of dashes.
- # }
- #
- # \examples{
- # # Create a PhyloSim object.
- # # Provide the phylo object
- # # and the root sequence.
- # sim<-PhyloSim(
- # name="TinySim",
- # phylo=rcoal(3),
- # root.seq=NucleotideSequence(string="ATG",processes=list(list(JC69())))
- # );
- # # run the simulation
- # Simulate(sim)
- # # get the resulting aligment
- # getAlignment(sim)
- # # via virtual fileld:
- # sim$alignment
- # }
- #
- # @author
- #
- # \seealso{
- # @seeclass
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "getAlignment",
- class="PhyloSim",
- function(
- this,
- ...
- ){
- this$.alignment;
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: setAlignment
- ##
- ###########################################################################/**
- #
- # @RdocMethod setAlignment
- #
- # @title "Forbidden action: setting the alignment stored in a PhyloSim object"
- #
- # \description{
- # @get "title".
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{An object.}
- # \item{value}{Not used.}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # Throws an error.
- # }
- #
- # @author
- #
- # \seealso{
- # @seeclass
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "setAlignment",
- class="PhyloSim",
- function(
- this,
- value,
- ...
- ){
- this$alignment <- value;
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: .recoverAlignment
- ##
- setMethodS3(
- ".recoverAlignment",
- class="PhyloSim",
- function(
- this,
- paranoid=FALSE,
- ...
- ){
- # Refuse to build alignment if at least one of the sequences is NA:
- for (seq in this$.sequences){
- if(!is.Sequence(seq)){
- throw("Cannot build alignment because the simulation is incomplete!\n");
- }
- }
- # The list holding all the partial alignment matrices:
- aln.mat<-list();
-
- # Assigning NA-s here to prevent creation of these variables in the global
- # environment.
- row.names<-NA;
- from.node<-NA;
- to.node<-NA;
- from.seq<-NA;
- to.seq<-NA;
- edge<-NA;
- from.name<-NA;
- to.name<-NA;
- from.mat<-NA;
- to.mat<-NA;
- # Initialize the variables:
- init.vars<-function(){
- # Getting the edge:
- edge<<-getEdge(this, edge.number);
-
- # Getting the nodes:
- from.node<<-edge[[1,"from"]];
- to.node<<-edge[[1,"to"]];
- # Getting the sequence objects:
- from.seq<<-getSeqFromNode(this, from.node)
- to.seq<<-getSeqFromNode(this, to.node)
-
- # Getting sequence names:
- from.name<<-from.seq$name;
- to.name<<-to.seq$name;
- }
- # Initialize the aligment matrices:
- init.aln.mats<-function(){
- # Initialize "from" element in aln.mat if necessary:
- if( is.null(aln.mat[[from.name]] )){
- # Create a row of the states:
- tmp<-rbind(as.character(lapply(from.seq$.sites, getState)));
- # Label the columns by the site position:
- colnames(tmp)<-seq(along.with=from.seq$.sites);
- # Label the row with the sequence name:
- rownames(tmp)<-from.name;
- # Set the corresponding list element in aln.mat:
- aln.mat[[ from.name ]]<-tmp;
- }
- # Set from.mat
- from.mat<<-aln.mat[[ from.name ]];
-
- # Initialize "to" element int aln.mat if necessary
- if( is.null(aln.mat[[to.name]]) ){
- # Create a new entry if we are dealing with a tip:
- if(is.tip(this, to.node)){
- # Create a vector of states:
- tmp<-rbind(as.character(lapply(to.seq$.sites, getState)));
- # Label columns by position:
- colnames(tmp)<-seq(along.with=to.seq$.sites);
- # Label row by sequence name:
- rownames(tmp)<-to.name;
- aln.mat[[ to.name ]]<-tmp;
- }
- else {
- # A "to" element can be null only if its a tip:
- throw("aln.mat inconsistency!\n");
- }
- }
- # Set to.mat:
- to.mat<<-aln.mat[[ to.name ]];
- # Save row names:
- # The order is important! First "from", than "to"!
- row.names<<-c(rownames(from.mat), rownames(to.mat));
- }
- # Get the sequence position of a given alignment column from
- # the column labels:
- get.seq.pos<-function(mat=NA, col=NA){
- # Column number cannot be larger than length:
- if(col > dim(mat)[[2]]){
- throw("Invalid column number!\n");
- }
- # Get the corresponding column name:
- tmp<-colnames(mat)[[col]];
- # Return if NA:
- if(is.na(tmp)){
- return(NA);
- }
- else{
- return(as.numeric(tmp));
- }
- }
-
- # Check if two positions from the two *sequences* are homologous.
- is.homologous<-function(from.pos=NA, to.pos=NA){
- # Check position validity:
- if(to.pos > to.seq$length){
- throw("to.pos too big ",to.pos);
- }
- if(from.pos > from.seq$length){
- throw("from.pos too big ",from.pos);
-
- }
- # Check if the ancestral from to.seq/to.pos is from.seq/from.pos:
- return(equals(to.seq$.sites[[ to.pos ]]$.ancestral, from.seq$.sites[[ from.pos ]]));
- }
- # Get the symbol length from "to.seq" at position to.pos:
- get.to.symlen<-function(pos=NA){
- len<-stringLength(to.mat[to.name , pos]);
- if( is.na(len) | (len < 1) ){
- throw("Trouble in getting to.symlen!");
- } else {
- return(len);
- }
- }
- # Get the symbol length from "from.seq" at position from.pos:
- get.from.symlen<-function(pos=NA){
- len<-stringLength(from.mat[from.name , pos]);
- if( is.na(len) | (len < 1) ){
- throw("Trouble in getting from.symlen!");
- } else {
- return(len);
- }
- }
- make.gap.in.from<-function(label=NA,symlen=NA){
- # Create the gap symbol:
- gap<-paste(rep("-",times=symlen),collapse="");
- # Create the vector with gaps:
- gaps<-cbind(rep(gap, times=dim(from.mat)[[1]] ));
- # Label the column:
- colnames(gaps)<-c(label);
- # Bind the gaps with the corresponding column from to.mat,
- # and than bind with res.mat:
- res.mat<<-cbind(res.mat, rbind( gaps, cbind(to.mat[,j]) ) );
- # Restore rownames:
- rownames(res.mat)<<-row.names;
- # Increment counter for to.mat:
- j<<-j+1;
- }
- make.gap.in.to<-function(label=NA,symlen=NA){
- # See above.
- gap<-paste(rep("-",times=symlen),collapse="");
- gaps<-cbind(rep(gap, times=dim(to.mat)[[1]] ));
- colnames(gaps)<-c(label);
- res.mat<<-cbind(res.mat, rbind( cbind(from.mat[,i]), gaps ) );
- rownames(res.mat)<<-row.names;
- i<<-i+1;
- }
- emmit.homologous<-function(){
- # Bind the two columns into one column:
- tmp<-cbind(rbind( cbind(from.mat[,i]), cbind(to.mat[,j]) ) );
- # Label the column by from.pos:
- colnames(tmp)<-c(from.pos);
- # Set res.mat
- res.mat<<-cbind(res.mat, tmp );
- # resotre rownames:
- rownames(res.mat)<<-row.names;
-
- i<<-i+1;
- j<<-j+1;
- }
- # Iterate over the reverse of the edge matrix:
- for (edge.number in rev(seq(from=1, to=this$nedges))){
- # Call variable initialization:
- init.vars();
- # Initialize partial alignment matrices:
- init.aln.mats();
-
- # The matrix holding the resulting partial alignment:
- res.mat<-c();
- # Column counter for from.mat:
- i<-1;
- # Column counter for to.mat:
- j<-1;
-
- while(i <=dim(from.mat)[[2]] | j <=dim(to.mat)[[2]] ){
- # First of all, check for counter overflow:
- if(i > dim(from.mat)[[2]]){
- # If i is greater than the length of from.mat, but we
- # are still iterating, that means that we have parts left from
- # to.mat, so we have to create gaps in from.mat, and increment j.
- make.gap.in.from(symlen=get.to.symlen(pos=j));
- next();
- }
- else if (j > dim(to.mat)[[2]]){
- # If j is greater than the length of to.mat and we still iterating,
- # that means that we have still some columns in from.mat, so we create
- # gaps in to.mat and increment i. We label the new column with from.pos.
- from.pos<-get.seq.pos(mat=from.mat, col=i);
- make.gap.in.to(label=from.pos,get.from.symlen(pos=i));
- next();
- }
- # Now figure out the positions:
- from.pos<-get.seq.pos(mat=from.mat, col=i);
- to.pos<-get.seq.pos(mat=to.mat, col=j);
- # Now check for the gaps wich have been introduced before:
- if(is.na(from.pos)){
- # If we have a gap in from.mat,
- # than emmit the columnt with a gap in "to":
- make.gap.in.to(symlen=get.from.symlen(pos=i));
- next();
- }
- if(is.na(to.pos)){
- # Existent gap in to.mat:
- make.gap.in.from(symlen=get.to.symlen(pos=j));
- next();
- }
- # Now we have some real alignment to do here:
- if(is.homologous(from.pos=from.pos, to.pos=to.pos)){
- # We have to homologous columns, bind them, and emmit:
- emmit.homologous();
- next();
- }
- else if(is.Process(to.seq$.sites[[to.pos]]$.ancestral)){
- # The two columns are not homologous. The column in "to"
- # was inserted by a process. Make gap in from:
- make.gap.in.from(symlen=get.from.symlen(pos=i));
- next();
- }
- else {
- # The only possibility left is a deletion in the child sequence.
- # Make gaps in "to", label new column by from.pos:
- make.gap.in.to(label=from.pos,symlen=get.from.symlen(pos=i));
- next();
- }
- } # while i | j
-
- # Replace the "from" element in aln.mat with the resulting partial
- # alignment matrix:
- aln.mat [[ from.name ]]<-res.mat;
- } # for edge.number
- alignment <-aln.mat[[ this$rootSeq$name ]];
- # Check the correcteness of the alignment if paranoid:
- if(paranoid){
- .checkAlignmentConsistency(this, alignment);
- }
- # Call garbage collection:
- gc();
- gc();
- # The whole alignment is associated with the root node:
- return(alignment);
-
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: .checkAlignmentConsistency
- ##
- setMethodS3(
- ".checkAlignmentConsistency",
- class="PhyloSim",
- function(
- this,
- aln,
- ...
- ){
-
- # First check if the sequences are intact:
- for(node in this$nodes){
- seq.node<-getSeqFromNode(this, node);
- seq.aln<-aln[seq.node$name, ];
- seq.aln<-seq.aln[ grep("^[^-]+$",seq.aln) ];
- seq.aln<-paste(seq.aln, collapse="");
- if(seq.aln != seq.node$string){
- throw("The alignment is inconsistent with the sequence objects!\n Blaming ",seq.node$name, ".\n");
- }
- }
- for(edge.number in rev(seq(from=1, to=this$nedges))){
- # Getting the edge:
- edge<-getEdge(this, edge.number);
-
- # Getting the nodes:
- from.node<-edge[[1,"from"]];
- to.node<-edge[[1,"to"]];
- # Getting the sequence objects:
- from.seq<-getSeqFromNode(this, from.node)
- to.seq<-getSeqFromNode(this, to.node)
-
- # Getting sequence names:
- from.name<-from.seq$name;
- to.name<-to.seq$name;
- # Initializing positions:
- from.pos<-1;
- to.pos<-1;
- is.gap<-function(string){
-
- res<-length(grep("^-+$",string));
- if( res == 1 ){
- return(TRUE);
- }
- else if (res > 1){
- throw("is.gap: argument vector too long!\n");
- }
- else {
- return(FALSE);
- }
-
- }
-
- # Iterate over edges:
- for (i in 1:dim(aln)[[2]]){
- # Overflow in "from" counter,
- if(from.pos > from.seq$length){
- to.char<-aln[to.name,i];
- if(!is.gap(to.char)){
- # we have a final insertion:
- if(!is.Process(to.seq$.sites[[to.pos]]$.ancestral)){
- throw("Alignment insertion inconsistency!\n");
- }
- to.pos<-to.pos+1;
- }
- next();
- }
-
- # Overflow in "to" counter (final deletion):
- if(to.pos > to.seq$length){
- break();
- }
-
- # Get the symbols from alignment:
- from.char<-aln[from.name,i];
- to.char<-aln[to.name,i];
- is.gap.to<-is.gap(to.char);
- is.gap.from<-is.gap(from.char);
-
- # Skip if we have to gap symbols:
- if( is.gap.from & is.gap.to ){
- next();
- }
- # Deletion in to.seq:
- else if(is.gap.to & !is.gap.from ){
- from.pos<-(from.pos+1);
- }
- # Insertion in to.seq:
- else if(!is.gap.to & is.gap.from ){
- # Check ancestral pointer for inserted sites:
- if(!is.Process(to.seq$.sites[[to.pos]]$.ancestral)){
- throw("Alignment insertion inconsistency!\n");
- }
- to.pos<-(to.pos+1);
- } else {
- # We must have a homology here:
- if(!equals(to.seq$.sites[[ to.pos ]]$.ancestral, from.seq$.sites[[ from.pos ]])){
- throw("Non-homologous sites aligned! Alignment is inconsistent!\n");
- }
- from.pos<-(from.pos+1);
- to.pos<-(to.pos+1);
- }
-
- }
- }
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: saveAlignment
- ##
- ###########################################################################/**
- #
- # @RdocMethod saveAlignment
- #
- # @title "Save the alignment stored in a PhyloSim object in a Fasta file"
- #
- # \description{
- # @get "title".
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{A PhyloSim object.}
- # \item{file}{The name of the output file.}
- # \item{skip.internal}{Do not save sequences corresponding to internal nodes.}
- # \item{paranoid}{Check the consistency of the alignment.}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # The PhyloSim object (invisible).
- # }
- #
- # \examples{
- # # Create a PhyloSim object.
- # # Provide the phylo object
- # # and the root sequence.
- # sim<-PhyloSim(
- # name="TinySim",
- # phylo=rcoal(3),
- # root.seq=NucleotideSequence(string="ATG",processes=list(list(JC69())))
- # );
- # # run the simulation
- # Simulate(sim)
- # # save the alignment
- # file<-paste("PhyloSim_dummy_fasta_",Sys.getpid(),".fas",sep="");
- # saveAlignment(sim,file=file,paranoid=TRUE);
- # # print out the Fasta file
- # cat(paste(scan(file=file,what=character(),sep="\n"),collapse="\n"));cat("\n");
- # # delete Fasta file
- # unlink(file);
- # }
- #
- # @author
- #
- # \seealso{
- # @seeclass
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "saveAlignment",
- class="PhyloSim",
- function(
- this,
- file="phylosim.fas",
- skip.internal=FALSE,
- paranoid=FALSE,
- ...
- ){
- if(any(is.na(this$.alignment))){
- warning("Alignment is undefined, nothin to save!\n");
- return();
- }
- else {
- if(paranoid){
- .checkAlignmentConsistency(this, this$.alignment);
- }
- sink(file);
- if(!skip.internal){
- for(i in 1:dim(this$.alignment)[[1]]){
- cat(">",rownames(this$.alignment)[[i]],"\n");
- cat(paste(this$.alignment[i,],collapse=""),"\n");
- }
- } else {
- for(i in 1:dim(this$.alignment)[[1]]){
- name<-rownames(this$.alignment)[[i]];
- if(!any((length(grep("^Node \\d+$",name,perl=TRUE,value=FALSE)) > 0),(length(grep("^Root node \\d+$",name,perl=TRUE,value=FALSE)) > 0))){
- cat(">",name,"\n");
- cat(paste(this$.alignment[i,],collapse=""),"\n");
- }
- }
- }
- sink(NULL);
- }
- return(invisible(this));
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: plot.PhyloSim
- ##
- ###########################################################################/**
- #
- # @RdocMethod plot
- #
- # @title "Plot a PhyloSim object"
- #
- # \description{
- # @get "title".
- #
- # This method plots the aggregated alignment alongside the tree used for simulation. Various options
- # allow for control over the plot style.
- #
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{x}{A PhyloSim object.}
- # \item{plot.tree}{Whether to plot the tree alongside the alignment. TRUE or FALSE; defaults to TRUE.}
- # \item{plot.ancestors}{Whether to plot the ancestral sequences. TRUE or FALSE; defaults to TRUE.}
- # \item{plot.chars}{Whether to plot the actual text of the characters.}
- # \item{plot.legend}{Whether to plot the legend showing the character-to-color mapping.}
- # \item{plot.labels}{Whether to plot the sequence labels along the y-axis}
- # \item{aspect.ratio}{(Experimental; when set, this option forces the num.pages value to 1) Constrains the alignment residues to have a certain aspect ratio; values above 1 cause vertically-stretched blocks. FALSE disables aspect ratio control, numerical values set the aspect ratio; defaults to FALSE.}
- # \item{num.pages}{Optionally split the alignment over a number of vertically-stacked pages. This is useful for long alignments. 'auto' chooses a sensible number of pages, numerical values specify a number; defaults to 'auto'.}
- # \item{char.text.size}{Text size for the aligned characters. This may require tweaking depending on the DPI and output format. Defaults to 'auto'.}
- # \item{axis.text.size}{Text size for the sequence labels along the y-axis. This may require tweaking depending on the DPI and output format. Defaults to 'auto'.}
- # \item{color.scheme}{Color scheme to use ("auto", "binary", "dna", "protein", "codon", "combined", "combined_codon"). Defaults to 'auto'. When set to 'auto', the function will choose an appropriate coloring scheme based on the alignment content.}
- # \item{color.branches}{The event count used to color the branches ("substitutions" by default). See \code{\link{getBranchEvents.PhyloSim}}.}
- # \item{tree.xlim}{The x-axis limits of the tree panel.}
- # \item{aln.xlim}{The x-axis limits of the alignment panel (in alignment column coordinates).}
- # \item{tracks}{Tracks to display above or below the alignment as colored blocks.
- #
- # The input format for tracks is a list of data frames with the following possible fields, all of which are optional and can be omitted:
- # \itemize{
- # \item pos - the sequence position (starting with 1) of the feature. Defaults to NULL.
- # \item score - the score (between 0 and 1) of the feature. Scores above 1 or below zero will be truncated.
- # Defaults to 1.
- # \item y_lo - the lower Y offset (between 0 and 1) of the feature. Defaults to 0. Use a y_lo and y_hi
- # value for each row in the track data frame to create a wiggle plot like effect.
- # \item y_hi - the upper Y offset (between 0 and 1) of the feature. Defaults to 1. Use just a y_hi value
- # for each row in the track data frame to create a bar plot like effect.
- # \item [the fields below are considered unique per track; the values from the first row in the track
- # data frame are used.]
- # \item id - the display ID for the track. Defaults to 'Track'.
- # \item layout - set to 'above' to put the track above the alignment, 'below' for below.
- # \item height - the number of alignment rows for the track to span in height. Defaults to 3.
- # \item color.gradient - a comma-separated list of colors to interpolate between when coloring
- # the blocks. Examples: 'white,red' 'blue,gray,red' '#FF00FF,#FFFFFF'. Defaults to 'white,black'.
- # \item color - a single color to use when coloring the blocks. Mutually exclusive with color.gradient,
- # and if a color.gradient value exists then this value will be ignored. Defaults to black.
- # \item background - a color for the background of the track. Defaults to white.
- # }}
- # \item{aln.length.tolerance}{The desired alignment/sequence length ratio (A/S ratio) to trim the alignment to.
- # The A/S ratio is defined as the ratio between the alignment length and the mean ungapped sequence length, and
- # the alignment trimming procedure will remove blocks of indel-containing columns (in a sensible order) until
- # either (a) the desired indel tolerance is reached, or (b) no more columns can be removed without yielding an empty
- # alignment. A track is added below the alignment to indicate how many indels each resulting alignment column used
- # used to harbor, and black squares are overlaid onto the alignment where extant sequence data has been trimmed.
- # Defaults to NULL (no trimming); values in the range of 0.9 to 1.3 tend to work well at improving the
- # legibility of very gappy alignments.}
- # \item{plot.nongap.bl}{If set to TRUE, plots the non-gap branch length (defined as the branch length of the subtree of non-gapped sequences) as a track below the alignment. Defaults to FALSE.}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # The PhyloSim object (invisible).
- # }
- #
- # \examples{
- # # Create a PhyloSim object.
- # # Provide the phylo object
- # # and the root sequence.
- # sim<-PhyloSim(
- # name="TinySim",
- # phylo=rcoal(3),
- # root.seq=NucleotideSequence(string="ATGCTAGCTAGG",processes=list(list(JC69())))
- # );
- # # plot the aggregated phylo object
- # plot(sim)
- # # run simulation
- # Simulate(sim)
- # # Plot the alignment without the tree or ancestral sequences.
- # plot(sim, plot.ancestors=FALSE, plot.tree=FALSE)
- # # Force a DNA-based color scheme
- # # (default is 'auto' to auto-detect based on the sequence composition)
- # plot(sim, color.scheme='dna', plot.legend=TRUE)
- # }
- #
- # @author
- #
- # \seealso{
- # @seeclass
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "plot",
- class="PhyloSim",
- function(
- x,
- plot.tree,
- plot.ancestors,
- plot.chars,
- plot.legend,
- plot.labels,
- aspect.ratio,
- num.pages,
- char.text.size,
- axis.text.size,
- color.scheme,
- color.branches,
- tree.xlim,
- aln.xlim,
- tracks,
- aln.length.tolerance,
- plot.nongap.bl,
- ...
- ){
-
- if(missing(char.text.size)){
- char.text.size <- 'auto'
- }
- if(missing(axis.text.size)){
- axis.text.size <- 'auto'
- }
- if(missing(color.scheme)){
- color.scheme <- 'auto'
- }
- if(missing(color.branches)){
- color.branches <- 'substitutions'
- }
- if(missing(plot.tree)){
- plot.tree <- TRUE
- }
- if(missing(plot.ancestors)){
- plot.ancestors <- TRUE
- }
- if(missing(plot.chars)){
- plot.chars <- TRUE
- }
- if(missing(plot.legend)){
- plot.legend <- FALSE
- }
- if(missing(plot.labels)){
- plot.labels <- TRUE
- }
- if(missing(aspect.ratio)){
- aspect.ratio <- FALSE
- }
- if(missing(num.pages)){
- num.pages='auto'
- }
- if(missing(color.scheme)){
- color.scheme='auto'
- }
- if(missing(tree.xlim)){
- tree.xlim=NULL
- }
- if(missing(aln.xlim)){
- aln.xlim=NULL
- }
- if(missing(tracks)){
- tracks=NULL
- }
- if(any(is.na(x$.phylo))) {
- plot.tree <- FALSE
- }
- if(missing(aln.length.tolerance)){
- aln.length.tolerance=NULL
- }
- if(missing(plot.nongap.bl)){
- plot.nongap.bl=FALSE
- }
-
- if(all(!is.na(x$.alignment), is.matrix(x$.alignment))){
- .plotWithAlignment(x,
- plot.tree=plot.tree,
- plot.ancestors=plot.ancestors,
- plot.chars=plot.chars,
- plot.legend=plot.legend,
- plot.labels=plot.labels,
- aspect.ratio=aspect.ratio,
- num.pages=num.pages,
- char.text.size=char.text.size,
- axis.text.size=axis.text.size,
- color.scheme=color.scheme,
- color.branches=color.branches,
- tree.xlim=tree.xlim,
- aln.xlim=aln.xlim,
- tracks=tracks,
- aln.length.tolerance=aln.length.tolerance,
- plot.nongap.bl=plot.nongap.bl
- );
- return(invisible(x));
- }
- plot(x$.phylo);
- nodelabels();
- return(invisible(x));
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: .plotWithAlignment
- ##
- setMethodS3(
- ".plotWithAlignment",
- class="PhyloSim",
- function(
- x,
- plot.tree,
- plot.ancestors,
- plot.chars,
- plot.legend,
- plot.labels,
- aspect.ratio,
- num.pages,
- char.text.size,
- axis.text.size,
- color.scheme,
- color.branches,
- tree.xlim,
- aln.xlim,
- tracks,
- aln.length.tolerance,
- plot.nongap.bl,
- ...
- ){
- # ugly empirical fix of some R CMD check warnings:
- id <-NA;
- pos <-NA;
- char <-NA;
- xend <-NA;
- yend <-NA;
- y <-NA;
- substitutions<-NA;
- event.count <-NA;
- type <-NA;
- track_index <-NA;
- xx <-NA;
- yy <-NA;
- xmin <-NA;
- xmax <-NA;
- ymin <-NA;
- ymax <-NA;
- ### First, we need to define a bunch of sparsely-documented utility functions. ###
- # Re-orders the alignment rows by matching the tree's tips.
- sort.aln.by.tree = function(aln,tree) {
- names <- dimnames(aln)[[1]]
- if (length(names) > length(tree$tip.label)) {
- return(aln)
- }
- newPositions <- rev(match(tree$tip.label,names))
- aln <- aln[newPositions,]
- dimnames(aln) <- list(names[newPositions])
- return(aln)
- }
-
- # Extracts a list of child node IDs for the given node. Returns (-1,-1) if the node is a leaf.
- child.nodes <- function(phylo,node) {
- edge.indices <- which(phylo$edge[,1]==node)
- nodes <- phylo$edge[edge.indices,2]
- if (length(nodes)==0) {
- nodes <- list(c(-1,-1))
- } else {
- nodes <- list(nodes)
- }
- return(list(nodes))
- }
- # Extracts the parent node ID for the given node. Returns -1 if the node is root.
- parent.node <- function(phylo,node) {
- edge.index <- which(phylo$edge[,2]==node)
- node <- phylo$edge[edge.index,1]
- if (length(node)==0) {
- node <- -1
- }
- return(node)
- }
- generic.count <- function(sim,node,type) {
- if (type == 'insertions' || type == 'ins') {
- type <- 'insertion'
- }
- if (type == 'deletions' || type == 'del') {
- type <- 'deletion'
- }
- if (type == 'substitutions' || type == 'subst' || type == 'sub') {
- type <- 'substitution'
- }
- if (type == 'syn') {
- type <- 'synonymous'
- }
- if (type == 'nsyn') {
- type <- 'non-synonymous'
- }
- # Default value of 0.
- cnt <- 0
- if(type=='none' || type == '') {
- return(as.numeric(cnt))
- }
- # Find the edge which points to the given node.
- edge.index <- which(phylo$edge[,2]==node)
- if (length(edge.index) > 0) {
- bs <- sim$.branch.stats
- cnt <- bs[[paste(edge.index)]][[type]]
- if (is.null(cnt) || is.na(cnt)) {
- cnt <- 0
- }
- }
- return(as.numeric(cnt))
- }
-
- subst.count <- function(sim,node) {
- return(generic.count(sim,node,'substitution'))
- }
- ins.count <- function(sim,node) {
- return(generic.count(sim,node,'insertion'))
- }
- del.count <- function(sim,node) {
- return(generic.count(sim,node,'deletion'))
- }
- syn.count <- function(sim,node) {
- return(generic.count(sim,node,'synonymous'))
- }
- nsyn.count <- function(sim,node) {
- return(generic.count(sim,node,'non-synonymous'))
- }
- # Finds the node with a given label.
- node.with.label <- function(tree,label) {
- return(which(tree$tip.label %in% label))
- }
-
- # Extracts the length of the branch above the given node. Returns 0 if the node is root.
- branch.length <- function(phylo,node) {
- edge.index <- which(phylo$edge[,2]==node)
- bl <- phylo$edge.length[edge.index]
- if (length(bl)==0) {
- bl <- 0
- }
- return(bl)
- }
- # The maximum root-to-tip length in the tree.
- max.length.to.root <- function(phylo) {
- max.length <- 0
- for (i in 1:length(phylo$tip.label)) {
- cur.length <- length.to.root(phylo,i)
- max.length <- max(max.length,cur.length)
- }
- return(max.length)
- }
- # The length from the root to the given node. Can be given either as a node ID or a tip label.
- length.to.root <- function(phylo,node) {
- tip.index <- node
- if (is.character(node)) {
- tip.index <- which(phylo$tip.label==node)
- }
-
- cur.node.b <- tip.index
- p.edges <- phylo$edge
- p.lengths <- phylo$edge.length
-
- length <- 0
- while(length(which(p.edges[,2]==cur.node.b)) > 0) {
- cur.edge.index <- which(p.edges[,2]==cur.node.b)
- cur.edge.length <- p.lengths[cur.edge.index]
- length <- length + cur.edge.length
- cur.node.a <- p.edges[cur.edge.index,1]
- cur.node.b <- cur.node.a # Move up to the next edge
- }
- return(length)
- }
- # Returns a data frame defining segments to draw the phylogenetic tree.
- phylo.layout.df <- function(phylo,layout.ancestors=FALSE,align.seq.names=NULL) {
- # Number of nodes and leaves.
- n.nodes <- length(phylo$tip.label)+phylo$Nnode
- n.leaves <- length(phylo$tip.label)
- # Create the skeleton data frame.
- df <- data.frame(
- node=c(1:n.nodes), # Nodes with IDs 1 to N.
- x=0, # These will contain the x and y coordinates after the layout procedure below.
- y=0,
- label=c(phylo$tip.label,((n.leaves+1):n.nodes)), # The first n.leaves nodes are the labeled tips.
- is.leaf=c(rep(TRUE,n.leaves),rep(FALSE,n.nodes-n.leaves)), # Just for convenience, store a boolean whether it's a leaf or not.
- parent=0, # Will contain the ID of the current node's parent
- children=0, # Will contain a list of IDs of the current node's children
- branch.length=0 # Will contain the branch lengths
- )
- # Collect the parents, children, and branch lengths for each node
- parent <- c()
- bl <- list()
- children <- list()
- event.count <- list()
- for (i in 1:nrow(df)) {
- node <- df[i,]$node
- parent <- c(parent,parent.node(phylo,node))
- bl <- c(bl,branch.length(phylo,node))
- children <- c(children,child.nodes(phylo,node))
- event.count <- c(event.count,generic.count(x,node,color.branches))
- }
- df$parent <- parent
- df$children <- children
- df$branch.length <- bl
- df$event.count <- as.numeric(event.count)
- # Start the layout procedure by equally spacing the leaves in the y-dimension.
- df[df$is.leaf==TRUE,]$y = c(1:n.leaves)
- found.any.internal.node.sequences <- FALSE
- # For each leaf: travel up towards the root, laying out each internal node along the way.
- for (i in 1:n.leaves) {
- cur.node <- i
- while (length(cur.node) > 0 && cur.node != -1) {
- # We always use branch lengths: x-position is simply the length to the root.
- df[cur.node,]$x <- length.to.root(phylo,cur.node)
- # The y-position for internal nodes is the mean of the y-position of the two children.
- children <- unlist(df[cur.node,]$children)
- if (length(children) > 0 && children[1] != -1) {
- child.a <- children[1]
- child.b <- children[2]
- child.a.y <- df[child.a,]$y
- child.b.y <- df[child.b,]$y
- df[cur.node,]$y <- (child.a.y+child.b.y)/2
- }
- # Try to find the index of this node in the alignment names.
- if (!is.null(align.seq.names)) {
- lbl <- df[cur.node,]$label
- index.in.names <- which(align.seq.names == lbl | align.seq.names %in% c(paste('Node',lbl),paste('Root node',lbl)))
- if (length(index.in.names)>0) {
- df[cur.node,]$y <- index.in.names
- if (!df[cur.node,]$is.leaf) {
- found.any.internal.node.sequences <- TRUE
- }
- }
- }
-
- cur.node <- unlist(df[cur.node,]$parent)
- }
- }
-
- # We have a data frame with each node positioned.
- # Now we go through and make two line segments for each node (for a 'square corner' type tree plot).
- line.df <- data.frame()
- for (i in 1:nrow(df)) {
- row <- df[i,] # Data frame row for the current node.
- if (row$parent == -1) {
- next; # Root node!
- }
- p.row <- df[row$parent,] # Data frame row for the parent node.
- if (layout.ancestors && found.any.internal.node.sequences) {
- horiz.line <- data.frame(
- x=row$x,
- xend=p.row$x,
- y=row$y,
- yend=p.row$y,
- lbl=row$label,
- event.count=row$event.count
- )
- line.df <- rbind(line.df,horiz.line)
- } else {
- horiz.line <- data.frame(
- x=row$x,
- xend=p.row$x,
- y=row$y,
- yend=row$y,
- lbl=row$label,
- event.count=row$event.count
- ) # First a line from row.x to parent.
- vert.line <- data.frame(
- x=p.row$x,
- xend=p.row$x,
- y=row$y,
- yend=p.row$y,
- lbl=row$label,
- event.count=row$event.count
- ) # Now a line from row.y to parent.y
-
- #horiz.line <- data.frame(x=row$x,xend=(p.row$x+row$x)/2,y=row$y,yend=row$y,lbl=row$label) # First a line from row.x to parent.
- #vert.line <- data.frame(x=(p.row$x+row$x)/2,xend=p.row$x,y=row$y,yend=p.row$y,lbl=row$label) # Now a line from row.y to parent.y
- line.df <- rbind(line.df,horiz.line,vert.line)
- }
- }
- return(line.df)
- }
- # Call this to put a ggplot panel into a specified layout position [for example: print(p,vp=subplot(1,2)) ]
- subplot <- function(x, y) viewport(layout.pos.col=x, layout.pos.row=y)
- # Call this to create a layout with x and y rows and columns, respectively
- vplayout <- function(x, y) {
- grid.newpage()
- pushViewport(viewport(layout=grid.layout(y,x)))
- }
- # Creates a color aesthetic for alignments
- alignment.colors <- function(scheme,darken=F) {
- scheme <- tolower(scheme)
- if (scheme == 'binary') {
- cols <- c(
- '0' = "#000000",
- '1' = '#FFFFFF'
- )
- } else if (scheme == 'dna') {
- cols <- c(
- 'G' = "#FFFF00",
- 'C' = "#00FF00",
- 'T' = "#FF0000",
- 'A' = "#0000FF"
- )
- } else if (scheme == 'numeric') {
- #cols <- heat.colors(10)
- cols <- colorRampPalette(c("red","yellow","green"))(10)
- cols <- c(
- '0' = cols[1],
- '1' = cols[2],
- '2' = cols[3],
- '3' = cols[4],
- '4' = cols[5],
- '5' = cols[6],
- '6' = cols[7],
- '7' = cols[8],
- '8' = cols[9],
- '9' = cols[10]
- )
- } else if (scheme == 'taylor' || scheme == 'protein') {
- cols <- c(
- 'A' = "#CCFF00", 'a' = "#CCFF00",
- 'C' = "#FFFF00", 'c' = "#FFFF00",
- 'D' = "#FF0000", 'd' = "#FF0000",
- 'E' = "#FF0066", 'e' = "#FF0066",
- 'F' = "#00FF66", 'f' = "#00FF66",
- 'G' = "#FF9900", 'g' = "#FF9900",
- 'H' = "#0066FF", 'h' = "#0066FF",
- 'I' = "#66FF00", 'i' = "#66FF00",
- 'K' = "#6600FF", 'k' = "#6600FF",
- 'L' = "#33FF00", 'l' = "#33FF00",
- 'M' = "#00FF00", 'm' = "#00FF00",
- 'N' = "#CC00FF", 'n' = "#CC00FF",
- 'P' = "#FFCC00", 'p' = "#FFCC00",
- 'Q' = "#FF00CC", 'q' = "#FF00CC",
- 'R' = "#0000FF", 'r' = "#0000FF",
- 'S' = "#FF3300", 's' = "#FF3300",
- 'T' = "#FF6600", 't' = "#FF6600",
- 'V' = "#99FF00", 'v' = "#99FF00",
- 'W' = "#00CCFF", 'w' = "#00CCFF",
- 'Y' = "#00FFCC", 'y' = "#00FFCC",
- '2' = "#888888", '2' = "#888888",
- 'O' = "#424242", 'o' = "#424242",
- 'B' = "#7D7D7D", 'b' = "#7D7D7D",
- 'Z' = "#EEEEEE", 'z' = "#EEEEEE",
- 'X' = "#000000", 'x' = "#000000"
- )
- } else if (scheme == 'codon') {
- # Get the protein colors.
- protein.colors <- alignment.colors('protein')
- # Create a list of codons.
- ca <- CodonAlphabet()
- nucs <- c('G','A','C','T')
- codons <- expand.grid(nucs,nucs,nucs,stringsAsFactors=F)
- # For each codon give the protein color.
- codon.colors <- c()
- for (i in 1:nrow(codons)) {
- codon = paste(codons[i,],collapse='')
- aa <- translateCodon(ca,codon)
- codon.colors[codon] = protein.colors[aa]
- }
- cols <- codon.colors
- } else if (scheme == 'combined') {
- dna.colors <- alignment.colors('dna')
- binary.colors <- alignment.colors('binary')
- protein.colors <- alignment.colors('protein')
- cols <- c(dna.colors,protein.colors,binary.colors)
- } else if (scheme == 'combined_codon') {
- dna.colors <- alignment.colors('dna')
- # Make the DNA stand out here by being a little darker
- for (i in 1:length(dna.colors)) {
- color <- dna.colors[i]
- darker.color <- darker(color)
- dna.colors[i] <- darker.color
- }
- binary.colors <- alignment.colors('binary')
- protein.colors <- alignment.colors('protein')
- codon.colors <- alignment.colors('codon')
- # Put them all together. (One remaining issue: the protein G,A,T,C will be colored as DNA!)
- cols <- c(dna.colors,protein.colors,binary.colors,codon.colors)
- }
- if (darken) {
- for (i in 1:length(cols)) {
- color <- cols[i]
- darker.color <- darker(color,0.85)
- cols[i] <- darker.color
- }
- }
-
- return(cols)
- }
- darker <- function(color,factor=0.7) {
- x <- col2rgb(color)
- x <- round(x * factor)
- y <- rgb(x[1],x[2],x[3],maxColorValue=255)
- return(y)
- }
- lighter <- function(color) {
- x <- col2rgb(color)
- x <- round(x * 1.2)
- y <- rgb(min(x[1],255),min(x[2],255),min(x[3],255),maxColorValue=255)
- return(y)
- }
- # Scores each column according to the branch length of the subtree created by
- # non-gap residues (we use nongap.bl as the var name), and stores the 'nongap.str'
- # which is the pasted list of labels for non-gapped sequences at this site.
- # This information is used by the 'remove.gaps' function to remove columns to
- # reach a certain alignment length threshold.
- score.aln.columns <- function(tree,aln) {
- aln.length <- length(aln[1,])
- score.df <- data.frame()
- for (i in 1:aln.length) {
- aln.column <- aln[,i]
- nongap.seqs <- names(aln.column[aln.column != '-'])
- gap.seqs <- names(aln.column[aln.column == '-'])
- # Get the non-gap branch length.
- if (length(nongap.seqs) == 1) {
- nongap.node <- node.with.label(tree,nongap.seqs[1])
- nongap.bl <- branch.length(tree,nongap.node)
- } else {
- nongap.tree <- drop.tip(tree,gap.seqs)
- nongap.bl <- sum(nongap.tree$edge.length)
- }
- nongap.str <- paste(nongap.seqs,collapse=';')
- cur.df <- data.frame(
- pos=i,
- score=nongap.bl,
- nongap.str=nongap.str,
- stringsAsFactors=F
- )
- score.df <- rbind(score.df,cur.df)
- }
- score.df <- score.df[order(score.df$score,score.df$pos),]
- return(score.df)
- }
- # Goes through the alignment of the sim object and removes columns until either the desired
- # tolerance is reached or there are no more low-scoring columns to remove (whichever comes
- # first). Tolerance is defined as the ratio of the alignment length to the mean sequence
- # length. So, an alignment with no indels at all is exactly 1; an alignment with lots
- # of deletions (but no insertion) is less than 1; an alignment with lots of insertion is
- # greater than 1.
- #
- # Values of 0.9 - 1.3 tend to give good results in a variety of situations.
- #
- remove.gaps <- function(sim,tolerance) {
- aln <- sim$.alignment
- tree <- sim$.phylo
- col.scores <- score.aln.columns(tree,aln)
- # Store the deletion markers in a separate data frame.
- deletion.df <- NULL
- # Get the mean sequence length
- lengths <- apply(aln,1,function(x) {
- x <- x[x != '-']
- return(stringLength(paste(x,collapse='')))
- })
- mean.seq.length <- mean(lengths)
- repeat {
- aln.length <- length(aln[1,])
- ratio <- aln.length / mean.seq.length
- if (ratio < tolerance) {
- break;
- }
- # Take the next site from the sorted scores
- lowest.scores <- col.scores[1,]
- col.scores <- col.scores[-c(1),]
- cur.pos <- lowest.scores$pos # Current position of lowest-scoring column.
- cur.score <- lowest.scores$score # The current column score.
- cur.str <- lowest.scores$nongap.str # The current column's nongap pattern.
- # Grab the entire 'current chunk' of alignment which has the same
- # score and non-gap string.
- repeat {
- first.pos <- col.scores[1,'pos']
- first.score <- col.scores[1,'score']
- first.str <- col.scores[1,'nongap.str']
- if (cur.score == max(col.scores$score)) {
- lowest.scores <- NULL
- cur.ratio <- length(aln[1,]) / mean.seq.length
- print(sprintf("Nothing left to remove at ratio %.2f!",cur.ratio))
- break;
- }
- if (first.pos == cur.pos + 1 && first.score == cur.score && first.str == cur.str) {
- cur.pos <- col.scores[1,]$pos
- lowest.scores <- rbind(lowest.scores,col.scores[1,])
- col.scores <- col.scores[-1,]
- } else {
- #print("Done!")
- break;
- }
- }
- if (is.null(lowest.scores)) {
- break;
- }
- # remove.us should be a contiguous vector of integers,
- # representing the set of columns to remove.
- remove.us <- lowest.scores$pos
- if (any(diff(remove.us) > 1)) {
- print("ERROR: Removing a non-consecutive set of columns!")
- }
- #print(paste("Removing at ",paste(remove.us[1])))
- # Go through columns from right to left, making sure to update
- # the new positions of columns on the right side of the splice.
- rev.pos <- rev(remove.us)
- for (i in 1:length(rev.pos)) {
- cur.pos <- rev.pos[i]
- aln <- splice.column(aln, cur.pos)
- # Update new positions of column scores
- above <- which(col.scores$pos > cur.pos)
- col.scores[above,]$pos <- col.scores[above,]$pos - 1
- # Update new positions of deletion locations.
- if (!is.null(deletion.df) && nrow(deletion.df) > 0) {
- above <- which(deletion.df$pos > cur.pos)
- deletion.df[above,]$pos <- deletion.df[above,]$pos - 1
- }
- }
- # Add a single entry to the data frame of deletions, to be used
- # by the plot function to indicate deletion points.
- cur.deletion <- data.frame(
- pos = cur.pos, # The first position AFTER the deletion splice.
- length = length(rev.pos), # The length of the block removed.
- nongap.str = as.character(cur.str), # nongap sequence IDs for this deletion
- stringsAsFactors=F
- )
- #print(cur.deletion)
- deletion.df <- rbind(deletion.df,cur.deletion)
- }
- # Create a new PhyloSim object, assign the tree & aln, and return.
- sim.temp <- PhyloSim();
- sim.temp$.alignment <- aln
- sim.temp$.phylo <- tree
- sim.temp$.indels <- deletion.df
- return(sim.temp)
- }
- splice.column <- function(aln,pos) {
- return(aln[,-pos])
- }
-
- ####################################
- ### Let the real plotting begin! ###
- ####################################
-
- # Apply the deletion tolerance if needed.
- if (!is.null(aln.length.tolerance)) {
- x <- remove.gaps(x, tolerance=aln.length.tolerance)
- indels <- x$.indels
- } else {
- indels <- NULL
- }
- df <- data.frame()
- aln <- x$.alignment
- phylo <- x$.phylo
- # Do some reordering of alignment & tree.
- if (!any(is.na(phylo))) {
- x$.phylo <- reorder(x$.phylo, order="cladewise");
- phylo <- x$.phylo
- aln <- sort.aln.by.tree(aln,phylo)
- }
- names <- dimnames(aln)[[1]]
-
- #print(paste("Aln length:",length(aln[1,])))
- #print(paste("Num seqs:",length(names)))
- # Create a factor of all the characters in the alignment.
- char.levels <- sort(unique(as.vector(aln)))
- names.levels <- names
-
- for (i in 1:length(names)) {
- char.list <- aln[i,]
- name <- names[i]
- # Store the indices of where the gaps are -- we won't plot the gaps.
- gaps <- char.list == '-'
- seq.pos <- seq(1,length(char.list))
- # Get the position and character of each non-gap residue.
- pos.nogaps <- seq.pos[gaps==FALSE]
- char.nogaps <- as.character(char.list[gaps==FALSE])
- # Create a data frame with 1 row per residue to plot.
- df <- rbind(df,data.frame(
- id=factor(x=rep(name,length(pos.nogaps)),levels=names.levels), # Sequence ID to which this residue belongs
- seq_index=rep(i,length(pos.nogaps)), # Index of the containing sequence
- pos=pos.nogaps, # Alignment position
- char=factor(x=char.nogaps,levels=char.levels) # Character of the residue
- ))
- }
-
- # Turn the IDs into a factor to plot along the y axis.
-
- if (!plot.ancestors) {
- # Remove the ancestral nodes from the plot.
- tip.name.indices <- grep("node",names,ignore.case=TRUE,invert=TRUE)
- names <- names[tip.name.indices]
- df$id <- factor(df$id,levels=names)
- df <- subset(df,id %in% names)
- }
- df$type <- 'aln'
-
- ### Add indels to the data frame.
- if (!is.null(indels)) {
- # For each non-gap sequence of each chunk deleted, add a row
- # to the data frame.
- del.df <- data.frame()
- # For each position, store the count of indels in a track.
- max.pos <- max(c(df$pos,indels$pos))
- indel.histogram <- data.frame(
- pos=1:max.pos - 0.5,
- count=0,
- length=0
- )
- for (i in 1:nrow(indels)) {
- row <- indels[i,]
- seqs <- strsplit(row$nongap.str,";")[[1]]
- if (length(seqs) == 0) {
- # Edge case: an indel row came from a column with no aligned sequence,
- # so the 'nongap.str' is empty. Just continue on...
- next()
- }
- cur.df <- data.frame(
- id=seqs,
- pos=row$pos,
- length=row$length,
- type='indel'
- )
- del.df <- rbind(del.df,cur.df)
- # Tick up the histogram.
- indel.histogram[row$pos,]$count <- indel.histogram[row$pos,]$count + 1
- indel.histogram[row$pos,]$length <- indel.histogram[row$pos,]$length + row$length
- }
- # Sync the two data frame's columns, fill empty stuff with NAs.
- columns.from.df <- colnames(df)[!(colnames(df) %in% colnames(del.df))]
- columns.from.del <- colnames(del.df)[!(colnames(del.df) %in% colnames(df))]
- del.df[,columns.from.df] <- NA
- df[,columns.from.del] <- NA
- df <- rbind(df,del.df)
- # Transform the histogram and add it to our tracks.
- max.count <- max(indel.histogram$count)
- max.length <- max(indel.histogram$length)
- indel.histogram$y_lo <- 0
- indel.histogram$score <- indel.histogram$count / (max.count+1)
- indel.histogram$y_hi <- indel.histogram$score
- indel.histogram$id <- 'Hidden Indel Count'
- indel.histogram$height <- 5
- indel.histogram$layout <- 'below'
- indel.histogram$color.gradient <- 'darkblue,darkblue'
- indel.histogram$type <- 'track'
- #indel.len <- indel.histogram
- #indel.len$id <- 'Hidden Indel Length'
- #indel.len$score <- indel.len$length / (max.length+1)
- #indel.len$y_hi <- indel.len$score
- #indel.len$color.gradient <- 'darkgreen,darkgreen'
-
- if (!is.null(tracks)) {
- tracks <- c(tracks,list(indel.histogram))
- } else {
- tracks <- list(indel.histogram)
- }
- }
- if (plot.nongap.bl && is.phylo(phylo)) {
- score.df <- score.aln.columns(phylo,aln)
- max.bl <- max(score.df$score)
- bl.track <- data.frame(
- id = 'Non-gap Branch Length',
- layout = 'below',
- pos = score.df$pos,
- score = score.df$score / max.bl * 0.9,
- y_hi = score.df$score / max.bl * 0.9,
- color.gradient = 'red,black,black'
- )
- if (!is.null(tracks)) {
- tracks <- c(tracks,list(bl.track))
- } else {
- tracks <- list(bl.track)
- }
- }
-
- # Track input format is a list of data frames with one row per feature to be
- # displayed, with the following columns. The only mandatory column is 'pos';
- # all others have sensible default values.
- # ---
- # [fields below are unique per row]
- # pos: The sequence position (starting with 1) of the feature
- # score: The score (between 0 and 1) of the feature. Scores above 1 or below
- # zero will be truncated.
- # y_lo: The lower Y offset (between 0 and 1) of the feature block.
- # y_hi: The upper Y offset (between 0 and 1) of the feature block. These two
- # values allow bars to be positioned along the y-axis.
- # [fields below are unique per track; the value from the first row is used.]
- # id: The display ID for the track.
- # layout: Set to 'above' to put the track above the alignment, 'below' for below.
- # height: The number of alignment rows for the track to span in height.
- # color.gradient: A comma-separated list of colors to interpolate between when coloring
- # the blocks. Examples: 'white,red' 'blue,gray,red' '#FF00FF,#FFFFFF'
- # color: A single color to use for the track, no matter what the scores. Overridden
- # color.gradient if both exist.
- #
- # ---
- #
-
- # What we do is add the track data to the alignment data frame (so it gets
- # paged and scaled properly along with the alignment data) and then separate
- # it out before plotting, so it can be plotted separately from the alignment.
- df$track_index <- -1
- if (!is.null(tracks)) {
- df$score <- NA
- i <- 0
- for (track in tracks) {
- i <- i + 1
- track$track_index <- i
- track$type <- 'track'
-
- # Add default values.
- if (is.null(track$background)) {
- track$background <- 'white'
- }
- if (length(track$pos) == 0) {
- # If no positions are included in the data frame, set the first position
- # to 1 and put the foreground color same as the background.
- # This has the effect of creating a 'spacer' row.
- track$pos <- 1
- track$color <- track[1,]$background
- }
- if (is.null(track$layout)) {
- track$layout <- 'above'
- }
- if (is.null(track$color.gradient)) {
- if (!is.null(track$color)) {
- # No color gradient exists, but we have 'color' instead. Use that.
- track$color.gradient <- paste(track[1,]$color, track[1,]$color, sep=',')
- } else {
- # No color.gradient OR color value exists, so use the default white-to-black.
- track$color.gradient <- 'white,black'
- }
- }
- if (is.null(track$score)) {
- track$score <- 1
- }
- if (is.null(track$y_lo)) {
- track$y_lo = 0
- }
- if (is.null(track$y_hi)) {
- track$y_hi = 1
- }
- if (is.null(track$id)) {
- track$id <- paste('Track',i)
- }
- if (is.null(track$height) || is.na(track$height)) {
- track$height <- 4
- }
- # Ensure that we don't have positions at zero or below.
- track[track$pos <= 0,'pos'] <- 1
- # Limit score range
- track$score <- pmin(track$score,1)
- track$score <- pmax(track$score,0)
-
- # Sync the two data frame's columns, fill empty stuff with NAs.
- columns.from.aln <- colnames(df)[!(colnames(df) %in% colnames(track))]
- columns.from.track <- colnames(track)[!(colnames(track) %in% colnames(df))]
- track[,columns.from.aln] <- NA
- df[,columns.from.track] <- NA
- track$type <- 'track'
- df <- rbind(df,track)
- }
- }
- if (num.pages != 'auto') {
- num.pages <- as.numeric(num.pages)
- }
- aln.length <- max(df$pos)
- chars.per.page <- aln.length
- num.seqs <- length(names)
-
- if (tolower(num.pages) == 'auto' || num.pages > 1) {
- # Plotting a multi-page alignment.
- if (tolower(num.pages) == 'auto') {
- # If we have tracks, add the total track height to the 'num.seqs' variable.
- track.rows <- subset(df,type == 'track')
- if (nrow(track.rows) > 0) {
- track.indices <- sort(unique(track.rows$track_index))
- for (track.index in track.indices) {
- sub.track <- subset(track.rows,track_index==track.index)
- if (nrow(sub.track) > 0) {
- num.seqs <- num.seqs + sub.track[1,]$height
- }
- }
- }
- # Formula to get a square-ish total plot.
- num.pages <- sqrt(aln.length/num.seqs)
- num.pages <- ceiling(num.pages)+1
- # One-page rule for short alignments.
- if (aln.length < 30) {num.pages <- 1}
- }
- # Add a 'page' factor to the data frame and use facet_grid.
- chars.per.page <- ceiling(aln.length / num.pages)
- df$page <- floor((df$pos-1) / chars.per.page) + 1
- if (nrow(subset(df,page <= 0)) > 0) {
- df[df$page <= 0,]$page <- 1 # Fix errors where pos=0.5 goes to page 0
- }
- df$pos <- df$pos - (chars.per.page*(df$page-1))
- page.labels <- paste((0:(num.pages-1))*chars.per.page+1,(1:num.pages)*chars.per.page,sep="-")
- page.numbers <- sort(unique(df$page))
- page.labels <- page.labels[page.numbers]
- df$page <- factor(df$page,levels=page.numbers,labels=page.labels)
- num.pages <- length(page.labels)
- #print(paste("Num pages:",num.pages))
- } else {
- # We've only got one page. Create a factor...
- df$page <- 1
- df$page <- factor(df$page,levels=c(1),labels=c('1'))
- # Store some values which will be used later.
- aln.length <- max(df$pos)
- chars.per.page <- aln.length
- }
- if (color.scheme == 'auto') {
- all.chars <- unlist(aln)
- all.chars <- all.chars[all.chars != '-']
- n.chars <- length(unique(toupper(all.chars)))
- dna <- c('a','t','g','c')
- dna <- c(dna,toupper(dna))
- protein <- letters
- protein <- protein[!(protein %in% c('b','j','o','u','x','z'))]
- protein <- c(protein,toupper(protein))
- ca <- CodonAlphabet()
- nucs <- c('G','A','C','T')
- codon.grid <- expand.grid(nucs,nucs,nucs,stringsAsFactors=F)
- codons <- c()
- for (i in 1:nrow(codon.grid)) {
- codons <- c(codons,paste(codon.grid[i,],collapse=''))
- }
- if(any(all.chars %in% codons)) {
- # If we see any codon alphabets, use the combined_codon
- color.scheme <- 'combined_codon'
- } else {
- # Else just use the combined color scheme. It's good enough!
- color.scheme <- 'combined'
- }
- }
- legend.title <- color.scheme
- # Remove any tracks from the main aln data frame.
- tracks <- subset(df,type=='track')
- indels <- subset(df,type=='indel')
- df <- subset(df,type=='aln')
- # Set positional values in the aln data frame.
- df$xx <- df$pos
- df$yy <- as.numeric(df$id)
- df$xmin <- df$xx - .5
- df$xmax <- df$xx + .5
- df$ymin <- df$yy - .5
- df$ymax <- df$yy + .5
- darken.colors <- FALSE
- if (nrow(indels) > 0) {
- darken.colors <- TRUE
- }
-
- color.map <- alignment.colors(color.scheme,darken=darken.colors)
- df$colors <- color.map[as.character(df$char)]
-
- # Set the base for y-axis configurations.
- y.lim <- c(0,length(names)+1)
- y.breaks <- 1:length(names)
- y.labels <- names
-
- # Create the ggplot panel.
- p <- ggplot(df,aes(x=xx,y=yy,xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax))
- p <- p + geom_rect(aes(fill=colors))
- if (!is.null(aln.xlim)) {
- aln.limits <- aln.xlim
- } else {
- aln.limits <- c(0,max(df$pos)+1)
- }
-
- if (nrow(tracks) > 0) {
- # Handle tracks.
- # These variables hold the current above and below offsets.
- cur.y.above <- length(names) + 1 - .5
- y.lim[2] <- length(names) + 1 - .5
- cur.y.below <- 0.5
- bg.out <- data.frame()
- track.out <- data.frame()
- track.indices <- sort(unique(tracks$track_index))
- for (track.index in track.indices) {
- sub.track <- subset(tracks,track_index==track.index)
- # Collect the track-wide features.
- track.id <- sub.track[1,]$id
- track.height <- sub.track[1,]$height
- track.layout <- sub.track[1,]$layout
- color.gradient <- strsplit(sub.track[1,]$color.gradient,',')[[1]]
- track.bg <- sub.track[1,]$background
- track.ramp <- colorRamp(colors=color.gradient)
- sub.track$colors <- rgb(track.ramp(sub.track$score),maxColorValue=255)
- sub.track$xx <- sub.track$pos
- sub.track$xmin <- sub.track$pos-.5
- sub.track$xmax <- sub.track$pos+.5
- if (track.layout == 'below') {
- sub.track$yy <- cur.y.below
- sub.track$ymin <- cur.y.below - track.height
- sub.track$ymax <- cur.y.below
- y.lim[1] <- cur.y.below - track.height
- y.breaks <- c(y.lim[1] + track.height/2,y.breaks)
- y.labels <- c(as.character(track.id),y.labels)
- # Temporarily store the current min and max for y.
- cur.y.min <- cur.y.below - track.height
- cur.y.max <- cur.y.below
- # Add our track height to the state variable.
- cur.y.below <- cur.y.below - track.height
- } else {
- sub.track$yy <- cur.y.above
- sub.track$ymin <- cur.y.above
- sub.track$ymax <- cur.y.above + track.height
-
- y.lim[2] <- cur.y.above + track.height
- y.breaks <- c(y.breaks,y.lim[2] - track.height/2)
- y.labels <- c(y.labels,as.character(track.id))
- # Temporarily store the current min and max for y.
- cur.y.min <- cur.y.above
- cur.y.max <- cur.y.above + track.height
- # Add our track height to the state variable.
- cur.y.above <- cur.y.above + track.height
- }
- # Adjust bar position if we have y_lo and y_hi values.
- if(length(sub.track$y_lo) > 0) {
- sub.track$ymin <- cur.y.min + track.height*sub.track$y_lo
- }
- if (length(sub.track$y_hi) > 0) {
- sub.track$ymax <- cur.y.min + track.height*sub.track$y_hi
- }
- track.out <- rbind(track.out,sub.track)
- # Create a background rectangle for each page, to place behind the current track.
- pages <- sort(unique(df$page))
- page.indices <- sort(unique(as.integer(df$page)))
- for (page.index in 1:length(pages)) {
- pg <- pages[page.index]
- current.track.bg <- data.frame(
- page=pg,
- colors=track.bg,
- xx=0.5,
- yy=cur.y.min,
- xmin=0.5,
- xmax=chars.per.page + 0.5,
- ymin=cur.y.min,
- ymax=cur.y.max
- )
- bg.out <- rbind(bg.out,current.track.bg)
- }
- }
- # Add layers for the backgrounds and tracks.
- p <- p + geom_rect(aes(fill=colors),data=bg.out)
- p <- p + geom_rect(aes(fill=colors),data=track.out)
- }
- if (nrow(indels) > 0) {
- # Plot indels as small bars on top of characters.
- indel.width <- 0.25
- indels$pos <- indels$pos - .5 # Position indel on seq boundary.
- indels$xx <- indels$pos
- indels$yy <- as.numeric(indels$id)
- indels$xmin <- indels$xx - indel.width
- indels$xmax <- indels$xx + indel.width
- indels$ymin <- indels$yy - .5
- indels$ymax <- indels$yy + .5
- indels$colors <- 'black'
- p <- p + geom_rect(aes(fill=colors),data=indels)
- }
-
- #color.map <- alignment.colors(color.scheme)
- p <- p + scale_fill_identity()
- p <- p + scale_x_continuous(limits=aln.limits,expand=c(0,0))
-
- if (plot.labels) {
- p <- p + scale_y_continuous(limits=y.lim,breaks=y.breaks,labels=y.labels,expand=c(0,0))
- } else {
- p <- p + scale_y_continuous(limits=y.lim,breaks=y.breaks,labels=rep('',length(y.labels)),expand=c(0,0))
- }
- if (aspect.ratio) {
- axis.text.size <- 5
- char.text.size <- 2
- }
-
- if (num.pages > 1) {
- p <- p + facet_grid(page ~ .)
- }
- if (char.text.size == 'auto') {
- char.text.size <- 125 / (num.pages * (num.seqs+1))
- char.text.size <- min(char.text.size,10)
- char.text.size <- max(char.text.size,1)
- }
-
- if (plot.chars) {
- p <- p + geom_text(aes(label=char),colour='black',size=char.text.size)
- }
- if (aspect.ratio) {
- p <- p + coord_equal(ratio=aspect.ratio)
- }
- if (axis.text.size == 'auto') {
- axis.text.size <- 500 / (num.pages * (num.seqs+1))
- axis.text.size <- min(axis.text.size,10)
- axis.text.size <- max(axis.text.size,1)
- }
-
- plot.theme <- theme(
- axis.text.y = element_text(size=axis.text.size,hjust=1),
- axis.title.x = element_blank(),
- axis.title.y = element_blank(),
- axis.ticks = element_blank(),
- panel.grid.minor = element_blank(),
- panel.grid.major = element_blank(),
-
- plot.margin = unit(c(0,0,0,0),'npc')
- )
- p <- p + plot.theme
- if (!plot.legend) {
- p <- p + theme(legend.position='none')
- }
-
- if (plot.tree) {
- if (plot.ancestors) {
- tree.df <- phylo.layout.df(phylo,layout.ancestors=TRUE,align.seq.names=names)
- } else {
- tree.df <- phylo.layout.df(phylo,align.seq.names=names)
- }
- if (num.pages > 1) {
- tree.df$page <- 1
- df.copy <- tree.df
- for (i in 2 : num.pages) {
- df.copy$page <- i
- tree.df <- rbind(tree.df,df.copy)
- }
- }
- max.length <- max.length.to.root(phylo)
- aln.length <- length(aln[1,])
- n.leaves <- length(names)
-
- if (max(tree.df$event.count) > 0) {
- #print(paste("Coloring tree by",color.branches))
- q <- ggplot(tree.df,aes(colour=event.count))
- q <- q + scale_colour_gradient()
- } else {
- q <- ggplot(tree.df)
- }
- q <- q + geom_segment(aes(x=x,y=y,xend=xend,yend=yend))
- if (plot.labels) {
- q <- q + scale_y_continuous(limits=y.lim,breaks=y.breaks,labels=y.labels,expand=c(0,0))
- } else {
- q <- q + scale_y_continuous(limits=y.lim,breaks=y.breaks,labels=rep('',length(y.labels)),expand=c(0,0))
- }
- if (!is.null(tree.xlim)) {
- tree.limits <- tree.xlim
- } else {
- tree.limits <- c(0,max.length)
- }
- q <- q + scale_x_continuous(limits=tree.limits,expand=c(0.05,0))
- q <- q + plot.theme
- # q <- q + theme(plot.margin = unit(c(0,0,0,0),'npc'))
- if (num.pages > 1) {
- q <- q + facet_grid(page ~ .)
- q <- q + theme(strip.text.y=element_blank())
- }
- if (!plot.legend) {
- q <- q + theme(legend.position='none')
- }
- if (aspect.ratio) {
- warning("Aspect ratio set while plotting tree -- the alignment and tree won't line up!")
- }
-
- vplayout(4,1)
- print(p,vp=subplot(2:4,1))
- print(q,vp=subplot(1,1))
- } else {
- p <- p + plot.theme
- print(p)
- }
-
- # this method has no meaningful return value
- },
- private=TRUE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: summary.PhyloSim
- ##
- ###########################################################################/**
- #
- # @RdocMethod summary
- #
- # @title "Summarize the properties of an object"
- #
- # \description{
- # @get "title".
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{object}{An object}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # Returns a PSRootSummary object.
- # }
- #
- # \examples{
- # # Create a PhyloSim object.
- # # Provide the phylo object
- # # and the root sequence.
- # sim<-PhyloSim(
- # name="TinySim",
- # phylo=rcoal(3),
- # root.seq=NucleotideSequence(string="ATG",processes=list(list(JC69())))
- # );
- # # get a summary
- # summary(sim)
- # }
- #
- # @author
- #
- # \seealso{
- # @seeclass
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "summary",
- class="PhyloSim",
- function(
- object,
- ...
- ){
- this<-object;
- this$.summary$"Name"<-this$name;
- this$.summary$"Id"<-this$id;
- if(is.Sequence(this$rootSeq)){
- root.seq<-this$rootSeq$id;
- } else {
- this$.summary$"Root Sequence"<-"undefined";
- }
- if(is.Sequence(this$rootSeq)){
- this$.summary$"Root Sequence big rate"<-this$rootSeq$bigRate;
- }
- if(is.phylo(this$.phylo)){
- this$.summary$"Tree length"<-this$treeLength;
- phylo.details<-grep(pattern="[[:alnum:]]+",x=capture.output(print(this$.phylo)),perl=TRUE,value=TRUE);
- phylo.details<-paste("\n",phylo.details,collapse="",sep="\t");
- this$.summary$"Phylo object details"<-phylo.details;
- } else {
- this$.summary$"Phylo object details"<-"undefined";
- }
- aln<-"undefined";
- if(is.matrix(this$alignment)){
- aln<-"defined";
- }
- this$.summary$"Alignment"<-aln;
- this$.summary$"Log file"<-this$.log.file;
- this$.summary$"Log level"<-this$.log.level;
- NextMethod();
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##### Logging Methods #####
- ##
- ## Method: getLogFile
- ##
- ###########################################################################/**
- #
- # @RdocMethod getLogFile
- #
- # @title "Get the name of the file used for logging"
- #
- # \description{
- # @get "title".
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{A PhyloSim object.}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # A character vector of length one.
- # }
- #
- # \examples{
- # # Create a PhyloSim object
- # sim<-PhyloSim();
- # # get the name of the log file
- # getLogFile(sim)
- # # modify log file name
- # setLogFile(sim,"OldLog.txt")
- # # get/set log file name via virtual field
- # sim$logFile
- # sim$logFile<-"NewLog"
- # sim$logFile
- # }
- #
- # @author
- #
- # \seealso{
- # @seeclass
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "getLogFile",
- class="PhyloSim",
- function(
- this,
- ...
- ){
- this$.log.file;
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: setLogFile
- ##
- ###########################################################################/**
- #
- # @RdocMethod setLogFile
- #
- # @title "Set the name of the file used for logging"
- #
- # \description{
- # @get "title".
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{A PhyloSim object.}
- # \item{value}{The name of the file used for logging.}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # The new logfile.
- # }
- #
- # \examples{
- # # Create a PhyloSim object
- # sim<-PhyloSim();
- # # get the name of the log file
- # getLogFile(sim)
- # # modify log file name
- # setLogFile(sim,"OldLog.txt")
- # # get/set log file name via virtual field
- # sim$logFile
- # sim$logFile<-"NewLog"
- # sim$logFile
- # }
- #
- # @author
- #
- # \seealso{
- # @seeclass
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "setLogFile",
- class="PhyloSim",
- function(
- this,
- value,
- ...
- ){
- if(missing(value)){
- throw("No value provided!\n");
- }
- value<-as.character(value);
- if( length(value) != 1 ){
- throw("The new value must be a character vector of length 1!\n");
- }
- else{
- if( file.access(value,mode=0) == c(0) ){
- warning("The specified file already exists and it will be overwritten during simulation!\n");
- }
- this$.log.file<-value;
- }
- return(this$.log.file);
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: getLogLevel
- ##
- ###########################################################################/**
- #
- # @RdocMethod getLogLevel
- #
- # @title "Get the log level from a PhyloSim object"
- #
- # \description{
- # @get "title".
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{A PhyloSim object.}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # The log level as an integer vector of length one.
- # }
- #
- # \examples{
- # # Create a PhyloSim object
- # sim<-PhyloSim();
- # # get/set log level
- # getLogLevel(sim)
- # setLogLevel(sim,0)
- # # set/get log level via virtual field
- # sim$logLevel<- -1
- # sim$logLevel
- # # clean up
- # unlink(sim$logFile)
- # }
- #
- # @author
- #
- # \seealso{
- # setLogLevel PhyloSim
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "getLogLevel",
- class="PhyloSim",
- function(
- this,
- ...
- ){
- this$.log.level;
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: setLogLevel
- ##
- ###########################################################################/**
- #
- # @RdocMethod setLogLevel
- #
- # @title "Set the log level for a given PhyloSim object"
- #
- # \description{
- # @get "title".
- #
- # No logging is performed if the log level is negative. If the log level is zero, the messages passed to
- # the \code{Log} method will be writen in the log file. If the log level is positive, the messages passed to
- # the \code{Debug} method are saved as well.
- #
- # The default log level is -1. The specified file will be truncated in the case it already exists.
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{A PhyloSim object.}
- # \item{value}{The new log level as an integer.}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # The new level as an integer vector of length one.
- # }
- #
- # \examples{
- # # Create a PhyloSim object
- # sim<-PhyloSim();
- # # get/set log level
- # getLogLevel(sim)
- # setLogLevel(sim,0)
- # # set/get log level via virtual field
- # sim$logLevel<- -1
- # sim$logLevel
- # }
- #
- # @author
- #
- # \seealso{
- # getLogLevel PhyloSim
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "setLogLevel",
- class="PhyloSim",
- function(
- this,
- value,
- ...
- ){
- if(missing(value)){
- throw("No value provided!\n");
- }
- if((!is.numeric(value)) | length(value) != 1 ){
- throw("The new value must be a numeric vector of length 1!\n");
- }
- else{
- # Create/wipe out log file.
- if(value >= 0 ){
- if(file.access(this$.log.file,mode=0) == c(0)){
- warning("The log file already existed and it was wiped out!\n");
- }
- # Creating the associated connection:
- this$.log.connection<-file(paste(this$.log.file),"w+");
- }
- this$.log.level<-value;
- }
- return(this$.log.level);
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: .getMessageTemplate
- ##
- setMethodS3(
- ".getMessageTemplate",
- class="PhyloSim",
- function(
- this,
- ...
- ){
- template<-list(
- time=paste("[",Sys.time(),"]",sep=""),
- level="Info",
- event=""
- );
- return(template);
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: .logMessage
- ##
- setMethodS3(
- ".logMessage",
- class="PhyloSim",
- function(
- this,
- message,
- ...
- ){
- if(missing(message)){
- throw("No message given!\n");
- }
- else if (!is.list(message)){
- throw("The message should be a list");
- }
- else if( length(intersect(names(message),c("time","level","event"))) != 3){
- throw("The \"time\", \"level\" and \"event\" elements are mandatory in the message list!\n");
- }
- else {
- writeLines(paste(message[["time"]]," "),con=this$.log.connection,sep="");
- message[["time"]]<-NULL;
- writeLines(paste(message[["level"]]," "),con=this$.log.connection,sep="");
- message[["level"]]<-NULL;
- writeLines(paste(message[["event"]]," "),con=this$.log.connection,sep="");
- message[["event"]]<-NULL;
- writeLines(paste(message,collapse=", "),con=this$.log.connection,sep="");
- writeLines("\n",con=this$.log.connection,sep="");
- return(TRUE);
- }
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: Log
- ##
- ###########################################################################/**
- #
- # @RdocMethod Log
- #
- # @title "Save a message in the PhyloSim log file"
- #
- # \description{
- # @get "title".
- #
- # The message is written to the log file only if the log level is non-negative. You can use this method for logging
- # in the case you write classes for PhyloSim.
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{A PhyloSim object.}
- # \item{message}{A character vector of length one.}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # The message (invisible).
- # }
- #
- # \examples{
- # # create a PhyloSim object,
- # # with logLevel set to zero
- # sim<-PhyloSim(log.level=0);
- # # log a message
- # Log(sim,"Hiya there!");
- # # close log connection
- # close(sim$.log.connection)
- # # print out the log file
- # cat(paste(scan(file=sim$LogFile,what=character(),sep="\n"),collapse="\n"));cat("\n");
- # # clean up
- # unlink(sim$logFile)
- # }
- #
- # @author
- #
- # \seealso{
- # @seeclass
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "Log",
- class="PhyloSim",
- function(
- this,
- message,
- ...
- ){
-
- if(this$.log.level < 0){
- return(invisible(FALSE))
- }
- if(missing(message)){
- throw("No message given!\n");
- } else {
- template<-.getMessageTemplate(this);
- template$level<-"Info";
- message<-c(template,as.list(message));
- .logMessage(this, message);
- return(invisible(message));
- }
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: Debug
- ##
- ###########################################################################/**
- #
- # @RdocMethod Debug
- #
- # @title "Save a debug message in the PhyloSim log file"
- #
- # \description{
- # @get "title".
- #
- # The debug message is written to the log file only if the log level is non-negative. You can use this method for logging
- # debug messages in the case you write classes for PhyloSim.
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{A PhyloSim object.}
- # \item{message}{A character vector of length one.}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # The message (invisible).
- # }
- #
- # \examples{
- # # create a PhyloSim object,
- # # with logLevel set to zero
- # sim<-PhyloSim(log.level=0);
- # # log a debug message
- # Debug(sim,"Some useful detail...");
- # # close log connection
- # close(sim$.log.connection)
- # # print out the log file
- # cat(paste(scan(file=sim$LogFile,what=character(),sep="\n"),collapse="\n"));cat("\n");
- # # clean up
- # unlink(sim$logFile)
- # }
- #
- # @author
- #
- # \seealso{
- # @seeclass
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "Debug",
- class="PhyloSim",
- function(
- this,
- message,
- ...
- ){
-
- if(missing(message)){
- throw("No message given!\n");
- }
- else if( this$.log.level <= 0){
- return(invisible(FALSE))
- }
- else {
- template<-.getMessageTemplate(this);
- template$level<-"DEBUG";
- message<-c(template,as.list(message));
- .logMessage(this, message);
- return(invisible(message));
- }
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: .UpdateBranchStats
- ##
- setMethodS3(
- ".UpdateBranchStats",
- class="PhyloSim",
- function(
- this,
- event,
- details,
- branch.number,
- ...
- ){
- if(details$type == "substitution"){
-
- if(details$accepted) {
-
- if(is.null(this$.branch.stats[[as.character(branch.number)]]$substitution)){
- this$.branch.stats[[as.character(branch.number)]]$substitution<-1;
- } else {
- this$.branch.stats[[as.character(branch.number)]]$substitution<-(this$.branch.stats[[as.character(branch.number)]]$substitution + 1);
- }
- name<-event$name;
-
- if(is.null(this$.branch.stats[[as.character(branch.number)]][[name]])){
- this$.branch.stats[[as.character(branch.number)]][[name]]<-1;
- }
- else {
- this$.branch.stats[[as.character(branch.number)]][[name]]<-(this$.branch.stats[[as.character(branch.number)]][[name]] + 1);
- }
- }
-
- # Special stuff for the GY94 codon model:
- if(is.GY94(event$.process)){
- # Increment synonymous counter:
- if(event$.type == "synonymous"){
- if(is.null(this$.branch.stats[[as.character(branch.number)]][["nr.syn.subst"]])){
- # First event of this type on this branch, initialize the list element to 1.
- this$.branch.stats[[as.character(branch.number)]][["nr.syn.subst"]]<-1;
- }
- else {
- this$.branch.stats[[as.character(branch.number)]][["nr.syn.subst"]]<-(this$.branch.stats[[as.character(branch.number)]][["nr.syn.subst"]] + 1);
- }
- }
- # Increment non-synonymous counter:
- else if(event$.type == "non-synonymous"){
- if(is.null(this$.branch.stats[[as.character(branch.number)]][["nr.nsyn.subst"]])){
- # First event of this type on this branch, initialize the list element to 1.
- this$.branch.stats[[as.character(branch.number)]][["nr.nsyn.subst"]]<-1;
- }
- else {
- this$.branch.stats[[as.character(branch.number)]][["nr.nsyn.subst"]]<-(this$.branch.stats[[as.character(branch.number)]][["nr.nsyn.subst"]] + 1);
- }
- } else {
- throw("The event generated by the GY94 has no type!\n");
- }
- }
- }
- else if(details$type == "deletion"){
- if(is.null(this$.branch.stats[[as.character(branch.number)]]$deletion)){
- this$.branch.stats[[as.character(branch.number)]]$deletion<-1;
- }
- else {
- this$.branch.stats[[as.character(branch.number)]]$deletion<-(this$.branch.stats[[as.character(branch.number)]]$deletion + 1);
- }
- }
- else if(details$type == "insertion"){
- if(is.null(this$.branch.stats[[as.character(branch.number)]]$insertion)){
- this$.branch.stats[[as.character(branch.number)]]$insertion<-1;
- }
- else {
- this$.branch.stats[[as.character(branch.number)]]$insertion<-(this$.branch.stats[[as.character(branch.number)]]$insertion + 1);
- }
- }
- else {
- throw("Invalid event type!\n");
- }
-
- },
- private=TRUE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: getBranchEvents
- ##
- ###########################################################################/**
- #
- # @RdocMethod getBranchEvents
- #
- # @title "Get the list of events having per-branch statistics recorded"
- #
- # \description{
- # @get "title".
- #
- # During simulation the number of events performed on every branch is recorded. The recorded events can be "basic"
- # events, like "insertion", "deletion" and "A->T" or events which are sums of basic events, like "substituion". The
- # \code{getBranchEvents} method returns a character vector with the names of the events having per-branch
- # statistics recorded. The method should be called after the simulation is finished.
- #
- # The per-branch statistics can be exported as phylo objects by using the \code{exportStatTree} method.
- # The branch lengths of the exported phylo objects are set to the value of the respective per-branch event count.
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{A PhyloSim object.}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # A character vector.
- # }
- #
- # \examples{
- # # Create a PhyloSim object.
- # # Provide the phylo object
- # # and the root sequence.
- #
- # # NOTE: this will be a little bit slow
- # sim<-PhyloSim(
- # phylo=rcoal(3),
- # root.seq=CodonSequence(
- # string="ATGATTATT",
- # processes=list(list(GY94(kappa=2,omega.default=0.5))))
- # );
- # # make the tree longer to have more events
- # scaleTree(sim,5)
- # # plot the tree
- # plot(sim)
- # # run simulation
- # Simulate(sim)
- # # get the list of recorded per-branch event counts
- # getBranchEvents(sim)
- # # export the number of subtitions as a phylo object
- # subst<-exportStatTree(sim,"substitution")
- # # plot the exported phylo object
- # plot(subst)
- # #export the number of synonymous substitutions as a phylo object
- # subst<-exportStatTree(sim,"nr.syn.subst")
- # # plot the exported phylo object
- # plot(subst)
- # }
- #
- # @author
- #
- # \seealso{
- # @seeclass
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "getBranchEvents",
- class="PhyloSim",
- function(
- this,
- ...
- ){
- tmp<-character();
- for(branch in this$.branch.stats){
- tmp<-c(tmp,names(branch));
- }
- return(unique(sort(tmp)));
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: setBranchEvents
- ##
- ###########################################################################/**
- #
- # @RdocMethod setBranchEvents
- #
- # @title "Forbidden action: setting the list of events having per-branch statistics recorded"
- #
- # \description{
- # @get "title".
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{An object.}
- # \item{value}{Not used.}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # Throws an error.
- # }
- #
- # @author
- #
- # \seealso{
- # @seeclass
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "setBranchEvents",
- class="PhyloSim",
- function(
- this,
- value,
- ...
- ){
- virtualAssignmentForbidden(this);
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: exportStatTree
- ##
- ###########################################################################/**
- #
- # @RdocMethod exportStatTree
- #
- # @title "Export the per-branch counts of an event as a phylo object"
- #
- # \description{
- # @get "title".
- #
- # During simulation the number of events performed on every branch is recorded. The recorded events can be "basic"
- # events, like "insertion", "deletion" and "A->T" or events which are sums of basic events, like "substituion". The
- # \code{getBranchEvents} method returns a character vector with the names of the events having per-branch
- # statistics recorded. The method should be called after the simulation is finished.
- #
- # The per-branch statistics can be exported as phylo objects by using the \code{exportStatTree} method.
- # The branch lengths of the exported phylo objects are set to the value of the respective per-branch event count.
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{A PhyloSim object.}
- # \item{event}{The name of the event as returned by the \code{getBranchEvents} method.}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # A phylo object.
- # }
- #
- # \examples{
- # # Create a PhyloSim object.
- # # Provide the phylo object
- # # and the root sequence.
- #
- # # NOTE: this will be a little bit slow
- # sim<-PhyloSim(
- # phylo=rcoal(3),
- # root.seq=CodonSequence(
- # string="ATGATTATT",
- # processes=list(list(GY94(kappa=2,omega.default=0.5)))
- # )
- # );
- # # make the tree longer to have more events
- # scaleTree(sim,5)
- # # plot the tree
- # plot(sim)
- # # run simulation
- # Simulate(sim)
- # # get the list of recorded per-branch event counts
- # getBranchEvents(sim)
- # # export the number of substitutions as a phylo object
- # subst<-exportStatTree(sim,"substitution")
- # # plot the exported phylo object
- # plot(subst)
- # #export the number of synonymous substitutions as a phylo object
- # subst<-exportStatTree(sim,"nr.syn.subst")
- # # plot the exported phylo object
- # plot(subst)
- # }
- #
- # @author
- #
- # \seealso{
- # @seeclass
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "exportStatTree",
- class="PhyloSim",
- function(
- this,
- event,
- ...
- ){
- if(!is.matrix(this$.alignment)){
- throw("Simulation is not complete, cannot export statistics!\n");
- }
- else if(missing(event)){
- throw("No event name specified!\n");
- }
- else if(length(intersect(event, this$branchEvents)) != 1 ){
- throw("Invalid event name!");
- }
- else {
- phylo.copy<-this$phylo;
- phylo.copy$edge.length<-.getStatBrlen(this, event);
- return(phylo.copy);
- }
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: .getStatBrlen
- ##
- setMethodS3(
- ".getStatBrlen",
- class="PhyloSim",
- function(
- this,
- event,
- ...
- ){
- tmp<-numeric();
- for(i in dimnames(this$edges)[[1]]){
- if(is.null(this$.branch.stats[[i]][[event]])){
- tmp[[as.numeric(i)]]<-0;
- }
- else {
- tmp[[as.numeric(i)]]<-this$.branch.stats[[i]][[event]];
- }
- }
- return(tmp);
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##### Phylo object interface methods #####
- ##
- ## Method: getEdges
- ##
- ###########################################################################/**
- #
- # @RdocMethod getEdges
- #
- # @title "Get the edge matrix from a phylo object aggregated by a PhyloSim object"
- #
- # \description{
- # @get "title".
- #
- # The rows of the edge matrix contain the nodes connected by the edge and the edge length.
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{A PhyloSim object.}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # A matrix.
- # }
- #
- # \examples{
- # # create a PhyloSim object
- # sim<-PhyloSim(phylo=rcoal(5));
- # # get the edge matrix
- # getEdges(sim)
- # # get the edge matrix via virtual field
- # sim$edges
- # }
- #
- # @author
- #
- # \seealso{
- # @seeclass
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "getEdges",
- class="PhyloSim",
- function(
- this,
- ...
- ){
- if(!all(is.na(this$.phylo))){
- if(is.phylo(this$.phylo)){
- if(length(this$.phylo$edge.length) > 2){
- if(attr(this$.phylo, "order") != "cladewise"){
- throw("The order of the phylo object is not cladewise! Someone must have been messing with that!\n");
- }
- }
- tmp<-cbind(this$.phylo$edge,this$.phylo$edge.length);
- colnames(tmp)<-c("from","to","length");
- rownames(tmp)<-1:dim(tmp)[[1]];
- return(tmp);
- }
- else{
- throw("The phylo object is invalid!\n");
- }
- }
- else{
- throw("The phylo object is not set!\n");
- }
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: setEdges
- ##
- ###########################################################################/**
- #
- # @RdocMethod setEdges
- #
- # @title "Forbidden action: setting the edge matrix for a phylo object aggregated by a PhyloSim object"
- #
- # \description{
- # @get "title".
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{An object.}
- # \item{value}{Not used.}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # Throws an error.
- # }
- #
- # @author
- #
- # \seealso{
- # @seeclass
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "setEdges",
- class="PhyloSim",
- function(
- this,
- value,
- ...
- ){
- virtualAssignmentForbidden(this);
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: getNtips
- ##
- ###########################################################################/**
- #
- # @RdocMethod getNtips
- #
- # @title "Get the number of the tips form a phylo object aggregated by a PhyloSim object"
- #
- # \description{
- # @get "title".
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{A PhyloSim object}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # A numeric vector of length one.
- # }
- #
- # \examples{
- # # create a PhyloSim object
- # sim<-PhyloSim(phylo=rcoal(5));
- # # get the number of tips
- # getNtips(sim)
- # # get the number of tips via virtual field
- # sim$ntips
- # }
- #
- # @author
- #
- # \seealso{
- # @seeclass
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "getNtips",
- class="PhyloSim",
- function(
- this,
- ...
- ){
- if(!all(is.na(this$.phylo))){
- if(is.phylo(this$.phylo)){
- return(length(this$.phylo$tip.label));
- }
- else{
- throw("The phylo object is invalid!\n");
- }
- }
- else{
- throw("The phylo object is not set!\n");
- }
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: setNtips
- ##
- ###########################################################################/**
- #
- # @RdocMethod setNtips
- #
- # @title "Forbidden action: setting the number of the tips for a phylo object aggregated by a PhyloSim object"
- #
- # \description{
- # @get "title".
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{An object.}
- # \item{value}{Not used.}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # Throws an error.
- # }
- #
- # @author
- #
- # \seealso{
- # @seeclass
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "setNtips",
- class="PhyloSim",
- function(
- this,
- value,
- ...
- ){
- virtualAssignmentForbidden(this);
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: getTipLabels
- ##
- ###########################################################################/**
- #
- # @RdocMethod getTipLabels
- #
- # @title "Get the tip labels from a phylo object aggregated by a PhyloSim object"
- #
- # \description{
- # @get "title".
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{A PhyloSim object.}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # A matrix containing the tip labels.
- # }
- #
- # \examples{
- # # create a PhyloSim object
- # sim<-PhyloSim(phylo=rcoal(5));
- # # get the tip labels
- # getTipLabels(sim)
- # # get the tip lables via virtual field
- # sim$tipLabels
- # }
- #
- # @author
- #
- # \seealso{
- # @seeclass
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "getTipLabels",
- class="PhyloSim",
- function(
- this,
- ...
- ){
- if(!all(is.na(this$.phylo))){
- if(is.phylo(this$.phylo)){
- tmp<-rbind(this$.phylo$tip.label);
- rownames(tmp)<-c("Labels:");
- colnames(tmp)<-c(1:length(tmp));
- return(tmp);
- }
- else{
- throw("The phylo object is invalid!\n");
- }
- }
- else{
- throw("The phylo object is not set!\n");
- }
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: setTipLabels
- ##
- ###########################################################################/**
- #
- # @RdocMethod setTipLabels
- #
- # @title "Forbidden action: setting the tip labels for a phylo object aggregated by a PhyloSim object"
- #
- # \description{
- # @get "title".
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{An object.}
- # \item{value}{Not used.}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # Throws an error.
- # }
- #
- # @author
- #
- # \seealso{
- # @seeclass
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "setTipLabels",
- class="PhyloSim",
- function(
- this,
- value,
- ...
- ){
- virtualAssignmentForbidden(this);
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: getNodes
- ##
- ###########################################################################/**
- #
- # @RdocMethod getNodes
- #
- # @title "Get the node identifiers from a PhyloSim object"
- #
- # \description{
- # @get "title".
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{A PhyloSim object.}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # A numeric vector.
- # }
- #
- # \examples{
- # # create a PhyloSim object
- # sim<-PhyloSim(phylo=rcoal(5));
- # # get the node IDs
- # getNodes(sim)
- # # get the node IDs via virtual field
- # sim$nodes
- # }
- #
- # @author
- #
- # \seealso{
- # @seeclass
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "getNodes",
- class="PhyloSim",
- function(
- this,
- ...
- ){
- if(!all(is.na(this$.phylo))){
- if(is.phylo(this$.phylo)){
- # This is dumb but safe:
- return(sort(unique(as.vector(this$.phylo$edge))));
- }
- else{
- throw("The phylo object is invalid!\n");
- }
- }
- else{
- throw("The phylo object is not set!\n");
- }
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: getNedges
- ##
- ###########################################################################/**
- #
- # @RdocMethod getNedges
- #
- # @title "Get the number of edges from phylo object aggregated by a PhyloSim object"
- #
- # \description{
- # @get "title".
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{A PhyloSim object.}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # A numeric vector of length one.
- # }
- #
- # \examples{
- # # create a PhyloSim object
- # sim<-PhyloSim(phylo=rcoal(5));
- # # get the number of the edges
- # getNedges(sim)
- # # get the number of the edges via virtual field
- # sim$nedges
- # }
- #
- # @author
- #
- # \seealso{
- # @seeclass
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "getNedges",
- class="PhyloSim",
- function(
- this,
- ...
- ){
- if(!all(is.na(this$.phylo))){
- if(is.phylo(this$.phylo)){
- return(dim(this$.phylo$edge)[[1]]);
- }
- else{
- throw("The phylo object is invalid!\n");
- }
- }
- else{
- throw("The phylo object is not set!\n");
- }
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: setNedges
- ##
- ###########################################################################/**
- #
- # @RdocMethod setNedges
- #
- # @title "Forbidden action: setting the number of edges for phylo object aggregated by a PhyloSim object"
- #
- # \description{
- # @get "title".
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{An object.}
- # \item{value}{Not used.}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # Throws an error.
- # }
- #
- # @author
- #
- # \seealso{
- # @seeclass
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "setNedges",
- class="PhyloSim",
- function(
- this,
- value,
- ...
- ){
- virtualAssignmentForbidden(this);
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: setNodes
- ##
- ###########################################################################/**
- #
- # @RdocMethod setNodes
- #
- # @title "Forbidden action: setting the node identifiers for a PhyloSim object"
- #
- # \description{
- # @get "title".
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{An object.}
- # \item{value}{Not used.}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # Throws an error.
- # }
- #
- # @author
- #
- # \seealso{
- # @seeclass
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "setNodes",
- class="PhyloSim",
- function(
- this,
- value,
- ...
- ){
- virtualAssignmentForbidden(this);
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: getTips
- ##
- ###########################################################################/**
- #
- # @RdocMethod getTips
- #
- # @title "Get the node identifiers of the tip nodes from a PhyloSim object"
- #
- # \description{
- # @get "title".
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{A PhyloSim object}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # A numeric vector.
- # }
- #
- # \examples{
- # # create a PhyloSim object
- # sim<-PhyloSim(phylo=rcoal(5));
- # # get the tip IDs
- # getTips(sim)
- # # get the tip IDs via virtual field
- # sim$tips
- # }
- #
- # @author
- #
- # \seealso{
- # @seeclass
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "getTips",
- class="PhyloSim",
- function(
- this,
- ...
- ){
- if(!all(is.na(this$.phylo))){
- if(is.phylo(this$.phylo)){
- # This is dumb but safe:
- #return(sort(unique(as.vector(this$.phylo$edge))));
- return(1:(getNtips(this)));
- }
- else{
- throw("The phylo object is invalid!\n");
- }
- }
- else{
- throw("The phylo object is not set!\n");
- }
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: setTips
- ##
- ###########################################################################/**
- #
- # @RdocMethod setTips
- #
- # @title "Forbidden action: setting the node identifiers of the tip nodes for a PhyloSim object"
- #
- # \description{
- # @get "title".
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{An object.}
- # \item{value}{Not used.}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # Throws an error.
- # }
- #
- # @author
- #
- # \seealso{
- # @seeclass
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "setTips",
- class="PhyloSim",
- function(
- this,
- value,
- ...
- ){
- virtualAssignmentForbidden(this);
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: getRootNode
- ##
- ###########################################################################/**
- #
- # @RdocMethod getRootNode
- #
- # @title "Get the identifier of the root node from a PhyloSim object"
- #
- # \description{
- # @get "title".
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{A PhyloSim object.}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # A numeric vector of length one.
- # }
- #
- # \examples{
- # # create a PhyloSim object
- # sim<-PhyloSim(phylo=rcoal(5));
- # # get the root node ID
- # getRootNode(sim)
- # # get the root node ID via virtual field
- # sim$rootNode
- # }
- #
- # @author
- #
- # \seealso{
- # @seeclass
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "getRootNode",
- class="PhyloSim",
- function(
- this,
- ...
- ){
- if(!all(is.na(this$.phylo))){
- if(is.phylo(this$.phylo)){
- # Relying on cladewise order:
- return(this$.phylo$edge[1,1]);
- }
- else{
- throw("The phylo object is invalid!\n");
- }
- }
- else{
- throw("The phylo object is not set!\n");
- }
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: setRootNode
- ##
- ###########################################################################/**
- #
- # @RdocMethod setRootNode
- #
- # @title "Forbidden action: setting the identifier of the root node for a PhyloSim object"
- #
- # \description{
- # @get "title".
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{An object.}
- # \item{value}{Not used.}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # Throws an error.
- # }
- #
- # @author
- #
- # \seealso{
- # @seeclass
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "setRootNode",
- class="PhyloSim",
- function(
- this,
- value,
- ...
- ){
- virtualAssignmentForbidden(this);
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: is.tip
- ##
- ###########################################################################/**
- #
- # @RdocMethod is.tip
- #
- # @title "Check if a node is a tip"
- #
- # \description{
- # @get "title".
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{A PhyloSim object.}
- # \item{node}{A node identifier (integer vector of length one).}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # TRUE or FALSE
- # }
- #
- # \examples{
- # # create a PhyloSim object
- # sim<-PhyloSim(phylo=rcoal(5));
- # # check if node 4 is a tip
- # is.tip(sim,4)
- # # check if node 6 is a tip
- # is.tip(sim,6)
- # }
- #
- # @author
- #
- # \seealso{
- # @seeclass
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "is.tip",
- class="PhyloSim",
- function(
- this,
- node=NA,
- ...
- ){
- if(missing(node)){
- throw("No node number specified!\n");
- }
- else if(!is.numeric(node)){
- throw("The node number must be numeric!\n");
- }
- else {
- return(round(node) <= this$ntips);
- }
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: getEdge
- ##
- ###########################################################################/**
- #
- # @RdocMethod getEdge
- #
- # @title "Get and edge from the edge matrix"
- #
- # \description{
- # @get "title".
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{A PhyloSim object.}
- # \item{number}{The edge number.}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # The edge as a matrix with a single row.
- # }
- #
- # \examples{
- # # create a PhyloSim object
- # sim<-PhyloSim(phylo=rcoal(5));
- # # get edge number 3
- # getEdge(sim,3)
- # }
- #
- # @author
- #
- # \seealso{
- # @seeclass
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "getEdge",
- class="PhyloSim",
- function(
- this,
- number=NA,
- ...
- ){
- if(missing(number)){
- throw("No object provided!\n");
- }
- else if(!is.numeric(number)){
- throw("The edge number must be numeric!\n");
- }
- else {
- number<-round(number);
- tmp<-rbind(c(this$.phylo$edge[number,],this$.phylo$edge.length[number]));
- colnames(tmp)<-c("from","to","length");
- rownames(tmp)<-c("Edge:");
- return(tmp);
- }
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: getTreeLength
- ##
- ###########################################################################/**
- #
- # @RdocMethod getTreeLength
- #
- # @title "Get the tree length from a PhyloSim object"
- #
- # \description{
- # @get "title".
- #
- # This method retruns the sum of the edge lengths stored in the aggregated phylo object.
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{A PhyloSim object.}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # A numeric vector of length one.
- # }
- #
- # \examples{
- # # create a PhyloSim object
- # sim<-PhyloSim(phylo=rcoal(5));
- # # get the tree length
- # getTreeLength(sim)
- # # get tree length via virtual field
- # sim$treeLength
- # }
- #
- # @author
- #
- # \seealso{
- # @seeclass
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "getTreeLength",
- class="PhyloSim",
- function(
- this,
- ...
- ){
-
- if(!all(is.na(this$.phylo))){
- if(is.phylo(this$.phylo)){
- return(sum(this$.phylo$edge.length));
- }
- else{
- throw("The phylo object is invalid!\n");
- }
- }
- else{
- throw("The phylo object is not set!\n");
- }
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: setTreeLength
- ##
- ###########################################################################/**
- #
- # @RdocMethod setTreeLength
- #
- # @title "Forbidden action: setting the tree length for a PhyloSim object"
- #
- # \description{
- # @get "title".
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{An object.}
- # \item{value}{Not used.}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # Throws an error.
- # }
- #
- # @author
- #
- # \seealso{
- # @seeclass
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "setTreeLength",
- class="PhyloSim",
- function(
- this,
- value,
- ...
- ){
-
- virtualAssignmentForbidden(this);
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: scaleTree
- ##
- ###########################################################################/**
- #
- # @RdocMethod scaleTree
- #
- # @title "Scale the branch lengths of a phylo object aggragted by a PhyloSim object"
- #
- # \description{
- # @get "title".
- # This method multiples all the edge lengths by the specified factor.
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{A PhyloSim object.}
- # \item{factor}{A numeric vector of length one.}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # The PhyloSim object (invisible).
- # }
- #
- # \examples{
- # # create a PhyloSim object
- # sim<-PhyloSim(phylo=rcoal(5));
- # # get the tree length
- # sim$treeLength
- # # scale tree
- # scaleTree(sim,10)
- # # get the scaled tree length
- # sim$treeLength
- # }
- #
- # @author
- #
- # \seealso{
- # @seeclass
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "scaleTree",
- class="PhyloSim",
- function(
- this,
- factor,
- ...
- ){
- if(missing(factor)){
- throw("No branch length scaling factor specified!\n");
- } else if((!is.numeric(factor)) | (length(factor) != 1)){
- throw("The scaling factor must be a numeric vector of length 1!\n");
- } else if(!is.phylo(this$.phylo)){
- throw("The phylo object is not set or it is invalid!\n");
- } else {
- this$.phylo$edge.length<-(this$.phylo$edge.length * factor);
- return(invisible(this));
-
- }
-
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ###########################################################################/**
- #
- # @RdocMethod readAlignment
- #
- # @title "Read alignment from file"
- #
- # \description{
- # @get "title".
- #
- # This method reads an alignment by using the \code{read.dna} function from the \code{\link{ape}}
- # package and stores in the \code{PhyloSim} object. If a tree is already attached to the \code{PhyloSim}
- # object, the alignment must at least contain the sequences corresponding to tip nodes (but it
- # may also contain additional ancestral sequences).
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{A PhyloSim object.}
- # \item{file}{A file name specified by either a variable of mode character, or a double-quoted string.}
- # \item{format}{a character string specifying the format of the DNA sequences. Four choices are possible: "interleaved", "sequential", "clustal", or "fasta", or any unambiguous abbreviation of these.}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # The PhyloSim object (invisible).
- # }
- #
- # \examples{
- # # get a safe file name
- # fname<-paste("PhyloSim_dummy_fas_",Sys.getpid(),sep="")
- # # write out a fasta alignment
- # cat("> t3\nGTCTTT-CG-\n",file=fname);
- # cat("> t4\nG--TC-TCGG\n",file=fname,append=TRUE);
- # cat("> t2\nG--TC-TCGG\n",file=fname,append=TRUE);
- # cat("> t1\nGTC-G-TCGG",file=fname,append=TRUE);
- # # construct a PhyloSim object,
- # # set the phylo object
- # sim<-PhyloSim(phylo=rcoal(4))
- # # read the alignment
- # readAlignment(sim,fname)
- # # remove alignment file
- # unlink(fname)
- # # plot the tree & alignment
- # plot(sim)
- # }
- #
- # @author
- #
- # \seealso{
- # @seeclass
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "readAlignment",
- class="PhyloSim",
- function(
- this,
- file,
- format="fasta",
- ...
- ){
- aln<-toupper(read.dna(file=file,format=format,as.matrix=TRUE,as.character=TRUE));
- aln.names<-dimnames(aln)[[1]];
- if (!all(is.na(this$.phylo))) {
- tip.labels<-this$tipLabels;
- length.overlap <- length(intersect(tip.labels,aln.names))
- if(length.overlap != length(tip.labels)){
- throw("The alignment must contain all sequences corresponding to tip nodes!");
- }
- if (length(aln.names) > length(tip.labels)) {
- warning("Alignment has more sequences than the tree's leaf count -- either it contains ancestral sequences or something is wrong!")
- }
- }
- this$.alignment<-aln;
- return(invisible(this));
-
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ##
- ## Method: getAlignmentLength
- ##
- ###########################################################################/**
- #
- # @RdocMethod getAlignmentLength
- #
- # @title "Get the alignment length from a PhyloSim object"
- #
- # \description{
- # @get "title".
- #
- # This method retruns the number of columns in the alignment stored in the PhyloSim object.
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{A PhyloSim object.}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # A numeric vector of length one.
- # }
- #
- # \examples{
- # # create a PhyloSim object and run a simulation:
- # sim<-Simulate(
- # PhyloSim(phy=rcoal(3),
- # root=NucleotideSequence(string="ATGC", proc=list(list(JC69())) ) )
- # )
- # # get the alignment length
- # getAlignmentLength(sim)
- # }
- #
- # @author
- #
- # \seealso{
- # @seeclass
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "getAlignmentLength",
- class="PhyloSim",
- function(
- this,
- ...
- ){
-
- if(!all(is.na(this$.alignment))){
- return(dim(this$.alignment)[2]);
- }
- else{
- throw("The alignment object is not set!\n");
- }
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ###########################################################################/**
- #
- # @RdocMethod readTree
- #
- # @title "Read tree from file"
- #
- # \description{
- # @get "title".
- #
- # This method reads a tree by using the \code{read.tree} function from the \code{\link{ape}}
- # package and stores in the \code{PhyloSim} object. If an alignment is already attached
- # to the \code{PhyloSim} object, it must contain all sequences corresponding to tip nodes.
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{this}{A PhyloSim object.}
- # \item{file}{A file name specified by either a variable of mode character, or a double-quoted string.}
- # \item{...}{Not used.}
- # }
- #
- # \value{
- # The PhyloSim object (invisible).
- # }
- #
- # \examples{
- # # get a safe file name
- # fname<-paste("PhyloSim_dummy_fas_",Sys.getpid(),sep="")
- # # write out a fasta alignment
- # cat("(a,(b,c));",file=fname);
- # # construct a PhyloSim object:
- # sim<-PhyloSim()
- # # read the alignment
- # readTree(sim,fname)
- # # remove alignment file
- # unlink(fname)
- # # plot the tree
- # plot(sim)
- # }
- #
- # @author
- #
- # \seealso{
- # @seeclass
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "readTree",
- class="PhyloSim",
- function(
- this,
- file,
- ...
- ){
- tree <- read.tree(file)
- if (!any(is.na(this$.alignment))) {
- # Check for overlap between leaves and seqs.
- aln <- this$.alignment;
- aln.names <- dimnames(aln)[[1]];
- tip.labels <- tree$tip.label;
- aln.tree.overlap <- length(intersect(tip.labels,aln.names))
- if(aln.tree.overlap != length(tip.labels)){
- throw("The alignment must contain all sequences corresponding to tip nodes!");
- }
- if (length(aln.names) > length(tip.labels)) {
- warning("Alignment has more sequences than the tree's leaf count -- either it contains ancestral sequences or something is wrong!")
- }
- }
-
- this$.phylo <- tree
- return(invisible(this));
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );
- ###########################################################################/**
- #
- # @RdocMethod Undocumented
- # \alias{Undocumented}
- # \alias{newMatrix}
- # \alias{setEquDist.CodonSubst}
- # \alias{BrownianPath}
- # \alias{buildFromPAML}
- # \alias{checkConsistency}
- # \alias{clearStates}
- # \alias{areSynonymous}
- # \alias{attachHookToNode}
- # \alias{attachProcess}
- # \alias{attachSeqToNode}
- # \alias{copySubSequence}
- # \alias{Debug}
- # \alias{deleteSubSequence}
- # \alias{detachHookFromNode}
- # \alias{detachProcess}
- # \alias{detachSeqFromNode}
- # \alias{enableVirtual}
- # \alias{exportStatTree}
- # \alias{flagTotalRate}
- # \alias{generateInsert}
- # \alias{getAcceptBy}
- # \alias{getAcceptWin}
- # \alias{getAlignment}
- # \alias{getAlphabet}
- # \alias{getAlphabets}
- # \alias{getAlphabet}
- # \alias{getAncestral}
- # \alias{getAncestral}
- # \alias{getBaseFreqs}
- # \alias{getBigRate}
- # \alias{getBranchEvents}
- # \alias{getCodonFreqs}
- # \alias{getComments}
- # \alias{getCumulativeRates}
- # \alias{getCumulativeRatesFromRange}
- # \alias{getDeletionTolerance}
- # \alias{getDist}
- # \alias{getEdge}
- # \alias{getEdges}
- # \alias{getEquDist}
- # \alias{getEventRate}
- # \alias{getEventRateAtSite}
- # \alias{getEvents}
- # \alias{getEventsAtSite}
- # \alias{getGenerateBy}
- # \alias{getHandler}
- # \alias{getId}
- # \alias{getInsertHook}
- # \alias{getInsertionTolerance}
- # \alias{getInsertionTolerance}
- # \alias{getKappa}
- # \alias{getLength}
- # \alias{getLengthParam1}
- # \alias{getLengthParam2}
- # \alias{getLogFile}
- # \alias{getLogLevel}
- # \alias{getMatrix}
- # \alias{getMaxLength}
- # \alias{getMethodsList}
- # \alias{getNedges}
- # \alias{getNodes}
- # \alias{getNtips}
- # \alias{getOmegas}
- # \alias{getParameterAtSite}
- # \alias{getParameterAtSites}
- # \alias{getPhylo}
- # \alias{getProbs}
- # \alias{getProcess}
- # \alias{getProcesses}
- # \alias{getProposeBy}
- # \alias{getQMatrix}
- # \alias{getRate}
- # \alias{getRateList}
- # \alias{getRateMultipliers}
- # \alias{getRateParam}
- # \alias{getRateParamList}
- # \alias{getRootNode}
- # \alias{getRootSeq}
- # \alias{getScale}
- # \alias{getScaledMatrix}
- # \alias{getSeqFromNode}
- # \alias{getSequence}
- # \alias{getSequences}
- # \alias{getSite}
- # \alias{getSites}
- # \alias{getSiteSpecificParamIds}
- # \alias{getSiteSpecificParamList}
- # \alias{getSize}
- # \alias{getSizes}
- # \alias{getState}
- # \alias{getStates}
- # \alias{getString}
- # \alias{getSymbolFreqs}
- # \alias{getSymbolLength}
- # \alias{getSymbols}
- # \alias{getTableId}
- # \alias{getTargetState}
- # \alias{getTemplateSeq}
- # \alias{getTheta}
- # \alias{getTipLabels}
- # \alias{getTips}
- # \alias{getToleranceMargin}
- # \alias{getTotalRate}
- # \alias{getTotalRates}
- # \alias{getTotalRatesFromRange}
- # \alias{getTransTable}
- # \alias{getTreeLength}
- # \alias{getType}
- # \alias{getUniqueAlphabets}
- # \alias{getUniqueProcesses}
- # \alias{getWriteProtected}
- # \alias{globalConsistencyCheck}
- # \alias{hasSiteSpecificParameter}
- # \alias{hasSymbols}
- # \alias{hasUndefinedRate}
- # \alias{insertSequence}
- # \alias{intersect}
- # \alias{list}
- # \alias{isAttached}
- # \alias{isEmpty}
- # \alias{is}
- # \alias{default}
- # \alias{isStartCodon}
- # \alias{isStopCodon}
- # \alias{is}
- # \alias{tip}
- # \alias{Log}
- # \alias{my}
- # \alias{all}
- # \alias{equal}
- # \alias{newAAMatrix}
- # \alias{omegaHist}
- # \alias{omegaVarM0}
- # \alias{omegaVarM1}
- # \alias{omegaVarM10Cont}
- # \alias{omegaVarM10Cont}
- # \alias{omegaVarM2}
- # \alias{omegaVarM3}
- # \alias{omegaVarM4}
- # \alias{Perform}
- # \alias{plotParametersAtSites}
- # \alias{plot}
- # \alias{plusGamma}
- # \alias{plusInvGamma}
- # \alias{proposeLength}
- # \alias{rescaleQMatrix}
- # \alias{sampleState}
- # \alias{sampleStates}
- # \alias{saveAlignment}
- # \alias{Scale}
- # \alias{scaleTree}
- # \alias{setAcceptBy}
- # \alias{setAcceptWin}
- # \alias{setAlignment}
- # \alias{setAlphabet}
- # \alias{setAlphabets}
- # \alias{setAncestral}
- # \alias{setBaseFreqs}
- # \alias{setBigRate}
- # \alias{setBranchEvents}
- # \alias{setCodonFreqs}
- # \alias{setCodonFreqs}
- # \alias{setComments}
- # \alias{setCumulativeRates}
- # \alias{setDeletionTolerance}
- # \alias{setDist}
- # \alias{setEdges}
- # \alias{setEquDist}
- # \alias{setEquDist}
- # \alias{setEvents}
- # \alias{setGenerateBy}
- # \alias{setHandler}
- # \alias{setId}
- # \alias{setInsertHook}
- # \alias{setInsertionTolerance}
- # \alias{setInsertionTolerance}
- # \alias{setKappa}
- # \alias{setKappa}
- # \alias{setLength}
- # \alias{setLengthParam1}
- # \alias{setLengthParam2}
- # \alias{setLogFile}
- # \alias{setLogLevel}
- # \alias{setMatrix}
- # \alias{setMaxLength}
- # \alias{setMethodsList}
- # \alias{setName}
- # \alias{setNedges}
- # \alias{setNodes}
- # \alias{setNtips}
- # \alias{setOmegas}
- # \alias{setParameterAtSite}
- # \alias{setParameterAtSites}
- # \alias{setPhylo}
- # \alias{setPosition}
- # \alias{setProbs}
- # \alias{setProcess}
- # \alias{setProcesses}
- # \alias{setProposeBy}
- # \alias{setQMatrix}
- # \alias{setRate}
- # \alias{setRateList}
- # \alias{setRateMultipliers}
- # \alias{setRateParam}
- # \alias{setRateParamList}
- # \alias{setRootNode}
- # \alias{setRootSeq}
- # \alias{setScale}
- # \alias{setScaledMatrix}
- # \alias{setSequence}
- # \alias{setSequences}
- # \alias{setSite}
- # \alias{setSiteSpecificParamIds}
- # \alias{setSiteSpecificParamList}
- # \alias{setSize}
- # \alias{setSizes}
- # \alias{setState}
- # \alias{setStates}
- # \alias{setString}
- # \alias{setSymbolLength}
- # \alias{setSymbols}
- # \alias{setTableId}
- # \alias{setTargetState}
- # \alias{setTemplateSeq}
- # \alias{setTheta}
- # \alias{setTipLabels}
- # \alias{setTips}
- # \alias{setToleranceMargin}
- # \alias{setTotalRate}
- # \alias{setTotalRates}
- # \alias{setTransTable}
- # \alias{setTreeLength}
- # \alias{setType}
- # \alias{setUniqueAlphabets}
- # \alias{setUniqueProcesses}
- # \alias{setWriteProtected}
- # \alias{Simulate}
- # \alias{Translate}
- # \alias{translateCodon}
- # \alias{virtualAssignmentForbidden}
- # \alias{intersect.list}
- # \alias{is.tip}
- # \alias{my.all.equal}
- # \alias{plot.PSRoot}
- # \alias{revComp}
- # \alias{readAlignment}
- # \alias{readTree}
- # \alias{getOmegaScalingFactor}
- # \alias{saveLoadReference}
- # \alias{getAlignmentLength}
- #
- # @title "Undocumented object (PhyloSim package)"
- #
- # \description{
- # @get "title".
- #
- # See the corresponding specific methods if applicable.
- # }
- #
- # @synopsis
- #
- # \arguments{
- # \item{...}{Not used.}
- # }
- #
- # @author
- #
- # \seealso{
- # @seeclass
- # }
- #
- #*/###########################################################################
- setMethodS3(
- "Undocumented",
- class="PhyloSim",
- function(
- ...
- ){
- cat("This method has no documentation!\n");
- },
- private=FALSE,
- protected=FALSE,
- overwrite=FALSE,
- conflict="warning",
- validators=getOption("R.methodsS3:validators:setMethodS3")
- );