From de2fc8672f4162245dd60866e3cdc546badd598b Mon Sep 17 00:00:00 2001 From: Abhirupa Ghosh <100681585+AbhirupaGhosh@users.noreply.github.com> Date: Tue, 20 Jan 2026 11:24:24 -0700 Subject: [PATCH 01/17] Enhance plotting functions and documentation Updated documentation for plotPRC, plotROC, plotCM, plotDensity, and plotTopFeatsVI functions. Added new functions for plotting ROC curves, confusion matrices, and density of predicted class probabilities. --- R/plot_ml.R | 258 ++++++++++++++++++++++++++++++++-------------------- 1 file changed, 157 insertions(+), 101 deletions(-) diff --git a/R/plot_ml.R b/R/plot_ml.R index 6e3eb40..088c481 100644 --- a/R/plot_ml.R +++ b/R/plot_ml.R @@ -19,20 +19,36 @@ #' @importFrom tune extract_fit_parsnip #' @importFrom vip vip #' @importFrom yardstick pr_curve -#' @importFrom graphics barplot NULL -#' plotPRC() -#' -#' Plots the precision-recall curve given a set of test data plus predicted AMR -#' phenotypes. -#' -#' @param test_data_plus_predictions Test data (tibble) with an added column for -#' predicted phenotype labels, such as the output of `predict()`. -#' @return A precision-recall curve as a `ggplot2` object -#' @export +#' Plot a Precision-Recall Curve +#' +#' Generates a precision-recall curve (PRC) for AMR phenotype prediction results. +#' @param test_data_plus_predictions A tibble containing test data with added +#' prediction columns, typically the output of `runMLmodels()`. +#' +#' @return A `ggplot2` object showing the precision-recall curve. +#' +#' @details +#' The function uses `yardstick::pr_curve()` to compute the PR curve and then +#' visualizes it using `ggplot2`. +#' +#' @examples +#' \dontrun{ +#' test_data_plus_predictions <- readr::read_tsv(results/ML_pred/Sfl_drug_AMP_domains_binary_prediction.tsv) +#' plotPRC(test_data_plus_predictions) +#' } +#' +#' @export plotPRC <- function(test_data_plus_predictions) { .checkArgTestDataPlusPredictions(test_data_plus_predictions) +test_data_plus_predictions <- test_data_plus_predictions |> +dplyr::mutate( +genome_drug.resistant_phenotype = factor( +genome_drug.resistant_phenotype, +levels = c("Resistant", "Susceptible") +) +) prc <- yardstick::pr_curve( test_data_plus_predictions, @@ -46,113 +62,154 @@ plotPRC <- function(test_data_plus_predictions) { return(prc) } -#' plotTopFeatsVI() -#' -#' Generates a plot showing the top features and their variable importance -#' scores. -#' -#' @param fit Best model fit, such as the output of `fitBestModel()` -#' @param n_top_feats [num] Number of top features to plot -#' @return Variable importance plot (a `ggplot2` object) +#' Plot a Receiver Operating Characteristic (ROC) Curve +#' +#' Generates a ROC curve for AMR phenotype prediction results. +#' +#' @param test_data_plus_predictions A tibble with test data and prediction +#' columns (output of `runMLmodels()`). +#' +#' @return A ROC curve plotted using `ggplot2::autoplot()`. +#' #' @export -plotTopFeatsVI <- function(fit, n_top_feats = 10) { - .checkArgWflow(fit) - .checkArgNTopFeats(n_top_feats) - - vip <- fit |> - tune::extract_fit_parsnip() |> - vip::vip(num_features = n_top_feats) + - ggplot2::xlab("Top Features") + +plotROC <- function(test_data_plus_predictions) { + .checkArgTestDataPlusPredictions(test_data_plus_predictions) + test_data_plus_predictions <- test_data_plus_predictions |> +dplyr::mutate( +genome_drug.resistant_phenotype = factor( +genome_drug.resistant_phenotype, +levels = c("Resistant", "Susceptible") +) +) + + roc <- yardstick::roc_curve( + test_data_plus_predictions, + genome_drug.resistant_phenotype, .pred_Resistant + ) |> + ggplot2::autoplot(type = "se") + ggplot2::theme(panel.grid = ggplot2::element_blank()) - - return(vip) + + return(roc) } -#' plotDefaultEval() -#' -#' Plots performance metric or runtime vs. training data proportion or number -#' of cross-validation folds, colored by model. -#' -#' @param default_eval_tibble Output of `findOptimalMLDefaults()` -#' @param x_default_eval [chr] x value of default evaluation plot: "train_prop" -#' or "n_fold" -#' @param y_default_eval [chr] y value of default evaluation plot. It can be -#' "avg_runtime_sec" or one of the following performance metrics: -#' "avg_f1_score", "avg_log2_apop", "avg_bal_acc", or "avg_nmcc" -#' @param xlab [chr] Label for x axis -#' @param ylab [chr] Label for y axis -#' @return A `ggplot2` scatterplot (performance metric or runtime vs. -#' `train_prop` or `n_fold`), colored by model +#' Plot a Confusion Matrix Heatmap +#' +#' Produces a heatmap visualization of the confusion matrix for AMR predictions. +#' +#' @param test_data_plus_predictions A tibble containing true and predicted +#' phenotype labels. +#' +#' @return A heatmap (`ggplot2` object) showing the confusion matrix. +#' #' @export -plotDefaultEval <- function( - default_eval_tibble, x_default_eval = "train_prop", - y_default_eval = "avg_f1_score", xlab = "Train Data Proportion", - ylab = "Average F1 Score" -) { - .checkArgTibble(default_eval_tibble) - .checkArgXDefaultEval(x_default_eval) - .checkArgYDefaultEval(y_default_eval) - .checkArgXYLabs(xlab = xlab, ylab = ylab) - - if (x_default_eval == "n_fold") { - default_eval_tibble <- default_eval_tibble |> - dplyr::filter(train_prop == 0.8) - } else { - default_eval_tibble <- default_eval_tibble |> - dplyr::filter(train_prop != 0.8) - } - - default_eval_plot <- ggplot2::ggplot( - default_eval_tibble, - ggplot2::aes( - x = unlist(default_eval_tibble[x_default_eval]), - y = unlist(default_eval_tibble[y_default_eval]), color = model +plotCM <- function(test_data_plus_predictions) { + .checkArgTestDataPlusPredictions(test_data_plus_predictions) + test_data_plus_predictions <- test_data_plus_predictions |> + dplyr::mutate( + genome_drug.resistant_phenotype = factor( + genome_drug.resistant_phenotype, + levels = c("Resistant", "Susceptible") + ), + .pred_class = factor( + .pred_class, + levels = c("Resistant", "Susceptible") + ) ) - ) + - ggplot2::geom_line(size = 1.5) + - ggplot2::geom_point(size = 3) + - ggplot2::theme( - axis.line = ggplot2::element_line(linewidth = 1.5), - axis.ticks = ggplot2::element_line(linewidth = 1.5, colour = "black"), - axis.text = ggplot2::element_text(size = 16, colour = "black"), - axis.title = ggplot2::element_text(size = 16, face = "bold"), - panel.grid = ggplot2::element_blank(), - panel.background = ggplot2::element_blank(), - legend.text = ggplot2::element_text(size = 16), - legend.title = ggplot2::element_text(size = 16, face = "bold") - ) + - ggplot2::labs(x = xlab, y = ylab, color = "Model") + test_data_plus_predictions |> +yardstick::conf_mat(truth = genome_drug.resistant_phenotype, + estimate = .pred_class) |> +ggplot2::autoplot(type = "heatmap") +} - return(default_eval_plot) +#' Plot Density of Predicted Class Probabilities +#' +#' Visualizes how predicted class probabilities differ between resistant and +#' susceptible genome-drug combinations. +#' +#' @param test_data_plus_predictions Tibble with prediction probabilities and +#' true labels. +#' +#' @return A ggplot2 density plot. +#' +#' @export +plotDensity <- function(test_data_plus_predictions) { + test_data_plus_predictions |> +ggplot2::ggplot(ggplot2::aes(x = .pred_Resistant, +fill = genome_drug.resistant_phenotype)) + +ggplot2::geom_density(alpha = 0.5) +} + +#' Plot Top Feature Importances +#' +#' Creates a bar plot showing the most important features affecting +#' AMR phenotype predictions. +#' +#' @param topfeat A tibble containing feature importance scores +#' (output of `runMLmodels()`). +#' @param n_top_feats Number of top features to display (default: 10). +#' +#' @return A bar plot of variable importance (`ggplot2` object). +#' +#' @examples +#' \dontrun{ +#' topfeat <- readr::read_tsv(results/ML_top_features/Sfl_drug_AMP_domains_binary_top_features.tsv) +#' plotTopFeatsVI(topfeat) +#' } +#' +#' @export +plotTopFeatsVI <- function(topfeat, n_top_feats = 10) { + .checkArgNTopFeats(n_top_feats) + + vip <- topfeat |> + dplyr::slice_max(order_by = Importance, n = n_top_feats) |> + dplyr::mutate( + Variable = factor(Variable, levels = rev(Variable)), # preserve order as shown in table + Sign = factor(Sign, levels = c("POS", "NEG")) + ) |> + ggplot2::ggplot(ggplot2::aes(x = Importance, y = Variable, fill = Sign)) + + ggplot2::geom_col() + + ggplot2::scale_fill_manual( + values = c( + "POS" = "#c6d8d3", + "NEG" = "#f6c9a1" + ) + ) + + ggplot2::labs( + x = "Importance", + y = "Features" + ) + + ggplot2::theme_minimal(base_size = 14) + + ggplot2::theme( + panel.grid.minor = ggplot2::element_blank(), + axis.text.y = ggplot2::element_text(size = 10) + ) + + return(vip) } -#' getBaselineComparisonBarplot() -#' -#' Generates a bar plot that compares model performance with and without -#' randomly shuffled AMR phenotype labels. -#' -#' @param non_shuffled_label_results Output of `runMLPipeline()` -#' (`shuffle_labels = FALSE`) -#' @param shuffled_label_results Output of `runMLPipeline()` -#' (`shuffle_labels = TRUE`) -#' @return A bar plot with balanced accuracy comparisons per antibiotic +#' Compare Baseline Performance With and Without Shuffled Labels +#' +#' Produces a bar plot comparing balanced accuracy for each antibiotic using +#' true AMR labels vs. randomly shuffled labels. +#' +#' @param non_shuffled_label_results Output of `runMLPipeline(shuffle_labels = FALSE)` +#' @param shuffled_label_results Output of `runMLPipeline(shuffle_labels = TRUE)` +#' +#' @return A base R barplot comparing balanced accuracy across models. +#' #' @export getBaselineComparisonBarplot <- function( non_shuffled_label_results, shuffled_label_results ) { - .checkArgTibble(non_shuffled_label_results) - .checkArgTibble(shuffled_label_results) - - drugs <- non_shuffled_label_results |> - dplyr::select(antibiotic) |> - dplyr::pull() + .checkArgTibble(non_shuffled_label_results$performance_tibble) + .checkArgTibble(shuffled_label_results$performance_tibble) - non_shuffled_bal_acc <- non_shuffled_label_results |> + non_shuffled_bal_acc <- non_shuffled_label_results$performance_tibble |> dplyr::select(bal_acc) |> dplyr::pull() - shuffled_bal_acc <- shuffled_label_results |> + shuffled_bal_acc <- shuffled_label_results$performance_tibble |> dplyr::select(bal_acc) |> dplyr::pull() @@ -160,13 +217,12 @@ getBaselineComparisonBarplot <- function( nrow = 2, byrow = TRUE ) - colnames(bal_acc_matrix) <- drugs rownames(bal_acc_matrix) <- c("Non-Shuffled Labels", "Shuffled Labels") baseline_comparison_barplot <- barplot(bal_acc_matrix, beside = TRUE, legend.text = TRUE, col = c("skyblue", "lightpink"), - ylab = "Balanced Accuracy", xlab = "Antibiotic" + ylab = "Balanced Accuracy" ) return(baseline_comparison_barplot) From 9486c555d0933cbeb340b2ece9dce4dc637e9b63 Mon Sep 17 00:00:00 2001 From: AbhirupaGhosh Date: Wed, 21 Jan 2026 06:33:43 +0000 Subject: [PATCH 02/17] Style code (GHA) --- R/core_ml.R | 217 ++++++++++++------- R/generate_matrices_ml.R | 197 ++++++++++------- R/globals.R | 2 - R/plot_ml.R | 121 ++++++----- R/prep_ml.R | 6 +- R/run_ML.R | 450 +++++++++++++++++++++------------------ R/run_ml_pipeline.R | 41 ++-- vignettes/intro.Rmd | 30 +-- 8 files changed, 605 insertions(+), 459 deletions(-) diff --git a/R/core_ml.R b/R/core_ml.R index db874db..1fb0c1b 100644 --- a/R/core_ml.R +++ b/R/core_ml.R @@ -73,7 +73,8 @@ NULL #' @return An `rsplit` object #' @export splitMLInputTibble <- function(ml_input_tibble, split = c(0.6, 0.2), seed = 5280) { - .checkArgTibble(ml_input_tibble, ml = TRUE); .checkArgSplit(split) + .checkArgTibble(ml_input_tibble, ml = TRUE) + .checkArgSplit(split) .checkArgSeed(seed) set.seed(seed) @@ -85,7 +86,7 @@ splitMLInputTibble <- function(ml_input_tibble, split = c(0.6, 0.2), seed = 5280 # If in CV mode: # Still retain a stratified testing holdout purely for final reporting metrics; # CV is only performed on the training portion. - prop_train_for_holdout <- 0.8 # 80 percent train, 20 percent reserved test + prop_train_for_holdout <- 0.8 # 80 percent train, 20 percent reserved test data_split <- rsample::initial_split( ml_input_tibble, prop = prop_train_for_holdout, @@ -115,7 +116,8 @@ splitMLInputTibble <- function(ml_input_tibble, split = c(0.6, 0.2), seed = 5280 #' @return A `recipe` object #' @export buildRecipe <- function(train_data, use_pca = FALSE, pca_threshold = 0.95) { - .checkArgTibble(train_data, ml = TRUE); .checkArgUsePCA(use_pca) + .checkArgTibble(train_data, ml = TRUE) + .checkArgUsePCA(use_pca) .checkArgPCAThreshold(pca_threshold) target_var <- .getTargetVarName(train_data) |> as.character() @@ -124,8 +126,10 @@ buildRecipe <- function(train_data, use_pca = FALSE, pca_threshold = 0.95) { nm <- names(train_data) id_cols <- setdiff(nm[grepl("^genome", nm)], target_var) - rec <- recipes::recipe(formula = stats::reformulate(".", response = target_var), - data = train_data) + rec <- recipes::recipe( + formula = stats::reformulate(".", response = target_var), + data = train_data + ) # Only update roles if we actually have ID columns to mark as metadata if (length(id_cols) > 0) { @@ -146,7 +150,6 @@ buildRecipe <- function(train_data, use_pca = FALSE, pca_threshold = 0.95) { } - #' buildLRModel() #' #' Builds a logistic regression model. @@ -158,13 +161,17 @@ buildRecipe <- function(train_data, use_pca = FALSE, pca_threshold = 0.95) { buildLRModel <- function(multi_class = FALSE) { .checkArgMultiClass(multi_class) - if(!multi_class) { - lr_mod <- parsnip::logistic_reg(penalty = hardhat::tune(), - mixture = hardhat::tune()) |> + if (!multi_class) { + lr_mod <- parsnip::logistic_reg( + penalty = hardhat::tune(), + mixture = hardhat::tune() + ) |> parsnip::set_engine(engine = "glmnet") - } else if(multi_class) { - lr_mod <- parsnip::multinom_reg(penalty = hardhat::tune(), - mixture = hardhat::tune()) |> + } else if (multi_class) { + lr_mod <- parsnip::multinom_reg( + penalty = hardhat::tune(), + mixture = hardhat::tune() + ) |> parsnip::set_engine(engine = "glmnet") } @@ -181,9 +188,11 @@ buildLRModel <- function(multi_class = FALSE) { #' @return A `workflow` object #' @export buildWflow <- function(parsnip_mod, recipe) { - .checkArgParsnipMod(parsnip_mod); .checkArgRecipe(recipe) + .checkArgParsnipMod(parsnip_mod) + .checkArgRecipe(recipe) - wflow <- workflows::workflow() |> workflows::add_model(parsnip_mod) |> + wflow <- workflows::workflow() |> + workflows::add_model(parsnip_mod) |> workflows::add_recipe(recipe) return(wflow) @@ -203,21 +212,21 @@ buildWflow <- function(parsnip_mod, recipe) { #' @return A logistic regression tuning grid as a tibble #' @export buildTuningGrid <- function( - model = "LR", - penalty_vec = 10^seq(-4, -1, length.out = 10), - mix_vec = 0:5 / 5 + model = "LR", + penalty_vec = 10^seq(-4, -1, length.out = 10), + mix_vec = 0:5 / 5 ) { .checkArgModel(model) - + if (model == "LR") { .checkArgPenaltyVec(penalty_vec) .checkArgMixVec(mix_vec) - + penalty <- rep(penalty_vec, each = length(mix_vec)) mixture <- rep(mix_vec, length(penalty_vec)) grid <- tibble::tibble(penalty, mixture) } - + return(grid) } @@ -237,13 +246,14 @@ buildTuningGrid <- function( #' @export tuneGrid <- function(wflow, data_split, grid = buildTuningGrid(model = "LR"), n_fold = 5) { - .checkArgTibble(grid); .checkArgWflow(wflow) + .checkArgTibble(grid) + .checkArgWflow(wflow) .checkArgDataSplit(data_split) split_class <- class(data_split)[1] # Always do CV on the training portion of the split - train_df <- rsample::training(data_split) + train_df <- rsample::training(data_split) target_var <- .getTargetVarName(train_df) if (identical(split_class, "initial_split")) { @@ -259,9 +269,9 @@ tuneGrid <- function(wflow, data_split, grid = buildTuningGrid(model = "LR"), tune_res <- tune::tune_grid( wflow, resamples = resamples, - grid = grid, - control = tune::control_grid(save_pred = TRUE), - metrics = yardstick::metric_set( + grid = grid, + control = tune::control_grid(save_pred = TRUE), + metrics = yardstick::metric_set( yardstick::f_meas, yardstick::pr_auc, yardstick::spec, @@ -286,7 +296,8 @@ tuneGrid <- function(wflow, data_split, grid = buildTuningGrid(model = "LR"), #' @return Best model workflow #' @export selectBestModel <- function(tune_res, wflow, select_best_metric = "mcc") { - .checkArgTuneRes(tune_res); .checkArgWflow(wflow) + .checkArgTuneRes(tune_res) + .checkArgWflow(wflow) .checkArgSelectBestMetric(select_best_metric) best_mod <- tune::select_best(tune_res, metric = select_best_metric) @@ -306,7 +317,8 @@ selectBestModel <- function(tune_res, wflow, select_best_metric = "mcc") { #' @return Best model fit #' @export fitBestModel <- function(final_mod, train_data) { - .checkArgWflow(final_mod); .checkArgTibble(train_data, ml = TRUE) + .checkArgWflow(final_mod) + .checkArgTibble(train_data, ml = TRUE) fit <- final_mod |> parsnip::fit(data = train_data) @@ -324,8 +336,7 @@ fitBestModel <- function(final_mod, train_data) { model <- class(fit$fit$actions$model$spec)[1] - if(model %in% c("logistic_reg", "multinom_reg")) { - + if (model %in% c("logistic_reg", "multinom_reg")) { penalty <- fit$fit$fit$spec$args$penalty mixture <- tryCatch( @@ -334,7 +345,6 @@ fitBestModel <- function(final_mod, train_data) { ) tibble::tibble(penalty = penalty, mixture = mixture) - } else { stop("The `fit` object provided must correspond to 'logistic_reg' or 'multinom_reg'.") } @@ -353,7 +363,8 @@ fitBestModel <- function(final_mod, train_data) { #' labels #' @export predictML <- function(fit, test_data) { - .checkArgWflow(fit); .checkArgTibble(test_data, ml = TRUE) + .checkArgWflow(fit) + .checkArgTibble(test_data, ml = TRUE) test_data_plus_predictions <- parsnip::augment(fit, test_data) @@ -396,7 +407,8 @@ getConfusionMatrix <- function(test_data_plus_predictions) { mcc <- test_data_plus_predictions |> yardstick::mcc(truth = !!target_var, estimate = .pred_class) |> - dplyr::select(.estimate) |> as.numeric() + dplyr::select(.estimate) |> + as.numeric() nmcc <- (mcc + 1) / 2 @@ -413,15 +425,21 @@ getConfusionMatrix <- function(test_data_plus_predictions) { .calculateF1 <- function(test_data_plus_predictions) { .checkArgTestDataPlusPredictions(test_data_plus_predictions) - if(!("genome_drug.resistant_phenotype" %in% + if (!("genome_drug.resistant_phenotype" %in% colnames(test_data_plus_predictions))) { - stop(paste("`test_data_plus_predictions` does not have a column for", - "`genome_drug.resistant_phenotype`.")) + stop(paste( + "`test_data_plus_predictions` does not have a column for", + "`genome_drug.resistant_phenotype`." + )) } f1 <- test_data_plus_predictions |> - yardstick::f_meas(truth = genome_drug.resistant_phenotype, - estimate = .pred_class) |> dplyr::select(.estimate) |> as.numeric() |> + yardstick::f_meas( + truth = genome_drug.resistant_phenotype, + estimate = .pred_class + ) |> + dplyr::select(.estimate) |> + as.numeric() |> round(2) return(f1) @@ -437,16 +455,21 @@ getConfusionMatrix <- function(test_data_plus_predictions) { .calculateAUPRC <- function(test_data_plus_predictions) { .checkArgTestDataPlusPredictions(test_data_plus_predictions) - if(!("genome_drug.resistant_phenotype" %in% + if (!("genome_drug.resistant_phenotype" %in% colnames(test_data_plus_predictions))) { - stop(paste("`test_data_plus_predictions` does not have a column for", - "`genome_drug.resistant_phenotype`.")) + stop(paste( + "`test_data_plus_predictions` does not have a column for", + "`genome_drug.resistant_phenotype`." + )) } auprc <- test_data_plus_predictions |> yardstick::pr_auc( - truth = genome_drug.resistant_phenotype, .pred_Resistant) |> - dplyr::select(.estimate) |> as.numeric() |> round(2) + truth = genome_drug.resistant_phenotype, .pred_Resistant + ) |> + dplyr::select(.estimate) |> + as.numeric() |> + round(2) return(auprc) } @@ -461,26 +484,33 @@ getConfusionMatrix <- function(test_data_plus_predictions) { .calculateLog2APOP <- function(test_data_plus_predictions) { .checkArgTestDataPlusPredictions(test_data_plus_predictions) - if(!("genome_drug.resistant_phenotype" %in% + if (!("genome_drug.resistant_phenotype" %in% colnames(test_data_plus_predictions))) { - stop(paste("`test_data_plus_predictions` does not have a column for", - "`genome_drug.resistant_phenotype`.")) + stop(paste( + "`test_data_plus_predictions` does not have a column for", + "`genome_drug.resistant_phenotype`." + )) } auprc <- .calculateAUPRC(test_data_plus_predictions) prior <- sum( - test_data_plus_predictions$genome_drug.resistant_phenotype == "Resistant") / + test_data_plus_predictions$genome_drug.resistant_phenotype == "Resistant" + ) / nrow(test_data_plus_predictions) - if(prior > 0.3 && prior < 0.7) { - warning(paste("Classes are roughly balanced.", - "Calculation of log2(AUPRC/prior) may be inappropriate.")) - } else if(prior >= 0.7) { - warning(paste("Classes are imbalanced toward the resistant phenotype.", - "Calculation of log2(AUPRC/prior) may be inappropriate.")) + if (prior > 0.3 && prior < 0.7) { + warning(paste( + "Classes are roughly balanced.", + "Calculation of log2(AUPRC/prior) may be inappropriate." + )) + } else if (prior >= 0.7) { + warning(paste( + "Classes are imbalanced toward the resistant phenotype.", + "Calculation of log2(AUPRC/prior) may be inappropriate." + )) } - log2_apop <- log2(auprc/prior) |> round(2) + log2_apop <- log2(auprc / prior) |> round(2) return(log2_apop) } @@ -495,16 +525,21 @@ getConfusionMatrix <- function(test_data_plus_predictions) { .calculateBalAcc <- function(test_data_plus_predictions) { .checkArgTestDataPlusPredictions(test_data_plus_predictions) - if(!("genome_drug.resistant_phenotype" %in% + if (!("genome_drug.resistant_phenotype" %in% colnames(test_data_plus_predictions))) { - stop(paste("`test_data_plus_predictions` does not have a column for", - "`genome_drug.resistant_phenotype`.")) + stop(paste( + "`test_data_plus_predictions` does not have a column for", + "`genome_drug.resistant_phenotype`." + )) } bal_acc <- test_data_plus_predictions |> yardstick::bal_accuracy( - truth = genome_drug.resistant_phenotype, estimate = .pred_class) |> - dplyr::select(.estimate) |> as.numeric() |> round(2) + truth = genome_drug.resistant_phenotype, estimate = .pred_class + ) |> + dplyr::select(.estimate) |> + as.numeric() |> + round(2) return(bal_acc) } @@ -519,15 +554,21 @@ getConfusionMatrix <- function(test_data_plus_predictions) { .calculateSensitivity <- function(test_data_plus_predictions) { .checkArgTestDataPlusPredictions(test_data_plus_predictions) - if(!("genome_drug.resistant_phenotype" %in% + if (!("genome_drug.resistant_phenotype" %in% colnames(test_data_plus_predictions))) { - stop(paste("`test_data_plus_predictions` does not have a column for", - "`genome_drug.resistant_phenotype`.")) + stop(paste( + "`test_data_plus_predictions` does not have a column for", + "`genome_drug.resistant_phenotype`." + )) } sens <- test_data_plus_predictions |> - yardstick::sens(truth = genome_drug.resistant_phenotype, - estimate = .pred_class) |> dplyr::select(.estimate) |> as.numeric() |> + yardstick::sens( + truth = genome_drug.resistant_phenotype, + estimate = .pred_class + ) |> + dplyr::select(.estimate) |> + as.numeric() |> round(2) return(sens) @@ -543,15 +584,21 @@ getConfusionMatrix <- function(test_data_plus_predictions) { .calculateSpecificity <- function(test_data_plus_predictions) { .checkArgTestDataPlusPredictions(test_data_plus_predictions) - if(!("genome_drug.resistant_phenotype" %in% + if (!("genome_drug.resistant_phenotype" %in% colnames(test_data_plus_predictions))) { - stop(paste("`test_data_plus_predictions` does not have a column for", - "`genome_drug.resistant_phenotype`.")) + stop(paste( + "`test_data_plus_predictions` does not have a column for", + "`genome_drug.resistant_phenotype`." + )) } spec <- test_data_plus_predictions |> - yardstick::spec(truth = genome_drug.resistant_phenotype, - estimate = .pred_class) |> dplyr::select(.estimate) |> as.numeric() |> + yardstick::spec( + truth = genome_drug.resistant_phenotype, + estimate = .pred_class + ) |> + dplyr::select(.estimate) |> + as.numeric() |> round(2) return(spec) @@ -598,30 +645,36 @@ calculateEvalMets <- function(test_data_plus_predictions) { #' `Importance`, and a column for `Sign` (or, for multi-class, a tibble with #' per-class columns of importance scores for each `Variable`) #' @export -extractTopFeats <- function(fit, prop_vi_top_feats = c(0, 1), - n_top_feats = NA) { +extractTopFeats <- function( + fit, prop_vi_top_feats = c(0, 1), + n_top_feats = NA +) { .checkArgWflow(fit) - if(!is.na(n_top_feats)) {prop_vi_top_feats <- NA} + if (!is.na(n_top_feats)) { + prop_vi_top_feats <- NA + } # Arg checking for every permutation of `prop_vi_top_feats` and `n_top_feats` - if(is.na(n_top_feats) & any(!is.na(prop_vi_top_feats))) { + if (is.na(n_top_feats) & any(!is.na(prop_vi_top_feats))) { .checkArgPropVITopFeats(prop_vi_top_feats) - } else if(any(is.na(prop_vi_top_feats)) & !is.na(n_top_feats)) { + } else if (any(is.na(prop_vi_top_feats)) & !is.na(n_top_feats)) { .checkArgNTopFeats(n_top_feats) - } else if(any(!is.na(prop_vi_top_feats)) & !is.na(n_top_feats)) { + } else if (any(!is.na(prop_vi_top_feats)) & !is.na(n_top_feats)) { stop("Set either `n_top_feats` or `prop_vi_top_feats` to `NA` but not both.") - } else if(any(is.na(prop_vi_top_feats)) & is.na(n_top_feats)) { + } else if (any(is.na(prop_vi_top_feats)) & is.na(n_top_feats)) { stop("Please specify either `n_top_feats` or `prop_vi_top_feats`.") } - feats_arranged <- fit |> workflowsets::extract_fit_parsnip() |> vip::vi() |> + feats_arranged <- fit |> + workflowsets::extract_fit_parsnip() |> + vip::vi() |> dplyr::arrange(dplyr::desc(Importance)) - if(!is.na(n_top_feats)) { + if (!is.na(n_top_feats)) { top_feats_and_VIs <- feats_arranged |> dplyr::slice(1:n_top_feats) - } else if(any(!is.na(prop_vi_top_feats))) { + } else if (any(!is.na(prop_vi_top_feats))) { cum_vi_lower <- prop_vi_top_feats[1] * sum(feats_arranged$Importance) cum_vi_upper <- prop_vi_top_feats[2] * sum(feats_arranged$Importance) @@ -638,9 +691,11 @@ extractTopFeats <- function(fit, prop_vi_top_feats = c(0, 1), # Take a different approach if using multi-class (the previous code would give # a less meaningful result). - if(class(fit$fit$actions$model$spec)[1] == "multinom_reg") { - warning(paste("Extracting top features from a multi-class model.", - "The `prop_vi_top_feats` and `n_top_feats` arguments do not apply.")) + if (class(fit$fit$actions$model$spec)[1] == "multinom_reg") { + warning(paste( + "Extracting top features from a multi-class model.", + "The `prop_vi_top_feats` and `n_top_feats` arguments do not apply." + )) fit_penalty <- .getFitHps(fit)["penalty"] |> as.numeric() glmnet_fit <- parsnip::extract_fit_engine(fit) diff --git a/R/generate_matrices_ml.R b/R/generate_matrices_ml.R index bb1dc68..b19beef 100644 --- a/R/generate_matrices_ml.R +++ b/R/generate_matrices_ml.R @@ -156,7 +156,6 @@ skipImbalancedMatrix <- function(genome_ids, split, stratify_by = NULL, verbosity = c("minimal", "debug")) { - verbosity <- match.arg(verbosity) log <- .make_logger(verbosity) @@ -197,8 +196,10 @@ skipImbalancedMatrix <- function(genome_ids, if (!dir.exists(matrix_path)) dir.create(matrix_path, recursive = TRUE) log("info", paste0("Matrix output directory: ", matrix_path)) - log("debug", paste0("Stratification: ", - ifelse(is.null(stratify_column), "None", stratify_column))) + log("debug", paste0( + "Stratification: ", + ifelse(is.null(stratify_column), "None", stratify_column) + )) # Feature and matrix types feature_types <- list( @@ -220,9 +221,11 @@ skipImbalancedMatrix <- function(genome_ids, # Safe DBI-quoting quote_condition <- function(group_cols, group_values, con) { - ids <- vapply(group_cols, - function(col) DBI::dbQuoteIdentifier(con, col), - character(1)) + ids <- vapply( + group_cols, + function(col) DBI::dbQuoteIdentifier(con, col), + character(1) + ) vals <- vapply( group_cols, function(col) { @@ -256,7 +259,6 @@ skipImbalancedMatrix <- function(genome_ids, log("debug", paste0("Found ", nrow(all_groups), " groups for type: ", group_type)) for (i in seq_len(nrow(all_groups))) { - # New connection for this group con <- DBI::dbConnect(duckdb::duckdb(), parquet_duckdb_path) @@ -268,13 +270,14 @@ skipImbalancedMatrix <- function(genome_ids, condition_string <- quote_condition(group_cols, group_values, con) # Strat filter - strat_filter <- if (!is.null(stratify_column)) + strat_filter <- if (!is.null(stratify_column)) { sprintf("AND \"%s\" IS NOT NULL AND \"%s\" != ''", stratify_column, stratify_column) - else "" + } else { + "" + } # Genome selection logic if (group_type %in% c("drug_class", "drug_class_year", "drug_class_country")) { - genome_ids <- DBI::dbGetQuery(con, sprintf(" WITH class_phenotypes AS ( SELECT \"genome_drug.genome_id\" AS genome_id, @@ -290,7 +293,6 @@ skipImbalancedMatrix <- function(genome_ids, FROM class_phenotypes WHERE any_resistant = 1 OR all_susceptible = 1 ", condition_string))[[1]] - } else { genome_ids <- DBI::dbGetQuery(con, sprintf(" SELECT DISTINCT \"genome_drug.genome_id\" @@ -310,19 +312,24 @@ skipImbalancedMatrix <- function(genome_ids, ", condition_string)) phenotype_summary <- paste( - apply(phenotype_counts_all, 1, - function(row) paste0(row["phenotype"], "=", row["count"])), + apply( + phenotype_counts_all, 1, + function(row) paste0(row["phenotype"], "=", row["count"]) + ), collapse = "; " ) # Apply skip logic if (skipImbalancedMatrix(genome_ids, phenotype_counts_all, n_fold, split, - verbosity = verbosity)) { - + verbosity = verbosity + )) { readr::write_lines( - sprintf("%s\tToo few samples for CV/split\t%d\t%s", - group_label, length(genome_ids), phenotype_summary), - log_path, append = TRUE + sprintf( + "%s\tToo few samples for CV/split\t%d\t%s", + group_label, length(genome_ids), phenotype_summary + ), + log_path, + append = TRUE ) DBI::dbDisconnect(con, shutdown = FALSE) @@ -331,9 +338,12 @@ skipImbalancedMatrix <- function(genome_ids, if (length(genome_ids) < 40) { readr::write_lines( - sprintf("%s\tToo few observations\t%d\t%s", - group_label, length(genome_ids), phenotype_summary), - log_path, append = TRUE + sprintf( + "%s\tToo few observations\t%d\t%s", + group_label, length(genome_ids), phenotype_summary + ), + log_path, + append = TRUE ) DBI::dbDisconnect(con, shutdown = FALSE) @@ -351,9 +361,12 @@ skipImbalancedMatrix <- function(genome_ids, if (nrow(phen2) < 2) { readr::write_lines( - sprintf("%s\tOnly one phenotype class\t%d\t%s", - group_label, length(genome_ids), phenotype_summary), - log_path, append = TRUE + sprintf( + "%s\tOnly one phenotype class\t%d\t%s", + group_label, length(genome_ids), phenotype_summary + ), + log_path, + append = TRUE ) DBI::dbDisconnect(con, shutdown = FALSE) @@ -363,13 +376,14 @@ skipImbalancedMatrix <- function(genome_ids, # Create selected_genomes DBI::dbExecute(con, "CREATE OR REPLACE TEMP TABLE selected_genomes (genome_id VARCHAR)") DBI::dbWriteTable(con, "selected_genomes", - data.frame(genome_id = genome_ids), append = TRUE) + data.frame(genome_id = genome_ids), + append = TRUE + ) # Feature and matrix generation steps for (ftype in names(feature_types)) { - fview <- feature_types[[ftype]]$view - fid <- feature_types[[ftype]]$id_col + fid <- feature_types[[ftype]]$id_col # binary view DBI::dbExecute(con, sprintf(" @@ -389,13 +403,14 @@ skipImbalancedMatrix <- function(genome_ids, } for (mtype in names(matrix_types)) { - binary_only <- matrix_types[[mtype]]$binary_only if (ftype == "struct" && !binary_only) next - mview <- sprintf("%s_%s", ftype, - ifelse(grepl("binary", mtype), "binary", "counts")) - value_col <- matrix_types[[mtype]]$value_col + mview <- sprintf( + "%s_%s", ftype, + ifelse(grepl("binary", mtype), "binary", "counts") + ) + value_col <- matrix_types[[mtype]]$value_col filter_clause <- matrix_types[[mtype]]$filter # select features with non-zero variance @@ -409,29 +424,38 @@ skipImbalancedMatrix <- function(genome_ids, keep_features <- DBI::dbGetQuery(con, keep_query)[["feature_id"]] if (length(keep_features) == 0) { - log("info", paste0("All features filtered for ", - ftype, " - ", mtype, " - ", group_label)) + log("info", paste0( + "All features filtered for ", + ftype, " - ", mtype, " - ", group_label + )) next } - DBI::dbExecute(con, - "CREATE OR REPLACE TEMP TABLE keep_features (feature_id VARCHAR)") + DBI::dbExecute( + con, + "CREATE OR REPLACE TEMP TABLE keep_features (feature_id VARCHAR)" + ) DBI::dbWriteTable(con, - "keep_features", - data.frame(feature_id = keep_features), - append = TRUE) + "keep_features", + data.frame(feature_id = keep_features), + append = TRUE + ) mtype_label <- matrix_types[[mtype]]$label - long_out_path <- file.path(matrix_path, - sprintf("%s_%s_%s_%s_%s_sparse.parquet", - bug, group_type, group_label, ftype, mtype_label)) + long_out_path <- file.path( + matrix_path, + sprintf( + "%s_%s_%s_%s_%s_sparse.parquet", + bug, group_type, group_label, ftype, mtype_label + ) + ) long_out_path_sql <- gsub("\\\\", "/", long_out_path) # phenotype case phenotype_case <- if (group_type %in% - c("drug_class", "drug_class_year", "drug_class_country")) { + c("drug_class", "drug_class_year", "drug_class_country")) { " CASE WHEN MAX(CASE WHEN f.\"genome_drug.resistant_phenotype\"='Resistant' @@ -451,13 +475,20 @@ skipImbalancedMatrix <- function(genome_ids, " } - strat_col_select <- if (!is.null(stratify_by)) - sprintf(", f.\"%s\"", stratify_column) else "" + strat_col_select <- if (!is.null(stratify_by)) { + sprintf(", f.\"%s\"", stratify_column) + } else { + "" + } - strat_col_group <- if (!is.null(stratify_by)) - sprintf(", f.\"%s\"", stratify_column) else "" + strat_col_group <- if (!is.null(stratify_by)) { + sprintf(", f.\"%s\"", stratify_column) + } else { + "" + } - copy_sql <- sprintf(" + copy_sql <- sprintf( + " COPY ( SELECT f.\"genome_drug.genome_id\" AS genome_id, @@ -478,18 +509,21 @@ skipImbalancedMatrix <- function(genome_ids, TO '%s' (FORMAT 'parquet', COMPRESSION 'zstd') ", - fid, value_col, phenotype_case, strat_col_select, - mview, fid, condition_string, - strat_filter, fid, strat_col_group, fid, - long_out_path_sql) + fid, value_col, phenotype_case, strat_col_select, + mview, fid, condition_string, + strat_filter, fid, strat_col_group, fid, + long_out_path_sql + ) ok <- try(DBI::dbExecute(con, copy_sql), silent = TRUE) # On copy failure, log + continue without stopping entire pipeline if (inherits(ok, "try-error")) { readr::write_lines( - sprintf("%s\tCOPY_failed\t%d\t%s", - group_label, length(genome_ids), phenotype_summary), + sprintf( + "%s\tCOPY_failed\t%d\t%s", + group_label, length(genome_ids), phenotype_summary + ), log_path, append = TRUE ) @@ -530,7 +564,7 @@ skipImbalancedMatrix <- function(genome_ids, # Normalize paths to forward slashes for consistency matrix_path <- gsub("\\\\", "/", file.path(path, paste0("matrix_", stratify_by))) - LOO_path <- gsub("\\\\", "/", file.path(path, paste0("LOO_matrix_", stratify_by))) + LOO_path <- gsub("\\\\", "/", file.path(path, paste0("LOO_matrix_", stratify_by))) if (!dir.exists(matrix_path)) { log("info", paste0("The matrix directory ", matrix_path, " does not exist.")) @@ -626,9 +660,11 @@ skipImbalancedMatrix <- function(genome_ids, out_file <- gsub("\\\\", "/", file.path( LOO_path, - paste0(sub_prefix, "_", stratify_by, "_", - drug_class, "_leaveout_", leave_one_out, "_", - sub_feature, "_sparse.parquet") + paste0( + sub_prefix, "_", stratify_by, "_", + drug_class, "_leaveout_", leave_one_out, "_", + sub_feature, "_sparse.parquet" + ) )) arrow::write_parquet(combined, out_file) created <<- c(created, out_file) @@ -702,7 +738,7 @@ skipImbalancedMatrix <- function(genome_ids, # Build one matrix per feature type and matrix type for (ftype in names(feature_types)) { fview <- feature_types[[ftype]]$view - fid <- feature_types[[ftype]]$id_col + fid <- feature_types[[ftype]]$id_col for (mtype in names(matrix_types)) { binary_only <- matrix_types[[mtype]]$binary_only @@ -722,8 +758,9 @@ skipImbalancedMatrix <- function(genome_ids, # Selected genomes DBI::dbExecute(con, "CREATE OR REPLACE TEMP TABLE selected_genomes (genome_id VARCHAR)") DBI::dbWriteTable(con, "selected_genomes", - data.frame(genome_id = genomes_to_keep), - append = TRUE) + data.frame(genome_id = genomes_to_keep), + append = TRUE + ) # Binary view DBI::dbExecute(con, sprintf(" @@ -763,13 +800,15 @@ skipImbalancedMatrix <- function(genome_ids, DBI::dbExecute(con, "CREATE OR REPLACE TEMP TABLE keep_features (feature_id VARCHAR)") DBI::dbWriteTable(con, "keep_features", - data.frame(feature_id = keep_features), - append = TRUE) + data.frame(feature_id = keep_features), + append = TRUE + ) + - - copy_sql <- sprintf(" + copy_sql <- sprintf( + " COPY ( - SELECT + SELECT f.\"genome_drug.genome_id\" AS genome_id, %s AS feature_id, MAX(CAST(%s AS DOUBLE)) AS value, @@ -779,26 +818,26 @@ skipImbalancedMatrix <- function(genome_ids, JOIN keep_features kf ON %s = kf.feature_id JOIN metadata f ON genome_id = f.\"genome_drug.genome_id\" WHERE resistant_classes <> 'Intermediate' - GROUP BY - f.\"genome_drug.genome_id\", - %s, + GROUP BY + f.\"genome_drug.genome_id\", + %s, resistant_classes - ORDER BY - f.\"genome_drug.genome_id\", + ORDER BY + f.\"genome_drug.genome_id\", %s ) TO '%s' (FORMAT 'parquet', COMPRESSION 'zstd') - ", - fid, # %s -> feature_id expression column name - value_col, # %s -> value column to CAST - mview, # %s -> source view (binary or counts) - fid, # %s -> join to keep_features - fid, # %s -> group by feature id - fid, # %s -> order by feature id - out_file_sql # %s -> destination parquet file + ", + fid, # %s -> feature_id expression column name + value_col, # %s -> value column to CAST + mview, # %s -> source view (binary or counts) + fid, # %s -> join to keep_features + fid, # %s -> group by feature id + fid, # %s -> order by feature id + out_file_sql # %s -> destination parquet file ) - + ok <- try(DBI::dbExecute(con, copy_sql), silent = TRUE) if (inherits(ok, "try-error")) { log("info", paste0("COPY failed for MDR matrix: ", out_file)) diff --git a/R/globals.R b/R/globals.R index a6595d2..131d016 100644 --- a/R/globals.R +++ b/R/globals.R @@ -8,7 +8,6 @@ "_PACKAGE" utils::globalVariables(c( - # Prediction columns from tidymodels ".estimate", ".pred_Resistant", @@ -52,7 +51,6 @@ utils::globalVariables(c( "pair_id", "parts", "phenotype", - "precision", "prefix", "prefix_key", diff --git a/R/plot_ml.R b/R/plot_ml.R index b41efb5..42d51ec 100644 --- a/R/plot_ml.R +++ b/R/plot_ml.R @@ -22,33 +22,33 @@ NULL #' Plot a Precision-Recall Curve -#' +#' #' Generates a precision-recall curve (PRC) for AMR phenotype prediction results. #' @param test_data_plus_predictions A tibble containing test data with added #' prediction columns, typically the output of `runMLmodels()`. -#' +#' #' @return A `ggplot2` object showing the precision-recall curve. -#' -#' @details +#' +#' @details #' The function uses `yardstick::pr_curve()` to compute the PR curve and then #' visualizes it using `ggplot2`. -#' +#' #' @examples #' \dontrun{ -#' test_data_plus_predictions <- readr::read_tsv(results/ML_pred/Sfl_drug_AMP_domains_binary_prediction.tsv) +#' test_data_plus_predictions <- readr::read_tsv(results / ML_pred / Sfl_drug_AMP_domains_binary_prediction.tsv) #' plotPRC(test_data_plus_predictions) -#' } -#' +#' } +#' #' @export plotPRC <- function(test_data_plus_predictions) { .checkArgTestDataPlusPredictions(test_data_plus_predictions) -test_data_plus_predictions <- test_data_plus_predictions |> -dplyr::mutate( -genome_drug.resistant_phenotype = factor( -genome_drug.resistant_phenotype, -levels = c("Resistant", "Susceptible") -) -) + test_data_plus_predictions <- test_data_plus_predictions |> + dplyr::mutate( + genome_drug.resistant_phenotype = factor( + genome_drug.resistant_phenotype, + levels = c("Resistant", "Susceptible") + ) + ) prc <- yardstick::pr_curve( test_data_plus_predictions, @@ -63,24 +63,24 @@ levels = c("Resistant", "Susceptible") } #' Plot a Receiver Operating Characteristic (ROC) Curve -#' +#' #' Generates a ROC curve for AMR phenotype prediction results. -#' +#' #' @param test_data_plus_predictions A tibble with test data and prediction #' columns (output of `runMLmodels()`). -#' +#' #' @return A ROC curve plotted using `ggplot2::autoplot()`. -#' +#' #' @export plotROC <- function(test_data_plus_predictions) { .checkArgTestDataPlusPredictions(test_data_plus_predictions) test_data_plus_predictions <- test_data_plus_predictions |> -dplyr::mutate( -genome_drug.resistant_phenotype = factor( -genome_drug.resistant_phenotype, -levels = c("Resistant", "Susceptible") -) -) + dplyr::mutate( + genome_drug.resistant_phenotype = factor( + genome_drug.resistant_phenotype, + levels = c("Resistant", "Susceptible") + ) + ) roc <- yardstick::roc_curve( test_data_plus_predictions, @@ -88,19 +88,19 @@ levels = c("Resistant", "Susceptible") ) |> ggplot2::autoplot(type = "se") + ggplot2::theme(panel.grid = ggplot2::element_blank()) - + return(roc) } #' Plot a Confusion Matrix Heatmap -#' +#' #' Produces a heatmap visualization of the confusion matrix for AMR predictions. -#' +#' #' @param test_data_plus_predictions A tibble containing true and predicted #' phenotype labels. -#' +#' #' @return A heatmap (`ggplot2` object) showing the confusion matrix. -#' +#' #' @export plotCM <- function(test_data_plus_predictions) { .checkArgTestDataPlusPredictions(test_data_plus_predictions) @@ -116,62 +116,66 @@ plotCM <- function(test_data_plus_predictions) { ) ) test_data_plus_predictions |> -yardstick::conf_mat(truth = genome_drug.resistant_phenotype, - estimate = .pred_class) |> -ggplot2::autoplot(type = "heatmap") + yardstick::conf_mat( + truth = genome_drug.resistant_phenotype, + estimate = .pred_class + ) |> + ggplot2::autoplot(type = "heatmap") } #' Plot Density of Predicted Class Probabilities -#' +#' #' Visualizes how predicted class probabilities differ between resistant and #' susceptible genome-drug combinations. -#' +#' #' @param test_data_plus_predictions Tibble with prediction probabilities and #' true labels. -#' +#' #' @return A ggplot2 density plot. -#' +#' #' @export plotDensity <- function(test_data_plus_predictions) { test_data_plus_predictions |> -ggplot2::ggplot(ggplot2::aes(x = .pred_Resistant, -fill = genome_drug.resistant_phenotype)) + -ggplot2::geom_density(alpha = 0.5) + ggplot2::ggplot(ggplot2::aes( + x = .pred_Resistant, + fill = genome_drug.resistant_phenotype + )) + + ggplot2::geom_density(alpha = 0.5) } #' Plot Top Feature Importances -#' +#' #' Creates a bar plot showing the most important features affecting #' AMR phenotype predictions. -#' +#' #' @param topfeat A tibble containing feature importance scores #' (output of `runMLmodels()`). #' @param n_top_feats Number of top features to display (default: 10). -#' +#' #' @return A bar plot of variable importance (`ggplot2` object). -#' +#' #' @examples #' \dontrun{ -#' topfeat <- readr::read_tsv(results/ML_top_features/Sfl_drug_AMP_domains_binary_top_features.tsv) +#' topfeat <- readr::read_tsv(results / ML_top_features / Sfl_drug_AMP_domains_binary_top_features.tsv) #' plotTopFeatsVI(topfeat) -#' } -#' +#' } +#' #' @export plotTopFeatsVI <- function(topfeat, n_top_feats = 10) { .checkArgNTopFeats(n_top_feats) - vip <- topfeat |> + vip <- topfeat |> dplyr::slice_max(order_by = Importance, n = n_top_feats) |> dplyr::mutate( - Variable = factor(Variable, levels = rev(Variable)), # preserve order as shown in table + Variable = factor(Variable, levels = rev(Variable)), # preserve order as shown in table Sign = factor(Sign, levels = c("POS", "NEG")) - ) |> + ) |> ggplot2::ggplot(ggplot2::aes(x = Importance, y = Variable, fill = Sign)) + ggplot2::geom_col() + ggplot2::scale_fill_manual( values = c( - "POS" = "#c6d8d3", - "NEG" = "#f6c9a1" + "POS" = "#c6d8d3", + "NEG" = "#f6c9a1" ) ) + ggplot2::labs( @@ -183,20 +187,20 @@ plotTopFeatsVI <- function(topfeat, n_top_feats = 10) { panel.grid.minor = ggplot2::element_blank(), axis.text.y = ggplot2::element_text(size = 10) ) - - return(vip) + + return(vip) } #' Compare Baseline Performance With and Without Shuffled Labels -#' +#' #' Produces a bar plot comparing balanced accuracy for each antibiotic using #' true AMR labels vs. randomly shuffled labels. -#' +#' #' @param non_shuffled_label_results Output of `runMLPipeline(shuffle_labels = FALSE)` #' @param shuffled_label_results Output of `runMLPipeline(shuffle_labels = TRUE)` -#' +#' #' @return A base R barplot comparing balanced accuracy across models. -#' +#' #' @export plotBaselineComparison <- function( non_shuffled_label_results, @@ -268,7 +272,6 @@ plotFishers <- function( alpha = 0.05, label_top_n = 5 ) { - required_cols <- c("gene", "adj_p_value", "sig_after_bh") missing_cols <- setdiff(required_cols, colnames(fisher_df)) diff --git a/R/prep_ml.R b/R/prep_ml.R index d47c160..4a5954e 100644 --- a/R/prep_ml.R +++ b/R/prep_ml.R @@ -111,8 +111,10 @@ loadMLInputTibble <- function(parquet_path) { if (exists(".ml_logger")) { log <- .ml_logger("minimal") - log("debug", paste0("ML tibble constructed: ", nrow(ml_input_tibble), - " genomes x ", getNumFeat(ml_input_tibble), " features")) + log("debug", paste0( + "ML tibble constructed: ", nrow(ml_input_tibble), + " genomes x ", getNumFeat(ml_input_tibble), " features" + )) } if (anyDuplicated(dplyr::pull(ml_input_tibble, genome_id)) != 0) { diff --git a/R/run_ML.R b/R/run_ML.R index eba37f8..2ed07e7 100644 --- a/R/run_ML.R +++ b/R/run_ML.R @@ -4,9 +4,11 @@ #' the ML matrices with these new split/CV values instead. #' @noRd .resolveSplitParams <- function(parquet_path, - defaults = list(split = c(0.8, 0), - seed = 5280, - n_fold = 5)) { + defaults = list( + split = c(0.8, 0), + seed = 5280, + n_fold = 5 + )) { # matrix_dir is the directory that contains the parquet files matrix_dir <- normalizePath(dirname(parquet_path)) params_json <- .readMLParameters(matrix_dir) @@ -16,8 +18,8 @@ } list( - split = if (!is.null(params_json$split)) params_json$split else defaults$split, - seed = if (!is.null(params_json$seed)) params_json$seed else defaults$seed, + split = if (!is.null(params_json$split)) params_json$split else defaults$split, + seed = if (!is.null(params_json$seed)) params_json$seed else defaults$seed, n_fold = if (!is.null(params_json$n_fold)) params_json$n_fold else defaults$n_fold ) } @@ -53,8 +55,9 @@ #' #' # LOO analysis stratified by year #' paths_loo <- createMLResultDir("/path/to/results", -#' stratify_by = "year", -#' LOO = TRUE) +#' stratify_by = "year", +#' LOO = TRUE +#' ) #' #' # MDR analysis #' paths_mdr <- createMLResultDir("/path/to/results", MDR = TRUE) @@ -90,16 +93,17 @@ createMLResultDir <- function(path, ) } else { # Determine prefixes (only in non-MDR mode) - full_prefix <- paste0(ifelse(isTRUE(LOO), "LOO_", ""), - ifelse(isTRUE(cross_test), "cross_test_", "")) + full_prefix <- paste0( + ifelse(isTRUE(LOO), "LOO_", ""), + ifelse(isTRUE(cross_test), "cross_test_", "") + ) half_prefix <- ifelse(isTRUE(LOO), "LOO_", "") # Determine suffix suffix <- if (is.null(stratify_by) || identical(stratify_by, "")) { "" } else { - switch( - stratify_by, + switch(stratify_by, "country" = "_country", "year" = "_year", stop("`stratify_by` must be NULL, 'country', or 'year'.") @@ -127,20 +131,20 @@ createMLResultDir <- function(path, return(paths) } - # createAllMLResultDir <- function(path) { - # createMLResultDir(path, stratify_by = NULL, LOO = FALSE, cross_test = FALSE, MDR = FALSE) - # createMLResultDir(path, stratify_by = NULL, LOO = FALSE, cross_test = TRUE, MDR = FALSE) - # createMLResultDir(path, stratify_by = NULL, LOO = FALSE, cross_test = FALSE, MDR = TRUE) - # createMLResultDir(path, stratify_by = "year", LOO = FALSE, cross_test = FALSE, MDR = FALSE) - # createMLResultDir(path, stratify_by = "year", LOO = FALSE, cross_test = TRUE, MDR = FALSE) - # createMLResultDir(path, stratify_by = "year", LOO = TRUE, cross_test = FALSE, MDR = FALSE) - # createMLResultDir(path, stratify_by = "year", LOO = TRUE, cross_test = TRUE, MDR = FALSE) - # createMLResultDir(path, stratify_by = "country", LOO = FALSE, cross_test = FALSE, MDR = FALSE) - # createMLResultDir(path, stratify_by = "country", LOO = FALSE, cross_test = TRUE, MDR = FALSE) - # createMLResultDir(path, stratify_by = "country", LOO = TRUE, cross_test = FALSE, MDR = FALSE) - # createMLResultDir(path, stratify_by = "country", LOO = TRUE, cross_test = TRUE, MDR = FALSE) - # } - # +# createAllMLResultDir <- function(path) { +# createMLResultDir(path, stratify_by = NULL, LOO = FALSE, cross_test = FALSE, MDR = FALSE) +# createMLResultDir(path, stratify_by = NULL, LOO = FALSE, cross_test = TRUE, MDR = FALSE) +# createMLResultDir(path, stratify_by = NULL, LOO = FALSE, cross_test = FALSE, MDR = TRUE) +# createMLResultDir(path, stratify_by = "year", LOO = FALSE, cross_test = FALSE, MDR = FALSE) +# createMLResultDir(path, stratify_by = "year", LOO = FALSE, cross_test = TRUE, MDR = FALSE) +# createMLResultDir(path, stratify_by = "year", LOO = TRUE, cross_test = FALSE, MDR = FALSE) +# createMLResultDir(path, stratify_by = "year", LOO = TRUE, cross_test = TRUE, MDR = FALSE) +# createMLResultDir(path, stratify_by = "country", LOO = FALSE, cross_test = FALSE, MDR = FALSE) +# createMLResultDir(path, stratify_by = "country", LOO = FALSE, cross_test = TRUE, MDR = FALSE) +# createMLResultDir(path, stratify_by = "country", LOO = TRUE, cross_test = FALSE, MDR = FALSE) +# createMLResultDir(path, stratify_by = "country", LOO = TRUE, cross_test = TRUE, MDR = FALSE) +# } +# #' Create machine learning input list #' @@ -174,8 +178,9 @@ createMLResultDir <- function(path, #' #' # Cross-test with year stratification #' inputs_ct <- createMLinputList("/path/to/results", -#' stratify_by = "year", -#' cross_test = TRUE) +#' stratify_by = "year", +#' cross_test = TRUE +#' ) #' #' # MDR analysis #' inputs_mdr <- createMLinputList("/path/to/results", MDR = TRUE) @@ -187,10 +192,10 @@ createMLinputList <- function(path, LOO = FALSE, MDR = FALSE, cross_test = FALSE) { - # Validate inputs - if (!is.character(path) || length(path) != 1 || is.na(path)) + if (!is.character(path) || length(path) != 1 || is.na(path)) { stop("`path` must be a valid file path string.") + } path <- normalizePath(path) @@ -225,21 +230,17 @@ createMLinputList <- function(path, # Multi-drug resistance models # ============================ if (MDR) { - parsed <- tibble::tibble(ref_file = files_vec) |> dplyr::mutate( parts = stringr::str_split(basename(ref_file), "_"), - species = purrr::map_chr(parts, ~ .x[1]), - mdr_tag = purrr::map_chr(parts, ~ .x[2]), # always "MDR" + mdr_tag = purrr::map_chr(parts, ~ .x[2]), # always "MDR" phenotype = purrr::map_chr(parts, ~ paste(.x[3:4], collapse = "_")), # Feature is 5th + 6th tokens feature_type = purrr::map_chr(parts, ~ .x[5]), feature_subtype = purrr::map_chr(parts, ~ stringr::str_remove(.x[6], "_sparse.parquet")), - feature = purrr::map2_chr(feature_type, feature_subtype, paste, sep = "_"), - output_prefix = paste0("MDR_", phenotype, "_", feature) ) @@ -247,38 +248,43 @@ createMLinputList <- function(path, dplyr::mutate( test_file = NA_character_, matrix_path = paths$matrix_path, - out_perf = paths$ML_performance, - out_top = paths$ML_top_features, - out_models= paths$ML_models, - out_pred = paths$ML_prediction + out_perf = paths$ML_performance, + out_top = paths$ML_top_features, + out_models = paths$ML_models, + out_pred = paths$ML_prediction ) return(out) - # ============================ - # For all other modeling types - # ============================ + # ============================ + # For all other modeling types + # ============================ } else { - parsed <- tibble::tibble(ref_file = files_vec) |> dplyr::mutate( - parts = stringr::str_split(basename(ref_file), "_"), + parts = stringr::str_split(basename(ref_file), "_"), i_sparse = purrr::map_int(parts, ~ .get_idx(.x, "sparse.parquet")), - i_strat = purrr::map_int(parts, ~ { - if (is.null(stratify_by)) return(NA_integer_) + i_strat = purrr::map_int(parts, ~ { + if (is.null(stratify_by)) { + return(NA_integer_) + } .get_idx(.x, stratify_by) }), # Feature = last two tokens before sparse.parquet feature = purrr::map2_chr(parts, i_sparse, ~ { - i <- .y; x <- .x - if (is.na(i) || i < 3) return(NA_character_) + i <- .y + x <- .x + if (is.na(i) || i < 3) { + return(NA_character_) + } paste(x[(i - 2):(i - 1)], collapse = "_") }), # Drug or drug class extraction drug_or_class = purrr::map2_chr(parts, i_strat, ~ { - i <- .y; x <- .x + i <- .y + x <- .x # Stratified models if (!is.na(i)) { @@ -304,32 +310,40 @@ createMLinputList <- function(path, # Stratification value (if present) strat_value = purrr::map2_chr(parts, i_strat, ~ { - i <- .y; x <- .x - if (is.na(i)) return("") + i <- .y + x <- .x + if (is.na(i)) { + return("") + } # default position is two tokens after the strat label j <- i + 2 # if there's an intervening 'leaveout', skip over it if (j <= length(x) && identical(x[j], "leaveout")) j <- j + 1 - if (j <= length(x)) return(x[j]) - "" # no stratification + if (j <= length(x)) { + return(x[j]) + } + "" # no stratification }), # Prefix key for grouping prefix_key = purrr::map2_chr(parts, i_strat, ~ { - i <- .y; x <- .x + i <- .y + x <- .x # Case A: stratified -> prefix before the stratify label if (!is.na(i)) { - if (i - 1 >= 1) return(paste(x[1:(i - 1)], collapse = "_")) + if (i - 1 >= 1) { + return(paste(x[1:(i - 1)], collapse = "_")) + } return("") } # Case B: unstratified -> prefix is first two tokens - if (x[2] == "drug" && x[3] != "class"){ + if (x[2] == "drug" && x[3] != "class") { # Case A: Cje_drug_X return(paste(x[1:2], collapse = "_")) } - if (x[2] == "drug" && x[3] == "class"){ + if (x[2] == "drug" && x[3] == "class") { # Case A: Cje_drug_X return(paste(x[1:3], collapse = "_")) } @@ -345,18 +359,17 @@ createMLinputList <- function(path, test_file = NA_character_, output_prefix = gsub("_sparse\\.parquet$", "", basename(ref_file)), matrix_path = paths$matrix_path, - out_perf = paths$ML_performance, - out_top = paths$ML_top_features, - out_models= paths$ML_models, - out_pred = paths$ML_prediction + out_perf = paths$ML_performance, + out_top = paths$ML_top_features, + out_models = paths$ML_models, + out_pred = paths$ML_prediction ) return(out) - # ============================ - # Cross-test modeling, no LOO - # ============================ + # ============================ + # Cross-test modeling, no LOO + # ============================ } else if (cross_test && !LOO) { - if (is.null(stratify_by)) { # Case A: stratify_by = NULL, pair across abx within same feature + prefix pairs <- parsed |> @@ -366,8 +379,10 @@ createMLinputList <- function(path, dplyr::select(test_file = ref_file, feature, prefix_key, strat_value, test_drug = drug_or_class), by = c("feature", "prefix_key", "strat_value") ) |> - dplyr::filter(ref_file != test_file, - ref_drug != test_drug) |> + dplyr::filter( + ref_file != test_file, + ref_drug != test_drug + ) |> dplyr::distinct() |> dplyr::mutate( output_prefix = paste0( @@ -380,10 +395,10 @@ createMLinputList <- function(path, out <- pairs |> dplyr::mutate( matrix_path = paths$matrix_path, - out_perf = paths$ML_performance, - out_top = paths$ML_top_features, - out_models= paths$ML_models, - out_pred = paths$ML_prediction + out_perf = paths$ML_performance, + out_top = paths$ML_top_features, + out_models = paths$ML_models, + out_pred = paths$ML_prediction ) return(out) @@ -392,30 +407,29 @@ createMLinputList <- function(path, # Case B: stratify_by != NULL, pair same drug/class, prefix, feature, # but across different stratification groups pairs <- parsed |> - dplyr::select(ref_file, feature, prefix_key, strat_value, - drug_or_class) |> - + dplyr::select( + ref_file, feature, prefix_key, strat_value, + drug_or_class + ) |> # self-join ONLY on prefix_key, drug/class, feature dplyr::inner_join( parsed |> - dplyr::select(test_file = ref_file, - feature, prefix_key, strat_value_test = strat_value, - drug_or_class), + dplyr::select( + test_file = ref_file, + feature, prefix_key, strat_value_test = strat_value, + drug_or_class + ), by = c("prefix_key", "feature", "drug_or_class") ) |> - # do NOT test file against itself dplyr::filter(ref_file != test_file) |> - # enforce different stratification group dplyr::filter(strat_value != strat_value_test) |> - # remove symmetric duplicates (A,B == B,A) dplyr::rowwise() |> dplyr::mutate(pair_id = paste(sort(c(ref_file, test_file)), collapse = "||")) |> dplyr::ungroup() |> dplyr::distinct(pair_id, .keep_all = TRUE) |> - dplyr::mutate( output_prefix = paste0( prefix_key, "_", @@ -429,19 +443,18 @@ createMLinputList <- function(path, out <- pairs |> dplyr::mutate( matrix_path = paths$matrix_path, - out_perf = paths$ML_performance, - out_top = paths$ML_top_features, - out_models= paths$ML_models, - out_pred = paths$ML_prediction + out_perf = paths$ML_performance, + out_top = paths$ML_top_features, + out_models = paths$ML_models, + out_pred = paths$ML_prediction ) return(out) - # ============================ - # Cross-test + LOO modeling - # ============================ + # ============================ + # Cross-test + LOO modeling + # ============================ } else if (cross_test && LOO) { - # LOO requires special directory structure resolution test_path <- file.path(path, stringr::str_remove(basename(paths$matrix_path), "^LOO_")) test_path <- normalizePath(test_path) @@ -461,10 +474,10 @@ createMLinputList <- function(path, out <- loo_pairs |> dplyr::mutate( matrix_path = paths$matrix_path, - out_perf = paths$ML_performance, - out_top = paths$ML_top_features, - out_models= paths$ML_models, - out_pred = paths$ML_prediction + out_perf = paths$ML_performance, + out_top = paths$ML_top_features, + out_models = paths$ML_models, + out_pred = paths$ML_prediction ) return(out) @@ -472,9 +485,11 @@ createMLinputList <- function(path, } # If we ever get here, something wasn't covered - stop("Unhandled combination of arguments: ", - "MDR=", MDR, ", cross_test=", cross_test, ", LOO=", LOO, - ", stratify_by=", if (is.null(stratify_by)) "NULL" else stratify_by) + stop( + "Unhandled combination of arguments: ", + "MDR=", MDR, ", cross_test=", cross_test, ", LOO=", LOO, + ", stratify_by=", if (is.null(stratify_by)) "NULL" else stratify_by + ) } @@ -544,13 +559,15 @@ createMLinputList <- function(path, #' #' # Run with more threads and minimal output #' runMDRmodels("/path/to/results", -#' threads = 32, -#' verbose = FALSE) +#' threads = 32, +#' verbose = FALSE +#' ) #' #' # Run without saving model fits (save disk space) #' runMDRmodels("/path/to/results", -#' threads = 16, -#' return_fit = FALSE) +#' threads = 16, +#' return_fit = FALSE +#' ) #' } #' #' @seealso @@ -571,12 +588,12 @@ runMDRmodels <- function(path, use_saved_split = TRUE, shuffle_labels = FALSE, use_pca = FALSE) { - files <- createMLinputList(path, - stratify_by = NULL, - LOO = FALSE, - cross_test = FALSE, - MDR = TRUE) + stratify_by = NULL, + LOO = FALSE, + cross_test = FALSE, + MDR = TRUE + ) if (nrow(files) == 0) { message("No MDR files found to process. Exiting.") @@ -594,18 +611,19 @@ runMDRmodels <- function(path, # Auto tags for shuffled and PCA shuffle_tag <- if (isTRUE(shuffle_labels)) "shuffled_" else "" - pca_tag <- if (isTRUE(use_pca)) paste0("_pca", as.character(pca_threshold)) else "" + pca_tag <- if (isTRUE(use_pca)) paste0("_pca", as.character(pca_threshold)) else "" results_list <- future.apply::future_lapply( seq_len(nrow(files)), FUN = function(i) { - - ref_parquet <- files$ref_file[i] + ref_parquet <- files$ref_file[i] output_prefix <- files$output_prefix[i] if (interactive()) { - message(sprintf("[runMDRmodels] %d/%d: %s", - i, nrow(files), basename(ref_parquet))) + message(sprintf( + "[runMDRmodels] %d/%d: %s", + i, nrow(files), basename(ref_parquet) + )) } ml_input <- loadMLInputTibble(ref_parquet) @@ -619,32 +637,37 @@ runMDRmodels <- function(path, list(split = split, seed = 5280, n_fold = n_fold) } - res <- try({ - runMLPipeline( - ml_input_tibble = ml_input, - test_data = NA, - model = "LR", - split = sp$split, - n_fold = sp$n_fold, - prop_vi_top_feats = prop_vi_top_feats, - n_top_feats = NA, - use_pca = use_pca, - pca_threshold = pca_threshold, - shuffle_labels = shuffle_labels, - penalty_vec = 10^seq(-4, -1, length.out = 10), - mix_vec = 0:5 / 5, - select_best_metric = "mcc", - seed = sp$seed, - verbose = verbose, - return_tune_res = return_tune_res, - return_fit = return_fit, - return_pred = return_pred - ) - }, silent = TRUE) + res <- try( + { + runMLPipeline( + ml_input_tibble = ml_input, + test_data = NA, + model = "LR", + split = sp$split, + n_fold = sp$n_fold, + prop_vi_top_feats = prop_vi_top_feats, + n_top_feats = NA, + use_pca = use_pca, + pca_threshold = pca_threshold, + shuffle_labels = shuffle_labels, + penalty_vec = 10^seq(-4, -1, length.out = 10), + mix_vec = 0:5 / 5, + select_best_metric = "mcc", + seed = sp$seed, + verbose = verbose, + return_tune_res = return_tune_res, + return_fit = return_fit, + return_pred = return_pred + ) + }, + silent = TRUE + ) if (inherits(res, "try-error")) { - warning("Model failed for: ", output_prefix, - "\n Error: ", attr(res, "condition")$message) + warning( + "Model failed for: ", output_prefix, + "\n Error: ", attr(res, "condition")$message + ) return(NULL) } @@ -652,19 +675,25 @@ runMDRmodels <- function(path, base <- paste0(shuffle_tag, output_prefix, pca_tag) if (!is.null(res$performance_tibble)) { - readr::write_tsv(res$performance_tibble, - file.path(files$out_perf[i], paste0(base, "_performance.tsv"))) + readr::write_tsv( + res$performance_tibble, + file.path(files$out_perf[i], paste0(base, "_performance.tsv")) + ) } if (!is.null(res$top_feat_tibble)) { - readr::write_tsv(res$top_feat_tibble, - file.path(files$out_top[i], paste0(base, "_top_features.tsv"))) + readr::write_tsv( + res$top_feat_tibble, + file.path(files$out_top[i], paste0(base, "_top_features.tsv")) + ) } if (!is.null(res$fit)) { saveRDS(res$fit, file.path(files$out_models[i], paste0(base, "_model_fit.rds"))) } if (!is.null(res$pred)) { - readr::write_tsv(res$pred, - file.path(files$out_pred[i], paste0(base, "_prediction.tsv"))) + readr::write_tsv( + res$pred, + file.path(files$out_pred[i], paste0(base, "_prediction.tsv")) + ) } NULL @@ -783,21 +812,24 @@ runMDRmodels <- function(path, #' #' # Cross-test with year stratification #' runMLmodels("/path/to/results", -#' stratify_by = "year", -#' cross_test = TRUE, -#' threads = 32) +#' stratify_by = "year", +#' cross_test = TRUE, +#' threads = 32 +#' ) #' #' # LOO analysis stratified by country with cross-testing #' runMLmodels("/path/to/results", -#' stratify_by = "country", -#' LOO = TRUE, -#' cross_test = TRUE, -#' verbose = TRUE) +#' stratify_by = "country", +#' LOO = TRUE, +#' cross_test = TRUE, +#' verbose = TRUE +#' ) #' #' # Run without saving model fits (save disk space) #' runMLmodels("/path/to/results", -#' stratify_by = "year", -#' return_fit = FALSE) +#' stratify_by = "year", +#' return_fit = FALSE +#' ) #' } #' #' @seealso @@ -823,19 +855,21 @@ runMLmodels <- function(path, use_saved_split = TRUE, shuffle_labels = FALSE, use_pca = FALSE) { - if (!is.null(stratify_by)) { - if (!is.character(stratify_by) || length(stratify_by) != 1L) + if (!is.character(stratify_by) || length(stratify_by) != 1L) { stop("`stratify_by` must be NULL or a single string: 'year' or 'country'.") - if (!stratify_by %in% c("year", "country")) + } + if (!stratify_by %in% c("year", "country")) { stop("`stratify_by` must be NULL, 'year', or 'country'.") + } } files <- createMLinputList(path, - stratify_by = stratify_by, - LOO = LOO, - MDR = FALSE, - cross_test = cross_test) + stratify_by = stratify_by, + LOO = LOO, + MDR = FALSE, + cross_test = cross_test + ) if (nrow(files) == 0) { message("No files found to process. Exiting.") @@ -864,8 +898,7 @@ runMLmodels <- function(path, strat_suffix <- if (is.null(stratify_by) || identical(stratify_by, "")) { "" } else { - switch( - stratify_by, + switch(stratify_by, "country" = "_country", "year" = "_year", stop("`stratify_by` must be NULL, 'year', or 'country'.") @@ -874,18 +907,19 @@ runMLmodels <- function(path, # Auto naming for shuffled and PCA shuffle_tag <- if (isTRUE(shuffle_labels)) "shuffled_" else "" - pca_tag <- if (isTRUE(use_pca)) paste0("_pca", as.character(pca_threshold)) else "" + pca_tag <- if (isTRUE(use_pca)) paste0("_pca", as.character(pca_threshold)) else "" results_list <- future.apply::future_lapply( seq_len(nrow(files)), FUN = function(i) { - - ref_parquet <- files$ref_file[i] + ref_parquet <- files$ref_file[i] output_prefix <- files$output_prefix[i] if (interactive()) { - message(sprintf("[runMLmodels] %d/%d: %s", - i, nrow(files), basename(ref_parquet))) + message(sprintf( + "[runMLmodels] %d/%d: %s", + i, nrow(files), basename(ref_parquet) + )) } ml_input <- loadMLInputTibble(ref_parquet) @@ -910,32 +944,37 @@ runMLmodels <- function(path, list(split = split, seed = 5280, n_fold = n_fold) } - res <- try({ - runMLPipeline( - ml_input_tibble = ml_input, - test_data = test_data, - model = "LR", - split = sp$split, - n_fold = sp$n_fold, - prop_vi_top_feats = prop_vi_top_feats, - n_top_feats = NA, - use_pca = use_pca, - pca_threshold = pca_threshold, - shuffle_labels = shuffle_labels, - penalty_vec = 10^seq(-4, -1, length.out = 10), - mix_vec = 0:5 / 5, - select_best_metric = "mcc", - seed = sp$seed, - verbose = verbose, - return_tune_res = return_tune_res, - return_fit = return_fit, - return_pred = return_pred - ) - }, silent = TRUE) + res <- try( + { + runMLPipeline( + ml_input_tibble = ml_input, + test_data = test_data, + model = "LR", + split = sp$split, + n_fold = sp$n_fold, + prop_vi_top_feats = prop_vi_top_feats, + n_top_feats = NA, + use_pca = use_pca, + pca_threshold = pca_threshold, + shuffle_labels = shuffle_labels, + penalty_vec = 10^seq(-4, -1, length.out = 10), + mix_vec = 0:5 / 5, + select_best_metric = "mcc", + seed = sp$seed, + verbose = verbose, + return_tune_res = return_tune_res, + return_fit = return_fit, + return_pred = return_pred + ) + }, + silent = TRUE + ) if (inherits(res, "try-error")) { - warning("Model failed for: ", output_prefix, - "\n Error: ", attr(res, "condition")$message) + warning( + "Model failed for: ", output_prefix, + "\n Error: ", attr(res, "condition")$message + ) return(NULL) } @@ -943,19 +982,25 @@ runMLmodels <- function(path, base <- paste0(shuffle_tag, config_prefix, output_prefix, pca_tag, strat_suffix) if (!is.null(res$performance_tibble)) { - readr::write_tsv(res$performance_tibble, - file.path(files$out_perf[i], paste0(base, "_performance.tsv"))) + readr::write_tsv( + res$performance_tibble, + file.path(files$out_perf[i], paste0(base, "_performance.tsv")) + ) } if (!is.null(res$top_feat_tibble)) { - readr::write_tsv(res$top_feat_tibble, - file.path(files$out_top[i], paste0(base, "_top_features.tsv"))) + readr::write_tsv( + res$top_feat_tibble, + file.path(files$out_top[i], paste0(base, "_top_features.tsv")) + ) } if (!is.null(res$fit)) { saveRDS(res$fit, file.path(files$out_models[i], paste0(base, "_model_fit.rds"))) } if (!is.null(res$pred)) { - readr::write_tsv(res$pred, - file.path(files$out_pred[i], paste0(base, "_prediction.tsv"))) + readr::write_tsv( + res$pred, + file.path(files$out_pred[i], paste0(base, "_prediction.tsv")) + ) } NULL @@ -973,7 +1018,6 @@ runMLmodels <- function(path, } - #' Run the entire AMR ML pipeline from a parquet-backed DuckDB #' #' This function provides a complete end-to-end AMR machine learning workflow. @@ -1006,11 +1050,12 @@ runModelingPipeline <- function(parquet_duckdb_path, pca_threshold = 0.99, verbose = TRUE, use_saved_split = TRUE) { - parquet_duckdb_path <- normalizePath(parquet_duckdb_path) if (!file.exists(parquet_duckdb_path)) { - stop("Parquet-backed DuckDB at ", parquet_duckdb_path, " not found.\n", - "Are you using `{Bug}.duckdb` instead of `{Bug}_parquet.duckdb?`") + stop( + "Parquet-backed DuckDB at ", parquet_duckdb_path, " not found.\n", + "Are you using `{Bug}.duckdb` instead of `{Bug}_parquet.duckdb?`" + ) } out_root <- dirname(parquet_duckdb_path) @@ -1024,9 +1069,9 @@ runModelingPipeline <- function(parquet_duckdb_path, generateMLInputs( parquet_duckdb_path = parquet_duckdb_path, out_path = out_root, - n_fold = n_fold, - split = split, - min_n = min_n, + n_fold = n_fold, + split = split, + min_n = min_n, verbosity = if (verbose) "minimal" else "debug" ) @@ -1089,12 +1134,13 @@ runModelingPipeline <- function(parquet_duckdb_path, # All done! if (verbose) { message("\n=== AMR-ML Pipeline Complete ===") - message("All matrices, models, top feature lists, and performance metrics saved under:\n ", - out_root) + message( + "All matrices, models, top feature lists, and performance metrics saved under:\n ", + out_root + ) message("\nTo inspect model outputs, see directories such as:") message(" ML_performance/, ML_models/, ML_prediction/, ML_top_features/") } invisible(out_root) } - diff --git a/R/run_ml_pipeline.R b/R/run_ml_pipeline.R index 2a97c00..ad9f691 100644 --- a/R/run_ml_pipeline.R +++ b/R/run_ml_pipeline.R @@ -93,20 +93,21 @@ runMLPipeline <- function( .checkArgReturnPred(return_pred) - # Set `n_fold` to `NA` if not using cross-validation. if (split[2] != 0) { n_fold <- NA } # Confirm resolved split params - if (verbose) { - mode <- if (split[2] == 0) "cv" else "splits" - message(sprintf("ML split mode: %s | split = c(%.2f, %.2f) | n_fold = %s | seed = %s", - mode, split[1], split[2], - ifelse(is.na(n_fold), "NA", as.character(n_fold)), - as.character(seed))) - } + if (verbose) { + mode <- if (split[2] == 0) "cv" else "splits" + message(sprintf( + "ML split mode: %s | split = c(%.2f, %.2f) | n_fold = %s | seed = %s", + mode, split[1], split[2], + ifelse(is.na(n_fold), "NA", as.character(n_fold)), + as.character(seed) + )) + } # Create a variable indicating whether external `test_data` was provided. This # will be set to `TRUE` later if the `test_data` argument is not `NA`. @@ -116,10 +117,10 @@ runMLPipeline <- function( # Determine whether multi-class classification is to be performed. if (as.character(.getTargetVarName(ml_input_tibble)) == "resistant_classes") { - multi_class <- TRUE - } else { - multi_class <- FALSE - } + multi_class <- TRUE + } else { + multi_class <- FALSE + } if (model != "LR" & multi_class) { stop(paste( @@ -262,7 +263,7 @@ runMLPipeline <- function( mix_vec = mix_vec ) } - + recipe <- buildRecipe(train_data, use_pca = use_pca, pca_threshold = pca_threshold @@ -421,14 +422,16 @@ runMLPipeline <- function( all_results[["fit"]] <- fit } - if(return_pred) { - if(!multi_class){ + if (return_pred) { + if (!multi_class) { all_results[["pred"]] <- test_data_plus_predictions |> - dplyr::select(c(genome_id, .pred_class, .pred_Resistant, - .pred_Susceptible, genome_drug.resistant_phenotype)) - } - all_results[["pred"]] <- test_data_plus_predictions + dplyr::select(c( + genome_id, .pred_class, .pred_Resistant, + .pred_Susceptible, genome_drug.resistant_phenotype + )) } + all_results[["pred"]] <- test_data_plus_predictions + } return(all_results) } diff --git a/vignettes/intro.Rmd b/vignettes/intro.Rmd index 996eb6b..af5bc8e 100644 --- a/vignettes/intro.Rmd +++ b/vignettes/intro.Rmd @@ -264,19 +264,19 @@ ml_tibble_reduced <- removeTopFeats(ml_tibble, top_features) ### Precision-recall curve ```{r plot-prc} -test_data_plus_predictions <- readr::read_tsv(results/ML_pred/Sfl_drug_AMP_domains_binary_prediction.tsv) +test_data_plus_predictions <- readr::read_tsv(results / ML_pred / Sfl_drug_AMP_domains_binary_prediction.tsv) plotPRC(test_data_plus_predictions) ``` ### ROC curve ```{r plot-roc} -test_data_plus_predictions <- readr::read_tsv(results/ML_pred/Sfl_drug_AMP_domains_binary_prediction.tsv) +test_data_plus_predictions <- readr::read_tsv(results / ML_pred / Sfl_drug_AMP_domains_binary_prediction.tsv) plotROC(test_data_plus_predictions) ``` ### Variable importance plot ```{r plot-vi} -topfeat <- readr::read_tsv(results/ML_top_features/Sfl_drug_AMP_domains_binary_top_features.tsv) +topfeat <- readr::read_tsv(results / ML_top_features / Sfl_drug_AMP_domains_binary_top_features.tsv) plotTopFeatsVI(topfeat) ``` ### Baseline comparison barplot @@ -326,7 +326,6 @@ You can label the top N features to highlight the strongest hits (default is 5) ```{r} plotFishers(fisher_results) plotFishers(fisher_results, alpha = 0.01, label_top_n = 5) - ``` ## Wrapper to run all models @@ -338,14 +337,15 @@ Given a DuckDB file produced by `runDataProcessing()`, it: 5. saves performance metrics, fitted models, predictions, and top feature rankings ``` {r} runModelingPipeline(parquet_duckdb_path, - threads = 16, - n_fold = 5, - split = c(1, 0), - min_n = 25, - prop_vi_top_feats = c(0, 1), - pca_threshold = 0.99, - verbose = TRUE, - use_saved_split = TRUE) + threads = 16, + n_fold = 5, + split = c(1, 0), + min_n = 25, + prop_vi_top_feats = c(0, 1), + pca_threshold = 0.99, + verbose = TRUE, + use_saved_split = TRUE +) ``` Merge the performance and top features of each kind of models into a parquet that will serve as starting data for `amRshiny` package @@ -357,7 +357,7 @@ buildPerformancePq( LOO = FALSE, MDR = FALSE, cross_test = FALSE, - out_parquet = NULL, + out_parquet = NULL, compression = "zstd", verbose = TRUE ) @@ -367,8 +367,8 @@ buildTopFeatsPq( LOO = FALSE, MDR = FALSE, cross_test = FALSE, - out_parquet = NULL, + out_parquet = NULL, compression = "zstd", verbose = TRUE -) +) ``` From d3a97bfff946b1be6abb2cf1b9e3067084bea6f4 Mon Sep 17 00:00:00 2001 From: Abhirupa Ghosh <100681585+AbhirupaGhosh@users.noreply.github.com> Date: Fri, 29 May 2026 16:39:23 -0600 Subject: [PATCH 03/17] Refactor plot functions to accept file inputs --- R/plot_ml.R | 58 +++++++++++++++++++---------------------------------- 1 file changed, 21 insertions(+), 37 deletions(-) diff --git a/R/plot_ml.R b/R/plot_ml.R index 42d51ec..a7b459a 100644 --- a/R/plot_ml.R +++ b/R/plot_ml.R @@ -24,8 +24,8 @@ NULL #' Plot a Precision-Recall Curve #' #' Generates a precision-recall curve (PRC) for AMR phenotype prediction results. -#' @param test_data_plus_predictions A tibble containing test data with added -#' prediction columns, typically the output of `runMLmodels()`. +#' @param test_data_plus_predictions_file A file containing test data with added +#' prediction columns, typically the output of `runMLmodels(return_pred=TRUE)`. #' #' @return A `ggplot2` object showing the precision-recall curve. #' @@ -35,12 +35,13 @@ NULL #' #' @examples #' \dontrun{ -#' test_data_plus_predictions <- readr::read_tsv(results / ML_pred / Sfl_drug_AMP_domains_binary_prediction.tsv) -#' plotPRC(test_data_plus_predictions) +#' test_data_plus_predictions_file <- "results / ML_pred / Sfl_drug_AMP_domains_binary_prediction.tsv" +#' plotPRC(test_data_plus_predictions_file) #' } #' #' @export -plotPRC <- function(test_data_plus_predictions) { +plotPRC <- function(test_data_plus_predictions_file) { + test_data_plus_predictions <- readr::read_tsv(test_data_plus_predictions_file) .checkArgTestDataPlusPredictions(test_data_plus_predictions) test_data_plus_predictions <- test_data_plus_predictions |> dplyr::mutate( @@ -66,13 +67,14 @@ plotPRC <- function(test_data_plus_predictions) { #' #' Generates a ROC curve for AMR phenotype prediction results. #' -#' @param test_data_plus_predictions A tibble with test data and prediction -#' columns (output of `runMLmodels()`). +#' @param test_data_plus_predictions_file A file with test data and prediction +#' columns (output of `runMLmodels(return_pred=TRUE)`). #' #' @return A ROC curve plotted using `ggplot2::autoplot()`. #' #' @export -plotROC <- function(test_data_plus_predictions) { +plotROC <- function(test_data_plus_predictions_file) { + test_data_plus_predictions <- readr::read_tsv(test_data_plus_predictions_file) .checkArgTestDataPlusPredictions(test_data_plus_predictions) test_data_plus_predictions <- test_data_plus_predictions |> dplyr::mutate( @@ -96,13 +98,14 @@ plotROC <- function(test_data_plus_predictions) { #' #' Produces a heatmap visualization of the confusion matrix for AMR predictions. #' -#' @param test_data_plus_predictions A tibble containing true and predicted +#' @param test_data_plus_predictions_file A file containing true and predicted #' phenotype labels. #' #' @return A heatmap (`ggplot2` object) showing the confusion matrix. #' #' @export -plotCM <- function(test_data_plus_predictions) { +plotCM <- function(test_data_plus_predictions_file) { + test_data_plus_predictions <- readr::read_tsv(test_data_plus_predictions_file) .checkArgTestDataPlusPredictions(test_data_plus_predictions) test_data_plus_predictions <- test_data_plus_predictions |> dplyr::mutate( @@ -123,32 +126,12 @@ plotCM <- function(test_data_plus_predictions) { ggplot2::autoplot(type = "heatmap") } -#' Plot Density of Predicted Class Probabilities -#' -#' Visualizes how predicted class probabilities differ between resistant and -#' susceptible genome-drug combinations. -#' -#' @param test_data_plus_predictions Tibble with prediction probabilities and -#' true labels. -#' -#' @return A ggplot2 density plot. -#' -#' @export -plotDensity <- function(test_data_plus_predictions) { - test_data_plus_predictions |> - ggplot2::ggplot(ggplot2::aes( - x = .pred_Resistant, - fill = genome_drug.resistant_phenotype - )) + - ggplot2::geom_density(alpha = 0.5) -} - #' Plot Top Feature Importances #' #' Creates a bar plot showing the most important features affecting #' AMR phenotype predictions. #' -#' @param topfeat A tibble containing feature importance scores +#' @param topfeat_file A file containing feature importance scores #' (output of `runMLmodels()`). #' @param n_top_feats Number of top features to display (default: 10). #' @@ -156,12 +139,13 @@ plotDensity <- function(test_data_plus_predictions) { #' #' @examples #' \dontrun{ -#' topfeat <- readr::read_tsv(results / ML_top_features / Sfl_drug_AMP_domains_binary_top_features.tsv) -#' plotTopFeatsVI(topfeat) +#' topfeat_file <- "results / ML_top_features / Sfl_drug_AMP_domains_binary_top_features.tsv" +#' plotTopFeatsVI(topfeat_file) #' } #' #' @export -plotTopFeatsVI <- function(topfeat, n_top_feats = 10) { +plotTopFeatsVI <- function(topfeat_file, n_top_feats = 10) { + topfeat <- readr::read_tsv(topfeat_file) .checkArgNTopFeats(n_top_feats) vip <- topfeat |> @@ -279,8 +263,8 @@ plotFishers <- function( stop("Missing required columns: ", paste(missing_cols, collapse = ", ")) } - plot_df <- fisher_df %>% - dplyr::arrange(adj_p_value) %>% + plot_df <- fisher_df |> + dplyr::arrange(adj_p_value) |> dplyr::mutate( rank = dplyr::row_number(), neg_log10_adj_p = -log10(adj_p_value), @@ -312,7 +296,7 @@ plotFishers <- function( ) if (label_top_n > 0) { - label_df <- plot_df %>% + label_df <- plot_df |> dplyr::slice_head(n = label_top_n) p <- p + From 241836639dc00c8d4e904242d6428d35064a9f7b Mon Sep 17 00:00:00 2001 From: AbhirupaGhosh Date: Fri, 29 May 2026 22:42:41 +0000 Subject: [PATCH 04/17] Style code (GHA) --- R/plot_ml.R | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/R/plot_ml.R b/R/plot_ml.R index a7b459a..6ee976d 100644 --- a/R/plot_ml.R +++ b/R/plot_ml.R @@ -41,7 +41,7 @@ NULL #' #' @export plotPRC <- function(test_data_plus_predictions_file) { - test_data_plus_predictions <- readr::read_tsv(test_data_plus_predictions_file) + test_data_plus_predictions <- readr::read_tsv(test_data_plus_predictions_file) .checkArgTestDataPlusPredictions(test_data_plus_predictions) test_data_plus_predictions <- test_data_plus_predictions |> dplyr::mutate( @@ -74,7 +74,7 @@ plotPRC <- function(test_data_plus_predictions_file) { #' #' @export plotROC <- function(test_data_plus_predictions_file) { - test_data_plus_predictions <- readr::read_tsv(test_data_plus_predictions_file) + test_data_plus_predictions <- readr::read_tsv(test_data_plus_predictions_file) .checkArgTestDataPlusPredictions(test_data_plus_predictions) test_data_plus_predictions <- test_data_plus_predictions |> dplyr::mutate( @@ -105,7 +105,7 @@ plotROC <- function(test_data_plus_predictions_file) { #' #' @export plotCM <- function(test_data_plus_predictions_file) { - test_data_plus_predictions <- readr::read_tsv(test_data_plus_predictions_file) + test_data_plus_predictions <- readr::read_tsv(test_data_plus_predictions_file) .checkArgTestDataPlusPredictions(test_data_plus_predictions) test_data_plus_predictions <- test_data_plus_predictions |> dplyr::mutate( From f2eb43c427d87e3f732505aa6113665f4f885713 Mon Sep 17 00:00:00 2001 From: Abhirupa Ghosh <100681585+AbhirupaGhosh@users.noreply.github.com> Date: Mon, 1 Jun 2026 17:45:56 -0600 Subject: [PATCH 05/17] Implement drug performance plotting functions in plot_ml.R Add functions to plot drug phenotype distribution, model performance, and cross-drug generalization heatmaps. --- R/plot_ml.R | 798 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 798 insertions(+) diff --git a/R/plot_ml.R b/R/plot_ml.R index 6ee976d..fab0759 100644 --- a/R/plot_ml.R +++ b/R/plot_ml.R @@ -313,3 +313,801 @@ plotFishers <- function( return(p) } + +#' Plot drug phenotype distribution +#' +#' Reads metadata and generates a stacked bar plot showing counts of resistant +#' and susceptible phenotypes per antibiotic. +#' +#' @param metadata_path Character. Path to directory containing `metadata.parquet`. +#' +#' @return A ggplot object. +#' @export +#' +#' @examples +#' plotDrugDist(metadata_path = "data/Campylobacter/") +plotDrugDist <- function(metadata_path = "."){ + +metadata <- arrow::read_parquet(file.path(metadata_path,"metadata.parquet")) + +##################### phenotype distribution (drugs) ######################### +drug_dist <- metadata |> + dplyr::distinct(genome.genome_id, + genome_drug.antibiotic, + drug_abbr, + genome_drug.resistant_phenotype) |> + dplyr::count(genome_drug.antibiotic, + drug_abbr, + genome_drug.resistant_phenotype) |> + dplyr::group_by(genome_drug.antibiotic, drug_abbr) |> + dplyr::mutate(total = sum(n)) |> + dplyr::ungroup() |> + dplyr::mutate( + label = paste0(genome_drug.antibiotic, " (", drug_abbr, ")"), + label = forcats::fct_reorder(label, total) + ) + +p <- ggplot2::ggplot( + drug_dist, + ggplot2::aes( + x = label, + y = n, + fill = genome_drug.resistant_phenotype + ) +) + + ggplot2::geom_col(color = "black", width = 0.8) + + ggplot2::coord_flip() + + ggplot2::scale_fill_manual( + values = c( + "Resistant" = "#d4872a", + "Susceptible" = "#5b8db8" + ), + name = "Phenotype" + ) + + ggplot2::labs( + x = "Antibiotic", + y = "Number of unique genomes" + ) + + ggplot2::theme_classic(base_size = 14) + +p +} + +#' Plot drug-level model performance +#' +#' Generates heatmaps and ridge plots summarizing model performance (MCC) +#' across drugs and feature types. +#' +#' @param metadata_path Character. Path to `metadata.parquet`. +#' @param performance_path Character. Path to `all_performance.parquet`. +#' +#' @return A patchwork ggplot object combining multiple panels. +#' @export +#' +#' @examples +#' plotDrugPerf(metadata_path = "data/Campylobacter/", performance_path = "data/Campylobacter/ML_performance/") +plotDrugPerf <- function(metadata_path = ".", performance_path = ".") { + +metadata <- arrow::read_parquet(file.path(metadata_path,"metadata.parquet")) + +performance <- arrow::read_parquet(file.path(performance_path,"all_performance.parquet")) + +######################## drug performances ################################# +median_drug <- performance |> + dplyr::filter( + drug_label == "drug", + !shuffled # keep real models; remove if you want both + ) |> + dplyr::group_by(drug_or_class, feature_type, feature_subtype) |> + dplyr::summarise(median_mcc = median(mcc, na.rm = TRUE), .groups = "drop") |> + dplyr::left_join(plot_df, by = c("drug_or_class" = "drug_abbr")) |> + dplyr::mutate(drug_or_class = reorder(drug_or_class, total)) + +drug_p1 <- ggplot2::ggplot(median_drug, + ggplot2::aes(x = feature_type, + y = drug_or_class, + fill = median_mcc)) + + ggplot2::geom_tile(color = "grey90", width = 0.9) + + + ggplot2::scale_fill_gradientn( + colors = c( + "#C4B8A8", # low + "#FAFAF7", # around 0 + "#5F84C9", # medium/high (~0.7–0.9) + "#0F2A5A" # very dark for ~1 + ), + values = scales::rescale(c(-1, 0, 0.85, 1)), + name = "Best MCC" + ) + + + ggplot2::labs(x = "Feature type") + + ggplot2::theme_minimal(base_size = 12) + + ggplot2::theme( + axis.text = ggplot2::element_text(size = 10, colour = "black"), + axis.title = ggplot2::element_text(size = 12), + axis.title.y = ggplot2::element_blank(), + axis.text.x = ggplot2::element_text(angle = 45, hjust = 1), + legend.position = "bottom" + ) + + ggplot2::coord_fixed() + +drug_p1 + +median_feature <- performance |> + dplyr::filter( + drug_label == "drug", + !shuffled + ) |> + dplyr::group_by(drug_or_class, feature_type) |> + dplyr::summarise(median_mcc = median(mcc, na.rm = TRUE), .groups = "drop") |> + dplyr::left_join(plot_df, by = c("drug_or_class" = "drug_abbr")) |> + dplyr::mutate(drug_or_class = reorder(drug_or_class, total)) + + +feat_pal <- c( + "args" = "#56B4E9", # sky blue + "cogs" = "#E69F00", # orange + "genes" = "#009E73", # bluish green + "domains" = "#F0E442", # yellow + "proteins" = "#CC79A7", # reddish purple + "struct" = "#D55E00" # vermillion +) + +rc_perf <- ggplot2::ggplot(median_feature |> + dplyr::distinct(drug_or_class, + feature_type, median_mcc), + ggplot2::aes(x = median_mcc, y = drug_or_class)) + + + ggridges::geom_density_ridges( + scale = 0.75, + rel_min_height = 0.01, + alpha = 0.4, + fill = "grey90", + colour = "grey70" + ) + + + ggplot2::geom_point( + position = position_jitter(height = 0.1), + size = 2, + alpha = 0.8, + aes(color = feature_type) + ) + + ggplot2::scale_color_manual(values = feat_pal, name = "Feature type") + + ggplot2::stat_summary( + fun = median, + geom = "point", + size = 2, + color = "black" + ) + + ggplot2::theme_minimal(base_size = 14) + + ggplot2::theme( + axis.text = ggplot2::element_text(size = 10, colour = "black"), + axis.title = ggplot2::element_text(size = 12), + axis.title.y = ggplot2::element_blank(), + axis.text.x = ggplot2::element_text(angle = 45, hjust = 1), + legend.position = "right", + panel.grid.minor = ggplot2::element_blank(), + axis.line = ggplot2::element_line(color = "black") + ) + +rc_perf + +final_plot = drug_p1 + + rc_perf + + patchwork::plot_layout( + widths = c(2, 2), # adjust proportions + guides = "collect" + ) & + ggplot2::theme( + legend.position = "bottom" + ) + +final_plot + +} + +#' Plot cross-drug generalization heatmap +#' +#' Creates a heatmap showing cross-drug model performance (MCC), where models +#' trained on one drug are evaluated on another. +#' +#' @param cross_test_performance_path Character. Path to `cross_drug_perf.parquet`. +#' @param drug_performance_path Character. Path to `all_performance.parquet`. +#' +#' @return A ComplexHeatmap object. +#' @export +#' +#' @examples +#' plotCrossDrug(cross_test_performance_path = "data/Campylobacter/cross_test_ML_performance", drug_performance_path = "data/Campylobacter/ML_performance/") +plotCrossDrug <- function(cross_test_performance_path = ".", drug_performance_path = ".") { + + cross_drug <- arrow::read_parquet(file.path(cross_test_performance_path,"cross_drug_perf.parquet")) + performance <- arrow::read_parquet(file.path(drug_performance_path,"all_performance.parquet")) + +###################### CROSS DRUG Testing ############################# +heatmap_df <- cross_drug |> + # dplyr::filter(tested_on %in% (cross_drug |> dplyr::pull(drug_or_class))) |> + dplyr::group_by(drug_or_class, tested_on) |> + dplyr::summarise(median_mcc = median(mcc, na.rm = TRUE), .groups = "drop") + +same_drugs <- performance |> + dplyr::filter(drug_label == "drug", + drug_or_class %in% (cross_drug |> + dplyr::distinct(drug_or_class) |> + dplyr::pull())) |> + dplyr::group_by(drug_or_class) |> + dplyr::summarise(median_mcc = median(mcc, na.rm = TRUE), .groups = "drop") |> + dplyr::mutate(tested_on = drug_or_class) |> + dplyr::distinct(drug_or_class, tested_on, median_mcc) + +heatmap_df <- heatmap_df |> + dplyr::add_row(same_drugs) |> + dplyr::left_join(metadata |> + dplyr::distinct(drug_abbr, class_abbr), + by = c("drug_or_class" = "drug_abbr")) |> + dplyr::rename("drug_class" = "class_abbr") |> + dplyr::left_join(metadata |> + dplyr::distinct(drug_abbr, class_abbr), + by = c("tested_on" = "drug_abbr")) + +# Row annotation (already similar to what you did) +annotation_row <- heatmap_df |> + dplyr::distinct(drug_or_class, drug_class) |> + tibble::column_to_rownames("drug_or_class") + +# Column annotation +annotation_col <- heatmap_df |> + dplyr::distinct(tested_on, class_abbr) |> + tibble::column_to_rownames("tested_on") + +mat <- heatmap_df |> + dplyr::select(drug_or_class, tested_on, median_mcc) |> + tidyr::pivot_wider(names_from = tested_on, values_from = median_mcc) |> + tibble::column_to_rownames("drug_or_class") |> + as.matrix() + +row_order <- heatmap_df |> + dplyr::distinct(drug_or_class, drug_class) |> + dplyr::arrange(drug_class, drug_or_class) |> + dplyr::pull(drug_or_class) + +col_order <- heatmap_df |> + dplyr::distinct(tested_on, class_abbr) |> + dplyr::arrange(class_abbr, tested_on) |> + dplyr::pull(tested_on) + +# mat[is.na(mat)] <- 0 +mat <- mat[row_order, col_order] + +# Align annotations +annotation_row <- annotation_row[row_order, , drop = FALSE] +annotation_col <- annotation_col[col_order, , drop = FALSE] + + +# Collect all classes from both row and column +classes <- base::union( + annotation_row$drug_class, + annotation_col$class_abbr +) + +# Create ONE named color vector +class_colors <- stats::setNames( + scales::hue_pal() (length(classes)), + classes +) + +heat_colors <- colorRampPalette(RColorBrewer::brewer.pal(11, "RdBu"))(100) + +# ---- Convert annotations ---- +ha_row <- ComplexHeatmap::rowAnnotation( + drug_class = annotation_row$drug_class, + col = list(drug_class = class_colors), + show_annotation_name = FALSE, + show_legend = FALSE +) + +ha_col <- ComplexHeatmap::HeatmapAnnotation( + class_abbr = annotation_col$class_abbr, + col = list(class_abbr = class_colors), + show_annotation_name = FALSE, na_col = "grey3" +) + +# ---- Color function (instead of breaks + palette) ---- +col_fun <- circlize::colorRamp2( + seq(-max_val, max_val, length.out = length(heat_colors)), + heat_colors +) +# ---- Heatmap ---- +cross_drug_hm <- ComplexHeatmap::Heatmap( + mat, + name = "median_mcc", + col = col_fun, + cluster_rows = FALSE, + cluster_columns = FALSE, + row_order = row_order, + column_order = col_order, + left_annotation = ha_row, + top_annotation = ha_col, + show_row_names = TRUE, + show_column_names = TRUE, + column_title = "tested on", + row_title = "trained on", + column_title_side = "bottom", + row_title_side = "right", + row_names_gp = grid::gpar(fontsize = 14), + column_names_gp = grid::gpar(fontsize = 14), + column_names_rot = 0, + + # remove borders like pheatmap + rect_gp = grid::gpar(col = NA), + + # legends + show_heatmap_legend = TRUE +) + +cross_drug_hm +} + +#' Plot stratified model performance +#' +#' Visualizes model performance (MCC) stratified by year or country, +#' comparing within-group vs cross-group evaluation. +#' +#' @param year_or_country Character. Either "year" or "country". +#' @param stratified_performance_path Character. Path to stratified performance files. +#' @param stratified_cross_performance_path Character. Path to cross-stratified performance files. +#' +#' @return A ggplot object. +#' @export +#' +#' @examples +#' plotStratifiedPerf("year", stratified_performance_path = "data/Campylobacter/ML_year_performance", +#' stratified_cross_performance_path = "data/Campylobacter/cross_test_ML_year_performance") +plotStratifiedPerf <- function(year_or_country = "year", + stratified_performance_path = ".", + stratified_cross_performance_path = ".") { + + perf <- arrow::read_parquet(file.path(stratified_performance_path, + paste0(year_or_country,"_perf.parquet"))) + + cross_test <- arrow::read_parquet(file.path(stratified_cross_performance_path, + paste0("cross_", + year_or_country, + "_perf.parquet"))) +if(year_or_country == "year") { +all <- perf |> + dplyr::rename("train_year" = "strat_value") |> + dplyr::mutate(test_year = train_year) |> + dplyr::select(drug_label, drug_or_class, + train_year, test_year, feature_type, feature_subtype, mcc) |> + dplyr::bind_rows(cross_test |> + dplyr::select(drug_label, drug_or_class, + train_year, test_year, feature_type, + feature_subtype, mcc)) |> + dplyr::mutate(category = dplyr::if_else( + train_year == test_year, "same year bin", "different year bin")) +} +else { +all <- perf |> + dplyr::rename("train_country" = "strat_value") |> + dplyr::mutate(test_country = train_country) |> + dplyr::select(drug_label, drug_or_class, + train_country, test_country, + feature_type, feature_subtype, mcc) |> + dplyr::bind_rows(cross_test |> + dplyr::select(drug_label, drug_or_class, + train_country, test_country, + feature_type, feature_subtype, mcc)) |> + dplyr::mutate(category = dplyr::if_else( + train_country == test_country, "same country", "different country")) +} + + fill_vals <- if (year_or_country == "year") { + c( + "same year bin" = "#b3cde3", + "different year bin" = "#fbb4ae" + ) + } else { + c( + "same country" = "#b3cde3", + "different country" = "#fbb4ae" + ) + } + + plot <- ggplot2::ggplot( + all |> + dplyr::filter(drug_label == "drug", !is.na(mcc)), + ggplot2::aes(x = mcc, y = drug_or_class, fill = category) + ) + + ggridges::geom_density_ridges( + alpha = 0.5, + scale = 1, + rel_min_height = 0.01, + position = "identity" + ) + + ggplot2::geom_vline(xintercept = 0, linetype = "dashed", color = "black") + + ggplot2::scale_fill_manual(values = fill_vals) + + ggplot2::theme_minimal(base_size = 14) + + ggplot2::labs( + title = if (year_or_country == "year") { + "Temporal performance by drug" + } else { + "Geographical performance by drug" + }, + x = "MCC", + y = "Drug", + fill = "Tested on" + ) + + ggplot2::theme( + axis.title = ggplot2::element_text(colour = "black", size = 10), + axis.text.x = ggplot2::element_text(angle = 45, hjust = 1, size = 10, colour = "black"), + axis.text.y = ggplot2::element_text(size = 10, colour = "black"), + axis.title.y = ggplot2::element_blank(), + legend.title = ggplot2::element_text(size = 12), + legend.text = ggplot2::element_text(size = 10), + legend.position = "bottom", + plot.title = ggplot2::element_text(face = "bold"), + panel.grid.minor = ggplot2::element_blank(), + plot.margin = margin(0, 0, 0, 0) + ) +plot +} + +#' Plot multi-drug resistance (MDR) model performance +#' +#' Generates violin plots of performance, feature importance summaries, +#' and prediction confusion-style visualizations for MDR models. +#' +#' @param MDR_performance_path Character. Path to `MDR_perf.parquet`. +#' @param MDR_top_feature_path Character. Path to `MDR_top_features.parquet`. +#' @param MDR_pred_path Character. Path to `MDR_pred.parquet`. +#' +#' @return A list of ggplot objects. +#' @export +#' +#' @examples +#' plotMDR(MDR_performance_path = "data/Campylobacter/MDR_ML_performance", MDR_top_feature_path = "data/Campylobacter/MDR_ML_top_features", +#' MDR_pred_path = "data/Campylobacter/MDR_ML_pred") +plotMDR <- function(MDR_performance_path = ".", MDR_top_feature_path = ".", + MDR_pred_path = ".") { + +MDR_perf <- arrow::read_parquet(file.path(MDR_performance_path, "MDR_perf.parquet")) + +# ---- Violin plot ---- +perf_plot <- ggplot2::ggplot(MDR_perf, + ggplot2::aes(x = feature_type, y = mcc)) + + + # violins (overall distribution per feature type) + ggplot2::geom_violin(fill = "grey85", color = NA, alpha = 0.8) + + + # points (colored by binary vs counts) + ggplot2::geom_jitter( + ggplot2::aes(color = feature_subtype), + width = 0.12, size = 2, alpha = 0.8 + ) + + + ggplot2::scale_color_manual(values = c( + "binary" = "#7B9CB5", + "counts" = "#CC8644" + )) + + + ggplot2::theme_minimal(base_size = 12) + + ggplot2::labs( + # title = "MDR model performances", + # subtitle = "Violin = distribution per feature type; points = binary vs counts", + x = "Feature type", + y = "MCC", + color = "Feature\nsubtype" + ) + + + ggplot2::theme( + legend.position = "right", + plot.title = ggplot2::element_text(face = "bold") + ) + + ggplot2::theme( + axis.title = ggplot2::element_text(colour = "black", size = 10), + axis.text.x = ggplot2::element_text(angle = 45, hjust = 1, size = 10, colour = "black"), + axis.text.y = ggplot2::element_text(size = 14, colour = "black"), + legend.title = ggplot2::element_text(size = 12), + legend.text = ggplot2::element_text(size = 10), + legend.position = "none", + title = ggplot2::element_text(face = "bold"), + + panel.background = ggplot2::element_blank(), + panel.grid.minor = ggplot2::element_blank(), + panel.grid.major.x = ggplot2::element_blank(), # remove vertical lines + panel.grid.major.y = ggplot2::element_line(color = "grey80"), # keep horizontal lines + + axis.line = ggplot2::element_line(color = "black") + ) + + ggplot2::scale_y_continuous(limits = c(0, 1)) + +perf_plot + +MDR_pred <- arrow::read_parquet(file.path(MDR_pred_path, "MDR_pred.parquet")) |> + dplyr::mutate( + diff_top2 = purrr::pmap_dbl(dplyr::across(dplyr::contains(".pred") & dplyr::where(is.numeric)), function(...) { + x <- c(...) + sx <- sort(x, decreasing = TRUE) + sx[1] - sx[2] + }) # Difference between prediction probabilities of top two classes + ) |> + dplyr::select(genome_id, resistant_classes, .pred_class, diff_top2, + feature_type, feature_subtype, seed) |> + dplyr::group_by(resistant_classes, .pred_class, feature_type) |> + dplyr::summarise(mean_margin = mean(diff_top2), n = n(), .groups = "drop") |> + dplyr::group_by(resistant_classes, feature_type) |> # normalize within true class + dplyr::mutate(sum = sum(n), prop = n / sum(n)) |> + dplyr::ungroup() + +MDR_pred_plot <- ggplot2::ggplot(MDR_pred, + ggplot2::aes(x = resistant_classes, + y = .pred_class)) + + ggplot2::geom_tile(ggplot2::aes(fill = prop)) + + ggplot2::geom_point(ggplot2::aes(size = mean_margin), color = "black") + + ggplot2::facet_wrap(~ feature_type) + + ggplot2::scale_fill_distiller( + palette = "RdBu", + direction = 1, # flip with -1 if needed + name = "Prediction proportion" + ) + + ggplot2::labs(x = "true class", y = "predicted class") + + ggplot2::scale_size(range = c(1, 6), name = "Mean margin") + + ggplot2::coord_equal() + + ggplot2::theme_minimal() + + ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45, hjust = 1)) + + ggplot2::theme( + axis.title = ggplot2::element_text(colour = "black", size = 10), + axis.text.x = ggplot2::element_text(angle = 45, hjust = 1, size = 10, colour = "black"), + axis.text.y = ggplot2::element_text(size = 10, colour = "black"), + legend.title = ggplot2::element_text(size = 12), + legend.text = ggplot2::element_text(size = 10), + legend.position = "right", + title = ggplot2::element_text(face = "bold") + ) + +MDR_pred_plot + +# MDR_feat <- arrow::read_parquet(file.path( +# MDR_top_feature_path,"MDR_top_features.parquet")) |> +# pivot_longer(-c(Variable, feature_type, feature_subtype, seed), +# values_to = "Importance", +# names_to = "Resistant_classes") |> +# filter(!Importance == 0) +# +# MDR_feat_clean <- MDR_feat |> +# dplyr::filter(feature_type != "struct") |> +# dplyr::group_by(Resistant_classes, feature_type, feature_subtype, seed) |> +# dplyr::slice_max(Importance, n = top_n, with_ties = FALSE) |> +# dplyr::ungroup() |> +# dplyr::mutate(Variable = gsub( ".NCBIFAM", "", Variable)) |> +# dplyr::mutate(Variable = gsub("^X", "", Variable)) |> +# dplyr::mutate(Variable = dplyr::if_else( +# feature_type == "domains", gsub("_.*", "", Variable), Variable)) |> +# dplyr::mutate(Variable = dplyr::if_else( +# feature_type == "proteins", gsub("fig.", "fig|", Variable), Variable)) |> +# dplyr::left_join(cluster_feature, by = c("Variable" = "feature")) |> +# dplyr::mutate( +# cluster = dplyr::coalesce(cluster, Variable) +# ) +# +# cluster_df <- MDR_feat_clean |> +# dplyr::group_by(Resistant_classes, cluster) |> +# dplyr::summarise( +# Importance = median(Importance, na.rm = TRUE), +# .groups = "drop" +# ) +# +# top_clusters <- cluster_df |> +# group_by(Resistant_classes) |> +# group_modify(~{ +# +# df <- .x +# +# top_pos <- df |> +# arrange(desc(Importance)) |> +# slice_head(n = 10) +# +# top_neg <- df |> +# arrange(Importance) |> +# slice_head(n = 10) +# +# bind_rows(top_pos, top_neg) +# }) |> +# ungroup() +# +# top_clusters <- top_clusters |> +# dplyr::left_join(protein_names, by = c("cluster" = "proteinID")) |> +# dplyr::mutate( +# proteinName = dplyr::coalesce(proteinName, cluster), # fallback +# proteinName = stringr::str_trunc(proteinName, 50) +# ) |> +# dplyr::distinct(Resistant_classes, proteinName, Importance) |> +# # ✅ reorder AFTER naming +# dplyr::group_by(Resistant_classes) |> +# dplyr::mutate( +# proteinName = forcats::fct_reorder(proteinName, Importance) +# ) |> +# dplyr::ungroup() +# +# ggplot(top_clusters, +# aes(x = Importance, y = proteinName)) + +# +# # line (lollipop stem) +# geom_segment( +# aes(x = 0, xend = Importance, +# y = proteinName, yend = proteinName), +# color = "grey60" +# ) + +# +# # dot +# geom_point( +# aes(color = Importance > 0), +# size = 3 +# ) + +# +# facet_wrap(~ Resistant_classes, scales = "free_y") + +# +# scale_color_manual( +# values = c("TRUE" = "#5b8db8", # positive +# "FALSE" = "#d4872a"), # negative +# guide = "none" +# ) + +# +# theme_minimal(base_size = 13) + +# labs( +# x = "Median importance", +# y = "Cluster" +# ) + +# theme( +# panel.grid.minor = element_blank(), +# strip.text = element_text(face = "bold") +# ) + + +} + +#' Compare shuffled vs real model performance +#' +#' Creates boxplots comparing performance (MCC) between real and shuffled labels +#' across feature types. +#' +#' @param metadata_path Character. Path to `metadata.parquet`. +#' @param performance_path Character. Path to `all_performance.parquet`. +#' +#' @return A ggplot object. +#' @export +#' +#' @examples +#' plotShuffleVsReal(metadata_path = "data/Campylobacter/", performance_path = "data/Campylobacter/ML_performance") +plotShuffleVsReal <- function(metadata_path = ".", performance_path = ".") { + + metadata <- arrow::read_parquet(file.path(metadata_path,"metadata.parquet")) + performance <- arrow::read_parquet(file.path(performance_path,"all_performance.parquet")) + + performance |> + dplyr::mutate( + shuffled_label = dplyr::if_else(shuffled, "shuffled", "real") + ) |> + ggplot2::ggplot(ggplot2::aes(x = feature_subtype, y = mcc, fill = shuffled_label)) + + ggplot2::geom_boxplot( + width = 0.55, outlier.size = 0.8, outlier.alpha = 0.4, + outlier.color = "grey50", linewidth = 0.4 + ) + + ggplot2::geom_hline(yintercept = 0, linetype = "dashed", color = "grey60", linewidth = 0.4) + + ggplot2::scale_fill_manual( + values = c("real" = "#7B9CB5", "shuffled" = "#C4B8A8"), + name = NULL + ) + + ggplot2::scale_y_continuous(limits = c(-0.2, 1), breaks = seq(-0.2, 1, 0.2)) + + ggplot2::facet_wrap(~ feature_type, nrow = 1) + + ggplot2::theme_minimal(base_size = 12) + + theme( + panel.grid.major.x = element_blank(), + panel.grid.minor = element_blank(), + panel.grid.major.y = element_line(color = "#E5E2D9", linewidth = 0.4), + strip.text = element_text(color = "grey30", face = "bold", size = 10), + # strip.background = element_rect(fill = "#EEEAE0", color = NA), + axis.text = element_text(color = "grey45"), + legend.position = "top", + legend.text = element_text(color = "grey40", size = 10) + ) + + labs( + x = NULL, y = "MCC" + ) + +} + +#' Plot top contributing feature clusters +#' +#' Identifies top contributing clusters across feature types and drugs, +#' and visualizes their relative contributions. +#' +#' @param top_feat_path Character. Path to `all_top_features.parquet`. +#' @param cluster_feature_path Character. Path to `cluster_feature.parquet`. +#' @param protein_names_path Character. Path to `protein_names.parquet`. +#' @param top_n Integer. Number of top features to retain per model. +#' +#' @return A ggplot object. +#' @export +#' +#' @examples +#' plotTopClusters(top_feat_path = "data/Campylobacter/ML_top_features", cluster_feature_path = "data/Campylobacter/", +#' protein_names_path = "data/Campylobacter/", top_n = 10) +plotTopClusters <- function(top_feat_path = ".", cluster_feature_path = ".", + protein_names_path = ".", top_n = 10) { + + ################### Top features ######################### + + top_feat <- arrow::read_parquet(file.path(top_feat_path, "all_top_features.parquet")) + cluster_feature <- arrow::read_parquet(file.path(cluster_feature_path, "cluster_feature.parquet")) + protein_names <- arrow::read_parquet(file.path(protein_names_path, "protein_names.parquet")) + + # which clusters appear in top n across feature types per drug + # join top features with cluster mapping, filter out struct and shuffled + top_feat_clean <- top_feat |> + dplyr::filter(!shuffled, feature_type != "struct", drug_label == "drug") |> + dplyr::group_by(drug_or_class, feature_type, feature_subtype, seed) |> + dplyr::slice_max(Importance, n = top_n, with_ties = FALSE) |> + dplyr::ungroup() |> + dplyr::mutate(Variable = gsub( ".NCBIFAM", "", Variable)) |> + dplyr::mutate(Variable = gsub("^X", "", Variable)) |> + dplyr::mutate(Variable = dplyr::if_else( + feature_type == "domains", gsub("_.*", "", Variable), Variable)) |> + dplyr::mutate(Variable = dplyr::if_else( + feature_type == "proteins", gsub("fig.", "fig|", Variable), Variable)) |> + dplyr::left_join(cluster_feature, by = c("Variable" = "feature")) |> + dplyr::mutate( + cluster = dplyr::coalesce(cluster, Variable), # fallback to Variable if no match + Importance_signed = dplyr::if_else(Sign == "NEG", -Importance, Importance) + ) + + shared_mat <- top_feat_clean |> + dplyr::group_by(drug_or_class, feature_type, cluster) |> + dplyr::summarise(abs_imp = median(Importance, na.rm = TRUE), .groups = "drop") |> + + # convert to contribution within each feature_type + dplyr::group_by(drug_or_class, feature_type) |> + dplyr::mutate(contribution = abs_imp / sum(abs_imp, na.rm = TRUE)) |> + + # pick top n contributors + dplyr::slice_max(contribution, n = top_n, with_ties = FALSE) |> + dplyr::ungroup() |> + + dplyr::add_count(drug_or_class, cluster, name = "n_feat_types") |> + dplyr::left_join(protein_names, by = c("cluster" = "proteinID")) |> + dplyr::mutate( + proteinName = stringr::str_trunc(proteinName, 50), + proteinName = forcats::fct_reorder(proteinName, n_feat_types) + ) + + feat_plot <- ggplot2::ggplot(shared_mat, + ggplot2::aes(x = feature_type, + y = proteinName, + fill = contribution)) + + ggplot2::geom_tile(color = "#FAFAF7", linewidth = 0.5, width = 0.9, height = 0.9) + + # coord_fixed() + + ggplot2::scale_fill_distiller( + palette = "RdPu", + direction = 1, + name = "contribution", + na.value = "#EEEAE0" + ) + + ggplot2::facet_wrap(~ drug_or_class, scales = "free_y") + + ggplot2::theme_minimal(base_size = 12) + + ggplot2::theme( + panel.grid.major.x = ggplot2::element_blank(), + panel.grid.minor = ggplot2::element_blank(), + panel.grid.major.y = ggplot2::element_line(color = "#E5E2D9", linewidth = 0.4), + strip.text = ggplot2::element_text(color = "grey30", face = "bold", size = 10), + strip.background = ggplot2::element_rect(fill = "#EEEAE0", color = NA), + axis.title.y = ggplot2::element_blank(), + axis.text.x = ggplot2::element_text(color = "black", angle = 30, hjust = 1, size = 6), + axis.text.y = ggplot2::element_text(color = "black", size = 6), + legend.position = "bottom", + legend.text = ggplot2::element_text(color = "grey40", size = 10), + legend.title = ggplot2::element_text(color = "grey40", size = 10) + ) + + feat_plot +} From 7ece64468621c326cfc65b75973397381bf8c815 Mon Sep 17 00:00:00 2001 From: AbhirupaGhosh Date: Mon, 1 Jun 2026 23:47:29 +0000 Subject: [PATCH 06/17] Style code (GHA) --- R/plot_ml.R | 1163 ++++++++++++++++++++++++++------------------------- 1 file changed, 597 insertions(+), 566 deletions(-) diff --git a/R/plot_ml.R b/R/plot_ml.R index fab0759..825b488 100644 --- a/R/plot_ml.R +++ b/R/plot_ml.R @@ -326,51 +326,54 @@ plotFishers <- function( #' #' @examples #' plotDrugDist(metadata_path = "data/Campylobacter/") -plotDrugDist <- function(metadata_path = "."){ - -metadata <- arrow::read_parquet(file.path(metadata_path,"metadata.parquet")) - -##################### phenotype distribution (drugs) ######################### -drug_dist <- metadata |> - dplyr::distinct(genome.genome_id, - genome_drug.antibiotic, - drug_abbr, - genome_drug.resistant_phenotype) |> - dplyr::count(genome_drug.antibiotic, - drug_abbr, - genome_drug.resistant_phenotype) |> - dplyr::group_by(genome_drug.antibiotic, drug_abbr) |> - dplyr::mutate(total = sum(n)) |> - dplyr::ungroup() |> - dplyr::mutate( - label = paste0(genome_drug.antibiotic, " (", drug_abbr, ")"), - label = forcats::fct_reorder(label, total) - ) +plotDrugDist <- function(metadata_path = ".") { + metadata <- arrow::read_parquet(file.path(metadata_path, "metadata.parquet")) + + ##################### phenotype distribution (drugs) ######################### + drug_dist <- metadata |> + dplyr::distinct( + genome.genome_id, + genome_drug.antibiotic, + drug_abbr, + genome_drug.resistant_phenotype + ) |> + dplyr::count( + genome_drug.antibiotic, + drug_abbr, + genome_drug.resistant_phenotype + ) |> + dplyr::group_by(genome_drug.antibiotic, drug_abbr) |> + dplyr::mutate(total = sum(n)) |> + dplyr::ungroup() |> + dplyr::mutate( + label = paste0(genome_drug.antibiotic, " (", drug_abbr, ")"), + label = forcats::fct_reorder(label, total) + ) -p <- ggplot2::ggplot( - drug_dist, - ggplot2::aes( - x = label, - y = n, - fill = genome_drug.resistant_phenotype - ) -) + - ggplot2::geom_col(color = "black", width = 0.8) + - ggplot2::coord_flip() + - ggplot2::scale_fill_manual( - values = c( - "Resistant" = "#d4872a", - "Susceptible" = "#5b8db8" - ), - name = "Phenotype" - ) + - ggplot2::labs( - x = "Antibiotic", - y = "Number of unique genomes" + p <- ggplot2::ggplot( + drug_dist, + ggplot2::aes( + x = label, + y = n, + fill = genome_drug.resistant_phenotype + ) ) + - ggplot2::theme_classic(base_size = 14) + ggplot2::geom_col(color = "black", width = 0.8) + + ggplot2::coord_flip() + + ggplot2::scale_fill_manual( + values = c( + "Resistant" = "#d4872a", + "Susceptible" = "#5b8db8" + ), + name = "Phenotype" + ) + + ggplot2::labs( + x = "Antibiotic", + y = "Number of unique genomes" + ) + + ggplot2::theme_classic(base_size = 14) -p + p } #' Plot drug-level model performance @@ -387,123 +390,125 @@ p #' @examples #' plotDrugPerf(metadata_path = "data/Campylobacter/", performance_path = "data/Campylobacter/ML_performance/") plotDrugPerf <- function(metadata_path = ".", performance_path = ".") { - -metadata <- arrow::read_parquet(file.path(metadata_path,"metadata.parquet")) - -performance <- arrow::read_parquet(file.path(performance_path,"all_performance.parquet")) - -######################## drug performances ################################# -median_drug <- performance |> - dplyr::filter( - drug_label == "drug", - !shuffled # keep real models; remove if you want both - ) |> - dplyr::group_by(drug_or_class, feature_type, feature_subtype) |> - dplyr::summarise(median_mcc = median(mcc, na.rm = TRUE), .groups = "drop") |> - dplyr::left_join(plot_df, by = c("drug_or_class" = "drug_abbr")) |> - dplyr::mutate(drug_or_class = reorder(drug_or_class, total)) - -drug_p1 <- ggplot2::ggplot(median_drug, - ggplot2::aes(x = feature_type, - y = drug_or_class, - fill = median_mcc)) + - ggplot2::geom_tile(color = "grey90", width = 0.9) + - - ggplot2::scale_fill_gradientn( - colors = c( - "#C4B8A8", # low - "#FAFAF7", # around 0 - "#5F84C9", # medium/high (~0.7–0.9) - "#0F2A5A" # very dark for ~1 - ), - values = scales::rescale(c(-1, 0, 0.85, 1)), - name = "Best MCC" - ) + - - ggplot2::labs(x = "Feature type") + - ggplot2::theme_minimal(base_size = 12) + - ggplot2::theme( - axis.text = ggplot2::element_text(size = 10, colour = "black"), - axis.title = ggplot2::element_text(size = 12), - axis.title.y = ggplot2::element_blank(), - axis.text.x = ggplot2::element_text(angle = 45, hjust = 1), - legend.position = "bottom" - ) + - ggplot2::coord_fixed() + metadata <- arrow::read_parquet(file.path(metadata_path, "metadata.parquet")) -drug_p1 + performance <- arrow::read_parquet(file.path(performance_path, "all_performance.parquet")) -median_feature <- performance |> - dplyr::filter( - drug_label == "drug", - !shuffled - ) |> - dplyr::group_by(drug_or_class, feature_type) |> - dplyr::summarise(median_mcc = median(mcc, na.rm = TRUE), .groups = "drop") |> - dplyr::left_join(plot_df, by = c("drug_or_class" = "drug_abbr")) |> - dplyr::mutate(drug_or_class = reorder(drug_or_class, total)) - - -feat_pal <- c( - "args" = "#56B4E9", # sky blue - "cogs" = "#E69F00", # orange - "genes" = "#009E73", # bluish green - "domains" = "#F0E442", # yellow - "proteins" = "#CC79A7", # reddish purple - "struct" = "#D55E00" # vermillion -) - -rc_perf <- ggplot2::ggplot(median_feature |> - dplyr::distinct(drug_or_class, - feature_type, median_mcc), - ggplot2::aes(x = median_mcc, y = drug_or_class)) + - - ggridges::geom_density_ridges( - scale = 0.75, - rel_min_height = 0.01, - alpha = 0.4, - fill = "grey90", - colour = "grey70" - ) + - - ggplot2::geom_point( - position = position_jitter(height = 0.1), - size = 2, - alpha = 0.8, - aes(color = feature_type) - ) + - ggplot2::scale_color_manual(values = feat_pal, name = "Feature type") + - ggplot2::stat_summary( - fun = median, - geom = "point", - size = 2, - color = "black" + ######################## drug performances ################################# + median_drug <- performance |> + dplyr::filter( + drug_label == "drug", + !shuffled # keep real models; remove if you want both + ) |> + dplyr::group_by(drug_or_class, feature_type, feature_subtype) |> + dplyr::summarise(median_mcc = median(mcc, na.rm = TRUE), .groups = "drop") |> + dplyr::left_join(plot_df, by = c("drug_or_class" = "drug_abbr")) |> + dplyr::mutate(drug_or_class = reorder(drug_or_class, total)) + + drug_p1 <- ggplot2::ggplot( + median_drug, + ggplot2::aes( + x = feature_type, + y = drug_or_class, + fill = median_mcc + ) ) + - ggplot2::theme_minimal(base_size = 14) + - ggplot2::theme( - axis.text = ggplot2::element_text(size = 10, colour = "black"), - axis.title = ggplot2::element_text(size = 12), - axis.title.y = ggplot2::element_blank(), - axis.text.x = ggplot2::element_text(angle = 45, hjust = 1), - legend.position = "right", - panel.grid.minor = ggplot2::element_blank(), - axis.line = ggplot2::element_line(color = "black") - ) + ggplot2::geom_tile(color = "grey90", width = 0.9) + + ggplot2::scale_fill_gradientn( + colors = c( + "#C4B8A8", # low + "#FAFAF7", # around 0 + "#5F84C9", # medium/high (~0.7–0.9) + "#0F2A5A" # very dark for ~1 + ), + values = scales::rescale(c(-1, 0, 0.85, 1)), + name = "Best MCC" + ) + + ggplot2::labs(x = "Feature type") + + ggplot2::theme_minimal(base_size = 12) + + ggplot2::theme( + axis.text = ggplot2::element_text(size = 10, colour = "black"), + axis.title = ggplot2::element_text(size = 12), + axis.title.y = ggplot2::element_blank(), + axis.text.x = ggplot2::element_text(angle = 45, hjust = 1), + legend.position = "bottom" + ) + + ggplot2::coord_fixed() -rc_perf + drug_p1 -final_plot = drug_p1 + - rc_perf + - patchwork::plot_layout( - widths = c(2, 2), # adjust proportions - guides = "collect" - ) & - ggplot2::theme( - legend.position = "bottom" + median_feature <- performance |> + dplyr::filter( + drug_label == "drug", + !shuffled + ) |> + dplyr::group_by(drug_or_class, feature_type) |> + dplyr::summarise(median_mcc = median(mcc, na.rm = TRUE), .groups = "drop") |> + dplyr::left_join(plot_df, by = c("drug_or_class" = "drug_abbr")) |> + dplyr::mutate(drug_or_class = reorder(drug_or_class, total)) + + + feat_pal <- c( + "args" = "#56B4E9", # sky blue + "cogs" = "#E69F00", # orange + "genes" = "#009E73", # bluish green + "domains" = "#F0E442", # yellow + "proteins" = "#CC79A7", # reddish purple + "struct" = "#D55E00" # vermillion ) -final_plot + rc_perf <- ggplot2::ggplot( + median_feature |> + dplyr::distinct( + drug_or_class, + feature_type, median_mcc + ), + ggplot2::aes(x = median_mcc, y = drug_or_class) + ) + + ggridges::geom_density_ridges( + scale = 0.75, + rel_min_height = 0.01, + alpha = 0.4, + fill = "grey90", + colour = "grey70" + ) + + ggplot2::geom_point( + position = position_jitter(height = 0.1), + size = 2, + alpha = 0.8, + aes(color = feature_type) + ) + + ggplot2::scale_color_manual(values = feat_pal, name = "Feature type") + + ggplot2::stat_summary( + fun = median, + geom = "point", + size = 2, + color = "black" + ) + + ggplot2::theme_minimal(base_size = 14) + + ggplot2::theme( + axis.text = ggplot2::element_text(size = 10, colour = "black"), + axis.title = ggplot2::element_text(size = 12), + axis.title.y = ggplot2::element_blank(), + axis.text.x = ggplot2::element_text(angle = 45, hjust = 1), + legend.position = "right", + panel.grid.minor = ggplot2::element_blank(), + axis.line = ggplot2::element_line(color = "black") + ) + rc_perf + + final_plot <- drug_p1 + + rc_perf + + patchwork::plot_layout( + widths = c(2, 2), # adjust proportions + guides = "collect" + ) & + ggplot2::theme( + legend.position = "bottom" + ) + + final_plot } #' Plot cross-drug generalization heatmap @@ -520,132 +525,137 @@ final_plot #' @examples #' plotCrossDrug(cross_test_performance_path = "data/Campylobacter/cross_test_ML_performance", drug_performance_path = "data/Campylobacter/ML_performance/") plotCrossDrug <- function(cross_test_performance_path = ".", drug_performance_path = ".") { - - cross_drug <- arrow::read_parquet(file.path(cross_test_performance_path,"cross_drug_perf.parquet")) - performance <- arrow::read_parquet(file.path(drug_performance_path,"all_performance.parquet")) - -###################### CROSS DRUG Testing ############################# -heatmap_df <- cross_drug |> - # dplyr::filter(tested_on %in% (cross_drug |> dplyr::pull(drug_or_class))) |> - dplyr::group_by(drug_or_class, tested_on) |> - dplyr::summarise(median_mcc = median(mcc, na.rm = TRUE), .groups = "drop") - -same_drugs <- performance |> - dplyr::filter(drug_label == "drug", - drug_or_class %in% (cross_drug |> - dplyr::distinct(drug_or_class) |> - dplyr::pull())) |> - dplyr::group_by(drug_or_class) |> - dplyr::summarise(median_mcc = median(mcc, na.rm = TRUE), .groups = "drop") |> - dplyr::mutate(tested_on = drug_or_class) |> - dplyr::distinct(drug_or_class, tested_on, median_mcc) - -heatmap_df <- heatmap_df |> - dplyr::add_row(same_drugs) |> - dplyr::left_join(metadata |> - dplyr::distinct(drug_abbr, class_abbr), - by = c("drug_or_class" = "drug_abbr")) |> - dplyr::rename("drug_class" = "class_abbr") |> - dplyr::left_join(metadata |> - dplyr::distinct(drug_abbr, class_abbr), - by = c("tested_on" = "drug_abbr")) - -# Row annotation (already similar to what you did) -annotation_row <- heatmap_df |> - dplyr::distinct(drug_or_class, drug_class) |> - tibble::column_to_rownames("drug_or_class") - -# Column annotation -annotation_col <- heatmap_df |> - dplyr::distinct(tested_on, class_abbr) |> - tibble::column_to_rownames("tested_on") - -mat <- heatmap_df |> - dplyr::select(drug_or_class, tested_on, median_mcc) |> - tidyr::pivot_wider(names_from = tested_on, values_from = median_mcc) |> - tibble::column_to_rownames("drug_or_class") |> - as.matrix() - -row_order <- heatmap_df |> - dplyr::distinct(drug_or_class, drug_class) |> - dplyr::arrange(drug_class, drug_or_class) |> - dplyr::pull(drug_or_class) - -col_order <- heatmap_df |> - dplyr::distinct(tested_on, class_abbr) |> - dplyr::arrange(class_abbr, tested_on) |> - dplyr::pull(tested_on) - -# mat[is.na(mat)] <- 0 -mat <- mat[row_order, col_order] - -# Align annotations -annotation_row <- annotation_row[row_order, , drop = FALSE] -annotation_col <- annotation_col[col_order, , drop = FALSE] - - -# Collect all classes from both row and column -classes <- base::union( - annotation_row$drug_class, - annotation_col$class_abbr -) - -# Create ONE named color vector -class_colors <- stats::setNames( - scales::hue_pal() (length(classes)), - classes -) - -heat_colors <- colorRampPalette(RColorBrewer::brewer.pal(11, "RdBu"))(100) - -# ---- Convert annotations ---- -ha_row <- ComplexHeatmap::rowAnnotation( - drug_class = annotation_row$drug_class, - col = list(drug_class = class_colors), - show_annotation_name = FALSE, - show_legend = FALSE -) - -ha_col <- ComplexHeatmap::HeatmapAnnotation( - class_abbr = annotation_col$class_abbr, - col = list(class_abbr = class_colors), - show_annotation_name = FALSE, na_col = "grey3" -) - -# ---- Color function (instead of breaks + palette) ---- -col_fun <- circlize::colorRamp2( - seq(-max_val, max_val, length.out = length(heat_colors)), - heat_colors -) -# ---- Heatmap ---- -cross_drug_hm <- ComplexHeatmap::Heatmap( - mat, - name = "median_mcc", - col = col_fun, - cluster_rows = FALSE, - cluster_columns = FALSE, - row_order = row_order, - column_order = col_order, - left_annotation = ha_row, - top_annotation = ha_col, - show_row_names = TRUE, - show_column_names = TRUE, - column_title = "tested on", - row_title = "trained on", - column_title_side = "bottom", - row_title_side = "right", - row_names_gp = grid::gpar(fontsize = 14), - column_names_gp = grid::gpar(fontsize = 14), - column_names_rot = 0, - - # remove borders like pheatmap - rect_gp = grid::gpar(col = NA), - - # legends - show_heatmap_legend = TRUE -) - -cross_drug_hm + cross_drug <- arrow::read_parquet(file.path(cross_test_performance_path, "cross_drug_perf.parquet")) + performance <- arrow::read_parquet(file.path(drug_performance_path, "all_performance.parquet")) + + ###################### CROSS DRUG Testing ############################# + heatmap_df <- cross_drug |> + # dplyr::filter(tested_on %in% (cross_drug |> dplyr::pull(drug_or_class))) |> + dplyr::group_by(drug_or_class, tested_on) |> + dplyr::summarise(median_mcc = median(mcc, na.rm = TRUE), .groups = "drop") + + same_drugs <- performance |> + dplyr::filter( + drug_label == "drug", + drug_or_class %in% (cross_drug |> + dplyr::distinct(drug_or_class) |> + dplyr::pull()) + ) |> + dplyr::group_by(drug_or_class) |> + dplyr::summarise(median_mcc = median(mcc, na.rm = TRUE), .groups = "drop") |> + dplyr::mutate(tested_on = drug_or_class) |> + dplyr::distinct(drug_or_class, tested_on, median_mcc) + + heatmap_df <- heatmap_df |> + dplyr::add_row(same_drugs) |> + dplyr::left_join( + metadata |> + dplyr::distinct(drug_abbr, class_abbr), + by = c("drug_or_class" = "drug_abbr") + ) |> + dplyr::rename("drug_class" = "class_abbr") |> + dplyr::left_join( + metadata |> + dplyr::distinct(drug_abbr, class_abbr), + by = c("tested_on" = "drug_abbr") + ) + + # Row annotation (already similar to what you did) + annotation_row <- heatmap_df |> + dplyr::distinct(drug_or_class, drug_class) |> + tibble::column_to_rownames("drug_or_class") + + # Column annotation + annotation_col <- heatmap_df |> + dplyr::distinct(tested_on, class_abbr) |> + tibble::column_to_rownames("tested_on") + + mat <- heatmap_df |> + dplyr::select(drug_or_class, tested_on, median_mcc) |> + tidyr::pivot_wider(names_from = tested_on, values_from = median_mcc) |> + tibble::column_to_rownames("drug_or_class") |> + as.matrix() + + row_order <- heatmap_df |> + dplyr::distinct(drug_or_class, drug_class) |> + dplyr::arrange(drug_class, drug_or_class) |> + dplyr::pull(drug_or_class) + + col_order <- heatmap_df |> + dplyr::distinct(tested_on, class_abbr) |> + dplyr::arrange(class_abbr, tested_on) |> + dplyr::pull(tested_on) + + # mat[is.na(mat)] <- 0 + mat <- mat[row_order, col_order] + + # Align annotations + annotation_row <- annotation_row[row_order, , drop = FALSE] + annotation_col <- annotation_col[col_order, , drop = FALSE] + + + # Collect all classes from both row and column + classes <- base::union( + annotation_row$drug_class, + annotation_col$class_abbr + ) + + # Create ONE named color vector + class_colors <- stats::setNames( + scales::hue_pal()(length(classes)), + classes + ) + + heat_colors <- colorRampPalette(RColorBrewer::brewer.pal(11, "RdBu"))(100) + + # ---- Convert annotations ---- + ha_row <- ComplexHeatmap::rowAnnotation( + drug_class = annotation_row$drug_class, + col = list(drug_class = class_colors), + show_annotation_name = FALSE, + show_legend = FALSE + ) + + ha_col <- ComplexHeatmap::HeatmapAnnotation( + class_abbr = annotation_col$class_abbr, + col = list(class_abbr = class_colors), + show_annotation_name = FALSE, na_col = "grey3" + ) + + # ---- Color function (instead of breaks + palette) ---- + col_fun <- circlize::colorRamp2( + seq(-max_val, max_val, length.out = length(heat_colors)), + heat_colors + ) + # ---- Heatmap ---- + cross_drug_hm <- ComplexHeatmap::Heatmap( + mat, + name = "median_mcc", + col = col_fun, + cluster_rows = FALSE, + cluster_columns = FALSE, + row_order = row_order, + column_order = col_order, + left_annotation = ha_row, + top_annotation = ha_col, + show_row_names = TRUE, + show_column_names = TRUE, + column_title = "tested on", + row_title = "trained on", + column_title_side = "bottom", + row_title_side = "right", + row_names_gp = grid::gpar(fontsize = 14), + column_names_gp = grid::gpar(fontsize = 14), + column_names_rot = 0, + + # remove borders like pheatmap + rect_gp = grid::gpar(col = NA), + + # legends + show_heatmap_legend = TRUE + ) + + cross_drug_hm } #' Plot stratified model performance @@ -661,47 +671,63 @@ cross_drug_hm #' @export #' #' @examples -#' plotStratifiedPerf("year", stratified_performance_path = "data/Campylobacter/ML_year_performance", -#' stratified_cross_performance_path = "data/Campylobacter/cross_test_ML_year_performance") -plotStratifiedPerf <- function(year_or_country = "year", - stratified_performance_path = ".", +#' plotStratifiedPerf("year", +#' stratified_performance_path = "data/Campylobacter/ML_year_performance", +#' stratified_cross_performance_path = "data/Campylobacter/cross_test_ML_year_performance" +#' ) +plotStratifiedPerf <- function(year_or_country = "year", + stratified_performance_path = ".", stratified_cross_performance_path = ".") { - - perf <- arrow::read_parquet(file.path(stratified_performance_path, - paste0(year_or_country,"_perf.parquet"))) - - cross_test <- arrow::read_parquet(file.path(stratified_cross_performance_path, - paste0("cross_", - year_or_country, - "_perf.parquet"))) -if(year_or_country == "year") { -all <- perf |> - dplyr::rename("train_year" = "strat_value") |> - dplyr::mutate(test_year = train_year) |> - dplyr::select(drug_label, drug_or_class, - train_year, test_year, feature_type, feature_subtype, mcc) |> - dplyr::bind_rows(cross_test |> - dplyr::select(drug_label, drug_or_class, - train_year, test_year, feature_type, - feature_subtype, mcc)) |> - dplyr::mutate(category = dplyr::if_else( - train_year == test_year, "same year bin", "different year bin")) -} -else { -all <- perf |> - dplyr::rename("train_country" = "strat_value") |> - dplyr::mutate(test_country = train_country) |> - dplyr::select(drug_label, drug_or_class, - train_country, test_country, - feature_type, feature_subtype, mcc) |> - dplyr::bind_rows(cross_test |> - dplyr::select(drug_label, drug_or_class, - train_country, test_country, - feature_type, feature_subtype, mcc)) |> - dplyr::mutate(category = dplyr::if_else( - train_country == test_country, "same country", "different country")) -} - + perf <- arrow::read_parquet(file.path( + stratified_performance_path, + paste0(year_or_country, "_perf.parquet") + )) + + cross_test <- arrow::read_parquet(file.path( + stratified_cross_performance_path, + paste0( + "cross_", + year_or_country, + "_perf.parquet" + ) + )) + if (year_or_country == "year") { + all <- perf |> + dplyr::rename("train_year" = "strat_value") |> + dplyr::mutate(test_year = train_year) |> + dplyr::select( + drug_label, drug_or_class, + train_year, test_year, feature_type, feature_subtype, mcc + ) |> + dplyr::bind_rows(cross_test |> + dplyr::select( + drug_label, drug_or_class, + train_year, test_year, feature_type, + feature_subtype, mcc + )) |> + dplyr::mutate(category = dplyr::if_else( + train_year == test_year, "same year bin", "different year bin" + )) + } else { + all <- perf |> + dplyr::rename("train_country" = "strat_value") |> + dplyr::mutate(test_country = train_country) |> + dplyr::select( + drug_label, drug_or_class, + train_country, test_country, + feature_type, feature_subtype, mcc + ) |> + dplyr::bind_rows(cross_test |> + dplyr::select( + drug_label, drug_or_class, + train_country, test_country, + feature_type, feature_subtype, mcc + )) |> + dplyr::mutate(category = dplyr::if_else( + train_country == test_country, "same country", "different country" + )) + } + fill_vals <- if (year_or_country == "year") { c( "same year bin" = "#b3cde3", @@ -713,7 +739,7 @@ all <- perf |> "different country" = "#fbb4ae" ) } - + plot <- ggplot2::ggplot( all |> dplyr::filter(drug_label == "drug", !is.na(mcc)), @@ -750,7 +776,7 @@ all <- perf |> panel.grid.minor = ggplot2::element_blank(), plot.margin = margin(0, 0, 0, 0) ) -plot + plot } #' Plot multi-drug resistance (MDR) model performance @@ -766,205 +792,208 @@ plot #' @export #' #' @examples -#' plotMDR(MDR_performance_path = "data/Campylobacter/MDR_ML_performance", MDR_top_feature_path = "data/Campylobacter/MDR_ML_top_features", -#' MDR_pred_path = "data/Campylobacter/MDR_ML_pred") -plotMDR <- function(MDR_performance_path = ".", MDR_top_feature_path = ".", +#' plotMDR( +#' MDR_performance_path = "data/Campylobacter/MDR_ML_performance", MDR_top_feature_path = "data/Campylobacter/MDR_ML_top_features", +#' MDR_pred_path = "data/Campylobacter/MDR_ML_pred" +#' ) +plotMDR <- function(MDR_performance_path = ".", MDR_top_feature_path = ".", MDR_pred_path = ".") { - -MDR_perf <- arrow::read_parquet(file.path(MDR_performance_path, "MDR_perf.parquet")) - -# ---- Violin plot ---- -perf_plot <- ggplot2::ggplot(MDR_perf, - ggplot2::aes(x = feature_type, y = mcc)) + - - # violins (overall distribution per feature type) - ggplot2::geom_violin(fill = "grey85", color = NA, alpha = 0.8) + - - # points (colored by binary vs counts) - ggplot2::geom_jitter( - ggplot2::aes(color = feature_subtype), - width = 0.12, size = 2, alpha = 0.8 - ) + - - ggplot2::scale_color_manual(values = c( - "binary" = "#7B9CB5", - "counts" = "#CC8644" - )) + - - ggplot2::theme_minimal(base_size = 12) + - ggplot2::labs( - # title = "MDR model performances", - # subtitle = "Violin = distribution per feature type; points = binary vs counts", - x = "Feature type", - y = "MCC", - color = "Feature\nsubtype" - ) + - - ggplot2::theme( - legend.position = "right", - plot.title = ggplot2::element_text(face = "bold") - ) + - ggplot2::theme( - axis.title = ggplot2::element_text(colour = "black", size = 10), - axis.text.x = ggplot2::element_text(angle = 45, hjust = 1, size = 10, colour = "black"), - axis.text.y = ggplot2::element_text(size = 14, colour = "black"), - legend.title = ggplot2::element_text(size = 12), - legend.text = ggplot2::element_text(size = 10), - legend.position = "none", - title = ggplot2::element_text(face = "bold"), - - panel.background = ggplot2::element_blank(), - panel.grid.minor = ggplot2::element_blank(), - panel.grid.major.x = ggplot2::element_blank(), # remove vertical lines - panel.grid.major.y = ggplot2::element_line(color = "grey80"), # keep horizontal lines - - axis.line = ggplot2::element_line(color = "black") - ) + - ggplot2::scale_y_continuous(limits = c(0, 1)) - -perf_plot - -MDR_pred <- arrow::read_parquet(file.path(MDR_pred_path, "MDR_pred.parquet")) |> - dplyr::mutate( - diff_top2 = purrr::pmap_dbl(dplyr::across(dplyr::contains(".pred") & dplyr::where(is.numeric)), function(...) { - x <- c(...) - sx <- sort(x, decreasing = TRUE) - sx[1] - sx[2] - }) # Difference between prediction probabilities of top two classes - ) |> - dplyr::select(genome_id, resistant_classes, .pred_class, diff_top2, - feature_type, feature_subtype, seed) |> - dplyr::group_by(resistant_classes, .pred_class, feature_type) |> - dplyr::summarise(mean_margin = mean(diff_top2), n = n(), .groups = "drop") |> - dplyr::group_by(resistant_classes, feature_type) |> # normalize within true class - dplyr::mutate(sum = sum(n), prop = n / sum(n)) |> - dplyr::ungroup() - -MDR_pred_plot <- ggplot2::ggplot(MDR_pred, - ggplot2::aes(x = resistant_classes, - y = .pred_class)) + - ggplot2::geom_tile(ggplot2::aes(fill = prop)) + - ggplot2::geom_point(ggplot2::aes(size = mean_margin), color = "black") + - ggplot2::facet_wrap(~ feature_type) + - ggplot2::scale_fill_distiller( - palette = "RdBu", - direction = 1, # flip with -1 if needed - name = "Prediction proportion" + MDR_perf <- arrow::read_parquet(file.path(MDR_performance_path, "MDR_perf.parquet")) + + # ---- Violin plot ---- + perf_plot <- ggplot2::ggplot( + MDR_perf, + ggplot2::aes(x = feature_type, y = mcc) ) + - ggplot2::labs(x = "true class", y = "predicted class") + - ggplot2::scale_size(range = c(1, 6), name = "Mean margin") + - ggplot2::coord_equal() + - ggplot2::theme_minimal() + - ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45, hjust = 1)) + - ggplot2::theme( - axis.title = ggplot2::element_text(colour = "black", size = 10), - axis.text.x = ggplot2::element_text(angle = 45, hjust = 1, size = 10, colour = "black"), - axis.text.y = ggplot2::element_text(size = 10, colour = "black"), - legend.title = ggplot2::element_text(size = 12), - legend.text = ggplot2::element_text(size = 10), - legend.position = "right", - title = ggplot2::element_text(face = "bold") - ) -MDR_pred_plot - -# MDR_feat <- arrow::read_parquet(file.path( -# MDR_top_feature_path,"MDR_top_features.parquet")) |> -# pivot_longer(-c(Variable, feature_type, feature_subtype, seed), -# values_to = "Importance", -# names_to = "Resistant_classes") |> -# filter(!Importance == 0) -# -# MDR_feat_clean <- MDR_feat |> -# dplyr::filter(feature_type != "struct") |> -# dplyr::group_by(Resistant_classes, feature_type, feature_subtype, seed) |> -# dplyr::slice_max(Importance, n = top_n, with_ties = FALSE) |> -# dplyr::ungroup() |> -# dplyr::mutate(Variable = gsub( ".NCBIFAM", "", Variable)) |> -# dplyr::mutate(Variable = gsub("^X", "", Variable)) |> -# dplyr::mutate(Variable = dplyr::if_else( -# feature_type == "domains", gsub("_.*", "", Variable), Variable)) |> -# dplyr::mutate(Variable = dplyr::if_else( -# feature_type == "proteins", gsub("fig.", "fig|", Variable), Variable)) |> -# dplyr::left_join(cluster_feature, by = c("Variable" = "feature")) |> -# dplyr::mutate( -# cluster = dplyr::coalesce(cluster, Variable) -# ) -# -# cluster_df <- MDR_feat_clean |> -# dplyr::group_by(Resistant_classes, cluster) |> -# dplyr::summarise( -# Importance = median(Importance, na.rm = TRUE), -# .groups = "drop" -# ) -# -# top_clusters <- cluster_df |> -# group_by(Resistant_classes) |> -# group_modify(~{ -# -# df <- .x -# -# top_pos <- df |> -# arrange(desc(Importance)) |> -# slice_head(n = 10) -# -# top_neg <- df |> -# arrange(Importance) |> -# slice_head(n = 10) -# -# bind_rows(top_pos, top_neg) -# }) |> -# ungroup() -# -# top_clusters <- top_clusters |> -# dplyr::left_join(protein_names, by = c("cluster" = "proteinID")) |> -# dplyr::mutate( -# proteinName = dplyr::coalesce(proteinName, cluster), # fallback -# proteinName = stringr::str_trunc(proteinName, 50) -# ) |> -# dplyr::distinct(Resistant_classes, proteinName, Importance) |> -# # ✅ reorder AFTER naming -# dplyr::group_by(Resistant_classes) |> -# dplyr::mutate( -# proteinName = forcats::fct_reorder(proteinName, Importance) -# ) |> -# dplyr::ungroup() -# -# ggplot(top_clusters, -# aes(x = Importance, y = proteinName)) + -# -# # line (lollipop stem) -# geom_segment( -# aes(x = 0, xend = Importance, -# y = proteinName, yend = proteinName), -# color = "grey60" -# ) + -# -# # dot -# geom_point( -# aes(color = Importance > 0), -# size = 3 -# ) + -# -# facet_wrap(~ Resistant_classes, scales = "free_y") + -# -# scale_color_manual( -# values = c("TRUE" = "#5b8db8", # positive -# "FALSE" = "#d4872a"), # negative -# guide = "none" -# ) + -# -# theme_minimal(base_size = 13) + -# labs( -# x = "Median importance", -# y = "Cluster" -# ) + -# theme( -# panel.grid.minor = element_blank(), -# strip.text = element_text(face = "bold") -# ) + # violins (overall distribution per feature type) + ggplot2::geom_violin(fill = "grey85", color = NA, alpha = 0.8) + + + # points (colored by binary vs counts) + ggplot2::geom_jitter( + ggplot2::aes(color = feature_subtype), + width = 0.12, size = 2, alpha = 0.8 + ) + + ggplot2::scale_color_manual(values = c( + "binary" = "#7B9CB5", + "counts" = "#CC8644" + )) + + ggplot2::theme_minimal(base_size = 12) + + ggplot2::labs( + # title = "MDR model performances", + # subtitle = "Violin = distribution per feature type; points = binary vs counts", + x = "Feature type", + y = "MCC", + color = "Feature\nsubtype" + ) + + ggplot2::theme( + legend.position = "right", + plot.title = ggplot2::element_text(face = "bold") + ) + + ggplot2::theme( + axis.title = ggplot2::element_text(colour = "black", size = 10), + axis.text.x = ggplot2::element_text(angle = 45, hjust = 1, size = 10, colour = "black"), + axis.text.y = ggplot2::element_text(size = 14, colour = "black"), + legend.title = ggplot2::element_text(size = 12), + legend.text = ggplot2::element_text(size = 10), + legend.position = "none", + title = ggplot2::element_text(face = "bold"), + panel.background = ggplot2::element_blank(), + panel.grid.minor = ggplot2::element_blank(), + panel.grid.major.x = ggplot2::element_blank(), # remove vertical lines + panel.grid.major.y = ggplot2::element_line(color = "grey80"), # keep horizontal lines + + axis.line = ggplot2::element_line(color = "black") + ) + + ggplot2::scale_y_continuous(limits = c(0, 1)) + + perf_plot + MDR_pred <- arrow::read_parquet(file.path(MDR_pred_path, "MDR_pred.parquet")) |> + dplyr::mutate( + diff_top2 = purrr::pmap_dbl(dplyr::across(dplyr::contains(".pred") & dplyr::where(is.numeric)), function(...) { + x <- c(...) + sx <- sort(x, decreasing = TRUE) + sx[1] - sx[2] + }) # Difference between prediction probabilities of top two classes + ) |> + dplyr::select( + genome_id, resistant_classes, .pred_class, diff_top2, + feature_type, feature_subtype, seed + ) |> + dplyr::group_by(resistant_classes, .pred_class, feature_type) |> + dplyr::summarise(mean_margin = mean(diff_top2), n = n(), .groups = "drop") |> + dplyr::group_by(resistant_classes, feature_type) |> # normalize within true class + dplyr::mutate(sum = sum(n), prop = n / sum(n)) |> + dplyr::ungroup() + + MDR_pred_plot <- ggplot2::ggplot( + MDR_pred, + ggplot2::aes( + x = resistant_classes, + y = .pred_class + ) + ) + + ggplot2::geom_tile(ggplot2::aes(fill = prop)) + + ggplot2::geom_point(ggplot2::aes(size = mean_margin), color = "black") + + ggplot2::facet_wrap(~feature_type) + + ggplot2::scale_fill_distiller( + palette = "RdBu", + direction = 1, # flip with -1 if needed + name = "Prediction proportion" + ) + + ggplot2::labs(x = "true class", y = "predicted class") + + ggplot2::scale_size(range = c(1, 6), name = "Mean margin") + + ggplot2::coord_equal() + + ggplot2::theme_minimal() + + ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45, hjust = 1)) + + ggplot2::theme( + axis.title = ggplot2::element_text(colour = "black", size = 10), + axis.text.x = ggplot2::element_text(angle = 45, hjust = 1, size = 10, colour = "black"), + axis.text.y = ggplot2::element_text(size = 10, colour = "black"), + legend.title = ggplot2::element_text(size = 12), + legend.text = ggplot2::element_text(size = 10), + legend.position = "right", + title = ggplot2::element_text(face = "bold") + ) + MDR_pred_plot + + # MDR_feat <- arrow::read_parquet(file.path( + # MDR_top_feature_path,"MDR_top_features.parquet")) |> + # pivot_longer(-c(Variable, feature_type, feature_subtype, seed), + # values_to = "Importance", + # names_to = "Resistant_classes") |> + # filter(!Importance == 0) + # + # MDR_feat_clean <- MDR_feat |> + # dplyr::filter(feature_type != "struct") |> + # dplyr::group_by(Resistant_classes, feature_type, feature_subtype, seed) |> + # dplyr::slice_max(Importance, n = top_n, with_ties = FALSE) |> + # dplyr::ungroup() |> + # dplyr::mutate(Variable = gsub( ".NCBIFAM", "", Variable)) |> + # dplyr::mutate(Variable = gsub("^X", "", Variable)) |> + # dplyr::mutate(Variable = dplyr::if_else( + # feature_type == "domains", gsub("_.*", "", Variable), Variable)) |> + # dplyr::mutate(Variable = dplyr::if_else( + # feature_type == "proteins", gsub("fig.", "fig|", Variable), Variable)) |> + # dplyr::left_join(cluster_feature, by = c("Variable" = "feature")) |> + # dplyr::mutate( + # cluster = dplyr::coalesce(cluster, Variable) + # ) + # + # cluster_df <- MDR_feat_clean |> + # dplyr::group_by(Resistant_classes, cluster) |> + # dplyr::summarise( + # Importance = median(Importance, na.rm = TRUE), + # .groups = "drop" + # ) + # + # top_clusters <- cluster_df |> + # group_by(Resistant_classes) |> + # group_modify(~{ + # + # df <- .x + # + # top_pos <- df |> + # arrange(desc(Importance)) |> + # slice_head(n = 10) + # + # top_neg <- df |> + # arrange(Importance) |> + # slice_head(n = 10) + # + # bind_rows(top_pos, top_neg) + # }) |> + # ungroup() + # + # top_clusters <- top_clusters |> + # dplyr::left_join(protein_names, by = c("cluster" = "proteinID")) |> + # dplyr::mutate( + # proteinName = dplyr::coalesce(proteinName, cluster), # fallback + # proteinName = stringr::str_trunc(proteinName, 50) + # ) |> + # dplyr::distinct(Resistant_classes, proteinName, Importance) |> + # # ✅ reorder AFTER naming + # dplyr::group_by(Resistant_classes) |> + # dplyr::mutate( + # proteinName = forcats::fct_reorder(proteinName, Importance) + # ) |> + # dplyr::ungroup() + # + # ggplot(top_clusters, + # aes(x = Importance, y = proteinName)) + + # + # # line (lollipop stem) + # geom_segment( + # aes(x = 0, xend = Importance, + # y = proteinName, yend = proteinName), + # color = "grey60" + # ) + + # + # # dot + # geom_point( + # aes(color = Importance > 0), + # size = 3 + # ) + + # + # facet_wrap(~ Resistant_classes, scales = "free_y") + + # + # scale_color_manual( + # values = c("TRUE" = "#5b8db8", # positive + # "FALSE" = "#d4872a"), # negative + # guide = "none" + # ) + + # + # theme_minimal(base_size = 13) + + # labs( + # x = "Median importance", + # y = "Cluster" + # ) + + # theme( + # panel.grid.minor = element_blank(), + # strip.text = element_text(face = "bold") + # ) } #' Compare shuffled vs real model performance @@ -981,41 +1010,39 @@ MDR_pred_plot #' @examples #' plotShuffleVsReal(metadata_path = "data/Campylobacter/", performance_path = "data/Campylobacter/ML_performance") plotShuffleVsReal <- function(metadata_path = ".", performance_path = ".") { - - metadata <- arrow::read_parquet(file.path(metadata_path,"metadata.parquet")) - performance <- arrow::read_parquet(file.path(performance_path,"all_performance.parquet")) - - performance |> - dplyr::mutate( - shuffled_label = dplyr::if_else(shuffled, "shuffled", "real") - ) |> - ggplot2::ggplot(ggplot2::aes(x = feature_subtype, y = mcc, fill = shuffled_label)) + + metadata <- arrow::read_parquet(file.path(metadata_path, "metadata.parquet")) + performance <- arrow::read_parquet(file.path(performance_path, "all_performance.parquet")) + + performance |> + dplyr::mutate( + shuffled_label = dplyr::if_else(shuffled, "shuffled", "real") + ) |> + ggplot2::ggplot(ggplot2::aes(x = feature_subtype, y = mcc, fill = shuffled_label)) + ggplot2::geom_boxplot( - width = 0.55, outlier.size = 0.8, outlier.alpha = 0.4, - outlier.color = "grey50", linewidth = 0.4 - ) + + width = 0.55, outlier.size = 0.8, outlier.alpha = 0.4, + outlier.color = "grey50", linewidth = 0.4 + ) + ggplot2::geom_hline(yintercept = 0, linetype = "dashed", color = "grey60", linewidth = 0.4) + ggplot2::scale_fill_manual( - values = c("real" = "#7B9CB5", "shuffled" = "#C4B8A8"), - name = NULL - ) + + values = c("real" = "#7B9CB5", "shuffled" = "#C4B8A8"), + name = NULL + ) + ggplot2::scale_y_continuous(limits = c(-0.2, 1), breaks = seq(-0.2, 1, 0.2)) + - ggplot2::facet_wrap(~ feature_type, nrow = 1) + + ggplot2::facet_wrap(~feature_type, nrow = 1) + ggplot2::theme_minimal(base_size = 12) + - theme( - panel.grid.major.x = element_blank(), - panel.grid.minor = element_blank(), - panel.grid.major.y = element_line(color = "#E5E2D9", linewidth = 0.4), - strip.text = element_text(color = "grey30", face = "bold", size = 10), - # strip.background = element_rect(fill = "#EEEAE0", color = NA), - axis.text = element_text(color = "grey45"), - legend.position = "top", - legend.text = element_text(color = "grey40", size = 10) - ) + - labs( - x = NULL, y = "MCC" - ) - + theme( + panel.grid.major.x = element_blank(), + panel.grid.minor = element_blank(), + panel.grid.major.y = element_line(color = "#E5E2D9", linewidth = 0.4), + strip.text = element_text(color = "grey30", face = "bold", size = 10), + # strip.background = element_rect(fill = "#EEEAE0", color = NA), + axis.text = element_text(color = "grey45"), + legend.position = "top", + legend.text = element_text(color = "grey40", size = 10) + ) + + labs( + x = NULL, y = "MCC" + ) } #' Plot top contributing feature clusters @@ -1032,59 +1059,63 @@ plotShuffleVsReal <- function(metadata_path = ".", performance_path = ".") { #' @export #' #' @examples -#' plotTopClusters(top_feat_path = "data/Campylobacter/ML_top_features", cluster_feature_path = "data/Campylobacter/", -#' protein_names_path = "data/Campylobacter/", top_n = 10) -plotTopClusters <- function(top_feat_path = ".", cluster_feature_path = ".", +#' plotTopClusters( +#' top_feat_path = "data/Campylobacter/ML_top_features", cluster_feature_path = "data/Campylobacter/", +#' protein_names_path = "data/Campylobacter/", top_n = 10 +#' ) +plotTopClusters <- function(top_feat_path = ".", cluster_feature_path = ".", protein_names_path = ".", top_n = 10) { - ################### Top features ######################### - + top_feat <- arrow::read_parquet(file.path(top_feat_path, "all_top_features.parquet")) cluster_feature <- arrow::read_parquet(file.path(cluster_feature_path, "cluster_feature.parquet")) protein_names <- arrow::read_parquet(file.path(protein_names_path, "protein_names.parquet")) - + # which clusters appear in top n across feature types per drug # join top features with cluster mapping, filter out struct and shuffled - top_feat_clean <- top_feat |> + top_feat_clean <- top_feat |> dplyr::filter(!shuffled, feature_type != "struct", drug_label == "drug") |> dplyr::group_by(drug_or_class, feature_type, feature_subtype, seed) |> dplyr::slice_max(Importance, n = top_n, with_ties = FALSE) |> - dplyr::ungroup() |> - dplyr::mutate(Variable = gsub( ".NCBIFAM", "", Variable)) |> + dplyr::ungroup() |> + dplyr::mutate(Variable = gsub(".NCBIFAM", "", Variable)) |> dplyr::mutate(Variable = gsub("^X", "", Variable)) |> dplyr::mutate(Variable = dplyr::if_else( - feature_type == "domains", gsub("_.*", "", Variable), Variable)) |> + feature_type == "domains", gsub("_.*", "", Variable), Variable + )) |> dplyr::mutate(Variable = dplyr::if_else( - feature_type == "proteins", gsub("fig.", "fig|", Variable), Variable)) |> + feature_type == "proteins", gsub("fig.", "fig|", Variable), Variable + )) |> dplyr::left_join(cluster_feature, by = c("Variable" = "feature")) |> dplyr::mutate( - cluster = dplyr::coalesce(cluster, Variable), # fallback to Variable if no match + cluster = dplyr::coalesce(cluster, Variable), # fallback to Variable if no match Importance_signed = dplyr::if_else(Sign == "NEG", -Importance, Importance) ) - + shared_mat <- top_feat_clean |> dplyr::group_by(drug_or_class, feature_type, cluster) |> dplyr::summarise(abs_imp = median(Importance, na.rm = TRUE), .groups = "drop") |> - # convert to contribution within each feature_type dplyr::group_by(drug_or_class, feature_type) |> dplyr::mutate(contribution = abs_imp / sum(abs_imp, na.rm = TRUE)) |> - # pick top n contributors dplyr::slice_max(contribution, n = top_n, with_ties = FALSE) |> dplyr::ungroup() |> - dplyr::add_count(drug_or_class, cluster, name = "n_feat_types") |> dplyr::left_join(protein_names, by = c("cluster" = "proteinID")) |> dplyr::mutate( proteinName = stringr::str_trunc(proteinName, 50), proteinName = forcats::fct_reorder(proteinName, n_feat_types) ) - - feat_plot <- ggplot2::ggplot(shared_mat, - ggplot2::aes(x = feature_type, - y = proteinName, - fill = contribution)) + + + feat_plot <- ggplot2::ggplot( + shared_mat, + ggplot2::aes( + x = feature_type, + y = proteinName, + fill = contribution + ) + ) + ggplot2::geom_tile(color = "#FAFAF7", linewidth = 0.5, width = 0.9, height = 0.9) + # coord_fixed() + ggplot2::scale_fill_distiller( @@ -1093,7 +1124,7 @@ plotTopClusters <- function(top_feat_path = ".", cluster_feature_path = ".", name = "contribution", na.value = "#EEEAE0" ) + - ggplot2::facet_wrap(~ drug_or_class, scales = "free_y") + + ggplot2::facet_wrap(~drug_or_class, scales = "free_y") + ggplot2::theme_minimal(base_size = 12) + ggplot2::theme( panel.grid.major.x = ggplot2::element_blank(), @@ -1108,6 +1139,6 @@ plotTopClusters <- function(top_feat_path = ".", cluster_feature_path = ".", legend.text = ggplot2::element_text(color = "grey40", size = 10), legend.title = ggplot2::element_text(color = "grey40", size = 10) ) - + feat_plot } From 6a2787ecb9711ac03513f4a12123f241f355bb1c Mon Sep 17 00:00:00 2001 From: Alexander McKim Date: Tue, 30 Jun 2026 12:41:23 -0600 Subject: [PATCH 07/17] adding back pr 23 changes and more runnable examples --- DESCRIPTION | 10 ++ NAMESPACE | 13 +- R/arg_check_ml.R | 114 ++++++++-------- R/core_ml.R | 56 ++++---- R/data.R | 6 + R/generate_matrices_ml.R | 19 ++- R/globals.R | 26 ++++ R/ife_ml.R | 32 ++--- R/plot_ml.R | 217 ++++++++++++++++++++++--------- R/prep_ml.R | 14 +- R/run_ML.R | 6 +- R/run_ml_pipeline.R | 44 ++++--- man/calculateEvalMets.Rd | 4 +- man/demo_fit.Rd | 4 + man/demo_ml_tibble.Rd | 4 + man/dot-calculateMCC.Rd | 20 +++ man/dot-calculatenMCC.Rd | 9 +- man/dot-parquet2LOODrugMatrix.Rd | 20 +++ man/plotBaselineComparison.Rd | 30 +++-- man/plotCM.Rd | 26 ++++ man/plotCrossDrug.Rd | 47 +++++++ man/plotDefaultEval.Rd | 48 ------- man/plotDrugDist.Rd | 21 +++ man/plotDrugPerf.Rd | 38 ++++++ man/plotMDR.Rd | 34 +++++ man/plotPRC.Rd | 36 ++--- man/plotROC.Rd | 26 ++++ man/plotShuffleVsReal.Rd | 31 +++++ man/plotStratifiedPerf.Rd | 34 +++++ man/plotTopClusters.Rd | 37 ++++++ man/plotTopFeatsVI.Rd | 20 +-- man/runIFE.Rd | 4 +- vignettes/intro.Rmd | 2 +- 33 files changed, 754 insertions(+), 298 deletions(-) create mode 100644 man/dot-calculateMCC.Rd create mode 100644 man/dot-parquet2LOODrugMatrix.Rd create mode 100644 man/plotCM.Rd create mode 100644 man/plotCrossDrug.Rd delete mode 100644 man/plotDefaultEval.Rd create mode 100644 man/plotDrugDist.Rd create mode 100644 man/plotDrugPerf.Rd create mode 100644 man/plotMDR.Rd create mode 100644 man/plotROC.Rd create mode 100644 man/plotShuffleVsReal.Rd create mode 100644 man/plotStratifiedPerf.Rd create mode 100644 man/plotTopClusters.Rd diff --git a/DESCRIPTION b/DESCRIPTION index 7956cdb..62b06a0 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -38,24 +38,34 @@ Imports: DBI, dplyr, duckdb, + circlize, + forcats, ggplot2, ggrepel, + ggridges, glmnet, glue, + graphics, + grDevices, jsonlite, hardhat, methods, parsnip, purrr, + patchwork, + RColorBrewer, readr, recipes, rlang, rsample, + scales, sgof, + stats, stringr, tibble, tidyr, tune, + utils, vip, withr, workflows, diff --git a/NAMESPACE b/NAMESPACE index 0b9d60e..789117a 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -19,9 +19,17 @@ export(getConfusionMatrix) export(getNumFeat) export(loadMLInputTibble) export(plotBaselineComparison) -export(plotDefaultEval) +export(plotCM) +export(plotCrossDrug) +export(plotDrugDist) +export(plotDrugPerf) export(plotFishers) +export(plotMDR) export(plotPRC) +export(plotROC) +export(plotShuffleVsReal) +export(plotStratifiedPerf) +export(plotTopClusters) export(plotTopFeatsVI) export(predictML) export(removeTopFeats) @@ -87,6 +95,7 @@ importFrom(ggplot2,xlab) importFrom(ggplot2,ylim) importFrom(glmnet,glmnet) importFrom(glue,glue) +importFrom(grDevices,colorRampPalette) importFrom(graphics,barplot) importFrom(hardhat,tune) importFrom(jsonlite,fromJSON) @@ -124,7 +133,9 @@ 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(stringr,str_remove) importFrom(stringr,str_split) importFrom(tibble,add_column) diff --git a/R/arg_check_ml.R b/R/arg_check_ml.R index 2056323..a0b8af8 100644 --- a/R/arg_check_ml.R +++ b/R/arg_check_ml.R @@ -41,11 +41,11 @@ NULL "genome_drug.resistant_phenotype" %in% colnames(tibble), "resistant_classes" %in% colnames(tibble) )) { - stop(paste( - "The input tibble must have a target variable column named", - "either `genome_drug.resistant_phenotype` or `resistant_classes`, but", + stop( + "The input tibble must have a target variable column named ", + "either `genome_drug.resistant_phenotype` or `resistant_classes`, but ", "not both." - )) + ) } } } @@ -107,11 +107,11 @@ NULL is_train_val_test <- (train_prop > 0 && val_prop > 0 && test_prop > 0) if (!(is_pure_cv || is_train_val_test)) { - stop(paste0( + stop( "Unsupported split proportions. Use either:\n", " - pure CV: split = c(1, 0)\n", " - classical splits (train/val/test): c(train_prop, val_prop) with train > 0, val > 0, and test = 1 - train - val > 0." - )) + ) } # Sum must be <= 1 for all valid cases @@ -192,10 +192,10 @@ NULL } if (smallest_n_obs_rs < 1 || smallest_n_obs_rs %% 1 != 0) { - stop(paste( - "The `smallest_n_obs_rs` argument must be a whole, positive", + stop( + "The `smallest_n_obs_rs` argument must be a whole, positive ", "number greater than 0." - )) + ) } } @@ -207,10 +207,10 @@ NULL #' .checkArgUsePCA <- function(use_pca) { if (!is.logical(use_pca)) { - stop(paste( - "The `use_pca` argument can only take logical values:", + stop( + "The `use_pca` argument can only take logical values: ", "`TRUE` or `FALSE`." - )) + ) } } @@ -240,10 +240,10 @@ NULL #' .checkArgMultiClass <- function(multi_class) { if (!is.logical(multi_class)) { - stop(paste( - "The `multi_class` argument can only take logical values:", + stop( + "The `multi_class` argument can only take logical values: ", "`TRUE` or `FALSE`." - )) + ) } } @@ -281,10 +281,10 @@ NULL #' .checkArgModel <- function(model) { if (!(model %in% c("LR", "RF", "BT"))) { - stop(paste( - "The `model` argument must be one of:", + stop( + "The `model` argument must be one of: ", "'LR' (logistic regression), 'RF' (random forest), 'BT' (boosted tree)." - )) + ) } } @@ -434,10 +434,10 @@ NULL #' .checkArgSelectBestMetric <- function(select_best_metric) { if (!(select_best_metric %in% c("f_meas", "pr_auc", "mcc", "bal_accuracy"))) { - stop(paste( - "The `select_best_metric` argument must be one of: 'f_meas',", + stop( + "The `select_best_metric` argument must be one of: 'f_meas', ", "'pr_auc', 'mcc', 'bal_accuracy'." - )) + ) } } @@ -453,10 +453,10 @@ NULL .checkArgTibble(test_data_plus_predictions, ml = TRUE) if (!(".pred_class" %in% colnames(test_data_plus_predictions))) { - stop(paste( - "`test_data_plus_predictions` is missing a `.pred_class`", + stop( + "`test_data_plus_predictions` is missing a `.pred_class` ", "column. Try using the output of `predictML()`." - )) + ) } } @@ -489,10 +489,10 @@ NULL #' .checkArgPropVITopFeats <- function(prop_vi_top_feats) { if (!is.numeric(prop_vi_top_feats) || length(prop_vi_top_feats) != 2) { - stop(paste( - "The `prop_vi_top_feats` argument must be a vector of length", + stop( + "The `prop_vi_top_feats` argument must be a vector of length ", "2 with numeric values." - )) + ) } if (any(prop_vi_top_feats > 1) || any(prop_vi_top_feats < 0)) { @@ -515,10 +515,10 @@ NULL #' .checkArgShuffleLabels <- function(shuffle_labels) { if (!is.logical(shuffle_labels)) { - stop(paste( - "The `shuffle_labels` argument can only take logical values:", + stop( + "The `shuffle_labels` argument can only take logical values: ", "`TRUE` or `FALSE`." - )) + ) } } @@ -530,10 +530,10 @@ NULL #' .checkArgReturnTuneRes <- function(return_tune_res) { if (!is.logical(return_tune_res)) { - stop(paste( - "The `return_tune_res` argument can only take logical values:", + stop( + "The `return_tune_res` argument can only take logical values: ", "`TRUE` or `FALSE`." - )) + ) } } @@ -545,10 +545,10 @@ NULL #' .checkArgReturnFit <- function(return_fit) { if (!is.logical(return_fit)) { - stop(paste( - "The `return_fit` argument can only take logical values:", + stop( + "The `return_fit` argument can only take logical values: ", "`TRUE` or `FALSE`." - )) + ) } } @@ -561,10 +561,10 @@ NULL #' .checkArgReturnPred <- function(return_pred) { if (!is.logical(return_pred)) { - stop(paste( - "The `return_pred` argument can only take logical values:", + stop( + "The `return_pred` argument can only take logical values: ", "`TRUE` or `FALSE`." - )) + ) } } @@ -576,10 +576,10 @@ NULL #' .checkArgVerbose <- function(verbose) { if (!is.logical(verbose)) { - stop(paste( - "The `verbose` argument can only take logical values:", + stop( + "The `verbose` argument can only take logical values: ", "`TRUE` or `FALSE`." - )) + ) } } @@ -613,10 +613,10 @@ NULL #' .checkArgByVI <- function(by_vi) { if (!is.logical(by_vi)) { - stop(paste( - "The `by_vi` argument can only take logical values:", + stop( + "The `by_vi` argument can only take logical values: ", "`TRUE` or `FALSE`." - )) + ) } } @@ -629,10 +629,10 @@ NULL #' .checkArgByNum <- function(by_num) { if (!is.logical(by_num)) { - stop(paste( - "The `by_num` argument can only take logical values:", + stop( + "The `by_num` argument can only take logical values: ", "`TRUE` or `FALSE`." - )) + ) } } @@ -657,10 +657,10 @@ NULL } if (any(diff(percent_removal_vec) <= 0)) { - stop(paste( - "The elements of `percent_removal_vec` must be in ascending", + stop( + "The elements of `percent_removal_vec` must be in ascending ", "order with no repeats." - )) + ) } } @@ -673,10 +673,10 @@ NULL #' .checkArgReturnFeats <- function(return_feats) { if (!is.logical(return_feats)) { - stop(paste( - "The `return_feats` argument can only take logical values:", + stop( + "The `return_feats` argument can only take logical values: ", "`TRUE` or `FALSE`." - )) + ) } } @@ -713,10 +713,10 @@ NULL if (!(y_default_eval %in% c("avg_f1_score", "avg_log2_apop", "avg_bal_acc", "avg_mcc", "avg_nmcc")) ) { - stop(paste( - "`y_default_eval` must be one of:", + stop( + "`y_default_eval` must be one of: ", "'avg_f1_score', 'avg_log2_apop', 'avg_bal_acc', 'avg_mcc', 'avg_nmcc'." - )) + ) } } diff --git a/R/core_ml.R b/R/core_ml.R index ec29d41..7e398f2 100644 --- a/R/core_ml.R +++ b/R/core_ml.R @@ -531,10 +531,10 @@ getConfusionMatrix <- function(test_data_plus_predictions) { if (!("genome_drug.resistant_phenotype" %in% colnames(test_data_plus_predictions))) { - stop(paste( - "`test_data_plus_predictions` does not have a column for", + stop( + "`test_data_plus_predictions` does not have a column for ", "`genome_drug.resistant_phenotype`." - )) + ) } f1 <- test_data_plus_predictions |> @@ -561,10 +561,10 @@ getConfusionMatrix <- function(test_data_plus_predictions) { if (!("genome_drug.resistant_phenotype" %in% colnames(test_data_plus_predictions))) { - stop(paste( - "`test_data_plus_predictions` does not have a column for", + stop( + "`test_data_plus_predictions` does not have a column for ", "`genome_drug.resistant_phenotype`." - )) + ) } auprc <- test_data_plus_predictions |> @@ -590,10 +590,10 @@ getConfusionMatrix <- function(test_data_plus_predictions) { if (!("genome_drug.resistant_phenotype" %in% colnames(test_data_plus_predictions))) { - stop(paste( - "`test_data_plus_predictions` does not have a column for", + stop( + "`test_data_plus_predictions` does not have a column for ", "`genome_drug.resistant_phenotype`." - )) + ) } auprc <- .calculateAUPRC(test_data_plus_predictions) @@ -603,15 +603,15 @@ getConfusionMatrix <- function(test_data_plus_predictions) { nrow(test_data_plus_predictions) if (prior > 0.3 && prior < 0.7) { - warning(paste( - "Classes are roughly balanced.", + warning( + "Classes are roughly balanced. ", "Calculation of log2(AUPRC/prior) may be inappropriate." - )) + ) } else if (prior >= 0.7) { - warning(paste( - "Classes are imbalanced for this model.", + warning( + "Classes are imbalanced for this model. ", "The use of the log2(AUPRC/prior) metric may be more informative in this imbalanced model." - )) + ) } log2_apop <- log2(auprc / prior) |> round(2) @@ -631,10 +631,10 @@ getConfusionMatrix <- function(test_data_plus_predictions) { if (!("genome_drug.resistant_phenotype" %in% colnames(test_data_plus_predictions))) { - stop(paste( - "`test_data_plus_predictions` does not have a column for", + stop( + "`test_data_plus_predictions` does not have a column for ", "`genome_drug.resistant_phenotype`." - )) + ) } bal_acc <- test_data_plus_predictions |> @@ -660,10 +660,10 @@ getConfusionMatrix <- function(test_data_plus_predictions) { if (!("genome_drug.resistant_phenotype" %in% colnames(test_data_plus_predictions))) { - stop(paste( - "`test_data_plus_predictions` does not have a column for", + stop( + "`test_data_plus_predictions` does not have a column for ", "`genome_drug.resistant_phenotype`." - )) + ) } sens <- test_data_plus_predictions |> @@ -690,10 +690,10 @@ getConfusionMatrix <- function(test_data_plus_predictions) { if (!("genome_drug.resistant_phenotype" %in% colnames(test_data_plus_predictions))) { - stop(paste( - "`test_data_plus_predictions` does not have a column for", + stop( + "`test_data_plus_predictions` does not have a column for ", "`genome_drug.resistant_phenotype`." - )) + ) } spec <- test_data_plus_predictions |> @@ -799,7 +799,7 @@ extractTopFeats <- function( dplyr::arrange(dplyr::desc(Importance)) if (!is.na(n_top_feats)) { - top_feats_and_VIs <- feats_arranged |> dplyr::slice(1:n_top_feats) + top_feats_and_VIs <- feats_arranged |> dplyr::slice(seq_len(n_top_feats)) } else if (any(!is.na(prop_vi_top_feats))) { cum_vi_lower <- prop_vi_top_feats[1] * sum(feats_arranged$Importance) cum_vi_upper <- prop_vi_top_feats[2] * sum(feats_arranged$Importance) @@ -818,10 +818,10 @@ extractTopFeats <- function( # Take a different approach if using multi-class (the previous code would give # a less meaningful result). if (is(fit$fit$actions$model$spec, "multinom_reg")) { - warning(paste( - "Extracting top features from a multi-class model.", + warning( + "Extracting top features from a multi-class model. ", "The `prop_vi_top_feats` and `n_top_feats` arguments do not apply." - )) + ) fit_penalty <- .getFitHps(fit)["penalty"] |> as.numeric() glmnet_fit <- parsnip::extract_fit_engine(fit) diff --git a/R/data.R b/R/data.R index 0f06f74..07cbed9 100644 --- a/R/data.R +++ b/R/data.R @@ -7,6 +7,9 @@ #' @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 @@ -15,4 +18,7 @@ #' #' @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/generate_matrices_ml.R b/R/generate_matrices_ml.R index 55e04c6..1f4de04 100644 --- a/R/generate_matrices_ml.R +++ b/R/generate_matrices_ml.R @@ -944,10 +944,10 @@ skipImbalancedMatrix <- function(genome_ids, } }), prefix = purrr::pmap_chr(list(parts, idx_sparse), \(x, i) { - if (is.na(i) || i < 3) { + if (is.na(i) || i < 5) { NA_character_ } else { - paste(x[1:(i - 4)], collapse = "_") + paste(x[seq_len(i - 4)], collapse = "_") } }), drug = purrr::pmap_chr(list(parts, idx_sparse), \(x, i) { @@ -986,7 +986,7 @@ skipImbalancedMatrix <- function(genome_ids, ## -------------------------------------------------------------- ## Leave-one-drug-out with genome blocking ## -------------------------------------------------------------- - purrr::walk(names(grouped), function(leave_one_out) { + new_files <- purrr::map(names(grouped), function(leave_one_out) { ## Read held-out drug to get its genome IDs heldout_tbl <- arrow::read_parquet(grouped[[leave_one_out]]$file) @@ -1036,7 +1036,6 @@ skipImbalancedMatrix <- function(genome_ids, )) arrow::write_parquet(combined, out_file) - created <<- c(created, out_file) log( "debug", @@ -1045,7 +1044,13 @@ skipImbalancedMatrix <- function(genome_ids, " | removed ", length(heldout_genomes), " genomes" ) ) - }) + + out_file + }) |> + purrr::compact() |> + unlist() + + created <- c(created, new_files) } log("info", "All LOO-drug matrices generated with genome-level blocking") @@ -1105,10 +1110,10 @@ skipImbalancedMatrix <- function(genome_ids, } }), prefix = purrr::pmap_chr(list(parts, idx_sparse), \(x, i) { - if (is.na(i) || i < 3) { + if (is.na(i) || i < 5) { NA_character_ } else { - paste(x[1:(i - 4)], collapse = "_") + paste(x[seq_len(i - 4)], collapse = "_") } }), drug = purrr::pmap_chr(list(parts, idx_sparse), \(x, i) { diff --git a/R/globals.R b/R/globals.R index ef6b85e..3c440fe 100644 --- a/R/globals.R +++ b/R/globals.R @@ -4,6 +4,9 @@ # Package-level imports for packages used with :: notation #' @importFrom jsonlite fromJSON write_json #' @importFrom glmnet glmnet +#' @importFrom grDevices colorRampPalette +#' @importFrom stats median +#' @importFrom stats reorder #' @keywords internal "_PACKAGE" @@ -22,10 +25,19 @@ utils::globalVariables(c( "cum_imp", # Data columns + "abs_imp", "adj_p_value", "antibiotic", "bal_acc", + "category", + "class_abbr", + "cluster", + "contribution", + "diff_top2", "drug", + "drug_abbr", + "drug_class", + "drug_label", "drug_or_class", "feature", "feature_id", @@ -36,6 +48,7 @@ utils::globalVariables(c( "fit_mtry", "fit_trees", "gene", + "genome_drug.antibiotic", "genome_drug.genome_id", "genome_drug.resistant_phenotype", "genome_id", @@ -43,9 +56,13 @@ utils::globalVariables(c( "i_strat", "idx_sparse", "idx_strat", + "label", "model", "neg_log10_adj_p", "mcc", + "mean_margin", + "median_mcc", + "n_feat_types", "nmcc", "num_obs", "seed", @@ -56,6 +73,7 @@ utils::globalVariables(c( "pair_id", "parts", "phenotype", + "proteinName", "precision", "prefix", "prefix_key", @@ -66,14 +84,22 @@ utils::globalVariables(c( "res_prop", "resistant_classes", "run_time_sec", + "shuffled", + "shuffled_label", "sig_after_bh", "significance", "strat_value", "strat_value_test", "stratification", + "test_country", "test_drug", "test_file", + "test_year", + "tested_on", + "total", + "train_country", "train_prop", + "train_year", "value", # Base R functions that need explicit import diff --git a/R/ife_ml.R b/R/ife_ml.R index e6e38c0..d4fae1e 100644 --- a/R/ife_ml.R +++ b/R/ife_ml.R @@ -34,10 +34,10 @@ removeTopFeats <- function(ml_input_tibble, top_feat_tibble) { feats_to_remove <- top_feat_tibble |> dplyr::pull(Variable) if (length(feats_to_remove) == 0) { - stop(paste( - "No more feats to remove; try adjusting `percent_removal_vec`", + stop( + "No more feats to remove; try adjusting `percent_removal_vec` ", "in `runIFE()` to include a wider range between percentages." - )) + ) } ml_input_tibble_top_feats_removed <- ml_input_tibble |> @@ -92,10 +92,10 @@ runIFE <- function( .checkArgVerbose(verbose) if (sum(c(by_num, by_vi)) != 1) { - stop(paste( - "Set either `by_num` or `by_vi` to `TRUE` and the other to", + stop( + "Set either `by_num` or `by_vi` to `TRUE` and the other to ", "`FALSE` (not both `TRUE` or both `FALSE`)." - )) + ) } # Create empty objects to store results. @@ -119,7 +119,7 @@ runIFE <- function( message("Training with all features") } - for (i in 1:(length(percent_removal_vec) + 1)) { + for (i in seq_len(length(percent_removal_vec) + 1)) { if (i == 1) { n_feats_removed <- 0 } @@ -227,19 +227,19 @@ runIFE <- function( # Removal of top features if (by_vi) { if (verbose) { - message(paste("Cumulatively removed top ", - percent_removal_vec[i], "% of total variable importance", - sep = "" - )) + message( + "Cumulatively removed top ", + percent_removal_vec[i], "% of total variable importance" + ) } ml_input_tibble <- ml_input_tibble |> removeTopFeats(feats_to_remove_tibble) } else if (by_num) { if (verbose) { - message(paste("Cumulatively removed top ", - percent_removal_vec[i], "% of all features", - sep = "" - )) + message( + "Cumulatively removed top ", + percent_removal_vec[i], "% of all features" + ) } feats_to_remove_tibble <- ml_res$top_feat_tibble @@ -264,7 +264,7 @@ runIFE <- function( # Adjust `percent_removal_vec` to appropriate length (if IFE stopped running # early due to encountering an error). - percent_removal_vec <- percent_removal_vec[1:(length(num_obs_vec) - 1)] + percent_removal_vec <- percent_removal_vec[seq_len(length(num_obs_vec) - 1)] ife_performance_tibble <- tibble::tibble( percent_removed = c(0, percent_removal_vec), diff --git a/R/plot_ml.R b/R/plot_ml.R index bc1b9b9..7ab9f8c 100644 --- a/R/plot_ml.R +++ b/R/plot_ml.R @@ -16,6 +16,7 @@ #' @importFrom ggplot2 theme #' @importFrom ggplot2 xlab #' @importFrom ggplot2 ylim +#' @importFrom graphics barplot #' @importFrom tune extract_fit_parsnip #' @importFrom vip vip #' @importFrom yardstick pr_curve @@ -24,8 +25,9 @@ NULL #' Plot a Precision-Recall Curve #' #' Generates a precision-recall curve (PRC) for AMR phenotype prediction results. -#' @param test_data_plus_predictions_file A file containing test data with added -#' prediction columns, typically the output of `runMLmodels(return_pred=TRUE)`. +#' @param test_data_plus_predictions A tibble of test data with added prediction +#' columns (e.g. the output of `predictML()` or `runMLmodels(return_pred=TRUE)`), +#' or a path to a TSV file containing the same. #' #' @return A `ggplot2` object showing the precision-recall curve. #' @@ -34,14 +36,16 @@ NULL #' visualizes it using `ggplot2`. #' #' @examples -#' \dontrun{ -#' test_data_plus_predictions_file <- "results / ML_pred / Sfl_drug_AMP_domains_binary_prediction.tsv" -#' plotPRC(test_data_plus_predictions_file) -#' } +#' data(demo_fit) +#' data(demo_ml_tibble) +#' preds <- predictML(demo_fit, demo_ml_tibble) +#' plotPRC(preds) #' -#' @export -plotPRC <- function(test_data_plus_predictions_file) { - test_data_plus_predictions <- readr::read_tsv(test_data_plus_predictions_file) +#' @export +plotPRC <- function(test_data_plus_predictions) { + if (is.character(test_data_plus_predictions)) { + test_data_plus_predictions <- readr::read_tsv(test_data_plus_predictions) + } .checkArgTestDataPlusPredictions(test_data_plus_predictions) test_data_plus_predictions <- test_data_plus_predictions |> dplyr::mutate( @@ -67,14 +71,23 @@ plotPRC <- function(test_data_plus_predictions_file) { #' #' Generates a ROC curve for AMR phenotype prediction results. #' -#' @param test_data_plus_predictions_file A file with test data and prediction -#' columns (output of `runMLmodels(return_pred=TRUE)`). +#' @param test_data_plus_predictions A tibble of test data with prediction +#' columns (e.g. the output of `predictML()` or `runMLmodels(return_pred=TRUE)`), +#' or a path to a TSV file containing the same. #' #' @return A ROC curve plotted using `ggplot2::autoplot()`. #' +#' @examples +#' data(demo_fit) +#' data(demo_ml_tibble) +#' preds <- predictML(demo_fit, demo_ml_tibble) +#' plotROC(preds) +#' #' @export -plotROC <- function(test_data_plus_predictions_file) { - test_data_plus_predictions <- readr::read_tsv(test_data_plus_predictions_file) +plotROC <- function(test_data_plus_predictions) { + if (is.character(test_data_plus_predictions)) { + test_data_plus_predictions <- readr::read_tsv(test_data_plus_predictions) + } .checkArgTestDataPlusPredictions(test_data_plus_predictions) test_data_plus_predictions <- test_data_plus_predictions |> dplyr::mutate( @@ -98,14 +111,23 @@ plotROC <- function(test_data_plus_predictions_file) { #' #' Produces a heatmap visualization of the confusion matrix for AMR predictions. #' -#' @param test_data_plus_predictions_file A file containing true and predicted -#' phenotype labels. +#' @param test_data_plus_predictions A tibble containing true and predicted +#' phenotype labels (e.g. the output of `predictML()`), or a path to a TSV file +#' containing the same. #' #' @return A heatmap (`ggplot2` object) showing the confusion matrix. #' +#' @examples +#' data(demo_fit) +#' data(demo_ml_tibble) +#' preds <- predictML(demo_fit, demo_ml_tibble) +#' plotCM(preds) +#' #' @export -plotCM <- function(test_data_plus_predictions_file) { - test_data_plus_predictions <- readr::read_tsv(test_data_plus_predictions_file) +plotCM <- function(test_data_plus_predictions) { + if (is.character(test_data_plus_predictions)) { + test_data_plus_predictions <- readr::read_tsv(test_data_plus_predictions) + } .checkArgTestDataPlusPredictions(test_data_plus_predictions) test_data_plus_predictions <- test_data_plus_predictions |> dplyr::mutate( @@ -131,21 +153,23 @@ plotCM <- function(test_data_plus_predictions_file) { #' Creates a bar plot showing the most important features affecting #' AMR phenotype predictions. #' -#' @param topfeat_file A file containing feature importance scores -#' (output of `runMLmodels()`). +#' @param topfeat A tibble of feature importance scores (e.g. the output of +#' `extractTopFeats()` or `runMLmodels()`), with `Variable`, `Importance`, and +#' `Sign` columns, or a path to a TSV file containing the same. #' @param n_top_feats Number of top features to display (default: 10). #' #' @return A bar plot of variable importance (`ggplot2` object). #' #' @examples -#' \dontrun{ -#' topfeat_file <- "results / ML_top_features / Sfl_drug_AMP_domains_binary_top_features.tsv" -#' plotTopFeatsVI(topfeat_file) -#' } +#' data(demo_fit) +#' top_feats <- extractTopFeats(demo_fit, n_top_feats = 10) +#' plotTopFeatsVI(top_feats, n_top_feats = 10) #' #' @export -plotTopFeatsVI <- function(topfeat_file, n_top_feats = 10) { - topfeat <- readr::read_tsv(topfeat_file) +plotTopFeatsVI <- function(topfeat, n_top_feats = 10) { + if (is.character(topfeat)) { + topfeat <- readr::read_tsv(topfeat) + } .checkArgNTopFeats(n_top_feats) vip <- topfeat |> @@ -185,19 +209,18 @@ plotTopFeatsVI <- function(topfeat_file, n_top_feats = 10) { #' #' @return A base R barplot comparing balanced accuracy across models. #' -#' @param non_shuffled_label_results Output of `runMLPipeline()` -#' (`shuffle_labels = FALSE`) -#' @param shuffled_label_results Output of `runMLPipeline()` -#' (`shuffle_labels = TRUE`) -#' @return A bar plot with balanced accuracy comparisons per antibiotic #' @examples -#' non_shuffled <- tibble::tibble( -#' antibiotic = c("AMP", "CIP", "CRO"), -#' bal_acc = c(0.88, 0.81, 0.92) +#' non_shuffled <- list( +#' performance_tibble = tibble::tibble( +#' antibiotic = c("AMP", "CIP", "CRO"), +#' bal_acc = c(0.88, 0.81, 0.92) +#' ) #' ) -#' shuffled <- tibble::tibble( -#' antibiotic = c("AMP", "CIP", "CRO"), -#' bal_acc = c(0.52, 0.49, 0.55) +#' shuffled <- list( +#' performance_tibble = tibble::tibble( +#' antibiotic = c("AMP", "CIP", "CRO"), +#' bal_acc = c(0.52, 0.49, 0.55) +#' ) #' ) #' plotBaselineComparison(non_shuffled, shuffled) #' @export @@ -352,14 +375,14 @@ plotFishers <- function( #' @export #' #' @examples -#' plotDrugDist(metadata_path = "data/Campylobacter/") +#' plotDrugDist(metadata_path = system.file("extdata", package = "amRml")) plotDrugDist <- function(metadata_path = ".") { metadata <- arrow::read_parquet(file.path(metadata_path, "metadata.parquet")) ##################### phenotype distribution (drugs) ######################### drug_dist <- metadata |> dplyr::distinct( - genome.genome_id, + genome_drug.genome_id, genome_drug.antibiotic, drug_abbr, genome_drug.resistant_phenotype @@ -408,18 +431,42 @@ plotDrugDist <- function(metadata_path = ".") { #' Generates heatmaps and ridge plots summarizing model performance (MCC) #' across drugs and feature types. #' -#' @param metadata_path Character. Path to `metadata.parquet`. -#' @param performance_path Character. Path to `all_performance.parquet`. +#' @param metadata_path Character. Path to a directory containing `metadata.parquet`. +#' @param performance_path A performance tibble (with `drug_label`, `shuffled`, +#' `drug_or_class`, `feature_type`, `feature_subtype`, and `mcc` columns), or a +#' directory path containing `all_performance.parquet`. #' #' @return A patchwork ggplot object combining multiple panels. #' @export #' #' @examples -#' plotDrugPerf(metadata_path = "data/Campylobacter/", performance_path = "data/Campylobacter/ML_performance/") +#' # Several models per drug x feature type so the ridge densities can be drawn. +#' performance <- tidyr::expand_grid( +#' drug_or_class = c("AMP", "CIP", "NAL"), +#' feature_type = c("genes", "proteins"), +#' feature_subtype = c("binary", "counts"), +#' replicate = 1:8 +#' ) +#' performance$drug_label <- "drug" +#' performance$shuffled <- FALSE +#' performance$mcc <- 0.6 + 0.3 * sin(seq_len(nrow(performance))) +#' plotDrugPerf( +#' metadata_path = system.file("extdata", package = "amRml"), +#' performance_path = performance +#' ) plotDrugPerf <- function(metadata_path = ".", performance_path = ".") { metadata <- arrow::read_parquet(file.path(metadata_path, "metadata.parquet")) - performance <- arrow::read_parquet(file.path(performance_path, "all_performance.parquet")) + if (!is.data.frame(performance_path)) { + performance_path <- arrow::read_parquet( + file.path(performance_path, "all_performance.parquet") + ) + } + performance <- performance_path + + plot_df <- metadata |> + dplyr::distinct(genome_drug.genome_id, genome_drug.antibiotic, drug_abbr) |> + dplyr::count(genome_drug.antibiotic, drug_abbr, name = "total") ######################## drug performances ################################# median_drug <- performance |> @@ -525,12 +572,11 @@ plotDrugPerf <- function(metadata_path = ".", performance_path = ".") { rc_perf - final_plot <- drug_p1 + - rc_perf + - patchwork::plot_layout( - widths = c(2, 2), # adjust proportions - guides = "collect" - ) & + final_plot <- patchwork::wrap_plots( + drug_p1, rc_perf, + widths = c(2, 2), # adjust proportions + guides = "collect" + ) & ggplot2::theme( legend.position = "bottom" ) @@ -543,17 +589,51 @@ plotDrugPerf <- function(metadata_path = ".", performance_path = ".") { #' Creates a heatmap showing cross-drug model performance (MCC), where models #' trained on one drug are evaluated on another. #' -#' @param cross_test_performance_path Character. Path to `cross_drug_perf.parquet`. -#' @param drug_performance_path Character. Path to `all_performance.parquet`. +#' @param cross_test_performance_path A cross-drug performance tibble (with +#' `drug_or_class`, `tested_on`, and `mcc` columns), or a directory path +#' containing `cross_drug_perf.parquet`. +#' @param drug_performance_path A performance tibble (with `drug_label`, +#' `drug_or_class`, and `mcc` columns), or a directory path containing +#' `all_performance.parquet`. +#' @param metadata_path Character. Path to a directory containing `metadata.parquet`. #' #' @return A ComplexHeatmap object. #' @export #' #' @examples -#' plotCrossDrug(cross_test_performance_path = "data/Campylobacter/cross_test_ML_performance", drug_performance_path = "data/Campylobacter/ML_performance/") -plotCrossDrug <- function(cross_test_performance_path = ".", drug_performance_path = ".") { - cross_drug <- arrow::read_parquet(file.path(cross_test_performance_path, "cross_drug_perf.parquet")) - performance <- arrow::read_parquet(file.path(drug_performance_path, "all_performance.parquet")) +#' cross_drug <- tibble::tibble( +#' drug_or_class = c("AMP", "AMP", "CIP", "CIP", "NAL", "NAL"), +#' tested_on = c("CIP", "NAL", "AMP", "NAL", "AMP", "CIP"), +#' mcc = c(0.3, 0.2, 0.4, 0.25, 0.15, 0.35) +#' ) +#' performance <- tibble::tibble( +#' drug_label = "drug", +#' drug_or_class = c("AMP", "CIP", "NAL"), +#' mcc = c(0.8, 0.7, 0.6) +#' ) +#' plotCrossDrug( +#' cross_test_performance_path = cross_drug, +#' drug_performance_path = performance, +#' metadata_path = system.file("extdata", package = "amRml") +#' ) +plotCrossDrug <- function( + cross_test_performance_path = ".", + drug_performance_path = ".", + metadata_path = "." +) { + if (!is.data.frame(cross_test_performance_path)) { + cross_test_performance_path <- arrow::read_parquet( + file.path(cross_test_performance_path, "cross_drug_perf.parquet") + ) + } + cross_drug <- cross_test_performance_path + if (!is.data.frame(drug_performance_path)) { + drug_performance_path <- arrow::read_parquet( + file.path(drug_performance_path, "all_performance.parquet") + ) + } + performance <- drug_performance_path + metadata <- arrow::read_parquet(file.path(metadata_path, "metadata.parquet")) ###################### CROSS DRUG Testing ############################# heatmap_df <- cross_drug |> @@ -634,6 +714,7 @@ plotCrossDrug <- function(cross_test_performance_path = ".", drug_performance_pa ) heat_colors <- colorRampPalette(RColorBrewer::brewer.pal(11, "RdBu"))(100) + max_val <- max(abs(mat), na.rm = TRUE) # ---- Convert annotations ---- ha_row <- ComplexHeatmap::rowAnnotation( @@ -698,10 +779,12 @@ plotCrossDrug <- function(cross_test_performance_path = ".", drug_performance_pa #' @export #' #' @examples +#' \dontrun{ #' plotStratifiedPerf("year", #' stratified_performance_path = "data/Campylobacter/ML_year_performance", #' stratified_cross_performance_path = "data/Campylobacter/cross_test_ML_year_performance" #' ) +#' } plotStratifiedPerf <- function(year_or_country = "year", stratified_performance_path = ".", stratified_cross_performance_path = ".") { @@ -819,10 +902,12 @@ plotStratifiedPerf <- function(year_or_country = "year", #' @export #' #' @examples +#' \dontrun{ #' plotMDR( #' MDR_performance_path = "data/Campylobacter/MDR_ML_performance", MDR_top_feature_path = "data/Campylobacter/MDR_ML_top_features", #' MDR_pred_path = "data/Campylobacter/MDR_ML_pred" #' ) +#' } plotMDR <- function(MDR_performance_path = ".", MDR_top_feature_path = ".", MDR_pred_path = ".") { MDR_perf <- arrow::read_parquet(file.path(MDR_performance_path, "MDR_perf.parquet")) @@ -1028,17 +1113,29 @@ plotMDR <- function(MDR_performance_path = ".", MDR_top_feature_path = ".", #' Creates boxplots comparing performance (MCC) between real and shuffled labels #' across feature types. #' -#' @param metadata_path Character. Path to `metadata.parquet`. -#' @param performance_path Character. Path to `all_performance.parquet`. +#' @param metadata_path Character. Unused; retained for backward compatibility. +#' @param performance_path A performance tibble (with `feature_type`, +#' `feature_subtype`, `mcc`, and `shuffled` columns), or a directory path +#' containing `all_performance.parquet`. #' #' @return A ggplot object. #' @export #' #' @examples -#' plotShuffleVsReal(metadata_path = "data/Campylobacter/", performance_path = "data/Campylobacter/ML_performance") +#' performance <- tibble::tibble( +#' feature_type = rep(c("genes", "proteins"), each = 4), +#' feature_subtype = rep(c("binary", "counts"), times = 4), +#' mcc = c(0.7, 0.6, 0.65, 0.55, 0.1, 0.05, 0.08, 0.02), +#' shuffled = rep(c(FALSE, TRUE), each = 4) +#' ) +#' plotShuffleVsReal(performance_path = performance) plotShuffleVsReal <- function(metadata_path = ".", performance_path = ".") { - metadata <- arrow::read_parquet(file.path(metadata_path, "metadata.parquet")) - performance <- arrow::read_parquet(file.path(performance_path, "all_performance.parquet")) + if (!is.data.frame(performance_path)) { + performance_path <- arrow::read_parquet( + file.path(performance_path, "all_performance.parquet") + ) + } + performance <- performance_path performance |> dplyr::mutate( @@ -1086,10 +1183,12 @@ plotShuffleVsReal <- function(metadata_path = ".", performance_path = ".") { #' @export #' #' @examples +#' \dontrun{ #' plotTopClusters( #' top_feat_path = "data/Campylobacter/ML_top_features", cluster_feature_path = "data/Campylobacter/", #' protein_names_path = "data/Campylobacter/", top_n = 10 #' ) +#' } plotTopClusters <- function(top_feat_path = ".", cluster_feature_path = ".", protein_names_path = ".", top_n = 10) { ################### Top features ######################### diff --git a/R/prep_ml.R b/R/prep_ml.R index 9ce2ef2..efd4a97 100644 --- a/R/prep_ml.R +++ b/R/prep_ml.R @@ -96,18 +96,18 @@ loadMLInputTibble <- function(parquet_path) { "genome_drug.resistant_phenotype" %in% colnames(long_tibble), "resistant_classes" %in% colnames(long_tibble) )) { - stop(paste( - "The tibble loaded from `parquet_path` must have a target", - "variable column named either `genome_drug.resistant_phenotype` or", + stop( + "The tibble loaded from `parquet_path` must have a target ", + "variable column named either `genome_drug.resistant_phenotype` or ", "`resistant_classes`, but not both." - )) + ) } if (!all(c("genome_id", "feature_id", "value") %in% colnames(long_tibble))) { - stop(paste( - "The data is missing one or more of the following required", + stop( + "The data is missing one or more of the following required ", "columns: 'genome_id', 'feature_id', 'value'." - )) + ) } target_var <- .getTargetVarName(long_tibble) diff --git a/R/run_ML.R b/R/run_ML.R index e85d42e..a28d821 100644 --- a/R/run_ML.R +++ b/R/run_ML.R @@ -326,7 +326,7 @@ createMLinputList <- function(path, # Case A: stratified -> prefix before the stratify label if (!is.na(i)) { if (i - 1 >= 1) { - return(paste(x[1:(i - 1)], collapse = "_")) + return(paste(x[seq_len(i - 1)], collapse = "_")) } return("") } @@ -659,7 +659,7 @@ runMDRmodels <- function(path, if (inherits(res, "try-error")) { warning( "Model failed for: ", output_prefix, - "\n Error: ", attr(res, "condition")$message + "\n ", attr(res, "condition")$message ) return(NULL) } @@ -964,7 +964,7 @@ runMLmodels <- function(path, if (inherits(res, "try-error")) { warning( "Model failed for: ", output_prefix, - "\n Error: ", attr(res, "condition")$message + "\n ", attr(res, "condition")$message ) return(NULL) } diff --git a/R/run_ml_pipeline.R b/R/run_ml_pipeline.R index 5c47271..cab507d 100644 --- a/R/run_ml_pipeline.R +++ b/R/run_ml_pipeline.R @@ -136,10 +136,10 @@ runMLPipeline <- function( } if (model != "LR" & multi_class) { - stop(paste( - "Multi-class classification can be performed only with", + stop( + "Multi-class classification can be performed only with ", "logistic regression in our current implementation." - )) + ) } if (verbose) { @@ -158,7 +158,11 @@ runMLPipeline <- function( dplyr::mutate(prop = n / sum(n)) if (verbose) { - print(rs_props_original_ml_input_tibble) + rs_props_summary <- paste( + utils::capture.output(rs_props_original_ml_input_tibble), + collapse = "\n" + ) + message(rs_props_summary) } res_prop_original_ml_input_tibble <- rs_props_original_ml_input_tibble |> @@ -183,19 +187,19 @@ runMLPipeline <- function( external_test_data <- TRUE if (multi_class) { - stop(paste( - "Multi-class classification cannot be performed with", + stop( + "Multi-class classification cannot be performed with ", "external test data with the current implementation." - )) + ) } num_obs_test_data <- nrow(test_data) if (verbose) { - message(paste( - "`test_data` has", num_obs_test_data, - "observations of", getNumFeat(test_data), "features." - )) + message( + "`test_data` has ", num_obs_test_data, + " observations of ", getNumFeat(test_data), " features." + ) } rs_props_test_data <- test_data |> @@ -203,7 +207,11 @@ runMLPipeline <- function( dplyr::mutate(prop = n / sum(n)) if (verbose) { - print(rs_props_test_data) + test_data_props_summary <- paste( + utils::capture.output(rs_props_test_data), + collapse = "\n" + ) + message(test_data_props_summary) } # Get the proportion of resistant genomes in the external `test_data`. @@ -248,12 +256,12 @@ runMLPipeline <- function( # Stratify by R/S phenotype res_train_indices <- sample(res_indices, n_train_res) - sus_train_indices <- (1:n_train)[!(1:n_train %in% res_indices)] + sus_train_indices <- seq_len(n_train)[!(seq_len(n_train) %in% res_indices)] train_indices <- c(res_train_indices, sus_train_indices) - val_indices <- (1:nrow(original_ml_input_tibble))[!(1:nrow( + val_indices <- seq_len(nrow(original_ml_input_tibble))[!(seq_len(nrow( original_ml_input_tibble - ) %in% train_indices)] + )) %in% train_indices)] data_split[[2]] <- as.integer(train_indices) data_split[[3]] <- as.integer(val_indices) @@ -261,7 +269,7 @@ runMLPipeline <- function( n_train <- nrow(original_ml_input_tibble) n_val <- 0 - data_split[[2]] <- 1:nrow(original_ml_input_tibble) + data_split[[2]] <- seq_len(nrow(original_ml_input_tibble)) data_split[[3]] <- NA data_split[[4]] <- tibble::tibble(id = "Resample1") } @@ -314,7 +322,7 @@ runMLPipeline <- function( nmcc <- .calculatenMCC(test_data_plus_predictions) if (verbose) { - message(paste("Matthews correlation coefficient:", mcc, "| nMCC:", nmcc)) + message("Matthews correlation coefficient: ", mcc, " | nMCC: ", nmcc) } top_feat_tibble <- extractTopFeats(fit, @@ -342,7 +350,7 @@ runMLPipeline <- function( # `top_feats`. n_random_feats <- n_top_feats - length(top_feats) - random_indices <- sample(1:length(zero_vi_feats), n_random_feats) + random_indices <- sample(seq_along(zero_vi_feats), n_random_feats) random_feats <- zero_vi_feats[random_indices] top_feats <- c(top_feats, random_feats) diff --git a/man/calculateEvalMets.Rd b/man/calculateEvalMets.Rd index 0506645..5b61252 100644 --- a/man/calculateEvalMets.Rd +++ b/man/calculateEvalMets.Rd @@ -15,8 +15,8 @@ F1 score, AUPRC, balanced accuracy, nMCC, and log2(AUPRC/prior) } \description{ Returns the F1 score, area under the precision-recall curve (AUPRC), balanced -accuracy, normalized (to a 0 to 1 scale instead of -1 to 1) Matthews -correlation coefficient (nMCC), and log2(AUPRC/prior) based on the AMR +accuracy, Matthews correlation coefficient (MCC), normalized MCC (nMCC), +and log2(AUPRC/prior) based on the AMR phenotype predictions by an ML model compared against the actual values. } \examples{ diff --git a/man/demo_fit.Rd b/man/demo_fit.Rd index 0810296..9f29307 100644 --- a/man/demo_fit.Rd +++ b/man/demo_fit.Rd @@ -16,4 +16,8 @@ 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 index 021834d..79f685a 100644 --- a/man/demo_ml_tibble.Rd +++ b/man/demo_ml_tibble.Rd @@ -19,4 +19,8 @@ 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/dot-calculateMCC.Rd b/man/dot-calculateMCC.Rd new file mode 100644 index 0000000..edcc550 --- /dev/null +++ b/man/dot-calculateMCC.Rd @@ -0,0 +1,20 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/core_ml.R +\name{.calculateMCC} +\alias{.calculateMCC} +\title{.calculateMCC()} +\usage{ +.calculateMCC(test_data_plus_predictions) +} +\arguments{ +\item{test_data_plus_predictions}{Test data (tibble) with an added column for +predicted phenotype labels, such as the output of \code{predictML()}} +} +\value{ +Matthews correlation coefficient (MCC), range -1 to 1 +} +\description{ +Returns the Matthews correlation coefficient (MCC) +based on the AMR phenotype predictions by an +ML model compared against the actual values. +} diff --git a/man/dot-calculatenMCC.Rd b/man/dot-calculatenMCC.Rd index de13978..8763081 100644 --- a/man/dot-calculatenMCC.Rd +++ b/man/dot-calculatenMCC.Rd @@ -11,11 +11,10 @@ predicted phenotype labels, such as the output of \code{predictML()}} } \value{ -Normalized (to a 0 to 1 scale instead of -1 to 1) Matthews -correlation coefficient (nMCC) +Normalized Matthews correlation coefficient (nMCC), range 0 to 1 } \description{ -Returns the normalized (to a 0 to 1 scale instead of -1 to 1) Matthews -correlation coefficient (nMCC) based on the AMR phenotype predictions by an -ML model compared against the actual values. +Returns the normalized (0 to 1) Matthews correlation coefficient (nMCC) +based on the AMR phenotype predictions by an ML model compared against +the actual values. } diff --git a/man/dot-parquet2LOODrugMatrix.Rd b/man/dot-parquet2LOODrugMatrix.Rd new file mode 100644 index 0000000..9ad18ce --- /dev/null +++ b/man/dot-parquet2LOODrugMatrix.Rd @@ -0,0 +1,20 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/generate_matrices_ml.R +\name{.parquet2LOODrugMatrix} +\alias{.parquet2LOODrugMatrix} +\title{Build leave-one-out (LOO) merged parquet matrices from drug parquet files.} +\usage{ +.parquet2LOODrugMatrix(path, verbosity = c("minimal", "debug")) +} +\arguments{ +\item{path}{Character. Base directory containing stratified parquet matrices. +Expected subdirs: matrix/..} + +\item{verbosity}{Character. "minimal" or "debug"; when "debug", prints detailed steps.} +} +\value{ +Invisibly returns a tibble with paths of created LOO parquet files. +} +\description{ +Build leave-one-out (LOO) merged parquet matrices from drug parquet files. +} diff --git a/man/plotBaselineComparison.Rd b/man/plotBaselineComparison.Rd index b212ab2..8eab040 100644 --- a/man/plotBaselineComparison.Rd +++ b/man/plotBaselineComparison.Rd @@ -2,32 +2,34 @@ % Please edit documentation in R/plot_ml.R \name{plotBaselineComparison} \alias{plotBaselineComparison} -\title{plotBaselineComparison()} +\title{Compare Baseline Performance With and Without Shuffled Labels} \usage{ plotBaselineComparison(non_shuffled_label_results, shuffled_label_results) } \arguments{ -\item{non_shuffled_label_results}{Output of \code{runMLPipeline()} -(\code{shuffle_labels = FALSE})} +\item{non_shuffled_label_results}{Output of \code{runMLPipeline(shuffle_labels = FALSE)}} -\item{shuffled_label_results}{Output of \code{runMLPipeline()} -(\code{shuffle_labels = TRUE})} +\item{shuffled_label_results}{Output of \code{runMLPipeline(shuffle_labels = TRUE)}} } \value{ -A bar plot with balanced accuracy comparisons per antibiotic +A base R barplot comparing balanced accuracy across models. } \description{ -Generates a bar plot that compares model performance with and without -randomly shuffled AMR phenotype labels. +Produces a bar plot comparing balanced accuracy for each antibiotic using +true AMR labels vs. randomly shuffled labels. } \examples{ -non_shuffled <- tibble::tibble( - antibiotic = c("AMP", "CIP", "CRO"), - bal_acc = c(0.88, 0.81, 0.92) +non_shuffled <- list( + performance_tibble = tibble::tibble( + antibiotic = c("AMP", "CIP", "CRO"), + bal_acc = c(0.88, 0.81, 0.92) + ) ) -shuffled <- tibble::tibble( - antibiotic = c("AMP", "CIP", "CRO"), - bal_acc = c(0.52, 0.49, 0.55) +shuffled <- list( + performance_tibble = tibble::tibble( + antibiotic = c("AMP", "CIP", "CRO"), + bal_acc = c(0.52, 0.49, 0.55) + ) ) plotBaselineComparison(non_shuffled, shuffled) } diff --git a/man/plotCM.Rd b/man/plotCM.Rd new file mode 100644 index 0000000..673bca0 --- /dev/null +++ b/man/plotCM.Rd @@ -0,0 +1,26 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plot_ml.R +\name{plotCM} +\alias{plotCM} +\title{Plot a Confusion Matrix Heatmap} +\usage{ +plotCM(test_data_plus_predictions) +} +\arguments{ +\item{test_data_plus_predictions}{A tibble containing true and predicted +phenotype labels (e.g. the output of \code{predictML()}), or a path to a TSV file +containing the same.} +} +\value{ +A heatmap (\code{ggplot2} object) showing the confusion matrix. +} +\description{ +Produces a heatmap visualization of the confusion matrix for AMR predictions. +} +\examples{ +data(demo_fit) +data(demo_ml_tibble) +preds <- predictML(demo_fit, demo_ml_tibble) +plotCM(preds) + +} diff --git a/man/plotCrossDrug.Rd b/man/plotCrossDrug.Rd new file mode 100644 index 0000000..1415f1e --- /dev/null +++ b/man/plotCrossDrug.Rd @@ -0,0 +1,47 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plot_ml.R +\name{plotCrossDrug} +\alias{plotCrossDrug} +\title{Plot cross-drug generalization heatmap} +\usage{ +plotCrossDrug( + cross_test_performance_path = ".", + drug_performance_path = ".", + metadata_path = "." +) +} +\arguments{ +\item{cross_test_performance_path}{A cross-drug performance tibble (with +\code{drug_or_class}, \code{tested_on}, and \code{mcc} columns), or a directory path +containing \code{cross_drug_perf.parquet}.} + +\item{drug_performance_path}{A performance tibble (with \code{drug_label}, +\code{drug_or_class}, and \code{mcc} columns), or a directory path containing +\code{all_performance.parquet}.} + +\item{metadata_path}{Character. Path to a directory containing \code{metadata.parquet}.} +} +\value{ +A ComplexHeatmap object. +} +\description{ +Creates a heatmap showing cross-drug model performance (MCC), where models +trained on one drug are evaluated on another. +} +\examples{ +cross_drug <- tibble::tibble( + drug_or_class = c("AMP", "AMP", "CIP", "CIP", "NAL", "NAL"), + tested_on = c("CIP", "NAL", "AMP", "NAL", "AMP", "CIP"), + mcc = c(0.3, 0.2, 0.4, 0.25, 0.15, 0.35) +) +performance <- tibble::tibble( + drug_label = "drug", + drug_or_class = c("AMP", "CIP", "NAL"), + mcc = c(0.8, 0.7, 0.6) +) +plotCrossDrug( + cross_test_performance_path = cross_drug, + drug_performance_path = performance, + metadata_path = system.file("extdata", package = "amRml") +) +} diff --git a/man/plotDefaultEval.Rd b/man/plotDefaultEval.Rd deleted file mode 100644 index ca2dcfe..0000000 --- a/man/plotDefaultEval.Rd +++ /dev/null @@ -1,48 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/plot_ml.R -\name{plotDefaultEval} -\alias{plotDefaultEval} -\title{plotDefaultEval()} -\usage{ -plotDefaultEval( - default_eval_tibble, - x_default_eval = "train_prop", - y_default_eval = "avg_f1_score", - xlab = "Train data proportion", - ylab = "Average F1 score" -) -} -\arguments{ -\item{default_eval_tibble}{Output of \code{findOptimalMLDefaults()}} - -\item{x_default_eval}{\link[rlang:chr]{chr} x value of default evaluation plot: "train_prop" -or "n_fold"} - -\item{y_default_eval}{\link[rlang:chr]{chr} y value of default evaluation plot. It can be -"avg_runtime_sec" or one of the following performance metrics: -"avg_f1_score", "avg_log2_apop", "avg_bal_acc", or "avg_nmcc"} - -\item{xlab}{\link[rlang:chr]{chr} Label for x axis} - -\item{ylab}{\link[rlang:chr]{chr} Label for y axis} -} -\value{ -A \code{ggplot2} scatterplot (performance metric or runtime vs. -\code{train_prop} or \code{n_fold}), colored by model -} -\description{ -Plots performance metric or runtime vs. training data proportion or number -of cross-validation folds, colored by model. -} -\examples{ -default_eval <- tibble::tibble( - train_prop = c(0.5, 0.6, 0.7, 0.5, 0.6, 0.7), - n_fold = rep(5, 6), - model = rep(c("LR", "RF"), each = 3), - avg_f1_score = c(0.72, 0.78, 0.83, 0.70, 0.75, 0.80) -) -plotDefaultEval(default_eval, - x_default_eval = "train_prop", - y_default_eval = "avg_f1_score" -) -} diff --git a/man/plotDrugDist.Rd b/man/plotDrugDist.Rd new file mode 100644 index 0000000..9188dfe --- /dev/null +++ b/man/plotDrugDist.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plot_ml.R +\name{plotDrugDist} +\alias{plotDrugDist} +\title{Plot drug phenotype distribution} +\usage{ +plotDrugDist(metadata_path = ".") +} +\arguments{ +\item{metadata_path}{Character. Path to directory containing \code{metadata.parquet}.} +} +\value{ +A ggplot object. +} +\description{ +Reads metadata and generates a stacked bar plot showing counts of resistant +and susceptible phenotypes per antibiotic. +} +\examples{ +plotDrugDist(metadata_path = system.file("extdata", package = "amRml")) +} diff --git a/man/plotDrugPerf.Rd b/man/plotDrugPerf.Rd new file mode 100644 index 0000000..6b83cdb --- /dev/null +++ b/man/plotDrugPerf.Rd @@ -0,0 +1,38 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plot_ml.R +\name{plotDrugPerf} +\alias{plotDrugPerf} +\title{Plot drug-level model performance} +\usage{ +plotDrugPerf(metadata_path = ".", performance_path = ".") +} +\arguments{ +\item{metadata_path}{Character. Path to a directory containing \code{metadata.parquet}.} + +\item{performance_path}{A performance tibble (with \code{drug_label}, \code{shuffled}, +\code{drug_or_class}, \code{feature_type}, \code{feature_subtype}, and \code{mcc} columns), or a +directory path containing \code{all_performance.parquet}.} +} +\value{ +A patchwork ggplot object combining multiple panels. +} +\description{ +Generates heatmaps and ridge plots summarizing model performance (MCC) +across drugs and feature types. +} +\examples{ +# Several models per drug x feature type so the ridge densities can be drawn. +performance <- tidyr::expand_grid( + drug_or_class = c("AMP", "CIP", "NAL"), + feature_type = c("genes", "proteins"), + feature_subtype = c("binary", "counts"), + replicate = 1:8 +) +performance$drug_label <- "drug" +performance$shuffled <- FALSE +performance$mcc <- 0.6 + 0.3 * sin(seq_len(nrow(performance))) +plotDrugPerf( + metadata_path = system.file("extdata", package = "amRml"), + performance_path = performance +) +} diff --git a/man/plotMDR.Rd b/man/plotMDR.Rd new file mode 100644 index 0000000..6573e4d --- /dev/null +++ b/man/plotMDR.Rd @@ -0,0 +1,34 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plot_ml.R +\name{plotMDR} +\alias{plotMDR} +\title{Plot multi-drug resistance (MDR) model performance} +\usage{ +plotMDR( + MDR_performance_path = ".", + MDR_top_feature_path = ".", + MDR_pred_path = "." +) +} +\arguments{ +\item{MDR_performance_path}{Character. Path to \code{MDR_perf.parquet}.} + +\item{MDR_top_feature_path}{Character. Path to \code{MDR_top_features.parquet}.} + +\item{MDR_pred_path}{Character. Path to \code{MDR_pred.parquet}.} +} +\value{ +A list of ggplot objects. +} +\description{ +Generates violin plots of performance, feature importance summaries, +and prediction confusion-style visualizations for MDR models. +} +\examples{ +\dontrun{ +plotMDR( + MDR_performance_path = "data/Campylobacter/MDR_ML_performance", MDR_top_feature_path = "data/Campylobacter/MDR_ML_top_features", + MDR_pred_path = "data/Campylobacter/MDR_ML_pred" +) +} +} diff --git a/man/plotPRC.Rd b/man/plotPRC.Rd index 3bad222..7581ca3 100644 --- a/man/plotPRC.Rd +++ b/man/plotPRC.Rd @@ -2,37 +2,29 @@ % Please edit documentation in R/plot_ml.R \name{plotPRC} \alias{plotPRC} -\title{plotPRC()} +\title{Plot a Precision-Recall Curve} \usage{ plotPRC(test_data_plus_predictions) } \arguments{ -\item{test_data_plus_predictions}{Test data (tibble) with an added column for -predicted phenotype labels, such as the output of \code{predict()}.} +\item{test_data_plus_predictions}{A tibble of test data with added prediction +columns (e.g. the output of \code{predictML()} or \code{runMLmodels(return_pred=TRUE)}), +or a path to a TSV file containing the same.} } \value{ -A precision-recall curve as a \code{ggplot2} object +A \code{ggplot2} object showing the precision-recall curve. } \description{ -Plots the precision-recall curve given a set of test data plus predicted AMR -phenotypes. +Generates a precision-recall curve (PRC) for AMR phenotype prediction results. +} +\details{ +The function uses \code{yardstick::pr_curve()} to compute the PR curve and then +visualizes it using \code{ggplot2}. } \examples{ -preds <- tibble::tibble( - genome_id = paste0("g", 1:10), - genome_drug.resistant_phenotype = factor( - rep(c("Resistant", "Susceptible"), each = 5), - levels = c("Resistant", "Susceptible") - ), - .pred_class = factor( - c( - "Resistant", "Resistant", "Susceptible", "Resistant", "Susceptible", - "Susceptible", "Resistant", "Susceptible", "Susceptible", "Resistant" - ), - levels = c("Resistant", "Susceptible") - ), - .pred_Resistant = c(0.9, 0.8, 0.4, 0.7, 0.3, 0.2, 0.6, 0.1, 0.2, 0.55), - .pred_Susceptible = c(0.1, 0.2, 0.6, 0.3, 0.7, 0.8, 0.4, 0.9, 0.8, 0.45) -) +data(demo_fit) +data(demo_ml_tibble) +preds <- predictML(demo_fit, demo_ml_tibble) plotPRC(preds) + } diff --git a/man/plotROC.Rd b/man/plotROC.Rd new file mode 100644 index 0000000..45e0f54 --- /dev/null +++ b/man/plotROC.Rd @@ -0,0 +1,26 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plot_ml.R +\name{plotROC} +\alias{plotROC} +\title{Plot a Receiver Operating Characteristic (ROC) Curve} +\usage{ +plotROC(test_data_plus_predictions) +} +\arguments{ +\item{test_data_plus_predictions}{A tibble of test data with prediction +columns (e.g. the output of \code{predictML()} or \code{runMLmodels(return_pred=TRUE)}), +or a path to a TSV file containing the same.} +} +\value{ +A ROC curve plotted using \code{ggplot2::autoplot()}. +} +\description{ +Generates a ROC curve for AMR phenotype prediction results. +} +\examples{ +data(demo_fit) +data(demo_ml_tibble) +preds <- predictML(demo_fit, demo_ml_tibble) +plotROC(preds) + +} diff --git a/man/plotShuffleVsReal.Rd b/man/plotShuffleVsReal.Rd new file mode 100644 index 0000000..d000e89 --- /dev/null +++ b/man/plotShuffleVsReal.Rd @@ -0,0 +1,31 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plot_ml.R +\name{plotShuffleVsReal} +\alias{plotShuffleVsReal} +\title{Compare shuffled vs real model performance} +\usage{ +plotShuffleVsReal(metadata_path = ".", performance_path = ".") +} +\arguments{ +\item{metadata_path}{Character. Unused; retained for backward compatibility.} + +\item{performance_path}{A performance tibble (with \code{feature_type}, +\code{feature_subtype}, \code{mcc}, and \code{shuffled} columns), or a directory path +containing \code{all_performance.parquet}.} +} +\value{ +A ggplot object. +} +\description{ +Creates boxplots comparing performance (MCC) between real and shuffled labels +across feature types. +} +\examples{ +performance <- tibble::tibble( + feature_type = rep(c("genes", "proteins"), each = 4), + feature_subtype = rep(c("binary", "counts"), times = 4), + mcc = c(0.7, 0.6, 0.65, 0.55, 0.1, 0.05, 0.08, 0.02), + shuffled = rep(c(FALSE, TRUE), each = 4) +) +plotShuffleVsReal(performance_path = performance) +} diff --git a/man/plotStratifiedPerf.Rd b/man/plotStratifiedPerf.Rd new file mode 100644 index 0000000..1fe5854 --- /dev/null +++ b/man/plotStratifiedPerf.Rd @@ -0,0 +1,34 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plot_ml.R +\name{plotStratifiedPerf} +\alias{plotStratifiedPerf} +\title{Plot stratified model performance} +\usage{ +plotStratifiedPerf( + year_or_country = "year", + stratified_performance_path = ".", + stratified_cross_performance_path = "." +) +} +\arguments{ +\item{year_or_country}{Character. Either "year" or "country".} + +\item{stratified_performance_path}{Character. Path to stratified performance files.} + +\item{stratified_cross_performance_path}{Character. Path to cross-stratified performance files.} +} +\value{ +A ggplot object. +} +\description{ +Visualizes model performance (MCC) stratified by year or country, +comparing within-group vs cross-group evaluation. +} +\examples{ +\dontrun{ +plotStratifiedPerf("year", + stratified_performance_path = "data/Campylobacter/ML_year_performance", + stratified_cross_performance_path = "data/Campylobacter/cross_test_ML_year_performance" +) +} +} diff --git a/man/plotTopClusters.Rd b/man/plotTopClusters.Rd new file mode 100644 index 0000000..0420d67 --- /dev/null +++ b/man/plotTopClusters.Rd @@ -0,0 +1,37 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plot_ml.R +\name{plotTopClusters} +\alias{plotTopClusters} +\title{Plot top contributing feature clusters} +\usage{ +plotTopClusters( + top_feat_path = ".", + cluster_feature_path = ".", + protein_names_path = ".", + top_n = 10 +) +} +\arguments{ +\item{top_feat_path}{Character. Path to \code{all_top_features.parquet}.} + +\item{cluster_feature_path}{Character. Path to \code{cluster_feature.parquet}.} + +\item{protein_names_path}{Character. Path to \code{protein_names.parquet}.} + +\item{top_n}{Integer. Number of top features to retain per model.} +} +\value{ +A ggplot object. +} +\description{ +Identifies top contributing clusters across feature types and drugs, +and visualizes their relative contributions. +} +\examples{ +\dontrun{ +plotTopClusters( + top_feat_path = "data/Campylobacter/ML_top_features", cluster_feature_path = "data/Campylobacter/", + protein_names_path = "data/Campylobacter/", top_n = 10 +) +} +} diff --git a/man/plotTopFeatsVI.Rd b/man/plotTopFeatsVI.Rd index 9bd920b..21fbb48 100644 --- a/man/plotTopFeatsVI.Rd +++ b/man/plotTopFeatsVI.Rd @@ -2,23 +2,27 @@ % Please edit documentation in R/plot_ml.R \name{plotTopFeatsVI} \alias{plotTopFeatsVI} -\title{plotTopFeatsVI()} +\title{Plot Top Feature Importances} \usage{ -plotTopFeatsVI(fit, n_top_feats = 10) +plotTopFeatsVI(topfeat, n_top_feats = 10) } \arguments{ -\item{fit}{Best model fit, such as the output of \code{fitBestModel()}} +\item{topfeat}{A tibble of feature importance scores (e.g. the output of +\code{extractTopFeats()} or \code{runMLmodels()}), with \code{Variable}, \code{Importance}, and +\code{Sign} columns, or a path to a TSV file containing the same.} -\item{n_top_feats}{\link[pillar:num]{num} Number of top features to plot} +\item{n_top_feats}{Number of top features to display (default: 10).} } \value{ -Variable importance plot (a \code{ggplot2} object) +A bar plot of variable importance (\code{ggplot2} object). } \description{ -Generates a plot showing the top features and their variable importance -scores. +Creates a bar plot showing the most important features affecting +AMR phenotype predictions. } \examples{ data(demo_fit) -plotTopFeatsVI(demo_fit, n_top_feats = 10) +top_feats <- extractTopFeats(demo_fit, n_top_feats = 10) +plotTopFeatsVI(top_feats, n_top_feats = 10) + } diff --git a/man/runIFE.Rd b/man/runIFE.Rd index 2e39c9e..f817ba9 100644 --- a/man/runIFE.Rd +++ b/man/runIFE.Rd @@ -4,7 +4,7 @@ \alias{runIFE} \title{runIFE Removes top features identified by ML models and retrains iteratively; -returns nMCC at each iteration.} +returns MCC at each iteration.} \usage{ runIFE( ml_input_tibble, @@ -49,7 +49,7 @@ list along with top features removed per iteration if \code{return_feats = TRUE} \description{ runIFE Removes top features identified by ML models and retrains iteratively; -returns nMCC at each iteration. +returns MCC at each iteration. } \examples{ data(demo_ml_tibble) diff --git a/vignettes/intro.Rmd b/vignettes/intro.Rmd index cbe4f2b..b0c5182 100644 --- a/vignettes/intro.Rmd +++ b/vignettes/intro.Rmd @@ -196,7 +196,7 @@ plotPRC(results$pred) ``` ```{r plot-vi, fig.width = 5, fig.height = 4} -plotTopFeatsVI(results$fit, n_top_feats = 10) +plotTopFeatsVI(results$top_feat_tibble, n_top_feats = 10) ``` For a baseline comparison against random labels, fit a shuffled-label pipeline and compare: From 7177ab2829ecb92f8d7fa8375a841d6b683ed473 Mon Sep 17 00:00:00 2001 From: Abhirupa Ghosh <100681585+AbhirupaGhosh@users.noreply.github.com> Date: Wed, 1 Jul 2026 17:15:29 -0600 Subject: [PATCH 08/17] plot enhancements --- R/plot_ml.R | 178 +++++++++++++++++++++++++++++----------------------- 1 file changed, 100 insertions(+), 78 deletions(-) diff --git a/R/plot_ml.R b/R/plot_ml.R index 7ab9f8c..569ff88 100644 --- a/R/plot_ml.R +++ b/R/plot_ml.R @@ -248,7 +248,7 @@ plotBaselineComparison <- function( baseline_comparison_barplot <- barplot(bal_acc_matrix, beside = TRUE, legend.text = TRUE, col = c("skyblue", "lightpink"), - ylab = "Balanced Accuracy" + ylab = "Balanced accuracy" ) return(baseline_comparison_barplot) @@ -382,7 +382,7 @@ plotDrugDist <- function(metadata_path = ".") { ##################### phenotype distribution (drugs) ######################### drug_dist <- metadata |> dplyr::distinct( - genome_drug.genome_id, + genome.genome_id, genome_drug.antibiotic, drug_abbr, genome_drug.resistant_phenotype @@ -434,7 +434,7 @@ plotDrugDist <- function(metadata_path = ".") { #' @param metadata_path Character. Path to a directory containing `metadata.parquet`. #' @param performance_path A performance tibble (with `drug_label`, `shuffled`, #' `drug_or_class`, `feature_type`, `feature_subtype`, and `mcc` columns), or a -#' directory path containing `all_performance.parquet`. +#' directory path containing `all_perf.parquet`. #' #' @return A patchwork ggplot object combining multiple panels. #' @export @@ -459,13 +459,13 @@ plotDrugPerf <- function(metadata_path = ".", performance_path = ".") { if (!is.data.frame(performance_path)) { performance_path <- arrow::read_parquet( - file.path(performance_path, "all_performance.parquet") + file.path(performance_path, "all_perf.parquet") ) } performance <- performance_path plot_df <- metadata |> - dplyr::distinct(genome_drug.genome_id, genome_drug.antibiotic, drug_abbr) |> + dplyr::distinct(genome.genome_id, genome_drug.antibiotic, drug_abbr) |> dplyr::count(genome_drug.antibiotic, drug_abbr, name = "total") ######################## drug performances ################################# @@ -548,7 +548,7 @@ plotDrugPerf <- function(metadata_path = ".", performance_path = ".") { ) + ggplot2::geom_point( position = position_jitter(height = 0.1), - size = 2, + size = 1.5, alpha = 0.8, aes(color = feature_type) ) + @@ -590,11 +590,11 @@ plotDrugPerf <- function(metadata_path = ".", performance_path = ".") { #' trained on one drug are evaluated on another. #' #' @param cross_test_performance_path A cross-drug performance tibble (with -#' `drug_or_class`, `tested_on`, and `mcc` columns), or a directory path +#' `drug_or_class`, `test_drug`, and `mcc` columns), or a directory path #' containing `cross_drug_perf.parquet`. #' @param drug_performance_path A performance tibble (with `drug_label`, #' `drug_or_class`, and `mcc` columns), or a directory path containing -#' `all_performance.parquet`. +#' `all_perf.parquet`. #' @param metadata_path Character. Path to a directory containing `metadata.parquet`. #' #' @return A ComplexHeatmap object. @@ -603,7 +603,7 @@ plotDrugPerf <- function(metadata_path = ".", performance_path = ".") { #' @examples #' cross_drug <- tibble::tibble( #' drug_or_class = c("AMP", "AMP", "CIP", "CIP", "NAL", "NAL"), -#' tested_on = c("CIP", "NAL", "AMP", "NAL", "AMP", "CIP"), +#' test_drug = c("CIP", "NAL", "AMP", "NAL", "AMP", "CIP"), #' mcc = c(0.3, 0.2, 0.4, 0.25, 0.15, 0.35) #' ) #' performance <- tibble::tibble( @@ -624,12 +624,12 @@ plotCrossDrug <- function( if (!is.data.frame(cross_test_performance_path)) { cross_test_performance_path <- arrow::read_parquet( file.path(cross_test_performance_path, "cross_drug_perf.parquet") - ) + ) } cross_drug <- cross_test_performance_path if (!is.data.frame(drug_performance_path)) { drug_performance_path <- arrow::read_parquet( - file.path(drug_performance_path, "all_performance.parquet") + file.path(drug_performance_path, "all_perf.parquet") ) } performance <- drug_performance_path @@ -637,61 +637,63 @@ plotCrossDrug <- function( ###################### CROSS DRUG Testing ############################# heatmap_df <- cross_drug |> - # dplyr::filter(tested_on %in% (cross_drug |> dplyr::pull(drug_or_class))) |> - dplyr::group_by(drug_or_class, tested_on) |> + dplyr::filter(!is.na(drug), !is.na(test_drug)) |> + # dplyr::filter(test_drug %in% (cross_drug |> dplyr::pull(drug))) |> + dplyr::group_by(drug, test_drug) |> dplyr::summarise(median_mcc = median(mcc, na.rm = TRUE), .groups = "drop") same_drugs <- performance |> dplyr::filter( drug_label == "drug", drug_or_class %in% (cross_drug |> - dplyr::distinct(drug_or_class) |> + dplyr::distinct(drug) |> dplyr::pull()) ) |> dplyr::group_by(drug_or_class) |> dplyr::summarise(median_mcc = median(mcc, na.rm = TRUE), .groups = "drop") |> - dplyr::mutate(tested_on = drug_or_class) |> - dplyr::distinct(drug_or_class, tested_on, median_mcc) + dplyr::mutate(test_drug = drug_or_class) |> + dplyr::distinct(drug_or_class, test_drug, median_mcc) |> + dplyr::rename(drug = drug_or_class) heatmap_df <- heatmap_df |> dplyr::add_row(same_drugs) |> dplyr::left_join( metadata |> dplyr::distinct(drug_abbr, class_abbr), - by = c("drug_or_class" = "drug_abbr") + by = c("drug" = "drug_abbr") ) |> dplyr::rename("drug_class" = "class_abbr") |> dplyr::left_join( metadata |> dplyr::distinct(drug_abbr, class_abbr), - by = c("tested_on" = "drug_abbr") + by = c("test_drug" = "drug_abbr") ) # Row annotation (already similar to what you did) annotation_row <- heatmap_df |> - dplyr::distinct(drug_or_class, drug_class) |> - tibble::column_to_rownames("drug_or_class") + dplyr::distinct(drug, drug_class) |> + tibble::column_to_rownames("drug") # Column annotation annotation_col <- heatmap_df |> - dplyr::distinct(tested_on, class_abbr) |> - tibble::column_to_rownames("tested_on") + dplyr::distinct(test_drug, class_abbr) |> + tibble::column_to_rownames("test_drug") mat <- heatmap_df |> - dplyr::select(drug_or_class, tested_on, median_mcc) |> - tidyr::pivot_wider(names_from = tested_on, values_from = median_mcc) |> - tibble::column_to_rownames("drug_or_class") |> + dplyr::select(drug, test_drug, median_mcc) |> + tidyr::pivot_wider(names_from = test_drug, values_from = median_mcc) |> + tibble::column_to_rownames("drug") |> as.matrix() row_order <- heatmap_df |> - dplyr::distinct(drug_or_class, drug_class) |> - dplyr::arrange(drug_class, drug_or_class) |> - dplyr::pull(drug_or_class) + dplyr::distinct(drug, drug_class) |> + dplyr::arrange(drug_class, drug) |> + dplyr::pull(drug) col_order <- heatmap_df |> - dplyr::distinct(tested_on, class_abbr) |> - dplyr::arrange(class_abbr, tested_on) |> - dplyr::pull(tested_on) + dplyr::distinct(test_drug, class_abbr) |> + dplyr::arrange(class_abbr, test_drug) |> + dplyr::pull(test_drug) # mat[is.na(mat)] <- 0 mat <- mat[row_order, col_order] @@ -708,10 +710,16 @@ plotCrossDrug <- function( ) # Create ONE named color vector + # class_colors <- stats::setNames( + # scales::hue_pal()(length(classes)), + # classes + # ) + class_colors <- stats::setNames( - scales::hue_pal()(length(classes)), - classes - ) + colorRampPalette(RColorBrewer::brewer.pal(9, "Pastel1"))(length(classes)), + classes +) + heat_colors <- colorRampPalette(RColorBrewer::brewer.pal(11, "RdBu"))(100) max_val <- max(abs(mat), na.rm = TRUE) @@ -727,7 +735,8 @@ plotCrossDrug <- function( ha_col <- ComplexHeatmap::HeatmapAnnotation( class_abbr = annotation_col$class_abbr, col = list(class_abbr = class_colors), - show_annotation_name = FALSE, na_col = "grey3" + show_annotation_name = FALSE, na_col = "grey3", + annotation_legend_param = list(labels_gp = grid::gpar(fontsize = 14)) ) # ---- Color function (instead of breaks + palette) ---- @@ -746,6 +755,12 @@ plotCrossDrug <- function( column_order = col_order, left_annotation = ha_row, top_annotation = ha_col, + + + width = grid::unit(ncol(mat), "in"), + height = grid::unit(nrow(mat), "in"), + + show_row_names = TRUE, show_column_names = TRUE, column_title = "tested on", @@ -754,13 +769,15 @@ plotCrossDrug <- function( row_title_side = "right", row_names_gp = grid::gpar(fontsize = 14), column_names_gp = grid::gpar(fontsize = 14), - column_names_rot = 0, + column_names_rot = 45, # remove borders like pheatmap rect_gp = grid::gpar(col = NA), # legends - show_heatmap_legend = TRUE + show_heatmap_legend = TRUE, + + use_raster = FALSE ) cross_drug_hm @@ -801,54 +818,59 @@ plotStratifiedPerf <- function(year_or_country = "year", "_perf.parquet" ) )) - if (year_or_country == "year") { + # if (year_or_country == "year") { all <- perf |> - dplyr::rename("train_year" = "strat_value") |> - dplyr::mutate(test_year = train_year) |> + # dplyr::rename("train_year" = "strat_value") |> + dplyr::mutate(strat_value_test = strat_value) |> dplyr::select( drug_label, drug_or_class, - train_year, test_year, feature_type, feature_subtype, mcc + strat_value, strat_value_test, feature_type, feature_subtype, mcc ) |> dplyr::bind_rows(cross_test |> dplyr::select( drug_label, drug_or_class, - train_year, test_year, feature_type, + strat_value, strat_value_test, feature_type, feature_subtype, mcc )) |> dplyr::mutate(category = dplyr::if_else( - train_year == test_year, "same year bin", "different year bin" - )) - } else { - all <- perf |> - dplyr::rename("train_country" = "strat_value") |> - dplyr::mutate(test_country = train_country) |> - dplyr::select( - drug_label, drug_or_class, - train_country, test_country, - feature_type, feature_subtype, mcc - ) |> - dplyr::bind_rows(cross_test |> - dplyr::select( - drug_label, drug_or_class, - train_country, test_country, - feature_type, feature_subtype, mcc - )) |> - dplyr::mutate(category = dplyr::if_else( - train_country == test_country, "same country", "different country" + strat_value == strat_value_test, "same", "different" )) - } + # } else { + # all <- perf |> + # dplyr::rename("train_country" = "strat_value") |> + # dplyr::mutate(test_country = train_country) |> + # dplyr::select( + # drug_label, drug_or_class, + # train_country, test_country, + # feature_type, feature_subtype, mcc + # ) |> + # dplyr::bind_rows(cross_test |> + # dplyr::select( + # drug_label, drug_or_class, + # train_country, test_country, + # feature_type, feature_subtype, mcc + # )) |> + # dplyr::mutate(category = dplyr::if_else( + # train_country == test_country, "same country", "different country" + # )) + # } + + # fill_vals <- if (year_or_country == "year") { + # c( + # "same year bin" = "#b3cde3", + # "different year bin" = "#fbb4ae" + # ) + # } else { + # c( + # "same country" = "#b3cde3", + # "different country" = "#fbb4ae" + # ) + # } - fill_vals <- if (year_or_country == "year") { - c( - "same year bin" = "#b3cde3", - "different year bin" = "#fbb4ae" - ) - } else { - c( - "same country" = "#b3cde3", - "different country" = "#fbb4ae" + fill_vals <- c( + "same" = "#b3cde3", + "different" = "#fbb4ae" ) - } plot <- ggplot2::ggplot( all |> @@ -875,16 +897,16 @@ plotStratifiedPerf <- function(year_or_country = "year", fill = "Tested on" ) + ggplot2::theme( - axis.title = ggplot2::element_text(colour = "black", size = 10), - axis.text.x = ggplot2::element_text(angle = 45, hjust = 1, size = 10, colour = "black"), - axis.text.y = ggplot2::element_text(size = 10, colour = "black"), + axis.title = ggplot2::element_text(colour = "black", size = 14), + axis.text.x = ggplot2::element_text(angle = 45, hjust = 1, size = 14, colour = "black"), + axis.text.y = ggplot2::element_text(size = 14, colour = "black"), axis.title.y = ggplot2::element_blank(), legend.title = ggplot2::element_text(size = 12), legend.text = ggplot2::element_text(size = 10), legend.position = "bottom", plot.title = ggplot2::element_text(face = "bold"), panel.grid.minor = ggplot2::element_blank(), - plot.margin = margin(0, 0, 0, 0) + plot.margin = ggplot2::margin(0, 0, 0, 0) ) plot } @@ -1116,7 +1138,7 @@ plotMDR <- function(MDR_performance_path = ".", MDR_top_feature_path = ".", #' @param metadata_path Character. Unused; retained for backward compatibility. #' @param performance_path A performance tibble (with `feature_type`, #' `feature_subtype`, `mcc`, and `shuffled` columns), or a directory path -#' containing `all_performance.parquet`. +#' containing `all_perf.parquet`. #' #' @return A ggplot object. #' @export @@ -1132,7 +1154,7 @@ plotMDR <- function(MDR_performance_path = ".", MDR_top_feature_path = ".", plotShuffleVsReal <- function(metadata_path = ".", performance_path = ".") { if (!is.data.frame(performance_path)) { performance_path <- arrow::read_parquet( - file.path(performance_path, "all_performance.parquet") + file.path(performance_path, "all_perf.parquet") ) } performance <- performance_path From 49fc7e000e4f97bc0e25f5fa617c5892466cfc9a Mon Sep 17 00:00:00 2001 From: AbhirupaGhosh Date: Wed, 1 Jul 2026 23:17:14 +0000 Subject: [PATCH 09/17] Style code (GHA) --- R/plot_ml.R | 55 ++++++++++++++++++++++++----------------------------- 1 file changed, 25 insertions(+), 30 deletions(-) diff --git a/R/plot_ml.R b/R/plot_ml.R index 569ff88..8a89296 100644 --- a/R/plot_ml.R +++ b/R/plot_ml.R @@ -624,7 +624,7 @@ plotCrossDrug <- function( if (!is.data.frame(cross_test_performance_path)) { cross_test_performance_path <- arrow::read_parquet( file.path(cross_test_performance_path, "cross_drug_perf.parquet") - ) + ) } cross_drug <- cross_test_performance_path if (!is.data.frame(drug_performance_path)) { @@ -637,7 +637,7 @@ plotCrossDrug <- function( ###################### CROSS DRUG Testing ############################# heatmap_df <- cross_drug |> - dplyr::filter(!is.na(drug), !is.na(test_drug)) |> + dplyr::filter(!is.na(drug), !is.na(test_drug)) |> # dplyr::filter(test_drug %in% (cross_drug |> dplyr::pull(drug))) |> dplyr::group_by(drug, test_drug) |> dplyr::summarise(median_mcc = median(mcc, na.rm = TRUE), .groups = "drop") @@ -716,10 +716,10 @@ plotCrossDrug <- function( # ) class_colors <- stats::setNames( - colorRampPalette(RColorBrewer::brewer.pal(9, "Pastel1"))(length(classes)), - classes -) - + colorRampPalette(RColorBrewer::brewer.pal(9, "Pastel1"))(length(classes)), + classes + ) + heat_colors <- colorRampPalette(RColorBrewer::brewer.pal(11, "RdBu"))(100) max_val <- max(abs(mat), na.rm = TRUE) @@ -735,7 +735,7 @@ plotCrossDrug <- function( ha_col <- ComplexHeatmap::HeatmapAnnotation( class_abbr = annotation_col$class_abbr, col = list(class_abbr = class_colors), - show_annotation_name = FALSE, na_col = "grey3", + show_annotation_name = FALSE, na_col = "grey3", annotation_legend_param = list(labels_gp = grid::gpar(fontsize = 14)) ) @@ -755,12 +755,8 @@ plotCrossDrug <- function( column_order = col_order, left_annotation = ha_row, top_annotation = ha_col, - - - width = grid::unit(ncol(mat), "in"), + width = grid::unit(ncol(mat), "in"), height = grid::unit(nrow(mat), "in"), - - show_row_names = TRUE, show_column_names = TRUE, column_title = "tested on", @@ -776,7 +772,6 @@ plotCrossDrug <- function( # legends show_heatmap_legend = TRUE, - use_raster = FALSE ) @@ -819,22 +814,22 @@ plotStratifiedPerf <- function(year_or_country = "year", ) )) # if (year_or_country == "year") { - all <- perf |> - # dplyr::rename("train_year" = "strat_value") |> - dplyr::mutate(strat_value_test = strat_value) |> + all <- perf |> + # dplyr::rename("train_year" = "strat_value") |> + dplyr::mutate(strat_value_test = strat_value) |> + dplyr::select( + drug_label, drug_or_class, + strat_value, strat_value_test, feature_type, feature_subtype, mcc + ) |> + dplyr::bind_rows(cross_test |> dplyr::select( drug_label, drug_or_class, - strat_value, strat_value_test, feature_type, feature_subtype, mcc - ) |> - dplyr::bind_rows(cross_test |> - dplyr::select( - drug_label, drug_or_class, - strat_value, strat_value_test, feature_type, - feature_subtype, mcc - )) |> - dplyr::mutate(category = dplyr::if_else( - strat_value == strat_value_test, "same", "different" - )) + strat_value, strat_value_test, feature_type, + feature_subtype, mcc + )) |> + dplyr::mutate(category = dplyr::if_else( + strat_value == strat_value_test, "same", "different" + )) # } else { # all <- perf |> # dplyr::rename("train_country" = "strat_value") |> @@ -868,9 +863,9 @@ plotStratifiedPerf <- function(year_or_country = "year", # } fill_vals <- c( - "same" = "#b3cde3", - "different" = "#fbb4ae" - ) + "same" = "#b3cde3", + "different" = "#fbb4ae" + ) plot <- ggplot2::ggplot( all |> From 2cb0e6604cc03f7df6a39dfc6fa1c2bd0bd2dd5c Mon Sep 17 00:00:00 2001 From: Abhirupa Ghosh <100681585+AbhirupaGhosh@users.noreply.github.com> Date: Mon, 6 Jul 2026 14:12:41 -0600 Subject: [PATCH 10/17] Update plot_ml.R --- R/plot_ml.R | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/R/plot_ml.R b/R/plot_ml.R index 8a89296..694e62b 100644 --- a/R/plot_ml.R +++ b/R/plot_ml.R @@ -736,7 +736,7 @@ plotCrossDrug <- function( class_abbr = annotation_col$class_abbr, col = list(class_abbr = class_colors), show_annotation_name = FALSE, na_col = "grey3", - annotation_legend_param = list(labels_gp = grid::gpar(fontsize = 14)) + annotation_legend_param = list(labels_gp = grid::gpar(fontsize = 10)) ) # ---- Color function (instead of breaks + palette) ---- @@ -755,8 +755,8 @@ plotCrossDrug <- function( column_order = col_order, left_annotation = ha_row, top_annotation = ha_col, - width = grid::unit(ncol(mat), "in"), - height = grid::unit(nrow(mat), "in"), + width = grid::unit(ncol(mat) * 1.5, "cm"), + height = grid::unit(nrow(mat)* 1.5, "cm"), show_row_names = TRUE, show_column_names = TRUE, column_title = "tested on", From 6d150d9d18791e11037b4c6312584c19c3bec796 Mon Sep 17 00:00:00 2001 From: AbhirupaGhosh Date: Mon, 6 Jul 2026 20:14:52 +0000 Subject: [PATCH 11/17] Style code (GHA) --- R/plot_ml.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/plot_ml.R b/R/plot_ml.R index 694e62b..551df87 100644 --- a/R/plot_ml.R +++ b/R/plot_ml.R @@ -756,7 +756,7 @@ plotCrossDrug <- function( left_annotation = ha_row, top_annotation = ha_col, width = grid::unit(ncol(mat) * 1.5, "cm"), - height = grid::unit(nrow(mat)* 1.5, "cm"), + height = grid::unit(nrow(mat) * 1.5, "cm"), show_row_names = TRUE, show_column_names = TRUE, column_title = "tested on", From 1e42c53c883e78720fc699bae48a68947f552bf1 Mon Sep 17 00:00:00 2001 From: Abhirupa Ghosh <100681585+AbhirupaGhosh@users.noreply.github.com> Date: Mon, 6 Jul 2026 18:27:18 -0600 Subject: [PATCH 12/17] Update plot_ml.R --- R/plot_ml.R | 36 +++++++++++++++++++++++++----------- 1 file changed, 25 insertions(+), 11 deletions(-) diff --git a/R/plot_ml.R b/R/plot_ml.R index 551df87..5be662e 100644 --- a/R/plot_ml.R +++ b/R/plot_ml.R @@ -489,15 +489,29 @@ plotDrugPerf <- function(metadata_path = ".", performance_path = ".") { ) + ggplot2::geom_tile(color = "grey90", width = 0.9) + ggplot2::scale_fill_gradientn( - colors = c( - "#C4B8A8", # low - "#FAFAF7", # around 0 - "#5F84C9", # medium/high (~0.7–0.9) - "#0F2A5A" # very dark for ~1 - ), - values = scales::rescale(c(-1, 0, 0.85, 1)), - name = "Best MCC" - ) + + colors = c( + "#C4B8A8", + "#FAFAF7", + "#5F84C9", + "#0F2A5A" + ), + values = scales::rescale(c(-1, 0, 0.85, 1)), + breaks = c(-1, -0.5, 0, 0.5, 1), + labels = scales::label_number(accuracy = 0.01), + name = "Median MCC" +) + +# ggplot2::scale_fill_gradientn( +# colors = c( +# "white", +# "#BDD7E7", +# "#6BAED6", +# "#2171B5", +# "#08306B" +# ), +# values = scales::rescale(c(0, 0.4, 0.6, 0.8, 1)), +# labels = scales::label_number(accuracy = 0.1), +# name = "Median MCC" +# ) + ggplot2::labs(x = "Feature type") + ggplot2::theme_minimal(base_size = 12) + ggplot2::theme( @@ -547,10 +561,10 @@ plotDrugPerf <- function(metadata_path = ".", performance_path = ".") { colour = "grey70" ) + ggplot2::geom_point( - position = position_jitter(height = 0.1), + position = ggplot2::position_jitter(height = 0.1), size = 1.5, alpha = 0.8, - aes(color = feature_type) + ggplot2::aes(color = feature_type) ) + ggplot2::scale_color_manual(values = feat_pal, name = "Feature type") + ggplot2::stat_summary( From 4601ada2ac21fde4c1e30e669bd1c6f57d384500 Mon Sep 17 00:00:00 2001 From: AbhirupaGhosh Date: Tue, 7 Jul 2026 00:29:30 +0000 Subject: [PATCH 13/17] Style code (GHA) --- R/plot_ml.R | 46 +++++++++++++++++++++++----------------------- 1 file changed, 23 insertions(+), 23 deletions(-) diff --git a/R/plot_ml.R b/R/plot_ml.R index 5be662e..a2fe496 100644 --- a/R/plot_ml.R +++ b/R/plot_ml.R @@ -489,29 +489,29 @@ plotDrugPerf <- function(metadata_path = ".", performance_path = ".") { ) + ggplot2::geom_tile(color = "grey90", width = 0.9) + ggplot2::scale_fill_gradientn( - colors = c( - "#C4B8A8", - "#FAFAF7", - "#5F84C9", - "#0F2A5A" - ), - values = scales::rescale(c(-1, 0, 0.85, 1)), - breaks = c(-1, -0.5, 0, 0.5, 1), - labels = scales::label_number(accuracy = 0.01), - name = "Median MCC" -) + -# ggplot2::scale_fill_gradientn( -# colors = c( -# "white", -# "#BDD7E7", -# "#6BAED6", -# "#2171B5", -# "#08306B" -# ), -# values = scales::rescale(c(0, 0.4, 0.6, 0.8, 1)), -# labels = scales::label_number(accuracy = 0.1), -# name = "Median MCC" -# ) + + colors = c( + "#C4B8A8", + "#FAFAF7", + "#5F84C9", + "#0F2A5A" + ), + values = scales::rescale(c(-1, 0, 0.85, 1)), + breaks = c(-1, -0.5, 0, 0.5, 1), + labels = scales::label_number(accuracy = 0.01), + name = "Median MCC" + ) + + # ggplot2::scale_fill_gradientn( + # colors = c( + # "white", + # "#BDD7E7", + # "#6BAED6", + # "#2171B5", + # "#08306B" + # ), + # values = scales::rescale(c(0, 0.4, 0.6, 0.8, 1)), + # labels = scales::label_number(accuracy = 0.1), + # name = "Median MCC" + # ) + ggplot2::labs(x = "Feature type") + ggplot2::theme_minimal(base_size = 12) + ggplot2::theme( From 34c691a0041bad377866ac023aa19223a89f124c Mon Sep 17 00:00:00 2001 From: Abhirupa Ghosh <100681585+AbhirupaGhosh@users.noreply.github.com> Date: Tue, 7 Jul 2026 09:54:57 -0600 Subject: [PATCH 14/17] resolve comments and suggestions --- R/plot_ml.R | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/R/plot_ml.R b/R/plot_ml.R index a2fe496..ac5b1c8 100644 --- a/R/plot_ml.R +++ b/R/plot_ml.R @@ -207,7 +207,7 @@ plotTopFeatsVI <- function(topfeat, n_top_feats = 10) { #' @param non_shuffled_label_results Output of `runMLPipeline(shuffle_labels = FALSE)` #' @param shuffled_label_results Output of `runMLPipeline(shuffle_labels = TRUE)` #' -#' @return A base R barplot comparing balanced accuracy across models. +#' @return A barplot comparing balanced accuracy across models. #' #' @examples #' non_shuffled <- list( @@ -598,7 +598,7 @@ plotDrugPerf <- function(metadata_path = ".", performance_path = ".") { final_plot } -#' Plot cross-drug generalization heatmap +#' Plot cross-drug prediction as a heatmap #' #' Creates a heatmap showing cross-drug model performance (MCC), where models #' trained on one drug are evaluated on another. @@ -649,7 +649,7 @@ plotCrossDrug <- function( performance <- drug_performance_path metadata <- arrow::read_parquet(file.path(metadata_path, "metadata.parquet")) - ###################### CROSS DRUG Testing ############################# + heatmap_df <- cross_drug |> dplyr::filter(!is.na(drug), !is.na(test_drug)) |> # dplyr::filter(test_drug %in% (cross_drug |> dplyr::pull(drug))) |> @@ -897,9 +897,9 @@ plotStratifiedPerf <- function(year_or_country = "year", ggplot2::theme_minimal(base_size = 14) + ggplot2::labs( title = if (year_or_country == "year") { - "Temporal performance by drug" + "Performance of AMR drug models with temporal holdouts" } else { - "Geographical performance by drug" + "Performance of AMR drug models with geographical holdouts" }, x = "MCC", y = "Drug", @@ -1222,7 +1222,7 @@ plotShuffleVsReal <- function(metadata_path = ".", performance_path = ".") { #' } plotTopClusters <- function(top_feat_path = ".", cluster_feature_path = ".", protein_names_path = ".", top_n = 10) { - ################### Top features ######################### + top_feat <- arrow::read_parquet(file.path(top_feat_path, "all_top_features.parquet")) cluster_feature <- arrow::read_parquet(file.path(cluster_feature_path, "cluster_feature.parquet")) From aa2abd26df7419ae16af8135c2cb9dd935404a2f Mon Sep 17 00:00:00 2001 From: AbhirupaGhosh Date: Tue, 7 Jul 2026 15:56:38 +0000 Subject: [PATCH 15/17] Style code (GHA) --- R/plot_ml.R | 2 -- 1 file changed, 2 deletions(-) diff --git a/R/plot_ml.R b/R/plot_ml.R index ac5b1c8..6cf4910 100644 --- a/R/plot_ml.R +++ b/R/plot_ml.R @@ -1222,8 +1222,6 @@ plotShuffleVsReal <- function(metadata_path = ".", performance_path = ".") { #' } plotTopClusters <- function(top_feat_path = ".", cluster_feature_path = ".", protein_names_path = ".", top_n = 10) { - - top_feat <- arrow::read_parquet(file.path(top_feat_path, "all_top_features.parquet")) cluster_feature <- arrow::read_parquet(file.path(cluster_feature_path, "cluster_feature.parquet")) protein_names <- arrow::read_parquet(file.path(protein_names_path, "protein_names.parquet")) From 769e69f73bfde9bf768e634b4c224837691aa086 Mon Sep 17 00:00:00 2001 From: Alexander McKim Date: Thu, 9 Jul 2026 15:21:33 -0600 Subject: [PATCH 16/17] CI fixes and genome_drug metadata --- R/ife_ml.R | 2 +- R/plot_ml.R | 8 ++++---- R/run_ML.R | 4 ++-- R/run_ml_pipeline.R | 6 ++---- man/plotBaselineComparison.Rd | 2 +- man/plotCrossDrug.Rd | 10 +++++----- man/plotDrugPerf.Rd | 2 +- man/plotShuffleVsReal.Rd | 2 +- man/runIFE.Rd | 2 +- 9 files changed, 18 insertions(+), 20 deletions(-) diff --git a/R/ife_ml.R b/R/ife_ml.R index d4fae1e..f94a2a3 100644 --- a/R/ife_ml.R +++ b/R/ife_ml.R @@ -82,7 +82,7 @@ removeTopFeats <- function(ml_input_tibble, top_feat_tibble) { #' @export runIFE <- function( ml_input_tibble, by_num = TRUE, by_vi = FALSE, - percent_removal_vec = 10 * 1:9, mix_vec = 0, return_feats = FALSE, + percent_removal_vec = 10 * seq_len(9), mix_vec = 0, return_feats = FALSE, verbose = TRUE ) { .checkArgByNum(by_num) diff --git a/R/plot_ml.R b/R/plot_ml.R index 6cf4910..d0203e1 100644 --- a/R/plot_ml.R +++ b/R/plot_ml.R @@ -382,7 +382,7 @@ plotDrugDist <- function(metadata_path = ".") { ##################### phenotype distribution (drugs) ######################### drug_dist <- metadata |> dplyr::distinct( - genome.genome_id, + genome_drug.genome_id, genome_drug.antibiotic, drug_abbr, genome_drug.resistant_phenotype @@ -465,7 +465,7 @@ plotDrugPerf <- function(metadata_path = ".", performance_path = ".") { performance <- performance_path plot_df <- metadata |> - dplyr::distinct(genome.genome_id, genome_drug.antibiotic, drug_abbr) |> + dplyr::distinct(genome_drug.genome_id, genome_drug.antibiotic, drug_abbr) |> dplyr::count(genome_drug.antibiotic, drug_abbr, name = "total") ######################## drug performances ################################# @@ -604,7 +604,7 @@ plotDrugPerf <- function(metadata_path = ".", performance_path = ".") { #' trained on one drug are evaluated on another. #' #' @param cross_test_performance_path A cross-drug performance tibble (with -#' `drug_or_class`, `test_drug`, and `mcc` columns), or a directory path +#' `drug`, `test_drug`, and `mcc` columns), or a directory path #' containing `cross_drug_perf.parquet`. #' @param drug_performance_path A performance tibble (with `drug_label`, #' `drug_or_class`, and `mcc` columns), or a directory path containing @@ -616,7 +616,7 @@ plotDrugPerf <- function(metadata_path = ".", performance_path = ".") { #' #' @examples #' cross_drug <- tibble::tibble( -#' drug_or_class = c("AMP", "AMP", "CIP", "CIP", "NAL", "NAL"), +#' drug = c("AMP", "AMP", "CIP", "CIP", "NAL", "NAL"), #' test_drug = c("CIP", "NAL", "AMP", "NAL", "AMP", "CIP"), #' mcc = c(0.3, 0.2, 0.4, 0.25, 0.15, 0.35) #' ) diff --git a/R/run_ML.R b/R/run_ML.R index a28d821..bd0fe0a 100644 --- a/R/run_ML.R +++ b/R/run_ML.R @@ -334,11 +334,11 @@ createMLinputList <- function(path, # Case B: unstratified -> prefix is first two tokens if (x[2] == "drug" && x[3] != "class") { # Case A: Cje_drug_X - return(paste(x[1:2], collapse = "_")) + return(paste(x[seq_len(2)], collapse = "_")) } if (x[2] == "drug" && x[3] == "class") { # Case A: Cje_drug_X - return(paste(x[1:3], collapse = "_")) + return(paste(x[seq_len(3)], collapse = "_")) } }) ) diff --git a/R/run_ml_pipeline.R b/R/run_ml_pipeline.R index cab507d..db0fc64 100644 --- a/R/run_ml_pipeline.R +++ b/R/run_ml_pipeline.R @@ -144,10 +144,8 @@ runMLPipeline <- function( if (verbose) { message( - paste( - "`ml_input_tibble` has", num_obs_ml_input_tibble, "observations of", - getNumFeat(ml_input_tibble), "features." - ) + "`ml_input_tibble` has ", num_obs_ml_input_tibble, " observations of ", + getNumFeat(ml_input_tibble), " features." ) } diff --git a/man/plotBaselineComparison.Rd b/man/plotBaselineComparison.Rd index 8eab040..8d6403d 100644 --- a/man/plotBaselineComparison.Rd +++ b/man/plotBaselineComparison.Rd @@ -12,7 +12,7 @@ plotBaselineComparison(non_shuffled_label_results, shuffled_label_results) \item{shuffled_label_results}{Output of \code{runMLPipeline(shuffle_labels = TRUE)}} } \value{ -A base R barplot comparing balanced accuracy across models. +A barplot comparing balanced accuracy across models. } \description{ Produces a bar plot comparing balanced accuracy for each antibiotic using diff --git a/man/plotCrossDrug.Rd b/man/plotCrossDrug.Rd index 1415f1e..4c9163d 100644 --- a/man/plotCrossDrug.Rd +++ b/man/plotCrossDrug.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/plot_ml.R \name{plotCrossDrug} \alias{plotCrossDrug} -\title{Plot cross-drug generalization heatmap} +\title{Plot cross-drug prediction as a heatmap} \usage{ plotCrossDrug( cross_test_performance_path = ".", @@ -12,12 +12,12 @@ plotCrossDrug( } \arguments{ \item{cross_test_performance_path}{A cross-drug performance tibble (with -\code{drug_or_class}, \code{tested_on}, and \code{mcc} columns), or a directory path +\code{drug}, \code{test_drug}, and \code{mcc} columns), or a directory path containing \code{cross_drug_perf.parquet}.} \item{drug_performance_path}{A performance tibble (with \code{drug_label}, \code{drug_or_class}, and \code{mcc} columns), or a directory path containing -\code{all_performance.parquet}.} +\code{all_perf.parquet}.} \item{metadata_path}{Character. Path to a directory containing \code{metadata.parquet}.} } @@ -30,8 +30,8 @@ trained on one drug are evaluated on another. } \examples{ cross_drug <- tibble::tibble( - drug_or_class = c("AMP", "AMP", "CIP", "CIP", "NAL", "NAL"), - tested_on = c("CIP", "NAL", "AMP", "NAL", "AMP", "CIP"), + drug = c("AMP", "AMP", "CIP", "CIP", "NAL", "NAL"), + test_drug = c("CIP", "NAL", "AMP", "NAL", "AMP", "CIP"), mcc = c(0.3, 0.2, 0.4, 0.25, 0.15, 0.35) ) performance <- tibble::tibble( diff --git a/man/plotDrugPerf.Rd b/man/plotDrugPerf.Rd index 6b83cdb..8fd65df 100644 --- a/man/plotDrugPerf.Rd +++ b/man/plotDrugPerf.Rd @@ -11,7 +11,7 @@ plotDrugPerf(metadata_path = ".", performance_path = ".") \item{performance_path}{A performance tibble (with \code{drug_label}, \code{shuffled}, \code{drug_or_class}, \code{feature_type}, \code{feature_subtype}, and \code{mcc} columns), or a -directory path containing \code{all_performance.parquet}.} +directory path containing \code{all_perf.parquet}.} } \value{ A patchwork ggplot object combining multiple panels. diff --git a/man/plotShuffleVsReal.Rd b/man/plotShuffleVsReal.Rd index d000e89..3ed40f0 100644 --- a/man/plotShuffleVsReal.Rd +++ b/man/plotShuffleVsReal.Rd @@ -11,7 +11,7 @@ plotShuffleVsReal(metadata_path = ".", performance_path = ".") \item{performance_path}{A performance tibble (with \code{feature_type}, \code{feature_subtype}, \code{mcc}, and \code{shuffled} columns), or a directory path -containing \code{all_performance.parquet}.} +containing \code{all_perf.parquet}.} } \value{ A ggplot object. diff --git a/man/runIFE.Rd b/man/runIFE.Rd index f817ba9..31afb6b 100644 --- a/man/runIFE.Rd +++ b/man/runIFE.Rd @@ -10,7 +10,7 @@ runIFE( ml_input_tibble, by_num = TRUE, by_vi = FALSE, - percent_removal_vec = 10 * 1:9, + percent_removal_vec = 10 * seq_len(9), mix_vec = 0, return_feats = FALSE, verbose = TRUE From 98f416ce9810c749a69a77ad0fe55a271dd6fea1 Mon Sep 17 00:00:00 2001 From: Abhirupa Ghosh <100681585+AbhirupaGhosh@users.noreply.github.com> Date: Thu, 9 Jul 2026 15:54:25 -0600 Subject: [PATCH 17/17] Color changes in the plots --- R/plot_ml.R | 25 ++++++++++++------------- 1 file changed, 12 insertions(+), 13 deletions(-) diff --git a/R/plot_ml.R b/R/plot_ml.R index d0203e1..3fde725 100644 --- a/R/plot_ml.R +++ b/R/plot_ml.R @@ -382,7 +382,7 @@ plotDrugDist <- function(metadata_path = ".") { ##################### phenotype distribution (drugs) ######################### drug_dist <- metadata |> dplyr::distinct( - genome_drug.genome_id, + genome.genome_id, genome_drug.antibiotic, drug_abbr, genome_drug.resistant_phenotype @@ -465,7 +465,7 @@ plotDrugPerf <- function(metadata_path = ".", performance_path = ".") { performance <- performance_path plot_df <- metadata |> - dplyr::distinct(genome_drug.genome_id, genome_drug.antibiotic, drug_abbr) |> + dplyr::distinct(genome.genome_id, genome_drug.antibiotic, drug_abbr) |> dplyr::count(genome_drug.antibiotic, drug_abbr, name = "total") ######################## drug performances ################################# @@ -489,17 +489,16 @@ plotDrugPerf <- function(metadata_path = ".", performance_path = ".") { ) + ggplot2::geom_tile(color = "grey90", width = 0.9) + ggplot2::scale_fill_gradientn( - colors = c( - "#C4B8A8", - "#FAFAF7", - "#5F84C9", - "#0F2A5A" - ), - values = scales::rescale(c(-1, 0, 0.85, 1)), - breaks = c(-1, -0.5, 0, 0.5, 1), - labels = scales::label_number(accuracy = 0.01), - name = "Median MCC" - ) + + colors = c( + "#EAF2FF", # very light blue + "#5F84C9", # medium blue + "#0F2A5A" # dark blue + ), + values = scales::rescale(c(-1, 0, 1)), + breaks = c(-1, -0.5, 0, 0.5, 1), + labels = scales::label_number(accuracy = 0.01), + name = "Median MCC" +) + # ggplot2::scale_fill_gradientn( # colors = c( # "white",