diff --git a/DESCRIPTION b/DESCRIPTION index 4430f2a5..2197253b 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -45,7 +45,8 @@ Suggests: purrr, apeglm, utils, - DEXSeq + DEXSeq, + Seurat Enhances: parallel SystemRequirements: biocViews: diff --git a/NAMESPACE b/NAMESPACE index e9c86d79..ccf252fe 100755 --- a/NAMESPACE +++ b/NAMESPACE @@ -6,7 +6,7 @@ export(plotBambu) export(prepareAnnotations) export(readFromGTF) export(transcriptToGeneExpression) -export(generateUniqueCountsSEFromQuantData) +export(getUniqueCountsSe) export(writeBambuOutput) export(writeToGTF) export(writeAnnotationsToGTF) diff --git a/R/bambu.R b/R/bambu.R index fb130080..91ea73fb 100644 --- a/R/bambu.R +++ b/R/bambu.R @@ -86,7 +86,9 @@ #' read class files will be saved. #' @param discovery A logical variable indicating whether annotations #' are to be extended. Defaults to TRUE -#' @param quant A logical variable indicating whether quantification will +#' @param assignDist A logical variable indicating whether read-class-to-transcript +#' assignment will be performed, defaults to TRUE +#' @param quant A logical variable indicating whether quantification will #' be performed. If false the output type will change. Defaults to TRUE #' @param stranded A boolean for strandedness, defaults to FALSE. #' @param ncore specifying number of cores used when parallel processing @@ -109,15 +111,9 @@ #' @param fusionMode A logical variable indicating whether run in fusion mode #' @param verbose A logical variable indicating whether processing messages will #' be printed. -#' @param opt.singlecell A list of single-cell specific parameters: -#' \describe{ -#' \item{extractBarcodeUMI}{Logical, whether to extract cell barcodes and -#' UMIs from BAM tags or read names. Defaults to FALSE} -#' \item{dedupUMI}{Logical, whether to perform UMI-based deduplication per -#' barcode. Defaults to FALSE} -#' \item{clusters}{A named list mapping cluster names to barcode vectors, -#' used for cluster-level quantification. Defaults to NULL} -#' } +#' @param opt.singlecell A list of parameters for bambu's single-cell module +#' (\code{extractBarcodeUMI}, \code{dedupUMI}, and \code{clusters}). See +#' \code{\link{bambu.singlecell}} for details. #' @details #' @return \code{bambu} will output different results depending on whether #' \emph{quant} mode is on. By default, \emph{quant} is set to TRUE, so @@ -198,7 +194,8 @@ bambu <- function(reads, annotations = NULL, genome = NULL, NDR = NULL, annotations <- checkInputs(annotations, reads, readClass.outputDir = rcOutDir, genomeSequence = genome, discovery = discovery, - sampleData = sampleData, quantData = quantData) + sampleData = sampleData, quantData = quantData, + clusters = opt.singlecell$clusters) } opt.discovery <- setDiscoveryParameters(discoveryParameters = opt.discovery) #below line is to be compatible with earlier version of running bambu @@ -285,37 +282,44 @@ bambu <- function(reads, annotations = NULL, genome = NULL, NDR = NULL, countsSeCompressed.all <- NULL ColNames <- c() colData.all <- list() + if (!is.null(clusters)) { + clusterDf <- if (is.atomic(clusters) && !is.null(names(clusters))) { + data.frame(id = names(clusters), cluster = as.character(clusters)) + } else if (is.data.frame(clusters)) { + clusters + } else { + fread(clusters) + } + } for(i in seq_along(quantData)){ quantData_i <- quantData[[i]] - #load in the barcode clustering from file if provided - # single-cell mode: iter is an integer vector of column indices, one per barcode - iter <- seq_len(nrow(getSampleData(quantData_i))) if(!is.null(clusters)){ - if(class(clusters[[i]])!="CompressedCharacterList"){ # !is.list(clusters) is FALSE for CompressedCharacterList - clusterMaps <- NULL - for(j in seq_along(getSampleData(quantData_i)$sampleName)){ #load in a file per sample name provided - clusterMap <- fread(clusters[[j]], header = FALSE, - data.table = FALSE) - # read.table(clusters[[j]], - # sep = ifelse(grepl(".tsv$",clusters[[j]]), "\t", ","), - # header = FALSE) - clusterMap[,1] <- paste0(getSampleData(quantData_i)$sampleName[j], - "_",clusterMap[,1]) - clusterMaps <- rbind(clusterMaps, clusterMap) - } - clustering <- splitAsList(clusterMaps[,1], clusterMaps[,2]) - rm(clusterMaps) - rm(clusterMap) - iter <- clustering - - } else{ #if clusters is a list - iter <- clusters[[i]] - } - # cluster mode: convert barcode strings to integer column indices; - # iter becomes a named list of integer vectors, one per cluster - iter <- lapply(iter, match, getSampleData(quantData_i)$id) + # cluster mode: one column per cluster, pooling that cluster's cells/spots + sampleName <- names(quantData)[i] + # find each clustered cell/spot's column in this sample; those from other + # samples don't match here and drop out as NA + idx <- match(clusterDf$id, getSampleData(quantData_i)$id) + inSample <- !is.na(idx) # cells/spots belonging to sample i + # columnGroups is a named list, one entry per cluster, each holding + # that cluster's columns, e.g. + # cluster0: 1 4 5 6 ... + # cluster1: 13 19 20 ... + columnGroups <- split(idx[inSample], clusterDf$cluster[inSample]) + # each cluster's barcodes, in the same index order as in columnGroups + clusterBarcodes <- split(getSampleData(quantData_i)$barcode[idx[inSample]], + clusterDf$cluster[inSample]) + clusterLabels <- names(columnGroups) + # two samples can both have a cluster, so prefix the sample name back + # so the columns don't clash when samples are combined into a single SE + # later + names(columnGroups) <- paste0(sampleName, "_", clusterLabels) + } else { + # no clusters: columnGroups is an integer vector with one column per + # input column, i.e. one per cell/spot for single-cell/spatial or one + # per sample for bulk (1, 2, 3, ...) + columnGroups <- seq_len(nrow(getSampleData(quantData_i))) } - countsSeCompressed <- bplapply(iter, FUN = function(columnIdx){ # previous i changed to j to avoid duplicated assignment + countsSeCompressed <- bplapply(columnGroups, FUN = function(columnIdx){ # previous i changed to j to avoid duplicated assignment incompatibleCounts_i <- getIncompatibleCounts(quantData_i)[, columnIdx, drop = FALSE] @@ -328,13 +332,14 @@ bambu <- function(reads, annotations = NULL, genome = NULL, NDR = NULL, end.ptm <- proc.time() message("Total Time ", round((end.ptm - start.ptm)[3] / 60, 3), " mins.") if(!is.null(clusters)){ - ColNames <- c(ColNames, names(iter)) + ColNames <- c(ColNames, names(columnGroups)) colData.all[[i]] <- data.frame( - id = names(countsSeCompressed), - sampleName = sub("_[^_]+$", "", names(countsSeCompressed)), - cluster = sub(".*_", "", names(countsSeCompressed)), - row.names = names(countsSeCompressed) + id = names(columnGroups), + sampleName = sampleName, + cluster = clusterLabels, + row.names = names(columnGroups) ) + colData.all[[i]]$barcodes <- unname(clusterBarcodes) } else{ ColNames <- c(ColNames, rownames(getSampleData(quantData_i))) colData.all[[i]] <- data.frame(getSampleData(quantData_i)) @@ -364,56 +369,249 @@ bambu <- function(reads, annotations = NULL, genome = NULL, NDR = NULL, } } -#' Single-cell isoform reconstruction and quantification -#' @title Single-cell isoform reconstruction and quantification with Bambu -#' @description Analyse single-cell long-read RNA-seq data with Bambu, -#' performing isoform discovery and quantification at single-cell resolution. -#' This function calls the main \code{\link{bambu}} function with cell -#' barcode/UMI extraction and UMI-based deduplication enabled by default, -#' and returns a \emph{SummarizedExperiment} object with per-cell transcript -#' expression estimates. +#' @title Single-cell / Spatial transcript discovery and quantification with Bambu +#' @description A function wrapper for \code{\link{bambu}} that runs on long-read RNA-seq data to +#' perform transcript discovery and quantification at the single-cell/spatial level. +#' The function calls \code{\link{bambu}} with cell +#' barcode/UMI extraction and UMI-based deduplication enabled by default, with an optional +#' parameter for clustered-level transcript quantification. #' -#' We recommend processing single-cell data using the Bambu Nextflow pipeline, -#' which handles preprocessing, barcode demultiplexing, and alignment prior to -#' running this function. See \url{https://github.com/GoekeLab/bambu-singlecell-spatial}. -#' @param reads A string or vector of strings specifying paths to BAM files. -#' BAM files must contain cell barcode and UMI information, either as BAM tags -#' (\code{CB} for cell barcode, \code{UB} for UMI) or encoded in the read name -#' using the format \code{CB_UMI#READNAME} (note: \code{CB} and \code{UMI} -#' must not contain underscores). -#' @param annotations A path to a .gtf file or a \code{TxDb} object for -#' transcript annotations. Defaults to NULL. -#' @param genome A path to a fasta file or a \code{BSGenome} object. -#' Defaults to NULL. -#' @param NDR Numeric specifying the maximum NDR rate for novel transcript -#' discovery. Defaults to NULL. -#' @param discovery Logical, whether transcript discovery is performed. -#' Defaults to TRUE. -#' @param assignDist Logical, whether to assign reads to transcripts. -#' Defaults to TRUE. -#' @param quant Logical, whether quantification is performed. Defaults to TRUE. -#' @param clusters A named list mapping cluster names to barcode vectors, -#' used for cluster-level quantification. Defaults to NULL. -#' @param stranded Logical, whether reads are stranded. Defaults to FALSE. -#' @param ncore Integer specifying the number of cores for parallel processing. -#' Defaults to 1. +#' For general users, we recommend to use bambu-pipe, an end-to-end single-cell/spatial pipeline that +#' handles preprocessing, barcode demultiplexing, and alignment prior to +#' transcript discovery and quantification. See \url{https://github.com/GoekeLab/bambu-pipe}. +#' @inheritParams bambu +#' @param reads For \code{output} in \code{"readClasses"}, +#' \code{"extendedAnnotations"}, or \code{"quantData"}, the input reads as in +#' \code{\link{bambu}}: paths to BAM files, a \code{BamFileList}, or the read-class +#' objects returned by the \code{"readClasses"} stage. For \code{output} in +#' \code{"uniqueCounts"}, \code{"EM"}, or \code{"clusteredEM"}, the per-sample +#' \code{quantData} list returned by the \code{"quantData"} stage. +#' @param output Bambu performs transcript discovery and quantification on single-cell +#' data in distinct steps, that are each specified using the \code{output} option. Must +#' be one of \code{"readClasses"}, \code{"extendedAnnotations"}, \code{"quantData"}, +#' \code{"uniqueCounts"}, \code{"EM"}, or \code{"clusteredEM"}. See the \strong{Details} +#' section for what each stage does and the \strong{Value} section for what it returns. +#' @param clusters Assignment of cells or spatial spots to clusters, supplied as one of: +#' a named vector whose names are the cell/spot identifiers and whose values are the +#' cluster labels; a \code{data.frame} object with an \code{id} and a \code{cluster} +#' column; or a path to a \code{.csv}/\code{.tsv}/\code{.txt} file holding +#' that \code{data.frame}. In every case the identifier follows the +#' \code{sampleName_barcode} format and the cluster label is the cluster that cell/spot +#' belongs to. Must be used with \code{output = "clusteredEM"}, where expression is +#' aggregated to the cluster level instead of per individual cell/spot. #' @param ... Additional arguments passed to \code{\link{bambu}}, such as -#' \code{verbose}, \code{lowMemory}, \code{opt.discovery}, etc. +#' \code{sampleData}, \code{opt.em}, \code{trackReads}, \code{returnDistTable}, +#' \code{yieldSize}, and \code{verbose}. See \code{\link{bambu}} for the full set. +#' @details +#' Single-cell and spatial long read RNA-Seq protocols tag every read with a barcode +#' and a unique molecular identifier (UMI). The barcode identifies which cell +#' (single-cell data) or spatial spot (spatial data) a read came from; below, "cell" +#' and "single-cell" refer to either. \code{\link{bambu.singlecell}} reads the barcode +#' and UMI for each read and collapses duplicate molecules, so that transcript +#' quantification can be performed per barcode. +#' +#' \code{\link{bambu.singlecell}} runs one stage at a time, each specified with the +#' \code{output} argument for the returned object. The step-by-step workflow is +#' described below and in the examples: +#' \enumerate{ +#' \item \strong{Read Class Construction} (\code{output = "readClasses"}): build +#' per-sample read classes from the demultiplexed reads. Requires \code{reads} +#' (paths to BAM files), \code{genome} (reference genome sequence), and +#' \code{annotations} (optional, but recommended). +#' \item \strong{Transcript Discovery} (\code{output = "extendedAnnotations"}): +#' discover novel transcripts and extend the reference annotations. Requires +#' \code{reads} (the read classes from step 1) and \code{annotations} (optional, +#' but recommended; use the same annotations as in step 1, or omit them in both). +#' \item \strong{Read to Transcript Assignment} (\code{output = "quantData"}): +#' assign reads to the transcripts, giving the per-sample read-to-transcript +#' assignments. Requires \code{reads} (the read classes from step 1) and +#' \code{annotations} (the extended annotations from step 2 or the original +#' annotations, which cannot be NULL if transcript discovery was skipped). +#' \item \strong{Transcript Quantification} with the Expectation-Maximization +#' (EM) algorithm, or without it (EM-free): +#' \itemize{ +#' \item \strong{EM-free unique counts} (\code{output = "uniqueCounts"}): +#' transcript-level unique counts, ready for downstream analysis. Apply +#' \code{\link{transcriptToGeneExpression}} for gene counts, or cluster cells +#' by gene expression to run the clustered EM. Requires \code{reads} (the +#' quantData list from step 3) and the same \code{annotations} used in step 3. +#' \item \strong{Single-Cell EM} (\code{output = "EM"}): estimates +#' single-cell transcript expression with the EM. Requires \code{reads} (the +#' quantData list from step 3) and the same \code{annotations} used in step 3. +#' Single-cell EM is resource intensive; for higher efficiency use \code{uniqueCounts}, +#' or use \code{clusteredEM} to run the EM directly on cell clusters for more +#' accurate, stable estimates. +#' \item \strong{Clustered EM} (\code{output = "clusteredEM"}): runs the EM +#' directly on user-defined cell clusters, giving more stable estimates than +#' quantifying single cells and clustering afterwards. Requires \code{reads} +#' (the quantData list from step 3), the same \code{annotations} used in +#' step 3, and \code{clusters}. +#' } +#' } +#' This makes it practical to run different stages in different R sessions to minimise +#' memory usage, and also lets users reuse intermediate results and explore different +#' downstream settings without reprocessing large datasets from the raw alignments. See +#' the \strong{Value} section for the returned output from each stage and \strong{Examples} +#' for a detailed step-by-step run. +#' @return The object returned depends on the specified \code{output} argument: +#' \describe{ +#' \item{\code{output = "readClasses"}}{returns a list of read-class objects for +#' each sample that represent a summarised representation of similar reads with +#' cell-barcode-to-read mapping information. ReadClass objects can be used as input +#' for transcript discovery and read to transcript assignment.} +#' \item{\code{output = "extendedAnnotations"}}{returns the extended annotations +#' after transcript discovery as a \code{GRangesList} object. The extended +#' annotation object can be used as input for transcript to read assignment and +#' quantification.} +#' \item{\code{output = "quantData"}}{returns a list of \code{quantData} objects +#' that contain the read-to-transcript assignments for each sample. The quantData +#' object is used to calculate unique counts per transcript and gene-level counts, +#' and also can be used as input for EM quantification.} +#' \item{\code{output = "uniqueCounts"}}{returns a \code{SummarizedExperiment} of +#' transcript-level unique counts (reads compatible with a single transcript), with +#' one column per cell. Unique counts are ready for downstream analysis as is. +#' Collapse them to gene-level counts with \code{\link{transcriptToGeneExpression}}, +#' or use them to cluster cells by gene expression profile (see \code{output = "clusteredEM"}).} +#' \item{\code{output = "EM"}}{returns a \code{SummarizedExperiment} with transcript +#' counts estimated using Bambu's EM for each barcode/single cell. This step is +#' resource intensive, consider using unique counts or the clustered EM.} +#' \item{\code{output = "clusteredEM"}}{returns a \code{SummarizedExperiment} with +#' one column per cluster as specified in the cluster argument. \code{colData} +#' describes the \code{cluster} label and the set of cell barcodes belonging to the +#' cluster. Quantifying transcript expression for cell clusters provides more stable +#' EM estimates.} +#' } +#' @seealso \code{\link{bambu}} for the underlying function and the full +#' parameter set; \code{\link{transcriptToGeneExpression}} for preparing gene-level counts +#' for cell clustering and clustered quantification, or for direct use in downstream analysis. #' @examples -#' sc.bam <- system.file("extdata", "demultiplexed.bam", package = "bambu") -#' fa.file <- system.file("extdata", +#' ## single-cell example data: two demultiplexed samples (PacBio and ONT) +#' sce.dir <- system.file("extdata", "single_cell", package = "bambu") +#' reads <- file.path(sce.dir, c( +#' "GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000.bam", +#' "GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000.bam")) +#' annotations <- system.file("extdata", +#' "Homo_sapiens.GRCh38.91_chr9_1_1000000.gtf", package = "bambu") +#' genome <- system.file("extdata", #' "Homo_sapiens.GRCh38.dna_sm.primary_assembly_chr9_1_1000000.fa", #' package = "bambu") -#' gr <- readRDS(system.file("extdata", -#' "annotationGranges_txdbGrch38_91_chr9_1_1000000.rds", -#' package = "bambu")) -#' se <- bambu.singlecell(reads = sc.bam, annotations = gr, genome = fa.file) +#' +#' ## Each stage is selected with an 'output' preset; its output feeds the next. +#' ## 1. read-class construction from the demultiplexed BAMs +#' readClasses <- bambu.singlecell(reads = reads, output = "readClasses", +#' annotations = annotations, genome = genome) +#' +#' ## 2. extend the annotations with novel transcripts discovered by Bambu +#' extendedAnnotations <- bambu.singlecell(reads = readClasses, +#' output = "extendedAnnotations", annotations = annotations) +#' +#' ## 3. assign reads to transcripts (per-sample quantData list) +#' quantData <- bambu.singlecell(reads = readClasses, +#' output = "quantData", annotations = extendedAnnotations) +#' +#' ## 4. EM-free quantification. Steps 4a-4b give unique counts and gene counts ready +#' ## for downstream analysis; stop here, or continue to step 5 for EM-based +#' ## quantification. +#' +#' ## 4a. transcript-level unique counts (reads compatible with a single transcript). +#' ## These EM-free counts are ready to use downstream; collapse them to gene counts, +#' ## or cluster cells by gene expression to run the clusteredEM step. +#' uniqueCountsSe <- bambu.singlecell(reads = quantData, +#' output = "uniqueCounts", annotations = extendedAnnotations) +#' +#' ## 4b. gene-level counts: transcriptToGeneExpression() sums the transcript counts +#' ## within each gene to give a gene-level SummarizedExperiment. +#' geneCountsSe <- transcriptToGeneExpression(uniqueCountsSe) +#' +#' ## 5. (optional) EM-based quantification. Steps 5a-5b estimate transcript expression +#' ## with the EM, either at single-cell resolution (5a) or per user-defined cluster +#' ## for more stable estimates (5b). +#' +#' ## 5a. single-cell transcript counts estimated with the EM. +#' EMCountsSe <- bambu.singlecell(reads = quantData, +#' output = "EM", annotations = extendedAnnotations) +#' +#' ## 5b. per-cluster (pseudobulk) transcript counts estimated with the EM. +#' ## How the clusters are defined is up to the analysis. One common approach is to +#' ## cluster by cell-type using the gene-level counts from step 4 with Seurat (not a +#' ## bambu dependency, so this block runs only when Seurat is available): +#' if (requireNamespace("Seurat", quietly = TRUE)) { +#' seurat <- Seurat::CreateSeuratObject( +#' counts = assays(geneCountsSe)$counts, min.cells = 1) +#' seurat <- Seurat::NormalizeData(seurat) +#' ## this example is small, so cluster on all genes with exact PCA and a small +#' ## number of components (npcs must stay below the gene count) +#' seurat <- Seurat::ScaleData(seurat, features = rownames(seurat)) +#' seurat <- Seurat::RunPCA(seurat, features = rownames(seurat), +#' npcs = 10, approx = FALSE) +#' seurat <- Seurat::FindNeighbors(seurat, dims = 1:10) +#' seurat <- Seurat::FindClusters(seurat, resolution = 0.05) +#' +#' ## Seurat's active.ident is a named vector (names = ids in sampleName_barcode +#' ## format, values = cluster labels), one of the accepted input data type for +#' ## `clusters` argument in `bambu.singlecell()` +#' clusters <- seurat@active.ident +#' } +#' +#' ## In general, clusteredEM accepts the clusters argument as a named vector (e.g. +#' ## seurat@active.ident), a data.frame, or a path to a .csv/.tsv/.txt file. Here we +#' ## illustrate the input data type with a stored .csv: +#' clusters <- file.path(sce.dir, +#' "clusters_GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCell_multisample_chr9_1_1000000.csv") +#' clusteredEMCountsSe <- bambu.singlecell(reads = quantData, +#' output = "clusteredEM", annotations = extendedAnnotations, +#' clusters = clusters) #' @export -bambu.singlecell <- function(reads, annotations = NULL, genome = NULL, NDR = NULL, - clusters = NULL, discovery = TRUE, assignDist = TRUE, quant = TRUE, - stranded = FALSE, ncore = 1, ...) { +bambu.singlecell <- function(reads, output, annotations = NULL, genome = NULL, + clusters = NULL, NDR = NULL, stranded = FALSE, ncore = 1, ...) { + # 'output' is required: it names the single stage to run and return + if (missing(output) || is.null(output)) + stop("'output' must be one of 'readClasses', 'extendedAnnotations', ", + "'quantData', 'uniqueCounts', 'EM', or 'clusteredEM'.") + output <- match.arg(output, c("readClasses", "extendedAnnotations", + "quantData", "uniqueCounts", "EM", "clusteredEM")) + # For the EM/clusteredEM stages the quantData list is passed in through 'reads', + # but bambu() currently takes quantData as a separate argument, so it is rerouted + # (reads -> quantData) below and defaults to NULL otherwise. + # TODO: once bambu() accepts a quantData list directly via 'reads', this + # initializer and the reroute can be removed and 'reads' passed straight through. + quantData <- NULL + # validate the inputs each preset needs before setting the stage flags + # (a genome for BAM input is checked downstream in checkInputs()) + if (output %in% c("readClasses", "extendedAnnotations", "quantData") && + is.null(reads)) + stop("output = '", output, "' requires 'reads' to be provided.") + if (output %in% c("uniqueCounts", "EM", "clusteredEM") && is.null(reads)) + stop("output = '", output, "' requires the per-sample 'quantData' list ", + "(from output = 'quantData') to be supplied as 'reads'.") + # annotations is optional for extendedAnnotations (de novo discovery), but + # required for the assignment and quantification stages + if (output %in% c("quantData", "uniqueCounts", "EM", + "clusteredEM") && (is.null(annotations) || length(annotations) == 0)) + stop("output = '", output, "' requires 'annotations' to be provided.") + if (output %in% c("uniqueCounts", "EM") && !is.null(clusters)) + stop("output = '", output, "' ignores 'clusters'; use output = 'clusteredEM'.") + if (output == "clusteredEM" && is.null(clusters)) + stop("output = 'clusteredEM' requires 'clusters' to be provided.") + # each (named) preset returns exactly one object + if (output == "readClasses") { + discovery <- FALSE; assignDist <- FALSE; quant <- FALSE + } else if (output == "extendedAnnotations") { + discovery <- TRUE; assignDist <- FALSE; quant <- FALSE + } else if (output == "quantData") { + discovery <- FALSE; assignDist <- TRUE; quant <- FALSE + } else if (output == "uniqueCounts") { + return(getUniqueCountsSe(reads, annotations)) + } else { + discovery <- FALSE; assignDist <- FALSE; quant <- TRUE + # the quantification stages take the quantData list through 'reads' + quantData <- reads + reads <- NULL + } bambu(reads = reads, annotations = annotations, genome = genome, NDR = NDR, - opt.singlecell = list(extractBarcodeUMI = TRUE, dedupUMI = TRUE, clusters = clusters), + quantData = quantData, + opt.singlecell = list(extractBarcodeUMI = TRUE, dedupUMI = TRUE, + clusters = clusters), discovery = discovery, assignDist = assignDist, quant = quant, stranded = stranded, ncore = ncore, ...) } diff --git a/R/bambu_utilityFunctions.R b/R/bambu_utilityFunctions.R index c5056d20..75284eb5 100644 --- a/R/bambu_utilityFunctions.R +++ b/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) } diff --git a/R/transcriptToGeneExpression.R b/R/transcriptToGeneExpression.R index e39b9c8a..8e63e633 100644 --- a/R/transcriptToGeneExpression.R +++ b/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)) diff --git a/README.md b/README.md index c9fc2b67..7f4a3813 100755 --- a/README.md +++ b/README.md @@ -11,12 +11,13 @@ [![BioC status](http://bioconductor.org/shields/build/release/bioc/bambu.svg)](http://bioconductor.org/checkResults/release/bioc-LATEST/bambu/) [![BioC dev status](http://www.bioconductor.org/shields/build/devel/bioc/bambu.svg)](https://bioconductor.org/checkResults/devel/bioc-LATEST/bambu) [![License: GPL v3](https://img.shields.io/badge/License-GPLv3-blue.svg)](https://www.gnu.org/licenses/gpl-3.0) -[![CodeFactor](https://www.codefactor.io/repository/github/goekelab/bambu/badge)](https://www.codefactor.io/repository/github/goekelab/bambu) [![codecov](https://codecov.io/gh/GoekeLab/bambu/branch/master/graph/badge.svg?token=PMeRi0r1tj)](https://codecov.io/gh/GoekeLab/bambu) ***bambu*** is a R package for multi-sample transcript discovery and quantification using long read RNA-Seq data. You can use ***bambu*** after read alignment to obtain expression estimates for known and novel transcripts and genes. The output from ***bambu*** can directly be used for visualization and downstream analysis such as differential gene expression or transcript usage. +To use Bambu for single-cell and spatial long-read RNA-Seq data, we provide [bambu-pipe](https://github.com/GoekeLab/bambu-pipe), an end-to-end pipeline that handles raw data preprocessing, demultiplexing, read alignment, transcript discovery, and quantification. To use Bambu directly for single cell and spatial data, please refer to the documentation of `bambu.singlecell()`. + ### Content - [Installation](#installation) @@ -26,8 +27,7 @@ - [Modulating the sensitivity of discovery (pre and post analysis)](#Modulating-the-sensitivity-of-discovery-pre-and-post-analysis) - [Output](#Output) - [Visualization](#Visualization) - - [Single-Cell and Spatial](#Single-Cell-and-Spatial) - - [*bambu* Advanced Options](#Bambu-Advanced-Options) +- [*bambu* Advanced Options](#Bambu-Advanced-Options) - [Using a pretrained model](#Using-a-pretrained-model) - [De-novo transcript discovery](#De-novo-transcript-discovery) - [Storing and using preprocessed files (rcFiles)](#Storing-and-using-preprocessed-files-rcFiles) @@ -36,7 +36,6 @@ - [Quantification of gene expression](#Quantification-of-gene-expression) - [Including single exons](#Including-single-exons) - [Fusion gene/isoform detection](#Fusion-geneisoform-detection) - - [Custom single-cell and spatial analysis](#Custom-Single-Cell-and-Spatial) - [*bambu* Arguments](#Bambu-Arguments) - [Output Description](#Output-Description) - [Release History](#Release-History) @@ -291,12 +290,6 @@ plotBambu(se, type = "heatmap", group.var) # heatmap plotBambu(se, type = "pca", group.var) # PCA visualization ``` -### Single-Cell-and-Spatial - -There is a single-cell and spatial pipeline starting from fastq or demultiplexed bam files that include demultiplexing and aligning available here https://github.com/GoekeLab/bambu-singlecell-spatial. We recommend using this pipeline where possible. - -For advanced users see the #[Custom single-cell and spatial analysis](#Custom-Single-Cell-and-Spatial) section under advanced options - ### *Bambu* Advanced Options Below we include several advanced options and use-cases for *bambu*. We recommend reading and understanding the [paper](https://doi.org/10.1038/s41592-023-01908-w) before attempting to use these features. @@ -478,106 +471,6 @@ To use this feature, it is recommended to detect the fusion gene breakpoints usi se <- bambu(reads = fusionAligned.bam, annotations = fusionAnnotations, genome = fusionFasta, fusionMode = TRUE) ``` -### Custom single-cell and spatial - -If you want to run Bambu-Clump for single-cell or spatial analysis stand alone and not part of the Bambu-Pipe pipeline we recommend running it in 4 stages which we will describe seperately: Read Class Construction, Transcript Discovery, Read Class Assignment, and EM Quantification. Note that this section will only cover arguments that are different or unique to this analysis. - -#### Read Class Construction: - -**reads**: provided bam files should have barcodes in the read name or in the BC tag ( and UG tag for UMI identifiers). In the case where both tags and read names contain barcode information, tags will be used a prior. If not, a regular delimited headerless file that contain the demultiplexing information for each read should be provided to demultiplexed argument below. For exact requirements see https://github.com/GoekeLab/bambu-singlecell-spatial.
- -**demultiplexed**: should be either set to TRUE or the path to barcode mapping file. Otherwise, bambu will not look for barcodes and seperate reads by barcode rather than sample.
- -Optional: - -**cleanReads**: A logical TRUE/FALSE. Chimeric reads in samples can cause issues with barcode assignments. Setting this to TRUE will ensure only the first alignment per barcode is used (We recommend using this).
- -**sampleNames**: A vector of characters assigning names to each sample in the reads argument. By default the sample names are taken from the file names and appended to the barcodes in order to differentiate them. If your sample names are the same across multiple files, but matching barcodes between the samples should be counted seperately, provide them with different sample names using this argument. Similiarly if your samples have different names, but overlapping barcodes should be counted together, give them the same sample name with this argument.
- -**dedupUMI**: A logical TRUE/FALSE.
- -**barcodesToFilter**: A string vector indicating barcodes to be filtered out.
- -```rscript -readClassFile <- bambu(reads = samples, annotations = annotations, genome = fa.file, ncore = 1, discovery = FALSE, quant = FALSE, demultiplexed = barcode_maps, verbose = TRUE, assignDist = FALSE, lowMemory = as.logical("$params.lowMemory"), yieldSize = 10000000, sampleNames = ids, cleanReads = as.logical($cleanReads), dedupUMI = as.logical($deduplicateUMIs)) -``` - -#### Transcript Discovery: - -Transript discovery can be run as usual as typically bulk-level discovery is suitable. - -```rscript -extendedAnno <- bambu(reads = readClassFile, annotations = annotations, genome = fa.file, ncore = 1, discovery = TRUE, quant = FALSE, demultiplexed = TRUE, verbose = FALSE, assignDist = FALSE) -``` - -#### Read Class Assignment: - -This step was previously performed together with the quantification, but can be done seperately so that the arguments can be passed to the quantification seperately with different clustering. If you only want barcode level gene counts or unique transcript counts you can stop here and do not need to proceed to the EM quantification. - -**spatial**: This should be a path to your barcode whitelist that also contains the x and y coordinates as extra columns. If provided, the file should contain 3 columns with or without header, where the first column is the barcode, and the second and third column contains the x and y coordinates information accordingly. Compressed file format is accepted as well. - -```rscript -quantData <- bambu(reads = readClassFile, annotations = extendedAnno, genome = fa.file, ncore = 1, discovery = FALSE, quant = FALSE, demultiplexed = TRUE, verbose = FALSE, opt.em = list(degradationBias = FALSE), assignDist = TRUE, spatial = spatial) -``` - -#### EM quantification: - -If you plan to run this step with multiple processes we recommend restarting your R instance to ensure that environmental variables do not inflate the memory usage. - -**reads**: This argument is still mandatory but not needed when performing quantification alone as long as you provide the quantData argument
- -**quantData**: This is the summerized experiement output from the Read Class Assignment step
- -**clusters**: This is an optional argument which is either a path to a csv containing the barcode to cluster assignments or a CharacterList which can be produced using the code below.
- -**opt.em = list(degradationBias=FALSE)**: We recommend including this argument if you are doing barcode level EM quantification to greatly improve runtime with only a small reduction in quantification accuracy. - -```rscript -#use Seurat to generate clusters from gene counts -library(Seurat) - -clusterCells <- function(counts, resolution = 0.8, dim = 15){ - - cellMix <- CreateSeuratObject(counts = counts, - project = "cellMix", min.cells = 1)#, min.features = 200) - #cellMix <- subset(cellMix, subset = nFeature_RNA > nFeature_RNA_threshold & nFeature_RNA < nFeature_RNA_threshold_max) - #nFeature_RNA_threshold <- 1000, nFeature_RNA_threshold_max = 9000, - cellMix <- NormalizeData(cellMix, normalization.method = "LogNormalize", scale.factor = 10000) - cellMix <- FindVariableFeatures(cellMix, selection.method = "vst", nfeatures = 2500) - all.genes <- rownames(cellMix) - cellMix <- ScaleData(cellMix, features = all.genes) - npcs <- ifelse(ncol(counts)>50, 50, ncol(counts)-1) - cellMix <- RunPCA(cellMix, features = VariableFeatures(object = cellMix), npcs = npcs) - dim <- ifelse(dim >= dim(cellMix@reductions$pca)[2], dim(cellMix@reductions$pca)[2],dim) # if data dimension is small, otherwise, cap dimension at 15 - cellMix <- FindNeighbors(cellMix, dims = 1:dim) - cellMix <- FindClusters(cellMix, resolution = resolution) - cellMix <- RunUMAP(cellMix, dims = 1:dim) - - return(cellMix) -} - -quantData.gene <- transcriptToGeneExpression(quantData) -counts <- assays(quantData.gene)$counts #selecting first sample -cellMix <- clusterCells(counts) #resolution can be customized. For larger clusters: 0.2-0.6, for higher resolution: 0.8-2 -x <- setNames(names(cellMix@active.ident), cellMix@active.ident) -clusters_temp <- splitAsList(unname(x), paste)("cluster",names(x)))#make clusters names start with cluster, for better comprehension - - - -se <- bambu( reads = NULL, - annotations = rowRanges(quantDatas), - genome = "$genome", - quantData = quantDatas, - assignDist = FALSE, - ncore = $params.ncore, - discovery = FALSE, - quant = TRUE, - demultiplexed = TRUE, - verbose = FALSE, - opt.em = list(degradationBias = FALSE), - clusters = clusters_temp) -``` - ### *Bambu* Arguments |argument|description| @@ -592,21 +485,16 @@ se <- bambu( reads = NULL, | yieldSize | see Rsamtools. | | opt.discovery | A list of controlling parameters for isoform reconstruction process:
**prefix** specifying prefix for new gene Ids (genePrefix.number), defaults to empty
**remove.subsetTx** indicating whether filter to remove read classes which are a subset of known transcripts, defaults to TRUE
**min.readCount** specifying minimun read count to consider a read class valid in a sample, defaults to 2
**min.readFractionByGene** specifying minimum relative read count per gene, highly expressed genes will have many high read count low relative abundance transcripts that can be filtered, defaults to 0.05
**min.sampleNumber** specifying minimum sample number with minimum read count, gene read proportion, and TPS, defaults to 1
**min.exonDistance** specifying minimum distance to known transcript to be considered valid as new, defaults to 35bp
**min.exonOverlap** specifying minimum number of bases shared with annotation to be assigned to the same gene id, defaults to 10bp
**min.primarySecondaryDist** specifying the minimum number of distance threshold between a read class and the annotations internal exons. Read classes with distances less than the threshold are not annotated as novel and counted with the annotations for quantification, defaults to 5bp
**min.primarySecondaryDistStartEnd1** specifying the minimum number of distance threshold between a read class and the annotations start/end exons. Read classes with distances less than the threshold are not annotated as novel, defaults to 5bp
**min.primarySecondaryDistStartEnd2** specifying the minimum number of distance threshold between a read class and the annotations start/end exons. Read classes with distances less than the threshold are counted with the annotations, defaults to 5bp
**min.txScore.multiExon** specifying the minimum transcript probility score threshold for multi-exon transcripts for min.sampleNumber, defaults to 0
**min.txScore.singleExon** specifying the minimum transcript probability score threshold for single-exon transcripts for min.sampleNumber
**fitReadClassModel** a boolean specifying if bambu should train a model on each sample. If set to false bambu will use the default model for ranking novel transcripts. defaults to TRUE
**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
**returnModel** a boolean specifying if bambu will output the model it trained on the data, defaults to FALSE
**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.
**min.readFractionByEqClass** indicating the minimum relative read count of a subset transcript compared to all superset transcripts (ie the relative read count within the minimum equivalent class). This filter is applied on the set of annotations across all samples using the total read count, this is not a per-sample filter. Please use with caution. defaults to 0 | | opt.em | A list of controlling parameters for quantification algorithm estimation process:
**maxiter** specifying maximum number of run iterations, defaults to 10000
**degradationBias** correcting for degradation bias, defaults to TRUE
**conv** specifying the covergence threshold control, defaults to 0.0001
**minvalue** specifying the minvalue for convergence consideration, defaults to 0.00000001 | +| opt.singlecell | A list of single-cell specific parameters:
**extractBarcodeUMI** whether to extract cell barcodes and UMIs from BAM tags or read names, only applicable on single-cell and spatial samples, defaults to FALSE
**dedupUMI** whether to perform UMI-based deduplication per barcode, defaults to FALSE
**clusters** a path to a .csv/.tsv/.txt file assigning cells to clusters, with two columns: `id`, the cell identifier following the `sampleName_barcode` format, and `cluster`, the cluster label for that cell; .csv is comma-separated, .tsv/.txt tab-separated. When provided, transcript quantification is performed per cluster instead of per cell | | trackReads | When TRUE read names will be tracked and output as metadata in the final output as readToTranscriptMaps detailing the assignment of reads to transcripts.The output is a list with an entry for each sample. | | returnDistTable | When TRUE the calculated distance table between read classes and annotations will be output as metadata as distTables. The output is a list with an entry for each sample. | | discovery | A logical variable indicating whether annotations are to be extended for quantification, defaults to TRUE. | | quant | A logical variable indicating whether quantification will be performed, defaults to TRUE. | | verbose | A logical variable indicating whether processing messages will be printed. | -| mode | A string that will set other input arguments ['bulk', 'multiplexed', 'fusion', 'debug']
bulk -
    processByChromsome = FALSE
