
```{r setup, include=FALSE}
# load libraries
library(tidyverse)
library(conflicted)

# additional libraries for today
library(viridis)
library(magrittr)
# library(broom)
# library(pheatmap)

# resolve package conflicts
filter <- dplyr::filter
select <- dplyr::select

# configure knitr settings
knitr::opts_chunk$set(echo = TRUE, fig.width = 6, fig.height = 4)
```

## 2020-07-15 Clustering Methods in R 

### Setup

Install the `broom` package. This package summarizes model information, so it's going to be very useful to make the clustering data easier to work with, but it'll be especially useful when we review statistics.

```{r}
# install.packages('broom')
```

Install the `pheatmap` package. `ggplot2` doesn't really do heatmaps, so you need a different plotting package.

```{r}
# install.packages('pheatmap')
```

And load the libraries

```{r}
library(broom)
library(pheatmap)
```

---

To demonstrate all of the clustering methods, we'll use the `iris` dataset again. There's an underlying variable, Species, that explains some of the variation in the data. We'll see if the clustering methods can recover the three species.

```{r}
ggplot(iris, aes(x = Sepal.Width, y = Sepal.Length, color = Species)) +
  geom_point(size = 3) +
  labs(x = 'Sepal Width (cm)', y = 'Sepal Length (cm)') +
  theme_classic() +
  theme(axis.title = element_text(size = 14), legend.title = element_text(size = 14))
```

<br>

### kmeans

For kmeans clustering, we have to pick the number of clusters we want. All of clustering methods only work on numeric data, so you have to remove any categorical variables before running the clustering function.

```{r}
iris %>% select(-Species) %>% kmeans(4)
```

Notice that when we used the `kmeans()` function, we didn't get a nice table as the result. If we check what type of object it is with the `class()` function, it tells us it's class "kmeans". R unfortunately has a lot of special types that aren't data tables that you have to deal with

```{r}
iris %>% select(-Species) %>% kmeans(4) %>% class()

### more typical types
# data.frame for tables
iris %>% class()
# or a tibble, the tidyverse version of a data.frame
who %>% class()
# There's different classes for plots, which makes sense because they're pictures
class(ggplot(iris, aes(x = Sepal.Width)) + geom_density())
```

The `broom` package functions transforms special model objects like "kmeans" objects into tables that you can do further work with. The `itdy()` function transforms the model results into a table.

```{r}
iris %>% select(-Species) %>% kmeans(4) %>% tidy()
```

But we don't really care about the model parameters, what we want is the cluster assignments with the original data so we can plot the clusters. The `augment()` function does this

```{r}
iris %>% select(-Species) %>% kmeans(4) %>% augment(iris) -> iris_kmeans4
```

And let's visualize the clusters we just created!

```{r}
ggplot(iris_kmeans4, aes(x = Sepal.Width, 
                         y = Sepal.Length, 
                         color = .cluster)) +
  geom_point(size = 3) +
  labs(x = 'Sepal Width (cm)', 
       y = 'Sepal Length (cm)',
       color = 'cluster') +
  theme_classic() +
  theme(axis.title = element_text(size = 14), 
        legend.title = element_text(size = 14))
```

Notice that clusters can overlap. Also it looks like from the visualization that four is probably not the correct number of clusters.

#### Test Different Kmeans Cluster Numbers

Because you don't know the correct number of kmeans clusters, we test a bunch of them. By running many kmeans and picking the kmeans with the fewest number of clusters with the best fit, we can do a little better than just guessing.

```{r}
# drop the categorial columns from the data
iris %>% select(-Species) -> iris_num

### do a bunch of kmeans
# make a table of numbers of clusters
tibble(k = 2:15) %>% 
# group by those numbers so the next command is applied to all of them
  group_by(k) %>% 
# do() is a helper function to get things to work with dplyr; do kmeans for all those cluster numbers
  do(kclust = kmeans(iris_num, .$k)) %>% 
# remember you use glance to get the model parameters
  glance(kclust) -> kmeans_params  

# plot to see the inflection point and pick number of clusters 
kmeans_params %>%
  mutate(group = 1) %>%   # just do this (add a grouping variable) to make geom_line() happy
  ggplot(aes(x = as.factor(k), y = tot.withinss, group = group)) + 
    geom_point(size = 3) + 
    geom_line(size = 1) + 
    labs(x = 'Number of Clusters', y = 'Goodness of Fit \n (within cluster sum of squares)') +
    theme_classic() +
    theme(axis.title = element_text(size = 14))
```

Let's try some of the numbers of clusters where we see an inflection point from the plot above.

