Update devel_pre_v4_beta from private release

DESCRIPTION

@@ -45,7 +45,8 @@ Suggests:
     purrr,
     apeglm,
     utils,
-    DEXSeq
+    DEXSeq,
+    Seurat
 Enhances: parallel
 SystemRequirements:
 biocViews: 

NAMESPACE

@@ -6,7 +6,7 @@ export(plotBambu)
 export(prepareAnnotations)
 export(readFromGTF)
 export(transcriptToGeneExpression)
-export(generateUniqueCountsSEFromQuantData)
+export(getUniqueCountsSe)
 export(writeBambuOutput)
 export(writeToGTF)
 export(writeAnnotationsToGTF)

R/bambu_utilityFunctions.R

@@ -81,10 +81,11 @@ updateParameters <- function(Parameters, Parameters.default) {
 #' @param reads path to BAM file(s)
 #' @param readClass.file path to readClass file(s)
 #' @param readClass.outputDir path to readClass output directory
+#' @param clusters NULL, a named vector (names are cell/spot ids in sampleName_barcode format, values are cluster labels), a data.frame with id (sampleName_barcode) and cluster columns, or a path to a .csv/.tsv/.txt file holding that data.frame
 #' @importFrom methods is
 #' @noRd
 checkInputs <- function(annotations, reads, readClass.outputDir, genomeSequence,
-                        discovery, sampleData, quantData){
+                        discovery, sampleData, quantData, clusters = NULL){
     # ===# Check annotation inputs #===#
     if (!is.null(annotations)) {
         if (is(annotations, "CompressedGRangesList")) {
@@ -159,6 +160,7 @@ checkInputs <- function(annotations, reads, readClass.outputDir, genomeSequence,
             use of Rsamtools for opening.")
     }
 
+    # ===# Check sampleData inputs #===#
     if(!is.null(sampleData)){
         if (!all(grepl("\\.(csv|tsv|txt)$", na.omit(sampleData), ignore.case = TRUE))){
             stop("Not all paths for sample metadata files are .csv/.tsv/.txt files")
@@ -173,6 +175,31 @@ checkInputs <- function(annotations, reads, readClass.outputDir, genomeSequence,
             )
         }
     }
+
+    # ===# Check clusters inputs #===#
+    if(!is.null(clusters)){
+        if(is.atomic(clusters) && !is.null(names(clusters))){ # named vector check
+            if(length(clusters) == 0){
+                stop("clusters named vector must not be empty")
+            }
+        } else if(is.data.frame(clusters)){ # data.frame object check
+            if(!all(c("id", "cluster") %in% colnames(clusters))){
+                stop("clusters data frame must have columns named 'id' and 'cluster'")
+            }
+        } else if(is.character(clusters) && length(clusters) == 1){ # path check
+            if(!file.exists(clusters)){
+                stop("clusters file does not exist: ", clusters)
+            }
+            if(!grepl("\\.(csv|tsv|txt)$", clusters, ignore.case = TRUE)){
+                stop("clusters must be a .csv/.tsv/.txt file")
+            }
+            if(!all(c("id", "cluster") %in% colnames(fread(clusters, nrows = 0)))){
+                stop("clusters file must be a data frame with columns named 'id' and 'cluster'")
+            }
+        } else {
+            stop("clusters must be a named vector, a data.frame with 'id' and 'cluster' columns, or a .csv/.tsv/.txt file path")
+        }
+    }
     return(annotations)
 }
 

R/transcriptToGeneExpression.R

@@ -24,7 +24,7 @@ transcriptToGeneExpression <- function(se) {
     }
     counts.total = colSums(counts)
     counts.total[counts.total==0] = 1
-    counts.CPM = counts/counts.total * 10^6
+    counts.CPM = counts %*% Diagonal(x = 1 / counts.total) * 10^6
 
     ## geneRanges
     exByGene <- reducedRangesByGenes(rowRanges(se))
@@ -50,20 +50,54 @@ transcriptToGeneExpression <- function(se) {
     return(seOutput)
 }
 
-#' Generate a SummarizedExperiment of unique counts from quantData
+#' @title Unique-count SummarizedExperiment from Bambu read-to-transcript assignments
 #' @description This function is intended to be used after the transcript
-#'   discovery and \code{assignDist} steps in \code{\link{bambu}}. It builds a
-#'   transcript-level SummarizedExperiment containing raw unique counts (reads
-#'   uniquely assigned to a single transcript) without EM estimation, which can
-#'   be passed directly to \code{\link{transcriptToGeneExpression}} to obtain
-#'   gene-level counts as uniqueCounts + nonuniqueCounts + incompatibleCounts.
-#' @param quantData a list of quantData objects produced by the assignDist step
-#' @param annotations a GRangesList of transcript annotations
-#' @return A SummarizedExperiment object with \code{assays$uniqueCounts},
-#'   \code{metadata$incompatibleCounts}, and \code{metadata$nonuniqueCounts}
+#' discovery and read-to-transcript assignment steps in \code{\link{bambu}} /
+#' \code{\link{bambu.singlecell}}. It generates a transcript-level SummarizedExperiment
+#' containing raw unique counts (reads uniquely assigned to a single transcript) without
+#' EM estimation. This function is useful for highly multiplexed, sparse data (such as
+#' single cell and spatial data) where the EM does not have sufficient information to
+#' provide accurate transcript expression estimates.
+#' @param quantData A list of \code{quantData} objects, one per sample, produced by the
+#' read-to-transcript assignment step of \code{\link{bambu.singlecell}}
+#' (\code{output = "quantData"}) or the equivalent \code{assignDist = TRUE} run of
+#' \code{\link{bambu}}.
+#' @param annotations A \code{GRangesList} of transcript annotations matching the ones used to
+#' produce \code{quantData}, typically the extended annotations from transcript discovery.
+#' @return A \code{SummarizedExperiment} with one row per transcript and one column per cell
+#' (or sample). Pass it to \code{\link{transcriptToGeneExpression}} to collapse the unique
+#' counts to the gene level. It contains:
+#' \describe{
+#'     \item{\code{assays(se)$counts}}{a sparse matrix of unique counts, i.e. reads uniquely
+#'     assigned to a single transcript.}
+#'     \item{\code{rowRanges(se)}}{the transcript \code{annotations} provided in the argument.}
+#'     \item{\code{colData(se)}}{per-cell (or per-sample) metadata carried over from
+#'     \code{quantData}, such as \code{id}, \code{sampleName}, and \code{barcode}.}
+#'     \item{\code{metadata(se)$incompatibleCounts}}{per-gene counts of reads not compatible
+#'     with any annotated transcript. \code{\link{transcriptToGeneExpression}} adds these back
+#'     into the gene-level counts to give more accurate gene expression estimates, so reads
+#'     that cannot be pinned to one transcript still count toward their gene.}
+#'     \item{\code{metadata(se)$nonuniqueCounts}}{per-gene counts of reads compatible with more
+#'     than one transcript (ambiguous assignments). Like \code{incompatibleCounts}, these are
+#'     added back into the gene-level counts by \code{\link{transcriptToGeneExpression}} for
+#'     more accurate gene expression estimates, and also indicate how many reads were
+#'     ambiguously assigned.}
+#'     \item{\code{metadata(se)$seType}}{a label identifying the \code{SummarizedExperiment} object as \code{"uniqueCounts"} type}
+#' }
+#' @seealso \code{\link{bambu.singlecell}} and \code{\link{bambu}} for producing \code{quantData};
+#' \code{\link{transcriptToGeneExpression}} to collapse the result to gene-level counts.
+#' @examples
+#' ## This works on both bulk and single-cell quantData;
+#' ## here we demonstrate with the single-cell test data.
+#' rds.dir <- system.file("extdata", "single_cell", package = "bambu")
+#' quantData <- readRDS(file.path(rds.dir,
+#'     "quantData_GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCell_multisample_chr9_1_1000000.rds"))
+#' extendedAnnotations <- readRDS(file.path(rds.dir,
+#'     "extendedAnnotations_GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCell_multisample_chr9_1_1000000.rds"))
+#' uniqueCountsSe <- getUniqueCountsSe(quantData, extendedAnnotations)
 #' @import data.table
 #' @export
-generateUniqueCountsSEFromQuantData <- function(quantData, annotations) {
+getUniqueCountsSe <- function(quantData, annotations) {
     uniqueCountsList <- lapply(quantData, function(x) {
         readClassDt <- getReadClassDt(x)
         x_filtered <- readClassDt %>% filter(!multi_align & !is.na(eqClass.match))