% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/combineMarkers.R
\name{combineMarkers}
\alias{combineMarkers}
\title{Combine pairwise DE results into a marker list}
\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{
\item{de.lists}{A list-like object where each element is a data.frame or \linkS4class{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.}

\item{pairs}{A matrix, data.frame or \linkS4class{DataFrame} with two columns and number of rows equal to the length of \code{de.lists}.
Each row should specify the pair of clusters being compared for the corresponding element of \code{de.lists}.}

\item{pval.field}{A string specifying the column name of each element of \code{de.lists} that contains the p-value.}

\item{effect.field}{A string specifying the column name of each element of \code{de.lists} that contains the effect size.
If \code{NULL}, effect sizes are not reported in the output.}

\item{pval.type}{A string specifying how p-values are to be combined across pairwise comparisons for a given group/cluster.}

\item{min.prop}{Numeric scalar specifying the minimum proportion of significant comparisons per gene,
Defaults to 0.5 when \code{pval.type="some"}, otherwise defaults to zero.}

\item{log.p.in}{A logical scalar indicating if the p-values in \code{de.lists} were log-transformed.}

\item{log.p.out}{A logical scalar indicating if log-transformed p-values/FDRs should be returned.}

\item{output.field}{A string specifying the prefix of the field names containing the effect sizes.
Defaults to \code{"stats"} if \code{full.stats=TRUE}, otherwise it is set to \code{effect.field}.}

\item{full.stats}{A logical scalar indicating whether all statistics in \code{de.lists} should be stored in the output for each pairwise comparison.}

\item{sorted}{Logical scalar indicating whether each output DataFrame should be sorted by a statistic relevant to \code{pval.type}.}

\item{flatten}{Logical scalar indicating whether the individual effect sizes should be flattened in the output DataFrame.
If \code{FALSE}, effect sizes are reported as a nested matrix for easier programmatic use.}

\item{BPPARAM}{A \linkS4class{BiocParallelParam} object indicating whether and how parallelization should be performed across genes.}
}
\value{
A named \linkS4class{List} of \linkS4class{DataFrame}s where each DataFrame contains the consolidated marker statistics for each gene (row) for the cluster of the same name.
The DataFrame for cluster \eqn{X} contains the fields:
\describe{
\item{\code{Top}:}{Integer, the minimum rank across all pairwise comparisons.
This is only reported if \code{pval.type="any"}.}
\item{\code{p.value}:}{Numeric, the combined p-value across all comparisons if \code{log.p.out=FALSE}.}
\item{\code{FDR}:}{Numeric, the BH-adjusted p-value for each gene if \code{log.p.out=FALSE}.}
\item{\code{log.p.value}:}{Numeric, the (natural) log-transformed version p-value.
Replaces the \code{p.value} field if \code{log.p.out=TRUE}.}
\item{\code{log.FDR}:}{Numeric, the (natural) log-transformed adjusted p-value.
Replaces the \code{FDR} field if \code{log.p.out=TRUE}.}
\item{\code{summary.<OUTPUT>}:}{Numeric, named by replacing \code{<OUTPUT>} with \code{output.field}.
This contains the summary effect size, obtained by combining effect sizes from all pairwise comparison into a single value.
Only reported when \code{effect.field} is not \code{NULL}.}
\item{\code{<OUTPUT>.Y}:}{Comparison-specific statistics, named by replacing \code{<OUTPUT>} with \code{output.field}.
One of these fields is present for every other cluster \eqn{Y} in \code{clusters} and contains statistics for the comparison of \eqn{X} to \eqn{Y}.
If \code{full.stats=FALSE}, each field is numeric and contains the effect size of the comparison of \eqn{X} over \eqn{Y}.
Otherwise, each field is a nested DataFrame containing the full statistics for that comparison (i.e., the same asthe corresponding entry of \code{de.lists}).
Only reported if \code{flatten=FALSE} and (for \code{full.stats=FALSE}) if \code{effect.field} is not \code{NULL}.
}
\item{\code{each.<OUTPUT>}:}{A nested DataFrame of comparison-specific statistics, named by replacing \code{<OUTPUT>} with \code{output.field}.
If \code{full.stats=FALSE}, one column is present for every other cluster \eqn{Y} in \code{clusters} and contains the effect size of the comparison of \eqn{X} to \eqn{Y}.
Otherwise, each column contains another nested DataFrame containing the full set of statistics for that comparison.
Only reported if \code{flatten=FALSE} and (for \code{full.stats=FALSE}) if \code{effect.field} is not \code{NULL}.
}
}
}
\description{
Combine multiple pairwise differential expression comparisons between groups or clusters into a single ranked list of markers for each cluster.
}
\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 \code{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 \code{\link{pairwiseTTests}}.
The effect size field in the output will vary according to the type of input statistics, for example:
\itemize{
\item \code{logFC.Y} from \code{\link{pairwiseTTests}}, containing log-fold changes in mean expression (usually in base 2).
\item \code{AUC.Y} from \code{\link{pairwiseWilcox}}, containing the area under the curve, i.e., the concordance probability. 
\item \code{logFC.Y} from \code{\link{pairwiseBinom}}, containing log2-fold changes in the expressing proportion.
}
}
\section{Consolidating with DE against any other cluster}{

By default, each DataFrame is sorted by the \code{Top} value when \code{pval.type="any"}.
Taking all rows with \code{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.

To demonstrate, let us define a marker set with an T of 1 for a given cluster.
The set of genes with \code{Top <= 1} will contain the top gene from each pairwise comparison to every other cluster.
If T is instead, say, 5, the set will consist of the \emph{union} of the top 5 genes from each pairwise comparison.
Obviously, multiple genes can have the same \code{Top} as different genes may have the same rank across different pairwise comparisons.
Conversely, the marker set may be smaller than the product of \code{Top} and the number of other clusters, as the same gene may be shared across different comparisons.

This approach does not explicitly favour genes that are uniquely expressed in a cluster.
Rather, it focuses on combinations of genes that - together - drive separation of a cluster from the others.
This is more general and robust but tends to yield a less focused marker set compared to the other \code{pval.type} settings.

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.
}

\section{Consolidating with DE against all other clusters}{

If \code{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 \emph{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 \code{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 \emph{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.
}

\section{Consolidating with DE against some other clusters}{

The \code{pval.type="some"} setting serves as a compromise between \code{"all"} and \code{"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 \code{min.prop} - see \code{?\link{combinePValues}}.)
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 \code{pval.type="any"}.

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

\section{Consolidating against some other clusters, rank-style}{

A slightly different flavor of the \dQuote{some cluster} approach is achieved by setting \code{method="any"} with \code{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 \code{min.prop} of the pairwise comparisons.
For example, if \code{min.prop=0.3}, any gene with a value of \code{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 \code{"any"} setting in a safer manner than \code{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 \code{min.prop} does not affect the combined p-value and summary effect size calculations for \code{pval.type="any"}.
}

\section{Correcting for multiple testing}{

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

If \code{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.
}

\section{Ordering of the output}{

\itemize{
\item Within each DataFrame, if \code{sorted=TRUE}, genes are ranked by the \code{Top} column if available and the \code{p.value} (or \code{log.p.value}) if not.
Otherwise, the input order of the genes is preserved.
\item For the DataFrame corresponding to cluster \eqn{X}, the \code{<OUTPUT>.Y} columns are sorted according to the order of cluster IDs in \code{pairs[,2]} for all rows where \code{pairs[,1]} is \eqn{X}.
\item In the output List, the DataFrames themselves are sorted according to the order of cluster IDs in \code{pairs[,1]}.
Note that DataFrames are only created for clusters present in \code{pairs[,1]}.
Clusters unique to \code{pairs[,2]} will only be present within a DataFrame as \eqn{Y}.
}
}

\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"]]

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

Berger RL and Hsu JC (1996). 
Bioequivalence trials, intersection-union tests and equivalence confidence sets.
\emph{Statist. Sci.} 11, 283-319.
}
\seealso{
\code{\link{pairwiseTTests}} and \code{\link{pairwiseWilcox}}, for functions that can generate \code{de.lists} and \code{pairs}.

\code{\link{findMarkers}}, which automatically performs \code{combineMarkers} on the t-test or Wilcoxon test results.
}
\author{
Aaron Lun
}