diff --git a/DESCRIPTION b/DESCRIPTION index 4191ecb..1ed16f7 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -53,9 +53,11 @@ Imports: glue, graphics, grDevices, - jsonlite, hardhat, + igraph, + jsonlite, methods, + networkD3, parsnip, purrr, patchwork, @@ -89,4 +91,4 @@ biocViews: Visualization URL: https://github.com/JRaviLab/amRml BugReports: https://github.com/JRaviLab/amRml/issues -Config/roxygen2/version: 8.0.0 +Config/roxygen2/version: 8.1.0 diff --git a/NAMESPACE b/NAMESPACE index e4fb286..717843b 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -2,8 +2,25 @@ export(.getTargetVarName) export(applyBenjaminiHochberg) +export(buildFeatureNetwork) export(buildLRModel) +export(buildPerfPq) +export(buildPerfPqCrossCountry) +export(buildPerfPqCrossDrug) +export(buildPerfPqCrossYear) +export(buildPerfPqLOOCountry) +export(buildPerfPqLOODrug) +export(buildPerfPqLOOYear) +export(buildPerfPqMDR) +export(buildPerfPqYearCountry) +export(buildPredPqMDR) export(buildRecipe) +export(buildTopFeatsPq) +export(buildTopFeatsPqLOOCountry) +export(buildTopFeatsPqLOODrug) +export(buildTopFeatsPqLOOYear) +export(buildTopFeatsPqMDR) +export(buildTopFeatsPqYearCountry) export(buildTuningGrid) export(buildWflow) export(calculateEvalMets) @@ -13,16 +30,20 @@ export(createMLResultDir) export(createMLinputList) export(encodePhenotype) export(extractTopFeats) +export(findSharedDyads) +export(findUniqueDyads) export(fitBestModel) export(generateMLInputs) export(getConfusionMatrix) export(getNumFeat) export(loadMLInputTibble) +export(parse_ml_filename) export(plotBaselineComparison) export(plotCM) export(plotCrossDrug) export(plotDrugDist) export(plotDrugPerf) +export(plotFeatureNetworkD3) export(plotFishers) export(plotMDR) export(plotPRC) @@ -43,6 +64,7 @@ export(runModelingPipeline) export(selectBestModel) export(shuffleLabels) export(splitMLInputTibble) +export(topFeaturesPerDrugOrClass) export(tuneGrid) import(arrow) import(dplyr) @@ -51,116 +73,154 @@ import(ggrepel) import(purrr) import(stringr) import(tibble) -importFrom(DBI,dbConnect) -importFrom(DBI,dbDisconnect) -importFrom(DBI,dbExecute) -importFrom(DBI,dbGetQuery) -importFrom(DBI,dbWriteTable) -importFrom(arrow,read_parquet) -importFrom(arrow,write_parquet) -importFrom(dplyr,all_of) -importFrom(dplyr,arrange) -importFrom(dplyr,bind_rows) -importFrom(dplyr,case_when) -importFrom(dplyr,count) -importFrom(dplyr,desc) -importFrom(dplyr,distinct) -importFrom(dplyr,filter) -importFrom(dplyr,group_by) -importFrom(dplyr,if_any) -importFrom(dplyr,matches) -importFrom(dplyr,mutate) -importFrom(dplyr,n_distinct) -importFrom(dplyr,pull) -importFrom(dplyr,relocate) -importFrom(dplyr,rename) -importFrom(dplyr,row_number) -importFrom(dplyr,rowwise) -importFrom(dplyr,select) -importFrom(dplyr,slice) -importFrom(dplyr,summarise) -importFrom(dplyr,ungroup) +importFrom(DBI, + dbConnect, + dbDisconnect, + dbExecute, + dbGetQuery, + dbWriteTable +) +importFrom(arrow, + read_parquet, + write_parquet +) +importFrom(dplyr, + all_of, + arrange, + bind_rows, + case_when, + count, + desc, + distinct, + filter, + group_by, + if_any, + matches, + mutate, + n_distinct, + pull, + relocate, + rename, + row_number, + rowwise, + select, + slice, + summarise, + ungroup +) importFrom(duckdb,duckdb) -importFrom(ggplot2,aes) -importFrom(ggplot2,element_blank) -importFrom(ggplot2,element_line) -importFrom(ggplot2,element_text) -importFrom(ggplot2,geom_line) -importFrom(ggplot2,geom_path) -importFrom(ggplot2,geom_point) -importFrom(ggplot2,ggplot) -importFrom(ggplot2,labs) -importFrom(ggplot2,theme) -importFrom(ggplot2,xlab) -importFrom(ggplot2,ylim) +importFrom(ggplot2, + aes, + element_blank, + element_line, + element_text, + geom_line, + geom_path, + geom_point, + ggplot, + labs, + theme, + xlab, + ylim +) importFrom(glmnet,glmnet) importFrom(glue,glue) importFrom(grDevices,colorRampPalette) importFrom(graphics,barplot) importFrom(hardhat,tune) -importFrom(jsonlite,fromJSON) -importFrom(jsonlite,write_json) +importFrom(jsonlite, + fromJSON, + write_json +) importFrom(methods,is) -importFrom(parsnip,augment) -importFrom(parsnip,boost_tree) -importFrom(parsnip,extract_fit_engine) -importFrom(parsnip,fit) -importFrom(parsnip,logistic_reg) -importFrom(parsnip,multinom_reg) -importFrom(parsnip,rand_forest) -importFrom(parsnip,set_engine) -importFrom(parsnip,set_mode) -importFrom(purrr,imap_dfr) -importFrom(purrr,map) -importFrom(purrr,map_int) -importFrom(purrr,pmap_chr) -importFrom(purrr,walk) +importFrom(parsnip, + augment, + boost_tree, + extract_fit_engine, + fit, + logistic_reg, + multinom_reg, + rand_forest, + set_engine, + set_mode +) +importFrom(purrr, + imap_dfr, + map, + map_int, + pmap_chr, + walk +) importFrom(readr,write_lines) -importFrom(recipes,all_predictors) -importFrom(recipes,recipe) -importFrom(recipes,step_normalize) -importFrom(recipes,step_pca) -importFrom(recipes,step_zv) -importFrom(recipes,update_role) -importFrom(rlang,":=") -importFrom(rlang,eval_tidy) -importFrom(rlang,sym) -importFrom(rsample,initial_split) -importFrom(rsample,initial_validation_split) -importFrom(rsample,testing) -importFrom(rsample,training) -importFrom(rsample,validation_set) -importFrom(rsample,vfold_cv) -importFrom(stats,coef) -importFrom(stats,fisher.test) -importFrom(stats,median) -importFrom(stats,reformulate) -importFrom(stats,reorder) -importFrom(stats,sd) -importFrom(stringr,str_remove) -importFrom(stringr,str_split) -importFrom(tibble,add_column) -importFrom(tibble,is_tibble) -importFrom(tibble,tibble) -importFrom(tidyr,drop_na) -importFrom(tidyr,pivot_longer) -importFrom(tidyr,pivot_wider) -importFrom(tune,control_grid) -importFrom(tune,extract_fit_parsnip) -importFrom(tune,finalize_workflow) -importFrom(tune,select_best) -importFrom(tune,tune_grid) -importFrom(vip,vi) -importFrom(vip,vip) -importFrom(workflows,add_model) -importFrom(workflows,add_recipe) -importFrom(workflows,workflow) -importFrom(workflowsets,extract_fit_parsnip) -importFrom(workflowsets,extract_spec_parsnip) -importFrom(yardstick,bal_accuracy) -importFrom(yardstick,conf_mat) -importFrom(yardstick,f_meas) -importFrom(yardstick,mcc) -importFrom(yardstick,metric_set) -importFrom(yardstick,pr_auc) -importFrom(yardstick,pr_curve) +importFrom(recipes, + all_predictors, + recipe, + step_normalize, + step_pca, + step_zv, + update_role +) +importFrom(rlang, + ":=", + eval_tidy, + sym +) +importFrom(rsample, + initial_split, + initial_validation_split, + testing, + training, + validation_set, + vfold_cv +) +importFrom(stats, + coef, + fisher.test, + median, + reformulate, + reorder, + sd +) +importFrom(stringr, + str_remove, + str_split +) +importFrom(tibble, + add_column, + is_tibble, + tibble +) +importFrom(tidyr, + drop_na, + pivot_longer, + pivot_wider +) +importFrom(tune, + control_grid, + extract_fit_parsnip, + finalize_workflow, + select_best, + tune_grid +) +importFrom(vip, + vi, + vip +) +importFrom(workflows, + add_model, + add_recipe, + workflow +) +importFrom(workflowsets, + extract_fit_parsnip, + extract_spec_parsnip +) +importFrom(yardstick, + bal_accuracy, + conf_mat, + f_meas, + mcc, + metric_set, + pr_auc, + pr_curve +) 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..7ae0fd0 100644 --- a/R/feature_rescoring.R +++ b/R/feature_rescoring.R @@ -1,113 +1,167 @@ -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 - ) |> - 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 - ) + seed, feature_type, feature_subtype, + model, fit_penalty, fit_mixture, + shuffled, mcc ) |> - 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) + tidyr::pivot_wider( + names_from = shuffled, + values_from = mcc, + names_prefix = "shuffled_" ) |> - 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::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 (regardless of \code{filter_model} or \code{add_lasso_advtg}) to compute the per-fit sparsity score. +#' @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}) +#' @param add_lasso_advtg Logical indicating whether to weight each feature's contribution by a sparsity score that rewards model fits returning fewer features relative to the candidate feature space (default is \code{FALSE}) +#' +#' @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{sparsity_score}: \code{1 - feat_return_ratio} when \code{add_lasso_advtg = TRUE} (fits returning fewer features relative to the candidate space score closer to 1), else 1 for every row. +#' \item \code{adjusted_contribution}: \code{contribution * sparsity_score}. This, not raw \code{contribution}, is what every downstream column below is actually computed from. +#' \item \code{rank}: descending rank of \code{adjusted_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{adjusted_contribution} in descending order; features tied on \code{adjusted_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, + add_lasso_advtg = FALSE) + { + # 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 + ) + ) + } + + # The models are elastic net, fit with a mix of penalties from lasso to ridge. + # Lasso shrinks the feature space to fewer selected variables, while + # ridge retains variables (keeps correlated variables together). + # Calculate the sparsity score to give an advantage to lasso-like fits. + # + # NOTE: all_top_features_parquet has no model/fit_penalty/fit_mixture column, + # so this join is keyed at the (species, drug_label, drug_or_class, seed, + # feature_type, feature_subtype) level. That is only safe because there is + # currently never more than one surviving fit per group; the check below + # enforces that assumption instead of silently fanning out and + # double-counting feature rows if it were ever violated. + all_perf <- arrow::read_parquet(normalizePath(all_performance_parquet)) |> dplyr::filter(!shuffled) |> - dplyr::select( - species, drug_label, drug_or_class, - feature_type, feature_subtype, seed, Variable, - Importance, Sign - ) |> + 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( - 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 - ) + feat_return_ratio = n_feats_returned / n_feat, + sparsity_score = if (add_lasso_advtg) 1 - feat_return_ratio else 1 ) - cluster_feature <- arrow::read_parquet(normalizePath(cluster_feature_parquet)) - ## -------------------------------- - ## 1. Scale capability (data-driven) - ## -------------------------------- + dup_groups <- all_perf |> + dplyr::count(species, drug_label, drug_or_class, seed, + feature_type, feature_subtype) |> + dplyr::filter(n > 1) - 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" + if (nrow(dup_groups) > 0) { + stop( + "scoreFeaturesWithinSeed() assumes at most one fit per (species, ", + "drug_label, drug_or_class, seed, feature_type, feature_subtype) group, ", + "but ", nrow(dup_groups), " group(s) in all_performance_parquet have ", + "more than one. Re-check filterOptimalModel()'s thresholds, or update ", + "this join to key on model/fit_penalty/fit_mixture as well." ) + } - ## -------------------------------- - ## 2. Importance → contribution → rank - ## -------------------------------- - - ranked_features <- top_features |> + # 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, seed, + feature_type, feature_subtype, variable = Variable, + importance = Importance, sign = Sign + ) |> + dplyr::mutate( + variable = dplyr::case_when( + feature_type == "protein" ~ sub("fig.", "fig|", variable, fixed = TRUE), + feature_type == "AMRFinder" ~ sub( + "^X", "", + gsub("\\.NCBIFAM", "", variable) + ), + TRUE ~ variable + ) + ) |> + 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 +171,626 @@ 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), + adjusted_contribution = contribution * sparsity_score + ) |> + dplyr::arrange(dplyr::desc(adjusted_contribution), .by_group = TRUE) |> + dplyr::mutate( + rank = rank(dplyr::desc(adjusted_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(adjusted_contribution) + ) |> + dplyr::group_by(adjusted_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{adjusted_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(adjusted_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}. The quantile is computed once over every row of \code{feature_summary} (with \code{na.rm = TRUE}) — globally across all species/drugs/drug classes/feature types/subtypes, not per group — and is not restricted to rows that already pass the other conditions listed below: \code{dplyr::filter()} evaluates every condition passed to a single call against the same original, ungrouped 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 the full, ungrouped \code{feature_summary} (see \code{rank_score_quantile} for what that means for the last condition): +#' \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}. +#' } +#' 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. +#' +#' @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, + add_lasso_advtg = FALSE, + 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, + add_lasso_advtg = add_lasso_advtg + ) + + feature_summary <- summariseFeaturesAcrossSeeds(scored_features) + + top_filtered_features <- feature_summary |> + 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 >= quantile(median_rank_score, rank_score_quantile, na.rm = TRUE) + ) |> + 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 + ) - ## -------------------------------- - ## 3. Collapse to protein level - ## -------------------------------- + return(top_filtered_features) +} - feature_cluster <- arrow::read_parquet(normalizePath(feature_cluster_parquet)) +#' 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{target} column formatted as \code{":"} (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. +#' +#' @returns 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{":"}) 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{":"}) 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))}). +#' +#' @keywords internal +summariseDyads <- function(top_filtered_features, + dyad_feature_parquet + ) { + stopifnot(is.data.frame(top_filtered_features)) + stopifnot(file.exists(dyad_feature_parquet)) + + dyad_feature <- arrow::read_parquet(dyad_feature_parquet) + if (!all(c("source", "target") %in% names(dyad_feature))) { + stop( + "dyad_feature_parquet is expected to have 'source' (dyad id) and ", + "'target' (':') columns; found: ", + paste(names(dyad_feature), collapse = ", ") + ) + } - 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)) - - protein_scale_realization <- ranked_features |> - dplyr::filter(!is.na(cluster), shuffled == FALSE) |> - dplyr::group_by( - species, - drug_label, - drug_or_class, - cluster, - feature_type + tidyr::unite("target", feature_type, variable, sep = ":", remove = FALSE) |> + tidyr::unite("feature_type_subtype", feature_type, feature_subtype, sep = ":", remove = FALSE) |> + dplyr::left_join(dyad_feature, by = "target") |> + dplyr::select(species, drug_label, drug_or_class, feature_type_subtype, source, target, sign, + median_rank_score, median_contribution, subtype_csv) |> + dplyr::group_by(species, drug_label, drug_or_class, source) |> + dplyr::summarise( + frequency = dplyr::n(), + n_variables = dplyr::n_distinct(target), + variables_csv = paste(sort(unique(target)), 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(frequency), dplyr::desc(n_feature_types)) + + 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{target} column formatted as \code{":"} (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. +#' @param protein_names_parquet Path to the Parquet file with dyad name annotations. +#' +#' @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()}. +#' @export +buildFeatureNetwork <- function(top_features, + top_dyads, + dyad_feature_parquet, + protein_names_parquet + ) { + stopifnot(is.data.frame(top_features)) + stopifnot(is.data.frame(top_dyads)) + stopifnot(file.exists(dyad_feature_parquet)) + stopifnot(file.exists(protein_names_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", + "source", "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() + + if (!all(c("source", "target") %in% names(dyad_feature))) { + stop( + "dyad_feature_parquet is expected to have 'source' (dyad id) and ", + "'target' (':') columns, matching what ", + "summariseDyads() expects; found: ", paste(names(dyad_feature), collapse = ", ") + ) + } + + protein_names <- arrow::read_parquet(normalizePath(protein_names_parquet)) |> + dplyr::distinct() + + 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 `target`. + 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::mutate(model_id = make_model_id(drug_label, drug_or_class)) |> + dplyr::group_by(species, model_id, feature_type, variable) |> + dplyr::summarise( + # 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" + ) + + dyad_table <- top_dyads |> + dplyr::mutate(model_id = make_model_id(drug_label, drug_or_class)) |> + dplyr::group_by(species, model_id, source) |> + dplyr::summarise( + dyad_score = mean(dyad_median_rank_score, na.rm = TRUE), + .groups = "drop" + ) + + model_nodes <- dplyr::bind_rows( + feature_table |> + dplyr::distinct(species, model_id), + dyad_table |> + dplyr::distinct(species, model_id) + ) |> + dplyr::distinct(species, model_id) |> + dplyr::transmute( + name = model_id, + label = model_id, + node_type = "model", + species = species, + score = NA_real_, + breadth = NA_real_, + node_size = 4 + ) + + feature_nodes <- feature_table |> + dplyr::group_by(species, variable) |> dplyr::summarise( - observed_types = dplyr::n_distinct(feature_subtype), - observed_types_csv = paste(sort(unique(feature_subtype)), collapse = ","), + score = mean(feature_score, na.rm = TRUE), + breadth = dplyr::n_distinct(model_id), .groups = "drop" ) |> - dplyr::left_join(scale_capability, by = "feature_type") |> - dplyr::mutate( - scale_realization = dplyr::case_when( - observed_types == expected_types ~ "full_realization", - observed_types < expected_types ~ "partial_realization" - ) + dplyr::transmute( + name = variable, + label = variable, + node_type = "feature", + species = species, + score = score, + breadth = breadth, + node_size = pmax(3, pmin(10, breadth + 2)) ) - protein_scale_summary <- protein_scale_realization |> - dplyr::group_by(species, drug_label, drug_or_class, cluster) |> + dyad_nodes <- dyad_table |> + dplyr::group_by(species, source) |> 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 = "; " - ), + score = median(dyad_score, na.rm = TRUE), + breadth = dplyr::n_distinct(model_id), .groups = "drop" ) |> - dplyr::mutate( - realization_score = - (fully_realized_scales + - 0.5 * partially_realized_scales) / - (fully_realized_scales + partially_realized_scales) + dplyr::transmute( + name = source, + label = source, + node_type = "dyad", + species = species, + score = score, + breadth = breadth, + node_size = pmax(3, pmin(10, breadth + 2)) + ) + + nodes <- dplyr::bind_rows(model_nodes, feature_nodes, dyad_nodes) |> + dplyr::distinct(name, .keep_all = TRUE) + + feature_edges <- feature_table |> + dplyr::transmute( + from = model_id, + to = variable, + weight = feature_score, + edge_type = "model_feature" ) |> - dplyr::mutate( - coverage_boost = - n_scales / max(n_scales) + dplyr::distinct(from, to, edge_type, .keep_all = TRUE) + + dyad_edges <- dyad_table |> + dplyr::transmute( + from = model_id, + to = source, + weight = dyad_score, + edge_type = "model_dyad" ) |> - dplyr::mutate( - scale_factor = realization_score * coverage_boost + dplyr::distinct(from, to, edge_type, .keep_all = TRUE) + + feature_dyad_edges <- feature_table |> + dplyr::mutate(short_feature_type = shorten_feature_type(feature_type)) |> + tidyr::unite("target", short_feature_type, variable, sep = ":", remove = FALSE) |> + dplyr::left_join( + dyad_feature |> dplyr::add_count(target, name = "n_dyads"), + by = "target", + relationship = "many-to-many" + ) |> + dplyr::filter(!is.na(source)) |> + dplyr::transmute( + from = variable, + to = source, + weight = 1 / n_dyads, + edge_type = "feature_dyad" + ) |> + dplyr::distinct(from, to, edge_type, .keep_all = TRUE) + + edges <- dplyr::bind_rows(feature_edges, feature_dyad_edges, dyad_edges) + + missing_vertices <- setdiff(unique(c(edges$from, edges$to)), nodes$name) + if (length(missing_vertices) > 0) { + extra_nodes <- tibble::tibble(name = missing_vertices) |> + dplyr::mutate( + label = name, + node_type = dplyr::case_when( + grepl("^drug\\.|^drug_class\\.", name) ~ "model", + grepl("^fig\\||^dyad", name) ~ "dyad", + TRUE ~ "feature" + ), + species = NA_character_, + score = NA_real_, + breadth = NA_real_, + node_size = 4 + ) + + nodes <- dplyr::bind_rows(nodes, extra_nodes) |> + dplyr::distinct(name, .keep_all = TRUE) + } + + graph <- if (nrow(edges) > 0) { + igraph::graph_from_data_frame( + d = edges, + directed = FALSE, + vertices = nodes + ) + } else { + igraph::graph_from_data_frame( + d = data.frame(from = character(), to = character()), + directed = FALSE, + vertices = nodes ) + } - ## -------------------------------- - ## 4. Shuffle vs non-shuffle - ## -------------------------------- + feature_network <- list( + feature_table = feature_table, + dyad_table = dyad_table, + nodes = nodes, + edges = edges, + graph = graph + ) - shuffle_delta <- ranked_features |> - dplyr::filter(!is.na(cluster)) |> - dplyr::group_by( - species, drug_label, - drug_or_class, - cluster - ) |> - 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_ - ), + return(feature_network) +} - # "Good" if non-shuffled exists AND (non-shuffled < shuffled OR shuffled is missing) - good_feature = !is.na(mean_rank_nonshuffle) & improvement > 0, - .groups = "drop" +#' 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( + "", label, "", + ifelse(is.na(species), "", paste0("
Species: ", species)), + ifelse(is.na(score), "", paste0("
Score: ", signif(score, 3))), + ifelse(is.na(breadth), "", paste0("
Breadth: ", breadth)) + ) ) - ## -------------------------------- - ## 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" - ) |> + links <- feature_network$edges |> + dplyr::filter(!is.na(from), !is.na(to)) |> 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") + nodes |> dplyr::select(name, id), + by = c("from" = "name") ) |> + dplyr::rename(source = id) |> 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") + nodes |> dplyr::select(name, id), + by = c("to" = "name") ) |> - dplyr::group_by(species, drug_label, drug_or_class) |> + dplyr::rename(target = id) |> + dplyr::filter(!is.na(source), !is.na(target)) |> 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 - ), + value = dplyr::if_else(is.na(weight), 1, weight) + ) |> + dplyr::select(source, target, value, edge_type) + + stopifnot(nrow(nodes) > 0) + stopifnot(nrow(links) > 0) - # already bounded [0,1] - scale_score = scale_factor, - robustness_score = - contrib_score * - stability_score * - scale_factor * - delta_score + colour_scale <- networkD3::JS( + "d3.scaleOrdinal() + .domain(['model', 'feature', 'dyad']) + .range(['#4C78A8', '#F58518', '#54A24B'])" + ) + + networkD3::forceNetwork( + Links = links, + Nodes = nodes, + Source = "source", + Target = "target", + Value = "value", + NodeID = "label", + 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; + }" ) + ) +} - return(robustness) +#' Find dyads that appear across multiple drugs/classes +#' +#' @param top_dyads The tibble of summarized dyads generated from `summarisedyads()` +#' @param label The \code{drug_label} value to filter dyads by, either \code{"drug"} or \code{"drug_class"} (default is \code{"drug"}) +#' @param min_drugs_or_classes The minimum number of distinct drugs or classes required for a dyad to be considered shared (default is \code{2}) +#' +#' @returns 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. +#' +#' @export +findSharedDyads <- function(top_dyads = summariseDyads(top_features, dyad_feature_parquet), + label = "drug", + min_drugs_or_classes = 2 + ) { + shared_dyads <- top_dyads |> + dplyr::filter(!is.na(source), drug_label == label) |> + dplyr::group_by(source) |> + dplyr::mutate( + n_drug_or_class = dplyr::n_distinct(drug_or_class), + drug_or_class_csv = paste(sort(unique(drug_or_class)), collapse = ", ") + ) |> + dplyr::filter(n_drug_or_class >= min_drugs_or_classes) |> + dplyr::ungroup() |> + dplyr::select(source, n_drug_or_class, drug_or_class_csv) |> + dplyr::arrange(dplyr::desc(n_drug_or_class)) + + return(shared_dyads) } + +#' Find the dyads that are unique to a single drug/class +#' +#' @param top_dyads The tibble of summarized dyads generated from `summarisedyads()` +#' @param label The \code{drug_label} value to filter dyads by, either \code{"drug"} or \code{"drug_class"} (default is \code{"drug"}) +#' @param protein_names_parquet The path to the Parquet file containing the annotations to protein dyad names +#' +#' @returns 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. +#' +#' @export +#' @examples +#' findUniquedyads(summarisedyads(top_features, dyad_feature_parquet), label = "drug", protein_names_parquet) +findUniqueDyads <- function(top_dyads = summariseDyads(top_features, dyad_feature_parquet), + label = "drug", + protein_names_parquet +) { + + protein_names <- arrow::read_parquet(normalizePath(protein_names_parquet)) |> + dplyr::distinct() + + unique_dyads <- top_dyads |> + dplyr::filter(!is.na(source), drug_label == label) |> + dplyr::group_by(source) |> + dplyr::mutate( + n_drug_or_class = dplyr::n_distinct(drug_or_class), + drug_or_class_csv = paste(sort(unique(drug_or_class)), collapse = ", ") + ) |> + dplyr::filter(n_drug_or_class == 1) |> + dplyr::ungroup() |> + dplyr::select(source, drug_or_class_csv, dyad_median_rank_score) |> + dplyr::arrange(dplyr::desc(dyad_median_rank_score)) |> +dplyr::rename(drug_or_class = drug_or_class_csv) + + return(unique_dyads) +} + +# 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, protein_names_parquet = protein_names_parquet) + # plotFeatureNetworkD3(feature_network) 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/Sfl_parquet.duckdb b/inst/extdata/Sfl_parquet.duckdb deleted file mode 100644 index 8a0ba86..0000000 Binary files a/inst/extdata/Sfl_parquet.duckdb 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/metadata.parquet b/inst/extdata/metadata.parquet deleted file mode 100644 index f5aa64e..0000000 Binary files a/inst/extdata/metadata.parquet and /dev/null 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/amRml-package.Rd b/man/amRml-package.Rd index 50cd609..b4e00b0 100644 --- a/man/amRml-package.Rd +++ b/man/amRml-package.Rd @@ -4,9 +4,9 @@ \name{amRml-package} \alias{amRml} \alias{amRml-package} -\title{amRml: Machile learning framework to predict AMR across molecular scales | pkg 2 of 3-pkg amR suite} +\title{amRml: Machine learning framework to predict AMR across molecular scales | pkg 2 of 3-pkg amR suite} \description{ -Comprehensive machine learning (ML) pipeline for predicting antimicrobial resistance (AMR) in ESKAPE pathogens across molecular scales (gene, protein, domain, structure). Implements logistic regression using tidymodels framework. Supports multi-drug resistance (MDR) and cross-resistance predictions with multiple holdout strategies (geographic, temporal, species). Includes model evaluation, feature importance analysis, cross-scale performance comparison, and functional characterization of top features. Designed to work with data from amRdata package and generate results for amRshiny visualization. Can work independently with example ESKAPE pathogen datasets, too. +Comprehensive machine learning (ML) pipeline for predicting antimicrobial resistance (AMR) in ESKAPE pathogens across molecular scales (gene, protein, domain, structure). Implements logistic regression using tidymodels framework. Supports multi-drug resistance (MDR) and cross-resistance predictions with multiple holdout strategies (geographic, temporal, species). Includes model evaluation, feature importance analysis, cross-scale performance comparison, and functional characterization of top features. Designed to work with data from amRdata package and generate results for amRviz visualization. Can work independently with example ESKAPE pathogen datasets, too. } \seealso{ Useful links: @@ -30,6 +30,9 @@ Other contributors: \item Evan Brenner \email{evan.brenner@cuanschutz.edu} [contributor] \item Abhirupa Ghosh \email{abhirupa.ghosh@cuanschutz.edu} [contributor] \item David Mayer \email{david.mayer@cuanschutz.edu} [contributor] + \item Emily Boyer \email{emily.boyer@cuanschutz.edu} [contributor] + \item Charmie Vang \email{charmie.vang@cuanschutz.edu} [contributor] + \item Alexander McKim \email{alexander.mckim@cuanschutz.edu} (\href{https://orcid.org/0000-0002-7802-7591}{ORCID}) [contributor] } } 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{":"} (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/buildPerfPq.Rd b/man/buildPerfPq.Rd new file mode 100644 index 0000000..2fb7278 --- /dev/null +++ b/man/buildPerfPq.Rd @@ -0,0 +1,45 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/merge_ml_results.R +\name{buildPerfPq} +\alias{buildPerfPq} +\title{Build Parquet file of ML performance results} +\usage{ +buildPerfPq( + path, + stratify_by = NULL, + LOO = FALSE, + MDR = FALSE, + cross_test = FALSE, + out_parquet = "all_perf.parquet", + compression = "zstd", + verbose = TRUE +) +} +\arguments{ +\item{path}{Base directory containing ML results} + +\item{stratify_by}{NULL, "year", or "country"} + +\item{LOO}{Logical; leave-one-out analysis} + +\item{MDR}{Logical; MDR mode} + +\item{cross_test}{Logical; cross-testing mode} + +\item{out_parquet}{Output file name} + +\item{compression}{Compression method (default: "zstd")} + +\item{verbose}{Logical; print progress} +} +\value{ +A tibble with metadata columns + performance metrics +} +\description{ +Reads all \verb{_performance.tsv} files, parses metadata from filenames, +combines them into a single table, and writes a Parquet output. +} +\examples{ +buildPerfPq("data/Campylobacter") + +} diff --git a/man/buildPerfPqCrossCountry.Rd b/man/buildPerfPqCrossCountry.Rd new file mode 100644 index 0000000..bf8ce85 --- /dev/null +++ b/man/buildPerfPqCrossCountry.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/merge_ml_results.R +\name{buildPerfPqCrossCountry} +\alias{buildPerfPqCrossCountry} +\title{Build Parquet file from cross country testing ML performances} +\usage{ +buildPerfPqCrossCountry(perf_dir_path) +} +\arguments{ +\item{perf_dir_path}{Directory containing performance TSV files} +} +\value{ +Writes a Parquet file to the same directory +} +\description{ +Build Parquet file from cross country testing ML performances +} +\examples{ +buildPerfPqCrossCountry(perf_dir_path = "data/Campylobacter/cross_test_ML_country_performance/") + +} diff --git a/man/buildPerfPqCrossDrug.Rd b/man/buildPerfPqCrossDrug.Rd new file mode 100644 index 0000000..b648834 --- /dev/null +++ b/man/buildPerfPqCrossDrug.Rd @@ -0,0 +1,22 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/merge_ml_results.R +\name{buildPerfPqCrossDrug} +\alias{buildPerfPqCrossDrug} +\title{Build Parquet file from cross drug testing ML performances} +\usage{ +buildPerfPqCrossDrug(perf_dir_path) +} +\arguments{ +\item{perf_dir_path}{Directory containing performance TSV files} +} +\value{ +Writes a Parquet file to the same directory +} +\description{ +Build Parquet file from cross drug testing ML performances +} +\examples{ +# Cross testing +buildPerfPqCrossDrug(perf_dir_path = "data/Campylobacter/cross_test_ML_performance/") + +} diff --git a/man/buildPerfPqCrossYear.Rd b/man/buildPerfPqCrossYear.Rd new file mode 100644 index 0000000..c679886 --- /dev/null +++ b/man/buildPerfPqCrossYear.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/merge_ml_results.R +\name{buildPerfPqCrossYear} +\alias{buildPerfPqCrossYear} +\title{Build Parquet file from cross year testing ML performances} +\usage{ +buildPerfPqCrossYear(perf_dir_path) +} +\arguments{ +\item{perf_dir_path}{Directory containing performance TSV files} +} +\value{ +Writes a Parquet file to the same directory +} +\description{ +Build Parquet file from cross year testing ML performances +} +\examples{ +buildPerfPqCrossYear(perf_dir_path = "data/Campylobacter/cross_test_ML_year_performance/") + +} diff --git a/man/buildPerfPqLOOCountry.Rd b/man/buildPerfPqLOOCountry.Rd new file mode 100644 index 0000000..869df95 --- /dev/null +++ b/man/buildPerfPqLOOCountry.Rd @@ -0,0 +1,22 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/merge_ml_results.R +\name{buildPerfPqLOOCountry} +\alias{buildPerfPqLOOCountry} +\title{Build Parquet file from LOO country ML performances} +\usage{ +buildPerfPqLOOCountry(perf_dir_path) +} +\arguments{ +\item{perf_dir_path}{Directory containing performance TSV files} +} +\value{ +Writes a Parquet file to the same directory +} +\description{ +Build Parquet file from LOO country ML performances +} +\examples{ +# LOO country +buildPerfPqLOOCountry(perf_dir_path = "data/Campylobacter/LOO_ML_country_performance/") + +} diff --git a/man/buildPerfPqLOODrug.Rd b/man/buildPerfPqLOODrug.Rd new file mode 100644 index 0000000..a7943bd --- /dev/null +++ b/man/buildPerfPqLOODrug.Rd @@ -0,0 +1,22 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/merge_ml_results.R +\name{buildPerfPqLOODrug} +\alias{buildPerfPqLOODrug} +\title{Build Parquet file from LOO drug ML performances} +\usage{ +buildPerfPqLOODrug(perf_dir_path) +} +\arguments{ +\item{perf_dir_path}{Directory containing performance TSV files} +} +\value{ +Writes a Parquet file to the same directory +} +\description{ +Build Parquet file from LOO drug ML performances +} +\examples{ +# LOO drug +buildPerfPqLOODrug(perf_dir_path = "data/Campylobacter/LOO_ML_performance/") + +} diff --git a/man/buildPerfPqLOOYear.Rd b/man/buildPerfPqLOOYear.Rd new file mode 100644 index 0000000..b9abc6b --- /dev/null +++ b/man/buildPerfPqLOOYear.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/merge_ml_results.R +\name{buildPerfPqLOOYear} +\alias{buildPerfPqLOOYear} +\title{Build Parquet file from LOO year ML performances} +\usage{ +buildPerfPqLOOYear(perf_dir_path) +} +\arguments{ +\item{perf_dir_path}{Directory containing performance TSV files} + +\item{top_feat_dir_path}{Directory containing top feature TSV files} +} +\value{ +Writes a Parquet file to the same directory +} +\description{ +Build Parquet file from LOO year ML performances +} +\examples{ +# LOO year +buildPerfPqLOOYear(perf_dir_path = "data/Campylobacter/LOO_ML_year_performance/") + +} diff --git a/man/buildPerfPqMDR.Rd b/man/buildPerfPqMDR.Rd new file mode 100644 index 0000000..42b14c8 --- /dev/null +++ b/man/buildPerfPqMDR.Rd @@ -0,0 +1,22 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/merge_ml_results.R +\name{buildPerfPqMDR} +\alias{buildPerfPqMDR} +\title{Build Parquet file from MDR ML performances} +\usage{ +buildPerfPqMDR(perf_dir_path) +} +\arguments{ +\item{perf_dir_path}{Directory containing performance TSV files} +} +\value{ +Writes a Parquet file to the same directory +} +\description{ +Build Parquet file from MDR ML performances +} +\examples{ +# MDR +buildPerfPqMDR(perf_dir_path = "data/Campylobacter/MDR_ML_performance/") + +} diff --git a/man/buildPerfPqYearCountry.Rd b/man/buildPerfPqYearCountry.Rd new file mode 100644 index 0000000..b0444e3 --- /dev/null +++ b/man/buildPerfPqYearCountry.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/merge_ml_results.R +\name{buildPerfPqYearCountry} +\alias{buildPerfPqYearCountry} +\title{Build Parquet file from year/country ML performances} +\usage{ +buildPerfPqYearCountry(perf_dir_path) +} +\arguments{ +\item{perf_dir_path}{Directory containing performance TSV files} +} +\value{ +Writes a Parquet file to the same directoryYear stratified + +buildPerfPqYearCountry(perf_dir_path = "data/Campylobacter/ML_year_performance/")Country stratified + +buildPerfPqYearCountry(perf_dir_path = "data/Campylobacter/ML_country_performance/") +} +\description{ +Build Parquet file from year/country ML performances +} diff --git a/man/buildPredPqMDR.Rd b/man/buildPredPqMDR.Rd new file mode 100644 index 0000000..a7c3103 --- /dev/null +++ b/man/buildPredPqMDR.Rd @@ -0,0 +1,22 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/merge_ml_results.R +\name{buildPredPqMDR} +\alias{buildPredPqMDR} +\title{Build Parquet file from MDR ML predictions} +\usage{ +buildPredPqMDR(pred_dir_path) +} +\arguments{ +\item{perf_dir_path}{Directory containing prediction TSV files} +} +\value{ +Writes a Parquet file to the same directory +} +\description{ +Build Parquet file from MDR ML predictions +} +\examples{ +# MDR +buildPredPqMDR(pred_dir_path = "data/Campylobacter/MDR_ML_pred/") + +} diff --git a/man/buildTopFeatsPq.Rd b/man/buildTopFeatsPq.Rd new file mode 100644 index 0000000..1d99a7c --- /dev/null +++ b/man/buildTopFeatsPq.Rd @@ -0,0 +1,45 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/merge_ml_results.R +\name{buildTopFeatsPq} +\alias{buildTopFeatsPq} +\title{Build Parquet file of top ML features} +\usage{ +buildTopFeatsPq( + path, + stratify_by = NULL, + LOO = FALSE, + MDR = FALSE, + cross_test = FALSE, + out_parquet = "all_top_features.parquet", + compression = "zstd", + verbose = TRUE +) +} +\arguments{ +\item{path}{Base directory containing ML results} + +\item{stratify_by}{NULL, "year", or "country"} + +\item{LOO}{Logical; leave-one-out analysis} + +\item{MDR}{Logical; MDR mode} + +\item{cross_test}{Logical; cross-testing mode} + +\item{out_parquet}{Output file name} + +\item{compression}{Compression method (default: "zstd")} + +\item{verbose}{Logical; print progress} +} +\value{ +A tibble with metadata columns + feature importance data +} +\description{ +Reads all \verb{_top_features.tsv} files, parses metadata from filenames, +combines them into a single table, and writes a Parquet file. +} +\examples{ +buildTopFeatsPq("data/Campylobacter") + +} diff --git a/man/buildTopFeatsPqLOOCountry.Rd b/man/buildTopFeatsPqLOOCountry.Rd new file mode 100644 index 0000000..882b469 --- /dev/null +++ b/man/buildTopFeatsPqLOOCountry.Rd @@ -0,0 +1,23 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/merge_ml_results.R +\name{buildTopFeatsPqLOOCountry} +\alias{buildTopFeatsPqLOOCountry} +\title{Build Parquet file from LOO country ML top features} +\usage{ +buildTopFeatsPqLOOCountry(top_feat_dir_path) +} +\arguments{ +\item{top_feat_dir_path}{Directory containing top feature TSV files} +} +\value{ +Writes a Parquet file to the same directory +} +\description{ +Build Parquet file from LOO country ML top features +} +\examples{ + +# LOO country +buildTopFeatsPqLOOCountry(top_feat_dir_path = "data/Campylobacter/LOO_ML_country_top_features/") + +} diff --git a/man/buildTopFeatsPqLOODrug.Rd b/man/buildTopFeatsPqLOODrug.Rd new file mode 100644 index 0000000..9b62f5c --- /dev/null +++ b/man/buildTopFeatsPqLOODrug.Rd @@ -0,0 +1,22 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/merge_ml_results.R +\name{buildTopFeatsPqLOODrug} +\alias{buildTopFeatsPqLOODrug} +\title{Build Parquet file from LOO drug ML top features} +\usage{ +buildTopFeatsPqLOODrug(top_feat_dir_path) +} +\arguments{ +\item{top_feat_dir_path}{Directory containing top feature TSV files} +} +\value{ +Writes a Parquet file to the same directory +} +\description{ +Build Parquet file from LOO drug ML top features +} +\examples{ +# LOO drug +buildTopFeatsPqLOODrug(top_feat_dir_path = "data/Campylobacter/LOO_ML_top_features/") + +} diff --git a/man/buildTopFeatsPqLOOYear.Rd b/man/buildTopFeatsPqLOOYear.Rd new file mode 100644 index 0000000..e3fca09 --- /dev/null +++ b/man/buildTopFeatsPqLOOYear.Rd @@ -0,0 +1,23 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/merge_ml_results.R +\name{buildTopFeatsPqLOOYear} +\alias{buildTopFeatsPqLOOYear} +\title{Build Parquet file from LOO year ML top features} +\usage{ +buildTopFeatsPqLOOYear(top_feat_dir_path) +} +\arguments{ +\item{top_feat_dir_path}{Directory containing top feature TSV files} +} +\value{ +Writes a Parquet file to the same directory +} +\description{ +Build Parquet file from LOO year ML top features +} +\examples{ + +# LOO year +buildTopFeatsPqLOOYear(top_feat_dir_path = "data/Campylobacter/LOO_ML_year_top_features/") + +} diff --git a/man/buildTopFeatsPqMDR.Rd b/man/buildTopFeatsPqMDR.Rd new file mode 100644 index 0000000..69d946c --- /dev/null +++ b/man/buildTopFeatsPqMDR.Rd @@ -0,0 +1,22 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/merge_ml_results.R +\name{buildTopFeatsPqMDR} +\alias{buildTopFeatsPqMDR} +\title{Build Parquet file from MDR ML top features} +\usage{ +buildTopFeatsPqMDR(top_feat_dir_path) +} +\arguments{ +\item{top_feat_dir_path}{Directory containing top feature TSV files} +} +\value{ +Writes a Parquet file to the same directory +} +\description{ +Build Parquet file from MDR ML top features +} +\examples{ +# MDR +buildTopFeatsPqMDR(top_feat_dir_path = "data/Campylobacter/MDR_ML_top_features/") + +} diff --git a/man/buildTopFeatsPqYearCountry.Rd b/man/buildTopFeatsPqYearCountry.Rd new file mode 100644 index 0000000..0b8d4ed --- /dev/null +++ b/man/buildTopFeatsPqYearCountry.Rd @@ -0,0 +1,25 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/merge_ml_results.R +\name{buildTopFeatsPqYearCountry} +\alias{buildTopFeatsPqYearCountry} +\title{Build Parquet file from year/country ML top features} +\usage{ +buildTopFeatsPqYearCountry(top_feat_dir_path) +} +\arguments{ +\item{top_feat_dir_path}{Directory containing top feature TSV files} +} +\value{ +Writes a Parquet file to the same directory +} +\description{ +Build Parquet file from year/country ML top features +} +\examples{ +# Year stratified +buildTopFeatsPqYearCountry(top_feat_dir_path = "data/Campylobacter/ML_year_top_features/") + +# Country stratified +buildTopFeatsPqYearCountry(top_feat_dir_path = "data/Campylobacter/ML_country_top_features/") + +} 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/parse_ml_filename.Rd b/man/parse_ml_filename.Rd new file mode 100644 index 0000000..8bd6ec2 --- /dev/null +++ b/man/parse_ml_filename.Rd @@ -0,0 +1,34 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/merge_ml_results.R +\name{parse_ml_filename} +\alias{parse_ml_filename} +\title{Parse ML result filename into structured metadata} +\usage{ +parse_ml_filename(filename) +} +\arguments{ +\item{filename}{Character. Filename (not full path) to parse.} +} +\value{ +A named list with elements: +\describe{ +\item{shuffled}{Logical; whether dataset was shuffled} +\item{species}{Species name} +\item{drug_label}{"drug" or "drug_class"} +\item{drug_or_class}{Drug name or class} +\item{strat_label}{Stratification type (year/country) or NA} +\item{strat_value}{Stratification value or NA} +\item{feature_type}{Feature category} +\item{feature_subtype}{binary/counts} +\item{seed}{Integer seed} +} +} +\description{ +Extracts metadata fields encoded in machine learning result filenames. +Supports shuffled runs, drug vs drug_class, optional stratification, +feature types, and seed information. +} +\examples{ +parse_ml_filename("Csp_drug_AMX_genes_binary_42_top_features.tsv") + +} 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..c1b13b1 --- /dev/null +++ b/man/scoreFeaturesWithinSeed.Rd @@ -0,0 +1,57 @@ +% 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, + add_lasso_advtg = 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 (regardless of \code{filter_model} or \code{add_lasso_advtg}) to compute the per-fit sparsity score.} + +\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{add_lasso_advtg}{Logical indicating whether to weight each feature's contribution by a sparsity score that rewards model fits returning fewer features relative to the candidate feature space (default is \code{FALSE})} +} +\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{sparsity_score}: \code{1 - feat_return_ratio} when \code{add_lasso_advtg = TRUE} (fits returning fewer features relative to the candidate space score closer to 1), else 1 for every row. +\item \code{adjusted_contribution}: \code{contribution * sparsity_score}. This, not raw \code{contribution}, is what every downstream column below is actually computed from. +\item \code{rank}: descending rank of \code{adjusted_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{adjusted_contribution} in descending order; features tied on \code{adjusted_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{":"} (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{":"}) 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{":"}) 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..c77e2f6 --- /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{adjusted_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..7727cbe --- /dev/null +++ b/man/topFeaturesPerDrugOrClass.Rd @@ -0,0 +1,72 @@ +% 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, + add_lasso_advtg = FALSE, + 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 (regardless of \code{filter_model} or \code{add_lasso_advtg}) to compute the per-fit sparsity score.} + +\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{add_lasso_advtg}{Logical indicating whether to weight each feature's contribution by a sparsity score that rewards model fits returning fewer features relative to the candidate feature space (default is \code{FALSE})} + +\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}. The quantile is computed once over every row of \code{feature_summary} (with \code{na.rm = TRUE}) — globally across all species/drugs/drug classes/feature types/subtypes, not per group — and is not restricted to rows that already pass the other conditions listed below: \code{dplyr::filter()} evaluates every condition passed to a single call against the same original, ungrouped 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 the full, ungrouped \code{feature_summary} (see \code{rank_score_quantile} for what that means for the last condition): +\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}. +} +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. +} +\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. +}