/BrownianInsertor.R
R | 660 lines | 230 code | 62 blank | 368 comment | 33 complexity | 652347048cfdbbb71074e00555530cbf MD5 | raw file
1 2## Copyright 2009 Botond Sipos 3## See the package description for licensing information. 4## 5##########################################################################/** 6# 7# @RdocClass BrownianInsertor 8# 9# @title "The BrownianInsertor class" 10# 11# \description{ 12# The \code{BrownianInsertor} class inherits from the \code{DiscreteInsertor} 13# or \code{ContinuousInsertor} class depending on the \code{type} constructor argument 14# ("discrete" or "continuous"). 15# 16# This process generates the insert sequence based on the sites flanking the insertions as follows: 17# \itemize{ 18# \item An insert length is sampled by calling the function stored in the \code{proposeBy} virtual field. 19# \item A sequence object is constructed. 20# \item The processes attached to both flanking sites are attached to the insert sequence. If there are no common processes, the processes from a randomly chosen site will be attached to the insert. 21# \item The site-process specific parameters are sampled from Brownian paths with linear trends having the values from the flanking sites as endpoints. 22# } 23# 24# The "noisiness" of the Brownian path can be controlled through the \code{scale} virtual field/constructor parameter. 25# 26# @classhierarchy 27# } 28# 29# @synopsis 30# 31# \arguments{ 32# \item{name}{Object name.} 33# \item{type}{Process type (see above).} 34# \item{scale}{Brownian path scale parameter.} 35# \item{...}{Additional arguments.} 36# } 37# 38# \section{Fields and Methods}{ 39# @allmethods 40# } 41# 42# \examples{ 43# # create a BrownianInsertor process, discrete type 44# p<-BrownianInsertor( 45# type="discrete", 46# scale=0.05, 47# sizes=1:4, 48# probs=c(3/6,1/6,1/6,1/6), 49# rate=0.05 50# ) 51# # get object summary 52# summary(p) 53# # plot insert length distribution 54# plot(p) 55# # create a nucleotide sequence, attach processes 56# s<-NucleotideSequence(string="AAAAAAAAAAA",processes=list(list(p,JC69()))) 57# # create simulation object 58# sim<-PhyloSim(root.seq=s, phylo=rcoal(2)) 59# # simulate and show alignment 60# Simulate(sim) 61# sim$alignment 62# # check the rate multipliers and insertion tolerances in one of the sequences 63# res<-sim$sequences[[2]] 64# getRateMultipliers(res,p) 65# getInsertionTolerance(res,p) 66# } 67# 68# @author 69# 70# \seealso{ 71# DiscreteInsertor ContinuousInsertor GeneralInsertor GeneralInDel 72# } 73# 74#*/########################################################################### 75setConstructorS3( 76 "BrownianInsertor", 77 function( 78 name="Anonymous", 79 type="discrete", 80 scale=0.001, 81 ... 82 ) { 83 84 if(type == "continuous"){ 85 this<-ContinuousInsertor( 86 ... 87 ); 88 } 89 else if (type == "discrete") { 90 this<-DiscreteInsertor( 91 ... 92 ); 93 } 94 else { 95 throw("Invalid insertor process type!\n"); 96 } 97 98 this<-extend( 99 this, 100 "BrownianInsertor", 101 .scale = NA, 102 .type = type 103 ); 104 105 # Using virtual field to clear Id cache: 106 this$name<-name; 107 # setting scale 108 this$scale<-scale; 109 110 this$generateBy<-function(process=NA,length=NA,target.seq=NA,event.pos=NA,insert.pos=NA){ 111 112 if(is.na(target.seq)){ 113 return(NA); 114 } 115 116 if(is.na(length) | (length(length) == 0) | length == 0){ 117 throw("Invalid insert length!\n"); 118 } 119 120 # The start and end of the Brownian path: 121 122 start; 123 end; 124 proc<-list(); 125 126 if( (event.pos == 1) || (event.pos == target.seq$.length) ){ 127 128 start<-clone(target.seq$.sites[[event.pos]]); 129 start$.state=NA; 130 end<-clone(start); 131 proc<-getProcesses(start); 132 133 } else { 134 start<-clone(target.seq$.sites[[insert.pos]]); 135 start$.state=NA; 136 end<-clone(target.seq$.sites[[insert.pos + 1]]); 137 end$.state=NA; 138 139 proc.start<-getProcesses(start); 140 proc.end<-getProcesses(end); 141 142 # Calculate the intersection of process list: 143 144 proc<-PSRoot$intersect.list(proc.start,proc.end); 145 146 # No common processes: 147 148 if(length(proc) == 0){ 149 coin.flip<-sample(c(0,1),1); 150 if (coin.flip) { 151 proc <- getProcesses(start) 152 end <- clone(start) 153 } else { 154 proc <- getProcesses(end) 155 start <- clone(end) 156 } 157 } 158 } 159 160 161 # Create the insert sequence: 162 163 class.seq<-class(target.seq)[[1]]; 164 insert<-do.call(class.seq,list(length=length)); 165 setProcesses(this=insert,value=list(proc)); 166 167 # For every process... 168 169 for (p in proc){ 170 171 # ... and site specific parameter: 172 for(param in getSiteSpecificParamIds(p)){ 173 174 start.value<-getParameterAtSite(p,site=start,id=param)$value; 175 end.value<-getParameterAtSite(p,site=end,id=param)$value; 176 177 path<-seq(from=start.value,to=end.value,length.out=(insert$.length + 2)); 178 path<-path[2:(length(path)-1)]; 179 brownian.path<-abs(BrownianInsertor$BrownianPath(p=path, a=this$.scale)); 180 181 # Tolerance values are probabilities, do not alow them to wander beyond 1: 182 if(param == "insertion.tolerance" || param == "deletion.tolerance"){ 183 brownian.path[which(brownian.path > 1)]<-1; 184 } 185 186 setParameterAtSites( 187 insert, 188 process = p, 189 id = param, 190 value = brownian.path 191 ); 192 193 } 194 } 195 196 return(insert); 197 198 199 } 200 201 return(this); 202 }, 203 enforceRCC=TRUE 204); 205 206## 207## Method: checkConsistency 208## 209###########################################################################/** 210# 211# @RdocMethod checkConsistency 212# 213# @title "Check object consistency" 214# 215# \description{ 216# @get "title". 217# } 218# 219# @synopsis 220# 221# \arguments{ 222# \item{this}{An object.} 223# \item{...}{Not used.} 224# } 225# 226# 227# \value{ 228# Returns an invisible TRUE if no inconsistencies found in the object, throws 229# an error otherwise. 230# } 231# 232# @author 233# 234# \seealso{ 235# @seeclass 236# } 237# 238#*/########################################################################### 239setMethodS3( 240 "checkConsistency", 241 class="BrownianInsertor", 242 function( 243 this, 244 ... 245 ){ 246 247 wp<-this$writeProtected; 248 if (wp) { 249 this$writeProtected<-FALSE; 250 } 251 252 may.fail<-function(this) { 253 254 this$scale<-this$scale; 255 if( (this$.type != "discrete") && (this$.type != "continuous") ){ 256 throw("BrownianInsertor type is invalid!\n"); 257 } 258 259 } 260 tryCatch(may.fail(this),finally=this$writeProtected<-wp); 261 NextMethod(); 262 263 }, 264 private=FALSE, 265 protected=FALSE, 266 overwrite=FALSE, 267 conflict="warning", 268 validators=getOption("R.methodsS3:validators:setMethodS3") 269); 270 271## 272## Method: BrownianPath 273## 274###########################################################################/** 275# 276# @RdocMethod BrownianPath 277# 278# @title "Generate a Brownian path" 279# 280# \description{ 281# @get "title". 282# 283# This method generates a Brownian path given the scale parameter a (determining "noisiness") 284# and the vector p describing the trends. More useful as a static method. 285# } 286# 287# @synopsis 288# 289# \arguments{ 290# \item{this}{A BrownianInsertor object.} 291# \item{p}{Path parameter (a numeric vector).} 292# \item{a}{Scale paramater (a numeric vector of length one).} 293# \item{...}{Not used.} 294# } 295# 296# \value{ 297# A numeric vector. 298# } 299# 300# \examples{ 301# path<-BrownianInsertor$BrownianPath(a=2,p=1:10); 302# path 303# plot(path) 304# } 305# 306# @author 307# 308# \seealso{ 309# @seeclass 310# } 311# 312#*/########################################################################### 313setMethodS3( 314 "BrownianPath", 315 class="BrownianInsertor", 316 function( 317 this, 318 p=NA, 319 a=NA, 320 ... 321 ){ 322 323 generate_brownian<-function(length, a){ 324 cumsum(rnorm(length,0,sd=a)); 325 } 326 327 generate_bridge <- function (length,a){ 328 b <- generate_brownian(length,a) 329 b - (1:length)/(length) * b[length] 330 } 331 332 generate_path <- function (p,a){ 333 n <- length(p); 334 b <- generate_bridge (n+1,a); 335 p + b[1:n]; 336 } 337 338 return(generate_path(p,a)); 339 340 341 }, 342 private=FALSE, 343 protected=FALSE, 344 overwrite=FALSE, 345 conflict="warning" 346); 347 348 349## 350## Method: getType 351## 352###########################################################################/** 353# 354# @RdocMethod getType 355# 356# @title "Get the type of the BrownianInsertor process" 357# 358# \description{ 359# @get "title". 360# 361# If type is \code{discrete}, than the process inherits from \code{DiscreteInsertor}. 362# If the type is \code{continuous}, then the object inherits from \code{ContinuousInsertor}. 363# } 364# 365# @synopsis 366# 367# \arguments{ 368# \item{this}{A BrownianInsertor object.} 369# \item{...}{Not used.} 370# } 371# 372# \value{ 373# A character vector of length one. 374# } 375# 376# \examples{ 377# p<-BrownianInsertor(type="discrete") 378# # get type 379# getType(p) 380# # get upstream classes 381# class(p) 382# p<-BrownianInsertor(type="continuous") 383# # get type via virtual field 384# p$type 385# # get upstream classes 386# class(p) 387# } 388# 389# @author 390# 391# \seealso{ 392# @seeclass 393# } 394# 395#*/########################################################################### 396setMethodS3( 397 "getType", 398 class="BrownianInsertor", 399 function( 400 this, 401 ... 402 ){ 403 404 this$.type; 405 406 }, 407 private=FALSE, 408 protected=FALSE, 409 overwrite=FALSE, 410 conflict="warning", 411 validators=getOption("R.methodsS3:validators:setMethodS3") 412); 413 414## 415## Method: setType 416## 417###########################################################################/** 418# 419# @RdocMethod setType 420# 421# @title "Forbidden action: setting the type of a BrownianInsertor object" 422# 423# \description{ 424# @get "title". 425# 426# The type can be set only from the \code{type} constructor argument and cannot be changed later. 427# } 428# 429# @synopsis 430# 431# \arguments{ 432# \item{this}{An object.} 433# \item{value}{Not used.} 434# \item{...}{Not used.} 435# } 436# 437# \value{ 438# Throws an error. 439# } 440# 441# @author 442# 443# \seealso{ 444# @seeclass 445# } 446# 447#*/########################################################################### 448setMethodS3( 449 "setType", 450 class="BrownianInsertor", 451 function( 452 this, 453 value, 454 ... 455 ){ 456 457 throw("The type of the BrownianInsertor objects can be set only from the constructor!\n"); 458 459 }, 460 private=FALSE, 461 protected=FALSE, 462 overwrite=FALSE, 463 conflict="warning", 464 validators=getOption("R.methodsS3:validators:setMethodS3") 465); 466 467## 468## Method: getScale 469## 470###########################################################################/** 471# 472# @RdocMethod getScale 473# 474# @title "Get scale parameter" 475# 476# \description{ 477# @get "title". 478# } 479# 480# @synopsis 481# 482# \arguments{ 483# \item{this}{A BrownianInsertor object.} 484# \item{...}{Not used.} 485# } 486# 487# \value{ 488# A numeric vector of length one. 489# } 490# 491# \examples{ 492# # create a BrownianInsertor object 493# p<-BrownianInsertor(scale=0.002) 494# # set/get scale parameter 495# setScale(p,0.1) 496# getScale(p) 497# # set/get scale parameter via virtual field 498# p$scale<-0.1 499# p$scale 500# } 501# 502# @author 503# 504# \seealso{ 505# @seeclass 506# } 507# 508#*/########################################################################### 509setMethodS3( 510 "getScale", 511 class="BrownianInsertor", 512 function( 513 this, 514 ... 515 ){ 516 517 this$.scale; 518 519 }, 520 private=FALSE, 521 protected=FALSE, 522 overwrite=FALSE, 523 conflict="warning", 524 validators=getOption("R.methodsS3:validators:setMethodS3") 525); 526 527## 528## Method: setScale 529## 530###########################################################################/** 531# 532# @RdocMethod setScale 533# 534# @title "Set scale parameter" 535# 536# \description{ 537# @get "title". 538# } 539# 540# @synopsis 541# 542# \arguments{ 543# \item{this}{A BrownianInsertor object.} 544# \item{value}{A numeric vector of length one.} 545# \item{...}{Not used.} 546# } 547# 548# \value{ 549# value (invisible). 550# } 551# 552# \examples{ 553# # create a BrownianInsertor object 554# p<-BrownianInsertor(scale=0.002) 555# # set/get scale parameter 556# setScale(p,0.1) 557# getScale(p) 558# # set/get scale parameter via virtual field 559# p$scale<-0.1 560# p$scale 561# } 562# 563# @author 564# 565# \seealso{ 566# @seeclass 567# } 568# 569#*/########################################################################### 570setMethodS3( 571 "setScale", 572 class="BrownianInsertor", 573 function( 574 this, 575 value, 576 ... 577 ){ 578 579 580 .checkWriteProtection(this); 581 if(!is.numeric(value) || (length(value) != 1)){ 582 throw("The value of the scale paramter must be a numeric vector of length 1!\n"); 583 } 584 this$.scale<-value; 585 return(invisible(this$.scale)); 586 }, 587 private=FALSE, 588 protected=FALSE, 589 overwrite=FALSE, 590 conflict="warning", 591 validators=getOption("R.methodsS3:validators:setMethodS3") 592); 593 594## 595## Method: summary 596## 597###########################################################################/** 598# 599# @RdocMethod summary 600# 601# @title "Summarize the properties of an object" 602# 603# \description{ 604# @get "title". 605# } 606# 607# @synopsis 608# 609# \arguments{ 610# \item{object}{An object} 611# \item{...}{Not used.} 612# } 613# 614# \value{ 615# Returns a PSRootSummary object. 616# } 617# 618# \examples{ 619# 620# # create an object 621# p<-BrownianInsertor( 622# type="discrete", 623# scale=0.05, 624# sizes=1:4, 625# probs=c(3/6,1/6,1/6,1/6), 626# rate=0.05 627# ) 628# # get a summary 629# summary(p) 630# } 631# 632# @author 633# 634# \seealso{ 635# @seeclass 636# } 637# 638#*/########################################################################### 639setMethodS3( 640 "summary", 641 class="BrownianInsertor", 642 function( 643 object, 644 ... 645 ){ 646 this<-object; 647 .addSummaryNameId(this); 648 this$.summary$"Type"<-this$.type; 649 this$.summary$"Brownian path scale parameter"<-this$.scale; 650 651 NextMethod(); 652 653 }, 654 private=FALSE, 655 protected=FALSE, 656 overwrite=FALSE, 657 conflict="warning", 658 validators=getOption("R.methodsS3:validators:setMethodS3") 659); 660