CRUK CI Summer School 2020 - introduction to single-cell RNA-seq analysis

Cluster marker genes

1 Cluster marker genes

projDir <- "/mnt/scratcha/bioinformatics/baller01/20200511_FernandesM_ME_crukBiSs2020"
outDirBit <- "AnaWiSce/Attempt1"
nbPcToComp <- 50

library(ggplot2)
library(scater)

library(RColorBrewer)
library(pheatmap)
fontsize <- theme(axis.text=element_text(size=12), axis.title=element_text(size=16))

Source: we will follow the OSCA chapter on marker detection (with some of its text copied here with little modification). See also the Hemberg group chapter on differential analysis section.

To interpret our clustering results, 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 most obvious case, the marker genes for each cluster are a priori 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.

Identification of marker genes is usually based around the retrospective detection of differential expression between clusters. Genes that are more strongly DE are more likely to have caused separate clustering of cells in the first place. Several different statistical tests are available to quantify the differences in expression profiles, and different approaches can be used to consolidate test results into a single ranking of genes for each cluster.

1.1 Load data

We will load the R file keeping the SCE (SingleCellExperiment) object with the normalised counts for 500 cells per sample used for feature selection and dimensionality reduction then clustering.

setName <- "caron"
setSuf <- "_5hCellPerSpl"

# Read object in:
tmpFn <- sprintf("%s/%s/Robjects/%s_sce_nz_postDeconv%s_clustered.Rds", projDir, outDirBit, setName, setSuf)
print(tmpFn)

[1] "/mnt/scratcha/bioinformatics/baller01/20200511_FernandesM_ME_crukBiSs2020/AnaWiSce/Attempt1/Robjects/caron_sce_nz_postDeconv_5hCellPerSpl_clustered.Rds"

if(!file.exists(tmpFn))
{
    knitr::knit_exit()
}
sce <- readRDS(tmpFn)
sce

class: SingleCellExperiment 
dim: 18372 5500 
metadata(0):
assays(2): counts logcounts
rownames(18372): ENSG00000238009 ENSG00000237491 ... ENSG00000275063
  ENSG00000271254
rowData names(11): ensembl_gene_id external_gene_name ... detected
  gene_sparsity
colnames: NULL
colData names(22): Sample Barcode ... cluster kmeans10
reducedDimNames(3): PCA TSNE UMAP
altExpNames(0):

#head(rowData(sce))

#any(duplicated(rowData(sce)$ensembl_gene_id))
# some function(s) used below complain about 'strand' already being used in row data,
# so rename that column now:
#colnames(rowData(sce))[colnames(rowData(sce)) == "strand"] <- "strandNum"
#assayNames(sce)
#reducedDimNames(sce)

1.1.1 Detecting genes differentially expressed between clusters

1.1.1.1 Differential expression analysis

We will identify genes for each cluster whose expression differ to that of other clusters, using findMarkers(). It fits a linear model to the log-expression values for each gene using limma [@doi:10.1093/nar/gkv007] and allows testing for differential expression in each cluster compared to the others while accounting for known, uninteresting factors.

sce$clusterStg <- factor(paste0("c", sce$cluster),
             levels = paste0("c", levels( sce$cluster)) )
markers <- findMarkers(sce, sce$clusterStg)

Results are compiled in a single table per cluster that stores the outcome of comparisons against the other clusters. One can then select differentially expressed genes from each pairwise comparison between clusters.

We will define a set of genes for cluster 1 by selecting the top 10 genes of each comparison, and check the test output, eg adjusted p-values and log-fold changes.

# get output table for clsuter 1:
marker.set <- markers[["c1"]]
head(marker.set, 10)

DataFrame with 10 rows and 25 columns
                      Top      p.value          FDR summary.logFC   logFC.c2
                <integer>    <numeric>    <numeric>     <numeric>  <numeric>
ENSG00000188643         1 9.60012e-203 4.40933e-200      1.329499  1.3372126
ENSG00000115977         1  1.21980e-81  7.70111e-80      0.521398  0.4081095
ENSG00000081189         1 2.99575e-215 1.83460e-212      1.551063  1.2176288
ENSG00000019582         1  0.00000e+00  0.00000e+00      3.889273 -0.0441531
ENSG00000204287         1  0.00000e+00  0.00000e+00      3.305402  0.2934284
ENSG00000198034         1 5.45114e-202 2.44264e-199      1.856998  2.0279499
ENSG00000244734         1  0.00000e+00  0.00000e+00    -13.510067 -1.2905633
ENSG00000100721         1  0.00000e+00  0.00000e+00      3.271716  1.0963887
ENSG00000108298         1 7.60633e-160 2.11733e-157     -1.447336  0.2700361
ENSG00000169575         1  0.00000e+00  0.00000e+00      2.745276  1.6570717
                 logFC.c3  logFC.c4  logFC.c5   logFC.c6  logFC.c7  logFC.c8
                <numeric> <numeric> <numeric>  <numeric> <numeric> <numeric>
