diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 10cf5483..5bc3a6b2 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -156,7 +156,10 @@ jobs: sudo chmod 777 "$RENV_PATHS_CACHE" sudo env R_CLI_DYNAMIC=false TERM=dumb Rscript --no-init-file -e 'install.packages(c("renv", "BiocManager"), repos="https://cloud.r-project.org", lib=.Library)' sudo env RENV_PATHS_CACHE="$RENV_PATHS_CACHE" R_CLI_DYNAMIC=false TERM=dumb Rscript --no-init-file -e 'renv::restore(project=".", library=.Library, prompt=FALSE)' - Rscript --no-init-file -e 'for (pkg in c("dada2","Biostrings","ShortRead","vegan","dplyr","MASS","zCompositions")) if (!require(pkg,character.only=TRUE,quietly=TRUE)) stop(pkg, " failed to install")' + Rscript --no-init-file -e 'for (pkg in c("dada2","Biostrings","ShortRead","vegan","MASS","zCompositions")) if (!require(pkg,character.only=TRUE,quietly=TRUE)) stop(pkg, " failed to install")' + + - name: Check DADA2 CSV assembly + run: Rscript --vanilla test/unit/test_dada2_tables.R # Saved on failure too, since the packages built before a failure are the ones # the next attempt should reuse. Root writes the cache, so it is made readable first. @@ -247,7 +250,7 @@ jobs: julia --project=. -e ' using RCall R"cat(.libPaths(), sep=\"\n\")" - for pkg in ["dada2", "dplyr", "vegan", "tibble", "MASS", "zCompositions"] + for pkg in ["dada2", "vegan", "tibble", "MASS", "zCompositions"] R"if (!require($pkg, character.only=TRUE, quietly=TRUE)) stop($pkg, \" not found\")" end ' diff --git a/R/_renv_dependencies.R b/R/_renv_dependencies.R index d9b4c570..20394ab6 100644 --- a/R/_renv_dependencies.R +++ b/R/_renv_dependencies.R @@ -3,7 +3,7 @@ if (FALSE) { library(dada2) library(vegan) - library(dplyr) + # Retain while installer and CI explicitly verify this package. library(tibble) library(MASS) } diff --git a/install.jl b/install.jl index 7e835e62..eb8e9e47 100644 --- a/install.jl +++ b/install.jl @@ -869,14 +869,14 @@ function setup_r_packages(; rebuild::Bool=false)::Tuple{Symbol,String} # renv's exit code: restore can fail on one leaf package and still leave a # working library, or "succeed" against a stale cache. probe = raw""" - pkgs <- c("dada2", "vegan", "dplyr", "tibble") + pkgs <- c("dada2", "vegan", "tibble") cat(paste(pkgs[!vapply(pkgs, requireNamespace, logical(1), quietly = TRUE)], collapse = ",")) """ bad = try strip(read(Cmd(`Rscript -e $probe`; dir = PROJECT_ROOT), String)) catch - "dada2,vegan,dplyr,tibble" + "dada2,vegan,tibble" end isempty(bad) ? (:ok, "") : (:incomplete, String(bad)) diff --git a/renv.lock b/renv.lock index dfac546e..a487543c 100644 --- a/renv.lock +++ b/renv.lock @@ -1668,65 +1668,6 @@ "NeedsCompilation": "yes", "Repository": "CRAN" }, - "dplyr": { - "Package": "dplyr", - "Version": "1.2.1", - "Source": "Repository", - "Type": "Package", - "Title": "A Grammar of Data Manipulation", - "Authors@R": "c( person(\"Hadley\", \"Wickham\", , \"hadley@posit.co\", role = c(\"aut\", \"cre\"), comment = c(ORCID = \"0000-0003-4757-117X\")), person(\"Romain\", \"François\", role = \"aut\", comment = c(ORCID = \"0000-0002-2444-4226\")), person(\"Lionel\", \"Henry\", role = \"aut\"), person(\"Kirill\", \"Müller\", role = \"aut\", comment = c(ORCID = \"0000-0002-1416-3412\")), person(\"Davis\", \"Vaughan\", , \"davis@posit.co\", role = \"aut\", comment = c(ORCID = \"0000-0003-4777-038X\")), person(\"Posit Software, PBC\", role = c(\"cph\", \"fnd\")) )", - "Description": "A fast, consistent tool for working with data frame like objects, both in memory and out of memory.", - "License": "MIT + file LICENSE", - "URL": "https://dplyr.tidyverse.org, https://github.com/tidyverse/dplyr", - "BugReports": "https://github.com/tidyverse/dplyr/issues", - "Depends": [ - "R (>= 4.1.0)" - ], - "Imports": [ - "cli (>= 3.6.2)", - "generics", - "glue (>= 1.3.2)", - "lifecycle (>= 1.0.5)", - "magrittr (>= 1.5)", - "methods", - "pillar (>= 1.9.0)", - "R6", - "rlang (>= 1.1.7)", - "tibble (>= 3.2.0)", - "tidyselect (>= 1.2.0)", - "utils", - "vctrs (>= 0.7.1)" - ], - "Suggests": [ - "broom", - "covr", - "DBI", - "dbplyr (>= 2.2.1)", - "ggplot2", - "knitr", - "Lahman", - "lobstr", - "nycflights13", - "purrr", - "rmarkdown", - "RSQLite", - "stringi (>= 1.7.6)", - "testthat (>= 3.1.5)", - "tidyr (>= 1.3.0)", - "withr" - ], - "VignetteBuilder": "knitr", - "Config/build/compilation-database": "true", - "Config/Needs/website": "tidyverse/tidytemplate", - "Config/testthat/edition": "3", - "Encoding": "UTF-8", - "LazyData": "true", - "RoxygenNote": "7.3.3", - "NeedsCompilation": "yes", - "Author": "Hadley Wickham [aut, cre] (ORCID: ), Romain François [aut] (ORCID: ), Lionel Henry [aut], Kirill Müller [aut] (ORCID: ), Davis Vaughan [aut] (ORCID: ), Posit Software, PBC [cph, fnd]", - "Maintainer": "Hadley Wickham ", - "Repository": "CRAN" - }, "farver": { "Package": "farver", "Version": "2.1.2", diff --git a/src/pipeline/dada2/dada2_functions.r b/src/pipeline/dada2/dada2_functions.r index 8792317a..8fd9fd40 100644 --- a/src/pipeline/dada2/dada2_functions.r +++ b/src/pipeline/dada2/dada2_functions.r @@ -15,8 +15,6 @@ # https://creativecommons.org/licenses/by/4.0/. library(dada2) -library(dplyr) -library(tibble) # library(yaml) # BiocParallel defaults to MulticoreParam, which forks. ShortRead::qa(), reached @@ -214,10 +212,50 @@ write_fasta <- function(seq_table, tables_dir, prefix) { # Writes taxonomy assignments to CSV. Bootstrap confidence values (0-100 per # taxonomic rank) # Returns taxa_df (SeqName + Sequence + taxonomy), used by write_combined_table(). +# +# Preserve left-row order, right-row order within matches, and the original +# column order when assembling the taxonomy and count tables. +base_left_join <- function(x, y, by) { + x_cols <- names(x) + y_cols <- names(y) + stopifnot(!anyDuplicated(x_cols), !anyDuplicated(y_cols), + all(by %in% x_cols), all(by %in% y_cols)) + suffix_names <- function(values, reserved, suffix) { + out <- c(reserved, values) + while (anyDuplicated(out)) { + duplicate <- which(duplicated(out)) + out[duplicate] <- paste0(out[duplicate], suffix) + } + tail(out, length(values)) + } + # Match left_join's repeated suffixing, including names such as Species.x. + x_out <- x_cols + x_aux <- !x_cols %in% by + y_aux <- !y_cols %in% by + x_out[x_aux] <- suffix_names(x_cols[x_aux], c(by, y_cols[y_aux]), ".x") + y_out <- suffix_names(y_cols, x_cols, ".y")[y_aux] + + # Merge only unique internal names, so character selection can never confuse + # a taxonomy column with a sample whose name resembles a generated suffix. + x_ids <- paste0(".x", seq_along(x_cols)) + y_ids <- paste0(".y", seq_along(y_cols)) + names(x) <- x_ids + names(y) <- y_ids + x$.left_position <- seq_len(nrow(x)) + y$.right_position <- seq_len(nrow(y)) + joined <- merge(x, y, by.x = x_ids[match(by, x_cols)], + by.y = y_ids[match(by, y_cols)], all.x = TRUE, sort = FALSE) + joined <- joined[order(joined$.left_position, joined$.right_position), + c(x_ids, y_ids[y_aux]), drop = FALSE] + names(joined) <- c(x_out, y_out) + rownames(joined) <- NULL + joined +} + write_taxa_table <- function(tax_matrix, boot_matrix, index, tables_dir, prefix) { taxa_df <- as.data.frame(tax_matrix, stringsAsFactors = FALSE) taxa_df$Sequence <- rownames(taxa_df) - taxa_df <- dplyr::left_join(index, taxa_df, by = "Sequence") + taxa_df <- base_left_join(index, taxa_df, by = "Sequence") write.csv(taxa_df, file.path(tables_dir, paste0(prefix, ".csv")), quote = FALSE, row.names = FALSE) @@ -225,11 +263,11 @@ write_taxa_table <- function(tax_matrix, boot_matrix, index, tables_dir, prefix) boot_df <- as.data.frame(boot_matrix, stringsAsFactors = FALSE) colnames(boot_df) <- paste0(colnames(boot_df), "_boot") boot_df$Sequence <- rownames(boot_matrix) - boot_df <- dplyr::left_join(index, boot_df, by = "Sequence") + boot_df <- base_left_join(index, boot_df, by = "Sequence") write.csv(boot_df, file.path(tables_dir, paste0(prefix, "_bootstraps.csv")), quote = FALSE, row.names = FALSE) - combined_df <- dplyr::left_join(taxa_df, boot_df, by = c("SeqName", "Sequence")) + combined_df <- base_left_join(taxa_df, boot_df, by = c("SeqName", "Sequence")) write.csv(combined_df, file.path(tables_dir, paste0(prefix, "_combined.csv")), quote = FALSE, row.names = FALSE) @@ -239,10 +277,13 @@ write_taxa_table <- function(tax_matrix, boot_matrix, index, tables_dir, prefix) # Joins taxonomy with per-sample counts write_combined_table <- function(taxa_df, index, seq_table, tables_dir, tax_filename, asv_filename) { seq_t <- as.data.frame(t(seq_table), stringsAsFactors = FALSE) - seq_t <- tibble::rownames_to_column(seq_t, "Sequence") - write.csv(dplyr::left_join(taxa_df, seq_t, by = "Sequence"), + # tibble::rownames_to_column(seq_t, "Sequence") put the row names in a leading + # column and dropped them from the row names; data.frame() does the same. + seq_t <- data.frame(Sequence = rownames(seq_t), seq_t, row.names = NULL, + check.names = FALSE, stringsAsFactors = FALSE) + write.csv(base_left_join(taxa_df, seq_t, by = "Sequence"), file.path(tables_dir, tax_filename), quote = FALSE, row.names = FALSE) - write.csv(dplyr::left_join(index, seq_t, by = "Sequence"), + write.csv(base_left_join(index, seq_t, by = "Sequence"), file.path(tables_dir, asv_filename), quote = FALSE, row.names = FALSE) invisible(NULL) } diff --git a/test/unit/test_dada2_tables.R b/test/unit/test_dada2_tables.R new file mode 100644 index 00000000..228bae1b --- /dev/null +++ b/test/unit/test_dada2_tables.R @@ -0,0 +1,74 @@ +# Regression checks for CSV assembly, without loading DADA2 or dplyr. +# Run from the repository root: Rscript --vanilla test/unit/test_dada2_tables.R +expressions <- parse("src/pipeline/dada2/dada2_functions.r") +wanted <- c("base_left_join", "write_taxa_table", "write_combined_table") +loaded <- character() +for (expr in expressions) { + if (is.call(expr) && identical(expr[[1]], as.name("<-")) && + is.symbol(expr[[2]]) && as.character(expr[[2]]) %in% wanted) { + eval(expr) + loaded <- c(loaded, as.character(expr[[2]])) + } +} +stopifnot(setequal(loaded, wanted)) + +# Unsorted keys, repeated matches, unmatched rows, overlapping column names, +# NA keys, and helper-column name collisions must preserve the intended order. +x <- data.frame(value = 1:4, key = c("b", "a", "z", NA), + .left_join_position = 11:14) +y <- data.frame(key = c("a", "b", "b", NA), value = 21:24, + .right_join_position = 31:34) +got <- base_left_join(x, y, "key") +stopifnot(identical(got$key, c("b", "b", "a", "z", NA_character_)), + identical(got$value.x, c(1L, 1L, 2L, 3L, 4L)), + identical(got$value.y, c(22L, 23L, 21L, NA_integer_, 24L)), + identical(names(got), c("value.x", "key", ".left_join_position", + "value.y", ".right_join_position"))) +stopifnot(nrow(base_left_join(x[FALSE, ], y, "key")) == 0L, + nrow(base_left_join(x, y[FALSE, ], "key")) == nrow(x)) +# Multiple keys: matching only one component must not produce a match. +a <- data.frame(id = c("a", "a"), key = c("x", "y"), value = 1:2) +b <- data.frame(id = c("a", "b"), key = c("y", "x"), other = 3:4) +stopifnot(identical(base_left_join(a, b, c("id", "key"))$other, + c(NA_integer_, 3L))) + +# Rank names and their suffixed variants are also valid sample IDs. +# Repeated suffixing must preserve both the taxonomy strings and count columns. +tax_collision <- data.frame(Sequence = c("AA", "GG"), Species = c("taxA", "taxG")) +count_collision <- data.frame(Sequence = c("GG", "AA"), Species = c(10L, 20L), + Species.x = c(30L, 40L), Species.x.x = c(50L, 60L)) +collision <- base_left_join(tax_collision, count_collision, "Sequence") +stopifnot(identical(names(collision), + c("Sequence", "Species.x.x.x", "Species.y", "Species.x", "Species.x.x")), + identical(collision[[2]], c("taxA", "taxG")), + identical(collision[[3]], c(20L, 10L)), + identical(collision[[4]], c(40L, 30L)), + identical(collision[[5]], c(60L, 50L))) + +# Fixed expected CSV bytes cover every writer output, including non-syntactic +# sample names and an unmatched taxonomy row. No sequencing dependencies needed. +out <- tempfile("dada2-tables-") +dir.create(out) +index <- data.frame(SeqName = c("ASV3", "ASV1", "ASV2"), + Sequence = c("TT", "GG", "AA")) +tax <- matrix(c("G", "A"), ncol = 1, + dimnames = list(c("GG", "AA"), "Kingdom")) +boot <- matrix(c(90L, 80L), ncol = 1, dimnames = dimnames(tax)) +taxa <- write_taxa_table(tax, boot, index, out, "taxonomy") +counts <- matrix(c(1L, 2L, 3L, 4L, 5L, 6L), nrow = 2, + dimnames = list(c("sample-1", "sample 2"), c("AA", "TT", "GG"))) +write_combined_table(taxa, index, counts, out, "tax_counts.csv", "asv_counts.csv") +expected <- list( + taxonomy.csv = c("SeqName,Sequence,Kingdom", "ASV3,TT,NA", "ASV1,GG,G", "ASV2,AA,A"), + taxonomy_bootstraps.csv = c("SeqName,Sequence,Kingdom_boot", "ASV3,TT,NA", "ASV1,GG,90", "ASV2,AA,80"), + taxonomy_combined.csv = c("SeqName,Sequence,Kingdom,Kingdom_boot", "ASV3,TT,NA,NA", "ASV1,GG,G,90", "ASV2,AA,A,80"), + tax_counts.csv = c("SeqName,Sequence,Kingdom,sample-1,sample 2", "ASV3,TT,NA,3,4", "ASV1,GG,G,5,6", "ASV2,AA,A,1,2"), + asv_counts.csv = c("SeqName,Sequence,sample-1,sample 2", "ASV3,TT,3,4", "ASV1,GG,5,6", "ASV2,AA,1,2") +) +for (name in names(expected)) { + path <- file.path(out, name) + bytes <- readBin(path, "raw", n = file.info(path)$size) + stopifnot(identical(bytes, charToRaw(paste0(paste(expected[[name]], collapse = "\n"), "\n")))) +} +unlink(out, recursive = TRUE) +cat("DADA2 join and five CSV fixtures passed\n")