PageRenderTime 49ms CodeModel.GetById 33ms app.highlight 12ms RepoModel.GetById 0ms app.codeStats 1ms

/BrownianInsertor.R

http://github.com/sbotond/phylosim
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