---
title: "Total proteomics IGFBP3"
author: "Cihan Erkut"
date: "20240305"
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(tidyverse) # CRAN
library(glue) # CRAN
library(ggrepel) # CRAN
library(ggthemes) # CRAN
library(UpSetR) # CRAN
library(readxl) # 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

The following analysis were used:

CU0780-6: Data independent analysis (DIA) / sgNTC +/- D4
CU0780-7: DIA / sgIGFBP3 +/- D4

Also UniProt mappings.

```{r}
table_sgNTC <- read_excel(
  path = here("raw/CU0780-6.xlsx"),
  sheet = "filtered_proteinGroups"
) |>
  mutate(
    Potential.contaminant = if_else(PG.Organisms != "Homo sapiens", "+", NA)
  ) |>
  rename(
    Majority.protein.IDs = PG.ProteinGroups
  ) |>
  select(Majority.protein.IDs, Potential.contaminant, contains("Quantity")) |>
  rename_with(
    .fn = \(x) str_remove_all(x, "(\\[\\d*\\] oecf3-)|(.raw.PG.Quantity)")
  ) |>
  mutate(
    n = row_number(),
    experiment = "TP2"
  ) |>
  unite(
    col = "protein_group",
    experiment, n,
    sep = "_"
  )

table_samples_sgNTC <- read_excel(
  path = here("raw/CU0780-6.xlsx"),
  sheet = "Samples"
)

se_sgNTC <- import_MaxQuant(
  proteins = table_sgNTC,
  expdesign = table_samples_sgNTC,
  filter = "Potential.contaminant",
  intensities = "CU0780",
  names = "protein_group",
  ids = "protein_group"
)

table_sgIGFBP3 <- read_excel(
  path = here("raw/CU0780-7.xlsx"),
  sheet = "filtered_proteinGroups"
) |>
  mutate(
    Potential.contaminant = if_else(PG.Organisms != "Homo sapiens", "+", NA)
  ) |>
  rename(
    Majority.protein.IDs = PG.ProteinGroups
  ) |>
  select(Majority.protein.IDs, Potential.contaminant, contains("Quantity")) |>
  rename_with(
    .fn = \(x) str_remove_all(x, "(\\[\\d*\\] oecf3-)|(.raw.PG.Quantity)")
  ) |>
  mutate(
    n = row_number(),
    experiment = "TP3"
  ) |>
  unite(
    col = "protein_group",
    experiment, n,
    sep = "_"
  )

table_samples_sgIGFBP3 <- read_excel(
  path = here("raw/CU0780-7.xlsx"),
  sheet = "Samples"
)

se_sgIGFBP3 <- import_MaxQuant(
  proteins = table_sgIGFBP3,
  expdesign = table_samples_sgIGFBP3,
  filter = "Potential.contaminant",
  intensities = "CU0780",
  names = "protein_group",
  ids = "protein_group"
)

table_annotation <- bind_rows(
  map_uniprot_to_ensembl(idmapping_file = here("data/HUMAN_9606_idmapping.dat"), se = se_sgNTC),
  map_uniprot_to_ensembl(idmapping_file = here("data/HUMAN_9606_idmapping.dat"), se = se_sgIGFBP3)
)

```

## Batch effects

No batch effects attributable to a known factor could be found.

```{r}
list_data <- list(
  sgNTC = prepare_data(se_sgNTC, "dNTC", "dD4", dia = TRUE),
  sgIGFBP3 = prepare_data(se_sgIGFBP3, "dNTC", "dD4", dia = TRUE)
)

list_data_withoutBatch <- list(
  sgNTC = prepare_data(se_sgNTC, "dNTC", "dD4", dia = TRUE, batchCorrect = TRUE),
  sgIGFBP3 = prepare_data(se_sgIGFBP3, "dNTC", "dD4", dia = TRUE, batchCorrect = TRUE)
)

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

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

## Call results

We apply differential abundance analysis separately to each dataset but do not use it further, because we mainly focus on the interaction between IGFBP3 KO and TBXT inhibition.

```{r}
table_results <- map_dfr(
  .x = list_data,
  .f = \(x) protein_dea(se = x, LFC_cutoff = log2(1.15)),
  .id = "Condition"
) |>
  mutate(
    across(Condition, \(x) as_factor(x))
  ) |> 
  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() |> 
  filter(
    !is.na(Genes)
  )

table_results |> 
  ggplot(aes(
    x = P.Value
  )) +
  geom_histogram(
    bins = 100
  ) +
  facet_wrap(~ Condition)
