finally fixing #30, #79

NAMESPACE

@@ -42,6 +42,7 @@ importFrom(methods,is)
 importFrom(parallel,mclapply)
 importFrom(scistreer,annotate_tree)
 importFrom(scistreer,get_mut_graph)
+importFrom(scistreer,ladderize)
 importFrom(scistreer,mut_to_tree)
 importFrom(scistreer,perform_nni)
 importFrom(scistreer,score_tree)

R/RcppExports.R

@@ -25,10 +25,6 @@ viterbi_compute <- function(log_delta, logprob, logPi, n, m, nu, z) {
     .Call('_numbat_viterbi_compute', PACKAGE = 'numbat', log_delta, logprob, logPi, n, m, nu, z)
 }
 
-node_depth <- function(ntip, e1, e2, nedge, xx, method) {
-    .Call('_numbat_node_depth', PACKAGE = 'numbat', ntip, e1, e2, nedge, xx, method)
-}
-
 roman2int_internal <- function(letters, nchar) {
     .Call('_numbat_roman2int_internal', PACKAGE = 'numbat', letters, nchar)
 }

R/main.R

@@ -8,7 +8,7 @@
 #' @import tidygraph
 #' @import ggplot2
 #' @import ggraph
-#' @importFrom scistreer perform_nni get_mut_graph score_tree annotate_tree mut_to_tree to_phylo
+#' @importFrom scistreer perform_nni get_mut_graph score_tree annotate_tree mut_to_tree ladderize to_phylo
 #' @importFrom ggtree %<+%
 #' @importFrom methods is as
 #' @importFrom igraph vcount ecount E V V<- E<-

R/trees.R

@@ -14,60 +14,6 @@ upgma <- function(D, method = "average", ...) {
   result
 }
 
-#' from ape
-#' @keywords internal
-ladderize <- function(phy, right = TRUE) {
-    desc_fun <- function(x) {
-        parent <- x[, 1]
-        children <- x[, 2]
-        res <- vector("list", max(x))
-        for (i in seq_along(parent)) res[[parent[i]]] <- c(res[[parent[i]]], children[i])
-        return(res)
-    }
-
-    if(!is.null(phy$edge.length)){
-        el <- numeric(max(phy$edge))
-        el[phy$edge[, 2]] <- phy$edge.length
-    }
-
-    nb.tip <- length(phy$tip.label)
-    nb.node <- phy$Nnode
-    nb.edge <- dim(phy$edge)[1]
-
-    phy <- reorder(phy, "postorder")
-    N <- node_depth(as.integer(nb.tip), as.integer(phy$edge[, 1]), as.integer(phy$edge[, 2]),
-            as.integer(nb.edge), double(nb.tip + nb.node), 1L)
-
-    ii <- order(x <- phy$edge[,1], y <- N[phy$edge[,2]], decreasing = right)
-    desc <- desc_fun(phy$edge[ii,])
-
-    tmp <- integer(nb.node)
-    new_anc <- integer(nb.node)
-    new_anc[1] <- tmp[1] <- nb.tip + 1L
-    k <- nb.node
-    pos <- 1L
-
-    while(pos > 0L && k > 0){
-        current <- tmp[pos]
-        new_anc[k] <- current
-        k <- k - 1L
-        dc <- desc[[current]]
-        ind <- (dc > nb.tip)
-        if(any(ind)){
-            l <- sum(ind)
-            tmp[pos -1L + seq_len(l)] <-  dc[ind]
-            pos <- pos + l - 1L
-        }
-        else pos <- pos - 1L
-    }
-    edge <- cbind(rep(new_anc, lengths(desc[new_anc])), unlist(desc[new_anc]))
-    phy$edge <- edge
-    if(!is.null(phy$edge.length)) phy$edge.length <- el[edge[,2]]
-    attr(phy, "order") <- "postorder"
-    phy <- reorder(phy, "cladewise")
-    phy
-}
-
 #' Mark the tumor lineage of a phylogeny
 #' @param gtree tbl_graph Single-cell phylogeny
 #' @return tbl_graph Phylogeny annotated with tumor versus normal compartment