ENSG00000188643  1.337213  1.337213  1.327932   1.337213  1.329499  1.337213
ENSG00000115977  0.237450  0.432485 -0.154086   0.492783  0.417690  0.428059
ENSG00000081189  1.572606  1.570229  1.551063   1.495518  0.970556  1.572606
ENSG00000019582  3.889273  3.959482  3.028084   2.143768  1.082220  3.984027
ENSG00000204287  3.513484  3.545370  3.305402   2.039801  1.269541  3.533434
ENSG00000198034  0.194784  0.767665  0.251410   2.672021  0.431460  0.688086
ENSG00000244734 -1.320825  0.326670 -2.386501 -13.510067 -1.373549  0.520068
ENSG00000100721  3.420225  3.420225  3.271716   2.923737  3.317503  3.418893
ENSG00000108298 -0.519027 -0.198044 -0.927478   1.326099 -0.644297 -0.270588
ENSG00000169575  2.827075  2.784402  2.745276   2.541961  2.705674  2.756574
                  logFC.c9 logFC.c10 logFC.c11  logFC.c12 logFC.c13 logFC.c14
                 <numeric> <numeric> <numeric>  <numeric> <numeric> <numeric>
ENSG00000188643  1.3104313  1.268088 -0.150252  1.3152143  1.328237  1.319371
ENSG00000115977  0.4922509  0.397145  0.188344  0.4337223  0.382133  0.328245
ENSG00000081189  0.0375915  0.254803  0.381641  0.0594653  0.332399  1.140873
ENSG00000019582 -1.4604987 -0.672980 -0.336633 -0.7821347 -0.309313  0.961137
ENSG00000204287 -0.6838107 -0.134829  0.209194 -1.1630891 -0.759056  0.899975
ENSG00000198034  0.9339553  1.452751 -0.900094  0.5857575  0.879309  1.145308
ENSG00000244734 -2.3195297 -2.706101  0.754421  0.3265927  0.402572 -0.525518
ENSG00000100721  2.0020345  1.045096 -0.187818  0.3844558  0.990486  1.617084
ENSG00000108298 -0.4605559  0.298537 -1.447336 -0.5890394 -0.217389 -0.192680
ENSG00000169575  2.7098921  1.641334 -0.225350  2.8145924  2.815142  0.929480
                 logFC.c15  logFC.c16 logFC.c17 logFC.c18 logFC.c19 logFC.c20
                 <numeric>  <numeric> <numeric> <numeric> <numeric> <numeric>
ENSG00000188643  0.8893410   1.337213  1.337213  1.337213  1.337213  1.337213
ENSG00000115977  0.4119520   0.411171  0.446193 -0.068711  0.337583  0.363816
ENSG00000081189  0.0734917   1.550773  1.374269  1.531552  1.031078  1.572606
ENSG00000019582 -0.9752284   3.461483  3.737241  2.461530  1.019503  3.986356
ENSG00000204287 -0.7117804   3.042733  3.457508  3.066186  1.078753  3.546942
ENSG00000198034  1.8569982   2.470039  0.858512  1.746264  1.646451  2.670146
ENSG00000244734 -0.2882371 -10.599858 -8.045918 -2.401120 -2.840030  0.234350
ENSG00000100721  0.6751496   3.195810  3.374781  3.291750  3.343118  3.420225
ENSG00000108298  0.7091609   0.909686 -0.352559  0.392223  0.115234  1.714652
ENSG00000169575  1.8686088   2.723215  2.844063  2.800898  2.775603  2.844063
                logFC.c21  logFC.c22
                <numeric>  <numeric>
ENSG00000188643  1.337213   1.117229
ENSG00000115977 -0.353520   0.521398
ENSG00000081189  1.572606   0.779136
ENSG00000019582  3.745102   0.501519
ENSG00000204287  3.517578   0.863246
ENSG00000198034 -0.174714   1.155415
ENSG00000244734  0.651856 -10.751479
ENSG00000100721  3.404903   1.940488
ENSG00000108298 -1.456636   0.453718
ENSG00000169575  2.827143   1.994102

# add gene annotation:
tmpDf <- marker.set
tmpDf$ensembl_gene_id <- rownames(tmpDf)
tmpDf2 <- base::merge(tmpDf, rowData(sce), by="ensembl_gene_id", all.x=TRUE, all.y=F, sort=F)

Write Table to file:

tmpFn <- sprintf("%s/%s/Tables/%s_sce_nz_postDeconv%s_c1.tsv", projDir, outDirBit, setName, setSuf)
print(tmpFn)

[1] "/mnt/scratcha/bioinformatics/baller01/20200511_FernandesM_ME_crukBiSs2020/AnaWiSce/Attempt1/Tables/caron_sce_nz_postDeconv_5hCellPerSpl_c1.tsv"

write.table(tmpDf2, file=tmpFn, sep="\t", quote=FALSE, row.names=FALSE)

Show expression of merker on t-SNE and UMAP:

tsne1 <- plotTSNE(sce, colour_by=tmpDf2[1,"ensembl_gene_id"]) + fontsize
umap1 <- plotUMAP(sce, colour_by=tmpDf2[1,"ensembl_gene_id"]) + fontsize
tsne1

umap1

Gene set enrichment analyses used for bulk RNA-seq may be used to characterise clusters further.

1.1.1.2 Heatmap

As for bulk RNA, differences in expression profiles of the top genes can be visualised with a heatmap.

# select some top genes:
top.markers <- rownames(marker.set)[marker.set$Top <= 10]

