diff --git a/DESCRIPTION b/DESCRIPTION
index 3f59bac..1e2e237 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -51,9 +51,11 @@ Imports:
     glue,
     graphics,
     grDevices,
-    jsonlite,
     hardhat,
+    igraph,
+    jsonlite,
     methods,
+    networkD3,
     parsnip,
     purrr,
     patchwork,
diff --git a/NAMESPACE b/NAMESPACE
index 19d42ab..1107888 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -2,6 +2,7 @@
 
 export(.getTargetVarName)
 export(applyBenjaminiHochberg)
+export(buildFeatureNetwork)
 export(buildLRModel)
 export(buildPerfPq)
 export(buildPerfPqCrossCountry)
@@ -29,6 +30,8 @@ export(createMLResultDir)
 export(createMLinputList)
 export(encodePhenotype)
 export(extractTopFeats)
+export(findSharedDyads)
+export(findUniqueDyads)
 export(fitBestModel)
 export(generateMLInputs)
 export(getConfusionMatrix)
@@ -41,6 +44,7 @@ export(plotCM)
 export(plotCrossDrug)
 export(plotDrugDist)
 export(plotDrugPerf)
+export(plotFeatureNetworkD3)
 export(plotFishers)
 export(plotMDR)
 export(plotPRC)
@@ -61,6 +65,7 @@ export(runModelingPipeline)
 export(selectBestModel)
 export(shuffleLabels)
 export(splitMLInputTibble)
+export(topFeaturesPerDrugOrClass)
 export(tuneGrid)
 import(arrow)
 import(dplyr)
diff --git a/R/data.R b/R/data.R
index 07cbed9..9d6cba8 100644
--- a/R/data.R
+++ b/R/data.R
@@ -1,24 +1,38 @@
-#' Demo ML input tibble
-#'
-#' Stratified subset (30 Resistant + 30 Susceptible) of the AMP-genes-binary
-#' matrix from the bundled `Sfl_parquet.duckdb`, restricted to 80 feature
-#' columns.
-#'
-#' @format A tibble with 60 rows and 82 columns: `genome_id`,
-#'   `genome_drug.resistant_phenotype`, and 80 binary feature columns.
-#' @source `inst/scripts/make_demo_data.R`.
-#' @examples
-#' data(demo_ml_tibble)
-#' dim(demo_ml_tibble)
-"demo_ml_tibble"
+# FLAG for Abhirupa: data/demo_fit.rda and data/demo_ml_tibble.rda were
+# deleted in b8f45b1 ("Add example data and remove old, unused data"), but
+# the roxygen docs below (and @examples in plot_ml.R, run_ml_pipeline.R,
+# ife_ml.R, core_ml.R that call data(demo_ml_tibble)) still reference them.
+# Confirmed via `exists("demo_ml_tibble")` after devtools::load_all() -> FALSE,
+# the object genuinely isn't in the package anymore, so devtools::document()
+# fatally errors trying to build .Rd files for data that doesn't exist ('demo_ml_tibble'
+# is not an exported object from 'namespace:amRml'). Commented out both blocks
+# below as the minimal unblock so NAMESPACE regeneration can proceed for the
+# rest of the package. The real question -- restore the .rda files via
+# inst/scripts/make_demo_data.R, or retire these two blocks + the four
+# @examples blocks elsewhere that still call data(demo_ml_tibble) -- is still
+# open and needs her call; nothing below was deleted, just disabled.
 
