---
title: "Pulldown proteomics DARPin"
author: "Cihan Erkut"
date: "20240116"
date-modified: "`r Sys.time()`"
format:
  html:
    theme: materia
    toc-title: Contents
    highlight-style: github
    df-print: tibble
    toc: true
    number-sections: true
    embed-resources: true
    code-line-numbers: true
    code-block-border-left: true
    code-fold: true
    code-tools: true
knitr:
  opts_chunk: 
    tidy: "styler"
execute:
  cache: true
---

## Setup

Load libraries and functions.

```{r}
#| label: setup

library(DEP) # BioConductor
library(SummarizedExperiment) # BioConductor
library(limma) # Bioconductor
library(ComplexHeatmap) # BioConductor
library(tidyverse) # CRAN
library(ggrepel) # CRAN
library(ggthemes) # CRAN
library(ggbreak) # CRAN
library(readxl) # CRAN
library(glue) # CRAN
library(here) # CRAN
```

```{r}
#| label: functions

import::here(map_uniprot_to_ensembl, .from = "functions.R")
import::here(prepare_data, .from = "functions.R")
import::here(plot_PCA_se, .from = "functions.R")
import::here(protein_dea, .from = "functions.R")
import::here(plot_PCA, .from = "functions.R")
```

## Load data

We load the data from the following analysis:

MG0075-3: T-DARPin pulldown in U-CH2
MG0075-4: T-DARPin pulldown in UM-Chor1
CU1045-4: Anti-TBXT pulldown in RIPA buffer in UM-Chor1

We also import ID mappings from UniProt.

```{r}
#| label: load data

table_samples_pulldown_UCH2 <- read_excel(
  path = here("raw/MG0075-3.xlsx"),
  sheet = "sampleNames"
) |>
  mutate(
    across(label, ~ str_remove(.x, "Sample")),
    cell = "UCH2"
  ) |>
  extract(
    col = condition,
    into = c("condition", "replicate"),
    regex = "(.*)_rep(\\d?)"
  ) |>
  unite(
    col = "condition",
    cell, condition
  )

table_pulldown_UCH2 <- read_excel(
  path = here("raw/MG0075-3.xlsx"),
  sheet = "filtered_proteinGroups"
) |>
  filter(
    !str_detect(`Protein IDs`, "DARPin")
  ) |>
  mutate(
    n = row_number(),
    experiment = "PD1"
  ) |>
  unite(
    col = "protein_group",
    experiment, n,
    sep = "_"
  )

table_samples_pulldown_UMChor1 <- read_excel(
  path = here("raw/MG0075-4.xlsx"),
  sheet = "sampleNames"
) |>
  mutate(
    across(label, ~ str_remove(.x, "Sample")),
    cell = "UMChor1"
  ) |>
  extract(
    col = condition,
    into = c("condition", "replicate"),
    regex = "(.*)_rep(\\d?)"
  ) |>
  unite(
    col = "condition",
    cell, condition
  )

table_pulldown_UMChor1 <- read_excel(
  path = here("raw/MG0075-4.xlsx"),
  sheet = "filtered_proteinGroups"
) |>
  filter(
    !str_detect(`Protein IDs`, "DARPin")
  ) |>
  mutate(
    n = row_number(),
    experiment = "PD2"
  ) |>
  unite(
    col = "protein_group",
    experiment, n,
    sep = "_"
  )

table_samples_pulldown_RIPA <- read_excel(
  path = here("raw/CU1045-4.xlsx"),
  sheet = "sampleNames"
) |>
  mutate(
    across(label, ~ str_remove(.x, "Sample")),
    buffer = "RIPA"
  ) |>
  extract(
    col = condition,
    into = c("condition", "replicate"),
    regex = "(.*)_rep(\\d?)"
  ) |>
  unite(
    col = "condition",
    buffer, condition
  )

table_pulldown_RIPA <- read_excel(
  path = here("raw/CU1045-4.xlsx"),
  sheet = "filtered_proteinGroups"
) |>
  mutate(
    n = row_number(),
    experiment = "PD4"
  ) |>
  unite(
    col = "protein_group",
    experiment, n,
    sep = "_"
  )

se_pulldown_UCH2 <- import_MaxQuant(
  proteins = table_pulldown_UCH2,
  expdesign = table_samples_pulldown_UCH2,
  filter = "Potential contaminant",
  names = "protein_group",
  ids = "protein_group"
)

se_pulldown_UMChor1 <- import_MaxQuant(
  proteins = table_pulldown_UMChor1,
  expdesign = table_samples_pulldown_UMChor1,
  filter = "Potential contaminant",
  names = "protein_group",
  ids = "protein_group"
)

se_pulldown_RIPA <- import_MaxQuant(
  proteins = table_pulldown_RIPA,
  expdesign = table_samples_pulldown_RIPA,
  filter = "Potential contaminant",
  names = "protein_group",
  ids = "protein_group"
)

table_annotation <- bind_rows(
  map_uniprot_to_ensembl(idmapping_file = here("data/HUMAN_9606_idmapping.dat"), se = se_pulldown_UCH2),
  map_uniprot_to_ensembl(idmapping_file = here("data/HUMAN_9606_idmapping.dat"), se = se_pulldown_UMChor1),
  map_uniprot_to_ensembl(idmapping_file = here("data/HUMAN_9606_idmapping.dat"), se = se_pulldown_RIPA)
)
```

## Batch effects

