diff --git a/bin/workflow_glue_r/NAMESPACE b/bin/workflow_glue_r/NAMESPACE index e22afa7..f4b73bb 100644 --- a/bin/workflow_glue_r/NAMESPACE +++ b/bin/workflow_glue_r/NAMESPACE @@ -4,7 +4,6 @@ export(bambu_discovery_enabled) export(bambu_filter_transcripts) export(bambu_normalise_tsv_df) export(bambu_resolve_inputs) -export(bambu_strip_alias) export(bambu_write_outputs) export(de_analysis_arg_parser) export(de_validate_inputs) diff --git a/bin/workflow_glue_r/R/bambu.R b/bin/workflow_glue_r/R/bambu.R index 1e4c5b8..8cb87b0 100644 --- a/bin/workflow_glue_r/R/bambu.R +++ b/bin/workflow_glue_r/R/bambu.R @@ -192,10 +192,19 @@ bambu_effective_threads <- function(args, bam_count) { } bambu_filter_transcripts <- function(se) { - counts_mat <- SummarizedExperiment::assays(se)$counts - full_length_mat <- SummarizedExperiment::assays(se)$fullLengthCounts + assays <- SummarizedExperiment::assays(se) + counts_mat <- assays$counts + full_length_mat <- assays$fullLengthCounts - gene_ids <- SummarizedExperiment::rowData(se)$GENEID + row_data <- SummarizedExperiment::rowData(se) + if (!"GENEID" %in% names(row_data)) { + stop("rowData(se) does not contain required column 'GENEID'.") + } + if (is.null(full_length_mat) && is.null(counts_mat)) { + stop("Neither counts nor fullLengthCounts assay found in bambu output.") + } + + gene_ids <- row_data$GENEID qc_stats <- list( total_transcripts_before_filter = nrow(se), total_genes_before_filter = length(unique(gene_ids)), @@ -217,6 +226,53 @@ bambu_filter_transcripts <- function(se) { qc_stats$transcripts_filtered <- sum(!keep_idx) se <- se[keep_idx, ] + + # Keep incompatible gene-level counts in sync with the filtered transcript set. + # transcriptToGeneExpression() expects incompatibleCounts GENEIDs to be a + # subset of rowData(se)$GENEID after filtering. + sample_names <- colnames(se) + empty_incompatible_counts <- function() { + cols <- c( + list(GENEID = character(0)), + stats::setNames(rep(list(numeric(0)), length(sample_names)), sample_names) + ) + data.table::as.data.table(cols) + } + + incompatible_counts <- S4Vectors::metadata(se)$incompatibleCounts + if (is.null(incompatible_counts)) { + incompatible_counts <- empty_incompatible_counts() + } else { + if (!"GENEID" %in% names(incompatible_counts)) { + incompatible_counts <- empty_incompatible_counts() + } else { + kept_genes <- unique(SummarizedExperiment::rowData(se)$GENEID) + incompatible_gene_ids <- incompatible_counts$GENEID + + rows_to_keep <- incompatible_gene_ids %in% kept_genes + incompatible_counts <- incompatible_counts[rows_to_keep, , drop = FALSE] + data.table::set( + incompatible_counts, + j = "GENEID", + value = incompatible_gene_ids[rows_to_keep] + ) + + for (sample_name in sample_names) { + if (!sample_name %in% names(incompatible_counts)) { + data.table::set( + incompatible_counts, + j = sample_name, + value = numeric(nrow(incompatible_counts)) + ) + } + } + + keep_cols <- c("GENEID", sample_names) + incompatible_counts <- incompatible_counts[, keep_cols, with = FALSE] + } + } + S4Vectors::metadata(se)$incompatibleCounts <- incompatible_counts + qc_stats$total_transcripts_after_filter <- nrow(se) qc_stats$total_genes_after_filter <- length(unique(SummarizedExperiment::rowData(se)$GENEID)) @@ -388,12 +444,7 @@ bambu_write_outputs <- function(se, gene_se, sample_df, args, qc_stats) { writeLines(capture.output(sessionInfo()), file.path(args$out_dir, "session_info.txt")) } -main_run_bambu <- function( - args, - analysis_fn = bambu::bambu, - prepare_annotations_fn = bambu::prepareAnnotations, - gene_expression_fn = bambu::transcriptToGeneExpression -) { +main_run_bambu <- function(args) { set.seed(42) # bambu's parallel worker code may rely on these being attached for generics # such as seqlengths(). @@ -404,7 +455,7 @@ main_run_bambu <- function( inputs <- bambu_resolve_inputs(args) args$threads <- bambu_effective_threads(args, length(inputs$bam_paths)) - annotation_obj <- prepare_annotations_fn(args$annotation) + annotation_obj <- bambu::prepareAnnotations(args$annotation) if (!is.null(args$ndr)) { message(sprintf("Using user-specified NDR = %.3f", args$ndr)) @@ -428,13 +479,15 @@ main_run_bambu <- function( message(sprintf("Running bambu with threads = %d", args$threads)) message("Running bambu...") - se <- do.call(analysis_fn, bambu_build_args(args, inputs$reads, annotation_obj)) + se <- do.call(bambu::bambu, bambu_build_args(args, inputs$reads, annotation_obj)) message("Bambu completed successfully") colnames(se) <- inputs$aliases filtered <- bambu_filter_transcripts(se) se <- filtered$se qc_stats <- filtered$qc_stats + gene_se <- bambu::transcriptToGeneExpression(se) + colnames(gene_se) <- inputs$aliases message( sprintf( "Filtering: keeping %d / %d transcripts", @@ -473,9 +526,6 @@ main_run_bambu <- function( qc_stats$transcripts_detected_per_sample <- as.list(detected_per_sample) qc_stats$median_transcripts_detected <- stats::median(detected_per_sample) - gene_se <- gene_expression_fn(se) - colnames(gene_se) <- inputs$aliases - bambu_write_outputs( se, gene_se, diff --git a/bin/workflow_glue_r/R/cli.R b/bin/workflow_glue_r/R/cli.R index abc8de4..d8bc8ac 100644 --- a/bin/workflow_glue_r/R/cli.R +++ b/bin/workflow_glue_r/R/cli.R @@ -1,4 +1,4 @@ -workflow_glue_r_components <- function(env = globalenv()) { +workflow_glue_r_components <- function(env = environment(workflow_glue_r_components)) { parser_suffix <- "_arg_parser" parser_names <- grep( paste0(parser_suffix, "$"), @@ -48,7 +48,7 @@ workflow_glue_r_usage <- function(components = workflow_glue_r_components()) { paste(lines, collapse = "\n") } -workflow_glue_r_cli <- function(argv = commandArgs(trailingOnly = TRUE), env = globalenv()) { +workflow_glue_r_cli <- function(argv = commandArgs(trailingOnly = TRUE), env = environment(workflow_glue_r_components)) { components <- workflow_glue_r_components(env = env) if (length(argv) < 1 || argv[[1]] %in% c("-h", "--help", "help")) { diff --git a/bin/workflow_glue_r/tests/testthat.R b/bin/workflow_glue_r/tests/testthat.R index f4a710b..a60e493 100644 --- a/bin/workflow_glue_r/tests/testthat.R +++ b/bin/workflow_glue_r/tests/testthat.R @@ -14,7 +14,10 @@ workflow_glue_r_load(pkg_dir) Sys.setenv( WORKFLOW_GLUE_R_PACKAGE_DIR = pkg_dir, WORKFLOW_GLUE_R_REPO_ROOT = repo_root, - TEST_DATA = Sys.getenv("TEST_DATA", unset = file.path(repo_root, "test_data")) + TEST_DATA = Sys.getenv("TEST_DATA", unset = file.path(repo_root, "test_data")), + TESTTHAT_PKG = "workflowGlueR", + TESTTHAT_PARALLEL = "true", + TESTTHAT_CPUS = "8" ) testthat_dir <- file.path(pkg_dir, "tests", "testthat") diff --git a/bin/workflow_glue_r/tests/testthat/test_bambu.R b/bin/workflow_glue_r/tests/testthat/test_bambu.R index 17b2d6a..417b316 100644 --- a/bin/workflow_glue_r/tests/testthat/test_bambu.R +++ b/bin/workflow_glue_r/tests/testthat/test_bambu.R @@ -231,94 +231,6 @@ testthat::test_that("transcriptome mode mapped to bambu args", { testthat::expect_equal(fixed$yieldSize, 250000L) }) -# End-to-end unit test with mocked bambu analysis function. -# Verifies --bams input resolution, sample sheet reordering, and discovery settings. -testthat::test_that("bams input with discovery mode", { - fixture_dir <- tempfile("bambu-discover-") - dir.create(fixture_dir) - bam_dir <- file.path(fixture_dir, "bams") - dir.create(bam_dir) - sample_a <- file.path(bam_dir, "sampleA.aligned.sorted.bam") - sample_b <- file.path(bam_dir, "sampleB.bam") - file.create(sample_b) - file.create(sample_a) - - sample_sheet <- file.path(fixture_dir, "sample_sheet.csv") - writeLines( - paste( - "alias,condition", - "sampleB,treated", - "sampleA,control", - sep = "\n" - ), - sample_sheet - ) - - captured <- new.env(parent = emptyenv()) - fake_analysis <- function( - reads, - annotations, - genome, - ncore, - discovery, - lowMemory, - NDR = NULL, - yieldSize = NULL, - ... - ) { - captured$reads <- reads - captured$annotations <- annotations - captured$genome <- genome - captured$ncore <- ncore - captured$discovery <- discovery - captured$NDR <- NDR - captured$low_memory <- lowMemory - captured$arg_yield_size <- yieldSize - captured$yield_size <- Rsamtools::yieldSize(reads[[1]]) - - base_se <- make_test_tx_se(sample_names = c("sampleA", "sampleB")) - SummarizedExperiment::SummarizedExperiment( - assays = SummarizedExperiment::assays(base_se), - rowRanges = make_test_bambu_row_ranges(fixture_dir) - ) - } - - argv <- list( - annotation = "annotation.gtf", - genome = "genome.fa", - out_dir = file.path(fixture_dir, "out"), - bams = paste(c(sample_a, sample_b), collapse = ","), - aliases = "sampleA,sampleB", - sample_sheet = sample_sheet, - transcriptome_mode = "discover", - ndr = 0.25, - threads = 2 - ) - - argv <- workflow_glue_r_normalise_args(argv, bambu_arg_spec()) - result <- suppressMessages(main_run_bambu( - argv, - analysis_fn = fake_analysis, - prepare_annotations_fn = function(annotation) { - captured$annotation_path <- annotation - structure(list(path = annotation), class = "mockAnnotation") - }, - gene_expression_fn = function(se) make_test_gene_se(sample_names = colnames(se)) - )) - - testthat::expect_equal(captured$annotation_path, "annotation.gtf") - testthat::expect_equal(captured$genome, "genome.fa") - testthat::expect_equal(captured$ncore, 1L) - testthat::expect_true(captured$discovery) - testthat::expect_equal(captured$NDR, 0.25) - testthat::expect_true(captured$low_memory) - testthat::expect_equal(captured$arg_yield_size, 250000L) - testthat::expect_equal(captured$yield_size, 250000) - testthat::expect_equal(result$sample_df$alias, c("sampleA", "sampleB")) - testthat::expect_equal(result$sample_df$condition, c("control", "treated")) - testthat::expect_true(file.exists(file.path(argv$out_dir, "bambu_qc_stats.json"))) -}) - ### # Output serialization # @@ -465,6 +377,41 @@ testthat::test_that("QC stats match filtered results", { testthat::expect_equal(result$qc_stats$transcripts_filtered, 3) }) +# Contract test: after transcript filtering, the filtered object must remain +# acceptable input for gene-level aggregation. +testthat::test_that("bambu_filter_transcripts contract with transcriptToGeneExpression", { + fixture_dir <- tempfile("bambu-filter-contract-") + dir.create(fixture_dir) + + base_se <- make_test_tx_se(sample_names = "sampleA") + se <- SummarizedExperiment::SummarizedExperiment( + assays = SummarizedExperiment::assays(base_se), + rowRanges = make_test_bambu_row_ranges(fixture_dir) + ) + + full_length <- matrix( + c(5, 4, 0, 0), + nrow = nrow(se), + ncol = ncol(se), + dimnames = dimnames(SummarizedExperiment::assays(se)$counts) + ) + SummarizedExperiment::assays(se, withDimnames = FALSE)[["fullLengthCounts"]] <- full_length + + S4Vectors::metadata(se)$incompatibleCounts <- data.table::data.table( + GENEID = "gene2", + sampleA = 7L + ) + + filtered <- bambu_filter_transcripts(se) + testthat::expect_true(all(SummarizedExperiment::rowData(filtered$se)$GENEID == "gene1")) + testthat::expect_equal( + unique(S4Vectors::metadata(filtered$se)$incompatibleCounts$GENEID), + character(0) + ) + + testthat::expect_error(bambu::transcriptToGeneExpression(filtered$se), NA) +}) + ### # CLI integration tests #