<img src="images/logo-sm.png" style="position:absolute;top:40px;right:10px;" width="200" />


# Introduction

In this workshop we will be using the dataset:-

> Caron M, St-Onge P, Sontag T, Wang YC et al. Single-cell analysis of childhood leukemia reveals a link between developmental states and ribosomal protein expression as a source of intra-individual heterogeneity. Sci Rep 2020 May 15;10(1):8079

The materials are based on materials kindly-shared by Harvard Bioinformatics Core and Cancer Research Uk Cambridge Institute

- [Harvard Bioinformatics Core](https://hbctraining.github.io/scRNA-seq_online/schedule/links-to-lessons.html)
- [Cancer Research Uk](https://bioinformatics-core-shared-training.github.io/Bitesize_SingleCell/)

The *processed* single-cell data for this study have been made available on Gene Expression Omnibus (GEO), and we have downloaded them for you and re-organised the data into folders. For more information, the GEO entry is [here](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE132509) and the code to download is given in the R script [download_GSE132509.R](download_GSE132509.R). **This code is for your information only, and you do not need to re-run**.

*We will only analyse a subset of the dataset during the practical*. However, the other samples are provided if you want extra practice in your own time.

![](images/flowchart1.png)


<div class="information">
Much of the code and information presented in this workshop is *exploratory* rather than describing a strict workflow that must always be followed. In particular, the some of the cut-offs and parameters used **may not be appropriate for your own dataset**. 
</div>

The two packages we are going to use are `Seurat` and `tidyverse`; which you should already have installed.

```{r}
library(Seurat)
library(tidyverse)
```

You will need to install these on your own computer after the course

```{r eval=FALSE}
install.packages("tidyverse")
install.packages("Seurat")
```

# Data Import

We start by importing the sample information for the dataset. The `sampleSheet` data frame contains information about each of the samples in the dataset; which GEO Accession ID (GSE) they correspond to, and the disease status of the sample.

```{r}
sampleSheet <- read.delim("GSE132509/meta_data.tsv")
sampleSheet
```

## Code to import a single biological sample

We are going to import the data for the first healthy sample - `GSM3872442` The code to read a single sample is given below. This assumes that the output from the *cellranger* tool can be found in the `GSM3872442` directory. This output consists of the files `barcodes.tsv.gz`, `features.tsv.gz` and `matrix.mtx.gz`. As part of the course setup ([`download_GSE132509.R`]("download_GSE132509.R")) we have already renamed the outputs to the convention that `Seurat` is expecting. If your data have been processed using `cellranger` they should already be in the correct location and format (presumably the files for this dataset had to be re-organised as part of the upload to GEO?).

```{r}
## check to see that the expected files are present
list.files("GSE132509/GSM3872442/")
```

It is then a two-step process to create a object that `Seurat` can use for analysis.

```{r}
## Read all the files into a 'sparse matrix'
seurat_data <- Read10X(data.dir = "GSE132509/GSM3872442")
## convert into Seurat object type for further analysis
seurat_obj <- CreateSeuratObject(counts = seurat_data, 
                                 min.features = 100, ### Cells must have this many genes measured
                                 min.cells = 10, ### Genes must be present in this many cells
                                 project = "PBMMC_1")
seurat_obj
```

Note that the "samples" referred to in the output are individual cells that have been sequenced. We can inspect the data that we have just imported and display the counts;, although most of them will be empty due to the *sparse* nature of single-cell data.

Note that we tend to look at summaries or plots of the data rather than the raw counts themselves. This is just to explain what the object that we have just created contains.

```{r}
## each row is a gene, each column is a sample
GetAssayData(seurat_obj, "counts")[1:10,1:2]
seurat_obj@meta.data
```

The "meta-data" in the object is used to record any QA data that we generate or other grouping information related to the samples (cells). We will see how to add to this data frame later. For now it contains the number of UMI (`nCount_RNA`) and genes (`nFeature_RNA`) detected in each cell. The `orig.ident` is used to identify which biological sample each cell originates from. It will be possible to store multiple samples in this object, as we will see in the following exercise.



******
******
******

#### Exercise
<div class="exercise">
- Use to sampleSheet to determine the directory containing data from the second healthy individual (`PBMMC_2`).
- Use the `Read10X` and `CreateSeuratObject` functions to import the data from this sample, and save as a new variable (`seurat_obj2` for example)
</div>

```{r}


```

We have now loaded the data from two separate samples into RStudio. For QC purposes it will be convenient to *merge* the two samples into one object. In practice, this means combining the counts matrix and meta information for both samples. If you have created the `seurat_obj2` from the exercise we can use the `merge` function to create a single object. At this point, it will also be useful to save the merged object to disk so we can start from the merged data at a later point.

```{r eval=FALSE}
merged_data <- merge(seurat_obj,seurat_obj2)
## create a directory to store the Robjects
dir.create("Robjects",showWarnings = FALSE)
saveRDS(merged_data, file = "Robjects/merged_data.RDS")
```

If you didn't manage to complete this exercise, the output is available in the course folder.

```{r }
merged_data <- readRDS("Robjects/merged_data_BACKUP.RDS")
```


<div class="information">
If you have more than two samples in your dataset, you might want to take a more efficient approach to reading the data. An example is given at the end of this section.
</div>

### Quality control

As with other high-throughput datasets, **QC is an important part of the workflow and should not be overlooked**. Whilst it is difficult to give hard cut-offs that will work for all experiments, the `ggplot2` and `dplyr` packages provide convenient infrastructure for generating QC summaries and plots. You can re-use the code presented here on your own data, but you will have to interpret the results and use your judgement. 

Furthermore, the QC process is a delicate balance between retaining as many high quality cells as possible and removing as many truly poor quality cells as possible. Ultimately we want the clusters in our dataset to reflect biological variability and not technical or uninteresting variation. If we are not careful we could mistake poor quality cells for cells that have low complexity but are still biologically-meaningful.

In order to remove or reduce the impact of poor-quality data on our downstream analysis we will attempt to filter using some QC metrics. The three principle means of doing this are to apply thresholds for inclusion on three characteristics:

- Cells with small library sizes are considered to be of low quality as the RNA has not been efficiently captured
- The number of expressed genes in each cell defined as the number of genes with non-zero counts for that cell; any cell with very few expressed genes is likely to be of poor quality as the diverse transcript population has not been successfully captured. 
-  The proportion of UMIs mapped to genes in the mitochondrial genome; high proportions are indicative of poor-quality cells, possibly because of loss of cytoplasmic RNA from perforated cells (the reasoning is that mitochondria are larger than individual transcript molecules and less likely to escape through tears in the cell membrane).

Prior to performing QC, you should have a good idea what type of cells you expect to find; especially with regards to mitochondrial genes and complexity. We must then take this biology into account when assessing the quality of our data.

The meta data stored in the merged data already contain some useful information about QC. We can add further to this by calculating the percentage of a particular set of genes that are expressed. A common use-case for this is Mitochondrial genes that can be indicative of cells under stress. The `PercentageFeatureSet` function can be used and makes use of a *regular expression* to identify the genes we are interested in. 


<div class="information">
The regular expression used here assumes that our Mitochondrial genes start with "MT". If this is not the case for your own data then you will need to modify the `pattern` argument in the below code.
</div>

```{r}
## we look for genes that start with MT
## this will work for human genes
merged_data[["mito.ratio"]] <- PercentageFeatureSet(merged_data, pattern="^MT")/100
merged_data@meta.data
```

The meta data can also include information about our sample grouping (which we already have in our `sampleSheet` object). The `left_join` function can be used to create a single data frame, compatible with `tidyverse`, to hold all the information about our samples. At this point we also do some re-naming which will helping with data interpretation.

```{r}
qc_data <- merged_data@meta.data %>% 
  data.frame %>% 
  left_join(sampleSheet, by = c("orig.ident"="SampleName")) %>% 
  dplyr::rename(SampleName = orig.ident) %>% 
  dplyr::rename(nUMI = nCount_RNA, nGene=nFeature_RNA)
```


```{r eval=FALSE}
View(qc_data)
```


### Number of cells per-sample

The number of cells for each sample can be informative - especially if we have some idea about how many cells to expect as part of the experimental protocol. We can make a barplot using a `geom_bar`. Using the `SampleName` as our `x` aesthetic will automatically put the number of occurrences of that name on the `y` aesthetic (hence we do not need to specify a `y` value inside `aes`. Since the sample group information is contained in the `qc_data` data frame we can also colour by the `SampleGroup` to see if there is any underlying difference.

```{r}
qc_data %>% 
  ggplot(aes(x = SampleName, fill=SampleName)) + geom_bar() + theme(axis.text.x = element_text(angle=45,vjust=1,hjust=1))
```

There is clearly some variation between samples, and  the number of cells in *PBMMC_1* sample is a bit lower. The number of cells does not tell the whole story. Even though *PBMMC_2* has more cells they might not all be usable. To explore this further we can look at the number of UMI and genes (and the relationship between these).

### UMI and gene counts per cell

To recap, from our `qc_data` we have the `uMI` for each cell and the `nGene` metric. We also know which sample each cell originated from. Using the `ggplot2` interface we can translate these quantities into a plot to visualise the distributions of `nUMI` and `nGene` on a per-sample basis. Thus we can see what a typical number of *UMI* and *genes* is for a sample.

A typical way to representing distributions is using a violin plot. We can choose to apply a log$_{10}$ transformation to the `nUMI` values to make them easier to compare.

```{r}
qc_data %>% 
  ggplot(aes(x = nUMI,y = SampleName,fill = SampleName)) + geom_violin() + scale_x_log10() 
```

An equivalent plot for `nGene` can be produced by changing the `x` aesthetic.

```{r}
qc_data %>% 
  ggplot(aes(x = nGene,y = SampleName,fill = SampleName)) + geom_violin() + scale_x_log10() 
```

Along with these graphical representations we can also produce numeric summaries:-

```{r}
qc_data %>% 
  group_by(SampleName) %>% 
  summarise_if(is.numeric, median)
```


### Complexity

We can now look at the relationship between the number of UMI and number of genes detected; which we would hope to be somewhat correlated.

<div class="exercise">
**Exercise:** Make a plot to compare the `nUMI` and `nGene` values for each sample. See below for an example
</div>


```{r}

```

![](images/uMI_vs_gene.png)


We would ideally like each cell to be measuring as many genes as possible. We can assess this by defining a *novelty* score as the ratio between the `nGene` and `nUMI` quantities. If this score is low for a particular cell it means that there is not a wide variety of genes that have been sequenced in the cell. Even though the `nUMI` might be high for such cells, the same genes have been sequenced repeatedly. Using `dplyr` and `ggplot2` we can add this new score and visualise it's distribution. **The score is expected to be around 0.8 for good quality cells**

```{r}
qc_data <- mutate(qc_data, Novelty = log10(nGene) / log10(nUMI))
ggplot(qc_data, aes(x = Novelty, y = SampleName)) + geom_violin()

```


By colouring the plot of UMI vs nGene by this new score we seem to highlight any outlier points. However, exactly where to drawn the line and exclude cells is not completely clear. 


```{r}

ggplot(qc_data, aes(x = log10(nUMI), y= log10(nGene),col=Novelty)) + geom_point() + facet_wrap(~SampleName)

```

We still have another piece of information to consider, which is the number of mitochondrial genes. Recall that we want to remove any cells where this percentage is quite high. 

## Mitochondrial counts

The code to visualise the distribution of mitochondrial counts is highly similar to that to visualise the `nUMI` and `nGene`.

```{r}
qc_data %>% 
  ggplot(aes(x=mito.ratio,y = SampleName)) + geom_violin() 
```

## Considering multiple QC metrics


Based on the QC we have generated, we decide on the following criteria for *keeping* a cell in the analysis

- `nUMI > 500`
- `nGene > 250`
- `Novelty > 0.8`
- `mito.ratio < 0.2`

We can now add an extra column to the `qc_data` to indicate whether a cell would be retained given these criteria.

```{r}
qc_data %>% 
  mutate(kept_in_analysis = nUMI >= 500 & nGene >=250 & Novelty > 0.8 & mito.ratio < 0.2) %>%
  count(kept_in_analysis)
```

We can also overlay this information on the plot of `nUMI` vs `nGene` to see if we have been successful in removing most of the outlier points.

```{r}
qc_data %>% 
  mutate(kept_in_analysis = nUMI >= 500 & nGene >=250 & Novelty > 0.8 & mito.ratio < 0.2) %>%
  ggplot(aes(x = log10(nUMI), y = log10(nGene), col=kept_in_analysis)) + geom_point() + facet_wrap(~SampleName)
```
As we would expect, most of the removed cells are from *PBMMC_2*. Exactly how many are removed can be calculated as follows:-

```{r}
qc_data %>% 
  mutate(kept_in_analysis = nUMI >= 500 & nGene >=250 & Novelty > 0.8 & mito.ratio < 0.2) %>%
  count(kept_in_analysis,SampleName)
```

The plot and table above we useful for visualisation purposes, and once we are happy we can proceed to filter our cells. This requires us to modify the data stored in the Seurat object and then use the `subset` function.

```{r}
### Put the new meta data we've created into the Seurat object

qc_data <- data.frame(qc_data)
rownames(qc_data) <- rownames(merged_data@meta.data)

merged_data <- AddMetaData(merged_data, qc_data)

filtered_data <- subset(x = merged_data, 
                         subset= (nUMI >= 500) & 
                           (nGene >= 250) & 
                           (Novelty > 0.80) & 
                           (mito.ratio < 0.20))

```

We can now save the Seurat data as an R object again and remove the merged data as we no longer need it

```{r}
saveRDS(filtered_data, file = "Robjects/seurat_filtered.RDS")
rm(merged_data)
```



# Appendix 

<div class="information">
The following section will be useful if you are required to import a dataset with many samples.

</div>

## Importing many samples using a `for` loop

Since our dataset comprises multiple samples, we need a way to import all the data and create a single Seurat object for analysis. The `Read10X` and `CreateSeuratObject` functions are designed to import the data from a single directory and create a corresponding Seurat object. To deal with multiple directories and samples we have to read them separately and use the `merge` function to join the datasets together.

```{r eval=FALSE}
## DON'T RUN THIS CODE
merged_data <- merge(x = seurat_obj,y=seurat_obj2)
```

**But how do we read data from multiple directories?** We will employ a "for loop" to process the remainder of the samples. In programming, looping is a strategy for repeating same lines of code but for a set of different inputs. To see why this is beneficial, lets look at the code for importing the second and third samples in our dataset:-

```{r eval=FALSE}
## No need to run this code!
## It is just for an example

seurat_data <- Read10X(data.dir = "GSE132509/GSM3872435")
seurat_obj <- CreateSeuratObject(counts = seurat_data, 
                                 min.features = 100, 
                                 project = "ETV6-RUNX1_2")
seurat_obj2

merged_data <- merge(x = seurat_obj,y=seurat_obj2)

seurat_data <- Read10X(data.dir = "GSE132509/GSM3872436")
seurat_obj2 <- CreateSeuratObject(counts = seurat_data, 
                                 min.features = 100, 
                                 project = "ETV6-RUNX1_3")
seurat_obj <- merged_data <- merge(x = seurat_obj,y=seurat_obj2)

```

The only changes required to import the data from a different biological sample are; 1) the directory name containing the raw data 2) the name assigned to the object.

We could manually create code chunks to read the entire set of samples in the same manner. This might be feasible for a few samples, but probably involve some copying and pasting and be prone to errors.

We can start constructing the `for` loop re-writing the code chunk using *variables* in place of the the directory name and project name. This is the only part of the code chunk that is specific to a particular sample.

```{r eval=FALSE}
seurat_data <- Read10X(data.dir = the_data_dir)
seurat_obj <- CreateSeuratObject(counts = seurat_data, 
                                 min.features = 100, 
                                 project = project_name)
```

If we want to read the data for a particular sample, all that is required is to set the values of `the_data_dir` and `project_name` appropriately, and run the rest of the code.

```{r eval=FALSE}
the_data_dir <- "GSE132509/GSM3872435"
project_name <- "ETV6-RUNX1_2"

seurat_data <- Read10X(data.dir = the_data_dir)
seurat_obj <- CreateSeuratObject(counts = seurat_data, 
                                 min.features = 100, 
                                 project = project_name)

```


The next stage is to define which samples we want to use as possible inputs to the code chunk above. To keep the data manageable we will concentrate on either `ETV6-RUNX1` or `PBMMC` samples.

```{r eval=FALSE}
## slice can be used to select particular rows from a data frame
dirs <- slice(sampleSheet, c(2:4,9:11)) %>% 
        pull(GSE)

samples <- slice(sampleSheet, c(2:4,9:11)) %>% 
        pull(SampleName)
```

```{r eval=FALSE}
dirs
samples
```
Both `dirs` and `samples` are examples of a vector data type in R. If we want to know the item in the 1st of second element of `dirs` we can do

```{r eval=FALSE}
dirs[1]
dirs[2]
```

The value inside the `[]` can itself be a variable rather than a number

```{r eval=FALSE}
i <- 1
dirs[i]
```


We "loop" over the `dirs` variable. This will use a indexed value of `dirs` starting at 1 until all values of `dirs` have been processed. At each *iteration* of the loop we can find out the current value of `dirs` and `samples` using `dirs[i]` and `samples[i]` respectively.

However, the previous code would over-write the Seurat object at each iteration. The solution is to use a `list` 

```{r warning=FALSE,eval=FALSE}
## create an empty list to store the seurat objects.
seurat_list <- NULL

# Create a Seurat object for each sample
# Assume we have already read the first sample, and loop over the remainder

dirs <- dirs[-1]
samples <- samples[-1]

for (i in 1:length(dirs)){
        the_data_dir <- dirs[i]
        project_name <- samples[i]
        
        message(paste("Reading data from ", dirs[i]))
        ## paste0 used to create the relative path to the data directory
        seurat_data <- Read10X(data.dir = paste0("GSE132509/",the_data_dir))
        seurat_list[[i]] <- CreateSeuratObject(counts = seurat_data, 
                                         min.features = 100, 
                                         project = project_name)
        
}


```

The `merge` function that we saw previously can now be used to combine all the datasets into a single object. The second argument, `y`, can be a `list` of `Seurat` objects.

```{r eval=FALSE}
merged_data <- merge(x = seurat_obj,y=seurat_list)
```

We will now work with the merged dataset, so no longer require the `seurat_obj` and `seurat_data` corresponding to particular samples. These objects are using some of the memory available to RStudio, so we can delete them.

```{r eval=FALSE}
rm(seurat_list)
rm(seurat_data)
rm(seurat_obj)
gc()
```
We can save this object to disk so we don't have to repeat this step

```{r eval=FALSE}
dir.create("Robjects",showWarnings = FALSE)
saveRDS(merged_data, file = "Robjects/merged_data.RDS")
```

```{r eval=FALSE}
merged_data <- readRDS("Robjects/merged_data.RDS")
```