We test for batch effects attributable to a known factor and find replicate as the batch factor. Then we correct.

```{r}
#| label: test batch effects

table_pmap_args <- read_tsv(
  file = here("data/pulldown_dea_pmap_pars.tsv")
)

list_data_withBatch <- pmap(
  .l = table_pmap_args,
  .f = ~ prepare_data(get(..1), ..2, ..3, quantile_cutoff = 0.025)
) |>
  set_names(
    table_pmap_args$test_condition
  )

list_data <- pmap(
  .l = table_pmap_args,
  .f = ~ prepare_data(get(..1), ..2, ..3, batchCorrect = TRUE, quantile_cutoff = 0.025)
) |>
  set_names(
    table_pmap_args$test_condition
  )

map2(
  .x = list_data_withBatch,
  .y = names(list_data_withBatch),
  .f = \(x, y) plot_PCA_se(se = x, title = paste0(y, " without batch correction"))
)

map2(
  .x = list_data,
  .y = names(list_data),
  .f = \(x, y) plot_PCA_se(se = x, paste0(y, " with batch correction"))
)
```

## Results

Call results for each dataset and make a combined table.

```{r}
#| label: call results

table_results <- map_dfr(
  .x = list_data,
  .f = \(x) protein_dea(se = x, LFC_cutoff = log2(10), FDR = 0.05)
) |>
  rename(
    protein_group = Protein
  ) |>
  left_join(table_annotation) |>
  mutate(
    Condition = case_match(
      Condition,
      "UCH2_A2" ~ "U-CH2_A2",
      "UCH2_B1" ~ "U-CH2_B1",
      "UCH2_D4" ~ "U-CH2_D4",
      "UMChor1_A2" ~ "UM-Chor1_A2",
      "UMChor1_B1" ~ "UM-Chor1_B1",
      "UMChor1_D4" ~ "UM-Chor1_D4",
      "RIPA_TBXT" ~ "UM-Chor1_TBXT_RIPA"
    ),
    across(Condition, \(x) as_factor(x)),
    sig_bait = Significant & logFC > 0,
    ind = case_when(
      gene_names == "TBXT" ~ "TBXT",
      sig_bait ~ "Bait",
      !sig_bait ~ "Control"
    ),
    ind = fct_relevel(ind, "Control")
  ) |>
  select(Condition, gene_names, logFC, CI.L, CI.R, t, P.Value, adj.P.Val, Significant, ind, UniProtKB_AC, protein_group) |>
  rename(
    Genes = gene_names,
  ) |>
  filter(
    !is.na(Genes)
  ) |>
  distinct()
```

## Visualizations

```{r}
plot_volcano_all <- table_results |>
  ggplot(aes(
    x = logFC,
    y = -log10(P.Value),
    color = ind
  )) +
  geom_point(
    data = table_results[table_results$ind == "Control", ]
  ) +
  geom_point(
    data = table_results[table_results$ind == "Bait", ]
  ) +
  geom_point(
    data = table_results[table_results$ind == "TBXT", ]
  ) +
  geom_vline(
    xintercept = log2(10),
    linetype = "dotted"
  ) +
  scale_x_continuous(
    name = "log-fold change",
    limits = c(-11, 11)
  ) +
  scale_y_continuous(
    name = "-log10 p-value",
    limits = c(0, 20)
  ) +
  scale_color_manual(
    values = c(
      TBXT = "#FF0000AA",
      Bait = "#00000044",
      Control = "#AAAAAA11"
    )
  ) +
  facet_wrap(~Condition) +
  theme_few() +
  theme(
    aspect.ratio = 1,
    legend.position = "bottom",
    legend.title = element_blank()
  )

plot_bubble <- table_results |>
  filter(
    ind != "Control"
  ) |>
  mutate(
    across(Genes, ~ fct_relevel(fct_reorder(.x, logFC, max), "TBXT", after = Inf)),
    across(Condition, ~ fct_relevel(
      .x,
      "U-CH2_A2",
      "UM-Chor1_A2",
      "U-CH2_B1",
      "UM-Chor1_B1",
      "U-CH2_D4",
      "UM-Chor1_D4",
    ))
  ) |>
  ggplot(aes(
    x = 2^logFC,
    y = Genes,
    size = t,
    color = Condition
  )) +
  geom_point() +
  scale_x_continuous(
    name = "Fold change",
    limits = c(0, 2100),
    breaks = c(0, 100, 200, 2000, 2100)
  ) +
  scale_x_break(
    breaks = c(250, 2000),
    scales = c(0.4, 1)
  ) +
  scale_color_brewer(
    type = "qual",
    palette = "Paired"
  ) +
  scale_size(
    name = "moderated\nt-statistic"
  ) +
  guides(
    color = guide_legend(
      override.aes = list(
        size = 3
      )
    )
  ) +
  theme(
    axis.title.y = element_blank(),
    panel.background = element_rect(
      fill = "gray95"
    )
  )

plot_volcano_all

plot_bubble
```

## Save output

```{r}
write_tsv(
  x = table_results,
  file = here("data/table_results_apms_full.tsv")
)

ggsave(
  filename = here("plots/20250430_plot_volcano_all.pdf"),
  plot = plot_volcano_all,
  width = 18,
  height = 22,
  units = "cm"
)

ggsave(
  filename = here("plots/20250514_plot_bubble.pdf"),
  plot = plot_bubble,
  width = 20,
  height = 10,
  units = "cm"
)
```

