diff --git a/NAMESPACE b/NAMESPACE index c8bbd0a4..266aba08 100644 --- a/NAMESPACE +++ b/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) diff --git a/R/RcppExports.R b/R/RcppExports.R index eb4cb156..de4ce30d 100644 --- a/R/RcppExports.R +++ b/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) } diff --git a/R/main.R b/R/main.R index 278de604..07ff6467 100644 --- a/R/main.R +++ b/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<- diff --git a/R/trees.R b/R/trees.R index 02691f2c..b816cb20 100644 --- a/R/trees.R +++ b/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) -} \ No newline at end of file +} + +#' 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] +} diff --git a/R/vis.R b/R/vis.R index b068c17b..c2cc36a7 100644 --- a/R/vis.R +++ b/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 %>% diff --git a/man/get_ordered_tips.Rd b/man/get_ordered_tips.Rd new file mode 100644 index 00000000..b9d42f6d --- /dev/null +++ b/man/get_ordered_tips.Rd @@ -0,0 +1,12 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/trees.R +\name{get_ordered_tips} +\alias{get_ordered_tips} +\title{Get ordered tips from a tree} +\usage{ +get_ordered_tips(tree) +} +\description{ +Get ordered tips from a tree +} +\keyword{internal} diff --git a/man/ladderize.Rd b/man/ladderize.Rd deleted file mode 100644 index c3113740..00000000 --- a/man/ladderize.Rd +++ /dev/null @@ -1,12 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/trees.R -\name{ladderize} -\alias{ladderize} -\title{from ape} -\usage{ -ladderize(phy, right = TRUE) -} -\description{ -from ape -} -\keyword{internal} diff --git a/man/plot_phylo_heatmap.Rd b/man/plot_phylo_heatmap.Rd index 831c4826..9f1cae44 100644 --- a/man/plot_phylo_heatmap.Rd +++ b/man/plot_phylo_heatmap.Rd @@ -27,12 +27,11 @@ plot_phylo_heatmap( annot_bar_width = 0.25, clone_bar_width = 0.25, bar_label_size = 7, + tvn_line = TRUE, clone_line = FALSE, - superclone = FALSE, exclude_gap = FALSE, root_edge = TRUE, - raster = FALSE, - show_phylo = TRUE + raster = FALSE ) } \arguments{ @@ -76,17 +75,15 @@ plot_phylo_heatmap( \item{bar_label_size}{numeric Size of sidebar text labels} -\item{clone_line}{logical Whether to display borders for clones in the heatmap} +\item{tvn_line}{logical Whether to draw line separating tumor and normal cells} -\item{superclone}{logical Wehether to display superclones in the clone bar} +\item{clone_line}{logical Whether to display borders for clones in the heatmap} \item{exclude_gap}{logical Whether to mark gap regions} \item{root_edge}{logical Whether to plot root edge} \item{raster}{logical Whether to raster images} - -\item{show_phylo}{logical Whether to display phylogeny on y axis} } \value{ ggplot panel diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 3c0a2e91..2db09143 100644 --- a/src/RcppExports.cpp +++ b/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}, diff --git a/src/misc.cpp b/src/misc.cpp index 64f36ef5..92334904 100644 --- a/src/misc.cpp +++ b/src/misc.cpp @@ -1,38 +1,6 @@ #include 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) {