/PhyloSim.R
R | 6741 lines | 2944 code | 500 blank | 3297 comment | 491 complexity | 599a6b1f6c5fd46221d5d716d177dc82 MD5 | raw file
Large files files are truncated, but you can click here to view the full file
1## 2## Copyright 2009 Botond Sipos 3## See the package description for licensing information. 4## 5##########################################################################/** 6# 7# @RdocClass PhyloSim 8# \alias{phylosim} 9# 10# @title "The PhyloSim class" 11# 12# \description{ 13# 14# PhyloSim is an extensible object-oriented framework for the Monte Carlo simulation 15# of sequence evolution written in 100 percent \code{R}. 16# It is built on the top of the \code{\link{R.oo}} and \code{\link{ape}} packages and uses 17# Gillespie's direct method to simulate substitutions, insertions and deletions. 18# 19# Key features offered by the framework: 20# \itemize{ 21# \item Simulation of the evolution of a set of discrete characters with arbitrary states evolving 22# by a continuous-time Markov process with an arbitrary rate matrix. 23# \item Explicit implementations of the most popular substitution models (for nucleotides, amino acids and codons). 24# \item Simulation under the popular models of among-sites rate variation, like the gamma (+G) and invariants plus gamma (+I+G) models. 25# \item The possibility to simulate with arbitrarily complex patterns of among-sites rate variation by setting the site specific rates according to any \code{R} expression. 26# \item Simulation with one or more separate insertion and/or deletion processes acting on the sequences, which sample the insertion/deletion length from an arbitrary discrete distribution or an \code{R} expression (so all the probability distributions implemented in \code{R} are readily available for this purpose). 27# \item Simulation of the effects of variable functional constraints over the sites by site-process-specific insertion and deletion tolerance parameters, which determine the rejection probability of a proposed insertion/deletion. 28# \item The possibility of having a different set of processes and site-process-specific parameters for every site, which allow for an arbitrary number of partitions in the simulated data. 29# \item Simulation of heterotachy and other cases of non-homogeneous evolution by allowing the user to set "node hook" functions altering the site properties at internal nodes of the phylogeny. 30# \item The possibility to export the counts of various events ("branch statistics") as phylo objects (see \code{\link{exportStatTree.PhyloSim}}). 31# } 32# 33# General notes: 34# \itemize{ 35# \item The \code{Sequence} objects have no "immortal links". The simulation 36# is aborted if the sequence length shrinks to zero. It is up to the user 37# to choose sensible indel rates and sequence lengths to prevent that. 38# \item The sites near the beginning and end of the sequences have less sites proposing 39# insertion and deletion events around the so the insertion and deletion processes 40# have an "edge effect". The user can simulate 41# realistic flanking sequences to alleviate the edge effect in the simulation settings where 42# it may be an issue. 43# } 44# 45# Notes on performance: 46# \itemize{ 47# \item The pure \code{R} implementation offers felxibility, but also comes 48# with a slower simulation speed. If the \code{PSIM_FAST} object is present in the environment, a "fast & careless" 49# mode is enabled. In this mode most of the error checking is skipped, increasing the speed. 50# It is recomended that simulations are only run in fast mode if you are sure that the simulation 51# settings are free from errors. It is probably a good practice to set up the simulations in normal mode 52# with short sequences and enable fast mode when running the actual simulation with long sequences. 53# \item Please note, that no "branch statistics" are saved in fast mode. 54# \item Logging also has a negative impact on performance, so it's not a good idea to run 55# large simulations with the logging enabled. 56# \item The time needed to run a simulation depends not only on the number of the sites, 57# but also on the length of the tree. 58# \item Constructing \code{Sequence} objects with large number of sites is expensive. Avoid doing 59# that inside a cycle. 60# \item In the case of \code{Sequence} objects with a large number of sites (more than 10 000) the 61# amount of available memory can be limiting as well. 62# } 63# 64# The examples below demonstrate only some more common simulation settings, 65# the framework offers much more flexibility. See the package 66# vignette (\code{vignette("PhyloSim",package="phylosim")}) and the 67# examples directory (\url{http://github.com/sbotond/phylosim/tree/master/examples/}) 68# for additional examples. 69# 70# @classhierarchy 71# } 72# 73# 74# \references{ 75# Gillespie, DT (1977) Exact stochastic simulation of coupled chemical reactions - 76# J. Phys. Chem. 81 (25):2340-2361 \url{http://dx.doi.org/10.1021/j100540a008} 77# } 78# 79# @synopsis 80# 81# \arguments{ 82# \item{phylo}{A rooted phylo object, constructed by the APE package.} 83# \item{root.seq}{A valid Sequence object with Process objects attached. Used as the starting sequence during simulation.} 84# \item{name}{The name of the object (a character vector of length one).} 85# \item{log.file}{Name of the file used for logging.} 86# \item{log.level}{An integer specifying the verbosity of logging (see \code{\link{setLogLevel.PhyloSim}}).} 87# \item{...}{Not used.} 88# } 89# 90# \section{Fields and Methods}{ 91# @allmethods 92# } 93# 94# \examples{ 95# set.seed(1) 96# ## The following examples demonstrate 97# ## the typical use of the framework. 98# ## See the package vignette and 99# ## \url{http://github.com/sbotond/phylosim/tree/master/examples/} 100# 101# ## The ll() method gives information about the methods defined 102# ## in the immediate class of an object. 103# ## Useful when exploring the framework. 104# 105# s<-Sequence() 106# ll(s) 107# ll(PhyloSim()) 108# ll(GTR()) 109# 110# ## Example 1 - A short simulation: 111# ## simulate nucleotide seqeunces and display 112# ## the resulting alignment matrix. 113# 114# Simulate( 115# PhyloSim(phy=rcoal(3), 116# root=NucleotideSequence(string="ATGC", proc=list(list(JC69())) ) ) 117# )$alignment 118# 119# # Construct a phylo object for the following 120# # simulations, scale total tree length to 1: 121# 122# tmp<-PhyloSim(phylo=rcoal(3)) 123# scaleTree(tmp,1/tmp$treeLength) 124# tmp$treeLength 125# t<-tmp$phylo 126# 127# ## Example 3 - simulating rate variation, 128# ## insertions and deletions. 129# ## See the examples/example_3_clean.R file 130# ## in the phylosim GitHub repository. 131# 132# # construct a GTR process object 133# gtr<-GTR( 134# name="MyGTR", 135# rate.params=list( 136# "a"=1, "b"=2, "c"=3, 137# "d"=1, "e"=2, "f"=3 138# ), 139# base.freqs=c(2,2,1,1)/6 140# ) 141# # get object summary 142# summary(gtr) 143# # display a bubble plot 144# plot(gtr) 145# 146# # construct root sequence object 147# s<-NucleotideSequence(length=20) 148# # attach process via virtual field 149# s$processes<-list(list(gtr)) 150# # sample states from the equilibrium 151# # distribution of the attached processes 152# sampleStates(s) 153# # create among-site rate variation by sampling 154# # the "rate.multiplier" site-process-specific parameter 155# # from a discrete gamma distribution (GTR+G). 156# plusGamma(s,gtr,shape=0.1) 157# # make the range 11:12 invariable 158# setRateMultipliers(s,gtr,0,11:12) 159# # get the rate multipliers for s and gtr 160# getRateMultipliers(s,gtr) 161# 162 # construct a deletion process object 163# # proposing lengths in the range 1:3 164# d<-DiscreteDeletor( 165# rate=0.1, 166# name="MyDel", 167# sizes=c(1:3), 168# probs=c(3/6,2/6,1/6) 169# ) 170# # get object 171# summary(d) 172# # plot deletion length distribution 173# plot(d) 174# # attach deletion process d to sequence s 175# attachProcess(s,d) 176# # create a region rejecting all deletions 177# setDeletionTolerance(s,d,0,11:12) 178# 179# # construct an insertion process object 180# # proposing lengths in the range 1:3 181# i<-DiscreteInsertor( 182# rate=0.1, 183# name="MyDel", 184# sizes=c(1:2), 185# probs=c(1/2,1/2), 186# template.seq=NucleotideSequence(length=1,processes=list(list(JC69()))) 187# ) 188# # states will be sampled from the JC69 equilibrium distribution 189# # get object 190# summary(i) 191# # plot insertion length distribution 192# plot(i) 193# # attach insertion process i to sequence s 194# attachProcess(s,i) 195# # create a region rejecting all insertions 196# setInsertionTolerance(s,i,0,11:12) 197# 198# # plot total site rates 199# plot(s) 200# # construct simulation object 201# sim<-PhyloSim(root.seq=s, phylo=t) 202# # get object summary 203# summary(sim) 204# # plot tree 205# plot(sim) 206# # run simulation 207# Simulate(sim) 208# # get the list of recorded per-branch event counts 209# getBranchEvents(sim) 210# # export the number of substitutions as a phylo object 211# subst<-exportStatTree(sim,"substitution") 212# # plot the exported phylo object 213# plot(subst) 214# # plot tree and alignment 215# plot(sim) 216# # save and display alingment 217# file<-paste("PhyloSim_dummy_fasta_",Sys.getpid(),".fas",sep=""); 218# saveAlignment(sim,file=file); 219# # print out the Fasta file 220# cat(paste(scan(file=file,what=character(),sep="\n"),collapse="\n"));cat("\n"); 221# # delete Fasta file 222# unlink(file); 223# 224# # See \url{http://github.com/sbotond/phylosim/tree/master/examples/examples_class.R} 225# # for the full list of PhyloSim constructor examples. 226# 227# ## See the package vignette and 228# ## the GitHub repository for even more examples. 229# } 230# 231# @author 232# 233# \seealso{ 234# \code{\link{PSRoot} \link{Alphabet} \link{AminoAcidAlphabet} 235# \link{AminoAcidSequence} \link{AminoAcidSubst} 236# \link{AnyAlphabet} \link{BinaryAlphabet} \link{BinarySequence} 237# \link{BinarySubst} \link{BrownianInsertor} \link{CodonAlphabet} 238# \link{CodonSequence} \link{CodonUNREST} \link{ContinuousDeletor} 239# \link{ContinuousInsertor} \link{cpREV} \link{DiscreteDeletor} 240# \link{DiscreteInsertor} \link{Event} \link{F81} \link{F84} 241# \link{FastFieldDeletor} \link{GeneralDeletor} 242# \link{GeneralInDel} \link{GeneralInsertor} \link{GeneralSubstitution} 243# \link{GTR} \link{GY94} \link{HKY} \link{JC69} \link{JTT} \link{JTT.dcmut} 244# \link{K80} \link{K81} \link{LG} \link{mtArt} \link{mtMam} \link{mtREV24} 245# \link{MtZoa} \link{NucleotideAlphabet} \link{NucleotideSequence} \link{PAM} 246# \link{PAM.dcmut} \link{PhyloSim} \link{Process} \link{QMatrix} \link{Sequence} 247# \link{Site} \link{T92} \link{TN93} \link{UNREST} \link{WAG}} 248# } 249# 250#*/########################################################################### 251setConstructorS3( 252"PhyloSim", 253 function( 254 phylo=NA, 255 root.seq=NA, 256 name=NA, 257 log.file=NA, 258 log.level=-1, # no loggin is performed by default 259 ... 260 ) { 261 262 this<-PSRoot(); 263 this<-extend(this, 264 "PhyloSim", 265 .name="Anonymous", 266 .phylo=NA, 267 .root.sequence=NA, 268 .sequences=list(), # references to the sequence objects 269 .node.hooks=list(), # references to the node hook functions. 270 .branch.stats=list(), # branch statistics. 271 .alignment=NA, # the resulting alignment in fasat format. 272 .log.file=NA, # the name of the log file. 273 .log.connection=NA, # connection for the log file. 274 .log.level=NA # log level 275 ); 276 277 if(!all(is.na(phylo))){ 278 this$phylo<-phylo; 279 } 280 281 if(!all(is.na(root.seq))){ 282 this$rootSeq<-root.seq; 283 } 284 285 if(!missing(name)){ 286 this$name<-name; 287 } 288 289 if(!missing(log.file)){ 290 this$logFile<-log.file; 291 } else { 292 # Setting default log file: 293 tmp<-this$id; 294 tmp<-gsub(":","_",tmp); 295 this$logFile<-paste(tmp,".log",sep=""); 296 } 297 298 # Setting log level: 299 this$logLevel<-log.level; 300 301 return(this); 302 303 }, 304 enforceRCC=TRUE 305); 306 307## 308## Method: checkConsistency 309## 310###########################################################################/** 311# 312# @RdocMethod checkConsistency 313# 314# @title "Check object consistency" 315# 316# \description{ 317# @get "title". 318# } 319# 320# @synopsis 321# 322# \arguments{ 323# \item{this}{An object.} 324# \item{...}{Not used.} 325# } 326# 327# 328# \value{ 329# Returns an invisible TRUE if no inconsistencies found in the object, throws 330# an error otherwise. 331# } 332# 333# @author 334# 335# \seealso{ 336# @seeclass 337# } 338# 339#*/########################################################################### 340setMethodS3( 341 "checkConsistency", 342 class="PhyloSim", 343 function( 344 this, 345 ... 346 ){ 347 348 may.fail<-function(this) { 349 350 # Checking the name: 351 this$name<-this$name; 352 # Checking the phylo object: 353 if (!any(is.na(this$.phylo)) & !is.phylo(this$.phylo) ){ 354 throw("The phylo object is invalid!\n"); 355 } 356 # Checking the log level: 357 if(!is.numeric(this$.log.level) | (length(this$.log.level) != 1) ){ 358 throw("The log level must be numeric vector of length 1!\n"); 359 } 360 # Checking lof file: 361 if(!is.character(this$.log.file) | (length(this$.log.level) != 1) ){ 362 throw("The log file must be charcter vector of length 1!\n"); 363 } 364 # Checking the sequences: 365 for (seq in this$.sequences){ 366 if(is.Sequence(seq)){ 367 checkConsistency(seq); 368 } 369 } 370 # Checking node hooks: 371 for (hook in this$.node.hooks){ 372 if(!is.null(hook) & !is.function(hook)){ 373 throw("Invalid node hook found!\n"); 374 } 375 } 376 # Checking the alignment: 377 if(!any(is.na(this$.alignment))){ 378 .checkAlignmentConsistency(this, this$.alignment); 379 } 380 381 } 382 tryCatch(may.fail(this)); 383 return(invisible(TRUE)); 384 385 }, 386 private=FALSE, 387 protected=FALSE, 388 overwrite=FALSE, 389 conflict="warning", 390 validators=getOption("R.methodsS3:validators:setMethodS3") 391) 392 393## 394## Method: is.phylo.default 395## 396###########################################################################/** 397# 398# @RdocDefault is.phylo 399# 400# @title "Check if an object is an instance of the phylo class" 401# 402# \description{ 403# @get "title". 404# Phylo objects are created by the \pkg{APE} package. This method just return the value of \code{inherits(this,"phylo")}. 405# } 406# 407# @synopsis 408# 409# \arguments{ 410# \item{this}{An object.} 411# \item{...}{Not used.} 412# } 413# 414# \value{ 415# TRUE or FALSE. 416# } 417# 418# \examples{ 419# # load APE 420# library(ape); 421# # create some objects 422# o1<-Object(); 423# o2<-rcoal(3); 424# # check if they are phylo objects 425# is.phylo(o1); 426# is.phylo(o2); 427# 428# } 429# 430# @author 431# 432# \seealso{ 433# The \pkg{ape} package. 434# } 435# 436#*/########################################################################### 437setMethodS3( 438 "is.phylo", 439 class="default", 440 function( 441 this, 442 ... 443 ){ 444 445 inherits(this,"phylo"); 446 447 }, 448 private=FALSE, 449 protected=FALSE, 450 overwrite=FALSE, 451 conflict="warning", 452 validators=getOption("R.methodsS3:validators:setMethodS3") 453) 454 455## 456## Method: setPhylo 457## 458###########################################################################/** 459# 460# @RdocMethod setPhylo 461# 462# @title "Set the phylo object for a PhyloSim object" 463# 464# \description{ 465# @get "title". 466# 467# The internal structure of the provided phylo object is reordered in a cladeweise fashion. 468# } 469# 470# @synopsis 471# 472# \arguments{ 473# \item{this}{A PhyloSim object.} 474# \item{value}{A phylo object created by the \pkg{ape} package.} 475# \item{...}{Not used.} 476# } 477# 478# \value{ 479# A phylo object or FALSE. 480# } 481# 482# \examples{ 483# #create a PhyloSim object 484# sim<-PhyloSim(); 485# # creat a phylo object 486# tree<-rcoal(3); 487# # get/set phylo object 488# setPhylo(sim,tree); 489# getPhylo(sim,tree); 490# # get/set phylo object via virtual field 491# sim$tree<-rcoal(5); 492# sim$tree; 493# } 494# 495# @author 496# 497# \seealso{ 498# The PhyloSim class, the \pkg{ape} package. 499# } 500# 501#*/########################################################################### 502setMethodS3( 503 "setPhylo", 504 class="PhyloSim", 505 function( 506 this, 507 value, 508 ... 509 ){ 510 511 if(missing(value)){ 512 throw("No object provided!\n"); 513 } 514 else if(!is.phylo(value)){ 515 throw("The new value must be a \"phylo\" object!\n"); 516 } 517 else if(!is.rooted(value)){ 518 throw("The new value must be a rooted \"phylo\" object!\n"); 519 } 520 else { 521 522 .checkTipLabels(value); 523 this$.phylo<-value; 524 this$.phylo<-reorder(this$.phylo, order="cladewise"); 525 for (i in this$nodes){ 526 this$.sequences[[i]]<-NA; 527 } 528 return(this$.phylo); 529 530 } 531 return(FALSE); 532 533 }, 534 private=FALSE, 535 protected=FALSE, 536 overwrite=FALSE, 537 conflict="warning", 538 validators=getOption("R.methodsS3:validators:setMethodS3") 539); 540 541## 542## Method: .checkTipLabels 543## 544setMethodS3( 545 ".checkTipLabels", 546 class="phylo", 547 function( 548 this, 549 ... 550 ){ 551 552 for(label in this$tip.label){ 553 if(length(grep("^Node \\d+$",label,perl=TRUE,value=FALSE)) > 0){ 554 throw("Sorry, but the node labels matching \"Node \\d+\" are reserved for internal nodes! Blaming label: ",label,".\n"); 555 } 556 else if(length(grep("^Root node \\d+$",label,perl=TRUE,value=FALSE)) > 0){ 557 throw("Sorry, but the node labels matching \"Root node \\d+\" are reserved for the root node! Blaming label: ",label,".\n"); 558 } 559 560 } 561 562 }, 563 private=TRUE, 564 protected=FALSE, 565 overwrite=FALSE, 566 conflict="warning", 567 validators=getOption("R.methodsS3:validators:setMethodS3") 568) 569 570## 571## Method: getPhylo 572## 573###########################################################################/** 574# 575# @RdocMethod getPhylo 576# 577# @title "Get the phylo object aggregated in a PhyloSim object" 578# 579# \description{ 580# @get "title". 581# 582# } 583# 584# @synopsis 585# 586# \arguments{ 587# \item{this}{A PhyloSim object.} 588# \item{...}{Not used.} 589# } 590# 591# \value{ 592# A phylo object or NA. 593# } 594# 595# \examples{ 596# #create a PhyloSim object 597# sim<-PhyloSim(); 598# # creat a phylo object 599# tree<-rcoal(3); 600# # get/set phylo object 601# setPhylo(sim,tree); 602# getPhylo(sim,tree); 603# # get/set phylo object via virtual field 604# sim$tree<-rcoal(5); 605# sim$tree; 606# } 607# 608# @author 609# 610# \seealso{ 611# The PhyloSim class, the \pkg{ape} package. 612# } 613# 614#*/########################################################################### 615setMethodS3( 616 "getPhylo", 617 class="PhyloSim", 618 function( 619 this, 620 ... 621 ){ 622 623 this$.phylo; 624 625 }, 626 private=FALSE, 627 protected=FALSE, 628 overwrite=FALSE, 629 conflict="warning", 630 validators=getOption("R.methodsS3:validators:setMethodS3") 631) 632 633## 634## Method: setRootSeq 635## 636###########################################################################/** 637# 638# @RdocMethod setRootSeq 639# 640# @title "Set the root sequence for a PhyloSim object" 641# 642# \description{ 643# @get "title". 644# 645# The root sequence will be used as a starting point for the simulation. The phylo object must be set before 646# trying to set the root sequence object. 647# } 648# 649# @synopsis 650# 651# \arguments{ 652# \item{this}{A PhyloSim object.} 653# \item{value}{A valid Sequence object.} 654# \item{...}{Not used.} 655# } 656# 657# \value{ 658# The root Sequence object if succesfull, FALSE otherwise. 659# } 660# 661# \examples{ 662# # create some objects 663# sim<-PhyloSim(phylo=rcoal(3)); 664# seq<-NucleotideSequence(string="ATGCC"); 665# # set/get root sequence 666# setRootSeq(sim, seq); 667# getRootSeq(sim, seq); 668# # set/get root sequence via virtual field 669# sim$rootSeq<-BinarySequence(string="111000111000"); 670# sim$rootSeq; 671# 672# } 673# 674# @author 675# 676# \seealso{ 677# @seeclass Sequence Process 678# } 679# 680#*/########################################################################### 681setMethodS3( 682 "setRootSeq", 683 class="PhyloSim", 684 function( 685 this, 686 value, 687 ... 688 ){ 689 690 if(missing(value)){ 691 throw("No object provided!\n"); 692 } 693 else if(!is.Sequence(value)){ 694 throw("The new value must be a sequence object!\n"); 695 } 696 else { 697 698 this$.root.sequence<-clone(value); 699 this$.root.sequence$name<-paste("Root node",this$rootNode); 700 701 # Call garbage collection: 702 gc(); 703 gc(); 704 705 return(this$.root.sequence); 706 707 } 708 return(FALSE); 709 710 711 }, 712 private=FALSE, 713 protected=FALSE, 714 overwrite=FALSE, 715 conflict="warning", 716 validators=getOption("R.methodsS3:validators:setMethodS3") 717) 718 719## 720## Method: getRootSeq 721## 722###########################################################################/** 723# 724# @RdocMethod getRootSeq 725# 726# @title "Get the root sequence aggregated by a PhyloSim object" 727# 728# \description{ 729# @get "title". 730# 731# } 732# 733# @synopsis 734# 735# \arguments{ 736# \item{this}{A PhyloSim object.} 737# \item{...}{Not used.} 738# } 739# 740# \value{ 741# The root Sequence object or NA. 742# } 743# 744# \examples{ 745# # create some objects 746# sim<-PhyloSim(phylo=rcoal(3)); 747# seq<-NucleotideSequence(string="ATGCC"); 748# # set/get root sequence 749# setRootSeq(sim, seq); 750# getRootSeq(sim, seq); 751# # set/get root sequence via virtual field 752# sim$rootSeq<-BinarySequence(string="111000111000"); 753# sim$rootSeq; 754# 755# } 756# 757# @author 758# 759# \seealso{ 760# @seeclass Sequence Process 761# } 762# 763#*/########################################################################### 764setMethodS3( 765 "getRootSeq", 766 class="PhyloSim", 767 function( 768 this, 769 ... 770 ){ 771 772 this$.root.sequence; 773 774 775 }, 776 private=FALSE, 777 protected=FALSE, 778 overwrite=FALSE, 779 conflict="warning", 780 validators=getOption("R.methodsS3:validators:setMethodS3") 781) 782 783## 784## Method: as.character.PhyloSim 785## 786###########################################################################/** 787# 788# @RdocMethod as.character 789# 790# @title "Return the character representation of a PhyloSim object" 791# 792# \description{ 793# @get "title". 794# 795# The character representation is the identifier of the PhyloSim object as returned by the \code{getId} method. 796# } 797# 798# @synopsis 799# 800# \arguments{ 801# \item{x}{A PhyloSim object.} 802# \item{...}{Not used.} 803# } 804# 805# \value{ 806# A character vector of length one. 807# } 808# 809# \examples{ 810# # create a PhyloSim object 811# o<-PhyloSim(name="MySim"); 812# # get character representation 813# as.character(o) 814# } 815# 816# @author 817# 818# \seealso{ 819# @seeclass 820# } 821# 822#*/########################################################################### 823setMethodS3( 824 "as.character", 825 class="PhyloSim", 826 function( 827 x, 828 ... 829 ){ 830 831 return(getId(x)); 832 833 }, 834 private=FALSE, 835 protected=FALSE, 836 overwrite=FALSE, 837 conflict="warning", 838 validators=getOption("R.methodsS3:validators:setMethodS3") 839) 840 841## 842## Method: getName 843## 844###########################################################################/** 845# 846# @RdocMethod getName 847# 848# @title "Get the name of a PhyloSim object" 849# 850# \description{ 851# @get "title". 852# } 853# 854# @synopsis 855# 856# \arguments{ 857# \item{this}{A PhyloSim object.} 858# \item{...}{Not used.} 859# } 860# 861# \value{ 862# A character vector of length one. 863# } 864# 865# \examples{ 866# # create a PhyloSim object 867# o<-PhyloSim(); 868# # set/get name 869# setName(o,"MySim"); 870# getName(o,"MySim"); 871# # set/get name via virtual field 872# o$name<-"George"; 873# o$name 874# } 875# 876# @author 877# 878# \seealso{ 879# @seeclass 880# } 881# 882#*/########################################################################### 883setMethodS3( 884 "getName", 885 class="PhyloSim", 886 function( 887 this, 888 ... 889 ){ 890 891 this$.name; 892 }, 893 private=FALSE, 894 protected=FALSE, 895 overwrite=FALSE, 896 conflict="warning", 897 validators=getOption("R.methodsS3:validators:setMethodS3") 898); 899 900## 901## Method: setName 902## 903###########################################################################/** 904# 905# @RdocMethod setName 906# 907# @title "Set the name of a PhyloSim object" 908# 909# \description{ 910# @get "title". 911# } 912# 913# @synopsis 914# 915# \arguments{ 916# \item{this}{A PhyloSim object.} 917# \item{new.name}{A character vector of length one.} 918# \item{...}{Not used.} 919# } 920# 921# \value{ 922# The new name. 923# } 924# 925# \examples{ 926# # create a PhyloSim object 927# o<-PhyloSim(); 928# # set/get name 929# setName(o,"MySim"); 930# getName(o,"MySim"); 931# # set/get name via virtual field 932# o$name<-"George"; 933# o$name 934# } 935# 936# @author 937# 938# \seealso{ 939# @seeclass 940# } 941# 942#*/########################################################################### 943setMethodS3( 944 "setName", 945 class="PhyloSim", 946 function( 947 this, 948 new.name, 949 ... 950 ){ 951 952 this$.name<-as.character(new.name); 953 }, 954 private=FALSE, 955 protected=FALSE, 956 overwrite=FALSE, 957 conflict="warning", 958 validators=getOption("R.methodsS3:validators:setMethodS3") 959); 960 961## 962## Method: getId 963## 964###########################################################################/** 965# 966# @RdocMethod getId 967# 968# @title "Get the unique identifier of a PhyloSim object" 969# 970# \description{ 971# @get "title". 972# The unique identifier is the concatenation of the class, the object name as returned by getName() and the object hash 973# as returned by hashCode(). 974# } 975# 976# @synopsis 977# 978# \arguments{ 979# \item{this}{A PhyloSim object.} 980# \item{...}{Not used.} 981# } 982# 983# \value{ 984# A character vector of length one. 985# } 986# 987# \examples{ 988# # create a PhyloSim object 989# o<-PhyloSim(name="MySim"); 990# # get id 991# getId(o); 992# # get id via virtual field 993# o$id; 994# } 995# 996# @author 997# 998# \seealso{ 999# @seeclass 1000# } 1001# 1002#*/########################################################################### 1003setMethodS3( 1004 "getId", 1005 class="PhyloSim", 1006 function( 1007 this, 1008 ... 1009 ){ 1010 1011 this.class<-class(this)[1]; 1012 id<-paste(this.class,this$.name,hashCode(this),sep=":"); 1013 return(id); 1014 1015 }, 1016 private=FALSE, 1017 protected=FALSE, 1018 overwrite=FALSE, 1019 conflict="warning", 1020 validators=getOption("R.methodsS3:validators:setMethodS3") 1021); 1022 1023## 1024## Method: setId 1025## 1026###########################################################################/** 1027# 1028# @RdocMethod setId 1029# 1030# @title "Forbidden action: setting the unique identifier of a PhyloSim object" 1031# 1032# \description{ 1033# @get "title". 1034# 1035# } 1036# 1037# @synopsis 1038# 1039# \arguments{ 1040# \item{this}{An object.} 1041# \item{value}{Not used.} 1042# \item{...}{Not used.} 1043# } 1044# 1045# \value{ 1046# Throws an error. 1047# } 1048# 1049# @author 1050# 1051# \seealso{ 1052# @seeclass 1053# } 1054# 1055#*/########################################################################### 1056setMethodS3( 1057 "setId", 1058 class="PhyloSim", 1059 function( 1060 this, 1061 value, 1062 ... 1063 ){ 1064 1065 throw("Id is generated automatically and it cannot be set!\n"); 1066 1067 }, 1068 private=FALSE, 1069 protected=FALSE, 1070 overwrite=FALSE, 1071 conflict="warning", 1072 validators=getOption("R.methodsS3:validators:setMethodS3") 1073); 1074 1075## 1076## Method: Simulate 1077## 1078###########################################################################/** 1079# 1080# @RdocMethod Simulate 1081# 1082# @title "Run a simulation according to a PhyloSim object" 1083# 1084# \description{ 1085# @get "title". 1086# 1087# The phylo object and the root sequence must be set before attempting to run a simulation. 1088# Also the bigRate of the root sequence must not be NA or zero, so at least one sane 1089# Process object must be attached to the root sequence object. 1090# } 1091# 1092# @synopsis 1093# 1094# \arguments{ 1095# \item{this}{A PhyloSim object.} 1096# \item{quiet}{TRUE or FALSE (default).} 1097# \item{...}{Not used.} 1098# } 1099# 1100# \value{ 1101# The PhyloSim object (invisible). 1102# } 1103# 1104# \examples{ 1105# # Create a PhyloSim object. 1106# # Provide the phylo object 1107# # and the root sequence. 1108# sim<-PhyloSim( 1109# name="TinySim", 1110# phylo=rcoal(3), 1111# root.seq=NucleotideSequence(string="ATGC",processes=list(list(JC69()))) 1112# ); 1113# # Run the simulation 1114# Simulate(sim); 1115# # Print the resulting sequences 1116# sim$sequences 1117# # Print the resulting alignment 1118# sim$alignment 1119# } 1120# 1121# @author 1122# 1123# \seealso{ 1124# @seeclass 1125# } 1126# 1127#*/########################################################################### 1128setMethodS3( 1129 "Simulate", 1130 class="PhyloSim", 1131 function( 1132 this, 1133 quiet=FALSE, 1134 ... 1135 ){ 1136 1137 1138 if(!is.phylo(this$.phylo)){ 1139 throw("Cannot simulate because the phylo object is not set or it is invalid!\n"); 1140 } 1141 # Check for the root sequence: 1142 else if(!is.Sequence(this$.root.sequence)){ 1143 throw("Cannot simulate because the root sequence is not set or it is invalid!\n"); 1144 } 1145 # Check bigRate validity: 1146 else if(is.na(this$.root.sequence$bigRate)){ 1147 throw("Cannot simulate because the bigRate of the root sequence is NA!\n"); 1148 } 1149 else{ 1150 1151 # Warn for zero bigRate: 1152 if(this$.root.sequence$bigRate == 0){ 1153 warning("The bigRate of the root sequence is zero! You are running a pointless simulation!\n"); 1154 } 1155 1156 if(exists(x="PSIM_FAST")){ 1157 if(!quiet){ 1158cat("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n"); 1159cat("!! WARNING: fast & careless mode is on, most of the error checking is omitted! !!\n"); 1160cat("!! Please note that this also disables the saving of branch statistics. !!\n"); 1161cat("!! You can go back to normal mode by deleting the PSIM_FAST object. !!\n"); 1162cat("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n"); 1163 } 1164 Log(this,"WARNING: fast & careless mode is on, most of the error checking is omitted!"); 1165 } 1166 1167 # Attach root sequence to root node: 1168 Log(this,paste("Attaching root sequence ",this$.root.sequence$id,sep="")); 1169 attachSeqToNode(this, node=getRootNode(this),seq=this$.root.sequence); 1170 1171 # Write protecting the root sequence: 1172 Log(this,paste("Write protecting root sequence ",this$.root.sequence$id,sep="")); 1173 this$.root.sequence$writeProtected<-TRUE; 1174 1175 # Traverse the tree and simulate: 1176 Log(this,paste("Starting simulation on the object",this$id)); 1177 edge.counter<-1; 1178 n.edges<-this$nedges; 1179 for(edge in 1:n.edges){ 1180 if(!quiet){ cat("Simulating edge",edge,"of", n.edges,"\n");} 1181 Log(this,paste("Starting to simulate edge",edge,"of",n.edges)); 1182 .simulateEdge(this,number=edge); 1183 edge.counter<-edge.counter+1; 1184 } 1185 } 1186 Log(this, "Simulation finished, building alignment!\n"); 1187 this$.alignment<-.recoverAlignment(this); 1188 # Flush the log connection: 1189 if(!is.na(this$.log.connection)){ 1190 close(this$.log.connection); 1191 } 1192 1193 # Call the garbage collector: 1194 gc(); 1195 gc(); 1196 1197 return(invisible(this)); 1198 }, 1199 private=FALSE, 1200 protected=FALSE, 1201 overwrite=FALSE, 1202 conflict="warning", 1203 validators=getOption("R.methodsS3:validators:setMethodS3") 1204); 1205 1206## 1207## Method: .simulateEdge 1208## 1209setMethodS3( 1210 ".simulateEdge", 1211 class="PhyloSim", 1212 function( 1213 this, 1214 number=NA, 1215 ... 1216 ){ 1217 1218 # Get edge: 1219 edge<-getEdge(this, number); 1220 # Get parent node: 1221 start.seq<-getSeqFromNode(this, edge[[1,"from"]]); 1222 # Evolve sequence: 1223 new.seq<-.evolveBranch(this, start.seq=start.seq, branch.length=edge[1,"length"], old.node=edge[[1,"from"]],new.node=edge[[1,"to"]], branch.number=number); 1224 # Write protect the sequence: 1225 new.seq$writeProtected<-TRUE; 1226 # Attach sequence to children node: 1227 attachSeqToNode(this, node=edge[1,"to"], seq=new.seq); 1228 1229 }, 1230 private=TRUE, 1231 protected=FALSE, 1232 overwrite=FALSE, 1233 conflict="warning", 1234 validators=getOption("R.methodsS3:validators:setMethodS3") 1235); 1236 1237## 1238## Method: .evolveBranch 1239## 1240setMethodS3( 1241 ".evolveBranch", 1242 class="PhyloSim", 1243 function( 1244 this, 1245 start.seq=NA, 1246 branch.length=NA, 1247 old.node=NA, 1248 new.node=NA, 1249 branch.number=NA, 1250 ... 1251 ){ 1252 1253 if(!exists(x="PSIM_FAST")){ 1254 1255 if(missing(start.seq)){ 1256 throw("No starting sequence provided!\n"); 1257 } 1258 else if(missing(branch.length)){ 1259 throw("No branch length provided!\n"); 1260 } 1261 else if(!is.numeric(branch.length)){ 1262 throw("The branch length must be numeric!\n"); 1263 } 1264 } 1265 if(.checkSeq(this, start.seq) ){ 1266 1267 # Cloning the starting sequence: 1268 seq<-clone(start.seq); 1269 1270 # Set the name of the sequence object: 1271 if(is.tip(this, new.node)){ 1272 seq$name<-this$tipLabels[[new.node]]; 1273 } 1274 else { 1275 seq$name<-paste("Node",new.node); 1276 } 1277 .GillespieDirect(this, seq=seq, branch.length=branch.length, branch.number=branch.number); 1278 1279 # Call the node hook if exists: 1280 hook<-this$.node.hooks[[as.character(new.node)]]; 1281 if(!is.null(hook) & is.function(hook)){ 1282 Log(this,paste("Calling node hook for node",new.node)); 1283 seq<-hook(seq=seq); 1284 if(!is.Sequence(seq)){ 1285 throw("Node hook returned an invalid sequence object!\n"); 1286 } 1287 else if(is.na(seq$bigRate)){ 1288 throw("Node hook returned sequence with NA bigRate!\n"); 1289 } 1290 else if(seq$bigRate == 0.0){ 1291 throw("Node hook returned sequence with zero bigRate!\n"); 1292 } 1293 else{ 1294 checkConsistency(seq, omit.sites=TRUE); 1295 } 1296 } 1297 1298 # Return the resulting sequence object: 1299 return(seq); 1300 1301 } 1302 }, 1303 private=FALSE, 1304 protected=FALSE, 1305 overwrite=FALSE, 1306 conflict="warning", 1307 validators=getOption("R.methodsS3:validators:setMethodS3") 1308); 1309 1310## 1311## Method: .GillespieDirect 1312## 1313setMethodS3( 1314 ".GillespieDirect", 1315 class="PhyloSim", 1316 function( 1317 this, 1318 seq=NA, 1319 branch.length=NA, 1320 branch.number=NA, 1321 ... 1322 ){ 1323 1324 Debug(this, paste("Branch length is",branch.length)); 1325 1326 # Initialize time: 1327 1328 time<-0.0; 1329 1330 # Sample the next waiting time until 1331 # the branch length is consumed: 1332 1333 while( (time<-time + rexp(n=1, rate=(big.rate<-getBigRate(seq)))) <= branch.length){ 1334 1335 # Generate a random number between zero and the bigRate: 1336 1337 E<-runif(n=1,min=0,max=big.rate); 1338 1339 # Identify the target site: 1340 1341 site.number<-which(.getCumulativeRatesFast(seq) >= E)[[1]]; 1342 1343 # Get the events from the target site: 1344 1345 site<-seq$.sites[[site.number]]; 1346 site$.position<-site.number; 1347 events<-getEvents(site); 1348 site$.position<-NULL; 1349 1350 # Get the rates: 1351 rates<-double(); 1352 for(e in events){ 1353 rates<-c(rates,e$.rate); 1354 } 1355 1356 # Calculate the corresponding cumulative rates: 1357 if(site.number > 1){ 1358 rates<-cumsum(c(seq$.cumulative.rates[[site.number - 1]], rates)); 1359 } 1360 else { 1361 rates<-cumsum(c(0.0, rates)); 1362 } 1363 1364 # Pick the event: 1365 1366 event.number<-which(rates >= E)[[1]] - 1; 1367 event<-events[[event.number]]; 1368 1369 # Log the event: 1370 1371 Log(this,paste("Performing event [",event$.name,"] at position",event$.position,"generated by the process",event$.process$.id)); 1372 1373 # Perform the event: 1374 1375 event.details<-Perform(event); 1376 Debug(this,paste("Remaining branch length is",(branch.length-time) )); 1377 1378 # Log event details: 1379 1380 # Log deletion event details: 1381 if(event$.name == "Deletion"){ 1382 Log(this,paste("The process",event$.process,"proposed to delete range",paste(event.details$range,collapse="--"),". Accepted:",event.details$accepted)); 1383 } 1384 # Log insertion event details: 1385 else if(event$.name == "Insertion"){ 1386 message<-paste("The process ",event$.process," proposed insertion at position ",event.details$position,". Accepted: ",event.details$accepted,sep=""); 1387 if(event.details$accepted == TRUE){ 1388 message<-paste(message,"."," Insert length was ",event.details$length,sep=""); 1389 } 1390 Log(this, message); 1391 } 1392 1393 # Update branch statistics if not in fast mode: 1394 if(!exists(x="PSIM_FAST")){ 1395 .UpdateBranchStats(this,event,event.details, branch.number); 1396 } 1397 1398 # Abort if sequence length shrunk to zero: 1399 1400 if(seq$.length == 0){ 1401 message<-paste("Terminating the simulation because the length of the sequence ",seq$name," shrunk to zero! Please be more careful when tuning the indel rates!\n"); 1402 Log(this, message); 1403 throw(message); 1404 } 1405 1406 } #/while 1407 1408 # Calling garbage collection: 1409 gc(); 1410 gc(); 1411 return(seq); 1412 1413 }, 1414 private=FALSE, 1415 protected=FALSE, 1416 overwrite=FALSE, 1417 conflict="warning", 1418 validators=getOption("R.methodsS3:validators:setMethodS3") 1419); 1420 1421## 1422## Method: attachSeqToNode 1423## 1424###########################################################################/** 1425# 1426# @RdocMethod attachSeqToNode 1427# 1428# @title "Assotiate a Sequence object with a given node of a phylo object aggregated by a PhyloSim object" 1429# 1430# \description{ 1431# @get "title". 1432# 1433# This method is mainly used internally. 1434# 1435# } 1436# 1437# @synopsis 1438# 1439# \arguments{ 1440# \item{this}{A PhyloSim object.} 1441# \item{node}{Node identifier.} 1442# \item{seq}{A Sequence object.} 1443# \item{...}{Not used.} 1444# } 1445# 1446# \value{ 1447# The PhyloSim object (invisible). 1448# } 1449# 1450# 1451# @author 1452# 1453# \seealso{ 1454# @seeclass 1455# } 1456# 1457#*/########################################################################### 1458setMethodS3( 1459 "attachSeqToNode", 1460 class="PhyloSim", 1461 function( 1462 this, 1463 node=NA, 1464 seq=NA, 1465 ... 1466 ){ 1467 1468 if(!is.phylo(this$.phylo)){ 1469 throw("The phylo object is not set, sequence to node is not possible!\n"); 1470 } 1471 if(missing(node)){ 1472 throw("No node specified!\n"); 1473 } 1474 else if(missing(seq)){ 1475 throw("No sequence object given"); 1476 } 1477 else if(.checkNode(this,node) & .checkSeq(this, seq)){ 1478 1479 if(is.Sequence(this$.sequences[[node]])){ 1480 throw("The node has already an attached sequence. Detach that before trying to attach a new one!\n"); 1481 } 1482 else { 1483 this$.sequences[[as.numeric(node)]]<-seq; 1484 return(invisible(this)); 1485 } 1486 1487 } 1488 1489 return(invisible(this)); 1490 }, 1491 private=FALSE, 1492 protected=FALSE, 1493 overwrite=FALSE, 1494 conflict="warning", 1495 validators=getOption("R.methodsS3:validators:setMethodS3") 1496); 1497 1498## 1499## Method: attachHookToNode 1500## 1501###########################################################################/** 1502# 1503# @RdocMethod attachHookToNode 1504# 1505# @title "Attach a callback function to a given node of a phylo object aggregated by a PhyloSim object" 1506# 1507# \description{ 1508# @get "title". 1509# 1510# A "node hook" is a function which accepts a Sequence object through the named argument "seq" and returns a 1511# Sequence object. The node hook function must accept any object which inherits from the \code{Sequence} class! 1512# 1513# After simulating the branch leading to the node, the resulting Sequence object is passed 1514# to the node hook and the returned object is used to simulate the downstream branches. 1515# 1516# By using node hooks the attached processes can be replaced during simulation, hence enabling the simulation of 1517# non-homogeneous sequence evolution. 1518# } 1519# 1520# @synopsis 1521# 1522# \arguments{ 1523# \item{this}{A PhyloSim object.} 1524# \item{node}{Node identifier.} 1525# \item{fun}{A function (see above).} 1526# \item{...}{Not used.} 1527# } 1528# 1529# \value{ 1530# The PhyloSim object (invisible). 1531# } 1532# 1533# \examples{ 1534# # Create a PhyloSim object. 1535# # Provide the phylo object 1536# # and the root sequence. 1537# sim<-PhyloSim( 1538# name="TinySim", 1539# phylo=rcoal(3), 1540# root.seq=NucleotideSequence(string="ATGC",processes=list(list(JC69()))) 1541# ); 1542# # create a node hook function 1543# hook<-function(seq=NA){ 1544# # replace the substitution process with F84 1545# if(inherits(seq,"NucleotideSequence")){ 1546# cat("Replacing JC69 with F84.\n"); 1547# seq$processes<-list(list(F84(rate.params=list("Kappa" = 2)))); 1548# } 1549# return(seq); 1550# } 1551# # attach hook function to node 5 1552# attachHookToNode(sim,5,hook); 1553# # Run the simulation 1554# Simulate(sim); 1555# # Check if the processes have been truly replaced 1556# lapply(sim$sequences, getUniqueProcesses.Sequence) 1557# # Print the resulting alignment 1558# sim$alignment 1559# } 1560# 1561# @author 1562# 1563# \seealso{ 1564# @seeclass 1565# } 1566# 1567#*/########################################################################### 1568setMethodS3( 1569 "attachHookToNode", 1570 class="PhyloSim", 1571 function( 1572 this, 1573 node=NA, 1574 fun=NA, 1575 ... 1576 ){ 1577 1578 if(!is.phylo(this$.phylo)){ 1579 throw("The phylo object is not set, attaching node hook is not possible!\n"); 1580 } 1581 if(missing(node)){ 1582 throw("No node specified!\n"); 1583 } 1584 else if(missing(fun)){ 1585 throw("No function given!"); 1586 } 1587 else if(!is.function(fun)){ 1588 throw("The argument \"fun\" must be a function!\n"); 1589 } 1590 else if( length(intersect(names(formals(fun)), "seq")) == 0 ){ 1591 throw("The function argument must have a an argument named \"seq\""); 1592 } 1593 else if(!is.Sequence(fun(Sequence(length=1)))){ 1594 throw("The insert hook function must return a Sequence object!\n"); 1595 } 1596 else if( .checkNode(this,node) ){ 1597 if(is.function(this$.node.hooks[[as.character(node)]])){ 1598 throw("The node has already an attached node hook. Detach that before trying to attach a new one!\n"); 1599 } 1600 else { 1601 this$.node.hooks[[as.character(node)]]<-fun; 1602 return(invisible(this)); 1603 } 1604 1605 } 1606 1607 return(invisible(this)); 1608 }, 1609 private=FALSE, 1610 protected=FALSE, 1611 overwrite=FALSE, 1612 conflict="warning", 1613 validators=getOption("R.methodsS3:validators:setMethodS3") 1614); 1615 1616 1617## 1618## Method: .checkNode 1619## 1620setMethodS3( 1621 ".checkNode", 1622 class="PhyloSim", 1623 function( 1624 this, 1625 node=NA, 1626 ... 1627 ){ 1628 1629 if(missing(node)){ 1630 throw("No node specified!\n"); 1631 } else if( length(intersect(node, getNodes(this))) != 1){ 1632 throw("The specified node is invalid!\n"); 1633 } 1634 else { 1635 return(TRUE); 1636 } 1637 }, 1638 private=FALSE, 1639 protected=FALSE, 1640 overwrite=FALSE, 1641 conflict="warning", 1642 validators=getOption("R.methodsS3:validators:setMethodS3") 1643); 1644 1645## 1646## Method: .checkSeq 1647## 1648setMethodS3( 1649 ".checkSeq", 1650 class="PhyloSim", 1651 function( 1652 this, 1653 seq=NA, 1654 ... 1655 ){ 1656 1657 if(missing(seq)){ 1658 throw("No sequence specified!\n"); 1659 } else if(!is.Sequence(seq)){ 1660 throw("The sequence object is invalid!\n"); 1661 } 1662 else { 1663 return(TRUE); 1664 } 1665 }, 1666 private=FALSE, 1667 protected=FALSE, 1668 overwrite=FALSE, 1669 conflict="warning", 1670 validators=getOption("R.methodsS3:validators:setMethodS3") 1671); 1672 1673## 1674## Method: detachSeqFromNode 1675## 1676###########################################################################/** 1677# 1678# @RdocMethod detachSeqFromNode 1679# 1680# @title "Detach a Sequence object from a given node of a phylo object aggregated by a PhyloSim object" 1681# 1682# \description{ 1683# @get "title". 1684# 1685# This method is mainly used internally. 1686# 1687# } 1688# 1689# @synopsis 1690# 1691# \arguments{ 1692# \item{this}{A PhyloSim object.} 1693# \item{node}{Node identifier.} 1694# \item{...}{Not used.} 1695# } 1696# 1697# \value{ 1698# The PhyloSim object (invisible). 1699# } 1700# 1701# @author 1702# 1703# \seealso{ 1704# @seeclass 1705# } 1706# 1707#*/########################################################################### 1708setMethodS3( 1709 "detachSeqFromNode", 1710 class="PhyloSim", 1711 function( 1712 this, 1713 node=NA, 1714 ... 1715 ){ 1716 1717 if(missing(node)){ 1718 throw("No node specified!\n"); 1719 } 1720 else if( .checkNode(this,node) ){ 1721 1722 this$.sequences[[as.numeric(node)]]<-NA; 1723 } 1724 1725 return(invisible(this)); 1726 }, 1727 private=FALSE, 1728 protected=FALSE, 1729 overwrite=FALSE, 1730 conflict="warning", 1731 validators=getOption("R.methodsS3:validators:setMethodS3") 1732); 1733 1734## 1735## Method: detachHookFromNode 1736## 1737###########################################################################/** 1738# 1739# @RdocMethod detachHookFromNode 1740# 1741# @title "Detach a node hook function from a given node of a phylo object aggregated by a PhyloSim object" 1742# 1743# \description{ 1744# @get "title". 1745# 1746# } 1747# 1748# @synopsis 1749# 1750# \arguments{ 1751# \item{this}{A PhyloSim object.} 1752# \item{node}{Node identifier.} 1753# \item{...}{Not used.} 1754# } 1755# 1756# \value{ 1757# The PhyloSim object (invisible). 1758# } 1759# 1760# \examples{ 1761# # Create a PhyloSim object. 1762# # Provide the phylo object 1763# # and the root sequence. 1764# sim<-PhyloSim( 1765# name="TinySim", 1766# phylo=rcoal(3), 1767# root.seq=NucleotideSequence(string="ATGC",processes=list(list(JC69()))) 1768# ); 1769# # create a node hook function 1770# hook<-function(seq=NA){ 1771# # replace the substitution process with F84 1772# if(inherits(seq,"NucleotideSequence")){ 1773# cat("Replacing JC69 with F84.\n"); 1774# seq$processes<-list(list(F84(rate.params=list("Kappa" = 2)))); 1775# } 1776# return(seq); 1777# } 1778# # attach hook function to node 5 1779# attachHookToNode(sim,5,hook); 1780# # detach hook from node 5 1781# detachHookFromNode(sim,5); 1782# # Run the simulation again 1783# Simulate(sim); # You should not see the message printed out by the "hook" function. 1784# 1785# } 1786# 1787# @author 1788# 1789# \seealso{ 1790# attachHookToNode PhyloSim Simulate.PhyloSim 1791# } 1792# 1793#*/########################################################################### 1794setMethodS3( 1795 "detachHookFromNode", 1796 class="PhyloSim", 1797 function( 1798 this, 1799 node=NA, 1800 ... 1801 ){ 1802 1803 if(missing(node)){ 1804 throw("No node specified!\n"); 1805 } 1806 else if( .checkNode(this,node) ){ 1807 1808 this$.node.hooks[[as.character(node)]]<-NA; 1809 } 1810 1811 return(invisible(this)); 1812 }, 1813 private=FALSE, 1814 protected=FALSE, 1815 overwrite=FALSE, 1816 conflict="warning", 1817 validators=getOption("R.methodsS3:validators:setMethodS3") 1818); 1819 1820 1821## 1822## Method: getSeqFromNode 1823## 1824###########################################################################/** 1825# 1826# @RdocMethod getSeqFromNode 1827# 1828# @title "Get the Sequence object associated with a given node of a phylo object aggregated by a PhyloSim object" 1829# 1830# \description{ 1831# @get "title". 1832# } 1833# 1834# @synopsis 1835# 1836# \arguments{ 1837# \item{this}{A PhyloSim object.} 1838# \item{node}{Node identifier.} 1839# \item{...}{Not used.} 1840# } 1841# 1842# \value{ 1843# A Sequence object. 1844# } 1845# 1846# \examples{ 1847# # Create a PhyloSim object. 1848# # Provide the phylo object 1849# # and the root sequence. 1850# sim<-PhyloSim( 1851# name="TinySim", 1852# phylo=rcoal(3), 1853# root.seq=NucleotideSequence(string="ATG",processes=list(list(JC69()))) 1854# ); 1855# # get the sequence associated with node 5 1856# getSeqFromNode(sim,5) # Should be NA 1857# # Run the simulation 1858# Simulate(sim) 1859# # try again 1860# getSeqFromNode(sim,5) 1861# } 1862# 1863# @author 1864# 1865# \seealso{ 1866# @seeclass 1867# } 1868# 1869#*/########################################################################### 1870setMethodS3( 1871 "getSeqFromNode", 1872 class="PhyloSim", 1873 function( 1874 this, 1875 node=NA, 1876 ... 1877 ){ 1878 1879 if(missing(node)){ 1880 throw("No node specified!\n"); 1881 } 1882 else if( .checkNode(this,node) ){ 1883 1884 return(this$.sequences[[as.numeric(node)]]); 1885 } 1886 1887 }, 1888 private=FALSE, 1889 protected=FALSE, 1890 overwrite=FALSE, 1891 conflict="warning", 1892 validators=getOption("R.methodsS3:validators:setMethodS3") 1893); 1894 1895## 1896## Method: getSequences 1897## 1898###########################################################################/** 1899# 1900# @RdocMethod getSequences 1901# 1902# @title "Gets all the Sequence objects associated with the nodes of a phylo object aggregated by a PhyloSim object" 1903# 1904# \description{ 1905# @get "title". 1906# 1907# The order of the Sequence objects in the returned list reflects the identifiers of the associated nodes. 1908# } 1909# 1910# @synopsis 1911# 1912# \arguments{ 1913# \item{this}{A PhyloSim object.} 1914# \item{...}{Not used.} 1915# } 1916# 1917# \value{ 1918# A list of sequence objects. 1919# } 1920# 1921# \examples{ 1922# # Create a PhyloSim object. 1923# # Provide the phylo object 1924# # and the root sequence. 1925# sim<-PhyloSim( 1926# name="TinySim", 1927# phylo=rcoal(3), 1928# root.seq=NucleotideSequence(string="ATG",processes=list(list(JC69()))) 1929# ); 1930# # run the simulation 1931# Simulate(sim) 1932# # get all the associated sequence objects 1933# getSequences(sim) 1934# # get the sequence associated with node 3 1935# # via virtual field 1936# sim$sequences[[3]] 1937# } 1938# 1939# @author 1940# 1941# \seealso{ 1942# @seeclass 1943# } 1944# 1945#*/########################################################################### 1946setMethodS3( 1947 "getSequences", 1948 class="PhyloSim", 1949 function( 1950 this, 1951 ... 1952 ){ 1953 1954 slist<-list(); 1955 for (node in getNodes(this)){ 1956 slist[[node]]<-getSeqFromNode(this, node=node); 1957 } 1958 return(slist); 1959 1960 }, 1961 private=FALSE, 1962 protected=FALSE, 1963 overwrite=FALSE, 1964 conflict="warning", 1965 validators=getOption("R.methodsS3:validators:setMethodS3") 1966); 1967 1968## 1969## Method: setSequences 1970## 1971###########################################################################/** 1972# 1973# @RdocMethod setSequences 1974# 1975# @title "Forbidden action: setting the Sequence objects associated with the nodes of a phylo object aggregated by a PhyloSim object" 1976# 1977# \description{ 1978# @get "title". 1979# } 1980# 1981# @synopsis 1982# 1983# \arguments{ 1984# \item{this}{An object.} 1985# \item{value}{Not used.} 1986# \item{...}{Not used.} 1987# } 1988# 1989# \value{ 1990# Throws an error. 1991# } 1992# 1993# @author 1994# 1995# \seealso{ 1996# @seeclass 1997# } 1998# 1999#*/########################################################################### 2000setMethodS3( 2001 "setSequences", 2002 class="PhyloSim", 2003 function( 2004 this, 2005 value, 2006 ... 2007 ){ 2008 2009 virtualAssignmentForbidden(this); 2010 2011 }, 2012 private=FALSE, 2013 protected=FALSE, 2014 overwrite=FALSE, 2015 conflict="warning", 2016 validators=getOption("R.methodsS3:validators:setMethodS3") 2017); 2018 2019## 2020## Method: getAlignment 2021## 2022###########################################################################/** 2023# 2024# @RdocMethod getAlignment 2025# 2026# @title "Get the alignment stored in a PhyloSim object" 2027# 2028# \description{ 2029# @get "title". 2030# } 2031# 2032# @synopsis 2033# 2034# \arguments{ 2035# \item{this}{A PhyloSim object.} 2036# \item{...}{Not used.} 2037# } 2038# 2039# \value{ 2040# The alignment as a matrix. Gap are represented by strings composed of dashes. 2041# } 2042# 2043# \examples{ 2044# # Create a PhyloSim object. 2045# # Provide the phylo object 2046# # and the root sequence. 2047# sim<-PhyloSim( 2048# name="TinySim", 2049# phylo=rcoal(3), 2050# root.seq=NucleotideSequence(string="ATG",processes=list(list(JC69()))) 2051# ); 2052# # run the simulation 2053# Simulate(sim) 2054# # get the resulting aligment 2055# getAlignment(sim) 2056# # via virtual fileld: 2057# sim$alignment 2058# } 2059# 2060# @author 2061# 2062# \seealso{ 2063# @seeclass 2064# } 2065# 2066#*/########################################################################### 2067setMethodS3( 2068 "getAlignment", 2069 class="PhyloSim", 2070 function( 2071 this, 2072 ... 2073 ){ 2074 2075 this$.alignment; 2076 2077 }, 2078 private=FALSE, 2079 protected=FALSE, 2080 overwrite=FALSE, 2081 conflict="warning", 2082 validators=getOption("R.methodsS3:validators:setMethodS3") 2083); 2084 2085## 2086## Method: setAlignment 2087## 2088###########################################################################/** 2089# 2090# @RdocMethod setAlignment 2091# 2092# @title "Forbidden action: setting the alignment stored in a PhyloSim object" 2093# 2094# \description{ 2095# @get "title". 2096# } 2097# 2098# @synopsis 2099# 2100# \arguments{ 2101# \item{this}{An object.} 2102# \item{value}{Not used.} 2103# \item{...}{Not used.} 2104# } 2105# 2106# \value{ 2107# Throws an error. 2108# } 2109# 2110# @author 2111# 2112# \seealso{ 2113# @seeclass 2114# } 2115# 2116#*/########################################################################### 2117setMethodS3( 2118 "setAlignment", 2119 class="PhyloSim", 2120 function( 2121 this, 2122 value, 2123 ... 2124 ){ 2125 2126 this$alignment <- value; 2127 2128 }, 2129 private=FALSE, 2130 protected=FALSE, 2131 overwrite=FALSE, 2132 conflict="warning", 2133 validators=getOption("R.methodsS3:validators:setMethodS3") 2134); 2135 2136## 2137## Method: .recoverAlignment 2138## 2139setMethodS3( 2140 ".recoverAlignment", 2141 class="PhyloSim", 2142 function( 2143 this, 2144 paranoid=FALSE, 2145 ... 2146 ){ 2147 2148 # Refuse to build alignment if at least one of the sequences is NA: 2149 for (seq in this$.sequences){ 2150 if(!is.Sequence(seq)){ 2151 throw("Cannot build alignment because the simulation is incomplete!\n"); 2152 } 2153 } 2154 2155 # The list holding all the partial alignment matrices: 2156 aln.mat<-list(); 2157 2158 # Assigning NA-s here to prevent creation of these variables in the global 2159 # environment. 2160 row.names<-NA; 2161 from.node<-NA; 2162 to.node<-NA; 2163 from.seq<-NA; 2164 to.seq<-NA; 2165 edge<-NA; 2166 from.name<-NA; 2167 to.name<-NA; 2168 from.mat<-NA; 2169 to.mat<-NA; 2170 2171 # Initialize the variables: 2172 2173 init.vars<-function(){ 2174 2175 # Getting the edge: 2176 edge<<-getEdge(this, edge.number); 2177 2178 # Getting the nodes: 2179 from.node<<-edge[[1,"from"]]; 2180 to.node<<-edge[[1,"to"]]; 2181 2182 # Getting the sequence objects: 2183 from.seq<<-getSeqFromNode(this, from.node) 2184 to.seq<<-getSeqFromNode(this, to.node) 2185 2186 # Getting sequence names: 2187 from.name<<-from.seq$name; 2188 to.name<<-to.seq$name; 2189 2190 } 2191 2192 # Initialize the aligment matrices: 2193 init.aln.mats…
Large files files are truncated, but you can click here to view the full file