Merge branch 'cw-7337' into 'dev'
Stream collate RDS to reduce max mem [CW-7337] See merge request epi2melabs/workflows/wf-transcriptomes!338
This commit is contained in:
Sam Nicholls
commit
cca505faae
@ -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.
|
||||
|
||||
@ -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"))
|
||||
|
||||
|
||||
@ -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"))
|
||||
|
||||
@ -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/*"
|
||||
|
||||
Loading…
Reference in New Issue
Block a user