recap_lab08.qmd

---
title: "Recap Lab 8"
author: "Leon Eyrich Jessen"
format:
  revealjs:
    embed-resources: true
    theme: moon
    slide-number: c/t
    width: 1600
    height: 900
    mainfont: avenir
    logo: images/r4bds_logo_small.png
    footer: "R for Bio Data Science"
---


<!--# ---------------------------------------------------------------------- -->
<!--# SLIDE ---------------------------------------------------------------- -->
<!--# ---------------------------------------------------------------------- -->
# Tidy PCA Assignment

```{r}
#| echo: false
#| eval: true
#| message: false

library("tidyverse")
library("broom")
library("ggrepel")
```



<!--# ---------------------------------------------------------------------- -->
<!--# SLIDE ---------------------------------------------------------------- -->
<!--# ---------------------------------------------------------------------- -->
## About Data - I Choose

### The BLOSUM (BLOcks SUbstitution Matrix)

- A substitution matrix used extensively in protein science e.g. sequence alignment

- Contains log-odds scores indicating how likely a substitution is

- Derived from very conserved regions of protein families in the BLOCKS database

- Using blocks of aligned sequence segments

- Blocks are chosen such that they have less than x% similarity

- E.g. BLOSUM62 is derived from blocks with less than 62% similarity

**Note, only numbers computed from observed substitutions here! No information on physico-chemical features of amino acids!**

<!--# ---------------------------------------------------------------------- -->
<!--# SLIDE ---------------------------------------------------------------- -->
<!--# ---------------------------------------------------------------------- -->
## Get Data

### Download

```{r}
#| echo: true
#| eval: true
#| message: false
#| warning: false

base_url <- "https://raw.githubusercontent.com/"
file_url <- "rdpstaff/AlignmentTools/master/src/data/blosum62.txt"

bl62 <- read_table(
  file = str_c(base_url, file_url),
  comment = "#")
bl62
```



<!--# ---------------------------------------------------------------------- -->
<!--# SLIDE ---------------------------------------------------------------- -->
<!--# ---------------------------------------------------------------------- -->
## Wrangle Data

### Adjust Column Names

```{r}
#| echo: true
#| eval: true
#| message: false
#| warning: false

colnames(bl62) <- c("aa", colnames(bl62))
```

### Subset to Data on the 20 Proteinogenic AA

```{r}
#| echo: true
#| eval: true
#| message: false
#| warning: false

std_aa = c("A", "R", "N", "D", "C", "Q", "E", "G", "H", "I",
           "L", "K", "M", "F", "P", "S", "T", "W", "Y", "V")
bl62 <- bl62 |> 
  filter(aa %in% std_aa) |> # Work on rows
  select(aa, std_aa)        # Work on columns
bl62
```



<!--# ---------------------------------------------------------------------- -->
<!--# SLIDE ---------------------------------------------------------------- -->
<!--# ---------------------------------------------------------------------- -->
## Wrangle Data

### Save to File

```{r}
#| echo: true
#| eval: true
#| message: false
#| warning: false

write_tsv(x = bl62,
          file = "data/blosum62.tsv") # Very small file, so skipping gz
```



<!--# ---------------------------------------------------------------------- -->
<!--# SLIDE ---------------------------------------------------------------- -->
<!--# ---------------------------------------------------------------------- -->
## Visualise Data

:::: {.columns}

::: {.column width="40%"}

```{r}
#| echo: true
#| eval: true
#| message: false
#| warning: false

pl <- bl62 |>
  pivot_longer(
    cols = -aa,
    names_to = "aa2",
    values_to = "log_odds_score") |> 
  ggplot(aes(
    x = factor(aa,
               levels = std_aa),
    y = factor(aa2,
               levels = rev(std_aa)),
    fill = log_odds_score,
    label = log_odds_score)) +
  geom_tile() +
  geom_text(colour = "darkgrey",
            size = 7) +
  scale_fill_gradient2(low = "blue",
                       mid = "white",
                       high = "red",
                       midpoint = 0) +
  scale_x_discrete(position = "top") +
  coord_fixed() +
  theme_minimal(base_size = 22,
                base_family = "Avenir") +
  theme(legend.position = "none",
        plot.title = element_text(
          hjust = 0.5,
          vjust = -4)) +
  labs(
    x = "",
    y = "",
    title = str_c(
    "The BLOSUM62 Visualised ",
    "as a Heatmap"))
```

