From 33785df8ef7244060a3d99471c69f06e5dcb2a51 Mon Sep 17 00:00:00 2001 From: GabrielNakamura Date: Tue, 26 Dec 2023 12:32:30 -0300 Subject: [PATCH] adding new PD_deficit function with tibbles --- R/PD_deficit_2.R | 81 ++++++++++++++++++++++++------------------------ _pkgdown.yml | 0 2 files changed, 40 insertions(+), 41 deletions(-) mode change 100644 => 100755 _pkgdown.yml diff --git a/R/PD_deficit_2.R b/R/PD_deficit_2.R index 0c49783..2b87f3c 100644 --- a/R/PD_deficit_2.R +++ b/R/PD_deficit_2.R @@ -1,43 +1,42 @@ +#' Calculate PD deficit at species and insertion levels +#' +#' @param phylo a phylo object +#' @param data data frame with species and insertion levels. This is the output of FishPhyloMaker function +#' +#' @return a list with two tibbles. One with results of deficit for all species and another with summary results for each category of insertion +#' @export +#' +#' @examples PD_deficit_2 <- -function (phylo, data, level = "Congeneric_insertion") -{ - if (is.null(phylo) == TRUE) { - stop("/n A phylogenetic tree must be provided") - } - if (class(phylo) != "phylo") { - stop("/n A phylo object must be provided") - } - names_exclude <- phylo$tip.label[na.omit(match(data[which(data$insertion == - "present in tree"), "s"], phylo$tip.label))] - if (all(!is.na(match(phylo$tip.label, names_exclude))) == - TRUE) { - PD_present <- sum(phylo$edge.length) - PD_level <- 0 - PD_total <- sum(phylo$edge.length) - Darwinian_deficit <- PD_level/PD_total - res <- c(PD_present, PD_level, PD_total, Darwinian_deficit) - names(res) <- c("PDintree", "PDdeficit", "PDtotal", "Darwinian_deficit") - return(res) - } - else { - exclude <- ape::drop.tip(phylo, tip = names_exclude)$tip.label - phylo_present <- ape::drop.tip(phylo, tip = exclude) - PD_present <- sum(phylo_present$edge.length) - if (length(level) == 1) { - level_exclude <- ape::drop.tip(phylo, tip = phylo$tip.label[na.omit(match(data[which(data$insertions == - level), "s"], phylo$tip.label))])$tip.label - } - if (length(level) > 1) { - level_exclude <- ape::drop.tip(phylo, data[!is.na(match(data$insertion, - level)), "s"])$tip.label - } - phylo_level <- ape::drop.tip(phylo, level_exclude) - PD_level <- sum(phylo_level$edge.length) - PD_total <- sum(phylo$edge.length) - Darwinian_deficit <- PD_level/PD_total - res <- c(PD_present, PD_level, PD_total, Darwinian_deficit) - names(res) <- c("PDintree", "PDdeficit", "PDtotal", "Darwinian_deficit") - return(res) - } -} \ No newline at end of file + function(phylo, data){ + tree_tib <- tidytree::as_tibble(phylo) # change to a tibble object + + # join with insertion data + tree_tib2 <- + tree_tib %>% + right_join(data, by = c("label" = "s")) + + # Summary stats for whole table + tibble_res_all <- + tree_tib2 %>% + group_by(insertions) %>% + add_count(name = "n.insertion.category") %>% + mutate(pd.total.category = sum(branch.length), + pd.relative.category = sum(branch.length)/sum(phylo$edge.length), + pd.relative.spp = branch.length/sum(phylo$edge.length)) %>% + rename(species = "label", + order = "o", + family = "f") %>% + select(species, order, family, insertions, branch.length, n.insertion.category, pd.total.category, pd.relative.category, pd.relative.spp) + + # preparing the output with summary stats + summary_table <- + tibble_res_all %>% + distinct(insertions, n.insertion.category, pd.total.category, pd.relative.category) + + list_res <- vector(mode = "list") + list_res$pd_deficit_long <- tibble_res_all + list_res$summary_deficit <- summary_table + return(list_res) + } \ No newline at end of file diff --git a/_pkgdown.yml b/_pkgdown.yml old mode 100644 new mode 100755