# have matrix to annotate sample with cluster and sample:
tmpData <- logcounts(sce)[top.markers,]
# concat sample and barcode names to make unique name across the whole data set
tmpCellNames <- paste(colData(sce)$Sample.Name, colData(sce)$Barcode, sep="_")
# use these to namecolumn of matrix the show as heatmap:
colnames(tmpData) <- tmpCellNames # colData(sce)$Barcode                    

# columns annotation with cell name:
mat_col <- data.frame(cluster = sce$cluster,
              sample = sce$Sample.Name,
              type = sce$source_name
        )
rownames(mat_col) <- colnames(tmpData)
rownames(mat_col) <- tmpCellNames # colData(sce)$Barcode

# Prepare colours for clusters:
colourCount = length(unique(sce$cluster))
getPalette = colorRampPalette(brewer.pal(9, "Set1"))

mat_colors <- list(group = getPalette(colourCount))
names(mat_colors$group) <- unique(sce$cluster)

# plot heatmap:
pheatmap(tmpData,
           border_color      = NA,
           show_colnames     = FALSE,
           show_rownames     = FALSE,
           drop_levels       = TRUE,
           annotation_col    = mat_col,
           annotation_colors = mat_colors
           )

One can sort both the gene and sample dendrograms to improve the heatmap.

library(dendsort)

mat <- tmpData
mat_cluster_cols <- hclust(dist(t(mat)))

sort_hclust <- function(...) as.hclust(dendsort(as.dendrogram(...)))

mat_cluster_cols <- sort_hclust(mat_cluster_cols)
#plot(mat_cluster_cols, main = "Sorted Dendrogram", xlab = "", sub = "")

mat_cluster_rows <- sort_hclust(hclust(dist(mat)))

pheatmap(tmpData,
           border_color      = NA,
           show_colnames     = FALSE,
           show_rownames     = FALSE,
           drop_levels       = TRUE,
           annotation_col    = mat_col,
           annotation_colors = mat_colors,
           cluster_cols      = mat_cluster_cols,
           cluster_rows      = mat_cluster_rows
         )

To demonstrate how to interpret the results, we will use cluster 9 as our cluster of interest. The relevant DataFrame contains log2-fold changes of expression in cluster 9 over each other cluster, along with several statistics obtained by combining p-values (Simes 1986) across the pairwise comparisons involving cluster 9.

chosen <- "c9"
interesting <- markers[[chosen]]
print(colnames(interesting))

 [1] "Top"           "p.value"       "FDR"           "summary.logFC"
 [5] "logFC.c1"      "logFC.c2"      "logFC.c3"      "logFC.c4"     
 [9] "logFC.c5"      "logFC.c6"      "logFC.c7"      "logFC.c8"     
[13] "logFC.c10"     "logFC.c11"     "logFC.c12"     "logFC.c13"    
[17] "logFC.c14"     "logFC.c15"     "logFC.c16"     "logFC.c17"    
[21] "logFC.c18"     "logFC.c19"     "logFC.c20"     "logFC.c21"    
[25] "logFC.c22"    

Of particular interest is the Top field. The set of genes with Top ≤X is the union of the top X genes (ranked by p-value) from each pairwise comparison involving cluster 9. For example, the set of all genes with Top values of 1 contains the gene with the lowest p-value from each comparison. Similarly, the set of genes with Top values less than or equal to 10 contains the top 10 genes from each comparison. The Top field represents findMarkers()’s approach to consolidating multiple pairwise comparisons into a single ranking for each cluster; each DataFrame produced by findMarkers() will order genes based on the Top value by default.

interesting[1:10,1:4]

DataFrame with 10 rows and 4 columns
                      Top      p.value          FDR summary.logFC
                <integer>    <numeric>    <numeric>     <numeric>
ENSG00000117632         1 1.39119e-173 1.06496e-170      -2.68610
ENSG00000197956         1 2.23405e-155 1.36813e-152      -4.00260
ENSG00000177954         1 4.17810e-185 4.26445e-182       1.96167
ENSG00000019582         1 2.14971e-315 3.94945e-311       4.48858
ENSG00000227507         1  3.44902e-71  3.31756e-69       2.16907
ENSG00000204472         1 2.18155e-175 1.82179e-172      -1.74721
ENSG00000204287         1 7.85309e-223 1.10982e-219       4.20139
ENSG00000211751         1 1.49359e-231 2.49457e-228      -3.07451
ENSG00000244734         1 3.12618e-243 7.17927e-240     -11.19054
ENSG00000223609         1 1.39274e-305 1.27937e-301      -8.44988

We use the Top field to identify a set of genes that is guaranteed to distinguish cluster 9 from any other cluster. Here, we examine the top 6 genes from each pairwise comparison.

best.set <- interesting[interesting$Top <= 6,]
logFCs <- getMarkerEffects(best.set)
logFCs.ens <- rownames(logFCs)
rownames(logFCs) <- rowData(sce)[rownames(logFCs), "Symbol"]

library(pheatmap)
pheatmap(logFCs, breaks=seq(-5, 5, length.out=101))

1.1.2 Using the log-fold change

Our previous findMarkers() call considers both up- and downregulated genes to be potential markers. However, downregulated genes are less appealing as markers as it is more difficult to interpret and experimentally validate an absence of expression. To focus on up-regulated markers, we can instead perform a one-sided t-test to identify genes that are upregulated in each cluster compared to the others. This is achieved by setting direction=“up” in the findMarkers() call.

