In hbc/bcbioSinglecell: Bcbio Single-Cell RNA-Seq

## nolint start
suppressPackageStartupMessages({
    library(goalie)
    library(basejump)
    library(ggplot2)
    library(bcbioSingleCell)
})
prepareTemplate()
source("_setup.R")
## nolint end

Load bcbioSingleCell object

object <- import(params[["bcb_file"]])
assert(
    is(object, "bcbioSingleCell"),
    validObject(object)
)
print(object)

[bcbio][] run data was imported from:
r metadata(object)[["uploadDir"]].

Sample metadata

sampleData(object)

Reads per cell {.tabset}

These are counts of how many reads are assigned to a given cellular barcode. It is normal for single cell RNA-seq data to contain a large number of low complexity barcodes. The bcbio pipeline filters out most of these barcodes, and here we have applied a threshold cutoff of a minimum of r metadata(object)[["cellularBarcodeCutoff"]] reads per cell. The unfiltered read count distributions are shown here.

Histogram

For high quality data, the proportional histogram should contain a single large peak that represents cells that were encapsulated. If we see a strong shoulder, or a bimodal distribution of the cells, that can indicate a couple problems. It might be that there is free floating RNA, which happens when cells are dying. It could also be that there are a set of cells that failed for some reason. Finally, it could also be that there are biologically different types of cells, and one type is much smaller than the other. If this is the case we would expect to see less RNA being sequenced from the smaller cells.

plotReadsPerCell(
    object = object,
    geom = "histogram",
    interestingGroups = "sampleName"
)

ECDF

An empirical distribution function (ECDF) plot will show the frequency distribution of the reads per cell. You can see that the vast majority of low complexity barcodes plateau at a read depth below 1000 reads per cell.

plotReadsPerCell(object = object, geom = "ecdf")

UMI counts per cell {.tabset}

Now let's assess the distribution of unique molecular identifier (UMI)-deconvoluted counts per cell. In general, the distributions should be relatively uniform per sample. Here we are also including violin and ridgeline plots, with the average number of genes per cell labeled.

markdownHeader("Violin", level = 2L)
plotCountsPerCell(
    object = object,
    geom = "violin",
    min = params[["min_counts"]],
    max = params[["max_counts"]]
)

markdownHeader("Ridgeline", level = 2L)
plotCountsPerCell(
    object = object,
    geom = "ridgeline",
    min = params[["min_counts"]],
    max = params[["max_counts"]]
)

markdownHeader("Histogram", level = 2L)
plotCountsPerCell(
    object = object,
    geom = "histogram",
    min = params[["min_counts"]],
    max = params[["max_counts"]]
)

markdownHeader("ECDF", level = 2L)
plotCountsPerCell(
    object = object,
    geom = "ecdf",
    interestingGroups = "sampleName",
    min = params[["min_counts"]],
    max = params[["max_counts"]]
)

Filter cells by UMI count

Let's apply this step first and then proceed to evaluating gene detection, mitocondrial transcript abundance, and novelty scores.

object <- filterCells(
    object = object,
    minCounts = params[["min_counts"]],
    maxCounts = params[["max_counts"]]
)

Let's take a look at the UMI per cell distributions after this filtering step. Note that we haven't applied very strict filtering here -- we're going to cut off the "low quality" cells based on the gene detection rate, novelty score, and mitochondrial abundance.

plotCountsPerCell(
    object = object,
    geom = "histogram",
    min = params[["min_counts"]],
    max = params[["max_counts"]]
)

Genes detected per cell {.tabset}

Here by "detected", we mean genes with a non-zero count measurement per cell. Seeing gene detection in the range of 500-5000 is normal for most single-cell experiments.

markdownHeader("Violin", level = 2L)
plotFeaturesPerCell(
    object = object,
    geom = "violin",
    min = min(params[["min_features"]]),
    max = max(params[["max_features"]])
)

markdownHeader("Ridgeline", level = 2L)
plotFeaturesPerCell(
    object = object,
    geom = "ridgeline",
    min = min(params[["min_features"]]),
    max = max(params[["max_features"]])
)

markdownHeader("Histogram", level = 2L)
plotFeaturesPerCell(
    object = object,
    geom = "histogram",
    min = min(params[["min_features"]]),
    max = max(params[["max_features"]])
)

markdownHeader("ECDF", level = 2L)
plotFeaturesPerCell(
    object = object,
    geom = "ecdf",
    min = min(params[["min_features"]]),
    max = max(params[["max_features"]])
)

UMIs vs. features detected

If we graph out the total number of UMI counts per cell vs. the genes detected per cell, we can assess whether there is a large population of low quality cells with low counts and/or gene detection.

plotCountsVsFeatures(object)

Novelty score {.tabset}

Another way to QC the data is to look for less novelty, that is cells that have less genes detected per count than other cells. We can see the samples where we sequenced each cell less have a higher overall novelty, that is because we have not started saturated the sequencing for any given gene for these samples. Outlier cells in these samples might be cells that we have a less complex RNA species than other cells. Sometimes we can detect contamination with low complexity cell types like red blood cells via this metric.

markdownHeader("Violin", level = 2L)
plotNovelty(
    object = object,
    geom = "violin",
    min = min(params[["min_novelty"]])
)

markdownHeader("Ridgeline", level = 2L)
plotNovelty(
    object = object,
    geom = "ridgeline",
    min = min(params[["min_novelty"]])
)

markdownHeader("Histogram", level = 2L)
plotNovelty(
    object = object,
    geom = "histogram",
    min = min(params[["min_novelty"]])
)

markdownHeader("ECDF", level = 2L)
plotNovelty(
    object = object,
    geom = "ecdf",
    min = min(params[["min_novelty"]])
)

Mitochondrial abundance {.tabset}

We evaluate overall mitochondrial gene expression as a biomarker of cellular stress during sample preparation.

markdownHeader("Violin", level = 2L)
plotMitoRatio(
    object = object,
    geom = "violin",
    max = max(params[["max_mito_ratio"]])
)

markdownHeader("Ridgeline", level = 2L)
plotMitoRatio(
    object = object,
    geom = "ridgeline",
    max = max(params[["max_mito_ratio"]])
)

markdownHeader("Histogram", level = 2L)
plotMitoRatio(
    object = object,
    geom = "histogram",
    max = max(params[["max_mito_ratio"]])
)

markdownHeader("ECDF", level = 2L)
plotMitoRatio(
    object = object,
    geom = "ecdf",
    max = max(params[["max_mito_ratio"]])
)

Filter cells

object <- filterCells(
    object = object,
    minCounts = params[["min_counts"]],
    maxCounts = params[["max_counts"]],
    minFeatures = params[["min_features"]],
    maxFeatures = params[["max_features"]],
    maxMitoRatio = params[["max_mito_ratio"]],
    minNovelty = params[["min_novelty"]],
    nCells = params[["n_cells"]],
    minCellsPerFeature = params[["min_cells_per_feature"]]
)

plotQc(object, geom = "violin")

Save filtered data

name <- basenameSansExt(params[["bcb_file"]])
assignAndSaveData(
    name = paste(name, "filtered", sep = "_"),
    object = object,
    dir = params[["data_dir"]]
)

export(
    object = object,
    con = params[["output_dir"]],
    compress = TRUE
)

hbc/bcbioSinglecell documentation built on Jan. 14, 2024, 3:41 a.m.