How to extract row/column order after clustering in a heatmap?

[CaoWZ254]

The ggalign package has been very beneficial for my actual work, but there's a confusion that I haven't found a good solution for.

In bioinformatics analysis, there is a type of heatmap. Its main body uses relative abundance for representation, while the side annotations display the absolute abundance of each gene to express richer gene abundance information. My question is: After clustering using align_dendro(), the column order of the main heatmap changes accordingly. At this point, how can one conveniently extract the changed column order and apply it to the absolute abundance plot?

First, the ComplexHeatmap package provides simple parameters to achieve this purpose, but I am not familiar with base R graphics grammar:

library(ComplexHeatmap)

set.seed(100)
expr_matrix <- matrix(sample(c(0,1), 400, replace = TRUE),
                     nrow = 20, ncol = 20)

rownames(expr_matrix) <- paste0("Gene", 1:20)
colnames(expr_matrix) <- paste0("Sample", 1:20)

# Calculate z-score
zscore_matrix <- t(scale(t(expr_matrix)))

# Calculate total expression per row
expr_counts <- rowSums(expr_matrix)

# Create heatmap
zscore_matrix |>
 Heatmap(
   name = " ",
   rect_gp = gpar(col = "white", lwd = 2),
   row_names_gp = gpar(fontsize = 8),
   column_names_rot = 45,
   column_names_gp = gpar(fontsize = 8),
   cluster_rows = TRUE,
   cluster_columns = TRUE,
   show_row_dend = TRUE,
   show_column_dend = TRUE,
   right_annotation = rowAnnotation(
     ExprCount = anno_barplot(
       expr_counts,
       bar_width = 0.8,
       gp = gpar(fill = "steelblue"),
       axis_param = list(
         gp = gpar(fontsize = 6)
       )
     )
   )
 ) -> ht

Secondly, I wrote my own set of code to apply to the ggalign workflow:

library(tidyverse)
library(ggalign)

get_gene_order <- function(
    zscore_matrix,
    method = "complete",
    distance = "euclidean"
) {
  dist_mat <- dist(zscore_matrix, method = distance)
  hc_row <- hclust(dist_mat, method = method)
  hc_row$labels[hc_row$order]
}

create_row_barplot <- function(
    expr_matrix,
    gene_order,
    fill_color = "steelblue"
) {
  Gene <- rownames(expr_matrix)
  Count <- rowSums(expr_matrix)

  tibble::tibble(Gene, Count) |>
    mutate(Gene = forcats::fct_relevel(Gene, gene_order)) |>
    ggplot(aes(x = Count, y = Gene)) +
    geom_bar(stat = "identity", fill = fill_color, color = "black") +
    theme_minimal() +
    scale_x_continuous(labels = scales::label_number(accuracy = 1)) +
    theme(
      axis.ticks.y = element_blank(),
      axis.text.y = element_blank(),
      axis.text.x = element_text(angle = 90, hjust = 1),
      axis.title = element_blank(),
      panel.grid = element_blank(),
      panel.border = element_rect(color = "black")
    )
}


ggene_heatmap <- function(
    expr_matrix,
    bar_fill = "steelblue",
    cluster_palette_row = "Set1",
    cluster_palette_col = "Dark2",
    dendro_size = 0.2,
    tile_color = "black",
    tile_width = 0.9,
    tile_height = 0.9,
    axis_text_x_angle = 45,
    cluster_method = "complete",
    cluster_distance = "euclidean"
) {
  # Step 1: Calculate z-score
  zscore_matrix <- t(scale(t(expr_matrix)))

  # Step 2: Get row clustering order
  gene_order <- get_gene_order(
    zscore_matrix,
    method = cluster_method,
    distance = cluster_distance
  )

  # Step 3: Draw right-side barplot
  p_bar <- create_row_barplot(expr_matrix, gene_order, fill_color = bar_fill)

  # Step 4: Assemble ggalign heatmap
  zscore_matrix |>
    ggheatmap(filling = NULL) +
    geom_tile(
      aes(fill = value),
      color = tile_color,
      width = tile_width,
      height = tile_height
    ) +
    scale_fill_viridis_c(name = "") +
    theme_minimal() +
    theme(
      axis.text.x = element_text(angle = axis_text_x_angle, hjust = 1),
      axis.title = element_blank()
    ) +
    # Top dendro
    anno_top(size = dendro_size) +
    align_dendro() +
    scale_color_brewer(palette = cluster_palette_col, guide = "none") +
    theme(
      axis.title = element_blank(),
      axis.ticks = element_blank(),
      axis.text = element_blank()
    ) +
    # Left dendro
    anno_left(size = dendro_size) +
    align_dendro() +
    scale_color_brewer(palette = cluster_palette_row, guide = "none") +
    theme(
      axis.title = element_blank(),
      axis.ticks = element_blank(),
      axis.text = element_blank()
    ) +
    # Right bar chart
    quad_switch(position = "right", size = dendro_size) +
    p_bar
}

