---
output:
  html_document
bibliography: ref.bib
---

# Marker gene detection {#marker-detection}

<script>
document.addEventListener("click", function (event) {
    if (event.target.classList.contains("rebook-collapse")) {
        event.target.classList.toggle("active");
        var content = event.target.nextElementSibling;
        if (content.style.display === "block") {
            content.style.display = "none";
        } else {
            content.style.display = "block";
        }
    }
})
</script>

<style>
.rebook-collapse {
  background-color: #eee;
  color: #444;
  cursor: pointer;
  padding: 18px;
  width: 100%;
  border: none;
  text-align: left;
  outline: none;
  font-size: 15px;
}

.rebook-content {
  padding: 0 18px;
  display: none;
  overflow: hidden;
  background-color: #f1f1f1;
}
</style>

## Motivation

To interpret our clustering results from Chapter \@ref(clustering), we identify the genes that drive separation between clusters.
These marker genes allow us to assign biological meaning to each cluster based on their functional annotation.
In the simplest case, we have _a priori_ knowledge of the marker genes associated with particular cell types, allowing us to treat the clustering as a proxy for cell type identity.
The same principle can be applied to discover more subtle differences between clusters (e.g., changes in activation or differentiation state) based on the behavior of genes in the affected pathways.

The most straightforward approach to marker gene detection involves testing for differential expression between clusters.
If a gene is strongly DE between clusters, it is likely to have driven the separation of cells in the clustering algorithm.
Several methods are available to quantify the differences in expression profiles between clusters and obtain a single ranking of genes for each cluster.
We will demonstrate some of these choices in this chapter using the 10X PBMC dataset:

