adding new PD_deficit function with tibbles

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)
-  }
-}
+  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)
+  }