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

```
# Session 3: Learning Objectives
- Perform normalisation and select for highly variable genes

# Session 4: 
- Explain the purpose of dimension reduction
- Understand the difference between linear and non-linear dimension reduction methods (pros/ cons)
- Perform different dimension reduction techniques

```
# Session 3
## Libraries
```{r libraries}
library(DropletUtils) 
library(SingleCellExperiment)
library(scater)
library(scuttle)
library(GenomicRanges)
library(scran)
library(dittoSeq) #BiocManager::install("dittoSeq")
library(uwot) #install.packages("uwot")
```

## Import data
```{r import}
dat_dir <- file.path("Data", "10x_chromium")

sce <- read10xCounts(file.path(dat_dir, "10k_PBMC_3p_nextgem_Chromium_X_intron_filtered_feature_bc_matrix.h5"))
```

## Setup (from Day 1)

```{r day1_sce}
# Rename sce$Sample to "10k_PBMC"
sce$Sample <- "10k_PBMC"

# Change rownames to gene symbols
rownames(sce) <- uniquifyFeatureNames(rowData(sce)$ID, rowData(sce)$Symbol)

# Identify mitochondria genes (hint: all human mitochondria genes start with "MT-"). Useful function: grep()
is.mito <- grep("^MT-", rownames(rowData(sce)))

stats <- perCellQCMetrics(sce, subsets=list(mito=is.mito))

# What other QC metrics could be useful? 
## Ribosomal genes

#Convert to stats object to dataframe for easy plotting
stats <- data.frame(stats)

# Utilise isOutlier to determine the threshold for high MT percentage, and what cells have high mito percentage.
high.mito <- isOutlier(stats$subsets_mito_percent, type="higher") #Threshold is 9.323
attr(high.mito, "thresholds")


# Calculate outliers for all the fields using function perCellQCFilters()
discard_reasons <- perCellQCFilters(stats, sub.fields=c("subsets_mito_percent"))

# Summarise how many cells are discarded for each reason
colSums(as.matrix(discard_reasons))

# Add QC data to colData of unfiltered using addPerCellQC()
sce <- addPerCellQC(sce, subsets=list(mito=is.mito))
sce$discard <- discard_reasons$discard


unfiltered <- sce

#Discared the unwanted cells
sce <- sce[,sce$discard == FALSE]

# Add colnames
colnames(sce) <- sce$Barcode
```

## Normalisation
There are many types of normalisation strategies. The overarching goal of normalisation is to remove systematic differences in sequencing coverage and prevent these technical differences from interfering with comparisons of cell expression profiles. Today we will focus on scaling normalisation. Other methods include: normalisation by spike-ins, normalisation by deconvolution (built on removal of composition biases in bulk RNA seq). 

References: 
- https://genomebiology.biomedcentral.com/articles/10.1186/gb-2010-11-10-r106 
- Spike-ins: https://genome.cshlp.org/content/27/11/1795.long
- Deconvolution (scran): https://genomebiology.biomedcentral.com/articles/10.1186/s13059-016-0947-7
- ScTransform: https://genomebiology.biomedcentral.com/articles/10.1186/s13059-019-1874-1
- Benchmarking paper: https://genomebiology.biomedcentral.com/articles/10.1186/s13059-020-02136-7

Scaling normalisation: 
Library size is defined as the total sum of counts across all genes for each cell, expected value assumed to scale with any cell-specific biases. 
The “library size factor” for each cell is then directly proportional to its library size where the proportionality constant is defined such that the mean size factor across all cells is equal to 1.

```{r normalisation}
# Calculate the library size factor using the function librarySizeFactors()
lib.sf.sce <- librarySizeFactors(sce)
summary(lib.sf.sce)

# Graph the distribution of size factors derived from the library size as a histogram.
hist(log10(lib.sf.sce), xlab="Log10[Size factor]", col='grey80')

# Question: What does the function logNormCounts do? 
?logNormCounts()
## Answer: 
sce <- logNormCounts(sce)
```

## Log-transformation
Question: What is more interesting - a gene that is expressed at an average count of 50 in cell type A and 10 in cell type B, or a gene that is expressed at an average count of 1100 in A and 1000 in B? How can log-transformation be useful in looking at these types of differences?

