Harden abundance heatmaps

Author: wolski

NEWS.md

@@ -9,6 +9,8 @@
 - `setup_analysis()` now stops with an informative error (listing the offending keys) when a hierarchy-key/sample combination has more than one observation, instead of silently returning a different-schema count table that crashed downstream. Pass `debug = TRUE` to recover the old behaviour and return the count table for inspection.
 - Removed the unused `impute_with_zcomp()`, `estimate_lod_global()`, and `function_lod_quantile()` exports (and the `zCompositions` dependency). For missing-value imputation use `AggregateLimpa$new(lfqdata, impute_only = TRUE)$aggregate()`.
 - Hardened `plot_pca()`: errors early on duplicated sample names, an all-missing matrix, or too few samples instead of returning `NULL` (which broke `pca_plotly()`); joins scores to annotation with an explicit `by`; makes `prcomp(center = TRUE, scale. = FALSE)` explicit.
+- Hardened abundance heatmaps for sparse significant-feature subsets: when row or column distances are non-finite because of missing values, `plot_heatmap()` now falls back to the input order instead of returning a `ComplexHeatmap` object that fails during drawing.
+- `LFQDataPlotter$heatmap()` now shows only the `top_n` most variable features (default 1000), ranked by the prolfqua per-feature statistic (CV for untransformed data, sd for transformed, via `LFQDataStats`). Row clustering uses `stats::hclust`, which errors above 65536 features, so peptide-list / entrapment searches with tens of thousands of degenerate protein groups no longer crash the QC heatmap. Pass `top_n = NULL` (or `Inf`) to keep every feature.
 
 # prolfqua 1.6.1
 

R/LFQDataPlotter.R

@@ -1,3 +1,31 @@
+# Keep the `top_n` most variable features for the intensity heatmap.
+#
+# Ranks features by their prolfqua variability statistic (CV for untransformed
+# data, sd for transformed data; see LFQDataStats) and keeps the most variable.
+# The heatmap clusters rows with stats::hclust, which errors above 65536
+# features, so this keeps clustering feasible (and the plot legible) for
+# peptide-list / entrapment searches with tens of thousands of degenerate
+# protein groups. `top_n = NULL` / non-finite keeps all features.
+# Returns a subset of wide$data in the original row order.
+.select_most_variable_features <- function(lfqdata, wide, top_n) {
+  mat <- wide$data
+  if (is.null(top_n) || !is.finite(top_n) || nrow(mat) <= top_n) {
+    return(mat)
+  }
+  lfqstats <- lfqdata$get_Stats(stats = "all")
+  metric <- lfqstats$stat # "CV" (untransformed) or "sd" (transformed)
+  keys <- c(lfqdata$hierarchy_keys(), lfqdata$isotope_label())
+  # wide$rowdata is in the same row order as the matrix, so joining the
+  # per-feature statistic onto it aligns the ranking with the matrix rows.
+  ranked <- dplyr::left_join(
+    wide$rowdata,
+    dplyr::select(lfqstats$stats(), dplyr::all_of(c(keys, metric))),
+    by = keys
+  )
+  keep <- utils::head(order(ranked[[metric]], decreasing = TRUE, na.last = TRUE), top_n)
+  mat[sort(keep), , drop = FALSE]
+}
+
 #' LFQDataPlotter ----
 #' Create various visualization of the LFQdata
 #' @return An R6 class generator.
