Merge branch 'devel_pre_v4' into devel_pre_v4_jg

Resolve conflicts across 4 files by taking the devel_pre_v4 refactor
(isoreParameters -> opt.discovery/rcAssignmentParameters, demultiplexed
-> extractBarcodeUMI; removed sampleNames/cleanReads/barcodesToFilter
handling) and porting the applicable JG TODO review comments. Drop the
stale setIsoreParameters Pass-2 header and comments that referred to
code devel_pre_v4 deleted.

Also fix a pre-existing brace mismatch in checkInputs() introduced by
2047190 that prevented the package from parsing.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>

Author: jonathangoeke

NAMESPACE

@@ -1,6 +1,7 @@
 # Generated by roxygen2: do not edit by hand
 
 export(bambu)
+export(bambu.singlecell)
 export(plotBambu)
 export(prepareAnnotations)
 export(readFromGTF)

R/bambu-assignDist.R

@@ -6,11 +6,11 @@
 #' @inheritParams bambu
 #' @import data.table
 #' @noRd
-assignReadClasstoTranscripts <- function(readClassList, annotations, isoreParameters, 
-                                        verbose, sampleMetadata, demultiplexed,
+assignReadClasstoTranscripts <- function(readClassList, annotations, rcAssignmentParameters,
+                                        verbose, sampleMetadata, extractBarcodeUMI,
                                         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, rcAssignmentParameters, verbose, returnDistTable)
     readClassList <- splitReadClassFiles(readClassList)
     readClassDt <- genEquiRCs(metadata(readClassList)$readClassDist, annotations, verbose) 
     readClassDt$eqClass.match = match(readClassDt$eqClassById,metadata(readClassList)$eqClassById)
@@ -21,7 +21,7 @@ assignReadClasstoTranscripts <- function(readClassList, annotations, isoreParame
         mutate(aval = 1) %>%
         data.table()
     #return non-em counts
-    ColData <- generateColData(readClassList, sampleMetadata, demultiplexed)
+    ColData <- generateColData(readClassList, sampleMetadata, extractBarcodeUMI)
     quantData <- SummarizedExperiment(assays = SimpleList(
         counts = generateUniqueCounts(readClassDt, metadata(readClassList)$countMatrix, annotations)),
         rowRanges = annotations,

R/bambu-extendAnnotations-utilityExtend.R

@@ -124,7 +124,7 @@ filterTranscriptsByAnnotation <- function(rowDataCombined, annotationGrangesList
   } else if(is.null(NDR)) {
           NDR <- 0.5
   }
-  filterSet <- (rowDataCombined$NDR <= NDR | rowDataCombined$readClassType == "equal:compatible")
+  filterSet <- ((!is.na(rowDataCombined$NDR) & rowDataCombined$NDR <= NDR) | rowDataCombined$readClassType == "equal:compatible")
   lowConfidenceTranscripts <- combindRowDataWithRanges(
         rowDataCombined[!filterSet,], 
         exonRangesCombined[!filterSet])
@@ -248,7 +248,7 @@ calculateNDROnTranscripts <- function(combinedTranscripts, useTxScore = FALSE){
     } else {
         combinedTranscripts$NDR <- calculateNDR(combinedTranscripts$maxTxScore, equal)
     }
-    combinedTranscripts$NDR[combinedTranscripts$maxTxScore==-1] <- 1
+    combinedTranscripts$NDR[combinedTranscripts$maxTxScore==-1] <- NA
     return(combinedTranscripts)
 }
 
@@ -916,15 +916,14 @@ addGeneIdsToReadClassTable <- function(readClassTable, distTable,
 #' @description This function train a model for use on other data
 #' @param extendedAnnotations A GRangesList object produced from bambu(quant = FALSE) or rowRanges(se)
 #' @param NDR The maximum NDR for novel transcripts to be in extendedAnnotations (0-1). If not provided a recommended NDR is calculated.
-#' @param includeRef A boolean which if TRUE will also filter out reference annotations based on their NDR
 #' @param prefix A string which determines which transcripts are considered novel by bambu and will be filtered (by default = 'Bambu')
 #' @param baselineFDR a value between 0-1. Bambu uses this FDR on the trained model to recommend an equivilent NDR threshold to be used for the sample. By default, a baseline FDR of 0.1 is used. This does not impact the analysis if an NDR is set.
 #' @param defaultModels a bambu trained model object that bambu will use when fitReadClassModel==FALSE or the data is not suitable for training, defaults to the pretrained model in the bambu package
 #' Output - returns a similiar GRangesList object with entries swapped into or out of metadata(extendedAnnotations)$lowConfidenceTranscripts
 #' @details 
 #' @return extendedAnnotations with a new NDR threshold
 #' @export
-setNDR <- function(extendedAnnotations, NDR = NULL, includeRef = FALSE, prefix = 'Bambu', baselineFDR = 0.1, defaultModels2 = defaultModels){
+setNDR <- function(extendedAnnotations, NDR = NULL, prefix = 'Bambu', baselineFDR = 0.1, defaultModels2 = defaultModels){
     #Check to see if the annotations/gtf are dervived from Bambu
     if(is.null(mcols(extendedAnnotations)$NDR)){
         warning("Annotations were not extended by Bambu (or the wrong prefix was provided). NDR can not be set")
@@ -941,17 +940,10 @@ setNDR <- function(extendedAnnotations, NDR = NULL, includeRef = FALSE, prefix =
         message("Recommending a novel discovery rate (NDR) of: ", NDR)
     }
 
-    #If reference annotations should be filtered too (note that reference annotations with no read support arn't filtered)
-    if(includeRef){
-        toRemove <- (!is.na(mcols(extendedAnnotations)$NDR) & mcols(extendedAnnotations)$NDR > NDR)
-        toAdd <- !is.na(mcols(metadata(extendedAnnotations)$lowConfidenceTranscripts)$NDR) & 
-            mcols(metadata(extendedAnnotations)$lowConfidenceTranscripts)$NDR <= NDR  
-    } else {
-        toRemove <- (mcols(extendedAnnotations)$NDR > NDR & 
-            grepl(prefix, mcols(extendedAnnotations)$TXNAME))
-        toAdd <- (mcols(metadata(extendedAnnotations)$lowConfidenceTranscripts)$NDR <= NDR & 
-            grepl(prefix, mcols(metadata(extendedAnnotations)$lowConfidenceTranscripts)$TXNAME))     
-    }
+    toRemove <- (mcols(extendedAnnotations)$NDR > NDR &
+        grepl(prefix, mcols(extendedAnnotations)$TXNAME))
+    toAdd <- (mcols(metadata(extendedAnnotations)$lowConfidenceTranscripts)$NDR <= NDR &
+        grepl(prefix, mcols(metadata(extendedAnnotations)$lowConfidenceTranscripts)$TXNAME))
 
   # TODO: [POOR NAMING] temp stores annotations being swapped between low/high confidence sets; rename to swappedAnnotations or similar
   temp <- c(metadata(extendedAnnotations)$lowConfidenceTranscripts[!toAdd], extendedAnnotations[toRemove])
@@ -974,7 +966,7 @@ setNDR <- function(extendedAnnotations, NDR = NULL, includeRef = FALSE, prefix =
 # Call count: 1 call, 1 file
 #' Extend annotations by clusters
 #' @noRd
-isore.extendAnnotations.clusters <- function(readClassList, annotations, clusters, NDR, isoreParameters, stranded, bpParameters, fusionMode, verbose = FALSE){
+isore.extendAnnotations.clusters <- function(readClassList, annotations, clusters, NDR, discoveryParameters, stranded, bpParameters, fusionMode, verbose = FALSE){
     message("--- Start extending annotations for clusters ---")
     #if clustering is a csv, create a list with the barcodes for each cluster
     #csv must have two cols with heading barcode, cluster
@@ -992,7 +984,7 @@ isore.extendAnnotations.clusters <- function(readClassList, annotations, cluster
     for(i in seq_along(clusters)){
         print(names(clusters)[i])
         ###TODO need to account for the sample name here which is added to the barcode
-        index <- match(clusters[[i]],gsub('demultiplexed','',metadata(readClassList[[1]])$samples)) 
+        index <- match(clusters[[i]],gsub('demultiplexed','',metadata(readClassList[[1]])$samples))
         index <- index[!is.na(index)]
         print(length(index))
         if(length(index)<20) next
@@ -1006,21 +998,21 @@ isore.extendAnnotations.clusters <- function(readClassList, annotations, cluster
         rowData(rcf.filt)$startSD <- 0
         rowData(rcf.filt)$endSD <- 0
         rowData(rcf.filt)$readCount.posStrand <- 0
-        thresholdIndex <- which(rowData(rcf.filt)$readCount>=isoreParameters$min.readCount)
-        model <- trainBambu(rcf.filt, verbose = verbose, min.readCount = isoreParameters$min.readCount)
+        thresholdIndex <- which(rowData(rcf.filt)$readCount>=discoveryParameters$min.readCount)
+        model <- trainBambu(rcf.filt, verbose = verbose, min.readCount = discoveryParameters$min.readCount)
         txScore <- getTranscriptScore(rowData(rcf.filt)[thresholdIndex,], model,
                                 defaultModels)
         rowData(rcf.filt)$txScore <- rep(NA,nrow(rcf.filt))
         rowData(rcf.filt)$txScore[thresholdIndex] <- txScore
         #txScores = cbind(txScores, rowData(rcf.filt)$txScore)
         rcfs.clusters[[names(clusters)[i]]] <- rcf.filt
         annotations.clusters[[names(clusters)[i]]] <- bambu.extendAnnotations(list(rcf.filt), annotations, NDR,
-                                isoreParameters, stranded, bpParameters, fusionMode, verbose)
+                                discoveryParameters, stranded, bpParameters, fusionMode, verbose)
     }
     if(length(rcfs.clusters)>0){
         print("--- Merging all individual clusters ---")
         annotations.clusters[["merged"]] <- bambu.extendAnnotations(rcfs.clusters, annotations, NDR,
-            isoreParameters, stranded, bpParameters, fusionMode, verbose)
+            discoveryParameters, stranded, bpParameters, fusionMode, verbose)
     }
 
     return(annotations.clusters)

R/bambu-extendAnnotations.R

@@ -7,14 +7,14 @@
 #' @inheritParams bambu
 #' @noRd
 bambu.extendAnnotations <- function(readClassList, annotations, NDR,
-    isoreParameters, stranded, bpParameters, fusionMode = FALSE, verbose = FALSE) {
+    discoveryParameters, stranded, bpParameters, fusionMode = FALSE, verbose = FALSE) {
     start.ptm_all <- proc.time()
     combinedTxCandidates <- isore.combineTranscriptCandidates(readClassList,
         stranded, ## stranded used for unspliced reduce  
-        min.readCount = isoreParameters[["min.readCount"]], 
-        min.readFractionByGene = isoreParameters[["min.readFractionByGene"]],
-        min.txScore.multiExon = isoreParameters[["min.txScore.multiExon"]],
-        min.txScore.singleExon = isoreParameters[["min.txScore.singleExon"]],
+        min.readCount = discoveryParameters[["min.readCount"]], 
+        min.readFractionByGene = discoveryParameters[["min.readFractionByGene"]],
+        min.txScore.multiExon = discoveryParameters[["min.txScore.multiExon"]],
+        min.txScore.singleExon = discoveryParameters[["min.txScore.singleExon"]],
         bpParameters,
         verbose)
     end.ptm_all <- proc.time()
@@ -24,21 +24,21 @@ bambu.extendAnnotations <- function(readClassList, annotations, NDR,
     annotations <- isore.extendAnnotations(
         combinedTranscripts = combinedTxCandidates,
         annotationGrangesList = annotations,
-        remove.subsetTx = isoreParameters[["remove.subsetTx"]],
-        min.sampleNumber = isoreParameters[["min.sampleNumber"]],
+        remove.subsetTx = discoveryParameters[["remove.subsetTx"]],
+        min.sampleNumber = discoveryParameters[["min.sampleNumber"]],
         NDR = NDR,
-        min.exonDistance = isoreParameters[["min.exonDistance"]],
-        min.exonOverlap = isoreParameters[["min.exonOverlap"]],
+        min.exonDistance = discoveryParameters[["min.exonDistance"]],
+        min.exonOverlap = discoveryParameters[["min.exonOverlap"]],
         min.primarySecondaryDist = 
-        isoreParameters[['min.primarySecondaryDist']], 
+        discoveryParameters[['min.primarySecondaryDist']], 
         min.primarySecondaryDistStartEnd = 
-        isoreParameters[['min.primarySecondaryDistStartEnd1']],
+        discoveryParameters[['min.primarySecondaryDistStartEnd1']],
         min.readFractionByEqClass =  
-        isoreParameters[['min.readFractionByEqClass']],
+        discoveryParameters[['min.readFractionByEqClass']],
         fusionMode = fusionMode,
-        prefix = isoreParameters[["prefix"]],
-        baselineFDR = isoreParameters[["baselineFDR"]],
-        defaultModels = isoreParameters[["defaultModels"]],
+        prefix = discoveryParameters[["prefix"]],
+        baselineFDR = discoveryParameters[["baselineFDR"]],
+        defaultModels = discoveryParameters[["defaultModels"]],
         verbose = verbose)
     end.ptm_all <- proc.time()
     if (verbose) message("extend annotations in ",