## Feature selection
In this step, we want to select genes that contain useful biological information and remove the genes that contain random noise. The easiest approach is taking the most variable genes based on their overall expression--assuming that genuine biological differences will have increased variation compared to genes that are affected by technical noise.  
```{r feature_selection}
# Question: What does the function modelGeneVar() do? 
?modelGeneVar()
# Answer: 
allfeatures <- modelGeneVar(sce)

# Visualizing the fit:
fit.dat <- metadata(allfeatures)

# Question: What does this graph tell you?
plot(fit.dat$mean, fit.dat$var, xlab="Mean of log-expression",
    ylab="Variance of log-expression")
curve(fit.dat$trend(x), col="dodgerblue", add=TRUE, lwd=2)

# Look at the top most variable genes in allfeatures
allfeatures[order(allfeatures$bio, decreasing=TRUE),] 

# Choose the top highly variable genes using getTopHVGs()
# Question: What are different ways to choose the "top" genes?
topn2000_features <- getTopHVGs(allfeatures, n=2000) #top 2000 genes
topprop10_features <- getTopHVGs(allfeatures, prop=0.1) #top 10% genes
```

## Dimension Reduction
Question: How would you isolate CD3+ CD4+ CCR7+ T cells in flow cytometry? 
Answer: 

Question: What if you are looking at 200 genes, and you don't know what genes define a particular subset? How would you define that subset?

Question? What does supervised/ manual gating vs unsupervised techniques mean?
```{r}
#To learn more about dittoseq package: https://bioconductor.org/packages/devel/bioc/vignettes/dittoSeq/inst/doc/dittoSeq.html

# Plot scatter plot with CD3 and CD4 gene expression
dittoScatterPlot(
    object = sce, assay.x = "logcounts", assay.y = "logcounts",
    x.var = "CD3E", y.var = "CD4")

# Plot any other genes of your choice
dittoScatterPlot(
    object = sce, assay.x = "logcounts", assay.y = "logcounts",
    x.var = "??", y.var = "??")
```

## Plotting a 2D graph into 1D
Question: How would you best represent a 2D graph in only one dimension?

## PCA example with non-biological data
Video regarding PCA: https://www.youtube.com/embed/FgakZw6K1QQ

```{r}
library("factoextra")

# Import data
data(mtcars)

# Question: What is mtcars dataset showing?
```

```{r}
# Run PCA
res.pca <- prcomp(mtcars, scale = TRUE)
```

```{r}
# Scree plot
## Question: What does this graph show?
## Question: How much explained variance does the first two PCA dimensions capture? 
## Answer: 
fviz_eig(res.pca)
```

```{r}
## Question: What does this graph show?
## Answer:
fviz_pca_ind(res.pca,
             repel = TRUE     # Avoid text overlapping
             )
```

```{r}
# Add what country the car manufacturer is based
mtcars.country <- c(rep("Japan", 3), rep("US",4), rep("Europe", 7),rep("US",3), "Europe", rep("Japan", 3), rep("US",4), rep("Europe", 3), "US", rep("Europe", 3))

# Colour and group the indivdual car types by country. 
# Question: What does this show? 
fviz_pca_ind(res.pca, 
             col.ind = mtcars.country, # Color by the quality of representation
             addEllipses = TRUE, ellipse.type = "confidence",
             repel = TRUE     # Avoid text overlapping
             )
```

```{r}
# PCA contribution plot
## Question: What does this graph show? 
fviz_contrib(res.pca, choice = "var", axes = 1, top = 10)

## Question: What is the top contributing factor to PCA 2?
```

Variable correlation plots

## Question: What does this graph show? 
- Positively correlated variables are grouped together.
- Negatively correlated variables are positioned on opposite sides of the plot origin (opposed quadrants).
- The distance between variables and the origin measures the quality of the variables on the factor map. Variables that are away from the origin are well represented on the factor map.

## Question: What do the colours mean? 
The cos2 values are used to estimate the quality of the representation
The closer a variable is to the circle of correlations, the better its representation on the factor map (and the more important it is to interpret these components)
Variables that are closed to the center of the plot are less important for the first components.

```{r}

fviz_pca_var(res.pca,
             col.var = "cos2", # Color by contributions to the PC
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel = TRUE     # Avoid text overlapping
             )
```