@@ -97,16 +125,32 @@ LFQDataPlotter <- R6::R6Class(
     #' Without the z-scoring, the proteins would group according
     #' to their abundance, e.g., high abundant proteins would be one cluster.
     #'
+    #' Only the \code{top_n} most variable features are shown. Row clustering
+    #' uses \code{stats::hclust}, which errors above 65536 features, so the rows
+    #' are ranked by their variability statistic (CV for untransformed data, sd
+    #' for transformed data; see \code{\link{LFQDataStats}}) and the most
+    #' variable are kept. This keeps the heatmap feasible and legible for
+    #' peptide-list / entrapment searches with tens of thousands of features.
+    #'
     #' @param na_fraction max fraction of NA's per row
     #' @param rownames show rownames (default FALSE - do not show.)
     #' @param max_rownames_chars maximum displayed row label length
     #' @param max_sample_label_chars maximum displayed sample label length.
     #'   Labels keep their suffix because sample prefixes are often shared.
+    #' @param top_n keep the \code{top_n} most variable features (default 1000);
+    #'   \code{NULL} or \code{Inf} keeps all features.
     #' @return ComplexHeatmap::Heatmap
-    heatmap = function(na_fraction = 0.3, rownames = FALSE, max_rownames_chars = 60, max_sample_label_chars = 20) {
+    heatmap = function(
+      na_fraction = 0.3,
+      rownames = FALSE,
+      max_rownames_chars = 60,
+      max_sample_label_chars = 20,
+      top_n = 1000
+    ) {
       wide <- self$lfq$data_wide(as.matrix = TRUE)
+      data <- .select_most_variable_features(self$lfq, wide, top_n)
       fig <- prolfqua::plot_heatmap(
-        wide$data,
+        data,
         wide$annotation,
         self$lfq$factor_keys(),
         self$lfq$sample_name(),

R/tidyMS_plotting.R

@@ -115,6 +115,37 @@
   circlize::colorRamp2(c(-m, 0, m), c("green", "black", "red"))
 }
 
+.finite_heatmap_dist <- function(matrix) {
+  tryCatch(
+    {
+      distance <- stats::dist(matrix)
+      if (any(!is.finite(distance))) {
+        return(NULL)
+      }
+      distance
+    },
+    error = function(e) NULL
+  )
+}
+
+.can_cluster_heatmap_columns <- function(matrix) {
+  if (ncol(matrix) < 2) {
+    return(FALSE)
+  }
+  !is.null(.finite_heatmap_dist(t(matrix)))
+}
+
+.cluster_heatmap_rows <- function(matrix) {
+  if (nrow(matrix) < 3) {
+    return(matrix)
+  }
+  distance <- .finite_heatmap_dist(matrix)
+  if (is.null(distance)) {
+    return(matrix)
+  }
+  matrix[stats::hclust(distance)$order, , drop = FALSE]
+}
+
 # White-to-red color function for correlation heatmaps, ranging over the
 # observed values (R-squared lives in [0, 1]).
 .cor_col_fun <- function(cres, R2 = FALSE) {
@@ -500,50 +531,33 @@ plot_heatmap <- function(
   }
 
   resdata <- t(scale(t(matrix)))
-  resdataf <- prolfqua::remove_na_rows(resdata, floor(ncol(resdata) * na_fraction))
-
-  if (nrow(resdataf) >= 3) {
-    gg <- stats::hclust(stats::dist(resdataf))
-    plot_data <- resdataf[gg$order, ]
-    res <- ComplexHeatmap::Heatmap(
-      plot_data,
-      name = "row z-score",
-      col = .abundance_col_fun(plot_data),
-      na_col = .HEATMAP_NA_COL,
-      cluster_rows = FALSE,
-      cluster_columns = TRUE,
-      top_annotation = .heatmap_top_annotation(annotation, factor_keys, sample_name, colnames(plot_data)),
-      show_row_names = show_rownames,
-      show_column_names = TRUE,
-      row_labels = .truncate_plot_labels(rownames(plot_data), max_rownames_chars),
-      column_labels = .suffix_plot_labels(colnames(plot_data), max_sample_label_chars),
-      border = FALSE,
-      use_raster = FALSE,
-      heatmap_legend_param = list(title = "row z-score"),
-      ... = ...
-    )
+  na_threshold <- floor(ncol(resdata) * na_fraction)
+  keep_rows <- rowSums(is.na(resdata)) <= na_threshold
+  resdataf <- resdata[keep_rows, , drop = FALSE]
+  plot_data <- if (nrow(resdataf) >= 3) {
+    .cluster_heatmap_rows(resdataf)
   } else {
-    res <- tryCatch(
-      ComplexHeatmap::Heatmap(
-        resdata,
-        name = "row z-score",
-        col = .abundance_col_fun(resdata),
-        na_col = .HEATMAP_NA_COL,
-        cluster_rows = FALSE,
-        cluster_columns = TRUE,
-        top_annotation = .heatmap_top_annotation(annotation, factor_keys, sample_name, colnames(resdata)),
-        show_row_names = show_rownames,
-        show_column_names = TRUE,
-        row_labels = .truncate_plot_labels(rownames(resdata), max_rownames_chars),
-        column_labels = .suffix_plot_labels(colnames(resdata), max_sample_label_chars),
-        border = FALSE,
-        use_raster = FALSE,
-        heatmap_legend_param = list(title = "row z-score"),
-        ... = ...
-      ),
-      error = .error_handler # nolint object_usage_linter. defined in utilities.R
-    )
+    resdata
   }
+  cluster_columns <- .can_cluster_heatmap_columns(plot_data)
+
+  res <- ComplexHeatmap::Heatmap(
+    plot_data,
+    name = "row z-score",
+    col = .abundance_col_fun(plot_data),
+    na_col = .HEATMAP_NA_COL,
+    cluster_rows = FALSE,
+    cluster_columns = cluster_columns,
+    top_annotation = .heatmap_top_annotation(annotation, factor_keys, sample_name, colnames(plot_data)),
+    show_row_names = show_rownames,
+    show_column_names = TRUE,
+    row_labels = .truncate_plot_labels(rownames(plot_data), max_rownames_chars),
+    column_labels = .suffix_plot_labels(colnames(plot_data), max_sample_label_chars),
+    border = FALSE,
+    use_raster = FALSE,
+    heatmap_legend_param = list(title = "row z-score"),
+    ... = ...
+  )
   invisible(res)
 }
 

man/LFQDataPlotter.Rd

@@ -268,6 +268,29 @@ test_that("abundance heatmap shows missing values in light gray", {
   )
 })
 
+test_that("abundance heatmap draws sparse one-row matrices with all-NA columns", {
+  matrix <- matrix(
+    c(NA, 12, 13, NA, 12.5, 13.5),
+    nrow = 1,
+    dimnames = list("P1", c("C1", "T1", "T2", "C2", "T3", "T4"))
+  )
+  annotation <- data.frame(
+    sample = colnames(matrix),
+    group = c("control", "treated", "treated", "control", "treated", "treated")
+  )
+
+  p <- plot_heatmap(
+    matrix,
+    annotation,
+    factor_keys = "group",
+    sample_name = "sample",
+    show_rownames = TRUE
+  )
+
+  expect_s4_class(p, "Heatmap")
+  expect_no_error(ComplexHeatmap::draw(p))
+})
+
 test_that("abundance color mapping spans green, black, and red", {
   col_fun <- prolfqua:::.abundance_col_fun(matrix(c(-2, 0, 2), nrow = 1))
   cols <- grDevices::col2rgb(col_fun(c(-2, 0, 2)))
@@ -345,3 +368,38 @@ test_that("volcano_plotly does not cap small non-zero FDR values by default", {
 
   expect_equal(max(capped_y), 4)
 })
+
+test_that("heatmap keeps only the top_n most variable features", {
+  istar <- sim_lfq_data_protein_config()
+  lfq <- LFQData$new(istar$data, istar$config)
+  wide <- lfq$data_wide(as.matrix = TRUE)
+  n_all <- nrow(wide$data)
+  expect_gt(n_all, 3)
+
+  # top_n >= feature count (or NULL / Inf) keeps every feature.
+  expect_equal(nrow(prolfqua:::.select_most_variable_features(lfq, wide, top_n = n_all + 10)), n_all)
+  expect_equal(nrow(prolfqua:::.select_most_variable_features(lfq, wide, top_n = NULL)), n_all)
+  expect_equal(nrow(prolfqua:::.select_most_variable_features(lfq, wide, top_n = Inf)), n_all)
+
+  # top_n < feature count keeps exactly top_n rows, and they are the ones with
+  # the highest variability statistic (CV here — untransformed data).
+  k <- 3L
+  sub <- prolfqua:::.select_most_variable_features(lfq, wide, top_n = k)
+  expect_equal(nrow(sub), k)
+
+  st <- lfq$get_Stats(stats = "all")
+  keys <- c(lfq$hierarchy_keys(), lfq$isotope_label())
+  ranked <- dplyr::left_join(
+    wide$rowdata,
+    dplyr::select(st$stats(), dplyr::all_of(c(keys, st$stat))),
+    by = keys
+  )
+  metric <- stats::setNames(ranked[[st$stat]], rownames(wide$data))
+  kept <- metric[rownames(sub)]
+  dropped <- metric[setdiff(rownames(wide$data), rownames(sub))]
+  expect_gte(min(kept, na.rm = TRUE), max(dropped, na.rm = TRUE))
+
+  # end-to-end through the plotter: no hclust blow-up, returns a Heatmap.
+  p <- lfq$get_Plotter()$heatmap(top_n = k)
+  expect_s4_class(p, "Heatmap")
+})