Store combined transcript models in long form to bound peak memory

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:

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,

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"
     ))
 }

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))
+})