Contents

1 Overview

scrapper implements bindings to C++ code for analyzing single-cell data, mostly from the libscran libraries. Each function performs an individual analysis step ranging from quality control to clustering and marker detection. scrapper is mostly intended for other Bioconductor package developers to build more user-friendly end-to-end workflows.

2 Quick start

Let’s fetch a small single-cell RNA-seq dataset for demonstration purposes:

library(scRNAseq)
sce.z <- ZeiselBrainData()
sce.z
## class: SingleCellExperiment 
## dim: 20006 3005 
## metadata(0):
## assays(1): counts
## rownames(20006): Tspan12 Tshz1 ... mt-Rnr1 mt-Nd4l
## rowData names(1): featureType
## colnames(3005): 1772071015_C02 1772071017_G12 ... 1772066098_A12
##   1772058148_F03
## colData names(9): tissue group # ... level1class level2class
## reducedDimNames(0):
## mainExpName: gene
## altExpNames(2): repeat ERCC

We run it through the scrapper analysis pipeline:

# Finding the mitochondrial genes for QC purposes.
is.mito <- grepl("^mt-", rownames(sce.z))

# We can set the number of threads higher in real applications, but
# we want to avoid stressing out Bioconductor's build system.
nthreads <- 2

library(scrapper)
res <- analyze(
    sce.z,
    rna.subsets=list(mito=is.mito),
    num.threads=nthreads
)

Now we can have a look at some of the results. For example, we can make a t-SNE plot:

tabulate(res$clusters)
##  [1] 289 257 199 392 192 321 205 177 152 318 124  34  77 137
plot(res$tsne[,1], res$tsne[,2], pch=res$clusters)

We can have a look at the markers defining each cluster:

# Top markers for each cluster, say, based on the median AUC:
top.markers <- lapply(res$rna.markers$auc, function(x) {
    head(rownames(x)[order(x$median, decreasing=TRUE)], 10)
})
head(top.markers)
## $`1`
##  [1] "Gad1"   "Gad2"   "Ndrg4"  "Stmn3"  "Tspyl4" "Scg5"   "Syp"    "Snap25"
##  [9] "Nap1l5" "Rab3a" 
## 
## $`2`
##  [1] "Plp1"   "Mbp"    "Cnp"    "Mobp"   "Mog"    "Sept4"  "Ugt8a"  "Trf"   
##  [9] "Gsn"    "Cldn11"
## 
## $`3`
##  [1] "Stmn3"  "Calm1"  "Scg5"   "Chn1"   "Snap25" "Pcsk2"  "Calm2"  "Mdh1"  
##  [9] "Nme1"   "Stmn2" 
## 
## $`4`
##  [1] "Chn1"          "Stmn3"         "Hpca"          "Calm2"        
##  [5] "Neurod6"       "3110035E14Rik" "Crym"          "Calm1"        
##  [9] "Ppp3r1"        "Ppp3ca"       
## 
## $`5`
##  [1] "Ywhah"  "Syp"    "Snap25" "Ndrg4"  "Vsnl1"  "Rab3a"  "Chn1"   "Stmn3" 
##  [9] "Uchl1"  "Ppp3r1"
## 
## $`6`
##  [1] "Atp1b1"        "Ppp3ca"        "Nsg2"          "Tspan13"      
##  [5] "Syp"           "Crym"          "Chn1"          "3110035E14Rik"
##  [9] "Thy1"          "Gria1"
# Aggregating all marker statistics for the first cluster
# (delta.mean is the same as the log-fold change in this case).
df1 <- reportGroupMarkerStatistics(res$rna.markers, 1)
df1 <- df1[order(df1$auc.median, decreasing=TRUE),]
head(df1[,c("mean", "detected", "auc.median", "delta.mean.median")])
##            mean  detected auc.median delta.mean.median
## Gad1   4.783640 1.0000000  0.9976632          4.587905
## Gad2   4.433167 0.9965398  0.9951791          4.293656
## Ndrg4  4.395829 0.9965398  0.9897291          3.755445
## Stmn3  4.713631 0.9930796  0.9869674          4.226543
## Tspyl4 3.364068 1.0000000  0.9831016          2.447678
## Scg5   4.750306 1.0000000  0.9765381          3.095245

Finally, we can convert the results into a SingleCellExperiment for interoperability with other Bioconductor packages like scater:

sce <- convertAnalyzeResults(res)
sce
## class: SingleCellExperiment 
## dim: 20006 2874 
## metadata(0):
## assays(2): filtered normalized
## rownames(20006): Tspan12 Tshz1 ... mt-Rnr1 mt-Nd4l
## rowData names(5): means variances fitted residuals is.highly.variable
## colnames(2874): 1772071015_C02 1772071017_G12 ... 1772063068_D01
##   1772066098_A12
## colData names(5): sum detected subsets.mito sizeFactor clusters
## reducedDimNames(3): pca tsne umap
## mainExpName: rna
## altExpNames(0):

3 Step-by-step

The analyze() function just calls a series of lower-level functions in scrapper. First, some quality control:

counts <- assay(sce.z)
rna.qc.metrics <- computeRnaQcMetrics(counts, subsets=list(mt=is.mito), num.threads=nthreads)
rna.qc.thresholds <- suggestRnaQcThresholds(rna.qc.metrics)
rna.qc.filter <- filterRnaQcMetrics(rna.qc.thresholds, rna.qc.metrics)
filtered <- counts[,rna.qc.filter,drop=FALSE]

Then normalization:

size.factors <- centerSizeFactors(rna.qc.metrics$sum[rna.qc.filter])
normalized <- normalizeCounts(filtered, size.factors)

Feature selection and PCA:

gene.var <- modelGeneVariances(normalized, num.threads=nthreads)
top.hvgs <- chooseHighlyVariableGenes(gene.var$statistics$residuals)
pca <- runPca(normalized[top.hvgs,], num.threads=nthreads)

Some visualization and clustering:

umap.out <- runUmap(pca$components, num.threads=nthreads)
tsne.out <- runTsne(pca$components, num.threads=nthreads)
snn.graph <- buildSnnGraph(pca$components, num.threads=nthreads)
clust.out <- clusterGraph(snn.graph)

And finally marker detection:

markers <- scoreMarkers(normalized, groups=clust.out$membership, num.threads=nthreads)

Readers are referred to the OSCA book for more details on the theory behind each step.

4 Dealing with multiple batches

Let’s fetch another brain dataset and combine it with the previous one.

sce.t <- TasicBrainData()
common <- intersect(rownames(sce.z), rownames(sce.t))
combined <- cbind(assay(sce.t)[common,], assay(sce.z)[common,])
block <- rep(c("tasic", "zeisel"), c(ncol(sce.t), ncol(sce.z)))

We call analyze() with block= to account for the batch structure. This ensures that the various calculations are not affected by inter-block differences. It also uses MNN correction to batch effects in the low-dimensional space prior to further analyses.

# No mitochondrial genes in the combined set, so we'll just skip it.
blocked_res <- analyze(combined, block=block, num.threads=nthreads)

# Visually check whether the datasets were suitably merged together.
# Note that these two datasets don't have the same cell types, so
# complete overlap should not be expected.
plot(blocked_res$tsne[,1], blocked_res$tsne[,2], pch=16, col=factor(blocked_res$block))