ggene_heatmap(expr_matrix)

Although the above code can run, I believe this is not the best solution for the ggalign workflow. Therefore, I sincerely ask if you have a method to achieve either of the following goals:

1. Extract the row/column order of the heatmap after clustering.
2. Precisely align the row/column order of two plots when their axis labels are completely identical.

I look forward to your reply!

[Yunuuuu]

Sorry for the delay, I missed your message.

How to extract row/column order after clustering in a heatmap?

ggalign uses ggalign_stat to extract statistics from the layout. For heatmap layouts, the first argument specifies the position of the annotation where the statistics are applied, and the second argument indicates which box (note that ggalign uses box to refer to entities in the layout, and some entities may not contain any plots) the statistics belong to.

library(ggalign)
#> Loading required package: ggplot2
#> ========================================
#> ggalign version 1.2.0
#>
#> If you use it in published research, please cite:
#> Peng, Y.; Jiang, S.; Song, Y.; et al. ggalign: Bridging the Grammar of Graphics and Biological Multilayered Complexity. Advanced Science. 2025. doi:10.1002/advs.202507799
#> ========================================
set.seed(100)
expr_matrix <- matrix(
    sample(c(0, 1), 400, replace = TRUE),
    nrow = 20, ncol = 20
)
rownames(expr_matrix) <- paste0("Gene", 1:20)
colnames(expr_matrix) <- paste0("Sample", 1:20)

# Calculate z-score
zscore_matrix <- t(scale(t(expr_matrix)))

myheatmap <- ggheatmap(zscore_matrix) +
    scale_fill_viridis_c(NULL) +
    theme(axis.text.x = element_text(angle = -60, hjust = 0)) +
    anno_top(size = 0.2) +
    align_dendro() +
    anno_left(size = 0.2) +
    align_dendro()
ggalign_stat(myheatmap, "top", 1)
#>
#> Call:
#> stats::hclust(d = d, method = method)
#>
#> Cluster method   : complete
#> Distance         : euclidean
#> Number of objects: 20
ggalign_stat(myheatmap, "left", 1)
#>
#> Call:
#> stats::hclust(d = d, method = method)
#>
#> Cluster method   : complete
#> Distance         : euclidean
#> Number of objects: 20

^{Created on 2026-01-15 with reprex v2.1.0}
~

However, in your example, you shouldn't rely on the extracted ordering. ggalign can add a barplot, and as long as the initial data ordering is matched, ggalign will reorder them to match the layout ordering or groups. Additionally, ggalign allows setting data in each box independently, so you can apply the absolute abundance data directly instead of inheriting from the layout.

library(ggalign)
#> Loading required package: ggplot2
#> ========================================
#> ggalign version 1.2.0
#> 
#> If you use it in published research, please cite: 
#> Peng, Y.; Jiang, S.; Song, Y.; et al. ggalign: Bridging the Grammar of Graphics and Biological Multilayered Complexity. Advanced Science. 2025. doi:10.1002/advs.202507799
#> ========================================
set.seed(100)
expr_matrix <- matrix(
    sample(c(0, 1), 400, replace = TRUE),
    nrow = 20, ncol = 20
)
rownames(expr_matrix) <- paste0("Gene", 1:20)
colnames(expr_matrix) <- paste0("Sample", 1:20)

# Calculate z-score
zscore_matrix <- t(scale(t(expr_matrix)))