multiplexed -
    demultiplex = TRUE
    cleanReads = TRUE
    opt.em = list(degradationBias = FALSE)
    quant = FALSE
    processByChromosome = TRUE
fusion -
    NDR = 1
    fusionMode = TRUE
debug -
    verbose = TRUE
    trackReads = TRUE
    returnDistTable = TRUE | -| demultiplexed | A logical variable indicating whether the input bam file is demultiplexed. The barcode and umi either need to be present in the read name or the $BC and $UG tags, defaults to FALSE. Alternatively a path to a csv file can be provided where column 1 is read names, column 2 is barcodes, and column 3 is UMI. | -| spatial | A path to the barcode whitelist containing X and Y coordinates, defaults to null. If provided, the file should contain 3 columns with or without header, where the first column is the barcode, and the second and third column contains the x and y coordinates information accordingly. Compressed file format is accepted as well.| +| mode | A string that will set other input arguments ['bulk', 'multiplexed', 'fusion', 'debug']
bulk -
    processByBam = TRUE
    processByChromosome = FALSE
multiplexed -
    opt.singlecell = list(extractBarcodeUMI = TRUE)
    opt.em = list(degradationBias = FALSE)
    quant = FALSE
    processByChromosome = TRUE
fusion -
    NDR = 1
    fusionMode = TRUE