```

Here, we create a new LFQ matrix from proteins that were detected in both analysis, so that it could be modelled for two factors.

`design_factorial` is the model matrix here.

```{r}
table_paired_proteins <- table_annotation |> 
  filter(
    protein_group %in% union(rownames(list_data$sgNTC), rownames(list_data$sgIGFBP3))
  ) |> 
  filter(
    n() == 2,
    .by = gene_names
  ) |> 
  mutate(
    tmp = if_else(str_detect(protein_group, "TP2"), "sgNTC", "sgIGFBP3")
  ) |> 
  pivot_wider(
    names_from = tmp,
    values_from = protein_group
  )

matrix_sgNTC <- assay(list_data$sgNTC[table_paired_proteins$sgNTC,])
matrix_sgIGFBP3 <- assay(list_data$sgIGFBP3[table_paired_proteins$sgIGFBP3,])

rownames(matrix_sgNTC) <- table_paired_proteins$gene_names
rownames(matrix_sgIGFBP3) <- table_paired_proteins$gene_names

colnames(matrix_sgNTC) <- paste0("sgNTC_", colnames(matrix_sgNTC))
colnames(matrix_sgIGFBP3) <- paste0("sgIGFBP3_", colnames(matrix_sgIGFBP3))

matrix_factorial <- cbind(
  matrix_sgNTC,
  matrix_sgIGFBP3
)

plot_PCA(
  source_matrix = matrix_factorial,
  title = "Combined"
)

table_samples_factorial <- tibble(
  Sample = colnames(matrix_factorial)
) |> 
  separate_wider_regex(
    cols = Sample,
    patterns = c(
      KD = "^.*?",
      "_",
      DARPin = ".*?",
      "_",
      replicate = "\\d*$"
    ),
    cols_remove = FALSE
  ) |> 
  mutate(
    across(KD, \(x) fct_relevel(as_factor(x), "sgNTC")),
    across(DARPin, \(x) fct_relevel(as_factor(x), "dNTC"))
  ) |> 
  column_to_rownames("Sample")

design_factorial <- model.matrix(~KD * DARPin, data = table_samples_factorial)

table_results_factorial <- lmFit(
  object = matrix_factorial,
  design = design_factorial
) |> 
  eBayes(
    trend = TRUE,
    robust = TRUE
  ) |> 
  topTable(
    coef = 4,
    number = Inf,
    sort.by = "t",
    confint = TRUE
  ) |> 
  as_tibble(
    rownames = "Genes"
  ) |> 
  mutate(
    Significant = adj.P.Val < 0.05
  )

table_hits_factorial <- table_results_factorial |> 
  filter(
    Significant
  ) |> 
  mutate(
    Direction = if_else(logFC > 0, "Up", "Down")
  )

table_results_combined <- table_results |> 
  left_join(
    table_hits_factorial,
    suffix = c("_Single", "_Factorial"),
    by = "Genes"
  )
```

And some quality control, some visualization.

```{r}
table_results_combined |> 
  ggplot(aes(
    x = logFC_Single,
    y = -log10(P.Value_Single),
    color = Direction
  )) +
  geom_point() +
  facet_wrap(~ Condition) +
  theme_few()


table_results_factorial |> 
  ggplot(aes(
    x = logFC,
    y = -log10(P.Value),
    color = Significant
  )) +
  geom_point()

table_results_factorial |> 
  ggplot(aes(
    x = P.Value
  )) +
  geom_histogram(
    binwidth = 0.01
  )

table_results_nonsig <- table_results |> 
  filter(
    !Significant
  )

table_results_sig <- table_results |> 
  filter(
    Significant
  )

table_results_nonsig |> 
  ggplot(aes(
    x = logFC,
    y = -log10(P.Value)
  )) +
  geom_point(
    color = "gray50",
    size = 0.3
  ) +
  geom_point(
    data = table_results_sig,
    color = "black",
    size = 0.3
  ) +
  scale_x_continuous(
    name = "Fold change (log2)",
    limits = c(-6.1, 6.1)
  ) +
  scale_y_continuous(
    name = "-log10(p-value)"
  ) +
  facet_wrap(~ Condition) +
  theme_few()
```

This is the main results table for the multifactorial linear regression.

```{r}
table_results_hits <- table_results |> 
  filter(
    Significant
  ) |> 
  mutate(
    Regulation = if_else(logFC > 0, "Up", "Down")
  )

list_upset <- map2(
  .x = rep(c("sgNTC", "sgIGFBP3"), each = 2),
  .y = rep(c("Up", "Down"), times = 2),
  .f = \(x, y) unique(pull(filter(table_results_hits, Condition == x, Regulation == y), Genes))
)

