```{r setup, include=FALSE} 
knitr::opts_chunk$set(warning = FALSE, message = FALSE) 
``` 
<img src="images/logo-sm.png" style="position:absolute;top:40px;right:10px;" width="200" />

# Introduction

In Session 1, we read in our single cell RNA-seq data from 10X and created a Seurat object in R. We then performed some quality control to assess for poor quality/empty cells, and subset the data according to the QC plots we generated. 

In this Session, we are going to normalise the data, select the genes/features that are most variable between cells and scale the data, regressing out unwanted sources of variation. We are also going to visualise out data using t-SNE and UMAP plots and cluster our cells

# Learning Objectives

* Identify unwanted sources of variation in the data

* Learn how to normalise, find variable features and scale our data

* Learn how to perform PCA on the data and why we do this

* Use SCTransform to transform data and regress out unwanted variation 

* Integrate multiple datasets to remove batch effects 

* t-SNE and UMAP plots

* Cluster our cells into meaningful groups (by cell type)

![](images/flowchart2.png)

# Pre-processing
## Sources of variation

Firstly, lets load in our data from the first session if we haven't got it already loaded.


```{r}
#Re-load libraries if needed
library(Seurat)
library(tidyverse)
```

```{r}
#Reload R objects if needed
seurat <- readRDS("Robjects/seurat_filtered_BACKUP.RDS")
```

Our end goal is to identify groups of cells that show interesting biological variation. However at the moment our data contains lots of potentially **unwanted variation** and **noise**. We have expression of many genes across many cells in multiple samples! This is an example of *high dimensional data*.

Potential sources of variation -

* Differences in sequencing depth (This is also common to bulk seq)
* Influential cell cycle genes can drive a lot of variation
* Poor quality cells (Hopefully we have removed these in the first session!)
* Differences between samples

## Normalisation

First of all we want to be able to compare gene expression between cells, do do this we need to normalise for sequencing depth. Currently each cell has a different number for `nCount_RNA` (the number of UMI).


```{r}
# We use the $ operator to show us some of the information stored in the 'seurat data' object
#seurat$ tab will show you the list
  
summary(seurat$nCount_RNA)
```

We can use the `NormalizeData` function to use log normalisation to take the sequencing depth into account. 

Briefly, this normalises the gene expression measurements for each cell by the total expression, multiplies this by a scale factor, and log-transforms the result. Now we can directly compare expression between cells.

```{r}
seurat <- NormalizeData(seurat)

```

## Find Variable Features

At this point we want to define the features or genes that are the **most variable** in our dataset. ie, the genes that are showing most differences in expression within our dataset and are likely to be the most informative for defining clusters of cell types. 

We find features that define clusters, using the `FindVariableFeatures` function. Seurat does this by calculating a normalised intensity for each gene, then selecting the top 'n' most variable features (usually 2000). The samples in our dataset derive from human peripheral blood mononuclear cells (PBMCs), so seeing immune-related genes in the top 20 is in-line with what we would expect.

```{r}
# Identify the most variable genes
seurat <- FindVariableFeatures(seurat)

# Identify the top 20 most highly variable genes. We need the VariableFeatures function to pull the names
top20 <- head(VariableFeatures(seurat), 20)
top20
```

```{r}
# Plot the variable features
VarPlot<- VariableFeaturePlot(seurat, )
LabelPoints(plot = VarPlot, points = top20, repel = TRUE)

```
## Scale

The next step is to scale the data - which applies a linear transformation to the data that is a standard pre-processing step prior to dimensional reduction techniques like PCA. It stops highly expressed genes dominating the analysis.

```{r}
seurat <- ScaleData(seurat)
```

Now the data is pre-processed and ready for us to move onto the dimensionality reduction techniques.

# PCA

As we mentioned earlier, the data are *high dimensional data*. This basically means we have many observations (cells) and many features (genes). This makes it impossible to discern patterns in expression to a human viewer! We need to use dimensional reduction techniques to increase the **interpretability** of the data - *without* losing too much meaningful information.

To do this we use a Principal Component Analysis or **PCA**. This breaks down multivariate data table as smaller set of variables or summary indices in order to observe patterns and outliers. It has many uses for interpreting multivariate data.

<div class="information">
If you are interested in learning more about PCA and want to understand the mathematical principles behind the method, there are plenty of online resources to explore further:

- https://builtin.com/data-science/step-step-explanation-principal-component-analysis

