lab_pca_hmap.Rmd

---
title: "PCA and Heatmap"
subtitle: "Workshop on ggplot"
author: "`r paste0('<b>Lokesh Mano</b> • ',format(Sys.time(), '%d-%b-%Y'))`"
output:
  bookdown::html_document2:
          toc: true
          toc_float: true
          toc_depth: 4
          number_sections: true
          theme: flatly
          highlight: tango
          df_print: default
          code_folding: "none"
          self_contained: false
          keep_md: false
          encoding: 'UTF-8'
          css: "assets/lab.css"
---

```{r, include=FALSE}
hooks = knitr::knit_hooks$get()
hook_foldable = function(type) {
  force(type)
  function(x, options) {
    res = hooks[[type]](x, options)
    
    if (isFALSE(options[[paste0("fold.", type)]])) return(res)
    
    paste0(
      "<details><summary>", type, "</summary>\n\n",
      res,
      "\n\n</details>"
    )
  }
}
knitr::knit_hooks$set(
  output = hook_foldable("output"),
  plot = hook_foldable("plot")
)
```

```{r,child="assets/header-lab.Rmd"}
```

```{r, echo=FALSE, warning=FALSE, message=FALSE}
library(dplyr)
library(ggplot2)
library(reshape2)
data("iris")
md <- read.table("data/metadata_raw.csv", header = T, sep = ";")
rownames(md) <- md$Sample_ID

library(tidyverse)
```

# PCA

Let us first make a PCA object. For this, we will use the `VST` data, because it makes sense to use the normalized data for building the PCA. To run PCA, we use the R function `prcomp()`. It takes in a matrix where samples are rows and variables are columns. Therefore, we transpose our count matrix using the function `t()`. If we do not transpose, then PCA is run on the genes rather than the samples. 

```{r, warning=FALSE, message=FALSE}
gc_vst <- read.table("data/counts_vst.txt", header = T, row.names = 1, sep = "\t")
vst_pca <- prcomp(t(gc_vst))
```

After you computer the PCA, if you type the object `vst_pca$` and press `TAB`, you will notice that this R object has multiple vecors and data.frames within it. Some of the important ones are

* `sdev:` the standard deviations of the principal components
* `x:` the coordinates of the samples (observations) on the principal components.
* `rotation:` the matrix of variable loadings (columns are eigenvectors).

<i class="fas fa-exclamation-circle"></i> Note &nbsp; There are quite a few functions in R from different packages that can run PCA. So, one should look into the structure of the PCA object and import it into `ggplot` accordingly!

## Variance of components (Scree plot)

First, let us look into plotting the variance explained by the top PCs.

```{r, warning=FALSE, message=FALSE}
frac_var <- function(x) x^2/sum(x^2)
library(scales)

vst_pca$sdev %>% 
  as_tibble() %>% 
  frac_var() %>% 
  mutate(Comp = colnames(vst_pca$x)) %>% 
  slice(1:9) %>% 
  ggplot(aes(x=Comp, y = value)) + 
  geom_bar(stat = "identity", fill = "#4DC5F9") +
  geom_hline(yintercept = 0.03, linetype=2) +
  xlab("Principal Components") +
  scale_y_continuous(name = "Variance Explained", breaks = seq(0,0.8,0.1), labels = percent_format(accuracy = 5L)) +
  theme_classic(base_size = 14)

```

## PCA plot

So, looks like the first two components explain almost 85% of the data. Now, let us look into building a plot out of these components. From the above object, to get the scatter plot for the samples, we need to look into `vst_pca$x`. Then, we combine this data (as shown below) with the metadata to use different aesthetics and colors on the plot.

```{r, warning=FALSE, message=FALSE}
vst_pca$x
```

And, if you check the `class()` of this object, you will realize that this is a `matrix`. To be able to comfortably use `tidyverse` on this object, we must first convert this to a `data.frame`. 

```{r, warning=FALSE, message=FALSE}
vst_pca_all <- vst_pca$x %>%
  as.data.frame() %>%
  rownames_to_column(var = "Sample_ID") %>%
  full_join(md, by = "Sample_ID")

# Just to keep the order the right way.
vst_pca_all$Sample_Name <- factor(vst_pca_all$Sample_Name, levels = c("t0_A","t0_B","t0_C","t2_A","t2_B","t2_C","t6_A","t6_B","t6_C","t24_A","t24_B","t24_C"))
vst_pca_all$Time <- factor(vst_pca_all$Time, levels = c("t0","t2","t6","t24"))
vst_pca_all$Replicate <- factor(vst_pca_all$Replicate, levels = c("A","B","C"))

ggplot(vst_pca_all, aes(x=PC1, y=PC2, color = Time)) +
  geom_point(size = 3, aes(shape = Replicate)) +
  geom_vline(xintercept = 0, linetype=2) +
  geom_hline(yintercept = 0, linetype=2) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) 
```



