---
title: "Total 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(ComplexHeatmap) # BioConductor
library(DEP) # BioConductor
library(SummarizedExperiment) # BioConductor
library(limma) # Bioconductor
library(tidyverse) # CRAN
library(ggrepel) # CRAN
library(ggthemes) # CRAN
library(readxl) # CRAN
library(here) # CRAN
library(UpSetR) # 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")
```

## Load data

Load the data from the analysis CU0519-2. Also load the ID mappings that was downloaded from UniProt.

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

table_samples_total <- read_excel(
  path = here("raw/CU0519-2.xlsx"),
  sheet = "sampleNames"
) |>
  separate(
    col = condition,
    into = c("condition", "replicate"),
    sep = "_"
  ) |>
  mutate(across(label, ~ str_replace(sprintf("%02s", .x), " ", "0")))

table_total <- read_excel(
  path = here("raw/CU0519-2.xlsx"),
  sheet = "filtered_proteinGroups"
) |>
  filter(
    !str_detect(`Protein IDs`, "Darpin")
  ) |>
  mutate(
    n = row_number(),
    experiment = "TP1"
  ) |>
  unite(
    col = "protein_group",
    experiment, n,
    sep = "_"
  )

se_total <- import_MaxQuant(
  proteins = table_total,
  expdesign = table_samples_total,
  filter = "Potential contaminant",
  names = "protein_group",
  ids = "protein_group"
)

table_annotation <- map_uniprot_to_ensembl(
  idmapping_file = here("data/HUMAN_9606_idmapping.dat"),
  se = se_total
)
```

## Batch effects

There are no batch effects attributable to a known factor. We would rather not attempt to correct.

```{r}
#| label: check batch effects
#| fig-width: 6

prepare_data(
  se = se_total[, 1:12],
  batchDetect = TRUE,
  batchCorrect = FALSE
) |>
  plot_PCA_se(
    title = "Before batch correction"
  )

prepare_data(
  se = se_total[, 1:12],
  batchDetect = TRUE,
  batchCorrect = TRUE
) |>
  plot_PCA_se(
    title = "After batch correction"
  )

prepare_data(
  se = se_total[, 13:18],
  batchDetect = TRUE,
  batchCorrect = FALSE
) |>
  plot_PCA_se(
    title = "Before batch correction"
  )

prepare_data(
  se = se_total[, 13:18],
  batchDetect = TRUE,
  batchCorrect = TRUE,
) |>
  plot_PCA_se(
    title = "After batch correction"
  )
```

## Prepare data

Prepare the data for differential abundance analysis.

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

list_data <- list(
  total_A2 = prepare_data(se_total, "dNTC", "dA2"),
  total_B1 = prepare_data(se_total, "dNTC", "dB1"),
  total_D4 = prepare_data(se_total, "dNTC", "dD4"),
  total_KO = prepare_data(se_total, "sgNTC", "sgBrachyury")
)
```

## Differential abundance analysis

Mathematically, we define the int_reg function as such:

$$
\text{int_reg}_{Protein}(adj\_p, \text{logFC}) = \begin{cases} 
\text{sign}(\text{logFC}) & \text{if } adj\_p < 0.05 \text{ and } |\text{logFC}| > \log_2(1.15) \\
0 & \text{otherwise} 
\end{cases}

$$
For each protein, we sum up the int_reg across DARPin samples. If the absolute value of the sum is 2 or 3, we denote that protein as "commonly deregulated".

$$
\text{common}_{Protein}(\text{int_reg}) = \begin{cases} 
\text{TRUE} & \text{if } \sum_{Sample} |\text{int_reg}_{Protein,Sample}| > 1 \\
\text{FALSE} & \text{otherwise} 
\end{cases}

$$
Similarly, we define a consensus regulation direction for each protein. If the sum is greater than 1, it is "up", if less than -1, it is "down". This eliminates the 3 discordant cases that would normally be counted as "common" and receive a consensus regulation tag.

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

table_results <- map_dfr(
  .x = list_data,
  .f = \(x) protein_dea(se = x, LFC_cutoff = log2(1.15))
) |>
  mutate(
    across(Condition, \(x) str_remove(x, "d")),
    across(Condition, \(x) recode(x, "sgBrachyury" = "KO")),
    across(Condition, \(x) as_factor(x)),
    across(Protein, \(x) str_remove_all(x, "\\.\\d*$"))
  ) |>
  rename(
    protein_group = Protein
  ) |>
  left_join(table_annotation) |>
  select(Condition, gene_names, logFC, CI.L, CI.R, t, P.Value, adj.P.Val, Significant, UniProtKB_AC, protein_group) |>
  rename(
    Genes = gene_names,
  ) |>
  distinct() |>
  mutate(
    int_reg = Significant * sign(logFC),
    Regulation = case_match(
      int_reg,
      -1 ~ "Down",
      1 ~ "Up"
    )
  )

table_common_proteins <- table_results |>
  filter(
    Condition != "KO"
  ) |>
  summarise(
    median_DARPin_logFC = median(logFC, na.rm = TRUE),
    abs_sum_int_reg = abs(sum(int_reg)),
    DARPin_regulation = case_when(
      sum(int_reg) > 1 ~ "Up",
      sum(int_reg) < -1 ~ "Down"
    ),
    .by = protein_group
  ) |>
  mutate(
    Common = abs_sum_int_reg > 1
  )

table_results_common <- table_results |>
  left_join(
    table_common_proteins
  )
```

