hashedDrops {DropletUtils} | R Documentation

## Demultiplex cell hashing data

### Description

Demultiplex cell barcodes into their samples of origin based on the most significant hash tag oligo (HTO). Also identify potential doublets based on the presence of multiple significant HTOs.

### Usage

    hashedDrops(
      x,
      ambient = NULL,
      min.prop = 0.05,
      pseudo.count = 5,
      doublet.nmads = 3,
      doublet.min = 2,
      doublet.mixture = FALSE,
      confident.nmads = 3,
      confident.min = 2
    )

### Arguments

`x` | A numeric/integer matrix-like object containing UMI counts. Rows correspond to HTOs and columns correspond to cell barcodes. Each barcode is assumed to correspond to a cell, i.e., cell calling is assumed to have already been performed.
`ambient` | A numeric vector of length equal to `nrow(x)`, specifying the relative abundance of each HTO in the ambient solution - see Details.
`min.prop` | Numeric scalar in (0, 1) specifying the expected minimum proportion of barcodes contributed by each sample. Only used for estimating the ambient profile when `ambient=NULL`.
`pseudo.count` | A numeric scalar specifying the minimum pseudo-count when computing log-fold changes.
`doublet.nmads` | A numeric scalar specifying the number of median absolute deviations (MADs) to use to identify doublets.
`doublet.min` | A numeric scalar specifying the minimum threshold on the log-fold change to use to identify doublets.
`doublet.mixture` | Logical scalar indicating whether to use a 2-component mixture model to identify doublets.
`confident.nmads` | A numeric scalar specifying the number of MADs to use to identify confidently assigned singlets.
`confident.min` | A numeric scalar specifying the minimum threshold on the log-fold change to use to identify singlets.

### Details

Note that this function is still experimental; feedback is welcome.

The idea behind cell hashing is that cells from the same sample are stained with reagent conjugated with a single HTO. Cells are mixed across multiple samples and subjected to droplet-based single-cell sequencing. Cell barcode libraries can then be demultiplexed into individual samples based on whether their unique HTO is detected.

We identify the sample of origin for each cell barcode as that corresponding to the most abundant HTO. (See below for more details on exactly how “most abundant” is defined.) The log-fold change between the largest and second-largest abundances is also reported for each barcode, with large log-fold changes representing confident assignment to a single sample. We also report the log-fold change of the second-most abundant HTO over the estimated level of ambient contamination. Large log-fold changes indicate that the second HTO has greater abundance than expected, consistent with a doublet.

To facilitate quality control, we explicitly identify problematic barcodes as outliers on the relevant metrics.

- By default, we identify putative doublets as those with `LogFC2` values that are (i) `doublet.nmads` MADs above the median and (ii) greater than `doublet.min`. The hard threshold is more-or-less arbitrary and aims to avoid overly aggressive detection of large outliers in a naturally right-skewed distribution (given that the log-fold changes are positive by definition, and most of the distribution is located near zero).

- Alternatively, if `doublet.mixture=TRUE`, we fit a two-component mixture model to the `LogFC2` distribution. Doublets are identified as all members of the component with the larger mean. This avoids the need for the arbitrary parameters mentioned above but only works when there are many doublets, otherwise both components will be fitted to the non-doublet values. (Initialization of the model assumes at least 5% doublets.)

Of the non-doublet libraries, we consider them to be confidently assigned to a single sample if their `LogFC` values are (i) not less than `confident.nmads` MADs below the median and (ii) greater than `confident.min`. The hard threshold is again arbitrary, but this time it aims to avoid insufficiently aggressive outlier detection - typically from an inflation of the MAD when the `LogFC` values are large, positive and highly variable.

### Value

A DataFrame with one row per column of `x`, containing the following fields:

- `Total`, integer specifying the total count across all HTOs for each barcode.

- `Best`, integer specifying the HTO with the highest abundance for each barcode.

- `Second`, integer specifying the HTO with the second-highest abundance.

- `LogFC`, numeric containing the log-fold change between the abundances of the best and second-best HTO.

- `LogFC2`, numeric containing the log-fold change in the second-best HTO over the ambient contamination.

- `Doublet`, logical specifying whether a barcode is a doublet.

- `Confident`, logical specifying whether a barcode is a confidently assigned singlet.

In addition, the metadata contains `ambient`, a numeric vector containing the (estimate of the) ambient profile; `doublet.threshold`, the threshold applied to `LogFC2` to identify doublets; and `confident.threshold`, the threshold applied to non-doublet `LogFC` values to identify confident singlets.

### Adjusting abundances for ambient contamination

HTO abundances require some care to compute due to the presence of ambient contamination in each library. Ideally, the experiment would be performed in such a manner that the concentration of each HTO is the same. However, if one HTO is present at higher concentration in the ambient solution, this might incorrectly cause us to assign all barcodes to the corresponding sample.

