---
title: "Supplementary code for: A BASiCS workflow for expression variability analysis using scRNA-Seq data"
author: 
  - Alan O'Callaghan
  - Nils Eling
  - John C. Marioni
  - Catalina A. Vallejos
bibliography: Workflow.bib
output:
  bookdown::pdf_document2:
    fig_width: 6
    fig_height: 4
---


```{r setup_knitr, include = FALSE, cache = FALSE}
library("BiocStyle")
library("knitr")
knitr::opts_chunk$set(
  message = FALSE, error = FALSE, warning = FALSE,
  cache = 0, cache.path = "cache_supp/",
  fig.path = "figure_supp/"
)
```

# Introduction

This file provides the code required to pre-process the scRNAseq datasets 
analysed in the [F1000Research BASiCS Workflow](https://f1000research.com/articles/11-59#ref33). 
In particular, we provide the code used to perform quality control for the 
presomitic and somitic mesoderm cells extracted from @Ibarra-Soria2018. 

First, we load all necessary libraries:

```{r libs}
library("SingleCellExperiment")
library("scater")
library("scran")
library("BASiCS")
library("ggplot2")
library("ggpointdensity")
library("biomaRt")
library("patchwork")
library("viridis")
library("readxl")
## set default theme for plots
theme_set(theme_bw())
```

# Download and load the data

The full dataset, described by @Ibarra-Soria2018, is stored under the accession 
number E-MTAB-6153 on ArrayExpress and can be obtained from 
[here](https://www.ebi.ac.uk/biostudies/arrayexpress/studies/E-MTAB-6153#). 
This workflows uses the *processed* data which consists of 2 files: 

- `rawCounts.tsv`: matrix of raw expression counts (UMI-based) for all genes
and cells.
- `cellAnnotation.tsv`: cell type labels identified by [@Ibarra-Soria2018].

These will be downloaded into a *downloads* folder. 

First, we download and load the `rawCounts.tsv` file (832.9 MB in size).

```{r download-droplet-data}
options(timeout = 10000)

# Check whether the 'downloads' directory exists. If not, create it. 
if (!dir.exists("downloads")) {
  dir.create("downloads", showWarnings = FALSE)   
}

# Check whether the data was downloaded already
if (!file.exists("downloads/rawCounts.tsv")) {
  website <- "https://www.ebi.ac.uk/"
  folder <- "biostudies/files/E-MTAB-6153/"
  file <- "rawCounts.tsv"
  download.file(
    paste0(website, folder, file),
    destfile = "downloads/rawCounts.tsv"
  )
}
```

```{r read-droplet-data}
# Load the downloaded file (this will take some time!)
raw_counts <- read.delim("downloads/rawCounts.tsv", header = TRUE)
```

These data consist of `r nrow(raw_counts)` genes and `r ncol(raw_counts)` cells,
i.e. this is the raw data prior to the QC steps applied by @Ibarra-Soria2018. 

Prior to using BASiCS, it is important to perform quality control to remove 
poor quality cells and to control for potential heterogeneous substructure 
within the population of cells under study. In this case, [@Ibarra-Soria2018]
has already characterised population structure using clustering techniques. 
Here, we will use such information in order to define the cell populations to
be use as the input for BASiCS. 

The cluster labels identified in the original publication are then loaded (
`cellAnnotation.tsv` file; 883 KB in size):

```{r cluster-labels-droplet}
# Check whether the data was downloaded already
if (!file.exists("downloads/cellAnnotation.tsv")) {
  website <- "https://www.ebi.ac.uk/"
  folder <- "biostudies/files/E-MTAB-6153/"
  file <- "cellAnnotation.tsv"
  download.file(
    paste0(website, folder, file),
    destfile = "downloads/cellAnnotation.tsv"
  )
}
# Load the downloaded file
cluster_labels <- read.table("downloads/cellAnnotation.tsv",
  sep = "\t", header = TRUE, stringsAsFactors = FALSE)
```

The table above contains annotation information for all `r nrow(cluster_labels)` 
cells. 

The data contains information about the three samples from which the cells were 
extracted. This information can help us to explore potential batch effects.

```{r table-batch}
cluster_labels$sample <- as.factor(cluster_labels$sample)
table(cluster_labels$sample)
```

The authors integrated the data from these three samples using the `cellranger`
`aggr` program, which uses downsampling to generate samples with similar depth. 

Furthermore, the data contains the 20 major cell types identified by the authors. 
The number of cells per cluster is displayed below:

```{r table-cluster-labels-droplet}
table(cluster_labels$cellType, useNA = "ifany")
```

Note that `cellType` is missing for `r sum(is.na(cluster_labels$cellType))` cells. 
The latter correspond to cells discarded during QC by [@Ibarra-Soria2018]. Such
cells are labelled with a `badQuality` tag in `cluster_label$cell`. The QC used
by the authors consisted of filtering out cells for which less than 1,000 genes
were detected or that had more than 3\% of reads mapped to mitochondrial genes. 
They also removed potential doublets defined as cells for which expression of 
*Xist* and any of *Kdm5d*, *Eif2s3y*, *Gm29650*, *Uty* or *Ddx3y* (genes located 
in the Y chromosome) was detected. `r sum(!is.na(cluster_labels$cellType))` were
used by the authors in their analysis after applying this exclusion criteria. 

Note that the authors also identified some sub-structure within these main clusters.
For example, they identified two sub-types of presomitic mesoderm cells:

```{r table-subclusters}
# Extract sub-type cell labels (sample sub-index is removed)
cluster_labels$subCellType <- sub("_.*", "", cluster_labels$cell)
with(cluster_labels[cluster_labels$cellType == "presomiticMesoderm", ], 
     table(subCellType, cellType))
```

# Select populations of interest

We selected the somitic 
($n = $ `r sum(cluster_labels$cellType == "somiticMesoderm", na.rm = TRUE)`) and 
pre-somitic mesoderm 
($n = $ `r sum(cluster_labels$cellType == "presomiticMesoderm", na.rm = TRUE)`)
cells. This is an arbitrary choice for illustration purposes only. The steps
of this workflow could be applied to any other pair of cell types. 

```{r droplet-cell-selection}
# Identifies which cells are in the groups of interest
ind_som <- which(
  cluster_labels$cellType == "somiticMesoderm" |
  cluster_labels$cellType == "presomiticMesoderm"
)
# Selects the cells of interest from the raw counts and cluster-label metadata
raw_counts <- raw_counts[, ind_som]
cluster_labels <- cluster_labels[ind_som, ]
```

After selecting these cells, the distribution of the number of cells across
samples is displayed below:

```{r table-batch-aftersel}
table(cluster_labels$cellType, cluster_labels$sample)
```

# Generating a SingleCellExperiment object 

For pre-processing and visualisation purposes, we load the data into a 
`SingleCellExperiment` object. The metadata will be stored in the `colData` slot.

```{r droplet-sce}
droplet_sce <- SingleCellExperiment(
  assays = list(counts = as(as.matrix(raw_counts), "dgCMatrix")),
  colData = cluster_labels
)
# Delete the original table of raw expression counts as no longer needed
rm(raw_counts)
```

# Gene annotation using BioMart

Input data was annotated using Ensembl gene identifiers. 
To facilitate interpretation, it is often useful to obtain a mapping from 
Ensembl gene IDs to gene symbols using the BioMart suite
([http://www.biomart.org](http://www.biomart.org)) via the 
Bioconductor package `r Biocpkg("biomaRt")` [@Durinck2009].
This can also be used to obtain gene-pathways mappings and other
metadata (e.g. gene length), useful for performing functional analysis
of gene sets identified in downstream analyses.

```{r obtain-gene-symbols}
# Check if data storage folder exists and create it if needed
if (!dir.exists("rds")) {
  dir.create("rds", showWarnings = FALSE)
}
if (!file.exists("rds/genenames.Rds")) {
  # Initialize mart and dataset
  ensembl <- useEnsembl(
    biomart = "genes",
    version = 104,
    dataset = "mmusculus_gene_ensembl"
  )
  # Select gene ID and gene name
  genenames <- getBM(
    attributes = c("ensembl_gene_id", "external_gene_name", "gene_biotype", 
                   "chromosome_name", "description"),
    mart = ensembl
  )
  rownames(genenames) <- genenames$ensembl_gene_id
  # Store information for future use
  saveRDS(genenames, "rds/genenames.Rds")
}
# Load pre-downloaded gene annotation information
genenames <- readRDS("rds/genenames.Rds")
```

We add this information as `rowData` within the `SingleCellExperiment` 
object created above. 

```{r}
## Merge biomaRt annotation
rowdata <- data.frame(ensembl_gene_id = rownames(droplet_sce))
rowdata <- merge(rowdata, genenames, by = "ensembl_gene_id", all.x = TRUE)
rownames(rowdata) <- rowdata$ensembl_gene_id
# Sort to match the order in the original data
rowdata <- rowdata[rownames(droplet_sce),]
## Check if  order is correct after merge;
stopifnot(all(rownames(rowdata) == rownames(droplet_sce)))
## add to the SingleCellExperiment object
rowData(droplet_sce) <- rowdata
```

For the remaining analysis, we will only focus on the 
`r format(sum(rowData(droplet_sce)$gene_biotype == "protein_coding", na.rm=TRUE), big.mark=",")`
protein coding genes that are contained in the data. These are selected below. 

```{r}
ind_pc <- which(rowData(droplet_sce)$gene_biotype == "protein_coding")
droplet_sce <- droplet_sce[ind_pc, ]
```

# QC and exploratory data analysis

In order to perform quality control (QC), we calculate some cell-specific
summary statistics: the total number of counts per cell (`sum`) and the total 
number of detected genes per cell (`detected`). These will be added in the
`colData` slot using the `addPerCellQC()` function from the `r Biocpkg("scater")`
package. Note that additional QC metrics are often calculated (e.g. the 
percentage of reads mapped to mitochondrial reads). These are not calculated 
here as the authors already applied relevant QC. 

```{r add-cellqc-metrics}
droplet_sce <- addPerCellQC(droplet_sce)
```

The distribution of these metrics per each cluster can be visualised as follows:

```{r cellqc-metrics, fig.cap="Cell-level QC metrics stratified by cell sub-type and sample (before QC). For each cell cluster: distribution of the total number of UMIs per cell (left) and the total number of detected genes per cell (right)."}
p_cell_qc1 <- plotColData(droplet_sce, y = "sum", x = "subCellType") +
  geom_hline(yintercept = 25000, lty = 2, col = "red") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
p_cell_qc2 <- plotColData(droplet_sce, y = "detected", x = "subCellType") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
p_cell_qc3 <- plotColData(droplet_sce, y = "sum", x = "sample") +
  geom_hline(yintercept = 25000, lty = 2, col = "red")
p_cell_qc4 <- plotColData(droplet_sce, y = "detected", x = "sample")

(p_cell_qc1 + p_cell_qc2) / (p_cell_qc3 + p_cell_qc4)
```

As seen in Figure \@ref(fig:cellqc-metrics), the  `presomiticMesoderm.b` 
sub-cluster is characterised by a substantially larger number of UMIs and 
detected genes. These cells have been highlighted by the authors as potential 
doublets. Therefore, as in the original publication, the entire sub-cluster will 
be removed from subsequent analyses. For this workflow, we also apply an 
additional exclusion criteria to remove cells with more 25,000 UMIs as potential 
doublets. 

```{r drop-filt}
ind_retain <- colData(droplet_sce)$subCellType != "presomiticMesoderm.b" &
  colData(droplet_sce)$sum <= 25000
droplet_sce <- droplet_sce[, ind_retain]
```

After applying this filter, the dataset contains `r ncol(droplet_sce)` cells:
`r sum(colData(droplet_sce)$subCellType == "presomiticMesoderm.a")` in the 
presomitic Mesoderm group and 
`r sum(colData(droplet_sce)$subCellType == "somiticMesoderm")` in the somitic
Mesoderm group. 

The distribution of the total number of UMIs and detected genes after QC is
shown in Figure \@ref(fig:cellqc-metrics-postqc).

```{r cellqc-metrics-postqc, fig.cap="Cell-level QC metrics stratified by cell sub-type and sample (after QC). For each cell cluster: distribution of the total number of UMIs per cell (left) and the total number of detected genes per cell (right)."}
p_cell_qc1 <- plotColData(droplet_sce, y = "sum", x = "subCellType")  +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
p_cell_qc2 <- plotColData(droplet_sce, y = "detected", x = "subCellType") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
p_cell_qc3 <- plotColData(droplet_sce, y = "sum", x = "sample") +
  geom_hline(yintercept = 25000, lty = 2, col = "red")
p_cell_qc4 <- plotColData(droplet_sce, y = "detected", x = "sample")

(p_cell_qc1 + p_cell_qc2) / (p_cell_qc3 + p_cell_qc4)
```

<!--- We removed all cells that expressed <1,000 genes or that had >3% of their 
transcripts mapped to mitochondrial genes. We further removed any cells that 
expressed both Xist and any of Kdm5d, Eif2s3y, Gm29650, Uty or Ddx3y (genes in 
the Y chromosome), as these are probably doublets. ---> 

To further explore the underlying structure of the data, we perform global
scaling normalisation using `r Biocpkg("scran")` and
principal component analysis (PCA) of log-transformed normalised expression 
counts using `r Biocpkg("scater")`.
As seen in Figure \@ref(fig:pca-visualisation-mesoderm-batch), this analysis 
suggests the absence of strong batch effects.
It should be noted that the estimation of global scaling normalisation factors 
using `r Biocpkg("scran")` is not strictly necessary in the 
`r Biocpkg("BASiCS")` workflow. 
Here, we only use it as part of the exploratory data analysis. 
Moreover, count-based models for dimensionality reduction 
(e.g. [@Townes2019;@Lopez2018]) could be used as an alternative to PCA, 
removing the need for log normalisation.

```{r pca-visualisation-mesoderm-batch, fig.cap="First two principal components of log-transformed expression counts after scran normalisation. Colour indicates the experimental condition (left), the number of detected genes (middle) and animal of origin (right) for each cell."}
## Global scaling normalisation 
droplet_sce <- computeSumFactors(droplet_sce, colData(droplet_sce)$subCellType)
## Obtain log-normalised expression counts
droplet_sce <- logNormCounts(droplet_sce)
## Run PCA
droplet_sce <- runPCA(droplet_sce)
## Plot first two PCs colour-coded according to cell-level metadata
p_type <- plotPCA(
    droplet_sce,
    ncomponents = c(2, 1),
    colour_by = "subCellType"
  ) +
  theme(legend.position = "bottom")
p_detected <- plotPCA(
    droplet_sce,
    ncomponents = c(2, 1),
    colour_by = "detected"
  ) +
  theme(legend.position = "bottom")
p_batch <- plotPCA(
    droplet_sce,
    ncomponents = c(2, 1),
    colour_by = "sample"
  ) +
  theme(legend.position = "bottom")
p_type + p_detected  + p_batch 
```


As seen inn Figure \@ref(fig:pca-visualisation-mesoderm-batch), the first 
principal component (which explains approximately 9\%
of the total variability) captures the presomitic versus somitic Mesoderm
structure of the data. The second principal component (which explains 
approximately 3\% of the total variability) is correlated with the number of
detected features per cell. We do not observe substantial batch effects
driven by the animal of origin from which cells were collected (although some
minor differences are seen, due to differences in the number of detected genes). 

In addition to cell-specific QC, we also recommend a gene filtering 
step prior to using `r Biocpkg("BASiCS")`. 
The purpose of this filter is to remove lowly expressed genes that were largely
undetected through sequencing, making reliable variability estimates difficult
to obtain.
Here, we remove genes that are not detected in at least 20 cells across all 
cells. This is to ensure a reliable estimate of variability, roughly in line 
with the sample size requirements for the negative binomial distribution 
outlined in @Lloyd-Smith2007.


```{r gene-selection, fig.cap="Average UMI-count for each gene is plotted against the number of cells in which that gene was detected. Dashed grey lines are shown at the thresholds below which genes are removed."}
## Calculate gene-level QC metrics using scater
droplet_sce <- addPerFeatureQC(droplet_sce, exprs_values = "counts")
## Remove genes with zero total counts across all cells
droplet_sce <- droplet_sce[rowData(droplet_sce)$detected != 0, ]

## Transform "detected" into number of cells and define inclusion criteria
rowData(droplet_sce)$detected_cells <-
  rowData(droplet_sce)$detected * ncol(droplet_sce) / 100

## Set detection threshold to detecting expression in at least 20 cells
detected_threshold <- 20
rowData(droplet_sce)$include_gene <- 
  rowData(droplet_sce)$detected_cells >= detected_threshold

## Plot gene-level metrics, highlighting those included for our analysis
plotRowData(droplet_sce,
    x = "detected_cells",
    y = "mean",
    colour_by = "include_gene"
  ) +
  xlab("Number of cells in which expression was detected") +
  ylab("Average number of read counts across all cells") +
  scale_x_log10() +
  scale_y_log10() +
  theme(legend.position = "bottom") +
  geom_vline(
    xintercept = detected_threshold,
    linetype = "dashed",
    col = "grey60"
  )

## Apply gene filter
droplet_sce <- droplet_sce[rowData(droplet_sce)$include_gene, ]
```

The final dataset used in subsequent analyses contains 
`r ncol(droplet_sce)` cells and `r nrow(droplet_sce)` genes.
As the analysis will be carried out separately for pre-somitic and somitic
mesoderm cells, we create separate `SingleCellExperiment` object for each 
group of cells. These datasets contain 
`r sum(colData(droplet_sce)$cellType == "presomiticMesoderm")` and 
`r sum(colData(droplet_sce)$cellType == "somiticMesoderm")`, respectively. 

```{r separate-SCE}
sce_psm <- droplet_sce[, colData(droplet_sce)$cellType == "presomiticMesoderm"]
sce_sm <- droplet_sce[, colData(droplet_sce)$cellType == "somiticMesoderm"]
```

These are then stored for their use in the main manuscript.

```{r}
saveRDS(sce_psm, "rds/sce_psm.Rds")
saveRDS(sce_sm, "rds/sce_sm.Rds")
```

# Session information

```{r SessionInfo}
sessionInfo()
```