debug -
    verbose = TRUE
    trackReads = TRUE
    returnDistTable = TRUE | +| sampleData | A character vector of paths to metadata CSV files (or NA if unavailable for specific samples), defaults to NULL. Files must contain a "sampleName" column for bulk data or a "barcode" column for single-cell/spatial data. For bulk data, one metadata CSV file for all samples is sufficient, whereas single-cell/spatial data requires one metadata CSV file per sample. | | assignDist | A logical variable indicating whether read class to transcript assignment will be performed, defaults to TRUE. | | quantData | Advanced use only. A list of se outputs from the assignDist step. Used only to run quantification | -| sampleNames | A vector of strings representing the sample name associated with each input bam. bam files with the same sample name will be combined | -| cleanReads | A logical variable indicating whether only the first sequenced alignment in a read should be kept. This helps to remove chimeric reads, but will remove alignments from fusion genes, defaults to FALSE. | -| dedupUMI | A logical variable indicating whether UMI deduplication is performed. The longest read per UMI will be used and the rest discarded, defaults to FALSE.| -| barcodesToFilter | A vector of strings indicating the barcodes to be filtered out in reads.| -| clusters | Either a list containing the barcodes for each cluster, or a path to a csv file containg the barcode to cluster mapping. When provided, clusters will be used during discovery and EM quant steps, defaults to null. | | processByChromosome | A logical variable indicating if read classes will be constructed with all reads together (FALSE), or done by chromsome which uses less memory, but provides less information for the junction correction model (TRUE), defaults to FALSE | ### setNDR() arguments diff --git a/inst/extdata/demultiplexed.bam b/inst/extdata/demultiplexed.bam deleted file mode 100644 index ae85fd59..00000000 Binary files a/inst/extdata/demultiplexed.bam and /dev/null differ diff --git a/inst/extdata/single_cell/GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000.bam b/inst/extdata/single_cell/GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000.bam new file mode 100644 index 00000000..ee6f7f68 Binary files /dev/null and b/inst/extdata/single_cell/GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000.bam differ diff --git a/inst/extdata/single_cell/GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000.bam.bai b/inst/extdata/single_cell/GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000.bam.bai new file mode 100644 index 00000000..b2661a33 Binary files /dev/null and b/inst/extdata/single_cell/GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000.bam.bai differ diff --git a/inst/extdata/single_cell/GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000.bam b/inst/extdata/single_cell/GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000.bam new file mode 100644 index 00000000..d8b3cb4b Binary files /dev/null and b/inst/extdata/single_cell/GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000.bam differ diff --git a/inst/extdata/single_cell/GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000.bam.bai b/inst/extdata/single_cell/GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000.bam.bai new file mode 100644 index 00000000..93b79e4e Binary files /dev/null and b/inst/extdata/single_cell/GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000.bam.bai differ diff --git a/inst/extdata/single_cell/clusters_GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCell_multisample_chr9_1_1000000.csv b/inst/extdata/single_cell/clusters_GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCell_multisample_chr9_1_1000000.csv new file mode 100644 index 00000000..73d24ef4 --- /dev/null +++ b/inst/extdata/single_cell/clusters_GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCell_multisample_chr9_1_1000000.csv @@ -0,0 +1,248 @@ +id,cluster +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_AAACCGCTCACTTCCG,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_AACCCCTGTATCGCTC,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_AACGGGTCAGTGAAGC,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_AAGCATTTCATCCGCC,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_AAGGCCCGTCGCTTGT,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_AAGGGCCAGGCCCCAA,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_AAGTAGCAGGTCGCAT,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_AATCTAGGTCCAACAG,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_AATTACCAGAACAGAC,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_ACAATTGGTGCTCCTG,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_ACACGGAAGCTGTGGT,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_ACAGAGCCAGTCTCGA,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_ACCCAACAGGGCTGAA,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_ACGTGGGCAGTATGGA,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_ACTCCCGGTATCCGCT,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_AGAGTAGGTAGGCATT,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_AGCCGCATCGCCACTC,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_AGCGGTCGTAAGGTAG,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_AGCGTGTTCACAACAT,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_AGGGGTCAGAGCAACG,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_AGGTTCTCAGCAGGAT,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_AGTTGCGGTCGGTATG,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_AGTTTACAGCTGTATG,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_ATAGCGGAGGGCCATT,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_ATCGAGACAGTTAGCT,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_ATGCAATTCCTATCCG,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_ATGCAGCTCACAACAT,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_ATGCTCGCACTTAAAT,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CAAATCGTCGCCATAA,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CAATCACCAAAACTAA,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CAATTCTGTAATCGCG,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CCAACTCGTACGTTGG,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CCACCGTTCGTAGTCA,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CCACGGGGTAACGAGA,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CCAGTAGGTCGAGGGA,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CCGTTGGAGGGAATGT,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CCTAGACAGCAGCTGT,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CCTTATCAGCTAGGTC,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CGAGTAACACTTAAAT,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CGCAGTAGTAACATGC,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CGCATCACACCCTTAC,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CGTACAGGTCCGGCCT,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CGTACATCACGATGTT,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CGTGTCTGTACATATG,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CGTTTAGAGCCGATAC,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CTAAACCCAATGCTAC,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CTACAATAGACGCGAA,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CTCGATAAGGTCCCGA,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CTCTCCGCACATAACC,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CTGACGTAGCCTTAAG,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CTGGTGAGTGGAAGTA,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CTGTAACCAGTACGCC,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CTGTATTGTATGTACA,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_GCACGTTAGGGTTCCG,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_GCTTACTTCGCTACTT,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_GGAAATGGTCACGCTT,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_GGACACAAGATTGTGT,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_GGCATTGGTGCCTCGC,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_GGGAGATGTGGAAGTA,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_GGGATCAAGGCATACC,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_GGGTTGAAGCAAGAGT,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_GGTACTATCCGTTACT,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_GTAATTCAGACCTATG,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_GTCTCGCTCCGACTCC,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_AATTACCAGCGATGAG,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_ACATCTTGTAGGCCGT,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_ACATGCTCAATCGGGG,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_ACGACACTCATGGCCG,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_ACGAGTTGTCAGCCTT,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_ACTCAGGGTGCTTCGT,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_AGACTCGGTATACTTC,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_AGGTTTAGTAACTGCG,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_AGTAACGCACTGACGC,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_ATACTTCGTCACTAGG,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_ATGGGGCAGGCCATAC,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_ATGTGTCCAAATTCGC,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CAAGGGCAGCGCCCAA,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CAATAGGCAACCTGCA,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CAATTAATCCACCCGT,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CAGTATTAGCAAGACG,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CCACGTTGTGCTCTGT,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CCATAGCAGCGCATTA,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CGATGAGAGCTACAAT,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CGGTCATAGCACCTAA,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CGTCAAATCACTCGCC,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CTAACCCGTAAGCATC,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CTAGTCAAGGCATTCA,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CTCGTTTAGAAGGGCC,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CTGGATGGTAGGGGCA,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_GCAATGTGTACTCGAG,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_GGACGTTTCGTTATGG,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_GGAGTGGGTAATGTCA,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_GGCTGAAAGCTACATG,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_GGGTTCAGTCATTTCG,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_GTCGGTCCACCCTTCG,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_AAACCTCCAAGACAGC,cluster2 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_AAACGCACAGTTCGCA,cluster2 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_AAGCCGGCACGCTCAA,cluster2 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_AATTGGTTCCGGCACA,cluster2 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_ACAGGCTTCAGCTAAT,cluster2 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_ACTCTATGTCCGTAAG,cluster2 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_ACTGGGGCAAGACAGC,cluster2 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_AGACATAAGCTCGAAC,cluster2 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_AGGAAACTCCGACTAA,cluster2 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_AGGTGATAGGCTTCAA,cluster2 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_ATTCACGAGCGTTTGA,cluster2 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CAATCGGGTACAAGCT,cluster2 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CACCCTGGTGGAAATG,cluster2 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CACCTAAAGCGATGGA,cluster2 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CGTACATCACTCACCT,cluster2 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CTGGTAATCCTGTTCG,cluster2 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CTTAGTGAGGTCCGAT,cluster2 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_GCATAATTCCCGTGTC,cluster2 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_GCATACATCACTCATC,cluster2 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_GGGATCAAGATCTGCA,cluster2 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_AGTCTGGGTGGTATGT,cluster3 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_ATGGACACAGCAAGAC,cluster3 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CAAGTCTCACAAGACA,cluster3 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CAGTAGGAGGTGATTT,cluster3 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CCTGCAGAGAGGGCGT,cluster3 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CTGACTCGTGGTATTA,cluster3 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CTTATAGCAACACGAA,cluster3 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CTTGGCGAGGCCACAT,cluster3 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_CTTGTGCGTCTTAATC,cluster3 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_GCGACAGCACATGCGC,cluster3 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_GGGGTAATCCCACGCT,cluster3 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000_GGGGTTAAGGGGCGAT,cluster3 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_AAACCGCTCAGCCCTC,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_AAACTCCCACCATAAC,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_AAATCTGGTATGTTAG,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_AACATAGAGCAGTCTA,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_AAGCCTATCCCGAAAC,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_AATGGGAGTAGTTTAG,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_AATGGGTAGAAAGCCG,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_AATGTGCGTGCTTCGT,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_AATTGGTTCCGGACAT,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_ACAGTTTAGCTACAGC,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_ACCCGATAGGGTACGA,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_ACGAGTTGTCAGCCTT,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_ACGTCAATCCCTAATC,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_ACGTCAGGTATGTACA,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_ACGTGGGCAGTATGGA,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_ACTCGATCACACTCAG,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_AGAAGCTAGACGAGGC,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_AGCAGATAGGTTCATG,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_AGCGTACTCAATGTGT,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_AGCTGCCCAGTCAACT,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_AGGACCAAGACGCCCA,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_AGGCTCAGTGTCATGG,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_AGGCTCCCAAAGCAAT,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_AGTCCCGAGATGGTTT,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_AGTCGGCGTAGGCCCG,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_ATCAGCTCAGTAGGCG,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_ATCGAAGCACACTGCA,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_ATGATGCGTACGCATG,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_ATGCTACCAGTTAACC,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_ATTCTAGCAAGGGGAA,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_CAAAGCTGTATATGTG,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_CAATATATCAGGGACA,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_CACAATCCACGCAGGA,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_CACGAGTTCGGCCATG,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_CAGTCAGGTAAACTGC,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_CAGTCGCCACGACACC,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_CAGTGCCCAAAAGCCA,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_CATACAAGTCCCTGTG,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_CATCTTCAGGGAGGTC,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_CCAACATCACGTTAGT,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_CCACGTTGTGCTCTGT,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_CCATAGCAGCGCATTA,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_CCATATCTCCGACCCT,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_CCGACTACAAAAGTTA,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_CCGATACAGCCCTAGG,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_CCGCAATAGCTGATAT,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_CCTTAGCGTATTCGTT,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_CGACGTTAGACCAGAG,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_CGATGAGAGCTACAAT,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_CGCATCACACCCTTAC,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_CGCGTATAGAATGTAG,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_CTACAGGGTATGAGCT,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_CTCACGCTCCCAAAGC,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_GCACAGCGTACGTTGG,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_GCATAGACAGGCACTC,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_GCCTTCGGTACTAATG,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_GCGAAGCAGAGTTAGA,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_GCGGATAAGACTCGAC,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_GCGTAGCGTCTAAGGT,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_GCTGGCATCCCTATTA,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_GGAGTGGGTAATGTCA,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_GGTACTATCCGTTACT,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_GGTCAGCGTCAGTTGG,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_GTAATGAAGCGCTCGA,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_GTCGTAGGTAGCATCG,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_GTGTTCGCACAGGTCC,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_GTTCACTTCGCTAATG,cluster0 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_ACAGGACCAACATGGA,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_AGCGGACAGATGAGTA,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_AGCGTGTTCACAACAT,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_AGGCACAAGGCTCAAC,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_AGGGCTTTCGTTCGGG,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_AGTAATCTCAACCCAA,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_ATAGTGCGTAATGTTG,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_ATCTTGGAGCCACTAA,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_CAGTATTAGCAAGACG,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_CCACTCTCACGCAACC,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_CTTAACCTCCTGGTAA,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_GCATTCGCACCAGGCT,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_GCCAACCCAAAACGGC,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_GCTTACTTCGCTACTT,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_GGACATGAGACAGAAC,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_GGCCTGGGTGTGCATG,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_GGCGCAACAAAGCTCT,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_GGGCAATCACCGACAC,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_GTATCCCGTAAAGCCT,cluster1 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_AAACCTCCAAGACAGC,cluster2 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_AACGGCGCAACTTCAA,cluster2 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_ACAAGTAAGCACAGGA,cluster2 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_ACCACTCGTACCTGCT,cluster2 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_ACTCATACAACCAAGA,cluster2 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_ACTGGGGCAAGACAGC,cluster2 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_AGCATCGGTCAATGCG,cluster2 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_AGCCCCTTCACCAACA,cluster2 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_AGCTGATGTACCTATT,cluster2 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_CAATCGGGTACAAGCT,cluster2 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_CATACCGCAATGGCGG,cluster2 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_CCCATTGTCACGCGCA,cluster2 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_CGCAATATCGGCGTAT,cluster2 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_CGCATAGGTGTAATTG,cluster2 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_CGCCAACCACCCACAA,cluster2 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_CGGTCATAGCACCTAA,cluster2 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_CGTAATTAGGTGTCCG,cluster2 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_CTATAAGGTACGGATT,cluster2 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_CTATAGCCACAATGTC,cluster2 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_GCAAAGTAGCGTGAAG,cluster2 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_GGAGTTGAGGTGCCCA,cluster2 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_GTCCATGCACCACTCT,cluster2 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_AACCGCAGTGCGCCTC,cluster3 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_AACTAAGAGCCAAAGG,cluster3 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_AAGTACGCACAGTTTA,cluster3 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_AATCAAGAGCATCTAG,cluster3 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_ACCTGCTGTAATTGGC,cluster3 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_ACTGCATCACTCGACA,cluster3 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_AGCCTTCGTAGGAAAT,cluster3 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_AGGGTATAGCTTAAGT,cluster3 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_CTCTCCGCACATAACC,cluster3 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_GCATACATCATGACCA,cluster3 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_GCTACATCACAATCAC,cluster3 +GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000_GGGAGATGTAGTCTCC,cluster3 diff --git a/inst/extdata/single_cell/extendedAnnotations_GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCell_multisample_chr9_1_1000000.rds b/inst/extdata/single_cell/extendedAnnotations_GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCell_multisample_chr9_1_1000000.rds new file mode 100644 index 00000000..5d11b9e9 Binary files /dev/null and b/inst/extdata/single_cell/extendedAnnotations_GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCell_multisample_chr9_1_1000000.rds differ diff --git a/inst/extdata/single_cell/quantData_GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCell_multisample_chr9_1_1000000.rds b/inst/extdata/single_cell/quantData_GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCell_multisample_chr9_1_1000000.rds new file mode 100644 index 00000000..063ef1e7 Binary files /dev/null and b/inst/extdata/single_cell/quantData_GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCell_multisample_chr9_1_1000000.rds differ diff --git a/inst/extdata/single_cell/seGene_GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCell_multisample_chr9_1_1000000.rds b/inst/extdata/single_cell/seGene_GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCell_multisample_chr9_1_1000000.rds new file mode 100644 index 00000000..0af6e2fb Binary files /dev/null and b/inst/extdata/single_cell/seGene_GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCell_multisample_chr9_1_1000000.rds differ diff --git a/inst/extdata/single_cell/seReadClass_GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCell_multisample_chr9_1_1000000.rds b/inst/extdata/single_cell/seReadClass_GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCell_multisample_chr9_1_1000000.rds new file mode 100644 index 00000000..d2bd26e1 Binary files /dev/null and b/inst/extdata/single_cell/seReadClass_GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCell_multisample_chr9_1_1000000.rds differ diff --git a/inst/extdata/single_cell/seTxClusteredEM_GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCell_multisample_chr9_1_1000000.rds b/inst/extdata/single_cell/seTxClusteredEM_GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCell_multisample_chr9_1_1000000.rds new file mode 100644 index 00000000..aae988be Binary files /dev/null and b/inst/extdata/single_cell/seTxClusteredEM_GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCell_multisample_chr9_1_1000000.rds differ diff --git a/inst/extdata/single_cell/seTxEM_GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCell_multisample_chr9_1_1000000.rds b/inst/extdata/single_cell/seTxEM_GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCell_multisample_chr9_1_1000000.rds new file mode 100644 index 00000000..b4c48a91 Binary files /dev/null and b/inst/extdata/single_cell/seTxEM_GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCell_multisample_chr9_1_1000000.rds differ diff --git a/inst/extdata/single_cell/seTxUnique_GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCell_multisample_chr9_1_1000000.rds b/inst/extdata/single_cell/seTxUnique_GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCell_multisample_chr9_1_1000000.rds new file mode 100644 index 00000000..2422cd6b Binary files /dev/null and b/inst/extdata/single_cell/seTxUnique_GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCell_multisample_chr9_1_1000000.rds differ diff --git a/man/bambu.Rd b/man/bambu.Rd index e3f4eca4..a32aacf7 100644 --- a/man/bambu.Rd +++ b/man/bambu.Rd @@ -9,10 +9,14 @@ bambu( annotations = NULL, genome = NULL, NDR = NULL, + mode = NULL, opt.discovery = NULL, + opt.rcAssignment = NULL, opt.em = NULL, + opt.singlecell = NULL, rcOutDir = NULL, discovery = TRUE, + assignDist = TRUE, quant = TRUE, stranded = FALSE, ncore = 1, @@ -20,8 +24,11 @@ bambu( trackReads = FALSE, returnDistTable = FALSE, lowMemory = FALSE, + sampleData = NULL, fusionMode = FALSE, - verbose = FALSE + verbose = FALSE, + quantData = NULL, + processByChromosome = FALSE ) } \arguments{ @@ -57,11 +64,8 @@ reconstruction process: annotation to be assigned to the same gene id, defaults to 10bp} \item{min.primarySecondaryDist}{specifying the minimum number of distance threshold, defaults to 5bp} - \item{min.primarySecondaryDistStartEnd1}{specifying the minimum number + \item{min.primarySecondaryDistStartEnd1}{specifying the minimum number of distance threshold, used for extending annotation, defaults to 5bp} - \item{min.primarySecondaryDistStartEnd2}{specifying the minimum number - of distance threshold, used for estimating distance to annotation, - defaults to 5bp} \item{min.txScore.multiExon}{specifying the minimum transcript level threshold for multi-exon transcripts during sample combining, defaults to 0} @@ -86,6 +90,18 @@ reconstruction process: defaults to "Bambu"} }} +\item{opt.rcAssignment}{A list of controlling parameters for the read class +to transcript assignment process: +\describe{ + \item{min.exonDistance}{specifying minimum distance to known transcript + to be considered a valid match, defaults to 35bp} + \item{min.primarySecondaryDist}{specifying the minimum distance + threshold between primary and secondary assignments, defaults to 5bp} + \item{min.primarySecondaryDistStartEnd2}{specifying the minimum + distance threshold for start/end positions used for read assignment, + defaults to 5bp} +}} + \item{opt.em}{A list of controlling parameters for quantification algorithm estimation process: \describe{ @@ -96,16 +112,23 @@ algorithm estimation process: defaults to 0.0001} \item{minvalue}{specifying the minvalue for convergence consideration, defaults to 0.00000001} - \item{sig.digit}{specifying the maximum significant digits of the reported estimates} + \item{sig.digit}{specifying the maximum significant digits of the reported estimates} }} +\item{opt.singlecell}{A list of parameters for bambu's single-cell module +(\code{extractBarcodeUMI}, \code{dedupUMI}, and \code{clusters}). See +\code{\link{bambu.singlecell}} for details.} + \item{rcOutDir}{A string variable specifying the path to where read class files will be saved.} \item{discovery}{A logical variable indicating whether annotations are to be extended. Defaults to TRUE} -\item{quant}{A logical variable indicating whether quantification will +\item{assignDist}{A logical variable indicating whether read-class-to-transcript +assignment will be performed, defaults to TRUE} + +\item{quant}{A logical variable indicating whether quantification will be performed. If false the output type will change. Defaults to TRUE} \item{stranded}{A boolean for strandedness, defaults to FALSE.} @@ -127,6 +150,12 @@ distTables. The output is a list with an entry for each sample.} \item{lowMemory}{Read classes will be processed by chromosomes when lowMemory is specified. This option provides an efficient way to process big samples.} +\item{sampleData}{A character vector of paths to metadata CSV files (or \code{NA} if +unavailable for specific samples); defaults to \code{NULL}. Files must contain a +"sampleName" column for bulk data or a "barcode" column for single-cell/spatial data. +For bulk data, one metadata CSV file for all samples is sufficient, whereas single-cell/spatial +data requires one metadata CSV file per sample.} + \item{fusionMode}{A logical variable indicating whether run in fusion mode} \item{verbose}{A logical variable indicating whether processing messages will diff --git a/man/bambu.singlecell.Rd b/man/bambu.singlecell.Rd new file mode 100644 index 00000000..06cc6d9d --- /dev/null +++ b/man/bambu.singlecell.Rd @@ -0,0 +1,234 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/bambu.R +\name{bambu.singlecell} +\alias{bambu.singlecell} +\title{Single-cell / Spatial transcript discovery and quantification with Bambu} +\usage{ +bambu.singlecell( + reads, + output, + annotations = NULL, + genome = NULL, + clusters = NULL, + NDR = NULL, + stranded = FALSE, + ncore = 1, + ... +) +} +\arguments{ +\item{reads}{For \code{output} in \code{"readClasses"}, +\code{"extendedAnnotations"}, or \code{"quantData"}, the input reads as in +\code{\link{bambu}}: paths to BAM files, a \code{BamFileList}, or the read-class +objects returned by the \code{"readClasses"} stage. For \code{output} in +\code{"uniqueCounts"}, \code{"EM"}, or \code{"clusteredEM"}, the per-sample +\code{quantData} list returned by the \code{"quantData"} stage.} + +\item{output}{Bambu performs transcript discovery and quantification on single-cell +data in distinct steps, that are each specified using the \code{output} option. Must +be one of \code{"readClasses"}, \code{"extendedAnnotations"}, \code{"quantData"}, +\code{"uniqueCounts"}, \code{"EM"}, or \code{"clusteredEM"}. See the \strong{Details} +section for what each stage does and the \strong{Value} section for what it returns.} + +\item{annotations}{A path to a .gtf file or a \code{TxDb} object +or a GRangesList object obtained by \code{\link{prepareAnnotations}}.} + +\item{genome}{A path to a fasta file or a BSGenome object.} + +\item{clusters}{Assignment of cells or spatial spots to clusters, supplied as one of: +a named vector whose names are the cell/spot identifiers and whose values are the +cluster labels; a \code{data.frame} object with an \code{id} and a \code{cluster} +column; or a path to a \code{.csv}/\code{.tsv}/\code{.txt} file holding +that \code{data.frame}. In every case the identifier follows the +\code{sampleName_barcode} format and the cluster label is the cluster that cell/spot +belongs to. Must be used with \code{output = "clusteredEM"}, where expression is +aggregated to the cluster level instead of per individual cell/spot.} + +\item{NDR}{specifying the maximum NDR rate to novel transcript +output from detected read classes, defaults to an automatic recommendation} + +\item{stranded}{A boolean for strandedness, defaults to FALSE.} + +\item{ncore}{specifying number of cores used when parallel processing +is used, defaults to 1.} + +\item{...}{Additional arguments passed to \code{\link{bambu}}, such as +\code{sampleData}, \code{opt.em}, \code{trackReads}, \code{returnDistTable}, +\code{yieldSize}, and \code{verbose}. See \code{\link{bambu}} for the full set.} +} +\value{ +The object returned depends on the specified \code{output} argument: +\describe{ + \item{\code{output = "readClasses"}}{returns a list of read-class objects for + each sample that represent a summarised representation of similar reads with + cell-barcode-to-read mapping information. ReadClass objects can be used as input + for transcript discovery and read to transcript assignment.} + \item{\code{output = "extendedAnnotations"}}{returns the extended annotations + after transcript discovery as a \code{GRangesList} object. The extended + annotation object can be used as input for transcript to read assignment and + quantification.} + \item{\code{output = "quantData"}}{returns a list of \code{quantData} objects + that contain the read-to-transcript assignments for each sample. The quantData + object is used to calculate unique counts per transcript and gene-level counts, + and also can be used as input for EM quantification.} + \item{\code{output = "uniqueCounts"}}{returns a \code{SummarizedExperiment} of + transcript-level unique counts (reads compatible with a single transcript), with + one column per cell. Unique counts are ready for downstream analysis as is. + Collapse them to gene-level counts with \code{\link{transcriptToGeneExpression}}, + or use them to cluster cells by gene expression profile (see \code{output = "clusteredEM"}).} + \item{\code{output = "EM"}}{returns a \code{SummarizedExperiment} with transcript + counts estimated using Bambu's EM for each barcode/single cell. This step is + resource intensive, consider using unique counts or the clustered EM.} + \item{\code{output = "clusteredEM"}}{returns a \code{SummarizedExperiment} with + one column per cluster as specified in the cluster argument. \code{colData} + describes the \code{cluster} label and the set of cell barcodes belonging to the + cluster. Quantifying transcript expression for cell clusters provides more stable + EM estimates.} +} +} +\description{ +A function wrapper for \code{\link{bambu}} that runs on long-read RNA-seq data to +perform transcript discovery and quantification at the single-cell/spatial level. +The function calls \code{\link{bambu}} with cell +barcode/UMI extraction and UMI-based deduplication enabled by default, with an optional +parameter for clustered-level transcript quantification. + +For general users, we recommend to use bambu-pipe, an end-to-end single-cell/spatial pipeline that +handles preprocessing, barcode demultiplexing, and alignment prior to +transcript discovery and quantification. See \url{https://github.com/GoekeLab/bambu-pipe}. +} +\details{ +Single-cell and spatial long read RNA-Seq protocols tag every read with a barcode +and a unique molecular identifier (UMI). The barcode identifies which cell +(single-cell data) or spatial spot (spatial data) a read came from; below, "cell" +and "single-cell" refer to either. \code{\link{bambu.singlecell}} reads the barcode +and UMI for each read and collapses duplicate molecules, so that transcript +quantification can be performed per barcode. + +\code{\link{bambu.singlecell}} runs one stage at a time, each specified with the +\code{output} argument for the returned object. The step-by-step workflow is +described below and in the examples: +\enumerate{ + \item \strong{Read Class Construction} (\code{output = "readClasses"}): build + per-sample read classes from the demultiplexed reads. Requires \code{reads} + (paths to BAM files), \code{genome} (reference genome sequence), and + \code{annotations} (optional, but recommended). + \item \strong{Transcript Discovery} (\code{output = "extendedAnnotations"}): + discover novel transcripts and extend the reference annotations. Requires + \code{reads} (the read classes from step 1) and \code{annotations} (optional, + but recommended; use the same annotations as in step 1, or omit them in both). + \item \strong{Read to Transcript Assignment} (\code{output = "quantData"}): + assign reads to the transcripts, giving the per-sample read-to-transcript + assignments. Requires \code{reads} (the read classes from step 1) and + \code{annotations} (the extended annotations from step 2 or the original + annotations, which cannot be NULL if transcript discovery was skipped). + \item \strong{Transcript Quantification} with the Expectation-Maximization + (EM) algorithm, or without it (EM-free): + \itemize{ + \item \strong{EM-free unique counts} (\code{output = "uniqueCounts"}): + transcript-level unique counts, ready for downstream analysis. Apply + \code{\link{transcriptToGeneExpression}} for gene counts, or cluster cells + by gene expression to run the clustered EM. Requires \code{reads} (the + quantData list from step 3) and the same \code{annotations} used in step 3. + \item \strong{Single-Cell EM} (\code{output = "EM"}): estimates + single-cell transcript expression with the EM. Requires \code{reads} (the + quantData list from step 3) and the same \code{annotations} used in step 3. + Single-cell EM is resource intensive; for higher efficiency use \code{uniqueCounts}, + or use \code{clusteredEM} to run the EM directly on cell clusters for more + accurate, stable estimates. + \item \strong{Clustered EM} (\code{output = "clusteredEM"}): runs the EM + directly on user-defined cell clusters, giving more stable estimates than + quantifying single cells and clustering afterwards. Requires \code{reads} + (the quantData list from step 3), the same \code{annotations} used in + step 3, and \code{clusters}. + } +} +This makes it practical to run different stages in different R sessions to minimise +memory usage, and also lets users reuse intermediate results and explore different +downstream settings without reprocessing large datasets from the raw alignments. See +the \strong{Value} section for the returned output from each stage and \strong{Examples} +for a detailed step-by-step run. +} +\examples{ +## single-cell example data: two demultiplexed samples (PacBio and ONT) +sce.dir <- system.file("extdata", "single_cell", package = "bambu") +reads <- file.path(sce.dir, c( + "GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNAPacBio_Rep1_Run2_demultiplexed_chr9_1_1000000.bam", + "GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCellcDNA_Rep1_Run1_demultiplexed_chr9_1_1000000.bam")) +annotations <- system.file("extdata", + "Homo_sapiens.GRCh38.91_chr9_1_1000000.gtf", package = "bambu") +genome <- system.file("extdata", + "Homo_sapiens.GRCh38.dna_sm.primary_assembly_chr9_1_1000000.fa", + package = "bambu") + +## Each stage is selected with an 'output' preset; its output feeds the next. +## 1. read-class construction from the demultiplexed BAMs +readClasses <- bambu.singlecell(reads = reads, output = "readClasses", + annotations = annotations, genome = genome) + +## 2. extend the annotations with novel transcripts discovered by Bambu +extendedAnnotations <- bambu.singlecell(reads = readClasses, + output = "extendedAnnotations", annotations = annotations) + +## 3. assign reads to transcripts (per-sample quantData list) +quantData <- bambu.singlecell(reads = readClasses, + output = "quantData", annotations = extendedAnnotations) + +## 4. EM-free quantification. Steps 4a-4b give unique counts and gene counts ready +## for downstream analysis; stop here, or continue to step 5 for EM-based +## quantification. + +## 4a. transcript-level unique counts (reads compatible with a single transcript). +## These EM-free counts are ready to use downstream; collapse them to gene counts, +## or cluster cells by gene expression to run the clusteredEM step. +uniqueCountsSe <- bambu.singlecell(reads = quantData, + output = "uniqueCounts", annotations = extendedAnnotations) + +## 4b. gene-level counts: transcriptToGeneExpression() sums the transcript counts +## within each gene to give a gene-level SummarizedExperiment. +geneCountsSe <- transcriptToGeneExpression(uniqueCountsSe) + +## 5. (optional) EM-based quantification. Steps 5a-5b estimate transcript expression +## with the EM, either at single-cell resolution (5a) or per user-defined cluster +## for more stable estimates (5b). + +## 5a. single-cell transcript counts estimated with the EM. +EMCountsSe <- bambu.singlecell(reads = quantData, + output = "EM", annotations = extendedAnnotations) + +## 5b. per-cluster (pseudobulk) transcript counts estimated with the EM. +## How the clusters are defined is up to the analysis. One common approach is to +## cluster by cell-type using the gene-level counts from step 4 with Seurat (not a +## bambu dependency, so this block runs only when Seurat is available): +if (requireNamespace("Seurat", quietly = TRUE)) { + seurat <- Seurat::CreateSeuratObject( + counts = assays(geneCountsSe)$counts, min.cells = 1) + seurat <- Seurat::NormalizeData(seurat) + ## this example is small, so cluster on all genes with exact PCA and a small + ## number of components (npcs must stay below the gene count) + seurat <- Seurat::ScaleData(seurat, features = rownames(seurat)) + seurat <- Seurat::RunPCA(seurat, features = rownames(seurat), + npcs = 10, approx = FALSE) + seurat <- Seurat::FindNeighbors(seurat, dims = 1:10) + seurat <- Seurat::FindClusters(seurat, resolution = 0.05) + + ## Seurat's active.ident is a named vector (names = ids in sampleName_barcode + ## format, values = cluster labels), one of the accepted input data type for + ## `clusters` argument in `bambu.singlecell()` + clusters <- seurat@active.ident +} + +## In general, clusteredEM accepts the clusters argument as a named vector (e.g. +## seurat@active.ident), a data.frame, or a path to a .csv/.tsv/.txt file. Here we +## illustrate the input data type with a stored .csv: +clusters <- file.path(sce.dir, + "clusters_GIS_cellMix_HepG2-A549-H9-HEYA8_5primeSingleCell_multisample_chr9_1_1000000.csv") +clusteredEMCountsSe <- bambu.singlecell(reads = quantData, + output = "clusteredEM", annotations = extendedAnnotations, + clusters = clusters) +} +\seealso{ +\code{\link{bambu}} for the underlying function and the full +parameter set; \code{\link{transcriptToGeneExpression}} for preparing gene-level counts +for cell clustering and clustered quantification, or for direct use in downstream analysis. +} diff --git a/man/getUniqueCountsSe.Rd b/man/getUniqueCountsSe.Rd new file mode 100644 index 00000000..51b0166d --- /dev/null +++ b/man/getUniqueCountsSe.Rd @@ -0,0 +1,62 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/transcriptToGeneExpression.R +\name{getUniqueCountsSe} +\alias{getUniqueCountsSe} +\title{Unique-count SummarizedExperiment from Bambu read-to-transcript assignments} +\usage{ +getUniqueCountsSe(quantData, annotations) +} +\arguments{ +\item{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}}.} + +\item{annotations}{A \code{GRangesList} of transcript annotations matching the ones used to +produce \code{quantData}, typically the extended annotations from transcript discovery.} +} +\value{ +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} +} +} +\description{ +This function is intended to be used after the transcript +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. +} +\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) +} +\seealso{ +\code{\link{bambu.singlecell}} and \code{\link{bambu}} for producing \code{quantData}; +\code{\link{transcriptToGeneExpression}} to collapse the result to gene-level counts. +}