okl_qc_processing.Rmd

---
title: "Descriptive analysis of OKL affinity data"
output: html_document
editor_options:
  chunk_output_type: console
---

```{r setup}
library(tidyverse)
library(synExtra)
library(here)
library(fst)
library(data.table)
library(powerjoin)
# library(dtplyr)

theme_set(theme_minimal())

synapser::synLogin()
syn <- synDownloader("~/data/", .cache = TRUE)
```

```{r loading}
inputs <- list(
  okl_single_dose = "syn26486828",
  okl_single_dose_datasets = "syn52504516",
  okl_single_dose_no_agg = "syn52505898",
  kinase_info = "syn51286743"
)

input_files <- map(
  inputs, syn
)

single_dose <- input_files[["okl_single_dose"]] %>%
  fread()

single_dose_datasets <- input_files[["okl_single_dose_datasets"]] %>%
  fread()

single_dose_no_agg <- input_files[["okl_single_dose_no_agg"]] %>%
  fread()

kinase_info <- input_files[["kinase_info"]] %>%
  read_csv()
```

Do duplicate measurements agree on average?

Can't take into account dupes that only exist because we repeated certain
problematic measurements.

```{r}
single_dose_dupes <- single_dose_no_agg %>%
  anti_join(
    single_dose_no_agg %>%
      filter(experiment_group == "repeat") %>%
      mutate(experiment_group = "original"),
    by = c("lspci_id", "DiscoveRx Gene Symbol", "Compound Concentration (nM)", "experiment_group")
  ) %>%
  group_by(lspci_id, `DiscoveRx Gene Symbol`, `Compound Concentration (nM)`) %>%
  filter(n() > 1) %>%
  mutate(dup_id = 1:n()) %>%
  ungroup()

single_dose_dupes %>%
  count(name, `Compound Concentration (nM)`)

single_dose_dupes %>%
  count(dup_id)

p <- single_dose_dupes %>%
  filter(!qpcr_miss) %>%
  select(
    lspci_id, name, `DiscoveRx Gene Symbol`,
    `Compound Concentration (nM)`, dup_id, `Percent Control`, experiment
  ) %>%
  pivot_wider(names_from = dup_id, values_from = c(experiment, `Percent Control`)) %>%
  ggplot(
    aes(`Percent Control_1`, `Percent Control_2`, color = fct_inseq(as.character(`Compound Concentration (nM)`)))
  ) +
    geom_point(alpha = 0.5, shape = 16) +
    labs(x = "Percent Control Experiment 1", y = "Percent Control Experiment 2", title = "Agreement between duplicate experiments")

ggsave(
  "plots/agreement_duplicate_experiments.pdf", p,
  width = 4, height = 4
)

single_dose_dupes %>%
  filter(!qpcr_miss) %>%
  select(
    lspci_id, name, `DiscoveRx Gene Symbol`,
    `Compound Concentration (nM)`, dup_id, `Percent Control`, experiment
  ) %>%
  pivot_wider(names_from = dup_id, values_from = c(experiment, `Percent Control`)) %>%
  filter(`Percent Control_1` < 10, `Percent Control_2` > 30) %>% View()
```

```{r}
p <- single_dose %>%
  filter(`Compound Concentration (nM)` %in% c(12.5, 100, 1000, 10000)) %>%
  ggplot(aes(`Percent Control`, color = as.ordered(`Compound Concentration (nM)`))) +
    geom_freqpoly() +
    scale_color_viridis_d(direction = -1) +
    labs(color = "Concentration")
ggsave("all_combs_percent_control_hist.pdf", p, width = 5, height = 3)
```

```{r}
per_conc <- single_dose %>%
  select(hmsl_id, `DiscoveRx Gene Symbol`, `Compound Concentration (nM)`, `Percent Control`) %>%
  pivot_wider(
    id_cols = c(hmsl_id, `DiscoveRx Gene Symbol`),
    names_from = `Compound Concentration (nM)`,
    values_from = `Percent Control`,
    values_fn = mean
  )

p <- per_conc %>%
  drop_na(`100`, `1000`) %>%
  ggplot(aes(`100`, `1000`)) +
  geom_point()
```

```{r}
library(rasterly)
p <- per_conc %>%
  drop_na(`100`, `1000`) %>%
  ggRasterly(
    mapping = aes(x = `100`, y = `1000`), color = rev(viridis_map),
    shape = 16
  ) +
  theme_minimal() +
  geom_segment(aes(x = 0, y = 0, xend = 100, yend = 100), size = 0.2) +
  labs(
    color = "Density"
  )
ggsave(
  here("percent_control_100_vs_1000.pdf"),
  p, width = 4, height = 3
)
```

```{r}
library(rasterly)
p <- per_conc %>%
  drop_na(`100`, `10000`) %>%
  mutate(
    across(where(is.numeric), ~log10(.x + 0.01))
  ) %>%
  ggRasterly(
    mapping = aes(x = `100`, y = `10000`), color = rev(viridis_map),
    shape = 16
  ) +
  theme_minimal() +
  labs(
    color = "Density"
  )
ggsave(
  here("percent_control_100_vs_10000_log.pdf"),
  p, width = 4, height = 3
)
```

