plotting anno_tax_density with 2 groups?

[ctekellogg]

Hi, I was wondering if it is possible to have a density taxAnnotation where you plot the histograms for two subsets of the data (e.g. the two sample categories) on a single anno_tax_density plot. I've managed to plot them side by side:

tax_anno = taxAnnotation(
      `Derm Prop.Abun. Log2` = anno_tax_density(data = derm, taxa = rownames(tax_table(derm)),
        joyplot_scale = 1, gp = grid::gpar(fill = "#F65A6D", alpha = 0.2),
        trans = function(ps) {
          ps %>%
            tax_transform("log2", zero_replace = "halfmin")
        }),
      `Pis Prop.Abun. Log2` = anno_tax_density(data = pis, taxa = rownames(tax_table(pis)),
        joyplot_scale = 1, gp = grid::gpar(fill = "#6a51a3", alpha = 0.2),
        trans = function(ps) {
          ps %>%
            tax_transform("log2", zero_replace = "halfmin")
        }))

But ideally these would overlap to show the different abundance distributions of the microbial genera in this dataset.

Thank you!!

[david-barnett]

These heatmaps and annotations are made with ComplexHeatmap, there doesn't seem to be a native way to group the density plots.

https://jokergoo.github.io/ComplexHeatmap-reference/book/heatmap-annotations.html#density-annotation

I had an AI have a go at rolling something custom that meets your needs, and this solution does seem to work and the approach makes sense but it's not a straightforward piece of code, so please check carefully that the densities match your side-by-side version :)

• sadly i don't have time to give this further human attention myself, good luck!

library(microViz)
library(phyloseq)
library(dplyr)
library(ComplexHeatmap)
library(grid)

data("ibd", package = "microViz")

# Prepare an example object similar to the microViz heatmaps docs
psq <- ibd %>%
  tax_filter(min_prevalence = 5) %>%
  tax_mutate(Species = NULL) %>%
  tax_fix() %>%
  tax_agg("Genus")

# Pick the two most frequent DiseaseState levels
sd <- data.frame(sample_data(psq))
grps <- names(sort(table(sd$DiseaseState), decreasing = TRUE))[1:2]
psq_a <- subset_samples(psq, DiseaseState == grps[1])
psq_b <- subset_samples(psq, DiseaseState == grps[2])

# Choose taxa to display
taxa <- tax_top(psq, n = 20, rank = "Genus")

# Transformation for the density annotation (matches microViz heatmaps style)
trans_fun <- function(ps) {
  ps %>%
    tax_transform("compositional", zero_replace = 1) %>%
    tax_transform("log2", zero_replace = "halfmin", chain = TRUE)
}

# Compute shared x-axis limits after applying the same transform and "only_detected" rule
calc_x <- function(ps, taxa, undetected = 0, only_detected = TRUE) {
  otu_counts <- as(otu_table(ps), "matrix")
  if (taxa_are_rows(ps)) otu_counts <- t(otu_counts)
  keep <- if (only_detected) otu_counts[, taxa, drop = FALSE] > undetected else TRUE

  ps_t <- trans_fun(ps)
  otu_t <- as(otu_table(ps_t), "matrix")
  if (taxa_are_rows(ps_t)) otu_t <- t(otu_t)

  x <- otu_t[, taxa, drop = FALSE]
  x[!keep] <- NaN
  x
}

x_a <- calc_x(psq_a, taxa)
x_b <- calc_x(psq_b, taxa)
xlim_shared <- range(c(x_a, x_b), finite = TRUE)

# Optional: enforce a shared vertical scale across the two layers (recommended for fair overlays)
max_dens <- function(xmat) {
  apply(xmat, 2, function(v) {
    v <- v[is.finite(v)]
    if (length(v) < 2) return(NA_real_)
    max(stats::density(v)$y)
  })
}
max_density_shared <- max(c(max_dens(x_a), max_dens(x_b)), finite = TRUE)

# Build two density layers with identical scales, then overlay them
layer_a <- anno_tax_density(
  data = psq_a, taxa = taxa, which = "row",
  trans = trans_fun,
  xlim = xlim_shared,
  max_density = max_density_shared,
  gp = gpar(fill = "#F65A6D", alpha = 0.25),
  border = TRUE, axis = TRUE
)

layer_b <- anno_tax_density(
  data = psq_b, taxa = taxa, which = "row",
  trans = trans_fun,
  xlim = xlim_shared,
  max_density = max_density_shared,
  gp = gpar(fill = "#6a51a3", alpha = 0.25),
  border = FALSE, axis = FALSE
)

overlap_anno <- AnnotationFunction(
  fun = function(index, ...) {
    layer_a@fun(index, ...)
    layer_b@fun(index, ...)
  },
  var_import = list(layer_a = layer_a, layer_b = layer_b),
  which = "row",
  n = length(taxa),
  subsettable = TRUE,
  width = unit(30, "mm")
)

# Draw a composition heatmap with the overlaid density tax annotation
psq %>%
  tax_transform("compositional", zero_replace = 1) %>%
  comp_heatmap(
    taxa = taxa,
    tax_anno = taxAnnotation(
      "test" = overlap_anno
    )
  )

[ctekellogg]

Thanks! I tried to use AI to tackle this as well but the annotation didn't quite fit as nicely (it was smaller in height than the height of the heatmap) so your prompt worked better than mine. I'll give this a shot and fiddle with it a bit. Thanks!

[david-barnett]

great :) and if you're curious, I used claude code with opus 4.6 - i'm finding it works quite nicely with R