To adjust for ambient contamination, we assume that the ambient contamination in each library follows the same profile as `ambient`. We further assume that a minority of HTOs in a library are actually driven by the presence of cell(s), the rest coming from the ambient solution. We estimate the level of ambient contamination in each barcode by scaling `ambient`, using a DESeq-like normalization algorithm to compute the scaling factor. (The requisite assumption of a non-DE majority follows from the two assumptions above.) We then subtract the scaled ambient proportions from the HTO count profile to remove the effect of contamination. Abundances that would otherwise be negative are set to zero.

For experiments with 3-4 HTOs, we assume that higher-order multiplets are negligible and define the scaling factor as the third-largest ratio between the HTO counts and `ambient`. For experiments with only 2 HTOs, doublet detection is not possible as the second-most abundant HTO is always used to estimate the ambient contamination.

### Getting the ambient proportions

Ideally, `ambient` would be obtained from libraries that do not correspond to cell-containing droplets. For example, we could get this information from the `metadata` of the `emptyDrops` output, had we run `emptyDrops` on the HTO count matrix (see below).

Unfortunately, in some cases (e.g., public data), counts are provided for only the cell-containing barcodes. If `ambient=NULL`, the function will fit a two-component mixture model to each HTO's count distribution. All barcodes assigned to the lower component are considered to have background counts for that HTO, and the mean of those counts is used as an estimate of the ambient contribution.

The initialization of the mixture model is controlled by `min.prop`, which starts the means of the lower and upper components at the `min.prop` and `1-min.prop` quantiles, respectively. This means that each sample is expected to contribute somewhere between `[min.prop, 1-min.prop]` barcodes. Larger values improve convergence but require stronger assumptions about the relative proportions of multiplexed samples.

### Computing the log-fold changes

After subtraction of the ambient noise but before calculation of the log-fold changes, we need to add a pseudo-count to ensure that the log-fold changes are well-defined. We set the pseudo-count to the average ambient HTO count (i.e., the average of the scaled `ambient`), effectively exploiting the ambient contamination as a natural pseudo-count that scales with barcode-specific capture efficiency and sequencing depth. (In libraries with low sequencing depth, we still enforce a minimum pseudo-count of `pseudo.count`.)

This scaling behavior is useful as it ensures that shrinkage of the log-fold changes is not more severe for libraries that have not been sequenced as deeply. We thus avoid excessive variability in the log-fold change distribution that would otherwise reduce the precision of outlier detection. The implicit assumption is that the number of contaminating transcript molecules is roughly the same in each droplet, meaning that any differences in ambient coverage between libraries reflect technical biases that would also affect cell-derived HTO counts.

Another nice aspect of this entire procedure (subtraction and re-addition) is that it collapses to a no-op if the experiment is well-executed with identical concentrations of all HTOs in the ambient solution.

### Use only on non-empty droplets

This function assumes that cell calling has already been performed, e.g., with `emptyDrops`. Specifically, `x` should only contain columns corresponding to non-empty droplets. If empty droplets are included, their log-fold changes will simply reflect stochastic sampling in the ambient solution and violate the assumptions involved in outlier detection.

If `x` contains columns for both empty and non-empty droplets, it is straightforward to simply run `emptyDrops` on the HTO count matrix to identify the latter. Note that some fiddling with the `lower=` argument may be required, depending on the sequencing depth of the HTO libraries.

### Author(s)

Aaron Lun

### References

Stoeckius M, Zheng S, Houck-Loomis B et al. (2018) Cell Hashing with barcoded antibodies enables multiplexing and doublet detection for single cell genomics. Genome Biol. 19, 1:224

### See Also

`emptyDrops`, to identify which barcodes are likely to contain cells.

### Examples

    # Mocking up an example dataset with 10 HTOs and 10% doublets.
    ncells <- 1000
    nhto <- 10
    y <- matrix(rpois(ncells*nhto, 50), nrow=nhto)
    true.sample <- sample(nhto, ncells, replace=TRUE)
    y[cbind(true.sample, seq_len(ncells))] <- 1000

    ndoub <- ncells/10
    next.sample <- (true.sample[1:ndoub]  + 1) %% nrow(y)
    next.sample[next.sample==0] <- nrow(y)
    y[cbind(next.sample, seq_len(ndoub))] <- 500

    # Computing hashing statistics.
    stats <- hashedDrops(y)

    # Doublets show up in the top-left,
    # singlets in the bottom right.
    plot(stats$LogFC, stats$LogFC2)

    # Most cells should be singlets with low NMAD.
    hist(stats$LogFC2, breaks=50)

    # Identify confident singlets or doublets at the given NMAD threshold.
    summary(stats$Confident)
    summary(stats$Doublet)

    # Chcecking against the known truth, in this case
    # 'Best' contains the putative sample of origin.
    table(stats$Best, true.sample)

[Package DropletUtils version 1.8.0 Index]