## Loading plot

Now, let us say you want to plot the four genes that contribute the most to the four directions in the PCA plot. We could obtain them by looking at the `vst_pca$rotation` matrix. We could get those genes and their respective coordinates as follows.

```{r, message=FALSE, warning=FALSE}
genes.selected=vst_pca$rotation[c(which.max(vst_pca$rotation[,"PC1"]), which.min(vst_pca$rotation[,"PC1"]), which.max(vst_pca$rotation[,"PC2"]), which.min(vst_pca$rotation[,"PC2"])),c("PC1","PC2")]
genes.selected <- genes.selected %>%
  as.data.frame() %>%
  rownames_to_column(var = "Gene_ID")
genes.selected
```

A loading plot shows how strongly each variable (gene) influences a principal component. As an example, we could plot the four genes we selected.

```{r, message=FALSE, warning=TRUE}
ggplot(genes.selected, aes(x=PC1, y=PC2)) +
  geom_point() +
  geom_segment(aes(xend=PC1, yend=PC2), x=0, y=0, color="Grey") +
  geom_label(aes(x=PC1, y=PC2, label=Gene_ID), size=2, vjust="outward") +
  theme_bw() +
  theme(panel.grid.major = element_blank(),  panel.grid.minor = element_blank()) 
```


# PCA bi-plot

By merging the PCA plot with the loadings plot one can create a so-called PCA bi-plot.

```{r, message=FALSE, warning=FALSE}
scale=500
ggplot(data=vst_pca_all, mapping=aes(x=PC1, y=PC2)) +
     geom_point(size = 3, aes(shape = Replicate, color = Time)) +
     geom_vline(xintercept = 0, linetype=2) +
     geom_hline(yintercept = 0, linetype=2) +
     geom_segment(data=genes.selected, mapping=aes(xend=scale*PC1,yend=scale*PC2), x=0, y=0, arrow=arrow(), color="grey") +
     geom_label(data=genes.selected, mapping=aes(x=scale*PC1,y=scale*PC2, label=Gene_ID), size=2, hjust="outward", vjust="outward") +
     theme_bw() +
     theme(panel.grid.major = element_blank(),  panel.grid.minor = element_blank())
     
```


<i class="fas fa-exclamation-circle"></i> Note &nbsp; Similarly, let us say, you have environmental variables (continous variables like pH and so on) from the same samples and you would like to see how they would fit this bi-plot. One can use `envfit()` function from the `vegan` package, this function would return both the `p-value` and the `coordinates` of each of the variables in your environmental matrix. Then you could subset the significant variables and plot them in the same way as above. 

# Heatmap

## pheatmap

For heatmap, let us look into `pheatmap` library which is not part of `ggplot`, but it is a well known package for building heatmaps. It contains a lot of internal aesthetics that you can add that are very informative and intuitive. Let us first start with making a correlation matrix and plot it.

```{r, warning=FALSE, message=FALSE}
vst_cor <- as.matrix(cor(gc_vst, method="spearman"))
library(pheatmap)

pheatmap(vst_cor,border_color=NA,annotation_col=md[,"Time",drop=F],
         annotation_row=md[,"Time",drop=F],annotation_legend=T)
```

<i class="fas fa-exclamation-circle"></i> Note &nbsp; If you get a `gpar()` error, this is due to the reason that the metadata table you provide here does not have the sample ids as `rownames`. `pheatmap` maps the meta-data to the count data using `rownames`  

Now, this is based on a correlation matrix where you have a very simple square matrix with samples and their correlations to each to each other. Now let us look into how we can make a heatmap from an expression dataset. 

For this we will use the dataset: `Time_t2_vs_t0.txt` here that basically has the list of genes that are differentially expressed in `t2` vs `t0` and we would like to visualize this in an heatmap. To be more precise, these are the top and bottom 200 (2 X 100) genes with `adjusted p-value < 0.01`. We use just these 200, for the  sake of visualization. Let us take the `t24 vs t0` here as these are most different based on our PCA.

<i class="fas fa-exclamation-circle"></i> Note &nbsp; The genes are already sorted in the file. First 100 genes with `negative log2FoldChange` and the next 100 with `positive log2FoldChange`

