Contents

1 Introduction

DecontPro assess and decontaminate single-cell protein expression data, such as those generated from CITE-seq or Total-seq. The count matrix is decomposed into three matrices, the native, the ambient and the background that represent the contribution from the true protein expression on cells, the ambient material and other non-specific background contamination.

2 Installation

DecontX Package can be installed from Bioconductor:

if (!requireNamespace("BiocManager", quietly = TRUE)) {
    install.packages("BiocManager")
}
BiocManager::install("decontX")

Then the package can be loaded in R using the following command:

library(decontX)

To see the latest updates and releases or to post a bug, see our GitHub page at https://github.com/campbio/decontX.

3 Importing data

Here we use an example dataset from SingleCellMultiModal package.

library(SingleCellMultiModal)
dat <- CITEseq("cord_blood", dry.run = FALSE)
#> Warning: 'ExperimentList' contains 'data.frame' or 'DataFrame',
#>   potential for errors with mixed data types
counts <- experiments(dat)$scADT

For this tutorial, we sample only 1000 droplets from the dataset to demonstrate the use of functions. When analyzing your dataset, sub-sampling should be done with caution, as decontPro approximates contamination profile using the dataset. A biased sampling may introduce bias to the contamination profile approximation.

set.seed(42)
sample_id <- sample(dim(counts)[2], 1000, replace = FALSE)
counts_sample <- counts[, sample_id]

4 Generate cell clusters

decontPro requires a vector indicating the cell types of each droplet. Here we use Seurat for clustering.

library(Seurat)
library(dplyr)
adt_seurat <- CreateSeuratObject(counts_sample, assay = "ADT")
#> Warning: Data is of class matrix. Coercing to dgCMatrix.
adt_seurat <- NormalizeData(adt_seurat, normalization.method = "CLR", margin = 2) %>%
  ScaleData(assay = "ADT") %>%
  RunPCA(assay = "ADT", features = rownames(adt_seurat), npcs = 10,
  reduction.name = "pca_adt") %>%
  FindNeighbors(dims = 1:10, assay = "ADT", reduction = "pca_adt") %>%
  FindClusters(resolution = 0.5)
#> Warning in irlba(A = t(x = object), nv = npcs, ...): You're computing too large
#> a percentage of total singular values, use a standard svd instead.
#> Warning: Requested number is larger than the number of available items (13).
#> Setting to 13.
#> Warning: Requested number is larger than the number of available items (13).
#> Setting to 13.
#> Warning: Requested number is larger than the number of available items (13).
#> Setting to 13.
#> Warning: Requested number is larger than the number of available items (13).
#> Setting to 13.
#> Warning: Requested number is larger than the number of available items (13).
#> Setting to 13.
#> Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
#> 
#> Number of nodes: 1000
#> Number of edges: 32498
#> 
#> Running Louvain algorithm...
#> Maximum modularity in 10 random starts: 0.8567
#> Number of communities: 9
#> Elapsed time: 0 seconds
adt_seurat <- RunUMAP(adt_seurat,
                      dims = 1:10,
                      assay = "ADT",
                      reduction = "pca_adt",
                      reduction.name = "adtUMAP",
                      verbose = FALSE)
#> Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
#> To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
#> This message will be shown once per session
DimPlot(adt_seurat, reduction = "adtUMAP", label = TRUE)

FeaturePlot(adt_seurat, 
            features = c("CD3", "CD4", "CD8", "CD19", "CD14", "CD16", "CD56"))


clusters <- as.integer(Idents(adt_seurat))

5 Run DecontPro

You can run DecontPro by simply:

counts <- as.matrix(counts_sample)
out <- deconPro(counts,
                clusters)

Priors (delta_sd and background_sd) may need tuning with the help of plotting the decontaminated results. The two priors encode belief if a more spread-out count should be considered contamination vs. native. We tested the default values on datasets ranging 5k to 10k droplets and 10 to 30 ADTs and the results are reasonable. A more contaminated or a smaller dataset may need larger priors. In this tutorial, since we sampled only 1k droplets from the original 10k droplets, we use slightly larger priors:

