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

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  tidy = TRUE,
  tidy.opts = list(width.cutoff = 95),
  message = FALSE,
  warning = FALSE,
  fig.width = 10,
  time_it = TRUE
)

```

# Introduction

To re-cap -  we have obtained our raw counts matrices from GEO (or cellranger). We have performed QC and subset the data according to the QC plots we generated. We have normalised the data via `SCTransform` - regressing out unwanted sources of variation, and integrated across our datasets. Lastly, we have clustered the cells into separate groups based on the similarity of their expression profiles.

In the next part of the course we are going to explore the genes that show the most differences in their expression between clusters. The top genes showing the biggest change in their expression compared to the other clusters are often called *marker genes*. We can then use these marker genes to assign an *identity class* to our clusters of interest.

![](images/flowchart3.png)

## Learning Objectives

* Find differentially expressed genes using `FindMarkers()` and `FindAllMarkers()`.
* Use the bioconductor human annotation package `org.Hs.eg.db` to add some more useful annotation to our marker genes.
* Assign cluster identity using cell type canonical marker genes
* Carry out a simple pathway analysis

Firstly, we will load up the packages we have been using so far for this part of the course.

```{r}
library(Seurat, quietly = T)
library(tidyverse, quietly = T)
```
And then we can load up the R objects we have used in the previous sessions.
```{r}
seurat <- readRDS("Robjects/seurat_integrated.RDS")
```

We can simply run the name of the Seurat object (in our case `seurat`) to give us a summary of what the object contains and what calculations have been performed on it.

```{r}
seurat
```
We can recap our clusters by taking a look at the dimensionality reduction plot - `DimPlot`

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


# Find differentially expressed genes

Our next objective is to start assigning some *biological meaning* to our data. We have some nice clusters, but no clue yet as to what cell types our clusters describe.

In Seurat we can find the top differentially-expressed genes, or *markers*, within each cluster. The function `FindAllMarkers()` handles most of the differential expression testing within the Seurat package and finds marker genes across all clusters. As a default, Seurat performs differential expression based on the non-parameteric Wilcoxon rank sum test, but other methods are available within the function. 


Usually, we are wanting to find the most differentially-expressed genes in each cluster across the whole dataset. For this we use `FindAllMarkers`. We will only report the positively expressed markers with a minimum log fold change of 0.25.

```{r}
markers <- FindAllMarkers(seurat, 
                          only.pos = TRUE, 
                          min.pct = 0.25, 
                          logfc.threshold = 0.25, 
                          verbose = F) # keeps it quiet!
markers %>% head()
```

A bit of manipulation in `dplyr` enables us to get the top 5 DE genes in each cluster.

```{r}
markers %>%
  select(gene, cluster, avg_log2FC, p_val, p_val_adj) %>%
    group_by(cluster) %>% # important to look in each cluster!
    slice_max(n = 5, order_by = avg_log2FC)
```


# Annotating marker genes

We may want to share our findings with collaborators or perform further bioinformatic analysis on our DE genes in clusters of interest. In this dataset, the gene level information is given as a gene name. Whilst this is useful we may want to embellish these results with some gene description information or some other gene IDs (eg Ensembl, Entrez).

The `org.Hs.eg.db` database is a genome wide annotation package for human that has many different gene mappings based on Entrez Gene identifiers. These are up-to-date databases of gene annotation that can be installed like any other bioconductor package

```{r}
# Dont run this if you have run already!
#install.packages("BiocManager")
#BiocManager::install("org.Hs.eg.db")

library(org.Hs.eg.db)
```

In our case we are going to use the gene symbols (rownames in our Seurat object) to extract the gene description and Entrez ID from the database.

We are going to filter the database by a key or set of keys in order to extract the information we want. Valid names for the key can be retrieved with the `keytypes` function.

```{r}
keytypes(org.Hs.eg.db)
```

We are going to use the keytype 'SYMBOL' which designates the gene names that we have in our rownames of the Seurat object.

Unfortunately, the authors of `dplyr` and `AnnotationDbi` have both decided to use the name `select` in their packages. To avoid confusion and problems from packages attempting to use the wrong tool, the following code is sometimes used:-

`AnnotationDbi::select` which tells us to use the 'select' function from the `AnnotationDBI` package.
```{r}
anno <- AnnotationDbi::select(org.Hs.eg.db,
                              keys=rownames(seurat),
                              columns=c("SYMBOL","GENENAME", "ENTREZID"),  # info we want to extract
                              keytype = "SYMBOL") # information we have

# Check it worked!
head(anno)  
```

Now we can combine our annotation table to our marker genes using `dplyr` and `left_join`

```{r}
anno_markers <- markers %>% 
  dplyr::select(gene, cluster, avg_log2FC, p_val, p_val_adj) %>%
  left_join(anno, by = c("gene" = "SYMBOL"))

