qc.Rmd

---
title: "QC of single cell data"
author: "Sydney Precision Bioinformatics Group"
date: "04 December 2019"
output:
  html_document:
    css: style.css
    code_folding: show
    fig_height: 10
    fig_width: 10
    toc: yes
    number_sections: true
    toc_depth: 3
    toc_float: yes
---

# Introduction

One of the first steps in single-cell RNA-seq analysis is to perform a quality check of the data we have sequenced in an experiment. This process is a necessary step before we proceed to data merging and other downstream data analyses. There are multiple quality control (QC) tools for sequencing data that we can use for both bulk and single-cell RNA-seq data. Here, we will look at some QC measures and generate a report on the quality of our single-cell data.

## Load in packages

Load in packages required to generate the QC plots.

```{r load qc pkg, warning=FALSE, message=FALSE}
library(DropletUtils)
library(dplyr)
library(ggplot2)
library(scales)
library(ggpubr)
library(stringr)
library(forcats)
library(SingleCellExperiment)

theme_set(theme_classic(16))
```

# The mouse liver dataset (Su et al. 2017)

For demonstration purpose, we will perform some QC checks on a single-cell mouse liver dataset generated by Su et al. (2017). This liver dataset contains 507 cells at seven developmental stages between embryonic day 11.5 and postnatal day 2.5. The cells were sequenced using the Fluidigm C1 platform and the reads are paired end. The file containing the raw count matrix is located in the `data` folder in the zip file you have downloaded as `GSE87795_counts.csv`.

The dataset of read counts is stored in a `csv` file (which is simply an Excel file). We can load the data into R using the `read.csv()` function. Here, we will use two options:

  + `header = TRUE` means that we will take the first row of the file as the column names.
  + `row.names = 1` means that we will take the first column of the file as the row names.


```{r}
datapath = "./data/"
# datapath = "/home/data/"
liver = read.csv(file = paste0(datapath, "GSE87795_counts.csv"),
                 header = TRUE,
                 row.names = 1)
```

First, let's check some basics of this data. This dataset contains 51918 rows and 507 columns. Each row corresponds to a gene, and each column corresponds to a cell. The first five rows and five columns of the dataset can be shown using `liver[1:5, 1:5]`.

```{r}
dim(liver)
liver[1:5, 1:5]
```


Notice that the column names of this data contain information regarding the developmental stage of each cell when it was sequenced duing the experiment. We can clean up this information using a short line of code.

With some code, we can produce a frequency table of the number of cells at each developmental stage.

```{r}
stage = str_split(colnames(liver), "_") %>%
  sapply("[[", 1)

table(stage)
```




## A quick pre-processing of the matrix

A simple way to make single-cell datasets manageable for analysis is to remove genes that are not expressed in any cells (i.e. genes that have zero counts across all cells). Since these genes do not contain any measurable signals, removing them does not affect the downstream analysis. In addition, a smaller dataset actually improves the performance and speed of most algorithms.

```{r}
liver = liver[rowSums(liver) != 0, ]
dim(liver)
```


Here, we first used `rowSums(liver)` to calculate the total counts for each gene across all cells. `rowSums(liver) != 0` then looks for those genes that are *not equal* (represented by the `!=` symbol) zeroes. And then we will subset only those rows in the data (remember, rows of the data are genes) that are not purely zeroes using the square brackets `[,]`.

<h4> <span class="glyphicon glyphicon-education" aria-hidden="true"></span> Quick quiz </h4>

1. How many genes do we have now?
1. Did we lose any cells (columns)?





# Waterfall plot

The waterfall plot shows the log-count against the log-rank of each barcode. The barcodes are ranked based on the number of counts each barcode has. We will compute these statistics using the `barcodeRanks` function from the `DropletUtils` package. And then we will extract the statistics and perform some plotting using the `ggplot2` package.

The waterfall plots show two points of interest:

  + **inflection point** is the point on the curve where the first derivative is minimised
  + **knee point** is the point where the signed curvature is minimised


These two points potentially indicate empty droplets (background barcodes). To elaborate on this concept, it is helpful to know that for microfluidic sequencing systems each droplet contains a bead and each bead is marked uniquely with a barcode. Intuitively, we can consider a droplet and the barcode associated with it as being a unique cell. Due to limitation of the technology, sometimes the droplet doesn't always contain a cell, or it might contain multiples cells (which we will not explore here). We can now use the waterfall plot to detect those empty droplets. The barcodes with the number of counts dropping rapidly beyond the knee point may suggest empty reads.



```{r, fig.height=6, fig.width=8}
barcode = DropletUtils::barcodeRanks(liver)
barcode_data = as.data.frame(barcode)

barcode_points = data.frame(
  type = c("inflection", "knee"),
  value = c(barcode@metadata$inflection, barcode@metadata$knee))


ggplot(data = barcode_data, aes(x = rank, y = total)) +
  geom_point() +
  geom_hline(data = barcode_points,
             aes(yintercept = value,
                 colour = type), linetype = 2) +
  scale_x_log10(breaks = trans_breaks("log10", function(x) 10^x),
                labels = trans_format("log10", math_format(10^.x))) +
  scale_y_log10(breaks = trans_breaks("log10", function(x) 10^x),
                labels = trans_format("log10", math_format(10^.x))) +
  labs(title = "Waterfall plot of read counts (log)",
       x = "Log Rank",
       y = "Log counts")
```



<h4> <span class="glyphicon glyphicon-education" aria-hidden="true"></span> Quick quiz </h4>

Based on the waterfall plot, should we filter out any cells?

# Library size

Next, we can evaluate the total library size, that is the total read counts per cell. The library size may be used as an indicator of sequencing depth and it is a typical diagnostic for RNA-Seq data (for both the bulk and singl-cell type).

A small library size for a cell may suggest abnormalies and lack of captured signals. In such a case, the cells with unusually small library size should be removed.

We will begin with a `data.frame` capturing the characteristics of each column `cell_plotdf_full` and then perform plotting with `ggplot2`. This `cell_plotdf_full` has four columns:

1. `cell_name`: The name of each cell in the data
1. `stage`: The developmental time stage of each cell.
1. `library`: The total number of read counts for a cell.
1. `num_genes`: The total number of non-zero measurements for a cell. i.e. Number of genes expressed in a cell.


For the ease of visualisation, we will subset our `cell_plotdf_full` to only those cells in the 11.5, 12.5 and 13.5 development stage. Then, we will also calculate the mean of `library` and `num_genes` for *each* development stage.



```{r, fig.height=6, fig.width=8}
cell_plotdf_full = data.frame(
  cell_name = colnames(liver),
  stage = stage,
  library = colSums(liver),
  num_genes = colSums(liver != 0)
)

head(cell_plotdf_full)

cell_plotdf = cell_plotdf_full %>%
 dplyr::filter(stage %in% c("E11.5", "E12.5", "E13.5"))


cell_plotdf_mean = cell_plotdf %>%
  dplyr::group_by(stage) %>%
  dplyr::summarise(library_mean = mean(library),
                   num_genes_mean = mean(num_genes))


cell_plotdf_mean

ggplot(cell_plotdf,
       aes(x = library/1e6, colour = stage, fill = stage)) +
  geom_histogram(aes(y = ..density..), colour = "black", fill = "white", bins = 30) +
  geom_density(alpha = 0.4) +
  geom_vline(data = cell_plotdf_mean,
             aes(xintercept = library_mean/1e6),
             colour = "blue", linetype = "dashed", size = 1.5) +
  # scale_x_continuous(breaks = 1:6) +
  labs(title = "Histogram of library size",
       x = "library size (in millions)",
       y = "density") +
  facet_wrap(~ stage, ncol = 3)
```




When we compare the distribution of library size for all cells at three developmental stages, we can see that mean library size is lowest in E11.5 cells (mean = 2.838642 million) and highest in E13.5 cells (mean = 3.639297 million). E12.5 cells (mean = 3.094663 million) have a greater distribution of cells with few reads. Cells with very few reads are likely to have failed to capture the transcriptome of a cell, and thus should be filtered.



# Number of uniquely expressed genes

Closely related to the idea of a library size is the number of uniquely expressed genes.
In addition to sequencing depth, we can also filter for cells that express genes with sufficiently good coverage of the transcriptome. This is to ensure that we capture cells with not only sufficient reads but also those that have reads that are relatively evenly distributed across the transcriptome. To this end, we measure the number of uniquely expressed genes in individual cells.

```{r, fig.height=6, fig.width=8, message=FALSE, warning=FALSE}
ggplot(cell_plotdf,
       aes(x = num_genes/1000, colour = stage, fill = stage)) +
  geom_histogram(aes(y = ..density..), colour = "black", fill = "white", bins = 30) +
  geom_density(alpha = 0.4) +
  geom_vline(data = cell_plotdf_mean,
             aes(xintercept = num_genes_mean/1000),
             colour = "blue", linetype = "dashed", size = 1.5) +
  labs(title = "Histogram of number of unique genes expressed",
       x = "Number of unique genes expressed (in thousands)",
       y = "density") +
  facet_wrap(~ stage, ncol = 3)
```

From above, we can see that most cells have between 9000-10,000 detected genes. This is roughly the amount expected for high-depth scRNA-seq; however, the expected value may vary between experimental protocols, sequencing depth, and cell type used. Note that, unlike library size, a "heavy tail" can be observed for each of the distribution, indicating an unequal detection of genes across populations. Filtering out cells with low number of genes identified is another way to improve the quality of the dataset.

# Mitochondrial gene expression

Unhealthy or dying cells are associated with high expression of mitochondrial genes, and thus the proportion of mitochondrial expression is commonly used as a quality metric.

```{r, message = FALSE, warning = FALSE}
mito_genes <- c("ENSMUSG00000064336", "ENSMUSG00000064337", "ENSMUSG00000064338", "ENSMUSG00000064339", "ENSMUSG00000064340", "ENSMUSG00000064341", "ENSMUSG00000064342", "ENSMUSG00000064343", "ENSMUSG00000064344", "ENSMUSG00000064345", "ENSMUSG00000064346", "ENSMUSG00000064347", "ENSMUSG00000064348", "ENSMUSG00000064349", "ENSMUSG00000064350", "ENSMUSG00000064351", "ENSMUSG00000064352", "ENSMUSG00000064353", "ENSMUSG00000064354", "ENSMUSG00000064355", "ENSMUSG00000064356", "ENSMUSG00000064357", "ENSMUSG00000064358", "ENSMUSG00000064359", "ENSMUSG00000064360", "ENSMUSG00000064361", "ENSMUSG00000065947", "ENSMUSG00000064363", "ENSMUSG00000064364", "ENSMUSG00000064365", "ENSMUSG00000064366", "ENSMUSG00000064367", "ENSMUSG00000064368", "ENSMUSG00000064369", "ENSMUSG00000064370", "ENSMUSG00000064371", "ENSMUSG00000064372", "ENSMUSG00000096105")


liver_mito = liver[mito_genes, ]

mito_plotdf_full = data.frame(
  cell_name = colnames(liver_mito),
  stage = stage,
  mito_percent = colSums(liver_mito)/colSums(liver) * 100
)

mito_plotdf = mito_plotdf_full %>% dplyr::filter(stage %in% c("E11.5", "E12.5", "E13.5"))

mito_plotdf_mean = mito_plotdf %>%
  group_by(stage) %>%
  dplyr::summarise(mito_percent_mean = mean(mito_percent))

mito_plotdf_mean

g1 = ggplot(mito_plotdf,
       aes(x = mito_percent, colour = stage, fill = stage)) +
  geom_histogram(aes(y = ..density..), colour = "black", fill = "white", bins = 30) +
  geom_density(alpha = 0.4) +
  geom_vline(data = mito_plotdf_mean,
             aes(xintercept = mito_percent_mean),
             colour = "blue", linetype = "dashed", size = 1.5) +
  labs(x = "Mitochondrial percentage",
       y = "Density") +
  facet_wrap(~ stage, nrow = 3)


g2 = ggplot(mito_plotdf,
            aes(x = stage,
                y = mito_percent,
                fill = stage)) +
  geom_violin(alpha = 0.8, draw_quantiles = 0.5) +
  labs(y = "Mitochondrial percentage")

ggpubr::ggarrange(g1, g2, ncol=  2, common.legend = TRUE, legend = "right")
```


