Nothing
R/Sequence.R
In phylosim: Flexible Simulations of Biological Sequence Evolution
##
## Copyright 2009 Botond Sipos
## See the file ../COPYING for licensing issues.
##
##########################################################################/**
#
# @RdocClass Sequence
#
# @title "The Sequence class"
#
# \description{
#
# This is the class representing a sequence. The backbone of the Sequence objects are
# lists aggregating Site objects. The class has fields for keeping track of cumulative
# site rates, the sum of all active event rates and methods for performing actions
# on a collection of sites (positions).
#
# The Sequence objects have a field specifying an ancestral object, which can be a Sequence
# object (when the object is obtained through clone() ) or the "Root insertion process" object
# (for newly created objects).
#
# @classhierarchy
# }
#
# @synopsis
#
# \arguments{
# \item{name}{The name of the Sequence object.}
# \item{string}{A string containing symbols belonging to the associated Alphabet object.
# It can be used to set the initial states of the aggregated Site objects. It also specifies the length of the sequence}
# \item{length}{The length of the sequence. It cannot be used when 'string' is specified.}
# \item{alphabets}{A list of Alphabet objects to be associated with the Site objects.
# The list is recycled in the case it is shorter than the sequence length.}
# \item{processes}{A list of lists of Process objects to be attached
# (recycled if shorter than sequence length). }
# \item{ancestral.obj}{The ancestral object (Sequence or Process).}
# \item{...}{Not used.}
# }
#
# \section{Fields and Methods}{
# @allmethods
# }
#
# \examples{
# # create a sequence object by
# # providng alphabets, processes and states
# s.one<-Sequence(
# name="Seq",
# string="AATTGGCCTTAAGGCCTTAA",
# alphabets=list(NucleotideAlphabet()),
# processes=list(list(JC69()))
# )
# s.one
# # check if inherits from Sequence
# is.Sequence(s.one)
# # get object summary
# summary(s.one)
# # create a sequence object,
# # specifying length, alphabets
# # and ancestral object
# s<-Sequence(
# name="Seq",
# length=20,
# ancestral.obj=s.one
# )
# # get sequence string
# s$string
# # get the list of site objects
# s$sites
# # get object id
# s$id
# # set and get name
# s$name<-"SeqSeq"
# s$seq
# # get length
# s$length
# # get and set ancestral object
# s$ancestral
# s$ancestral<-Sequence();
# # set alphabets
# setAlphabets(s,list(NucleotideAlphabet()))
# # set states
# # "A"-s in the range 1:10
# setStates(s,"A",1:10)
# # a pattern of "ATGC"-s in the range 11:20
# setStates(s,c("A","T","G","C"),11:20)
# # get states from range 10:12
# getStates(s,10:12)
# # attach a JC69 process to range 1:10
# jc69<-JC69()
# attachProcess(s,jc69,1:10)
# # set the rate multiplier site-process specific parameter for jc69
# setParameterAtSites(s,jc69,"rate.multiplier",2,1:10)
# # get "rate.multiplier" for jc69 from range 1:2
# getParameterAtSites(s, jc69, "rate.multiplier",1:2)
# # attach a GTR process to range 11:20
# gtr<-GTR()
# attachProcess(s,gtr,11:20)
# # set and get rate multipliers for gtr
# setRateMultipliers(s, gtr, c(0,5,1), 11:20)
# getRateMultipliers(s, gtr, 11:20)
# # get processes from range 1:5
# getProcesses(s,1:5)
# # replace the processes with a (GTR, JC69), JC69 pattern
# s$processes<-list(list(GTR(), JC69()), list(JC69()))
# # get processes from range 1:2
# getProcesses(s,1:2)
# # get unique processes
# s$uniqueProcesses
# # get unique alphabets
# s$uniqueAlphabets
# # get symbol frequencies
# getSymbolFreqs(s)
# # get the big rate
# s$bigRate
# # get a vector of total rates from range 1:4
# getTotalRatesFromRange(s,1:4)
# # get a vector of cumulative rates from range 1:4
# getCumulativeRatesFromRange(s,1:4)
# # reset all site states
# clearStates(s)
# s
# # sample states from the equilibrium distributions
# # of the attached substitution processes
# sampleStates(s)
# s
# # clone s
# s.clone<-clone(s)
# # insert a sequence in s.clone after position 2
# insertSequence(s.clone,NucleotideSequence(string="AAAAAA"),2)
# s.clone
# # delete positions 1,2,3 and 10
# deleteSubSequence(s.clone, c(1:3,10))
# s.clone
# # copy positions 7:10 into a new sequence object
# sc<-copySubSequence(s.clone, 7:10)
# sc
# # get events from sc in the range 1:2
# getEvents(sc,1:2)
# }
#
# @author
#
# \seealso{
# Alphabet Site Process Event
# }
#
#*/###########################################################################
setConstructorS3(
"Sequence",
function(
name=NA,
string=NA,
length=NA,
alphabets=NA,
processes=NA,
ancestral.obj=NA,
...
){
# Marking the instance as static by default:
STATIC<-TRUE;
# Extending the PSRoot class:
this<-extend(
PSRoot(),
"Sequence",
.name="Anonymous",
.length=NA,
.sites=list(),
.ancestral.obj=NA,
.cumulative.rates=NA,
.total.rates=NA,
.cumulative.rate.flag=TRUE,
.flagged.sites=integer(0),
.write.protected=FALSE,
.is.sequence=TRUE,
.root.ins=NA
);
# Initializing the variables for length and states:
len<-0;
str<-list();
# Optional argument: name
if(!missing(name)) {
this$name<-name;
STATIC<-FALSE;
}
# The user can specify a sequence
# or the sequence length, but not both.
if (!missing(string) & !missing(length)) {
throw("You can specify the sequence, or the sequence length, but not both!\n");}
else if (!missing(string)){
STATIC<-FALSE;
# An alphabet list must be specified along the sequence!
if(missing(alphabets)) {throw("A list of valid alphabets must be specified when a string is given!\n");}
}
# Deal with the string or length argument:
if (!missing(length)) {
STATIC<-FALSE;
len<-length;
}
else if( !missing(string) ) {
str<-strsplit(string,split="",fixed=TRUE)[[1]];
len<-length(str);
}
this$.length<-len;
root.ins<-NA;
if (!is.Process(Sequence$.root.ins)) {
# Dummy proces used as ancestral object for sites.
this$.root.ins<-Process(name="Root insertion process");
this$.root.ins$comments<-"This is just a dummy process object serving as ancestral for newly created site and sequence objects.";
root.ins<-this$.root.ins;
} else {
root.ins<-Sequence$.root.ins
}
# Site template object:
site.template<-Site(
ancestral=root.ins,
sequence=this
);
# Clone and store the site objects:
if(!STATIC) {
if ( len > 0 ) {
for(position in 1:len) {
this$.sites[[position]]<-clone.Object(site.template);
}
}
}
# Optional argument: ancestral object
if (!missing(ancestral.obj)) {
STATIC<-FALSE;
this$ancestral<-ancestral.obj;
} else {
this$ancestral<-root.ins;
}
# Setting the alphabets:
if(!missing(alphabets)) {
STATIC<-FALSE;
# setAlphabets will check the arguments
setAlphabets(this, alphabets);
}
# Setting the processes:
if (!missing(processes)) {
STATIC<-FALSE;
# setProcesses will take care about the arguments
# and deal with alphabet mismatch.
setProcesses(this,processes)
}
# Initializing these vectors properly is
# importtant for the insertion method!
this$.total.rates<-double(len);
this$.cumulative.rates<-double(len);
if(!STATIC){
# Now we are prepared to set the states:
if (!missing(string)) {
setStates(this, str);
}
# Calculate cumulative rates for the first time, but only if
# states are defined. This is expensive, as total rates are calculated.
if (!missing(string) & (length(str) > 0) ) {
.recalculateCumulativeRates(this);
}
}
return(this);
},
enforceRCC=TRUE
);
##
## Method: is.Sequence
##
###########################################################################/**
#
# @RdocDefault is.Sequence
#
# @title "Check whether an object inherits from the Sequence class"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{An object.}
# \item{...}{Not used.}
# }
#
# \value{
# TRUE or FALSE
# }
#
# \examples{
# # create some objects
# seq<-Sequence(length=10)
# a<-Alphabet()
# # check if they inherit from Sequence
# is.Sequence(seq)
# is.Sequence(a)
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"is.Sequence",
class="default",
function(
this,
...
){
if(!is.PSRoot(this)) {return(FALSE)}
if(!is.null(this$.is.sequence)){return(TRUE)}
if ( inherits(this, "Sequence")) {
this$.is.sequence<-TRUE;
return(TRUE);
} else {
return(FALSE)
}
},
private=FALSE,
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{omit.sites}{Do not check aggregated site objects.}
# \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="Sequence",
function(
this,
omit.sites=FALSE,
...
){
if(this$length != length(this$.sites)) {
throw("Sequence length inconsistency detected");
} else if (length(this$.cumulative.rates) != this$.length) {
throw("Cumulative rates vector length mismatch!\n");
} else if (length(this$.total.rates) != this$.length) {
throw("Total rates vector length mismatch!\n");
} else if (!identical(this$.cumulative.rates, cumsum(this$.total.rates))) {
throw("Cumulative rates vector is not consistent with total rates vector!\n");
}
if(!is.numeric(this$.flagged.sites)) {
throw("Flagged sites vector is not numeric!\n");
} else if (length(this$.flagged.sites) > 0) {
if ( (min(this$.flagged.sites) < 1) | ( max(this$.flagged.sites) > this$.length) ) {
throw("Inconsistency in the flagged sites vector!\n");
}
}
if(!is.character(this$.name)) {
throw("Sequence name is invalid!\n");
} else if(stringLength(this$.name) == 0) {
throw("Sequence name is of length zero!\n");
}
if(!is.Sequence(this$.ancestral.obj) & !is.Process(this$.ancestral.obj)) {
throw("The ancestral object is invalid!\n");
}
if(!is.logical(this$.cumulative.rate.flag)) {
throw("Cumulative rate flag is not logical!\n");
}
# Calling consistency check on sites.
# This will be painfully slow!
if(!omit.sites){
for(site in this$.sites){
checkConsistency(site);
}
}
return(invisible(TRUE));
},
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 Sequence object"
#
# \description{
# @get "title".
# The identifier is the concatenation of the object name and the object hash code as returned
# by hashCode().
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{...}{Not used.}
# }
#
# \value{
# A charcter vector of length one.
# }
#
# \examples{
# # create a Sequence object.
# s<-Sequence(length=5)
# # get id
# getId(s)
# # get id via virtual field
# s$id
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"getId",
class="Sequence",
function(
this,
...
){
this.class<-class(this)[1];
id<-paste(this.class,this$.name,hashCode(this),sep=":");
},
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 Sequence 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="Sequence",
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: getName
##
###########################################################################/**
#
# @RdocMethod getName
#
# @title "Get the name of a Sequence object"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{...}{Not used.}
# }
#
# \value{
# A character vector of length one.
# }
#
# \examples{
# # create object
# s<-Sequence(length=10);
# # get object name
# getName(s)
# # get name via virtual field
# s$name
#
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"getName",
class="Sequence",
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 Sequence object"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{new.name}{A character vector of length one.}
# \item{...}{Not used.}
# }
#
# \value{
# The new object name (invisible).
# }
#
# \examples{
# # create a Sequence object
# s<-Sequence(name="MySeq",length=4)
# # get sequence name
# s$name
# # rename object
# setName(s,"SO")
# s$name
# # rename via virtual field
# s$name<-"SeqSeq"
# s$name
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"setName",
class="Sequence",
function(
this,
new.name,
...
){
.checkWriteProtection(this);
this$.name<-as.character(new.name);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getString
##
###########################################################################/**
#
# @RdocMethod getString
#
# @title "Get the string representation of a Sequence object"
#
# \description{
# @get "title".
# The string representation is the concatenation of the states of the
# aggregated Site object. Undefined states (NA-s) are represented by question marks.
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{...}{Not used.}
# }
#
# \value{
# A character vector of length one.
# }
#
# \examples{
# # create object
# s<-Sequence(length=10)
# # get character representation
# getString(s) # a bunch of '?'-s
# # get string reperesentation via virtual field
# s$string
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"getString",
class="Sequence",
function(
this,
...
){
str<-character();
for(site in this$.sites){
if(is.na(site$.state)){
if(is.Alphabet(site$.alphabet)){
str<-c(str,rep("?",site$.alphabet$.symbolLength));
}
else {
str<-c(str,"?");
}
}
else {
str<-c(str,site$.state);
}
}
return(paste(str,collapse=""));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: setString
##
###########################################################################/**
#
# @RdocMethod setString
#
# @title "Forbidden action: setting the string representation of a Sequence 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(
"setString",
class="Sequence",
function(
this,
value,
...
){
virtualAssignmentForbidden(this);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getLength
##
###########################################################################/**
#
# @RdocMethod getLength
#
# @title "Get the number of Site objects aggregated in a Sequence object"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{...}{Not used.}
# }
#
# \value{
# An integer vector of length one.
# }
#
# \examples{
# # create a Sequence object
# s<-Sequence(length=5)
# # get sequence length
# getLength(s)
# # get length via virtual field
# s$length
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"getLength",
class="Sequence",
function(
this,
...
){
this$.length;
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: setLength
##
###########################################################################/**
#
# @RdocMethod setLength
#
# @title "Forbidden action: setting the length of a Sequence object"
#
# \description{
# @get "title".
# The length of the Sequence objects can be specified when the object is constructed,
# or modified later by the "insertSequence" and "deleteSubSequence" methods.
# }
#
# @synopsis
#
# \arguments{
# \item{this}{An object.}
# \item{value}{Not used.}
# \item{...}{Not used.}
# }
#
# \value{
# Throws an error.
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"setLength",
class="Sequence",
function(
this,
value,
...
){
virtualAssignmentForbidden(this);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getSites
##
###########################################################################/**
#
# @RdocMethod getSites
#
# @title "Get the list of the Site object aggregated in a Sequence object"
#
# \description{
# @get "title".
# Warning: there is no setSites method!
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{...}{Not used.}
# }
#
# \value{
# A list of Site objects.
# }
#
# \examples{
# # create a sequence object
# s<-Sequence(alphabets=list(NucleotideAlphabet()),string="AATTGCCC")
# # get the list of aggregated Site objects
# getSites(s)
# # get Site objects via virtual field
# s$sites
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"getSites",
class="Sequence",
function(
this,
...
){
this$.sites;
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getStates
##
###########################################################################/**
#
# @RdocMethod getStates
#
# @title "Get the states of a set of Site objects aggregated by a Sequence object"
#
# \description{
# @get "title".
#
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{index}{An integer vector specifying a set of positions.
# It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# A character vector.
# }
#
# \examples{
# # create a sequence object
# s<-Sequence(alphabets=list(NucleotideAlphabet()),string="AATTGCCCCCTTGG")
# # get all Site states
# getStates(s)
# # get the states for a collection of sites
# getStates(s,c(1:3,5,8))
# # get states via virtual field
# s$states
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"getStates",
class="Sequence",
function(
this,
index,
...
){
if (missing(index)) {
index<-seq(along.with=this$.sites);
} else {
index<-.checkIndexSanity(this, index);
}
lapply(
this$.sites[index],
function(site) {
site$.state;
}
);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: setStates
##
###########################################################################/**
#
# @RdocMethod setStates
#
# @title "Set the states for a set of Site objects aggregated by a Sequence object"
#
# \description{
# @get "title".
# The value vector is recycled, which is useful when creating repeated patterns.
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{value}{A character vector containg the states (recycled if shorter than the index vector). The states must be compatible with the corresponding Alphabet object.}
# \item{index}{An integer vector specifying a set of positions. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# The Sequence object (invisible).
# }
#
# \examples{
# # create a sequence object of length 10
# s<-Sequence(alphabets=list(NucleotideAlphabet()),length=10)
# # set the states in some ranges
# setStates(s,c("A","T","A"),index=1:5)
# setStates(s,c("G","C"),index=6:10)
# # display sequence
# s
# # set the states for the whole Sequence object
# setStates(s,c("A","T","T","A"))
# s
# # set states via virtual field
# s$states<-c("A","T")
# s
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"setStates",
class="Sequence",
function(
this,
value,
index,
...
){
.checkWriteProtection(this);
if(missing(value)) {
throw("No new values specified!\n");
}
else if (missing(index)) {
index<-seq(along.with=this$.sites);
} else {
index<-.checkIndexSanity(this, index);
}
# Recycling value vector by using rep().
if (length(value) < length(this$.sites) ) {
value<-rep(as.character(value),length.out=length(index))
}
for (i in 1:length(index)) {
this$.sites[[ index[[i]] ]]$state<-value[i];
}
# Flagging the changed sites:
this$.cumulative.rate.flag<-TRUE;
this$.flagged.sites<-c(this$.flagged.sites, index);
.recalculateCumulativeRates(this);
invisible(this);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getAlphabets
##
###########################################################################/**
#
# @RdocMethod getAlphabets
#
# @title "Get the list of the Alphabet objects attached to the Site objects aggregated by a Sequence object"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{index}{An integer vector specifying a set of positions. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# A list of Alphabet objects.
# }
#
# \examples{
# # create a Sequence object with NucleotideAlphabet
# #and BinaryAlphabet objects attached
# s<-Sequence(alphabets=list(NucleotideAlphabet(),BinaryAlphabet()),length=5)
# setStates(s,c("A","0"))
# # get the list of attached Alphabet objects
# getAlphabets(s)
# # get Alphabets from a range
# getAlphabets(s,c(2:3,5))
# # get alphabets via virtual field
# s$alphabets
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"getAlphabets",
class="Sequence",
function(
this,
index,
...
){
if (missing(index)) {
index<-seq(along.with=this$.sites);
} else {
index<-.checkIndexSanity(this, index);
}
lapply(
this$.sites[index],
function(site) {
site$alphabet;
}
);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: setAlphabets
##
###########################################################################/**
#
# @RdocMethod setAlphabets
#
# @title "Assotiate Alphabet objects to a set of Site objects aggregated by a Sequence object"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{value}{A list of Alphabet objects, recycled if shorter than the index vector.}
# \item{index}{An integer vector specifying a set of positions. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# The Sequence object (invisible).
# }
#
# \examples{
# # create a Sequence object
# s<-Sequence(length=10)
# # set the alphabets for range 1:5
# setAlphabets(s,list(NucleotideAlphabet(),BinaryAlphabet()),1:5)
# # set the alphabets for range 6:10
# setAlphabets(s,list(AminoAcidAlphabet()),6:10)
# # get the list of attached Alphabet objects
# getAlphabets(s)
# # get Alphabets from a range
# getAlphabets(s,c(2:3,5))
# # set alphabets via virtual field
# s$alphabets<-list(BinaryAlphabet(),NucleotideAlphabet())
# # get alphabets via virtual field
# s$alphabets
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"setAlphabets",
class="Sequence",
function(
this,
value,
index,
...
){
.checkWriteProtection(this);
if(missing(value)) {
throw("No new values specified!\n");
}
if(!is.list(value)) {
throw("The value parameter must be a list!\n");
} else {
for(a in value) {
if(!is.Alphabet(a)) {
throw("The value parameter must be a list of valid alphabet objects!\n");
}
}
}
if (missing(index)) {
index<-seq(along.with=this$.sites);
} else {
index<-.checkIndexSanity(this, index);
}
# Recycling value vector. rep() cannot be used here,
# because we loose the object references!
value.counter<-1;
for (i in index) {
if(value.counter > length(value)) {
value.counter<-1;
}
setAlphabet(this$.sites[[i]], value[[value.counter]]);
value.counter<-(value.counter + 1);
}
invisible(this);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getUniqueAlphabets
##
###########################################################################/**
#
# @RdocMethod getUniqueAlphabets
#
# @title "Get the list of unique Alphabet objects associated to Site objects aggaregated by a Sequence object"
#
# \description{
# @get "title".
# The returned list contains unique instances of the Alphabet class. The symbol sets are not compared, so
# two instances of the same class are considered to be different.
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{...}{Not used.}
# }
#
# \value{
# A list of Alphabet objects.
# }
#
# \examples{
# # create a Sequence object with some Alphabet objects attached
# s<-Sequence(
# alphabets=list(NucleotideAlphabet(),
# BinaryAlphabet(),
# NucleotideAlphabet()),
# length=10
# )
# # get the list of attached alphabets
# s$alphabets
# # get the unique list of attahced Alphabet objects
# getUniqueAlphabets(s)
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"getUniqueAlphabets",
class="Sequence",
function(
this,
...
){
tmp<-list();
lapply(
this$.sites,
function(site) {
tmp<<-c(tmp,list(site$.alphabet))
}
);
return(unique(tmp));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: setUniqueAlphabets
##
###########################################################################/**
#
# @RdocMethod setUniqueAlphabets
#
# @title "Forbidden action: setting the list of unique Alphabet objects attached to the Site object aggregated by a Sequence 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(
"setUniqueAlphabets",
class="Sequence",
function(
this,
value,
...
){
virtualAssignmentForbidden(this);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: attachProcess
##
###########################################################################/**
#
# @RdocMethod attachProcess
#
# @title "Attach a Process object to a set of Site objects aggregated by a Sequence object"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{process}{A Process object.}
# \item{index}{An integer vector specifying a set of positions. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# The Sequence object (invisible).
# }
#
# \examples{
# # create a Sequence object of length 6
# s<-Sequence(length=10,alphabets=list(NucleotideAlphabet()))
# # attach a JC69 substitution process
# attachProcess(s,JC69())
# # get the list of attached processes
# s$processes
# # attach the GTR substitution process to range 3:6
# attachProcess(s,GTR(),3:6)
# # get the list of attached processes
# s$processes
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"attachProcess",
class="Sequence",
function(
this,
process,
index,
...
){
.checkWriteProtection(this);
if(!exists(x="PSIM_FAST")){
if(!is.Process(process)){
throw("Process object invalid!\n");
}
}
if (missing(index)) {
index<-seq(along.with=this$.sites);
} else {
index<-.checkIndexSanity(this, index);
}
for(i in index){
attachProcess(this$.sites[[i]],process);
}
return(invisible(this));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: detachProcess
##
###########################################################################/**
#
# @RdocMethod detachProcess
#
# @title "Detach a Process object from a set of Site objects aggregated by a Sequence object"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{process}{A Process object.}
# \item{index}{An integer vector specifying a set of positions. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# The Sequence object (invisible).
# }
#
# \examples{
# # create a sequence object with two processes attached
# p<-JC69()
# s<-Sequence(length=4,alphabets=list(NucleotideAlphabet()),processes=list(list(p,K80())))
# # get the list of attached processes
# s$processes
# # detach JC69 from range c(1,4)
# detachProcess(s,p,c(1,4))
# s$processes
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"detachProcess",
class="Sequence",
function(
this,
process,
index,
...
){
.checkWriteProtection(this);
if(!exists(x="PSIM_FAST")){
if(!is.Process(process)){
throw("Process object invalid!\n");
}
}
if (missing(index)) {
index<-seq(along.with=this$.sites);
} else {
index<-.checkIndexSanity(this, index);
}
lapply(
this$.sites[index],
function(site) {
detachProcess(site,process);
}
);
return(invisible(this));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getProcesses
##
###########################################################################/**
#
# @RdocMethod getProcesses
#
# @title "Get the Process objects attached to the Site objects aggregated by a Sequence object"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{index}{An integer vector specifying a set of positions. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# A list of lists of Process objects.
# }
#
# \examples{
# # create a sequence object with some processes attached
# s<-Sequence(
# length=4,
# alphabets=list(NucleotideAlphabet()),
# processes=list(list(JC69(),K80()),list(GTR()))
# )
# # get the list of lists of attached processes from positions 1 and 3
# getProcesses(s,c(1,3))
# # get processes via virtual field
# s$processes
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"getProcesses",
class="Sequence",
function(
this,
index,
...
){
if (missing(index)) {
index<-seq(along.with=this$.sites);
} else {
index<-.checkIndexSanity(this, index);
}
lapply(
this$.sites[index],
function(site) {
site$processes;
}
);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getUniqueProcesses
##
###########################################################################/**
#
# @RdocMethod getUniqueProcesses
#
# @title "Get the list of unique Process instances attached to the Site objects aggregated by a Sequence object"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{...}{Not used.}
# }
#
# \value{
# A list of Process objects.
# }
#
# \examples{
# # create a sequence object and attach processes
# p<-JC69()
# s<-Sequence(
# length=4,
# alphabets=list(NucleotideAlphabet()),
# processes=list(list(p,K80()),list(p))
# )
# # get the unique list of attached Process instances
# getUniqueProcesses(s)
# # via virtual field
# s$uniqueProcesses
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"getUniqueProcesses",
class="Sequence",
function(
this,
...
){
tmp<-list();
lapply(
this$.sites,
function(site) {
tmp<<-c(tmp,site$processes)
}
);
return(unique(tmp));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: setUniqueProcesses
##
###########################################################################/**
#
# @RdocMethod setUniqueProcesses
#
# @title "Forbidden action: setting the list of unique Process instances attached to the sites of a Sequence 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(
"setUniqueProcesses",
class="Sequence",
function(
this,
value,
...
){
virtualAssignmentForbidden(this);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: setProcesses
##
###########################################################################/**
#
# @RdocMethod setProcesses
#
# @title "Specify a set of Process objects to be attached to a set of Site objects aggregated by a Sequence object"
#
# \description{
# @get "title".
# The Process objects in a given inner list correspond to the set of processes to be attached to one Site object.
# Process objects already attached to a given Site are skipped. Attached processes which are not memebers of the list
# are detached, so specifying an empty list will detach all processes.
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{value}{A list of list of Process objects, recycled if shorter than the index vector.}
# \item{index}{An integer vector specifying a set of positions. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# The Sequence object (invisible).
# }
#
# \examples{
# # create a sequence of length 5
# s<-Sequence(length=5,alphabets=list(NucleotideAlphabet()));
# # set a pattern of processes
# setProcesses(s,list(list(JC69(),K81())))
# # get attached processes
# s$processes
# # detach all processes from range 1:3
# setProcesses(s,list(list()),1:3)
# s$processes
# # detach all processes via virtual field
# s$processes<-list(list())
# # create a process pattern in the full sequence via virtual field
# s$processes<-list(list(JC69()),list(GTR(),K80()))
# s$processes
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"setProcesses",
class="Sequence",
function(
this,
value,
index,
...
){
.checkWriteProtection(this);
if(!is.list(value)) {
throw("The value parameter must be a list!\n");
} else {
lapply(
value,
function(element) {
if(!is.list(element)){
throw("The value parameter must be a list of lists containing process objects!\n");
}
}
);
}
if (missing(index)) {
index<-seq(along.with=this$.sites);
} else {
index<-.checkIndexSanity(this, index);
}
value.counter<-1;
# Recycling value vector. rep() cannot be used here,
# because we loose the object references!
for (i in index) {
if(value.counter > length(value)) {
value.counter<-1;
}
setProcesses(this$.sites[[i]], value[[value.counter]]);
value.counter<-(value.counter + 1);
}
invisible(this);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: setParameterAtSites
##
###########################################################################/**
#
# @RdocMethod setParameterAtSites
#
# @title "Set the values of a site-process specific paramater for a process and a collection of Site objects aggregated by a Sequence object"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{process}{A valid Process object.}
# \item{id}{The identifier of the site-process specific parameter.}
# \item{value}{A vector containing the new values of the site-process specific parameter, recycled if shorter than the index vector. It should be consistent with the type of the parameter.}
# \item{index}{An integer vector specifying a set of positions. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# The Sequence object (invisible).
# }
#
# \examples{
# # create a sequence, attach a process
# p<-K80()
# s<-Sequence(length=6,alphabets=list(NucleotideAlphabet()),processes=list(list(p)))
# # set a new pattern of rate multipliers in the range 1:3,
# # the default value is 1.0 by the way
# setParameterAtSites(s,p,"rate.multiplier",c(2,3),1:3)
# # get rate multipliers
# getParameterAtSites(s,p,"rate.multiplier")
# # set a new value for the whole sequence
# setParameterAtSites(s,p,"rate.multiplier",0.5)
# # get rate multipliers
# getParameterAtSites(s,p,"rate.multiplier")
# }
#
# @author
#
# \seealso{
# Site Process @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"setParameterAtSites",
class="Sequence",
function(
this,
process,
id,
value,
index,
...
){
.checkWriteProtection(this);
if(!exists(x="PSIM_FAST")){
if(missing(process)) {
throw("No process given!\n");
}
else if(!is.Process(process)){
throw("Process object invalid!\n");}
else if (missing(id)) {
throw("No site-process specific parameter id given!\n");
} else if (!is.character(id)) {
throw("Parameter id must be character!\n");
} else if (missing(value)){
throw("No new value given!\n");
}
}
if (missing(index)) {
index<-seq(along.with=this$.sites);
} else {
index<-.checkIndexSanity(this, index);
}
if(length(value) == 1) {
lapply(
this$.sites[index],
function(site){
setParameterAtSite(process,site,id,value);
}
);
} else {
counter<-1;
lapply(
this$.sites[index],
function(site){
if( counter > length(value) ){
counter<<-1;
}
setParameterAtSite(process,site,id,value[[counter]]);
counter<<-(counter+1);
}
);
}
return(invisible(this));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: setRateMultipliers
##
###########################################################################/**
#
# @RdocMethod setRateMultipliers
#
# @title "Set the values of the rate multiplier parameters for a given Process object and a collection of Site object aggregated by a Sequence object"
#
# \description{
# @get "title".
# This method just calls \code{setParameterAtSites(this=this,process=process,id="rate.multiplier",value=value,index=index)} See setParameterAtSites.Sequence for details.
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{process}{A valid Process object.}
# \item{value}{A numeric vector containing the new values of the site-process specific parameter, recycled if shorter than the index vector.}
# \item{index}{An integer vector specifying a set of positions. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# The Sequence object (invisible).
# }
#
#
# @author
#
# \seealso{
# setParameterAtSites.Sequence
# }
#
#*/###########################################################################
setMethodS3(
"setRateMultipliers",
class="Sequence",
function(
this,
process,
value,
index,
...
){
if(missing(process)){
throw("No process given!\n");
}
else if(missing(value)){
throw("No value provided!\n");
}
#else if(!is.GeneralSubstitution(process)){
# throw("The specified process is not a substitution process!\n");
#}
setParameterAtSites(this=this,process=process,id="rate.multiplier",value=value,index=index);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getRateMultipliers
##
###########################################################################/**
#
# @RdocMethod getRateMultipliers
#
# @title "Get the values of the rate multiplier parameters for a given Process object and a collection of Site object aggregated by a Sequence object"
#
# \description{
# @get "title".
# This method just calls \code{getParameterAtSites(this=this,process=process,id="rate.multiplier",index=index)} See getParameterAtSites.Sequence for details.
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{process}{A valid Process object.}
# \item{index}{An integer vector specifying a set of positions. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# A numeric vector with the current values of the rate multiplier in the specified range.
# }
#
#
# @author
#
# \seealso{
# setParameterAtSites.Sequence
# }
#
#*/###########################################################################
setMethodS3(
"getRateMultipliers",
class="Sequence",
function(
this,
process,
index,
...
){
if(missing(process)){
throw("No process given!\n");
}
rm<-getParameterAtSites(this=this,process=process,id="rate.multiplier",index=index);
return(as.numeric(lapply(rm,function(param){param$value})));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getParameterAtSites
##
###########################################################################/**
#
# @RdocMethod getParameterAtSites
#
# @title "Get the values of a site-process specific paramater for a process and a collection of Site objects aggregated by a Sequence object"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{process}{A valid Process object.}
# \item{id}{The identifier of the site-process specific parameter.}
# \item{index}{An integer vector specifying a set of positions. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# A list of site specific paramters. A site specific paramter is a list storing the id, the name, the value
# and the type of the parameter.
# }
#
# \examples{
# # create a sequence, attach a process
# p<-K80()
# s<-Sequence(length=6,alphabets=list(NucleotideAlphabet()),processes=list(list(p)))
# # set a new pattern of rate multipliers in the
# # range 1:3, the default value is 1.0 by the way
# setParameterAtSites(s,p,"rate.multiplier",c(2,3),1:3)
# # get rate multipliers
# getParameterAtSites(s,p,"rate.multiplier")
# # set a new value for the whole sequence
# setParameterAtSites(s,p,"rate.multiplier",0.5)
# # get rate multipliers
# getParameterAtSites(s,p,"rate.multiplier")
# }
#
# @author
#
# \seealso{
# Site Process @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"getParameterAtSites",
class="Sequence",
function(
this,
process,
id,
index,
...
){
if(!exists(x="PSIM_FAST")){
if(missing(process)) {
throw("No process given!\n");
}
else if(!is.Process(process)){
throw("Process object invalid!\n");}
else if (missing(id)) {
throw("No site-process specific parameter id given!\n");
} else if (!is.character(id)) {
throw("Parameter id must be character!\n");
}
}
if (missing(index)) {
index<-seq(along.with=this$.sites);
} else {
index<-.checkIndexSanity(this, index);
}
lapply(
this$.sites[index],
function(site){
getParameterAtSite(process,site,id);
}
);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getEvents
##
###########################################################################/**
#
# @RdocMethod getEvents
#
# @title "Get the list of active Event objects for a set of Site objects aggregated by a Sequence object"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{index}{An integer vector specifying a set of positions. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# A list of Event objects.
# }
#
# \examples{
# # create a sequence with a process attached
# s<-Sequence(
# string="ATGC",
# alphabets=list(NucleotideAlphabet()),
# processes=list(list(JC69()))
# )
# # get the active events from range 1:3
# getEvents(s,1:3)
# # get all active events via virtual field
# s$events
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"getEvents",
class="Sequence",
function(
this,
index,
...
){
if (missing(index)) {
index<-seq(along.with=this$.sites);
} else {
index<-.checkIndexSanity(this, index);
}
tmp<-list();
for (i in index){
# Setting the .positions field for then Events.
this$.sites[[i]]$.position<-i;
tmp<-c(tmp, getEvents(this$.sites[[i]]));
# Deleting the .position field;
this$.sites[[i]]$.position<-NULL;
}
tmp;
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getTotalRatesFromRange
##
###########################################################################/**
#
# @RdocMethod getTotalRatesFromRange
#
# @title "Get the vector of total site rates for a collection of Site objects aggregated by a Sequence object"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{index}{An integer vector specifying a set of positions. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# A numeric vector.
# }
#
# \examples{
# # create a sequence with some processes attached
# s<-Sequence(
# string="ATGC",
# alphabets=list(NucleotideAlphabet()),
# processes=list(list(JC69()),list(JC69(),GTR()))
# )
# # get total rates for positions 1 and 3
# getTotalRatesFromRange(s,c(1,3))
# # get all total rates via virtual field
# s$totalRates # via the "getTotalRates.Sequence" method
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"getTotalRatesFromRange",
class="Sequence",
function(
this,
index,
...
){
if (missing(index)) {
index<-seq(along.with=this$.sites);
} else {
index<-.checkIndexSanity(this, index);
}
if (this$.cumulative.rate.flag){
.recalculateCumulativeRates(this);
}
this$.total.rates[index];
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getTotalRates
##
###########################################################################/**
#
# @RdocMethod getTotalRates
#
# @title "Get the total site rates from a Sequence object"
#
# \description{
# @get "title".
# This method simply calls \code{getTotalRatesFromRange(this)}.
# See \code{getTotalRatesFromRange.Sequence} for more details.
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{...}{Not used.}
# }
#
# \value{
# A numeric vector containing the total site rates.
# }
#
# @author
#
# \seealso{
# getTotalRatesFromRange.Sequence
# }
#
#*/###########################################################################
setMethodS3(
"getTotalRates",
class="Sequence",
function(
this,
...
){
getTotalRatesFromRange(this);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: setTotalRates
##
###########################################################################/**
#
# @RdocMethod setTotalRates
#
# @title "Forbidden action: setting the list of total site rates for a Sequence 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(
"setTotalRates",
class="Sequence",
function(
this,
value,
...
){
virtualAssignmentForbidden(this);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getAncestral
##
###########################################################################/**
#
# @RdocMethod getAncestral
#
# @title "Get the ancestral object of a Sequence object"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{...}{Not used.}
# }
#
# \value{
# A Sequence object or a Process object.
# }
#
# \examples{
# # create a sequence object
# s<-Sequence(length=4)
# # get ancestral object
# getAncestral(s) # newly created sequences have the "Root insertion process" as ancestral
# # clone sequence
# cs<-clone(s)
# # get ancestral object id via virtual field
# cs$ancestral$id
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"getAncestral",
class="Sequence",
function(
this,
...
){
this$.ancestral.obj;
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getCumulativeRatesFromRange
##
###########################################################################/**
#
# @RdocMethod getCumulativeRatesFromRange
#
# @title "Get the cumulative site rates for a collection of Site objects aggregated by a Sequence object"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{index}{An integer vector specifying a set of positions. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# A numeric vector.
# }
#
# \examples{
# # create a sequence with some processes attached
# s<-Sequence(
# string="ATGC",
# alphabets=list(NucleotideAlphabet()),
# processes=list(list(JC69()),list(JC69(),GTR()))
# )
# # get cumulative rates for positions 1 and 3
# getCumulativeRatesFromRange(s,c(1,3))
# # get all cumulative rates via virtual field
# s$cumulativeRates # via the "getCumulativeRates.Sequence" method
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"getCumulativeRatesFromRange",
class="Sequence",
function(
this,
index,
...
){
if (missing(index)) {
index<-seq(along.with=this$.sites);
} else {
index<-.checkIndexSanity(this, index);
}
if (this$.cumulative.rate.flag){
.recalculateCumulativeRates(this);
}
this$.cumulative.rates[index];
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
setMethodS3(
".getCumulativeRatesFast",
class="Sequence",
function(
this,
...
){
if (this$.cumulative.rate.flag){
.recalculateCumulativeRates(this);
}
this$.cumulative.rates;
},
private=TRUE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: .recalculateCumulativeRates
##
setMethodS3(
".recalculateCumulativeRates",
class="Sequence",
function(
this,
target.site,
...
){
length<-this$.length;
if(length == 0){
return();
}
total.rates<-this$.total.rates;
sites<-this$.sites;
flagged.sites<-this$.flagged.sites;
if( length(flagged.sites) == 0 ) {
# Fresh start:
total.rates<-as.numeric(lapply(sites,getTotalRate));
} else {
# We have some flagged sites, recalculate just their total rates:
total.rates[flagged.sites]<-as.numeric(lapply(sites[flagged.sites],getTotalRate));
}
this$.total.rates<-total.rates;
this$.cumulative.rates<-cumsum(total.rates);
this$.flagged.sites<-integer(0);
this$.cumulative.rate.flag<-FALSE;
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getCumulativeRates
##
###########################################################################/**
#
# @RdocMethod getCumulativeRates
#
# @title "Get the total site rates from a Sequence object"
#
# \description{
# @get "title".
# This method simply calls \code{getCumulativeRatesFromRange(this)}.
# See \code{getCumulativeRates.Sequence} for more details.
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{...}{Not used.}
# }
#
# \value{
# A numeric vector containing the total site rates.
# }
#
# @author
#
# \seealso{
# getCumulativeRatesFromRange.Sequence
# }
#
#*/###########################################################################
setMethodS3(
"getCumulativeRates",
class="Sequence",
function(
this,
...
){
getCumulativeRatesFromRange(this);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: setCumulativeRates
##
###########################################################################/**
#
# @RdocMethod setCumulativeRates
#
# @title "Forbidden action: setting the cumulative rates for the sites aggregated by a Sequence 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(
"setCumulativeRates",
class="Sequence",
function(
this,
value,
...
){
virtualAssignmentForbidden(this);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getBigRate
##
###########################################################################/**
#
# @RdocMethod getBigRate
#
# @title "Get the sum of all active event rates from a Sequence object"
#
# \description{
# @get "title".
# The sum of active event rates depends on all Site object states and on the attached Process objects.
# It basically returns the last element of the cumulative site rates vector.
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{...}{Not used.}
# }
#
# \value{
# A numeric vector of length one.
# }
#
# \examples{
# # create a nucleotide sequence attach a process
# s<-NucleotideSequence(length=5);
# s$processes<-list(list(JC69()))
# # get the sum of active event rates
# getBigRate(s) # returns NA because site states are undefined
# # set site states
# s$states<-c("A","T")
# # get big rate via virtual field
# s$bigRate
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"getBigRate",
class="Sequence",
function(
this,
...
){
if (length(this$.sites) > 0) {
if(this$.cumulative.rate.flag) {
.recalculateCumulativeRates(this);
}
return(this$.cumulative.rates[this$.length]);
} else {
return(NA);
}
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: setBigRate
##
###########################################################################/**
#
# @RdocMethod setBigRate
#
# @title "Forbidden action: setting the sum of total active event rates for a Sequence 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(
"setBigRate",
class="Sequence",
function(
this,
value,
...
){
virtualAssignmentForbidden(this);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: setAncestral
##
###########################################################################/**
#
# @RdocMethod setAncestral
#
# @title "Set the ancestral object of a Sequence object"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{value}{A Sequence or a Process object.}
# \item{...}{Not used.}
# }
#
# \value{
# The new ancestral object (invisible).
# }
#
# \examples{
# # create a nucleotide sequence and a process object
# s<-NucleotideSequence(string="AGCT")
# p<-Process(name="MyProcess")
# # set the p as the ancestral of s
# setAncestral(s,p)
# s$ancestral
# # clone s
# cs<-clone(s)
# # set cs as ancestral of s via virtual field
# s$ancestral<-cs
# # get ancestral ids
# s$ancestral$id
# cs$ancestral$id
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"setAncestral",
class="Sequence",
function(
this,
value,
...
){
.checkWriteProtection(this);
if (!is.Sequence(value) & ! is.Process(value)) {
throw("Ancestral object must be a sequence or a process!\n");
} else {
this$.ancestral.obj<-value;
}
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: clone.Sequence
##
###########################################################################/**
#
# @RdocMethod clone
#
# @title "Clone a Sequence object"
#
# \description{
# @get "title".
#
# The cloning of Sequence objects involves the cloning of all aggregated Site objects. Because of that the
# cloning of long sequences is quite expensive.
# The cloned Site objects have the orginal Site objects as ancestral.
# The new Sequence objects has the original object as ancestral.
#
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{...}{Not used.}
# }
#
# \value{
# A Sequence object.
# }
#
# \examples{
# # cretate a nucleotide sequence
# s<-NucleotideSequence(string="ATG")
# # clone the sequence
# cs<-clone(s)
# # get some properties
# equals(s,s)
# equals(s,cs)
# cs$ancestral
# cs$sites[[1]]$ancestral
# }
#
# @author
#
# \seealso{
# Sequence clone.Object
# }
#
#*/###########################################################################
setMethodS3(
"clone",
class="Sequence",
function(
this,
...
){
# Cloning the whole sequence object:
that<-clone.Object(this)
# Disabling write protection:
if(that$writeProtected) {
that$writeProtected<-FALSE
}
# Setting the ancestral sequence:
that$.ancestral.obj<-this
# Resetting comments:
that$.comments<-list()
# Cloning sites:
clone.sites<-that$.sites
if(this$.length > 0) {
for (i in 1:this$.length) {
site<-this$.sites[[i]]
clone<-clone.Object(site)
clone$.ancestral<-site
clone$.sequence<-that
clone.sites[[i]]<-clone
}
}
that$.sites<-clone.sites
# Setting the name:
that$name<-paste("clone of",this$.name)
return(that)
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getWriteProtected
##
###########################################################################/**
#
# @RdocMethod getWriteProtected
#
# @title "Check if the object is write protected"
#
# \description{
# @get "title".
# Write protected objects cannot be modified through get/set methods and virtual fields.
# }
#
# @synopsis
#
# \arguments{
# \item{this}{An object.}
# \item{...}{Not used.}
# }
#
# \value{
# TRUE or FALSE
# }
#
# \examples{
#
# # create an object
# o<-Sequence()
# # toggle write protection
# o$writeProtected<-TRUE
# # check if it's write protected
# getWriteProtected(o)
# # check write protection via virtual field
# o$writeProtected
#
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"getWriteProtected",
class="Sequence",
function(
this,
...
){
this$.write.protected;
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: .checkIndexSanity
##
setMethodS3(
".checkIndexSanity",
class="Sequence",
function(
this,
index,
...
){
if (length(index) == 0 ) {
return(c());
}
if(!exists(x="PSIM_FAST")){
if( length(index) == 1 ) {
if(is.na(index)) {
warning("Index vector is NA! Coercing to empty vector!\n");
return(c());
}
if (is.nan(index)) {
warning("Index vector is NaN! Coercing to empty vector!\n");
return(c());
}
}
if(min(index) < 1 ) {
throw("Index vector element ",min(index)," too small!\n");
}
if( max(index) > this$.length ) {
throw("Index vector element ",max(index)," too big!\n");
}
}
return(index);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: setWriteProtected
##
###########################################################################/**
#
# @RdocMethod setWriteProtected
#
# @title "Set the write protection field for an object"
#
# \description{
# @get "title".
# Write protected objects cannot be modified through get/set methods and virtual fields.
# }
#
# @synopsis
#
# \arguments{
# \item{this}{An object.}
# \item{value}{A logical vector of size one.}
# \item{...}{Not used.}
# }
#
# \value{
# Invisible TRUE or FALSE.
# }
#
# \examples{
#
# # create an object
# o<-Sequence()
# # toggle write protection
# setWriteProtected(o,TRUE)
# # check write protection
# o$writeProtected
# # set write protection via virtual field
# o$writeProtected<-FALSE
# o$writeProtected
#
#
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"setWriteProtected",
class="Sequence",
function(
this,
value,
...
){
if(!is.logical(value)) {throw("The new value must be logical!\n")}
else {
this$.write.protected<-value;
}
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: .checkWriteProtection
##
setMethodS3(
".checkWriteProtection",
class="Sequence",
function(
this,
value,
...
){
if(getWriteProtected(this)) {throw("Cannot set value because the object is write protected!\n")}
else {return(FALSE)}
},
private=TRUE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: summary.Sequence
##
###########################################################################/**
#
# @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<-Sequence()
# # get a summary
# summary(a)
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"summary",
class="Sequence",
function(
object,
...
){
this<-object;
this$.summary$"Name"<-this$name;
this$.summary$"Id"<-this$id;
this$.summary$"Length"<-this$length;
this$.summary$"Big rate"<-this$bigRate;
this$.summary$"Ancestral object"<-this$ancestral$id;
if(this$.cumulative.rate.flag) {
this$.summary$"Cumulative rate flag"<-TRUE;
}
if(length(this$.flagged.sites) > 0 ) {
this$.summary$"Flagged sites"<-paste(this$.flagged.sites,collapse=" ");
}
if(this$writeProtected) {
this$.summary$"Write protected"<-TRUE;
}
NextMethod();
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: as.character.Sequence
##
###########################################################################/**
#
# @RdocMethod as.character
#
# @title "Get the string representation of a Sequence object"
#
# \description{
# @get "title".
# The string representation is the concatenation of the states of the
# aggregated Site object. Undefined states (NA-s) are represented by question marks.
# }
#
# @synopsis
#
# \arguments{
# \item{x}{A Sequence object.}
# \item{...}{Not used.}
# }
#
# \value{
# A character vector of length one.
# }
#
# \examples{
# # create object
# s<-Sequence(length=10)
# # get character representation
# as.character(s)
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"as.character",
class="Sequence",
function(
x,
...
){
getString(x);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: plot
##
###########################################################################/**
#
# @RdocMethod plot
#
# @title "Plot the total site rates for a collection of Site objects aggregated by a Sequence object"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{x}{A Sequence object.}
# \item{index}{An integer vector specifying a set of positions. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# Invisible TRUE or FALSE.
# }
#
# \examples{
# # create a nucleotide sequence with a process attached
# s<-NucleotideSequence(string="ATGGCCA",processes=list(list(JC69())))
# # plot total rates in range 1:4
# plot(s,1:4)
# # plot all total rates
# plot(s)
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"plot",
class="Sequence",
function(
x,
index=NA,
...
){
this<-x;
if(this$length == 0) {
warning("The sequence leght is zero, nothing to plot here!\n");
return(invisible(FALSE));
}
if( length(unique(this$totalRates)) == 1 & any(is.na(unique(this$totalRates))) ){
warning("The total rates are undefined, nothing to plot here!\n");
return(invisible(FALSE));
}
else {
if(missing(index)) {
index<-seq(along.with=1:this$length,by=1);
}
if(this$.cumulative.rate.flag){
.recalculateCumulativeRates(this);
}
what<-this$.total.rates[c(index)]
plot(
x=index,
y=what,
type="h",
lwd=1,
col="blue",
main=paste("Total rate plot for sequence", this$id),
xlab="Position",
ylab="Total rate",
ylim=c(0,max(what)),
xlim=c(min(index),max(index)),
xaxt="n"
);
axis(side=1, at=index, labels=index);
return(invisible(TRUE));
}
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: plotParameterAtSites
##
###########################################################################/**
#
# @RdocMethod plotParametersAtSites
#
# @title "Plot the value of a site-process specifc paramter for a collection of Site objects aggregated by a Sequence object"
#
# \description{
# @get "title".
# The type of the paramter must be numeric. The Process object must be attached to all positions specified
# in the index vector.
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{process}{A Process object.}
# \item{id}{The identifier of the site-process specific parameter.}
# \item{index}{An integer vector specifying a set of positions. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# Invisible TRUE or FALSE.
# }
#
# \examples{
# # create a nucleotide sequence with a process attached
# p<-JC69()
# s<-NucleotideSequence(string="ATGGCCA",processes=list(list(p)))
# # plot rate multipliers in range 1:4
# plotParametersAtSites(s,p,"rate.multiplier",1:4)
# # plot rate multiplier for the full sequence
# plotParametersAtSites(s,p,"rate.multiplier")
# }
#
# @author
#
# \seealso{
# Site Process Sequence
# }
#
#*/###########################################################################
setMethodS3(
"plotParametersAtSites",
class="Sequence",
function(
this,
process,
id,
index,
...
){
if(this$length == 0) {
warning("The sequence leght is zero, nothing to plot here!\n");
return(invisible(FALSE));
}
if(missing(index)) {
index<-seq(along.with=1:this$.length,by=1);
}
what<-apply(as.array(index),1,
function(pos){
tmp<-getParameterAtSites(this,process,id,pos)[[1]]$value;
if(!is.numeric(tmp)){
throw("Plot method failed becuase encountered non-numeric parameter value!\n");
}
return(tmp);
}
);
plot(
x=index,
y=what,
type="h",
lwd=1,
col="blue",
main=paste("Plot of parameter",id,"for process",process$id),
xlab="Position",
ylab="Value",
xlim=c(min(index),max(index)),
ylim=c(0,max(what)),
xaxt="n"
);
axis(side=1, at=index, labels=index);
invisible(TRUE);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: setDeletionTolerance
##
###########################################################################/**
#
# @RdocMethod setDeletionTolerance
#
# @title "Set the deletion tolerance site-process specific parameter for a collection of Site objects aggregated by a Sequence object"
#
# \description{
# @get "title".
# This method does some error checking and the calls \code{setParameterAtSites(this=this,process=process,id="deletion.tolerance",value=value,index=index)}.
# See \code{setParameterAtSites.Sequence} for more details.
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{process}{A Process object.}
# \item{value}{A numeric vector, recycled if shorter than the index vector.}
# \item{index}{An integer vector specifying a set of positions. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# The Sequence object (invisible).
# }
#
#
# @author
#
# \seealso{
# setParameterAtSites.Sequence
# }
#
#*/###########################################################################
setMethodS3(
"setDeletionTolerance",
class="Sequence",
function(
this,
process,
value,
index,
...
){
if(!exists(x="PSIM_FAST")){
if(missing(process)){
throw("No process given!\n");
}
else if (missing(value)){
throw("No value provided!\n");
}
else if(!inherits(process,"GeneralDeletor")){
throw("The specified process is not an insertion process!\n");
}
}
setParameterAtSites(this=this,process=process,id="deletion.tolerance",value=value,index=index);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getDeletionTolerance
##
###########################################################################/**
#
# @RdocMethod getDeletionTolerance
#
# @title "Get the deletion tolerance site-process specific parameter for a collection of Site objects aggregated by a Sequence object"
#
# \description{
# @get "title".
# This method does some error checking and the calls \code{getParameterAtSites(this=this,process=process,id="deletion.tolerance",index=index)}.
# See \code{getParameterAtSites.Sequence} for more details.
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{process}{A Process object.}
# \item{index}{An integer vector specifying a set of positions. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# A numeric vector.
# }
#
#
# @author
#
# \seealso{
# getParameterAtSites.Sequence
# }
#
#*/###########################################################################
setMethodS3(
"getDeletionTolerance",
class="Sequence",
function(
this,
process,
index,
...
){
if(!exists(x="PSIM_FAST")){
if(missing(process)){
throw("No process given!\n");
}
if(!inherits(process,"GeneralDeletor")){
throw("The specified process is not an insertion process!\n");
}
}
rm<-getParameterAtSites(this=this,process=process,id="deletion.tolerance",index=index);
return(as.numeric(lapply(rm,function(param){param$value})));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: setInsertionTolerance
##
###########################################################################/**
#
# @RdocMethod setInsertionTolerance
#
# @title "Set the insertion tolerance site-process specific parameter for a collection of Site objects aggregated by a Sequence object"
#
# \description{
# @get "title".
# This method does some error checking and the calls \code{setParameterAtSites(this=this,process=process,id="insertion.tolerance",value=value,index=index)}.
# See \code{setParameterAtSites.Sequence} for more details.
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{process}{A Process object.}
# \item{value}{A numeric vector, recycled if shorter than the index vector.}
# \item{index}{An integer vector specifying a set of positions. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# The Sequence object (invisible).
# }
#
#
# @author
#
# \seealso{
# setParameterAtSites.Sequence
# }
#
#*/###########################################################################
setMethodS3(
"setInsertionTolerance",
class="Sequence",
function(
this,
process,
value,
index,
...
){
if(!exists(x="PSIM_FAST")){
if(missing(process)){
throw("No process given!\n");
}
else if (missing(value)){
throw("No value provided!\n");
}
else if(!inherits(process,"GeneralInsertor")){
throw("The specified process is not an insertion process!\n");
}
}
setParameterAtSites(this=this,process=process,id="insertion.tolerance",value=value,index=index);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getInsertionTolerance
##
###########################################################################/**
#
# @RdocMethod getInsertionTolerance
#
# @title "Get the insertion tolerance site-process specific parameter for a collection of Site objects aggregated by a Sequence object"
#
# \description{
# @get "title".
# This method does some error checking and the calls \code{getParameterAtSites(this=this,process=process,id="insertion.tolerance",index=index)}.
# See \code{getParameterAtSites.Sequence} for more details.
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{process}{A Process object.}
# \item{index}{An integer vector specifying a set of positions. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# A numeric vector.
# }
#
#
# @author
#
# \seealso{
# getParameterAtSites.Sequence
# }
#
#*/###########################################################################
setMethodS3(
"getInsertionTolerance",
class="Sequence",
function(
this,
process,
index,
...
){
if(!exists(x="PSIM_FAST")){
if(missing(process)){
throw("No process given!\n");
}
if(!inherits(process,"GeneralInsertor")){
throw("The specified process is not an insertion process!\n");
}
}
rm<-getParameterAtSites(this=this,process=process,id="insertion.tolerance",index=index);
return(as.numeric(lapply(rm,function(param){param$value})));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: sampleStates
##
###########################################################################/**
#
# @RdocMethod sampleStates
#
# @title "Sample the states for a collection of Site objects aggregated by a Sequence object"
#
# \description{
# @get "title".
# This method samples new states from the equlibrium distribution of the attched process(es) for sites
# having undefined states (NA).
# If a site has more than one substitution process attached, then the method samples the new state from the
# mixture of equlibrium distributions. The weight of each equlibrium distribution is proportional to the
# site-process specific rate multiplier of the corresponding process at the given site.
#
# Sites having defined states are not touched. All sites with undefined states must have at least one
# substitution process (object inheriting from GeneralSubstitution) attached.
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{index}{An integer vector specifying a set of positions. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# The Sequence object (invisible).
# }
#
# \examples{
# # create a nucleotide sequence
# s<-NucleotideSequence(length=80)
# # create some processes
# jc69<-JC69(); # Jukes-Cantor
# hky<-HKY(base.freqs=c(0.1,0.4,0.1,0.4)) # "GC-rich" HKY
# # attach the processes
# s$processes<-list(list(jc69)) # jc is attached to all sites
# attachProcess(s,hky,60:80) # hky is attached to range 60:80
# # tweak rate multiplier for hky
# setRateMultipliers(s,hky,10,60:80)
# # set states in range 1:20
# setStates(s,"A",1:20)
# # sample remaining states
# sampleStates(s)
# # print sequence
# s
# }
#
# @author
#
# \seealso{
# sampleState.GeneralSubstitution GeneralSubstitution
# }
#
#*/###########################################################################
setMethodS3(
"sampleStates",
class="Sequence",
function(
this,
index,
...
){
if(!missing(index)){
index<-.checkIndexSanity(this, index);
}
else {
index<-seq(along.with=this$.sites);
}
for(site in this$.sites[index]){
# Sample states from the equlibrium distributions if
# the state is NA:
if(is.na(site$state)){
# Assemble the list of substitution processes:
subst.proc<-list();
for (proc in site$processes){
if(is.GeneralSubstitution(proc)){
subst.proc<-c(subst.proc, list(proc));
}
}
# Complain if we have no substitution processes to sample from:
if(length(subst.proc) == 0){
throw("Site state is NA and no substitution processes are attached. Cannot sample state!\n");
}
site.rates<-as.numeric(lapply(
subst.proc,
function(proc){
return(getParameterAtSite(proc, site,"rate.multiplier")$value);
}
));
# Normalizing site rates:
site.rates<-site.rates/sum(site.rates);
# Single subst process:
if(length(subst.proc) == 1){
site$state<-sampleState(subst.proc[[1]]);
}
else {
# Sample a substitution process according to the rate multipliers:
nproc<-sample(x=c(1:length(subst.proc)),size=1, replace=FALSE, prob=site.rates);
# Sample the state from the winner process:
site$state<-sampleState(subst.proc[[nproc]]);
}
} # if is.na...
}
return(invisible(this));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: clearStates
##
###########################################################################/**
#
# @RdocMethod clearStates
#
# @title "Set the states of a collection of Site objects aggregated by a Sequence object to undefined (NA)"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{index}{An integer vector specifying a set of positions. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# The Sequence object (invisible).
# }
#
# \examples{
# # create a nucleotide sequence
# s<-NucleotideSequence(string="ATGC")
# s
# # set states to NA in the range 2:3
# clearStates(s,2:3)
# s
# # set all states to NA
# clearStates(s)
# s
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"clearStates",
class="Sequence",
function(
this,
index,
...
){
if(!missing(index)){
index<-.checkIndexSanity(this, index);
}
else {
index<-seq(along.with=this$.sites);
}
for(site in this$.sites[index]){
site$.state<-NA;
site$.total.rate<-NA;
}
this$.cumulative.rate.flag<-TRUE;
return(invisible(this));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: setOmegas
##
setMethodS3(
"setOmegas",
class="Sequence",
function(
this,
process,
index,
...
){
# dummy method to force the creation of the generic function
setOmegas.CodonSequence(this,process,index,...);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getOmegas
##
setMethodS3(
"getOmegas",
class="Sequence",
function(
this,
process,
index,
...
){
# dummy method to force the creation of the generic function
getOmegas.CodonSequence(this,process,index,...);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getSymbolFreqs
##
###########################################################################/**
#
# @RdocMethod getSymbolFreqs
#
# @title "Get a table with the frequencies of the states of a collection of Site objects aggregated by a Sequence object"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{index}{An integer vector specifying a set of positions. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# A table.
# }
#
# \examples{
# # create a nucleotide sequence
# s<-NucleotideSequence(length=30,processes=list(list(JC69())))
# # sample states
# sampleStates(s)
# # get state frequencies from ranges 1:10 and 20:30
# getSymbolFreqs(s,c(1:10,20:30))
# # get symbol frequencies for the full sequence
# getSymbolFreqs(s)
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"getSymbolFreqs",
class="Sequence",
function(
this,
index,
...
){
if(!missing(index)){
index<-.checkIndexSanity(this, index);
} else {
index<-seq(along.with=this$.sites);
}
prop.table(table(as.character(lapply(this$.sites[index],getState))));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: .discretizeGamma
##
setMethodS3(
".discretizeGamma",
class="Sequence",
function(
this,
shape,
ncat,
...
){
# figure out cutting points
cut<-apply(rbind(1:(ncat-1)/ncat),1, function(x){qgamma(x,shape=shape,rate=shape)})
cut<-c(0,cut,Inf);
# incomplete gamma function
Igamma<-function(x,a){
pgamma(x,shape=a, scale=1)
}
# function to calculate a category mean
cm<-function(a,b,shape,ncat){
( Igamma(b * shape, shape+1) - Igamma(a * shape,shape+1)) * ncat;
}
# calculate category means
means<-c();
for (i in 1:(length(cut)-1)){
means<-c(means,cm(cut[i], cut[i+1], shape, ncat));
}
# return a vector with the means
return(means);
},
private=TRUE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: plusGamma
##
###########################################################################/**
#
# @RdocMethod plusGamma
#
# @title "Sample the rate multiplier parameters of a Process from a Gamma distribution for a collection of Site objects aggregated by a Sequence object"
#
# \description{
# @get "title".
#
# The gamma distribution is discretized by calculating the means of the portions corresponding
# to the categories having equal probabilities. If the \code{ncat} argument is not numeric, the
# rates are sampled from the continuous gamma distribution.
# }
#
# \references{
# Yang, Z. (1994) Maximum likelihood phylogenetic estimation from DNA sequences with variable
# rates over sites: approximate methods - Journal of Molecular Evolution 39:306-314 \url{http://dx.doi.org/10.1007/BF00160154}
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{process}{A Process object.}
# \item{shape}{The shape parameter of the gamma distribution.}
# \item{index}{A vector of positions.}
# \item{ncat}{Numer of categories in the discretized gamma distribution (4 by default).}
# \item{...}{Not used.}
# }
#
# \value{
# The Sequence object (invisible);
# }
#
# \examples{
# # create a sequence
# s<-NucleotideSequence(length=20)
# # attach a process
# p<-JC69()
# attachProcess(s,p)
# # get rate multipliers
# getRateMultipliers(s,p) # the default value is 1.0
# # sample rate multipliers in range 1:5 from a discrete
# #gamma distribution with shape parameter 0.5
# plusGamma(s,p,0.5,1:5)
# # get rate multipliers
# getRateMultipliers(s,p) # the default value is 1.0
# # sample rates from a continuous gamma distribution
# # with shape parameter 0.5
# plusGamma(s,p,0.5,ncat="cont")
# # get rate multipliers
# getRateMultipliers(s,p)
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"plusGamma",
class="Sequence",
function(
this,
process,
shape,
index,
ncat=4,
...
){
if(missing(process)){
throw("No process specified!\n");
}
else if(missing(shape)){
throw("No shape parameter specified!\n");
}
else if(!all(is.numeric(shape)) | length(shape) != 1){
throw("The shape parameter must be a numeric vector of lenght 1!\n");
}
if(missing(index)){
index<-seq(along.with=this$.sites);
}
else {
index<-.checkIndexSanity(this, index);
}
if(!is.numeric(ncat)){
# continuous gamma
setParameterAtSites(this, process=process, id="rate.multiplier",value=rgamma(length(index),shape=shape,rate=shape),index=index);
}
else{
setParameterAtSites(this, process=process, id="rate.multiplier",sample(.discretizeGamma(this,shape,ncat),size=length(index),replace=TRUE),index=index);
}
return(invisible(this));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: plusInvGamma
##
###########################################################################/**
#
# @RdocMethod plusInvGamma
#
# @title "Sample the rate multiplier parameters of a Process from an I+Gamma distribution for a collection of Site objects aggregated by a Sequence object"
#
# \description{
# @get "title".
#
# The gamma distribution is discretized by calculating the means of the portions corresponding
# to the categories having equal probabilities. If the \code{ncat} argument is not numeric, the
# rates are sampled from the continuous gamma distribution.
# }
#
# \references{
# Gu X, Fu, YX, Li, WH (1995) Maximum likelihood estimation of the heterogeneity of substitution
# rate among nucleotide sites - Mol. Biol. Evol. 12(4):546-57 \url{http://bit.ly/aE6xF0}
#
# Yang, Z (1994) Maximum likelihood phylogenetic estimation from DNA sequences with variable
# rates over sites: approximate methods - Journal of Molecular Evolution 39:306-314 \url{http://dx.doi.org/10.1007/BF00160154}
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{process}{A Process object.}
# \item{pinv}{The proportion of invariant sites.}
# \item{shape}{The shape parameter of the gamma distribution.}
# \item{index}{A vector of positions.}
# \item{ncat}{Numer of categories in the discretized gamma distribution (4 by default).}
# \item{...}{Not used.}
# }
#
# \value{
# The Sequence object (invisible);
# }
#
# \examples{
# # create a sequence
# s<-NucleotideSequence(length=20)
# # attach a process
# p<-JC69()
# attachProcess(s,p)
# # get rate multipliers
# getRateMultipliers(s,p) # the default value is 1.0
# # sample rate multipliers in range 1:5 from I+discrete Gamma
# plusInvGamma(s,p,pinv=0.5,shape=0.5,1:5)
# # get rate multipliers
# getRateMultipliers(s,p) # the default value is 1.0
# # sample rates from an I+continuos Gamma model
# plusInvGamma(s,p,pinv=0.5,shape=0.5,1:5,ncat="cont")
# # get rate multipliers
# getRateMultipliers(s,p)
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"plusInvGamma",
class="Sequence",
function(
this,
process,
pinv,
shape,
index,
ncat=4,
...
){
if(missing(process)){
throw("No process specified!\n");
}
else if(missing(pinv)){
throw("No proportion of invariant sites given!\n");
}
else if(!all(is.numeric(pinv)) | length(pinv) != 1){
throw("The pinv parameter must be a numeric vector of lenght 1!\n");
}
else if(pinv > 1){
throw("Tpe proportion of invariant sites cannot be larger than 1.!");
}
else if(missing(shape)){
throw("No shape parameter specified!\n");
}
else if(!all(is.numeric(shape)) | length(shape) != 1){
throw("The shape parameter must be a numeric vector of lenght 1!\n");
}
else {
if(missing(index)){
index<-seq(along.with=this$.sites);
}
else {
index<-.checkIndexSanity(this, index);
}
# discretize gamma distribution
dg<-c()
if(is.numeric(ncat)){
dg<-.discretizeGamma(this,shape,ncat);
}
# Iterating over the sites specified by the index vector:
for(site in index){
# Choose between invariant and gamma:
type<-sample(c("INV","GAMMA"),size=1, replace=FALSE, prob=c( pinv, (1-pinv) ) );
if(type == "INV"){
setParameterAtSites(this, process=process, id="rate.multiplier",value=0,index=c(site));
}
else {
if(!is.numeric(ncat)){
# continuous gamma
setParameterAtSites(this, process=process, id="rate.multiplier",value=rgamma(1,shape=shape,rate=shape),index=c(site));
}
else{
setParameterAtSites(this, process=process, id="rate.multiplier",sample(dg,size=1),index=c(site));
}
}
}
}
return(invisible(this));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: insertSequence
##
###########################################################################/**
#
# @RdocMethod insertSequence
#
# @title "Insert a Sequence object into another Sequence object after a specified position"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{The target Sequence object.}
# \item{insert}{The Sequence object to be inserted.}
# \item{position}{The position after the Sequence object will be inserted.}
# \item{process}{The Process object performing the insertion (optional).}
# \item{paranoid}{If TRUE, then the consistency of teh target objects is checked more rigurously after insertion.}
# \item{...}{Not used.}
# }
#
# \value{
# The Sequence object (invisible).
# }
#
# \examples{
# # create some sequence objects.
# target<-NucleotideSequence(string="AAAAAAAAAAA")
# insert<-NucleotideSequence(string="GGGGGGGGGGG")
# # insert after position 5
# insertSequence(target,insert,5)
# # print the target sequence
# target
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"insertSequence",
class="Sequence",
function(
this,
insert,
position,
process=NA,
paranoid=FALSE,
...
){
.checkWriteProtection(this);
if(!exists(x="PSIM_FAST")){
if(missing(insert)) {
throw("Insert sequence object is missing!\n");
}
else if (missing(position)) {
throw("Insertion position is missing!\n");
}
if(!is.Sequence(insert)) {
throw("Insert object not valid!\n");
}
else if (this$length == 0 & position != 0 ) {
throw("Acceptor sequence length is zero! The only valid insertion position is 0!\n");
}
else if ( !( position >= 0 & position <=(this$.length + 1))) {
throw("Insertion position ",position," is invalid!\n");
}
}
# Just return if insert has zero length:
if(insert$length == 0){
warning("The length of the sequence to be inserted is zero! Nothing to do here!\n");
return(invisible(FALSE));
}
# Clone insert object:
insert<-clone(insert);
# Set the generator process:
if(!missing(process)) {
if( (length(process) == 0) | !is.Process(process)){
throw("Process object invalid!\n");
}
} else {
process<-Sequence$.root.ins;
}
for(site in insert$.sites) {
site$.ancestral<-process;
site$.sequence<-this;
}
# Recalculate cumulative rates if the flag is on:
if(this$.cumulative.rate.flag) {
.recalculateCumulativeRates(this);
}
# Flagging cumulative rates:
this$.cumulative.rate.flag<-TRUE;
# Inserting new site objects:
if ( position == this$.length) {
# Insertion at the end of the sequence;
this$.sites<-c(this$.sites,insert$.sites);
this$.total.rates<-c(this$.total.rates,rep(c(NA),times=insert$.length) );
this$.cumulative.rates<-c(this$.cumulative.rates,rep(NA,times=insert$.length) );
} else if (position == 0) {
# Insertion in the sequence
this$.sites<-c(insert$.sites, this$.sites);
this$.total.rates<-c(rep(NA,times=insert$.length),this$.total.rates);
this$.cumulative.rates<-c(rep(NA,times=insert$.length),this$.cumulative.rates);
} else {
# Insertion at position 0
this$.sites<-c(this$.sites[1:position],insert$.sites,this$.sites[(position+1):this$.length]);
this$.total.rates<-c(this$.total.rates[1:position],rep(NA,times=insert$.length),this$.total.rates[(position+1):this$.length]);
this$.cumulative.rates<-c(this$.cumulative.rates[1:position],rep(NA,times=insert$.length),this$.cumulative.rates[(position+1):this$.length]);
}
# Checking if lengths are consistent:
if(!exists(x="PSIM_FAST")){
if(length(this$.sites) != (this$.length + insert$.length)) {
throw("Length inconsistency after insertion!\n");
}
}
# Setting new length:
this$.length<-(this$.length + insert$.length);
# Flagging the inserted sites:
this$.flagged.sites<-c(this$.flagged.sites,(position+1):(position+insert$.length));
if(!exists(x="PSIM_FAST")){
if(length(this$.total.rates) != this$.length) {
throw("Total rates vector inconsistency after insertion!\n");
}
if(length(this$.cumulative.rates) != this$.length) {
throw("Cumulative rates vector inconsistency after insertion!\n");
}
}
# Recalculating cumulative rates:
.recalculateCumulativeRates(this);
# Paranoid check of total rates:
if(paranoid) {
for (i in 1:this$.length) {
if(this$.sites[[i]]$totalRate != this$.total.rates[[i]]) {
throw("Object total rates inconsistent with total rates vector!\n");
}
}
}
# Deleting the insert:
rm(insert);
return(invisible(this));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: deleteSubSequence
##
###########################################################################/**
#
# @RdocMethod deleteSubSequence
#
# @title "Delete a collection of sites aggregated by a Sequence object"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{index}{An index vector specifying a collection of sites to be deleted. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# Invisible TRUE or errror.
# }
#
# \examples{
# # create a nucleotide sequence
# s<-NucleotideSequence(string="ATATATATATATATAT")
# # delete sites 2, 4 and 6
# deleteSubSequence(s,c(2,4,6))
# s
# # delete sites in the range 3:6
# deleteSubSequence(s,3:6)
# s
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"deleteSubSequence",
class="Sequence",
function(
this,
index,
...
){
if(length(index) == 0) {
return(FALSE);
}
.checkWriteProtection(this);
if(!exists(x="PSIM_FAST")){
if(missing(index)) {
throw("No index vector specified!\n");
}
if(length(index) == 0) {
return(FALSE);
}
}
index<-.checkIndexSanity(this, index);
# Avoid deletion on dirty sequence as
# that will cause havoc.
if(this$.cumulative.rate.flag){
.recalculateCumulativeRates(this);
}
# Flagging cumulative rates:
this$.cumulative.rate.flag<-TRUE;
min.index<-min(index);
# Deleting site objects:
this$.sites[index]<-NULL;
# Updating rate vectors:
this$.total.rates<-this$.total.rates[-index];
this$.cumulative.rates<-this$.cumulative.rates[-index];
# Flag the site before the deletion to
# to force cumulative rate recalculation:
if (min.index > 2 ) {
this$.flagged.sites<-c(this$.flagged.sites,(min.index - 1));
}
if(!exists(x="PSIM_FAST")){
if( length(this$.sites) != (this$.length - length(index) ) ) {
throw("Inconsistency after deleting sites!\n");
}
}
this$.length<-length(this$.sites);
.recalculateCumulativeRates(this);
return(invisible(TRUE));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: copySubSequence
##
###########################################################################/**
#
# @RdocMethod copySubSequence
#
# @title "Copy a collection of Site objects aggregated by a Sequence object into a new Sequence object"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{index}{An index vector specifying a collection of sites to be copied. It is set to 1:seq$length if omitted.}
# \item{process}{The Process object performing the copy (optional).}
# \item{...}{Not used.}
# }
#
# \value{
# A new Sequence object.
# }
#
# \examples{
# # create a nucleotide sequence
# s<-NucleotideSequence(string="ATATATATATATATATA")
# # copy sites in the range 3:8 in a new object
# s2<-copySubSequence(s,3:8)
# s2
# # copy sites 1,3 and 5 from s2
# s3<-copySubSequence(s2,c(1,3,5))
# s3
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"copySubSequence",
class="Sequence",
function(
this,
index,
process=NA,
...
){
if(missing(index)) {
index<-seq(along.with=this$.sites);
}
else {
index<-.checkIndexSanity(this, index);
}
if(!exists(x="PSIM_FAST")){
if(!is.na(process) & !is.Process(process)) {
throw("Process object invalid!\n");
}
}
# Avoid copying from dirty sequence:
if(this$.cumulative.rate.flag){
.recalculateCumulativeRates(this);
}
length<-length(index);
# Create an empty sequence object:
copy<-Sequence();
# Flag copy cumulative rates:
copy$.cumulative.rate.flag<-TRUE;
if(is.na(process)){
# Getting the root insertion process:
process<-Sequence$.root.ins;
}
# Setting the ancestral to sequence:
copy$.ancestral.obj<-process;
# Setting copy name:
copy$name<-paste("Copied from",this$name);
# Setting length:
copy$.length<-length;
# Clone the sites:
copy$.sites<-lapply(this$.sites[index],
function(site){
site.copy<-clone(site);
site.copy$.ancestral<-process;
return(site.copy);
}
);
# Copy total rates:
copy$.total.rates<-this$.total.rates[index];
# Create cumulative rates vector:
copy$.cumulative.rates<-cumsum(copy$.total.rates);
copy$.cumulative.rate.flag<-FALSE;
if(!exists(x="PSIM_FAST")){
if(length(copy$.sites) != length){
throw("Sites list length mismatch!\n")
}
else if(length(copy$.total.rates) != length){
throw("Total rates vector length mismatch!\n")
}
else if(length(copy$.cumulative.rates) != length){
throw("Cumulative rates vector length mismatch!\n")
}
}
return(copy);
},
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 file ../COPYING for licensing issues.
##
##########################################################################/**
#
# @RdocClass Sequence
#
# @title "The Sequence class"
#
# \description{
#
# This is the class representing a sequence. The backbone of the Sequence objects are
# lists aggregating Site objects. The class has fields for keeping track of cumulative
# site rates, the sum of all active event rates and methods for performing actions
# on a collection of sites (positions).
#
# The Sequence objects have a field specifying an ancestral object, which can be a Sequence
# object (when the object is obtained through clone() ) or the "Root insertion process" object
# (for newly created objects).
#
# @classhierarchy
# }
#
# @synopsis
#
# \arguments{
# \item{name}{The name of the Sequence object.}
# \item{string}{A string containing symbols belonging to the associated Alphabet object.
# It can be used to set the initial states of the aggregated Site objects. It also specifies the length of the sequence}
# \item{length}{The length of the sequence. It cannot be used when 'string' is specified.}
# \item{alphabets}{A list of Alphabet objects to be associated with the Site objects.
# The list is recycled in the case it is shorter than the sequence length.}
# \item{processes}{A list of lists of Process objects to be attached
# (recycled if shorter than sequence length). }
# \item{ancestral.obj}{The ancestral object (Sequence or Process).}
# \item{...}{Not used.}
# }
#
# \section{Fields and Methods}{
# @allmethods
# }
#
# \examples{
# # create a sequence object by
# # providng alphabets, processes and states
# s.one<-Sequence(
# name="Seq",
# string="AATTGGCCTTAAGGCCTTAA",
# alphabets=list(NucleotideAlphabet()),
# processes=list(list(JC69()))
# )
# s.one
# # check if inherits from Sequence
# is.Sequence(s.one)
# # get object summary
# summary(s.one)
# # create a sequence object,
# # specifying length, alphabets
# # and ancestral object
# s<-Sequence(
# name="Seq",
# length=20,
# ancestral.obj=s.one
# )
# # get sequence string
# s$string
# # get the list of site objects
# s$sites
# # get object id
# s$id
# # set and get name
# s$name<-"SeqSeq"
# s$seq
# # get length
# s$length
# # get and set ancestral object
# s$ancestral
# s$ancestral<-Sequence();
# # set alphabets
# setAlphabets(s,list(NucleotideAlphabet()))
# # set states
# # "A"-s in the range 1:10
# setStates(s,"A",1:10)
# # a pattern of "ATGC"-s in the range 11:20
# setStates(s,c("A","T","G","C"),11:20)
# # get states from range 10:12
# getStates(s,10:12)
# # attach a JC69 process to range 1:10
# jc69<-JC69()
# attachProcess(s,jc69,1:10)
# # set the rate multiplier site-process specific parameter for jc69
# setParameterAtSites(s,jc69,"rate.multiplier",2,1:10)
# # get "rate.multiplier" for jc69 from range 1:2
# getParameterAtSites(s, jc69, "rate.multiplier",1:2)
# # attach a GTR process to range 11:20
# gtr<-GTR()
# attachProcess(s,gtr,11:20)
# # set and get rate multipliers for gtr
# setRateMultipliers(s, gtr, c(0,5,1), 11:20)
# getRateMultipliers(s, gtr, 11:20)
# # get processes from range 1:5
# getProcesses(s,1:5)
# # replace the processes with a (GTR, JC69), JC69 pattern
# s$processes<-list(list(GTR(), JC69()), list(JC69()))
# # get processes from range 1:2
# getProcesses(s,1:2)
# # get unique processes
# s$uniqueProcesses
# # get unique alphabets
# s$uniqueAlphabets
# # get symbol frequencies
# getSymbolFreqs(s)
# # get the big rate
# s$bigRate
# # get a vector of total rates from range 1:4
# getTotalRatesFromRange(s,1:4)
# # get a vector of cumulative rates from range 1:4
# getCumulativeRatesFromRange(s,1:4)
# # reset all site states
# clearStates(s)
# s
# # sample states from the equilibrium distributions
# # of the attached substitution processes
# sampleStates(s)
# s
# # clone s
# s.clone<-clone(s)
# # insert a sequence in s.clone after position 2
# insertSequence(s.clone,NucleotideSequence(string="AAAAAA"),2)
# s.clone
# # delete positions 1,2,3 and 10
# deleteSubSequence(s.clone, c(1:3,10))
# s.clone
# # copy positions 7:10 into a new sequence object
# sc<-copySubSequence(s.clone, 7:10)
# sc
# # get events from sc in the range 1:2
# getEvents(sc,1:2)
# }
#
# @author
#
# \seealso{
# Alphabet Site Process Event
# }
#
#*/###########################################################################
setConstructorS3(
"Sequence",
function(
name=NA,
string=NA,
length=NA,
alphabets=NA,
processes=NA,
ancestral.obj=NA,
...
){
# Marking the instance as static by default:
STATIC<-TRUE;
# Extending the PSRoot class:
this<-extend(
PSRoot(),
"Sequence",
.name="Anonymous",
.length=NA,
.sites=list(),
.ancestral.obj=NA,
.cumulative.rates=NA,
.total.rates=NA,
.cumulative.rate.flag=TRUE,
.flagged.sites=integer(0),
.write.protected=FALSE,
.is.sequence=TRUE,
.root.ins=NA
);
# Initializing the variables for length and states:
len<-0;
str<-list();
# Optional argument: name
if(!missing(name)) {
this$name<-name;
STATIC<-FALSE;
}
# The user can specify a sequence
# or the sequence length, but not both.
if (!missing(string) & !missing(length)) {
throw("You can specify the sequence, or the sequence length, but not both!\n");}
else if (!missing(string)){
STATIC<-FALSE;
# An alphabet list must be specified along the sequence!
if(missing(alphabets)) {throw("A list of valid alphabets must be specified when a string is given!\n");}
}
# Deal with the string or length argument:
if (!missing(length)) {
STATIC<-FALSE;
len<-length;
}
else if( !missing(string) ) {
str<-strsplit(string,split="",fixed=TRUE)[[1]];
len<-length(str);
}
this$.length<-len;
root.ins<-NA;
if (!is.Process(Sequence$.root.ins)) {
# Dummy proces used as ancestral object for sites.
this$.root.ins<-Process(name="Root insertion process");
this$.root.ins$comments<-"This is just a dummy process object serving as ancestral for newly created site and sequence objects.";
root.ins<-this$.root.ins;
} else {
root.ins<-Sequence$.root.ins
}
# Site template object:
site.template<-Site(
ancestral=root.ins,
sequence=this
);
# Clone and store the site objects:
if(!STATIC) {
if ( len > 0 ) {
for(position in 1:len) {
this$.sites[[position]]<-clone.Object(site.template);
}
}
}
# Optional argument: ancestral object
if (!missing(ancestral.obj)) {
STATIC<-FALSE;
this$ancestral<-ancestral.obj;
} else {
this$ancestral<-root.ins;
}
# Setting the alphabets:
if(!missing(alphabets)) {
STATIC<-FALSE;
# setAlphabets will check the arguments
setAlphabets(this, alphabets);
}
# Setting the processes:
if (!missing(processes)) {
STATIC<-FALSE;
# setProcesses will take care about the arguments
# and deal with alphabet mismatch.
setProcesses(this,processes)
}
# Initializing these vectors properly is
# importtant for the insertion method!
this$.total.rates<-double(len);
this$.cumulative.rates<-double(len);
if(!STATIC){
# Now we are prepared to set the states:
if (!missing(string)) {
setStates(this, str);
}
# Calculate cumulative rates for the first time, but only if
# states are defined. This is expensive, as total rates are calculated.
if (!missing(string) & (length(str) > 0) ) {
.recalculateCumulativeRates(this);
}
}
return(this);
},
enforceRCC=TRUE
);
##
## Method: is.Sequence
##
###########################################################################/**
#
# @RdocDefault is.Sequence
#
# @title "Check whether an object inherits from the Sequence class"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{An object.}
# \item{...}{Not used.}
# }
#
# \value{
# TRUE or FALSE
# }
#
# \examples{
# # create some objects
# seq<-Sequence(length=10)
# a<-Alphabet()
# # check if they inherit from Sequence
# is.Sequence(seq)
# is.Sequence(a)
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"is.Sequence",
class="default",
function(
this,
...
){
if(!is.PSRoot(this)) {return(FALSE)}
if(!is.null(this$.is.sequence)){return(TRUE)}
if ( inherits(this, "Sequence")) {
this$.is.sequence<-TRUE;
return(TRUE);
} else {
return(FALSE)
}
},
private=FALSE,
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{omit.sites}{Do not check aggregated site objects.}
# \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="Sequence",
function(
this,
omit.sites=FALSE,
...
){
if(this$length != length(this$.sites)) {
throw("Sequence length inconsistency detected");
} else if (length(this$.cumulative.rates) != this$.length) {
throw("Cumulative rates vector length mismatch!\n");
} else if (length(this$.total.rates) != this$.length) {
throw("Total rates vector length mismatch!\n");
} else if (!identical(this$.cumulative.rates, cumsum(this$.total.rates))) {
throw("Cumulative rates vector is not consistent with total rates vector!\n");
}
if(!is.numeric(this$.flagged.sites)) {
throw("Flagged sites vector is not numeric!\n");
} else if (length(this$.flagged.sites) > 0) {
if ( (min(this$.flagged.sites) < 1) | ( max(this$.flagged.sites) > this$.length) ) {
throw("Inconsistency in the flagged sites vector!\n");
}
}
if(!is.character(this$.name)) {
throw("Sequence name is invalid!\n");
} else if(stringLength(this$.name) == 0) {
throw("Sequence name is of length zero!\n");
}
if(!is.Sequence(this$.ancestral.obj) & !is.Process(this$.ancestral.obj)) {
throw("The ancestral object is invalid!\n");
}
if(!is.logical(this$.cumulative.rate.flag)) {
throw("Cumulative rate flag is not logical!\n");
}
# Calling consistency check on sites.
# This will be painfully slow!
if(!omit.sites){
for(site in this$.sites){
checkConsistency(site);
}
}
return(invisible(TRUE));
},
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 Sequence object"
#
# \description{
# @get "title".
# The identifier is the concatenation of the object name and the object hash code as returned
# by hashCode().
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{...}{Not used.}
# }
#
# \value{
# A charcter vector of length one.
# }
#
# \examples{
# # create a Sequence object.
# s<-Sequence(length=5)
# # get id
# getId(s)
# # get id via virtual field
# s$id
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"getId",
class="Sequence",
function(
this,
...
){
this.class<-class(this)[1];
id<-paste(this.class,this$.name,hashCode(this),sep=":");
},
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 Sequence 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="Sequence",
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: getName
##
###########################################################################/**
#
# @RdocMethod getName
#
# @title "Get the name of a Sequence object"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{...}{Not used.}
# }
#
# \value{
# A character vector of length one.
# }
#
# \examples{
# # create object
# s<-Sequence(length=10);
# # get object name
# getName(s)
# # get name via virtual field
# s$name
#
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"getName",
class="Sequence",
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 Sequence object"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{new.name}{A character vector of length one.}
# \item{...}{Not used.}
# }
#
# \value{
# The new object name (invisible).
# }
#
# \examples{
# # create a Sequence object
# s<-Sequence(name="MySeq",length=4)
# # get sequence name
# s$name
# # rename object
# setName(s,"SO")
# s$name
# # rename via virtual field
# s$name<-"SeqSeq"
# s$name
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"setName",
class="Sequence",
function(
this,
new.name,
...
){
.checkWriteProtection(this);
this$.name<-as.character(new.name);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getString
##
###########################################################################/**
#
# @RdocMethod getString
#
# @title "Get the string representation of a Sequence object"
#
# \description{
# @get "title".
# The string representation is the concatenation of the states of the
# aggregated Site object. Undefined states (NA-s) are represented by question marks.
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{...}{Not used.}
# }
#
# \value{
# A character vector of length one.
# }
#
# \examples{
# # create object
# s<-Sequence(length=10)
# # get character representation
# getString(s) # a bunch of '?'-s
# # get string reperesentation via virtual field
# s$string
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"getString",
class="Sequence",
function(
this,
...
){
str<-character();
for(site in this$.sites){
if(is.na(site$.state)){
if(is.Alphabet(site$.alphabet)){
str<-c(str,rep("?",site$.alphabet$.symbolLength));
}
else {
str<-c(str,"?");
}
}
else {
str<-c(str,site$.state);
}
}
return(paste(str,collapse=""));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: setString
##
###########################################################################/**
#
# @RdocMethod setString
#
# @title "Forbidden action: setting the string representation of a Sequence 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(
"setString",
class="Sequence",
function(
this,
value,
...
){
virtualAssignmentForbidden(this);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getLength
##
###########################################################################/**
#
# @RdocMethod getLength
#
# @title "Get the number of Site objects aggregated in a Sequence object"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{...}{Not used.}
# }
#
# \value{
# An integer vector of length one.
# }
#
# \examples{
# # create a Sequence object
# s<-Sequence(length=5)
# # get sequence length
# getLength(s)
# # get length via virtual field
# s$length
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"getLength",
class="Sequence",
function(
this,
...
){
this$.length;
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: setLength
##
###########################################################################/**
#
# @RdocMethod setLength
#
# @title "Forbidden action: setting the length of a Sequence object"
#
# \description{
# @get "title".
# The length of the Sequence objects can be specified when the object is constructed,
# or modified later by the "insertSequence" and "deleteSubSequence" methods.
# }
#
# @synopsis
#
# \arguments{
# \item{this}{An object.}
# \item{value}{Not used.}
# \item{...}{Not used.}
# }
#
# \value{
# Throws an error.
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"setLength",
class="Sequence",
function(
this,
value,
...
){
virtualAssignmentForbidden(this);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getSites
##
###########################################################################/**
#
# @RdocMethod getSites
#
# @title "Get the list of the Site object aggregated in a Sequence object"
#
# \description{
# @get "title".
# Warning: there is no setSites method!
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{...}{Not used.}
# }
#
# \value{
# A list of Site objects.
# }
#
# \examples{
# # create a sequence object
# s<-Sequence(alphabets=list(NucleotideAlphabet()),string="AATTGCCC")
# # get the list of aggregated Site objects
# getSites(s)
# # get Site objects via virtual field
# s$sites
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"getSites",
class="Sequence",
function(
this,
...
){
this$.sites;
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getStates
##
###########################################################################/**
#
# @RdocMethod getStates
#
# @title "Get the states of a set of Site objects aggregated by a Sequence object"
#
# \description{
# @get "title".
#
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{index}{An integer vector specifying a set of positions.
# It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# A character vector.
# }
#
# \examples{
# # create a sequence object
# s<-Sequence(alphabets=list(NucleotideAlphabet()),string="AATTGCCCCCTTGG")
# # get all Site states
# getStates(s)
# # get the states for a collection of sites
# getStates(s,c(1:3,5,8))
# # get states via virtual field
# s$states
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"getStates",
class="Sequence",
function(
this,
index,
...
){
if (missing(index)) {
index<-seq(along.with=this$.sites);
} else {
index<-.checkIndexSanity(this, index);
}
lapply(
this$.sites[index],
function(site) {
site$.state;
}
);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: setStates
##
###########################################################################/**
#
# @RdocMethod setStates
#
# @title "Set the states for a set of Site objects aggregated by a Sequence object"
#
# \description{
# @get "title".
# The value vector is recycled, which is useful when creating repeated patterns.
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{value}{A character vector containg the states (recycled if shorter than the index vector). The states must be compatible with the corresponding Alphabet object.}
# \item{index}{An integer vector specifying a set of positions. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# The Sequence object (invisible).
# }
#
# \examples{
# # create a sequence object of length 10
# s<-Sequence(alphabets=list(NucleotideAlphabet()),length=10)
# # set the states in some ranges
# setStates(s,c("A","T","A"),index=1:5)
# setStates(s,c("G","C"),index=6:10)
# # display sequence
# s
# # set the states for the whole Sequence object
# setStates(s,c("A","T","T","A"))
# s
# # set states via virtual field
# s$states<-c("A","T")
# s
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"setStates",
class="Sequence",
function(
this,
value,
index,
...
){
.checkWriteProtection(this);
if(missing(value)) {
throw("No new values specified!\n");
}
else if (missing(index)) {
index<-seq(along.with=this$.sites);
} else {
index<-.checkIndexSanity(this, index);
}
# Recycling value vector by using rep().
if (length(value) < length(this$.sites) ) {
value<-rep(as.character(value),length.out=length(index))
}
for (i in 1:length(index)) {
this$.sites[[ index[[i]] ]]$state<-value[i];
}
# Flagging the changed sites:
this$.cumulative.rate.flag<-TRUE;
this$.flagged.sites<-c(this$.flagged.sites, index);
.recalculateCumulativeRates(this);
invisible(this);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getAlphabets
##
###########################################################################/**
#
# @RdocMethod getAlphabets
#
# @title "Get the list of the Alphabet objects attached to the Site objects aggregated by a Sequence object"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{index}{An integer vector specifying a set of positions. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# A list of Alphabet objects.
# }
#
# \examples{
# # create a Sequence object with NucleotideAlphabet
# #and BinaryAlphabet objects attached
# s<-Sequence(alphabets=list(NucleotideAlphabet(),BinaryAlphabet()),length=5)
# setStates(s,c("A","0"))
# # get the list of attached Alphabet objects
# getAlphabets(s)
# # get Alphabets from a range
# getAlphabets(s,c(2:3,5))
# # get alphabets via virtual field
# s$alphabets
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"getAlphabets",
class="Sequence",
function(
this,
index,
...
){
if (missing(index)) {
index<-seq(along.with=this$.sites);
} else {
index<-.checkIndexSanity(this, index);
}
lapply(
this$.sites[index],
function(site) {
site$alphabet;
}
);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: setAlphabets
##
###########################################################################/**
#
# @RdocMethod setAlphabets
#
# @title "Assotiate Alphabet objects to a set of Site objects aggregated by a Sequence object"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{value}{A list of Alphabet objects, recycled if shorter than the index vector.}
# \item{index}{An integer vector specifying a set of positions. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# The Sequence object (invisible).
# }
#
# \examples{
# # create a Sequence object
# s<-Sequence(length=10)
# # set the alphabets for range 1:5
# setAlphabets(s,list(NucleotideAlphabet(),BinaryAlphabet()),1:5)
# # set the alphabets for range 6:10
# setAlphabets(s,list(AminoAcidAlphabet()),6:10)
# # get the list of attached Alphabet objects
# getAlphabets(s)
# # get Alphabets from a range
# getAlphabets(s,c(2:3,5))
# # set alphabets via virtual field
# s$alphabets<-list(BinaryAlphabet(),NucleotideAlphabet())
# # get alphabets via virtual field
# s$alphabets
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"setAlphabets",
class="Sequence",
function(
this,
value,
index,
...
){
.checkWriteProtection(this);
if(missing(value)) {
throw("No new values specified!\n");
}
if(!is.list(value)) {
throw("The value parameter must be a list!\n");
} else {
for(a in value) {
if(!is.Alphabet(a)) {
throw("The value parameter must be a list of valid alphabet objects!\n");
}
}
}
if (missing(index)) {
index<-seq(along.with=this$.sites);
} else {
index<-.checkIndexSanity(this, index);
}
# Recycling value vector. rep() cannot be used here,
# because we loose the object references!
value.counter<-1;
for (i in index) {
if(value.counter > length(value)) {
value.counter<-1;
}
setAlphabet(this$.sites[[i]], value[[value.counter]]);
value.counter<-(value.counter + 1);
}
invisible(this);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getUniqueAlphabets
##
###########################################################################/**
#
# @RdocMethod getUniqueAlphabets
#
# @title "Get the list of unique Alphabet objects associated to Site objects aggaregated by a Sequence object"
#
# \description{
# @get "title".
# The returned list contains unique instances of the Alphabet class. The symbol sets are not compared, so
# two instances of the same class are considered to be different.
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{...}{Not used.}
# }
#
# \value{
# A list of Alphabet objects.
# }
#
# \examples{
# # create a Sequence object with some Alphabet objects attached
# s<-Sequence(
# alphabets=list(NucleotideAlphabet(),
# BinaryAlphabet(),
# NucleotideAlphabet()),
# length=10
# )
# # get the list of attached alphabets
# s$alphabets
# # get the unique list of attahced Alphabet objects
# getUniqueAlphabets(s)
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"getUniqueAlphabets",
class="Sequence",
function(
this,
...
){
tmp<-list();
lapply(
this$.sites,
function(site) {
tmp<<-c(tmp,list(site$.alphabet))
}
);
return(unique(tmp));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: setUniqueAlphabets
##
###########################################################################/**
#
# @RdocMethod setUniqueAlphabets
#
# @title "Forbidden action: setting the list of unique Alphabet objects attached to the Site object aggregated by a Sequence 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(
"setUniqueAlphabets",
class="Sequence",
function(
this,
value,
...
){
virtualAssignmentForbidden(this);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: attachProcess
##
###########################################################################/**
#
# @RdocMethod attachProcess
#
# @title "Attach a Process object to a set of Site objects aggregated by a Sequence object"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{process}{A Process object.}
# \item{index}{An integer vector specifying a set of positions. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# The Sequence object (invisible).
# }
#
# \examples{
# # create a Sequence object of length 6
# s<-Sequence(length=10,alphabets=list(NucleotideAlphabet()))
# # attach a JC69 substitution process
# attachProcess(s,JC69())
# # get the list of attached processes
# s$processes
# # attach the GTR substitution process to range 3:6
# attachProcess(s,GTR(),3:6)
# # get the list of attached processes
# s$processes
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"attachProcess",
class="Sequence",
function(
this,
process,
index,
...
){
.checkWriteProtection(this);
if(!exists(x="PSIM_FAST")){
if(!is.Process(process)){
throw("Process object invalid!\n");
}
}
if (missing(index)) {
index<-seq(along.with=this$.sites);
} else {
index<-.checkIndexSanity(this, index);
}
for(i in index){
attachProcess(this$.sites[[i]],process);
}
return(invisible(this));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: detachProcess
##
###########################################################################/**
#
# @RdocMethod detachProcess
#
# @title "Detach a Process object from a set of Site objects aggregated by a Sequence object"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{process}{A Process object.}
# \item{index}{An integer vector specifying a set of positions. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# The Sequence object (invisible).
# }
#
# \examples{
# # create a sequence object with two processes attached
# p<-JC69()
# s<-Sequence(length=4,alphabets=list(NucleotideAlphabet()),processes=list(list(p,K80())))
# # get the list of attached processes
# s$processes
# # detach JC69 from range c(1,4)
# detachProcess(s,p,c(1,4))
# s$processes
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"detachProcess",
class="Sequence",
function(
this,
process,
index,
...
){
.checkWriteProtection(this);
if(!exists(x="PSIM_FAST")){
if(!is.Process(process)){
throw("Process object invalid!\n");
}
}
if (missing(index)) {
index<-seq(along.with=this$.sites);
} else {
index<-.checkIndexSanity(this, index);
}
lapply(
this$.sites[index],
function(site) {
detachProcess(site,process);
}
);
return(invisible(this));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getProcesses
##
###########################################################################/**
#
# @RdocMethod getProcesses
#
# @title "Get the Process objects attached to the Site objects aggregated by a Sequence object"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{index}{An integer vector specifying a set of positions. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# A list of lists of Process objects.
# }
#
# \examples{
# # create a sequence object with some processes attached
# s<-Sequence(
# length=4,
# alphabets=list(NucleotideAlphabet()),
# processes=list(list(JC69(),K80()),list(GTR()))
# )
# # get the list of lists of attached processes from positions 1 and 3
# getProcesses(s,c(1,3))
# # get processes via virtual field
# s$processes
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"getProcesses",
class="Sequence",
function(
this,
index,
...
){
if (missing(index)) {
index<-seq(along.with=this$.sites);
} else {
index<-.checkIndexSanity(this, index);
}
lapply(
this$.sites[index],
function(site) {
site$processes;
}
);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getUniqueProcesses
##
###########################################################################/**
#
# @RdocMethod getUniqueProcesses
#
# @title "Get the list of unique Process instances attached to the Site objects aggregated by a Sequence object"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{...}{Not used.}
# }
#
# \value{
# A list of Process objects.
# }
#
# \examples{
# # create a sequence object and attach processes
# p<-JC69()
# s<-Sequence(
# length=4,
# alphabets=list(NucleotideAlphabet()),
# processes=list(list(p,K80()),list(p))
# )
# # get the unique list of attached Process instances
# getUniqueProcesses(s)
# # via virtual field
# s$uniqueProcesses
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"getUniqueProcesses",
class="Sequence",
function(
this,
...
){
tmp<-list();
lapply(
this$.sites,
function(site) {
tmp<<-c(tmp,site$processes)
}
);
return(unique(tmp));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: setUniqueProcesses
##
###########################################################################/**
#
# @RdocMethod setUniqueProcesses
#
# @title "Forbidden action: setting the list of unique Process instances attached to the sites of a Sequence 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(
"setUniqueProcesses",
class="Sequence",
function(
this,
value,
...
){
virtualAssignmentForbidden(this);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: setProcesses
##
###########################################################################/**
#
# @RdocMethod setProcesses
#
# @title "Specify a set of Process objects to be attached to a set of Site objects aggregated by a Sequence object"
#
# \description{
# @get "title".
# The Process objects in a given inner list correspond to the set of processes to be attached to one Site object.
# Process objects already attached to a given Site are skipped. Attached processes which are not memebers of the list
# are detached, so specifying an empty list will detach all processes.
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{value}{A list of list of Process objects, recycled if shorter than the index vector.}
# \item{index}{An integer vector specifying a set of positions. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# The Sequence object (invisible).
# }
#
# \examples{
# # create a sequence of length 5
# s<-Sequence(length=5,alphabets=list(NucleotideAlphabet()));
# # set a pattern of processes
# setProcesses(s,list(list(JC69(),K81())))
# # get attached processes
# s$processes
# # detach all processes from range 1:3
# setProcesses(s,list(list()),1:3)
# s$processes
# # detach all processes via virtual field
# s$processes<-list(list())
# # create a process pattern in the full sequence via virtual field
# s$processes<-list(list(JC69()),list(GTR(),K80()))
# s$processes
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"setProcesses",
class="Sequence",
function(
this,
value,
index,
...
){
.checkWriteProtection(this);
if(!is.list(value)) {
throw("The value parameter must be a list!\n");
} else {
lapply(
value,
function(element) {
if(!is.list(element)){
throw("The value parameter must be a list of lists containing process objects!\n");
}
}
);
}
if (missing(index)) {
index<-seq(along.with=this$.sites);
} else {
index<-.checkIndexSanity(this, index);
}
value.counter<-1;
# Recycling value vector. rep() cannot be used here,
# because we loose the object references!
for (i in index) {
if(value.counter > length(value)) {
value.counter<-1;
}
setProcesses(this$.sites[[i]], value[[value.counter]]);
value.counter<-(value.counter + 1);
}
invisible(this);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: setParameterAtSites
##
###########################################################################/**
#
# @RdocMethod setParameterAtSites
#
# @title "Set the values of a site-process specific paramater for a process and a collection of Site objects aggregated by a Sequence object"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{process}{A valid Process object.}
# \item{id}{The identifier of the site-process specific parameter.}
# \item{value}{A vector containing the new values of the site-process specific parameter, recycled if shorter than the index vector. It should be consistent with the type of the parameter.}
# \item{index}{An integer vector specifying a set of positions. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# The Sequence object (invisible).
# }
#
# \examples{
# # create a sequence, attach a process
# p<-K80()
# s<-Sequence(length=6,alphabets=list(NucleotideAlphabet()),processes=list(list(p)))
# # set a new pattern of rate multipliers in the range 1:3,
# # the default value is 1.0 by the way
# setParameterAtSites(s,p,"rate.multiplier",c(2,3),1:3)
# # get rate multipliers
# getParameterAtSites(s,p,"rate.multiplier")
# # set a new value for the whole sequence
# setParameterAtSites(s,p,"rate.multiplier",0.5)
# # get rate multipliers
# getParameterAtSites(s,p,"rate.multiplier")
# }
#
# @author
#
# \seealso{
# Site Process @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"setParameterAtSites",
class="Sequence",
function(
this,
process,
id,
value,
index,
...
){
.checkWriteProtection(this);
if(!exists(x="PSIM_FAST")){
if(missing(process)) {
throw("No process given!\n");
}
else if(!is.Process(process)){
throw("Process object invalid!\n");}
else if (missing(id)) {
throw("No site-process specific parameter id given!\n");
} else if (!is.character(id)) {
throw("Parameter id must be character!\n");
} else if (missing(value)){
throw("No new value given!\n");
}
}
if (missing(index)) {
index<-seq(along.with=this$.sites);
} else {
index<-.checkIndexSanity(this, index);
}
if(length(value) == 1) {
lapply(
this$.sites[index],
function(site){
setParameterAtSite(process,site,id,value);
}
);
} else {
counter<-1;
lapply(
this$.sites[index],
function(site){
if( counter > length(value) ){
counter<<-1;
}
setParameterAtSite(process,site,id,value[[counter]]);
counter<<-(counter+1);
}
);
}
return(invisible(this));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: setRateMultipliers
##
###########################################################################/**
#
# @RdocMethod setRateMultipliers
#
# @title "Set the values of the rate multiplier parameters for a given Process object and a collection of Site object aggregated by a Sequence object"
#
# \description{
# @get "title".
# This method just calls \code{setParameterAtSites(this=this,process=process,id="rate.multiplier",value=value,index=index)} See setParameterAtSites.Sequence for details.
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{process}{A valid Process object.}
# \item{value}{A numeric vector containing the new values of the site-process specific parameter, recycled if shorter than the index vector.}
# \item{index}{An integer vector specifying a set of positions. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# The Sequence object (invisible).
# }
#
#
# @author
#
# \seealso{
# setParameterAtSites.Sequence
# }
#
#*/###########################################################################
setMethodS3(
"setRateMultipliers",
class="Sequence",
function(
this,
process,
value,
index,
...
){
if(missing(process)){
throw("No process given!\n");
}
else if(missing(value)){
throw("No value provided!\n");
}
#else if(!is.GeneralSubstitution(process)){
# throw("The specified process is not a substitution process!\n");
#}
setParameterAtSites(this=this,process=process,id="rate.multiplier",value=value,index=index);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getRateMultipliers
##
###########################################################################/**
#
# @RdocMethod getRateMultipliers
#
# @title "Get the values of the rate multiplier parameters for a given Process object and a collection of Site object aggregated by a Sequence object"
#
# \description{
# @get "title".
# This method just calls \code{getParameterAtSites(this=this,process=process,id="rate.multiplier",index=index)} See getParameterAtSites.Sequence for details.
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{process}{A valid Process object.}
# \item{index}{An integer vector specifying a set of positions. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# A numeric vector with the current values of the rate multiplier in the specified range.
# }
#
#
# @author
#
# \seealso{
# setParameterAtSites.Sequence
# }
#
#*/###########################################################################
setMethodS3(
"getRateMultipliers",
class="Sequence",
function(
this,
process,
index,
...
){
if(missing(process)){
throw("No process given!\n");
}
rm<-getParameterAtSites(this=this,process=process,id="rate.multiplier",index=index);
return(as.numeric(lapply(rm,function(param){param$value})));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getParameterAtSites
##
###########################################################################/**
#
# @RdocMethod getParameterAtSites
#
# @title "Get the values of a site-process specific paramater for a process and a collection of Site objects aggregated by a Sequence object"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{process}{A valid Process object.}
# \item{id}{The identifier of the site-process specific parameter.}
# \item{index}{An integer vector specifying a set of positions. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# A list of site specific paramters. A site specific paramter is a list storing the id, the name, the value
# and the type of the parameter.
# }
#
# \examples{
# # create a sequence, attach a process
# p<-K80()
# s<-Sequence(length=6,alphabets=list(NucleotideAlphabet()),processes=list(list(p)))
# # set a new pattern of rate multipliers in the
# # range 1:3, the default value is 1.0 by the way
# setParameterAtSites(s,p,"rate.multiplier",c(2,3),1:3)
# # get rate multipliers
# getParameterAtSites(s,p,"rate.multiplier")
# # set a new value for the whole sequence
# setParameterAtSites(s,p,"rate.multiplier",0.5)
# # get rate multipliers
# getParameterAtSites(s,p,"rate.multiplier")
# }
#
# @author
#
# \seealso{
# Site Process @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"getParameterAtSites",
class="Sequence",
function(
this,
process,
id,
index,
...
){
if(!exists(x="PSIM_FAST")){
if(missing(process)) {
throw("No process given!\n");
}
else if(!is.Process(process)){
throw("Process object invalid!\n");}
else if (missing(id)) {
throw("No site-process specific parameter id given!\n");
} else if (!is.character(id)) {
throw("Parameter id must be character!\n");
}
}
if (missing(index)) {
index<-seq(along.with=this$.sites);
} else {
index<-.checkIndexSanity(this, index);
}
lapply(
this$.sites[index],
function(site){
getParameterAtSite(process,site,id);
}
);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getEvents
##
###########################################################################/**
#
# @RdocMethod getEvents
#
# @title "Get the list of active Event objects for a set of Site objects aggregated by a Sequence object"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{index}{An integer vector specifying a set of positions. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# A list of Event objects.
# }
#
# \examples{
# # create a sequence with a process attached
# s<-Sequence(
# string="ATGC",
# alphabets=list(NucleotideAlphabet()),
# processes=list(list(JC69()))
# )
# # get the active events from range 1:3
# getEvents(s,1:3)
# # get all active events via virtual field
# s$events
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"getEvents",
class="Sequence",
function(
this,
index,
...
){
if (missing(index)) {
index<-seq(along.with=this$.sites);
} else {
index<-.checkIndexSanity(this, index);
}
tmp<-list();
for (i in index){
# Setting the .positions field for then Events.
this$.sites[[i]]$.position<-i;
tmp<-c(tmp, getEvents(this$.sites[[i]]));
# Deleting the .position field;
this$.sites[[i]]$.position<-NULL;
}
tmp;
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getTotalRatesFromRange
##
###########################################################################/**
#
# @RdocMethod getTotalRatesFromRange
#
# @title "Get the vector of total site rates for a collection of Site objects aggregated by a Sequence object"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{index}{An integer vector specifying a set of positions. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# A numeric vector.
# }
#
# \examples{
# # create a sequence with some processes attached
# s<-Sequence(
# string="ATGC",
# alphabets=list(NucleotideAlphabet()),
# processes=list(list(JC69()),list(JC69(),GTR()))
# )
# # get total rates for positions 1 and 3
# getTotalRatesFromRange(s,c(1,3))
# # get all total rates via virtual field
# s$totalRates # via the "getTotalRates.Sequence" method
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"getTotalRatesFromRange",
class="Sequence",
function(
this,
index,
...
){
if (missing(index)) {
index<-seq(along.with=this$.sites);
} else {
index<-.checkIndexSanity(this, index);
}
if (this$.cumulative.rate.flag){
.recalculateCumulativeRates(this);
}
this$.total.rates[index];
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getTotalRates
##
###########################################################################/**
#
# @RdocMethod getTotalRates
#
# @title "Get the total site rates from a Sequence object"
#
# \description{
# @get "title".
# This method simply calls \code{getTotalRatesFromRange(this)}.
# See \code{getTotalRatesFromRange.Sequence} for more details.
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{...}{Not used.}
# }
#
# \value{
# A numeric vector containing the total site rates.
# }
#
# @author
#
# \seealso{
# getTotalRatesFromRange.Sequence
# }
#
#*/###########################################################################
setMethodS3(
"getTotalRates",
class="Sequence",
function(
this,
...
){
getTotalRatesFromRange(this);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: setTotalRates
##
###########################################################################/**
#
# @RdocMethod setTotalRates
#
# @title "Forbidden action: setting the list of total site rates for a Sequence 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(
"setTotalRates",
class="Sequence",
function(
this,
value,
...
){
virtualAssignmentForbidden(this);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getAncestral
##
###########################################################################/**
#
# @RdocMethod getAncestral
#
# @title "Get the ancestral object of a Sequence object"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{...}{Not used.}
# }
#
# \value{
# A Sequence object or a Process object.
# }
#
# \examples{
# # create a sequence object
# s<-Sequence(length=4)
# # get ancestral object
# getAncestral(s) # newly created sequences have the "Root insertion process" as ancestral
# # clone sequence
# cs<-clone(s)
# # get ancestral object id via virtual field
# cs$ancestral$id
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"getAncestral",
class="Sequence",
function(
this,
...
){
this$.ancestral.obj;
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getCumulativeRatesFromRange
##
###########################################################################/**
#
# @RdocMethod getCumulativeRatesFromRange
#
# @title "Get the cumulative site rates for a collection of Site objects aggregated by a Sequence object"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{index}{An integer vector specifying a set of positions. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# A numeric vector.
# }
#
# \examples{
# # create a sequence with some processes attached
# s<-Sequence(
# string="ATGC",
# alphabets=list(NucleotideAlphabet()),
# processes=list(list(JC69()),list(JC69(),GTR()))
# )
# # get cumulative rates for positions 1 and 3
# getCumulativeRatesFromRange(s,c(1,3))
# # get all cumulative rates via virtual field
# s$cumulativeRates # via the "getCumulativeRates.Sequence" method
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"getCumulativeRatesFromRange",
class="Sequence",
function(
this,
index,
...
){
if (missing(index)) {
index<-seq(along.with=this$.sites);
} else {
index<-.checkIndexSanity(this, index);
}
if (this$.cumulative.rate.flag){
.recalculateCumulativeRates(this);
}
this$.cumulative.rates[index];
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
setMethodS3(
".getCumulativeRatesFast",
class="Sequence",
function(
this,
...
){
if (this$.cumulative.rate.flag){
.recalculateCumulativeRates(this);
}
this$.cumulative.rates;
},
private=TRUE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: .recalculateCumulativeRates
##
setMethodS3(
".recalculateCumulativeRates",
class="Sequence",
function(
this,
target.site,
...
){
length<-this$.length;
if(length == 0){
return();
}
total.rates<-this$.total.rates;
sites<-this$.sites;
flagged.sites<-this$.flagged.sites;
if( length(flagged.sites) == 0 ) {
# Fresh start:
total.rates<-as.numeric(lapply(sites,getTotalRate));
} else {
# We have some flagged sites, recalculate just their total rates:
total.rates[flagged.sites]<-as.numeric(lapply(sites[flagged.sites],getTotalRate));
}
this$.total.rates<-total.rates;
this$.cumulative.rates<-cumsum(total.rates);
this$.flagged.sites<-integer(0);
this$.cumulative.rate.flag<-FALSE;
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getCumulativeRates
##
###########################################################################/**
#
# @RdocMethod getCumulativeRates
#
# @title "Get the total site rates from a Sequence object"
#
# \description{
# @get "title".
# This method simply calls \code{getCumulativeRatesFromRange(this)}.
# See \code{getCumulativeRates.Sequence} for more details.
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{...}{Not used.}
# }
#
# \value{
# A numeric vector containing the total site rates.
# }
#
# @author
#
# \seealso{
# getCumulativeRatesFromRange.Sequence
# }
#
#*/###########################################################################
setMethodS3(
"getCumulativeRates",
class="Sequence",
function(
this,
...
){
getCumulativeRatesFromRange(this);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: setCumulativeRates
##
###########################################################################/**
#
# @RdocMethod setCumulativeRates
#
# @title "Forbidden action: setting the cumulative rates for the sites aggregated by a Sequence 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(
"setCumulativeRates",
class="Sequence",
function(
this,
value,
...
){
virtualAssignmentForbidden(this);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getBigRate
##
###########################################################################/**
#
# @RdocMethod getBigRate
#
# @title "Get the sum of all active event rates from a Sequence object"
#
# \description{
# @get "title".
# The sum of active event rates depends on all Site object states and on the attached Process objects.
# It basically returns the last element of the cumulative site rates vector.
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{...}{Not used.}
# }
#
# \value{
# A numeric vector of length one.
# }
#
# \examples{
# # create a nucleotide sequence attach a process
# s<-NucleotideSequence(length=5);
# s$processes<-list(list(JC69()))
# # get the sum of active event rates
# getBigRate(s) # returns NA because site states are undefined
# # set site states
# s$states<-c("A","T")
# # get big rate via virtual field
# s$bigRate
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"getBigRate",
class="Sequence",
function(
this,
...
){
if (length(this$.sites) > 0) {
if(this$.cumulative.rate.flag) {
.recalculateCumulativeRates(this);
}
return(this$.cumulative.rates[this$.length]);
} else {
return(NA);
}
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: setBigRate
##
###########################################################################/**
#
# @RdocMethod setBigRate
#
# @title "Forbidden action: setting the sum of total active event rates for a Sequence 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(
"setBigRate",
class="Sequence",
function(
this,
value,
...
){
virtualAssignmentForbidden(this);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: setAncestral
##
###########################################################################/**
#
# @RdocMethod setAncestral
#
# @title "Set the ancestral object of a Sequence object"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{value}{A Sequence or a Process object.}
# \item{...}{Not used.}
# }
#
# \value{
# The new ancestral object (invisible).
# }
#
# \examples{
# # create a nucleotide sequence and a process object
# s<-NucleotideSequence(string="AGCT")
# p<-Process(name="MyProcess")
# # set the p as the ancestral of s
# setAncestral(s,p)
# s$ancestral
# # clone s
# cs<-clone(s)
# # set cs as ancestral of s via virtual field
# s$ancestral<-cs
# # get ancestral ids
# s$ancestral$id
# cs$ancestral$id
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"setAncestral",
class="Sequence",
function(
this,
value,
...
){
.checkWriteProtection(this);
if (!is.Sequence(value) & ! is.Process(value)) {
throw("Ancestral object must be a sequence or a process!\n");
} else {
this$.ancestral.obj<-value;
}
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: clone.Sequence
##
###########################################################################/**
#
# @RdocMethod clone
#
# @title "Clone a Sequence object"
#
# \description{
# @get "title".
#
# The cloning of Sequence objects involves the cloning of all aggregated Site objects. Because of that the
# cloning of long sequences is quite expensive.
# The cloned Site objects have the orginal Site objects as ancestral.
# The new Sequence objects has the original object as ancestral.
#
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{...}{Not used.}
# }
#
# \value{
# A Sequence object.
# }
#
# \examples{
# # cretate a nucleotide sequence
# s<-NucleotideSequence(string="ATG")
# # clone the sequence
# cs<-clone(s)
# # get some properties
# equals(s,s)
# equals(s,cs)
# cs$ancestral
# cs$sites[[1]]$ancestral
# }
#
# @author
#
# \seealso{
# Sequence clone.Object
# }
#
#*/###########################################################################
setMethodS3(
"clone",
class="Sequence",
function(
this,
...
){
# Cloning the whole sequence object:
that<-clone.Object(this)
# Disabling write protection:
if(that$writeProtected) {
that$writeProtected<-FALSE
}
# Setting the ancestral sequence:
that$.ancestral.obj<-this
# Resetting comments:
that$.comments<-list()
# Cloning sites:
clone.sites<-that$.sites
if(this$.length > 0) {
for (i in 1:this$.length) {
site<-this$.sites[[i]]
clone<-clone.Object(site)
clone$.ancestral<-site
clone$.sequence<-that
clone.sites[[i]]<-clone
}
}
that$.sites<-clone.sites
# Setting the name:
that$name<-paste("clone of",this$.name)
return(that)
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getWriteProtected
##
###########################################################################/**
#
# @RdocMethod getWriteProtected
#
# @title "Check if the object is write protected"
#
# \description{
# @get "title".
# Write protected objects cannot be modified through get/set methods and virtual fields.
# }
#
# @synopsis
#
# \arguments{
# \item{this}{An object.}
# \item{...}{Not used.}
# }
#
# \value{
# TRUE or FALSE
# }
#
# \examples{
#
# # create an object
# o<-Sequence()
# # toggle write protection
# o$writeProtected<-TRUE
# # check if it's write protected
# getWriteProtected(o)
# # check write protection via virtual field
# o$writeProtected
#
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"getWriteProtected",
class="Sequence",
function(
this,
...
){
this$.write.protected;
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: .checkIndexSanity
##
setMethodS3(
".checkIndexSanity",
class="Sequence",
function(
this,
index,
...
){
if (length(index) == 0 ) {
return(c());
}
if(!exists(x="PSIM_FAST")){
if( length(index) == 1 ) {
if(is.na(index)) {
warning("Index vector is NA! Coercing to empty vector!\n");
return(c());
}
if (is.nan(index)) {
warning("Index vector is NaN! Coercing to empty vector!\n");
return(c());
}
}
if(min(index) < 1 ) {
throw("Index vector element ",min(index)," too small!\n");
}
if( max(index) > this$.length ) {
throw("Index vector element ",max(index)," too big!\n");
}
}
return(index);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: setWriteProtected
##
###########################################################################/**
#
# @RdocMethod setWriteProtected
#
# @title "Set the write protection field for an object"
#
# \description{
# @get "title".
# Write protected objects cannot be modified through get/set methods and virtual fields.
# }
#
# @synopsis
#
# \arguments{
# \item{this}{An object.}
# \item{value}{A logical vector of size one.}
# \item{...}{Not used.}
# }
#
# \value{
# Invisible TRUE or FALSE.
# }
#
# \examples{
#
# # create an object
# o<-Sequence()
# # toggle write protection
# setWriteProtected(o,TRUE)
# # check write protection
# o$writeProtected
# # set write protection via virtual field
# o$writeProtected<-FALSE
# o$writeProtected
#
#
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"setWriteProtected",
class="Sequence",
function(
this,
value,
...
){
if(!is.logical(value)) {throw("The new value must be logical!\n")}
else {
this$.write.protected<-value;
}
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: .checkWriteProtection
##
setMethodS3(
".checkWriteProtection",
class="Sequence",
function(
this,
value,
...
){
if(getWriteProtected(this)) {throw("Cannot set value because the object is write protected!\n")}
else {return(FALSE)}
},
private=TRUE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: summary.Sequence
##
###########################################################################/**
#
# @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<-Sequence()
# # get a summary
# summary(a)
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"summary",
class="Sequence",
function(
object,
...
){
this<-object;
this$.summary$"Name"<-this$name;
this$.summary$"Id"<-this$id;
this$.summary$"Length"<-this$length;
this$.summary$"Big rate"<-this$bigRate;
this$.summary$"Ancestral object"<-this$ancestral$id;
if(this$.cumulative.rate.flag) {
this$.summary$"Cumulative rate flag"<-TRUE;
}
if(length(this$.flagged.sites) > 0 ) {
this$.summary$"Flagged sites"<-paste(this$.flagged.sites,collapse=" ");
}
if(this$writeProtected) {
this$.summary$"Write protected"<-TRUE;
}
NextMethod();
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: as.character.Sequence
##
###########################################################################/**
#
# @RdocMethod as.character
#
# @title "Get the string representation of a Sequence object"
#
# \description{
# @get "title".
# The string representation is the concatenation of the states of the
# aggregated Site object. Undefined states (NA-s) are represented by question marks.
# }
#
# @synopsis
#
# \arguments{
# \item{x}{A Sequence object.}
# \item{...}{Not used.}
# }
#
# \value{
# A character vector of length one.
# }
#
# \examples{
# # create object
# s<-Sequence(length=10)
# # get character representation
# as.character(s)
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"as.character",
class="Sequence",
function(
x,
...
){
getString(x);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: plot
##
###########################################################################/**
#
# @RdocMethod plot
#
# @title "Plot the total site rates for a collection of Site objects aggregated by a Sequence object"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{x}{A Sequence object.}
# \item{index}{An integer vector specifying a set of positions. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# Invisible TRUE or FALSE.
# }
#
# \examples{
# # create a nucleotide sequence with a process attached
# s<-NucleotideSequence(string="ATGGCCA",processes=list(list(JC69())))
# # plot total rates in range 1:4
# plot(s,1:4)
# # plot all total rates
# plot(s)
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"plot",
class="Sequence",
function(
x,
index=NA,
...
){
this<-x;
if(this$length == 0) {
warning("The sequence leght is zero, nothing to plot here!\n");
return(invisible(FALSE));
}
if( length(unique(this$totalRates)) == 1 & any(is.na(unique(this$totalRates))) ){
warning("The total rates are undefined, nothing to plot here!\n");
return(invisible(FALSE));
}
else {
if(missing(index)) {
index<-seq(along.with=1:this$length,by=1);
}
if(this$.cumulative.rate.flag){
.recalculateCumulativeRates(this);
}
what<-this$.total.rates[c(index)]
plot(
x=index,
y=what,
type="h",
lwd=1,
col="blue",
main=paste("Total rate plot for sequence", this$id),
xlab="Position",
ylab="Total rate",
ylim=c(0,max(what)),
xlim=c(min(index),max(index)),
xaxt="n"
);
axis(side=1, at=index, labels=index);
return(invisible(TRUE));
}
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: plotParameterAtSites
##
###########################################################################/**
#
# @RdocMethod plotParametersAtSites
#
# @title "Plot the value of a site-process specifc paramter for a collection of Site objects aggregated by a Sequence object"
#
# \description{
# @get "title".
# The type of the paramter must be numeric. The Process object must be attached to all positions specified
# in the index vector.
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{process}{A Process object.}
# \item{id}{The identifier of the site-process specific parameter.}
# \item{index}{An integer vector specifying a set of positions. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# Invisible TRUE or FALSE.
# }
#
# \examples{
# # create a nucleotide sequence with a process attached
# p<-JC69()
# s<-NucleotideSequence(string="ATGGCCA",processes=list(list(p)))
# # plot rate multipliers in range 1:4
# plotParametersAtSites(s,p,"rate.multiplier",1:4)
# # plot rate multiplier for the full sequence
# plotParametersAtSites(s,p,"rate.multiplier")
# }
#
# @author
#
# \seealso{
# Site Process Sequence
# }
#
#*/###########################################################################
setMethodS3(
"plotParametersAtSites",
class="Sequence",
function(
this,
process,
id,
index,
...
){
if(this$length == 0) {
warning("The sequence leght is zero, nothing to plot here!\n");
return(invisible(FALSE));
}
if(missing(index)) {
index<-seq(along.with=1:this$.length,by=1);
}
what<-apply(as.array(index),1,
function(pos){
tmp<-getParameterAtSites(this,process,id,pos)[[1]]$value;
if(!is.numeric(tmp)){
throw("Plot method failed becuase encountered non-numeric parameter value!\n");
}
return(tmp);
}
);
plot(
x=index,
y=what,
type="h",
lwd=1,
col="blue",
main=paste("Plot of parameter",id,"for process",process$id),
xlab="Position",
ylab="Value",
xlim=c(min(index),max(index)),
ylim=c(0,max(what)),
xaxt="n"
);
axis(side=1, at=index, labels=index);
invisible(TRUE);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: setDeletionTolerance
##
###########################################################################/**
#
# @RdocMethod setDeletionTolerance
#
# @title "Set the deletion tolerance site-process specific parameter for a collection of Site objects aggregated by a Sequence object"
#
# \description{
# @get "title".
# This method does some error checking and the calls \code{setParameterAtSites(this=this,process=process,id="deletion.tolerance",value=value,index=index)}.
# See \code{setParameterAtSites.Sequence} for more details.
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{process}{A Process object.}
# \item{value}{A numeric vector, recycled if shorter than the index vector.}
# \item{index}{An integer vector specifying a set of positions. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# The Sequence object (invisible).
# }
#
#
# @author
#
# \seealso{
# setParameterAtSites.Sequence
# }
#
#*/###########################################################################
setMethodS3(
"setDeletionTolerance",
class="Sequence",
function(
this,
process,
value,
index,
...
){
if(!exists(x="PSIM_FAST")){
if(missing(process)){
throw("No process given!\n");
}
else if (missing(value)){
throw("No value provided!\n");
}
else if(!inherits(process,"GeneralDeletor")){
throw("The specified process is not an insertion process!\n");
}
}
setParameterAtSites(this=this,process=process,id="deletion.tolerance",value=value,index=index);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getDeletionTolerance
##
###########################################################################/**
#
# @RdocMethod getDeletionTolerance
#
# @title "Get the deletion tolerance site-process specific parameter for a collection of Site objects aggregated by a Sequence object"
#
# \description{
# @get "title".
# This method does some error checking and the calls \code{getParameterAtSites(this=this,process=process,id="deletion.tolerance",index=index)}.
# See \code{getParameterAtSites.Sequence} for more details.
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{process}{A Process object.}
# \item{index}{An integer vector specifying a set of positions. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# A numeric vector.
# }
#
#
# @author
#
# \seealso{
# getParameterAtSites.Sequence
# }
#
#*/###########################################################################
setMethodS3(
"getDeletionTolerance",
class="Sequence",
function(
this,
process,
index,
...
){
if(!exists(x="PSIM_FAST")){
if(missing(process)){
throw("No process given!\n");
}
if(!inherits(process,"GeneralDeletor")){
throw("The specified process is not an insertion process!\n");
}
}
rm<-getParameterAtSites(this=this,process=process,id="deletion.tolerance",index=index);
return(as.numeric(lapply(rm,function(param){param$value})));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: setInsertionTolerance
##
###########################################################################/**
#
# @RdocMethod setInsertionTolerance
#
# @title "Set the insertion tolerance site-process specific parameter for a collection of Site objects aggregated by a Sequence object"
#
# \description{
# @get "title".
# This method does some error checking and the calls \code{setParameterAtSites(this=this,process=process,id="insertion.tolerance",value=value,index=index)}.
# See \code{setParameterAtSites.Sequence} for more details.
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{process}{A Process object.}
# \item{value}{A numeric vector, recycled if shorter than the index vector.}
# \item{index}{An integer vector specifying a set of positions. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# The Sequence object (invisible).
# }
#
#
# @author
#
# \seealso{
# setParameterAtSites.Sequence
# }
#
#*/###########################################################################
setMethodS3(
"setInsertionTolerance",
class="Sequence",
function(
this,
process,
value,
index,
...
){
if(!exists(x="PSIM_FAST")){
if(missing(process)){
throw("No process given!\n");
}
else if (missing(value)){
throw("No value provided!\n");
}
else if(!inherits(process,"GeneralInsertor")){
throw("The specified process is not an insertion process!\n");
}
}
setParameterAtSites(this=this,process=process,id="insertion.tolerance",value=value,index=index);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getInsertionTolerance
##
###########################################################################/**
#
# @RdocMethod getInsertionTolerance
#
# @title "Get the insertion tolerance site-process specific parameter for a collection of Site objects aggregated by a Sequence object"
#
# \description{
# @get "title".
# This method does some error checking and the calls \code{getParameterAtSites(this=this,process=process,id="insertion.tolerance",index=index)}.
# See \code{getParameterAtSites.Sequence} for more details.
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{process}{A Process object.}
# \item{index}{An integer vector specifying a set of positions. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# A numeric vector.
# }
#
#
# @author
#
# \seealso{
# getParameterAtSites.Sequence
# }
#
#*/###########################################################################
setMethodS3(
"getInsertionTolerance",
class="Sequence",
function(
this,
process,
index,
...
){
if(!exists(x="PSIM_FAST")){
if(missing(process)){
throw("No process given!\n");
}
if(!inherits(process,"GeneralInsertor")){
throw("The specified process is not an insertion process!\n");
}
}
rm<-getParameterAtSites(this=this,process=process,id="insertion.tolerance",index=index);
return(as.numeric(lapply(rm,function(param){param$value})));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: sampleStates
##
###########################################################################/**
#
# @RdocMethod sampleStates
#
# @title "Sample the states for a collection of Site objects aggregated by a Sequence object"
#
# \description{
# @get "title".
# This method samples new states from the equlibrium distribution of the attched process(es) for sites
# having undefined states (NA).
# If a site has more than one substitution process attached, then the method samples the new state from the
# mixture of equlibrium distributions. The weight of each equlibrium distribution is proportional to the
# site-process specific rate multiplier of the corresponding process at the given site.
#
# Sites having defined states are not touched. All sites with undefined states must have at least one
# substitution process (object inheriting from GeneralSubstitution) attached.
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{index}{An integer vector specifying a set of positions. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# The Sequence object (invisible).
# }
#
# \examples{
# # create a nucleotide sequence
# s<-NucleotideSequence(length=80)
# # create some processes
# jc69<-JC69(); # Jukes-Cantor
# hky<-HKY(base.freqs=c(0.1,0.4,0.1,0.4)) # "GC-rich" HKY
# # attach the processes
# s$processes<-list(list(jc69)) # jc is attached to all sites
# attachProcess(s,hky,60:80) # hky is attached to range 60:80
# # tweak rate multiplier for hky
# setRateMultipliers(s,hky,10,60:80)
# # set states in range 1:20
# setStates(s,"A",1:20)
# # sample remaining states
# sampleStates(s)
# # print sequence
# s
# }
#
# @author
#
# \seealso{
# sampleState.GeneralSubstitution GeneralSubstitution
# }
#
#*/###########################################################################
setMethodS3(
"sampleStates",
class="Sequence",
function(
this,
index,
...
){
if(!missing(index)){
index<-.checkIndexSanity(this, index);
}
else {
index<-seq(along.with=this$.sites);
}
for(site in this$.sites[index]){
# Sample states from the equlibrium distributions if
# the state is NA:
if(is.na(site$state)){
# Assemble the list of substitution processes:
subst.proc<-list();
for (proc in site$processes){
if(is.GeneralSubstitution(proc)){
subst.proc<-c(subst.proc, list(proc));
}
}
# Complain if we have no substitution processes to sample from:
if(length(subst.proc) == 0){
throw("Site state is NA and no substitution processes are attached. Cannot sample state!\n");
}
site.rates<-as.numeric(lapply(
subst.proc,
function(proc){
return(getParameterAtSite(proc, site,"rate.multiplier")$value);
}
));
# Normalizing site rates:
site.rates<-site.rates/sum(site.rates);
# Single subst process:
if(length(subst.proc) == 1){
site$state<-sampleState(subst.proc[[1]]);
}
else {
# Sample a substitution process according to the rate multipliers:
nproc<-sample(x=c(1:length(subst.proc)),size=1, replace=FALSE, prob=site.rates);
# Sample the state from the winner process:
site$state<-sampleState(subst.proc[[nproc]]);
}
} # if is.na...
}
return(invisible(this));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: clearStates
##
###########################################################################/**
#
# @RdocMethod clearStates
#
# @title "Set the states of a collection of Site objects aggregated by a Sequence object to undefined (NA)"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{index}{An integer vector specifying a set of positions. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# The Sequence object (invisible).
# }
#
# \examples{
# # create a nucleotide sequence
# s<-NucleotideSequence(string="ATGC")
# s
# # set states to NA in the range 2:3
# clearStates(s,2:3)
# s
# # set all states to NA
# clearStates(s)
# s
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"clearStates",
class="Sequence",
function(
this,
index,
...
){
if(!missing(index)){
index<-.checkIndexSanity(this, index);
}
else {
index<-seq(along.with=this$.sites);
}
for(site in this$.sites[index]){
site$.state<-NA;
site$.total.rate<-NA;
}
this$.cumulative.rate.flag<-TRUE;
return(invisible(this));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: setOmegas
##
setMethodS3(
"setOmegas",
class="Sequence",
function(
this,
process,
index,
...
){
# dummy method to force the creation of the generic function
setOmegas.CodonSequence(this,process,index,...);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getOmegas
##
setMethodS3(
"getOmegas",
class="Sequence",
function(
this,
process,
index,
...
){
# dummy method to force the creation of the generic function
getOmegas.CodonSequence(this,process,index,...);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: getSymbolFreqs
##
###########################################################################/**
#
# @RdocMethod getSymbolFreqs
#
# @title "Get a table with the frequencies of the states of a collection of Site objects aggregated by a Sequence object"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{index}{An integer vector specifying a set of positions. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# A table.
# }
#
# \examples{
# # create a nucleotide sequence
# s<-NucleotideSequence(length=30,processes=list(list(JC69())))
# # sample states
# sampleStates(s)
# # get state frequencies from ranges 1:10 and 20:30
# getSymbolFreqs(s,c(1:10,20:30))
# # get symbol frequencies for the full sequence
# getSymbolFreqs(s)
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"getSymbolFreqs",
class="Sequence",
function(
this,
index,
...
){
if(!missing(index)){
index<-.checkIndexSanity(this, index);
} else {
index<-seq(along.with=this$.sites);
}
prop.table(table(as.character(lapply(this$.sites[index],getState))));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: .discretizeGamma
##
setMethodS3(
".discretizeGamma",
class="Sequence",
function(
this,
shape,
ncat,
...
){
# figure out cutting points
cut<-apply(rbind(1:(ncat-1)/ncat),1, function(x){qgamma(x,shape=shape,rate=shape)})
cut<-c(0,cut,Inf);
# incomplete gamma function
Igamma<-function(x,a){
pgamma(x,shape=a, scale=1)
}
# function to calculate a category mean
cm<-function(a,b,shape,ncat){
( Igamma(b * shape, shape+1) - Igamma(a * shape,shape+1)) * ncat;
}
# calculate category means
means<-c();
for (i in 1:(length(cut)-1)){
means<-c(means,cm(cut[i], cut[i+1], shape, ncat));
}
# return a vector with the means
return(means);
},
private=TRUE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: plusGamma
##
###########################################################################/**
#
# @RdocMethod plusGamma
#
# @title "Sample the rate multiplier parameters of a Process from a Gamma distribution for a collection of Site objects aggregated by a Sequence object"
#
# \description{
# @get "title".
#
# The gamma distribution is discretized by calculating the means of the portions corresponding
# to the categories having equal probabilities. If the \code{ncat} argument is not numeric, the
# rates are sampled from the continuous gamma distribution.
# }
#
# \references{
# Yang, Z. (1994) Maximum likelihood phylogenetic estimation from DNA sequences with variable
# rates over sites: approximate methods - Journal of Molecular Evolution 39:306-314 \url{http://dx.doi.org/10.1007/BF00160154}
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{process}{A Process object.}
# \item{shape}{The shape parameter of the gamma distribution.}
# \item{index}{A vector of positions.}
# \item{ncat}{Numer of categories in the discretized gamma distribution (4 by default).}
# \item{...}{Not used.}
# }
#
# \value{
# The Sequence object (invisible);
# }
#
# \examples{
# # create a sequence
# s<-NucleotideSequence(length=20)
# # attach a process
# p<-JC69()
# attachProcess(s,p)
# # get rate multipliers
# getRateMultipliers(s,p) # the default value is 1.0
# # sample rate multipliers in range 1:5 from a discrete
# #gamma distribution with shape parameter 0.5
# plusGamma(s,p,0.5,1:5)
# # get rate multipliers
# getRateMultipliers(s,p) # the default value is 1.0
# # sample rates from a continuous gamma distribution
# # with shape parameter 0.5
# plusGamma(s,p,0.5,ncat="cont")
# # get rate multipliers
# getRateMultipliers(s,p)
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"plusGamma",
class="Sequence",
function(
this,
process,
shape,
index,
ncat=4,
...
){
if(missing(process)){
throw("No process specified!\n");
}
else if(missing(shape)){
throw("No shape parameter specified!\n");
}
else if(!all(is.numeric(shape)) | length(shape) != 1){
throw("The shape parameter must be a numeric vector of lenght 1!\n");
}
if(missing(index)){
index<-seq(along.with=this$.sites);
}
else {
index<-.checkIndexSanity(this, index);
}
if(!is.numeric(ncat)){
# continuous gamma
setParameterAtSites(this, process=process, id="rate.multiplier",value=rgamma(length(index),shape=shape,rate=shape),index=index);
}
else{
setParameterAtSites(this, process=process, id="rate.multiplier",sample(.discretizeGamma(this,shape,ncat),size=length(index),replace=TRUE),index=index);
}
return(invisible(this));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: plusInvGamma
##
###########################################################################/**
#
# @RdocMethod plusInvGamma
#
# @title "Sample the rate multiplier parameters of a Process from an I+Gamma distribution for a collection of Site objects aggregated by a Sequence object"
#
# \description{
# @get "title".
#
# The gamma distribution is discretized by calculating the means of the portions corresponding
# to the categories having equal probabilities. If the \code{ncat} argument is not numeric, the
# rates are sampled from the continuous gamma distribution.
# }
#
# \references{
# Gu X, Fu, YX, Li, WH (1995) Maximum likelihood estimation of the heterogeneity of substitution
# rate among nucleotide sites - Mol. Biol. Evol. 12(4):546-57 \url{http://bit.ly/aE6xF0}
#
# Yang, Z (1994) Maximum likelihood phylogenetic estimation from DNA sequences with variable
# rates over sites: approximate methods - Journal of Molecular Evolution 39:306-314 \url{http://dx.doi.org/10.1007/BF00160154}
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{process}{A Process object.}
# \item{pinv}{The proportion of invariant sites.}
# \item{shape}{The shape parameter of the gamma distribution.}
# \item{index}{A vector of positions.}
# \item{ncat}{Numer of categories in the discretized gamma distribution (4 by default).}
# \item{...}{Not used.}
# }
#
# \value{
# The Sequence object (invisible);
# }
#
# \examples{
# # create a sequence
# s<-NucleotideSequence(length=20)
# # attach a process
# p<-JC69()
# attachProcess(s,p)
# # get rate multipliers
# getRateMultipliers(s,p) # the default value is 1.0
# # sample rate multipliers in range 1:5 from I+discrete Gamma
# plusInvGamma(s,p,pinv=0.5,shape=0.5,1:5)
# # get rate multipliers
# getRateMultipliers(s,p) # the default value is 1.0
# # sample rates from an I+continuos Gamma model
# plusInvGamma(s,p,pinv=0.5,shape=0.5,1:5,ncat="cont")
# # get rate multipliers
# getRateMultipliers(s,p)
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"plusInvGamma",
class="Sequence",
function(
this,
process,
pinv,
shape,
index,
ncat=4,
...
){
if(missing(process)){
throw("No process specified!\n");
}
else if(missing(pinv)){
throw("No proportion of invariant sites given!\n");
}
else if(!all(is.numeric(pinv)) | length(pinv) != 1){
throw("The pinv parameter must be a numeric vector of lenght 1!\n");
}
else if(pinv > 1){
throw("Tpe proportion of invariant sites cannot be larger than 1.!");
}
else if(missing(shape)){
throw("No shape parameter specified!\n");
}
else if(!all(is.numeric(shape)) | length(shape) != 1){
throw("The shape parameter must be a numeric vector of lenght 1!\n");
}
else {
if(missing(index)){
index<-seq(along.with=this$.sites);
}
else {
index<-.checkIndexSanity(this, index);
}
# discretize gamma distribution
dg<-c()
if(is.numeric(ncat)){
dg<-.discretizeGamma(this,shape,ncat);
}
# Iterating over the sites specified by the index vector:
for(site in index){
# Choose between invariant and gamma:
type<-sample(c("INV","GAMMA"),size=1, replace=FALSE, prob=c( pinv, (1-pinv) ) );
if(type == "INV"){
setParameterAtSites(this, process=process, id="rate.multiplier",value=0,index=c(site));
}
else {
if(!is.numeric(ncat)){
# continuous gamma
setParameterAtSites(this, process=process, id="rate.multiplier",value=rgamma(1,shape=shape,rate=shape),index=c(site));
}
else{
setParameterAtSites(this, process=process, id="rate.multiplier",sample(dg,size=1),index=c(site));
}
}
}
}
return(invisible(this));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: insertSequence
##
###########################################################################/**
#
# @RdocMethod insertSequence
#
# @title "Insert a Sequence object into another Sequence object after a specified position"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{The target Sequence object.}
# \item{insert}{The Sequence object to be inserted.}
# \item{position}{The position after the Sequence object will be inserted.}
# \item{process}{The Process object performing the insertion (optional).}
# \item{paranoid}{If TRUE, then the consistency of teh target objects is checked more rigurously after insertion.}
# \item{...}{Not used.}
# }
#
# \value{
# The Sequence object (invisible).
# }
#
# \examples{
# # create some sequence objects.
# target<-NucleotideSequence(string="AAAAAAAAAAA")
# insert<-NucleotideSequence(string="GGGGGGGGGGG")
# # insert after position 5
# insertSequence(target,insert,5)
# # print the target sequence
# target
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"insertSequence",
class="Sequence",
function(
this,
insert,
position,
process=NA,
paranoid=FALSE,
...
){
.checkWriteProtection(this);
if(!exists(x="PSIM_FAST")){
if(missing(insert)) {
throw("Insert sequence object is missing!\n");
}
else if (missing(position)) {
throw("Insertion position is missing!\n");
}
if(!is.Sequence(insert)) {
throw("Insert object not valid!\n");
}
else if (this$length == 0 & position != 0 ) {
throw("Acceptor sequence length is zero! The only valid insertion position is 0!\n");
}
else if ( !( position >= 0 & position <=(this$.length + 1))) {
throw("Insertion position ",position," is invalid!\n");
}
}
# Just return if insert has zero length:
if(insert$length == 0){
warning("The length of the sequence to be inserted is zero! Nothing to do here!\n");
return(invisible(FALSE));
}
# Clone insert object:
insert<-clone(insert);
# Set the generator process:
if(!missing(process)) {
if( (length(process) == 0) | !is.Process(process)){
throw("Process object invalid!\n");
}
} else {
process<-Sequence$.root.ins;
}
for(site in insert$.sites) {
site$.ancestral<-process;
site$.sequence<-this;
}
# Recalculate cumulative rates if the flag is on:
if(this$.cumulative.rate.flag) {
.recalculateCumulativeRates(this);
}
# Flagging cumulative rates:
this$.cumulative.rate.flag<-TRUE;
# Inserting new site objects:
if ( position == this$.length) {
# Insertion at the end of the sequence;
this$.sites<-c(this$.sites,insert$.sites);
this$.total.rates<-c(this$.total.rates,rep(c(NA),times=insert$.length) );
this$.cumulative.rates<-c(this$.cumulative.rates,rep(NA,times=insert$.length) );
} else if (position == 0) {
# Insertion in the sequence
this$.sites<-c(insert$.sites, this$.sites);
this$.total.rates<-c(rep(NA,times=insert$.length),this$.total.rates);
this$.cumulative.rates<-c(rep(NA,times=insert$.length),this$.cumulative.rates);
} else {
# Insertion at position 0
this$.sites<-c(this$.sites[1:position],insert$.sites,this$.sites[(position+1):this$.length]);
this$.total.rates<-c(this$.total.rates[1:position],rep(NA,times=insert$.length),this$.total.rates[(position+1):this$.length]);
this$.cumulative.rates<-c(this$.cumulative.rates[1:position],rep(NA,times=insert$.length),this$.cumulative.rates[(position+1):this$.length]);
}
# Checking if lengths are consistent:
if(!exists(x="PSIM_FAST")){
if(length(this$.sites) != (this$.length + insert$.length)) {
throw("Length inconsistency after insertion!\n");
}
}
# Setting new length:
this$.length<-(this$.length + insert$.length);
# Flagging the inserted sites:
this$.flagged.sites<-c(this$.flagged.sites,(position+1):(position+insert$.length));
if(!exists(x="PSIM_FAST")){
if(length(this$.total.rates) != this$.length) {
throw("Total rates vector inconsistency after insertion!\n");
}
if(length(this$.cumulative.rates) != this$.length) {
throw("Cumulative rates vector inconsistency after insertion!\n");
}
}
# Recalculating cumulative rates:
.recalculateCumulativeRates(this);
# Paranoid check of total rates:
if(paranoid) {
for (i in 1:this$.length) {
if(this$.sites[[i]]$totalRate != this$.total.rates[[i]]) {
throw("Object total rates inconsistent with total rates vector!\n");
}
}
}
# Deleting the insert:
rm(insert);
return(invisible(this));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: deleteSubSequence
##
###########################################################################/**
#
# @RdocMethod deleteSubSequence
#
# @title "Delete a collection of sites aggregated by a Sequence object"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{index}{An index vector specifying a collection of sites to be deleted. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# Invisible TRUE or errror.
# }
#
# \examples{
# # create a nucleotide sequence
# s<-NucleotideSequence(string="ATATATATATATATAT")
# # delete sites 2, 4 and 6
# deleteSubSequence(s,c(2,4,6))
# s
# # delete sites in the range 3:6
# deleteSubSequence(s,3:6)
# s
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"deleteSubSequence",
class="Sequence",
function(
this,
index,
...
){
if(length(index) == 0) {
return(FALSE);
}
.checkWriteProtection(this);
if(!exists(x="PSIM_FAST")){
if(missing(index)) {
throw("No index vector specified!\n");
}
if(length(index) == 0) {
return(FALSE);
}
}
index<-.checkIndexSanity(this, index);
# Avoid deletion on dirty sequence as
# that will cause havoc.
if(this$.cumulative.rate.flag){
.recalculateCumulativeRates(this);
}
# Flagging cumulative rates:
this$.cumulative.rate.flag<-TRUE;
min.index<-min(index);
# Deleting site objects:
this$.sites[index]<-NULL;
# Updating rate vectors:
this$.total.rates<-this$.total.rates[-index];
this$.cumulative.rates<-this$.cumulative.rates[-index];
# Flag the site before the deletion to
# to force cumulative rate recalculation:
if (min.index > 2 ) {
this$.flagged.sites<-c(this$.flagged.sites,(min.index - 1));
}
if(!exists(x="PSIM_FAST")){
if( length(this$.sites) != (this$.length - length(index) ) ) {
throw("Inconsistency after deleting sites!\n");
}
}
this$.length<-length(this$.sites);
.recalculateCumulativeRates(this);
return(invisible(TRUE));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: copySubSequence
##
###########################################################################/**
#
# @RdocMethod copySubSequence
#
# @title "Copy a collection of Site objects aggregated by a Sequence object into a new Sequence object"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A Sequence object.}
# \item{index}{An index vector specifying a collection of sites to be copied. It is set to 1:seq$length if omitted.}
# \item{process}{The Process object performing the copy (optional).}
# \item{...}{Not used.}
# }
#
# \value{
# A new Sequence object.
# }
#
# \examples{
# # create a nucleotide sequence
# s<-NucleotideSequence(string="ATATATATATATATATA")
# # copy sites in the range 3:8 in a new object
# s2<-copySubSequence(s,3:8)
# s2
# # copy sites 1,3 and 5 from s2
# s3<-copySubSequence(s2,c(1,3,5))
# s3
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"copySubSequence",
class="Sequence",
function(
this,
index,
process=NA,
...
){
if(missing(index)) {
index<-seq(along.with=this$.sites);
}
else {
index<-.checkIndexSanity(this, index);
}
if(!exists(x="PSIM_FAST")){
if(!is.na(process) & !is.Process(process)) {
throw("Process object invalid!\n");
}
}
# Avoid copying from dirty sequence:
if(this$.cumulative.rate.flag){
.recalculateCumulativeRates(this);
}
length<-length(index);
# Create an empty sequence object:
copy<-Sequence();
# Flag copy cumulative rates:
copy$.cumulative.rate.flag<-TRUE;
if(is.na(process)){
# Getting the root insertion process:
process<-Sequence$.root.ins;
}
# Setting the ancestral to sequence:
copy$.ancestral.obj<-process;
# Setting copy name:
copy$name<-paste("Copied from",this$name);
# Setting length:
copy$.length<-length;
# Clone the sites:
copy$.sites<-lapply(this$.sites[index],
function(site){
site.copy<-clone(site);
site.copy$.ancestral<-process;
return(site.copy);
}
);
# Copy total rates:
copy$.total.rates<-this$.total.rates[index];
# Create cumulative rates vector:
copy$.cumulative.rates<-cumsum(copy$.total.rates);
copy$.cumulative.rate.flag<-FALSE;
if(!exists(x="PSIM_FAST")){
if(length(copy$.sites) != length){
throw("Sites list length mismatch!\n")
}
else if(length(copy$.total.rates) != length){
throw("Total rates vector length mismatch!\n")
}
else if(length(copy$.cumulative.rates) != length){
throw("Cumulative rates vector length mismatch!\n")
}
}
return(copy);
},
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.