```{r}
p <- per_conc %>%
  drop_na(`100`, `10000`) %>%
  plotly::plot_ly(x = ~`100`, y = ~`10000`) %>%®
  add_rasterly_heatmap(scaling = "log")
```


Number of targets per drug under threshold

```{r}
binding_stats <- single_dose %>%
  group_by(`Compound Concentration (nM)`, `hmsl_id`) %>%
  summarize(
    n_total = n(),
    n_decided = sum(tas %in% c(2, 10), na.rm = TRUE),
    n_binding = sum(tas == 2, na.rm = TRUE),
    n_nonbinding = sum(tas == 10, na.rm = TRUE),
    n_below_50 = sum(`Percent Control` < 50),
    .groups = "drop"
  ) %>%
  mutate(
    frac_binding = n_binding / n_decided,
    frac_nonbinding = n_nonbinding / n_decided,
    frac_below_50 = n_below_50 / n_total
  )
```

```{r}
p <- binding_stats %>%
  drop_na(n_binding) %>%
  filter(`Compound Concentration (nM)` %in% c(12.5, 100, 1000, 10000)) %>%
  ggplot(aes(n_binding, color = as.ordered(`Compound Concentration (nM)`))) +
    geom_freqpoly() +
    labs(color = "Concentration")

p <- binding_stats %>%
  filter(`Compound Concentration (nM)` %in% c(12.5, 100, 1000, 10000)) %>%
  ggplot(aes(n_below_50, color = as.ordered(`Compound Concentration (nM)`))) +
    geom_freqpoly() +
    labs(color = "Concentration", x = "N inhibition below 50%")
ggsave("n_below_50_histogram.pdf", p, width = 5, height = 3)
```




```{r}
p <- binding_stats %>%
  filter(`Compound Concentration (nM)` %in% c(12.5, 100, 1000, 10000)) %>%
  arrange(desc(n_binding)) %>%
  group_by(`Compound Concentration (nM)`) %>%
  mutate(
    rank_n_binding = seq_len(n())
  ) %>%
  ungroup() %>%
  ggplot(
    aes(
      rank_n_binding, n_binding,
      color = as.ordered(`Compound Concentration (nM)`)
    )
  ) +
    geom_step() +
    labs(
      title = "N of bound targets (TAS = 2) per drug",
      color = "Concentration"
    )

p <- binding_stats %>%
  filter(`Compound Concentration (nM)` %in% c(12.5, 100, 1000, 10000)) %>%
  arrange(desc(n_below_50)) %>%
  group_by(`Compound Concentration (nM)`) %>%
  mutate(
    rank_n_below_50 = seq_len(n())
  ) %>%
  ungroup() %>%
  ggplot(aes(rank_n_below_50, n_below_50, color = as.ordered(`Compound Concentration (nM)`))) +
    geom_step() +
    labs(
      title = "N of targets with >50% inhibition per drug",
      color = "Concentration"
    )

p <- binding_stats %>%
  filter(`Compound Concentration (nM)` %in% c(12.5, 100, 1000, 10000)) %>%
  arrange(desc(frac_binding)) %>%
  group_by(`Compound Concentration (nM)`) %>%
  mutate(
    rank_frac_binding = seq_len(n())
  ) %>%
  ungroup() %>%
  ggplot(
    aes(
      rank_frac_binding, frac_binding,
      color = as.ordered(`Compound Concentration (nM)`)
    )
  ) +
    geom_step(alpha = 0.8) +
    labs(
      title = "Fraction of targets that are bound by drug",
      color = "Concentration"
    )
ggsave("rank_frac_binding.pdf", p, width = 5, height = 3)

p <- binding_stats %>%
  filter(`Compound Concentration (nM)` %in% c(12.5, 100, 1000, 10000)) %>%
  arrange(desc(frac_below_50)) %>%
  group_by(`Compound Concentration (nM)`) %>%
  mutate(
    rank_frac_below_50 = seq_len(n())
  ) %>%
  ungroup() %>%
  ggplot(aes(rank_frac_below_50, frac_below_50, color = as.ordered(`Compound Concentration (nM)`))) +
    geom_step(alpha = 0.8) +
    labs(
      title = "Fraction of targets with >50% inhibition by drug",
      color = "Concentration"
    )
ggsave("rank_frac_below_50.pdf", p, width = 5, height = 3)
```

### Heatmaps