<button class="rebook-collapse">View set-up code ([Workflow Chapter 3](http://bioconductor.org/books/3.20/OSCA.workflows/unfiltered-human-pbmcs-10x-genomics.html#unfiltered-human-pbmcs-10x-genomics))</button>
<div class="rebook-content">

```r
#--- loading ---#
library(DropletTestFiles)
raw.path <- getTestFile("tenx-2.1.0-pbmc4k/1.0.0/raw.tar.gz")
out.path <- file.path(tempdir(), "pbmc4k")
untar(raw.path, exdir=out.path)

library(DropletUtils)
fname <- file.path(out.path, "raw_gene_bc_matrices/GRCh38")
sce.pbmc <- read10xCounts(fname, col.names=TRUE)

#--- gene-annotation ---#
library(scater)
rownames(sce.pbmc) <- uniquifyFeatureNames(
    rowData(sce.pbmc)$ID, rowData(sce.pbmc)$Symbol)

library(EnsDb.Hsapiens.v86)
location <- mapIds(EnsDb.Hsapiens.v86, keys=rowData(sce.pbmc)$ID, 
    column="SEQNAME", keytype="GENEID")

#--- cell-detection ---#
set.seed(100)
e.out <- emptyDrops(counts(sce.pbmc))
sce.pbmc <- sce.pbmc[,which(e.out$FDR <= 0.001)]

#--- quality-control ---#
stats <- perCellQCMetrics(sce.pbmc, subsets=list(Mito=which(location=="MT")))
high.mito <- isOutlier(stats$subsets_Mito_percent, type="higher")
sce.pbmc <- sce.pbmc[,!high.mito]

#--- normalization ---#
library(scran)
set.seed(1000)
clusters <- quickCluster(sce.pbmc)
sce.pbmc <- computeSumFactors(sce.pbmc, cluster=clusters)
sce.pbmc <- logNormCounts(sce.pbmc)

#--- variance-modelling ---#
set.seed(1001)
dec.pbmc <- modelGeneVarByPoisson(sce.pbmc)
top.pbmc <- getTopHVGs(dec.pbmc, prop=0.1)

#--- dimensionality-reduction ---#
set.seed(10000)
sce.pbmc <- denoisePCA(sce.pbmc, subset.row=top.pbmc, technical=dec.pbmc)

set.seed(100000)
sce.pbmc <- runTSNE(sce.pbmc, dimred="PCA")

set.seed(1000000)
sce.pbmc <- runUMAP(sce.pbmc, dimred="PCA")

#--- clustering ---#
g <- buildSNNGraph(sce.pbmc, k=10, use.dimred = 'PCA')
clust <- igraph::cluster_walktrap(g)$membership
colLabels(sce.pbmc) <- factor(clust)
```

</div>


``` r
sce.pbmc
```

```
## class: SingleCellExperiment 
## dim: 33694 3985 
## metadata(1): Samples
## assays(2): counts logcounts
## rownames(33694): RP11-34P13.3 FAM138A ... AC213203.1 FAM231B
## rowData names(2): ID Symbol
## colnames(3985): AAACCTGAGAAGGCCT-1 AAACCTGAGACAGACC-1 ...
##   TTTGTCAGTTAAGACA-1 TTTGTCATCCCAAGAT-1
## colData names(4): Sample Barcode sizeFactor label
## reducedDimNames(3): PCA TSNE UMAP
## mainExpName: NULL
## altExpNames(0):
```

## Scoring markers by pairwise comparisons



Our general strategy is to compare each pair of clusters and compute scores quantifying the differences in the expression distributions between clusters.
The scores for all pairwise comparisons involving a particular cluster are then consolidated into a single `DataFrame` for that cluster.
The `scoreMarkers()` function from *[scran](https://bioconductor.org/packages/3.20/scran)* returns a list of `DataFrame`s where each `DataFrame` corresponds to a cluster and each row of the `DataFrame` corresponds to a gene.
In the `DataFrame` for cluster $X$, the columns contain the `self.average`, the mean log-expression in $X$;
`other.average`, the grand mean across all other clusters;
`self.detected`, the proportion of cells with detected expression in $X$;
`other.detected`, the mean detected proportion across all other clusters;
and finally, a variety of effect size summaries generated from all pairwise comparisons involving $X$.


``` r
library(scran)
marker.info <- scoreMarkers(sce.pbmc, colLabels(sce.pbmc))
marker.info
```

```
## List of length 15
## names(15): 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
```

``` r
colnames(marker.info[["1"]]) # statistics for cluster 1.
```

```
##  [1] "self.average"          "other.average"         "self.detected"        
##  [4] "other.detected"        "mean.logFC.cohen"      "min.logFC.cohen"      
##  [7] "median.logFC.cohen"    "max.logFC.cohen"       "rank.logFC.cohen"     
## [10] "mean.AUC"              "min.AUC"               "median.AUC"           
## [13] "max.AUC"               "rank.AUC"              "mean.logFC.detected"  
## [16] "min.logFC.detected"    "median.logFC.detected" "max.logFC.detected"   
## [19] "rank.logFC.detected"
```

For each cluster, we can then rank candidate markers based on one of these effect size summaries.
We demonstrate below with the mean AUC for cluster 1, which probably contains NK cells based on the top genes in Figure \@ref(fig:pbmc-scored-markers-1) (and no _CD3E_ expression).
The next section will go into more detail on the differences between the various columns.


``` r
chosen <- marker.info[["1"]]
ordered <- chosen[order(chosen$mean.AUC, decreasing=TRUE),]
head(ordered[,1:4]) # showing basic stats only, for brevity.
```

```
## DataFrame with 6 rows and 4 columns
##       self.average other.average self.detected other.detected
##          <numeric>     <numeric>     <numeric>      <numeric>
## NKG7       4.39503      0.796798      0.985366      0.3885337
## GNLY       4.40275      0.365112      0.956098      0.2038427
## CTSW       2.60281      0.496249      0.960976      0.3072619
## HOPX       1.99060      0.172733      0.907317      0.1238411
## PRF1       2.22297      0.178523      0.887805      0.1344186
## KLRF1      1.60598      0.041239      0.858537      0.0374419
```

``` r
library(scater)
plotExpression(sce.pbmc, features=head(rownames(ordered)), 
    x="label", colour_by="label")
```

<div class="figure">
<img src="marker-detection_files/figure-html/pbmc-scored-markers-1-1.png" alt="Distribution of expression values across clusters for the top potential marker genes (as determined by the mean AUC) for cluster 1 in the PBMC dataset." width="768" />
<p class="caption">(\#fig:pbmc-scored-markers-1)Distribution of expression values across clusters for the top potential marker genes (as determined by the mean AUC) for cluster 1 in the PBMC dataset.</p>
</div>



We deliberately use pairwise comparisons rather than comparing each cluster to the average of all other cells.
The latter approach is sensitive to the population composition, which introduces an element of unpredictability to the marker sets due to variation in cell type abundances.
(In the worst case, the presence of one subpopulation containing a majority of the cells will drive the selection of top markers for every other cluster, pushing out useful genes that can distinguish between the smaller subpopulations.)
Moreover, pairwise comparisons naturally provide more information to interpret of the utility of a marker, e.g., by providing log-fold changes to indicate which clusters are distinguished by each gene (Section \@ref(obtaining-the-full-effects)).

Previous editions of this chapter used $p$-values from the tests corresponding to each effect size, e.g., Welch's $t$-test, the Wilcoxon ranked sum test.
While this is fine for ranking genes, the $p$-values themselves are statistically flawed and are of little use for inference -
see [Advanced Section 6.4](http://bioconductor.org/books/3.20/OSCA.advanced/marker-detection-redux.html#p-value-invalidity) for more details.
The `scoreMarkers()` function simplifies the marker detection procedure by omitting the $p$-values altogether, instead focusing on the underlying effect sizes.

## Effect sizes for pairwise comparisons

In the context of marker detection, the area under the curve (AUC) quantifies our ability to distinguish between two distributions in a pairwise comparison.
The AUC represents the probability that a randomly chosen observation from our cluster of interest is greater than a randomly chosen observation from the other cluster.
A value of 1 corresponds to upregulation, where all values of our cluster of interest are greater than any value from the other cluster;
a value of 0.5 means that there is no net difference in the location of the distributions;
and a value of 0 corresponds to downregulation.
The AUC is closely related to the $U$ statistic in the Wilcoxon ranked sum test (a.k.a., Mann-Whitney U-test).


``` r
auc.only <- chosen[,grepl("AUC", colnames(chosen))]
auc.only[order(auc.only$mean.AUC,decreasing=TRUE),]
```

```
## DataFrame with 33694 rows and 5 columns
##         mean.AUC    min.AUC median.AUC   max.AUC  rank.AUC
##        <numeric>  <numeric>  <numeric> <numeric> <integer>
## NKG7    0.960888   0.803539   0.988883  0.991803         1
## GNLY    0.957445   0.876231   0.972312  0.975021         2
## CTSW    0.928138   0.709581   0.958830  0.978412         2
## HOPX    0.926121   0.777746   0.948263  0.953659         3
## PRF1    0.923078   0.831292   0.938935  0.943902         5
## ...          ...        ...        ...       ...       ...
## RPS13   0.180687 0.00151257  0.0829163  0.695732       426
## RPL26   0.179083 0.00741161  0.0657345  0.813720       104
## RPL18A  0.175037 0.01090122  0.0844709  0.749735       264
## RPL39   0.160374 0.00086973  0.0649374  0.774848       185
## FTH1    0.157883 0.00134567  0.1055687  0.645732       766
```

Cohen's $d$ is a standardized log-fold change where the difference in the mean log-expression between groups is scaled by the average standard deviation across groups.
In other words, it is the number of standard deviations that separate the means of the two groups.
The interpretation is similar to the log-fold change; positive values indicate that the gene is upregulated in our cluster of interest, 
negative values indicate downregulation and values close to zero indicate that there is little difference.
Cohen's $d$ is roughly analogous to the $t$-statistic in various two-sample $t$-tests.


``` r
cohen.only <- chosen[,grepl("logFC.cohen", colnames(chosen))]
cohen.only[order(cohen.only$mean.logFC.cohen,decreasing=TRUE),]
```

```
## DataFrame with 33694 rows and 5 columns
##         mean.logFC.cohen min.logFC.cohen median.logFC.cohen max.logFC.cohen
##                <numeric>       <numeric>          <numeric>       <numeric>
## NKG7             4.75176        1.009867            5.74481         6.30987
## GNLY             3.67887        1.781382            4.02699         4.37410
## CTSW             2.88706        0.700977            3.18157         3.96973
## GZMA             2.66326        0.393043            3.18665         3.44040
## HOPX             2.64870        1.095214            2.88657         3.06690
## ...                  ...             ...                ...             ...
## FTH1            -2.14556        -5.15590          -1.949032       0.2995322
## FTL             -2.17697        -5.82525          -1.062434       0.2569453
## HLA-DRA         -2.21211        -7.13493          -1.291256       0.0234466
## CST3            -2.28988        -7.92584          -0.580258       0.0254270
## LYZ             -2.46892        -7.77904          -0.326158      -0.1797338
##         rank.logFC.cohen
##                <integer>
## NKG7                   1
## GNLY                   1
## CTSW                   3
## GZMA                   2
## HOPX                   4
## ...                  ...
## FTH1                4362
## FTL                 5319
## HLA-DRA             6161
## CST3                5867
## LYZ                30966
```

Finally, we also compute the log-fold change in the proportion of cells with detected expression between clusters.
This ignores any information about the magnitude of expression, only considering whether any expression is detected at all.
Again, positive values indicate that a greater proportion of cells express the gene in our cluster of interest compared to the other cluster.
Note that a pseudo-count is added to avoid undefined log-fold changes when no cells express the gene in either group.


``` r
detect.only <- chosen[,grepl("logFC.detected", colnames(chosen))]
detect.only[order(detect.only$mean.logFC.detected,decreasing=TRUE),]
```

```
## DataFrame with 33694 rows and 5 columns
##        mean.logFC.detected min.logFC.detected median.logFC.detected
##                  <numeric>          <numeric>             <numeric>
## KLRF1              4.68657           2.447509               5.19667
## PRSS23             4.36839           2.272651               4.40053
## XCL2               4.28046           1.249518               4.81949
## XCL1               4.27138           1.343521               4.70897
## SH2D1B             4.11136           0.853048               4.57000
## ...                    ...                ...                   ...
## MARCKS            -2.87363           -6.61517              -1.58496
## SLC7A7            -2.89504           -6.70323              -2.52746
## RAB31             -3.00009           -6.50367              -2.08496
## RAB32             -3.11607           -6.61471              -3.02746
## NCF2              -3.50919           -6.90913              -3.06253
##        max.logFC.detected rank.logFC.detected
##                 <numeric>           <integer>
## KLRF1             6.37216                   1
## PRSS23            5.96312                   2
## XCL2              6.18058                   1
## XCL1              6.16993                   1
## SH2D1B            5.73588                   1
## ...                   ...                 ...
## MARCKS        0.00000e+00               11805
## SLC7A7        1.60171e-16                7775
## RAB31        -1.00000e+00               31423
## RAB32         1.60171e-16                7775
## NCF2          0.00000e+00               11805
```

The AUC or Cohen's $d$ is usually the best choice for general purpose marker detection, as they are effective regardless of the magnitude of the expression values.
The log-fold change in the detected proportion is specifically useful for identifying binary changes in expression.
See [Advanced Section 6.2](http://bioconductor.org/books/3.20/OSCA.advanced/marker-detection-redux.html#properties-of-each-effect-size) for more information about the practical differences between the effect sizes.

## Summarizing pairwise effects 

In a dataset with $N$ clusters, each cluster is associated with $N-1$ values for each type of effect size described in the previous section.
To simplify interpretation, we summarize the effects for each cluster into some key statistics such as the mean and median.
Each summary statistic has a different interpretation when used for ranking:

- The most obvious summary statistic is the mean.
For cluster $X$, a large mean effect size (>0 for the log-fold changes, >0.5 for the AUCs) indicates that the gene is upregulated in $X$ compared to the average of the other groups.
- Another summary statistic is the median, where a large value indicates that the gene is upregulated in $X$ compared to most (>50%) other clusters.
The median provides greater robustness to outliers than the mean, which may or may not be desirable.
On one hand, the median avoids an inflated effect size if only a minority of comparisons have large effects;
on the other hand, it will also overstate the effect size by ignoring a minority of comparisons that have opposing effects.
- The minimum value (`min.*`) is the most stringent summary for identifying upregulated genes, as a large value indicates that the gene is upregulated in $X$ compared to _all_ other clusters.
Conversely, if the minimum is small (<0 for the log-fold changes, <0.5 for the AUCs), we can conclude that the gene is downregulated in $X$ compared to at least one other cluster.
- The maximum value (`max.*`) is the least stringent summary for identifying upregulated genes, as a large value can be obtained if there is strong upregulation in $X$ compared to _any_ other cluster.
Conversely, if the maximum is small, we can conclude that the gene is downregulated in $X$ compared to all other clusters.
- The minimum rank, a.k.a., "min-rank" (`rank.*`) is the smallest rank of each gene across all pairwise comparisons.
Specifically, genes are ranked _within_ each pairwise comparison based on decreasing effect size, and then the smallest rank _across_ all comparisons is reported for each gene.
If a gene has a small min-rank, we can conclude that it is one of the top upregulated genes in at least one comparison of $X$ to another cluster.



Each of these summaries is computed for each effect size, for each gene, and for each cluster.
Our next step is to choose one of these summary statistics for one of the effect sizes and to use it to rank the rows of the `DataFrame`.
The choice of summary determines the stringency of the marker selection strategy, i.e., how many other clusters must we differ from?
For identifying upregulated genes, ranking by the minimum is the most stringent and the maximum is the least stringent;
the mean and median fall somewhere in between and are reasonable defaults for most applications.
The example below uses the median Cohen's $d$ to obtain a ranking of upregulated markers for cluster 12 (Figure \@ref(fig:pbmc-scored-markers-again)), which probably contains monocytes.


``` r
chosen <- marker.info[["12"]] # using another cluster, for some variety.
ordered <- chosen[order(chosen$median.logFC.cohen,decreasing=TRUE),]
head(ordered[,1:4]) # showing basic stats only, for brevity.
```

```
## DataFrame with 6 rows and 4 columns
##               self.average other.average self.detected other.detected
##                  <numeric>     <numeric>     <numeric>      <numeric>
## LYZ                5.87025      1.879720      1.000000       0.631549
## S100A9             5.95360      1.767483      1.000000       0.685216
## S100A8             6.25358      1.504959      0.983607       0.611945
## FTL                6.33877      3.939846      1.000000       0.958558
## RP11-1143G9.4      3.81136      0.622745      0.983607       0.304197
## CTSS               3.79308      1.370146      0.950820       0.605402
```

``` r
plotExpression(sce.pbmc, features=head(rownames(ordered)), 
    x="label", colour_by="label")
```

<div class="figure">
<img src="marker-detection_files/figure-html/pbmc-scored-markers-again-1.png" alt="Distribution of expression values across clusters for the top potential marker genes (as determined by the median Cohen's $d$) for cluster 12 in the PBMC dataset." width="768" />
<p class="caption">(\#fig:pbmc-scored-markers-again)Distribution of expression values across clusters for the top potential marker genes (as determined by the median Cohen's $d$) for cluster 12 in the PBMC dataset.</p>
</div>



On some occasions, ranking by the minimum can be highly effective as it yields a concise set of highly cluster-specific markers. 
However, any gene that is expressed at the same level in two or more clusters will simply not be detected. 
This is likely to discard many interesting genes, especially if the clusters are finely resolved with weak separation.
To give a concrete example, consider a mixed population of CD4^+^-only, CD8^+^-only, double-positive and double-negative T cells.
Neither _Cd4_ or _Cd8_ would be detected as subpopulation-specific markers because each gene is expressed in two subpopulations such that the minimum effect would be small.
In practice, the minimum and maximum are most helpful for diagnosing discrepancies between the mean and median, rather than being used directly for ranking.

Ranking genes by the min-rank is similiar in stringency to ranking by the maximum effect size, in that both will respond to strong DE in a single comparison.
However, the min-rank is more useful as it ensures that a single comparison to another cluster with consistently large effects does not dominate the ranking.
If we select all genes with min-ranks less than or equal to $T$, the resulting set is the union of the top $T$ genes from all pairwise comparisons.
This guarantees that our set contains at least $T$ genes that can distinguish our cluster of interest from any other cluster, 
which permits a comprehensive determination of a cluster's identity.
We demonstrate below for cluster 12, taking the top $T=5$ genes with the largest Cohen's $d$ from each comparison to display in Figure \@ref(fig:pbmc-grouped-ranked-heat).


``` r
ordered <- chosen[order(chosen$rank.logFC.cohen),]
top.ranked <- ordered[ordered$rank.logFC.cohen <= 5,]
rownames(top.ranked)
```

```
##  [1] "S100A9"        "S100A8"        "S100A4"        "DUSP1"        
##  [5] "LYZ"           "FTL"           "CTSS"          "RP11-1143G9.4"
##  [9] "TYROBP"        "S100A6"        "FOS"           "S100A12"      
## [13] "H3F3B"
```

``` r
plotGroupedHeatmap(sce.pbmc, features=rownames(top.ranked), group="label", 
    center=TRUE, zlim=c(-3, 3))
```

<div class="figure">
<img src="marker-detection_files/figure-html/pbmc-grouped-ranked-heat-1.png" alt="Heatmap of the centered average log-expression values for the top potential marker genes for cluster 12 in the PBMC dataset. The set of markers was selected as those genes with Cohen's $d$-derived min-ranks less than or equal to 5." width="672" />
<p class="caption">(\#fig:pbmc-grouped-ranked-heat)Heatmap of the centered average log-expression values for the top potential marker genes for cluster 12 in the PBMC dataset. The set of markers was selected as those genes with Cohen's $d$-derived min-ranks less than or equal to 5.</p>
</div>

Our discussion above has focused mainly on potential markers that are upregulated in our cluster of interest, as these are the easiest to interpret and experimentally validate.
However, it also means that any cluster defined by downregulation of a marker will not contain that gene among the top features. 
This is occasionally relevant for subtypes or other states that are defined by low expression of particular genes. 
In such cases, focusing on upregulation may yield a disappointing set of markers,
and it may be worth examining some of the lowest-ranked genes to see if there is any consistent downregulation compared to other clusters.


``` r
# Omitting the decreasing=TRUE to focus on negative effects.
ordered <- chosen[order(chosen$median.logFC.cohen),1:4]
head(ordered)
```

```
## DataFrame with 6 rows and 4 columns
##        self.average other.average self.detected other.detected
##           <numeric>     <numeric>     <numeric>      <numeric>
## RPSA       0.830441       2.73850      0.377049       0.867700
## RPL23A     2.514675       3.86930      0.901639       0.905495
## EEF1B2     0.480169       2.01308      0.213115       0.825588
## RPL37A     1.942131       3.35917      0.721311       0.880085
## RPLP0      1.246229       2.62299      0.540984       0.865697
## RPS12      2.815977       4.33851      0.852459       0.919220
```

## Obtaining the full effects

For more complex questions, we may need to interrogate effect sizes from specific comparisons of interest.
To do so, we set `full.stats=TRUE` to obtain the effect sizes for all pairwise comparisons involving a particular cluster.
This is returned in the form of a nested `DataFrame` for each effect size type - 
in the example below, `full.AUC` contains the AUCs for the comparisons between cluster 12 and every other cluster.


``` r
marker.info <- scoreMarkers(sce.pbmc, colLabels(sce.pbmc), full.stats=TRUE)
chosen <- marker.info[["12"]]
chosen$full.AUC
```

```
## DataFrame with 33694 rows and 14 columns
##                      1         2         3         4         5         6
##              <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
## RP11-34P13.3       0.5  0.500000  0.500000       0.5  0.500000   0.50000
## FAM138A            0.5  0.500000  0.500000       0.5  0.500000   0.50000
## OR4F5              0.5  0.500000  0.500000       0.5  0.500000   0.50000
## RP11-34P13.7       0.5  0.497264  0.498379       0.5  0.499076   0.49858
## RP11-34P13.8       0.5  0.498632  0.500000       0.5  0.500000   0.50000
## ...                ...       ...       ...       ...       ...       ...
## AC233755.2    0.500000   0.50000  0.500000       0.5  0.500000  0.500000
## AC233755.1    0.500000   0.50000  0.500000       0.5  0.500000  0.500000
## AC240274.1    0.492683   0.49658  0.495948       0.5  0.494455  0.495739
## AC213203.1    0.500000   0.50000  0.500000       0.5  0.500000  0.500000
## FAM231B       0.500000   0.50000  0.500000       0.5  0.500000  0.500000
##                      7         8         9        10        11        13
##              <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
## RP11-34P13.3       0.5       0.5       0.5  0.500000       0.5       0.5
## FAM138A            0.5       0.5       0.5  0.500000       0.5       0.5
## OR4F5              0.5       0.5       0.5  0.500000       0.5       0.5
## RP11-34P13.7       0.5       0.5       0.5  0.489362       0.5       0.5
## RP11-34P13.8       0.5       0.5       0.5  0.500000       0.5       0.5
## ...                ...       ...       ...       ...       ...       ...
## AC233755.2       0.500       0.5  0.500000  0.500000  0.500000       0.5
## AC233755.1       0.500       0.5  0.500000  0.500000  0.500000       0.5
## AC240274.1       0.496       0.5  0.493895  0.489362  0.496732       0.5
## AC213203.1       0.500       0.5  0.500000  0.500000  0.500000       0.5
## FAM231B          0.500       0.5  0.500000  0.500000  0.500000       0.5
##                     14        15
##              <numeric> <numeric>
## RP11-34P13.3       0.5       0.5
## FAM138A            0.5       0.5
## OR4F5              0.5       0.5
## RP11-34P13.7       0.5       0.5
## RP11-34P13.8       0.5       0.5
## ...                ...       ...
## AC233755.2    0.500000       0.5
## AC233755.1    0.500000       0.5
## AC240274.1    0.482759       0.5
## AC213203.1    0.500000       0.5
## FAM231B       0.500000       0.5
```

Say we want to identify the genes that distinguish cluster 12 from other clusters with high _LYZ_ expression.
We subset `full.AUC` to the relevant comparisons and sort on our summary statistic of choice to obtain a ranking of markers within this subset.
This allows us to easily characterize subtle differences between closely related clusters.
To illustrate, we use the smallest rank from `computeMinRank()` to identify the top DE genes in cluster 12 compared to the other _LYZ_-high clusters (Figure \@ref(fig:pbmc-grouped-ranked-heat2)). 


``` r
lyz.high <- c("2", "4", "7", "12") # based on inspection of the previous Figure.
subset <- chosen$full.AUC[,colnames(chosen$full.AUC) %in% lyz.high]
to.show <- subset[computeMinRank(subset) <= 10,]
to.show
```

```
## DataFrame with 18 rows and 3 columns
##                       2         4         7
##               <numeric> <numeric> <numeric>
## CTSS           0.561257  0.791862  0.937902
## S100A9         0.588796  0.791569  0.959475
## S100A12        0.585836  0.712822  0.873311
## S100A8         0.747258  0.887295  0.976131
## S100A6         0.541432  0.810597  0.863344
## ...                 ...       ...       ...
## RP11-1143G9.4  0.746922  0.862851  0.869049
## FOS            0.784564  0.868267  0.943475
## H3F3B          0.689332  0.677693  0.734033
## JUND           0.627234  0.616803  0.733115
## FTL            0.574398  0.904567  0.993049
```

``` r
plotGroupedHeatmap(sce.pbmc[,colLabels(sce.pbmc) %in% lyz.high],
    features=rownames(to.show), group="label", center=TRUE, zlim=c(-3, 3))
```

<div class="figure">
<img src="marker-detection_files/figure-html/pbmc-grouped-ranked-heat2-1.png" alt="Heatmap of the centered average log-expression values for the top potential marker genes for cluster 12 relative to other _LYZ_-high clusters in the PBMC dataset. The set of markers was selected as those genes with AUC-derived min-ranks less than or equal to 10." width="672" />
<p class="caption">(\#fig:pbmc-grouped-ranked-heat2)Heatmap of the centered average log-expression values for the top potential marker genes for cluster 12 relative to other _LYZ_-high clusters in the PBMC dataset. The set of markers was selected as those genes with AUC-derived min-ranks less than or equal to 10.</p>
</div>



Similarly, we can use the full set of effect sizes to define our own summary statistic if the precomputed measures are too coarse.
For example, we may be interested in markers that are upregulated against some percentage - say, 80% - of other clusters.
This improves the cluster specificity of the ranking by being more stringent than the median yet not as strigent as the minimum.
We achieve this by computing and sorting on the 20^th^ percentile of effect sizes, as shown below.


``` r
stat <- rowQuantiles(as.matrix(chosen$full.AUC), p=0.2)
chosen[order(stat, decreasing=TRUE), 1:4] # just showing the basic stats for brevity.
```

```
## DataFrame with 33694 rows and 4 columns
##               self.average other.average self.detected other.detected
##                  <numeric>     <numeric>     <numeric>      <numeric>
## S100A8             6.25358      1.504959      0.983607       0.611945
## S100A9             5.95360      1.767483      1.000000       0.685216
## LYZ                5.87025      1.879720      1.000000       0.631549
## RP11-1143G9.4      3.81136      0.622745      0.983607       0.304197
## CTSS               3.79308      1.370146      0.950820       0.605402
## ...                    ...           ...           ...            ...
## RPL3               2.90331       4.13328      0.885246       0.927486
## RPS12              2.81598       4.33851      0.852459       0.919220
## RPS27              4.69602       5.30282      1.000000       0.949728
## RPS29              3.22834       4.40025      0.934426       0.918752
## RPL23A             2.51467       3.86930      0.901639       0.905495
```

## Using a log-fold change threshold

The Cohen's $d$ and AUC calculations consider both the magnitude of the difference between clusters as well as the variability within each cluster.
If the variability is lower, it is possible for a gene to have a large effect size even if the magnitude of the difference is small.
These genes tend to be somewhat uninformative for cell type identification despite their strong differential expression (e.g., ribosomal protein genes).
We would prefer genes with larger log-fold changes between clusters, even if they have higher variability.

To favor the detection of such genes, we can compute the effect sizes relative to a log-fold change threshold by setting `lfc=` in `scoreMarkers()`.
The definition of Cohen's $d$ is generalized to the standardized difference between the observed log-fold change and the specified `lfc` threshold.
Similarly, the AUC is redefined as the probability of randomly picking an expression value from one cluster that is greater than a random value from the other cluster plus `lfc`.
A large positive Cohen's $d$ and an AUC above 0.5 can only be obtained if the observed log-fold change between clusters is significantly greater than `lfc`.
We demonstrate below by obtaining the top markers for cluster 6 in the PBMC dataset with `lfc=2` (Figure \@ref(fig:pbmc-markers-5-lfc)).


``` r
marker.info.lfc <- scoreMarkers(sce.pbmc, colLabels(sce.pbmc), lfc=2)
chosen2 <- marker.info.lfc[["6"]] # another cluster for some variety.
chosen2 <- chosen2[order(chosen2$mean.AUC, decreasing=TRUE),]
chosen2[,c("self.average", "other.average", "mean.AUC")]
```

```
## DataFrame with 33694 rows and 3 columns
##            self.average other.average  mean.AUC
##               <numeric>     <numeric> <numeric>
## CCL5            3.86291      1.107731  0.744410
## NKG7            3.31704      0.873797  0.682214
## IL32            3.31048      1.126574  0.618120
## GZMA            1.98416      0.449038  0.469145
## TRAC            2.25625      0.892877  0.422834
## ...                 ...           ...       ...
## AC233755.2    0.0000000     0.0000000         0
## AC233755.1    0.0000000     0.0000000         0
## AC240274.1    0.0100134     0.0075477         0
## AC213203.1    0.0000000     0.0000000         0
## FAM231B       0.0000000     0.0000000         0
```

``` r
plotDots(sce.pbmc, rownames(chosen2)[1:10], group="label")
```

<div class="figure">
<img src="marker-detection_files/figure-html/pbmc-markers-6-lfc-1.png" alt="Dot plot of the top potential marker genes (as determined by the mean AUC) for cluster 6 in the PBMC dataset. Each row corrresponds to a marker gene and each column corresponds to a cluster. The size of each dot represents the proportion of cells with detected expression of the gene in the cluster, while the color is proportional to the average expression across all cells in that cluster." width="576" />
<p class="caption">(\#fig:pbmc-markers-6-lfc)Dot plot of the top potential marker genes (as determined by the mean AUC) for cluster 6 in the PBMC dataset. Each row corrresponds to a marker gene and each column corresponds to a cluster. The size of each dot represents the proportion of cells with detected expression of the gene in the cluster, while the color is proportional to the average expression across all cells in that cluster.</p>
</div>

Note that the interpretation of the AUC and Cohen's $d$ becomes slightly more complicated when `lfc` is non-zero.
If `lfc` is positive, a positive Cohen's $d$ and an AUC above 0.5 represents upregulation. 
However, a negative Cohen's $d$ or AUC below 0.5 may not represent downregulation; it may just indicate that the observed log-fold change is less than the specified `lfc`.
The converse applies when `lfc` is negative, where the only conclusive interpretation occurs for downregulated genes.
For the most part, this complication is not too problematic for routine marker detection, as we are mostly interested in upregulated genes with large positive Cohen's $d$ and AUCs above 0.5.

## Handling blocking factors {#marker-batch}

Large studies may contain factors of variation that are known and not interesting (e.g., batch effects, sex differences).
If these are not modelled, they can interfere with marker gene detection - most obviously by inflating the variance within each cluster, but also by distorting the log-fold changes if the cluster composition varies across levels of the blocking factor.
To avoid these issues, we specify the blocking factor via the `block=` argument, as demonstrated below for the 416B data set.

<button class="rebook-collapse">View set-up code ([Workflow Chapter 1](http://bioconductor.org/books/3.20/OSCA.workflows/lun-416b-cell-line-smart-seq2.html#lun-416b-cell-line-smart-seq2))</button>
<div class="rebook-content">

```r
#--- loading ---#
library(scRNAseq)
sce.416b <- LunSpikeInData(which="416b") 
sce.416b$block <- factor(sce.416b$block)

#--- gene-annotation ---#
library(AnnotationHub)
ens.mm.v97 <- AnnotationHub()[["AH73905"]]
rowData(sce.416b)$ENSEMBL <- rownames(sce.416b)
rowData(sce.416b)$SYMBOL <- mapIds(ens.mm.v97, keys=rownames(sce.416b),
    keytype="GENEID", column="SYMBOL")
rowData(sce.416b)$SEQNAME <- mapIds(ens.mm.v97, keys=rownames(sce.416b),
    keytype="GENEID", column="SEQNAME")

library(scater)
rownames(sce.416b) <- uniquifyFeatureNames(rowData(sce.416b)$ENSEMBL, 
    rowData(sce.416b)$SYMBOL)

#--- quality-control ---#
mito <- which(rowData(sce.416b)$SEQNAME=="MT")
stats <- perCellQCMetrics(sce.416b, subsets=list(Mt=mito))
qc <- quickPerCellQC(stats, percent_subsets=c("subsets_Mt_percent",
    "altexps_ERCC_percent"), batch=sce.416b$block)
sce.416b <- sce.416b[,!qc$discard]

#--- normalization ---#
library(scran)
sce.416b <- computeSumFactors(sce.416b)
sce.416b <- logNormCounts(sce.416b)

#--- variance-modelling ---#
dec.416b <- modelGeneVarWithSpikes(sce.416b, "ERCC", block=sce.416b$block)
chosen.hvgs <- getTopHVGs(dec.416b, prop=0.1)

#--- batch-correction ---#
library(limma)
assay(sce.416b, "corrected") <- removeBatchEffect(logcounts(sce.416b), 
    design=model.matrix(~sce.416b$phenotype), batch=sce.416b$block)

#--- dimensionality-reduction ---#
sce.416b <- runPCA(sce.416b, ncomponents=10, subset_row=chosen.hvgs,
    exprs_values="corrected", BSPARAM=BiocSingular::ExactParam())

set.seed(1010)
sce.416b <- runTSNE(sce.416b, dimred="PCA", perplexity=10)

#--- clustering ---#
my.dist <- dist(reducedDim(sce.416b, "PCA"))
my.tree <- hclust(my.dist, method="ward.D2")

library(dynamicTreeCut)
my.clusters <- unname(cutreeDynamic(my.tree, distM=as.matrix(my.dist),
    minClusterSize=10, verbose=0))
colLabels(sce.416b) <- factor(my.clusters)
```

</div>


``` r
m.out <- scoreMarkers(sce.416b, colLabels(sce.416b), block=sce.416b$block)
```

For each gene, each pairwise comparison between clusters is performed separately in each level of the blocking factor - in this case, the plate of origin.
By comparing within each batch, we cancel out any batch effects so that they are not conflated with the biological differences between subpopulations.
The effect sizes are then averaged across batches to obtain a single value per comparison, using a weighted mean that accounts for the number of cells involved in the comparison in each batch.
A similar correction is applied to the mean log-expression and proportion of detected cells inside and outside each cluster.


``` r
demo <- m.out[["1"]] 
ordered <- demo[order(demo$median.logFC.cohen, decreasing=TRUE),]
ordered[,1:4]
```

```
## DataFrame with 46604 rows and 4 columns
##         self.average other.average self.detected other.detected
##            <numeric>     <numeric>     <numeric>      <numeric>
## Myh11        4.03436      0.861019      0.988132       0.303097
## Cd200r3      7.97667      3.524762      0.977675       0.624507
## Pi16         6.27654      2.644421      0.957126       0.530395
## Actb        15.48533     14.808584      1.000000       1.000000
## Ctsd        11.61247      9.130141      1.000000       1.000000
## ...              ...           ...           ...            ...
## Spc24      0.4772577       5.03548      0.222281       0.862153
## Ska1       0.0787421       4.43426      0.118743       0.773950
## Pimreg     0.5263611       5.35494      0.258150       0.910706
## Birc5      1.5580536       7.07230      0.698746       0.976929
## Ccna2      0.9664521       6.55243      0.554104       0.948520
```

``` r
plotExpression(sce.416b, features=rownames(ordered)[1:6],
    x="label", colour_by="block")
```

<div class="figure">
<img src="marker-detection_files/figure-html/blocked-markers-416b-1.png" alt="Distribution of expression values across clusters for the top potential marker genes from cluster 1 in the 416B dataset. Each point represents a cell and is colored by the batch of origin." width="576" />
<p class="caption">(\#fig:blocked-markers-416b)Distribution of expression values across clusters for the top potential marker genes from cluster 1 in the 416B dataset. Each point represents a cell and is colored by the batch of origin.</p>
</div>

The `block=` argument works for all effect sizes shown above and is robust to differences in the log-fold changes or variance between batches.
However, it assumes that each pair of clusters is present in at least one batch.
In scenarios where cells from two clusters never co-occur in the same batch, the associated pairwise comparison will be impossible and is ignored during calculation of summary statistics.

<!--
### Using the `design=` argument

Another approach is to define a design matrix containing the batch of origin as the sole factor.
`findMarkers()` will then fit a linear model to the log-expression values, similar to the use of *[limma](https://bioconductor.org/packages/3.20/limma)* for bulk RNA sequencing data [@ritchie2015limma].
This handles situations where multiple batches contain unique clusters, as comparisons can be implicitly performed via shared cell types in each batch.
There is also a slight increase in power when information is shared across clusters for variance estimation.


``` r
# Setting up the design matrix (we remove intercept for full rank
# in the final design matrix with the cluster-specific terms).
design <- model.matrix(~sce.416b$block)
design <- design[,-1,drop=FALSE]

m.alt <- findMarkers(sce.416b, design=design, direction="up")
demo <- m.alt[["1"]]
demo[demo$Top <= 5,1:4]
```

```
## DataFrame with 12 rows and 4 columns
##                          Top     p.value         FDR summary.logFC
##                    <integer>   <numeric>   <numeric>     <numeric>
## Gm6977                     1 7.15187e-24 8.77120e-21      0.810553
## Myh11                      1 4.56882e-64 2.12925e-59      4.381806
## Tmsb4x                     2 9.48997e-46 2.21135e-41      1.478213
## Cd63                       2 1.80446e-15 7.85933e-13      0.813016
## Cd200r3                    2 2.40861e-45 3.74170e-41      6.684003
## ...                      ...         ...         ...           ...
## Actb                       4 5.61751e-36 2.90887e-32      0.961762
## Ctsd                       4 2.08646e-42 2.43094e-38      2.893014
## Fth1                       4 1.83949e-23 2.14319e-20      0.797407
## Ccl9                       5 1.75378e-30 3.71514e-27      5.396347
## CBFB-MYH11-mcherry         5 9.09026e-39 8.47285e-35      3.017758
```

The use of a linear model makes some strong assumptions, necessitating some caution when interpreting the results.
If the batch effect is not consistent across clusters, the variance will be inflated and the log-fold change estimates will be distorted.
Variances are also assumed to be equal across groups, which is not true in general.
In particular, the presence of clusters in which a gene is silent will shrink the residual variance towards zero, preventing the model from penalizing genes with high variance in other clusters.
Thus, we generally recommend the use of `block=` where possible.
-->

## Session Info {-}

<button class="rebook-collapse">View session info</button>
<div class="rebook-content">
```
R version 4.4.2 (2024-10-31)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 24.04.1 LTS

Matrix products: default
BLAS:   /home/biocbuild/bbs-3.20-bioc/R/lib/libRblas.so 
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_GB              LC_COLLATE=C              
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

time zone: America/New_York
tzcode source: system (glibc)

attached base packages:
[1] stats4    stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] scater_1.34.0               ggplot2_3.5.1              
 [3] scran_1.34.0                scuttle_1.16.0             
 [5] SingleCellExperiment_1.28.1 SummarizedExperiment_1.36.0
 [7] Biobase_2.66.0              GenomicRanges_1.58.0       
 [9] GenomeInfoDb_1.42.0         IRanges_2.40.0             
[11] S4Vectors_0.44.0            BiocGenerics_0.52.0        
[13] MatrixGenerics_1.18.0       matrixStats_1.4.1          
[15] BiocStyle_2.34.0            rebook_1.16.0              

loaded via a namespace (and not attached):
 [1] gridExtra_2.3           CodeDepends_0.6.6       rlang_1.1.4            
 [4] magrittr_2.0.3          compiler_4.4.2          dir.expiry_1.14.0      
 [7] vctrs_0.6.5             pkgconfig_2.0.3         crayon_1.5.3           
[10] fastmap_1.2.0           XVector_0.46.0          labeling_0.4.3         
[13] utf8_1.2.4              rmarkdown_2.29          graph_1.84.0           
[16] UCSC.utils_1.2.0        ggbeeswarm_0.7.2        xfun_0.49              
[19] bluster_1.16.0          zlibbioc_1.52.0         cachem_1.1.0           
[22] beachmat_2.22.0         jsonlite_1.8.9          DelayedArray_0.32.0    
[25] BiocParallel_1.40.0     irlba_2.3.5.1           parallel_4.4.2         
[28] cluster_2.1.6           R6_2.5.1                RColorBrewer_1.1-3     
[31] bslib_0.8.0             limma_3.62.1            jquerylib_0.1.4        
[34] Rcpp_1.0.13-1           bookdown_0.41           knitr_1.49             
[37] Matrix_1.7-1            igraph_2.1.1            tidyselect_1.2.1       
[40] abind_1.4-8             yaml_2.3.10             viridis_0.6.5          
[43] codetools_0.2-20        lattice_0.22-6          tibble_3.2.1           
[46] withr_3.0.2             evaluate_1.0.1          pillar_1.9.0           
[49] BiocManager_1.30.25     filelock_1.0.3          generics_0.1.3         
[52] munsell_0.5.1           scales_1.3.0            glue_1.8.0             
[55] metapod_1.14.0          pheatmap_1.0.12         tools_4.4.2            
[58] BiocNeighbors_2.0.0     ScaledMatrix_1.14.0     locfit_1.5-9.10        
[61] XML_3.99-0.17           cowplot_1.1.3           grid_4.4.2             
[64] edgeR_4.4.0             colorspace_2.1-1        GenomeInfoDbData_1.2.13
[67] beeswarm_0.4.0          BiocSingular_1.22.0     vipor_0.4.7            
[70] cli_3.6.3               rsvd_1.0.5              rappdirs_0.3.3         
[73] fansi_1.0.6             S4Arrays_1.6.0          viridisLite_0.4.2      
[76] dplyr_1.1.4             gtable_0.3.6            sass_0.4.9             
[79] digest_0.6.37           SparseArray_1.6.0       ggrepel_0.9.6          
[82] dqrng_0.4.1             farver_2.1.2            htmltools_0.5.8.1      
[85] lifecycle_1.0.4         httr_1.4.7              statmod_1.5.0          
```
</div>