-#' Demo LR fit
-#'
-#' A tuned logistic-regression workflow fitted on [demo_ml_tibble].
-#'
-#' @format A fitted `workflow` object (output of [fitBestModel()]).
-#' @source `inst/scripts/make_demo_data.R`.
-#' @examples
-#' data(demo_fit)
-#' class(demo_fit)
-"demo_fit"
+# #' Demo ML input tibble
+# #'
+# #' Stratified subset (30 Resistant + 30 Susceptible) of the AMP-genes-binary
+# #' matrix from the bundled `Sfl_parquet.duckdb`, restricted to 80 feature
+# #' columns.
+# #'
+# #' @format A tibble with 60 rows and 82 columns: `genome_id`,
+# #'   `genome_drug.resistant_phenotype`, and 80 binary feature columns.
+# #' @source `inst/scripts/make_demo_data.R`.
+# #' @examples
+# #' data(demo_ml_tibble)
+# #' dim(demo_ml_tibble)
+# "demo_ml_tibble"
+#
+# #' Demo LR fit
+# #'
+# #' A tuned logistic-regression workflow fitted on `demo_ml_tibble`.
+# #'
+# #' @format A fitted `workflow` object (output of [fitBestModel()]).
+# #' @source `inst/scripts/make_demo_data.R`.
+# #' @examples
+# #' data(demo_fit)
+# #' class(demo_fit)
+# "demo_fit"
diff --git a/R/feature_rescoring.R b/R/feature_rescoring.R
index b83a699..19ae859 100644
--- a/R/feature_rescoring.R
+++ b/R/feature_rescoring.R
@@ -1,113 +1,137 @@
-computeFeatureImprovement <- function(
-  all_top_features_parquet,
-  cluster_feature_parquet
-) {
-  stopifnot(file.exists(all_top_features_parquet))
+#' Filter model fits to those passing quality thresholds
+#'
+#' @param all_performance_parquet The path to the 'all performance parquet' file
+#' @param MCC_threshold The minimum non-shuffled MCC required to keep a model fit (default is \code{NULL}, no filtering)
+#' @param compare_to_shuffled Logical indicating whether to require the non-shuffled MCC to exceed (outperform) the shuffled-label MCC for the same fit (default is \code{TRUE}). A fit with no shuffled counterpart (\code{shuffled_MCC} is \code{NA} after pivoting) always passes this check regardless.
+#'
+#' @returns a tibble of model fits (one row per species, drug_label, drug_or_class, seed, feature_type, feature_subtype, model, fit_penalty, fit_mixture) that pass the requested quality thresholds, with \code{nonshuffled_MCC}, \code{shuffled_MCC}, and \code{MCC_diff} columns added. This does not pick a single "best" fit per group; it filters out fits that fail the MCC/shuffled-comparison criteria, so multiple passing fits per group can remain.
+#'
+#' @keywords internal
+#' @examples
+#' filterOptimalModel(all_performance_parquet = "inst/extdata/all_perf.parquet")
+filterOptimalModel <- function(all_performance_parquet,
+  MCC_threshold = NULL,
+  compare_to_shuffled = TRUE)
+  {
+  stopifnot(file.exists(all_performance_parquet))
+  all_perf <- arrow::read_parquet(normalizePath(all_performance_parquet))
 
-  # read the cluster-feature mapping
-  stopifnot(file.exists(cluster_feature_parquet))
-  cluster_feature <- arrow::read_parquet(normalizePath(cluster_feature_parquet))
+  if (!all(c(TRUE, FALSE) %in% unique(all_perf$shuffled))) {
+    stop("The 'shuffled' column must contain both TRUE and FALSE values. Run runModelingPipelineIntense() ")
+  }
 
-  features_rescored <- arrow::read_parquet(normalizePath(all_top_features_parquet)) |>
-    dplyr::filter(!shuffled) |>
+  all_perf |>
     dplyr::select(
       species, drug_label, drug_or_class,
-      feature_type, feature_subtype, seed, Variable,
-      Importance, Sign
+      seed, feature_type, feature_subtype,
+      model, fit_penalty, fit_mixture,
+      shuffled, mcc
     ) |>
-    dplyr::mutate(
-      Variable = dplyr::case_when(
-        feature_type == "domains" ~ sub("_.+$", "", Variable),
-        feature_type == "proteins" ~ sub("fig.", "fig|", Variable, fixed = TRUE),
-        feature_type == "args" ~ sub(
-          "^X", "",
-          gsub("\\.NCBIFAM", "", Variable)
-        ),
-        TRUE ~ Variable
-      )
+    tidyr::pivot_wider(
+      names_from = shuffled,
+      values_from = mcc,
+      names_prefix = "shuffled_"
     ) |>
-    dplyr::group_by(drug_label, drug_or_class, feature_type, feature_subtype, seed) |>
-    dplyr::mutate(
-      rank = dplyr::dense_rank(dplyr::desc(Importance)),
-      denom = sum(Importance, na.rm = TRUE),
-      contribution = dplyr::if_else(denom > 0, Importance / denom, 0),
-
-      # Safer rescaling within each output_prefix
-      min_imp = min(Importance, na.rm = TRUE),
-      max_imp = max(Importance, na.rm = TRUE),
-      range_imp = max_imp - min_imp,
-      rescaled = dplyr::if_else(range_imp > 0, (Importance - min_imp) / range_imp, 0)
-    ) |>
-    dplyr::ungroup() |>
-    dplyr::group_by(drug_label, drug_or_class, feature_type, feature_subtype, Variable) |>
-    dplyr::mutate(median_datatype = stats::median(rank, na.rm = TRUE)) |>
-    dplyr::ungroup() |>
-    dplyr::group_by(drug_label, drug_or_class, feature_type, shuffled, Variable) |>
-    dplyr::mutate(median_scale = stats::median(median_datatype, na.rm = TRUE)) |>
-    dplyr::ungroup() |>
-    # NOTE: correct join syntax (no quotes in join_by)
-    dplyr::left_join(cluster_feature, by = dplyr::join_by(Variable == feature)) |>
-    dplyr::group_by(drug_label, drug_or_class, shuffled, cluster) |>
-    dplyr::mutate(
-      median_drug_or_class       = stats::median(rank, na.rm = TRUE),
-      count_scales_for_cluster   = dplyr::n_distinct(feature_type),
-      feature_types_csv          = paste(sort(unique(feature_type)), collapse = ","),
-      lowest_contri              = min(contribution, na.rm = TRUE),
-      highest_contri             = max(contribution, na.rm = TRUE)
+    dplyr::rename(
+      nonshuffled_MCC = shuffled_FALSE,
+      shuffled_MCC = shuffled_TRUE
     ) |>
-    dplyr::ungroup() |>
-    dplyr::group_by(drug_label, drug_or_class, cluster) |>
+    dplyr::filter(!is.na(nonshuffled_MCC)) |>
     dplyr::mutate(
-      median_cluster = stats::median(median_drug_or_class, na.rm = TRUE),
-      count_features_for_cluster = dplyr::n_distinct(Variable)
+      MCC_diff = nonshuffled_MCC - shuffled_MCC
     ) |>
-    dplyr::ungroup()
-
-  return(features_rescored)
+    dplyr::filter(
+      if (!is.null(MCC_threshold)) (nonshuffled_MCC >= MCC_threshold) else TRUE,
+      if (compare_to_shuffled) (MCC_diff > 0 | is.na(shuffled_MCC)) else TRUE
+    )
 }
 
-computeFeatureScore <- function(
-  all_top_features_parquet,
-  cluster_feature_parquet
-) {
-  top_features <- arrow::read_parquet(normalizePath(all_top_features_parquet)) |>
+#' Score features within each seed
+#'
+#' @param all_top_features_parquet The path to the Parquet file containing all top features with their importance scores.
+#' @param core_contribution_threshold The cumulative-contribution cutoff, in \[0, 1\] (default is \code{0.75}, i.e. 75%), used to flag whether a feature falls within the "core" set of features that jointly account for that share of a seed's total importance.
+#' @param exclude_feature_types Feature types to drop before any scoring happens (default is \code{NULL}, i.e. no feature types are excluded and struct is currently included). struct variables are composite IDs (e.g. \code{polA.group_211.group_2176}, three dot-joined gene/domain identifiers) representing a co-occurrence/structural motif rather than a single molecular entity like the other five scales, and its candidate-variable count (tens of thousands per group) dwarfs the other scales by orders of magnitude — pooling it into this rank_score/contribution machinery would compare a compound signal against five primary ones on an incomparable scale. Pass \code{"struct"} here to exclude it once struct is meant to be reserved for post-hoc biological annotation rather than scoring/ranking/thresholding.
+#' @param filter_model Logical indicating whether to restrict scoring to (species, drug_label, drug_or_class, feature_type, feature_subtype, seed) groups that have at least one model fit passing \code{filterOptimalModel()}'s MCC/shuffled-comparison quality thresholds (default is \code{TRUE}). The join is not keyed on \code{model}/\code{fit_penalty}/\code{fit_mixture}; this is only safe because there is currently never more than one fit per (species, drug_label, drug_or_class, feature_type, feature_subtype, seed) group. If that assumption changes, this join (and the \code{all_perf} join below) would need a finer-grained key to avoid silently keeping/duplicating rows from multiple fits.
+#' @param all_performance_parquet The path to the all performance parquet file. Always required — it is read unconditionally 
+#' @param MCC_threshold The minimum MCC threshold passed through to \code{filterOptimalModel()} (default is \code{NULL}, no filtering)
+#' @param compare_to_shuffled Logical indicating whether to compare the model to shuffled data, passed through to \code{filterOptimalModel()} (default is \code{TRUE})
+#'
+#' @returns a tibble of scored top features, one row per feature within each species/drug_label/drug_or_class/feature_type/feature_subtype/seed group, with the following columns added:
+#' \itemize{
+#'   \item \code{contribution}: the feature's importance divided by the sum of importance across all features in the group.
+#'   \item \code{feat_return_ratio}: \code{n_feats_returned / n_feat} for the fit that produced this group. This assumes there is never more than one surviving fit per (species, drug_label, drug_or_class, seed, feature_type, feature_subtype) group; if that ever stops being true, this join would need a finer-grained key (\code{all_top_features_parquet} carries no \code{model}/\code{fit_penalty}/\code{fit_mixture} column to join on directly) or the ratio would need to be aggregated across fits before joining, to avoid fanning out and double-counting feature rows.
+#'   \item \code{rank}: descending rank of \code{contribution} within the group; ties receive the average of the ranks they span.
+#'   \item \code{n_features}: the number of rows (features) in the group.
+#'   \item \code{rank_score}: \code{(n_features - rank) / (n_features - 1)}; ranges 0-1 with higher values indicating higher importance (a single-feature group scores 1).
+#'   \item \code{cum_contrib}: the cumulative sum of \code{contribution} in descending order; features tied on \code{contribution} share the same \code{cum_contrib}, equal to the cumulative sum through the end of their tied block, so a tie is never split by arbitrary sort order.
+#'   \item \code{in_core}: TRUE when \code{cum_contrib <= core_contribution_threshold}; a tied block that would push the cumulative total past the threshold is excluded in its entirety (conservative: stays at-or-under the threshold rather than overshooting it).
+#' }
+#'
+#' @keywords internal
+#' @examples
+#' scoreFeaturesWithinSeed(all_top_features_parquet = "inst/extdata/all_top_features.parquet",
+#' all_performance_parquet = "inst/extdata/all_perf.parquet")
+#'
+scoreFeaturesWithinSeed <- function(all_top_features_parquet,
+  core_contribution_threshold = 0.75,
+  exclude_feature_types = NULL,
+  filter_model = TRUE,
+  all_performance_parquet,
+  MCC_threshold = NULL,
+  compare_to_shuffled = TRUE)
+  {
+  # check for the all_perf.parquet and all_top_features.parquet files
+  stopifnot(file.exists(all_top_features_parquet))
+  stopifnot(file.exists(all_performance_parquet))
+
+  all_top_features <- arrow::read_parquet(normalizePath(all_top_features_parquet)) |>
+    dplyr::filter(!shuffled, !feature_type %in% exclude_feature_types)
+
+  if (filter_model) {
+    filtered_model <- filterOptimalModel(
+      normalizePath(all_performance_parquet),
+      MCC_threshold = MCC_threshold,
+      compare_to_shuffled = compare_to_shuffled
+    )
+
+    all_top_features <- all_top_features |>
+      dplyr::semi_join(
+        filtered_model,
+        by = dplyr::join_by(
+          species, drug_label, drug_or_class,
+          feature_type, feature_subtype, seed
+        )
+      )
+  }
+
+  all_perf <- arrow::read_parquet(normalizePath(all_performance_parquet)) |>
     dplyr::filter(!shuffled) |>
+    dplyr::select(species, drug_label, drug_or_class, seed,
+      feature_type, feature_subtype, fit_penalty, fit_mixture,
+      mcc, n_feat, n_feats_returned) |>
+    dplyr::mutate(
+      feat_return_ratio = n_feats_returned / n_feat
+    )
+
+  # add different layers of scoring to the features within each seed
+  scored_top_features <- all_top_features |>
     dplyr::select(
-      species, drug_label, drug_or_class,
-      feature_type, feature_subtype, seed, Variable,
-      Importance, Sign
+      species, drug_label, drug_or_class, seed,
+      feature_type, feature_subtype, variable = Variable,
+      importance = Importance, sign = Sign
     ) |>
     dplyr::mutate(
-      Variable = dplyr::case_when(
-        feature_type == "domains" ~ sub("_.+$", "", Variable),
-        feature_type == "proteins" ~ sub("fig.", "fig|", Variable, fixed = TRUE),
-        feature_type == "args" ~ sub(
+      variable = dplyr::case_when(
+        feature_type == "protein" ~ sub("fig.", "fig|", variable, fixed = TRUE),
+        feature_type == "AMRFinder" ~ sub(
           "^X", "",
-          gsub("\\.NCBIFAM", "", Variable)
+          gsub("\\.NCBIFAM", "", variable)
         ),
-        TRUE ~ Variable
+        TRUE ~ variable
       )
-    )
-  cluster_feature <- arrow::read_parquet(normalizePath(cluster_feature_parquet))
-
-  ## --------------------------------
-  ## 1. Scale capability (data-driven)
-  ## --------------------------------
-
-  scale_capability <- top_features |>
-    dplyr::distinct(feature_type, feature_subtype) |>
-    dplyr::group_by(feature_type) |>
-    dplyr::summarise(
-      expected_types = dplyr::n_distinct(feature_subtype),
-      expected_types_csv = paste(sort(feature_subtype), collapse = ","),
-      .groups = "drop"
-    )
-
-  ## --------------------------------
-  ## 2. Importance → contribution → rank
-  ## --------------------------------
-
-  ranked_features <- top_features |>
+    ) |>
+    dplyr::left_join(all_perf, by = dplyr::join_by(species, drug_label, drug_or_class, seed,
+      feature_type, feature_subtype)) |>
     dplyr::group_by(
       species,
       drug_label,
@@ -117,173 +141,1274 @@ computeFeatureScore <- function(
       seed
     ) |>
     dplyr::mutate(
-      contribution = Importance / sum(Importance, na.rm = TRUE),
-      rank = dplyr::dense_rank(dplyr::desc(contribution)),
+      contribution = importance / sum(importance, na.rm = TRUE)
+    ) |>
+    dplyr::arrange(dplyr::desc(contribution), .by_group = TRUE) |>
+    dplyr::mutate(
+      rank = rank(dplyr::desc(contribution), ties.method = "average"),
       n_features = dplyr::n(),
-      rank_score = ifelse(n_features > 1, (n_features - rank) / (n_features - 1), 1)
+      rank_score = dplyr::if_else(
+        n_features > 1,
+        (n_features - rank) / (n_features - 1),
+        1
+      ),
+      running_contrib = cumsum(contribution)
+    ) |>
+    dplyr::group_by(contribution, .add = TRUE) |>
+    dplyr::mutate(
+      cum_contrib = max(running_contrib),
+      in_core = cum_contrib <= core_contribution_threshold
     ) |>
     dplyr::ungroup() |>
+    dplyr::select(-running_contrib)
+
+  return(scored_top_features)
+}
+
+#' Summarize a feature's scoring across seeds
+#'
+#' Answers "which molecular features are consistently important?" by
+#' collapsing the per-seed rows from `scoreFeaturesWithinSeed()` down to one
+#' row per feature.
+#'
+#' @param scored_top_features The tibble of scored top features with their contribution, rank, and rank score within each seed generated from `scoreFeaturesWithinSeed()`
+#'
+#' @returns a tibble with one row per species/drug_label/drug_or_class/feature_type/feature_subtype/variable, with:
+#' \itemize{
+#'   \item \code{seed_ratio}: the number of seeds the feature appears in, divided by the total number of distinct \code{seed} values present anywhere in \code{scored_top_features} (a single count computed once for the whole call, not per feature or per group — so this assumes every group was fit with the same set of seeds).
+#'   \item \code{mean_rank_score}, \code{median_rank_score}: mean/median of \code{rank_score} across seeds; ranges 0-1 with higher values indicating higher importance.
+#'   \item \code{median_rank}: median of \code{rank} across seeds. (\code{mean_rank} is not currently computed.)
+#'   \item \code{median_contribution}: median of \code{contribution} across seeds. (\code{mean_contribution} is not currently computed.)
+#'   \item \code{median_cum_contrib}: median of \code{cum_contrib} across seeds.
+#'   \item \code{best_rank}: the best (lowest) rank seen across seeds.
+#'   \item \code{rank_consistent}: TRUE if the feature's rank is identical in every seed.
+#'   \item \code{rank_score_sd}, \code{rank_score_cv}: standard deviation and coefficient of variation of \code{rank_score} across seeds, i.e. how much the normalized rank_score varies — computed on \code{rank_score} rather than raw rank so it is comparable across groups with different numbers of features.
+#'   \item \code{in_core_consistent}, \code{in_core}: whether the feature's \code{in_core} flag is identical across every seed; \code{in_core} is that shared value if consistent, else \code{FALSE}.
+#'   \item \code{sign_consistent}, \code{sign}: whether the feature's \code{sign} is identical across every seed; \code{sign} is that shared value if consistent, else \code{"MIXED"}.
+#' }
+#'
+#' @keywords internal
+#' @examples
+#' summariseFeaturesAcrossSeeds(scoreFeaturesWithinSeed(all_top_features.parquet))
+summariseFeaturesAcrossSeeds <- function(scored_top_features) {
+
+  # find max number of seeds
+  max_seeds <- scored_top_features |>
+    dplyr::summarise(n_seeds = dplyr::n_distinct(seed)) |>
+    dplyr::pull(n_seeds)
+
+  feature_summary <- scored_top_features |>
     dplyr::group_by(
       species,
       drug_label,
       drug_or_class,
       feature_type,
       feature_subtype,
-      Variable
+      variable
     ) |>
     dplyr::summarise(
-      n_seeds = dplyr::n_distinct(seed),
-      mean_rank = mean(rank, na.rm = TRUE),
-      median_rank = median(rank, na.rm = TRUE),
+      seed_ratio = dplyr::n_distinct(seed) / max_seeds,
+
       mean_rank_score = mean(rank_score, na.rm = TRUE),
+      median_rank_score = median(rank_score, na.rm = TRUE),
+      rank_score_sd = sd(rank_score, na.rm = TRUE),
+      # coefficient of variation: how large rank_score_sd is relative to mean_rank_score
+      rank_score_cv = dplyr::if_else(
+        mean_rank_score != 0,
+        rank_score_sd / mean_rank_score,
+        NA_real_
+      ),
+
+      median_rank = median(rank, na.rm = TRUE),
       best_rank = min(rank, na.rm = TRUE),
       rank_consistent = dplyr::n_distinct(rank) == 1,
-      rank_sd = sd(rank),
-      rank_score_cv = sd(rank_score) / mean(rank_score),
-      sign_consistent = dplyr::n_distinct(Sign) == 1,
-      sign = if (sign_consistent) dplyr::first(Sign) else "mixed",
+
+      median_contribution = median(contribution, na.rm = TRUE),
+      median_cum_contrib = median(cum_contrib, na.rm = TRUE),
+
+      in_core_consistent = dplyr::n_distinct(in_core) == 1,
+      in_core = if (in_core_consistent) dplyr::first(in_core) else FALSE,
+
+      sign_consistent = dplyr::n_distinct(sign) == 1,
+      sign = if (sign_consistent) dplyr::first(sign) else "MIXED",
       .groups = "drop"
+    )
+
+  return(feature_summary)
+}
+
+#' Build a per-drug top-feature table with cutoffs
+#'
+#' Runs `scoreFeaturesWithinSeed()` and `summariseFeaturesAcrossSeeds()` on
+#' `all_top_features_parquet`, then filters the resulting per-feature summary
+#' down to the top features for each drug/class.
+#'
+#' @inheritParams scoreFeaturesWithinSeed
+#' @param rank_score_quantile A value in \[0, 1\] (default is \code{0.95}). Keep only features whose \code{median_rank_score} is at or above this quantile of \code{median_rank_score}, where the quantile is computed separately within each (species, drug_label, drug_or_class) group (with \code{na.rm = TRUE}) — so a feature is compared only against other features from the same drug/class, not against the whole panel. This means a drug/class with weaker or more diffuse signal overall can still contribute its own top features, rather than being crowded out by drugs/classes with stronger or more concentrated signal. As with the other conditions listed below, this is not restricted to rows that already pass them: \code{dplyr::filter()} evaluates every condition passed to a single call against the same original data, so this threshold does not narrow as other conditions are applied
+#' @param cv_threshold The maximum allowed coefficient of variation (default is \code{1}). Keep only features with \code{rank_score_cv <= cv_threshold}, i.e. drop features whose rank_score is inconsistent across seeds relative to its mean. A feature present in every seed (\code{seed_ratio == 1}) always passes this check regardless of its \code{rank_score_cv} — including when \code{rank_score_cv} is \code{NA}, which happens whenever a group has only a single seed
+#' @param cumulative_contribution_threshold The cumulative-contribution cutoff, in \[0, 1\] (default is \code{0.75}, i.e. 75%). Keep only features with \code{median_cum_contrib <= cumulative_contribution_threshold}
+#' @param seed_ratio_threshold If not \code{NULL} (the default), keep only features whose \code{seed_ratio} exactly equals this value
+#' @param found_in_both_subtypes Logical indicating whether to additionally restrict to features that survive the filters above in both the binary and counts \code{feature_subtype} (default is \code{FALSE})
+#' @param compare_median_to_sd_rank_score Logical indicating whether to additionally require \code{median_rank_score > rank_score_sd} (default is \code{FALSE})
+#'
+#' @returns a tibble of top features for each drug/class: the `summariseFeaturesAcrossSeeds()` output (species, drug label, drug or class, feature type, feature subtype, variable, seed_ratio, mean/median rank score, median rank, median contribution, median cumulative contribution, best rank, rank/sign/in_core consistency flags, sign), filtered in two \code{dplyr::filter()} passes.
+#' The first pass keeps rows where all of the following hold, evaluated together against \code{feature_summary} (see \code{rank_score_quantile} for how the last condition is grouped):
+#' \itemize{
+#'   \item \code{seed_ratio == seed_ratio_threshold}, only applied when \code{seed_ratio_threshold} is not \code{NULL},
+#'   \item \code{seed_ratio == 1} OR \code{rank_score_cv <= cv_threshold} (see \code{cv_threshold}),
+#'   \item \code{in_core} is TRUE,
+#'   \item \code{sign_consistent} is TRUE (sign is the same in every seed; this does not require the sign to be negative),
+#'   \item \code{median_cum_contrib <= cumulative_contribution_threshold}, and
+#'   \item \code{median_rank_score} is at or above the \code{rank_score_quantile} quantile of \code{median_rank_score}, computed separately per (species, drug_label, drug_or_class) group.
+#' }
+#' Two columns are then added, grouped by (species, drug_label, drug_or_class, feature_type, variable): \code{n_subtype} and \code{subtype_csv}, recording how many/which \code{feature_subtype} values each combination has among the rows that survived the first pass.
+#' A second \code{dplyr::filter()} pass then optionally keeps only rows where \code{subtype_csv == "binary,counts"} (when \code{found_in_both_subtypes = TRUE}) and/or \code{median_rank_score > rank_score_sd} (when \code{compare_median_to_sd_rank_score = TRUE}).
+#' Every drug/class may not have variables from all feature types. Because the \code{rank_score_quantile} cutoff is now per drug/class rather than global, every drug/class with at least one feature surviving the other conditions will generally contribute something here, rather than potentially being excluded entirely by comparison to stronger-signal drugs/classes elsewhere in the panel.
+#'
+#' @export
+topFeaturesPerDrugOrClass <- function(
+  all_top_features_parquet,
+  core_contribution_threshold = 0.75,
+  exclude_feature_types = NULL,
+  filter_model = TRUE,
+  all_performance_parquet,
+  MCC_threshold = NULL,
+  compare_to_shuffled = TRUE,
+  rank_score_quantile = 0.95,
+  cv_threshold = 1,
+  cumulative_contribution_threshold = 0.75,
+  # additional filters
+  seed_ratio_threshold = NULL,
+  found_in_both_subtypes = FALSE,
+  compare_median_to_sd_rank_score = FALSE
+  )
+  {
+
+  scored_features <- scoreFeaturesWithinSeed(
+    all_top_features_parquet,
+    core_contribution_threshold = core_contribution_threshold,
+    exclude_feature_types = exclude_feature_types,
+    filter_model = filter_model,
+    all_performance_parquet,
+    MCC_threshold = MCC_threshold,
+    compare_to_shuffled = compare_to_shuffled
+  )
+
+  feature_summary <- summariseFeaturesAcrossSeeds(scored_features)
+
+  top_filtered_features <- feature_summary |>
+    dplyr::group_by(species, drug_label, drug_or_class) |>
+    dplyr::mutate(
+      rank_score_cutoff = quantile(median_rank_score, rank_score_quantile, na.rm = TRUE)
+    ) |>
+    dplyr::ungroup() |>
+    dplyr::filter(
+      if (!is.null(seed_ratio_threshold)) seed_ratio == seed_ratio_threshold else TRUE,
+      seed_ratio == 1 | rank_score_cv <= cv_threshold, # rank_score_cv can be NA if there is only one seed.
+      in_core,
+      sign_consistent,
+      median_cum_contrib <= cumulative_contribution_threshold,
+      median_rank_score >= rank_score_cutoff
+    ) |>
+    dplyr::select(-rank_score_cutoff) |>
+    dplyr::group_by(species, drug_label, drug_or_class, feature_type, variable) |>
+    dplyr::mutate(
+      n_subtype = dplyr::n_distinct(feature_subtype),
+      subtype_csv = paste(sort(unique(feature_subtype)), collapse = ",")
     ) |>
-    dplyr::arrange(dplyr::desc(mean_rank_score), mean_rank, best_rank)
+    dplyr::ungroup() |>
+    dplyr::filter(
+      if (found_in_both_subtypes) (subtype_csv == "binary,counts") else TRUE,
+      if (compare_median_to_sd_rank_score) (median_rank_score > rank_score_sd) else TRUE
+    )
+
+  return(top_filtered_features)
+}
 
-  ## --------------------------------
-  ## 3. Collapse to protein level
-  ## --------------------------------
+#' Aggregate mapped features to protein dyads
+#'
+#' Internal helper called on the output of `topFeaturesPerDrugOrClass()`.
+#'
+#' @param top_filtered_features The tibble of top features for each drug/class generated from `topFeaturesPerDrugOrClass()`.
+#' @param dyad_feature_parquet The path to the Parquet file containing the mapping of features to protein dyads. Must have a \code{feature} column formatted as \code{"<feature_type>:<variable>"} (using the short feature-type codes \code{amr}/\code{cog}/\code{defense}/\code{pfam}/\code{protein}) and a \code{dyad} column giving the dyad id.
+#'
+#' @returns a tibble with one row per species/drug_label/drug_or_class/dyad (\code{dyad} is the protein-dyad id), with:
+#' \itemize{
+#'   \item \code{frequency}: the number of top-feature rows mapped to this dyad.
+#'   \item \code{n_variables}, \code{variables_csv}: number of, and comma-separated list of, distinct \code{feature} (\code{"<feature_type>:<variable>"}) values mapped to this dyad.
+#'   \item \code{n_feature_types}, \code{feature_types_csv}: number of, and comma-separated list of, distinct \code{feature_type_subtype} (\code{"<feature_type>:<feature_subtype>"}) values mapped to this dyad.
+#'   \item \code{dyad_median_rank_score}: median of \code{median_rank_score} across the top-feature rows mapped to this dyad.
+#'   \item \code{dyad_median_contribution}: median of \code{median_contribution} across the top-feature rows mapped to this dyad.
+#'   \item \code{sign_consistent}, \code{sign}: whether \code{sign} is identical across every top-feature row mapped to this dyad; \code{sign} is that shared value if consistent, else \code{"MIXED"}.
+#' }
+#'
+#' @keywords internal
+summariseDyads <- function(top_filtered_features,
+                            dyad_feature_parquet
+                            ) {
+  stopifnot(is.data.frame(top_filtered_features))
+  stopifnot(file.exists(dyad_feature_parquet))
 
-  feature_cluster <- arrow::read_parquet(normalizePath(feature_cluster_parquet))
+  dyad_feature <- arrow::read_parquet(normalizePath(dyad_feature_parquet)) |> 
+    dplyr::distinct()
 
-  ranked_features <- ranked_features |>
+  top_dyads <- top_filtered_features |>
     dplyr::mutate(
-      Variable = dplyr::case_when(
-        feature_type == "domains" ~ sub("_.+$", "", Variable),
-        feature_type == "proteins" ~ sub("fig.", "fig|", Variable, fixed = TRUE),
-        TRUE ~ Variable
+      feature_type = dplyr::case_when(
+        feature_type == "AMRFinder"  ~ "amr",
+        feature_type == "COG"        ~ "cog",
+        feature_type == "DefenseCas" ~ "defense",
+        feature_type == "Pfam"       ~ "pfam",
+        TRUE                         ~ feature_type
       )
     ) |>
-    dplyr::left_join(feature_cluster, by = dplyr::join_by(Variable == feature))
+    tidyr::unite("feature", feature_type, variable, sep = ":", remove = FALSE) |>
+    tidyr::unite("feature_type_subtype", feature_type, feature_subtype, sep = ":", remove = FALSE) |>
+    dplyr::left_join(dyad_feature, by = "feature", relationship = "many-to-many") |>
+    dplyr::filter(!is.na(dyad)) |>
+    dplyr::select(species, drug_label, drug_or_class, feature_type_subtype, dyad, feature, sign,
+      median_rank_score, median_contribution, subtype_csv) |>
+    dplyr::group_by(species, drug_label, drug_or_class, dyad) |>
+    dplyr::summarise(
+      frequency = dplyr::n(),
+      n_variables = dplyr::n_distinct(feature),
+      variables_csv = paste(sort(unique(feature)), collapse = ","),
+      n_feature_types = dplyr::n_distinct(feature_type_subtype),
+      feature_types_csv = paste(sort(unique(feature_type_subtype)), collapse = ","),
+      dyad_median_rank_score = median(median_rank_score, na.rm = TRUE),
+      dyad_median_contribution = median(median_contribution, na.rm = TRUE),
+      sign_consistent = dplyr::n_distinct(sign) == 1,
+      sign = if (sign_consistent) dplyr::first(sign) else "MIXED",
+      .groups = "drop"
+    ) |>
+    dplyr::arrange(dplyr::desc(dyad_median_rank_score), dplyr::desc(n_feature_types))
 
-  protein_scale_realization <- ranked_features |>
-    dplyr::filter(!is.na(cluster), shuffled == FALSE) |>
-    dplyr::group_by(
-      species,
-      drug_label,
-      drug_or_class,
-      cluster,
-      feature_type
+  return(top_dyads)
+}
+
+#' Build a feature network from selected top features and top dyads
+#'
+#' @param top_features Output of \code{topFeaturesPerDrugOrClass()}.
+#' @param top_dyads Output of \code{summariseDyads()}.
+#' @param dyad_feature_parquet Path to the same Parquet file passed to \code{summariseDyads()}: a \code{feature} column formatted as \code{"<feature_type>:<variable>"} (using the short feature-type codes \code{amr}/\code{cog}/\code{defense}/\code{pfam}/\code{protein}) and a \code{dyad} column giving the dyad id. \code{feature} is reconstructed here from \code{top_features$feature_type}/\code{variable} the same way \code{summariseDyads()} builds it, so the feature-to-dyad edges use the same mapping as the dyad table itself.
+#'
+#' @returns A list with \code{feature_table}, \code{dyad_table}, \code{nodes}, \code{edges}, and \code{graph}.
+#' Node/edge weights are built from \code{median_rank_score} (features) and \code{dyad_median_rank_score} (dyads, via the \code{dyad_score} column, itself median-based per \code{summariseDyads()}), consistent with the seed-noise-robust selection made in \code{topFeaturesPerDrugOrClass()}.
+#' @example \dontrun{buildFeatureNetwork(top_features = topFeaturesPerDrugOrClass(...), 
+#' top_dyads = summariseDyads(...), 
+#' dyad_feature_parquet = "inst/extdata/dyad_feature.parquet",
+#' filtered_model = filterOptimalModel(...))}
+#' @export
+buildFeatureNetwork <- function(top_features,
+median_rank_score_quantile = 1,
+                               top_dyads,
+                               n_feature_types_threshold = 6,
+                               dyad_feature_parquet,
+                               filtered_model,
+                               parquet_dir_path = NULL
+                              ) {
+  stopifnot(is.data.frame(top_features))
+  stopifnot(is.data.frame(top_dyads))
+  stopifnot(file.exists(dyad_feature_parquet))
+
+  required_feature_cols <- c(
+    "species", "drug_label", "drug_or_class",
+    "feature_type", "variable",
+    "median_rank_score"
+  )
+  required_dyad_cols <- c(
+    "species", "drug_label", "drug_or_class",
+    "dyad", "dyad_median_rank_score"
+  )
+
+  missing_feature_cols <- setdiff(required_feature_cols, names(top_features))
+  missing_dyad_cols <- setdiff(required_dyad_cols, names(top_dyads))
+
+  if (length(missing_feature_cols) > 0) {
+    stop("top_features is missing required columns: ",
+         paste(missing_feature_cols, collapse = ", "))
+  }
+  if (length(missing_dyad_cols) > 0) {
+    stop("top_dyads is missing required columns: ",
+         paste(missing_dyad_cols, collapse = ", "))
+  }
+
+  dyad_feature <- arrow::read_parquet(normalizePath(dyad_feature_parquet)) |>
+    dplyr::distinct() |>
+    dplyr::filter(!is.na(dyad), !is.na(feature))
+
+  make_model_id <- function(drug_label, drug_or_class) {
+    paste(drug_label, drug_or_class, sep = ".")
+  }
+
+  # Same short feature-type codes used by summariseDyads() to build `feature`.
+  shorten_feature_type <- function(feature_type) {
+    dplyr::case_when(
+      feature_type == "AMRFinder"  ~ "amr",
+      feature_type == "COG"        ~ "cog",
+      feature_type == "DefenseCas" ~ "defense",
+      feature_type == "Pfam"       ~ "pfam",
+      TRUE                         ~ feature_type
+    )
+  }
+
+  feature_table <- top_features |>
+    dplyr::filter(median_rank_score >= quantile(median_rank_score, median_rank_score_quantile, na.rm = TRUE)) |>
+    dplyr::mutate(feature_type = shorten_feature_type(feature_type)) |>
+    dplyr::mutate(
+      feature = paste(feature_type, variable, sep = ":"),
+      model_id = make_model_id(drug_label, drug_or_class)
     ) |>
+    dplyr::group_by(species, model_id, feature_type, feature) |>
     dplyr::summarise(
-      observed_types = dplyr::n_distinct(feature_subtype),
-      observed_types_csv = paste(sort(unique(feature_subtype)), collapse = ","),
+      # mean of median_rank_score across subtype (bin/count) rows for this variable --
+      # this is where bin/count reconciliation currently happens (implicitly)
+      feature_score = mean(median_rank_score, na.rm = TRUE),
+      .groups = "drop"
+    )
+
+ feature_dyad_edges <- dyad_feature |>
+  dplyr::semi_join(
+    feature_table,
+    by = "feature"
+  ) 
+
+  feature_dyad_edges <- feature_dyad_edges |>
+  dplyr::semi_join(
+    top_dyads |> dplyr::filter(n_feature_types > n_feature_types_threshold),
+    by = "dyad"
+  )
+  
+  model_feature_edges <- feature_table |>
+  dplyr::transmute(
+    source = model_id,
+    target = feature,
+    weight = feature_score,
+    edge_type = "model_feature"
+  ) |>
+  dplyr::distinct()
+  
+  feature_dyad_edges <- feature_dyad_edges |>
+  dplyr::transmute(
+    source = feature,
+    target = dyad,
+    weight = 1,
+    edge_type = "feature_dyad"
+  ) |>
+  dplyr::distinct()
+  
+  edges <- dplyr::bind_rows(
+  model_feature_edges,
+  feature_dyad_edges
+)
+  
+ feature_nodes <- feature_table |>
+  # dplyr::filter(feature %in% feature_dyad_edges$source) |>
+  dplyr::mutate(
+    node_type = feature_type
+  )|> 
+  dplyr::group_by(node_type, feature) |>
+  dplyr::summarise(
+    score = mean(feature_score),
+    breadth = dplyr::n_distinct(model_id),
+    .groups = "drop"
+  ) |>
+   dplyr::rename(name = feature)
+  
+  dyad_nodes <- top_dyads |>
+    dplyr::filter(dyad %in% feature_dyad_edges$target) |>
+  dplyr::group_by(dyad) |>
+  dplyr::summarise(
+    score = mean(dyad_median_rank_score),
+    breadth = dplyr::n_distinct(
+      paste(drug_label, drug_or_class)
+    ),
+    .groups = "drop"
+  ) |>
+  dplyr::mutate(
+    node_type = "dyad"
+  ) |>
+  dplyr::rename(name = dyad)
+
+  model_nodes <- filtered_model |>
+dplyr::group_by(drug_label, drug_or_class) |>
+    dplyr::summarize(
+      score = median(nonshuffled_MCC, na.rm = TRUE),
+      breadth = dplyr::n(),
       .groups = "drop"
     ) |>
-    dplyr::left_join(scale_capability, by = "feature_type") |>
+  tidyr::unite("name", drug_label, drug_or_class, sep = ".") |>
+  dplyr::mutate(
+    node_type = "model")
+
+  nodes <- bind_rows(
+  model_nodes,
+  feature_nodes,
+  dyad_nodes
+)
+  cog_names <- arrow::read_parquet(file.path(parquet_dir_path, "protein_COG.parquet")) |> 
+    dplyr::distinct(name = query_name, label = description)
+  
+defense_names <- arrow::read_parquet(file.path(parquet_dir_path, "protein_DefenseCas.parquet")) |> 
+    dplyr::distinct(name = query_name, label = description)
+
+ amr_names <- arrow::read_parquet(file.path(parquet_dir_path, "protein_AMRFinder.parquet")) |> 
+    dplyr::distinct(name = query_name, label = description) 
+
+ pfam_names <- arrow::read_parquet(file.path(parquet_dir_path, "protein_Pfam.parquet")) |> 
+    dplyr::distinct(name = query_name, label = description)  
+
+   protein_names <- arrow::read_parquet(file.path(parquet_dir_path, "protein_names.parquet")) |>
+    dplyr::distinct(name = proteinID, label = proteinName)
+
+  gene_names <- arrow::read_parquet(file.path(parquet_dir_path, "gene_names.parquet")) |>
+    dplyr::distinct(name = Gene, label = Annotation)
+
+  nodes <- nodes |>
+    dplyr::left_join(cog_names, by = "name") |>
+    dplyr::left_join(defense_names, by = "name") |>
+    dplyr::left_join(amr_names, by = "name") |>
+    dplyr::left_join(pfam_names, by = "name") |>
+    dplyr::left_join(protein_names, by = "name") |>
+    dplyr::left_join(gene_names, by = "name") |>
+    dplyr::mutate(
+      label = dplyr::coalesce(label.x, label.y, label.x.x, label.y.y, label.x.x.x, label.y.y.y)
+    ) |>
+    dplyr::select(-dplyr::starts_with("label.")) |>
     dplyr::mutate(
-      scale_realization = dplyr::case_when(
-        observed_types == expected_types ~ "full_realization",
-        observed_types < expected_types ~ "partial_realization"
+      label = dplyr::if_else(is.na(label), name, label)
+    )
+  
+  graph <- igraph::graph_from_data_frame(
+  edges,
+  directed = FALSE,
+  vertices = nodes
+)
+
+  feature_network <- list(
+  feature_table = feature_table,
+  feature_dyad_edges = feature_dyad_edges,
+  model_feature_edges = model_feature_edges,
+  nodes = nodes,
+  edges = edges,
+  graph = graph
+)
+
+  return(feature_network)
+}
+
+#' Plot the feature network with networkD3
+#'
+#' @param feature_network Output of \code{buildFeatureNetwork()}.
+#' @param height Widget height in pixels (default is \code{800}).
+#' @param width Widget width (default is \code{"100\%"}).
+#'
+#' @returns A \code{networkD3} widget.
+#' @export
+plotFeatureNetworkD3 <- function(feature_network,
+                                height = 800,
+                                width = "100%"
+                              ) {
+
+  stopifnot(is.list(feature_network))
+  stopifnot(!is.null(feature_network$nodes))
+  stopifnot(!is.null(feature_network$edges))
+
+  nodes <- feature_network$nodes |>
+    dplyr::distinct(name, .keep_all = TRUE) |>
+    dplyr::mutate(
+      id = dplyr::row_number() - 1L,
+      group = node_type,
+      title = paste0(
+        "<b>", name, "</b>",
+        ifelse(is.na(score), "", paste0("<br>Score: ", signif(score, 2))),
+        ifelse(is.na(breadth), "", paste0("<br>Breadth: ", breadth))
       )
     )
 
-  protein_scale_summary <- protein_scale_realization |>
-    dplyr::group_by(species, drug_label, drug_or_class, cluster) |>
+  links <- feature_network$edges |>
+    dplyr::filter(!is.na(source), !is.na(target))|>
+  dplyr::left_join(
+    nodes |>
+      dplyr::select(source = name, source_id = id),
+    by = "source"
+  ) |>
+  dplyr::left_join(
+    nodes |>
+      dplyr::select(target = name, target_id = id),
+    by = "target"
+  )
+
+  stopifnot(nrow(nodes) > 0)
+  stopifnot(nrow(links) > 0)
+
+  colour_scale <- networkD3::JS(
+"
+d3.scaleOrdinal()
+  .domain([
+    'model',
+    'amr',
+    'cog',
+    'defense',
+    'gene',
+    'pfam',
+    'protein',
+    'struct',
+    'dyad'
+  ])
+  .range([
+    '#4a6b8a',
+    '#87ceeb',
+    '#e6ab80',
+    '#8b6b7a',
+    '#4e9a9a',
+    '#c4a35a',
+    '#9b7fba',
+    '#5b8db8',
+    '#d4735e'
+  ])
+"
+)
+  networkD3::forceNetwork(
+    Links = links,
+    Nodes = nodes,
+    Source = "source_id",
+    Target = "target_id",
+    Value = "weight",
+    NodeID = "name",
+    Group = "group",
+    opacity = 0.9,
+    zoom = TRUE,
+    fontSize = 14,
+    height = height,
+    width = width,
+    colourScale = colour_scale,
+    linkDistance = networkD3::JS(
+      "function(d) {
+         if (d.edge_type === 'feature_dyad') return 60;
+         if (d.edge_type === 'model_feature') return 120;
+         return 90;
+       }"
+    )
+  )
+}
+
+# final run would be:
+# top_features <- topFeaturesPerDrugOrClass(rank_score_quantile = 0.75)
+# top_dyads <- summarisedyads(top_features, dyad_feature_parquet = dyad_feature_parquet)
+# feature_network <- buildFeatureNetwork(top_features = top_features, top_dyads = top_dyads,
+#   dyad_feature_parquet = dyad_feature_parquet)
+  # plotFeatureNetworkD3(feature_network)
+
+#' Discover high-performing models, stable features, and protein dyads
+#'
+#' Runs the model-quality, feature-ranking, and protein-dyad aggregation
+#' workflow and organizes the results into question-oriented tables.
+#'
+#' The returned object can be used to identify:
+#' \itemize{
+#'   \item model fits with good non-shuffled performance;
+#'   \item model fits with good separation from shuffled-label models;
+#'   \item model fits satisfying both performance conditions;
+#'   \item top features with high seed coverage and median rank score;
+#'   \item top features with high median contribution;
+#'   \item top features present across every analyzed seed;
+#'   \item top features shared across drugs or drug classes;
+#'   \item top features unique to one drug or drug class;
+#'   \item high-scoring, shared, and unique protein dyads.
+#' }
+#'
+#' This is a high-level wrapper around \code{filterOptimalModel()},
+#' \code{topFeaturesPerDrugOrClass()}, and \code{summariseDyads()}.
+#'
+#' @param all_top_features_parquet Path to the Parquet file containing
+#' top features and their importance scores.
+#' @param all_performance_parquet Path to the Parquet file containing
+#' model-performance results.
+#' @param dyad_feature_parquet Path to the Parquet file mapping features
+#' to protein dyads.
+#'
+#' @param MCC_threshold Minimum non-shuffled MCC used to define good model
+#' performance. Default is \code{0.5}, meaning that only models with a non-shuffled MCC of 0.5 or higher are considered.
+#' @param compare_to_shuffled Logical indicating whether the models used for
+#' feature selection must outperform their shuffled-label counterparts.
+#' Default is \code{TRUE}.
+#'
+#' @param core_contribution_threshold Cumulative-contribution threshold
+#' passed to \code{scoreFeaturesWithinSeed()}. Default is \code{0.9}.
+#' @param exclude_feature_types Feature types to exclude before feature
+#' scoring. Default is \code{NULL}.
+#' @param filter_model Logical indicating whether feature scoring should be
+#' restricted to model groups passing the requested model-quality criteria.
+#' Default is \code{TRUE}.
+#' @param rank_score_quantile Quantile cutoff for median rank score, computed
+#' separately within each species/drug-label/drug-or-class group.
+#' Default is \code{0.95}.
+#' @param cv_threshold Maximum allowed rank-score coefficient of variation.
+#' Default is \code{1}.
+#' @param cumulative_contribution_threshold Maximum allowed median cumulative
+#' contribution. Default is \code{0.75}.
+#' @param seed_ratio_threshold Optional exact seed-ratio requirement passed
+#' to \code{topFeaturesPerDrugOrClass()}. Default is \code{1}.
+#' @param found_in_both_subtypes Logical indicating whether selected features
+#' must occur in both binary and counts feature subtypes.
+#' Default is \code{FALSE}.
+#' @param compare_median_to_sd_rank_score Logical indicating whether selected
+#' features must have \code{median_rank_score > rank_score_sd}.
+#' Default is \code{FALSE}.
+#'
+#' @param consistent_seed_ratio Minimum seed ratio used to define a
+#' consistently selected feature. Default is \code{0.8}.
+#' @param consistent_rank_score Minimum median rank score used to define a
+#' consistently highly ranked feature. Default is \code{0.75}.
+#' @param high_contribution_quantile Quantile of \code{median_contribution}
+#' used to define high-contribution features. It is computed separately
+#' within each species/drug-label/drug-or-class group. Default is
+#' \code{0.75}.
+#' @param high_dyad_score_quantile Quantile of
+#' \code{dyad_median_rank_score} used to define high-scoring dyads. It is
+#' computed separately within each species/drug-label/drug-or-class group.
+#' Default is \code{0.75}.
+#' @param min_models_shared Minimum number of distinct drugs or drug classes
+#' required for a feature or dyad to be considered shared.
+#' Default is \code{2}.
+#'
+#' @returns An object of class \code{feature_dyad_discovery}, implemented as
+#' a named list containing:
+#' \itemize{
+#'   \item \code{models$good_performance}: fit-level models satisfying
+#'   \code{MCC_threshold}, without requiring shuffled-model separation.
+#'   \item \code{models$good_shuffled_separation}: fit-level models whose
+#'   non-shuffled MCC exceeds the shuffled-label MCC.
+#'   \item \code{models$qualified}: fit-level models satisfying the model
+#'   criteria used for feature selection.
+#'   \item \code{models$summary}: model-level summary across seeds and feature
+#'   scales.
+#'   \item \code{features$top}: all selected top features.
+#'   \item \code{features$consistent}: selected features with high seed ratio
+#'   and high median rank score.
+#'   \item \code{features$high_contribution}: selected features with high
+#'   median contribution.
+#'   \item \code{features$all_seeds}: selected features present across all
+#'   analyzed seeds.
+#'   \item \code{features$shared}: selected features associated with at least
+#'   \code{min_models_shared} drugs or drug classes.
+#'   \item \code{features$unique}: selected features associated with exactly
+#'   one drug or drug class.
+#'   \item \code{dyads$top}: dyads mapped from the selected top features.
+#'   \item \code{dyads$high_score}: dyads with high median rank scores.
+#'   \item \code{dyads$shared}: dyads associated with at least
+#'   \code{min_models_shared} drugs or drug classes.
+#'   \item \code{dyads$unique}: dyads associated with exactly one drug or
+#'   drug class.
+#'   \item \code{parameters}: parameter values used for the analysis.
+#' }
+#'
+#' @export
+runFeatureDyadDiscovery <- function(
+    all_top_features_parquet,
+    all_performance_parquet,
+    dyad_feature_parquet,
+    MCC_threshold = 0.5,
+    compare_to_shuffled = TRUE,
+    core_contribution_threshold = 0.9,
+    exclude_feature_types = NULL,
+    filter_model = TRUE,
+    rank_score_quantile = 0.95,
+    cv_threshold = 1,
+    cumulative_contribution_threshold = 0.75,
+    seed_ratio_threshold = 1,
+    found_in_both_subtypes = FALSE,
+    compare_median_to_sd_rank_score = FALSE,
+    consistent_seed_ratio = 0.8,
+    consistent_rank_score = 0.75,
+    high_contribution_quantile = 0.75,
+    high_dyad_score_quantile = 0.75,
+    min_models_shared = 2
+) {
+
+  # -------------------------------------------------------------------------
+  # Validate files
+  # -------------------------------------------------------------------------
+
+  input_files <- c(
+    all_top_features_parquet = all_top_features_parquet,
+    all_performance_parquet = all_performance_parquet,
+    dyad_feature_parquet = dyad_feature_parquet
+  )
+
+  missing_files <- input_files[!file.exists(input_files)]
+
+  if (length(missing_files) > 0) {
+    stop(
+      "The following input file(s) do not exist: ",
+      paste(names(missing_files), collapse = ", ")
+    )
+  }
+
+  # -------------------------------------------------------------------------
+  # Validate analysis thresholds
+  # -------------------------------------------------------------------------
+
+  unit_interval_parameters <- c(
+    core_contribution_threshold = core_contribution_threshold,
+    rank_score_quantile = rank_score_quantile,
+    cumulative_contribution_threshold =
+      cumulative_contribution_threshold,
+    consistent_seed_ratio = consistent_seed_ratio,
+    consistent_rank_score = consistent_rank_score,
+    high_contribution_quantile = high_contribution_quantile,
+    high_dyad_score_quantile = high_dyad_score_quantile
+  )
+
+  invalid_parameters <- names(unit_interval_parameters)[
+    is.na(unit_interval_parameters) |
+      unit_interval_parameters < 0 |
+      unit_interval_parameters > 1
+  ]
+
+  if (length(invalid_parameters) > 0) {
+    stop(
+      "The following parameter(s) must be between 0 and 1: ",
+      paste(invalid_parameters, collapse = ", ")
+    )
+  }
+
+  if (!is.null(seed_ratio_threshold)) {
+    if (
+      length(seed_ratio_threshold) != 1 ||
+      is.na(seed_ratio_threshold) ||
+      seed_ratio_threshold < 0 ||
+      seed_ratio_threshold > 1
+    ) {
+      stop(
+        "seed_ratio_threshold must be NULL or a single value ",
+        "between 0 and 1."
+      )
+    }
+  }
+
+  if (
+    length(cv_threshold) != 1 ||
+    is.na(cv_threshold) ||
+    cv_threshold < 0
+  ) {
+    stop("cv_threshold must be a single non-negative value.")
+  }
+
+  if (
+    length(min_models_shared) != 1 ||
+    is.na(min_models_shared) ||
+    min_models_shared < 2 ||
+    min_models_shared != as.integer(min_models_shared)
+  ) {
+    stop(
+      "min_models_shared must be a single integer greater than ",
+      "or equal to 2."
+    )
+  }
+
+  # -------------------------------------------------------------------------
+  # 1. Find fits with good absolute performance
+  #
+  # This table answers:
+  # "Which model fits have sufficiently high non-shuffled MCC?"
+  # -------------------------------------------------------------------------
+
+  good_performance_models <- filterOptimalModel(
+    all_performance_parquet = all_performance_parquet,
+    MCC_threshold = MCC_threshold,
+    compare_to_shuffled = FALSE
+  )
+
+  # -------------------------------------------------------------------------
+  # 2. Find fits with good separation from shuffled models
+  #
+  # This table answers:
+  # "Which model fits outperform their shuffled-label counterparts?"
+  #
+  # MCC_threshold is deliberately NULL here so that separation can be
+  # examined independently of absolute MCC.
+  # -------------------------------------------------------------------------
+
+  good_shuffled_separation_models <- filterOptimalModel(
+    all_performance_parquet = all_performance_parquet,
+    MCC_threshold = NULL,
+    compare_to_shuffled = TRUE
+  )
+
+  # -------------------------------------------------------------------------
+  # 3. Find models satisfying the criteria used for feature selection
+  #
+  # If compare_to_shuffled is TRUE, these satisfy both the requested
+  # MCC threshold and shuffled-label comparison.
+  # -------------------------------------------------------------------------
+
+  qualified_models <- filterOptimalModel(
+    all_performance_parquet = all_performance_parquet,
+    MCC_threshold = MCC_threshold,
+    compare_to_shuffled = compare_to_shuffled
+  )
+
+  if (nrow(qualified_models) == 0) {
+    stop(
+      "No model fits passed the requested model-quality criteria."
+    )
+  }
+
+  # Summarize fit-level results to the species/drug or drug-class level.
+  model_summary <- qualified_models |>
+    dplyr::group_by(
+      species,
+      drug_label,
+      drug_or_class
+    ) |>
     dplyr::summarise(
-      n_scales = dplyr::n_distinct(feature_type),
-      fully_realized_scales =
-        sum(scale_realization == "full_realization"),
-      partially_realized_scales =
-        sum(scale_realization == "partial_realization"),
-      scale_support_csv =
-        paste(
-          feature_type,
-          "(", observed_types_csv, "/", expected_types, ")",
-          collapse = "; "
-        ),
+      n_passing_fits = dplyr::n(),
+      n_seeds = dplyr::n_distinct(seed),
+      n_feature_types = dplyr::n_distinct(feature_type),
+      n_feature_scales = dplyr::n_distinct(
+        paste(feature_type, feature_subtype, sep = ":")
+      ),
+      median_nonshuffled_MCC = median(
+        nonshuffled_MCC,
+        na.rm = TRUE
+      ),
+      best_nonshuffled_MCC = max(
+        nonshuffled_MCC,
+        na.rm = TRUE
+      ),
+      median_shuffled_MCC = median(
+        shuffled_MCC,
+        na.rm = TRUE
+      ),
+      median_MCC_diff = median(
+        MCC_diff,
+        na.rm = TRUE
+      ),
+      if(!is.na(median_MCC_diff)) {
+        minimum_MCC_diff = min(
+          MCC_diff,
+          na.rm = TRUE
+        )
+      } else {
+        minimum_MCC_diff = NULL
+      },
+      if(!is.na(median_MCC_diff)) {
+        maximum_MCC_diff = max(
+          MCC_diff,
+          na.rm = TRUE
+        )
+      } else {
+        maximum_MCC_diff = NULL
+      },
+      proportion_outperforming_shuffled = median(
+        MCC_diff > 0 | is.na(shuffled_MCC),
+        na.rm = TRUE
+      ),
       .groups = "drop"
     ) |>
-    dplyr::mutate(
-      realization_score =
-        (fully_realized_scales +
-          0.5 * partially_realized_scales) /
-          (fully_realized_scales + partially_realized_scales)
+    dplyr::arrange(
+      dplyr::desc(median_nonshuffled_MCC),
+      dplyr::desc(median_MCC_diff)
+    )
+
+  # -------------------------------------------------------------------------
+  # 4. Select top features
+  # -------------------------------------------------------------------------
+
+  top_features <- topFeaturesPerDrugOrClass(
+    all_top_features_parquet = all_top_features_parquet,
+    core_contribution_threshold = core_contribution_threshold,
+    exclude_feature_types = exclude_feature_types,
+    filter_model = filter_model,
+    all_performance_parquet = all_performance_parquet,
+    MCC_threshold = MCC_threshold,
+    compare_to_shuffled = compare_to_shuffled,
+    rank_score_quantile = rank_score_quantile,
+    cv_threshold = cv_threshold,
+    cumulative_contribution_threshold =
+      cumulative_contribution_threshold,
+    seed_ratio_threshold = seed_ratio_threshold,
+    found_in_both_subtypes = found_in_both_subtypes,
+    compare_median_to_sd_rank_score =
+      compare_median_to_sd_rank_score
+  )
+
+  if (nrow(top_features) == 0) {
+    stop(
+      "No features passed the requested feature-selection criteria."
+    )
+  }
+
+  # -------------------------------------------------------------------------
+  # 5. Consistently top-ranked features
+  #
+  # Answers:
+  # "Which features have both high seed coverage and high median rank?"
+  # -------------------------------------------------------------------------
+
+  consistent_features <- top_features |>
+    dplyr::filter(
+      seed_ratio >= consistent_seed_ratio,
+      median_rank_score >= consistent_rank_score
     ) |>
-    dplyr::mutate(
-      coverage_boost =
-        n_scales / max(n_scales)
+    dplyr::arrange(
+      species,
+      drug_label,
+      drug_or_class,
+      dplyr::desc(seed_ratio),
+      dplyr::desc(median_rank_score)
+    )
+
+  # -------------------------------------------------------------------------
+  # 6. High-contribution features
+  #
+  # The contribution cutoff is calculated separately for each species and
+  # drug/drug-class model so each model is evaluated relative to its own
+  # selected feature set.
+  # -------------------------------------------------------------------------
+
+  high_contribution_features <- top_features |>
+    dplyr::group_by(
+      species,
+      drug_label,
+      drug_or_class
     ) |>
     dplyr::mutate(
-      scale_factor = realization_score * coverage_boost
+      contribution_cutoff = stats::quantile(
+        median_contribution,
+        probs = high_contribution_quantile,
+        na.rm = TRUE
+      )
+    ) |>
+    dplyr::filter(
+      median_contribution >= contribution_cutoff
+    ) |>
+    dplyr::ungroup() |>
+    dplyr::select(-contribution_cutoff) |>
+    dplyr::arrange(
+      species,
+      drug_label,
+      drug_or_class,
+      dplyr::desc(median_contribution)
     )
 
-  ## --------------------------------
-  ## 4. Shuffle vs non-shuffle
-  ## --------------------------------
+  # -------------------------------------------------------------------------
+  # 7. Features present across all seeds
+  #
+  # near() avoids problems caused by floating-point representation.
+  # -------------------------------------------------------------------------
 
-  shuffle_delta <- ranked_features |>
-    dplyr::filter(!is.na(cluster)) |>
-    dplyr::group_by(
-      species, drug_label,
+  all_seed_features <- top_features |>
+    dplyr::filter(dplyr::near(seed_ratio, 1)) |>
+    dplyr::arrange(
+      species,
+      drug_label,
       drug_or_class,
-      cluster
+      dplyr::desc(median_rank_score),
+      dplyr::desc(median_contribution)
+    )
+
+  # -------------------------------------------------------------------------
+  # 8. Summarize sharedness of top features
+  #
+  # Sharedness is evaluated independently for drugs and drug classes because
+  # drug_label is part of the grouping.
+  #
+  # Feature subtype is not part of the feature identity here. Thus, binary
+  # and counts representations of the same variable are treated as the same
+  # molecular feature.
+  # -------------------------------------------------------------------------
+
+  feature_sharedness <- top_features |>
+    dplyr::group_by(
+      species,
+      drug_label,
+      feature_type,
+      variable
     ) |>
     dplyr::summarise(
-      mean_rank_nonshuffle = mean(mean_rank[shuffled == FALSE], na.rm = TRUE),
-      mean_rank_shuffle = mean(mean_rank[shuffled == TRUE], na.rm = TRUE),
-      delta_rank = mean_rank_shuffle - mean_rank_nonshuffle,
-      # If non-shuffled is missing -> NA (no evidence). If shuffled missing -> +Inf improvement.
-      improvement = dplyr::case_when(
-        !is.na(mean_rank_nonshuffle) ~ tidyr::replace_na(mean_rank_shuffle, Inf) - mean_rank_nonshuffle,
-        TRUE ~ NA_real_
+      n_drug_or_class = dplyr::n_distinct(drug_or_class),
+      drug_or_class_csv = paste(
+        sort(unique(drug_or_class)),
+        collapse = ", "
       ),
-
-      # "Good" if non-shuffled exists AND (non-shuffled < shuffled OR shuffled is missing)
-      good_feature = !is.na(mean_rank_nonshuffle) & improvement > 0,
+      maximum_seed_ratio = max(
+        seed_ratio,
+        na.rm = TRUE
+      ),
+      median_seed_ratio = median(
+        seed_ratio,
+        na.rm = TRUE
+      ),
+      maximum_rank_score = max(
+        median_rank_score,
+        na.rm = TRUE
+      ),
+      overall_median_rank_score = median(
+        median_rank_score,
+        na.rm = TRUE
+      ),
+      overall_median_contribution = median(
+        median_contribution,
+        na.rm = TRUE
+      ),
+      sign_consistent_across_models =
+        dplyr::n_distinct(sign) == 1,
+      sign = if (sign_consistent_across_models) {
+        dplyr::first(sign)
+      } else {
+        "MIXED"
+      },
       .groups = "drop"
     )
 
-  ## --------------------------------
-  ## 5. Cluster scoring
-  ## --------------------------------
-  robustness <- ranked_features |>
-    dplyr::filter(!is.na(cluster), shuffled == FALSE) |>
-    dplyr::group_by(species, drug_label, drug_or_class, cluster) |>
-    dplyr::summarize(
-      median_contribution = median(contribution),
-      median_rank = median(mean_rank), .groups = "drop"
+  # -------------------------------------------------------------------------
+  # 9. Features shared across drugs or drug classes
+  # -------------------------------------------------------------------------
+
+  shared_features <- feature_sharedness |>
+    dplyr::filter(
+      n_drug_or_class >= min_models_shared
     ) |>
-    dplyr::left_join(
-      protein_scale_summary |>
-        dplyr::distinct(species, drug_label, drug_or_class, cluster, n_scales, scale_factor),
-      by = c("species", "drug_label", "drug_or_class", "cluster")
+    dplyr::arrange(
+      species,
+      drug_label,
+      dplyr::desc(n_drug_or_class),
+      dplyr::desc(overall_median_rank_score),
+      dplyr::desc(overall_median_contribution)
+    )
+
+  # -------------------------------------------------------------------------
+  # 10. Features unique to one drug or drug class
+  #
+  # Join back to top_features so the output retains the associated
+  # drug_or_class and feature-level statistics.
+  # -------------------------------------------------------------------------
+
+  unique_feature_ids <- feature_sharedness |>
+    dplyr::filter(n_drug_or_class == 1) |>
+    dplyr::select(
+      species,
+      drug_label,
+      feature_type,
+      variable
+    )
+
+  unique_features <- top_features |>
+    dplyr::semi_join(
+      unique_feature_ids,
+      by = c(
+        "species",
+        "drug_label",
+        "feature_type",
+        "variable"
+      )
     ) |>
-    dplyr::left_join(
-      shuffle_delta |>
-        dplyr::distinct(species, drug_label, drug_or_class, cluster, delta_rank, good_feature),
-      by = c("species", "drug_label", "drug_or_class", "cluster")
+    dplyr::arrange(
+      species,
+      drug_label,
+      drug_or_class,
+      dplyr::desc(median_rank_score),
+      dplyr::desc(median_contribution)
+    )
+
+  # -------------------------------------------------------------------------
+  # 11. Aggregate top features to protein dyads
+  # -------------------------------------------------------------------------
+
+  top_dyads <- summariseDyads(
+    top_filtered_features = top_features,
+    dyad_feature_parquet = dyad_feature_parquet
+  )
+
+  # -------------------------------------------------------------------------
+  # 12. High-scoring dyads
+  #
+  # A dyad score is high relative to the other selected dyads for the same
+  # species and drug/drug-class model.
+  # -------------------------------------------------------------------------
+
+  high_score_dyads <- top_dyads |>
+    dplyr::group_by(
+      species,
+      drug_label,
+      drug_or_class
     ) |>
-    dplyr::group_by(species, drug_label, drug_or_class) |>
     dplyr::mutate(
-      # higher = better
-      contrib_score = dplyr::percent_rank(median_contribution),
-
-      # lower rank = better → invert
-      stability_score = 1 - dplyr::percent_rank(median_rank),
-
-      # robustness: handle NaN = strongest case (missing in shuffle)
-      delta_score = dplyr::case_when(
-        is.nan(delta_rank) ~ 1, # best possible signal
-        delta_rank > 0 ~ dplyr::percent_rank(delta_rank), # reward
-        delta_rank == 0 ~ 0, # neutral
-        delta_rank < 0 ~ -dplyr::percent_rank(abs(delta_rank)) # penalize
+      dyad_score_cutoff = stats::quantile(
+        dyad_median_rank_score,
+        probs = high_dyad_score_quantile,
+        na.rm = TRUE
+      )
+    ) |>
+    dplyr::filter(
+      dyad_median_rank_score >= dyad_score_cutoff
+    ) |>
+    dplyr::ungroup() |>
+    dplyr::select(-dyad_score_cutoff) |>
+    dplyr::arrange(
+      species,
+      drug_label,
+      drug_or_class,
+      dplyr::desc(dyad_median_rank_score),
+      dplyr::desc(dyad_median_contribution)
+    )
+
+  # -------------------------------------------------------------------------
+  # 13. Summarize dyad sharedness
+  #
+  # A dyad is identified independently within each species and drug-label
+  # level.
+  # -------------------------------------------------------------------------
+
+  dyad_sharedness <- top_dyads |>
+    dplyr::group_by(
+      species,
+      drug_label,
+      dyad
+    ) |>
+    dplyr::summarise(
+      n_drug_or_class = dplyr::n_distinct(drug_or_class),
+      drug_or_class_csv = paste(
+        sort(unique(drug_or_class)),
+        collapse = ", "
+      ),
+      total_frequency = sum(
+        frequency,
+        na.rm = TRUE
+      ),
+      n_distinct_variables = dplyr::n_distinct(
+        unlist(strsplit(variables_csv, ",", fixed = TRUE))
+      ),
+      maximum_dyad_rank_score = max(
+        dyad_median_rank_score,
+        na.rm = TRUE
+      ),
+      overall_dyad_median_rank_score = median(
+        dyad_median_rank_score,
+        na.rm = TRUE
+      ),
+      overall_dyad_median_contribution = median(
+        dyad_median_contribution,
+        na.rm = TRUE
       ),
+      sign_consistent_across_models =
+        dplyr::n_distinct(sign) == 1,
+      sign = if (sign_consistent_across_models) {
+        dplyr::first(sign)
+      } else {
+        "MIXED"
+      },
+      .groups = "drop"
+    )
 
-      # already bounded [0,1]
-      scale_score = scale_factor,
-      robustness_score =
-        contrib_score *
-          stability_score *
-          scale_factor *
-          delta_score
+  # -------------------------------------------------------------------------
+  # 14. Dyads shared across drugs or drug classes
+  # -------------------------------------------------------------------------
+
+  shared_dyads <- dyad_sharedness |>
+    dplyr::filter(
+      n_drug_or_class >= min_models_shared
+    ) |>
+    dplyr::arrange(
+      species,
+      drug_label,
+      dplyr::desc(n_drug_or_class),
+      dplyr::desc(overall_dyad_median_rank_score),
+      dplyr::desc(total_frequency)
     )
 
-  return(robustness)
-}
+  # -------------------------------------------------------------------------
+  # 15. Dyads unique to one drug or drug class
+  #
+  # Join back to top_dyads to retain the drug/class identity and the
+  # original dyad statistics.
+  # -------------------------------------------------------------------------
+
+  unique_dyad_ids <- dyad_sharedness |>
+    dplyr::filter(n_drug_or_class == 1) |>
+    dplyr::select(
+      species,
+      drug_label,
+      dyad
+    )
+
+  unique_dyads <- top_dyads |>
+    dplyr::semi_join(
+      unique_dyad_ids,
+      by = c(
+        "species",
+        "drug_label",
+        "dyad"
+      )
+    ) |>
+    dplyr::arrange(
+      species,
+      drug_label,
+      drug_or_class,
+      dplyr::desc(dyad_median_rank_score),
+      dplyr::desc(dyad_median_contribution)
+    )
+
+  # -------------------------------------------------------------------------
+  # 16. Return an organized analysis object
+  # -------------------------------------------------------------------------
+
+  result <- list(
+    models = list(
+      good_performance = good_performance_models,
+      good_shuffled_separation =
+        good_shuffled_separation_models,
+      qualified = qualified_models,
+      summary = model_summary
+    ),
+    features = list(
+      top = top_features,
+      consistent = consistent_features,
+      high_contribution = high_contribution_features,
+      all_seeds = all_seed_features,
+      sharedness = feature_sharedness,
+      shared = shared_features,
+      unique = unique_features
+    ),
+    dyads = list(
+      top = top_dyads,
+      high_score = high_score_dyads,
+      sharedness = dyad_sharedness,
+      shared = shared_dyads,
+      unique = unique_dyads
+    ),
+    parameters = list(
+      MCC_threshold = MCC_threshold,
+      compare_to_shuffled = compare_to_shuffled,
+      core_contribution_threshold =
+        core_contribution_threshold,
+      exclude_feature_types = exclude_feature_types,
+      filter_model = filter_model,
+      rank_score_quantile = rank_score_quantile,
+      cv_threshold = cv_threshold,
+      cumulative_contribution_threshold =
+        cumulative_contribution_threshold,
+      seed_ratio_threshold = seed_ratio_threshold,
+      found_in_both_subtypes = found_in_both_subtypes,
+      compare_median_to_sd_rank_score =
+        compare_median_to_sd_rank_score,
+      consistent_seed_ratio = consistent_seed_ratio,
+      consistent_rank_score = consistent_rank_score,
+      high_contribution_quantile =
+        high_contribution_quantile,
+      high_dyad_score_quantile =
+        high_dyad_score_quantile,
+      min_models_shared = min_models_shared
+    )
+  )
+
+  class(result) <- c(
+    "feature_dyad_discovery",
+    class(result)
+  )
+
+  return(result)
+}
\ No newline at end of file
diff --git a/data/demo_fit.rda b/data/demo_fit.rda
deleted file mode 100644
index b3fa4ea..0000000
Binary files a/data/demo_fit.rda and /dev/null differ
diff --git a/data/demo_ml_tibble.rda b/data/demo_ml_tibble.rda
deleted file mode 100644
index a09f699..0000000
Binary files a/data/demo_ml_tibble.rda and /dev/null differ
diff --git a/inst/extdata/all_perf.parquet b/inst/extdata/all_perf.parquet
new file mode 100644
index 0000000..808f6c1
Binary files /dev/null and b/inst/extdata/all_perf.parquet differ
diff --git a/inst/extdata/all_top_features.parquet b/inst/extdata/all_top_features.parquet
new file mode 100644
index 0000000..5ff13ed
Binary files /dev/null and b/inst/extdata/all_top_features.parquet differ
diff --git a/inst/extdata/cluster_feature.parquet b/inst/extdata/cluster_feature.parquet
index 1766468..6a4ca86 100644
Binary files a/inst/extdata/cluster_feature.parquet and b/inst/extdata/cluster_feature.parquet differ
diff --git a/inst/extdata/protein_names.parquet b/inst/extdata/protein_names.parquet
new file mode 100644
index 0000000..2d8e080
Binary files /dev/null and b/inst/extdata/protein_names.parquet differ
diff --git a/man/buildFeatureNetwork.Rd b/man/buildFeatureNetwork.Rd
new file mode 100644
index 0000000..3347d0d
--- /dev/null
+++ b/man/buildFeatureNetwork.Rd
@@ -0,0 +1,29 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/feature_rescoring.R
+\name{buildFeatureNetwork}
+\alias{buildFeatureNetwork}
+\title{Build a feature network from selected top features and top dyads}
+\usage{
+buildFeatureNetwork(
+  top_features,
+  top_dyads,
+  dyad_feature_parquet,
+  protein_names_parquet
+)
+}
+\arguments{
+\item{top_features}{Output of \code{topFeaturesPerDrugOrClass()}.}
+
+\item{top_dyads}{Output of \code{summariseDyads()}.}
+
+\item{dyad_feature_parquet}{Path to the same Parquet file passed to \code{summariseDyads()}: a \code{target} column formatted as \code{"<feature_type>:<variable>"} (using the short feature-type codes \code{amr}/\code{cog}/\code{defense}/\code{pfam}/\code{protein}) and a \code{source} column giving the dyad id. \code{target} is reconstructed here from \code{top_features$feature_type}/\code{variable} the same way \code{summariseDyads()} builds it, so the feature-to-dyad edges use the same mapping as the dyad table itself.}
+
+\item{protein_names_parquet}{Path to the Parquet file with dyad name annotations.}
+}
+\value{
+A list with \code{feature_table}, \code{dyad_table}, \code{nodes}, \code{edges}, and \code{graph}.
+Node/edge weights are built from \code{median_rank_score} (features) and \code{dyad_median_rank_score} (dyads, via the \code{dyad_score} column, itself median-based per \code{summariseDyads()}), consistent with the seed-noise-robust selection made in \code{topFeaturesPerDrugOrClass()}.
+}
+\description{
+Build a feature network from selected top features and top dyads
+}
diff --git a/man/demo_fit.Rd b/man/demo_fit.Rd
deleted file mode 100644
index 9f29307..0000000
--- a/man/demo_fit.Rd
+++ /dev/null
@@ -1,23 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/data.R
-\docType{data}
-\name{demo_fit}
-\alias{demo_fit}
-\title{Demo LR fit}
-\format{
-A fitted \code{workflow} object (output of \code{\link[=fitBestModel]{fitBestModel()}}).
-}
-\source{
-\code{inst/scripts/make_demo_data.R}.
-}
-\usage{
-data(demo_fit)
-}
-\description{
-A tuned logistic-regression workflow fitted on \link{demo_ml_tibble}.
-}
-\examples{
-data(demo_fit)
-class(demo_fit)
-}
-\keyword{datasets}
diff --git a/man/demo_ml_tibble.Rd b/man/demo_ml_tibble.Rd
deleted file mode 100644
index 79f685a..0000000
--- a/man/demo_ml_tibble.Rd
+++ /dev/null
@@ -1,26 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/data.R
-\docType{data}
-\name{demo_ml_tibble}
-\alias{demo_ml_tibble}
-\title{Demo ML input tibble}
-\format{
-A tibble with 60 rows and 82 columns: \code{genome_id},
-\code{genome_drug.resistant_phenotype}, and 80 binary feature columns.
-}
-\source{
-\code{inst/scripts/make_demo_data.R}.
-}
-\usage{
-data(demo_ml_tibble)
-}
-\description{
-Stratified subset (30 Resistant + 30 Susceptible) of the AMP-genes-binary
-matrix from the bundled \code{Sfl_parquet.duckdb}, restricted to 80 feature
-columns.
-}
-\examples{
-data(demo_ml_tibble)
-dim(demo_ml_tibble)
-}
-\keyword{datasets}
diff --git a/man/filterOptimalModel.Rd b/man/filterOptimalModel.Rd
new file mode 100644
index 0000000..d39e097
--- /dev/null
+++ b/man/filterOptimalModel.Rd
@@ -0,0 +1,29 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/feature_rescoring.R
+\name{filterOptimalModel}
+\alias{filterOptimalModel}
+\title{Filter model fits to those passing quality thresholds}
+\usage{
+filterOptimalModel(
+  all_performance_parquet,
+  MCC_threshold = NULL,
+  compare_to_shuffled = TRUE
+)
+}
+\arguments{
+\item{all_performance_parquet}{The path to the 'all performance parquet' file}
+
+\item{MCC_threshold}{The minimum non-shuffled MCC required to keep a model fit (default is \code{NULL}, no filtering)}
+
+\item{compare_to_shuffled}{Logical indicating whether to require the non-shuffled MCC to exceed (outperform) the shuffled-label MCC for the same fit (default is \code{TRUE}). A fit with no shuffled counterpart (\code{shuffled_MCC} is \code{NA} after pivoting) always passes this check regardless.}
+}
+\value{
+a tibble of model fits (one row per species, drug_label, drug_or_class, seed, feature_type, feature_subtype, model, fit_penalty, fit_mixture) that pass the requested quality thresholds, with \code{nonshuffled_MCC}, \code{shuffled_MCC}, and \code{MCC_diff} columns added. This does not pick a single "best" fit per group; it filters out fits that fail the MCC/shuffled-comparison criteria, so multiple passing fits per group can remain.
+}
+\description{
+Filter model fits to those passing quality thresholds
+}
+\examples{
+filterOptimalModel(all_performance_parquet = "inst/extdata/all_perf.parquet")
+}
+\keyword{internal}
diff --git a/man/findSharedDyads.Rd b/man/findSharedDyads.Rd
new file mode 100644
index 0000000..35293f3
--- /dev/null
+++ b/man/findSharedDyads.Rd
@@ -0,0 +1,25 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/feature_rescoring.R
+\name{findSharedDyads}
+\alias{findSharedDyads}
+\title{Find dyads that appear across multiple drugs/classes}
+\usage{
+findSharedDyads(
+  top_dyads = summariseDyads(top_features, dyad_feature_parquet),
+  label = "drug",
+  min_drugs_or_classes = 2
+)
+}
+\arguments{
+\item{top_dyads}{The tibble of summarized dyads generated from \code{summarisedyads()}}
+
+\item{label}{The \code{drug_label} value to filter dyads by, either \code{"drug"} or \code{"drug_class"} (default is \code{"drug"})}
+
+\item{min_drugs_or_classes}{The minimum number of distinct drugs or classes required for a dyad to be considered shared (default is \code{2})}
+}
+\value{
+a tibble with one row per shared \code{dyad}, with \code{n_drug_or_class} (the number of distinct \code{drug_or_class} values the dyad appears in) and \code{drug_or_class_csv} (a comma-separated string of those values), sorted by \code{n_drug_or_class} descending.
+}
+\description{
+Find dyads that appear across multiple drugs/classes
+}
diff --git a/man/findUniqueDyads.Rd b/man/findUniqueDyads.Rd
new file mode 100644
index 0000000..7f84464
--- /dev/null
+++ b/man/findUniqueDyads.Rd
@@ -0,0 +1,28 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/feature_rescoring.R
+\name{findUniqueDyads}
+\alias{findUniqueDyads}
+\title{Find the dyads that are unique to a single drug/class}
+\usage{
+findUniqueDyads(
+  top_dyads = summariseDyads(top_features, dyad_feature_parquet),
+  label = "drug",
+  protein_names_parquet
+)
+}
+\arguments{
+\item{top_dyads}{The tibble of summarized dyads generated from \code{summarisedyads()}}
+
+\item{label}{The \code{drug_label} value to filter dyads by, either \code{"drug"} or \code{"drug_class"} (default is \code{"drug"})}
+
+\item{protein_names_parquet}{The path to the Parquet file containing the annotations to protein dyad names}
+}
+\value{
+a tibble with one row per dyad unique to a single drug/class, with \code{drug_or_class}, \code{dyad}, \code{dyad_name} (from the protein name annotations), and \code{dyad_mean_rank_score}, sorted by \code{dyad_mean_rank_score} descending.
+}
+\description{
+Find the dyads that are unique to a single drug/class
+}
+\examples{
+findUniquedyads(summarisedyads(top_features, dyad_feature_parquet), label = "drug", protein_names_parquet)
+}
diff --git a/man/plotFeatureNetworkD3.Rd b/man/plotFeatureNetworkD3.Rd
new file mode 100644
index 0000000..bba4678
--- /dev/null
+++ b/man/plotFeatureNetworkD3.Rd
@@ -0,0 +1,21 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/feature_rescoring.R
+\name{plotFeatureNetworkD3}
+\alias{plotFeatureNetworkD3}
+\title{Plot the feature network with networkD3}
+\usage{
+plotFeatureNetworkD3(feature_network, height = 800, width = "100\%")
+}
+\arguments{
+\item{feature_network}{Output of \code{buildFeatureNetwork()}.}
+
+\item{height}{Widget height in pixels (default is \code{800}).}
+
+\item{width}{Widget width (default is \code{"100\%"}).}
+}
+\value{
+A \code{networkD3} widget.
+}
+\description{
+Plot the feature network with networkD3
+}
diff --git a/man/scoreFeaturesWithinSeed.Rd b/man/scoreFeaturesWithinSeed.Rd
new file mode 100644
index 0000000..c09186c
--- /dev/null
+++ b/man/scoreFeaturesWithinSeed.Rd
@@ -0,0 +1,52 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/feature_rescoring.R
+\name{scoreFeaturesWithinSeed}
+\alias{scoreFeaturesWithinSeed}
+\title{Score features within each seed}
+\usage{
+scoreFeaturesWithinSeed(
+  all_top_features_parquet,
+  core_contribution_threshold = 0.75,
+  exclude_feature_types = NULL,
+  filter_model = TRUE,
+  all_performance_parquet,
+  MCC_threshold = NULL,
+  compare_to_shuffled = TRUE
+)
+}
+\arguments{
+\item{all_top_features_parquet}{The path to the Parquet file containing all top features with their importance scores.}
+
+\item{core_contribution_threshold}{The cumulative-contribution cutoff, in [0, 1] (default is \code{0.75}, i.e. 75\%), used to flag whether a feature falls within the "core" set of features that jointly account for that share of a seed's total importance.}
+
+\item{exclude_feature_types}{Feature types to drop before any scoring happens (default is \code{NULL}, i.e. no feature types are excluded and struct is currently included). struct variables are composite IDs (e.g. \code{polA.group_211.group_2176}, three dot-joined gene/domain identifiers) representing a co-occurrence/structural motif rather than a single molecular entity like the other five scales, and its candidate-variable count (tens of thousands per group) dwarfs the other scales by orders of magnitude — pooling it into this rank_score/contribution machinery would compare a compound signal against five primary ones on an incomparable scale. Pass \code{"struct"} here to exclude it once struct is meant to be reserved for post-hoc biological annotation rather than scoring/ranking/thresholding.}
+
+\item{filter_model}{Logical indicating whether to restrict scoring to (species, drug_label, drug_or_class, feature_type, feature_subtype, seed) groups that have at least one model fit passing \code{filterOptimalModel()}'s MCC/shuffled-comparison quality thresholds (default is \code{TRUE}). The join is not keyed on \code{model}/\code{fit_penalty}/\code{fit_mixture}; this is only safe because there is currently never more than one fit per (species, drug_label, drug_or_class, feature_type, feature_subtype, seed) group. If that assumption changes, this join (and the \code{all_perf} join below) would need a finer-grained key to avoid silently keeping/duplicating rows from multiple fits.}
+
+\item{all_performance_parquet}{The path to the all performance parquet file. Always required — it is read unconditionally}
+
+\item{MCC_threshold}{The minimum MCC threshold passed through to \code{filterOptimalModel()} (default is \code{NULL}, no filtering)}
+
+\item{compare_to_shuffled}{Logical indicating whether to compare the model to shuffled data, passed through to \code{filterOptimalModel()} (default is \code{TRUE})}
+}
+\value{
+a tibble of scored top features, one row per feature within each species/drug_label/drug_or_class/feature_type/feature_subtype/seed group, with the following columns added:
+\itemize{
+\item \code{contribution}: the feature's importance divided by the sum of importance across all features in the group.
+\item \code{feat_return_ratio}: \code{n_feats_returned / n_feat} for the fit that produced this group. This assumes there is never more than one surviving fit per (species, drug_label, drug_or_class, seed, feature_type, feature_subtype) group; if that ever stops being true, this join would need a finer-grained key (\code{all_top_features_parquet} carries no \code{model}/\code{fit_penalty}/\code{fit_mixture} column to join on directly) or the ratio would need to be aggregated across fits before joining, to avoid fanning out and double-counting feature rows.
+\item \code{rank}: descending rank of \code{contribution} within the group; ties receive the average of the ranks they span.
+\item \code{n_features}: the number of rows (features) in the group.
+\item \code{rank_score}: \code{(n_features - rank) / (n_features - 1)}; ranges 0-1 with higher values indicating higher importance (a single-feature group scores 1).
+\item \code{cum_contrib}: the cumulative sum of \code{contribution} in descending order; features tied on \code{contribution} share the same \code{cum_contrib}, equal to the cumulative sum through the end of their tied block, so a tie is never split by arbitrary sort order.
+\item \code{in_core}: TRUE when \code{cum_contrib <= core_contribution_threshold}; a tied block that would push the cumulative total past the threshold is excluded in its entirety (conservative: stays at-or-under the threshold rather than overshooting it).
+}
+}
+\description{
+Score features within each seed
+}
+\examples{
+scoreFeaturesWithinSeed(all_top_features_parquet = "inst/extdata/all_top_features.parquet",
+all_performance_parquet = "inst/extdata/all_perf.parquet")
+
+}
+\keyword{internal}
diff --git a/man/summariseDyads.Rd b/man/summariseDyads.Rd
new file mode 100644
index 0000000..331f3fa
--- /dev/null
+++ b/man/summariseDyads.Rd
@@ -0,0 +1,29 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/feature_rescoring.R
+\name{summariseDyads}
+\alias{summariseDyads}
+\title{Aggregate mapped features to protein dyads}
+\usage{
+summariseDyads(top_filtered_features, dyad_feature_parquet)
+}
+\arguments{
+\item{top_filtered_features}{The tibble of top features for each drug/class generated from \code{topFeaturesPerDrugOrClass()}.}
+
+\item{dyad_feature_parquet}{The path to the Parquet file containing the mapping of features to protein dyads. Must have a \code{target} column formatted as \code{"<feature_type>:<variable>"} (using the short feature-type codes \code{amr}/\code{cog}/\code{defense}/\code{pfam}/\code{protein}) and a \code{source} column giving the dyad id.}
+}
+\value{
+a tibble with one row per species/drug_label/drug_or_class/source (\code{source} is the protein-dyad id), with:
+\itemize{
+\item \code{frequency}: the number of top-feature rows mapped to this dyad.
+\item \code{n_variables}, \code{variables_csv}: number of, and comma-separated list of, distinct \code{target} (\code{"<feature_type>:<variable>"}) values mapped to this dyad.
+\item \code{n_feature_types}, \code{feature_types_csv}: number of, and comma-separated list of, distinct \code{feature_type_subtype} (\code{"<feature_type>:<feature_subtype>"}) values mapped to this dyad.
+\item \code{dyad_median_rank_score}: median of \code{median_rank_score} across the top-feature rows mapped to this dyad.
+\item \code{dyad_median_contribution}: median of \code{median_contribution} across the top-feature rows mapped to this dyad.
+\item \code{sign_consistent}, \code{sign}: whether \code{sign} is identical across every top-feature row mapped to this dyad; \code{sign} is that shared value if consistent, else \code{"MIXED"}.
+}
+\code{top_filtered_features} rows whose \code{target} has no match in \code{dyad_feature_parquet} are not dropped: they collapse into one \code{source = NA} row per species/drug_label/drug_or_class, aggregating every unmapped feature for that group. Callers that want only real dyads must filter this out explicitly (e.g. \code{dplyr::filter(!is.na(source))}).
+}
+\description{
+Internal helper called on the output of \code{topFeaturesPerDrugOrClass()}.
+}
+\keyword{internal}
diff --git a/man/summariseFeaturesAcrossSeeds.Rd b/man/summariseFeaturesAcrossSeeds.Rd
new file mode 100644
index 0000000..5f428f5
--- /dev/null
+++ b/man/summariseFeaturesAcrossSeeds.Rd
@@ -0,0 +1,35 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/feature_rescoring.R
+\name{summariseFeaturesAcrossSeeds}
+\alias{summariseFeaturesAcrossSeeds}
+\title{Summarize a feature's scoring across seeds}
+\usage{
+summariseFeaturesAcrossSeeds(scored_top_features)
+}
+\arguments{
+\item{scored_top_features}{The tibble of scored top features with their contribution, rank, and rank score within each seed generated from \code{scoreFeaturesWithinSeed()}}
+}
+\value{
+a tibble with one row per species/drug_label/drug_or_class/feature_type/feature_subtype/variable, with:
+\itemize{
+\item \code{seed_ratio}: the number of seeds the feature appears in, divided by the total number of distinct \code{seed} values present anywhere in \code{scored_top_features} (a single count computed once for the whole call, not per feature or per group — so this assumes every group was fit with the same set of seeds).
+\item \code{mean_rank_score}, \code{median_rank_score}: mean/median of \code{rank_score} across seeds; ranges 0-1 with higher values indicating higher importance.
+\item \code{median_rank}: median of \code{rank} across seeds. (\code{mean_rank} is not currently computed.)
+\item \code{median_contribution}: median of \code{contribution} across seeds. (\code{mean_contribution} is not currently computed.)
+\item \code{median_cum_contrib}: median of \code{cum_contrib} across seeds.
+\item \code{best_rank}: the best (lowest) rank seen across seeds.
+\item \code{rank_consistent}: TRUE if the feature's rank is identical in every seed.
+\item \code{rank_score_sd}, \code{rank_score_cv}: standard deviation and coefficient of variation of \code{rank_score} across seeds, i.e. how much the normalized rank_score varies — computed on \code{rank_score} rather than raw rank so it is comparable across groups with different numbers of features.
+\item \code{in_core_consistent}, \code{in_core}: whether the feature's \code{in_core} flag is identical across every seed; \code{in_core} is that shared value if consistent, else \code{FALSE}.
+\item \code{sign_consistent}, \code{sign}: whether the feature's \code{sign} is identical across every seed; \code{sign} is that shared value if consistent, else \code{"MIXED"}.
+}
+}
+\description{
+Answers "which molecular features are consistently important?" by
+collapsing the per-seed rows from \code{scoreFeaturesWithinSeed()} down to one
+row per feature.
+}
+\examples{
+summariseFeaturesAcrossSeeds(scoreFeaturesWithinSeed(all_top_features.parquet))
+}
+\keyword{internal}
diff --git a/man/topFeaturesPerDrugOrClass.Rd b/man/topFeaturesPerDrugOrClass.Rd
new file mode 100644
index 0000000..fb52a01
--- /dev/null
+++ b/man/topFeaturesPerDrugOrClass.Rd
@@ -0,0 +1,69 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/feature_rescoring.R
+\name{topFeaturesPerDrugOrClass}
+\alias{topFeaturesPerDrugOrClass}
+\title{Build a per-drug top-feature table with cutoffs}
+\usage{
+topFeaturesPerDrugOrClass(
+  all_top_features_parquet,
+  core_contribution_threshold = 0.75,
+  exclude_feature_types = NULL,
+  filter_model = TRUE,
+  all_performance_parquet,
+  MCC_threshold = NULL,
+  compare_to_shuffled = TRUE,
+  rank_score_quantile = 0.95,
+  cv_threshold = 1,
+  cumulative_contribution_threshold = 0.75,
+  seed_ratio_threshold = NULL,
+  found_in_both_subtypes = FALSE,
+  compare_median_to_sd_rank_score = FALSE
+)
+}
+\arguments{
+\item{all_top_features_parquet}{The path to the Parquet file containing all top features with their importance scores.}
+
+\item{core_contribution_threshold}{The cumulative-contribution cutoff, in [0, 1] (default is \code{0.75}, i.e. 75\%), used to flag whether a feature falls within the "core" set of features that jointly account for that share of a seed's total importance.}
+
+\item{exclude_feature_types}{Feature types to drop before any scoring happens (default is \code{NULL}, i.e. no feature types are excluded and struct is currently included). struct variables are composite IDs (e.g. \code{polA.group_211.group_2176}, three dot-joined gene/domain identifiers) representing a co-occurrence/structural motif rather than a single molecular entity like the other five scales, and its candidate-variable count (tens of thousands per group) dwarfs the other scales by orders of magnitude — pooling it into this rank_score/contribution machinery would compare a compound signal against five primary ones on an incomparable scale. Pass \code{"struct"} here to exclude it once struct is meant to be reserved for post-hoc biological annotation rather than scoring/ranking/thresholding.}
+
+\item{filter_model}{Logical indicating whether to restrict scoring to (species, drug_label, drug_or_class, feature_type, feature_subtype, seed) groups that have at least one model fit passing \code{filterOptimalModel()}'s MCC/shuffled-comparison quality thresholds (default is \code{TRUE}). The join is not keyed on \code{model}/\code{fit_penalty}/\code{fit_mixture}; this is only safe because there is currently never more than one fit per (species, drug_label, drug_or_class, feature_type, feature_subtype, seed) group. If that assumption changes, this join (and the \code{all_perf} join below) would need a finer-grained key to avoid silently keeping/duplicating rows from multiple fits.}
+
+\item{all_performance_parquet}{The path to the all performance parquet file. Always required — it is read unconditionally}
+
+\item{MCC_threshold}{The minimum MCC threshold passed through to \code{filterOptimalModel()} (default is \code{NULL}, no filtering)}
+
+\item{compare_to_shuffled}{Logical indicating whether to compare the model to shuffled data, passed through to \code{filterOptimalModel()} (default is \code{TRUE})}
+
+\item{rank_score_quantile}{A value in [0, 1] (default is \code{0.95}). Keep only features whose \code{median_rank_score} is at or above this quantile of \code{median_rank_score}, where the quantile is computed separately within each (species, drug_label, drug_or_class) group (with \code{na.rm = TRUE}) — so a feature is compared only against other features from the same drug/class, not against the whole panel. This means a drug/class with weaker or more diffuse signal overall can still contribute its own top features, rather than being crowded out by drugs/classes with stronger or more concentrated signal. As with the other conditions listed below, this is not restricted to rows that already pass them: \code{dplyr::filter()} evaluates every condition passed to a single call against the same original data, so this threshold does not narrow as other conditions are applied}
+
+\item{cv_threshold}{The maximum allowed coefficient of variation (default is \code{1}). Keep only features with \code{rank_score_cv <= cv_threshold}, i.e. drop features whose rank_score is inconsistent across seeds relative to its mean. A feature present in every seed (\code{seed_ratio == 1}) always passes this check regardless of its \code{rank_score_cv} — including when \code{rank_score_cv} is \code{NA}, which happens whenever a group has only a single seed}
+
+\item{cumulative_contribution_threshold}{The cumulative-contribution cutoff, in [0, 1] (default is \code{0.75}, i.e. 75\%). Keep only features with \code{median_cum_contrib <= cumulative_contribution_threshold}}
+
+\item{seed_ratio_threshold}{If not \code{NULL} (the default), keep only features whose \code{seed_ratio} exactly equals this value}
+
+\item{found_in_both_subtypes}{Logical indicating whether to additionally restrict to features that survive the filters above in both the binary and counts \code{feature_subtype} (default is \code{FALSE})}
+
+\item{compare_median_to_sd_rank_score}{Logical indicating whether to additionally require \code{median_rank_score > rank_score_sd} (default is \code{FALSE})}
+}
+\value{
+a tibble of top features for each drug/class: the \code{summariseFeaturesAcrossSeeds()} output (species, drug label, drug or class, feature type, feature subtype, variable, seed_ratio, mean/median rank score, median rank, median contribution, median cumulative contribution, best rank, rank/sign/in_core consistency flags, sign), filtered in two \code{dplyr::filter()} passes.
+The first pass keeps rows where all of the following hold, evaluated together against \code{feature_summary} (see \code{rank_score_quantile} for how the last condition is grouped):
+\itemize{
+\item \code{seed_ratio == seed_ratio_threshold}, only applied when \code{seed_ratio_threshold} is not \code{NULL},
+\item \code{seed_ratio == 1} OR \code{rank_score_cv <= cv_threshold} (see \code{cv_threshold}),
+\item \code{in_core} is TRUE,
+\item \code{sign_consistent} is TRUE (sign is the same in every seed; this does not require the sign to be negative),
+\item \code{median_cum_contrib <= cumulative_contribution_threshold}, and
+\item \code{median_rank_score} is at or above the \code{rank_score_quantile} quantile of \code{median_rank_score}, computed separately per (species, drug_label, drug_or_class) group.
+}
+Two columns are then added, grouped by (species, drug_label, drug_or_class, feature_type, variable): \code{n_subtype} and \code{subtype_csv}, recording how many/which \code{feature_subtype} values each combination has among the rows that survived the first pass.
+A second \code{dplyr::filter()} pass then optionally keeps only rows where \code{subtype_csv == "binary,counts"} (when \code{found_in_both_subtypes = TRUE}) and/or \code{median_rank_score > rank_score_sd} (when \code{compare_median_to_sd_rank_score = TRUE}).
+Every drug/class may not have variables from all feature types. Because the \code{rank_score_quantile} cutoff is now per drug/class rather than global, every drug/class with at least one feature surviving the other conditions will generally contribute something here, rather than potentially being excluded entirely by comparison to stronger-signal drugs/classes elsewhere in the panel.
+}
+\description{
+Runs \code{scoreFeaturesWithinSeed()} and \code{summariseFeaturesAcrossSeeds()} on
+\code{all_top_features_parquet}, then filters the resulting per-feature summary
+down to the top features for each drug/class.
+}