@@ -356,4 +302,12 @@ get_move_opt = function(G, l_matrix) {
         head(1)
 
     return(move_opt)
-}
+}
+
+#' Get ordered tips from a tree
+#' @keywords internal
+get_ordered_tips = function(tree) {
+    is_tip <- tree$edge[,2] <= length(tree$tip.label)
+    ordered_tips <- tree$edge[is_tip, 2]
+    tree$tip.label[ordered_tips]
+}

R/vis.R

@@ -519,46 +519,46 @@ plot_phylo_heatmap = function(
     # inspect gtree object
     if (is(gtree, "tbl_graph")) {
         gtree = gtree %>% activate(edges) %>% mutate(length = ifelse(leaf, tip_length, length))
-        tree_obj = gtree %>% to_phylo()
+        tree_obj = gtree %>% to_phylo() %>% ladderize(right = FALSE)
     } else {
         tree_obj = gtree
         tvn_line = FALSE
         clone_bar = FALSE
     }
 
-    # plot phylogeny
-    if (show_phylo) {
+    cell_order = get_ordered_tips(tree_obj)
 
-    	p_tree = tree_obj %>%
-    		ggtree::ggtree(ladderize = TRUE, size = branch_width) +
-    		theme(
-    			plot.margin = margin(0,1,0,0, unit = 'mm'),
-    			axis.title.x = element_blank(),
-    			axis.ticks.x = element_blank(),
-    			axis.text.x = element_blank(),
-    			axis.line.y = element_blank(),
-    			axis.ticks.y = element_blank(),
-    			# axis.text.y = element_text(size = 5)
-    			axis.text.y = element_blank(),
-    			panel.background = element_rect(fill = "transparent",colour = NA),
-    			plot.background = element_rect(fill = "transparent", color = NA)
-    		) +
-    		guides(color = 'none')
-
-    	if (root_edge) {
-    		p_tree = p_tree + ggtree::geom_rootedge(size = branch_width)
-    	}
-
-    	# order the cells
-    	cell_order = p_tree$data %>% filter(isTip) %>% arrange(y) %>% pull(label)
-    } else {
-    	p_tree = gtree %>%
-    		ape::as.phylo() %>%
-    		ape::ladderize(right = FALSE)
+    if (show_phylo) {
+        p_tree = tryCatch(expr = {
+            p_tree = tree_obj %>%
+                ggtree::ggtree(ladderize = FALSE, size = branch_width) +
+                theme(
+                    plot.margin = margin(0,1,0,0, unit = 'mm'),
+                    axis.title.x = element_blank(),
+                    axis.ticks.x = element_blank(),
+                    axis.text.x = element_blank(),
+                    axis.line.y = element_blank(),
+                    axis.ticks.y = element_blank(),
+                    axis.text.y = element_blank(),
+                    panel.background = element_rect(fill = "transparent",colour = NA),
+                    plot.background = element_rect(fill = "transparent", color = NA)
+                ) +
+                guides(color = 'none')
+
+            if (root_edge) {
+                p_tree = p_tree + ggtree::geom_rootedge(size = branch_width)
+            }
+        },
+        error = function(e) {
+            print(e)
+            log_warn("Plotting phylogeny failed, continuing..")
 
-    	cell_order = rev(p_tree$tip.label)
+            p_tree = ggplot() + theme_void()
 
-    	p_tree = ggplot() + theme_void()
+            return(p_tree)
+        })
+    } else {
+        p_tree = ggplot() + theme_void()
     }
 
     joint_post = joint_post %>%

src/RcppExports.cpp

@@ -98,22 +98,6 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
-// node_depth
-NumericVector node_depth(int ntip, NumericVector e1, NumericVector e2, int nedge, NumericVector xx, int method);
-RcppExport SEXP _numbat_node_depth(SEXP ntipSEXP, SEXP e1SEXP, SEXP e2SEXP, SEXP nedgeSEXP, SEXP xxSEXP, SEXP methodSEXP) {
-BEGIN_RCPP
-    Rcpp::RObject rcpp_result_gen;
-    Rcpp::RNGScope rcpp_rngScope_gen;
-    Rcpp::traits::input_parameter< int >::type ntip(ntipSEXP);
-    Rcpp::traits::input_parameter< NumericVector >::type e1(e1SEXP);
-    Rcpp::traits::input_parameter< NumericVector >::type e2(e2SEXP);
-    Rcpp::traits::input_parameter< int >::type nedge(nedgeSEXP);
-    Rcpp::traits::input_parameter< NumericVector >::type xx(xxSEXP);
-    Rcpp::traits::input_parameter< int >::type method(methodSEXP);
-    rcpp_result_gen = Rcpp::wrap(node_depth(ntip, e1, e2, nedge, xx, method));
-    return rcpp_result_gen;
-END_RCPP
-}
 // roman2int_internal
 int roman2int_internal(Rcpp::StringVector letters, int nchar);
 RcppExport SEXP _numbat_roman2int_internal(SEXP lettersSEXP, SEXP ncharSEXP) {
@@ -176,7 +160,6 @@ static const R_CallMethodDef CallEntries[] = {
     {"_numbat_likelihood_compute", (DL_FUNC) &_numbat_likelihood_compute, 5},
     {"_numbat_forward_backward_compute", (DL_FUNC) &_numbat_forward_backward_compute, 5},
     {"_numbat_viterbi_compute", (DL_FUNC) &_numbat_viterbi_compute, 7},
-    {"_numbat_node_depth", (DL_FUNC) &_numbat_node_depth, 6},
     {"_numbat_roman2int_internal", (DL_FUNC) &_numbat_roman2int_internal, 2},
     {"_numbat_fit_lnpois_cpp", (DL_FUNC) &_numbat_fit_lnpois_cpp, 3},
     {"_numbat_poilog1", (DL_FUNC) &_numbat_poilog1, 3},

src/misc.cpp

@@ -1,38 +1,6 @@
 #include <Rcpp.h>
 using namespace Rcpp;
 
-// Based on https://github.com/cran/ape/blob/390386e67f9ff6cd8e6e523b7c43379a1551c565/src/plot_phylo.c
-// [[Rcpp::export]]
-NumericVector node_depth(int ntip, NumericVector e1, NumericVector e2,
-        int nedge, NumericVector xx, int method)
-/* method == 1: the node depths are proportional to the number of tips
-   method == 2: the node depths are evenly spaced */
-{
-
-    int i;
-
-    /* First set the coordinates for all tips */
-    for (i = 0; i < ntip; i++) xx[i] = 1;
-
-    /* Then compute recursively for the nodes; we assume `xx' has */
-    /* been initialized with 0's which is true if it has been */
-    /* created in R (the tree must be in pruningwise order) */
-    if (method == 1) {
-        for (i = 0; i < nedge; i++)
-            xx[e1[i] - 1] = xx[e1[i] - 1] + xx[e2[i] - 1];
-    } else { /* *method == 2 */
-        for (i = 0; i < nedge; i++) {
-            /* if a value > 0 has already been assigned to the ancestor
-               node of this edge, check that the descendant node is not
-               at the same level or more */
-            if (xx[e1[i] - 1])
-            if (xx[e1[i] - 1] >= xx[e2[i] - 1] + 1) continue;
-            xx[e1[i] - 1] = xx[e2[i] - 1] + 1;
-        }
-    }
-    return xx;
-}
-
 // gtools is orphaned
 // [[Rcpp::export]]
 int roman2int_internal(Rcpp::StringVector letters, int nchar) {