---
title: "scPipe: a flexible data preprocessing pipeline for 3' end scRNA-seq data"
author: "Luyi Tian"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{scPipe: flexible data preprocessing pipeline for 3' end scRNA-seq data}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

# Introduction

`scPipe` is a package designed to process single-cell RNA-sequencing (scRNA-seq) data generated by different protocols. It is designed for protocols with UMI, but can also adapt to non-UMI protocols. `scPipe` consist of two major components. The first is data preprocessing with raw fastq as input and a gene count matrix as output. The second component starts from the gene count matrix, includes quality control and a shiny app for clustering and visualization.

# Case Study: Preprocessing CEL-seq2 data

CEL-seq2 is a modified version of the CEL-seq method with higher sensitivity and lower cost. It uses a polyT primer with index sequences to capture a cell's mRNA. The index sequence consists of the a cell specific index, which is designed, and a molecular index, also called unique molecular identifier (UMI), which is a random sequence that should be different for each primer (the number of distinct UMIs will depend on its length). The ouput fastq files from a CEL-seq2 experiment is paired-ended where the first read contains the transcript and second read contains the cellindex+UMI with some polyA sequence.

## Getting started

We begin by creating a folder to store the results.

```{r, message=FALSE}
library(scPipe)
library(SingleCellExperiment)
data_dir = tempdir()
```

## Organising the files required

To process the data, we need the genome fasta file, gff3 exon annotation and a cell barcode annotation. The barcode annotation should be a `.csv` file with at least two columns, where the first column has the cell id and the second column contains the barcode sequence. We use ERCC spike-in genes for this demo. All files can be found in the `extdata` folder of the `scPipe` package:

```{r}
# file path:
ERCCfa_fn = system.file("extdata", "ERCC92.fa", package = "scPipe")
ERCCanno_fn = system.file("extdata", "ERCC92_anno.gff3", package = "scPipe")
barcode_annotation_fn = system.file("extdata", "barcode_anno.csv", package = "scPipe")
```

The read structure for this example is paired-ended, with the first longer read mapping to transcripts and second shorter read consisting of 6bp UMIs followed by 8bp cell barcodes. **NOTE**: by convention, `scPipe` always assumes `read1` refers to the read with the transcript sequence, which is usually the longer read. These data are also available in the in `extdata` folder:

```{r,eval=TRUE}
fq_R1 = system.file("extdata", "simu_R1.fastq.gz", package = "scPipe")
fq_R2 = system.file("extdata", "simu_R2.fastq.gz", package = "scPipe")
```

## Data Preprocessing

### Fastq reformatting

The pipeline starts with fastq file reformatting. We move the barcode and UMI sequences to the read name and leave the transcript sequence as is. This outputs a read name that looks like `@[barcode_sequence]*[UMI_sequence]#[readname]` ... The read structure of read 2 in our example dataset has the 8 bp long cell barcode starting at position 6 and the 6 bp long UMI sequence starting at the first position. So the read structure will be : `list(bs1=-1, bl1=0, bs2=6, bl2=8, us=0, ul=6)`. `bs1=-1, bl1=0` means we don't have an index in read 1 so we set a negative value as its start position and give it zero length. `bs2=6, bl2=8` means we have an index in read two which starts at position 6 in the read and is 8 bases in length. `us=0, ul=6` means we have a UMI at the start of read 2 which is 6 bases long. **NOTE**: we use a zero based index system, so the indexing of the sequence starts at zero.

```{r,eval=TRUE}
sc_trim_barcode(file.path(data_dir, "combined.fastq.gz"),
                fq_R1,
                fq_R2,
                read_structure = list(bs1=-1, bl1=0, bs2=6, bl2=8, us=0, ul=6))
```

### Aligning reads to a reference genome

Next we align reads to the genome. This example uses `Rsubread` but any aligner that support RNA-seq alignment and gives standard BAM output can be used here. The `Rsubread` package is not available on Windows, so the following alignment step will only run on Mac OSX or Linux.

```{r,eval=TRUE}
if(.Platform$OS.type != "windows"){
  Rsubread::buildindex(basename=file.path(data_dir, "ERCC_index"), reference=ERCCfa_fn)

  Rsubread::align(index=file.path(data_dir, "ERCC_index"),
      readfile1=file.path(data_dir, "combined.fastq.gz"),
      output_file=file.path(data_dir, "out.aln.bam"), phredOffset=64)
}
```

### Assigning reads to annotated exons

After the read alignment, we assign reads to exons based on the annotation provided using the `sc_exon_mapping` function. Currently `scPipe` only supports annotation in `gff3` or `bed` format.

```{r,eval=TRUE}
if(.Platform$OS.type != "windows"){
  sc_exon_mapping(file.path(data_dir, "out.aln.bam"),
                file.path(data_dir, "out.map.bam"),
                ERCCanno_fn)
}
```

### De-multiplexing data and counting genes

Next we use the `sc_demultiplex` function to split the reads into separate `.csv` files for each cell in the `/count` subfolder. Each file contains three columns and each row corresponds to a distinct read. The first column contains the gene ID that the read maps to, the second column gives the UMI sequence for that read and the third column is the distance to the transcript end position (TES) in bp. This file will be used for UMI deduplication and to generate a gene count matrix by calling the `sc_gene_counting` function.

```{r,eval=TRUE}
if(.Platform$OS.type != "windows"){
  sc_demultiplex(file.path(data_dir, "out.map.bam"), data_dir, barcode_annotation_fn, has_UMI=FALSE)

  sc_gene_counting(data_dir, barcode_annotation_fn)
}
```

We have now completed the preprocessing steps. The gene count matrix is available as a `.csv` file in `data_dir/gene_count.csv` and quality control statistics are saved in the `data_dir/stat` folder. These data are useful for later quality control (QC).


# Preprocessing data generated by other protocols

The greatest difference between scRNA-seq protocols is in their read structure, which is specified by the `read_structure` argument of the `sc_trim_barcode` function in `scPipe`. The major difference between CEL-seq and Drop-seq/Chromium 10x is that the cell barcode is unknown in advance for the latter protocols, so an extra step `sc_detect_bc` is required before running `sc_demultiplex`, to identify and generate the cell barcode annotation. `scPipe` is not optimized for the full-lengh protocols like SMART-seq, but it can process the data generated by such protocols by setting the `has_UMI` to `FALSE` in the `sc_demultiplex` function.

## Quality Control

The easiest way to create a *SingleCellExperiment* object from the output of `scPipe` preprocessing is using `create_sce_by_dir` function, that will read in the gene count matrix together with the QC information available in the `stat` folder.

```{r}
if(.Platform$OS.type != "windows"){
sce = create_sce_by_dir(data_dir)
dim(sce)
}
```

The dataset we analysed above used ERCC simulated reads from 10 cells of perfect quality. In order to demonstrate QC on a more reaslitic example, we have included an experimental dataset in the `scPipe` software generated by Dr Christine Biben. This experiment profiles gene expression in 383 blood cells for a subset of 1000 genes (to save space). We first create a *SingleCellExperiment* object directly without using the `create_sce_by_dir` function:

```{r}
data("sc_sample_data")
data("sc_sample_qc")
sce = SingleCellExperiment(assays = list(counts = as.matrix(sc_sample_data))) # generate new sce with gene count matrix
QC_metrics(sce) = sc_sample_qc
demultiplex_info(sce) = cell_barcode_matching
UMI_dup_info(sce) = UMI_duplication
```

There are several plots we can create to assess the overall quality of the experiment. The first is to look at the cell barcode demultiplexing statistics. To do this we generate a bar plot that shows the percentage of reads that uniquely match to the cell barcodes, as well as the unmatched proportion of reads and their alignment rates to introns and exons. If we observe a large proportion of unmatched reads that map to exons, this indicates a failure of the cell barcode demultiplexing.

```{r, fig.height=7, fig.width=7}
plot_demultiplex(sce)
```

A second plot shows the duplication rate which can be used to evaluate read deapth. UMIs are routinely used to mark individual molecules and after PCR amplification, the same molecule with have multiple copys, which can be identified and removed if they are observed to have the same UMI sequence. Therefore, the copy number of UMIs is an indication of the PCR amplification rate. 

```{r, fig.height=7, fig.width=7}
plot_UMI_dup(sce)
```

Next we calculate QC metrics and use the `detect_outlier` function to identify poor quality cells. The `detect_outlier` function has argument `comp` to define the maximum component of the gaussian mixture model. Using the default value of 1 is generally sufficient, but in cases where the data are heterogeneous in terms of quality control metrics, setting this value to 2 or 3 can give better results. This function will remove low quality cells if `type="low"`, large cells if `type="high"` or both when `type="both"`. The `conf` argument specifies the lower and upper confidence intervals for outliers and `detect_outlier` is insensistive to the interval values.

```{r, warning=FALSE, message=FALSE}
sce = calculate_QC_metrics(sce)
sce = detect_outlier(sce)
```

We can plot the alignment statistics per sample as a barplot using

```{r, fig.height=7, fig.width=7}
plot_mapping(sce, percentage = TRUE, dataname = "sc_sample_data")
```

and generate pairwise plots of the QC metrics as follows:

```{r, warning=FALSE, message=FALSE, fig.height=7, fig.width=7}
plot_QC_pairs(sce)
```

The final step will be to remove the low quality cells by the `remove_outliers` function. 

```{r}
sce = remove_outliers(sce)
dim(sce)
```

# Downstream analysis

Since the **scater** and **scran** packages both use the *SingleCellExperiment* class, it will be easy to further process this data using these packages for normalization and visualization. Other packages such as **SC3** may be useful for clustering and **MAST** and **edgeR** for differential expression analysis.

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