pairwiseBinom {scran} | R Documentation

## Perform pairwise binomial tests

### Description

Perform pairwise binomial tests between groups of cells, possibly after blocking on uninteresting factors of variation.

### Usage

    pairwiseBinom(x, groups, block = NULL, restrict = NULL,
      exclude = NULL, direction = c("any", "up", "down"),
      threshold = 1e-08, lfc = 0, log.p = FALSE,
      gene.names = rownames(x), clusters = NULL, subset.row = NULL,
      BPPARAM = SerialParam())

### Arguments

`x` | A numeric matrix-like object of counts, where each column corresponds to a cell and each row corresponds to a gene.
`groups` | A vector specifying the group assignment for each cell.
`block` | A factor specifying the blocking level for each cell.
`restrict` | A vector specifying the levels of `groups` for which to perform pairwise comparisons.
`exclude` | A vector specifying the levels of `groups` for which not to perform pairwise comparisons.
`direction` | A string specifying the direction of effects to be considered for the alternative hypothesis.
`threshold` | Numeric scalar specifying the value below which a gene is presumed to be not expressed.
`lfc` | Numeric scalar specifying the minimum absolute log-ratio in the proportion of expressing genes between groups.
`log.p` | A logical scalar indicating if log-transformed p-values/FDRs should be returned.
`gene.names` | A character vector of gene names with one value for each row of `x`.
`clusters` | Deprecated, same as `groups`.
`subset.row` | See `?"scran-gene-selection"`.
`BPPARAM` | A BiocParallelParam object indicating how parallelization should be performed across genes.

### Details

This function performs exact binomial tests to identify marker genes between pairs of groups of cells. Here, the null hypothesis is that the proportion of cells expressing a gene is the same between groups. A list of tables is returned where each table contains the statistics for all genes for a comparison between each pair of groups. This can be examined directly or used as input to `combineMarkers` for marker gene detection.

Effect sizes for each comparison are reported as log2-fold changes in the proportion of expressing cells in one group over the proportion in another group. We add a pseudo-count that squeezes the log-FCs towards zero to avoid undefined values when one proportion is zero. This is closely related to but somewhat more interpretable than the log-odds ratio, which would otherwise be the more natural statistic for a proportion-based test.

If `restrict` is specified, comparisons are only performed between pairs of groups in `restrict`. This can be used to focus on DEGs distinguishing between a subset of the groups (e.g., closely related cell subtypes). Similarly, if any entries of `groups` are `NA`, the corresponding cells are are ignored.

`x` can be a count matrix or any transformed counterpart where zeroes remain zero and non-zeroes remain non-zero. This is true of any scaling normalization and monotonic transformation like the log-transform. If the transformation breaks this rule, some adjustment of `threshold` is necessary.

A consequence of the transformation-agnostic behaviour of this function is that it will not respond to normalization. Differences in library size will not be considered by this function. However, this is not necessarily problematic for marker gene detection - users can treat this as retaining information about the total RNA content, analogous to spike-in normalization.

### Value

A list is returned containing `statistics` and `pairs`.

The `statistics` element is itself a list of DataFrames. Each DataFrame contains the statistics for a comparison between a pair of groups, including the overlap proportions, p-values and false discovery rates.

The `pairs` element is a DataFrame with one row corresponding to each entry of `statistics`. This contains the fields `first` and `second`, specifying the two groups under comparison in the corresponding DataFrame in `statistics`.

In each DataFrame in `statistics`, the log-fold change represents the log-ratio of the proportion of expressing cells in the `first` group compared to the expressing proportion in the `second` group.

### Direction and magnitude of the effect

If `direction="any"`, two-sided binomial tests will be performed for each pairwise comparisons between groups of cells. For other `direction`, one-sided tests in the specified direction will be used instead. This can be used to focus on genes that are upregulated in each group of interest, which is often easier to interpret.

In practice, the two-sided test is approximated by combining two one-sided tests using a Bonferroni correction. This is done for various logistical purposes; it is also the only way to combine p-values across blocks in a direction-aware manner. As a result, the two-sided p-value reported here will not be the same as that from `binom.test`. In practice, they are usually similar enough that this is not a major concern.

To interpret the setting of `direction`, consider the DataFrame for group X, in which we are comparing to another group Y. If `direction="up"`, genes will only be significant in this DataFrame if they are upregulated in group X compared to Y. If `direction="down"`, genes will only be significant if they are downregulated in group X compared to Y. See `?binom.test` for more details on the interpretation of one-sided Wilcoxon rank sum tests.

The magnitude of the log-fold change in the proportion of expressing cells can also be tested by setting `lfc`. By default, `lfc=0` meaning that we will reject the null upon detecting any difference in proportions. If this is set to some other positive value, the null hypothesis will change depending on `direction`:

- If `direction="any"`, the null hypothesis is that the true log-fold change in proportions is either `-lfc` or `lfc` with equal probability. A two-sided p-value is computed against this composite null.

- If `direction="up"`, the null hypothesis is that the true log-fold change is `lfc`, and a one-sided p-value is computed.

- If `direction="down"`, the null hypothesis is that the true log-fold change is `-lfc`, and a one-sided p-value is computed.

### Blocking on uninteresting factors

If `block` is specified, binomial tests are performed between groups of cells within each level of `block`. For each pair of groups, the p-values for each gene across all levels of `block` are combined using Stouffer's weighted Z-score method.

The weight for the p-value in a particular level of `block` is defined as N_x + N_y, where N_x and N_y are the number of cells in groups X and Y, respectively, for that level. This means that p-values from blocks with more cells will have a greater contribution to the combined p-value for each gene.

When combining across batches, one-sided p-values in the same direction are combined first. Then, if `direction="any"`, the two combined p-values from both directions are combined. This ensures that a gene only receives a low overall p-value if it changes in the same direction across batches.

When comparing two groups, blocking levels are ignored if no p-value was reported, e.g., if there were insufficient cells for a group in a particular level. If all levels are ignored in this manner, the entire comparison will only contain `NA` p-values and a warning will be emitted.

### Author(s)

Aaron Lun

### References

Whitlock MC (2005). Combining probability from independent tests: the weighted Z-method is superior to Fisher's approach. J. Evol. Biol. 18, 5:1368-73.

### See Also

`binom.test`, on which this function is based.

### Examples

    library(scater)
    sce <- mockSCE()
    sce <- logNormCounts(sce)

    # Any clustering method is okay.
    kout <- kmeans(t(logcounts(sce)), centers=2)

    # Vanilla application:
    out <- pairwiseBinom(logcounts(sce), groups=kout$cluster)
    out

    # Directional and with a minimum log-fold change:
    out <- pairwiseBinom(logcounts(sce), groups=kout$cluster,
        direction="up", lfc=1)
    out

[Package scran version 1.14.6 Index]