```{r}
percent_remaining_mats <- single_dose %>%
  inner_join(
    kinase_info %>%
      select(`DiscoveRx Gene Symbol`, `Kinase Form`)
  ) %>%
  # filter(`Kinase Form` == "Wild Type") %>%
  group_nest(`Compound Concentration (nM)`) %>%
  crossing(trans = c("log", "lin")) %>%
  mutate(
    mat = map2(
      data, trans,
      function(data, trans) {
        trans_fun <- if (trans == "log")
          function(x) log10(x + 0.01)
        else
          identity
        data %>%
          mutate(across(`Percent Control`, trans_fun)) %>%
          select(hmsl_id, `DiscoveRx Gene Symbol`, `Percent Control`) %>%
          pivot_wider(values_from = `Percent Control`, names_from = hmsl_id, values_fn = mean) %>%
          column_to_rownames("DiscoveRx Gene Symbol")
      }
    )
  ) %>%
  filter(`Compound Concentration (nM)` %in% c(12.5, 100, 1000, 10000))

library(pheatmap)
library(seriation)
viridis_map <- viridis::viridis(100)
percent_remaining_heatmaps <- percent_remaining_mats %>%
  mutate(
    heatmap = map(
      mat,
      function(mat) {
        set.seed(42)
        mat <- t(mat)
        dr <- dist(mat, method = "euclidean")
        dc <- dist(t(mat), method = "euclidean")
        cr <- hclust(dr, method = "average") %>%
          reorder(dr, method = "olo")
        cc <- hclust(dc, method = "average") %>%
          reorder(dc, method = "olo")
        pheatmap(
          mat,
          show_rownames = FALSE, show_colnames = FALSE,
          color = rev(viridis_map),
          cluster_rows = cr, cluster_cols = cc,
          treeheight_row = 20, treeheight_col = 20
        )
      }
    )
  )

dir.create("heatmaps", showWarnings = FALSE)
pwalk(
  percent_remaining_heatmaps,
  function(heatmap, `Compound Concentration (nM)`, trans, ...) {
    ggsave(
      paste0("heatmaps/heatmap_all_", `Compound Concentration (nM)`, "_", trans, ".pdf"),
      heatmap, width = 5, height = 3
    )
  }
)

```


```{r}
p <- single_dose %>%
  filter(hmsl_id == "HMSL10139", `DiscoveRx Gene Symbol` %in% c("TYK2(JH1domain-catalytic)", "JAK3(JH1domain-catalytic)", "NEK10", "PRKCE")) %>%
  ggplot(aes(`Compound Concentration (nM)`, `Percent Control`, color = `Entrez Gene Symbol`)) +
    # geom_line() +
    # geom_point(alpha = 0.8) +
    geom_beeswarm(alpha = 0.8) +
    scale_x_log10() +
    theme_minimal() +
    geom_line() +
    labs(title = "JAK inhibitor AZD1480")

ggsave(
  "dose_respons_scatter_AZD1480.pdf",
  p, width = 4, height = 2.5
)

p <- single_dose %>%
  filter(hmsl_id == "HMSL10139", `DiscoveRx Gene Symbol` %in% c("TYK2(JH1domain-catalytic)", "JAK3(JH1domain-catalytic)", "NEK10", "PRKCE")) %>%
  ggplot(aes(`Compound Concentration (nM)`, log10(`Percent Control` + 0.0005), color = `Entrez Gene Symbol`)) +
    # geom_line() +
    # geom_point(alpha = 0.8) +
    geom_beeswarm(alpha = 0.8) +
    scale_x_log10() +
    theme_minimal() +
    geom_line() +
    labs(title = "JAK inhibitor AZD1480", y = "log10(Percent Control)")

ggsave(
  "dose_response_scatter_log_AZD1480.pdf",
  p, width = 4, height = 2.5
)
```

```{r}
single_dose_trans <- single_dose %>%
  mutate(
    log_inhibition = log10(`Percent Control` + 0.01),
    log_conc = log10(`Compound Concentration (nM)`),
    inhib_hit = `Percent Control` <= 35 & !qpcr_miss
  ) %>%
  as_tibble()
```

