Merge Devel pre v4 to sc_spatial_wrapper

R/bambu-assignDist.R

@@ -2,11 +2,11 @@
 #' @inheritParams bambu
 #' @import data.table
 #' @noRd
-assignReadClasstoTranscripts <- function(readClassList, annotations, isoreParameters, 
-                                        verbose, sampleMetadata, demultiplexed,
+assignReadClasstoTranscripts <- function(readClassList, annotations,
+                                        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, 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-utilityExtend.R

@@ -112,7 +112,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])
@@ -224,7 +224,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)
 }
 
@@ -827,15 +827,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")
@@ -852,17 +851,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))
 
   temp <- c(metadata(extendedAnnotations)$lowConfidenceTranscripts[!toAdd], extendedAnnotations[toRemove])
   extendedAnnotations <- c(extendedAnnotations[!toRemove], metadata(extendedAnnotations)$lowConfidenceTranscripts[toAdd])
@@ -880,7 +872,7 @@ setNDR <- function(extendedAnnotations, NDR = NULL, includeRef = FALSE, prefix =
 
 #' 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
@@ -912,21 +904,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

@@ -3,14 +3,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()
@@ -20,21 +20,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 ",

R/bambu-processReads.R

@@ -14,7 +14,7 @@
 #' @noRd
 bambu.processReads <- function(reads, annotations, genomeSequence,
     readClass.outputDir=NULL, yieldSize=1000000, bpParameters, 
-    stranded=FALSE, verbose=FALSE, isoreParameters = setIsoreParameters(NULL),
+    stranded=FALSE, verbose=FALSE, discoveryParameters = setDiscoveryParameters(NULL),
     processByChromosome = FALSE, processByBam = TRUE, trackReads = trackReads, fusionMode = fusionMode, 
     demultiplexed = FALSE, cleanReads = FALSE, dedupUMI = FALSE, sampleNames = NULL, barcodesToFilter = NULL) {
     genomeSequence <- checkInputSequence(genomeSequence)
@@ -48,11 +48,11 @@ bambu.processReads <- function(reads, annotations, genomeSequence,
             names(reads)[seq_along(sampleNames)] <- sampleNames
         }
     }
-    min.readCount <- isoreParameters[["min.readCount"]]
-    fitReadClassModel <- isoreParameters[["fitReadClassModel"]]
-    defaultModels <- isoreParameters[["defaultModels"]]
-    returnModel <- isoreParameters[["returnModel"]]
-    min.exonOverlap <- isoreParameters[["min.exonOverlap"]]
+    min.readCount <- discoveryParameters[["min.readCount"]]
+    fitReadClassModel <- discoveryParameters[["fitReadClassModel"]]
+    defaultModels <- discoveryParameters[["defaultModels"]]
+    returnModel <- discoveryParameters[["returnModel"]]
+    min.exonOverlap <- discoveryParameters[["min.exonOverlap"]]
 
     if(processByBam){
         readClassList <- bplapply(seq_along(reads), function(i) {

R/bambu-quantify.R

@@ -2,9 +2,9 @@
 #' @inheritParams bambu
 #' @import data.table
 #' @noRd
-bambu.quantify <- function(readClassDt, countMatrix, incompatibleCountMatrix, txid.index, GENEIDs, emParameters, 
+bambu.quantify <- function(readClassDt, countMatrix, incompatibleCountMatrix, txid.index, GENEIDs, emParameters,
                            trackReads = FALSE, returnDistTable = FALSE,
-                           verbose = FALSE, isoreParameters = setIsoreParameters(NULL)) {
+                           verbose = FALSE) {
     start.ptm <- proc.time()
     readClassDt$nobs = countMatrix[readClassDt$eqClass.match]
     readClassDt$nobs[is.na(readClassDt$nobs)] = 0