diff --git a/NEWS b/NEWS index bdeb7127..5e8f8123 100644 --- a/NEWS +++ b/NEWS @@ -1,3 +1,8 @@ +bambu v3.13.2 (change date: 2026-08-13) +========================================== +Minor fixes: + * Reduced peak memory in combineSplicedTranscriptModels by storing the per-sample start, end and readCount values in long form rather than as one column trio per sample. Output is unchanged; the previous implementation is retained and used whenever the sparse path cannot guarantee an identical result. + bambu v3.2.0 (change date: 2023-Apr-26) ========================================== Minor fixes: diff --git a/R/bambu-extendAnnotations-utilityCombine.R b/R/bambu-extendAnnotations-utilityCombine.R index fdbd389a..6fe0189d 100644 --- a/R/bambu-extendAnnotations-utilityCombine.R +++ b/R/bambu-extendAnnotations-utilityCombine.R @@ -31,38 +31,237 @@ isore.combineTranscriptCandidates <- function(readClassList, #' combine spliced transcript models +#' +#' Dispatches to a sparse reduction, falling back to the original dense +#' implementation when the sparse path cannot reproduce it exactly (see +#' combineSplicedTranscriptModelsSparse). Both return the same 12 columns, +#' in the same order, with the same row order. #' @noRd -combineSplicedTranscriptModels <- function(readClassList, bpParameters, - min.readCount, min.readFractionByGene, min.txScore.multiExon, +combineSplicedTranscriptModels <- function(readClassList, bpParameters, + min.readCount, min.readFractionByGene, min.txScore.multiExon, min.txScore.singleExon, verbose){ bpParameters$progressbar = FALSE options(scipen = 999) #maintain numeric basepair locations not sci.notfi. start.ptm <- proc.time() n_sample <- length(readClassList) - nGroups = max(ceiling(n_sample/10),min(bpworkers(bpParameters), + nGroups = max(ceiling(n_sample/10),min(bpworkers(bpParameters), round(n_sample/2))) indexList <- sample(rep(seq_len(nGroups), length.out=n_sample)) indexList <- splitAsList(seq_len(n_sample), indexList) + combinedFeatureTibble <- NULL + ## With fewer than 3 samples nGroups is 1, so the dense path performs no + ## full_join at all and its NSample* columns stay logical while maxTxScore + ## keeps NAs. Reproducing that column-type quirk is not worth it; defer. + if (n_sample >= 3) + combinedFeatureTibble <- combineSplicedTranscriptModelsSparse( + readClassList, indexList, bpParameters, min.readCount, + min.readFractionByGene, min.txScore.multiExon, + min.txScore.singleExon) + if (is.null(combinedFeatureTibble)) + combinedFeatureTibble <- combineSplicedTranscriptModelsDense( + readClassList, indexList, bpParameters, min.readCount, + min.readFractionByGene, min.txScore.multiExon, + min.txScore.singleExon) + end.ptm <- proc.time() + if (verbose) message("combing spliced feature tibble objects across all ", + "samples in ", round((end.ptm - start.ptm)[3] / 60, 1)," mins.") + return(combinedFeatureTibble) +} + +#' combine spliced transcript models by joining per-sample columns +#' +#' The original implementation: builds one wide table carrying a start, end and +#' readCount column per sample, then reduces it row-wise. +#' @noRd +combineSplicedTranscriptModelsDense <- function(readClassList, indexList, + bpParameters, min.readCount, min.readFractionByGene, + min.txScore.multiExon, min.txScore.singleExon){ combinedFeatureTibbleList <- bplapply(seq_along(indexList), function(g){ indexVec <- indexList[[g]] return(sequentialCombineFeatureTibble(readClassList[indexVec], - indexVec, intraGroup = TRUE, - min.readCount = min.readCount, - min.readFractionByGene = min.readFractionByGene, + indexVec, intraGroup = TRUE, + min.readCount = min.readCount, + min.readFractionByGene = min.readFractionByGene, min.txScore.multiExon = min.txScore.multiExon, min.txScore.singleExon = min.txScore.singleExon)) }, BPPARAM = bpParameters) - combinedFeatureTibble <- - sequentialCombineFeatureTibble(combinedFeatureTibbleList, - indexList = NULL, intraGroup = FALSE) - combinedFeatureTibble <- updateStartEndReadCount(combinedFeatureTibble) - end.ptm <- proc.time() - if (verbose) message("combing spliced feature tibble objects across all ", - "samples in ", round((end.ptm - start.ptm)[3] / 60, 1)," mins.") - return(combinedFeatureTibble) + combinedFeatureTibble <- + sequentialCombineFeatureTibble(combinedFeatureTibbleList, + indexList = NULL, intraGroup = FALSE) + return(updateStartEndReadCount(combinedFeatureTibble)) +} + +#' combine spliced transcript models without materialising per-sample columns +#' +#' The dense path builds a table whose rows are distinct intron chains and whose +#' columns are 3n+9 for n samples. A read class is observed in only a few +#' samples, so that table is mostly NA, and the join fold holds a second copy of +#' it while the first is live. +#' +#' This stores only the populated (chain, sample) cells in long form and reduces +#' them by chain, so memory scales with observations rather than rows x samples. +#' +#' Returns NULL when it cannot guarantee an identical result, so the caller can +#' fall back to the dense implementation. +#' @noRd +combineSplicedTranscriptModelsSparse <- function(readClassList, indexList, + bpParameters, min.readCount, min.readFractionByGene, + min.txScore.multiExon, min.txScore.singleExon){ + res <- bplapply(seq_along(indexList), function(g){ + indexVec <- as.integer(indexList[[g]]) + return(sparseGroupFeatures(readClassList[indexVec], indexVec, + min.readCount = min.readCount, + min.readFractionByGene = min.readFractionByGene, + min.txScore.multiExon = min.txScore.multiExon, + min.txScore.singleExon = min.txScore.singleExon)) + }, BPPARAM = bpParameters) + if (any(vapply(res, is.null, logical(1)))) return(NULL) + ## Global first-appearance key ids, over the same group-major traversal the + ## dense path uses, so row order is preserved. + allk <- rbindlist(lapply(res, `[[`, "gk"), idcol = "g") + for (i in seq_along(res)) res[[i]]$gk <- NULL + allk[, gi := .GRP, by = c("intronStarts", "intronEnds", "chr", "strand")] + nkey <- max(allk$gi) + keyDT <- allk[!duplicated(gi), .(intronStarts, intronEnds, chr, strand)] + if (nrow(keyDT) != nkey) return(NULL) + agg <- allk[, .(nsrc = sum(nsrc), nsrp = sum(nsrp), nstx = sum(nstx), + mts = max(mts), mtsnf = max(mtsnf)), by = gi] + if (!identical(agg$gi, seq_len(nkey))) return(NULL) + cellCount <- vapply(seq_along(res), function(i) nrow(res[[i]]$cells), 0L) + gvec <- allk$g + giVec <- allk$gi + rm(allk) + GI <- integer(sum(cellCount)); ST <- integer(sum(cellCount)) + EN <- integer(sum(cellCount)); RC <- integer(sum(cellCount)) + pos <- 1L + for (i in seq_along(res)){ + cl <- res[[i]]$cells + keyMap <- giVec[gvec == i] # group-local id -> global id + if (nrow(cl)){ + sl <- pos:(pos + nrow(cl) - 1L) + GI[sl] <- keyMap[cl$li] + ST[sl] <- cl$start; EN[sl] <- cl$end; RC[sl] <- cl$readCount + pos <- pos + nrow(cl) + } + res[[i]]$cells <- NULL + } + rm(res, gvec, giVec) + if (anyNA(ST) || anyNA(EN) || anyNA(RC)) return(NULL) + ord <- order(GI, method = "radix") + GI <- GI[ord]; ST <- ST[ord]; EN <- EN[ord]; RC <- RC[ord] + rm(ord) + reduced <- reduceSparseCells(GI, ST, EN, RC, nkey) + if (is.null(reduced)) return(NULL) + ## A missing score is stored as -Inf while reducing so max() stays valid. + mts <- agg$mts; mts[is.infinite(mts) & mts < 0] <- NA_real_ + mtsnf <- agg$mtsnf; mtsnf[is.infinite(mtsnf) & mtsnf < 0] <- NA_real_ + ## pmax() on a logical returns double, so the dense path's NSample* columns + ## are double once any join has happened. Match that. + return(data.table(start = reduced$start, end = reduced$end, + readCount = reduced$readCount, intronStarts = keyDT$intronStarts, + intronEnds = keyDT$intronEnds, chr = keyDT$chr, + strand = keyDT$strand, maxTxScore = as.numeric(mts), + maxTxScore.noFit = as.numeric(mtsnf), + NSampleReadCount = as.numeric(agg$nsrc), + NSampleReadProp = as.numeric(agg$nsrp), + NSampleTxScore = as.numeric(agg$nstx))) +} + +#' summarise one group of samples into a key table and its populated cells +#' @noRd +sparseGroupFeatures <- function(readClassList, indexVec, min.readCount, + min.readFractionByGene, min.txScore.multiExon, min.txScore.singleExon){ + tabs <- vector("list", length(indexVec)) + for (s in seq_along(indexVec)){ + featureTibble <- extractFeaturesFromReadClassSE( + readClassSe = readClassList[[s]], sample_id = indexVec[s], + min.readCount = min.readCount, + min.readFractionByGene = min.readFractionByGene, + min.txScore.multiExon = min.txScore.multiExon, + min.txScore.singleExon = min.txScore.singleExon) + setDT(featureTibble) + ## One row per key per sample is what makes per-key counters equal to + ## per-row counters; a duplicate would fan out in the dense join. + if (anyDuplicated(featureTibble, + by = c("intronStarts", "intronEnds", "chr", "strand"))) + return(NULL) + tabs[[s]] <- featureTibble + } + big <- rbindlist(tabs) + rm(tabs) + big[, li := .GRP, by = c("intronStarts", "intronEnds", "chr", "strand")] + gk <- big[, .(intronStarts = intronStarts[1L], + intronEnds = intronEnds[1L], chr = chr[1L], strand = strand[1L], + nsrc = sum(NSampleReadCount), nsrp = sum(NSampleReadProp), + nstx = sum(NSampleTxScore, na.rm = TRUE), # NA counts as 0, as pmax0NA + mts = if (all(is.na(maxTxScore))) -Inf else + max(maxTxScore, na.rm = TRUE), + mtsnf = if (all(is.na(maxTxScore.noFit))) -Inf else + max(maxTxScore.noFit, na.rm = TRUE)), by = li] + if (!identical(gk$li, seq_len(nrow(gk)))) return(NULL) + gk[, li := NULL] + return(list(gk = gk, cells = big[, .(li, start, end, readCount)])) +} + +#' reduce long-form cells to one start, end and readCount per key +#' @noRd +reduceSparseCells <- function(GI, ST, EN, RC, nkey, cellChunk = 2e7){ + readCountSum <- numeric(nkey) + summed <- rowsum(as.numeric(RC), GI, reorder = FALSE) + readCountSum[as.integer(rownames(summed))] <- summed[, 1L] + rm(summed) + if (any(readCountSum > .Machine$integer.max)) return(NULL) + startOut <- rep(Inf, nkey) + endOut <- rep(Inf, nkey) + ## Chunk on key boundaries so no key is split across calls. + brk <- c(0L, which(GI[-1L] != GI[-length(GI)]), length(GI)) + step <- max(1L, as.integer(cellChunk)) + b0 <- 1L + while (b0 <= length(brk) - 1L){ + b1 <- b0 + while (b1 < length(brk) - 1L && (brk[b1 + 2L] - brk[b0]) <= step) + b1 <- b1 + 1L + idx <- (brk[b0] + 1L):brk[b1 + 1L] + keys <- GI[idx] + loK <- keys[1L]; hiK <- keys[length(keys)] + startOut[loK:hiK] <- upperMedianByGroup(keys - loK + 1L, ST[idx], + RC[idx], hiK - loK + 1L) + endOut[loK:hiK] <- upperMedianByGroup(keys - loK + 1L, EN[idx], + RC[idx], hiK - loK + 1L) + b0 <- b1 + 1L + } + if (any(is.infinite(startOut)) || any(is.infinite(endOut))) return(NULL) + return(list(start = as.integer(startOut), end = as.integer(endOut), + readCount = as.integer(readCountSum))) +} + +#' readCount-weighted median per group, without expanding the values +#' +#' readCountWeightedMedian() repeats each value by its read count, takes the +#' type 7 median and snaps up to the nearest observed value. Because the weights +#' are integers that is exactly the element at floor(N/2)+1 of the expanded +#' sorted vector, which can be located from the cumulative weights alone. +#' @noRd +upperMedianByGroup <- function(groupIndex, values, weights, ngroup){ + out <- rep(Inf, ngroup) + if (!length(groupIndex)) return(out) + ord <- order(groupIndex, values, method = "radix") + groupIndex <- groupIndex[ord] + values <- values[ord] + cumWeight <- cumsum(as.numeric(weights[ord])) + starts <- c(1L, which(groupIndex[-1L] != + groupIndex[-length(groupIndex)]) + 1L) + priorWeight <- c(0, cumWeight)[starts] + lastIdx <- c(starts[-1L] - 1L, length(cumWeight)) + groupTotal <- cumWeight[lastIdx] - priorWeight + target <- priorWeight + floor(groupTotal / 2) + 1 + hit <- findInterval(target - 0.5, cumWeight) + 1L + nonEmpty <- groupTotal > 0 + out[groupIndex[starts][nonEmpty]] <- values[hit[nonEmpty]] + return(out) } -#' Sequentially combine feature tibbles +#' Sequentially combine feature tibbles #' @noRd sequentialCombineFeatureTibble <- function(readClassList, indexList,intraGroup,min.readCount,min.readFractionByGene, diff --git a/R/globals.R b/R/globals.R index 962098d3..50b512b4 100644 --- a/R/globals.R +++ b/R/globals.R @@ -45,6 +45,7 @@ if (getRversion() >= "2.15.1") { "maxTxScore.noFit.combined","maxTxScore.noFit.new", "txScore.noFit","n.obs","nObs_list","K_list","txids_list", "anyEqual","txScore.noFit","txidTemp", - "subjectHits.y", "txNumberFiltered" + "subjectHits.y", "txNumberFiltered", + "li", "gi", "nsrc", "nsrp", "nstx", "mts", "mtsnf" )) } diff --git a/tests/testthat/test_combineTranscriptModels.R b/tests/testthat/test_combineTranscriptModels.R new file mode 100644 index 00000000..4506eb60 --- /dev/null +++ b/tests/testthat/test_combineTranscriptModels.R @@ -0,0 +1,100 @@ +context("Combine transcript models across samples") + +## Builds a minimal read class SE carrying exactly the columns +## extractFeaturesFromReadClassSE() consumes. Synthetic rather than a stored +## fixture so the sample set can be varied per test, and so the intron chains +## overlap across samples, which is what the combine has to reduce over. +makeReadClassSe <- function(nClass, seed, nShared = 6L) { + set.seed(seed) + ## The first nShared chains are the same in every sample; the rest are + ## sample specific. That gives both keys seen once and keys seen many times. + chainId <- c(seq_len(nShared), seed * 1000L + seq_len(nClass - nShared)) + starts <- 1000L * chainId + exonsByRc <- GenomicRanges::GRangesList(lapply(seq_len(nClass), function(i) { + GenomicRanges::GRanges(seqnames = "chr9", + ranges = IRanges::IRanges( + start = c(starts[i], starts[i] + 400L), + end = c(starts[i] + 200L, starts[i] + 600L)), + strand = "+") + })) + rd <- S4Vectors::DataFrame( + chr.rc = factor(rep("chr9", nClass), levels = "chr9"), + strand.rc = factor(rep("+", nClass), levels = c("+", "-", "*")), + intronStarts = as.character(starts + 201L), + intronEnds = as.character(starts + 399L), + confidenceType = rep("highConfidenceJunctionReads", nClass), + readCount = sample(seq_len(50L), nClass, replace = TRUE), + geneReadProp = runif(nClass), + txScore = runif(nClass), + txScore.noFit = runif(nClass), + numExons = rep(2L, nClass)) + ## A missing score has to survive the reduction as NA, so seed some. + rd$txScore[sample(nClass, max(1L, nClass %/% 10L))] <- NA_real_ + se <- SummarizedExperiment::SummarizedExperiment( + assays = list(counts = matrix(rd$readCount, ncol = 1)), + rowRanges = exonsByRc) + SummarizedExperiment::rowData(se) <- rd + se +} + +test_that("upperMedianByGroup reproduces readCountWeightedMedian", { + set.seed(1) + for (trial in seq_len(50)) { + ngroup <- sample(seq_len(6), 1) + n <- sample(seq(ngroup, 40), 1) + groupIndex <- sort(sample(seq_len(ngroup), n, replace = TRUE)) + groupIndex <- match(groupIndex, sort(unique(groupIndex))) + values <- sample(seq_len(500), n, replace = TRUE) + weights <- sample(seq_len(20), n, replace = TRUE) + got <- upperMedianByGroup(groupIndex, values, weights, + max(groupIndex)) + want <- vapply(seq_len(max(groupIndex)), function(g) { + keep <- groupIndex == g + dt <- data.table(v = values[keep], w = weights[keep]) + readCountWeightedMedian(dt, "v", "w") + }, numeric(1)) + expect_equal(got, want) + } +}) + +test_that("the sparse combine is identical to the dense combine", { + ## Three samples is the smallest set that takes the sparse path. + for (nSample in c(3L, 5L)) { + readClassList <- lapply(seq_len(nSample), function(i) + makeReadClassSe(20L + i, seed = i)) + bpParameters <- BiocParallel::SerialParam() + seed <- 42L + ## combineSplicedTranscriptModels draws from the RNG itself, so both + ## paths have to start from the same seed to see the same grouping. + set.seed(seed) + sparse <- combineSplicedTranscriptModels(readClassList, bpParameters, + min.readCount = 2, min.readFractionByGene = 0.05, + min.txScore.multiExon = 0, min.txScore.singleExon = 1, + verbose = FALSE) + set.seed(seed) + indexList <- sample(rep(seq_len(max(ceiling(nSample / 10), + min(BiocParallel::bpworkers(bpParameters), + round(nSample / 2)))), length.out = nSample)) + indexList <- splitAsList(seq_len(nSample), indexList) + dense <- combineSplicedTranscriptModelsDense(readClassList, indexList, + bpParameters, min.readCount = 2, min.readFractionByGene = 0.05, + min.txScore.multiExon = 0, min.txScore.singleExon = 1) + expect_identical(setDT(sparse), setDT(dense)) + } +}) + +test_that("the sparse combine defers when a sample repeats an intron chain", { + readClassList <- lapply(seq_len(3), function(i) + makeReadClassSe(20L, seed = i)) + ## Duplicating a key inside one sample would fan out in the dense join, so + ## the sparse path must decline rather than return a different answer. + rd <- SummarizedExperiment::rowData(readClassList[[1]]) + rd$intronStarts[2] <- rd$intronStarts[1] + rd$intronEnds[2] <- rd$intronEnds[1] + SummarizedExperiment::rowData(readClassList[[1]]) <- rd + indexList <- splitAsList(seq_len(3), c(1L, 1L, 2L)) + expect_null(combineSplicedTranscriptModelsSparse(readClassList, indexList, + BiocParallel::SerialParam(), min.readCount = 2, + min.readFractionByGene = 0.05, min.txScore.multiExon = 0, + min.txScore.singleExon = 1)) +})