Nothing
R/FieldDeletor.R
In phylosim: Flexible Simulations of Biological Sequence Evolution
##
## Copyright 2009 Botond Sipos
## See the package description for licensing information.
##
##########################################################################/**
#
# @RdocClass FastFieldDeletor
#
# @title "The FastFieldDeletor class"
#
# \description{
# The \code{\link{DiscreteDeletor}} and \code{\link{ContinuousDeletor}} processes
# generate deletion event objects with rates determined by the general rate of the
# process and the "rate.multiplier" parameter. The probability of rejecting an event
# is determined by the product of the "deletion.tolerance" parameters from the affected sites.
# If the majority of the sites have low deletion tolerance most of the events are rejected, which
# slows down the simulation without performing much events.
#
# The \code{FastFieldDeletor} process scales down the rate and length distribution of the proposed
# events based on the highest insertion tolerance value observed in the root sequence
# (see the package vignette for details), thus making the simulation more efficient.
#
# The available length distributions are (see also the package vignette):
# \itemize{
# \item Geometric (default) - \code{lengthParam1} is \emph{Lambda}.
# \item Poisson+1 - \code{lengthParam1} is \emph{Lambda}.
# \item Conway-Maxwell Poisson+1 - \code{lengthParam1} is \emph{Lambda}, \code{lengthParam2} is \emph{nu}
# \item Negative Binomial+1 - \code{lengthParam1} is \emph{Lambda}, \code{lengthParam2} is \emph{r}
# }
#
# Insertion proceses can insert sites with deletion tolerance higher than the largest
# deletion tolerance observed in the root sequence. The user can specify the largest expected
# tolerance value through the \code{toleranceMargin} virtual field. The process is then scaled by
# max(initial_highest_tolerance, \code{toleranceMargin}).
#
# @classhierarchy
# }
#
# @synopsis
#
# \arguments{
# \item{name}{Object name.}
# \item{type}{Type of the process (see above).}
# \item{length.param.1}{Object name.}
# \item{length.param.2}{Object name.}
# \item{tolerance.margin}{Object name.}
# \item{...}{Additional arguments.}
# }
#
# \section{Fields and Methods}{
# @allmethods
# }
#
# \examples{
# # create a FastFieldDeletor object, default (geometric) type
# # length.param.1 is "lambda"
# p<-FastFieldDeletor(rate=1,length.param.1=0.9, tolerance.margin=0.8)
# # get type
# p$type
# # get object summary
# summary(p)
# # set/get lambda
# p$lengthParam1<-0.8
# p$lengthParam1
# # set/get tolerance margin
# p$toleranceMargin<-0.5
# p$toleranceMargin
# # create a nucleotide sequence, attach process
# s<-NucleotideSequence(length=30,processes=list(list(p)))
# # set state pattern
# s$states<-c("A","A","T","T","G","G","C","C")
# # sample deletion tolerances
# setDeletionTolerance(s,p,sample(seq(from=0,to=0.8,by=0.01),30,replace=TRUE))
# # plot deletion tolerance
# plotParametersAtSites(s,p,"deletion.tolerance")
# # simulate
# sim<-PhyloSim(root.seq=s, phylo=rcoal(2))
# Simulate(sim)
# # show alignment
# sim$alignment
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setConstructorS3(
"FastFieldDeletor",
function(
name="Anonymous",
type="geometric",
length.param.1=NA, # mostly "Lambda"
length.param.2=NA,
tolerance.margin=0, # minimum tolerance value used for scaling.
...
) {
# supported types
ALLOWED.TYPES=c("geometric","poisson","neg.binomial","compoisson");
# Creating a GeneralDeletor Process.
this<-GeneralDeletor(
...
);
# Check if the specified type is valid:
if(length(intersect(ALLOWED.TYPES,type)) != 1){
throw("The specified field model type is invalid!\n");
}
# Extending as FastFieldDeletor:
this<-extend(
this,
"FastFieldDeletor",
.type=type, # field model flavour
.tolerance.margin=NA, # minimum tolerance used for scaling
.tolerance.max=NA, # maximum tolerance obseved at first call of
# getEventAtSites
.d=NA, # is max(.tolerance.max, .tolerance.margin)
.field.scaling.factor=NA,# the precalculated scaling factor
.length.param.1=NA, # mostly "Lambda"
.length.param.2=NA,
.ALLOWED.TYPES=ALLOWED.TYPES # supported types
);
# Set length parameter 1 if not missing:
if(!missing(length.param.1)){
this$lengthParam1<-length.param.1;
}
# Set length parameter 2 if not missing:
if(!missing(length.param.2)){
this$lengthParam2<-length.param.2;
}
# Set tolerance margin:
setToleranceMargin(this, tolerance.margin);
# Using virtual field to clear Id cache:
this$name<-name;
# Set the function proposing deletion lengths:
this$proposeBy<-function(process=NA,seq=NA,pos=NA){
# Check the length parameters:
.checkLengthParams(this);
# Check if this$.d is defined:
if(is.na(this$.d)){
throw("thid$.d is NA! This shouldn't happen!");;
}
d<-this$.d;
# Type specific length sampling expressions:
# Geometric + 1:
if(this$.type == "geometric"){
express<-expression(1 + rgeom(1,prob=( 1 - ( this$.length.param.1 * d) ) ) );
}
# Poisson+1:
else if(this$.type == "poisson"){
express<-expression( 1 + rpois(1,lambda=(this$.length.param.1 * d) ) );
}
# Negative Binomial + 1:
else if(this$.type == "neg.binomial"){
express<-expression(1 + rnbinom(1,this$.length.param.2,prob=( 1 - ( this$.length.param.1 * d))) );
}
# Conway-Maxwell Poisson + 1:
else if(this$.type == "compoisson"){
express<-expression(1 + rcom(1,lambda=( this$.length.param.1 * d), nu = this$.length.param.2));
}
return( round( eval(express) ) );
} # /proposeBy
# Set the function performing the accept/reject step:
this$acceptBy<-function(process=NA,sequence=NA,range=NA){
del.tol<-c();
for(site in sequence$.sites[range]){
# Reject if the range contains a site which is not attached to
# the process:
if(!isAttached(site, process)){
return(FALSE);
}
del.tol<-c(del.tol, getParameterAtSite(process, site, "deletion.tolerance")$value);
}
# Calculate acceptance probability:
accept.prob<-( prod(as.numeric(del.tol)) / this$.d );
# Accept/reject:
return( sample(c(TRUE,FALSE),replace=FALSE,prob=c(accept.prob,(1-accept.prob)),size=1) );
}
return(this);
},
enforceRCC=TRUE
);
##
## Method: .checkLengthParams
##
setMethodS3(
".checkLengthParams",
class="FastFieldDeletor",
function(
this,
...
){
# Check length parameter 1:
if(is.na(this$.length.param.1)){
throw("Length parameter 1 is NA! Cannot generate events or propose lengths!\n");
}
# Check length parameter 2:
if(length(intersect(c("neg.binomial","compoisson"),this$.type)) != 0) {
if(is.na(this$.length.param.2)){
throw("Length parameter 1 is NA! Cannot generate events or propose lengths!\n");
}
}
},
private=TRUE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## 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="FastFieldDeletor",
function(
this,
...
){
wp<-this$writeProtected;
if (wp) {
this$writeProtected<-FALSE;
}
may.fail<-function(this) {
# Check if the type is valid:
if(length(intersect(this$.ALLOWED.TYPES, this$.type)) != 1){
throw("The specified field model type is invalid!\n");
}
if((!is.numeric(this$.tolerance.margin) & !is.na(this$.tolerance.margin))){
throw("Tolerance margin is invalid!\n");
}
if((!is.numeric(this$.tolerance.max) & !is.na(this$.tolerance.max))){
throw(".tolerance.max is invalid!\n");
}
if((!is.numeric(this$.d) & !is.na(this$.d))){
throw(".d is invalid!\n");
} else if(!is.na(this$.tolerance.margin) & !is.na(this$.tolerance.max)) {
if(this$.d != max(this$.tolerance.margin, this$.tolerance.max)){
throw(".d is inconsistent!\n");
}
}
if((!is.numeric(this$.field.scaling.factor) & !is.na(this$.field.scaling.factor))){
throw(".field.scaling.factor is invalid!\n");
}
if((!is.numeric(this$.length.param.1) & !is.na(this$.length.param.1))){
throw("Length parameter 1 is invalid!\n");
}
if((!is.numeric(this$.length.param.2) & !is.na(this$.length.param.2))){
throw("Length parameter 2 is invalid!\n");
}
}
tryCatch(may.fail(this),finally=this$writeProtected<-wp);
NextMethod();
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: summary.FastFieldDeletor
##
###########################################################################/**
#
# @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 an object
# a<-FastFieldDeletor()
# # get a summary
# summary(a)
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"summary",
class="FastFieldDeletor",
function(
object,
...
){
this<-object;
.addSummaryNameId(this);
this$.summary$"Type"<-this$.type;
this$.summary$"Tolerance margin"<-this$.tolerance.margin;
this$.summary$"Length parameter 1"<-this$.length.param.1;
this$.summary$"Length parameter 2"<-this$.length.param.2;
this$.summary$"Scaling factor"<-this$.field.scaling.factor;
NextMethod();
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getEventsAtSite
##
###########################################################################/**
#
# @RdocMethod getEventsAtSite
#
# @title "Generate a deletion event object given the state of the target site"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A FastFieldDeletor object.}
# \item{target.site}{A Site object.}
# \item{...}{Not used.}
# }
#
# \value{
# A list of Event objects.
# }
#
#
# @author
#
# \seealso{
# GeneralDeletor getEventsAtSite.GeneralDeletor
# }
#
#*/###########################################################################
setMethodS3(
"getEventsAtSite",
class="FastFieldDeletor",
function(
this,
target.site,
...
){
if(!exists(x="PSIM_FAST")){
if(missing(target.site)) {
throw("No target site provided!\n");
}
if(!is.Site(target.site)) {
throw("Target site invalid!\n");
}
else if(!is.function(this$.propose.by)) {
throw("proposeBy is not set, cannot propose deletion!\n");
}
else if (!is.function(this$.accept.by)){
throw("acceptBy is not set, cannot generate deletion event deletion!\n");
}
# Complain if sequence has a zero length:
if(target.site$.sequence$.length == 0) {
throw("Sequence has zero length so there is nothing to delete! How did you get here anyway?\n");
}
}
# Clone the event template object:
deletion.event<-clone(this$.event.template);
# Set the target position passed in a temporary field:
deletion.event$.position<-target.site$.position;
# Set the target site:
deletion.event$.site<-target.site;
# Set event name:
deletion.event$.name<-"Deletion";
# Set the generator process:
deletion.event$.process<-this;
# Calculate the field model specific scaling factor if it is not yet calculated:
if(is.na(this$.field.scaling.factor)){
this$.field.scaling.factor<-.getScalingFactor(this,process=this,seq=target.site$.sequence);
}
# Event rate is the product of the general rate, the field model scaling factor and the
# site specific rate multiplier:
deletion.event$.rate<-(this$.rate * (target.site$.processes[[this$.id]]$site.params[["rate.multiplier"]]$value) * this$.field.scaling.factor);
# Set the handler for the deletion event:
deletion.event$.handler<-this$.handler.template;
# Return the event object in a list:
list(deletion.event);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getType
##
###########################################################################/**
#
# @RdocMethod getType
#
# @title "Get the type of a FastFieldDeletor object"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A FastFieldDeletor object.}
# \item{...}{Not used.}
# }
#
# \value{
# A character vector of length one.
# }
#
# \examples{
# # create a FastFieldDeletor, default type (geometric)
# p<-FastFieldDeletor()
# # get type
# getType(p)
# # create a FastFieldDeletor, poisson type
# p<-FastFieldDeletor(type="poisson")
# p$type
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"getType",
class="FastFieldDeletor",
function(
this,
...
){
this$.type;
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: setType
##
###########################################################################/**
#
# @RdocMethod setType
#
# @title "Forbidden action: setting the type of a FastFieldDeletor object"
#
# \description{
# @get "title".
#
# The type can be set only through the \code{type} constructor argument.
# }
#
# @synopsis
#
# \arguments{
# \item{this}{An object.}
# \item{value}{Not used.}
# \item{...}{Not used.}
# }
#
# \value{
# Throws an error.
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"setType",
class="FastFieldDeletor",
function(
this,
value,
...
){
throw("The type of the FastFieldDeletor process cannot be modified. Please set it by the constructor argument.");
},
private=TRUE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getLengthParam1
##
###########################################################################/**
#
# @RdocMethod getLengthParam1
#
# @title "Get the first length parameter"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A FastFieldDeletor object.}
# \item{...}{Not used.}
# }
#
# \value{
# A numeric vector of length one.
# }
#
# \examples{
# # create a geometric FastFieldDeletor
# p<-FastFieldDeletor()
# # set/get the first length parameter
# setLengthParam1(p,0.5)
# getLengthParam1(p)
# # set/get the first length parameter via virtual field
# p$lengthParam1<-0.2
# p$lengthParam1
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"getLengthParam1",
class="FastFieldDeletor",
function(
this,
...
){
this$.length.param.1;
;
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getLengthParam2
##
###########################################################################/**
#
# @RdocMethod getLengthParam2
#
# @title "Get the second length parameter"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A FastFieldDeletor object.}
# \item{...}{Not used.}
# }
#
# \value{
# A numeric vector of length one.
# }
#
# \examples{
# # create a geometric FastFieldDeletor
# p<-FastFieldDeletor()
# # set/get the second length parameter
# setLengthParam2(p,0.5)
# getLengthParam2(p)
# # set/get the second length parameter via virtual field
# p$lengthParam2<-0.2
# p$lengthParam2
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"getLengthParam2",
class="FastFieldDeletor",
function(
this,
...
){
this$.length.param.2;
;
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: setLengthParam1
##
###########################################################################/**
#
# @RdocMethod setLengthParam1
#
# @title "Set the first length parameter"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A FastFieldDeletor object.}
# \item{value}{A numeric vector of length one.}
# \item{...}{Not used.}
# }
#
# \value{
# value (invisible).
# }
#
# \examples{
# # create a geometric FastFieldDeletor
# p<-FastFieldDeletor()
# # set/get the first length parameter
# setLengthParam1(p,0.5)
# getLengthParam1(p)
# # set/get the first length parameter via virtual field
# p$lengthParam1<-0.2
# p$lengthParam1
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"setLengthParam1",
class="FastFieldDeletor",
function(
this,
value,
...
){
if(missing(value)){
throw("No new length parameter value specified!\n");
}
else if ((!is.numeric(value)) | (length(value) != 1 ) ) {
throw("The new value must be a numeric vector of length 1!\n");
}
else {
# First set the scaling factor to NA to force the recalculation:
this$.field.scaling.factor<-NA;
this$.length.param.1<-value;
}
return(invisible(value));
},
private=TRUE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: setLengthParam2
##
###########################################################################/**
#
# @RdocMethod setLengthParam2
#
# @title "Set the second length parameter"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A FastFieldDeletor object.}
# \item{value}{A numeric vector of length one.}
# \item{...}{Not used.}
# }
#
# \value{
# value (invisible).
# }
#
# \examples{
# # create a geometric FastFieldDeletor
# p<-FastFieldDeletor()
# # set/get the second length parameter
# setLengthParam2(p,0.5)
# getLengthParam2(p)
# # set/get the second length parameter via virtual field
# p$lengthParam2<-0.2
# p$lengthParam2
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"setLengthParam2",
class="FastFieldDeletor",
function(
this,
value,
...
){
if(missing(value)){
throw("No new length parameter value specified!\n");
}
else if ((!is.numeric(value)) | (length(value) != 1 ) ) {
throw("The new value must be a numeric vector of length 1!\n");
}
else {
# First set the scaling factor to NA to force the recalculation:
this$.field.scaling.factor<-NA;
this$.length.param.2<-value;
}
return(invisible(value));
},
private=TRUE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getToleranceMargin
##
###########################################################################/**
#
# @RdocMethod getToleranceMargin
#
# @title "Get the tolerance margin"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A FastFieldDeletor object.}
# \item{...}{Not used.}
# }
#
# \value{
# A numeric vector of length one.
# }
#
# \examples{
# # create a geometric FastFieldDeletor
# p<-FastFieldDeletor()
# # set/get tolerance margin
# setToleranceMargin(p,0.8)
# getToleranceMargin(p)
# # set/get tolerance margin via virtual field
# p$toleranceMargin<-0.75
# p$toleranceMargin
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"getToleranceMargin",
class="FastFieldDeletor",
function(
this,
...
){
this$.tolerance.margin;
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: setToleranceMargin
##
###########################################################################/**
#
# @RdocMethod setToleranceMargin
#
# @title "Set the tolerance margin"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A FastFieldDeletor object.}
# \item{value}{A numeric vector of length one.}
# \item{...}{Not used.}
# }
#
# \value{
# value (invisible).
# }
#
# \examples{
# # create a geometric FastFieldDeletor
# p<-FastFieldDeletor()
# # set/get tolerance margin
# setToleranceMargin(p,0.8)
# getToleranceMargin(p)
# # set/get tolerance margin via virtual field
# p$toleranceMargin<-0.75
# p$toleranceMargin
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"setToleranceMargin",
class="FastFieldDeletor",
function(
this,
value,
...
){
if(missing(value)){
throw("No new length parameter value specified!\n");
}
else if ((!is.numeric(value)) | (length(value) != 1 ) ) {
throw("The new value must be a numeric vector of length 1!\n");
}
else {
# First set the scaling factor to NA to force the recalculation:
this$.field.scaling.factor<-NA;
this$.tolerance.margin<-value;
}
return(invisible(value));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: .getScalingFactor
##
setMethodS3(
".getScalingFactor",
class="FastFieldDeletor",
function(
this,
process,
seq,
...
){
# Check if the length parameters needed for rate scaling are present:
.checkLengthParams(this);
if(is.na(this$.tolerance.max)){
# Get all deletion tolerance parameters for this process:
deletion.tolerance<-c();
for(site in seq$.sites){
if(isAttached(site, process)){
deletion.tolerance<-c(deletion.tolerance, getParameterAtSite(process, site, id="deletion.tolerance")$value);
}
} # for site
# Get the maximal tolerance value:
this$.tolerance.max<-max(deletion.tolerance);
}
d<-max(this$.tolerance.max, this$.tolerance.margin);
this$.d<-d;
# The type specific rate scaling factors:
exp<-expression();
# Geometric:
if(this$.type=="geometric"){
exp<-expression(d * (1 - this$.length.param.1) / (1 - (d * this$.length.param.1)) );
}
# Poisson+1:
else if(this$.type=="poisson"){
exp<-expression(d * exp( - ( this$.length.param.1 * (1 - d ) ) ) );
}
# Negative Binomial + 1:
else if(this$.type=="neg.binomial"){
exp<-expression( d * ( ( (1 - this$.length.param.1) / (1 - (d * this$.length.param.1))) ^ this$.length.param.2) );
}
# Conway-Maxwell Poisson + 1:
else if(this$.type=="compoisson"){
exp<-expression( d * (com.compute.z(lambda=this$.length.param.1,nu=this$.length.param.2 ) / com.compute.z(lambda=(d * this$.length.param.1),nu=this$.length.param.2 )) );
}
return(eval(exp));
},
private=TRUE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: setParameterAtSite.FastFieldDeletor
##
setMethodS3(
"setParameterAtSite",
class="FastFieldDeletor",
function(
this,
site,
id,
value,
...
){
if(value < 0 | value > 1) {
throw("The field deletion model accepts deletion tolerance only from the [0,1] interval!\n");
}
NextMethod();
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
Try the phylosim package in your browser
Any scripts or data that you put into this service are public.
phylosim documentation built on Nov. 22, 2019, 1:07 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
##
## Copyright 2009 Botond Sipos
## See the package description for licensing information.
##
##########################################################################/**
#
# @RdocClass FastFieldDeletor
#
# @title "The FastFieldDeletor class"
#
# \description{
# The \code{\link{DiscreteDeletor}} and \code{\link{ContinuousDeletor}} processes
# generate deletion event objects with rates determined by the general rate of the
# process and the "rate.multiplier" parameter. The probability of rejecting an event
# is determined by the product of the "deletion.tolerance" parameters from the affected sites.
# If the majority of the sites have low deletion tolerance most of the events are rejected, which
# slows down the simulation without performing much events.
#
# The \code{FastFieldDeletor} process scales down the rate and length distribution of the proposed
# events based on the highest insertion tolerance value observed in the root sequence
# (see the package vignette for details), thus making the simulation more efficient.
#
# The available length distributions are (see also the package vignette):
# \itemize{
# \item Geometric (default) - \code{lengthParam1} is \emph{Lambda}.
# \item Poisson+1 - \code{lengthParam1} is \emph{Lambda}.
# \item Conway-Maxwell Poisson+1 - \code{lengthParam1} is \emph{Lambda}, \code{lengthParam2} is \emph{nu}
# \item Negative Binomial+1 - \code{lengthParam1} is \emph{Lambda}, \code{lengthParam2} is \emph{r}
# }
#
# Insertion proceses can insert sites with deletion tolerance higher than the largest
# deletion tolerance observed in the root sequence. The user can specify the largest expected
# tolerance value through the \code{toleranceMargin} virtual field. The process is then scaled by
# max(initial_highest_tolerance, \code{toleranceMargin}).
#
# @classhierarchy
# }
#
# @synopsis
#
# \arguments{
# \item{name}{Object name.}
# \item{type}{Type of the process (see above).}
# \item{length.param.1}{Object name.}
# \item{length.param.2}{Object name.}
# \item{tolerance.margin}{Object name.}
# \item{...}{Additional arguments.}
# }
#
# \section{Fields and Methods}{
# @allmethods
# }
#
# \examples{
# # create a FastFieldDeletor object, default (geometric) type
# # length.param.1 is "lambda"
# p<-FastFieldDeletor(rate=1,length.param.1=0.9, tolerance.margin=0.8)
# # get type
# p$type
# # get object summary
# summary(p)
# # set/get lambda
# p$lengthParam1<-0.8
# p$lengthParam1
# # set/get tolerance margin
# p$toleranceMargin<-0.5
# p$toleranceMargin
# # create a nucleotide sequence, attach process
# s<-NucleotideSequence(length=30,processes=list(list(p)))
# # set state pattern
# s$states<-c("A","A","T","T","G","G","C","C")
# # sample deletion tolerances
# setDeletionTolerance(s,p,sample(seq(from=0,to=0.8,by=0.01),30,replace=TRUE))
# # plot deletion tolerance
# plotParametersAtSites(s,p,"deletion.tolerance")
# # simulate
# sim<-PhyloSim(root.seq=s, phylo=rcoal(2))
# Simulate(sim)
# # show alignment
# sim$alignment
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setConstructorS3(
"FastFieldDeletor",
function(
name="Anonymous",
type="geometric",
length.param.1=NA, # mostly "Lambda"
length.param.2=NA,
tolerance.margin=0, # minimum tolerance value used for scaling.
...
) {
# supported types
ALLOWED.TYPES=c("geometric","poisson","neg.binomial","compoisson");
# Creating a GeneralDeletor Process.
this<-GeneralDeletor(
...
);
# Check if the specified type is valid:
if(length(intersect(ALLOWED.TYPES,type)) != 1){
throw("The specified field model type is invalid!\n");
}
# Extending as FastFieldDeletor:
this<-extend(
this,
"FastFieldDeletor",
.type=type, # field model flavour
.tolerance.margin=NA, # minimum tolerance used for scaling
.tolerance.max=NA, # maximum tolerance obseved at first call of
# getEventAtSites
.d=NA, # is max(.tolerance.max, .tolerance.margin)
.field.scaling.factor=NA,# the precalculated scaling factor
.length.param.1=NA, # mostly "Lambda"
.length.param.2=NA,
.ALLOWED.TYPES=ALLOWED.TYPES # supported types
);
# Set length parameter 1 if not missing:
if(!missing(length.param.1)){
this$lengthParam1<-length.param.1;
}
# Set length parameter 2 if not missing:
if(!missing(length.param.2)){
this$lengthParam2<-length.param.2;
}
# Set tolerance margin:
setToleranceMargin(this, tolerance.margin);
# Using virtual field to clear Id cache:
this$name<-name;
# Set the function proposing deletion lengths:
this$proposeBy<-function(process=NA,seq=NA,pos=NA){
# Check the length parameters:
.checkLengthParams(this);
# Check if this$.d is defined:
if(is.na(this$.d)){
throw("thid$.d is NA! This shouldn't happen!");;
}
d<-this$.d;
# Type specific length sampling expressions:
# Geometric + 1:
if(this$.type == "geometric"){
express<-expression(1 + rgeom(1,prob=( 1 - ( this$.length.param.1 * d) ) ) );
}
# Poisson+1:
else if(this$.type == "poisson"){
express<-expression( 1 + rpois(1,lambda=(this$.length.param.1 * d) ) );
}
# Negative Binomial + 1:
else if(this$.type == "neg.binomial"){
express<-expression(1 + rnbinom(1,this$.length.param.2,prob=( 1 - ( this$.length.param.1 * d))) );
}
# Conway-Maxwell Poisson + 1:
else if(this$.type == "compoisson"){
express<-expression(1 + rcom(1,lambda=( this$.length.param.1 * d), nu = this$.length.param.2));
}
return( round( eval(express) ) );
} # /proposeBy
# Set the function performing the accept/reject step:
this$acceptBy<-function(process=NA,sequence=NA,range=NA){
del.tol<-c();
for(site in sequence$.sites[range]){
# Reject if the range contains a site which is not attached to
# the process:
if(!isAttached(site, process)){
return(FALSE);
}
del.tol<-c(del.tol, getParameterAtSite(process, site, "deletion.tolerance")$value);
}
# Calculate acceptance probability:
accept.prob<-( prod(as.numeric(del.tol)) / this$.d );
# Accept/reject:
return( sample(c(TRUE,FALSE),replace=FALSE,prob=c(accept.prob,(1-accept.prob)),size=1) );
}
return(this);
},
enforceRCC=TRUE
);
##
## Method: .checkLengthParams
##
setMethodS3(
".checkLengthParams",
class="FastFieldDeletor",
function(
this,
...
){
# Check length parameter 1:
if(is.na(this$.length.param.1)){
throw("Length parameter 1 is NA! Cannot generate events or propose lengths!\n");
}
# Check length parameter 2:
if(length(intersect(c("neg.binomial","compoisson"),this$.type)) != 0) {
if(is.na(this$.length.param.2)){
throw("Length parameter 1 is NA! Cannot generate events or propose lengths!\n");
}
}
},
private=TRUE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## 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="FastFieldDeletor",
function(
this,
...
){
wp<-this$writeProtected;
if (wp) {
this$writeProtected<-FALSE;
}
may.fail<-function(this) {
# Check if the type is valid:
if(length(intersect(this$.ALLOWED.TYPES, this$.type)) != 1){
throw("The specified field model type is invalid!\n");
}
if((!is.numeric(this$.tolerance.margin) & !is.na(this$.tolerance.margin))){
throw("Tolerance margin is invalid!\n");
}
if((!is.numeric(this$.tolerance.max) & !is.na(this$.tolerance.max))){
throw(".tolerance.max is invalid!\n");
}
if((!is.numeric(this$.d) & !is.na(this$.d))){
throw(".d is invalid!\n");
} else if(!is.na(this$.tolerance.margin) & !is.na(this$.tolerance.max)) {
if(this$.d != max(this$.tolerance.margin, this$.tolerance.max)){
throw(".d is inconsistent!\n");
}
}
if((!is.numeric(this$.field.scaling.factor) & !is.na(this$.field.scaling.factor))){
throw(".field.scaling.factor is invalid!\n");
}
if((!is.numeric(this$.length.param.1) & !is.na(this$.length.param.1))){
throw("Length parameter 1 is invalid!\n");
}
if((!is.numeric(this$.length.param.2) & !is.na(this$.length.param.2))){
throw("Length parameter 2 is invalid!\n");
}
}
tryCatch(may.fail(this),finally=this$writeProtected<-wp);
NextMethod();
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: summary.FastFieldDeletor
##
###########################################################################/**
#
# @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 an object
# a<-FastFieldDeletor()
# # get a summary
# summary(a)
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"summary",
class="FastFieldDeletor",
function(
object,
...
){
this<-object;
.addSummaryNameId(this);
this$.summary$"Type"<-this$.type;
this$.summary$"Tolerance margin"<-this$.tolerance.margin;
this$.summary$"Length parameter 1"<-this$.length.param.1;
this$.summary$"Length parameter 2"<-this$.length.param.2;
this$.summary$"Scaling factor"<-this$.field.scaling.factor;
NextMethod();
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getEventsAtSite
##
###########################################################################/**
#
# @RdocMethod getEventsAtSite
#
# @title "Generate a deletion event object given the state of the target site"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A FastFieldDeletor object.}
# \item{target.site}{A Site object.}
# \item{...}{Not used.}
# }
#
# \value{
# A list of Event objects.
# }
#
#
# @author
#
# \seealso{
# GeneralDeletor getEventsAtSite.GeneralDeletor
# }
#
#*/###########################################################################
setMethodS3(
"getEventsAtSite",
class="FastFieldDeletor",
function(
this,
target.site,
...
){
if(!exists(x="PSIM_FAST")){
if(missing(target.site)) {
throw("No target site provided!\n");
}
if(!is.Site(target.site)) {
throw("Target site invalid!\n");
}
else if(!is.function(this$.propose.by)) {
throw("proposeBy is not set, cannot propose deletion!\n");
}
else if (!is.function(this$.accept.by)){
throw("acceptBy is not set, cannot generate deletion event deletion!\n");
}
# Complain if sequence has a zero length:
if(target.site$.sequence$.length == 0) {
throw("Sequence has zero length so there is nothing to delete! How did you get here anyway?\n");
}
}
# Clone the event template object:
deletion.event<-clone(this$.event.template);
# Set the target position passed in a temporary field:
deletion.event$.position<-target.site$.position;
# Set the target site:
deletion.event$.site<-target.site;
# Set event name:
deletion.event$.name<-"Deletion";
# Set the generator process:
deletion.event$.process<-this;
# Calculate the field model specific scaling factor if it is not yet calculated:
if(is.na(this$.field.scaling.factor)){
this$.field.scaling.factor<-.getScalingFactor(this,process=this,seq=target.site$.sequence);
}
# Event rate is the product of the general rate, the field model scaling factor and the
# site specific rate multiplier:
deletion.event$.rate<-(this$.rate * (target.site$.processes[[this$.id]]$site.params[["rate.multiplier"]]$value) * this$.field.scaling.factor);
# Set the handler for the deletion event:
deletion.event$.handler<-this$.handler.template;
# Return the event object in a list:
list(deletion.event);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getType
##
###########################################################################/**
#
# @RdocMethod getType
#
# @title "Get the type of a FastFieldDeletor object"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A FastFieldDeletor object.}
# \item{...}{Not used.}
# }
#
# \value{
# A character vector of length one.
# }
#
# \examples{
# # create a FastFieldDeletor, default type (geometric)
# p<-FastFieldDeletor()
# # get type
# getType(p)
# # create a FastFieldDeletor, poisson type
# p<-FastFieldDeletor(type="poisson")
# p$type
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"getType",
class="FastFieldDeletor",
function(
this,
...
){
this$.type;
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: setType
##
###########################################################################/**
#
# @RdocMethod setType
#
# @title "Forbidden action: setting the type of a FastFieldDeletor object"
#
# \description{
# @get "title".
#
# The type can be set only through the \code{type} constructor argument.
# }
#
# @synopsis
#
# \arguments{
# \item{this}{An object.}
# \item{value}{Not used.}
# \item{...}{Not used.}
# }
#
# \value{
# Throws an error.
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"setType",
class="FastFieldDeletor",
function(
this,
value,
...
){
throw("The type of the FastFieldDeletor process cannot be modified. Please set it by the constructor argument.");
},
private=TRUE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getLengthParam1
##
###########################################################################/**
#
# @RdocMethod getLengthParam1
#
# @title "Get the first length parameter"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A FastFieldDeletor object.}
# \item{...}{Not used.}
# }
#
# \value{
# A numeric vector of length one.
# }
#
# \examples{
# # create a geometric FastFieldDeletor
# p<-FastFieldDeletor()
# # set/get the first length parameter
# setLengthParam1(p,0.5)
# getLengthParam1(p)
# # set/get the first length parameter via virtual field
# p$lengthParam1<-0.2
# p$lengthParam1
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"getLengthParam1",
class="FastFieldDeletor",
function(
this,
...
){
this$.length.param.1;
;
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getLengthParam2
##
###########################################################################/**
#
# @RdocMethod getLengthParam2
#
# @title "Get the second length parameter"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A FastFieldDeletor object.}
# \item{...}{Not used.}
# }
#
# \value{
# A numeric vector of length one.
# }
#
# \examples{
# # create a geometric FastFieldDeletor
# p<-FastFieldDeletor()
# # set/get the second length parameter
# setLengthParam2(p,0.5)
# getLengthParam2(p)
# # set/get the second length parameter via virtual field
# p$lengthParam2<-0.2
# p$lengthParam2
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"getLengthParam2",
class="FastFieldDeletor",
function(
this,
...
){
this$.length.param.2;
;
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: setLengthParam1
##
###########################################################################/**
#
# @RdocMethod setLengthParam1
#
# @title "Set the first length parameter"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A FastFieldDeletor object.}
# \item{value}{A numeric vector of length one.}
# \item{...}{Not used.}
# }
#
# \value{
# value (invisible).
# }
#
# \examples{
# # create a geometric FastFieldDeletor
# p<-FastFieldDeletor()
# # set/get the first length parameter
# setLengthParam1(p,0.5)
# getLengthParam1(p)
# # set/get the first length parameter via virtual field
# p$lengthParam1<-0.2
# p$lengthParam1
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"setLengthParam1",
class="FastFieldDeletor",
function(
this,
value,
...
){
if(missing(value)){
throw("No new length parameter value specified!\n");
}
else if ((!is.numeric(value)) | (length(value) != 1 ) ) {
throw("The new value must be a numeric vector of length 1!\n");
}
else {
# First set the scaling factor to NA to force the recalculation:
this$.field.scaling.factor<-NA;
this$.length.param.1<-value;
}
return(invisible(value));
},
private=TRUE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: setLengthParam2
##
###########################################################################/**
#
# @RdocMethod setLengthParam2
#
# @title "Set the second length parameter"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A FastFieldDeletor object.}
# \item{value}{A numeric vector of length one.}
# \item{...}{Not used.}
# }
#
# \value{
# value (invisible).
# }
#
# \examples{
# # create a geometric FastFieldDeletor
# p<-FastFieldDeletor()
# # set/get the second length parameter
# setLengthParam2(p,0.5)
# getLengthParam2(p)
# # set/get the second length parameter via virtual field
# p$lengthParam2<-0.2
# p$lengthParam2
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"setLengthParam2",
class="FastFieldDeletor",
function(
this,
value,
...
){
if(missing(value)){
throw("No new length parameter value specified!\n");
}
else if ((!is.numeric(value)) | (length(value) != 1 ) ) {
throw("The new value must be a numeric vector of length 1!\n");
}
else {
# First set the scaling factor to NA to force the recalculation:
this$.field.scaling.factor<-NA;
this$.length.param.2<-value;
}
return(invisible(value));
},
private=TRUE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getToleranceMargin
##
###########################################################################/**
#
# @RdocMethod getToleranceMargin
#
# @title "Get the tolerance margin"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A FastFieldDeletor object.}
# \item{...}{Not used.}
# }
#
# \value{
# A numeric vector of length one.
# }
#
# \examples{
# # create a geometric FastFieldDeletor
# p<-FastFieldDeletor()
# # set/get tolerance margin
# setToleranceMargin(p,0.8)
# getToleranceMargin(p)
# # set/get tolerance margin via virtual field
# p$toleranceMargin<-0.75
# p$toleranceMargin
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"getToleranceMargin",
class="FastFieldDeletor",
function(
this,
...
){
this$.tolerance.margin;
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: setToleranceMargin
##
###########################################################################/**
#
# @RdocMethod setToleranceMargin
#
# @title "Set the tolerance margin"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A FastFieldDeletor object.}
# \item{value}{A numeric vector of length one.}
# \item{...}{Not used.}
# }
#
# \value{
# value (invisible).
# }
#
# \examples{
# # create a geometric FastFieldDeletor
# p<-FastFieldDeletor()
# # set/get tolerance margin
# setToleranceMargin(p,0.8)
# getToleranceMargin(p)
# # set/get tolerance margin via virtual field
# p$toleranceMargin<-0.75
# p$toleranceMargin
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"setToleranceMargin",
class="FastFieldDeletor",
function(
this,
value,
...
){
if(missing(value)){
throw("No new length parameter value specified!\n");
}
else if ((!is.numeric(value)) | (length(value) != 1 ) ) {
throw("The new value must be a numeric vector of length 1!\n");
}
else {
# First set the scaling factor to NA to force the recalculation:
this$.field.scaling.factor<-NA;
this$.tolerance.margin<-value;
}
return(invisible(value));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: .getScalingFactor
##
setMethodS3(
".getScalingFactor",
class="FastFieldDeletor",
function(
this,
process,
seq,
...
){
# Check if the length parameters needed for rate scaling are present:
.checkLengthParams(this);
if(is.na(this$.tolerance.max)){
# Get all deletion tolerance parameters for this process:
deletion.tolerance<-c();
for(site in seq$.sites){
if(isAttached(site, process)){
deletion.tolerance<-c(deletion.tolerance, getParameterAtSite(process, site, id="deletion.tolerance")$value);
}
} # for site
# Get the maximal tolerance value:
this$.tolerance.max<-max(deletion.tolerance);
}
d<-max(this$.tolerance.max, this$.tolerance.margin);
this$.d<-d;
# The type specific rate scaling factors:
exp<-expression();
# Geometric:
if(this$.type=="geometric"){
exp<-expression(d * (1 - this$.length.param.1) / (1 - (d * this$.length.param.1)) );
}
# Poisson+1:
else if(this$.type=="poisson"){
exp<-expression(d * exp( - ( this$.length.param.1 * (1 - d ) ) ) );
}
# Negative Binomial + 1:
else if(this$.type=="neg.binomial"){
exp<-expression( d * ( ( (1 - this$.length.param.1) / (1 - (d * this$.length.param.1))) ^ this$.length.param.2) );
}
# Conway-Maxwell Poisson + 1:
else if(this$.type=="compoisson"){
exp<-expression( d * (com.compute.z(lambda=this$.length.param.1,nu=this$.length.param.2 ) / com.compute.z(lambda=(d * this$.length.param.1),nu=this$.length.param.2 )) );
}
return(eval(exp));
},
private=TRUE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: setParameterAtSite.FastFieldDeletor
##
setMethodS3(
"setParameterAtSite",
class="FastFieldDeletor",
function(
this,
site,
id,
value,
...
){
if(value < 0 | value > 1) {
throw("The field deletion model accepts deletion tolerance only from the [0,1] interval!\n");
}
NextMethod();
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
Try the phylosim package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.