Merge branch 'CW-7261-b' into 'dev'

Annotation quant issues test [CW-7261]

Closes CW-7261

See merge request epi2melabs/workflows/wf-transcriptomes!284

.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"

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,

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)