Created
September 12, 2022 02:04
-
-
Save roblanf/8a24cfb711ee6cf64b7c7861849d89b7 to your computer and use it in GitHub Desktop.
What is happening with is.monophyletic()
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
``` r | |
library(ape) | |
# make 1K trees to ensure we get some that seem odd... | |
trees = lapply(rep(4, 1000), rtree, rooted = FALSE) | |
# test monophyly of all six possible pairs | |
m1 = unlist(lapply(trees, is.monophyletic, tips = c("t1", "t2"))) | |
m2 = unlist(lapply(trees, is.monophyletic, tips = c("t1", "t3"))) | |
m3 = unlist(lapply(trees, is.monophyletic, tips = c("t1", "t4"))) | |
m4 = unlist(lapply(trees, is.monophyletic, tips = c("t2", "t3"))) | |
m5 = unlist(lapply(trees, is.monophyletic, tips = c("t2", "t4"))) | |
m6 = unlist(lapply(trees, is.monophyletic, tips = c("t3", "t4"))) | |
# each tree should have exactly two monophyletic pairs | |
totals = m1+m2+m3+m4+m5+m6 | |
# a data frame for digging deeper if you like | |
d = data.frame(t1_t2 = m1, t1_t3 = m2, t1_t4 = m3, t2_t3 = m4, t2_t4 = m5, t3_t4 = m6, total = totals) | |
# here's the histogram. I expected everything to get a score of two! | |
hist(totals) | |
``` | |
 | |
``` r | |
# and we can pick out an individual example like this | |
# here I'm choosing a random example which apparently has *no* monophyletic pairs... | |
tree_i = sample(which(totals==0), 1) | |
t = trees[[tree_i]] | |
plot(t, type="unrooted") | |
``` | |
 | |
``` r | |
t | |
#> | |
#> Phylogenetic tree with 4 tips and 2 internal nodes. | |
#> | |
#> Tip labels: | |
#> [1] "t4" "t2" "t3" "t1" | |
#> | |
#> Unrooted; includes branch lengths. | |
is.monophyletic(t, tips = c("t1", "t2")) | |
#> [1] FALSE | |
is.monophyletic(t, tips = c("t1", "t3")) | |
#> [1] FALSE | |
is.monophyletic(t, tips = c("t1", "t4")) | |
#> [1] FALSE | |
is.monophyletic(t, tips = c("t2", "t3")) | |
#> [1] FALSE | |
is.monophyletic(t, tips = c("t2", "t4")) | |
#> [1] FALSE | |
is.monophyletic(t, tips = c("t3", "t4")) | |
#> [1] FALSE | |
``` | |
<sup>Created on 2022-09-12 by the [reprex package](https://reprex.tidyverse.org) (v0.3.0)</sup> |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Can't reproduce.