Last active
August 26, 2020 13:48
-
-
Save josephwb/f3d35f8833a07f71002af7726b12652b to your computer and use it in GitHub Desktop.
Check the monophyletic status of every genus in a tree
This file contains 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
require(phangorn); | |
# check the monophyletic status of every genus in the tree `phy` | |
check_genus_monophyly <- function (phy) { | |
# get genera. assumes form: Genus_species_whatever | |
gens <- sort(unique(sub("_.*", "", phy$tip.label))); | |
n <- length(gens); | |
cat("Checking ", n, " genera for monophyly status.\n", sep=""); | |
# the number of descendant tips of the MRCA. preallocate for monotypic genera | |
ndec <- rep(1, n); | |
# get number of matches per genus. don't match shorter names to longer ones | |
yy <- match(gsub("_.*", "_", phy$tip.label), paste0(gens, "_")); | |
# the expected number of descendants == tips initially matched | |
nexp <- as.numeric(table(yy)); | |
# get descendants for all nodes (phangorn) | |
decs <- phangorn:::bip(phy); # almost instantaneous | |
# skip monotypic genera | |
idx <- which(nexp != 1); | |
if (length(idx) > 0) { # only enter if at least one genus has > 1 representative | |
for (i in 1:length(idx)) { | |
mrca <- getMRCA(phy, which(yy == idx[i])); | |
ndec[idx[i]] <- length(decs[[mrca]]); | |
} | |
} | |
cat("Found ", sum(ndec != nexp), " non-monophyletic genera.\n", sep=""); | |
return(data.frame(Genus=gens, Mono=(ndec == nexp), N.obs=ndec, N.exp=nexp, stringsAsFactors=FALSE)); | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment