Introduction to single-cell RNA-seq analysis

Cluster marker genes

Introduction

In order to aid the interpretation of the clustering results, it is helpful to identify genes that are contributing to the separation of cells into those clusters.

The main approach to achieve this, is to identify genes that are differently expressed between clusters. These may be, for example, exclusively expressed in a single cluster or perhaps differentiate between a few different clusters. There are different methods to identify differentially expressed genes, focusing on differences between the mean expression level, the differences in the ranking of genes by expression between clusters or the differences between clusters of the proportions of cells expressing a gene.

Our main objective in this section is to cover some of the methods that can be used to achieve this goal, and obtain a summary table of results.

As always, the OSCA chapter on marker detection contains additional detail and advice.

Setup

Before we start, let’s load our packages and read our data in.

library(scater)
library(scran)
library(pheatmap)
library(tidyverse)

We will load two SingleCellExperiment objects generated in previous sections:

• The SCE with logcounts normalised within each batch generated in the Dataset Integration and Batch Correction section
• The SCE with batch-corrected data integration using the batchelor::fastMNN() algorithm, including clusters covered in the Clustering section.

uncorrected <- readRDS("../Robjects/DataIntegration_all_sce_dimred.Rds")
rownames(uncorrected) <- uniquifyFeatureNames(rownames(uncorrected), rowData(uncorrected)$Symbol)

# data corrected using batch integration with mutual nearest neighbours
corrected <- readRDS("../Robjects/Caron_clustering_material.rds")
rownames(corrected) <- uniquifyFeatureNames(rownames(corrected), rowData(corrected)$Symbol)

Note that we also replace the gene names (rownames) of our objects to use common gene symbols instead of Ensembl IDs. Using the uniquifyFeatureNames() function makes this safely, avoiding duplicated gene symbols.

As a reminder, the corrected object contains the cell cluster assignments obtained in a colData column called “label”:

# visualise cluster assignments on the corrected data
plotTSNE(corrected, 
         colour_by = "label", 
         text_by = "label")

It will also be helpful for data visualisation purposes to copy the cluster assignment in the corrected data to our uncorrected (normalised) SCE object:

# first make sure that the cell names are in the same order
all(colnames(uncorrected) == colnames(corrected))

[1] TRUE

colLabels(uncorrected) <- colLabels(corrected)

Identifying cluster marker genes

Although we have defined our clusters based on the batch-corrected expression values, these should not be used for for gene-based analyses like marker gene detection.

Instead, we should use the uncorrected (normalised) expression values for differential expression between clusters.

The reason for this is that batch correction algorithms do not necessarily preserve the magnitude (or even direction) of the differences in gene expression when attempting to integrate multiple batches. For example, the cosine normalization used by the “Mutual Nearest Neighbours” algorithm (fastMNN()) shrinks the magnitude of the expression values so that the computed log-fold changes have no obvious interpretation. Of greater concern is the possibility that the correction introduces artificial agreement across batches.

In our case, the valid assay to use in differential analysis tests is the logcounts obtained from the deconvolution method (using computePooledFactors & logNormCounts).

Pairwise comparisons for scoring potential marker genes

The basic approach for marker gene identification across clusters is to perform statistical tests for each gene between every pair of clusters based on gene expression. We use scoreMarkers() function to do this. In addition to computing various scores, the function can take account of known batch effects. The function outputs a list of DataFrames, with one DataFrame of scores for each cluster. However, note that the blocking assumes that each pair of clusters is present in at least one of the blocks. If there are two clusters which are not both present in at least one block (in our cases Samples), then that pairwise comparison will by necessity be omitted.

markers <- scoreMarkers(uncorrected, 
                        groups = factor(uncorrected$label), 
                        block =uncorrected$SampleName)

head(markers[["8"]], n=2)

DataFrame with 2 rows and 19 columns
           self.average other.average self.detected other.detected mean.logFC.cohen min.logFC.cohen median.logFC.cohen
              <numeric>     <numeric>     <numeric>      <numeric>        <numeric>       <numeric>          <numeric>
AL627309.5 -0.000708272  -0.000842074     0.0000000    2.14173e-05        0.0419331      -0.0412962          0.0271537
LINC01409   0.017808816   0.025250072     0.0311542    2.70209e-02        0.0251980      -0.1367201          0.0221337
           max.logFC.cohen rank.logFC.cohen  mean.AUC   min.AUC median.AUC   max.AUC  rank.AUC mean.logFC.detected
                 <numeric>        <integer> <numeric> <numeric>  <numeric> <numeric> <integer>           <numeric>
AL627309.5        0.185322             4024  0.503696  0.499629   0.502047  0.513969      4750           0.1162902
LINC01409         0.258749             4376  0.504467  0.486847   0.505058  0.520581      6109           0.0536382
           min.logFC.detected median.logFC.detected max.logFC.detected rank.logFC.detected
                    <numeric>             <numeric>          <numeric>           <integer>
AL627309.5          -0.152179             0.0369594           0.614644                4438
LINC01409           -0.907950             0.0942543           0.538028                4999

