diff --git a/CHANGELOG.md b/CHANGELOG.md index 3364fb6..d0bebea 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,16 +7,18 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [v2.0.1] -This patch release of `wf-transcriptomes` handles an additional quantification failure edge case that was not observed before release, and fixes an issue encountered during joint discovery for many samples. +This patch release of `wf-transcriptomes` handles an additional quantification failure edge case that was not observed before release, and fixes issues encountered by users during joint discovery when providing many samples. Users of wf-transcriptomes v2.0.0 who have encountered issues during discovery and quantification should adopt this release. ### Fixed - "Error in full_join" encountered during `runPerSampleBambuQuant` when all read classes have no compatible transcript assignment. An empty quant table is correctly emitted instead. - "unable to find an inherited method for function 'rowData'" encountered during `runJointBambuDiscover` when providing many samples. The workflow now correctly handles data spilled to disk by bambu discover. +- Fatal memory error (exit 137) encountered during `collateBambuQuant` when providing many samples with a large reference. Collation now uses a two-pass approach to process the per-chunk RDS results to avoid exhausting memory limits. - Volcano plot class counts incorrect when `log2FoldChange` or `padj` columns contained NA values. - Adjusted p-values below 0.001 in the volcano selection table are now shown in scientific notation instead of being rounded to 0.000. - IGV track not correctly loading in EPI2ME Desktop when a sample consists of a single input BAM. + ## [v2.0.0] This release refreshes `wf-transcriptomes` around a new reference-guided transcriptomics workflow built on `bambu`, with `SQANTI3` transcript classification and QC, `DESeq2` for differential gene expression, `DEXSeq` for differential transcript usage, and per-sample modified base summarisation with `modkit` when modification tags are present in aligned BAMs. diff --git a/bin/workflow_glue_r/R/bambu.R b/bin/workflow_glue_r/R/bambu.R index 81f3c03..0db468e 100644 --- a/bin/workflow_glue_r/R/bambu.R +++ b/bin/workflow_glue_r/R/bambu.R @@ -959,6 +959,55 @@ bambu_sum_incompatible_counts <- function(tx_ses) { data.table::as.data.table(combined) } +bambu_add_incompatible_counts <- function(combined, incompatible_counts, sample_names) { + if (is.null(incompatible_counts)) { + return(combined) + } + + current_sample_cols <- setdiff(colnames(incompatible_counts), c("GENEID", "TXNAME")) + if (!identical(current_sample_cols, sample_names)) { + stop( + "Chunk quantification outputs have mismatched incompatible count columns.", + call. = FALSE + ) + } + + if (is.null(combined)) { + combined <- data.frame( + GENEID = as.character(incompatible_counts$GENEID), + stringsAsFactors = FALSE + ) + for (sample_col in sample_names) { + values <- incompatible_counts[[sample_col]] + values[is.na(values)] <- 0 + combined[[sample_col]] <- values + } + return(combined) + } + + idx <- match(as.character(incompatible_counts$GENEID), combined$GENEID) + if (anyNA(idx)) { + new_gene_ids <- as.character(incompatible_counts$GENEID[is.na(idx)]) + new_rows <- data.frame( + GENEID = new_gene_ids, + stringsAsFactors = FALSE + ) + for (sample_col in sample_names) { + new_rows[[sample_col]] <- numeric(length(new_gene_ids)) + } + combined <- rbind(combined, new_rows) + idx <- match(as.character(incompatible_counts$GENEID), combined$GENEID) + } + + for (sample_col in sample_names) { + values <- incompatible_counts[[sample_col]] + values[is.na(values)] <- 0 + combined[[sample_col]][idx] <- combined[[sample_col]][idx] + values + } + + combined +} + bambu_unique_row_ranges <- function(tx_ses, tx_names) { combined_row_ranges <- do.call(c, lapply(tx_ses, SummarizedExperiment::rowRanges)) unique_row_ranges <- combined_row_ranges[!duplicated(names(combined_row_ranges))] @@ -1030,6 +1079,124 @@ bambu_combine_transcript_chunks <- function(tx_ses) { ) } +bambu_read_chunk_transcript_se <- function(chunk_dir) { + readRDS(file.path(chunk_dir, "bambu_transcripts.rds")) +} + +bambu_validate_chunk_shape <- function(se, assay_names, sample_names) { + if (!identical(SummarizedExperiment::assayNames(se), assay_names)) { + stop("Chunk quantification outputs have mismatched assay sets.", call. = FALSE) + } + if (!identical(colnames(se), sample_names)) { + stop("Chunk quantification outputs have mismatched sample columns.", call. = FALSE) + } +} + +bambu_collect_chunk_metadata <- function(chunk_dirs) { + # First pass: discover the global row order/metadata before allocating + # combined assay matrices, without retaining all chunk SEs in memory. + first_se <- bambu_read_chunk_transcript_se(chunk_dirs[[1]]) + assay_names <- SummarizedExperiment::assayNames(first_se) + sample_names <- colnames(first_se) + col_data <- SummarizedExperiment::colData(first_se) + object_metadata <- S4Vectors::metadata(first_se) + empty_row_ranges <- SummarizedExperiment::rowRanges(first_se)[0] + + tx_names <- character(0) + combined_row_ranges <- empty_row_ranges + incompatible_counts <- NULL + + for (chunk_dir in chunk_dirs) { + se <- bambu_read_chunk_transcript_se(chunk_dir) + bambu_validate_chunk_shape(se, assay_names, sample_names) + + row_names <- rownames(se) + new_idx <- !row_names %in% tx_names + if (any(new_idx)) { + tx_names <- c(tx_names, row_names[new_idx]) + combined_row_ranges <- c( + combined_row_ranges, + SummarizedExperiment::rowRanges(se)[new_idx] + ) + } + + chunk_incompatible <- bambu_normalise_incompatible_counts( + S4Vectors::metadata(se)$incompatibleCounts, + sample_names = sample_names + ) + incompatible_counts <- bambu_add_incompatible_counts( + incompatible_counts, + chunk_incompatible, + sample_names + ) + } + + row_ranges <- combined_row_ranges[match(tx_names, names(combined_row_ranges))] + if (is.null(incompatible_counts)) { + incompatible_counts <- bambu_empty_incompatible_counts(sample_names) + } + object_metadata$incompatibleCounts <- data.table::as.data.table(incompatible_counts) + + list( + assay_names = assay_names, + sample_names = sample_names, + col_data = col_data, + object_metadata = object_metadata, + tx_names = tx_names, + row_ranges = row_ranges + ) +} + +bambu_fill_combined_assays_from_chunks <- function(chunk_dirs, chunk_metadata) { + combined_assays <- stats::setNames(lapply(chunk_metadata$assay_names, function(assay_name) { + matrix( + 0, + nrow = length(chunk_metadata$tx_names), + ncol = length(chunk_metadata$sample_names), + dimnames = list(chunk_metadata$tx_names, chunk_metadata$sample_names) + ) + }), chunk_metadata$assay_names) + + for (chunk_dir in chunk_dirs) { + se <- bambu_read_chunk_transcript_se(chunk_dir) + for (assay_name in chunk_metadata$assay_names) { + assay_mat <- SummarizedExperiment::assay(se, assay_name) + collapsed <- rowsum(assay_mat, group = rownames(se), reorder = FALSE) + idx <- match(rownames(collapsed), chunk_metadata$tx_names) + combined_assays[[assay_name]][idx, ] <- + combined_assays[[assay_name]][idx, , drop = FALSE] + collapsed + } + } + + if ("counts" %in% names(combined_assays) && "CPM" %in% names(combined_assays)) { + counts_mat <- combined_assays[["counts"]] + combined_assays[["CPM"]] <- t(t(counts_mat) / pmax(colSums(counts_mat), 1)) * 1e6 + } + + combined_assays +} + +# Directory-backed equivalent of bambu_combine_transcript_chunks(); production +# collation uses this to preserve combiner semantics without retaining all chunks. +bambu_combine_transcript_chunk_dirs <- function(chunk_dirs) { + if (length(chunk_dirs) < 1) { + stop("No chunk quantification results were provided for collation.", call. = FALSE) + } + if (length(chunk_dirs) == 1) { + return(bambu_read_chunk_transcript_se(chunk_dirs[[1]])) + } + + chunk_metadata <- bambu_collect_chunk_metadata(chunk_dirs) + combined_assays <- bambu_fill_combined_assays_from_chunks(chunk_dirs, chunk_metadata) + + SummarizedExperiment::SummarizedExperiment( + assays = combined_assays, + rowRanges = chunk_metadata$row_ranges, + colData = chunk_metadata$col_data, + metadata = chunk_metadata$object_metadata + ) +} + bambu_collate_chunk_outputs <- function( chunk_dirs, out_dir, @@ -1043,10 +1210,8 @@ bambu_collate_chunk_outputs <- function( stop("No chunk quantification directories were provided for collation.", call. = FALSE) } - tx_ses <- lapply(chunk_dirs, function(chunk_dir) { - readRDS(file.path(chunk_dir, "bambu_transcripts.rds")) - }) - raw_se <- bambu_combine_transcript_chunks(tx_ses) + # Avoid retaining every chunk SE in memory; large cohorts can otherwise exceed 25 GB here. + raw_se <- bambu_combine_transcript_chunk_dirs(chunk_dirs) sample_df <- workflow_glue_r_read_sample_sheet(file.path(chunk_dirs[[1]], "samples.csv")) diff --git a/bin/workflow_glue_r/tests/testthat/test_bambu.R b/bin/workflow_glue_r/tests/testthat/test_bambu.R index 8e04b41..48ae57c 100644 --- a/bin/workflow_glue_r/tests/testthat/test_bambu.R +++ b/bin/workflow_glue_r/tests/testthat/test_bambu.R @@ -1089,6 +1089,65 @@ testthat::test_that("chunk combiner sums duplicate transcript rows across chunks testthat::expect_equal(incompatible$sampleA, c(10, 20)) }) +testthat::test_that("streaming chunk directory combiner matches in-memory combiner", { + sample_df <- data.frame(alias = "sampleA", stringsAsFactors = FALSE) + base_se <- make_test_tx_se(sample_names = sample_df$alias) + fixture_dir <- tempfile("bambu-combine-dirs-") + dir.create(fixture_dir) + tx_se <- SummarizedExperiment::SummarizedExperiment( + assays = SummarizedExperiment::assays(base_se), + rowRanges = make_test_bambu_row_ranges(fixture_dir) + ) + + chunk1 <- tx_se + chunk2 <- tx_se + counts1 <- SummarizedExperiment::assay(chunk1, "counts") + counts2 <- SummarizedExperiment::assay(chunk2, "counts") + counts1[3:4, ] <- 0 + counts2[1:2, ] <- 0 + SummarizedExperiment::assay(chunk1, "counts", withDimnames = FALSE) <- counts1 + SummarizedExperiment::assay(chunk2, "counts", withDimnames = FALSE) <- counts2 + SummarizedExperiment::assay(chunk1, "CPM", withDimnames = FALSE) <- + t(t(counts1) / pmax(colSums(counts1), 1)) * 1e6 + SummarizedExperiment::assay(chunk2, "CPM", withDimnames = FALSE) <- + t(t(counts2) / pmax(colSums(counts2), 1)) * 1e6 + + S4Vectors::metadata(chunk1)$incompatibleCounts <- data.frame( + GENEID = c("gene1", "gene2"), + `01` = c(10, 0), + stringsAsFactors = FALSE + ) + S4Vectors::metadata(chunk2)$incompatibleCounts <- data.frame( + GENEID = c("gene1", "gene2"), + `01` = c(0, 20), + stringsAsFactors = FALSE + ) + + chunk_dirs <- file.path(fixture_dir, c("chunk1", "chunk2")) + dir.create(chunk_dirs[[1]]) + dir.create(chunk_dirs[[2]]) + saveRDS(chunk1, file.path(chunk_dirs[[1]], "bambu_transcripts.rds")) + saveRDS(chunk2, file.path(chunk_dirs[[2]], "bambu_transcripts.rds")) + + in_memory <- bambu_combine_transcript_chunks(list(chunk1, chunk2)) + streaming <- bambu_combine_transcript_chunk_dirs(chunk_dirs) + + testthat::expect_equal(rownames(streaming), rownames(in_memory)) + testthat::expect_equal(SummarizedExperiment::rowData(streaming), SummarizedExperiment::rowData(in_memory)) + testthat::expect_equal( + SummarizedExperiment::assay(streaming, "counts"), + SummarizedExperiment::assay(in_memory, "counts") + ) + testthat::expect_equal( + SummarizedExperiment::assay(streaming, "CPM"), + SummarizedExperiment::assay(in_memory, "CPM") + ) + testthat::expect_equal( + S4Vectors::metadata(streaming)$incompatibleCounts, + S4Vectors::metadata(in_memory)$incompatibleCounts + ) +}) + # Transcript filtering edge cases testthat::test_that("filters on counts when fullLengthCounts disagrees", { se <- make_test_tx_se(sample_names = c("sampleA", "sampleB")) diff --git a/modules/local/bambu_chunked.nf b/modules/local/bambu_chunked.nf index aa49e47..3b22155 100644 --- a/modules/local/bambu_chunked.nf +++ b/modules/local/bambu_chunked.nf @@ -85,7 +85,7 @@ process bambuEmpty { process collateBambuQuant { label "wf_transcriptomes" cpus 1 - memory "16 GB" + memory "8 GB" input: tuple val(meta), path(chunk_dirs, stageAs: "chunks/*") path annotation, stageAs: "annotation/*"