Filtering cells with greater than 10% of mitochondrial gene expression would remove cells that have abnormally high mitochondrial gene expression.

<h4> <span class="glyphicon glyphicon-education" aria-hidden="true"></span> Quick quiz </h4>

1. Based on the figure above, would you say most cells are unhealthy or healthy?
2. How does the percentage of mitochondrial gene expression change over time point?


# Contribution by the top 50 expressed genes

As well as the above metrics, it is often instructive to know the proportions of the reads consumed by the top 50 expressed genes. Moreover, an assessment of the top gene list before and after filtering can give you a good indication of whether the filtering has been successful, as the top genes are often occupied by mitochondrial genes from unhealthy cells.


We can calculate the top 50 expressed genes for all cells as shown below.

```{r,  fig.height=12, fig.width=8}
## We will go to each column of the matrix, and divide the column by its sum
## and convert this into a percentage.
## The result is a matrix with each column that sum up to 100%.
liver_percent = sweep(liver, 2, STATS = colSums(liver), FUN = "/") * 100

liver_percent %>% colSums %>% head

liver_percent_meanEachGene = rowMeans(liver_percent)

top50Genes = liver_percent_meanEachGene %>% sort(decreasing = TRUE) %>% head(50) %>% names


liver_percent_plotdf = liver_percent[top50Genes, ] %>%
  as.matrix %>%
  reshape2::melt(varnames = c("gene_name", "cell_name"),
                 value.name = "percent") %>%
  dplyr::mutate(gene_name = fct_reorder(gene_name, percent, mean))

liver_percent_plotdf %>%
  ggplot(aes(x = gene_name, y = percent, colour = gene_name)) +
  geom_violin() +
  geom_jitter(alpha = 0.2) +
  coord_flip() +
  labs(x = "genes",
       y = "percent of library",
       title = "Top 50 genes for all cell stages") +
  theme(legend.position = "none")
```



The plot above shows that there are three spike-in genes within the top 50 genes. The presence of multiple spike-ins at the top of the list may suggest that the concentration of spike-ins added to the experiment may need to be optimized for subsequent experiments. Overall, the relatively flat distributions indicate good coverage of the transcriptome.



Alternatively, we can calculate the top 50 expressed genes for each stage, but we will leave this as an extension exercise.


```{r,  fig.height=12, fig.width=8}
liver_percent_split = split.data.frame(x = t(liver_percent), f = stage) %>% lapply(t)

liver_percent_split_meanEachGene = liver_percent_split %>% lapply(rowMeans)


e115_top50 = liver_percent_split_meanEachGene$E11.5 %>% sort(decreasing = TRUE) %>% head(50) %>% names

e115_liver_percent = liver_percent[e115_top50, stage == "E11.5"] %>%
  as.matrix %>%
  reshape2::melt(varnames = c("gene_name", "cell_name"),
                 value.name = "percent") %>%
  dplyr::mutate(gene_name = fct_reorder(gene_name, percent, mean))

e115_liver_percent %>%
  ggplot(aes(x = gene_name, y = percent, colour = gene_name)) +
  geom_violin() +
  geom_jitter(alpha = 0.2) +
  coord_flip() +
  labs(x = "genes",
       y = "percent of library",
       title = "Top 50 expressed genes for E11.5 cells") +
  theme(legend.position = "none")
```



