---
title: "Day 3 - MSc Genomics Medicine (Practical)"
output: 
  html_document:
    df_print: paged
---

```
Learning Objectives
- Explain the benefits and limitations of clustering
- Perform different clustering techniques
- Plot and identify top markers for each cluster
- Compare and contrast different annotation methods
      - Manual gating (flow cytometry style)
      - Manual cluster annotation using different cluster resolutions
      - Annotation using a reference
- Discuss challenges with automated annotation methods 

```
#### CHARLOTTE THIS VERSION
# Libraries
```{r libraries, echo=FALSE}
library(DropletUtils) 
library(SingleCellExperiment)
library(scater)
library(scuttle)
library(GenomicRanges)
library(scran)
library(dittoSeq) #BiocManager::install("dittoSeq")
library(uwot) #install.packages("uwot")
library(bluster) #BiocManager::install("bluster")
library(SingleR) #BiocManager::install("SingleR")
library(celldex)#BiocManager::install("celldex")
```

# Import

## Version 1: Day 1/2 condensed
```{r}
dat_dir <- file.path("Paste you file path here")
sce <- read10xCounts(file.path(dat_dir, "10k_PBMC_3p_nextgem_Chromium_X_intron_filtered_feature_bc_matrix.h5"))

rownames(sce) <- uniquifyFeatureNames(rowData(sce)$ID, rowData(sce)$Symbol)
colnames(sce) <- sce$Barcode

#set.seed(100)
#e.out <- emptyDrops(counts(sce))
#sce <- sce[,which(e.out$FDR <= 0.001)]

unfiltered <- sce

sce <- addPerCellQC(sce, subsets = list(Mito = grep("^MT-", rowData(sce)$Symbol)))

high_mito <- isOutlier(sce$subsets_Mito_percent, nmads = 3, type="higher")
low_detection <- (sce$detected < 1000)
high_counts <- sce$sum > 45000

data.frame(ByLibSize=sum(high_counts), ByFeature=sum(low_detection), ByMT=sum(high_mito), Remaining=ncol(sce))

sce <- sce[,!high_mito & !low_detection & !high_counts]
sce

# Removing lowly expressed genes
num_reads <- 1
num_cells <- 0.01*ncol(sce)
keep <- which(DelayedArray::rowSums(counts(sce) >= num_reads ) >= num_cells)
sce <- sce[keep,]
sce


set.seed(100)
#clusters <- quickCluster(sce)
#sce <- computeSumFactors(sce, cluster=clusters)
sce <- computeSumFactors(sce)
sce <- logNormCounts(sce)

set.seed(100)
dec.pbmc <- modelGeneVarByPoisson(sce)
#allfeatures <- modelGeneVar(sce)
top.pbmc <- getTopHVGs(dec.pbmc, prop=0.1)
#top.pbmc <- getTopHVGs(allfeatures, n=2000) #top 2000 genes

set.seed(100)
sce <- runPCA(sce, subset_row=top.pbmc)

set.seed(100)
sce <- runUMAP(sce, dimred="PCA")

plotReducedDim(sce, dimred="UMAP", colour_by="CD3E") 


```
```{r}

saveRDS(sce, file.path("Your path here", "10kPBMC_day2_sce.rds"))
```

## Version 2: load QC'd dataset directly
```{r}
sce <- readRDS(file.path("Your path here", "10kPBMC_day2_sce.rds"))
colnames(sce) <- sce$Barcode

# Double check the UMAP is what you expect
plotReducedDim(sce, dimred="UMAP", colour_by="CD3E")
```

# Clustering
## Graph-based clustering (ie. Louvain clustering)
Graph-based clustering works by taking each cell as a node, and building a graph where cells are connected to its nearest-neighbours in high-dimensional space by edges. Based on the similarity between the cells, edges are weighted such that closely related cells have a higher weight. Different algorithms are then applied to determine "clusters" of cells. 

Popularised by the Seurat framework, the biggest advantage of graph-based clustering is arguably its scalability. It also requires no strong assumptions about the cluster shape or cell distribution within the clusters (ie. normality).

The biggest drawback of graph-based clustering is that no information beyond a cell's nearest-neighbour is retained. Information such as data density is lost. This makes it difficult for algorithms which determines what a cluster is since high-density regions are more likely to cause formation of subclusters. In other words, resolution of clustering becomes dependent on the density of cells which can be misleading or hard to interpret. 

clusterCells() uses the 10 nearest neighbors of each cell to construct a shared nearest neighbor graph. Two cells are connected by an edge if any of their nearest neighbors are shared, with the edge weight defined from the highest average rank of the shared neighbors (Xu and Su 2015). The Louvain method is then used to identify communities.
```{r}
clust.louvain <- clusterCells(sce, use.dimred="PCA", 
                              BLUSPARAM=NNGraphParam(cluster.fun="louvain", 
                                                     cluster.args = list(resolution =1)))
sce$clusters_res1 <- clust.louvain

# Question: How many cells are in each cluster?
table(sce$clusters_res1)

# Plot the clusters on a UMAP
dittoDimPlot(sce, "clusters_res1", do.label = TRUE)

# Question: What happens when you modify the resolution?
clust.louvain <- clusterCells(sce, use.dimred="PCA", BLUSPARAM=NNGraphParam(cluster.fun="louvain", cluster.args = list(resolution =??)))

sce$clusters_res?? <- clust.louvain

```
### Directly comparing two different cluster resolutions using a heatmap (and using hierarchical clustering)
Hierarchical clustering does not scale well at all. 