```{r}
### kmeans with three clusters
iris %>% select(-Species) %>% kmeans(3) %>% augment(iris) %>%
  ggplot(aes(x = Sepal.Width, y = Sepal.Length, 
             color = .cluster, shape = Species)) +
    geom_point(size = 3) +
    labs(x = 'Sepal Width (cm)', y = 'Sepal Length (cm)', color = 'cluster') +
    theme_classic() +
    theme(axis.title = element_text(size = 14), legend.title = element_text(size = 14))

### kmeans with six clusters
iris %>% select(-Species) %>% kmeans(6) %>% augment(iris) %>%
  ggplot(aes(x = Sepal.Width, y = Sepal.Length, 
             color = .cluster, shape = Species)) +
    geom_point(size = 3) +
    labs(x = 'Sepal Width (cm)', y = 'Sepal Length (cm)', color = 'cluster') +
    theme_classic() +
    theme(axis.title = element_text(size = 14), legend.title = element_text(size = 14))
```

Three clusters does a pretty good job of finding the species of the irises.

<br>

## Heatmaps

### The Basics

To plot a heatmap with `pheatmap()`, all that's needed is a completely numeric table. Row and column names will become labels in the heatmap.

```{r}
iris %>%
# drop the species column because it's not numeric
  select(-Species) %>%
# plot the heatmap
  pheatmap(.)
```

**NOTE:** The problem of rownames

```{r}
### rownames
iris %>% #class()
  as_tibble() %>% #class()
  mutate(sample_id = paste0('sample', 1:nrow(.))) %>%
  # as.data.frame() %>%
  column_to_rownames('sample_id') %>%
```

### Modifying Appearance

Before playing with the appearance, save the data as a modified table that's in the correct format for `pheatmap()`

```{r}
iris %>%
  select(-Species) -> iris_num
```

**Drop the row names** `pheatmap` tries to label each individual row. 99% of the time there are so many rows that this is unreadable and looks terrible. You can tell `pheatmap` not to label them.

```{r}
pheatmap(iris_num, show_rownames = F)
```

**Change the color palette** for the heatmap by supplying a different R color palette to the `color` argument.

```{r}
### use an existing color palette like viridis
pheatmap(iris_num, show_rownames = F, color = viridis(50))

### changing the number after the palette changes how many colors/breaks there are in the scale
pheatmap(iris_num, show_rownames = F, color = viridis(5))

### pick colors for high, medium, low
# colorRampPalette takes a list of colors and makes them into a continuous palette
pheatmap(iris_num, show_rownames = F, color = colorRampPalette(c("navy", "white", "firebrick3"))(50))
```

<br>

**Cut the heatmap** to emphasize regions of interest. The `cutree_*` argument cuts based on the hierarchical clustering at the level where the number of branches matches the number you supply.

```{r}
### cut by columns
pheatmap(iris_num, show_rownames = F, color = magma(50), cutree_cols = 3)

### cut by rows
pheatmap(iris_num, show_rownames = F, color = magma(50), cutree_rows = 2)

### both
pheatmap(iris_num, show_rownames = F, color = magma(50), cutree_cols = 2, cutree_rows = 2)
```

#### Add color annotation to row/columns

`pheatmap()` will take annotation in a separate table and add it to the heatmap. The table must be a dataframe because the rownames of the annotation table must match either the column names or the row names (whichever one you want to annotate) of the main heatmap table. We'll set up an annotation table for the columns in the chunk below.

```{r}
# make a data frame with columns listing what you want to annotate
data.frame(treat_type = rep(c('length', 'width'), 2), 
# make the row names of the table the same as the column names
           row.names = colnames(iris_num)) -> col_anno
```

Add the annotation onto the heatmap with the `annotation_col` argument.

```{r}
pheatmap(iris_num, show_rownames = F, color = magma(50), annotation_col = col_anno)
```

**Select your own annotation colors** by making a list where the annotation categories are assigned colors.

```{r}
anno_colors <- list(treat_type = c(length = 'cornsilk3', 
                                   width = 'midnightblue'))
```

Then supply the annotation color list to the `annotation_colors` argument in `pheatmap()`

```{r}
pheatmap(iris_num, 
         color = viridis(50), 
         annotation_col = col_anno, 
         annotation_colors = anno_colors,
         show_rownames = F)
```

<br>

#### Combine everything

```{r}
pheatmap(iris_num, 
         color = viridis(10),
         cutree_rows = 2,
         cutree_cols = 2,
         annotation_col = col_anno, 
         annotation_colors = anno_colors,
         show_rownames = F)
```

<br>

### Clustering options

When clustering is implemented, first a measure of similarity is calculated, then second clustering is applied to that measure of similarity. `pheatmap()` calls the functions `dist()` and `hclust()` under the hood and you can select arguments to them within `pheatmap()`

#### Change the distance calculation

You can change the distance calculated or calculate a correlation instead by supply different arguments to `clustering_distance_rows` and/or `clustering_distance_cols`. Check the documentation for both `pheatmap()` and `dist()` for all the options