anno_markers %>%  dplyr::select(gene, cluster, avg_log2FC, p_val, p_val_adj) %>%
    group_by(cluster) %>% # important to look in each cluster!
    slice_max(n = 5, order_by = avg_log2FC)

```

At this point we might want to save our result into a nicely formatted file where we can look up our clusters, top DE genes and their gene descriptions.

```{r}
anno_markers %>% write_csv("Annotated.results.csv")
```

# Assigning cluster identity

If you recall, our dataset is from the following paper:-

> 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 May15;10(1):8079

This paper uses the following *canonical marker genes* to assign cell cluster identity. This starts to give some *biological meaning* to our results.

These genes are known to be expressed in a particular cell type.

Gene    Cell type
----    ----
CD79A   B Cells
CD34    CD34+, 
MS4A1   CD20+ B cells, 
MZB1    BCMA,
CST3    Monocytes, 
SPN / CD43   Immature hematopoietic, 
HBA1    erythrocytes, 
CD3D    T cells, 
NKG7    NK cells


We can assign these marker genes as a variable
```{r}
cell_markers <- c("CD79A","CD34", "MS4A1", "MZB1","CST3" ,"LSN", "KIT", "HBA1", "CD3D", "NKG7")
```


Dot plots allow us to examine particular features or genes of interest across all clusters. Creating a dotplot of the cell type markers lets us assign identity based on their expression.

```{r}
DotPlot(seurat, features = cell_markers) + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
```
We can then manually create new identities for each cluster as follows.

```{r}
new.cluster.ids <- c("T Cells", #0
                     "Erythrocytes", #1  
                     "Monocytes",  #2
                     "B Cells", #3
                     "CD34+", #4
                     "NK Cells", #5
                     "B Cells", #6
                     "Erythrocytes", #7
                     "BCMA",#8 
                     "Monocytes", #9
                     "CD20+ B cells", #10
                     "BCMA" #11
)

names(new.cluster.ids) <- levels(seurat)
new.cluster.ids
seurat <- RenameIdents(seurat, new.cluster.ids)
```

```{r}
DimPlot(seurat, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
```
We can check our result to the UMAP plot produced from the full dataset within the paper - *not too bad!* Considering we have used a subset and not followed their exact method.
![](images/umaps.png)

# Further explore cell type clusters

At this point we may want to interrogate our newly defined cell type clusters futher. We can use `FindMarkers` to look for specific genes that are differentially-expressed between clusters that represent biological conditions we might want to compare

```{r}
T_cell_B_cell_markers <- FindMarkers(seurat, ident.1 = "T Cells"  , ident.2 = "B Cells")

T_cell_B_cell_markers

```

Seurat has some nice in-built visualization tools to help us look at genes and features of interest within our cell type clusters. `VlnPlot` plots a violin plot showing the expression of a gene of interest, in this case SCIMP.

```{r}
VlnPlot(seurat, features="SCIMP")
```

<div class="exercise">

**Exercise**

1. Use `FindMarkers` to find differentially expressed marker genes between 2 other clusters of interest
2. Pick a gene at the top of this list and use the following visualisation tools to look at its expression pattern
 * `VlnPlot`
 * `RidgePlot`
 * `FeaturePlot`
 
</div>

```{r}
# Put exercise work in here


```

<div class="information">
Seurat has a great deal of customisability and plotting options. For more information look at the following vignette
https://satijalab.org/seurat/articles/visualization_vignette
</div>



# GO term and pathway enrichment

We often want to carry out a *Gene Enrichment Analysis* on a list of genes of interest. This enables us to find *functional terms* that are statistically over-represented within a group of genes compared to a background list. 

Seurat has in-built functionality to carry out a gene enrichment analysis using the `enrichr` package, which must be installed before carrying out this step. 

`enrichr` allows users to check a list of genes, for example the ones found in a cluster, to gene set libraries of known biological function. For more information see https://maayanlab.cloud/Enrichr/#.

We can check for GO terms measuring BP (biological process), CC (cell component) and MF (molecular function).


```{r}
#install.packages(enrichr)
library(enrichR)

results <- DEenrichRPlot(seurat, ident.1 = "Monocytes", enrich.database = "GO_Biological_Process_2023", max.genes = 1000, num.pathway = 10)

results
```



<div class="exercise">
**Final Exercise!**

There are a number of different gene set libraries available in `enrichr` https://maayanlab.cloud/Enrichr/#libraries.

1. Have a look for over representation of functional terms in a different cluster of interest (for example Erythrocytes) and check against a different database

</div>
```{r, message=FALSE}

# Put exercise work in here

results <- DEenrichRPlot(seurat, ident.1 = ....., enrich.database = ......, max.genes = 1000, num.pathway = 10)


```


```{r}

```



