Skip to content

Instantly share code, notes, and snippets.

@roblanf
Created September 12, 2022 02:04
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save roblanf/8a24cfb711ee6cf64b7c7861849d89b7 to your computer and use it in GitHub Desktop.
Save roblanf/8a24cfb711ee6cf64b7c7861849d89b7 to your computer and use it in GitHub Desktop.
What is happening with is.monophyletic()
``` 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)
```
![](https://i.imgur.com/2EPJr4c.png)
``` 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")
```
![](https://i.imgur.com/ITCo2Rh.png)
``` 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>
@J-Moravec
Copy link

Can't reproduce.

utils::sessionInfo()
R version 4.2.1 (2022-06-23)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04.1 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] ape_5.6-2

loaded via a namespace (and not attached):
[1] compiler_4.2.1  parallel_4.2.1  Rcpp_1.0.9      nlme_3.1-159   
[5] grid_4.2.1      lattice_0.20-45

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment