diff --git a/DESCRIPTION b/DESCRIPTION index 6e7156cce..9c872a4af 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,5 +1,5 @@ Package: AMR -Version: 3.0.1.9095 +Version: 3.0.1.9096 Date: 2026-10-03 Title: Antimicrobial Resistance Data Analysis Description: Functions to simplify and standardise antimicrobial resistance (AMR) diff --git a/NEWS.md b/NEWS.md index 7791d0f90..7eb4521cc 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,4 +1,4 @@ -# AMR 3.0.1.9095 +# AMR 3.0.1.9096 Planned as v3.1.0, mid October 2026. @@ -24,6 +24,7 @@ Planned as v3.1.0, mid October 2026. Thus, `mo_domain()` was previously an alias of `mo_kingdom()`; it is now a distinct function returning the domain. Output of `mo_domain()` is therefore unchanged, while `mo_kingdom()` now returns the formal, new kingdom. * Faster parallel computing via the `future` package for `as.sir()` and `wisca()`: a non-sequential plan (e.g. `future::plan(future::multisession)`) must be active before using `parallel = TRUE`. +* `as.mo()` no longer guesses unknown MO codes (these become `NA` with a warning), and MO codes from versions before v2.0.0 are no longer supported ### New * EUCAST 2026 and CLSI 2026 breakpoints: over 5,700 new breakpoints added to the `clinical_breakpoints` data set; EUCAST 2026 is now the default for all MIC and disk diffusion interpretations @@ -41,6 +42,7 @@ Planned as v3.1.0, mid October 2026. * New `amr_course()` to download and unpack course or webinar materials from GitHub in one call * Typed missing value constants `NA_ab_` and `NA_mo_`, for use in pipelines that need missing values of a specific class * New `wisca_plot()` to assess the susceptibility and incidence distributions from the Monte Carlo simulations +* MO codes of all releases since v2.0.0 are now reliably translated to their current code, also for taxa that moved to another domain; a registry guarantees that an MO code never gets another meaning ### Fixed * Streptomycin breakpoints imported from WHONET were coded as streptoduocin (`STR`) instead of streptomycin (`STR1`) @@ -62,6 +64,7 @@ Planned as v3.1.0, mid October 2026. * Translation fixes for Italian CoNS/CoPS names (#256), Dutch antimicrobials, and `sir_df()` foreign-language output (#272) * Fixed some EUCAST Expert Rules, mostly on *S. pneumoniae* * Logical comparison (`==`, `!=`) of `sir` class with a `factor` is now possible on R >= 4.3 +* `mo_current()` and `as.mo()` could not find the current name of some synonyms with multiple sources (e.g. *Chaetomium abuense*) ### Updated * `top_n_microorganisms()`: new `property_for_each` argument for sub-grouping within top *n* groups; rank ordering enforced (only lower taxonomic ranks allowed); fixed `property = NULL` not being accepted; inner filter now tracks original row indices to prevent cross-group contamination diff --git a/R/mo.R b/R/mo.R index f98e17c7a..59cbbf74f 100755 --- a/R/mo.R +++ b/R/mo.R @@ -1254,68 +1254,66 @@ parse_and_convert <- function(x) { } replace_old_mo_codes <- function(x, property) { - # this function transform old MO codes to current codes, such as: - # B_ESCH_COL (AMR v0.5.0) -> B_ESCHR_COLI - ind <- x %like_case% "^[A-Z]_[A-Z_]+$" & !x %in% AMR_env$MO_lookup$mo - if (any(ind, na.rm = TRUE)) { - add_MO_lookup_to_AMR_env() - # get the ones that match - affected <- x[ind] - affected_unique <- unique(affected) - all_direct_matches <- TRUE - # find their new codes, once per code - solved_unique <- unlist(lapply( - strsplit(affected_unique, ""), - function(m) { - domain <- paste0("^", m[1]) - name <- m[3:length(m)] - name[name == "_"] <- " " - name <- tolower(paste0(name, ".*", collapse = "")) - name <- gsub(" .*", " ", name, fixed = TRUE) - name <- paste0("^", name) - results <- AMR_env$MO_lookup$mo[AMR_env$MO_lookup$domain %like_case% domain & - AMR_env$MO_lookup$fullname_lower %like_case% name] - if (length(results) > 1) { - all_direct_matches <<- FALSE - } - results[1L] - } - ), use.names = FALSE) - solved <- solved_unique[match(affected, affected_unique)] - # assign on places where a match was found - x[ind] <- solved - n_matched <- length(affected[!is.na(affected)]) - n_solved <- length(affected[!is.na(solved)]) - n_unsolved <- length(affected[is.na(solved)]) - n_unique <- length(affected_unique[!is.na(affected_unique)]) - if (n_unique < n_matched) { - n_unique <- paste0(n_unique, " unique, ") - } else { - n_unique <- "" + # MO codes of earlier releases (since v2.0.0) that are not in the current data set anymore are translated + # with the MO code registry (internal data `MO_RETIRED_CODES`, see data-raw/microorganisms_files/README.md), + # e.g. a taxon that moved to another domain. + # MO codes are never guessed: a wrong organism is worse than NA. Codes from releases before v2.0.0 are + # not supported (decision by Matthijs S. Berends, 3 October 2026). + add_MO_lookup_to_AMR_env() + not_current <- !is.na(x) & !x %in% AMR_env$MO_lookup$mo + if (!any(not_current)) { + return(x) + } + # list at most 10 items, and escape curly brackets of names such as "Kapabacteria {class}" for the cli package + listing <- function(v) { + v <- gsub("}", "}}", gsub("{", "{{", v, fixed = TRUE), fixed = TRUE) + out <- vector_and(utils::head(v, 10), quotes = FALSE, sort = FALSE) + if (length(v) > 10) { + out <- paste0(out, ", and ", length(v) - 10, " more") } - if (property != "mo") { - warning_( - "in {.help [{.fun mo_", property, "}](AMR::mo_", property, ")}: the input contained ", n_matched, - " old MO code", ifelse(n_matched == 1, "", "s"), - " (", n_unique, "from another AMR package version). ", - "Please update your MO codes with {.help [{.fun as.mo}](AMR::as.mo)} to increase speed." - ) - } else { - warning_( - "in {.help [{.fun as.mo}](AMR::as.mo)}: the input contained ", n_matched, - " old MO code", ifelse(n_matched == 1, "", "s"), - " (", n_unique, "from another AMR package version). ", - n_solved, " old MO code", ifelse(n_solved == 1, "", "s"), - ifelse(n_solved == 1, " was", " were"), - ifelse(all_direct_matches, " updated ", font_bold(" guessed ")), - "to ", ifelse(n_solved == 1, "a ", ""), - "currently used MO code", ifelse(n_solved == 1, "", "s"), - ifelse(n_unsolved > 0, - paste0(" and ", n_unsolved, " old MO code", ifelse(n_unsolved == 1, "", "s"), " could not be updated."), - "." - ) + out + } + in_function <- ifelse(property == "mo", + "in {.help [{.fun as.mo}](AMR::as.mo)}: ", + paste0("in {.help [{.fun mo_", property, "}](AMR::mo_", property, ")}: ") + ) + + # codes from earlier releases + is_retired <- not_current & x %in% MO_RETIRED_CODES$old_mo + if (any(is_retired)) { + retired <- unique(x[is_retired]) + current <- MO_RETIRED_CODES$mo[match(retired, MO_RETIRED_CODES$old_mo)] + x[is_retired] <- current[match(x[is_retired], retired)] + msg <- paste0( + in_function, "the input contained ", length(retired), " MO code", ifelse(length(retired) == 1, "", "s"), + " from an earlier AMR package version. " + ) + if (any(!is.na(current))) { + msg <- paste0( + msg, sum(!is.na(current)), ifelse(sum(!is.na(current)) == 1, " was", " were"), + " translated to the current MO code of the same taxon (", listing(paste0(retired[!is.na(current)], " -> ", current[!is.na(current)])), "). " ) } + if (any(is.na(current))) { + gone <- MO_RETIRED_CODES$fullname[match(retired[is.na(current)], MO_RETIRED_CODES$old_mo)] + msg <- paste0( + msg, sum(is.na(current)), ifelse(sum(is.na(current)) == 1, " was", " were"), + " set to NA, as the taxon is not in the current data set anymore (", listing(paste0(retired[is.na(current)], " = ", gone)), "). " + ) + } + warning_(msg, "Please update your MO codes with {.help [{.fun as.mo}](AMR::as.mo)}.") + } + + # anything else that looks like an MO code, but was never part of a release since v2.0.0 + is_unknown_code <- not_current & !is_retired & x %like_case% "^[A-Z]{1,2}_[][A-Z0-9_-]+$" + if (any(is_unknown_code)) { + unknown <- unique(x[is_unknown_code]) + x[is_unknown_code] <- NA_character_ + warning_( + in_function, "the input contained ", length(unknown), " unknown MO code", ifelse(length(unknown) == 1, "", "s"), + " (", listing(unknown), "), which ", ifelse(length(unknown) == 1, "was", "were"), + " set to NA. MO codes of AMR package versions before v2.0.0 are not supported." + ) } x } @@ -1395,15 +1393,24 @@ synonym_mo_to_accepted_mo <- function(x, fill_in_accepted = FALSE, dataset = AMR x_mycobank <- dataset$mycobank_renamed_to[match(out, dataset$mo)] x_lpsn <- dataset$lpsn_renamed_to[match(out, dataset$mo)] - out[must_be_corrected & !is.na(x_gbif)] <- dataset$mo[match(x_gbif[must_be_corrected & !is.na(x_gbif)], dataset$gbif)] - out[must_be_corrected & !is.na(x_mycobank)] <- dataset$mo[match(x_mycobank[must_be_corrected & !is.na(x_mycobank)], dataset$mycobank)] - out[must_be_corrected & !is.na(x_lpsn)] <- dataset$mo[match(x_lpsn[must_be_corrected & !is.na(x_lpsn)], dataset$lpsn)] + # priority LPSN > MycoBank > GBIF, but a source may only overwrite the result if its target exists in the + # data set (before 3.1.0, e.g. a MycoBank target that was not included replaced a valid GBIF target by NA) + by_gbif <- dataset$mo[match(x_gbif, dataset$gbif, incomparables = NA)] + by_mycobank <- dataset$mo[match(x_mycobank, dataset$mycobank, incomparables = NA)] + by_lpsn <- dataset$mo[match(x_lpsn, dataset$lpsn, incomparables = NA)] + out[must_be_corrected & !is.na(by_gbif)] <- by_gbif[must_be_corrected & !is.na(by_gbif)] + out[must_be_corrected & !is.na(by_mycobank)] <- by_mycobank[must_be_corrected & !is.na(by_mycobank)] + out[must_be_corrected & !is.na(by_lpsn)] <- by_lpsn[must_be_corrected & !is.na(by_lpsn)] is_still_synonym <- dataset$status[match(out, dataset$mo)] == "synonym" } x_no_synonym <- dataset$status[match(x, dataset$mo)] != "synonym" out[x_no_synonym] <- NA_character_ + if (isFALSE(fill_in_accepted)) { + # a synonym without any known current name is not 'renamed' to itself + out[!is.na(out) & out == x] <- NA_character_ + } if (isTRUE(fill_in_accepted)) { out[!is.na(x_no_synonym) & x_no_synonym] <- x[!is.na(x_no_synonym) & x_no_synonym] } diff --git a/R/sysdata.rda b/R/sysdata.rda index 4bb944625..899b93240 100755 Binary files a/R/sysdata.rda and b/R/sysdata.rda differ diff --git a/data-raw/_pre_commit_checks.R b/data-raw/_pre_commit_checks.R index fccb6b7c2..547623dce 100644 --- a/data-raw/_pre_commit_checks.R +++ b/data-raw/_pre_commit_checks.R @@ -624,6 +624,28 @@ usethis::use_data(antimicrobials, overwrite = TRUE, version = 2, compress = "xz" pre_commit_lst$AB_LOOKUP <- create_AB_AV_lookup(antimicrobials) pre_commit_lst$AV_LOOKUP <- create_AB_AV_lookup(antivirals) +# MO codes of earlier releases that are not in the current data set anymore, with the current code of the same +# taxon (see data-raw/microorganisms_files/README.md). as.mo() uses this to translate them, instead of guessing. +source("tests/testthat/helper-microorganisms.R") +mo_registry <- read_mo_registry(".") +mo_renames <- read_mo_renames(".") +retired <- mo_registry[!mo_registry$mo %in% as.character(microorganisms$mo), c("mo", "fullname", "rank"), drop = FALSE] +renamed <- match(retired$mo, mo_renames$mo) +retired$fullname[!is.na(renamed)] <- mo_renames$new_name[renamed[!is.na(renamed)]] +# the same name, or else the same name without a rank suffix (such as "Kapabacteria {class}") and the same rank +current_mo <- as.character(microorganisms$mo)[match(retired$fullname, microorganisms$fullname)] +by_rank <- as.character(microorganisms$mo)[match( + paste(mo_name_without_suffix(retired$fullname), retired$rank), + paste(mo_name_without_suffix(microorganisms$fullname), microorganisms$rank) +)] +pre_commit_lst$MO_RETIRED_CODES <- data.frame( + old_mo = retired$mo, + mo = ifelse(is.na(current_mo), by_rank, current_mo), + fullname = retired$fullname, + stringsAsFactors = FALSE +) +rm(mo_registry, mo_renames, retired, renamed, current_mo, by_rank) + # Export to package as internal data ---- # usethis::use_data() must receive unquoted object names, which is not flexible at all. # we'll use good old base::save() instead diff --git a/data-raw/_reproduction_scripts/reproduction_of_microorganisms.R b/data-raw/_reproduction_scripts/reproduction_of_microorganisms.R index 79fa96e30..9a6dbbec7 100644 --- a/data-raw/_reproduction_scripts/reproduction_of_microorganisms.R +++ b/data-raw/_reproduction_scripts/reproduction_of_microorganisms.R @@ -115,6 +115,9 @@ devtools::load_all(".") # to load the AMR package # keep the previous data set at hand, `AMR::microorganisms` will be replaced at the end of this script microorganisms_old <- AMR::microorganisms +# the integrity rules of the data set and the MO code registry, shared with the unit tests +source("tests/testthat/helper-microorganisms.R") + # Helper functions -------------------------------------------------------------------------------- @@ -143,6 +146,11 @@ review <- function(x, title = deparse(substitute(x))) { invisible(x) } +# taxonomic names in ASCII and without quotes +ascii_names <- function(x) { + gsub("[\"'`]", "", stringi::stri_trans_general(x, "Latin-ASCII")) +} + # the priority of sources, used everywhere in this script source_priority <- c("LPSN" = 1L, "MycoBank" = 2L, "GBIF" = 3L, "manually added" = 4L, "inferred" = 5L) source_prio <- function(source) { @@ -1052,8 +1060,11 @@ saveRDS(taxonomy_mycobank, "data-raw/taxonomy_mycobank.rds", version = 2) # Read GBIF data ---------------------------------------------------------------------------------- +# @resume-block gbif_raw after=taxonomy_gbif +# (needed by `current_gbif` when resuming from a checkpoint with run_microorganisms_build.R) taxonomy_gbif.bak <- vroom(file_gbif, guess_max = 5e5) colnames(taxonomy_gbif.bak) <- gsub(".*:(.*)", "\\1", colnames(taxonomy_gbif.bak)) +# @end-resume-block # include all fungal orders from the mycobank db include_fungal_orders <- unique(taxonomy_mycobank$order[ @@ -1282,10 +1293,7 @@ taxonomy_gbif <- readRDS("data-raw/taxonomy_gbif.rds") # Add full names ---------------------------------------------------------------------------------- # names must be ASCII before they are made unique, otherwise e.g. 'Fusarium aloës' and 'Fusarium aloes' -# both end up as 'Fusarium aloes' in the final data set (which happened to 12 names in 2026) -ascii_names <- function(x) { - gsub("[\"'`]", "", stringi::stri_trans_general(x, "Latin-ASCII")) -} +# both end up as 'Fusarium aloes' in the final data set (which happened to 12 names in 2026), see ascii_names() taxonomy_gbif <- taxonomy_gbif %>% mutate(across(domain:subspecies, ascii_names)) taxonomy_lpsn <- taxonomy_lpsn %>% @@ -1639,6 +1647,8 @@ saveRDS(taxonomy, "data-raw/taxonomy0.rds") # Ensure every referenced rank has its own row (domain through species). # Where possible, enrich with GBIF identifiers. +# @resume-block current_gbif after=taxonomy1 +# (needed by add_missing_parents() after filtering, when resuming from a checkpoint) current_gbif <- taxonomy_gbif.bak %>% filter(is.na(acceptedNameUsageID)) %>% mutate( @@ -1663,6 +1673,7 @@ current_gbif <- taxonomy_gbif.bak %>% phylum = if_else(taxonRank == "phylum", canonical, phylum) ) %>% select(-canonical) +# @end-resume-block taxonomy <- add_missing_parents(taxonomy, current_gbif) @@ -2035,6 +2046,8 @@ taxonomy$rank[is.na(taxonomy$rank)] <- "domain" # this part is required here, because it's needed for filtering on 'relevant' species to keep later on # (`pathogens` was read at the start of this script) +# @resume-block prevalence after=taxonomy2 +# (compute_prevalence() is used again later on, when resuming from a checkpoint) # get all established, both old and current taxonomic names established <- pathogens %>% filter(status == "established") %>% @@ -2147,6 +2160,7 @@ compute_prevalence <- function(taxonomy) { mutate(prevalence = if_else(rank == "genus", pmin(prevalence, gen_prevalence, na.rm = TRUE), prevalence)) %>% select(-gen_prevalence) } +# @end-resume-block taxonomy <- compute_prevalence(taxonomy) table(taxonomy$prevalence, useNA = "always") @@ -2513,6 +2527,77 @@ taxonomy <- taxonomy %>% distinct(fullname, .keep_all = TRUE) message("Removed ", n_before - nrow(taxonomy), " nonsense names") +# Break cycles of synonyms ------------------------------------------------------------------------ + +# Sources can contradict each other, e.g. in 2026 GBIF listed Capillidium denaeosporus as synonym of +# Conidiobolus denaeosporum and also the other way around, so that the 'current name' depended on the starting +# point. In a cycle, the statement of the most authoritative source counts: that record stays a synonym, and its +# target becomes accepted. +next_synonym_step <- function(df) { + by_gbif <- df$fullname[match(df$gbif_renamed_to, df$gbif, incomparables = NA)] + by_mycobank <- df$fullname[match(df$mycobank_renamed_to, df$mycobank, incomparables = NA)] + by_lpsn <- df$fullname[match(df$lpsn_renamed_to, df$lpsn, incomparables = NA)] + out <- coalesce(by_lpsn, by_mycobank, by_gbif) + out[df$status != "synonym"] <- NA_character_ + setNames(out, df$fullname) +} +synonym_next <- next_synonym_step(taxonomy) +synonym_cycles <- list() +for (start in names(synonym_next)[!is.na(synonym_next)]) { + path <- start + repeat { + nxt <- unname(synonym_next[path[length(path)]]) + if (is.na(nxt)) { + break + } + if (nxt %in% path) { + cycle <- sort(path[match(nxt, path):length(path)]) + synonym_cycles[[paste(cycle, collapse = " | ")]] <- cycle + break + } + path <- c(path, nxt) + if (length(path) > 25) { + break + } + } +} +# If multiple records in the cycle come from the most authoritative source, the sources contradict themselves and +# the direction cannot be determined automatically (in 2026, even the authors in `ref` of these GBIF records were +# unreliable). Then ALL records in the cycle become accepted, so that no wrong current name is given, and they are +# shown for a human decision. +cycle_fixes <- bind_rows(lapply(synonym_cycles, function(cycle_members) { + # (not named `cycle`, since tibble() below would then use its own new column `cycle`) + members <- taxonomy[match(cycle_members, taxonomy$fullname), , drop = FALSE] + best <- members[source_prio(members$source) == min(source_prio(members$source)), , drop = FALSE] + if (nrow(best) == 1) { + tibble( + cycle = paste(cycle_members, collapse = " | "), + made_accepted = unname(synonym_next[best$fullname]), + reason = paste("statement of", best$fullname, "from", best$source, "(most authoritative source)") + ) + } else { + tibble( + cycle = paste(cycle_members, collapse = " | "), + made_accepted = cycle_members, + reason = paste("contradicting statements of", toString(unique(best$source)), "- all accepted, NEEDS A HUMAN DECISION") + ) + } +})) +review(cycle_fixes, "Cycles of synonyms") +if (nrow(cycle_fixes) > 0) { + taxonomy <- taxonomy %>% + mutate( + fix = fullname %in% cycle_fixes$made_accepted, + status = if_else(fix, "accepted", status), + lpsn_renamed_to = if_else(fix, NA_character_, lpsn_renamed_to), + mycobank_renamed_to = if_else(fix, NA_character_, mycobank_renamed_to), + gbif_renamed_to = if_else(fix, NA_character_, gbif_renamed_to) + ) %>% + select(-fix) +} +rm(synonym_next, synonym_cycles, cycle_fixes) + + # Fix genera that are synonyms while they contain accepted species -------------------------------- # e.g. in 2026, MycoBank listed the genus Blastomyces as a synonym, while Blastomyces dermatitidis was accepted @@ -2560,151 +2645,59 @@ mo_domain <- taxonomy %>% ) mo_domain -# `domain` exists in the data set since 2026, so codes are kept per domain -# (in the 2026 build, this was `AMR::microorganisms %>% rename(domain = kingdom)`, which would now fail) -existing_mo_tbl <- microorganisms_old %>% - mutate(mo = as.character(mo)) - +# @resume-block mo_registry after=taxonomy2c +# (the MO code registry is used by the final checks, when resuming from a checkpoint) # parts of existing MO codes, split on the underscore (genus codes can contain digits or hyphens, such as # B_RBCLM1 or B_CND-P, so a regex with only [A-Z] does not work and caused repeated elements like # B_SCLLM_CNNM_LNSM_LNSM_LNSM in earlier versions) mo_part <- function(mo, i) { vapply(strsplit(as.character(mo), "_", fixed = TRUE), function(x) if (length(x) >= i) x[i] else NA_character_, character(1)) } +# parts of a name, e.g. the species of 'Escherichia coli' +mo_part_of_name <- function(x, i) { + vapply(strsplit(x, " ", fixed = TRUE), function(x) if (length(x) >= i) x[i] else "", character(1)) +} + +# Existing codes come from the MO code registry (data-raw/microorganisms_files/), which contains every code of +# every release since v2.0.0, and NOT from the previous data set: codes that only existed in a development +# version carry no weight (decision by Matthijs S. Berends, 3 October 2026). A registered code is never given +# to another taxon, also not when its taxon is not in the new data set anymore, since users store these codes. +mo_registry <- read_mo_registry(".") +mo_renames <- read_mo_renames(".") +if (is.null(mo_registry)) { + stop("MO code registry not found, see data-raw/microorganisms_files/README.md", call. = FALSE) +} +release_order <- function(release) { + as.integer(factor(release, levels = unique(release[order(numeric_version(sub("^v", "", release)))]))) +} +existing_mo_tbl <- mo_registry %>% + # approved renames, e.g. a corrected spelling: the code continues with the new name + left_join(mo_renames %>% select(mo, new_name), by = "mo") %>% + mutate( + fullname = coalesce(new_name, fullname), + species = if_else(!is.na(new_name) & rank %in% c("species", "subspecies"), mo_part_of_name(new_name, 2), species), + subspecies = if_else(!is.na(new_name) & rank == "subspecies", mo_part_of_name(new_name, 3), subspecies) + ) %>% + select(-new_name) %>% + # if a taxon had multiple codes over time, the code of the most recent release is used + arrange(desc(release_order(last_release)), desc(release_order(first_release))) %>% + distinct(domain, rank, fullname, .keep_all = TRUE) + # all codes that have ever been used must never be given to another taxon, also not when their taxon # is not in the new data set anymore (users may have stored old codes) -reserved_genus_codes <- existing_mo_tbl %>% +reserved_genus_codes <- mo_registry %>% filter(rank %in% c("genus", "species", "subspecies")) %>% transmute(domain, code = mo_part(mo, 2)) %>% distinct() -reserved_species_codes <- existing_mo_tbl %>% +reserved_species_codes <- mo_registry %>% filter(rank %in% c("species", "subspecies")) %>% transmute(domain, genus_code = mo_part(mo, 2), code = mo_part(mo, 3)) %>% distinct() +# @end-resume-block # phylum until family are abbreviated with 8 characters and prefixed with their rank -# Phylum - keep old and fill up for new ones -mo_phylum <- taxonomy %>% - filter(rank == "phylum") %>% - distinct(domain, phylum) %>% - left_join( - existing_mo_tbl %>% - filter(rank == "phylum") %>% - transmute( - domain, - phylum = fullname, - mo_old = gsub("^[A-Z]{1,2}_", "", as.character(mo)) - ), - by = c("domain", "phylum") - ) %>% - group_by(domain) %>% - mutate( - mo_phylum8 = AMR:::abbreviate_mo(phylum, minlength = 8, prefix = "[PHL]_"), - mo_phylum9 = AMR:::abbreviate_mo(phylum, minlength = 9, prefix = "[PHL]_"), - mo_phylum = if_else(!is.na(mo_old), mo_old, mo_phylum8), - mo_duplicated = duplicated(mo_phylum), - mo_phylum = if_else(mo_duplicated, mo_phylum9, mo_phylum), - mo_duplicated = duplicated(mo_phylum) - ) %>% - ungroup() -if (any(mo_phylum$mo_duplicated, na.rm = TRUE)) { - stop("Duplicate MO codes for phylum!") -} -mo_phylum <- mo_phylum %>% - select(domain, phylum, mo_phylum) - -# Class - keep old and fill up for new ones -mo_class <- taxonomy %>% - filter(rank == "class") %>% - distinct(domain, class) %>% - left_join( - existing_mo_tbl %>% - filter(rank == "class") %>% - transmute( - domain, - class = fullname, - mo_old = gsub("^[A-Z]{1,2}_", "", as.character(mo)) - ), - by = c("domain", "class") - ) %>% - group_by(domain) %>% - mutate( - mo_class8 = AMR:::abbreviate_mo(class, minlength = 8, prefix = "[CLS]_"), - mo_class9 = AMR:::abbreviate_mo(class, minlength = 9, prefix = "[CLS]_"), - mo_class = if_else(!is.na(mo_old), mo_old, mo_class8), - mo_duplicated = duplicated(mo_class), - mo_class = if_else(mo_duplicated, mo_class9, mo_class), - mo_duplicated = duplicated(mo_class) - ) %>% - ungroup() -if (any(mo_class$mo_duplicated, na.rm = TRUE)) { - stop("Duplicate MO codes for class!") -} -mo_class <- mo_class %>% - select(domain, class, mo_class) - -# Order - keep old and fill up for new ones -mo_order <- taxonomy %>% - filter(rank == "order") %>% - distinct(domain, order) %>% - left_join( - existing_mo_tbl %>% - filter(rank == "order") %>% - transmute( - domain, - order = fullname, - mo_old = gsub("^[A-Z]{1,2}_", "", as.character(mo)) - ), - by = c("domain", "order") - ) %>% - group_by(domain) %>% - mutate( - mo_order8 = AMR:::abbreviate_mo(order, minlength = 8, prefix = "[ORD]_"), - mo_order9 = AMR:::abbreviate_mo(order, minlength = 9, prefix = "[ORD]_"), - mo_order = if_else(!is.na(mo_old), mo_old, mo_order8), - mo_duplicated = duplicated(mo_order), - mo_order = if_else(mo_duplicated, mo_order9, mo_order), - mo_duplicated = duplicated(mo_order) - ) %>% - ungroup() -if (any(mo_order$mo_duplicated, na.rm = TRUE)) { - stop("Duplicate MO codes for order!") -} -mo_order <- mo_order %>% - select(domain, order, mo_order) - -# Family - keep old and fill up for new ones -mo_family <- taxonomy %>% - filter(rank == "family") %>% - distinct(domain, family) %>% - left_join( - existing_mo_tbl %>% - filter(rank == "family") %>% - transmute( - domain, - family = fullname, - mo_old = gsub("^[A-Z]{1,2}_", "", as.character(mo)) - ), - by = c("domain", "family") - ) %>% - group_by(domain) %>% - mutate( - mo_family8 = AMR:::abbreviate_mo(family, minlength = 8, prefix = "[FAM]_"), - mo_family9 = AMR:::abbreviate_mo(family, minlength = 9, prefix = "[FAM]_"), - mo_family = if_else(!is.na(mo_old), mo_old, mo_family8), - mo_duplicated = duplicated(mo_family), - mo_family = if_else(mo_duplicated, mo_family9, mo_family), - mo_duplicated = duplicated(mo_family) - ) %>% - ungroup() -if (any(mo_family$mo_duplicated, na.rm = TRUE)) { - stop("Duplicate MO codes for family!") -} -mo_family <- mo_family %>% - select(domain, family, mo_family) - # keep the old code where available, and give every new taxon the first candidate code that is not in use # and was never used before (`reserved`) - if all candidates are taken, a number is added to the first one assign_codes <- function(old, candidates, reserved = character(0)) { @@ -2726,6 +2719,47 @@ assign_codes <- function(old, candidates, reserved = character(0)) { out } +# Phylum until family - keep old (from the registry) and fill up for new ones, never with a registered code +higher_rank_codes <- function(rank_name, tag) { + reserved <- mo_registry %>% + filter(rank == rank_name) %>% + transmute(domain, code = sub("^[A-Z]{1,2}_", "", mo)) + out <- taxonomy %>% + filter(rank == rank_name) %>% + distinct(domain, name = .data[[rank_name]]) %>% + left_join( + existing_mo_tbl %>% + filter(rank == rank_name) %>% + transmute(domain, name = mo_name_without_suffix(fullname), mo_old = sub("^[A-Z]{1,2}_", "", mo)) %>% + distinct(domain, name, .keep_all = TRUE), + by = c("domain", "name") + ) %>% + group_by(domain) %>% + mutate( + code = assign_codes( + old = mo_old, + candidates = Map( + c, + AMR:::abbreviate_mo(name, minlength = 8, prefix = paste0("[", tag, "]_")), + AMR:::abbreviate_mo(name, minlength = 9, prefix = paste0("[", tag, "]_")) + ), + reserved = reserved$code[reserved$domain == cur_group()$domain] + ) + ) %>% + ungroup() + if (anyDuplicated(paste(out$domain, out$code)) > 0 || anyNA(out$code)) { + stop("Duplicate MO codes for ", rank_name, "!", call. = FALSE) + } + out %>% + select(domain, name, code) %>% + setNames(c("domain", rank_name, paste0("mo_", rank_name))) +} +mo_phylum <- higher_rank_codes("phylum", "PHL") +mo_class <- higher_rank_codes("class", "CLS") +mo_order <- higher_rank_codes("order", "ORD") +mo_family <- higher_rank_codes("family", "FAM") + + # construct code part for genus - keep old code where available and generate new ones where needed mo_genus <- taxonomy %>% filter(rank == "genus") %>% @@ -2863,6 +2897,13 @@ mo_subspecies <- mo_subspecies %>% mo_unknown <- existing_mo_tbl %>% filter(fullname %like% "unknown") %>% transmute(fullname, mo_unknown = as.character(mo)) +# unknowns that were added after the last release keep their code too, as long as it is not registered +mo_unknown <- mo_unknown %>% + bind_rows( + microorganisms_old %>% + filter(fullname %like% "unknown", !fullname %in% mo_unknown$fullname, !mo %in% mo_registry$mo) %>% + transmute(fullname, mo_unknown = as.character(mo)) + ) # apply the new codes! taxonomy1 <- taxonomy @@ -3588,23 +3629,33 @@ microorganisms_old %>% select(mo, fullname, rank, status, prevalence) %>% review("Previously manually added taxa that are not in the new data set") -# LPSN records must be prokaryotes and MycoBank records must be fungi, also in their MO code (issue #309) -if (!all(taxonomy$domain[taxonomy$source == "LPSN"] %in% c("Bacteria", "Archaea")) || - !all(taxonomy$mo[taxonomy$source == "LPSN"] %like% "^(B|A)_") || - !all(taxonomy$domain[taxonomy$source == "MycoBank"] == "Fungi") || - !all(taxonomy$mo[taxonomy$source == "MycoBank"] %like% "^F_")) { - warning("LPSN records outside prokaryotes or MycoBank records outside fungi, check 'Combine the datasets'", call. = FALSE) -} - -# these must never fail, as the package relies on them -stopifnot( - "MO codes must be unique" = !anyDuplicated(taxonomy$mo), - "full names must be unique" = !anyDuplicated(taxonomy$fullname), - "MO codes must not be empty" = !any(taxonomy$mo %in% c("", NA)), - "prevalence must not be empty" = !anyNA(taxonomy$prevalence), - "status must be accepted, synonym or unknown" = all(taxonomy$status %in% c("accepted", "synonym", "unknown")), - "sources must be known" = all(taxonomy$source %in% c("LPSN", "MycoBank", "GBIF", "manually added")) +# 'renamed to' identifiers of records that were removed after the earlier clean-up (e.g. by deduplication) lead +# nowhere, so remove them again: these synonyms then have no current name, instead of a broken one +n_dangling <- sum( + (!is.na(taxonomy$lpsn_renamed_to) & !taxonomy$lpsn_renamed_to %in% taxonomy$lpsn) | + (!is.na(taxonomy$mycobank_renamed_to) & !taxonomy$mycobank_renamed_to %in% taxonomy$mycobank) | + (!is.na(taxonomy$gbif_renamed_to) & !taxonomy$gbif_renamed_to %in% taxonomy$gbif) ) +taxonomy <- taxonomy %>% + mutate( + lpsn_renamed_to = if_else(lpsn_renamed_to %in% lpsn, lpsn_renamed_to, NA_character_), + mycobank_renamed_to = if_else(mycobank_renamed_to %in% mycobank, mycobank_renamed_to, NA_character_), + gbif_renamed_to = if_else(gbif_renamed_to %in% gbif, gbif_renamed_to, NA_character_) + ) +message("Removed dangling 'renamed to' identifiers of ", n_dangling, " records") + +# synonyms without a current name are no error, but every one of them is a name that as.mo() cannot update +review(mo_synonyms_without_current_name(taxonomy), "Synonyms without a current name") + +# the integrity rules of the data set and the MO code registry (shared with the unit tests in +# tests/testthat/helper-microorganisms.R): these must never fail, as the package and its users rely on them +integrity_issues <- mo_integrity_issues(taxonomy, registry = mo_registry, renames = mo_renames) +if (!is.null(mo_integrity_report(integrity_issues))) { + stop("The new data set breaks integrity rules, fix the cause in this script:\n", + mo_integrity_report(integrity_issues), + call. = FALSE + ) +} # Update other data sets -------------------------------------------------------------------------- @@ -3676,11 +3727,13 @@ taxonomy <- taxonomy %>% AMR:::dataset_UTF8_to_ASCII() # check again after the conversion to ASCII -stopifnot( - "full names must be unique after conversion to ASCII" = !anyDuplicated(taxonomy$fullname), - "MO codes must be unique after conversion to ASCII" = !anyDuplicated(taxonomy$mo), - "MO codes may only contain A-Z, 0-9, _, -, [ and ]" = all(taxonomy$mo %like_case% "^[][A-Z0-9_-]+$") -) +integrity_issues <- mo_integrity_issues(taxonomy, registry = mo_registry, renames = mo_renames) +if (!is.null(mo_integrity_report(integrity_issues))) { + stop("The data set breaks integrity rules after the conversion to ASCII:\n", + mo_integrity_report(integrity_issues), + call. = FALSE + ) +} microorganisms <- taxonomy diff --git a/data-raw/_reproduction_scripts/reproduction_of_microorganisms_CLAUDE.md b/data-raw/_reproduction_scripts/reproduction_of_microorganisms_CLAUDE.md index 1bdbb2897..de7b6b47b 100644 --- a/data-raw/_reproduction_scripts/reproduction_of_microorganisms_CLAUDE.md +++ b/data-raw/_reproduction_scripts/reproduction_of_microorganisms_CLAUDE.md @@ -65,20 +65,22 @@ These are needed for `TAXONOMY_VERSION` in `R/aa_globals.R` later. Write them do ## Step 2: run the build The script takes several hours (most of it is LPSN scraping, which is cached in -`data-raw/lpsn_scrape_cache.rds`, so a restart is fast). Run it detached from your tool calls, so that -it survives time limits, and write everything to a log: +`data-raw/lpsn_scrape_cache.rds`, so a restart is fast). Never `source()` the whole script. Use the runner +`data-raw/_reproduction_scripts/run_microorganisms_build.R`: it cuts the script into chunks at its +checkpoints, runs them one by one in a single R process, logs each chunk with its duration in +`/logs/progress.tsv`, stops with the name of the failing chunk, and can resume from any checkpoint. +See the chunks first with `--list`. Run it detached from your tool calls, so that it survives time limits: ```bash mkdir -p /logs cd -nohup Rscript -e 'options(AMR_build_view = FALSE, AMR_build_print_rows = 200, AMR_build_review_dir = "", warn = 1, width = 200); source("data-raw/_reproduction_scripts/reproduction_of_microorganisms.R", echo = TRUE, max.deparse.length = Inf, keep.source = TRUE)' > /logs/build_1.log 2>&1 & +nohup Rscript data-raw/_reproduction_scripts/run_microorganisms_build.R --review-dir= --log-dir=/logs > /logs/build_1.log 2>&1 & ``` -Monitor the log at sensible intervals (every 10 to 20 minutes, not continuously). Read every -`>>> REVIEW:` block as soon as it appears; you do not have to wait until the end to start thinking. -With `AMR_build_review_dir` set, every review block is also saved in full as a CSV file in -``, named after its title. These files are part of the PR. Set the same option in every -resumed run. +Monitor the log and `progress.tsv` at sensible intervals (every 10 to 20 minutes, not continuously). Read +every `>>> REVIEW:` block as soon as it appears; you do not have to wait until the end to start thinking. +Every review block is also saved in full as a CSV file in ``, named after its title. These +files are part of the PR. ### Intermediate results @@ -105,24 +107,18 @@ and must not be committed. ### If the run fails -1. Read the error and the 100 lines before it. Find the root cause in the script, fix it, and - check the syntax with `Rscript -e 'invisible(parse("