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", diff --git a/R/utilities.R b/R/utilities.R index 41f2972..e01db58 100644 --- a/R/utilities.R +++ b/R/utilities.R @@ -179,11 +179,42 @@ 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) %>% +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) + 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 + 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 <- raw_results %>% dplyr::select(ia_id, clade_name, relative_abundance) %>% tidyr::pivot_wider(names_from = ia_id, values_from = relative_abundance, values_fill = 0) diff --git a/tests/testthat/test-utilities.R b/tests/testthat/test-utilities.R index 6a0962e..9ef588d 100644 --- a/tests/testthat/test-utilities.R +++ b/tests/testthat/test-utilities.R @@ -200,3 +200,21 @@ test_that("test that sample order is retained", { }) +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, 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)) +})