```{r}
filter_drug <- function(name, data) {
  data_filtered <- if (is.numeric(name))
    data %>%
      filter(lspci_id == .env$name)
  else if (str_starts(name, "HMSL"))
    data %>%
      filter(hmsl_id == .env$name)
  else
    data %>%
      filter(name == .env$name)
  data_filtered
}

plot_beeswarm_single_drug <- function(name) {
  library(ggbeeswarm)
  lspci_id <- if (is.numeric(name))
    name
  else if (str_starts(name, "HMSL"))
    compound_mapping[external_id == name][["lspci_id"]][1]
  else
    compound_dictionary[pref_name == name][["lspci_id"]][1, 1]
  single_dose %>%
    filter(lspci_id == .env$lspci_id) %>%
    ggplot(aes(as.ordered(`Compound Concentration (nM)`), `Percent Control`)) +
      geom_quasirandom()
}

library(seriation)
cluster_fun <- function(x) {
  d <- dist(x)
  #   d <- as.dist(1 - cor(t(x)))
  # browser()
  clust <- hclust(d, method = "average")
  reorder(clust, dist = d, method = "olo")
}

viridis_map <- viridis::viridis(100)

plot_heatmap_single_drug <- function(name, data) {
  library(ggbeeswarm)
  data_filtered <- filter_drug(name, data)
  d <- data_filtered %>%
    mutate(across(`Compound Concentration (nM)`, ~fct_inseq(as.character(.x)))) %>%
    arrange(`Compound Concentration (nM)`) %>%
    select(`DiscoveRx Gene Symbol`, `Compound Concentration (nM)`, `Percent Control`) %>%
    pivot_wider(names_from = `Compound Concentration (nM)`, values_from = `Percent Control`) %>%
    column_to_rownames("DiscoveRx Gene Symbol")
  ComplexHeatmap::Heatmap(
    as.matrix(d),
    show_column_names = TRUE,
    cluster_columns = FALSE,
    cluster_rows = cluster_fun,
    col = rev(viridis_map)
  )
  # pheatmap(
  #   d, show_rownames = TRUE, show_colnames = FALSE,
  #   cluster_rows = FALSE,
  #   color = rev(viridis_map)
  # )
    # ggplot(
    #   aes(
    #     `DiscoveRx Gene Symbol`,
    #     as.ordered(`Compound Concentration (nM)`),
    #     fill = `Percent Control`
    #   )
    # ) +
    # geom_tile() +
    # scale_fill_viridis_c(direction = -1)
}

plot_heatmap_single_drug(
  "ERLOTINIB",
  single_dose_datasets %>%
  filter(`Compound Concentration (nM)` %in% c(12.5, 100, 1000, 10000), dataset == "only_original")
)

plot_single_drug <- function(drug, target, data, point_aes = NULL) {
  library(ggbeeswarm)
  data_filtered <- filter_drug(drug, data) %>%
    filter(str_detect(`DiscoveRx Gene Symbol`, target))
  p <- data_filtered %>%
    ggplot(aes(as.ordered(`Compound Concentration (nM)`), `Percent Control`)) +
      geom_quasirandom(
        mapping = point_aes,
        groupOnX = TRUE,
        data = ~filter(.x, !qpcr_miss)
      ) +
      geom_point(
        data = ~filter(.x, qpcr_miss),
        y = 100, color = "red", size = 3, shape = 8
      )
  if (length(unique(data_filtered$`DiscoveRx Gene Symbol`)) > 1)
    p <- p +
      facet_wrap(~`DiscoveRx Gene Symbol`)
  p
}

p <- plot_single_drug(
  "TOFACITINIB", "TYK2",
  single_dose_no_agg %>%
  filter(
    `Compound Concentration (nM)` %in% c(12.5, 100, 1000, 10000)
  ),
  point_aes = aes(color = experiment)
)
p

p <- plot_single_drug(
  "TOFACITINIB", "TYK2",
  single_dose_datasets %>%
    filter(
      `Compound Concentration (nM)` %in% c(12.5, 100, 1000, 10000),
      dataset == "original_repeat_replaced"
    )
)
p

p <- plot_heatmap_single_drug("HMSL10139")
withr::with_pdf(
  "AZD1480_dose_response_hm.pdf",
  ComplexHeatmap::draw(plot_heatmap_single_drug("HMSL10139")), width = 6, height = 2
)
ggsave(
  "AZD1480_dose_response_hm.pdf",
  p, width = 6, height = 2.5
)

```


```{r}
single_dose_datasets_trans <- single_dose_datasets %>%
  mutate(
    log_inhibition = log10(`Percent Control` + 0.01),
    log_conc = log10(`Compound Concentration (nM)`),
    inhib_hit = `Percent Control` <= 35 & !qpcr_miss
  ) %>%
  as_tibble()

single_dose_groups <- single_dose_datasets_trans %>%
  filter(!qpcr_miss) %>%
  arrange(`Compound Concentration (nM)`) %>%
  group_nest(
    dataset, lspci_id, hmsl_id, library, name, `DiscoveRx Gene Symbol`,
    `Entrez Gene Symbol`, entrezgene_id, ensembl_gene_id, hgnc_symbol,
    exclude_target
  )

single_dose_hits <- single_dose_groups %>%
  mutate(
    lin_mod = map(
      data,
      ~lm(`Percent Control` ~ log_conc, data = .x)
    ),
    inhib_concordant = map_lgl(
      data,
      \(df) {
        if (any(df$inhib_hit)) {
          inhib_idx <- which(df$inhib_hit)[1]
          if (inhib_idx == length(df$inhib_hit))
            TRUE
          else {
            after <- tail(df$`Percent Control`, -inhib_idx)
            x <- !any((after > 35) & (after > (2 * df$`Percent Control`[inhib_idx])))
            # if (x) browser()
            x
          }
        } else {
            TRUE
        }
      }
    )
  )

library(broom)
single_dose_classification_res <- single_dose_hits %>%
  mutate(
    lin_mod_res = map(lin_mod, tidy)
  )

single_dose_classification_res2 <- single_dose_classification_res %>%
  unnest(lin_mod_res) %>%
  filter(term == "log_conc") %>%
  select(-c(term, std.error, statistic)) %>%
  rename(
    slope_percent_control = estimate,
    p_value_slope = p.value
  )

single_dose_classification <- single_dose_classification_res2 %>%
  rowwise() %>%
  mutate(
    classification = {
      if (!inhib_concordant) {
        "discordant"
      } else if (all(!data$inhib_hit)) {
        "non-binding"
      } else if (inhib_concordant && sum(data$inhib_hit) >= 2) {
        "binding"
      } else if (
        !is.na(slope_percent_control) && !is.na(p_value_slope) &&
          slope_percent_control < 0 && p_value_slope < 0.1
      ) {
        "weakly-binding"
      } else {
        "undetermined"
      }
    } %>%
      factor(
        levels = c("binding", "weakly-binding", "non-binding", "undetermined", "discordant"),
        ordered = TRUE
      )
  ) %>%
  ungroup()

library(qs)
qsave(
  single_dose_classification,
  file.path("classification", "okl_eurofins_classification_raw.qs")
)
# single_dose_classification <- qread(file.path("classification", "okl_eurofins_classification_raw.qs"))

write_csv(
  single_dose_classification %>%
    select(where(negate(is.list))),
  "classification/okl_eurofins_classification.csv"
)
# single_dose_classification <- read_csv("classification/okl_eurofins_classification.csv")
```

