combineMarkers: Combine pairwise DE results into a marker list

In MarioniLab/scran: Methods for Single-Cell RNA-Seq Data Analysis

Combine pairwise DE results into a marker list

Description

Combine multiple pairwise differential expression comparisons between groups or clusters into a single ranked list of markers for each cluster.

Usage

combineMarkers(
  de.lists,
  pairs,
  pval.field = "p.value",
  effect.field = "logFC",
  pval.type = c("any", "some", "all"),
  min.prop = NULL,
  log.p.in = FALSE,
  log.p.out = log.p.in,
  output.field = NULL,
  full.stats = FALSE,
  sorted = TRUE,
  flatten = TRUE,
  BPPARAM = SerialParam()
)

Arguments

de.lists	A list-like object where each element is a data.frame or DataFrame. Each element should represent the results of a pairwise comparison between two groups/clusters, in which each row should contain the statistics for a single gene/feature. Rows should be named by the feature name in the same order for all elements.
pairs	A matrix, data.frame or DataFrame with two columns and number of rows equal to the length of de.lists. Each row should specify the pair of clusters being compared for the corresponding element of de.lists.
pval.field	A string specifying the column name of each element of de.lists that contains the p-value.
effect.field	A string specifying the column name of each element of de.lists that contains the effect size. If NULL, effect sizes are not reported in the output.
pval.type	A string specifying how p-values are to be combined across pairwise comparisons for a given group/cluster.
min.prop	Numeric scalar specifying the minimum proportion of significant comparisons per gene, Defaults to 0.5 when pval.type="some", otherwise defaults to zero.
log.p.in	A logical scalar indicating if the p-values in de.lists were log-transformed.
log.p.out	A logical scalar indicating if log-transformed p-values/FDRs should be returned.
output.field	A string specifying the prefix of the field names containing the effect sizes. Defaults to "stats" if full.stats=TRUE, otherwise it is set to effect.field.
full.stats	A logical scalar indicating whether all statistics in de.lists should be stored in the output for each pairwise comparison.
sorted	Logical scalar indicating whether each output DataFrame should be sorted by a statistic relevant to pval.type.
flatten	Logical scalar indicating whether the individual effect sizes should be flattened in the output DataFrame. If FALSE, effect sizes are reported as a nested matrix for easier programmatic use.
BPPARAM	A BiocParallelParam object indicating whether and how parallelization should be performed across genes.

Details

An obvious strategy to characterizing differences between clusters is to look for genes that are differentially expressed (DE) between them. However, this entails a number of comparisons between all pairs of clusters to comprehensively identify genes that define each cluster. For all pairwise comparisons involving a single cluster, we would like to consolidate the DE results into a single list of candidate marker genes. Doing so is the purpose of the combineMarkers function.

DE statistics from any testing regime can be supplied to this function - see the Examples for how this is done with t-tests from pairwiseTTests. The effect size field in the output will vary according to the type of input statistics, for example:

• logFC.Y from pairwiseTTests, containing log-fold changes in mean expression (usually in base 2).

• AUC.Y from pairwiseWilcox, containing the area under the curve, i.e., the concordance probability.

• logFC.Y from pairwiseBinom, containing log2-fold changes in the expressing proportion.

Value

A named List of DataFrames where each DataFrame contains the consolidated marker statistics for each gene (row) for the cluster of the same name. The DataFrame for cluster X contains the fields:

Top:

Integer, the minimum rank across all pairwise comparisons - see ?computeMinRank for details. This is only reported if pval.type="any".

p.value:

Numeric, the combined p-value across all comparisons if log.p.out=FALSE.

FDR:

Numeric, the BH-adjusted p-value for each gene if log.p.out=FALSE.

log.p.value:

Numeric, the (natural) log-transformed version p-value. Replaces the p.value field if log.p.out=TRUE.

log.FDR:

Numeric, the (natural) log-transformed adjusted p-value. Replaces the FDR field if log.p.out=TRUE.

summary.<OUTPUT>:

Numeric, named by replacing <OUTPUT> with output.field. This contains the summary effect size, obtained by combining effect sizes from all pairwise comparison into a single value. Only reported when effect.field is not NULL.

<OUTPUT>.Y:

Comparison-specific statistics, named by replacing <OUTPUT> with output.field. One of these fields is present for every other cluster Y in clusters and contains statistics for the comparison of X to Y. If full.stats=FALSE, each field is numeric and contains the effect size of the comparison of X over Y. Otherwise, each field is a nested DataFrame containing the full statistics for that comparison (i.e., the same asthe corresponding entry of de.lists). Only reported if flatten=FALSE and (for full.stats=FALSE) if effect.field is not NULL.

each.<OUTPUT>:

A nested DataFrame of comparison-specific statistics, named by replacing <OUTPUT> with output.field. If full.stats=FALSE, one column is present for every other cluster Y in clusters and contains the effect size of the comparison of X to Y. Otherwise, each column contains another nested DataFrame containing the full set of statistics for that comparison. Only reported if flatten=FALSE and (for full.stats=FALSE) if effect.field is not NULL.

Consolidating with DE against any other cluster

When pval.type="any", each DataFrame is sorted by the min-rank in the Top column. Taking all rows with Top values less than or equal to T yields a marker set containing the top T genes (ranked by significance) from each pairwise comparison. This guarantees the inclusion of genes that can distinguish between any two clusters. Also see ?computeMinRank for more details on the rationale behind this metric.