<!-- ```{r} -->
<!-- liver.percent.melt = liver.percent %>% -->
<!--   as.matrix %>%  -->
<!--   reshape2::melt(varnames = c("gene.name", "cell.name"), value.name = "percent") %>%  -->
<!--   dplyr::mutate(stage = str_split(cell.name, "_") %>% sapply("[[", 1))  -->

<!-- e115.percent = liver.percent.melt %>%  -->
<!--   dplyr::filter(stage == "E11.5") %>%  -->
<!--   dplyr::mutate(gene.name = forcats::fct_reorder(gene.name, percent, mean, .desc = TRUE)) %>%  -->
<!--   dplyr::filter(as.integer(gene.name) <= 50) -->


<!-- e115.percent %>%  -->
<!--   ggplot(aes(x = fct_rev(gene.name), y = percent, colour = gene.name)) + -->
<!--   geom_violin() + -->
<!--   geom_jitter(alpha = 0.2) + -->
<!--   coord_flip() + -->
<!--   labs(x = "percent",  -->
<!--        y = "") -->
<!--   theme(legend.position = "none") -->
<!-- ``` -->



# Cell Filtering

Here we will give a demonstration of how you can remove some of the cells from the dataset.

In the `DropletUtils` package, it was recommended that we should only retain those cells above the inflection point in the plot. So we will perform subsetiting on the cells (columns) of this data.

```{r}
table(barcode_data$total >= barcode@metadata$inflection)
liver_new = liver[, barcode_data$total >= barcode@metadata$inflection]
dim(liver_new)
```


Recall in the mitochondiral plot, we have identified cells with a mitochondrial percentage greater than 10% as cells that we should remove. Thus, we will subset the data on the columns to retain those cells with a mitochondrial percentage less or equal to 10%.

```{r}
dim(liver)
liver_new2 = liver[, barcode_data$total >= barcode@metadata$inflection & mito_plotdf_full$mito_percent <= 10]
dim(liver_new2)
```

# Gene Filtering


Lastly in the quality control step, we will perform quality control from the gene persepctive, to filter out the genes that are expressed in less than 10 cells.


Finally we should check that the dimension of the dataset is correct after filtering.

```{r}
num_cell_expressed = rowSums(liver_new2 > 0)
liver_new3 = liver_new2[num_cell_expressed > 10, ]
dim(liver_new3)
```


# Create `SingleCellExperiment`

For further single cell data analysis, we will create a `SingleCellExperiment` object using our filtered count matrix.


```{r}
liver_sce = SingleCellExperiment::SingleCellExperiment(assay = list(counts = as.matrix(liver_new3)))
liver_sce
```



Noted that this experiment contains the spike-ins transcripts. We can use `splitAltExp` to remove the corresponding rows (Genes start with "ERCC-") from the main `SingleCellExperiment`, and stored these rows into a nested alternative Experiment (`altExp()`).


```{r}
is_spike = grepl("^ERCC-", rownames(liver_sce))
liver_sce = splitAltExps(liver_sce, ifelse(is_spike, "ERCC", "gene"))
altExpNames(liver_sce)
altExp(liver_sce)
```


Lastly in this preprocessing tutorial, we will perform log-transformation for the counts matrix and store the normalised matrix as `logcounts` in Assays.

```{r}
liver_counts <- counts(liver_sce)
lib_sizes <- colSums(liver_counts)
size_factors <- lib_sizes/mean(lib_sizes)
logcounts(liver_sce) <- log2(t(t(liver_counts)/size_factors) + 1)
assayNames(liver_sce)
```

```{r}
liver_sce
```

# Reference

Aaron Lun and Davide Risso (2019). SingleCellExperiment: S4 Classes for Single Cell Data. R package version 1.7.7.

Lun ATL, Riesenfeld S, Andrews T, Dao T, Gomes T,
participants in the 1st Human Cell Atlas Jamboree,
Marioni JC (2019). “EmptyDrops: distinguishing cells
from empty droplets in droplet-based single-cell RNA
sequencing data.” _Genome Biol._, *20*, 63. doi:
10.1186/s13059-019-1662-y (URL:
https://doi.org/10.1186/s13059-019-1662-y).

Griffiths JA, Richard AC, Bach K, Lun ATL, Marioni
JC (2018). “Detection and removal of barcode
swapping in single-cell RNA-seq data.” _Nat.
Commun._, *9*(1), 2667. doi:
10.1038/s41467-018-05083-x (URL:
https://doi.org/10.1038/s41467-018-05083-x).
# Session Info
```{r}
sessionInfo()
```