Troubleshoot erlotinib

```{r}
single_dose_classification %>% filter(name == "ERLOTINIB") %>% count(dataset, classification)

set.seed(42)
plot_single_drug(
  "ERLOTINIB",
  single_dose_classification %>%
    filter(name == "ERLOTINIB", dataset == "original_repeat_replaced", classification == "discordant") %>%
    pull(`DiscoveRx Gene Symbol`) %>%
    sample(20) %>%
    str_escape() %>%
    paste(collapse = "|"),
  single_dose_no_agg %>%
    filter(
      `Compound Concentration (nM)` %in% c(12.5, 100, 1000, 10000)
    ),
  aes(color = experiment)
) +
  theme_bw()

  single_dose_classification %>%
    filter(name == "ERLOTINIB", dataset == "original_repeat_replaced", classification == "discordant") %>%
    select(where(negate(is.list))) %>%
    View()

single_dose_classification %>%
    filter(name == "ERLOTINIB", dataset == "original_repeat_replaced", classification == "discordant") %>%
    chuck("data", 1)

set.seed(42)
plot_single_drug(
  "ERLOTINIB",
  single_dose_classification %>%
    filter(name == "ERLOTINIB", dataset == "original_repeat_replaced", classification == "discordant") %>%
    pull(`DiscoveRx Gene Symbol`) %>%
    sample(20) %>%
    str_escape() %>%
    paste(collapse = "|"),
  single_dose_datasets %>%
    filter(
      dataset == "only_original",
      `Compound Concentration (nM)` %in% c(12.5, 100, 1000, 10000)
    ),
  aes(color = experiment)
) +
  theme_bw()
```

   dataset                  classification     n
   <chr>                    <ord>          <int>
 1 only_original            binding           37
 2 only_original            weakly-binding    19
 3 only_original            non-binding      319
 4 only_original            undetermined      43
 5 only_original            discordant        50
 6 original_repeat_replaced binding           10
 7 original_repeat_replaced weakly-binding     1
 8 original_repeat_replaced non-binding      355

It was just an issue with the ordering of the data. The data has to be ordered
by concentration in the df fed into the function determining concordance.

```{r}
single_dose_classification %>%
  filter(`DiscoveRx Gene Symbol` == "LATS1") %>%
  mutate(
    dataset = if_else(superseded, "original", "repeat")
  ) %>%
  group_by(classification, dataset) %>%
  summarize(
    n = n(),
    frac_total = n() / nrow(.),
    frac_binding = n() / sum(map_lgl(.$data, ~any(.x$inhib_hit))),
    .groups = "drop"
  ) %>%
  group_by(dataset) %>%
  summarize(n_total = sum(n))

single_dose_classification %>%
  filter(`DiscoveRx Gene Symbol` == "LATS1") %>%
  mutate(
    dataset = if_else(superseded, "original", "repeat")
  ) %>%
  group_by(classification, dataset) %>%
  ggplot(aes(classification, fill = dataset)) +
    geom_bar(position = "dodge")
```


