diff --git a/README.md b/README.md index 4838372..77eac09 100644 --- a/README.md +++ b/README.md @@ -229,6 +229,9 @@ All aligned samples are analysed together with `bambu` to produce the primary cohort transcriptome, transcript counts, gene counts, and the `RDS` objects used for downstream differential analysis. This shared model is the main cohort-level result and is published under `cohort/`. +Before writing outputs, transcript filtering removes only transcripts with zero +total transcript counts across samples; it does not use `fullLengthCounts` for +this quantification filter. ### 5. Independent per-sample transcriptomes diff --git a/bin/workflow_glue_r/R/bambu.R b/bin/workflow_glue_r/R/bambu.R index 8cb87b0..9b4b9a5 100644 --- a/bin/workflow_glue_r/R/bambu.R +++ b/bin/workflow_glue_r/R/bambu.R @@ -194,14 +194,13 @@ bambu_effective_threads <- function(args, bam_count) { bambu_filter_transcripts <- function(se) { assays <- SummarizedExperiment::assays(se) counts_mat <- assays$counts - full_length_mat <- assays$fullLengthCounts 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.") + if (is.null(counts_mat)) { + stop("counts assay not found in bambu output.") } gene_ids <- row_data$GENEID @@ -211,11 +210,7 @@ bambu_filter_transcripts <- function(se) { samples = ncol(se) ) - if (is.null(full_length_mat)) { - keep_idx <- rowSums(counts_mat) > 0 - } else { - keep_idx <- rowSums(full_length_mat) > 0 - } + keep_idx <- rowSums(counts_mat) > 0 if (!any(keep_idx)) { stop( "All transcripts have zero counts after filtering. ", diff --git a/bin/workflow_glue_r/tests/testthat/test_bambu.R b/bin/workflow_glue_r/tests/testthat/test_bambu.R index 417b316..4620b7e 100644 --- a/bin/workflow_glue_r/tests/testthat/test_bambu.R +++ b/bin/workflow_glue_r/tests/testthat/test_bambu.R @@ -307,9 +307,10 @@ testthat::test_that("bambu outputs written correctly", { }) # Transcript filtering edge cases -testthat::test_that("zero full-length count transcripts filtered", { +testthat::test_that("filters on counts when fullLengthCounts disagrees", { se <- make_test_tx_se(sample_names = c("sampleA", "sampleB")) - # Set tx1 to have full-length counts, tx2 to have none (will be filtered) + # tx2 has counts but zero full-length counts: must be kept. + # tx3 has zero counts but non-zero full-length counts: must be filtered. counts <- SummarizedExperiment::assay(se, "counts") counts[1, ] <- c(10, 8) counts[2, ] <- c(5, 3) @@ -319,12 +320,14 @@ testthat::test_that("zero full-length count transcripts filtered", { full_length <- matrix(0, nrow = 4, ncol = 2) full_length[1, ] <- c(5, 4) + full_length[3, ] <- c(6, 6) SummarizedExperiment::assays(se, withDimnames = FALSE)[["fullLengthCounts"]] <- full_length result <- bambu_filter_transcripts(se) - testthat::expect_equal(nrow(result$se), 1) - testthat::expect_equal(result$qc_stats$transcripts_filtered, 3) + testthat::expect_equal(rownames(result$se), c("tx1", "tx2")) + testthat::expect_equal(nrow(result$se), 2) + testthat::expect_equal(result$qc_stats$transcripts_filtered, 2) }) testthat::test_that("filters on counts when no fullLengthCounts assay", { @@ -352,7 +355,6 @@ testthat::test_that("error when all transcripts filtered", { se <- make_test_tx_se(sample_names = c("sampleA", "sampleB")) # All transcripts have zero counts SummarizedExperiment::assay(se, "counts", withDimnames = FALSE) <- matrix(0, nrow = 4, ncol = 2) - SummarizedExperiment::assays(se, withDimnames = FALSE)[["fullLengthCounts"]] <- matrix(0, nrow = 4, ncol = 2) testthat::expect_error( bambu_filter_transcripts(se), @@ -389,8 +391,15 @@ testthat::test_that("bambu_filter_transcripts contract with transcriptToGeneExpr rowRanges = make_test_bambu_row_ranges(fixture_dir) ) + counts <- SummarizedExperiment::assays(se)$counts + counts["tx1", ] <- 10 + counts["tx2", ] <- 8 + counts["tx3", ] <- 0 + counts["tx4", ] <- 0 + SummarizedExperiment::assay(se, "counts", withDimnames = FALSE) <- counts + full_length <- matrix( - c(5, 4, 0, 0), + c(5, 4, 6, 0), nrow = nrow(se), ncol = ncol(se), dimnames = dimnames(SummarizedExperiment::assays(se)$counts) diff --git a/docs/06_pipeline_overview.md b/docs/06_pipeline_overview.md index 669ddaa..5cde13e 100644 --- a/docs/06_pipeline_overview.md +++ b/docs/06_pipeline_overview.md @@ -74,6 +74,9 @@ All aligned samples are analysed together with `bambu` to produce the primary cohort transcriptome, transcript counts, gene counts, and the `RDS` objects used for downstream differential analysis. This shared model is the main cohort-level result and is published under `cohort/`. +Before writing outputs, transcript filtering removes only transcripts with zero +total transcript counts across samples; it does not use `fullLengthCounts` for +this quantification filter. ### 5. Independent per-sample transcriptomes