```{r}
diff_t2_vs_t0 <- read.table("data/Time_t24_vs_t0.txt", sep = "\t", header = TRUE, row.names = 1)
hmap_t2_t0 <- subset(gc_vst, rownames(gc_vst) %in% diff_t2_vs_t0$gene)

pheatmap(hmap_t2_t0, border_color=NA, annotation_col=md[,"Time",drop=F], annotation_row=NULL, annotation_legend=T)
```

Further, we can also customize using the `pheatmap` and also use the continous color options we talked about in the earlier sections.

```{r}
library(wesanderson)
md$Time <- factor(md$Time, levels = c("t0","t2","t6","t24"))

pheatmap(hmap_t2_t0, color = wes_palette("Moonrise3", 100, type = "continuous"), border_color=NA, annotation_col=md[,"Time",drop=F], annotation_row=NULL, annotation_legend=T, show_rownames = FALSE, fontsize = 14)
```

## ggplot

We can use `geom_tile()` from `ggplot` to make heatmaps. The `ggplot` based heatmaps are, as we have seen for other cases, much more easier to customize.

```{r, warning=FALSE, message=FALSE}
hmap_t2_t0_long <- hmap_t2_t0 %>%
  rownames_to_column(var = "Gene") %>%
  gather(Sample_ID, VST, -Gene) %>% 
  full_join(md, by = "Sample_ID")  

hmap_t2_t0_long$Sample_Name <- factor(hmap_t2_t0_long$Sample_Name, levels = 
                                        c("t0_A","t0_B","t0_C","t2_A","t2_B","t2_C","t6_A","t6_B","t6_C","t24_A","t24_B","t24_C"))
hmap_t2_t0_long$Time <- factor(hmap_t2_t0_long$Time, levels = c("t0","t2","t6","t24"))
hmap_t2_t0_long$Replicate <- factor(hmap_t2_t0_long$Replicate, levels = c("A","B","C"))
hmap_t2_t0_long$Gene <- factor(hmap_t2_t0_long$Gene, levels = row.names(hmap_t2_t0))

ggplot(hmap_t2_t0_long) +
  geom_tile(aes(x = Sample_Name, y = Gene, fill = VST)) +
  scale_fill_gradientn(colours = rainbow(5)) +
  scale_x_discrete(limits = c("t0_A","t0_B","t0_C","t24_A","t24_B","t24_C")) +
  theme(axis.text.y = element_blank(), axis.ticks = element_blank()) 

```

# Exercise

Now, as I have mentioned earlier, building a plot similar to PCA really depends on how the object looks like. Now, let us try to make a `MDS` or `PCoA` plot from the same data as we have used. Here is how you get the MDS object in R.

```{r, message=FALSE, warning=FALSE}
gc_dist <- dist(t(gc_vst))
gc_mds <- cmdscale(gc_dist,eig=TRUE, k=2) 
```

<i class="fas fa-clipboard-list"></i> Task &nbsp; Now, try to replicate the example `MDS` plot below if you have enough time.

```{r, plot.fold=FALSE, echo=FALSE, warning=FALSE, message=FALSE}
Eigenvalues <- gc_mds$eig
Variance <- Eigenvalues / sum(Eigenvalues) 
Variance1 <- 100 * signif(Variance[1], 3)
Variance2 <- 100 * signif(Variance[2], 3)

gc_mds_long <- gc_mds$points %>%
  as.data.frame() %>%
  rownames_to_column("Sample_ID") %>%
  full_join(md, by = "Sample_ID")

gc_mds_long$Sample_Name <- factor(gc_mds_long$Sample_Name, levels = c("t0_A","t0_B","t0_C","t2_A","t2_B","t2_C","t6_A","t6_B","t6_C","t24_A","t24_B","t24_C"))
gc_mds_long$Time <- factor(gc_mds_long$Time, levels = c("t0","t2","t6","t24"))
gc_mds_long$Replicate <- factor(gc_mds_long$Replicate, levels = c("A","B","C"))

ggplot(gc_mds_long, aes(x=V1, y=V2, color = Time)) +
  geom_point(size = 3, aes(shape = Replicate)) +
  xlab(paste("PCO1: ", Variance1, "%")) +
  ylab(paste("PCO2: ", Variance2, "%")) +
  geom_vline(xintercept = 0, linetype=2) +
  geom_hline(yintercept = 0, linetype=2) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) 
```

<i class="fas fa-clipboard-list"></i> Hint &nbsp; The mds object `gc_mds` has "eigenvalues". You can calculate the variance by `Variance <- Eigenvalues / sum(Eigenvalues)`

# Session info

```{r, fold.output=FALSE, fold.plot=FALSE}
sessionInfo()
```

__End of document__