```{r}
single_dose_classification_by_compound <- single_dose_classification %>%
  count(dataset, lspci_id, name, classification) %>%
  complete(nesting(dataset, lspci_id, name), classification, fill = list(n = 0)) %>%
  group_by(dataset, lspci_id, name) %>%
  mutate(frac = n / sum(n)) %>%
  ungroup()

p <- single_dose_classification_by_compound %>%
  group_by(
    dataset
  ) %>%
  summarize(
    p = list(
      cur_data() %>%
        mutate(
          name = factor(
            name,
            levels = filter(., classification == "binding") %>%
              arrange(frac) %>%
              pull(name)
          )
        ) %>%
        arrange(desc(name)) %>%
        mutate(split_var = floor(seq(from = 0, to = 2.99999, length.out = n()))) %>%
        ggplot(aes(frac, name, color = classification)) +
          geom_point() +
          scale_color_manual(
            values = c(
              "binding" = "firebrick1",
              "weakly-binding" = "darkorange",
              "non-binding" = "black",
              "undetermined" = "gray",
              "discordant" = "blue"
            )
      ) +
      facet_wrap(~split_var, scales = "free_y", nrow = 1) +
      theme(strip.text = element_blank()) +
      labs(x = "Fraction of dose-response curves", y = "Compound")),
    .groups = "drop"
  )

dir.create("classification")
pwalk(
  p,
  \(dataset, p) {
    ggsave(
      file.path("classification", paste0("lm_classification_waterfall_", dataset, ".pdf")),
      p, width = 18, height = 10
    )
  }
)
single_dose_classification_by_kinase <- single_dose_classification %>%
  count(dataset, `DiscoveRx Gene Symbol`, hgnc_symbol, classification) %>%
  complete(nesting(dataset, `DiscoveRx Gene Symbol`, hgnc_symbol), classification, fill = list(n = 0)) %>%
  group_by(dataset, `DiscoveRx Gene Symbol`, hgnc_symbol) %>%
  mutate(frac = n / sum(n)) %>%
  ungroup()

p <- single_dose_classification_by_kinase %>%
  group_by(
    dataset
  ) %>%
  summarize(
    p = list(
      cur_data() %>%
        mutate(
          `DiscoveRx Gene Symbol` = factor(
            `DiscoveRx Gene Symbol`,
            levels = filter(., classification == "binding") %>%
              arrange(frac) %>%
              pull(`DiscoveRx Gene Symbol`)
          )
        ) %>%
        arrange(desc(`DiscoveRx Gene Symbol`)) %>%
        mutate(split_var = floor(seq(from = 0, to = 2.99999, length.out = n()))) %>%
        ggplot(aes(frac, `DiscoveRx Gene Symbol`, color = classification)) +
          geom_point() +
          scale_color_manual(
            values = c(
              "binding" = "firebrick1",
              "weakly-binding" = "darkorange",
              "non-binding" = "black",
              "undetermined" = "gray",
              "discordant" = "blue"
            )
      ) +
      facet_wrap(~split_var, scales = "free_y", nrow = 1) +
      theme(strip.text = element_blank()) +
      labs(x = "Fraction of dose-response curves", y = "Kinase")),
    .groups = "drop"
  )

dir.create("classification")
pwalk(
  p,
  \(dataset, p) {
    ggsave(
      file.path("classification", paste0("lm_classification_waterfall_kinases_", dataset, ".pdf")),
      p, width = 18, height = 14
    )
  }
)

library(ggbeeswarm)
p <- single_dose_classification %>%
  count(dataset, lspci_id, name, classification) %>%
  complete(nesting(dataset, lspci_id, name), classification, fill = list(n = 0)) %>%
  group_by(dataset, lspci_id, name) %>%
  mutate(frac = n / sum(n)) %>%
  ungroup() %>%
  group_by(dataset) %>%
  summarize(
    p = list(
      cur_data() %>%
      mutate(
        name = factor(
          name,
          levels = filter(., classification == "binding") %>%
            arrange(frac) %>%
            pull(name)
        )
      ) %>%
      ggplot(aes(classification, frac)) +
        geom_quasirandom()
    ),
    .groups = "drop"
  )

pwalk(
  p,
  \(dataset, p) {
    ggsave(
      file.path("classification", paste0("hist_classes_", dataset, ".pdf")), p,
      width = 8, height = 6
    )
  }
)

p <- single_dose_classification_by_compound %>%
  pivot_wider(names_from = dataset, values_from = c(n, frac)) %>%
  ggplot(
    aes(x = frac_only_original, y = frac_original_repeat_replaced)
  ) +
  geom_point() +
  facet_wrap(~classification, scales = "free") +
  geom_abline(slope = 1, intercept = 0)

library(ggalluvial)
single_dose_classification_change <- single_dose_classification %>%
  select(
    dataset, lspci_id, hmsl_id, library, name,
    `DiscoveRx Gene Symbol`, `Entrez Gene Symbol`,
    hgnc_symbol, inhib_concordant, classification
  ) %>%
  pivot_wider(names_from = dataset, values_from = c(inhib_concordant, classification))

single_dose_classification_change_agg <- single_dose_classification_change %>%
  count(classification_only_original, classification_original_repeat_replaced)


p <- single_dose_classification_change_agg %>%
  # filter(classification_only_original == "discordant" | classification_original_repeat_replaced == "discordant") %>%
  ggplot(
    aes(axis1 = classification_only_original, axis2 = classification_original_repeat_replaced, y = n)
  ) +
  geom_alluvium(
    aes(fill = classification_only_original == "discordant")
  ) +
  scale_fill_manual(
    values = c(`TRUE` = "red", `FALSE` = "white")
  ) +
  geom_stratum() +
  geom_text(stat = "stratum", aes(label = after_stat(stratum)))

library(egg)
set.seed(42)
p <- single_dose_classification %>%
  group_by(classification) %>%
  slice_sample(n = 10) %>%
  ungroup() %>%
  mutate(binding_name = paste(name, `Entrez Gene Symbol`, sep = " - ")) %>%
  unnest(data) %>%
  ggplot(aes(`Compound Concentration (nM)`, `Percent Control`)) +
    geom_rect(
      aes(fill = classification, xmin = 10, xmax = 10000,
      ymin = 0, ymax = 100),
      alpha = 0.5,
      data = ~.x %>%
        distinct(binding_name, classification),
      inherit.aes = FALSE
    ) +
    geom_point() +
    geom_line() +
    scale_x_log10() +
    geom_hline(yintercept = 35) +
    facet_wrap(~classification + binding_name) +
    # geom_custom(
    #   data = ~.x %>%
    #     distinct(binding_name, classification),
    #   mapping = aes(fill = classification),
    #   grob_fun = function(x) rectGrob(gp=gpar(fill=fill,col=NA, alpha=0.5)),
    #   inherit.aes = TRUE
    # ) +
    geom_text(
      aes(label = paste(round(estimate, digits = 1), signif(p.value, digits = 2))),
      data = ~.x %>%
        distinct(classification, binding_name, estimate, p.value),
      x = 20, y = 20
    )
ggsave(
  "classification/class_examples.pdf", p, width = 14, height = 12
)

single_dose_classification %>%
  filter(
    name == "SILMITASERTIB",
    `Entrez Gene Symbol` == "MAP2K2"
  ) %>%
  View()

library(ggrepel)
p <- single_dose_classification %>%
  count(dataset, lspci_id, name, classification) %>%
  complete(dataset, nesting(lspci_id, name), classification, fill = list(n = 0)) %>%
  group_by(lspci_id, name) %>%
  mutate(frac = n / sum(n)) %>%
  ungroup() %>%
  # mutate(
  #   name = factor(
  #     name,
  #     levels = filter(., classification == "binding") %>%
  #       arrange(frac) %>%
  #       pull(name)
  #   )
  # ) %>%
  select(dataset, lspci_id, name, classification, frac) %>%
  pivot_wider(names_from = classification, values_from = frac) %>%
  ggplot(aes(binding, discordant)) +
    geom_point() +
    geom_text_repel(
      aes(label = name),
      data = ~.x %>%
        mutate(name = if_else(discordant < 0.02, "", name)),
      min.segment.length = 0.2,
      max.overlaps = Inf,
      max.iter = 50000,
      max.time = 2,
      size = 3
    ) +
    labs(
      x = "Proportion of binding dose-response curves",
      y = "Proportion of discordant dose-response curves",
      title = "Discordant vs binding per compound"
    ) +
    scale_x_continuous(labels = scales::percent) +
    scale_y_continuous(labels = scales::percent) +
    facet_wrap(~dataset)

ggsave(
  "classification/discordant_vs_binding.pdf", p, width = 12, height = 7
)


p <- single_dose_classification %>%
  count(dataset, `DiscoveRx Gene Symbol`, classification) %>%
  complete(dataset, `DiscoveRx Gene Symbol`, classification, fill = list(n = 0)) %>%
  group_by(dataset, `DiscoveRx Gene Symbol`) %>%
  mutate(frac = n / sum(n)) %>%
  ungroup() %>%
  # mutate(
  #   name = factor(
  #     name,
  #     levels = filter(., classification == "binding") %>%
  #       arrange(frac) %>%
  #       pull(name)
  #   )
  # ) %>%
  select(dataset, `DiscoveRx Gene Symbol`, classification, frac) %>%
  pivot_wider(names_from = classification, values_from = frac) %>%
  ggplot(aes(binding, discordant)) +
    geom_point() +
    geom_text_repel(
      aes(label = `DiscoveRx Gene Symbol`),
      data = ~.x %>%
        mutate(`DiscoveRx Gene Symbol` = if_else(discordant < 0.08, "", `DiscoveRx Gene Symbol`)),
      min.segment.length = 0.2,
      max.overlaps = Inf,
      max.iter = 50000,
      max.time = 2,
      size = 3
    ) +
    labs(
      x = "Proportion of binding dose-response curves",
      y = "Proportion of discordant dose-response curves",
      title = "Discordant vs binding per kinase"
    ) +
    scale_x_continuous(labels = scales::percent) +
    scale_y_continuous(labels = scales::percent) +
    facet_wrap(~dataset)

ggsave(
  "classification/discordant_vs_binding_kinases.pdf", p, width = 12, height = 7
)

```

