From 03d28b2943dca5f1a419b754bcd41fa1ce5ae4ba Mon Sep 17 00:00:00 2001 From: Nick Waters Date: Tue, 16 Sep 2025 13:20:07 -0400 Subject: [PATCH 1/3] WIP handle all-unclassified mpa samples --- R/utilities.R | 15 +++++++++++++-- tests/testthat/test-utilities.R | 9 +++++++++ 2 files changed, 22 insertions(+), 2 deletions(-) diff --git a/R/utilities.R b/R/utilities.R index 41f2972..091e319 100644 --- a/R/utilities.R +++ b/R/utilities.R @@ -180,11 +180,22 @@ vdb_make_phylo <- function(metadata, sampleid_col = "sampleid", skip_seqs = TRUE get_metaphlan_analyses <- function(con, analysis_ids, schema="public") { - raw_results <- get_subset_pg_df("mgx_metaphlan", "ia_id", analysis_ids, schema=schema) %>% + raw_results <- get_subset_pg_df("mgx_metaphlan", "ia_id", analysis_ids, schema=schema) + # we have to deal with samples having no classified -- in those cases metaphlan reports + # estimated_number_of_reads_from_the_clade as 0, but to make it work with the rest of our code we modify that + # to nreads_input (eg all the input reads are from the unclassified clade) + if (sum(raw_results[raw_results$clade_name == "unclassified", "estimated_number_of_reads_from_the_clade"]) > 0){ + print(raw_results[raw_results$clade_name== "unclassified", "ia_id"]) + stop("some of the above analyses have metaphlan failures where estimated_number of reads from the clade is not zero as expected; please alert of the vdbR delevopers that the metaphlan output has changed") + } + # TODO: parametarize this to get sgb level results + results_with_zeros_handled <- raw_results %>% + dplyr::mutate(estimated_number_of_reads_from_the_clade = ifelse(clade_name == "unclassified", nreads_input, estimated_number_of_reads_from_the_clade)) %>% + dplyr::mutate(clade_name = ifelse(clade_name == "unclassified", "UNCLASSIFIED", clade_name)) %>% dplyr::filter(grepl("UNCLASSIFIED", clade_name) | grepl(".*\\|s__.*", clade_name)) %>% dplyr::filter(!grepl(".*t__.*", clade_name)) %>% dplyr::mutate(clade_name = ifelse(clade_name == "UNCLASSIFIED", "k__UNCLASSIFIED", clade_name)) - wide_results <- raw_results %>% + wide_results <- results_with_zeros_handled %>% dplyr::select(ia_id, clade_name, relative_abundance) %>% tidyr::pivot_wider(names_from = ia_id, values_from = relative_abundance, values_fill = 0) md <- raw_results %>% diff --git a/tests/testthat/test-utilities.R b/tests/testthat/test-utilities.R index 6a0962e..2e21547 100644 --- a/tests/testthat/test-utilities.R +++ b/tests/testthat/test-utilities.R @@ -200,3 +200,12 @@ test_that("test that sample order is retained", { }) +test_that("metaphlan failures are handled", { + # 32298 is a sample with no mapped reads + skip_if(Sys.getenv("GITHUB_ACTIONS") != "") + connect_database(bundled = FALSE) + analysis_ids = c(32298, 32299) + res = get_metaphlan_analyses(con = psql_con, analysis_ids) + # assert (a) that the column exists on output and (b) that we have full 100% abundance for that sample + expect_equal(sum(res[[1]]$`32298`), 100) +}) From dd6d378ccd2e891c20f6bdff9314e2ed9109de82 Mon Sep 17 00:00:00 2001 From: Nick Waters Date: Tue, 16 Sep 2025 14:02:40 -0400 Subject: [PATCH 2/3] default to failing when considering no-hit mpa samples --- R/utilities.R | 36 +++++++++++++++++++++++++-------- tests/testthat/test-utilities.R | 13 ++++++++++-- 2 files changed, 39 insertions(+), 10 deletions(-) diff --git a/R/utilities.R b/R/utilities.R index 091e319..e01db58 100644 --- a/R/utilities.R +++ b/R/utilities.R @@ -179,23 +179,43 @@ vdb_make_phylo <- function(metadata, sampleid_col = "sampleid", skip_seqs = TRUE -get_metaphlan_analyses <- function(con, analysis_ids, schema="public") { +get_metaphlan_analyses <- function(con, analysis_ids, schema="public", include_failures=FALSE) { raw_results <- get_subset_pg_df("mgx_metaphlan", "ia_id", analysis_ids, schema=schema) # we have to deal with samples having no classified -- in those cases metaphlan reports # estimated_number_of_reads_from_the_clade as 0, but to make it work with the rest of our code we modify that # to nreads_input (eg all the input reads are from the unclassified clade) - if (sum(raw_results[raw_results$clade_name == "unclassified", "estimated_number_of_reads_from_the_clade"]) > 0){ - print(raw_results[raw_results$clade_name== "unclassified", "ia_id"]) - stop("some of the above analyses have metaphlan failures where estimated_number of reads from the clade is not zero as expected; please alert of the vdbR delevopers that the metaphlan output has changed") + mpa_failures <- raw_results %>% dplyr::filter(clade_name == "unclassified") %>% dplyr::pull(ia_id) + if (include_failures){ + if (sum(raw_results[raw_results$clade_name == "unclassified", "estimated_number_of_reads_from_the_clade"]) > 0){ + print(raw_results[raw_results$clade_name== "unclassified", "ia_id"]) + stop("some of the above analyses have metaphlan failures where estimated_number of reads from the clade is not zero as expected; please alert of the vdbR delevopers that the metaphlan output has changed") + } + raw_results <- raw_results %>% + dplyr::mutate(estimated_number_of_reads_from_the_clade = ifelse(clade_name == "unclassified", nreads_input, estimated_number_of_reads_from_the_clade)) %>% + dplyr::mutate(clade_name = ifelse(clade_name == "unclassified", "UNCLASSIFIED", clade_name)) + } else{ + if(length(mpa_failures) != 0) { + # get corresponding sampleids for informative error message + for (tbl in c("isabl_api_analysis_targets", "isabl_api_experiment", "isabl_api_sample")){ + if (!tbl %in% ls()) { + get_table_from_database(tbl) + } + } + problem_sample_analyses <- data.frame(analysis_id = mpa_failures) %>% + dplyr::left_join( isabl_api_analysis_targets, by="analysis_id") %>% + dplyr::left_join(isabl_api_experiment %>% dplyr::select(experiment_id = id, sample_id), by="experiment_id") %>% + dplyr::left_join(isabl_api_sample %>% dplyr::select(sample_id = id, identifier), by="sample_id") %>% + dplyr::select(analysis_id, identifier) + stop(paste("Metaphlan detected no taxa in the following sample's analyses; these samples must be removed from your metadata to create a phyloseq object:\n", paste0(capture.output(print(problem_sample_analyses)), collapse = "\n"))) + } } # TODO: parametarize this to get sgb level results - results_with_zeros_handled <- raw_results %>% - dplyr::mutate(estimated_number_of_reads_from_the_clade = ifelse(clade_name == "unclassified", nreads_input, estimated_number_of_reads_from_the_clade)) %>% - dplyr::mutate(clade_name = ifelse(clade_name == "unclassified", "UNCLASSIFIED", clade_name)) %>% + raw_results <- raw_results %>% dplyr::filter(grepl("UNCLASSIFIED", clade_name) | grepl(".*\\|s__.*", clade_name)) %>% dplyr::filter(!grepl(".*t__.*", clade_name)) %>% dplyr::mutate(clade_name = ifelse(clade_name == "UNCLASSIFIED", "k__UNCLASSIFIED", clade_name)) - wide_results <- results_with_zeros_handled %>% + + wide_results <- raw_results %>% dplyr::select(ia_id, clade_name, relative_abundance) %>% tidyr::pivot_wider(names_from = ia_id, values_from = relative_abundance, values_fill = 0) md <- raw_results %>% diff --git a/tests/testthat/test-utilities.R b/tests/testthat/test-utilities.R index 2e21547..9ef588d 100644 --- a/tests/testthat/test-utilities.R +++ b/tests/testthat/test-utilities.R @@ -200,12 +200,21 @@ test_that("test that sample order is retained", { }) -test_that("metaphlan failures are handled", { +test_that("including metaphlan failures is handled", { # 32298 is a sample with no mapped reads skip_if(Sys.getenv("GITHUB_ACTIONS") != "") connect_database(bundled = FALSE) analysis_ids = c(32298, 32299) - res = get_metaphlan_analyses(con = psql_con, analysis_ids) + res = get_metaphlan_analyses(con = psql_con, analysis_ids, include_failures=TRUE) # assert (a) that the column exists on output and (b) that we have full 100% abundance for that sample expect_equal(sum(res[[1]]$`32298`), 100) }) + + +test_that("excluding metaphlan failures is handled", { + # 32298 is a sample with no mapped reads + skip_if(Sys.getenv("GITHUB_ACTIONS") != "") + connect_database(bundled = FALSE) + analysis_ids = c(32298, 32299) + expect_error(get_metaphlan_analyses(con = psql_con, analysis_ids, include_failures=FALSE)) +}) From 1981d8ede3ed23ed6083c9058e940581426d6d13 Mon Sep 17 00:00:00 2001 From: Nick Waters Date: Wed, 17 Sep 2025 16:17:10 -0400 Subject: [PATCH 3/3] Increment version number to 0.13.0 --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index e799a33..96264de 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: vdbR Title: A Package for Acessing and Processing Microbiome Data -Version: 0.12.0 +Version: 0.13.0 Authors@R: c(person(given = "Antonio", family = "Gomes",