R/skylineplot.R

   1 ## skylineplot.R (2004-07-4)
   2
   3 ##   Various methods to plot skyline objects (= skyline plots)
   4
   5 ## Copyright 2002-2004 Korbinian Strimmer
   6
   7 ## This file is part of the R-package `ape'.
   8 ## See the file ../COPYING for licensing issues.
   9
  10 # plot skyline
  11 plot.skyline <- function(x, show.years=FALSE, subst.rate, present.year, ...)
  12 {
  13   if (class(x) != "skyline")
  14     stop("object \"x\" is not of class \"skyline\"")
  15   t <- x$time
  16   m <- x$population.size
  17   lm <- length(m)
  18
  19   if (show.years)
  20   {
  21     plot((-c(0,t))/subst.rate+present.year,c(m,m[lm]),type="s",
  22      xlab="time (years)",ylab="effective population size",log="y", ...)
  23
  24   }
  25   else
  26   {
  27     plot(c(0,t),c(m,m[lm]),type="s", xlim=c(t[lm],0),
  28      xlab="time (past to present in units of substitutions)",ylab="effective population size",log="y", ...)
  29   }
  30
  31 }
  32
  33 # plot another skyline plot on top
  34 lines.skyline <- function(x, show.years=FALSE, subst.rate, present.year, ...)
  35 {
  36   if (class(x) != "skyline")
  37     stop("object \"x\" is not of class \"skyline\"")
  38   t <- x$time
  39   m <- x$population.size
  40   lm <- length(m)
  41
  42
  43   if (show.years)
  44   {
  45     lines((-c(0,t))/subst.rate+present.year,c(m,m[lm]),type="s", ...)
  46   }
  47   else
  48   {
  49     lines(c(0,t),c(m,m[lm]),type="s", ...)
  50   }
  51 }
  52
  53
  54 # convenience short cut (almost compatible with APE 0.1)
  55 skylineplot <- function(z, ...) plot(skyline(z, ...))
  56
  57
  58 #input: phylogenetic tree
  59 skylineplot.deluxe <- function(tree, ...)
  60 {
  61   if (class(tree) != "phylo")
  62     stop("object \"tree\" is not of class \"phylo\"")
  63
  64   ci <- coalescent.intervals(tree)
  65   classic <- skyline(ci)
  66   generalized <- skyline(ci, -1)
  67   plot(classic,col=grey(.8), ...)
  68   lines(generalized, ...)
  69   return(generalized)
  70 }
  71
  72
  73