markers.up <- findMarkers(sce, groups=sce$clusterStg, direction="up")
interesting.up <- markers.up[[chosen]]
interesting.up[1:10,1:4]

DataFrame with 10 rows and 4 columns
                      Top      p.value          FDR summary.logFC
                <integer>    <numeric>    <numeric>     <numeric>
ENSG00000177954         1 2.08905e-185 1.27934e-181      1.961669
ENSG00000019582         1 1.07486e-315 1.97473e-311      4.488583
ENSG00000227507         1  1.72451e-71  7.04061e-69      2.169067
ENSG00000204287         1 3.92654e-223 3.60692e-219      4.201388
ENSG00000023445         1  5.16434e-22  4.91603e-20      0.629716
ENSG00000133639         1 2.31676e-112 3.04026e-109      2.929176
ENSG00000211899         1 7.96547e-113 1.12570e-109      4.079709
ENSG00000104894         1 1.36061e-118 2.49970e-115      2.847488
ENSG00000169442         2 8.26415e-172 3.79572e-168      3.270898
ENSG00000196126         2 6.93749e-124 1.41617e-120      2.956772

The t-test also allows us to specify a non-zero log-fold change as the null hypothesis. This allows us to consider the magnitude of the log-fold change in our p-value calculations, in a manner that is more rigorous than simply filtering directly on the log-fold changes (McCarthy and Smyth 2009). (Specifically, a simple threshold does not consider the variance and can enrich for genes that have both large log-fold changes and large variances.) We perform this by setting lfc= in our findMarkers() call - when combined with direction=, this tests for genes with log-fold changes that are significantly greater than 1:

markers.up2 <- findMarkers(sce, groups=sce$clusterStg, direction="up", lfc=1)
interesting.up2 <- markers.up2[[chosen]]
interesting.up2[1:10,1:4]

DataFrame with 10 rows and 4 columns
                      Top      p.value          FDR summary.logFC
                <integer>    <numeric>    <numeric>     <numeric>
ENSG00000169442         1 1.25317e-121 7.67442e-118       3.27090
ENSG00000177954         1  2.67850e-72  6.15118e-69       1.96167
ENSG00000019582         1 1.25616e-250 2.30782e-246       4.48858
ENSG00000204287         1 3.06630e-181 2.81671e-177       4.20139
ENSG00000244734         1  1.11320e-44  1.27823e-41       3.07395
ENSG00000211899         1  7.01040e-84  3.21988e-80       4.07971
ENSG00000104894         1  8.98214e-74  2.75033e-70       2.86659
ENSG00000227507         2  1.23004e-29  8.36975e-27       2.20060
ENSG00000196126         2  1.38266e-80  5.08044e-77       2.95677
ENSG00000197958         2  7.42120e-51  1.04879e-47       2.19487

These two settings yield a more focused set of candidate marker genes that are upregulated in cluster 9.

best.set <- interesting.up2[interesting.up2$Top <= 5,]
logFCs <- getMarkerEffects(best.set)
logFCs.ens <- rownames(logFCs)
rownames(logFCs) <- rowData(sce)[rownames(logFCs), "Symbol"]

library(pheatmap)
pheatmap(logFCs, breaks=seq(-5, 5, length.out=101))

Of course, this increased stringency is not without cost. If only upregulated genes are requested from findMarkers(), any cluster defined by downregulation of a marker gene will not contain that gene among the top set of features in its DataFrame. This is occasionally relevant for subtypes or other states that are distinguished by high versus low expression of particular genes. Similarly, setting an excessively high log-fold change threshold may discard otherwise useful genes. For example, a gene upregulated in a small proportion of cells of a cluster will have a small log-fold change but can still be an effective marker if the focus is on specificity rather than sensitivity.

1.1.3 Finding cluster-specific markers

By default, findMarkers() will give a high ranking to genes that are differentially expressed in any pairwise comparison. This is because a gene only needs a very low p -value in a single pairwise comparison to achieve a low Top value. A more stringent approach would only consider genes that are differentially expressed in all pairwise comparisons involving the cluster of interest. To achieve this, we set pval.type=“all” in findMarkers() to use an intersection-union test (Berger and Hsu 1996) where the combined p-value for each gene is the maximum of the p-values from all pairwise comparisons. A gene will only achieve a low combined p-value if it is strongly DE in all comparisons to other clusters.

# We can combine this with 'direction='.
markers.up3 <- findMarkers(sce, groups=sce$clusterStg, pval.type="all", direction="up")
interesting.up3 <- markers.up3[[chosen]]
interesting.up3[1:10,1:3]

DataFrame with 10 rows and 3 columns
                    p.value         FDR summary.logFC
                  <numeric>   <numeric>     <numeric>
ENSG00000247982 2.40628e-17 4.42082e-13      1.095994
ENSG00000156738 4.71395e-16 4.33024e-12      2.032241
ENSG00000104921 1.21449e-13 7.43753e-10      0.556458
ENSG00000104894 4.58636e-12 2.10651e-08      2.243235
ENSG00000023445 3.40632e-11 1.25162e-07      0.493093
ENSG00000042980 4.57737e-10 1.40159e-06      0.389140
ENSG00000211898 4.71149e-09 1.23656e-05      1.255726
ENSG00000101017 4.96921e-08 1.14118e-04      0.284147
ENSG00000147535 9.18893e-08 1.87577e-04      0.599418
ENSG00000168081 2.10321e-07 3.70149e-04      0.225898