- Youtube video from StatQuest - https://www.youtube.com/watch?v=FgakZw6K1QQ

</div>

The good news is that, despite the complex mathematics going-on behind the scenes, this is a simple command to run.

```{r, message=TRUE}
# Perform PCA
seurat <- RunPCA(seurat)
```
The command prints out some information on the **principal components** ie PC_1, PC_1 and groups the genes that are behaving similarly. PC1 representing the most variation in the data and PC2 representing the second most variation in the data

## Elbowplot

We can see how much variance is explained by each principal component by using the`ElbowPlot` function. The point at which the plot starts to level out is usually taken as guidance for the number of principal components to use in the analysis.

```{r}
ElbowPlot(seurat, reduction = "pca", ndims = 50)
```

Looking at the above 'elbow plot', we can see that the top 9 PCs contain a lot of the information, while other PCs contain progressively less. However, it is still advisable to use more PCs in further analysis since they might contain information about rare cell types.


We can look at the most variable genes driving the PCA:

```{r}
print(x = seurat[["pca"]], 
      dims = 1:5, #Top 5 PCs 
      nfeatures = 5) #Top 5 Genes
```

In PC1 we can start to see known cell type markers such as CD3D, a T Cell marker. **Encouraging signs!**

## Dimensionality Reduction Plots (DimPlot)

Now that we have carried out a PCA on our dataset we can visualise our whole dataset for the first time!

To do that we need to use a Seurat function called `DimPlot` . This is short for Dimensional reduction plot and shows us a 2D scatter plot showing each cell according to the dimensionality reduction technique carried out on the data. 

In this case we have done a PCA and we can look at the first two principal components. By default, cells are colored by their identity class (this can be changed with the `group.by` parameter).

```{r}
DimPlot(seurat ,reduction = "pca")

DimPlot(seurat ,reduction = "pca", group.by = "SampleName")

```


# Feature Regression

## Cell Cycle Genes

**Cell cycle genes** can often have strong effects within groups of cells and cause a lot of the inter-cellular variation. This results with cells being clustered by their cell cycle stage instead of the cell type or condition.

To mitigate this, we can **identify** cell cycle stage by using known marker genes of the cell cycle and **regressing out** these effects.

![*We want to avoid this!*](images/cc_score-2.png)

Helpfully, the known cell phase marker genes from Tirosh et al, 2015 are included in Seurat by default.

```{r}
#Here they are!
cc.genes


s.genes <- cc.genes$s.genes
g2m.genes <- cc.genes$g2m.genes

# We also have this list saved if it hasn't loaded
# read.delim("regev_lab_cell_cycle_genes.txt")[,1]
```
<div class="information">
If you are interesting in cell cycle marker genes from organisms these can be found here https://github.com/hbc/tinyatlas/tree/master/cell_cycle
</div>

We can now use the `CellCycleScoring` function in Seurat to 'score' or identify each cell's cell cycle phase based on the expression of these genes within that cell.

```{r}
seurat <- CellCycleScoring(seurat, 
                                 g2m.features = g2m.genes, 
                                 s.features = s.genes)

                              
seurat@meta.data %>% select(S.Score, G2M.Score, Phase)
```

A scatter plot can be used to understand how the Phase categories are assigned based on the scores

```{r}
seurat@meta.data %>% select(S.Score, G2M.Score, Phase) %>% 
  ggplot(aes(x = S.Score, y = G2M.Score, col = Phase)) + geom_point()
```

We now have 3 new columns in our metadata. The highest score between S.Score and G2M.Score determines the phase. If both columns are negative, then the phase is G1. A bar plot shows us how many cells are in each phase. 

```{r}
as_tibble(seurat[[]]) %>% ggplot(aes(x = Phase, fill = Phase))+ geom_bar()
```

<div class="exercise">
**Exercise:** Now that we have categorized our cells by their cell cycle phase. We can visualise how the cell cycle drives variation in our data by looking at the PCA.

1. Have a go at visualizing the PCA data using a dimensionality reduction plot and colour the cells by their cell cycle phase

**hint** we have used the command for this earlier...
</div>

```{r}

# Exercise work in here

```

## Mitochondrial genes 
Mitochondrial genes expression can also be a big contributor of variation. To visualise the effects of mitochondrial gene expression on the cells expression profiles we can **categorise** the percent of mitochondrial genes into low, medium and high. `Dimplot` needs to work with **categorical data** to colour the cells by group (rather than continuous data like the raw percent scores).