ggene_heatmap <- function(expr_matrix,
                          bar_fill = "steelblue",
                          cluster_palette_row = "Set1",
                          cluster_palette_col = "Dark2",
                          dendro_size = 0.2,
                          tile_color = "black",
                          tile_width = 0.9,
                          tile_height = 0.9,
                          axis_text_x_angle = 45,
                          cluster_method = "complete",
                          cluster_distance = "euclidean") {
    ggheatmap(t(scale(t(expr_matrix))), filling = NULL) +
        geom_tile(
            aes(fill = value),
            color = tile_color,
            width = tile_width,
            height = tile_height
        ) +
        scale_fill_viridis_c(name = "") +
        theme_minimal() +
        theme(
            axis.text.x = element_text(
                angle = axis_text_x_angle,
                hjust = 1
            ),
            axis.title = element_blank()
        ) +
        # Top dendro
        anno_top(size = dendro_size) +
        align_dendro() +
        scale_color_brewer(palette = cluster_palette_col, guide = "none") +
        theme(
            axis.title = element_blank(),
            axis.ticks = element_blank(),
            axis.text = element_blank()
        ) +
        # Left dendro
        anno_left(size = dendro_size) +
        align_dendro() +
        scale_color_brewer(palette = cluster_palette_row, guide = "none") +
        theme(
            axis.title = element_blank(),
            axis.ticks = element_blank(),
            axis.text = element_blank()
        ) +
        # Right bar chart
        anno_right(size = dendro_size) +

        # I'll explain the `ggalign` and `scheme_data` deep in the ending
        ggalign(expr_matrix, aes(count, .discrete_y)) +
        scheme_data(function(data) {
            dplyr::summarise(data, count = sum(value), .by = .discrete_y)
        }) +

        geom_bar(stat = "identity", fill = bar_fill, color = "black") +
        theme_minimal() +
        scale_x_continuous(labels = scales::label_number(accuracy = 1)) +
        theme(
            axis.ticks.y = element_blank(),
            axis.text.y = element_blank(),
            axis.text.x = element_text(angle = 90, hjust = 1),
            axis.title = element_blank(),
            panel.grid = element_blank(),
            panel.border = element_rect(color = "black")
        )
}

ggene_heatmap(expr_matrix)

^{Created on 2026-01-15 with reprex v2.1.0}

ggalign() function accepts a matrix and treats the rows as observations. When building the plot, it will match the observations of the heatmap. For row annotations, this means it will match the heatmap data rows, and for column annotations, it will match the heatmap data columns.

When plotting, ggalign() translates the input matrix into a data frame using the fortify_data_frame() function and adds the following columns to the data frame:

• .panel: The panel for the aligned axis. Refers to the x-axis for vertical stack_layout() (including top and bottom annotations), and the y-axis for horizontal stack_layout() (including left and right annotations).

• .names (row names) and .index (row index): Character names (if available) and the integer index of the original data.

• .x/.y and .discrete_x/.discrete_y: Integer indices for x/y coordinates, and a factor of the data labels (only applicable when names exist). Depending on the axis to be aligned, for row annotations, it will be .y and .discrete_y; for column annotations, it will be .x and .discrete_x.

You should use .x/.y, or .discrete_x/.discrete_y as the x/y mapping to make the plot aligned.

scheme_data() is a function provided by ggalign that allows users to transform plot data into any desired format:

ggalign(expr_matrix, aes(count, .discrete_y)) +
    scheme_data(function(data) {
        dplyr::summarise(data, count = sum(value), .by = .discrete_y)
    })

Conceptually, the above code is equivalent to the following (although the actual internal implementation is more complex to ensure proper integration and the addition of necessary columns for alignment):

scheme_data_fn <- function(data) {
    dplyr::summarise(data, count = sum(value), .by = .discrete_y)
}

# `index` comes from the layout ordering (in your case, the phylogenetic tree)
data <- fortify_data_frame(data[index, , drop = FALSE])

# facet panel grouping information is added at this stage
# afterwards, the `scheme_data_fn` is applied
data <- scheme_data_fn(data)
ggplot(data, aes(count, .discrete_y))

[CaoWZ254]

A very intuitive explanation! To be honest, I spent a long time studying your tutorial book, but I just couldn't locate the scheme_data() function, because that chapter involves many explanations of low-level code and interactions between different parts of the codebase, which I couldn’t fully follow. I have a small suggestion: would it be possible to add a simplified figure to the “Heatmap” section of the gallery to illustrate the intuitive use of .discrete_x/.discrete_y and scheme_data?