This strategy will only report genes that are highly specific to the cluster of interest. When it works, it can be highly effective as it generates a small focused set of candidate 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. With pval.type=“all”, neither Cd4 or Cd8 would be detected as subpopulation-specific markers because each gene is expressed in two subpopulations. In comparison, pval.type=“any” will detect both of these genes as they will be DE between at least one pair of subpopulations.

If pval.type=“all” is too stringent yet pval.type=“any” is too generous, a compromise is to set pval.type=“some”. For each gene, we apply the Holm-Bonferroni correction across its p -values and take the middle-most value as the combined p-value. This effectively tests the global null hypothesis that at least 50% of the individual pairwise comparisons exhibit no DE. We then rank the genes by their combined p-values to obtain an ordered set of marker candidates. The aim is to improve the conciseness of the top markers for defining a cluster while mitigating the risk of discarding useful genes that are not DE to all other clusters. The downside is that taking this compromise position sacrifices the theoretical guarantees offered at the other two extremes.

markers.up4 <- findMarkers(sce, groups=sce$clusterStg, pval.type="some", direction="up")
interesting.up4 <- markers.up4[[chosen]]
interesting.up4[1:10,1:3]

DataFrame with 10 rows and 3 columns
                    p.value         FDR summary.logFC
                  <numeric>   <numeric>     <numeric>
ENSG00000211899 4.26299e-85 7.83196e-81       3.41992
ENSG00000156738 2.99846e-83 2.75438e-79       2.27727
ENSG00000019582 1.13419e-81 6.94579e-78       3.60427
ENSG00000105369 6.60277e-78 3.03265e-74       2.14283
ENSG00000007312 3.19141e-73 1.17265e-69       2.16430
ENSG00000104894 2.74322e-69 8.39974e-66       1.75596
ENSG00000196126 4.62457e-65 1.21375e-61       1.67076
ENSG00000177954 8.06719e-57 1.85263e-53       1.35080
ENSG00000231389 7.58851e-56 1.54907e-52       1.85604
ENSG00000169442 9.41890e-53 1.73044e-49       1.23526

In both cases, a different method is used to compute the summary effect size compared to pval.type=“any”. For pval.type=“all”, the summary log-fold change is defined as that corresponding to the pairwise comparison with the largest p-value, while for pval.type=“some”, it is defined as the log-fold change for the comparison with the middle-most p-value. This reflects the calculation of the combined p-value and avoids focusing on genes with strong changes in only one comparison.

1.1.4  Using the Wilcoxon rank sum test

The Wilcoxon rank sum test (also known as the Wilcoxon-Mann-Whitney test, or WMW test) is another widely used method for pairwise comparisons between groups of observations. Its strength lies in the fact that it directly assesses separation between the expression distributions of different clusters. The WMW test statistic is proportional to the area-under-the-curve (AUC), i.e., the concordance probability, which is the probability of a random cell from one cluster having higher expression than a random cell from another cluster. In a pairwise comparison, AUCs of 1 or 0 indicate that the two clusters have perfectly separated expression distributions. Thus, the WMW test directly addresses the most desirable property of a candidate marker gene, while the t-test only does so indirectly via the difference in the means and the intra-group variance.

We perform WMW tests by again using the findMarkers() function, this time with test=“wilcox”. This returns a list of DataFrames containing ranked candidate markers for each cluster. The direction=, lfc= and pval.type= arguments can be specified and have the same interpretation as described for t-tests. We demonstrate below by detecting upregulated genes in each cluster with direction=“up”.

markers.wmw <- findMarkers(sce, groups=sce$clusterStg, test="wilcox", direction="up")
print(names(markers.wmw))

 [1] "c1"  "c2"  "c3"  "c4"  "c5"  "c6"  "c7"  "c8"  "c9"  "c10" "c11" "c12"
[13] "c13" "c14" "c15" "c16" "c17" "c18" "c19" "c20" "c21" "c22"

To explore the results in more detail, we focus on the DataFrame for cluster 9. The interpretation of Top is the same as described for t-tests, and Simes’ method is again used to combine p-values across pairwise comparisons. If we want more focused sets, we can also change pval.type= as previously described.

interesting.wmw <- markers.wmw[[chosen]]
interesting.wmw[1:10,1:4]

DataFrame with 10 rows and 4 columns
                      Top      p.value          FDR summary.AUC
                <integer>    <numeric>    <numeric>   <numeric>
ENSG00000177954         1 9.93712e-116 2.60807e-112    0.982576
ENSG00000019582         1 6.21484e-122 1.90298e-118    0.997488
ENSG00000204287         1 1.45641e-144 2.67572e-140    0.994288
ENSG00000156738         1 1.21982e-130 7.47016e-127    0.914772
ENSG00000167526         1 1.24359e-100  1.75748e-97    0.972138
ENSG00000169442         2  2.38680e-81  1.21806e-78    0.981779
ENSG00000109475         2  4.76372e-83  2.65209e-80    0.928503
ENSG00000196126         2 6.11075e-138 5.61333e-134    0.959484
ENSG00000205542         2  2.43771e-37  3.29306e-35    0.978702
ENSG00000229117         2  1.07000e-96  1.31053e-93    0.940763

