Move min.exonDistance, min.primarySecondaryDist, min.primarySecondaryDistStartEnd1/2 from opt.discovery into a new opt.rcAssignment parameter with its own setRcAssignmentParameters().

Author: SuiYue-2308

R/bambu-assignDist.R

@@ -2,11 +2,11 @@
 #' @inheritParams bambu
 #' @import data.table
 #' @noRd
-assignReadClasstoTranscripts <- function(readClassList, annotations, isoreParameters, 
+assignReadClasstoTranscripts <- function(readClassList, annotations, isoreParameters, rcAssignmentParameters,
                                         verbose, sampleMetadata, demultiplexed,
                                         returnDistTable = FALSE, trackReads = TRUE) {
     if (is.character(readClassList)) readClassList <- readRDS(file = readClassList)
-    metadata(readClassList)$readClassDist <- calculateDistTable(readClassList, annotations, isoreParameters, verbose, returnDistTable)
+    metadata(readClassList)$readClassDist <- calculateDistTable(readClassList, annotations, isoreParameters, rcAssignmentParameters, verbose, returnDistTable)
     readClassList <- splitReadClassFiles(readClassList)
     readClassDt <- genEquiRCs(metadata(readClassList)$readClassDist, annotations, verbose) 
     readClassDt$eqClass.match = match(readClassDt$eqClassById,metadata(readClassList)$eqClassById)

R/bambu-extendAnnotations.R

@@ -3,7 +3,7 @@
 #' @inheritParams bambu
 #' @noRd
 bambu.extendAnnotations <- function(readClassList, annotations, NDR,
-    isoreParameters, stranded, bpParameters, fusionMode = FALSE, verbose = FALSE) {
+    isoreParameters, rcAssignmentParameters, stranded, bpParameters, fusionMode = FALSE, verbose = FALSE) {
     start.ptm_all <- proc.time()
     combinedTxCandidates <- isore.combineTranscriptCandidates(readClassList,
         stranded, ## stranded used for unspliced reduce  
@@ -23,12 +23,12 @@ bambu.extendAnnotations <- function(readClassList, annotations, NDR,
         remove.subsetTx = isoreParameters[["remove.subsetTx"]],
         min.sampleNumber = isoreParameters[["min.sampleNumber"]],
         NDR = NDR,
-        min.exonDistance = isoreParameters[["min.exonDistance"]],
+        min.exonDistance = rcAssignmentParameters[["min.exonDistance"]],
         min.exonOverlap = isoreParameters[["min.exonOverlap"]],
-        min.primarySecondaryDist = 
-        isoreParameters[['min.primarySecondaryDist']], 
-        min.primarySecondaryDistStartEnd = 
-        isoreParameters[['min.primarySecondaryDistStartEnd1']],
+        min.primarySecondaryDist =
+        rcAssignmentParameters[['min.primarySecondaryDist']],
+        min.primarySecondaryDistStartEnd =
+        rcAssignmentParameters[['min.primarySecondaryDistStartEnd1']],
         min.readFractionByEqClass =  
         isoreParameters[['min.readFractionByEqClass']],
         fusionMode = fusionMode,

R/bambu.R

@@ -141,7 +141,7 @@
 #'     genome = fa.file,  discovery = TRUE, quant = TRUE)
 #' @export
 bambu <- function(reads, annotations = NULL, genome = NULL, NDR = NULL,
-    mode = NULL, opt.discovery = NULL, opt.em = NULL, rcOutDir = NULL, discovery = TRUE, 
+    mode = NULL, opt.discovery = NULL, opt.rcAssignment = NULL, opt.em = NULL, rcOutDir = NULL, discovery = TRUE,
     assignDist = TRUE, quant = TRUE, stranded = FALSE,  ncore = 1, yieldSize = NULL,  
     trackReads = FALSE, returnDistTable = FALSE, lowMemory = FALSE, sampleData = NULL,
     fusionMode = FALSE, verbose = FALSE, demultiplexed = FALSE, quantData = NULL,
@@ -183,7 +183,7 @@ bambu <- function(reads, annotations = NULL, genome = NULL, NDR = NULL,
     isoreParameters <- setIsoreParameters(isoreParameters = opt.discovery)
     #below line is to be compatible with earlier version of running bambu
     if(!is.null(isoreParameters$max.txNDR)) NDR = isoreParameters$max.txNDR
-    
+    rcAssignmentParameters <- setRcAssignmentParameters(rcAssignmentParameters = opt.rcAssignment)
     emParameters <- setEmParameters(emParameters = opt.em)
     bpParameters <- setBiocParallelParameters(reads, ncore, verbose, demultiplexed)
 	xgb.set.config(nthread = 1)
@@ -223,7 +223,7 @@ bambu <- function(reads, annotations = NULL, genome = NULL, NDR = NULL,
         if (discovery) {
             message("--- Start extending annotations ---")
             extendedAnnotations <- bambu.extendAnnotations(readClassList, annotations, NDR,
-                                                           isoreParameters, stranded, bpParameters, fusionMode, verbose)
+                                                           isoreParameters, rcAssignmentParameters, stranded, bpParameters, fusionMode, verbose)
             metadata(extendedAnnotations)$warnings = warnings
 
             #### cluster based transcript discovery
@@ -242,9 +242,10 @@ bambu <- function(reads, annotations = NULL, genome = NULL, NDR = NULL,
             quantData <- bplapply(seq_along(readClassList), function(i){
               assignReadClasstoTranscripts(
                 readClassList = readClassList[[i]],
-                annotations = annotations, 
-                isoreParameters = isoreParameters, 
-                verbose = verbose, 
+                annotations = annotations,
+                isoreParameters = isoreParameters,
+                rcAssignmentParameters = rcAssignmentParameters,
+                verbose = verbose,
                 # for bulk data, there is one sampleData (keep sampleData[1]), for single-cell, there is one per sample
                 sampleMetadata = if(length(sampleData) == 1) sampleData[1] else sampleData[i],
                 demultiplexed = demultiplexed, 

R/bambu_utilityFunctions.R

@@ -22,11 +22,7 @@ setIsoreParameters <- function(isoreParameters){
         min.readCount = 2,
         min.readFractionByGene = 0.05,
         min.sampleNumber = 1,
-        min.exonDistance = 35,
-        min.exonOverlap = 10, 
-        min.primarySecondaryDist = 5,
-        min.primarySecondaryDistStartEnd1 = 5, # for creating new annotations
-        min.primarySecondaryDistStartEnd2 = 5, # for read assignment
+        min.exonOverlap = 10,
         min.txScore.multiExon = 0,
         min.txScore.singleExon = 1,
         fitReadClassModel = TRUE,
@@ -41,6 +37,20 @@ setIsoreParameters <- function(isoreParameters){
 }
 
 
+#' setRcAssignmentParameters
+#' @noRd
+setRcAssignmentParameters <- function(rcAssignmentParameters){
+    rcAssignmentParameters.default <- list(
+        min.exonDistance = 35,
+        min.primarySecondaryDist = 5,
+        min.primarySecondaryDistStartEnd1 = 5, # for extending annotations
+        min.primarySecondaryDistStartEnd2 = 5) # for read class assignment
+    rcAssignmentParameters <-
+        updateParameters(rcAssignmentParameters, rcAssignmentParameters.default)
+    return(rcAssignmentParameters)
+}
+
+
 #' setEmParameters
 #' @noRd
 setEmParameters <- function(emParameters){
@@ -242,11 +252,11 @@ handleWarnings <- function(readClassList, verbose){
 }
 
 #' Calculate the dist table used for Bambu Quantification
-calculateDistTable <- function(readClassList, annotations, isoreParameters, verbose, returnDistTable){
+calculateDistTable <- function(readClassList, annotations, isoreParameters, rcAssignmentParameters, verbose, returnDistTable){
     readClassDist <- isore.estimateDistanceToAnnotations(readClassList, annotations,
-                                                            min.exonDistance = isoreParameters[["min.exonDistance"]],
-                                                            min.primarySecondaryDist = isoreParameters[['min.primarySecondaryDist']],
-                                                            min.primarySecondaryDistStartEnd = isoreParameters[['min.primarySecondaryDistStartEnd2']],
+                                                            min.exonDistance = rcAssignmentParameters[["min.exonDistance"]],
+                                                            min.primarySecondaryDist = rcAssignmentParameters[['min.primarySecondaryDist']],
+                                                            min.primarySecondaryDistStartEnd = rcAssignmentParameters[['min.primarySecondaryDistStartEnd2']],
                                                             verbose = verbose)
         metadata(readClassDist)$distTable <- modifyIncompatibleAssignment(metadata(readClassDist)$distTable)
         if(returnDistTable) metadata(readClassDist)$distTableOld <- metadata(readClassDist)$distTable