amRml predicts antimicrobial resistance (AMR) from
bacterial genomic features. It consumes a DuckDB produced by
amRdata and produces ML matrices, tuned logistic regression
models, per-genome predictions, feature importances, and Fisher’s exact
tests as a non-ML baseline.
This vignette uses the Shigella flexneri (Sfl)
DuckDB bundled in inst/extdata.
fixture <- system.file("extdata", "Sfl_parquet.duckdb", package = "amRml")
out_dir <- file.path(tempdir(), "amRml_vignette")
dir.create(out_dir, showWarnings = FALSE, recursive = TRUE)generateMLInputs() reads the bug-level DuckDB (metadata
+ feature parquets) and writes one long-format sparse parquet per drug ×
feature × encoding combination into out_path/matrix/. With
stratify_by, it additionally writes year- or
country-stratified matrices into matrix_year/ or
matrix_country/.
generateMLInputs(
parquet_duckdb_path = fixture,
out_path = out_dir,
n_fold = 5,
split = c(1, 0),
min_n = 25,
verbosity = "minimal"
)
list.files(file.path(out_dir, "matrix"))[1:5]For a classical train/validation/test split instead of
cross-validation, set n_fold = NULL and pass a length-2
split of c(train_prop, val_prop); the test
proportion is the remainder.
generateMLInputs(
parquet_duckdb_path = fixture,
out_path = out_dir,
n_fold = NULL,
split = c(0.7, 0.15),
verbosity = "debug"
)The generated parquets are long-format sparse tibbles with these columns:
| Column | Description |
|---|---|
genome_id |
Unique identifier for each isolate |
feature_id |
Feature name (gene, protein, domain, or struct) |
value |
Binary presence/absence (0/1) or count |
genome_drug.resistant_phenotype |
"Resistant" or "Susceptible"
|
loadMLInputTibble() converts one of them to wide format
(one row per genome, one column per feature) ready for ML.
matrix_path <- file.path(
out_dir, "matrix", "Sfl_drug_AMP_genes_binary_sparse.parquet"
)
ml_tibble <- loadMLInputTibble(matrix_path)
n_features <- getNumFeat(ml_tibble)
target_var <- .getTargetVarName(ml_tibble)
c(n_features = n_features, target_var = target_var)runMLPipeline() runs the train/tune/fit/predict pipeline
on a single matrix in memory. Use it to iterate on one drug-feature
combo before scaling to all of them.
results <- runMLPipeline(
ml_input_tibble = ml_tibble,
model = "LR",
split = c(1, 0),
n_fold = 2,
n_top_feats = 20,
penalty_vec = 10^c(-3, -1),
mix_vec = c(0, 0.5, 1),
select_best_metric = "mcc",
return_fit = TRUE,
return_pred = TRUE,
verbose = FALSE
)
results$performance_tibble
head(results$top_feat_tibble)runMLPipeline() returns a named list:
performance_tibble — one row of model
performance metrics:
| Column | Description |
|---|---|
num_obs |
Number of observations |
res_prop |
Proportion of resistant samples |
n_feat |
Number of features |
model |
Model type ("LR") |
train_prop, val_prop
|
Train/validation split proportions |
fit_penalty, fit_mixture
|
Fitted hyperparameters |
mcc, nmcc, f1,
bal_acc, log2_apop
|
Performance metrics |
run_time_sec |
Runtime in seconds |
top_feat_tibble — ranked feature
importance:
| Column | Description |
|---|---|
Variable |
Feature name |
Importance |
Variable importance score |
Sign |
Direction of effect (POS = associated with resistance,
NEG = with susceptibility) |
Optional outputs (when return_* =
TRUE):
tune_res — tuning results from grid searchfit — the fitted workflow objectpred — predictions with .pred_class,
.pred_Resistant, .pred_Susceptible
The builders below are what runMLPipeline() chains
internally. Call them directly when you need control over any individual
step.
data_split <- splitMLInputTibble(ml_tibble, split = c(0.6, 0.2), seed = 123)
train_data <- rsample::training(data_split)
test_data <- rsample::testing(data_split)
recipe <- buildRecipe(train_data, use_pca = FALSE)
lr_mod <- buildLRModel(multi_class = FALSE)
wflow <- buildWflow(lr_mod, recipe)
grid <- buildTuningGrid(
model = "LR",
penalty_vec = 10^c(-3, -1),
mix_vec = c(0, 0.5, 1)
)
tune_res <- tuneGrid(wflow, data_split, grid, n_fold = 2)
best_wflow <- selectBestModel(tune_res, wflow, select_best_metric = "mcc")
fit <- fitBestModel(best_wflow, train_data)
preds <- predictML(fit, test_data)calculateEvalMets() returns all of nMCC, F1, balanced
accuracy, AUPRC, log2(AUPRC/prior), sensitivity, and specificity from a
tibble of predictions + truth.
calculateEvalMets(preds)
getConfusionMatrix(preds)extractTopFeats() ranks features by absolute coefficient
(for LR). Use n_top_feats for a fixed count or
prop_vi_top_feats for a percentile range.
top_features <- extractTopFeats(fit, n_top_feats = 20)
head(top_features)
plotPRC(results$pred)
plotTopFeatsVI(results$top_feat_tibble, n_top_feats = 10)For a baseline comparison against random labels, fit a shuffled-label pipeline and compare:
shuffled <- runMLPipeline(
ml_input_tibble = ml_tibble,
model = "LR",
split = c(1, 0),
n_fold = 2,
shuffle_labels = TRUE,
return_pred = TRUE
)
plotBaselineComparison(
non_shuffled_label_results = results$performance_tibble,
shuffled_label_results = shuffled$performance_tibble
)runIFE() retrains the model after iteratively removing
top-ranked features, helping identify the minimal predictive subset. It
runs the pipeline once per percentile in
percent_removal_vec
ife_results <- runIFE(
ml_tibble,
by_num = TRUE,
by_vi = FALSE,
percent_removal_vec = 10 * 1:9,
mix_vec = 0,
return_feats = TRUE,
verbose = FALSE
)
ife_results$ife_performance_tibble
ife_results$feats_removedremoveTopFeats() strips a given set of features from a
matrix tibble if you want to do this manually:
trimmed <- removeTopFeats(ml_tibble, head(top_features, 5))
ncol(ml_tibble) - ncol(trimmed)runFishers() runs a Fisher’s exact test of feature
presence vs. phenotype for each feature, applies Benjamini–Hochberg
correction, and computes per-class frequencies.
fisher_results <- runFishers(
matrix_path = matrix_path,
Q = 0.05,
alternative = "two.sided",
susceptible_label = "Susceptible",
resistant_label = "Resistant"
)
head(fisher_results)
plotFishers(fisher_results, alpha = 0.05, label_top_n = 5)runMLmodels
runMLmodels() trains a model on every matrix produced by
generateMLInputs() and writes performance TSVs into
out_path/ML_performance/, predictions into
ML_pred/, and top features into
ML_top_features/. Takes over an hour on Sfl with default
settings.
runMLmodels(
path = out_dir,
stratify_by = NULL,
LOO = FALSE,
cross_test = FALSE,
threads = max(1L, parallel::detectCores() - 1L),
split = c(1, 0),
n_fold = 5,
verbose = TRUE,
return_pred = TRUE,
use_saved_split = TRUE
)runModelingPipeline
For the full pipeline from a DuckDB to all outputs in one call:
runModelingPipeline(
parquet_duckdb_path = fixture,
threads = max(1L, parallel::detectCores() - 1L),
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
)
sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.4 LTS
#>
#> Matrix products: default
#> BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so; LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8
#> [4] LC_COLLATE=C.UTF-8 LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8
#> [7] LC_PAPER=C.UTF-8 LC_NAME=C LC_ADDRESS=C
#> [10] LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: UTC
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] amRml_0.99.0 BiocStyle_2.40.0
#>
#> loaded via a namespace (and not attached):
#> [1] DBI_1.3.0 rlang_1.3.0 magrittr_2.0.5
#> [4] furrr_0.4.0 otel_0.2.0 compiler_4.6.1
#> [7] systemfonts_1.3.2 vctrs_0.7.3 stringr_1.6.0
#> [10] tune_2.1.0 pkgconfig_2.0.3 shape_1.4.6.1
#> [13] fastmap_1.2.0 rmarkdown_2.31 prodlim_2026.03.11
#> [16] tzdb_0.5.0 ragg_1.5.2 purrr_1.2.2
#> [19] bit_4.6.0 xfun_0.60 glmnet_5.0
#> [22] cachem_1.1.0 jsonlite_2.0.0 recipes_1.3.3
#> [25] vip_0.4.6 parallel_4.6.1 R6_2.6.1
#> [28] bslib_0.11.0 stringi_1.8.7 rsample_1.3.2
#> [31] RColorBrewer_1.1-3 parallelly_1.48.0 rpart_4.1.27
#> [34] lubridate_1.9.5 jquerylib_0.1.4 Rcpp_1.1.2
#> [37] bookdown_0.47 assertthat_0.2.1 dials_1.4.4
#> [40] iterators_1.0.14 knitr_1.51 future.apply_1.20.2
#> [43] readr_2.2.0 Matrix_1.7-5 splines_4.6.1
#> [46] nnet_7.3-20 timechange_0.4.0 tidyselect_1.2.1
#> [49] yaml_2.3.12 timeDate_4052.112 codetools_0.2-20
#> [52] listenv_1.0.0 lattice_0.22-9 tibble_3.3.1
#> [55] withr_3.0.3 S7_0.2.2 evaluate_1.0.5
#> [58] future_1.75.0 desc_1.4.3 survival_3.8-6
#> [61] pillar_1.11.1 BiocManager_1.30.27 foreach_1.5.2
#> [64] generics_0.1.4 hms_1.1.4 ggplot2_4.0.3
#> [67] scales_1.4.0 globals_0.19.1 class_7.3-23
#> [70] glue_1.8.1 tools_4.6.1 data.table_1.18.4
#> [73] gower_1.0.2 fs_2.1.0 grid_4.6.1
#> [76] yardstick_1.4.0 tidyr_1.3.2 workflowsets_1.1.1
#> [79] ipred_0.9-15 duckdb_1.5.4.3 cli_3.6.6
#> [82] DiceDesign_1.10 textshaping_1.0.5 workflows_1.3.0
#> [85] parsnip_1.6.0 lava_1.9.2 arrow_25.0.0
#> [88] dplyr_1.2.1 gtable_0.3.6 sass_0.4.10
#> [91] digest_0.6.39 ggrepel_0.9.8 farver_2.1.2
#> [94] htmltools_0.5.9 pkgdown_2.2.1 lifecycle_1.0.5
#> [97] hardhat_1.4.3 bit64_4.8.2 MASS_7.3-65