The DataFrame contains the AUCs from comparing cluster 9 to every other cluster. A value greater than 0.5 indicates that the gene is upregulated in the current cluster compared to the other cluster, while values less than 0.5 correspond to downregulation. We would typically expect AUCs of 0.7-0.8 for a strongly upregulated candidate marker.

best.set <- interesting.wmw[interesting.wmw$Top <= 5,]
AUCs <- getMarkerEffects(best.set, prefix="AUC")
AUCs.ens <- rownames(AUCs)
rownames(AUCs) <- rowData(sce)[rownames(AUCs), "Symbol"]


library(pheatmap)
pheatmap(AUCs, breaks=seq(0, 1, length.out=21),
    color=viridis::viridis(21))

One practical advantage of the WMW test over the Welch t-test is that it is symmetric with respect to differences in the size of the groups being compared. This means that, all else being equal, the top-ranked genes on each side of a DE comparison will have similar expression profiles regardless of the number of cells in each group. In contrast, the t-test will favor genes where the larger group has the higher relative variance as this increases the estimated degrees of freedom and decreases the resulting p-value. This can lead to unappealing rankings when the aim is to identify genes upregulated in smaller groups. The WMW test is not completely immune to variance effects - for example, it will slightly favor detection of DEGs at low average abundance where the greater number of ties at zero deflates the approximate variance of the rank sum statistic - but this is relatively benign as the selected genes are still fairly interesting.

marker.t <- findMarkers(sce, groups=sce$source_name, 
    direction="up", restrict=c("PBMMC", "ETV6-RUNX1"))
marker.w <- findMarkers(sce, groups=sce$source_name, 
    direction="up", restrict=c("PBMMC", "ETV6-RUNX1"), test.type="wilcox")

# Upregulated in type 1:
type1 <- "PBMMC"
marker.type1.t <- marker.t[[type1]]
marker.type1.w <- marker.w[[type1]]
chosen.type1.t <- rownames(marker.type1.t)[1:30]
chosen.type1.w <- rownames(marker.type1.w)[1:30]
u.type1.t <- setdiff(chosen.type1.t, chosen.type1.w)
u.type1.w <- setdiff(chosen.type1.w, chosen.type1.t)

# Upregulated in gamma:
type2 <- "ETV6-RUNX1"
marker.type2.t <- marker.t[[type2]]
marker.type2.w <- marker.w[[type2]]
chosen.type2.t <- rownames(marker.type2.t)[1:30]
chosen.type2.w <- rownames(marker.type2.w)[1:30]
u.type2.t <- setdiff(chosen.type2.t, chosen.type2.w)
u.type2.w <- setdiff(chosen.type2.w, chosen.type2.t)

# Examining all uniquely detected markers in each direction.
library(scater)
subset <- sce[,sce$source_name %in% c(type1, type2)]
gridExtra::grid.arrange(
    plotExpression(subset, x="source_name", features=u.type1.t, ncol=2) +
        ggtitle(sprintf("Upregulated in %s, t-test-only", type1)),
    plotExpression(subset, x="source_name", features=u.type1.w, ncol=2) +
        ggtitle(sprintf("Upregulated in %s, WMW-test-only", type1)),
    plotExpression(subset, x="source_name", features=u.type2.t, ncol=2) +
        ggtitle(sprintf("Upregulated in %s, t-test-only", type2)),
    plotExpression(subset, x="source_name", features=u.type2.w, ncol=2) +
        ggtitle(sprintf("Upregulated in %s, WMW-test-only", type2)),
    ncol=2
)

The main disadvantage of the WMW test is that the AUCs are much slower to compute compared to t-statistics. This may be inconvenient for interactive analyses involving multiple iterations of marker detection. We can mitigate this to some extent by parallelizing these calculations using the BPPARAM= argument in findMarkers().

1.1.5 Using a binomial test

The binomial test identifies genes that differ in the proportion of expressing cells between clusters. (For the purposes of this section, a cell is considered to express a gene simply if it has non-zero expression for that gene.) This represents a much more stringent definition of marker genes compared to the other methods, as differences in expression between clusters are effectively ignored if both distributions of expression values are not near zero. The premise is that genes are more likely to contribute to important biological decisions if they were active in one cluster and silent in another, compared to more subtle “tuning” effects from changing the expression of an active gene. From a practical perspective, a binary measure of presence/absence is easier to validate.

We perform pairwise binomial tests between clusters using the findMarkers() function with test=“binom”. This returns a list of DataFrames containing marker statistics for each cluster such as the Top rank and its p

-value. Here, the effect size is reported as the log-fold change in this proportion between each pair of clusters. Large positive log-fold changes indicate that the gene is more frequently expressed in one cluster compared to the other. We focus on genes that are upregulated in each cluster compared to the others by setting direction=“up”.

markers.binom <- findMarkers(sce, test="binom", direction="up", groups=sce$clusterStg)
print(names(markers.binom))

 [1] "c1"  "c2"  "c3"  "c4"  "c5"  "c6"  "c7"  "c8"  "c9"  "c10" "c11" "c12"