This DataFrame contains the results for cluster 8. The first four columns contain summary statistics:

• self.average - the mean log-expression in the cluster of interest
  
• other.average - the grand mean across all other clusters
  
• self.detected - the proportion of cells with detected expression in the cluster of interest
  
• other.detected - the mean detected proportion across all other clusters.

The remaining columns contain summaries of three scores from the pairwise comparisons. The three scores are:

• logFC.cohen - “Cohen’s d” - this is the log fold change of mean gene expression that has been standardized by scaling by the average of the standard deviation across the groups. This can be interpreted in a similar way to log fold change in that a positive value indicates upregulation in the cluster of interest.
  
• AUC - “Area Under the Curve” - this quantifies the ability to distinguish between two gene expression distributions. It can be interpreted as the likelihood that any random cell in the cluster of interest will have a higher expression of the gene than any random cell in the other cluster. It ranges from 0 to 1, where 1 can be interpreted as upregulation, 0 downregulation, and 0.5 as no difference.
• logFC.detected - this is the log fold change in the proportion of cells in which the gene is detected in the cluster of interest versus the proportion of cells in which the gene is detected in the other cluster. Positive values indicate that the gene is detected in a more cells in the cluster of interest than the other cluster. Note this takes no account of the magnitude of the gene expression.

More detail on the differences between these effect size scores can be found in the “Advanced” Marker detection chapter of the OSCA book.

Whilst all the pairwise scores can be retrieved by adding the argument full.stats=TRUE to the scoreMarkers command, by default scoreMarkers returns 5 summary statistics for each score:

• mean.X - this is the mean of the score across all pairwise comparisons. It gives the relative expression of the gene versus the average of the other clusters.
  
• min.X - this is the minimum score obtained across all pairwise comparisons. This is the most stringent summary statistic for detecting upregulated genes, if the score is high, then the gene is upregulated in the cluster of interest relative to all other clusters. Conversely, if the minimum is low, then the gene is downregulated relative to at least one other cluster.
  
• median.X - this is the median of the score across all pairwise comparisons. It is more robust to outliers than the mean. Whilst this is beneficial in avoiding the situation where the effect size is very large in only a small number of comparisons, it may conversely be detrimental to marker gene selection if there are small number of clusters with very similar effect sizes, as these will be effectively ignored.
• max.X - this is the maximum score obtained across all pairwise comparisons. This is the least stringent summary statistic for detecting upregulated genes as a high score only indicates that the gene is upregulated in the cluster of interest relative to at least one other clusters. Conversely, if the maximum is low, then the gene is downregulated relative to all other clusters.
• rank.X - This is the minimum ranking (“min-rank”) of that gene by that score across all clusters. For each pairwise comparison the genes are ranked according to the score, this summary provides the lowest rank for that gene across all comparisons. Essentially, a gene with a low “min-rank” will be one of the best genes according to that score in at least one comparison.

The choice of the summary used for ranking will effect the stringency of the selection. See the the OSCA books “Basic” chapter on Marker gene detection for further discussion of these different summaries. In general the mean and median make reasonable defaults for most applications. In practice, the minimum and maximum are most helpful for diagnosing discrepancies between the mean and median, rather than being used directly for ranking.

Selecting genes based on a given min-rank, say 5, is useful as it will generate a list of genes that is the union of genes with a rank of 5 or less for any pairwise comparison.This will ensure we get at least 5 genes that distinguish the cluster of interest from all other clusters.

For example using the min-rank for Cohen’s d on cluster 11 yields 19 marker genes:

topRanked <- markers[["11"]] %>%
  as.data.frame() %>% 
  select(contains("cohen")) %>% 
  filter(rank.logFC.cohen <= 5) %>% 
  arrange(rank.logFC.cohen)
topRanked

We can then plot these on a heatmap to visualize the gene expression differences:

plotGroupedHeatmap(uncorrected, 
                   features=rownames(topRanked), 
                   group="label", 
                   block="SampleName",
                   center=TRUE, 
                   zlim=c(-3, 3))

We can plot the distribution of expression of specific genes as we have previously, e.g. FCGR3A (CD16), a monocyte marker:

plotTSNE(corrected, 
         by_exprs_values = "reconstructed", 
         colour_by = "FCGR3A",
         text_by = "label")

Adjusting the log-fold change

The AUC and Cohen’s d scores incorporate both the gene expression differences between the clusters and the variance in gene expression scores within each cluster. If a gene with low variance it is possible that it will be ranked highly even if the magnitude of the difference between the clusters is low. These genes will not necessarily make good marker genes. It may therefore be desirable to favour the detection of genes with larger log-fold changes. A log-fold change threshold can be set using the lfc= argument in score markers.

Consider for example cluster 6. The gene RPS18, a ribosomal protein gene, is identified as a good marker for this clusters based on a min-rank for Cohen’s d of 2.

topRanked <- markers[["6"]] %>%
  as.data.frame() %>% 
  select(contains("cohen")) %>% 
  filter(rank.logFC.cohen <= 5) %>% 
  arrange(rank.logFC.cohen)