```{r pca-swissroll}
#Reference for swissroll dataset:
#install.packages("remotes")
#remotes::install_github("jlmelville/snedata")
#library(snedata)
#swiss1000 <- swiss_roll(n = 1000)

# Import swiss roll dataset
swiss1000 <- read.csv(file = file.path("Data", "swiss1000.csv"))

# Plot dataset as scatterplot
## Question: Why is this dataset described as non-linear?
## Answer: 
plot(swiss1000$x, swiss1000$y, col = swiss1000$color, pch = 20)

## Run PCA
res.pca <- prcomp(swiss1000[,1:3], scale = TRUE)

## Plot PCA results
p <- fviz_pca_ind(res.pca, geom.ind = "point")
dat <- data.frame(p$data, "color" = swiss1000$color)

plot(dat$x, dat$y, col = dat$color, pch = 20)
```

# Reflection questions about PCA: 
- What are the benefits of PCA? 
- What are the limitations of PCA?

### Bonus: 
The package "factoextra" also has other methods of dimension reduction such as correspondance analysis (CA) and multiple factor analysis (MFA). Explore more: https://rpkgs.datanovia.com/factoextra/index.html 

## UMAP with non-biological data
https://jef.works//blog/2022/01/19/exploring-umap-parameters/
```{r}
#install.packages("uwot")
library(uwot)

set.seed(100) #this can be any number, we just set the seed so that the entire module has the same graph
swiss_umap <- umap(swiss1000[,1:3], n_neighbors = 50, learning_rate = 0.5, init = "random")

plot(swiss_umap[,1], swiss_umap[,2], col = swiss1000$color, pch = 20)

```

## PCA with biological data
```{r pca-sce}
# After feature selection (choosing top 2000 features), run PCA
sce <- runPCA(sce, subset_row=topn2000_features)

# Visualise the amount of variance captured per PC (scree plot)
percent.var <- attr(reducedDim(sce), "percentVar")
plot(percent.var, log="y", xlab="PC", ylab="Variance explained (%)")
```

```{r pca-sce-vis}
# Visualising PCAs
plotReducedDim(sce, dimred="PCA", colour_by="Sample") 

#Try visualizing the PCA and colouring by genes as well (ie. CD3E for T cells, CD19 for B cells)
plotReducedDim()

```
As demonstrated using the swiss roll dataset, PCA is a linear technique. This means that only variation along a line in the higher dimension space can be captured in each PC. Thus, PCA cannot capture all the variation in the n genes (n dimensions) using only the first two PCs. One workaround is to plot all of the top PCs together; however, this makes it very difficult to interpret and even then may not separate all the expected cellular sub populations. 
```{r pca2-sce-vis}
# Plot the top 4 PCs compared to the top 4 PCs 
## Try any genes of your choice
plotReducedDim(sce, dimred="PCA", ncomponents=4, colour_by="CD3E")
```
With any large datasets (biological or not), there are always a large number of biases/ confounders/ artifacts. One of the major challenges is identifying biological vs technical variability. It is typically near impossible to carry out true technical replicates, so computational methods are required to identify technical artifacts (and remove/ account for them). 
```{r}
# extract correlations between different variables and our PC scores. Typically, plotting individuals is a very helpful metric in a full dataset. 
explan_pcs <- getExplanatoryPCs(sce, variables = c("sum", "detected", "subsets_mito_percent")) 

plotExplanatoryPCs(explan_pcs/100)
```

When fitting a linear model regressing expression values for each gene against just that variable, we can calculate the marginal R2 for each variable. By plotting a density plot of the gene-wise marginal R2 values for the variables, we can determine if technical covariates (ie. which hospital the sample was collected, sequencing batch) have high correlations or not. If they do, we may need to repeat the normalisation step while accounting for these covariates or adjust down in downstream analysis. 

```{r}
plotExplanatoryVariables(sce, variables = c("sum", "detected", "subsets_mito_percent"))
```

## UMAP with biological data
UMAP is much faster than tSNE (see below); however, we still run it on the top PCs from PCA for consistency. 

UMAP has multiple parameters that will affect the visualisation. For more discussion details regarding the parameters, please look at: https://umap-learn.readthedocs.io/en/latest/parameters.html. For more details on the R implementation of UMAP, please look at: https://github.com/jlmelville/uwot.