names(list_upset) <- paste(rep(c("sgNTC", "sgIGFBP3"), each = 2), rep(c("Up", "Down"), times = 2))
  
```

## Visualization

Firstly, the upset plot

```{r}
pdf(
  file = here("plots/20240514_plot_upset_dia.pdf"),
  width = 5,
  height = 5
)

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

dev.off()
```


The code below creates the lollipop plot that shows the effect size for significant genes. Some of them are hand-picked for labeling in the context of the manuscript.

```{r}
labels_hits_IGFBP3 <- c("DDX60", "ISG15", "MYD88", "OAS3", "PLA2G4A", "STAT1", "IFRD1", "IFIT5", "COL2A1")

table_hits_IGFBP3 <- table_results |> 
  select(Condition, Genes, logFC) |> 
  inner_join(
    table_results_factorial,
    by = "Genes"
  ) |> 
  filter(
    Significant
  ) |> 
  mutate(
    mean_logFC = logFC.x[1],
    .by = Genes
  ) |> 
  mutate(
    `Knock out` = if_else(logFC.y > 0, "logFC KO > WT", "logFC KO < WT"),
    rank = dense_rank(mean_logFC),
    label = if_else(Genes %in% labels_hits_IGFBP3, Genes, NA)
  ) |> 
  select(Condition, Genes, logFC.x, P.Value, adj.P.Val, `Knock out`, rank, label)

plot_hits_IGFBP3 <- table_hits_IGFBP3 |>
  ggplot(aes(
    x = rank,
    y = logFC.x,
    group = Genes,
    color = `Knock out`,
    label = label
  )) +
  geom_hline(
    yintercept = 0,
    linetype = "dotted"
  ) +
  geom_label_repel(
    data = filter(table_hits_IGFBP3, Condition == "sgIGFBP3", label == "COL2A1"),
    color = "black",
    size = 2,
    nudge_y = -1,
    label.padding = 0.1,
    max.overlaps = Inf,
    min.segment.length = 0,
    point.padding = 0,
    box.padding = 0.1,
    segment.size = 0.2
  ) +
  geom_label_repel(
    data = filter(table_hits_IGFBP3, Condition == "sgIGFBP3", label != "COL2A1", !is.na(label)),
    color = "black",
    size = 2,
    nudge_y = 3,
    label.padding = 0.1,
    max.overlaps = Inf,
    min.segment.length = 0,
    point.padding = 0,
    box.padding = 0.1,
    segment.size = 0.2
  ) +
  geom_point(
    data = filter(table_hits_IGFBP3, Condition == "sgIGFBP3"),
    size = 0.5,
  ) +
  geom_line(
    linewidth = 0.3
  ) +
  geom_point(
    data = filter(table_hits_IGFBP3, Condition == "sgNTC"),
    size = 0.5,
    color = "black"
  ) +
  scale_color_manual(
    values = c(
      `logFC KO > WT` = "blue",
      `logFC KO < WT` = "red"
    )
  ) +
  scale_x_continuous(
    name = "Genes"
  ) +
  scale_y_continuous(
    name = "Fold change D4 vs. NTC (log2)",
    limits = c(-4.5, 4.5),
    breaks = seq(-4, 4, by = 2)
  ) +
  theme_few() +
  theme(
    axis.ticks.x = element_blank(),
    legend.position = "bottom",
    legend.title = element_blank(),
    axis.text.x.bottom = element_blank()
  )


plot_hits_IGFBP3

```

```{r}
ggsave(
  filename = here("plots/20240522_plot_hits_IGFBP3.pdf"),
  plot = plot_hits_IGFBP3,
  width = 10,
  height = 9,
  units = "cm"
)
```

## Save results

```{r}
table_hits_combined_wide <- table_results_combined |> 
  filter(
    Significant_Factorial
  ) |> 
  select(-c(protein_group, B, AveExpr, Significant_Factorial)) |> 
  pivot_wider(
    names_from = Condition,
    values_from = c(logFC_Single, CI.L_Single, CI.R_Single, t_Single, P.Value_Single, adj.P.Val_Single, Significant_Single)
  ) |> 
  select(Genes, UniProtKB_AC, Direction, contains("Factorial"), contains("sgNTC"), contains("sgIGFBP3"))
```


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

write_tsv(
  x = table_results,
  file = here("data/table_results_proteomics_DIA_single.tsv")
)

write_tsv(
  x = table_results_factorial,
  file = here("data/table_results_proteomics_DIA_factorial.tsv")
)

writexl::write_xlsx(
  x = table_hits_combined_wide,
  path = here("tables/table_hits_combined_wide.xlsx"),
  format_headers = FALSE
)
```
