Shrink demo data for Bioconductor + fix protein-name resolution + fix cross-model drug selection

R/globals.R

@@ -43,5 +43,9 @@ utils::globalVariables(c(
   "trained_year",
 
   # data.table operator
-  ":="
+  ":=",
+
+  # rlang pronoun used to disambiguate an argument from a same-named data
+  # column inside dplyr filters (see makeCrossModelPerformancePlot).
+  ".env"
 ))

R/plots_crossmodel.R

@@ -193,32 +193,52 @@ makeCrossModelRidgePlot <- function(perf_data, bug, cross_model) {
 #' @param bug Species code(s) to include.
 #' @param drug Drug/class identifier(s) to include.
 #' @param cross_model "country" or "time".
-#' @return A plotly heatmap, or NULL when there is nothing to plot.
+#' @return A plotly heatmap, or an empty placeholder plot (with an explanatory
+#'   title) when there is nothing to compare.
 #' @keywords internal
 #' @noRd
 makeCrossModelPerformancePlot <- function(perf_data, bug, drug, cross_model) {
+  placeholder <- function(msg) {
+    plotly::plot_ly() |>
+      plotly::layout(title = list(text = msg, x = 0))
+  }
+
   if (is.null(perf_data) || !is.data.frame(perf_data) || !nrow(perf_data)) {
-    return(NULL)
+    return(placeholder("No data available"))
   }
 
   strat <- if (cross_model == "country") "country" else "year"
 
+  # `.env$drug`: perf_data has a column literally named `drug`, which would
+  # otherwise shadow the `drug` argument inside the data mask and make this
+  # filter ignore the selected drug entirely.
   df <- perf_data |>
     dplyr::filter(normalize_species(.data$species) %in% normalize_species(bug)) |>
-    dplyr::filter(.data$drug_or_class %in% drug) |>
+    dplyr::filter(.data$drug_or_class %in% .env$drug) |>
     dplyr::filter(.data$strat_label == strat)
 
   if (!nrow(df)) {
-    return(NULL)
+    return(placeholder("No data for selection"))
+  }
+
+  # A cross-stratum grid needs at least two trained-on strata to compare.
+  n_strata <- dplyr::n_distinct(df$strat_value)
+  if (n_strata < 2) {
+    unit <- if (cross_model == "country") "countries" else "time periods"
+    return(placeholder(paste0(
+      "Not enough data to compare: needs models trained on 2+ ", unit,
+      " (this selection has ", n_strata, ")."
+    )))
   }
 
-  # For self-evaluation rows (strat_value_test is NA), set tested = trained
+  # For self-evaluation rows (strat_value_test is NA), set tested = trained.
+  # coalesce(as.character(...)) rather than if_else() so an all-NA logical
+  # strat_value_test column (no cross-tested rows) does not raise a type error.
   df <- df |>
     dplyr::mutate(
-      strat_value_test = dplyr::if_else(
-        is.na(.data$strat_value_test),
-        .data$strat_value,
-        .data$strat_value_test
+      strat_value_test = dplyr::coalesce(
+        as.character(.data$strat_value_test),
+        .data$strat_value
       )
     )
 
@@ -234,7 +254,7 @@ makeCrossModelPerformancePlot <- function(perf_data, bug, drug, cross_model) {
     as.matrix()
 
   if (!length(models_performance)) {
-    return(NULL)
+    return(placeholder("No data for selection"))
   }
 
   plotly::plot_ly(

R/utils_annotation.R

@@ -76,7 +76,13 @@ load_feature_name_map <- function(species_code, model_scale,
   }
   if (scale == "protein" &&
     all(c("proteinID", "proteinName") %in% names(df))) {
-    return(tibble::tibble(Variable = df$proteinID, label = df$proteinName))
+    # proteinID uses "fig|624..." but the top-features Variable it is matched
+    # against uses "fig.624...", so normalise "|" -> "." (as the annotated-join
+    # path in makeFeatureImportancePlot() also does) or the lookup never hits.
+    return(tibble::tibble(
+      Variable = stringr::str_replace(df$proteinID, "\\|", "."),
+      label = df$proteinName
+    ))
   }
   if (scale == "domain" && all(c("DB.ID", "SignDesc") %in% names(df))) {
     # domain ids need deduping since one Pfam can occur many times

inst/scripts/subset_demo_name_maps.R

@@ -0,0 +1,81 @@
+# Regenerate the bundled feature-name maps ({gene,protein,domain}_names.parquet)
+# by keeping only the ids that a top-features Variable can actually reference.
+# Run from the package source root, AFTER subset_demo_top_features.R.
+#
+# amRml ships whole-proteome/-genome name maps (e.g. Shigella sonnei's
+# `protein_names.parquet` is ~9.6M rows / ~10 MB). The dashboard only ever uses
+# them to relabel features that appear in a top-features plot: Feature
+# Importance looks up `name_map$Variable` for the features it is about to draw
+# and ignores the rest. So any id that never appears as a top-features
+# `Variable` is dead weight and is dropped.
+#
+# Id formats differ by scale and must be reconciled with the top-features
+# `Variable` column:
+#   * protein -- name map `proteinID` uses "fig|624...", top-features `Variable`
+#     uses "fig.624..."; we normalise "|" -> "." before matching.
+#   * domain  -- top-features `Variable` is "PF21279_IPR..."; the name map key
+#     is the part before the first "_" ("PF21279").
+#   * gene    -- ids match directly.
+# Matching on the normalised form (rather than the raw proteinID) keeps the map
+# correct even though the current runtime lookup does not itself normalise "|".
+
+suppressPackageStartupMessages({
+  library(arrow)
+  library(dplyr)
+  library(stringr)
+})
+
+# name-map file stem -> (id column, top-features feature_type it labels)
+specs <- list(
+  gene    = list(id = "Gene",      feature_type = "genes"),
+  protein = list(id = "proteinID", feature_type = "proteins"),
+  domain  = list(id = "DB.ID",     feature_type = "domains")
+)
+
+# Normalise a top-features Variable vector into the name-map's id space.
+to_key <- function(variable, scale) {
+  if (scale == "domain") {
+    stringr::str_split_i(variable, "_", 1)
+  } else {
+    variable
+  }
+}
+# Normalise a name-map id vector into the top-features Variable space.
+from_id <- function(id, scale) {
+  if (scale == "protein") gsub("[|]", ".", id) else id
+}
+
+species_dirs <- list.dirs("inst/extdata", recursive = FALSE)
+
+for (sp_dir in species_dirs) {
+  # Variables actually present in this species' trimmed top-features, by scale.
+  tf_files <- list.files(
+    sp_dir,
+    pattern = "^(all|country|year)_top_features[.]parquet$", full.names = TRUE
+  )
+  if (!length(tf_files)) next
+  tf <- dplyr::bind_rows(lapply(tf_files, function(f) {
+    arrow::read_parquet(f, col_select = c("Variable", "feature_type"))
+  }))
+
+  for (scale in names(specs)) {
+    id_col <- specs[[scale]]$id
+    ftype <- specs[[scale]]$feature_type
+    fp <- file.path(sp_dir, paste0(scale, "_names.parquet"))
+    if (!file.exists(fp)) next
+
+    keys <- unique(to_key(tf$Variable[tf$feature_type == ftype], scale))
+    orig_size <- file.size(fp)
+    full <- arrow::read_parquet(fp)
+    if (!id_col %in% names(full)) next
+
+    kept <- full[from_id(full[[id_col]], scale) %in% keys, , drop = FALSE]
+    arrow::write_parquet(kept, fp)
+
+    message(sprintf(
+      "%s: %d -> %d rows (%.2f -> %.2f MB)",
+      sub("inst/extdata/", "", fp), nrow(full), nrow(kept),
+      orig_size / 1e6, file.size(fp) / 1e6
+    ))
+  }
+}

inst/scripts/subset_demo_top_features.R

@@ -0,0 +1,94 @@
+# Regenerate the bundled top-features parquets by trimming them to only the rows
+# the dashboard can actually surface. Run from the package source root.
+#
+# amRml emits two flavours of `*_top_features.parquet` per species, both far
+# larger than a Bioconductor demo needs:
+#
+#   * baseline  `all_top_features.parquet`   -- has a `shuffled` column and a
+#     `seed` column; `strat_label` is NA (no country/year stratification).
+#   * stratified `country`/`year_top_features.parquet` -- one row per
+#     resampling replicate (encoded in `filename`/`seed_from_name`), with
+#     `strat_label`/`strat_value` naming the trained-on stratum.
+#
+# Two facts let us shrink both losslessly from the dashboard's point of view:
+#
+#   1. `shuffled` rows are the permutation-null baseline. No amRviz code path
+#      reads the `shuffled` column, so shuffled == TRUE rows can never appear in
+#      a plot and are dropped. (The `filename` and `seed_from_name` columns on
+#      the stratified files are likewise read nowhere, so they are dropped too.)
+#
+#   2. Every view collapses replicates with max/mean(Importance) per feature and
+#      then keeps at most the top N. The Feature Importance and Cross-model
+#      sliders both cap N at 100, so no view can ask for more than the top 100
+#      features of a model. A "model" at the finest user-selectable granularity
+#      is species x drug_or_class x feature_type (scale) x feature_subtype
+#      (encoding), plus strat_label x strat_value for the stratified files.
+#
+# We rank features by their max-across-replicate importance within each model
+# and keep every replicate row of the top 100. Keeping all replicate rows (not
+# an aggregated summary) means the app's own max()/mean(Importance) still return
+# the exact values they would on the full file -- the trim is invisible to the
+# dashboard. Ranking per (scale, subtype[, stratum]) independently makes the
+# kept set a superset of what any single view slices (Feature Importance folds
+# `struct` in alongside the chosen scale; the Network and Cross-model views fold
+# scales/subtypes together), so no view loses a feature it would have shown.
+#
+# Abhirupa (upstream data producer) has confirmed the `shuffled` rows are not
+# intended as a dashboard input.
+
+suppressPackageStartupMessages({
+  library(arrow)
+  library(dplyr)
+})
+
+top_n <- 100
+# Columns present on the stratified files that no amRviz code path reads.
+drop_cols <- c("filename", "seed_from_name")
+
+# All baseline + stratified top-features files; MDR/LOO/cross are loaded
+# differently (or not at all) and are left untouched.
+tf_files <- list.files(
+  "inst/extdata",
+  pattern = "^(all|country|year)_top_features[.]parquet$",
+  recursive = TRUE, full.names = TRUE
+)
+
+for (tf_file in tf_files) {
+  orig_size <- file.size(tf_file)
+  full <- arrow::read_parquet(tf_file)
+
+  # Drop the permutation-null rows and the unread replicate-bookkeeping columns.
+  kept <- full
+  if ("shuffled" %in% names(kept)) {
+    kept <- dplyr::filter(kept, is.na(.data$shuffled) | !.data$shuffled)
+  }
+  kept <- dplyr::select(kept, -dplyr::any_of(drop_cols))
+
+  # Model granularity: add the stratification columns when present.
+  model_cols <- c("species", "drug_or_class", "feature_type", "feature_subtype")
+  strat_cols <- intersect(c("strat_label", "strat_value"), names(kept))
+  model_cols <- c(model_cols, strat_cols)
+
+  # Rank features by max-across-replicate importance within each model, then
+  # keep every row belonging to the top `top_n` features of that model.
+  top_features <- kept |>
+    dplyr::group_by(dplyr::across(dplyr::all_of(c(model_cols, "Variable")))) |>
+    dplyr::summarise(
+      .maxImp = max(.data$Importance, na.rm = TRUE), .groups = "drop"
+    ) |>
+    dplyr::group_by(dplyr::across(dplyr::all_of(model_cols))) |>
+    dplyr::slice_max(order_by = .data$.maxImp, n = top_n, with_ties = FALSE) |>
+    dplyr::ungroup() |>
+    dplyr::select(dplyr::all_of(c(model_cols, "Variable")))
+
+  kept <- kept |>
+    dplyr::semi_join(top_features, by = c(model_cols, "Variable"))
+
+  arrow::write_parquet(kept, tf_file)
+
+  message(sprintf(
+    "%s: %d -> %d rows (%.2f -> %.2f MB)",
+    sub("inst/extdata/", "", tf_file), nrow(full), nrow(kept),
+    orig_size / 1e6, file.size(tf_file) / 1e6
+  ))
+}

tests/testthat/test-plots.R

@@ -126,18 +126,18 @@ test_that("makeFeatureImportancePlot returns NULL for empty data", {
 
 # ── makeCrossModelPerformancePlot ────────────────────────────────────────────
 
-test_that("makeCrossModelPerformancePlot returns NULL for empty data", {
+test_that("makeCrossModelPerformancePlot returns placeholder for empty data", {
   result <- makeCrossModelPerformancePlot(
     perf_data = NULL, bug = "Sau",
     drug = "ampicillin", cross_model = "country"
   )
-  expect_null(result)
+  expect_s3_class(result, "plotly")
 })
 
-test_that("makeCrossModelPerformancePlot returns NULL for zero-row data", {
+test_that("makeCrossModelPerformancePlot returns placeholder for zero-row data", {
   result <- makeCrossModelPerformancePlot(
     perf_data = data.frame(), bug = "Sau",
     drug = "ampicillin", cross_model = "country"
   )
-  expect_null(result)
+  expect_s3_class(result, "plotly")
 })