counts <- as.matrix(counts_sample)
out <- decontPro(counts,
                 clusters,
                 delta_sd = 2e-4,
                 background_sd = 2e-5)
#> Chain 1: ------------------------------------------------------------
#> Chain 1: EXPERIMENTAL ALGORITHM:
#> Chain 1:   This procedure has not been thoroughly tested and may be unstable
#> Chain 1:   or buggy. The interface is subject to change.
#> Chain 1: ------------------------------------------------------------
#> Chain 1: 
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Gradient evaluation took 0.019048 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 190.48 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Begin eta adaptation.
#> Chain 1: Iteration:   1 / 250 [  0%]  (Adaptation)
#> Chain 1: Iteration:  50 / 250 [ 20%]  (Adaptation)
#> Chain 1: Iteration: 100 / 250 [ 40%]  (Adaptation)
#> Chain 1: Iteration: 150 / 250 [ 60%]  (Adaptation)
#> Chain 1: Iteration: 200 / 250 [ 80%]  (Adaptation)
#> Chain 1: Success! Found best value [eta = 1] earlier than expected.
#> Chain 1: 
#> Chain 1: Begin stochastic gradient ascent.
#> Chain 1:   iter             ELBO   delta_ELBO_mean   delta_ELBO_med   notes 
#> Chain 1:    100   -513803659.665             1.000            1.000
#> Chain 1:    200     -7177548.141            35.792           70.585
#> Chain 1:    300      -342758.348            30.508           19.941
#> Chain 1:    400      -253644.497            22.969           19.941
#> Chain 1:    500      -223785.346            18.402            1.000
#> Chain 1:    600      -212411.740            15.344            1.000
#> Chain 1:    700      -206476.954            13.156            0.351
#> Chain 1:    800      -202104.000            11.514            0.351
#> Chain 1:    900      -199781.145            10.236            0.133
#> Chain 1:   1000      -198758.154             9.213            0.133
#> Chain 1:   1100      -198020.264             8.376            0.054   MAY BE DIVERGING... INSPECT ELBO
#> Chain 1:   1200      -197397.014             7.678            0.054   MAY BE DIVERGING... INSPECT ELBO
#> Chain 1:   1300      -196906.436             7.088            0.029   MAY BE DIVERGING... INSPECT ELBO
#> Chain 1:   1400      -196442.570             6.582            0.029   MAY BE DIVERGING... INSPECT ELBO
#> Chain 1:   1500      -196230.683             6.143            0.022   MAY BE DIVERGING... INSPECT ELBO
#> Chain 1:   1600      -195981.145             5.759            0.022   MAY BE DIVERGING... INSPECT ELBO
#> Chain 1:   1700      -195830.350             5.420            0.012   MAY BE DIVERGING... INSPECT ELBO
#> Chain 1:   1800      -195732.312             5.119            0.012   MAY BE DIVERGING... INSPECT ELBO
#> Chain 1:   1900      -195651.552             4.850            0.005   MEDIAN ELBO CONVERGED   MAY BE DIVERGING... INSPECT ELBO
#> Chain 1: 
#> Chain 1: Drawing a sample of size 1000 from the approximate posterior... 
#> Chain 1: COMPLETED.
#> Warning: Pareto k diagnostic value is 41.49. Resampling is disabled. Decreasing
#> tol_rel_obj may help if variational algorithm has terminated prematurely.
#> Otherwise consider using sampling instead.

The output contains three matrices, and model parameters after inference. decontaminated_counts represent the true protein expression on cells.

decontaminated_counts <- out$decontaminated_counts

6 Plot Results

Plot ADT density before and after decontamination. For bimodal ADTs, the background peak should be removed. Note CD4 is tri-modal with the intermediate mode corresponding to monocytes. This can be used as a QC metric for decontamination as only the lowest mode should be removed.

plotDensity(counts,
            decontaminated_counts,
            c("CD3", "CD4", "CD8", "CD14", "CD16", "CD19"))