The two parameters that have the most affect on the granularity are: 
- n_neighbors: number of neighbours
- min_dist: minimum distance between two embedded points

Try varying these parameters and determine what affect they have on the UMAP. 
```{r umap-sce}
# Run UMAP on these default settings and visualise it
set.seed(100) #any random number works, we just set a seed to keep it consistent in the module
sce <- runUMAP(sce, dimred="PCA", n_neighbors = 15, min_dist = 0.01)
out1 <- plotReducedDim(sce, dimred="UMAP", colour_by="CD3E") + ggtitle("n_neighbour=15, min_dist=0.01")
out1

# Question: What happens when you make n_neighbors small?
set.seed(100) #any random number works, we just set a seed to keep it consistent in the module
sce <- runUMAP(sce, dimred="PCA", n_neighbors = ??, min_dist = 0.01)
out2 <- plotReducedDim(sce, dimred="UMAP", colour_by="CD3E") + ggtitle("n_neighbour=??, min_dist=0.01")
out2

# Question: What happens when you make n_neighbors big?
set.seed(100) #any random number works, we just set a seed to keep it consistent in the module
sce <- runUMAP(sce, dimred="PCA", n_neighbors = ??, min_dist = 0.01)
out3 <- plotReducedDim(sce, dimred="UMAP", colour_by="CD3E") + ggtitle("n_neighbour=??, min_dist=0.01")
out3

# Question: What happens when you make min_dist very very small?
set.seed(100) #any random number works, we just set a seed to keep it consistent in the module
sce <- runUMAP(sce, dimred="PCA", n_neighbors = 15, min_dist = ??)
out4 <- plotReducedDim(sce, dimred="UMAP", colour_by="CD3E") + ggtitle("n_neighbour=15, min_dist=??")
out4

# Question: What happens when you make min_dist bigger?
set.seed(100) #any random number works, we just set a seed to keep it consistent in the module
sce <- runUMAP(sce, dimred="PCA", n_neighbors = 15, min_dist = ??)
out5 <- plotReducedDim(sce, dimred="UMAP", colour_by="CD3E") + ggtitle("n_neighbour=15, min_dist=??")
out5

# Question: What happens when you change both n_neighbors and min_dist?
set.seed(100) #any random number works, we just set a seed to keep it consistent in the module
sce <- runUMAP(sce, dimred="PCA", n_neighbors = ??, min_dist = ??)
out6 <- plotReducedDim(sce, dimred="UMAP", colour_by="CD3E") + ggtitle("n_neighbour=??, min_dist=??")
out6


gridExtra::grid.arrange(out1, out2, out3, out4, out5, out6, ncol= 3)
```

# Save progress
```{r}
saveRDS(sce, file.path("Data", "10kPBMC_day2_sce.rds"))
```

# Bonus: More dimension reduction techniques
tSNE
- tSNE: Van der Maaten, L., and G. Hinton. 2008. “Visualizing Data Using T-SNE.” J. Mach. Learn. Res. 9 (2579-2605): 85. 

This method is computationally more intensive. There are multiple parameters you can change to affect how the dimension reduction works. For more details regarding these parameters, please take a look at the following article: https://distill.pub/2016/misread-tsne/ 
```{r bonus-tSNE}
# Example of tSNE on SCE data 
set.seed(100)
sce <- runTSNE(sce, dimred="PCA", perplexity=5)
out7 <- plotReducedDim(sce, dimred="TSNE", colour_by="CD3E") + ggtitle("perplexity = 5")

set.seed(100)
sce <- runTSNE(sce, dimred="PCA", perplexity=20)
out8 <- plotReducedDim(sce, dimred="TSNE", colour_by="CD3E") + ggtitle("perplexity = 20")

set.seed(100)
sce <- runTSNE(sce, dimred="PCA", perplexity=80)
out9 <- plotReducedDim(sce, dimred="TSNE", colour_by="CD3E") + ggtitle("perplexity = 80")

gridExtra::grid.arrange(out7, out8, out9, ncol=3)
```

Correspondance Analysis
Try CA using this package: https://www.bioconductor.org/packages/release/bioc/html/corral.html 
Vignette: https://rdrr.io/bioc/corral/f/vignettes/corralm_alignment.Rmd