topRanked

However, if we are interested in finding strong positive markers, a look at the expression shows that it would not really be a good positive marker as fold-change difference between it and other clusters is not great.

plotExpression(uncorrected, 
               features = "RPS18", 
               x = "label")

We can select for strong positive markers by setting an lfc threshold of 2:

markerslfc <- scoreMarkers(uncorrected, 
                        groups = factor(uncorrected$label), 
                        block =uncorrected$SampleName,
                        lfc=2)

RPS18 no longer appears in the candidate cluster marker genes list by min-rank of Cohen’s d.

topRanked <- markerslfc[["6"]] %>%
  as.data.frame() %>% 
  select(contains("cohen")) %>% 
  filter(rank.logFC.cohen <= 5) %>% 
  arrange(rank.logFC.cohen)
topRanked

In fact it’s min-rank for Cohen’s d has dropped to:

markerslfc[["6"]]["RPS18", "rank.logFC.cohen"]

[1] 57

From the OSCA book:

Given that scoreMarkers() already reports effect sizes, it is tempting to take the next step and obtain p-values for the pairwise comparisons. Unfortunately, the p-values from the relevant tests cannot be reliably used to reject the null hypothesis. This is because 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.

More details here

Session information

sessionInfo()

sessionInfo()

R version 4.1.1 (2021-08-10)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Big Sur 11.5.2

Matrix products: default
LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib

locale:
[1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8

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

other attached packages:
 [1] forcats_0.5.1               stringr_1.4.0               dplyr_1.0.8                 purrr_0.3.4                
 [5] readr_2.1.2                 tidyr_1.2.0                 tibble_3.1.6                tidyverse_1.3.1            
 [9] patchwork_1.1.1             pheatmap_1.0.12             igraph_1.2.11               cluster_2.1.2              
[13] bluster_1.2.1               scran_1.20.1                scater_1.20.1               ggplot2_3.3.5              
[17] scuttle_1.2.1               SingleCellExperiment_1.14.1 SummarizedExperiment_1.22.0 Biobase_2.52.0             
[21] GenomicRanges_1.44.0        GenomeInfoDb_1.28.4         IRanges_2.26.0              S4Vectors_0.30.2           
[25] BiocGenerics_0.38.0         MatrixGenerics_1.4.3        matrixStats_0.61.0          DT_0.21                    
[29] knitr_1.37                 

loaded via a namespace (and not attached):
 [1] ggbeeswarm_0.6.0          colorspace_2.0-3          ellipsis_0.3.2            XVector_0.32.0           
 [5] BiocNeighbors_1.10.0      fs_1.5.2                  rstudioapi_0.13           farver_2.1.0             
 [9] fansi_1.0.2               lubridate_1.8.0           xml2_1.3.3                sparseMatrixStats_1.4.2  
[13] jsonlite_1.8.0            broom_0.7.12              dbplyr_2.1.1              compiler_4.1.1           
[17] httr_1.4.2                dqrng_0.3.0               backports_1.4.1           assertthat_0.2.1         
[21] Matrix_1.4-0              fastmap_1.1.0             limma_3.48.3              cli_3.2.0                
[25] BiocSingular_1.8.1        htmltools_0.5.2           tools_4.1.1               rsvd_1.0.5               
[29] gtable_0.3.0              glue_1.6.2                GenomeInfoDbData_1.2.6    Rcpp_1.0.8               
[33] cellranger_1.1.0          jquerylib_0.1.4           vctrs_0.3.8               crosstalk_1.2.0          
[37] DelayedMatrixStats_1.14.3 xfun_0.30                 beachmat_2.8.1            rvest_1.0.2              
[41] lifecycle_1.0.1           irlba_2.3.5               statmod_1.4.36            edgeR_3.34.1             
[45] zlibbioc_1.38.0           scales_1.1.1              hms_1.1.1                 RColorBrewer_1.1-2       
[49] yaml_2.3.5                gridExtra_2.3             sass_0.4.0                stringi_1.7.6            
[53] ScaledMatrix_1.0.0        BiocParallel_1.26.2       rlang_1.0.2               pkgconfig_2.0.3          
[57] bitops_1.0-7              evaluate_0.15             lattice_0.20-45           labeling_0.4.2           
[61] htmlwidgets_1.5.4         cowplot_1.1.1             tidyselect_1.1.2          magrittr_2.0.2           
[65] R6_2.5.1                  generics_0.1.2            metapod_1.0.0             DelayedArray_0.18.0      
[69] DBI_1.1.2                 pillar_1.7.0              haven_2.4.3               withr_2.5.0              
[73] RCurl_1.98-1.6            modelr_0.1.8              crayon_1.5.0              utf8_1.2.2               
[77] rmarkdown_2.12            tzdb_0.2.0                viridis_0.6.2             locfit_1.5-9.5           
[81] grid_4.1.1                readxl_1.3.1              reprex_2.0.1              digest_0.6.29            
[85] munsell_0.5.0             beeswarm_0.4.0            viridisLite_0.4.0         vipor_0.4.5              
[89] bslib_0.3.1