find.lineages: Separate a paleobuddy simulation into monophyletic clades

Description Usage Arguments Value Author(s) Examples

View source: R/find.lineages.R

Description

Separates a sim object into sim objects each with a mother species and its descendants. Returns by default the list of sim objects descended from each species with an NA parent in the original input. Allows for the user to input a vector of species to be the mother of each resulting member of the returning list instead. Returns for each clade a vector with the original identity of member species as well.

Usage

1
find.lineages(sim, S = NULL)

Arguments

sim

A simulation containing lists for the speciation times, extinction times, parent identities and status (extant or extinct).

S

A vector of species in sim. If not supplied, S will be the starting species in the simulation (i.e. those for which the parent is NA).

Value

A list object with (named) sim objects corresponding to the clades descended from species in S. For each clade, an extra vector sim$LIN is included so the user can identify the order of species in the return with the order of species in the original simulation.

Author(s)

Bruno Petrucci and Matheus Januario.

Examples

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
# we will start with examples where S are the starting species

###
# first, let us try a simulation with 3 clades,
sim <- bd.sim(n0 = 3, pp = 0.1, qq = 0.1, tMax = 10, nFinal = c(20, Inf))

# using the functions
clades <- find.lineages(sim)

# set up for plotting side by syde
par(mfrow = c(1,length(clades)))

# for each clade
for (i in 1:length(clades)) {
  # change NA to 0 on the clade's TE
  clades[[i]]$TE[clades[[i]]$EXTANT] = 0
  
  # if there is only one lineage in the clade
  if (length(clades[[i]]$TE) < 2) {
    # placeholder plot
    plot(NA, xlim = c(-1, 1), ylim = c(-1, 1))
    text("simulation with \n just one lineage", x = 0, y = 0.5, cex = 2)
  }
  # if it is a proper phylogeny
  else {
    if (requireNamespace("ape", quietly = TRUE)) {
      p <- ape::plot.phylo(
        make.phylo(clades[[i]]),
        main = "red: extinction events \n blue: speciation events");
      ape::axisPhylo()
    }
    
    # checking speciation times:
    for (j in 2:length(clades[[i]]$TS)) {
      # the subtraction is just to adjust the wt with the plot scale
      lines(x = c(
        sort(clades[[i]]$TS, decreasing = TRUE)[2] - clades[[i]]$TS[j],
        sort(clades[[i]]$TS, decreasing = TRUE)[2] - clades[[i]]$TS[j]),
        y = c(p$y.lim[1], p$y.lim[2]), lwd = 2, col = "blue")
    }
    
    # checking extinction times:
    for (j in 1:length(sim$TE)) {
      # the subtraction is just to adjust the wt with the plot scale
      lines(x = c(
        sort(clades[[i]]$TS, decreasing = TRUE)[2] - clades[[i]]$TE[j],
        sort(clades[[i]]$TS, decreasing = TRUE)[2] - clades[[i]]$TE[j]),
        y = c(p$y.lim[1], p$y.lim[2]), lwd = 2, col = "red")
    }
  }
}

###
# it works with any number of clades, of course
sim <- bd.sim(n0 = 5, pp = 0.1, qq = 0.08, tMax = 10, nFinal = c(20, Inf))

# using the functions
clades <- find.lineages(sim)

# set up for plotting side by syde
par(mfrow = c(1,length(clades)))

# for each clade
for (i in 1:length(clades)) {
  # change NA to 0 on the clade's TE
  clades[[i]]$TE[clades[[i]]$EXTANT] = 0
  
  # if there is only one lineage in the clade
  if (length(clades[[i]]$TE) < 2) {
    # placeholder plot
    plot(NA, xlim = c(-1, 1), ylim = c(-1, 1))
    text("simulation with \n just one lineage", x = 0, y = 0.5, cex = 2)
  }
  # if it is a proper phylogeny
  else {
    if (requireNamespace("ape", quietly = TRUE)) {
      p <- ape::plot.phylo(
        make.phylo(clades[[i]]),
        main = "red: extinction events \n blue: speciation events");
      ape::axisPhylo()
    }
    
    # checking speciation times:
    for (j in 2:length(clades[[i]]$TS)) {
      # the subtraction is just to adjust the wt with the plot scale
      lines(x = c(
        sort(clades[[i]]$TS, decreasing = TRUE)[2] - clades[[i]]$TS[j],
        sort(clades[[i]]$TS, decreasing = TRUE)[2] - clades[[i]]$TS[j]),
        y = c(p$y.lim[1], p$y.lim[2]), lwd = 2, col = "blue")
    }
    
    # checking extinction times:
    for (j in 1:length(sim$TE)) {
      # the subtraction is just to adjust the wt with the plot scale
      lines(x = c(
        sort(clades[[i]]$TS, decreasing = TRUE)[2] - clades[[i]]$TE[j],
        sort(clades[[i]]$TS, decreasing = TRUE)[2] - clades[[i]]$TE[j]),
        y = c(p$y.lim[1], p$y.lim[2]), lwd = 2, col = "red")
    }
  }
}

###
# including one clade
sim <- bd.sim(n0 = 1, pp = 0.1, qq = 0.08, tMax = 10, nFinal = c(5, Inf))

par(mfrow = c(1, 2))

# plotting sim and find.lineages(sim) - should be equal
if (requireNamespace("ape", quietly = TRUE)) {
  ape::plot.phylo(make.phylo(sim), main="original")
  ape::axisPhylo()
  ape::plot.phylo(make.phylo(find.lineages(sim)[[1]]), main="after find.lineages()")
  ape::axisPhylo()
}

###
# now let us check that when S does not contain a starting species, we still
# get correct subsets of the simulation
sim <- bd.sim(1, 0.1, 0.05, 40, nFinal = c(5, Inf))

# making sure we have a couple of clades to explore
while ((length(which(sim$PAR == 1)) < 3) | (length(which(sim$PAR == 2)) < 3) |
       (length(which(sim$PAR == 3)) < 3)) {
  sim <- bd.sim(1, 0.2, 0.1, 10)
}

if (requireNamespace("ape", quietly = TRUE)) {
  # first we plot the clade started by 1
  ape::plot.phylo(make.phylo(sim), main="original")
  
  # this should look the same
  ape::plot.phylo(make.phylo(find.lineages(sim)[[1]]), main="after find.lineages()")
  
  # and these should be part of the previous phylogenies
  ape::plot.phylo(make.phylo(find.lineages(sim, c(2, 3))$clade_2),
                  main = "Clade_2")
  ape::plot.phylo(make.phylo(find.lineages(sim, c(2, 3))$clade_3),
                  main = "Clade_3")
}

brpetrucci/paleobuddy documentation built on Aug. 8, 2020, 2:03 a.m.