Repository navigation
[Feature]: speed up the updateStartEndReadCount function #595
Description
Feature Area
Performance / Memory Optimization
The Goal or Problem
updateStartEndReadCount is very slow on large datasets. This can be easily fixed
Current code:
updateStartEndReadCount <- function(combinedFeatureTibble){
setDT(combinedFeatureTibble)
combinedFeatureTibble[, rowID := .I]
colNames <- colnames(combinedFeatureTibble)
readCountCols <- sort(colNames[grep("^readCount", colNames)]) # to make sure it's ordered by sample name
startCols <- sort(colNames[grep("^start", colNames)])
endCols <- sort(colNames[grep("^end", colNames)])
startEndDt <- combinedFeatureTibble[,
.(start = readCountWeightedMedian(.SD,x,y),
end = readCountWeightedMedian(.SD,z,y),
readCount = sum(.SD[,y], na.rm = TRUE)),
by = rowID, env = I(list(x = startCols, y = readCountCols,z = endCols))]
combinedFeatureTibble <- startEndDt[combinedFeatureTibble[,.(intronStarts, intronEnds, chr, strand, maxTxScore,
maxTxScore.noFit, NSampleReadCount, NSampleReadProp,
NSampleTxScore, rowID)], on = "rowID"]
combinedFeatureTibble[, rowID := NULL]
return(combinedFeatureTibble)
}
#' Function to get median value without interpolation using certain column names
#' @noRd
readCountWeightedMedian <- function(dt, valuevar, timesvar){
sortVector <- rep(na.omit(unlist(dt[,..valuevar])),
times = as.integer(na.omit(unlist(dt[,..timesvar]))))
return(min(sortVector[sortVector>=quantile(sortVector, probs = 0.5)]))
}
The Proposal & Example Usage (Optional)
With the code below, run time is ~80x shorter.
My test data is an input to updateStartEndReadCount derived from my real data. You can download it at https://drive.google.com/file/d/1ek_6oPWqTvuTRBixSaUaZPMTyWZ_ZfiS/view?usp=sharing (92MB)
Ideally, such a parallel weighted median should be rewritten in C/C++ with Rcpp, but 50s is anyhow better than 1h.Alternatives Considered (Optional)
updateStartEndReadCount2 <- function(combinedFeatureTibble){
setDT(combinedFeatureTibble)
colNames <- colnames(combinedFeatureTibble)
readCountCols <- sort(colNames[grep("^readCount", colNames)]) # to make sure it's ordered by sample name
startCols <- sort(colNames[grep("^start", colNames)])
endCols <- sort(colNames[grep("^end", colNames)])
combinedFeatureTibble[, start:=apply(.SD[, c(startCols, readCountCols), with=FALSE], 1, function (x) {
matrixStats::weightedMedian(x=x[seq.int(from=1, to=length(x)/2)], w=x[seq.int(from=length(x)/2+1, to=length(x))], na.rm=TRUE, interpolate=FALSE, ties="min")
})]
combinedFeatureTibble[, end:=apply(.SD[, c(endCols, readCountCols), with=FALSE], 1, function (x) {
matrixStats::weightedMedian(x=x[seq.int(from=1, to=length(x)/2)], w=x[seq.int(from=length(x)/2+1, to=length(x))], na.rm=TRUE, interpolate=FALSE, ties="min")
})]
combinedFeatureTibble[, readCount:=rowSums(.SD, na.rm=TRUE), .SDcols=readCountCols]
combinedFeatureTibble[, c(startCols, readCountCols, endCols):=NULL]
return(combinedFeatureTibble)
}
RStudio output:
> library(data.table)
> cft <- fread("combinedFeatureTibble.tsv.gz")
> dim(cft)
[1] 3663960 15
> inputCFT1 <- data.table::copy(cft)
> inputCFT2 <- data.table::copy(cft)
> system.time(outCFT1 <- updateStartEndReadCount(inputCFT1))
Processed 3663960 groups out of 3663960. 100% done. Time elapsed: 4064s. ETA: 0s...
user system elapsed
4067.269 0.000 4065.019
> system.time(outCFT2 <- updateStartEndReadCount2(inputCFT2))
user system elapsed
51.130 0.000 51.109
> all.equal(outCFT1, outCFT2, ignore.col.order=TRUE)
[1] TRUE
> 4065/51
[1] 79.70588
Checklist
- I searched existing issues to avoid duplicates.
- I described a concrete use case and success criteria.
Activity
oleksii-nikolaienko commentedon Oct 12, 2026
should maybe add that the
updateStartEndReadCount(readCountWeightedMedianspecifically) takes ~75% of time of entire discovery step