Skip to content

Instantly share code, notes, and snippets.

@josephwb
Last active August 26, 2020 13:48
Show Gist options
  • Save josephwb/f3d35f8833a07f71002af7726b12652b to your computer and use it in GitHub Desktop.
Save josephwb/f3d35f8833a07f71002af7726b12652b to your computer and use it in GitHub Desktop.
Check the monophyletic status of every genus in a tree
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