Nothing
R/PopularSequences.R
In phylosim: Flexible Simulations of Biological Sequence Evolution
##
## Copyright 2009 Botond Sipos
## See the package description for licensing information.
##
##
## BinarySequence
##
##########################################################################/**
#
# @RdocClass BinarySequence
#
# @title "The BinarySequence class"
#
# \description{
# Sequence objects aggregating Site objects having a BinaryAlphabet attached by default.
#
# @classhierarchy
# }
#
# @synopsis
#
# \arguments{
# \item{name}{Name of the Sequence object.}
# \item{string}{A string specifying the length and the states of the Sequence object.}
# \item{length}{The length of the sequence. Mutually exclusive with "string".}
# \item{processes}{A list of lists of Process objects, to be attached to the aggregated Site objects. Recycled if shorter than the length of the sequence.}
# \item{ancestral.obj}{The ancestral object of the Sequence object (a valid Sequence or Process object).}
# \item{...}{Not used.}
# }
#
# \section{Fields and Methods}{
# @allmethods
# }
#
# \examples{
# # create an empty BinarySequence object
# s<-BinarySequence(length=50)
# s
# # set states
# s$states<-c(0,0,1,0,1,1)
# s
# # create a sequence object by specifying a string
# s<-BinarySequence(string="00000110010001111")
# s
# }
#
# @author
#
# \seealso{
# Sequence BinaryAlphabet
# }
#
#*/###########################################################################
setConstructorS3(
"BinarySequence",
function(
name=NA,
string=NA,
length=NA,
processes=NA,
ancestral.obj=NA,
...
){
if(!missing(string)){
this<-Sequence(string=string,alphabets=list(BinaryAlphabet()));
}
else if(!missing(length)) {
this<-Sequence(length=length,alphabets=list(BinaryAlphabet()));
}
else {
this<-Sequence(alphabets=list(BinaryAlphabet()));
}
this<-extend(this, "BinarySequence");
if(!missing(name)){
this$name<-name;
}
if(!missing(processes)){
setProcesses(this,processes);
}
if(!missing(ancestral.obj)){
this$.ancestral<-ancestral.obj;
}
return(this);
},
enforceRCC=TRUE
);
##
## NucleotideSequence
##
##########################################################################/**
#
# @RdocClass NucleotideSequence
#
# @title "The NucleotideSequence class"
#
# \description{
# Sequence objects aggregating Site objects having a NucleotideAlphabet attached by default.
#
# @classhierarchy
# }
#
# @synopsis
#
# \arguments{
# \item{name}{Name of the Sequence object.}
# \item{string}{A string specifying the length and the states of the Sequence object.}
# \item{length}{The length of the sequence. Mutually exclusive with "string".}
# \item{processes}{A list of lists of Process objects, to be attached to the aggregated Site objects. Recycled if shorter than the length of the sequence.}
# \item{ancestral.obj}{The ancestral object of the Sequence object (a valid Sequence or Process object).}
# \item{...}{Not used.}
# }
#
# \section{Fields and Methods}{
# @allmethods
# }
#
# \examples{
# # create an empty NucleotideSequence object
# s<-NucleotideSequence(length=50)
# s
# # set states
# s$states<-c("A","A","G","T")
# s
# # create a sequence object by specifying a string
# s<-NucleotideSequence(string="ATGCCGATTAGCAAA")
# s
# }
#
# @author
#
# \seealso{
# Sequence NucleotideAlphabet
# }
#
#*/###########################################################################
setConstructorS3(
"NucleotideSequence",
function(
name=NA,
string=NA,
length=NA,
processes=NA,
ancestral.obj=NA,
...
){
if(!missing(string)){
this<-Sequence(string=string,alphabets=list(NucleotideAlphabet()));
}
else if(!missing(length)) {
this<-Sequence(length=length,alphabets=list(NucleotideAlphabet()));
}
else {
this<-Sequence(alphabets=list(NucleotideAlphabet()));
}
this<-extend(this, "NucleotideSequence");
if(!missing(name)){
this$name<-name;
}
if(!missing(processes)){
setProcesses(this,processes);
}
if(!missing(ancestral.obj)){
this$.ancestral<-ancestral.obj;
}
return(this);
},
enforceRCC=TRUE
);
##
## revComp
##
###########################################################################/**
#
# @RdocMethod revComp
#
# @title "Reverse complmenet a NucleotideSequence object"
#
# \description{
# @get "title".
#
# The method reverse complements the sequence "in place", no object cloning is performed.
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A NucleotideSequence object}
# \item{...}{Not used.}
# }
#
# \value{
# Invisible TRUE.
# }
#
# \examples{
# s<-NucleotideSequence(string="ATGC")
# s
# revComp(s)
# s
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"revComp",
class="NucleotideSequence",
function(
this,
...
){
states<-rev(as.character(getStates(this)));
len<-length(states);
if(len == 0){return(invisible(TRUE))}
for(i in 1:len){
if(states[[i]] == "NA"){
state<-NA;
}
else if(states[i] == "A"){
state<-"T";
}
else if(states[i] == "T"){
state<-"A";
}
else if(states[i] == "G"){
state<-"C";
}
else if(states[i] == "C"){
state<-"G";
} else {
throw("Symbol not in NucleotideAlphabet!");
}
this$.sites[[i]]$.state<-state;
}
return(invisible(TRUE));
},
private=TRUE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## AminoAcidSequence
##
##########################################################################/**
#
# @RdocClass AminoAcidSequence
#
# @title "The AminoAcidSequence class"
#
# \description{
# Sequence objects aggregating Site objects having an AminoAcidAlphabet attached by default.
#
# @classhierarchy
# }
#
# @synopsis
#
# \arguments{
# \item{name}{Name of the Sequence object.}
# \item{string}{A string specifying the length and the states of the Sequence object.}
# \item{length}{The length of the sequence. Mutually exclusive with "string".}
# \item{processes}{A list of lists of Process objects, to be attached to the aggregated Site objects. Recycled if shorter than the length of the sequence.}
# \item{ancestral.obj}{The ancestral object of the Sequence object (a valid Sequence or Process object).}
# \item{...}{Not used.}
# }
#
# \section{Fields and Methods}{
# @allmethods
# }
#
# \examples{
# # create an empty AminoAcidSequence object
# s<-AminoAcidSequence(length=50)
# s
# # set states
# s$states<-c("C","C","G","Y")
# s
# # create a sequence object by specifying a string
# s<-AminoAcidSequence(string="CNGGYCCNGYYYY")
# s
# }
#
# @author
#
# \seealso{
# Sequence AminoAcidAlphabet
# }
#
#*/###########################################################################
setConstructorS3(
"AminoAcidSequence",
function(
name=NA,
string=NA,
length=NA,
processes=NA,
ancestral.obj=NA,
...
){
if(!missing(string)){
this<-Sequence(string=string,alphabets=list(AminoAcidAlphabet()));
}
else if(!missing(length)) {
this<-Sequence(length=length,alphabets=list(AminoAcidAlphabet()));
}
else {
this<-Sequence(alphabets=list(AminoAcidAlphabet()));
}
this<-extend(this, "AminoAcidSequence");
if(!missing(name)){
this$name<-name;
}
if(!missing(processes)){
setProcesses(this,processes);
}
if(!missing(ancestral.obj)){
this$.ancestral<-ancestral.obj;
}
return(this);
},
enforceRCC=TRUE
);
##
## CodonSequence
##
##########################################################################/**
#
# @RdocClass CodonSequence
#
# @title "The CodonSequence class"
#
# \description{
# Sequence objects aggregating Site objects having a CodonAlphabet attached by default.
#
# @classhierarchy
# }
#
# @synopsis
#
# \arguments{
# \item{name}{Name of the Sequence object.}
# \item{string}{A string specifying the length and the states of the Sequence object.}
# \item{length}{The length of the sequence. Mutually exclusive with "string".}
# \item{table.id}{The genetic code table to use in the attached CodonAlphabet object ("Standard" by default).}
# \item{processes}{A list of lists of Process objects, to be attached to the aggregated Site objects. Recycled if shorter than the length of the sequence.}
# \item{ancestral.obj}{The ancestral object of the Sequence object (a valid Sequence or Process object).}
# \item{...}{Not used.}
# }
#
# \section{Fields and Methods}{
# @allmethods
# }
#
# \examples{
# # create an empty CodonSequence object
# s<-CodonSequence(length=50)
# s
# # set states
# s$states<-c("ATG","CGA","TTT","CTA")
# s
# # create a sequence object by specifying a string
# s<-CodonSequence(string="ATCTTTCGAATG")
# s
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setConstructorS3(
"CodonSequence",
function(
name=NA,
string=NA,
length=NA,
table.id=1,
processes=NA,
ancestral.obj=NA,
...
){
if (!missing(string)){
len<-((stringLength(string)/3));
flen<-floor(len);
this<-Sequence(length=flen,alphabets=list(CodonAlphabet(table.id=table.id)));
if(len > flen){
warning("The length of the provided string was not multiple of 3. The incomplete codon was discarded!\n");
}
nuc<-strsplit(string,"",fixed=TRUE)[[1]];
j<-1; # counter for the codon position
for(i in seq(from=1,to=(flen * 3),by=3)){
# get the codon:
state<-paste(nuc[c(i,i+1,i+2)],collapse="");
# Check for stop codons:
if( isStopCodon(this$.sites[[j]]$.alphabet, state) ){
throw("The CodonSequence objects does not accept stop codons as valid states!\n");
}
# Set the state:
setStates(this,state,j);
j<-j+1;
}
}
else if(!missing(length)){
this<-Sequence(length=length,alphabets=list(CodonAlphabet(table.id=table.id)));
}
else {
this<-Sequence(alphabets=list(CodonAlphabet(table.id=table.id)));
}
this<-extend(this, "CodonSequence");
if(!missing(name)){
this$name<-name;
}
if(!missing(processes)){
setProcesses(this,processes);
}
if(!missing(ancestral.obj)){
this$.ancestral<-ancestral.obj;
}
return(this);
},
enforceRCC=TRUE
);
##
## Method: getOmegas
##
###########################################################################/**
#
# @RdocMethod getOmegas
# \alias{getOmegas.Sequence}
#
# @title "Get the omegas from a collection of sites"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A CodonSequence object.}
# \item{process}{A process object inheriting from GY94.}
# \item{index}{A vector of positions.}
# \item{...}{Not used.}
# }
#
# \value{
# A numeric vector.
# }
#
# \examples{
# # create a GY94 process
# p<-GY94()
# # create a CodonSequence object,
# # attach a process p
# s<-CodonSequence(length=20,processes=list(list(p)))
# # set omega values in range 1:5
# setOmegas(s,p,c(0.5,1,1.5),1:5)
# # get omega values from siutes 1,2,3,10, and 20
# getOmegas(s,p,c(1:3,10,20))
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"getOmegas",
class="CodonSequence",
function(
this,
process,
index,
...
){
if(missing(process)){
throw("No process given!\n");
}
else if(!is.GY94(process)){
throw("The specified process is not a GY94 codon substitution process!\n");
}
rm<-getParameterAtSites(this=this,process=process,id="omega",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: omegaHist
##
###########################################################################/**
#
# @RdocMethod omegaHist
#
# @title "Plot a histogram of omega values from a range"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A CodonSequence object.}
# \item{process}{A process object inheriting from GY94.}
# \item{breaks}{\code{breaks} parameter for \code{hist()}.}
# \item{index}{A vector of positions.}
# \item{...}{Not used.}
# }
#
# \value{
# The CodonSequence object (invisible).
# }
#
# \examples{
# # create a GY94 process
# p<-GY94()
# # create a CodonSequence object,
# # attach a process p
# s<-CodonSequence(length=20,processes=list(list(p)))
# # set omega values through omegaVarM2.CodonSequence
# omegaVarM2(s,p,p0=0.5,p1=0.2,omega=1.5)
# # get a histogram of omega values from the range 1:15
# omegaHist(s,p,breaks=10,1:15)
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"omegaHist",
class="CodonSequence",
function(
this,
process,
breaks,
index,
...
){
if(missing(process)){
throw("No process given!\n");
}
else if(!is.GY94(process)){
throw("The specified process is not a GY94 codon substitution process!\n");
}
if(missing(index)){
index<-seq(along.with=this$.sites);
}
if(missing(breaks)){
hist(getOmegas(this,process,index));
}
else {
omegas<-getOmegas(this,process,index);
hist(omegas,breaks=breaks,main="Histogram of omega values",xlab="Omega",freq=FALSE);
}
return(invisible(this));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: setOmegas
###########################################################################/**
#
# @RdocMethod setOmegas
# \alias{setOmegas.Sequence}
#
# @title "Set the omegas for a collection of sites"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A CodonSequence object.}
# \item{process}{A process object inheriting from GY94.}
# \item{value}{A vector containing the new values of the site-process specific parameter, recycled if shorter than the index vector.}
# \item{index}{A vector of positions. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# A numeric vector.
# }
#
# \examples{
# # create a GY94 process
# p<-GY94()
# # create a CodonSequence object,
# # attach a process p
# s<-CodonSequence(length=20,processes=list(list(p)))
# # set omega values in range 1:5
# setOmegas(s,p,c(0.5,1,1.5),1:5)
# # get omega values from siutes 1,2,3,10, and 20
# getOmegas(s,p,c(1:3,10,20))
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"setOmegas",
class="CodonSequence",
function(
this,
process,
value,
index,
...
){
if(missing(process)){
throw("No process specified!\n");
}
if(!is.GY94(process)){
throw("The sepcified process is not a GY94 codon substitution process!\n");
}
else if(missing(value)){
throw("No new value specified!\n");
}
else if(!all(is.numeric(value)) ){
throw("The new value must be a numeric vector!\n");
}
else {
if(missing(index)){
index<-seq(along.with=this$.sites);
}
else {
index<-.checkIndexSanity(this, index);
}
setParameterAtSites(this, process=process, id="omega",value=value,index=index);
return(invisible(TRUE));
}
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: omegaVarM0 - one ratio
##
###########################################################################/**
#
# @RdocMethod omegaVarM0
#
# @title "The M0 (one-ratio) model of variable omega ratios among sites"
#
# \description{
# @get "title".
#
# This method sets the \code{omega} site-process specific parameter
# in the specified range to values sampled from the M0 (one-ratio) model of
# variable omega ratios among sites.
#
# Distribution of omega values:
# \preformatted{
# CATEGORY PROBABILITY
#
# omega 1
# }
# }
#
# \references{
# Yang, Z., Nielsen, R., Goldman, N., Pedersen Krabbe, A-M. (2000) Codon-Substitution Models for Heterogeneous Selection Pressure at Amino Acid Sites - Genetics 155:431-449 \url{http://bit.ly/bvjucn}
#
# Goldman, N., Yang, Z. (1994) A codon-based model of nucleotide substitution for protein-coding DNA sequences - Mol Biol Evol 11(5):725-36 \url{http://bit.ly/aSVEoa}
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A CodonSequence object.}
# \item{process}{A process object inheriting from GY94.}
# \item{omega}{The fixed omega value.}
# \item{index}{A vector of positions.}
# \item{...}{Not used.}
# }
#
# \value{
# Invisible TRUE.
# }
#
# \examples{
# # create a GY94 object
# p<-GY94(kappa=2)
# # create a CodonSequence object, attach process p
# s<-CodonSequence(length=20, processes=list(list(p)))
# # sample states
# sampleStates(s)
# # sample omegas in range 1:5 from model M0
# omegaVarM0(s,p,omega=2,1:5)
# # get omega values
# getOmegas(s,p)
# # get a histogram of omega values in range 1:5
# omegaHist(s,p,breaks=50,1:5)
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"omegaVarM0",
class="CodonSequence",
function(
this,
process,
omega,
index,
...
){
if(missing(process)){
throw("No process specified!\n");
}
if(!is.GY94(process)){
throw("The sepcified process is not a GY94 codon substitution process!\n");
}
else if(missing(omega)){
throw("No new omega value specified!\n");
}
else if((!is.numeric(omega))| (length(omega) != 1)){
throw("The new value must be a numeric vector of length 1!\n");
}
else {
if(missing(index)){
index<-seq(along.with=this$.sites);
}
else {
index<-.checkIndexSanity(this, index);
}
setParameterAtSites(this, process=process, id="omega",value=omega,index=index);
return(invisible(TRUE));
}
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: omegaVarM1 - neutral
##
###########################################################################/**
#
# @RdocMethod omegaVarM1
#
# @title "The M1 (neutral) model of variable omega ratios among sites"
#
# \description{
# @get "title".
#
# This method sets the \code{omega} site-process specific parameter
# in the specified range to values sampled from the M1 (neutral) model of
# variable omega ratios among sites.
#
# Distribution of omega values:
# \preformatted{
# CATEGORY PROBABILITY
#
# omega_0 = 0 p0
# omega_1 = 1 1-p0
# }
# }
#
# \references{
# Yang, Z., Nielsen, R., Goldman, N., Pedersen Krabbe, A-M. (2000) Codon-Substitution Models for Heterogeneous Selection Pressure at Amino Acid Sites - Genetics 155:431-449 \url{http://bit.ly/bvjucn}
#
# Goldman, N., Yang, Z. (1994) A codon-based model of nucleotide substitution for protein-coding DNA sequences - Mol Biol Evol 11(5):725-36 \url{http://bit.ly/aSVEoa}
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A CodonSequence object.}
# \item{process}{A process object inheriting from GY94.}
# \item{p0}{See above.}
# \item{index}{A vector of positions.}
# \item{...}{Not used.}
# }
#
# \value{
# Invisible TRUE.
# }
#
# \examples{
# # create a GY94 object
# p<-GY94(kappa=2)
# # create a CodonSequence object, attach process p
# s<-CodonSequence(length=25, processes=list(list(p)))
# # sample states
# sampleStates(s)
# # sample omegas in range 1:20 from model M1
# omegaVarM1(s,p,p0=0.5,1:20)
# # get omega values
# getOmegas(s,p)
# # get a histogram of omega values in range 1:20
# omegaHist(s,p,breaks=50,1:20)
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"omegaVarM1",
class="CodonSequence",
function(
this,
process,
p0,
index,
...
){
if(missing(process)){
throw("No process specified!\n");
}
if(!is.GY94(process)){
throw("The sepcified process is not a GY94 codon substitution process!\n");
}
else if(missing(p0)){
throw("No p0 value specified!\n");
}
else if((!is.numeric(p0))| (length(p0) != 1)){
throw("The new value must be a numeric vector of length 1!\n");
}
else if(p0 < 0 | p0 > 1){
throw("The p0 parameter must be in the [0,1] interval!\n");
}
else {
if(missing(index)){
index<-seq(along.with=this$.sites);
}
else {
index<-.checkIndexSanity(this, index);
}
for(site in this$.sites[index]){
setParameterAtSite(this=process,site=site, id="omega", value=sample(c(0,1), size=1, replace=FALSE, prob=c(p0,(1-p0))));
}
return(invisible(TRUE));
}
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: omegaVarM2 - selection
##
###########################################################################/**
#
# @RdocMethod omegaVarM2
#
# @title "The M2 (selection) model of variable omega ratios among sites"
#
# \description{
# @get "title".
#
# This method sets the \code{omega} site-process specific parameter
# in the specified range to values sampled from the M2 (selection) model of
# variable omega ratios among sites.
#
# Distribution of omega values:
# \preformatted{
# CATEGORY PROBABILITY
#
# omega_0 = 0 p0
# omega_1 = 1 p1
# omega_2 1-p0-p1
# }
# }
#
# \references{
# Yang, Z., Nielsen, R., Goldman, N., Pedersen Krabbe, A-M. (2000) Codon-Substitution Models for Heterogeneous Selection Pressure at Amino Acid Sites - Genetics 155:431-449 \url{http://bit.ly/bvjucn}
#
# Goldman, N., Yang, Z. (1994) A codon-based model of nucleotide substitution for protein-coding DNA sequences - Mol Biol Evol 11(5):725-36 \url{http://bit.ly/aSVEoa}
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A CodonSequence object.}
# \item{process}{A process object inheriting from GY94.}
# \item{p0}{See above.}
# \item{p1}{See above.}
# \item{omega_2}{See above.}
# \item{index}{A vector of positions.}
# \item{...}{Not used.}
# }
#
# \value{
# Invisible TRUE.
# }
#
# \examples{
# # create a GY94 object
# p<-GY94(kappa=2)
# # create a CodonSequence object, attach process p
# s<-CodonSequence(length=25, processes=list(list(p)))
# # sample states
# sampleStates(s)
# # sample omegas in range 1:20 from model M2
# omegaVarM2(s,p,p0=0.2,p1=0.3,omega_2=4,1:20)
# # get omega values
# getOmegas(s,p)
# # get a histogram of omega values in range 1:20
# omegaHist(s,p,breaks=50,1:20)
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"omegaVarM2",
class="CodonSequence",
function(
this,
process,
p0,
p1,
omega_2,
index,
...
){
if(missing(process)){
throw("No process specified!\n");
}
if(!is.GY94(process)){
throw("The sepcified process is not a GY94 codon substitution process!\n");
}
else if(missing(p0)){
throw("No p0 value specified!\n");
}
else if((!is.numeric(p0))| (length(p0) != 1)){
throw("The p0 value must be a numeric vector of length 1!\n");
}
else if(p0 < 0 | p0 > 1){
throw("The p0 parameter must be in the [0,1] interval!\n");
}
else if(missing(p1)){
throw("No p1 value specified!\n");
}
else if((!is.numeric(p1))| (length(p1) != 1)){
throw("The p1 value must be a numeric vector of length 1!\n");
}
else if(p1 < 0 | p1 > 1){
throw("The p1 parameter must be in the [0,1] interval!\n");
}
else if(missing(omega_2)){
throw("No omega value specified!\n");
}
else if((!is.numeric(omega_2))| (length(omega_2) != 1)){
throw("The omega value must be a numeric vector of length 1!\n");
}
else {
if(missing(index)){
index<-seq(along.with=this$.sites);
}
else {
index<-.checkIndexSanity(this, index);
}
for(site in this$.sites[index]){
setParameterAtSite(this=process,site=site, id="omega", value=sample(c(0,1,omega_2), size=1, replace=FALSE, prob=c(p0,p1,(1-p0-p1))));
}
return(invisible(TRUE));
}
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: omegaVarM3 - discrete
##
###########################################################################/**
#
# @RdocMethod omegaVarM3
#
# @title "The M3 (discrete) model of variable omega ratios among sites"
#
# \description{
# @get "title".
#
# This method sets the \code{omega} site-process specific parameter
# in the specified range to values sampled from the M3 (discrete) model of
# variable omega ratios among sites.
#
# Distribution of omega values:
# \preformatted{
# CATEGORY PROBABILITY
#
# omega_0 p0
# omega_1 p1
# omega_2 p2
# ... ...
# omega_k pk
# }
# }
#
# \references{
# Yang, Z., Nielsen, R., Goldman, N., Pedersen Krabbe, A-M. (2000) Codon-Substitution Models for Heterogeneous Selection Pressure at Amino Acid Sites - Genetics 155:431-449 \url{http://bit.ly/bvjucn}
#
# Goldman, N., Yang, Z. (1994) A codon-based model of nucleotide substitution for protein-coding DNA sequences - Mol Biol Evol 11(5):725-36 \url{http://bit.ly/aSVEoa}
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A CodonSequence object.}
# \item{process}{A process object inheriting from GY94.}
# \item{omegas}{A vector of omega values (omega_0 ... omega_k).}
# \item{probs}{A vector of probabilities (p0 ... pk).}
# \item{index}{A vector of positions.}
# \item{...}{Not used.}
# }
#
# \value{
# Invisible TRUE.
# }
#
# \examples{
# # create a GY94 object
# p<-GY94(kappa=2)
# # create a CodonSequence object, attach process p
# s<-CodonSequence(length=25, processes=list(list(p)))
# # sample states
# sampleStates(s)
# # sample omegas in range 1:20 from model M3
# omegaVarM3(s,p,omegas=c(0,2,4),probs=c(1/3,1/3,1/3),1:20)
# # get omega values
# getOmegas(s,p)
# # get a histogram of omega values in range 1:20
# omegaHist(s,p,breaks=50,1:20)
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"omegaVarM3",
class="CodonSequence",
function(
this,
process,
omegas,
probs,
index,
...
){
if(missing(process)){
throw("No process specified!\n");
}
if(!is.GY94(process)){
throw("The sepcified process is not a GY94 codon substitution process!\n");
}
else if(missing(omegas)){
throw("No omega values specified!\n");
}
else if((!is.numeric(omegas))){
throw("The omegas must be numeric!\n");
}
else if(missing(probs)){
throw("No probabilities specified!\n");
}
else if(!is.numeric(probs)){
throw("The probabilities must be numeric!\n");
}
else if(length(omegas) != length(probs)){
throw("The length of the \"omegas\" and \"probs\" vector must be the same!\n");
}
else if(!PSRoot$my.all.equal(sum(probs),1.0)){
probs<-(probs/sum(probs));
warning("The provided probabilities were scaked in order to sum to one!\n");
}
if(missing(index)){
index<-seq(along.with=this$.sites);
}
else {
index<-.checkIndexSanity(this, index);
}
for(site in this$.sites[index]){
setParameterAtSite(this=process,site=site, id="omega", value=sample(omegas, size=1, replace=FALSE, prob=probs));
}
return(invisible(TRUE));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: omegaVarM4 - freqs
##
###########################################################################/**
#
# @RdocMethod omegaVarM4
#
# @title "The M4 (freqs) model of variable omega ratios among sites"
#
# \description{
# @get "title".
#
# This method sets the \code{omega} site-process specific parameter
# in the specified range to values sampled from the M4 (freqs) model of
# variable omega ratios among sites.
#
# Distribution of omega values:
# \preformatted{
# CATEGORY PROBABILITY
#
# omega_0 = 0 p0
# omega_1 = 1/3 p1
# omega_2 = 2/3 p2
# omega_3 = 1 p3
# omega_4 = 3 p4
# }
# }
#
# \references{
# Yang, Z., Nielsen, R., Goldman, N., Pedersen Krabbe, A-M. (2000) Codon-Substitution Models for Heterogeneous Selection Pressure at Amino Acid Sites - Genetics 155:431-449 \url{http://bit.ly/bvjucn}
#
# Goldman, N., Yang, Z. (1994) A codon-based model of nucleotide substitution for protein-coding DNA sequences - Mol Biol Evol 11(5):725-36 \url{http://bit.ly/aSVEoa}
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A CodonSequence object.}
# \item{process}{A process object inheriting from GY94.}
# \item{probs}{A vector of probabilities (p0 ... p4).}
# \item{index}{A vector of positions.}
# \item{...}{Not used.}
# }
#
# \value{
# Invisible TRUE.
# }
#
# \examples{
# # create a GY94 object
# p<-GY94(kappa=2)
# # create a CodonSequence object, attach process p
# s<-CodonSequence(length=25, processes=list(list(p)))
# # sample states
# sampleStates(s)
# # sample omegas in range 1:20 from model M4
# omegaVarM4(s,p,probs=c(2/5,1/5,1/5,1/10,1/10),1:20)
# # get omega values
# getOmegas(s,p)
# # get a histogram of omega values in range 1:20
# omegaHist(s,p,breaks=50,1:20)
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"omegaVarM4",
class="CodonSequence",
function(
this,
process,
probs,
index,
...
){
if(missing(process)){
throw("No process specified!\n");
}
if(!is.GY94(process)){
throw("The sepcified process is not a GY94 codon substitution process!\n");
}
else if(missing(probs)){
throw("No probabilities specified!\n");
}
else if(!is.numeric(probs)){
throw("The probs must be greater than zero!\n");
}
else if( length(probs) != 5){
throw("The length of the \"probs\" vector must be 5!\n");
}
else if(!PSRoot$my.all.equal(sum(probs),1.0)){
probs<-(probs/sum(probs));
warning("The provided probabilities were scaked in order to sum to one!\n");
}
if(missing(index)){
index<-seq(along.with=this$.sites);
}
else {
index<-.checkIndexSanity(this, index);
}
omegas<-c(0,(1/3),(2/3),1,3);
for(site in this$.sites[index]){
setParameterAtSite(this=process,site=site, id="omega", value=sample(omegas, size=1, replace=FALSE, prob=probs));
}
return(invisible(TRUE));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: .omegaVarM5Cont - gamma
##
setMethodS3(
".omegaVarM5Cont",
class="CodonSequence",
function(
this,
process,
alpha,
beta,
index,
...
){
if(missing(process)){
throw("No process specified!\n");
}
if(!is.GY94(process)){
throw("The sepcified process is not a GY94 codon substitution process!\n");
}
else if(missing(alpha)){
throw("No alpha (shape) value specified!\n");
}
else if((!is.numeric(alpha)) | (length(alpha) != 1)){
throw("The alpha (shape) parameter must be a numeric vector of length 1!\n");
}
else if(alpha < 0){
throw("The alpha (shape) must be greater than zero!\n");
}
else if(missing(beta)){
throw("No beta (scale) value specified!\n");
}
else if((!is.numeric(beta)) | (length(beta) != 1)){
throw("The beta (scale) parameter must be a numeric vector of length 1!\n");
}
else if(beta <= 0){
throw("The beta (scale) must be strictly positive!\n");
}
if(missing(index)){
index<-seq(along.with=this$.sites);
}
else {
index<-.checkIndexSanity(this, index);
}
for(site in this$.sites[index]){
setParameterAtSite(this=process,site=site, id="omega", value=rgamma(1,shape=alpha,scale=beta));
}
return(invisible(TRUE));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: .omegaVarM6Cont - 2gamma
##
setMethodS3(
".omegaVarM6Cont",
class="CodonSequence",
function(
this,
process,
p0,
alpha0,
beta0,
alpha1,
beta1,
index,
...
){
if(missing(process)){
throw("No process specified!\n");
}
if(!is.GY94(process)){
throw("The sepcified process is not a GY94 codon substitution process!\n");
}
else if(missing(p0)){
throw("No p0 value specified!\n");
}
else if((!is.numeric(p0)) | (length(p0) != 1)){
throw("The p0 parameter must be a numeric vector of length 1!\n");
}
else if( (p0 < 0) | (p0 > 1)){
throw("p0 must be in the [0,1] interval!\n");
}
else if(missing(alpha0)){
throw("No alpha0 (shape0) value specified!\n");
}
else if((!is.numeric(alpha0)) | (length(alpha0) != 1)){
throw("The alpha0 (shape0) parameter must be a numeric vector of length 1!\n");
}
else if(alpha0 < 0){
throw("The alpha0 (shape0) must be greater than zero!\n");
}
else if(missing(beta0)){
throw("No beta0 (scale0) value specified!\n");
}
else if((!is.numeric(beta0)) | (length(beta0) != 1)){
throw("The beta0 (scale0) parameter must be a numeric vector of length 1!\n");
}
else if(beta0 <= 0){
throw("The beta0 (scale0) must be strictly positive!\n");
}
else if(missing(alpha1)){
throw("No alpha1 (shape1) value specified!\n");
}
else if((!is.numeric(alpha1)) | (length(alpha1) != 1)){
throw("The alpha1 (shape1) parameter must be a numeric vector of length 1!\n");
}
else if(alpha1 < 0){
throw("The alpha1 (shape1) must be greater than zero!\n");
}
else if(missing(beta1)){
throw("No beta1 (scale1) value specified!\n");
}
else if((!is.numeric(beta1)) | (length(beta1) != 1)){
throw("The beta1 (scale1) parameter must be a numeric vector of length 1!\n");
}
else if(beta1 <= 0){
throw("The beta1 (scale1) must be strictly positive!\n");
}
if(missing(index)){
index<-seq(along.with=this$.sites);
}
else {
index<-.checkIndexSanity(this, index);
}
for(site in this$.sites[index]){
setParameterAtSite(this=process,site=site, id="omega", value=sample(c(rgamma(1,shape=alpha0,scale=beta0),rgamma(1,shape=alpha1,scale=beta1)),size=1,replace=FALSE,prob=c(p0,(1-p0)) ));
}
return(invisible(TRUE));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: .omegaVarM7Cont - beta
##
setMethodS3(
".omegaVarM7",
class="CodonSequence",
function(
this,
process,
p,
q,
index,
...
){
if(missing(process)){
throw("No process specified!\n");
}
if(!is.GY94(process)){
throw("The sepcified process is not a GY94 codon substitution process!\n");
}
else if(missing(p)){
throw("No p value specified!\n");
}
else if((!is.numeric(p)) | (length(p) != 1)){
throw("The p parameter must be a numeric vector of length 1!\n");
}
else if(p < 0){
throw("The p parameter must be greater than zero!\n");
}
else if(missing(q)){
throw("No q value specified!\n");
}
else if((!is.numeric(q)) | (length(q) != 1)){
throw("The q parameter must be a numeric vector of length 1!\n");
}
else if(q < 0){
throw("The q parameter must be positive!\n");
}
if(missing(index)){
index<-seq(along.with=this$.sites);
}
else {
index<-.checkIndexSanity(this, index);
}
for(site in this$.sites[index]){
setParameterAtSite(this=process,site=site, id="omega", value=rbeta(1,shape1=p,shape2=q));
}
return(invisible(TRUE));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: .omegaVarM8Cont - beta&omega
##
setMethodS3(
".omegaVarM8Cont",
class="CodonSequence",
function(
this,
process,
p0,
p,
q,
omega,
index,
...
){
if(missing(process)){
throw("No process specified!\n");
}
if(!is.GY94(process)){
throw("The sepcified process is not a GY94 codon substitution process!\n");
}
else if(missing(p0)){
throw("No p0 value specified!\n");
}
else if((!is.numeric(p0)) | (length(p0) != 1)){
throw("The p0 parameter must be a numeric vector of length 1!\n");
}
else if( (p0 < 0) | (p0 > 1)){
throw("The p0 parameter must be from the [0,1] interval!\n");
}
else if(missing(p)){
throw("No p value specified!\n");
}
else if((!is.numeric(p)) | (length(p) != 1)){
throw("The p parameter must be a numeric vector of length 1!\n");
}
else if(p < 0){
throw("The p parameter must be greater than zero!\n");
}
else if(missing(q)){
throw("No q value specified!\n");
}
else if((!is.numeric(q)) | (length(q) != 1)){
throw("The q parameter must be a numeric vector of length 1!\n");
}
else if(q < 0){
throw("The q parameter must be positive!\n");
}
else if(missing(omega)){
throw("No omega value specified!\n");
}
else if((!is.numeric(omega)) | (length(omega) != 1)){
throw("The omega parameter must be a numeric vector of length 1!\n");
}
else if(omega < 0){
throw("The omega parameter must be positive!\n");
}
if(missing(index)){
index<-seq(along.with=this$.sites);
}
else {
index<-.checkIndexSanity(this, index);
}
for(site in this$.sites[index]){
setParameterAtSite(this=process,site=site, id="omega", value=sample(c(rbeta(1,shape1=p,shape2=q),omega),size=1,replace=FALSE,prob=(c(p0,(1-p0)))));
}
return(invisible(TRUE));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: .omegaVarM9 - beta&gamma
##
setMethodS3(
".omegaVarM9Cont",
class="CodonSequence",
function(
this,
process,
p0,
p,
q,
alpha,
beta,
index,
...
){
if(missing(process)){
throw("No process specified!\n");
}
if(!is.GY94(process)){
throw("The sepcified process is not a GY94 codon substitution process!\n");
}
else if(missing(p0)){
throw("No p0 value specified!\n");
}
else if((!is.numeric(p0)) | (length(p0) != 1)){
throw("The p0 parameter must be a numeric vector of length 1!\n");
}
else if( (p0 < 0) | (p0 > 1)){
throw("p0 must be in the [0,1] interval!\n");
}
else if(missing(p)){
throw("No p value specified!\n");
}
else if((!is.numeric(p)) | (length(p) != 1)){
throw("The p parameter must be a numeric vector of length 1!\n");
}
else if(p < 0){
throw("The p must be greater than zero!\n");
}
else if(missing(q)){
throw("No q value specified!\n");
}
else if((!is.numeric(q)) | (length(q) != 1)){
throw("The q parameter must be a numeric vector of length 1!\n");
}
else if(q < 0){
throw("The q must be greater than zero!\n");
}
else if(missing(alpha)){
throw("No alpha (shape) value specified!\n");
}
else if((!is.numeric(alpha)) | (length(alpha) != 1)){
throw("The alpha (shape) parameter must be a numeric vector of length 1!\n");
}
else if(alpha < 0){
throw("The alpha (shape) must be greater than zero!\n");
}
else if(missing(beta)){
throw("No beta (scale) value specified!\n");
}
else if((!is.numeric(beta)) | (length(beta) != 1)){
throw("The beta (scale) parameter must be a numeric vector of length 1!\n");
}
else if(beta <= 0){
throw("The beta (scale) must be strictly positive!\n");
}
if(missing(index)){
index<-seq(along.with=this$.sites);
}
else {
index<-.checkIndexSanity(this, index);
}
for(site in this$.sites[index]){
setParameterAtSite(this=process,site=site, id="omega", value=sample(c(rbeta(1,shape1=p,shape2=q),rgamma(1,shape=alpha,scale=beta)),size=1,replace=FALSE,prob=c(p0,(1-p0)) ));
}
return(invisible(TRUE));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: .omegaVarM10Cont - beta&gamma+1
##
setMethodS3(
".omegaVarM10Cont",
class="CodonSequence",
function(
this,
process,
p0,
p,
q,
alpha,
beta,
index,
...
){
if(missing(process)){
throw("No process specified!\n");
}
if(!is.GY94(process)){
throw("The sepcified process is not a GY94 codon substitution process!\n");
}
else if(missing(p0)){
throw("No p0 value specified!\n");
}
else if((!is.numeric(p0)) | (length(p0) != 1)){
throw("The p0 parameter must be a numeric vector of length 1!\n");
}
else if( (p0 < 0) | (p0 > 1)){
throw("p0 must be in the [0,1] interval!\n");
}
else if(missing(p)){
throw("No p value specified!\n");
}
else if((!is.numeric(p)) | (length(p) != 1)){
throw("The p parameter must be a numeric vector of length 1!\n");
}
else if(p < 0){
throw("The p must be greater than zero!\n");
}
else if(missing(q)){
throw("No q value specified!\n");
}
else if((!is.numeric(q)) | (length(q) != 1)){
throw("The q parameter must be a numeric vector of length 1!\n");
}
else if(q < 0){
throw("The q must be greater than zero!\n");
}
else if(missing(alpha)){
throw("No alpha (shape) value specified!\n");
}
else if((!is.numeric(alpha)) | (length(alpha) != 1)){
throw("The alpha (shape) parameter must be a numeric vector of length 1!\n");
}
else if(alpha < 0){
throw("The alpha (shape) must be greater than zero!\n");
}
else if(missing(beta)){
throw("No beta (scale) value specified!\n");
}
else if((!is.numeric(beta)) | (length(beta) != 1)){
throw("The beta (scale) parameter must be a numeric vector of length 1!\n");
}
else if(beta <= 0){
throw("The beta (scale) must be strictly positive!\n");
}
if(missing(index)){
index<-seq(along.with=this$.sites);
}
else {
index<-.checkIndexSanity(this, index);
}
# It's not too elegant to fork the whole method for just shifting the gamma with 1...
for(site in this$.sites[index]){
setParameterAtSite(this=process,site=site, id="omega", value=sample(c(rbeta(1,shape1=p,shape2=q),(1 + rgamma(1,shape=alpha,scale=beta)) ),size=1,replace=FALSE,prob=c(p0,(1-p0)) ));
}
return(invisible(TRUE));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: .omegaVarM11Cont - beta&normal>1
##
setMethodS3(
".omegaVarM11Cont",
class="CodonSequence",
function(
this,
process,
p0,
p,
q,
mean,
sd,
index,
...
){
if(missing(process)){
throw("No process specified!\n");
}
if(!is.GY94(process)){
throw("The sepcified process is not a GY94 codon substitution process!\n");
}
else if(missing(p0)){
throw("No p0 value specified!\n");
}
else if((!is.numeric(p0)) | (length(p0) != 1)){
throw("The p0 parameter must be a numeric vector of length 1!\n");
}
else if( (p0 < 0) | (p0 > 1)){
throw("p0 must be in the [0,1] interval!\n");
}
else if(missing(p)){
throw("No p value specified!\n");
}
else if((!is.numeric(p)) | (length(p) != 1)){
throw("The p parameter must be a numeric vector of length 1!\n");
}
else if(p < 0){
throw("The p must be greater than zero!\n");
}
else if(missing(q)){
throw("No q value specified!\n");
}
else if((!is.numeric(q)) | (length(q) != 1)){
throw("The q parameter must be a numeric vector of length 1!\n");
}
else if(q < 0){
throw("The q must be greater than zero!\n");
}
else if(missing(mean)){
throw("No mean specified!\n");
}
else if((!is.numeric(mean)) | (length(mean) != 1)){
throw("The mean parameter must be a numeric vector of length 1!\n");
}
else if(missing(sd)){
throw("No sd value specified!\n");
}
else if((!is.numeric(sd)) | (length(sd) != 1)){
throw("The sd parameter must be a numeric vector of length 1!\n");
}
else if( sd < 0){
throw("The sd parameter must be positive!\n");
}
if(missing(index)){
index<-seq(along.with=this$.sites);
}
else {
index<-.checkIndexSanity(this, index);
}
rnorm.gt.1<-function(mean=NA,sd=NA){
# This is probably the most primitive way to truncate the distribution!
tmp<-rnorm(1,mean=mean,sd=sd);
while( tmp <= 1){
tmp<-rnorm(1,mean=mean,sd=sd);
};
return(tmp);
}
for(site in this$.sites[index]){
setParameterAtSite(this=process,site=site, id="omega", value=sample(c(rbeta(1,shape1=p,shape2=q), (rnorm.gt.1(mean=mean,sd=sd))),size=1,replace=FALSE,prob=c(p0,(1-p0)) ));
}
return(invisible(TRUE));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: .omegaVarM12Cont - 0&2normal>0
##
setMethodS3(
".omegaVarM12Cont",
class="CodonSequence",
function(
this,
process,
p0,
p1,
sd1,
mean2,
sd2,
index,
...
){
if(missing(process)){
throw("No process specified!\n");
}
if(!is.GY94(process)){
throw("The sepcified process is not a GY94 codon substitution process!\n");
}
else if(missing(p0)){
throw("No p0 value specified!\n");
}
else if((!is.numeric(p0)) | (length(p0) != 1)){
throw("The p0 parameter must be a numeric vector of length 1!\n");
}
else if( (p0 < 0) | (p0 > 1)){
throw("p0 must be in the [0,1] interval!\n");
}
else if(missing(p1)){
throw("No p1 value specified!\n");
}
else if((!is.numeric(p1)) | (length(p1) != 1)){
throw("The p1 parameter must be a numeric vector of length 1!\n");
}
else if( (p1 < 0) | (p1 > 1)){
throw("p1 must be in the [0,1] interval!\n");
}
else if(missing(sd1)){
throw("No sd1 value specified!\n");
}
else if((!is.numeric(sd1)) | (length(sd1) != 1)){
throw("The sd1 parameter must be a numeric vector of length 1!\n");
}
else if( sd1 < 0){
throw("The sd1 parameter must be positive!\n");
}
else if(missing(mean2)){
throw("No mean2 specified!\n");
}
else if((!is.numeric(mean2)) | (length(mean2) != 1)){
throw("The mean2 parameter must be a numeric vector of length 1!\n");
}
else if(missing(sd2)){
throw("No sd2 value specified!\n");
}
else if((!is.numeric(sd2)) | (length(sd2) != 1)){
throw("The sd2 parameter must be a numeric vector of length 1!\n");
}
else if( sd2 < 0){
throw("The sd2 parameter must be positive!\n");
}
if(missing(index)){
index<-seq(along.with=this$.sites);
}
else {
index<-.checkIndexSanity(this, index);
}
rnorm.gt.0<-function(mean=NA,sd=NA){
# This is probably the most primitive way to truncate the distribution!
tmp<-rnorm(1,mean=mean,sd=sd);
while( tmp <= 0){
tmp<-rnorm(1,mean=mean,sd=sd);
};
return(tmp);
}
for(site in this$.sites[index]){
setParameterAtSite(this=process,site=site, id="omega", value=sample(
c(0,rnorm.gt.0(mean=1,sd=sd1),rnorm.gt.0(mean=mean2,sd=sd2)),
size=1,
replace=FALSE,
prob=c(p0,((1-p0)*p1),((1-p0)*(1-p1)))
));
}
return(invisible(TRUE));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: .omegaVarM13Cont - 3normal>0
##
setMethodS3(
".omegaVarM13Cont",
class="CodonSequence",
function(
this,
process,
p0,
p1,
sd0,
sd1,
mean2,
sd2,
index,
...
){
if(missing(process)){
throw("No process specified!\n");
}
if(!is.GY94(process)){
throw("The sepcified process is not a GY94 codon substitution process!\n");
}
else if(missing(p0)){
throw("No p0 value specified!\n");
}
else if((!is.numeric(p0)) | (length(p0) != 1)){
throw("The p0 parameter must be a numeric vector of length 1!\n");
}
else if( (p0 < 0) | (p0 > 1)){
throw("p0 must be in the [0,1] interval!\n");
}
else if(missing(p1)){
throw("No p1 value specified!\n");
}
else if((!is.numeric(p1)) | (length(p1) != 1)){
throw("The p1 parameter must be a numeric vector of length 1!\n");
}
else if( (p1 < 0) | (p1 > 1)){
throw("p1 must be in the [0,1] interval!\n");
}
else if(missing(sd0)){
throw("No sd0 value specified!\n");
}
else if((!is.numeric(sd0)) | (length(sd0) != 1)){
throw("The sd0 parameter must be a numeric vector of length 1!\n");
}
else if( sd0 < 0){
throw("The sd0 parameter must be positive!\n");
}
else if(missing(sd1)){
throw("No sd1 value specified!\n");
}
else if((!is.numeric(sd1)) | (length(sd1) != 1)){
throw("The sd1 parameter must be a numeric vector of length 1!\n");
}
else if( sd1 < 0){
throw("The sd1 parameter must be positive!\n");
}
else if(missing(mean2)){
throw("No mean2 specified!\n");
}
else if((!is.numeric(mean2)) | (length(mean2) != 1)){
throw("The mean2 parameter must be a numeric vector of length 1!\n");
}
else if(missing(sd2)){
throw("No sd2 value specified!\n");
}
else if((!is.numeric(sd2)) | (length(sd2) != 1)){
throw("The sd2 parameter must be a numeric vector of length 1!\n");
}
else if( sd2 < 0){
throw("The sd2 parameter must be positive!\n");
}
if(missing(index)){
index<-seq(along.with=this$.sites);
}
else {
index<-.checkIndexSanity(this, index);
}
rnorm.gt.0<-function(mean=NA,sd=NA){
# This is probably the most primitive way to truncate the distribution!
tmp<-rnorm(1,mean=mean,sd=sd);
while( tmp <= 0){
tmp<-rnorm(1,mean=mean,sd=sd);
};
return(tmp);
}
for(site in this$.sites[index]){
setParameterAtSite(this=process,site=site, id="omega", value=sample(
c(rnorm.gt.0(mean=0,sd=sd0),rnorm.gt.0(mean=1,sd=sd1),rnorm.gt.0(mean=mean2,sd=sd2)),
size=1,
replace=FALSE,
prob=c(p0,((1-p0)*p1),((1-p0)*(1-p1)))
));
}
return(invisible(TRUE));
},
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{...}{Not used.}
# }
#
#
# \value{
# Returns an invisible TRUE if no inconsistencies found in the object, throws
# an error otherwise.
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"checkConsistency",
class="CodonSequence",
function(
this,
...
){
wp<-this$writeProtected;
if (wp) {
this$writeProtected<-FALSE;
}
may.fail<-function(this) {
for(pos in seq(along.with=this$.sites)){
if(!is.CodonAlphabet(this$.sites[[pos]]$.alphabet)){
throw("The alphabet attached to site ",pos," is not a codon alphabet!\n");
}
}
}
tryCatch(may.fail(this),finally=this$writeProtected<-wp);
NextMethod();
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: Translate
##
###########################################################################/**
#
# @RdocMethod Translate
#
# @title "Translate a CodonSequence object"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A CodonSequence object}
# \item{...}{Not used.}
# }
#
# \value{
# The translation as an AminoAcidSequence object.
# }
#
# \examples{
# # create a CodonSequence object
# s<-CodonSequence(string="ATCTTTCGAATGGGGCCCTCCCGA")
# # get the translation as an AminoAcidSequence object
# as<-Translate(s)
# as
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"Translate",
class="CodonSequence",
function(
this,
...
){
# Create an amino acid sequence of proper length:
that<-AminoAcidSequence(
name=paste("Translation of",this$name),
length=this$length
);
# Setting the states, using the site alphabet to translate codons for greater flexibility.
for (pos in seq(along.with=this$.sites)){
setState(that$.sites[[pos]], translateCodon(this$.sites[[pos]]$alphabet, this$.sites[[pos]]$state) );
}
return(that);
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
Try the phylosim package in your browser
Any scripts or data that you put into this service are public.
phylosim documentation built on Nov. 22, 2019, 1:07 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
##
## Copyright 2009 Botond Sipos
## See the package description for licensing information.
##
##
## BinarySequence
##
##########################################################################/**
#
# @RdocClass BinarySequence
#
# @title "The BinarySequence class"
#
# \description{
# Sequence objects aggregating Site objects having a BinaryAlphabet attached by default.
#
# @classhierarchy
# }
#
# @synopsis
#
# \arguments{
# \item{name}{Name of the Sequence object.}
# \item{string}{A string specifying the length and the states of the Sequence object.}
# \item{length}{The length of the sequence. Mutually exclusive with "string".}
# \item{processes}{A list of lists of Process objects, to be attached to the aggregated Site objects. Recycled if shorter than the length of the sequence.}
# \item{ancestral.obj}{The ancestral object of the Sequence object (a valid Sequence or Process object).}
# \item{...}{Not used.}
# }
#
# \section{Fields and Methods}{
# @allmethods
# }
#
# \examples{
# # create an empty BinarySequence object
# s<-BinarySequence(length=50)
# s
# # set states
# s$states<-c(0,0,1,0,1,1)
# s
# # create a sequence object by specifying a string
# s<-BinarySequence(string="00000110010001111")
# s
# }
#
# @author
#
# \seealso{
# Sequence BinaryAlphabet
# }
#
#*/###########################################################################
setConstructorS3(
"BinarySequence",
function(
name=NA,
string=NA,
length=NA,
processes=NA,
ancestral.obj=NA,
...
){
if(!missing(string)){
this<-Sequence(string=string,alphabets=list(BinaryAlphabet()));
}
else if(!missing(length)) {
this<-Sequence(length=length,alphabets=list(BinaryAlphabet()));
}
else {
this<-Sequence(alphabets=list(BinaryAlphabet()));
}
this<-extend(this, "BinarySequence");
if(!missing(name)){
this$name<-name;
}
if(!missing(processes)){
setProcesses(this,processes);
}
if(!missing(ancestral.obj)){
this$.ancestral<-ancestral.obj;
}
return(this);
},
enforceRCC=TRUE
);
##
## NucleotideSequence
##
##########################################################################/**
#
# @RdocClass NucleotideSequence
#
# @title "The NucleotideSequence class"
#
# \description{
# Sequence objects aggregating Site objects having a NucleotideAlphabet attached by default.
#
# @classhierarchy
# }
#
# @synopsis
#
# \arguments{
# \item{name}{Name of the Sequence object.}
# \item{string}{A string specifying the length and the states of the Sequence object.}
# \item{length}{The length of the sequence. Mutually exclusive with "string".}
# \item{processes}{A list of lists of Process objects, to be attached to the aggregated Site objects. Recycled if shorter than the length of the sequence.}
# \item{ancestral.obj}{The ancestral object of the Sequence object (a valid Sequence or Process object).}
# \item{...}{Not used.}
# }
#
# \section{Fields and Methods}{
# @allmethods
# }
#
# \examples{
# # create an empty NucleotideSequence object
# s<-NucleotideSequence(length=50)
# s
# # set states
# s$states<-c("A","A","G","T")
# s
# # create a sequence object by specifying a string
# s<-NucleotideSequence(string="ATGCCGATTAGCAAA")
# s
# }
#
# @author
#
# \seealso{
# Sequence NucleotideAlphabet
# }
#
#*/###########################################################################
setConstructorS3(
"NucleotideSequence",
function(
name=NA,
string=NA,
length=NA,
processes=NA,
ancestral.obj=NA,
...
){
if(!missing(string)){
this<-Sequence(string=string,alphabets=list(NucleotideAlphabet()));
}
else if(!missing(length)) {
this<-Sequence(length=length,alphabets=list(NucleotideAlphabet()));
}
else {
this<-Sequence(alphabets=list(NucleotideAlphabet()));
}
this<-extend(this, "NucleotideSequence");
if(!missing(name)){
this$name<-name;
}
if(!missing(processes)){
setProcesses(this,processes);
}
if(!missing(ancestral.obj)){
this$.ancestral<-ancestral.obj;
}
return(this);
},
enforceRCC=TRUE
);
##
## revComp
##
###########################################################################/**
#
# @RdocMethod revComp
#
# @title "Reverse complmenet a NucleotideSequence object"
#
# \description{
# @get "title".
#
# The method reverse complements the sequence "in place", no object cloning is performed.
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A NucleotideSequence object}
# \item{...}{Not used.}
# }
#
# \value{
# Invisible TRUE.
# }
#
# \examples{
# s<-NucleotideSequence(string="ATGC")
# s
# revComp(s)
# s
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"revComp",
class="NucleotideSequence",
function(
this,
...
){
states<-rev(as.character(getStates(this)));
len<-length(states);
if(len == 0){return(invisible(TRUE))}
for(i in 1:len){
if(states[[i]] == "NA"){
state<-NA;
}
else if(states[i] == "A"){
state<-"T";
}
else if(states[i] == "T"){
state<-"A";
}
else if(states[i] == "G"){
state<-"C";
}
else if(states[i] == "C"){
state<-"G";
} else {
throw("Symbol not in NucleotideAlphabet!");
}
this$.sites[[i]]$.state<-state;
}
return(invisible(TRUE));
},
private=TRUE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## AminoAcidSequence
##
##########################################################################/**
#
# @RdocClass AminoAcidSequence
#
# @title "The AminoAcidSequence class"
#
# \description{
# Sequence objects aggregating Site objects having an AminoAcidAlphabet attached by default.
#
# @classhierarchy
# }
#
# @synopsis
#
# \arguments{
# \item{name}{Name of the Sequence object.}
# \item{string}{A string specifying the length and the states of the Sequence object.}
# \item{length}{The length of the sequence. Mutually exclusive with "string".}
# \item{processes}{A list of lists of Process objects, to be attached to the aggregated Site objects. Recycled if shorter than the length of the sequence.}
# \item{ancestral.obj}{The ancestral object of the Sequence object (a valid Sequence or Process object).}
# \item{...}{Not used.}
# }
#
# \section{Fields and Methods}{
# @allmethods
# }
#
# \examples{
# # create an empty AminoAcidSequence object
# s<-AminoAcidSequence(length=50)
# s
# # set states
# s$states<-c("C","C","G","Y")
# s
# # create a sequence object by specifying a string
# s<-AminoAcidSequence(string="CNGGYCCNGYYYY")
# s
# }
#
# @author
#
# \seealso{
# Sequence AminoAcidAlphabet
# }
#
#*/###########################################################################
setConstructorS3(
"AminoAcidSequence",
function(
name=NA,
string=NA,
length=NA,
processes=NA,
ancestral.obj=NA,
...
){
if(!missing(string)){
this<-Sequence(string=string,alphabets=list(AminoAcidAlphabet()));
}
else if(!missing(length)) {
this<-Sequence(length=length,alphabets=list(AminoAcidAlphabet()));
}
else {
this<-Sequence(alphabets=list(AminoAcidAlphabet()));
}
this<-extend(this, "AminoAcidSequence");
if(!missing(name)){
this$name<-name;
}
if(!missing(processes)){
setProcesses(this,processes);
}
if(!missing(ancestral.obj)){
this$.ancestral<-ancestral.obj;
}
return(this);
},
enforceRCC=TRUE
);
##
## CodonSequence
##
##########################################################################/**
#
# @RdocClass CodonSequence
#
# @title "The CodonSequence class"
#
# \description{
# Sequence objects aggregating Site objects having a CodonAlphabet attached by default.
#
# @classhierarchy
# }
#
# @synopsis
#
# \arguments{
# \item{name}{Name of the Sequence object.}
# \item{string}{A string specifying the length and the states of the Sequence object.}
# \item{length}{The length of the sequence. Mutually exclusive with "string".}
# \item{table.id}{The genetic code table to use in the attached CodonAlphabet object ("Standard" by default).}
# \item{processes}{A list of lists of Process objects, to be attached to the aggregated Site objects. Recycled if shorter than the length of the sequence.}
# \item{ancestral.obj}{The ancestral object of the Sequence object (a valid Sequence or Process object).}
# \item{...}{Not used.}
# }
#
# \section{Fields and Methods}{
# @allmethods
# }
#
# \examples{
# # create an empty CodonSequence object
# s<-CodonSequence(length=50)
# s
# # set states
# s$states<-c("ATG","CGA","TTT","CTA")
# s
# # create a sequence object by specifying a string
# s<-CodonSequence(string="ATCTTTCGAATG")
# s
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setConstructorS3(
"CodonSequence",
function(
name=NA,
string=NA,
length=NA,
table.id=1,
processes=NA,
ancestral.obj=NA,
...
){
if (!missing(string)){
len<-((stringLength(string)/3));
flen<-floor(len);
this<-Sequence(length=flen,alphabets=list(CodonAlphabet(table.id=table.id)));
if(len > flen){
warning("The length of the provided string was not multiple of 3. The incomplete codon was discarded!\n");
}
nuc<-strsplit(string,"",fixed=TRUE)[[1]];
j<-1; # counter for the codon position
for(i in seq(from=1,to=(flen * 3),by=3)){
# get the codon:
state<-paste(nuc[c(i,i+1,i+2)],collapse="");
# Check for stop codons:
if( isStopCodon(this$.sites[[j]]$.alphabet, state) ){
throw("The CodonSequence objects does not accept stop codons as valid states!\n");
}
# Set the state:
setStates(this,state,j);
j<-j+1;
}
}
else if(!missing(length)){
this<-Sequence(length=length,alphabets=list(CodonAlphabet(table.id=table.id)));
}
else {
this<-Sequence(alphabets=list(CodonAlphabet(table.id=table.id)));
}
this<-extend(this, "CodonSequence");
if(!missing(name)){
this$name<-name;
}
if(!missing(processes)){
setProcesses(this,processes);
}
if(!missing(ancestral.obj)){
this$.ancestral<-ancestral.obj;
}
return(this);
},
enforceRCC=TRUE
);
##
## Method: getOmegas
##
###########################################################################/**
#
# @RdocMethod getOmegas
# \alias{getOmegas.Sequence}
#
# @title "Get the omegas from a collection of sites"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A CodonSequence object.}
# \item{process}{A process object inheriting from GY94.}
# \item{index}{A vector of positions.}
# \item{...}{Not used.}
# }
#
# \value{
# A numeric vector.
# }
#
# \examples{
# # create a GY94 process
# p<-GY94()
# # create a CodonSequence object,
# # attach a process p
# s<-CodonSequence(length=20,processes=list(list(p)))
# # set omega values in range 1:5
# setOmegas(s,p,c(0.5,1,1.5),1:5)
# # get omega values from siutes 1,2,3,10, and 20
# getOmegas(s,p,c(1:3,10,20))
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"getOmegas",
class="CodonSequence",
function(
this,
process,
index,
...
){
if(missing(process)){
throw("No process given!\n");
}
else if(!is.GY94(process)){
throw("The specified process is not a GY94 codon substitution process!\n");
}
rm<-getParameterAtSites(this=this,process=process,id="omega",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: omegaHist
##
###########################################################################/**
#
# @RdocMethod omegaHist
#
# @title "Plot a histogram of omega values from a range"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A CodonSequence object.}
# \item{process}{A process object inheriting from GY94.}
# \item{breaks}{\code{breaks} parameter for \code{hist()}.}
# \item{index}{A vector of positions.}
# \item{...}{Not used.}
# }
#
# \value{
# The CodonSequence object (invisible).
# }
#
# \examples{
# # create a GY94 process
# p<-GY94()
# # create a CodonSequence object,
# # attach a process p
# s<-CodonSequence(length=20,processes=list(list(p)))
# # set omega values through omegaVarM2.CodonSequence
# omegaVarM2(s,p,p0=0.5,p1=0.2,omega=1.5)
# # get a histogram of omega values from the range 1:15
# omegaHist(s,p,breaks=10,1:15)
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"omegaHist",
class="CodonSequence",
function(
this,
process,
breaks,
index,
...
){
if(missing(process)){
throw("No process given!\n");
}
else if(!is.GY94(process)){
throw("The specified process is not a GY94 codon substitution process!\n");
}
if(missing(index)){
index<-seq(along.with=this$.sites);
}
if(missing(breaks)){
hist(getOmegas(this,process,index));
}
else {
omegas<-getOmegas(this,process,index);
hist(omegas,breaks=breaks,main="Histogram of omega values",xlab="Omega",freq=FALSE);
}
return(invisible(this));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: setOmegas
###########################################################################/**
#
# @RdocMethod setOmegas
# \alias{setOmegas.Sequence}
#
# @title "Set the omegas for a collection of sites"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A CodonSequence object.}
# \item{process}{A process object inheriting from GY94.}
# \item{value}{A vector containing the new values of the site-process specific parameter, recycled if shorter than the index vector.}
# \item{index}{A vector of positions. It is set to 1:seq$length if omitted.}
# \item{...}{Not used.}
# }
#
# \value{
# A numeric vector.
# }
#
# \examples{
# # create a GY94 process
# p<-GY94()
# # create a CodonSequence object,
# # attach a process p
# s<-CodonSequence(length=20,processes=list(list(p)))
# # set omega values in range 1:5
# setOmegas(s,p,c(0.5,1,1.5),1:5)
# # get omega values from siutes 1,2,3,10, and 20
# getOmegas(s,p,c(1:3,10,20))
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"setOmegas",
class="CodonSequence",
function(
this,
process,
value,
index,
...
){
if(missing(process)){
throw("No process specified!\n");
}
if(!is.GY94(process)){
throw("The sepcified process is not a GY94 codon substitution process!\n");
}
else if(missing(value)){
throw("No new value specified!\n");
}
else if(!all(is.numeric(value)) ){
throw("The new value must be a numeric vector!\n");
}
else {
if(missing(index)){
index<-seq(along.with=this$.sites);
}
else {
index<-.checkIndexSanity(this, index);
}
setParameterAtSites(this, process=process, id="omega",value=value,index=index);
return(invisible(TRUE));
}
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: omegaVarM0 - one ratio
##
###########################################################################/**
#
# @RdocMethod omegaVarM0
#
# @title "The M0 (one-ratio) model of variable omega ratios among sites"
#
# \description{
# @get "title".
#
# This method sets the \code{omega} site-process specific parameter
# in the specified range to values sampled from the M0 (one-ratio) model of
# variable omega ratios among sites.
#
# Distribution of omega values:
# \preformatted{
# CATEGORY PROBABILITY
#
# omega 1
# }
# }
#
# \references{
# Yang, Z., Nielsen, R., Goldman, N., Pedersen Krabbe, A-M. (2000) Codon-Substitution Models for Heterogeneous Selection Pressure at Amino Acid Sites - Genetics 155:431-449 \url{http://bit.ly/bvjucn}
#
# Goldman, N., Yang, Z. (1994) A codon-based model of nucleotide substitution for protein-coding DNA sequences - Mol Biol Evol 11(5):725-36 \url{http://bit.ly/aSVEoa}
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A CodonSequence object.}
# \item{process}{A process object inheriting from GY94.}
# \item{omega}{The fixed omega value.}
# \item{index}{A vector of positions.}
# \item{...}{Not used.}
# }
#
# \value{
# Invisible TRUE.
# }
#
# \examples{
# # create a GY94 object
# p<-GY94(kappa=2)
# # create a CodonSequence object, attach process p
# s<-CodonSequence(length=20, processes=list(list(p)))
# # sample states
# sampleStates(s)
# # sample omegas in range 1:5 from model M0
# omegaVarM0(s,p,omega=2,1:5)
# # get omega values
# getOmegas(s,p)
# # get a histogram of omega values in range 1:5
# omegaHist(s,p,breaks=50,1:5)
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"omegaVarM0",
class="CodonSequence",
function(
this,
process,
omega,
index,
...
){
if(missing(process)){
throw("No process specified!\n");
}
if(!is.GY94(process)){
throw("The sepcified process is not a GY94 codon substitution process!\n");
}
else if(missing(omega)){
throw("No new omega value specified!\n");
}
else if((!is.numeric(omega))| (length(omega) != 1)){
throw("The new value must be a numeric vector of length 1!\n");
}
else {
if(missing(index)){
index<-seq(along.with=this$.sites);
}
else {
index<-.checkIndexSanity(this, index);
}
setParameterAtSites(this, process=process, id="omega",value=omega,index=index);
return(invisible(TRUE));
}
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: omegaVarM1 - neutral
##
###########################################################################/**
#
# @RdocMethod omegaVarM1
#
# @title "The M1 (neutral) model of variable omega ratios among sites"
#
# \description{
# @get "title".
#
# This method sets the \code{omega} site-process specific parameter
# in the specified range to values sampled from the M1 (neutral) model of
# variable omega ratios among sites.
#
# Distribution of omega values:
# \preformatted{
# CATEGORY PROBABILITY
#
# omega_0 = 0 p0
# omega_1 = 1 1-p0
# }
# }
#
# \references{
# Yang, Z., Nielsen, R., Goldman, N., Pedersen Krabbe, A-M. (2000) Codon-Substitution Models for Heterogeneous Selection Pressure at Amino Acid Sites - Genetics 155:431-449 \url{http://bit.ly/bvjucn}
#
# Goldman, N., Yang, Z. (1994) A codon-based model of nucleotide substitution for protein-coding DNA sequences - Mol Biol Evol 11(5):725-36 \url{http://bit.ly/aSVEoa}
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A CodonSequence object.}
# \item{process}{A process object inheriting from GY94.}
# \item{p0}{See above.}
# \item{index}{A vector of positions.}
# \item{...}{Not used.}
# }
#
# \value{
# Invisible TRUE.
# }
#
# \examples{
# # create a GY94 object
# p<-GY94(kappa=2)
# # create a CodonSequence object, attach process p
# s<-CodonSequence(length=25, processes=list(list(p)))
# # sample states
# sampleStates(s)
# # sample omegas in range 1:20 from model M1
# omegaVarM1(s,p,p0=0.5,1:20)
# # get omega values
# getOmegas(s,p)
# # get a histogram of omega values in range 1:20
# omegaHist(s,p,breaks=50,1:20)
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"omegaVarM1",
class="CodonSequence",
function(
this,
process,
p0,
index,
...
){
if(missing(process)){
throw("No process specified!\n");
}
if(!is.GY94(process)){
throw("The sepcified process is not a GY94 codon substitution process!\n");
}
else if(missing(p0)){
throw("No p0 value specified!\n");
}
else if((!is.numeric(p0))| (length(p0) != 1)){
throw("The new value must be a numeric vector of length 1!\n");
}
else if(p0 < 0 | p0 > 1){
throw("The p0 parameter must be in the [0,1] interval!\n");
}
else {
if(missing(index)){
index<-seq(along.with=this$.sites);
}
else {
index<-.checkIndexSanity(this, index);
}
for(site in this$.sites[index]){
setParameterAtSite(this=process,site=site, id="omega", value=sample(c(0,1), size=1, replace=FALSE, prob=c(p0,(1-p0))));
}
return(invisible(TRUE));
}
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: omegaVarM2 - selection
##
###########################################################################/**
#
# @RdocMethod omegaVarM2
#
# @title "The M2 (selection) model of variable omega ratios among sites"
#
# \description{
# @get "title".
#
# This method sets the \code{omega} site-process specific parameter
# in the specified range to values sampled from the M2 (selection) model of
# variable omega ratios among sites.
#
# Distribution of omega values:
# \preformatted{
# CATEGORY PROBABILITY
#
# omega_0 = 0 p0
# omega_1 = 1 p1
# omega_2 1-p0-p1
# }
# }
#
# \references{
# Yang, Z., Nielsen, R., Goldman, N., Pedersen Krabbe, A-M. (2000) Codon-Substitution Models for Heterogeneous Selection Pressure at Amino Acid Sites - Genetics 155:431-449 \url{http://bit.ly/bvjucn}
#
# Goldman, N., Yang, Z. (1994) A codon-based model of nucleotide substitution for protein-coding DNA sequences - Mol Biol Evol 11(5):725-36 \url{http://bit.ly/aSVEoa}
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A CodonSequence object.}
# \item{process}{A process object inheriting from GY94.}
# \item{p0}{See above.}
# \item{p1}{See above.}
# \item{omega_2}{See above.}
# \item{index}{A vector of positions.}
# \item{...}{Not used.}
# }
#
# \value{
# Invisible TRUE.
# }
#
# \examples{
# # create a GY94 object
# p<-GY94(kappa=2)
# # create a CodonSequence object, attach process p
# s<-CodonSequence(length=25, processes=list(list(p)))
# # sample states
# sampleStates(s)
# # sample omegas in range 1:20 from model M2
# omegaVarM2(s,p,p0=0.2,p1=0.3,omega_2=4,1:20)
# # get omega values
# getOmegas(s,p)
# # get a histogram of omega values in range 1:20
# omegaHist(s,p,breaks=50,1:20)
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"omegaVarM2",
class="CodonSequence",
function(
this,
process,
p0,
p1,
omega_2,
index,
...
){
if(missing(process)){
throw("No process specified!\n");
}
if(!is.GY94(process)){
throw("The sepcified process is not a GY94 codon substitution process!\n");
}
else if(missing(p0)){
throw("No p0 value specified!\n");
}
else if((!is.numeric(p0))| (length(p0) != 1)){
throw("The p0 value must be a numeric vector of length 1!\n");
}
else if(p0 < 0 | p0 > 1){
throw("The p0 parameter must be in the [0,1] interval!\n");
}
else if(missing(p1)){
throw("No p1 value specified!\n");
}
else if((!is.numeric(p1))| (length(p1) != 1)){
throw("The p1 value must be a numeric vector of length 1!\n");
}
else if(p1 < 0 | p1 > 1){
throw("The p1 parameter must be in the [0,1] interval!\n");
}
else if(missing(omega_2)){
throw("No omega value specified!\n");
}
else if((!is.numeric(omega_2))| (length(omega_2) != 1)){
throw("The omega value must be a numeric vector of length 1!\n");
}
else {
if(missing(index)){
index<-seq(along.with=this$.sites);
}
else {
index<-.checkIndexSanity(this, index);
}
for(site in this$.sites[index]){
setParameterAtSite(this=process,site=site, id="omega", value=sample(c(0,1,omega_2), size=1, replace=FALSE, prob=c(p0,p1,(1-p0-p1))));
}
return(invisible(TRUE));
}
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: omegaVarM3 - discrete
##
###########################################################################/**
#
# @RdocMethod omegaVarM3
#
# @title "The M3 (discrete) model of variable omega ratios among sites"
#
# \description{
# @get "title".
#
# This method sets the \code{omega} site-process specific parameter
# in the specified range to values sampled from the M3 (discrete) model of
# variable omega ratios among sites.
#
# Distribution of omega values:
# \preformatted{
# CATEGORY PROBABILITY
#
# omega_0 p0
# omega_1 p1
# omega_2 p2
# ... ...
# omega_k pk
# }
# }
#
# \references{
# Yang, Z., Nielsen, R., Goldman, N., Pedersen Krabbe, A-M. (2000) Codon-Substitution Models for Heterogeneous Selection Pressure at Amino Acid Sites - Genetics 155:431-449 \url{http://bit.ly/bvjucn}
#
# Goldman, N., Yang, Z. (1994) A codon-based model of nucleotide substitution for protein-coding DNA sequences - Mol Biol Evol 11(5):725-36 \url{http://bit.ly/aSVEoa}
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A CodonSequence object.}
# \item{process}{A process object inheriting from GY94.}
# \item{omegas}{A vector of omega values (omega_0 ... omega_k).}
# \item{probs}{A vector of probabilities (p0 ... pk).}
# \item{index}{A vector of positions.}
# \item{...}{Not used.}
# }
#
# \value{
# Invisible TRUE.
# }
#
# \examples{
# # create a GY94 object
# p<-GY94(kappa=2)
# # create a CodonSequence object, attach process p
# s<-CodonSequence(length=25, processes=list(list(p)))
# # sample states
# sampleStates(s)
# # sample omegas in range 1:20 from model M3
# omegaVarM3(s,p,omegas=c(0,2,4),probs=c(1/3,1/3,1/3),1:20)
# # get omega values
# getOmegas(s,p)
# # get a histogram of omega values in range 1:20
# omegaHist(s,p,breaks=50,1:20)
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"omegaVarM3",
class="CodonSequence",
function(
this,
process,
omegas,
probs,
index,
...
){
if(missing(process)){
throw("No process specified!\n");
}
if(!is.GY94(process)){
throw("The sepcified process is not a GY94 codon substitution process!\n");
}
else if(missing(omegas)){
throw("No omega values specified!\n");
}
else if((!is.numeric(omegas))){
throw("The omegas must be numeric!\n");
}
else if(missing(probs)){
throw("No probabilities specified!\n");
}
else if(!is.numeric(probs)){
throw("The probabilities must be numeric!\n");
}
else if(length(omegas) != length(probs)){
throw("The length of the \"omegas\" and \"probs\" vector must be the same!\n");
}
else if(!PSRoot$my.all.equal(sum(probs),1.0)){
probs<-(probs/sum(probs));
warning("The provided probabilities were scaked in order to sum to one!\n");
}
if(missing(index)){
index<-seq(along.with=this$.sites);
}
else {
index<-.checkIndexSanity(this, index);
}
for(site in this$.sites[index]){
setParameterAtSite(this=process,site=site, id="omega", value=sample(omegas, size=1, replace=FALSE, prob=probs));
}
return(invisible(TRUE));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: omegaVarM4 - freqs
##
###########################################################################/**
#
# @RdocMethod omegaVarM4
#
# @title "The M4 (freqs) model of variable omega ratios among sites"
#
# \description{
# @get "title".
#
# This method sets the \code{omega} site-process specific parameter
# in the specified range to values sampled from the M4 (freqs) model of
# variable omega ratios among sites.
#
# Distribution of omega values:
# \preformatted{
# CATEGORY PROBABILITY
#
# omega_0 = 0 p0
# omega_1 = 1/3 p1
# omega_2 = 2/3 p2
# omega_3 = 1 p3
# omega_4 = 3 p4
# }
# }
#
# \references{
# Yang, Z., Nielsen, R., Goldman, N., Pedersen Krabbe, A-M. (2000) Codon-Substitution Models for Heterogeneous Selection Pressure at Amino Acid Sites - Genetics 155:431-449 \url{http://bit.ly/bvjucn}
#
# Goldman, N., Yang, Z. (1994) A codon-based model of nucleotide substitution for protein-coding DNA sequences - Mol Biol Evol 11(5):725-36 \url{http://bit.ly/aSVEoa}
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A CodonSequence object.}
# \item{process}{A process object inheriting from GY94.}
# \item{probs}{A vector of probabilities (p0 ... p4).}
# \item{index}{A vector of positions.}
# \item{...}{Not used.}
# }
#
# \value{
# Invisible TRUE.
# }
#
# \examples{
# # create a GY94 object
# p<-GY94(kappa=2)
# # create a CodonSequence object, attach process p
# s<-CodonSequence(length=25, processes=list(list(p)))
# # sample states
# sampleStates(s)
# # sample omegas in range 1:20 from model M4
# omegaVarM4(s,p,probs=c(2/5,1/5,1/5,1/10,1/10),1:20)
# # get omega values
# getOmegas(s,p)
# # get a histogram of omega values in range 1:20
# omegaHist(s,p,breaks=50,1:20)
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"omegaVarM4",
class="CodonSequence",
function(
this,
process,
probs,
index,
...
){
if(missing(process)){
throw("No process specified!\n");
}
if(!is.GY94(process)){
throw("The sepcified process is not a GY94 codon substitution process!\n");
}
else if(missing(probs)){
throw("No probabilities specified!\n");
}
else if(!is.numeric(probs)){
throw("The probs must be greater than zero!\n");
}
else if( length(probs) != 5){
throw("The length of the \"probs\" vector must be 5!\n");
}
else if(!PSRoot$my.all.equal(sum(probs),1.0)){
probs<-(probs/sum(probs));
warning("The provided probabilities were scaked in order to sum to one!\n");
}
if(missing(index)){
index<-seq(along.with=this$.sites);
}
else {
index<-.checkIndexSanity(this, index);
}
omegas<-c(0,(1/3),(2/3),1,3);
for(site in this$.sites[index]){
setParameterAtSite(this=process,site=site, id="omega", value=sample(omegas, size=1, replace=FALSE, prob=probs));
}
return(invisible(TRUE));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: .omegaVarM5Cont - gamma
##
setMethodS3(
".omegaVarM5Cont",
class="CodonSequence",
function(
this,
process,
alpha,
beta,
index,
...
){
if(missing(process)){
throw("No process specified!\n");
}
if(!is.GY94(process)){
throw("The sepcified process is not a GY94 codon substitution process!\n");
}
else if(missing(alpha)){
throw("No alpha (shape) value specified!\n");
}
else if((!is.numeric(alpha)) | (length(alpha) != 1)){
throw("The alpha (shape) parameter must be a numeric vector of length 1!\n");
}
else if(alpha < 0){
throw("The alpha (shape) must be greater than zero!\n");
}
else if(missing(beta)){
throw("No beta (scale) value specified!\n");
}
else if((!is.numeric(beta)) | (length(beta) != 1)){
throw("The beta (scale) parameter must be a numeric vector of length 1!\n");
}
else if(beta <= 0){
throw("The beta (scale) must be strictly positive!\n");
}
if(missing(index)){
index<-seq(along.with=this$.sites);
}
else {
index<-.checkIndexSanity(this, index);
}
for(site in this$.sites[index]){
setParameterAtSite(this=process,site=site, id="omega", value=rgamma(1,shape=alpha,scale=beta));
}
return(invisible(TRUE));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: .omegaVarM6Cont - 2gamma
##
setMethodS3(
".omegaVarM6Cont",
class="CodonSequence",
function(
this,
process,
p0,
alpha0,
beta0,
alpha1,
beta1,
index,
...
){
if(missing(process)){
throw("No process specified!\n");
}
if(!is.GY94(process)){
throw("The sepcified process is not a GY94 codon substitution process!\n");
}
else if(missing(p0)){
throw("No p0 value specified!\n");
}
else if((!is.numeric(p0)) | (length(p0) != 1)){
throw("The p0 parameter must be a numeric vector of length 1!\n");
}
else if( (p0 < 0) | (p0 > 1)){
throw("p0 must be in the [0,1] interval!\n");
}
else if(missing(alpha0)){
throw("No alpha0 (shape0) value specified!\n");
}
else if((!is.numeric(alpha0)) | (length(alpha0) != 1)){
throw("The alpha0 (shape0) parameter must be a numeric vector of length 1!\n");
}
else if(alpha0 < 0){
throw("The alpha0 (shape0) must be greater than zero!\n");
}
else if(missing(beta0)){
throw("No beta0 (scale0) value specified!\n");
}
else if((!is.numeric(beta0)) | (length(beta0) != 1)){
throw("The beta0 (scale0) parameter must be a numeric vector of length 1!\n");
}
else if(beta0 <= 0){
throw("The beta0 (scale0) must be strictly positive!\n");
}
else if(missing(alpha1)){
throw("No alpha1 (shape1) value specified!\n");
}
else if((!is.numeric(alpha1)) | (length(alpha1) != 1)){
throw("The alpha1 (shape1) parameter must be a numeric vector of length 1!\n");
}
else if(alpha1 < 0){
throw("The alpha1 (shape1) must be greater than zero!\n");
}
else if(missing(beta1)){
throw("No beta1 (scale1) value specified!\n");
}
else if((!is.numeric(beta1)) | (length(beta1) != 1)){
throw("The beta1 (scale1) parameter must be a numeric vector of length 1!\n");
}
else if(beta1 <= 0){
throw("The beta1 (scale1) must be strictly positive!\n");
}
if(missing(index)){
index<-seq(along.with=this$.sites);
}
else {
index<-.checkIndexSanity(this, index);
}
for(site in this$.sites[index]){
setParameterAtSite(this=process,site=site, id="omega", value=sample(c(rgamma(1,shape=alpha0,scale=beta0),rgamma(1,shape=alpha1,scale=beta1)),size=1,replace=FALSE,prob=c(p0,(1-p0)) ));
}
return(invisible(TRUE));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: .omegaVarM7Cont - beta
##
setMethodS3(
".omegaVarM7",
class="CodonSequence",
function(
this,
process,
p,
q,
index,
...
){
if(missing(process)){
throw("No process specified!\n");
}
if(!is.GY94(process)){
throw("The sepcified process is not a GY94 codon substitution process!\n");
}
else if(missing(p)){
throw("No p value specified!\n");
}
else if((!is.numeric(p)) | (length(p) != 1)){
throw("The p parameter must be a numeric vector of length 1!\n");
}
else if(p < 0){
throw("The p parameter must be greater than zero!\n");
}
else if(missing(q)){
throw("No q value specified!\n");
}
else if((!is.numeric(q)) | (length(q) != 1)){
throw("The q parameter must be a numeric vector of length 1!\n");
}
else if(q < 0){
throw("The q parameter must be positive!\n");
}
if(missing(index)){
index<-seq(along.with=this$.sites);
}
else {
index<-.checkIndexSanity(this, index);
}
for(site in this$.sites[index]){
setParameterAtSite(this=process,site=site, id="omega", value=rbeta(1,shape1=p,shape2=q));
}
return(invisible(TRUE));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: .omegaVarM8Cont - beta&omega
##
setMethodS3(
".omegaVarM8Cont",
class="CodonSequence",
function(
this,
process,
p0,
p,
q,
omega,
index,
...
){
if(missing(process)){
throw("No process specified!\n");
}
if(!is.GY94(process)){
throw("The sepcified process is not a GY94 codon substitution process!\n");
}
else if(missing(p0)){
throw("No p0 value specified!\n");
}
else if((!is.numeric(p0)) | (length(p0) != 1)){
throw("The p0 parameter must be a numeric vector of length 1!\n");
}
else if( (p0 < 0) | (p0 > 1)){
throw("The p0 parameter must be from the [0,1] interval!\n");
}
else if(missing(p)){
throw("No p value specified!\n");
}
else if((!is.numeric(p)) | (length(p) != 1)){
throw("The p parameter must be a numeric vector of length 1!\n");
}
else if(p < 0){
throw("The p parameter must be greater than zero!\n");
}
else if(missing(q)){
throw("No q value specified!\n");
}
else if((!is.numeric(q)) | (length(q) != 1)){
throw("The q parameter must be a numeric vector of length 1!\n");
}
else if(q < 0){
throw("The q parameter must be positive!\n");
}
else if(missing(omega)){
throw("No omega value specified!\n");
}
else if((!is.numeric(omega)) | (length(omega) != 1)){
throw("The omega parameter must be a numeric vector of length 1!\n");
}
else if(omega < 0){
throw("The omega parameter must be positive!\n");
}
if(missing(index)){
index<-seq(along.with=this$.sites);
}
else {
index<-.checkIndexSanity(this, index);
}
for(site in this$.sites[index]){
setParameterAtSite(this=process,site=site, id="omega", value=sample(c(rbeta(1,shape1=p,shape2=q),omega),size=1,replace=FALSE,prob=(c(p0,(1-p0)))));
}
return(invisible(TRUE));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: .omegaVarM9 - beta&gamma
##
setMethodS3(
".omegaVarM9Cont",
class="CodonSequence",
function(
this,
process,
p0,
p,
q,
alpha,
beta,
index,
...
){
if(missing(process)){
throw("No process specified!\n");
}
if(!is.GY94(process)){
throw("The sepcified process is not a GY94 codon substitution process!\n");
}
else if(missing(p0)){
throw("No p0 value specified!\n");
}
else if((!is.numeric(p0)) | (length(p0) != 1)){
throw("The p0 parameter must be a numeric vector of length 1!\n");
}
else if( (p0 < 0) | (p0 > 1)){
throw("p0 must be in the [0,1] interval!\n");
}
else if(missing(p)){
throw("No p value specified!\n");
}
else if((!is.numeric(p)) | (length(p) != 1)){
throw("The p parameter must be a numeric vector of length 1!\n");
}
else if(p < 0){
throw("The p must be greater than zero!\n");
}
else if(missing(q)){
throw("No q value specified!\n");
}
else if((!is.numeric(q)) | (length(q) != 1)){
throw("The q parameter must be a numeric vector of length 1!\n");
}
else if(q < 0){
throw("The q must be greater than zero!\n");
}
else if(missing(alpha)){
throw("No alpha (shape) value specified!\n");
}
else if((!is.numeric(alpha)) | (length(alpha) != 1)){
throw("The alpha (shape) parameter must be a numeric vector of length 1!\n");
}
else if(alpha < 0){
throw("The alpha (shape) must be greater than zero!\n");
}
else if(missing(beta)){
throw("No beta (scale) value specified!\n");
}
else if((!is.numeric(beta)) | (length(beta) != 1)){
throw("The beta (scale) parameter must be a numeric vector of length 1!\n");
}
else if(beta <= 0){
throw("The beta (scale) must be strictly positive!\n");
}
if(missing(index)){
index<-seq(along.with=this$.sites);
}
else {
index<-.checkIndexSanity(this, index);
}
for(site in this$.sites[index]){
setParameterAtSite(this=process,site=site, id="omega", value=sample(c(rbeta(1,shape1=p,shape2=q),rgamma(1,shape=alpha,scale=beta)),size=1,replace=FALSE,prob=c(p0,(1-p0)) ));
}
return(invisible(TRUE));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: .omegaVarM10Cont - beta&gamma+1
##
setMethodS3(
".omegaVarM10Cont",
class="CodonSequence",
function(
this,
process,
p0,
p,
q,
alpha,
beta,
index,
...
){
if(missing(process)){
throw("No process specified!\n");
}
if(!is.GY94(process)){
throw("The sepcified process is not a GY94 codon substitution process!\n");
}
else if(missing(p0)){
throw("No p0 value specified!\n");
}
else if((!is.numeric(p0)) | (length(p0) != 1)){
throw("The p0 parameter must be a numeric vector of length 1!\n");
}
else if( (p0 < 0) | (p0 > 1)){
throw("p0 must be in the [0,1] interval!\n");
}
else if(missing(p)){
throw("No p value specified!\n");
}
else if((!is.numeric(p)) | (length(p) != 1)){
throw("The p parameter must be a numeric vector of length 1!\n");
}
else if(p < 0){
throw("The p must be greater than zero!\n");
}
else if(missing(q)){
throw("No q value specified!\n");
}
else if((!is.numeric(q)) | (length(q) != 1)){
throw("The q parameter must be a numeric vector of length 1!\n");
}
else if(q < 0){
throw("The q must be greater than zero!\n");
}
else if(missing(alpha)){
throw("No alpha (shape) value specified!\n");
}
else if((!is.numeric(alpha)) | (length(alpha) != 1)){
throw("The alpha (shape) parameter must be a numeric vector of length 1!\n");
}
else if(alpha < 0){
throw("The alpha (shape) must be greater than zero!\n");
}
else if(missing(beta)){
throw("No beta (scale) value specified!\n");
}
else if((!is.numeric(beta)) | (length(beta) != 1)){
throw("The beta (scale) parameter must be a numeric vector of length 1!\n");
}
else if(beta <= 0){
throw("The beta (scale) must be strictly positive!\n");
}
if(missing(index)){
index<-seq(along.with=this$.sites);
}
else {
index<-.checkIndexSanity(this, index);
}
# It's not too elegant to fork the whole method for just shifting the gamma with 1...
for(site in this$.sites[index]){
setParameterAtSite(this=process,site=site, id="omega", value=sample(c(rbeta(1,shape1=p,shape2=q),(1 + rgamma(1,shape=alpha,scale=beta)) ),size=1,replace=FALSE,prob=c(p0,(1-p0)) ));
}
return(invisible(TRUE));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: .omegaVarM11Cont - beta&normal>1
##
setMethodS3(
".omegaVarM11Cont",
class="CodonSequence",
function(
this,
process,
p0,
p,
q,
mean,
sd,
index,
...
){
if(missing(process)){
throw("No process specified!\n");
}
if(!is.GY94(process)){
throw("The sepcified process is not a GY94 codon substitution process!\n");
}
else if(missing(p0)){
throw("No p0 value specified!\n");
}
else if((!is.numeric(p0)) | (length(p0) != 1)){
throw("The p0 parameter must be a numeric vector of length 1!\n");
}
else if( (p0 < 0) | (p0 > 1)){
throw("p0 must be in the [0,1] interval!\n");
}
else if(missing(p)){
throw("No p value specified!\n");
}
else if((!is.numeric(p)) | (length(p) != 1)){
throw("The p parameter must be a numeric vector of length 1!\n");
}
else if(p < 0){
throw("The p must be greater than zero!\n");
}
else if(missing(q)){
throw("No q value specified!\n");
}
else if((!is.numeric(q)) | (length(q) != 1)){
throw("The q parameter must be a numeric vector of length 1!\n");
}
else if(q < 0){
throw("The q must be greater than zero!\n");
}
else if(missing(mean)){
throw("No mean specified!\n");
}
else if((!is.numeric(mean)) | (length(mean) != 1)){
throw("The mean parameter must be a numeric vector of length 1!\n");
}
else if(missing(sd)){
throw("No sd value specified!\n");
}
else if((!is.numeric(sd)) | (length(sd) != 1)){
throw("The sd parameter must be a numeric vector of length 1!\n");
}
else if( sd < 0){
throw("The sd parameter must be positive!\n");
}
if(missing(index)){
index<-seq(along.with=this$.sites);
}
else {
index<-.checkIndexSanity(this, index);
}
rnorm.gt.1<-function(mean=NA,sd=NA){
# This is probably the most primitive way to truncate the distribution!
tmp<-rnorm(1,mean=mean,sd=sd);
while( tmp <= 1){
tmp<-rnorm(1,mean=mean,sd=sd);
};
return(tmp);
}
for(site in this$.sites[index]){
setParameterAtSite(this=process,site=site, id="omega", value=sample(c(rbeta(1,shape1=p,shape2=q), (rnorm.gt.1(mean=mean,sd=sd))),size=1,replace=FALSE,prob=c(p0,(1-p0)) ));
}
return(invisible(TRUE));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: .omegaVarM12Cont - 0&2normal>0
##
setMethodS3(
".omegaVarM12Cont",
class="CodonSequence",
function(
this,
process,
p0,
p1,
sd1,
mean2,
sd2,
index,
...
){
if(missing(process)){
throw("No process specified!\n");
}
if(!is.GY94(process)){
throw("The sepcified process is not a GY94 codon substitution process!\n");
}
else if(missing(p0)){
throw("No p0 value specified!\n");
}
else if((!is.numeric(p0)) | (length(p0) != 1)){
throw("The p0 parameter must be a numeric vector of length 1!\n");
}
else if( (p0 < 0) | (p0 > 1)){
throw("p0 must be in the [0,1] interval!\n");
}
else if(missing(p1)){
throw("No p1 value specified!\n");
}
else if((!is.numeric(p1)) | (length(p1) != 1)){
throw("The p1 parameter must be a numeric vector of length 1!\n");
}
else if( (p1 < 0) | (p1 > 1)){
throw("p1 must be in the [0,1] interval!\n");
}
else if(missing(sd1)){
throw("No sd1 value specified!\n");
}
else if((!is.numeric(sd1)) | (length(sd1) != 1)){
throw("The sd1 parameter must be a numeric vector of length 1!\n");
}
else if( sd1 < 0){
throw("The sd1 parameter must be positive!\n");
}
else if(missing(mean2)){
throw("No mean2 specified!\n");
}
else if((!is.numeric(mean2)) | (length(mean2) != 1)){
throw("The mean2 parameter must be a numeric vector of length 1!\n");
}
else if(missing(sd2)){
throw("No sd2 value specified!\n");
}
else if((!is.numeric(sd2)) | (length(sd2) != 1)){
throw("The sd2 parameter must be a numeric vector of length 1!\n");
}
else if( sd2 < 0){
throw("The sd2 parameter must be positive!\n");
}
if(missing(index)){
index<-seq(along.with=this$.sites);
}
else {
index<-.checkIndexSanity(this, index);
}
rnorm.gt.0<-function(mean=NA,sd=NA){
# This is probably the most primitive way to truncate the distribution!
tmp<-rnorm(1,mean=mean,sd=sd);
while( tmp <= 0){
tmp<-rnorm(1,mean=mean,sd=sd);
};
return(tmp);
}
for(site in this$.sites[index]){
setParameterAtSite(this=process,site=site, id="omega", value=sample(
c(0,rnorm.gt.0(mean=1,sd=sd1),rnorm.gt.0(mean=mean2,sd=sd2)),
size=1,
replace=FALSE,
prob=c(p0,((1-p0)*p1),((1-p0)*(1-p1)))
));
}
return(invisible(TRUE));
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: .omegaVarM13Cont - 3normal>0
##
setMethodS3(
".omegaVarM13Cont",
class="CodonSequence",
function(
this,
process,
p0,
p1,
sd0,
sd1,
mean2,
sd2,
index,
...
){
if(missing(process)){
throw("No process specified!\n");
}
if(!is.GY94(process)){
throw("The sepcified process is not a GY94 codon substitution process!\n");
}
else if(missing(p0)){
throw("No p0 value specified!\n");
}
else if((!is.numeric(p0)) | (length(p0) != 1)){
throw("The p0 parameter must be a numeric vector of length 1!\n");
}
else if( (p0 < 0) | (p0 > 1)){
throw("p0 must be in the [0,1] interval!\n");
}
else if(missing(p1)){
throw("No p1 value specified!\n");
}
else if((!is.numeric(p1)) | (length(p1) != 1)){
throw("The p1 parameter must be a numeric vector of length 1!\n");
}
else if( (p1 < 0) | (p1 > 1)){
throw("p1 must be in the [0,1] interval!\n");
}
else if(missing(sd0)){
throw("No sd0 value specified!\n");
}
else if((!is.numeric(sd0)) | (length(sd0) != 1)){
throw("The sd0 parameter must be a numeric vector of length 1!\n");
}
else if( sd0 < 0){
throw("The sd0 parameter must be positive!\n");
}
else if(missing(sd1)){
throw("No sd1 value specified!\n");
}
else if((!is.numeric(sd1)) | (length(sd1) != 1)){
throw("The sd1 parameter must be a numeric vector of length 1!\n");
}
else if( sd1 < 0){
throw("The sd1 parameter must be positive!\n");
}
else if(missing(mean2)){
throw("No mean2 specified!\n");
}
else if((!is.numeric(mean2)) | (length(mean2) != 1)){
throw("The mean2 parameter must be a numeric vector of length 1!\n");
}
else if(missing(sd2)){
throw("No sd2 value specified!\n");
}
else if((!is.numeric(sd2)) | (length(sd2) != 1)){
throw("The sd2 parameter must be a numeric vector of length 1!\n");
}
else if( sd2 < 0){
throw("The sd2 parameter must be positive!\n");
}
if(missing(index)){
index<-seq(along.with=this$.sites);
}
else {
index<-.checkIndexSanity(this, index);
}
rnorm.gt.0<-function(mean=NA,sd=NA){
# This is probably the most primitive way to truncate the distribution!
tmp<-rnorm(1,mean=mean,sd=sd);
while( tmp <= 0){
tmp<-rnorm(1,mean=mean,sd=sd);
};
return(tmp);
}
for(site in this$.sites[index]){
setParameterAtSite(this=process,site=site, id="omega", value=sample(
c(rnorm.gt.0(mean=0,sd=sd0),rnorm.gt.0(mean=1,sd=sd1),rnorm.gt.0(mean=mean2,sd=sd2)),
size=1,
replace=FALSE,
prob=c(p0,((1-p0)*p1),((1-p0)*(1-p1)))
));
}
return(invisible(TRUE));
},
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{...}{Not used.}
# }
#
#
# \value{
# Returns an invisible TRUE if no inconsistencies found in the object, throws
# an error otherwise.
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"checkConsistency",
class="CodonSequence",
function(
this,
...
){
wp<-this$writeProtected;
if (wp) {
this$writeProtected<-FALSE;
}
may.fail<-function(this) {
for(pos in seq(along.with=this$.sites)){
if(!is.CodonAlphabet(this$.sites[[pos]]$.alphabet)){
throw("The alphabet attached to site ",pos," is not a codon alphabet!\n");
}
}
}
tryCatch(may.fail(this),finally=this$writeProtected<-wp);
NextMethod();
},
private=FALSE,
protected=FALSE,
overwrite=FALSE,
conflict="warning",
validators=getOption("R.methodsS3:validators:setMethodS3")
);
##
## Method: Translate
##
###########################################################################/**
#
# @RdocMethod Translate
#
# @title "Translate a CodonSequence object"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{this}{A CodonSequence object}
# \item{...}{Not used.}
# }
#
# \value{
# The translation as an AminoAcidSequence object.
# }
#
# \examples{
# # create a CodonSequence object
# s<-CodonSequence(string="ATCTTTCGAATGGGGCCCTCCCGA")
# # get the translation as an AminoAcidSequence object
# as<-Translate(s)
# as
# }
#
# @author
#
# \seealso{
# @seeclass
# }
#
#*/###########################################################################
setMethodS3(
"Translate",
class="CodonSequence",
function(
this,
...
){
# Create an amino acid sequence of proper length:
that<-AminoAcidSequence(
name=paste("Translation of",this$name),
length=this$length
);
# Setting the states, using the site alphabet to translate codons for greater flexibility.
for (pos in seq(along.with=this$.sites)){
setState(that$.sites[[pos]], translateCodon(this$.sites[[pos]]$alphabet, this$.sites[[pos]]$state) );
}
return(that);
},
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.