How to Stack Painted Phylogenetic Trees in R Like ggdensitree but with Colored Regimes

Question

I would like to plot a sample of phylogenetic trees where the branches are “painted” to represent regimes, similar to the result of plotSimmap from the phytools package:

However, instead of plotting a single tree, I want to stack multiple trees on top of each other in a single visualization, similar to the behavior of ggdensitree from the ggtree package:

The goal is to have a stacked plot where each tree retains its regime-specific branch colors. Is there a way to achieve this in R, possibly using ggtree, ggplot2, or another approach?

Edit

Reproducible Example: We can simulate, for example, a sample of 25 trees this way:

library(ape)
library(phytools)
library(ggtree)

set.seed(123)

trees_with_regimes <- list()

# Generate trees with mapped regimes
for (i in 1:25) {

  tree <- rtree(5)   # 5 tips
  
  # Define regimes, "A" and "B"
  regimes <- sample(c("A", "B"), size = 5, replace = TRUE)
  names(regimes) <- tree$tip.label # Assign regimes to tip labels
  
  # Map regimes onto the tree
  mapped_tree <- make.simmap(tree, regimes, model = "ER", nsim = 1)
  
  # Store the mapped tree
  trees_with_regimes[[i]] <- mapped_tree
}

If we plot each of these trees separately, they will show the regimes with different colors:

# works also with "plot" instead of "plotSimmap"
plotSimmap(trees_with_regimes[[1]], lwd  = 7, fsize = 2)

But when plotting the entire tree sample, the regime colors are lost:

# or any other function for plotting several stacked trees
ggdensitree(trees_with_regimes)

I would like the last plot to include the different colors for the two regimes.

Thanks in advance for any suggestions!

Answer 1

You can use the add parameter (a logical value indicating whether or not to add the plotted tree to the current plot (TRUE) or create a new plot (FALSE, the default)) in plotSimmap to plot your 25 Simmaps on top of each other whilst still retaining the colors.

If you want to create an ultrametric tree where all tips are equidistant from the root and branches are aligned vertically, use chronos() or force.ultrametric()

To force the alignment of the tips between tree1, tree2, and so on, you need to ensure that all trees in your list have the same tip labels in the same order.

library(ape)
#install.packages("phytools")
library(phytools)
library(ggtree)
library(ggplot2)
library(tidytree)
library(dplyr)

set.seed(123)
trees_with_regimes <- list()
# Generate trees with mapped regimes
for (i in 1:8) {
  # Generate basic tree
  tree <- rtree(5)
  
  # Make tree ultrametric with equal branch lengths
  # Option 1: Using chronos
  #tree <- chronos(tree, lambda = 0, model = "discrete")
  # OR Option 2: Using force.ultrametric
  tree <- force.ultrametric(tree)
  
  # Define regimes
  regimes <- sample(c("A", "B"), size = 5, replace = TRUE)
  names(regimes) <- tree$tip.label
  
  # Map regimes
  mapped_tree <- make.simmap(tree, regimes, model = "ER", nsim = 1)
  trees_with_regimes[[i]] <- mapped_tree
}

# Reorder all trees to align their tip labels
aligned_trees <- lapply(trees_with_regimes, function(tree) {
  reorderSimmap(tree, order = "cladewise")
})

# Function to reorder tip labels of a tree based on a reference tree
align_tips <- function(tree, reference_tip_labels) {
  # Match the tree to the reference tip labels
  tree$edge <- reorder(tree, order = "cladewise")$edge
  tree$tip.label <- reference_tip_labels
  tree
}

# Use the first tree's tip labels as the reference
reference_tip_labels <- aligned_trees[[1]]$tip.label

# Align all trees' tip labels to the reference
aligned_trees <- lapply(aligned_trees, align_tips, reference_tip_labels = reference_tip_labels)

# Function to plot multiple simmap trees stacked
plot_stacked_simmaps <- function(trees, 
                                 lwd = 7, # line width
                                 fsize = 2, # relative label font size
                                 alpha = 0.5, # alpha of additional plots
                                 colors = NULL, 
                                 showLabels = T) { # show Tip Labels or not
  # If no colors provided, use default blue and red
  if(is.null(colors)) {
    colors <- c("A" = "blue", "B" = "red")
  }
  
  # Plot the first tree normally
  plotSimmap(trees[[1]], lwd = lwd, fsize = fsize, colors = colors, ftype = if(!showLabels) "off" else "reg")
  
  # Add subsequent trees with lower opacity
  if(length(trees) > 1) {
    for(i in 2:length(trees)) {
      # Colors with transparency
      transparent_colors <- sapply(colors, function(x) {
        rgb_vals <- col2rgb(x)
        rgb(rgb_vals[1], rgb_vals[2], rgb_vals[3], alpha = alpha * 255, maxColorValue = 255)
      })
      
      # Add tree to existing plot
      plotSimmap(trees[[i]], 
                 add = TRUE, 
                 lwd = lwd, 
                 fsize = fsize,
                 colors = transparent_colors,
                 ftype = if(!showLabels) "off" else "reg")
    }
  }
}

# Plot the aligned stacked simmaps
plot_stacked_simmaps(aligned_trees, lwd = 7, fsize = 2, showLabels = T)

Res