```{r pseudo_kd}
compute_pseudo_kd <- function(df) {
  df <- df %>%
    mutate(below_50 = `Percent Control` < 50)
  if (length(unique(df$below_50)) == 2 && !df$below_50[1]) {
    over_under_50_idx <- which(df$below_50)[1] + c(-1, 0)
    mdl_data <- df[over_under_50_idx, ]
    if (nrow(mdl_data) == 1)
      browser()
    mdl <- lm(`Percent Control` ~ log_conc, data = mdl_data)
    list(
      pseudo_kd = 10^(chemCal::inverse.predict(mdl, 50)$Prediction),
      pseudo_kd_relation = "=",
      # Make sure all values beyond the kd value are below 50
      # otherwise label as discordant
      pseudo_kd_concordant = if (
        over_under_50_idx[2] == nrow(df) ||
        all(tail(df$`Percent Control`, -over_under_50_idx[2]) < 50)
      ) {
        "concordant"
      } else {
        "discordant"
      }
    )
  } else if (all(df$below_50)) {
    list(
      pseudo_kd = min(df$`Compound Concentration (nM)`),
      pseudo_kd_relation= "<",
      pseudo_kd_concordant = "concordant"
    )
  } else if (!any(df$below_50)) {
    list(
      pseudo_kd = max(df$`Compound Concentration (nM)`),
      pseudo_kd_relation= ">",
      pseudo_kd_concordant = "concordant"
    )
  } else {
    list(
      pseudo_kd = NA_real_,
      pseudo_kd_relation = NA_character_,
      pseudo_kd_concordant = NA_character_
    )
  }
}


library(chemCal)
pseudo_kd_raw <- single_dose_classification %>%
  mutate(
    pseudo_kd = map(
      data,
      compute_pseudo_kd
    )
  )

pseudo_kd <- pseudo_kd_raw %>%
  unnest_wider(
    pseudo_kd
  ) %>%
  mutate(
    pseudo_kd_n_doses_missing = map_int(
      data,
      \(x) {
        # Check that all doses are present
        length(setdiff(c(12.5, 100, 1000, 10000), x$`Compound Concentration (nM)`))
      }
    )
  )

qsave(
  pseudo_kd,
  "classification/okl_eurofins_pseudo_kd_raw.qs"
)
# pseudo_kd <- qread("classification/okl_eurofins_pseudo_kd_raw.qs")

write_csv(
  pseudo_kd %>%
    select(where(negate(is.list))),
  "classification/okl_eurofins_pseudo_kd.csv"
)
# pseudo_kd <- fread("classification/okl_eurofins_pseudo_kd.csv")
```