```{r}
# Categorize the % of mitochondrial genes expressed in cells into low, medium and high
seurat$mito.cat <- cut(seurat$mito.ratio, breaks = 3, labels=c("Low","Medium", "High"))

DimPlot(seurat,  reduction = "pca",  group.by= "mito.cat")

```


# SCTransform

We decide to **regress** for mitochondrial ratio and cell-cycle status. This prevents these variables from contributing so much to the PCA dimensional reduction.

We need to do this for each sample separately (i.e. not the merged dataset), so first we "split" the object into a list of subsetted objects.

`SCTransform` is a newer addition to the `Seurat` package. It allows us to do all the pre-processing - the data normalisation, data scaling & finding of variable features in one simple step. 

It is a separate and more statistically robust version of the steps we have done above. Importantly it allows us to regress out any sources of variation we have found in the above steps.


```{r}
seurat_split <- SplitObject(seurat, split.by = "SampleName")

seurat_split$PBMMC_1 <- SCTransform(seurat_split$PBMMC_1, vars.to.regress = c("mito.ratio","S.Score","G2M.Score"), vst.flavor = "v1" , verbose = F)

seurat_split$PBMMC_2 <- SCTransform(seurat_split$PBMMC_2, vars.to.regress = c("mito.ratio","S.Score","G2M.Score"), vst.flavor = "v1" , verbose = F)
  
```

<div class="information">

If you have a larger number of samples, you will probably find it more efficient to use a `for` loop to run this step. We briefly covered for loops in the Appendix of Session 1. The code would look like this:-

```
for(i in 1:length(seurat_split)){
  seurat_split[[i]] <- SCTransform(seurat_split[[i]], vars.to.regress = c("mito.ratio","S.Score","G2M.Score"), vst.flavor = "v1" , verbose = F)
}
```

</div>

<div class="information">

For your own dataset, you will need to consider what variables need to be regressed. There may be other sources of variation that we have not considered for this dataset. If in doubt - consult your local Bioinformatician.
</div>


# Integration

Before we can carry out a thorough analysis of the single cell clusters and visualise them we need to address the possibility that **batch effects** are influencing our datasets.

### What is a batch effect?

<div class="information">
A batch effect occurs when technical factors in an experiment cause changes in the data produced by the experiment. They have been observed throughout the entire history of high-throughput gene expression studies and are still as relevant in modern technologies.

In this case we have two separate samples. These could have been produced in the lab on separate days by different people using different techniques, equipment etc. They could even be from **entirely different experiments**.

This can produce variation **between** samples that is not informative. We need to remove this effect.
</div>

Seurat has developed a method specifically for dealing with this issue - integration. If our ultimate goal is to compare cell-type expression between groups, across multiple samples we'll need to integrate the data. 

If you go back to look at older single cell papers, you'll probably see that they haven't done this step, but it is now recognized to be more important. It is optional, but recommended. 

There are 4 steps for integration: 

1. Select variable features for integration by taking the highly variable genes from each sample.

2. Prepare the `SCTransform` object for integration

3. Identify 'anchors' between cells. These are mutual nearest neighbours across the data. For each cell, the closest neighbour is identified based on gene expression. If the neighbour's closest neighbour is the original cell, these will be marked as 'anchors' for the data. 

4. Integrate the conditions / samples by using the anchors along with their corresponding scores to transform cell expression values. 

```{r}
#The complete command list for this is below.

# Select variable features for integration
integ_features <- SelectIntegrationFeatures(object.list = seurat_split, nfeatures = 3000, verbose=F)

# Prepare the SCTransform object for integration
seurat_split <- PrepSCTIntegration(object.list = seurat_split, anchor.features = integ_features, verbose=F)

# Identify mutual nearest neighbours across datasets
integ_anchors <- FindIntegrationAnchors(object.list = seurat_split, 
                                        normalization.method = "SCT",
                                        anchor.features = integ_features, 
                                        verbose=F)
# Integrate across samples
seurat_integrated <- IntegrateData(anchorset = integ_anchors,normalization.method = "SCT", verbose=F)
```


We need to re-run PCA on our new integrated dataset. The data is now different as we have regressed out mitochondrial effects, cell cycle effects and removed 'batch' effects between samples by integration.

After PCA we should hopefully have the most informative genes for our experiment within the top principal components.

```{r}
# Run PCA
seurat_integrated <- RunPCA(object = seurat_integrated)

```

# t-SNE and UMAP plots