```{r}
library(pheatmap)

clust.louvainres1 <- clusterCells(sce, use.dimred="PCA", BLUSPARAM=NNGraphParam(cluster.fun="louvain", cluster.args = list(resolution =1)))
clust.louvainres2 <- clusterCells(sce, use.dimred="PCA", BLUSPARAM=NNGraphParam(cluster.fun="louvain", cluster.args = list(resolution =2)))

tab <- table(paste("Louvain_res1", clust.louvainres1), 
    paste("Louvain_res2", clust.louvainres2))
ivw <- pheatmap(log10(tab+10), main="Louvain Res 1 vs Res 2",
    color=viridis::viridis(100), silent=TRUE)

gridExtra::grid.arrange(ivw[[4]])
```
### Another way to compare cluster resolutions
```{r}
clust1 <- clusterCells(sce, use.dimred="PCA", BLUSPARAM=NNGraphParam(cluster.fun="louvain", cluster.args = list(resolution =1)))
clust12 <- clusterCells(sce, use.dimred="PCA", BLUSPARAM=NNGraphParam(cluster.fun="louvain", cluster.args = list(resolution =1.2)))
clust16 <- clusterCells(sce, use.dimred="PCA", BLUSPARAM=NNGraphParam(cluster.fun="louvain", cluster.args = list(resolution =1.6)))

combined <- cbind(res.1=clust1, res.12=clust12, res.16=clust16)

library(clustree) #install.packages("clustree")
set.seed(1111)
clustree(combined, prefix="res.", edge_arrow=FALSE)
```
# Characterising clusters
Having identified your clusters, the next step is to try and assign a biological phenotype to each cluster, in order to characterise what cell types are in your dataset

This is typically done by unsupervised differential gene expression (looking for the highly expressed genes in each cluster compared to thes rest of the dataset),
and supervised characterisation using known marker genes based on the literature. 

This is often an iterative process whereby you cluster at different resolutions in order to best capture the biology heterogeneity in the dataset.

## Own knowledge - using known markers!
```{r}
key_genes <- c("CD3E", "CD4", "CD8A", "CD14", "FCGR3A", "CD19", "MS4A1", "NCAM1", "PPBP")
multi_dittoDimPlot(sce, key_genes)

plotExpression(sce, features=key_genes, x="clusters_res1", colour_by="clusters_res??")
```


## Find top genes
Find candidate marker genes for groups of cells (e.g., clusters) by testing for differential expression between pairs of groups.

Review the `?findMarkers` documentation to look at the different types of statistical tests available.

```{r}
markers <- findMarkers(sce, sce$clusters_res??, pval.type="some", direction="up")

# pull out the top markers for cluster 1
marker.set <- markers[["1"]]
as.data.frame(marker.set[1:30,1:3])
```


### Name your clusters
```{r}
manual_annotation <- data.frame(cluster_id = seq(1, 14),
                                cluster_name = c("Non-classical monocyte", "Classical monocyte", "Intermediate monocyte", "Memory_1 CD4+ T cell", "Memory B cell", "Naive B cell", "Naive CD4 + T cell", "Naive CD8+ T cell", "Memory CD8+ T cell", "Memory_2 CD4+ T cell", "NK", "doublet1", "doublet2", "ILC"))

sce$cell_subset <- sce$clusters_res1

temp_dat <- data.frame( "cluster_id" = sce$cell_subset, "barcode" =  sce$Barcode)

temp_dat <- merge(temp_dat, manual_annotation, by.x = "cluster_id", by.y = "cluster_id")
row.names(temp_dat) <- temp_dat$barcode
temp_dat = temp_dat[colnames(sce),]
head(temp_dat)

sce$cluster_name<- temp_dat$cluster_name


```

# Visualisation
## Dot plots
```{r}
key_genes <- c("CD3E", "CD4", "CD8A", "CCR7", "CD14", "FCGR3A", "CD19", "MS4A1", "CD27", "IGHD", "IGHM", "NCAM1", "PPBP")

dittoDotPlot(sce, vars = key_genes, group.by = "clusters_res1")

dittoDotPlot(sce, vars = key_genes, group.by = "cluster_name")

```


You could also use heatmaps to visualise this.



# Automated Annotation (singleR)
Annotation is one of the most time consuming and difficult aspects of single-cell analysis. Often times, different groups will decide to use different cluster resolutions and give their subsets slightly different names in order to answer their questions of interest (even on the same publicly available dataset). 

Open-access annotated datasets are an incredibly valuable resource for both wet-lab and dry-lab researchers. In 2022, Tabula Sapiens dataset was released and it highlights the collaborative efforts made by a huge global community. https://tabula-sapiens-portal.ds.czbiohub.org/ For those who are interested in looking at this resource, please visit CZI's CellxGene portal for this dataset: https://cellxgene.cziscience.com/e/c5d88abe-f23a-45fa-a534-788985e93dad.cxg/  

Review regarding automation annotation: https://www.nature.com/articles/s41596-021-00534-0
Tabular Sapiens paper: https://www.science.org/stoken/author-tokens/ST-495/full 

```{r}
ref <- celldex::MonacoImmuneData()
unique(ref$label.main) # these are the major cell types in the dataset
unique(ref$label.fine) # these are the higher resolution cell types in the dataset
anno.single <- SingleR(sce, ref = ref, labels = ref$label.main)
```

```{r}
# incorporate the singleR results into your sce object. 
head(anno.single)
temp_dat = data.frame(anno.single$pruned.labels, row.names=row.names(anno.single))
head(temp_dat)
# make sure temp_dat is the same order as our sce object
sce$single_annot = temp_dat[colnames(sce),]

plotReducedDim(sce, dimred="UMAP", colour_by="single_annot")
```

use the heatmap code that we used earlier to compare the singleR results, to our annotation.

```{r}

```