[13] "c13" "c14" "c15" "c16" "c17" "c18" "c19" "c20" "c21" "c22"

interesting.binom <- markers.binom[[chosen]]
print(colnames(interesting.binom))

 [1] "Top"           "p.value"       "FDR"           "summary.logFC"
 [5] "logFC.c1"      "logFC.c2"      "logFC.c3"      "logFC.c4"     
 [9] "logFC.c5"      "logFC.c6"      "logFC.c7"      "logFC.c8"     
[13] "logFC.c10"     "logFC.c11"     "logFC.c12"     "logFC.c13"    
[17] "logFC.c14"     "logFC.c15"     "logFC.c16"     "logFC.c17"    
[21] "logFC.c18"     "logFC.c19"     "logFC.c20"     "logFC.c21"    
[25] "logFC.c22"    

The plot below confirms that the top genes exhibit strong differences in the proportion of expressing cells in cluster 9 compared to the others.

library(scater)
top.genes <- head(rownames(interesting.binom))
#plotExpression(sce, x="clusterStg", features=top.genes)
plotExpression(sce, x="clusterStg", features=top.genes[1])

plotExpression(sce, x="clusterStg", features=top.genes[2])

plotExpression(sce, x="clusterStg", features=top.genes[3])

plotExpression(sce, x="clusterStg", features=top.genes[4])

1.1.6 Combining multiple marker statistics

On occasion, we might want to combine marker statistics from several testing regimes into a single DataFrame. This allows us to easily inspect multiple statistics at once to verify that a particular gene is a strong candidate marker. For example, a large AUC from the WMW test indicates that the expression distributions are well-separated between clusters, while the log-fold change reported with the t-test provides a more interpretable measure of the magnitude of the change in expression. We use the multiMarkerStats() to merge the results of separate findMarkers() calls into one DataFrame per cluster, with statistics interleaved to facilitate a direct comparison between different test regimes.

combined <- multiMarkerStats(
    t=findMarkers(sce, groups=sce$clusterStg, direction="up"),
    wilcox=findMarkers(sce, groups=sce$clusterStg, test="wilcox", direction="up"),
    binom=findMarkers(sce, groups=sce$clusterStg, test="binom", direction="up")
)

# Interleaved marker statistics from both tests for each cluster.
print(colnames(combined[["c1"]]))

 [1] "Top"                 "p.value"             "FDR"                
 [4] "t.Top"               "wilcox.Top"          "binom.Top"          
 [7] "t.p.value"           "wilcox.p.value"      "binom.p.value"      
[10] "t.FDR"               "wilcox.FDR"          "binom.FDR"          
[13] "t.summary.logFC"     "wilcox.summary.AUC"  "binom.summary.logFC"
[16] "t.logFC.c2"          "wilcox.AUC.c2"       "binom.logFC.c2"     
[19] "t.logFC.c3"          "wilcox.AUC.c3"       "binom.logFC.c3"     
[22] "t.logFC.c4"          "wilcox.AUC.c4"       "binom.logFC.c4"     
[25] "t.logFC.c5"          "wilcox.AUC.c5"       "binom.logFC.c5"     
[28] "t.logFC.c6"          "wilcox.AUC.c6"       "binom.logFC.c6"     
[31] "t.logFC.c7"          "wilcox.AUC.c7"       "binom.logFC.c7"     
[34] "t.logFC.c8"          "wilcox.AUC.c8"       "binom.logFC.c8"     
[37] "t.logFC.c9"          "wilcox.AUC.c9"       "binom.logFC.c9"     
[40] "t.logFC.c10"         "wilcox.AUC.c10"      "binom.logFC.c10"    
[43] "t.logFC.c11"         "wilcox.AUC.c11"      "binom.logFC.c11"    
[46] "t.logFC.c12"         "wilcox.AUC.c12"      "binom.logFC.c12"    
[49] "t.logFC.c13"         "wilcox.AUC.c13"      "binom.logFC.c13"    
[52] "t.logFC.c14"         "wilcox.AUC.c14"      "binom.logFC.c14"    
[55] "t.logFC.c15"         "wilcox.AUC.c15"      "binom.logFC.c15"    
[58] "t.logFC.c16"         "wilcox.AUC.c16"      "binom.logFC.c16"    
[61] "t.logFC.c17"         "wilcox.AUC.c17"      "binom.logFC.c17"    
[64] "t.logFC.c18"         "wilcox.AUC.c18"      "binom.logFC.c18"    
[67] "t.logFC.c19"         "wilcox.AUC.c19"      "binom.logFC.c19"    
[70] "t.logFC.c20"         "wilcox.AUC.c20"      "binom.logFC.c20"    
[73] "t.logFC.c21"         "wilcox.AUC.c21"      "binom.logFC.c21"    
[76] "t.logFC.c22"         "wilcox.AUC.c22"      "binom.logFC.c22"    

#head(combined[["c1"]][,1:9])
combined[["c1"]]$Symbol <- rowData(sce)[rownames(combined[["c1"]]), "Symbol"]
tmpCol <- c("Symbol", colnames(combined[["c1"]])[1:9])
head(combined[["c1"]][,tmpCol])

