PageRenderTime 19ms CodeModel.GetById 8ms app.highlight 7ms RepoModel.GetById 1ms app.codeStats 0ms

/DiscreteDeletor.R

http://github.com/sbotond/phylosim
R | 642 lines | 235 code | 33 blank | 374 comment | 35 complexity | d6994f94aea6d6b6ab04641e545b1e48 MD5 | raw file
  1##	
  2## Copyright 2009 Botond Sipos	
  3## See the package description for licensing information.	
  4##	
  5##########################################################################/** 
  6#
  7# @RdocClass DiscreteDeletor
  8# 
  9# @title "The DiscreteDeletor class"
 10# 
 11# \description{ 
 12#	This class implements a process which performs deletions with
 13#	lengths sampled from a user-specified discrete distribution.
 14#	See \code{GeneralDeletor} for how the deletion process works.
 15#
 16#	@classhierarchy
 17# }
 18#	
 19# @synopsis
 20#	
 21# \arguments{
 22# 	\item{name}{The name of the object.}
 23# 	\item{rate}{The general rate.}
 24#	\item{sizes}{The deletion sizes to propose.}
 25#	\item{probs}{A vector with the probabilites of the deletion sizes.}
 26# 	\item{...}{Additional arguments.}
 27#	}
 28# 
 29# \section{Fields and Methods}{ 
 30# 	@allmethods
 31# }
 32# 
 33# \examples{ 
 34#	# create a DiscreteDeletor process
 35#	d<-DiscreteDeletor(
 36#		name="M.D.",
 37#		rate=0.25,
 38#		sizes=c(1,2),
 39#		probs=c(1/2,1/2)
 40# 	)
 41#	# get object summary
 42#	summary(d)
 43#	# set/get deletions sizes
 44#	d$sizes<-1:3
 45#	d$sizes
 46#	# set/get length probabilities
 47#	d$probs<-c(3,2,1)/6
 48#	d$probs
 49#	# plot length distribution
 50#	plot(d)
 51#	
 52#	# The following code illustrates how to use
 53#	# a DiscreteDeletor process in a simulation
 54#	
 55#	# create a sequence object, attach process d
 56#	s<-NucleotideSequence(string="AAAAAAAAAAGGGGAAAAAAAAAA",processes=list(list(d)))
 57#	# set the deletion tolerance to zero in the range 11:15
 58#	# creating a region rejecting all deletions
 59#	setDeletionTolerance(s,d,0,11:15)	
 60#	# get deletion tolerances
 61#	getDeletionTolerance(s,d)
 62#	# create a simulation object
 63#	sim<-PhyloSim(root.seq=s,phylo=rcoal(2))
 64#	# simulate
 65#	Simulate(sim)
 66#	# print resulting alignment
 67#	sim$alignment
 68# }
 69# 
 70# @author
 71#
 72# \seealso{ 
 73# 	GeneralDeletor ContinuousDeletor GeneralInDel
 74# }
 75# 
 76#*/###########################################################################
 77setConstructorS3(
 78  "DiscreteDeletor",
 79  function( 
 80		name="Anonymous",
 81		rate=NA,
 82		sizes=NA,
 83		probs=NA,
 84		... 
 85		)	{
 86
 87		this<-GeneralDeletor(
 88			 name=NA,
 89			 rate=rate,
 90			 propose.by=NA,
 91    	 accept.by=NA,
 92			...
 93		);
 94    this<-extend(
 95      this,
 96      "DiscreteDeletor",
 97			.sizes=NA,
 98			.probs=NA
 99    );
100		# Using virtual field to clear Id cache:
101		this$name<-name;
102
103		STATIC<-TRUE;
104		if(!missing(sizes)) {
105			this$sizes<-sizes;
106			STATIC<-FALSE;
107		}
108		
109		if(!missing(probs)) {
110			this$probs<-probs;
111			STATIC<-FALSE;
112		}
113
114		this$proposeBy<-function(process=NA,...){
115			if(!exists(x="PSIM_FAST")){
116			if( !is.numeric(process$.sizes) | !is.numeric(process$.probs) ){
117				throw("Cannot propose deletion length because the length distribution is not defined properly!\n");
118			}
119			}
120			if(length(process$.sizes) == 1){
121				return(process$.sizes[[1]]);
122			} else {
123				return(sample(x=process$.sizes,size=1,replace=FALSE,prob=process$.probs));
124			}
125		}
126
127    return(this);
128  },
129  enforceRCC=TRUE
130);
131
132##
133## Method: checkConsistency
134##
135###########################################################################/**
136#
137# @RdocMethod	checkConsistency
138# 
139# @title "Check object consistency"
140# 
141# \description{ 
142#		@get "title".
143# } 
144# 
145# @synopsis 
146#
147# \arguments{ 
148#       \item{this}{An object.} 
149#       \item{...}{Not used.} 
150# } 
151# 
152# 
153# \value{ 
154#		Returns an invisible TRUE if no inconsistencies found in the object, throws 
155#		an error otherwise. 
156# } 
157# 
158# @author 
159# 
160# \seealso{ 
161# 	@seeclass 
162# } 
163# 
164#*/###########################################################################
165setMethodS3(
166  "checkConsistency",
167  class="DiscreteDeletor",
168  function(
169    this,
170    ...
171  ){
172
173          wp<-this$writeProtected;
174      if (wp) {
175        this$writeProtected<-FALSE;
176      }
177
178      may.fail<-function(this) {
179				
180        if (is.numeric(this$sizes)) {
181          this$sizes<-this$sizes;
182        } 
183				else if (!is.na(this$sizes)){
184					throw("Deletion size vector is invalid!\n");
185				}
186
187        if (is.numeric(this$probs)) {
188          this$probs<-this$probs;
189        }
190				else if (!is.na(this$probs)){
191					throw("Deletion size probability vector is invalid!\n");
192				}
193
194      }
195      tryCatch(may.fail(this),finally=this$writeProtected<-wp);
196      NextMethod();
197
198  },
199  private=FALSE,
200  protected=FALSE,
201  overwrite=FALSE,
202  conflict="warning",
203  validators=getOption("R.methodsS3:validators:setMethodS3")
204);
205
206
207##	
208## Method: getSizes
209##	
210###########################################################################/**
211#
212# @RdocMethod getSizes
213# 
214# @title "Get the sizes of the proposed deletions" 
215# 
216# \description{ 
217#	@get "title".
218# } 
219# 
220# @synopsis 
221# 
222# \arguments{ 
223# 	\item{this}{A DiscreteDeletor object.} 
224# 	\item{...}{Not used.} 
225# } 
226# 
227# \value{ 
228# 	A vector of integers.
229# } 
230# 
231# \examples{
232#	# create a DiscreteDeletor object
233#	d<-DiscreteDeletor(rate=1)
234#	# set deletion sizes      
235#	setSizes(d,c(1,2,3))
236#	# get deletion sizes
237#	getSizes(d)
238#	# set/get sizes via virtual field
239#	d$sizes<-1:10
240#	d$sizes	
241# } 
242# 
243# @author 
244# 
245# \seealso{ 
246# 	@seeclass 
247# } 
248# 
249#*/###########################################################################
250setMethodS3(
251	"getSizes", 
252	class="DiscreteDeletor", 
253	function(
254		this,
255		...
256	){
257
258		this$.sizes;
259		
260	},
261	private=FALSE,
262	protected=FALSE,
263	overwrite=FALSE,
264	conflict="warning",
265	validators=getOption("R.methodsS3:validators:setMethodS3")
266);
267
268##	
269## Method: setSizes
270##	
271###########################################################################/**
272#
273# @RdocMethod setSizes
274# 
275# @title "Set the sizes of the proposed deletions" 
276# 
277# \description{ 
278#	@get "title".
279#	
280#	The provided numeric vector is rounded.
281# } 
282# 
283# @synopsis 
284# 
285# \arguments{ 
286# 	\item{this}{A DiscreteDeletor object.} 
287# 	\item{value}{A numeric vector.} 
288# 	\item{...}{Not used.} 
289# } 
290# 
291# \value{ 
292# 	A vector of integers (invisible).
293# } 
294# 
295# \examples{
296#	# create a DiscreteDeletor object
297#	d<-DiscreteDeletor(rate=1)
298#	# set deletion sizes      
299#	setSizes(d,c(1,2,3))
300#	# get deletion sizes
301#	getSizes(d)
302#	# set/get sizes via virtual field
303#	d$sizes<-1:10
304#	d$sizes	
305# } 
306# 
307# @author 
308# 
309# \seealso{ 
310# 	@seeclass 
311# } 
312# 
313#*/###########################################################################
314setMethodS3(
315	"setSizes", 
316	class="DiscreteDeletor", 
317	function(
318		this,
319		value,
320		...
321	){
322	
323		.checkWriteProtection(this);
324		if (missing(value)) {
325			throw("No new value provided!\n");
326		} else if (!is.numeric(value)) {
327			throw("The new value must be numeric vector!\n");
328		} else {
329			if(length(value) == 0 ) {
330				warning("Deletion size vector has zero length!\n");
331			}
332			this$.sizes<-round(value);	
333		}
334		return(invisible(this$.sizes));
335		
336	},
337	private=FALSE,
338	protected=FALSE,
339	overwrite=FALSE,
340	conflict="warning",
341	validators=getOption("R.methodsS3:validators:setMethodS3")
342);
343
344##	
345## Method: getProbs
346##	
347###########################################################################/**
348#
349# @RdocMethod getProbs
350# 
351# @title "Get the deletion length probabilities" 
352# 
353# \description{ 
354#	@get "title".
355# } 
356# 
357# @synopsis 
358# 
359# \arguments{ 
360# 	\item{this}{A DiscreteDeletor object.} 
361# 	\item{...}{Not used.} 
362# } 
363# 
364# \value{ 
365# 	A numeric vector with the deletion length probabilities.
366# } 
367# 
368# \examples{
369#	# create a DiscreteDeletor object
370#	d<-DiscreteDeletor(rate=1, sizes=1:3)
371#	# set/get length probabilities
372#	setProbs(d,c(1/3,1/3,1/3)) # equal probabilites
373#	getProbs(d)
374#	# set/get length probabilities via virtual field
375#	x<-c(2,2,1)
376#	# normalize x
377#	x<-x/sum(x)
378#	d$probs<-x
379#	d$probs
380# } 
381# 
382# @author 
383# 
384# \seealso{ 
385# 	@seeclass 
386# } 
387# 
388#*/###########################################################################
389setMethodS3(
390	"getProbs", 
391	class="DiscreteDeletor", 
392	function(
393		this,
394		...
395	){
396
397		this$.probs;
398		
399	},
400	private=FALSE,
401	protected=FALSE,
402	overwrite=FALSE,
403	conflict="warning",
404	validators=getOption("R.methodsS3:validators:setMethodS3")
405);
406
407##	
408## Method: setProbs
409##	
410###########################################################################/**
411#
412# @RdocMethod setProbs
413# 
414# @title "Set the deletion length probabilities" 
415# 
416# \description{ 
417#	@get "title".
418#
419#	The \code{sizes} virtual field must be set before setting the length probabilities.
420#	The length of the provided numeric vector must match with the length of the vector
421#	stored in the \code{sizes} virtual field. The vector is rescaled if the values do not
422#	sum to one and a warning is issued.
423# } 
424# 
425# @synopsis 
426# 
427# \arguments{ 
428# 	\item{this}{A DiscreteDeletor object.} 
429#	\item{value}{A numeric vector containg the length probabilities.}
430# 	\item{...}{Not used.} 
431# } 
432# 
433# \value{ 
434# 	The vector of probabilities.
435# } 
436# 
437# \examples{
438#	# create a DiscreteDeletor object
439#	d<-DiscreteDeletor(rate=1, sizes=1:3)
440#	# set/get length probabilities
441#	setProbs(d,c(1/3,1/3,1/3)) # equal probabilites
442#	getProbs(d)
443#	# set/get length probabilities via virtual field
444#	x<-c(2,2,1)
445#	# normalize x
446#	x<-x/sum(x)
447#	d$probs<-x
448#	d$probs
449# } 
450# 
451# @author 
452# 
453# \seealso{ 
454# 	@seeclass 
455# } 
456# 
457#*/###########################################################################
458setMethodS3(
459	"setProbs", 
460	class="DiscreteDeletor", 
461	function(
462		this,
463		value,
464		...
465	){
466	
467		.checkWriteProtection(this);
468		if (missing(value)) {
469			throw("No new value provided!\n");
470		} 
471		else if(!is.numeric(this$.sizes)) {
472			throw("Cannot set probabilities because indel size vector is not defined!\n");
473		}
474		else if (!is.numeric(value)) {
475			throw("The new value must be a numeric vector!\n");
476		}
477		else if(length(value) != length(this$.sizes)) {
478			throw("The length of the probabilities vector must be the same as the length of the deletion sizes vector");	
479		} 
480		else if( length(value[value < 0 ]) != 0 ) {
481			throw("The elements of the probability vector must not be negative!\n");
482		}
483		else {
484  		if(!isTRUE(all.equal(sum(value),1.0))){
485        		value<-value/sum(value);
486        		warning("The provided values were rescaled in order to sum to one!\n");
487    		}
488			this$.probs<-value;
489		}
490		
491	},
492	private=FALSE,
493	protected=FALSE,
494	overwrite=FALSE,
495	conflict="warning",
496	validators=getOption("R.methodsS3:validators:setMethodS3")
497);
498
499##	
500## Method: summary
501##	
502###########################################################################/**
503#
504# @RdocMethod summary
505#
506# @title "Summarize the properties of an object"
507#
508# \description{
509#       @get "title".
510# }
511#
512# @synopsis
513#
514# \arguments{
515#       \item{object}{An object}
516#       \item{...}{Not used.}
517# }
518#
519# \value{
520#  Returns a PSRootSummary object.
521# }
522#
523# \examples{
524#       # create an object
525#       a<-DiscreteDeletor(rate=1,sizes=c(1,2),probs=c(1/2,1/2))
526#       # get a summary
527#       summary(a)
528# }
529#
530# @author
531#
532# \seealso{
533#       @seeclass
534# }
535#
536#*/###########################################################################
537setMethodS3(
538	"summary", 
539	class="DiscreteDeletor", 
540	function(
541		object,
542		...
543	){
544
545		this<-object;
546		.addSummaryNameId(this);
547
548		expected.length<-NA;	
549		sd<-NA;
550		if( is.numeric(this$sizes) & is.numeric(this$probs)) {
551		expected.length<-weighted.mean(this$sizes, this$probs);		
552		sd<-sqrt(sum( (this$sizes - expected.length)^2 * this$probs ));
553		}
554
555		this$.summary$"Expected deletion length"<-expected.length;
556		this$.summary$"Standard deviation of deletion length"<-format(sd);
557		NextMethod();
558		
559	},
560	private=FALSE,
561	protected=FALSE,
562	overwrite=FALSE,
563	conflict="warning",
564	validators=getOption("R.methodsS3:validators:setMethodS3")
565);
566
567##	
568## Method: plot
569##	
570###########################################################################/**
571#
572# @RdocMethod plot
573# 
574# @title "Plot the deletion length distribution" 
575# 
576# \description{ 
577#	@get "title".
578# } 
579# 
580# @synopsis 
581# 
582# \arguments{ 
583# 	\item{x}{A DiscreteDeletor object.} 
584# 	\item{...}{Not used.} 
585# } 
586# 
587# \value{ 
588# 	The DiscreteDeletor object (invisible).
589# } 
590# 
591# \examples{
592#	d<-DiscreteDeletor(
593#		name="MyDiscDel",
594#		sizes=1:6,
595#		probs=c(0.25000000, 0.16666667, 0.08333333, 0.08333333, 0.16666667, 0.25000000)
596#	)
597#	# plot the deletion length distribution
598#	plot(d)
599# } 
600# 
601# @author 
602# 
603# \seealso{ 
604# 	@seeclass 
605# } 
606# 
607#*/###########################################################################
608setMethodS3(
609	"plot", 
610	class="DiscreteDeletor", 
611	function(
612		x,
613		...
614	){
615			this<-x;
616			if( !is.numeric(this$sizes) | !is.numeric(this$probs) ){
617				warning("Deletion length distribution is not defined properly! Nothing to plot here!\n");
618				return();
619			}
620			plot.default(	
621				x=this$sizes,
622				y=this$probs,
623				col=c("blue"),
624				lwd=2,
625				type="h",
626				main=paste("Deletion size distribution for:",this$id), 
627				xlab="Size",
628				ylab="Probability",
629				ylim=c(0,1),
630				xaxt="n"
631			);
632			axis(side=1, at=this$sizes, labels=this$sizes);
633			return(invisible(this));
634		
635	},
636	private=FALSE,
637	protected=FALSE,
638	overwrite=FALSE,
639	conflict="warning",
640	validators=getOption("R.methodsS3:validators:setMethodS3")
641);
642