:::

::: {.column width="60%"}

```{r}
#| echo: false
#| eval: true
#| message: false
#| warning: false
#| fig-width: 8
#| fig-height: 8
#| fig-align: center

pl
```

:::

::::



<!--# ---------------------------------------------------------------------- -->
<!--# SLIDE ---------------------------------------------------------------- -->
<!--# ---------------------------------------------------------------------- -->
## Create a PCA object

```{r}
#| echo: true
#| eval: true
#| message: false
#| warning: false

bl62_pca <- bl62 |>
  select_if(is.numeric) |>
  prcomp(center = TRUE,
         scale. = TRUE)
bl62_pca |>
  str()
```

- "Complicated" model object $\rightarrow$ `broom` to the rescue!

<!--# ---------------------------------------------------------------------- -->
<!--# SLIDE ---------------------------------------------------------------- -->
<!--# ---------------------------------------------------------------------- -->
## Make a scree plot using broom to tidy


:::: {.columns}

::: {.column width="40%"}

```{r}
#| echo: true
#| eval: true
#| message: false
#| warning: false

bl62_pca |>
  tidy("pcs")
```

:::

::: {.column width="60%"}

```{r}
#| echo: true
#| eval: true
#| message: false
#| warning: false

bl62_pca |>
  tidy("pcs") |> 
  mutate(percent = percent * 100) |> 
  ggplot(aes(x = PC,
             y = percent)) +
  geom_hline(yintercept = 0) +
  geom_col(colour = "black",
           alpha = 0.5) +
  theme_bw(base_size = 20) +
  theme(panel.grid.minor.x = element_blank(),
        panel.grid.major.x = element_blank()) +
  labs(title = "Scree Plot of PCA of BLOSUM62")
```

:::

::::



<!--# ---------------------------------------------------------------------- -->
<!--# SLIDE ---------------------------------------------------------------- -->
<!--# ---------------------------------------------------------------------- -->
## Augment using broom

```{r}
#| echo: true
#| eval: true
#| message: false
#| warning: false

bl62_pca_aug <- bl62_pca |>
  augment(bl62)
bl62_pca_aug
```



<!--# ---------------------------------------------------------------------- -->
<!--# SLIDE ---------------------------------------------------------------- -->
<!--# ---------------------------------------------------------------------- -->
## Add some chemical classes


:::: {.columns}

::: {.column width="40%"}

```{r}
#| echo: true
#| eval: true
#| message: false
#| warning: false

get_chem_class <- function(x){
  chem_cols <- c(
    "A" = "Hydrophobic",
    "R" = "Basic",
    "N" = "Neutral",
    "D" = "Acidic",
    "C" = "Sulphur",
    "Q" = "Neutral",
    "E" = "Acidic",
    "G" = "Polar",
    "H" = "Basic",
    "I" = "Hydrophobic",
    "L" = "Hydrophobic",
    "K" = "Basic",
    "M" = "Sulphur",
    "F" = "Hydrophobic",
    "P" = "Hydrophobic",
    "S" = "Polar",
    "T" = "Polar",
    "W" = "Hydrophobic",
    "Y" = "Polar",
    "V" = "Hydrophobic")
  return( chem_cols[x] ) # Example of avoiding dependencies => shareable!
}
```

- This is a named vector, think dictionary (Yes `R` can do that too!)

:::

::: {.column width="60%"}

![](images/chemistry_beaker.gif){width="50%" fig-align="center"}

:::

::::


<!--# ---------------------------------------------------------------------- -->
<!--# SLIDE ---------------------------------------------------------------- -->
<!--# ---------------------------------------------------------------------- -->
## Add some chemical classes