For each gene and cluster, the summary effect size is defined as the effect size from the pairwise comparison with the lowest p-value. The combined p-value is computed by applying Simes' method to all p-values. Neither of these values are directly used for ranking and are only reported for the sake of the user.

Consolidating with DE against all other clusters

If pval.type="all", the null hypothesis is that the gene is not DE in all contrasts. A combined p-value for each gene is computed using Berger's intersection union test (IUT). Ranking based on the IUT p-value will focus on genes that are DE in that cluster compared to all other clusters. This strategy is particularly effective when dealing with distinct clusters that have a unique expression profile. In such cases, it yields a highly focused marker set that concisely captures the differences between clusters.

However, it can be too stringent if the cluster's separation is driven by combinations of gene expression. For example, consider a situation involving four clusters expressing each combination of two marker genes A and B. With pval.type="all", neither A nor B would be detected as markers as it is not uniquely defined in any one cluster. This is especially detrimental with overclustering where an otherwise acceptable marker is discarded if it is not DE between two adjacent clusters.

For each gene and cluster, the summary effect size is defined as the effect size from the pairwise comparison with the largest p-value. This reflects the fact that, with this approach, a gene is only as significant as its weakest DE. Again, this value is not directly used for ranking and are only reported for the sake of the user.

Consolidating with DE against some other clusters

The pval.type="some" setting serves as a compromise between "all" and "any". A combined p-value is calculated by taking the middlemost value of the Holm-corrected p-values for each gene. (By default, this the median for odd numbers of contrasts and one-after-the-median for even numbers, but the exact proportion can be changed by setting min.prop - see ?combineParallelPValues.) Here, the null hypothesis is that the gene is not DE in at least half of the contrasts.

Genes are then ranked by the combined p-value. The aim is to provide a more focused marker set without being overly stringent, though obviously it loses the theoretical guarantees of the more extreme settings. For example, there is no guarantee that the top set contains genes that can distinguish a cluster from any other cluster, which would have been possible with pval.type="any".

For each gene and cluster, the summary effect size is defined as the effect size from the pairwise comparison with the min.prop-smallest p-value. This mirrors the p-value calculation but, again, is reported only for the benefit of the user.

Consolidating against some other clusters, rank-style

A slightly different flavor of the “some cluster” approach is achieved by setting method="any" with min.prop set to some positive value in (0, 1). A gene will only be high-ranked if it is among the top-ranked genes in at least min.prop of the pairwise comparisons. For example, if min.prop=0.3, any gene with a value of Top less than or equal to 5 will be in the top 5 DEGs of at least 30% of the comparisons.

This method increases the stringency of the "any" setting in a safer manner than pval.type="some". Specifically, we avoid comparing p-values across pairwise comparisons, which can be problematic if there are power differences across comparisons, e.g., due to differences in the number of cells across the other clusters.

Note that the value of min.prop does not affect the combined p-value and summary effect size calculations for pval.type="any".

Correcting for multiple testing

The BH method is then applied on the consolidated p-values across all genes to obtain the FDR field. The reported FDRs are intended only as a rough measure of significance. Properly correcting for multiple testing is not generally possible when clusters is determined from the same x used for DE testing.

If log.p=TRUE, log-transformed p-values and FDRs will be reported. This may be useful in over-powered studies with many cells, where directly reporting the raw p-values would result in many zeroes due to the limits of machine precision.

Ordering of the output

• Within each DataFrame, if sorted=TRUE, genes are ranked by the Top column if available and the p.value (or log.p.value) if not. Otherwise, the input order of the genes is preserved.

• For the DataFrame corresponding to cluster X, the <OUTPUT>.Y columns are sorted according to the order of cluster IDs in pairs[,2] for all rows where pairs[,1] is X.

• In the output List, the DataFrames themselves are sorted according to the order of cluster IDs in pairs[,1]. Note that DataFrames are only created for clusters present in pairs[,1]. Clusters unique to pairs[,2] will only be present within a DataFrame as Y.

Author(s)

Aaron Lun

References

Simes RJ (1986). An improved Bonferroni procedure for multiple tests of significance. Biometrika 73:751-754.

Berger RL and Hsu JC (1996). Bioequivalence trials, intersection-union tests and equivalence confidence sets. Statist. Sci. 11, 283-319.

See Also

pairwiseTTests and pairwiseWilcox, for functions that can generate de.lists and pairs.

findMarkers, which automatically performs combineMarkers on the t-test or Wilcoxon test results.

Examples

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

# Any clustering method is okay.
kout <- kmeans(t(logcounts(sce)), centers=3)
clusters <- paste0("Cluster", kout$cluster)

out <- pairwiseTTests(logcounts(sce), groups=clusters)
comb <- combineMarkers(out$statistics, out$pairs)
comb[["Cluster1"]]

out <- pairwiseWilcox(logcounts(sce), groups=clusters)
comb <- combineMarkers(out$statistics, out$pairs, effect.field="AUC")
comb[["Cluster2"]]

out <- pairwiseBinom(logcounts(sce), groups=clusters)
comb <- combineMarkers(out$statistics, out$pairs)
comb[["Cluster3"]]

MarioniLab/scran documentation built on Sept. 7, 2024, 6:25 a.m.