From 7b94eb2ca0407524364970526aaabfd769a64066 Mon Sep 17 00:00:00 2001 From: Sarah Griffiths Date: Fri, 22 May 2026 09:38:33 +0000 Subject: [PATCH] Annotation quant issues test [CW-7261] --- .gitlab-ci.yml | 9 ++ bin/workflow_glue_r/R/bambu.R | 43 ++++++- .../tests/testthat/test_bambu.R | 111 +++++++++++++++++- 3 files changed, 159 insertions(+), 4 deletions(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index c196ad6..47cb658 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -67,6 +67,7 @@ docker-run: - MATRIX_NAME: [ "int_discover_dna", "int_fixed_rna", "int_de_control_vs_control", "smoke_discover", "smoke_fixed", "smoke_direct_rna", "smoke_de", + "mouse_de_0countquant" ] rules: # NOTE As we're overriding the rules block for the included docker-run @@ -170,3 +171,11 @@ docker-run: test -f ${CI_PROJECT_NAME}/de_analysis/condition_treated_vs_control/results_dge.tsv && test -f ${CI_PROJECT_NAME}/de_analysis/condition_treated_vs_control/results_dtu_transcript.tsv && [ "$(find ${CI_PROJECT_NAME}/samples -type f -name 'gene_counts.tsv' | wc -l)" -eq 4 ] + # Tests a common error when there are 0 annotation counts for a chunk + - if: $MATRIX_NAME == "mouse_de_0countquant" + variables: + NF_BEFORE_SCRIPT: "mkdir -p ${CI_PROJECT_NAME}/data/ && wget https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-transcriptomes/mouse_subset_test.tar.gz -O ${CI_PROJECT_NAME}/data/mouse_subset_test.tar.gz && tar -xzvf ${CI_PROJECT_NAME}/data/mouse_subset_test.tar.gz -C ${CI_PROJECT_NAME}/data/" + NF_WORKFLOW_OPTS: "--fastq ${CI_PROJECT_NAME}/data/mouse_subset_test/samples \ + --de_analysis --ref_genome ${CI_PROJECT_NAME}/data/mouse_subset_test/mouse_subset.fa \ + --ref_annotation ${CI_PROJECT_NAME}/data/mouse_subset_test/mouse_subset.gtf.gz \ + --direct_rna --sample_sheet ${CI_PROJECT_NAME}/data/mouse_subset_test/sample_sheet.csv" \ No newline at end of file diff --git a/bin/workflow_glue_r/R/bambu.R b/bin/workflow_glue_r/R/bambu.R index 86f4788..8c22472 100644 --- a/bin/workflow_glue_r/R/bambu.R +++ b/bin/workflow_glue_r/R/bambu.R @@ -342,6 +342,17 @@ bambu_chunk_rc_files <- function(rc_files, aliases, sample_df) { }, seqnames, chunk_ids) } + +bambu_known_quant_edge_error_kind <- function(msg) { + if (grepl("Error in filter(., (uniqueStartLengthQuery <= primarySecondaryDistStartEnd", msg, fixed = TRUE)) { + return("uniqueStartLengthQuery_filter") + } + if (grepl("eqClassById` with `y$eqClassById` due to incompatible types.", msg, fixed = TRUE)) { + return("eqClassById_incompatible_types") + } + NA_character_ +} + bambu_write_discovery_outputs <- function(out_dir, rc_files, discovered_annotations, chunk_bundles, sample_df) { saveRDS(rc_files, file.path(out_dir, "bambu_rcfiles.rds")) saveRDS(discovered_annotations, file.path(out_dir, "bambu_discovered_annotations.rds")) @@ -440,6 +451,7 @@ bambu_run_discover_mode <- function(args, analysis_fn, prepare_annotations_fn, b ) } + bambu_run_quant_mode <- function(args, analysis_fn) { chunk_bundle <- readRDS(args$chunk_rds) discovered_annotations <- readRDS(args$discovered_annotation_rds) @@ -470,7 +482,11 @@ bambu_run_quant_mode <- function(args, analysis_fn) { sample_names = sample_names ) } else { - se <- bambu_call_analysis( + + # We attempted pre-quant validation (see CW-7261) but we do not seem to be able + # to reliably eliminate these chunk-specific empty/typing edge cases from the discover side. + se <- tryCatch( + bambu_call_analysis( analysis_fn, bambu_build_args( args, @@ -479,7 +495,28 @@ bambu_run_quant_mode <- function(args, analysis_fn) { discovery = FALSE, quant = TRUE ) - ) + ), + error = function(e) { + msg <- conditionMessage(e) + error_kind <- bambu_known_quant_edge_error_kind(msg) + if (is.na(error_kind)) { + stop(e) + } + warning( + sprintf( + "Skipping bambu quant for chunk '%s' (%s): known bambu chunk edge case (%s).", + chunk_bundle$chunk_id, + chunk_bundle$seqname, + error_kind + ), + call. = FALSE + ) + bambu_empty_quant_se( + discovered_annotations = discovered_annotations, + sample_names = sample_names + ) + } + ) } if (length(sample_names) > 0 && ncol(se) == length(sample_names)) { @@ -860,7 +897,7 @@ bambu_add_library_qc_stats <- function(se, qc_stats = list()) { sprintf( paste0( "Large library size variation detected (%.1fx difference).\n", - " Min: %d, Max: %d reads.\n", + " Min: %.0f, Max: %.0f reads.\n", " CPM normalization may not be appropriate for such variation." ), lib_size_ratio, diff --git a/bin/workflow_glue_r/tests/testthat/test_bambu.R b/bin/workflow_glue_r/tests/testthat/test_bambu.R index b765698..47edf4c 100644 --- a/bin/workflow_glue_r/tests/testthat/test_bambu.R +++ b/bin/workflow_glue_r/tests/testthat/test_bambu.R @@ -563,7 +563,7 @@ testthat::test_that("quant mode skips chunks with no discovered annotations on t suppressMessages(main_run_bambu( args, analysis_fn = fake_analysis - )), + )) , "contains no transcripts on this seqname" ) @@ -589,6 +589,115 @@ testthat::test_that("quant mode skips chunks with no discovered annotations on t testthat::expect_equal(qc$total_genes_before_filter, 0) }) +testthat::test_that("quant mode catches known uniqueStartLengthQuery/min-Inf edge case and writes empty outputs", { + fixture_dir <- tempfile("bambu-quant-known-edge-") + dir.create(fixture_dir) + + chunk_bundle <- list( + chunk_id = "chr1", + seqname = "chr1", + aliases = c("sampleA"), + sample_df = data.frame(alias = c("sampleA"), stringsAsFactors = FALSE), + rc_files = list(sampleA = make_test_tx_se(sample_names = "sampleA")), + annotation_tx_count = 1L + ) + chunk_rds <- file.path(fixture_dir, "chr1.rds") + saveRDS(chunk_bundle, chunk_rds) + + discovered_annotation_rds <- file.path(fixture_dir, "annotations.rds") + saveRDS(make_test_bambu_row_ranges(fixture_dir), discovered_annotation_rds) + + analysis_called <- FALSE + fake_analysis <- function(...) { + analysis_called <<- TRUE + stop( + paste( + "Error in filter(., (uniqueStartLengthQuery <= primarySecondaryDistStartEnd & : In argument: `==...`.", + "Caused by warning in `min()`: no non-missing arguments to min; returning Inf" + ), + call. = FALSE + ) + } + + args <- workflow_glue_r_normalise_args( + list( + mode = "quant", + genome = "genome.fa", + out_dir = file.path(fixture_dir, "out"), + chunk_rds = chunk_rds, + discovered_annotation_rds = discovered_annotation_rds, + transcriptome_mode = "discover", + ndr = NULL, + threads = 2 + ), + bambu_arg_spec() + ) + + testthat::expect_warning( + suppressMessages(main_run_bambu(args, analysis_fn = fake_analysis)), + "known bambu chunk edge case" + ) + + testthat::expect_true(analysis_called) + + se <- readRDS(file.path(args$out_dir, "bambu_transcripts.rds")) + testthat::expect_equal(nrow(se), 0) + testthat::expect_equal(colnames(se), c("sampleA")) +}) + +testthat::test_that("quant mode catches known eqClassById incompatible-type edge case and writes empty outputs", { + fixture_dir <- tempfile("bambu-quant-known-eqclass-edge-") + dir.create(fixture_dir) + + chunk_bundle <- list( + chunk_id = "chr2", + seqname = "chr2", + aliases = c("sampleA"), + sample_df = data.frame(alias = c("sampleA"), stringsAsFactors = FALSE), + rc_files = list(sampleA = make_test_tx_se(sample_names = "sampleA")), + annotation_tx_count = 1L + ) + chunk_rds <- file.path(fixture_dir, "chr2.rds") + saveRDS(chunk_bundle, chunk_rds) + + discovered_annotation_rds <- file.path(fixture_dir, "annotations.rds") + saveRDS(make_test_bambu_row_ranges(fixture_dir), discovered_annotation_rds) + + analysis_called <- FALSE + fake_analysis <- function(...) { + analysis_called <<- TRUE + stop( + "Can't join `x$eqClassById` with `y$eqClassById` due to incompatible types.", + call. = FALSE + ) + } + + args <- workflow_glue_r_normalise_args( + list( + mode = "quant", + genome = "genome.fa", + out_dir = file.path(fixture_dir, "out"), + chunk_rds = chunk_rds, + discovered_annotation_rds = discovered_annotation_rds, + transcriptome_mode = "discover", + ndr = NULL, + threads = 2 + ), + bambu_arg_spec() + ) + + testthat::expect_warning( + suppressMessages(main_run_bambu(args, analysis_fn = fake_analysis)), + "known bambu chunk edge case" + ) + + testthat::expect_true(analysis_called) + + se <- readRDS(file.path(args$out_dir, "bambu_transcripts.rds")) + testthat::expect_equal(nrow(se), 0) + testthat::expect_equal(colnames(se), c("sampleA")) +}) + testthat::test_that("empty mode writes valid empty outputs including bambu rds files", { fixture_dir <- tempfile("bambu-empty-mode-") dir.create(fixture_dir)