- Note how we are using our own custom function inside a `dplyr` pipeline here

```{r}
#| echo: true
#| eval: true
#| message: false
#| warning: false

get_chem_class(x = c("A", "R", "D"))
```

```{r}
#| echo: true
#| eval: true
#| message: false
#| warning: false

bl62_pca_aug <- bl62_pca_aug |> 
  mutate(chem_class = get_chem_class(aa))

bl62_pca_aug |>
  select(aa, chem_class)
```



<!--# ---------------------------------------------------------------------- -->
<!--# SLIDE ---------------------------------------------------------------- -->
<!--# ---------------------------------------------------------------------- -->
## Plot the PCA Plot

```{r}
#| echo: true
#| eval: true
#| message: false
#| warning: false

pca_plot_axes_labels <- bl62_pca |>
  tidy("eigenvalues") |>
  mutate(lbl = str_c("PC", PC, ", VE = ", round(percent*100,2), "%")) |> 
  pull("lbl")
pca_plot_axes_labels
```



<!--# ---------------------------------------------------------------------- -->
<!--# SLIDE ---------------------------------------------------------------- -->
<!--# ---------------------------------------------------------------------- -->
## Plot the PCA Plot

:::: {.columns}

::: {.column width="40%"}

```{r}
#| echo: true
#| eval: true
#| message: false
#| warning: false

pca_plot <- bl62_pca_aug |> 
  ggplot(aes(x = .fittedPC1,
             y = .fittedPC2,
             label = aa,
             colour = chem_class,
             fill = chem_class)) +
  geom_vline(xintercept = 0) +
  geom_hline(yintercept = 0) +
  geom_point(shape = 21,
             size = 6,
             alpha = 0.5) +
  geom_text(colour = "black") +
  scale_fill_manual(
    values = c("red", "blue", "black",
               "purple", "green", "yellow")) +
  scale_colour_manual(
    values = c("red", "blue", "black",
               "purple", "green", "yellow")) +
  coord_fixed() +
  theme_bw(base_size = 20,
           base_family = "avenir") +
  labs(
    title = "PCA: Scores Plot of BLOSUM62",
    x = pluck(pca_plot_axes_labels, 1),
    y = pluck(pca_plot_axes_labels, 2),
    fill = "Chemistry",
    colour = "Chemistry",
    caption = "Up: Aromatic, down: Aliphatic, Left: Charged, Right: Hydrophobic")
```

:::

::: {.column width="60%"}
```{r}
#| echo: false
#| eval: true
#| message: false
#| warning: false
#| #| fig-align: center
#| fig-width: 10
#| fig-height: 8

pca_plot
```

:::

::::



<!--# ---------------------------------------------------------------------- -->
<!--# SLIDE ---------------------------------------------------------------- -->
<!--# ---------------------------------------------------------------------- -->
## That's about it for tidyverse

- We have been over a ton of materials - Well done!

- Now we will have two great labs on R-packages and Shiny

- And then the mini-symposium, with talks on how R for Bio Data Science is being applied in industry

<!--# ---------------------------------------------------------------------- -->
<!--# SLIDE ---------------------------------------------------------------- -->
<!--# ---------------------------------------------------------------------- -->
## That's about it for tidyverse

Lastly, a quick case story:

- *A well seasoned Professor was curious about the tidyverse, but reluctant! Over perhaps 6 months, I would answer the occasional over-the-coffee questions and point to some materials, still feeling some push-back. Then suddenly, one day said colleague came into my office and proclaimed: "I get it now! The penny has dropped! Once you get your head around it, this thing is a complete game changer!*

- It was easy to for me to utter: "Could not agree more!"

- So, get your "head around it" and you WILL HAVE...



<!--# ---------------------------------------------------------------------- -->
<!--# SLIDE ---------------------------------------------------------------- -->
<!--# ---------------------------------------------------------------------- -->
## That's about it for tidyverse

![](images/palpatine_unlimited_power.gif){width="60%" fig-align="center"}