From 8827983a913c1bc87cb200487b36d2be52e142fb Mon Sep 17 00:00:00 2001 From: sdgamboa Date: Tue, 5 Dec 2023 22:13:40 -0500 Subject: [PATCH] add interna_nodes.Rmd vignette/artticle --- vignettes/articles/internal_nodes.Rmd | 154 ++++++++++++++++++++++++++ 1 file changed, 154 insertions(+) create mode 100644 vignettes/articles/internal_nodes.Rmd diff --git a/vignettes/articles/internal_nodes.Rmd b/vignettes/articles/internal_nodes.Rmd new file mode 100644 index 0000000..2d5f276 --- /dev/null +++ b/vignettes/articles/internal_nodes.Rmd @@ -0,0 +1,154 @@ +--- +title: "internal_nodes" +--- + +```{r, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>" +) +``` + +```{r setup, message=FALSE} +library(phytools) +library(dplyr) +library(tidyr) +library(tibble) +``` + +## Data + +```{r} +data('primate.tree') +data('primate.data') +tre <- primate.tree +dat <- primate.data +states_original <- dat |> + rownames_to_column(var = 'tip_label') |> + select(tip_label, Activity_pattern) |> + mutate(state = 1) |> + pivot_wider( + names_from = 'Activity_pattern', values_from = 'state', values_fill = 0 + ) |> + column_to_rownames(var = 'tip_label') |> + as.matrix() |> + magrittr::set_colnames(LETTERS[1:3]) +states_original <- states_original[tre$tip.label,] + + +macaca <- grep('^Macaca_', rownames(states_original)) # ? +galago <- grep('^Galago_', rownames(states_original)) # B +eulemur <- grep('^Eulemur', rownames(states_original)) # C + +states <- states_original +states[] <- 1 / ncol(states) +# m1[macaca, ] <- c(rep(1,3), rep(0, 6)) +states[galago, ] <- c(rep(0,4), rep(1, 4), rep(0, 4)) +states[eulemur, ] <- c(rep(0,8), rep(1, 4)) +myMod <- 'ARD' +myPi <- 'fitzjohn' +``` + + +## Annotate one genus with an additional tip of zero length + +This is the recommended way of adding an annotations to an internal node. + +```{r} +macaca <- grep('Macaca', rownames(states)) +macaca_mrc <- findMRCA(tree = tre, tips = rownames(states)[macaca]) +t1 <- bind.tip(tre, tip.label = 'Macaca_gn1', edge.length = 0, where = macaca_mrc) +m1 <- states_original +# m1[] <- 1 / 3 +m1 <- rbind(m1, matrix(c(A = 0, B = 0, C = 1), nrow = 1, byrow = TRUE)) +rownames(m1)[nrow(m1)] <- 'Macaca_gn1' +m1[macaca,] <- rep(1 / ncol(m1), length(macaca) * 3) + +fit1 <- fitMk( + tree = t1, x = m1, model = myMod, pi = myPi, lik.func = "pruning", + logscale = TRUE +) +ace1 <- ancr(fit1, tips=TRUE) +plot(ace1, args.plotTree = list(direction = "upwards")) +``` + +However, in some cases, one might need to add two or more tips carrying +conclicting annotations. If both are of length zero, this throws an error: + +```{r} +macaca_mrc2 <- findMRCA(t1, tips = rownames(states)[macaca]) +t2 <- bind.tip(t1, tip.label = 'Macaca_gn2', edge.length = 0, where = macaca_mrc2) +m2 <- states +m2 <- rbind(m2, matrix(c(A = 0, B = 0, C = 1), nrow = 1, byrow = TRUE)) +m2 <- rbind(m2, matrix(c(A = 0, B = 1, C = 0), nrow = 1, byrow = TRUE)) +rownames(m2)[nrow(m2) - 1] <- 'Macaca_gn1' +rownames(m2)[nrow(m2)] <- 'Macaca_gn2' +m2[grep('^Macaca_[^g][^n][^\\d].*', rownames(m2)),] <- rep(1 / ncol(m2), length(macaca) * 3) +fit2 <- tryCatch( + error = function(e) e, + { + fitMk( + tree = t2, x = m2, model = myMod, pi = myPi, + lik.func = "pruning", logscale = TRUE + ) + } +) +print(fit2$message) +``` + +## Annotate one genus with an additional tip of length 1e-06 + +A solultion to the issue described above could be adding small quantity to the +branch length of the additional tips instead of a hard zero, e.g., 1e-06. At +least with the example below, the output seems to be the same (compare with +the first very firt case in this document). + +```{r} +macaca <- grep('Macaca', rownames(states)) +macaca_mrc <- findMRCA(tree = tre, tips = rownames(states)[macaca]) +t1 <- bind.tip(tre, tip.label = 'Macaca_gn1', edge.length = 1e-06, where = macaca_mrc) +m1 <- states_original +# m1[] <- 1 / 3 +m1 <- rbind(m1, matrix(c(A = 0, B = 0, C = 1), nrow = 1, byrow = TRUE)) +rownames(m1)[nrow(m1)] <- 'Macaca_gn1' +m1[macaca,] <- rep(1 / ncol(m1), length(macaca) * 3) + +fit1 <- fitMk( + tree = t1, x = m1, model = myMod, pi = myPi, lik.func = "pruning", + logscale = TRUE +) +ace1 <- ancr(fit1, tips=TRUE) +plot(ace1, args.plotTree = list(direction = "upwards")) +``` + +## Addint two tips with branch length of 1e-06 + +Now let's see if this has the expected result if I have two annotations. +For this, I'll make all annotations uncertain, except for the tips of the genus +Eulermur (C/red) and Galago (B/white). I'll add two conflicting annotations +to the genus Macaca using two tips with branch length of 1e-06. +This seems to have given the expected result: half and half. + +```{r} +t3 <- bind.tip(tre, tip.label = 'Macaca_gn1', edge.length = 1e-06, where = macaca_mrc) +macaca_mrc3 <- findMRCA(t3, tips = rownames(states)[macaca]) +t4 <- bind.tip(t3, tip.label = 'Macaca_gn2', edge.length = 1e-06, where = macaca_mrc3) +m3 <- states +m3 <- rbind(m3, matrix(c(A = 1, B = 0, C = 0), nrow = 1, byrow = TRUE)) +m3 <- rbind(m3, matrix(c(A = 0, B = 1, C = 0), nrow = 1, byrow = TRUE)) +rownames(m3)[nrow(m3) - 1] <- 'Macaca_gn1' +rownames(m3)[nrow(m3)] <- 'Macaca_gn2' +m3[grep('^Macaca_[^g][^n][^\\d].*', rownames(m3)),] <- rep(1 / ncol(m3), length(macaca) * 3) +fit3 <- fitMk( + tree = t4, x = m3, model = myMod, pi = myPi, lik.func = "pruning", + logscale = TRUE +) +ace3 <- ancr(fit3, tips=TRUE) +plot(ace3, args.plotTree = list(direction = "upwards")) +``` + +# Session information + +```{r} +sessioninfo::session_info() +```