t–Stochastic Neighbourhood Embedding (t-SNE) and Uniform Manifold Approximation and Projection (UMAP) are dimensionality reduction techniques which allow us to **graphically** simplify a large dataset. They are similar to PCA however they are non linear and produce a 2D separation of the data, rather than producing many components. This allows us to view the data on a 2D plot.

Their main function is to calculate the 'distance' between a cell's expression profile compared to a number of its neighbours and represent that information in a lower dimension - i.e to simplify it.

### The core differences

t-SNE uses a Gaussian probability function to calculate the likelhood of two cells being neighbours (similar in expression), and repeats this step for all cells. Cells are rearranged according to these distances, creating the t-SNE plot. 

UMAP, creates a high dimensionsal graph of the data using 'weights' to determine the likelhood of cells being in close proximity then compresses this information into a low dimensional space.


<div class="information">
If you are interested in learning more about t-SNA and UMAP and want to understand the mathematical principles behind the methods then there are plenty of online resources available. Here are some.

https://umap-learn.readthedocs.io/en/latest/how_umap_works.html
https://www.datacamp.com/tutorial/introduction-t-sne

UMAP Youtube video from StatQuest - https://www.youtube.com/watch?v=eN0wFzBA4Sc

t-SNE Youtube video from StatQuest - https://www.youtube.com/watch?v=NEaUSP4YerM

</div>


### Which to use?

* t-SNE and UMAP are both non-linear, graph-based methods for dimensionality reduction.
* Both are both for data visualization only.
* Similar methodology
* UMAP is less computationally intensive
* Overall, UMAP is a newer technique but it does a better job of preserving the global data structure and representing it graphically. 

```{r}
#Lets run them both!
seurat_integrated <- RunTSNE(object = seurat_integrated)
seurat_integrated <- RunUMAP(object = seurat_integrated, dims = 1:20) # RunUMAP needs to know which PCA 
```



<div class="exercise">
**Exercise:** 

Lets look at are new UMAP and t-SNE reduction plots. We can use the dimensionality reduction plot `DimPlot` to visualise these. The `reduction=` parameter will need to be changed....


1. Make plots for the UMAP and t-SNE reduced data.
2. See if you can colour the cells by cell phase and sample group (`group.by=`)


</div>


```{r}
# Exercise work in here

# Remember our object is now called seurat_integrated

DimPlot(seurat_integrated, reduction = "tsne", group.by = "Phase")
DimPlot(seurat_integrated)


```

We should start to see an organised **structuring** or our data now. The cells are starting to form visual **clusters** and **groupings**. We will be able to formally define these clusters in the next steps.

## Feature Plots 

Now that we have our t-SNE and UMAP reductions we can explore additional metrics and see how they influence our groups of cells. 

Feature plots are similar to `DimPlots`, except they visualise continuous *features*. For example gene expression. We can also examine other continuous data such as number of UMI, number of genes etc...

```{r}
FeaturePlot(seurat_integrated, 
            reduction = "umap", 
            features = c("nUMI", "nGene", "S.Score", "G2M.Score"),
            label = TRUE)

```
We can also take a look here at a particular gene of interest...


```{r}

#CD3D is a T cell marker gene

FeaturePlot(seurat_integrated, 
            reduction = "umap", 
            features = c("CD3D"),
            label = TRUE)

## I wonder where our T cells  are....?

```

# Clusters

It is now time to *formally* categorize our cell groupings into **statistically defined clusters**.

Briefly, this is done in Seurat by using a shared nearest neighbor (SNN) modularity optimization based clustering algorithm. 

1. First calculate k-nearest neighbors and construct the SNN graph using `FindNeighbors`. 
2. Then optimize the modularity function to determine clusters using `FindClusters`.

<div class="information">

For a full description of the algorithms, see Waltman and van Eck (2013) 

There is also a "stat quest" video on the similar K-nearest neighbours here https://www.youtube.com/watch?v=HVXime0nQeI

</div>

```{r}
seurat_integrated <- FindNeighbors(seurat_integrated, dims = 1:20, verbose = FALSE)
seurat_integrated <- FindClusters(seurat_integrated, verbose = FALSE, resolution=0.3)
```

We can use the `DimPlot` function again to plot our new clusters. 

```{r}
DimPlot(seurat_integrated, label = TRUE)
```

Finally, we save the Seurat data as an R object again.

```{r}
saveRDS(seurat_integrated, file = "Robjects/seurat_integrated.RDS")
```