## Plots

Upset plot for set analysis.

```{r}
#| label: upset plot

treatments_darpin <- c("a2", "b1", "d4")
names(treatments_darpin) <- str_to_upper(treatments_darpin)

table_results_upset <- table_results |>
  filter(
    Significant,
    Condition != "KO"
  ) |>
  select(Condition, protein_group, Regulation) |>
  distinct()

list_deg <- map2(
  .x = rep(c("Up", "Down"), each = 3),
  .y = rep(names(treatments_darpin), 2),
  .f = \(x, y) unique(pull(filter(table_results_upset, Regulation == x & Condition == y), protein_group))
)

names(list_deg) <- paste(rep(names(treatments_darpin), 2), rep(c("Up", "Down"), each = 3))

pdf(
  file = here("plots/20240213_plot_upset_protein_darpin.pdf"),
  width = 7,
  height = 5
)

upset(
  data = fromList(list_deg),
  sets = rev(names(list_deg)),
  keep.order = TRUE,
  order.by = "freq",
  set_size.show = TRUE,
  set_size.scale_max = 1000
)

dev.off()
```

Volcano plots for the visualization of DEPs.

```{r}
#| label: volcano plots

table_highlight <- table_results |>
  filter(
    Genes %in% c("TBXT", "YAP1", "IGFBP3")
  )

table_p_cutoffs <- table_results |>
  mutate(
    tmp = adj.P.Val < 0.05
  ) |>
  summarise(
    p1 = max(P.Value),
    p2 = min(P.Value),
    .by = c(Condition, tmp)
  ) |>
  pivot_wider(
    names_from = tmp,
    values_from = c(p1, p2)
  ) |>
  reframe(
    x = -log10((p1_TRUE + p2_FALSE) / 2)
  )

plot_volcano_total <- table_results |>
  ggplot(aes(
    x = logFC,
    y = -log10(P.Value),
    color = Significant
  )) +
  geom_point(
    size = 0.2
  ) +
  geom_point(
    mapping = aes(
      color = Genes
    ),
    data = table_highlight,
  ) +
  geom_vline(
    xintercept = c(-1, 1) * log2(1.15),
    linetype = "dotted"
  ) +
  geom_hline(
    mapping = aes(
      yintercept = x
    ),
    data = table_p_cutoffs,
    linetype = "dotted"
  ) +
  scale_color_manual(
    values = c(
      `TRUE` = "black",
      `FALSE` = "gray70",
      "TBXT" = "red",
      "YAP1" = "blue",
      "IGFBP3" = "green"
    )
  ) +
  scale_x_continuous(
    name = "Fold change (log2)",
    limits = c(-1, 1) * max(abs(table_results$logFC)),
    breaks = seq(-6, 6, by = 3)
  ) +
  scale_y_continuous(
    name = "-log10(p-value)"
  ) +
  facet_wrap(~Condition, nrow = 2) +
  theme_few()

plot_volcano_total

ggsave(
  filename = here("plots/20240213_plot_volcano_tp1.pdf"),
  plot = plot_volcano_total,
  width = 15,
  height = 10,
  units = "cm"
)
```

## Save results

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

write_tsv(
  x = table_results_common,
  file = here("data/table_results_proteomics_tp1.tsv")
)
```

## Additional analysis for the revision

```{r}
table_results_KO_common <- table_results_common |>
  filter(
    Condition == "KO",
    Common,
    Significant
  )

with(
  table_results_KO_common,
  cor.test(
    x = logFC,
    y = median_DARPin_logFC,
    method = "spearman"
  )
)

plot_DEPs_in_KO <- table_results_KO_common |>
  ggplot(aes(
    x = median_DARPin_logFC,
    y = logFC
  )) +
  geom_hline(
    yintercept = 0,
    linetype = "dotted"
  ) +
  geom_vline(
    xintercept = 0,
    linetype = "dotted"
  ) +
  geom_point() +
  stat_smooth(
    method = "lm"
  ) +
  scale_x_continuous(
    name = expression("Median" ~ log[2] * "(fold change) (T-DARPins)")
  ) +
  scale_y_continuous(
    name = expression(Log[2] * "(fold change) (TBXT KO)")
  ) +
  annotate(
    geom = "text",
    x = 5,
    y = 9,
    label = "r = 0.87"
  ) +
  theme_few()

plot_DEPs_in_KO

ggsave(
  filename = here("plots/20250526_plot_DDA_DEPs_in_KO.pdf"),
  plot = plot_DEPs_in_KO,
  width = 12,
  height = 12,
  units = "cm"
)
```