```{r}
### default
pheatmap(iris_num, show_rownames = F, clustering_distance_rows = 'euclidean')

### correlation
pheatmap(iris_num, show_rownames = F, clustering_distance_rows = 'correlation')

### or any other option from dist(), like manhattan
pheatmap(iris_num, show_rownames = F, clustering_distance_rows = 'manhattan')
```

<br>

#### Change the clustering method

You can change the clustering method by supplying a different argument to `clustering_method` in `pheatmap()`. Check the documentation for both `pheatmap()` and `hclust()` for all the options

```{r}
### default
pheatmap(iris_num, show_rownames = F, clustering_method = 'complete')

### different populat method
pheatmap(iris_num, show_rownames = F, clustering_method = 'ward.D2')

### third option; see hclust() documentation for complete list
pheatmap(iris_num, show_rownames = F, clustering_method = 'average')
```

<br>

#### No clustering

If you don't want your heatmap clustered, you can set the `cluster_rows` and/or `cluster_cols` arguments to FALSE.

```{r}
pheatmap(iris_num, cluster_rows = FALSE)
```

#### Scaling

Not clustering, but in the same vein, data can be centered and scaled by either rows or columns. Scaling is helpful when the range of magnitudes in your data is very large, but it can also be used to make small differences look larger and more important.

```{r}
### default is no scaling
pheatmap(iris_num, scale = 'none')

### scale rows
pheatmap(iris_num, scale = 'row')

### scale columns
pheatmap(iris_num, scale = 'column')
```
<br>

#### A non-clustered heatmap with `geom_tile()`

`ggplot()` sort of has a heatmap geom, `geom_tile()`. However, it only creates blocks of color and doesn't cluster the rows and columns.

```{r}
iris %>%
  mutate(fake_gene = rep(paste0('gene', 1:50), 3)) %>%
  
ggplot(aes(x = Species, y = fake_gene, fill = Petal.Length)) + 
  geom_tile() +
  scale_fill_viridis() +
  theme_void() +
  theme(axis.text.x = element_text(angle = 45, vjust = 0.4),
        axis.text.y = element_blank())
```

#### Note of Caution

You can drastically change the appearance of a heatmap. The same data can appear totally different depending on the parameters selected as seen in the examples below.

```{r, fig.width = 5, fig.height = 3}
pheatmap(iris_num,
         annotation_col = col_anno,
         annotation_colors = anno_colors,
         cutree_rows = 2, 
         cutree_cols = 3,
         show_rownames = F)
```

```{r, fig.width = 5, fig.height = 3}
pheatmap(iris_num, 
         color = viridis(50), 
         scale = 'column',
         annotation_col = col_anno,
         annotation_colors = anno_colors,
         clustering_distance_rows = 'euclidean',
         clustering_distance_cols = 'canberra',
         cutree_rows = 3,
         show_rownames = F)
```

<br>

### Principal Component Analysis (PCA)

For PCA, no prep is necessary beyond removing any categorical variables. You calculate it using the `prcomp()` function

```{r}
iris %>% select(-Species) %>% prcomp()
```

There are multiple ways to view the output. Frequently when you look at code on the internet, people will use base R's `summary()` to sort of tidy up the PCA output. `summary()` does nicely show the percent of variation explained by each principal component (PC).

```{r}
iris %>% select(-Species) %>% prcomp() %>% summary()
```

All the `broom` functions work with PCAs as well! And just like with kmeans, you can conviently use `augment()` to add the information from the model back to the original table

```{r}
iris %>% select(-Species) %>% prcomp() %>% augment(iris) -> iris_pca
```

Let's plot the first two PCs and see how well we did at clustering by species.

```{r}
ggplot(iris_pca, aes(x = .fittedPC1, y = .fittedPC2, color = Species)) +
  geom_point(size = 3) +
  labs(x = 'PC1', y = 'PC2') +
  theme_classic() +
  theme(axis.title = element_text(size = 14), 
        legend.title = element_text(size = 14))
```

By default, `prcomp()` centers the PCA, but doesn't scale it. Scaling is good for when your measurements are in different scales for example centimeters and meters or meters and feet. However, if your data is already in the same scale, you should leave this set to `FALSE`. Let's try it with scaling in addition to centering.

```{r}
# with scaling
iris %>% select(-Species) %>% prcomp(center = TRUE, scale = TRUE) %>% augment(iris) -> iris_pca_scale
```

```{r}
# with scaling
ggplot(iris_pca_scale, aes(x = .fittedPC1, y = .fittedPC2,              color = Species)) +
  geom_point(size = 3) +
  labs(x = 'PC1', y = 'PC2') +
  theme_classic() +
  theme(axis.title = element_text(size = 14), legend.title = element_text(size = 14))
```


