DataFrame with 6 rows and 10 columns
                     Symbol       Top      p.value          FDR     t.Top
                <character> <integer>    <numeric>    <numeric> <integer>
ENSG00000019582        CD74         1  2.12393e-92  1.44522e-89         1
ENSG00000204287     HLA-DRA         1  1.16719e-94  8.93483e-92         1
ENSG00000100721       TCL1A         1 9.16292e-118 1.87046e-114         1
ENSG00000169575      VPREB1         1 4.16629e-130 1.53086e-126         1
ENSG00000188643     S100A16         2 5.92877e-144 5.44617e-140         1
ENSG00000244734         HBB         2  5.61294e-32  2.31732e-30         1
                wilcox.Top binom.Top    t.p.value wilcox.p.value binom.p.value
                 <integer> <integer>    <numeric>      <numeric>     <numeric>
ENSG00000019582          1         1  0.00000e+00   3.10099e-197   2.12393e-92
ENSG00000204287          1         1  0.00000e+00   6.54248e-209   1.16719e-94
ENSG00000100721          1         1  0.00000e+00   1.03986e-207  9.16292e-118
ENSG00000169575          1         1  0.00000e+00   3.90854e-206  4.16629e-130
ENSG00000188643          2         1 4.80006e-203   2.81709e-167  5.92877e-144
ENSG00000244734          2         2  4.81933e-86    1.76203e-36   5.61294e-32

In addition, multiMarkerStats() will compute a number of new statistics by combining the per-regime statistics. The combined Top value is obtained by simply taking the largest Top value across all tests for a given gene, while the reported p.value is obtained by taking the largest p-value. Ranking on either metric focuses on genes with robust differences that are highly ranked and detected by each of the individual testing regimes. Of course, this might be considered an overly conservative approach in practice, so it is entirely permissible to re-rank the DataFrame according to the Top or p.value for an individual regime (effectively limiting the use of the other regimes’ statistics to diagnostics only).

Write list to file:

tmpFn <- sprintf("%s/%s/Robjects/%s_sce_nz_postDeconv%s_clustMarkCombi.Rds", projDir, outDirBit, setName, setSuf)
print(tmpFn)

[1] "/mnt/scratcha/bioinformatics/baller01/20200511_FernandesM_ME_crukBiSs2020/AnaWiSce/Attempt1/Robjects/caron_sce_nz_postDeconv_5hCellPerSpl_clustMarkCombi.Rds"

saveRDS(combined, file=tmpFn)

1.2  Invalidity of p-values

1.2.1 11.5.1 From data snooping

All of our DE strategies for detecting marker genes between clusters are statistically flawed to some extent. The DE analysis is performed on the same data used to obtain the clusters, which represents “data dredging” (also known as fishing or data snooping). The hypothesis of interest - are there differences between clusters? - is formulated from the data, so we are more likely to get a positive result when we re-use the data set to test that hypothesis.

The practical effect of data dredging is best illustrated with a simple simulation. We simulate i.i.d. normal values, perform k-means clustering and test for DE between clusters of cells with findMarkers(). The resulting distribution of p-values is heavily skewed towards low values. Thus, we can detect “significant” differences between clusters even in the absence of any real substructure in the data. This effect arises from the fact that clustering, by definition, yields groups of cells that are separated in expression space. Testing for DE genes between clusters will inevitably yield some significant results as that is how the clusters were defined.

Distribution of \(p\)-values from a DE analysis between two clusters in a simulation with no true subpopulation structure:

library(scran)
set.seed(0)
y <- matrix(rnorm(100000), ncol=200)
clusters <- kmeans(t(y), centers=2)$cluster
out <- findMarkers(y, clusters)
hist(out[[1]]$p.value, col="grey80", xlab="p-value")

For marker gene detection, this effect is largely harmless as the p-values are used only for ranking. However, it becomes an issue when the p-values are used to define “significant differences” between clusters with respect to an error rate threshold. Meaningful interpretation of error rates require consideration of the long-run behavior, i.e., the rate of incorrect rejections if the experiment were repeated many times. The concept of statistical significance for differences between clusters is not applicable if clusters and their interpretations are not stably reproducible across (hypothetical) replicate experiments.

1.2.2  Nature of replication

The naive application of DE analysis methods will treat counts from the same cluster of cells as replicate observations. This is not the most relevant level of replication when cells are derived from the same biological sample (i.e., cell culture, animal or patient). DE analyses that treat cells as replicates fail to properly model the sample-to-sample variability (Lun and Marioni 2017). The latter is arguably the more important level of replication as different samples will necessarily be generated if the experiment is to be replicated. Indeed, the use of cells as replicates only masks the fact that the sample size is actually one in an experiment involving a single biological sample. This reinforces the inappropriateness of using the marker gene p-values to perform statistical inference.

“We strongly recommend selecting some markers for use in validation studies with an independent replicate population of cells. A typical strategy is to identify a corresponding subset of cells that express the upregulated markers and do not express the downregulated markers. Ideally, a different technique for quantifying expression would also be used during validation, e.g., fluorescent in situ hybridisation or quantitative PCR. This confirms that the subpopulation genuinely exists and is not an artifact of the scRNA-seq protocol or the computational analysis.”

See the OSCA chapter on Marker gene detection

Challenge Identify markers for a different cluster and try to identify the cell type.