```{r upload_synapse}
synStoreMany(
  c(
    "classification/okl_eurofins_pseudo_kd.csv",
    "classification/okl_eurofins_classification.csv",
    "classification/okl_eurofins_classification_raw.qs"
  ),
  parentId = "syn18508401",
  used = unname(inputs),
  executed = "https://github.com/labsyspharm/okl-analysis/blob/main/okl_descriptive.Rmd",
  forceVersion = FALSE
)

```


```{r}
of_interest <- c(
  "ERLOTINIB",
  "TIVOZANIB",
  "BARICITINIB"
)

single_dose_classification_of_interest <- single_dose_classification %>%
  filter(name %in% of_interest)

library(ComplexHeatmap)
library(seriation)
cluster_fun <- function(x) {
  d <- dist(x)
  #   d <- as.dist(1 - cor(t(x)))
  # browser()
  clust <- hclust(d, method = "average")
  reorder(clust, dist = d, method = "olo")
}

make_hm <- function(df_new, df_old) {
  df_old <- bind_rows(
    df_old,
    df_new %>%
      anti_join(
        distinct(df_old, `DiscoveRx Gene Symbol`, `Compound Concentration (nM)`)
      )
  )
  mat_old <- df_old %>%
    mutate(
      across(`Compound Concentration (nM)`, ~fct_inseq(as.character(.x)))
    ) %>%
    select(`DiscoveRx Gene Symbol`, `Compound Concentration (nM)`, `Percent Control`) %>%
    pivot_wider(names_from = `Compound Concentration (nM)`, values_from = `Percent Control`) %>%
    column_to_rownames("DiscoveRx Gene Symbol") %>%
    as.matrix()
  mat_new <- df_new %>%
    mutate(
      across(`Compound Concentration (nM)`, ~fct_inseq(as.character(.x)))
    ) %>%
    select(`DiscoveRx Gene Symbol`, `Compound Concentration (nM)`, `Percent Control`) %>%
    pivot_wider(names_from = `Compound Concentration (nM)`, values_from = `Percent Control`) %>%
    column_to_rownames("DiscoveRx Gene Symbol") %>%
    as.matrix()
  mat_new <- mat_new[rownames(mat_old), colnames(mat_old)]
  clust_new_rows <- cluster_fun(mat_new)
  # clust_new_cols <- cluster_fun(t(mat_new))
  list(
    hm_new = Heatmap(
      mat_new,
      cluster_rows = clust_new_rows,
      cluster_columns = FALSE,
      col = rev(viridisLite::plasma(100))
    ),
    hm_old = Heatmap(
      mat_old,
      cluster_rows = clust_new_rows,
      cluster_columns = FALSE,
      col = rev(viridisLite::plasma(100))
    )
  )
}

single_dose_classification_of_interest_hms <- single_dose_classification_of_interest %>%
  select(lspci_id, name, `DiscoveRx Gene Symbol`, superseded, data) %>%
  pivot_wider(names_from = superseded, values_from = data) %>%
  group_by(lspci_id, name) %>%
  summarize(
    across(c(`FALSE`, `TRUE`), ~list(bind_rows(.x))),
    .groups = "drop"
  ) %>%
  mutate(
    hms = map2(`FALSE`, `TRUE`, make_hm)
  )

pwalk(
  single_dose_classification_of_interest_hms,
  function(name, hms, ...) {
    withr::with_pdf(
      paste0("classification/", name, "_hm_new.pdf"),
      draw(hms$hm_new), width = 6, height = 10
    )
    withr::with_pdf(
      paste0("classification/", name, "_hm_old.pdf"),
      draw(hms$hm_old), width = 6, height = 10
    )
  }
)
```