chembl_comparison.Rmd

```{r}
library(tidyverse)
library(synExtra)
library(fst)
library(data.table)
library(powerjoin)

synapser::synLogin()

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

```{r}
kinomescan_pseudo_kds <- syn("syn51080578") %>%
  read_csv()

kinomescan_classes <- syn("syn51080579") %>%
  read_csv()

kinomescan_kinase_info <- syn("syn51286743") %>%
  read_csv()

single_dose <- syn("syn52504516") %>%
  read_csv()

compound_mapping <- syn("syn26260389") %>%
  read_fst(as.data.table = TRUE)

chembl_kd_raw <- syn("syn52854574") %>%
  read_csv()

# qpcr_misses <- syn("syn33531640") %>%
#   readxl::read_excel() %>%
#   left_join(
#     compound_mapping[, .(lspci_id, external_id)] %>%
#       unique(),
#     by = c("Compound Name" = "external_id")
#   )

# lsp_biochem <- syn("syn25173511") %>%
#   read_fst()
lsp_targets <- syn("syn25173506") %>%
  read_fst()

lsp_biochem_single <- syn("syn25173513") %>%
  read_fst()

nominal_targets <- syn("syn52947779") %>%
  read_csv() %>%
  power_inner_join(
    single_dose %>%
      distinct(lspci_id, hmsl_id),
    by = c("HMSLID" = "hmsl_id"),
    check = check_specs(
      duplicate_keys_right = "warn",
      unmatched_keys_left = "warn"
    )
  )

lsp_clinical <- syn("syn26260451") %>%
  read_fst() %>%
  as_tibble()
```

Aggregate ChEMBL IC50 values. Remove any values that
have relation that is not "=" unless the relation is ">" or ">=" and the value
is greater than 9000. This is to avoid underestimating IC50s for cases where
people just didn't assay higher concentrations.

```{r}
lsp_biochem <- lsp_biochem_single %>%
  filter(
    value_relation == "=" |
      (value_relation %in% c(">=", ">") & value > 9000)
  ) %>%
  group_by(lspci_id, lspci_target_id) %>%
  summarize(
    value = quantile(value, 0.25, na.rm = TRUE, names = FALSE),
    .groups = "drop"
  )
```

```{r pseudo_ic50_vs_kd}
class_colors <- c(
    "binding" = "firebrick1",
    "weakly-binding" = "darkorange",
    "undetermined" = "gray",
    "non-binding" = "black",
    "discordant" = "blue"
)

pseudo_kd_comp <- kinomescan_pseudo_kds %>%
  mutate(
    across(classification, \(x) factor(x, levels = names(class_colors)))
  ) %>%
  filter(dataset == "original_repeat_replaced") %>%
  power_left_join(
    lsp_targets %>%
      filter(organism == "Homo sapiens") %>%
      select(lspci_target_id, symbol),
    by = c("hgnc_symbol" = "symbol"),
    check = check_specs(
      duplicate_keys_left = "ignore",
      duplicate_keys_right = "ignore",
      unmatched_keys_left = "warn",
      unmatched_keys_right = "ignore"
    )
  ) %>%
  filter(
    !exclude_target
  ) %>%
  inner_join(
    lsp_biochem,
    by = c("lspci_id", "lspci_target_id")
  ) %>%
  left_join(
    lsp_clinical %>%
      filter(max_phase == 4) %>%
      transmute(lspci_id, fda_approved = TRUE) %>%
      distinct(),
    by = "lspci_id"
  ) %>%
  left_join(
    nominal_targets %>%
      mutate(nominal_target = TRUE),
    by = c("lspci_id", "DiscoveRx Gene Symbol")
  ) %>%
  replace_na(
    list(fda_approved = FALSE, nominal_target = FALSE)
  ) %>%
  mutate(
    fda_target = case_when(
      !fda_approved ~ "no",
      !nominal_target ~ "offtarget",
      TRUE ~ "ontarget"
    ) %>%
      factor(levels = c("ontarget", "offtarget", "no"))
  )


```

Check weird cases where we have pseudo Kd < 1000 uM but binding
is classified as non-binding. Probably results mostly of cases
where entire binding curve is between 50-35% control so the 35% threshold
is never triggered but the 50% threshold for Kd is.

```{r}
pseudo_kd_comp %>%
    filter(pseudo_kd < 1000, classification == "non-binding", pseudo_kd_concordant == "concordant")

```

Only 152 such cases. Not worth worrying about.

Make scatterplot of ChEMBL IC50 vs KINOMEscan pseudo-Kd
Also show distribution of pseudo-IC50 values for each binding class as
beeswarm plot.

```{r}
p_scatter <- pseudo_kd_comp %>%
  drop_na(pseudo_kd) %>%
  arrange(nominal_target) %>%
  filter(
    # pseudo_ic50_relation == "=",
    pseudo_kd_concordant == "concordant",
    classification != "discordant"
  ) %>%
  # Clamp values >10000 to 10000
  mutate(value = pmin(value, 10000)) %>%
  ggplot(
    # aes(value, pseudo_ic50, color = classification, shape = pseudo_ic50_relation)
    aes(
      value, pseudo_kd,
      fill = classification, color = fda_target,
      text = paste(
        name, `DiscoveRx Gene Symbol`,
        "<br>ChEMBL", signif(value, 2), "nM Eurofins", signif(pseudo_kd, 2), "nM"
      )
    )
  ) +
  geom_point(aes(size = fda_target, shape = pseudo_kd_relation)) +
  geom_text(
    aes(label = label),
    data = \(df) df %>%
      mutate(across(c(value, pseudo_kd), \(x) -log10(x))) %>%
      summarize(
        res = cor.test(
          value, pseudo_kd
        ) %>%
          broom::tidy() %>%
          list()
      ) %>%
      unnest(res) %>%
      mutate(label = paste("italic(R)^2 ==", signif(estimate**2, 2))),
    x = -Inf, y = Inf, hjust = 0, vjust = 1,
    inherit.aes = FALSE, parse = TRUE
  ) +
  scale_shape_manual(
    values = c(
      `<` = 25, `>` = 24, `=` = 21
    ),
    guide = "none"
  ) +
  scale_fill_manual(
    values = map_chr(class_colors, \(x) alpha(x, alpha = 0.5)),
    guide = guide_legend(override.aes = list(color = alpha("white", 0), shape = 21, alpha = 1))
  ) +
  scale_color_manual(
    values = c(
      `ontarget` = alpha("black", .9),
      `offtarget` = alpha("blue", 0),
      `no` = alpha("purple", 0)
    ),
    guide = "none"
  ) +
  scale_size_manual(
    values = c(`ontarget` = 1.5, `offtarget` = 1, `no` = 1),
    guide = "none"
  ) +
  scale_x_log10(labels = scales::label_number(big.mark = ""), limits = c(1, NA)) +
  scale_y_log10(labels = scales::label_number(big.mark = "")) +
  coord_equal() +
  labs(
    x = "IC50 (ChEMBL)",
    y = "Pseudo-Kd (KINOMEscan)"
  )

ggsave(
    "plots/chembl_ic50_vs_kinomescan_pseudo_kd_new.pdf", p_scatter,
    width = 8, height = 3
)

library(plotly)

p_scatter_plotly <- ggplotly(p_scatter, tooltip = "text") %>%
  toWebGL()

htmlwidgets::saveWidget(p_scatter_plotly, "plots/chembl_ic50_vs_kinomescan_pseudo_kd.html")

library(ggbeeswarm)
library(ggtext)
library(glue)

class_labels <- glue_data(
  distinct(pseudo_kd_comp, classification) %>%
    filter(classification != "discordant"),
  "<span style='font-weight: bold; color: {class_colors[classification]}'>{classification}</span>"
) %>%
  set_names(
    distinct(pseudo_kd_comp, classification) %>%
      filter(classification != "discordant") %>%
      pull(classification)
  )

p_beeswarm <- pseudo_kd_comp %>%
  filter(
    pseudo_kd_concordant == "concordant",
    # pseudo_ic50_relation == "=",
    classification != "discordant"
  ) %>%
  mutate(
    classification = factor(
      classification,
      levels = names(class_colors) %>%
        rev()
    ),
    value = pmin(value, 10000)
  ) %>%
  ggplot(
    aes(value, classification, color = classification)
  ) +
  geom_quasirandom(groupOnX = FALSE, shape = 16, size = 1, alpha = 0.5) +
  geom_boxplot(fill = NA, outlier.shape = NA) +
  scale_color_manual(values = class_colors) +
  guides(color = "none") +
  scale_x_log10(labels = scales::label_number(big.mark = ""), limits = c(1, NA)) +
  scale_y_discrete(labels = class_labels) +
  theme(axis.text.y = element_markdown()) +
  labs(x = "IC50 (ChEMBL)", y = NULL)

ggsave(
    "plots/chembl_ic50_by_kinomescan_class_beeswarm.pdf", p_beeswarm,
    width = 7, height = 3
)

library(patchwork)
p_combined <- egg::ggarrange(
    p_scatter +
        labs(x = NULL) +
        guides(color = "none") +
        theme(
            axis.text.x = element_blank(),
            plot.margin = unit(c(0.5, 0.5, 0, 0.5), "cm")
        ),
    p_beeswarm +
        coord_fixed(ratio = .32) +
        theme(
            plot.margin = unit(c(0, 0.5, 0.5, 0.5), "cm"),
            legend.position = "bottom"
        ),
    nrow = 2,
    heights = c(2, 1.2),
    draw = FALSE
)

ggsave(
    "plots/chembl_ic50_by_comparison_combined_new_with_nominal_fda.pdf", p_combined,
    width = 6, height = 4.5
)
```

Now load the ChEMBL Kd only data and compare to KINOMEscan pseudo-IC50

Remove any ChEMBL Kd values that have relation that is not "=" unless
the relation is ">" or ">=" and the value is greater than 9000.


```{r}
chembl_kd <- chembl_kd_raw %>%
  mutate(across(standard_value, as.numeric)) %>%
  filter(
    standard_relation == "=" |
      (standard_relation %in% c(">=", ">") & standard_value > 9000)
  ) %>%
  power_inner_join(
    compound_mapping %>%
      filter(source == "chembl") %>%
      distinct(lspci_id, chembl_id_compound = external_id),
    by = "chembl_id_compound",
    check = check_specs(
      duplicate_keys_right = "warn",
      unmatched_keys_left = "warn"
    )
  ) %>%
  group_by(gene_symbol, lspci_id) %>%
  summarize(
    chembl_kd = quantile(standard_value, 0.25, na.rm = TRUE, names = FALSE),
    chembl_kd_relation = paste(sort(unique(standard_relation)), collapse = ""),
    .groups = "drop"
  ) %>%
  power_inner_join(
    kinomescan_kinase_info %>%
      distinct(hgnc_symbol, `DiscoveRx Gene Symbol`),
    by = c("gene_symbol" = "hgnc_symbol"),
    check = check_specs(
      duplicate_keys_right = "ignore",
      unmatched_keys_left = "warn"
    )
  )

chembl_kd %>%
  count(chembl_kd_relation)

pseudo_kd_kd_comp <- kinomescan_pseudo_kds %>%
  mutate(
    across(classification, \(x) factor(x, levels = names(class_colors)))
  ) %>%
  filter(
    dataset == "original_repeat_replaced",
    !exclude_target
  ) %>%
  inner_join(
    chembl_kd,
    by = c("lspci_id", "DiscoveRx Gene Symbol")
  ) %>%
  left_join(
    lsp_clinical %>%
      filter(max_phase == 4) %>%
      transmute(lspci_id, fda_approved = TRUE) %>%
      distinct(),
    by = "lspci_id"
  ) %>%
  left_join(
    nominal_targets %>%
      mutate(nominal_target = TRUE),
    by = c("lspci_id", "DiscoveRx Gene Symbol")
  ) %>%
  replace_na(
    list(fda_approved = FALSE, nominal_target = FALSE)
  ) %>%
  mutate(
    fda_target = case_when(
      !fda_approved ~ "no",
      !nominal_target ~ "offtarget",
      TRUE ~ "ontarget"
    ) %>%
      factor(levels = c("ontarget", "offtarget", "no"))
  )

```


```{r}
p_scatter <- pseudo_kd_kd_comp %>%
  drop_na(pseudo_kd) %>%
  filter(
    # pseudo_ic50_relation == "=",
    pseudo_kd_concordant == "concordant",
    classification != "discordant"
  ) %>%
  mutate(chembl_kd = pmin(chembl_kd, 10000)) %>%
  arrange(desc(fda_target)) %>%
  ggplot(
    # aes(value, pseudo_ic50, color = classification, shape = pseudo_ic50_relation)
    aes(chembl_kd, pseudo_kd, fill = classification, color = fda_target)
  ) +
  geom_point(aes(size = fda_target, shape = pseudo_kd_relation)) +
  geom_text(
    aes(label = label),
    data = \(df) df %>%
      mutate(across(c(chembl_kd, pseudo_kd), \(x) -log10(x))) %>%
      summarize(
        res = cor.test(
          chembl_kd, pseudo_kd
        ) %>%
          broom::tidy() %>%
          list()
      ) %>%
      unnest(res) %>%
      mutate(label = paste("italic(R)^2 ==", signif(estimate**2, 2))),
    x = -Inf, y = Inf, hjust = 0, vjust = 1,
    inherit.aes = FALSE, parse = TRUE
  ) +
  scale_shape_manual(
    values = c(
      `<` = 25, `>` = 24, `=` = 21
    ),
    guide = "none"
  ) +
  scale_fill_manual(
    values = map_chr(class_colors, \(x) alpha(x, alpha = 0.5)),
    guide = guide_legend(override.aes = list(color = alpha("white", 0), shape = 21, alpha = 1))
  ) +
  scale_color_manual(
    values = c(
      `ontarget` = alpha("black", .9),
      `offtarget` = alpha("blue", 0),
      `no` = alpha("purple", 0)
    ),
    guide = "none"
  ) +
  scale_size_manual(
    values = c(`ontarget` = 1.5, `offtarget` = 1, `no` = 1),
    guide = "none"
  ) +
  # Limiting to between 1 and 10000 because barely any data outside this range
  scale_x_log10(labels = scales::label_number(big.mark = ""), limits = c(1, NA)) +
  scale_y_log10(labels = scales::label_number(big.mark = "")) +
  coord_equal() +
  labs(
    x = "Kd (ChEMBL)",
    y = "Pseudo-Kd (KINOMEscan)",
    fill = "Class"
  )
# NA warning comes from values outside of plotted range

ggsave(
    "plots/chembl_kd_vs_kinomescan_pseudo_kd_new2_with_nominal.pdf", p_scatter,
    width = 8, height = 4
)

p_beeswarm <- pseudo_kd_kd_comp %>%
  filter(
    pseudo_kd_concordant == "concordant",
    # pseudo_ic50_relation == "=",
    classification != "discordant"
  ) %>%
  mutate(
    classification = factor(
      classification,
      levels = names(class_colors) %>%
        rev()
    ),
    chembl_kd = pmin(chembl_kd, 10000)
  ) %>%
  ggplot(
    aes(chembl_kd, classification, color = classification)
  ) +
  geom_quasirandom(orientation = "y", shape = 16, size = 1, alpha = 0.5, width = 0.3) +
  geom_boxplot(fill = NA, outlier.shape = NA) +
  scale_color_manual(values = class_colors) +
  guides(color = "none") +
  scale_x_log10(labels = scales::label_number(big.mark = ""), limits = c(1, NA)) +
  scale_y_discrete(labels = class_labels) +
  theme(axis.text.y = element_markdown()) +
  labs(x = "Kd (ChEMBL)", y = NULL)

ggsave(
    "plots/chembl_kd_by_kinomescan_class_beeswarm.pdf", p_beeswarm,
    width = 7, height = 3
)

p_combined <- egg::ggarrange(
    p_scatter +
        labs(x = NULL) +
        guides(color = "none") +
        theme(
            axis.text.x = element_blank(),
            plot.margin = unit(c(0.5, 0.5, 0, 0.5), "cm")
        ),
    p_beeswarm +
        coord_fixed(ratio = .32) +
        theme(
            plot.margin = unit(c(0, 0.5, 0.5, 0.5), "cm"),
            legend.position = "bottom"
        ),
    nrow = 2,
    heights = c(2, 1.2),
    draw = FALSE
)

ggsave(
    "plots/chembl_kd_by_comparison_combined_new2_with_nominal_fda.pdf", p_combined,
    width = 6, height = 4.5
)
```

Compute correlation coefficients between the ChEMBL IC50s / Kds and the KINOMEscan pseudo-IC50s.

Using Fisher r-to-z transformation to test for significant difference between the correlations.

```{r}
two_cor_test <- function(r1, r2, n1, n2) {
  (.5*log((1+r1)/(1-r1)) - .5*log((1+r2)/(1-r2))) /
    sqrt(1/(n1-3) + 1/(n2-3))
}

pseudo_kd_comp_both <- bind_rows(
  ic50 = pseudo_kd_comp %>%
    filter(
      pseudo_kd_concordant == "concordant",
      # pseudo_ic50_relation == "=",
      classification != "discordant"
    ) %>%
    select(lspci_id, name, `DiscoveRx Gene Symbol`, gene_symbol = hgnc_symbol, value, pseudo_kd, classification, pseudo_kd_relation),
  kd = pseudo_kd_kd_comp %>%
    filter(
      pseudo_kd_concordant == "concordant",
      # pseudo_ic50_relation == "=",
      classification != "discordant"
    ) %>%
    select(lspci_id, name, `DiscoveRx Gene Symbol`, gene_symbol, value = chembl_kd, pseudo_kd, classification, pseudo_kd_relation),
  .id = "chembl_type"
)

pseudo_kd_comp_cor_res <- pseudo_kd_comp_both %>%
  group_by(chembl_type) %>%
  summarize(
    res = cor.test(
      ~value + pseudo_kd, data = cur_data()
    ) %>%
      broom::tidy() %>%
      list(),
    n = n(),
    .groups = "drop"
  ) %>%
  unnest(res)

pseudo_kd_comp_cor_res %>%
  select(chembl_type, estimate, n) %>%
  pivot_wider(names_from = chembl_type, values_from = c(estimate, n)) %>%
  with(
    two_cor_test(
      estimate_ic50, estimate_kd,
      n_ic50, n_kd
    )
  )

```

Check false negatives and false positives. Any particular compounds or
targets that are consistently discordant?

Use Kd data

```{r}
pseudo_kd_comp_false_class <- pseudo_kd_comp_both %>%
  filter(
    pseudo_kd_relation == "=" | (
      pseudo_kd_relation == "<" & pseudo_kd < 90
    ) | (
      pseudo_kd_relation == ">" & pseudo_kd > 9000
    )
  ) %>%
  mutate(
    false_class = case_when(
      (value > (10 * pseudo_kd)) & pseudo_kd != 10000 ~ "false positive",
      (value < (0.1 * pseudo_kd)) & pseudo_kd != 12.5 ~ "false negative",
      TRUE ~ "concordant"
    )
  )

write_csv(
  pseudo_kd_comp_false_class,
  "classification/pseudo_kd_chembl_comparison_false_class.csv"
)

pseudo_kd_comp_false_class %>%
  count(chembl_type, false_class)

pseudo_kd_comp_false_class_by_compound <- pseudo_kd_comp_false_class %>%
  count(chembl_type, lspci_id, name, false_class) %>%
  group_by(chembl_type, lspci_id, name) %>%
  filter(sum(n) > 50) %>%
  mutate(
    prop = n / sum(n)
  ) %>%
  ungroup()

pseudo_kd_comp_false_wide_by_compound <- pseudo_kd_comp_false_class_by_compound %>%
  pivot_wider(names_from = false_class, values_from = c(n, prop), values_fill = list(n = 0L, prop = 0))

View(
  pseudo_kd_comp_false_wide_by_compound %>%
    filter(chembl_type == "kd") %>%
    arrange(desc(`prop_false negative` + `prop_false positive`))
)

pseudo_kd_comp_false_class_by_kinase <- pseudo_kd_comp_false_class %>%
  count(chembl_type, `DiscoveRx Gene Symbol`, gene_symbol, false_class) %>%
  group_by(chembl_type, `DiscoveRx Gene Symbol`, gene_symbol) %>%
  filter(sum(n) > 50) %>%
  mutate(
    prop = n / sum(n)
  ) %>%
  ungroup()

pseudo_kd_comp_false_wide_by_kinase <- pseudo_kd_comp_false_class_by_kinase %>%
  pivot_wider(names_from = false_class, values_from = c(n, prop), values_fill = list(n = 0L, prop = 0))

View(
  pseudo_kd_comp_false_wide_by_kinase %>%
    filter(chembl_type == "kd") %>%
    arrange(desc(`prop_false negative` + `prop_false positive`))
)

pseudo_kd_comp_false_wide_by_kinase %>%
    filter(chembl_type == "kd") %>%
    arrange(desc(`prop_false negative`)) %>%
    head(n = 20) %>%
    select(-chembl_type, -gene_symbol) %>%
    print(n = Inf)

pseudo_kd_comp_false_wide_by_kinase %>%
    filter(chembl_type == "kd") %>%
    arrange(desc(`prop_false positive`)) %>%
    head(n = 20) %>%
    select(-chembl_type, -gene_symbol) %>%
    print(n = Inf)
```

Plot common false positives/negatives

```{r}
pseudo_kd_comp_false_cherrypicked <- pseudo_kd_comp_false_class_by_kinase %>%
  filter(chembl_type == "kd", false_class != "concordant") %>%
  group_by(false_class) %>%
  arrange(desc(n)) %>%
  slice_head(n = 5) %>%
  ungroup()

library(ggrepel)
p <- pseudo_kd_kd_comp %>%
  drop_na(pseudo_kd) %>%
  filter(
    # pseudo_ic50_relation == "=",
    pseudo_kd_concordant == "concordant",
    classification != "discordant"
  ) %>%
  ggplot(
    # aes(value, pseudo_ic50, color = classification, shape = pseudo_ic50_relation)
    aes(chembl_kd, pseudo_kd, fill = classification)
  ) +
  scale_fill_manual(values = map_chr(class_colors, \(x) alpha(x, alpha = 0.5))) +
  geom_point(
    aes(color = false_class, size = false_class),
    data = \(x) x %>%
      left_join(
        pseudo_kd_comp_false_class %>%
          filter(
            chembl_type == "kd"
          ) %>%
          semi_join(
            pseudo_kd_comp_false_cherrypicked,
            by = c("DiscoveRx Gene Symbol", "false_class")
          ) %>%
          transmute(
            lspci_id, `DiscoveRx Gene Symbol`,
            false_class
          ),
        by = c("lspci_id", "DiscoveRx Gene Symbol")
      ) %>%
      arrange(desc(is.na(false_class))),
    shape = 21
  ) +
  scale_color_manual(
    values = c(`false negative` = alpha("black", .8), `false positive` = alpha("black", .8)),
    na.value = alpha("purple", 0), guide = "none"
  ) +
  scale_size_manual(
    values = c(`false negative` = 2, `false positive` = 2),
    na.value = 1, guide = "none"
  ) +
  # geom_label_repel(
  #   aes(label = label),
  #   data = \(x) x %>%
  #     left_join(
  #       pseudo_kd_comp_false_class %>%
  #         filter(
  #           chembl_type == "kd"
  #         ) %>%
  #         semi_join(
  #           pseudo_kd_comp_false_cherrypicked,
  #           by = c("DiscoveRx Gene Symbol", "false_class")
  #         ) %>%
  #         transmute(
  #           lspci_id, `DiscoveRx Gene Symbol`,
  #           label = `DiscoveRx Gene Symbol`
  #         ),
  #       by = c("lspci_id", "DiscoveRx Gene Symbol")
  #     ) %>%
  #     mutate(label = replace_na(label, "")),
  #     show.legend = FALSE
  # ) +
  # Limiting to between 1 and 10000 because barely any data outside this range
  scale_x_log10(labels = scales::label_number(big.mark = ""), limits = c(1, 10000)) +
  scale_y_log10(labels = scales::label_number(big.mark = "")) +
  coord_equal() +
  labs(
    x = "Kd (ChEMBL)",
    y = "Pseudo-Kd (KINOMEscan)"
  )

ggsave(
    "plots/chembl_kd_vs_kinomescan_pseudo_kd_top_false.pdf", p,
    width = 10, height = 6
)

```


We can identify some kinases which seem to have disproportionately many false
negatives and positives. However, it turns out that most of them, e.g. KIT-1, FLT-3,
and JAK3, are measured in multiple variants in Kinomescan, such as autoinhibited,
phosphorylated, or multidomain variants. Therefore, we think that the Kinomescan
"false negatives" might actually be correct.

In order to pick false negatives and positives that don't have an easy explanation
we can exclude kinases which are measured in multiple variants in Kinomescan.

Manually adding EIF2AK4, JAK2, JAK3, because they only have one variant in Kinomescan
but they actually contain multiple kinase domains.

In the ChEMBL comparison, using the `exclude_target` column, we are already
removing any data points from mutants in Kinomescan.

```{r}
kinase_n_variants <- kinomescan_kinase_info %>%
  filter(`Kinase Form` == "Wild Type") %>%
  count(hgnc_symbol)

kinase_multiple_variants <- kinase_n_variants %>%
  filter(n > 1) %>%
  pull(hgnc_symbol) %>%
  c(c("EIF2AK4", "JAK2", "JAK3"))

kinomescan_kinase_info %>%
  filter(
    `Kinase Form` == "Wild Type",
    hgnc_symbol %in% kinase_multiple_variants
  )

pseudo_kd_comp_false_class_by_kinase_no_var <- pseudo_kd_comp_false_class %>%
  filter(!gene_symbol %in% kinase_multiple_variants) %>%
  count(chembl_type, `DiscoveRx Gene Symbol`, gene_symbol, false_class) %>%
  group_by(chembl_type, `DiscoveRx Gene Symbol`, gene_symbol) %>%
  filter(sum(n) > 50) %>%
  mutate(
    prop = n / sum(n)
  ) %>%
  ungroup()

pseudo_kd_comp_false_wide_by_kinase_no_var <- pseudo_kd_comp_false_class_by_kinase_no_var %>%
  pivot_wider(names_from = false_class, values_from = c(n, prop), values_fill = list(n = 0L, prop = 0))

View(
  pseudo_kd_comp_false_wide_by_kinase_no_var %>%
    filter(chembl_type == "kd") %>%
    arrange(desc(`prop_false negative` + `prop_false positive`))
)
```

Visualize which false positive/negative outliers are proteins with variants
in Kinomescan. Marking them large and with black outline

```{r}
p <- pseudo_kd_kd_comp %>%
  drop_na(pseudo_kd) %>%
  filter(
    # pseudo_ic50_relation == "=",
    pseudo_kd_concordant == "concordant",
    classification != "discordant"
  ) %>%
  ggplot(
    # aes(value, pseudo_ic50, color = classification, shape = pseudo_ic50_relation)
    aes(chembl_kd, pseudo_kd, fill = classification)
  ) +
  scale_fill_manual(values = map_chr(class_colors, \(x) alpha(x, alpha = 0.5))) +
  geom_point(
    aes(color = false_class, size = false_class),
    data = \(x) x %>%
      left_join(
        pseudo_kd_comp_false_class %>%
          filter(
            chembl_type == "kd",
            gene_symbol %in% kinase_multiple_variants
          ) %>%
          transmute(
            lspci_id, `DiscoveRx Gene Symbol`,
            false_class
          ),
        by = c("lspci_id", "DiscoveRx Gene Symbol")
      ) %>%
      arrange(desc(is.na(false_class))),
    shape = 21
  ) +
  scale_color_manual(
    values = c(`false negative` = alpha("black", .8), `false positive` = alpha("black", .8)),
    na.value = alpha("purple", 0), guide = "none"
  ) +
  scale_size_manual(
    values = c(`false negative` = 2, `false positive` = 2),
    na.value = 1, guide = "none"
  ) +
  # geom_label_repel(
  #   aes(label = label),
  #   data = \(x) x %>%
  #     left_join(
  #       pseudo_kd_comp_false_class %>%
  #         filter(
  #           chembl_type == "kd"
  #         ) %>%
  #         semi_join(
  #           pseudo_kd_comp_false_cherrypicked,
  #           by = c("DiscoveRx Gene Symbol", "false_class")
  #         ) %>%
  #         transmute(
  #           lspci_id, `DiscoveRx Gene Symbol`,
  #           label = `DiscoveRx Gene Symbol`
  #         ),
  #       by = c("lspci_id", "DiscoveRx Gene Symbol")
  #     ) %>%
  #     mutate(label = replace_na(label, "")),
  #     show.legend = FALSE
  # ) +
  # Limiting to between 1 and 10000 because barely any data outside this range
  scale_x_log10(labels = scales::label_number(big.mark = ""), limits = c(1, 10000)) +
  scale_y_log10(labels = scales::label_number(big.mark = "")) +
  coord_equal() +
  labs(
    x = "Kd (ChEMBL)",
    y = "Pseudo-Kd (KINOMEscan)"
  )

ggsave(
    "plots/chembl_kd_vs_kinomescan_pseudo_kd_false_multiple_variants_new.pdf", p,
    width = 10, height = 6
)

```

What proportion of false positives/negatives are due to multiple variants?

```{r}
pseudo_kd_comp_false_class_with_var <- pseudo_kd_comp_false_class %>%
  mutate(
    multiple_variants = gene_symbol %in% kinase_multiple_variants
  )

pseudo_kd_comp_false_class_with_var_xtab <- pseudo_kd_comp_false_class_with_var %>%
  count(false_class, multiple_variants) %>%
  pivot_wider(names_from = multiple_variants, values_from = n)

pseudo_kd_comp_false_class_with_var_xtab %>%
  column_to_rownames("false_class") %>%
  rstatix::pairwise_fisher_test(detailed = TRUE) %>%
  select(group1, group2, estimate, p)

library(janitor)
pseudo_kd_comp_false_class_with_var_xtab %>%
  as_tabyl(col_var_name = "multiple_variants", row_var_name = "false_class") %>%
  adorn_percentages("row") %>%
  adorn_pct_formatting(digits = 1) %>%
  adorn_ns()
```

  group1         group2         estimate         p
  <chr>          <chr>             <dbl>     <dbl>
1 concordant     false negative    5.50  2.87e-117
2 concordant     false positive    1.46  3.93e- 13
3 false negative false positive    0.266 2.54e- 58

False negative and positive group definitely overrepresented in kinases with
multiple variants.

Select some false positives and negatives that are not due to multiple variants
for further analysis, perhaps for collecting full IC50s from Reaction Biology.

```{r}
pseudo_kd_comp_for_verification <- pseudo_kd_comp_false_class_with_var %>%
  filter(!multiple_variants, chembl_type == "kd") %>%
  mutate(
    lfc = log2(pseudo_kd / value)
  ) %>%
  arrange(desc(abs(lfc)))

View(pseudo_kd_comp_for_verification)
```

Double check if the dose-resonse curves for these look geom_pointdensity

```{r}
kinomescan_classification_raw <- syn("syn52428590") %>%
  qs::qread()

pseudo_kd_comp_for_verification_hm <- pseudo_kd_comp_for_verification %>%
  filter(
    false_class == "false negative"
  ) %>%
  head(n = 100) %>%
  inner_join(
    kinomescan_classification_raw %>%
      filter(dataset == "original_repeat_replaced") %>%
      select(lspci_id, `DiscoveRx Gene Symbol`, data),
    by = c("lspci_id", "DiscoveRx Gene Symbol")
  ) %>%
  transmute(
    combo = paste(name, `DiscoveRx Gene Symbol`, sep = " "), data
  ) %>%
  unnest(data) %>%
  filter(`Compound Concentration (nM)` %in% c(12.5, 100, 1000, 10000)) %>%
  select(combo, `Compound Concentration (nM)`, `Percent Control`) %>%
  pivot_wider(names_from = `Compound Concentration (nM)`, values_from = `Percent Control`)

p <- pheatmap::pheatmap(
  pseudo_kd_comp_for_verification_hm %>%
    column_to_rownames("combo"),
  cluster_cols = FALSE,
  color = rev(viridisLite::viridis(100)),
  silent = TRUE
)

ggsave(
  "plots/chembl_kd_vs_kinomescan_pseudo_kd_false_negative_top100_heatmap.pdf",
  p, width = 8, height = 12
)

pseudo_kd_comp_for_verification_hm <- pseudo_kd_comp_for_verification %>%
  filter(
    false_class == "false positive"
  ) %>%
  head(n = 100) %>%
  inner_join(
    kinomescan_classification_raw %>%
      filter(dataset == "original_repeat_replaced") %>%
      select(lspci_id, `DiscoveRx Gene Symbol`, data),
    by = c("lspci_id", "DiscoveRx Gene Symbol")
  ) %>%
  transmute(
    combo = paste(name, `DiscoveRx Gene Symbol`, sep = " "), data
  ) %>%
  unnest(data) %>%
  filter(`Compound Concentration (nM)` %in% c(12.5, 100, 1000, 10000)) %>%
  select(combo, `Compound Concentration (nM)`, `Percent Control`) %>%
  pivot_wider(names_from = `Compound Concentration (nM)`, values_from = `Percent Control`)

p <- pheatmap::pheatmap(
  pseudo_kd_comp_for_verification_hm %>%
    column_to_rownames("combo"),
  cluster_cols = FALSE,
  color = rev(viridisLite::viridis(100)),
  silent = TRUE
)
ggsave(
  "plots/chembl_kd_vs_kinomescan_pseudo_kd_false_positive_top100_heatmap.pdf",
  p, width = 8, height = 12
)

```

```{r}
pseudo_ic50_comp_false_wide_brms_input <- pseudo_ic50_comp_false_wide_by_compound %>%
  mutate(
    n_total = concordant + `false positive` + `false negative`
  ) %>%
  transmute(
    chembl_type,
    lspci_id, name,
    n_total,
    y = data.frame(
      concordant = concordant,
      `false positive` = `false positive`,
      `false negative` = `false negative`
    ) %>%
      as.matrix()
  )

library(brms)
mdl <- brm(
  y | trials(n_total) ~ 1 + (1 | name),
  data = pseudo_ic50_comp_false_wide_brms_input %>%
    filter(chembl_type == "kd"),
  family = multinomial(refcat = "concordant"),
  cores = 4,
  chains = 4
)

mdl2 <- brm(
  y | trials(n_total) ~ 1 + (1 | name) + chembl_type,
  data = pseudo_ic50_comp_false_wide_brms_input,
  family = multinomial(refcat = "concordant"),
  cores = 4,
  chains = 4
)

pseudo_ic50_comp_false_fisher_by_compound <- pseudo_ic50_comp_false_wide_by_compound %>%
  mutate(
    n_total = concordant + `false positive` + `false negative`
  ) %>%
  pivot_longer(
    c(`false positive`, `false negative`),
    names_to = "false_class",
    values_to = "n_class"
  ) %>%
  mutate(
    n_other = n_total - n_class
  ) %>%
  group_by(
    chembl_type, false_class
  ) %>%
  mutate(
    n_rest_class = sum(n_class) - n_class,
    n_rest_other = sum(n_other) - n_other
  ) %>%
  ungroup()



pseudo_ic50_comp_false_fisher_res_by_compound <- pseudo_ic50_comp_false_fisher_by_compound %>%
  group_by(chembl_type, lspci_id, name) %>%
  summarize(
    res = fisher.test(
      matrix(
        c(n_class, n_rest_class, n_other, n_rest_other),
        ncol = 2
      ),
      simulate = TRUE
    ) %>%
      broom::tidy() %>%
      list(),
    .groups = "drop"
  ) %>%
  unnest(res)

library(brms)
data(HairEyeColor)
haireye <- as.data.frame(HairEyeColor)
names(haireye) <- tolower(names(haireye))
names(haireye)[names(haireye) == "freq"] <- "number"
haireye
he_formula <- brmsformula(number ~ sex + hair + eye + sex:hair + sex:eye + hair:eye)
get_prior(he_formula, data=haireye, family='poisson')


pseudo_ic50_comp_false_class_by_kinase <- pseudo_ic50_comp_false_class %>%
  count(chembl_type, `DiscoveRx Gene Symbol`, false_class) %>%
  group_by(chembl_type, `DiscoveRx Gene Symbol`) %>%
  mutate(
    prop = n / sum(n)
  )
```

```{r}
pseudo_ic50_comp %>%
  filter(
    pseudo_ic50_concordant == "concordant",
    # pseudo_ic50_relation == "=",
    classification != "discordant"
  ) %>%
  group_by(classification) %>%
  mutate(ic50 = value) %>%
  summarize(
    across(
      c(pseudo_ic50, ic50),
      list(
        # average = ~10^(mean(log10(.x))),
        # ci = ~list(10^DescTools::MeanCI(log10(.x)))
        ci = ~list(10^quantile(log10(.x), c(.05, .5, .95), na.rm = TRUE))
      )
    ),
    .groups = "drop"
  ) %>%
  unnest_wider(ends_with("ci"), names_sep = "_")


pseudo_ic50_comp_left <- pseudo_ic50_class %>%
  group_by(
    lspci_id, hmsl_id, `Entrez Gene Symbol`, lspci_target_id
  ) %>%
  arrange(classification == "discordant", pseudo_ic50_relation != "=", pseudo_ic50) %>%
  slice_head(n = 1) %>%
  ungroup() %>%
  left_join(
    lsp_biochem,
    by = c("lspci_id", "lspci_target_id")
  )

pseudo_ic50_comp_left %>%
  mutate(
    chembl_is_covered = !is.na(value)
  ) %>%
  count(chembl_is_covered)


```

```{r ic50_single_dose_comp}
single_dose_trans <- single_dose %>%
  mutate(
    log_inhibition = log10(`Percent Control` + 0.01),
  ) %>%
  as_tibble()


single_dose_comp <- single_dose_trans %>%
  # power_left_join(
  #   lsp_targets %>%
  #     filter(organism == "Homo sapiens") %>%
  #     select(lspci_target_id, symbol),
  #   by = c("Entrez Gene Symbol" = "symbol"),
  #   check = check_specs(
  #     duplicate_keys_left = "ignore",
  #     duplicate_keys_right = "ignore",
  #     unmatched_keys_left = "warn",
  #     unmatched_keys_right = "ignore"
  #   )
  # )  %>%
  filter(dataset == "original_repeat_replaced") %>%
  power_left_join(
    kinomescan_classes %>%
      filter(dataset == "original_repeat_replaced") %>%
      distinct(lspci_id, `DiscoveRx Gene Symbol`, classification) %>%
      mutate(across(classification, \(x) factor(x, levels = names(class_colors)))),
    by = c("lspci_id", "DiscoveRx Gene Symbol"),
    check = check_specs(
      duplicate_keys_left = "ignore",
      duplicate_keys_right = "warn",
      unmatched_keys_left = "warn",
      unmatched_keys_right = "warn"
    )
  ) %>%
  inner_join(
    chembl_kd %>%
      select(lspci_id, chembl_kd, chembl_kd_relation, `DiscoveRx Gene Symbol`),
    by = c("lspci_id", "DiscoveRx Gene Symbol")
  )

p <- single_dose_comp %>%
  mutate(
    percent_control_zero = `Percent Control` == 0,
    percent_control_trans = if_else(
      percent_control_zero,
      0.5 * min(`Percent Control`[`Percent Control` > 0]),
      `Percent Control`
    ),
    chembl_kd = pmin(chembl_kd, 10000)
  ) %>%
  group_by(`Compound Concentration (nM)`) %>%
  filter(n() > 1000) %>%
  ungroup() %>%
  filter(
    !qpcr_miss,
    !exclude_target,
    classification != "discordant"
  ) %>%
  ggplot(
    aes(chembl_kd, `Percent Control`)
  ) +
  # ggrastr::rasterize(geom_point(alpha = 0.5)) +
  ggrastr::rasterize(
    ggpointdensity::geom_pointdensity(alpha = 0.5, shape = 16, adjust = 2),
    dev = "ragg",
    dpi = 300
  ) +
  # scale_shape_manual(values = c(`TRUE` = 4, `FALSE` = 16)) +
  scale_x_log10(breaks = c(0.1, 1, 10, 100, 1000, 10000), labels = scales::label_number()) +
  scale_color_viridis_c(trans = scales::log10_trans(), option = "magma") +
  facet_grid(classification ~ `Compound Concentration (nM)`)

ggsave(
  "plots/single_dose_vs_chembl_kd_all_doses_new.pdf", p,
  width = 10, height = 10
)

density_2d <- function (x, y = NULL, nbin = 128, bandwidth) {
  xy <- xy.coords(x, y, setLab = FALSE)
  select <- is.finite(xy$x) & is.finite(xy$y)
  x <- cbind(xy$x, xy$y)[select, ]
  map <- grDevices:::.smoothScatterCalcDensity(x, nbin, bandwidth)
  mkBreaks <- function(u) u - diff(range(u))/(length(u) - 1)/2
  xbin <- cut(x[, 1], mkBreaks(map$x1), labels = FALSE)
  ybin <- cut(x[, 2], mkBreaks(map$x2), labels = FALSE)
  dens <- map$fhat[cbind(xbin, ybin)]
  dens[is.na(dens)] <- 0
  dens
}

ps <- single_dose_comp %>%
  mutate(
    percent_control_zero = `Percent Control` == 0,
    percent_control_trans = if_else(
      percent_control_zero,
      0.5 * min(`Percent Control`[`Percent Control` > 0]),
      `Percent Control`
    ),
    chembl_kd = pmin(chembl_kd, 10000)
  ) %>%
  filter(
    !qpcr_miss,
    !exclude_target,
    classification != "discordant"
  ) %>%
  group_by(`Compound Concentration (nM)`) %>%
  summarize(
    plot = list(
      cur_data() %>%
        mutate(
          density = density_2d(
            log10(chembl_kd), `Percent Control`
          )
        ) %>%
        arrange(density) %>%
      # ggplot(cur_data(), aes(value, percent_control_trans, shape = percent_control_zero)) +
      ggplot(aes(chembl_kd, `Percent Control`, color = density)) +
        geom_point(alpha = 0.8, shape = 16) +
        scale_color_viridis_c(
          trans = "log10",
          breaks = scales::log_breaks(base = 10),
          labels = scales::trans_format("log10", scales::math_format(10^.x))
        ) +
        # scale_color_manual(values = class_colors) +
        # scale_shape_manual(
        #   values = c(`TRUE` = 25, `FALSE` = 16)
        # ) +
        scale_x_log10(breaks = c(0.1, 1, 10, 100, 1000, 10000), labels = scales::label_number()) +
        # scale_y_log10() +
        guides(shape = "none") +
        labs(
          x = "Kd ChEMBL (nM)",
          y = "% control Eurofins",
          color = "Density",
          title = paste("Compound concentration", cur_group()$`Compound Concentration (nM)`, "nM")
        )
    )
  )

pwalk(
  ps,
  function(`Compound Concentration (nM)`, plot, ...) {
    ggsave(
      paste0("plots/single_dose_vs_chembl_kd_", `Compound Concentration (nM)`, ".pdf"),
      plot,
      width = 5, height = 3.5
    )
  }
)


ps <- single_dose_comp %>%
  mutate(
    percent_control_zero = `Percent Control` == 0,
    percent_control_trans = if_else(
      percent_control_zero,
      0.5 * min(`Percent Control`[`Percent Control` > 0]),
      `Percent Control`
    ),
    value = pmin(value, 30000)
  ) %>%
  filter(
    !qpcr_miss
  ) %>%
  group_by(`Compound Concentration (nM)`) %>%
  summarize(
    plot = list(
      # ggplot(cur_data(), aes(value, percent_control_trans, shape = percent_control_zero)) +
      ggplot(cur_data(), aes(value, `Percent Control`)) +
        geom_point(alpha = 0.5) +
        # scale_shape_manual(
        #   values = c(`TRUE` = 25, `FALSE` = 16)
        # ) +
        scale_x_log10(breaks = c(0.1, 1, 10, 100, 1000, 10000), labels = scales::label_number()) +
        # scale_y_log10() +
        guides(shape = "none") +
        labs(
          x = "IC50 ChEMBL (nM)",
          y = "% control Eurofins",
          title = paste("Compound concentration", cur_group()$`Compound Concentration (nM)`, "nM")
        )
    )
  )

pwalk(
  ps,
  function(`Compound Concentration (nM)`, plot, ...) {
    ggsave(
      paste0("plots/single_dose_vs_chembl_ic50_", `Compound Concentration (nM)`, "nM.pdf"),
      plot,
      width = 8, height = 8
    )
  }
)
```

Testing if there are any kinases or compounds that are particularly discordant
compared to the ChEBML data.

First, establishing the Percent Control thresholds necessary to ensure a 95%
probability that the ChEMBL IC50 is below 5000 nM.

To do that, first plot for each possible Percent Control threshold what proportion
of IC50s is below 5000 nM.


```{r}
percent_control_threshold_plotting_data <- single_dose_comp %>%
  group_by(`Compound Concentration (nM)`) %>%
  filter(n() > 1000) %>%
  ungroup() %>%
  filter(
    classification != "discordant",
    !qpcr_miss
  ) %>%
  # crossing(ic_50_threshold = 10**seq(2, 4, length.out = 6)) %>%
  crossing(ic_50_threshold = c(100, 200, 500, 1000, 5000, 10000)) %>%
  group_by(`Compound Concentration (nM)`, ic_50_threshold) %>%
  arrange(`Percent Control`) %>%
  mutate(
    cumsum_below_threshold = cumsum(value < ic_50_threshold),
    frac_below_threshold = cumsum_below_threshold / seq_len(n())
  ) %>%
  ungroup() %>%
  group_by(`Compound Concentration (nM)`, `Percent Control`, ic_50_threshold) %>%
  slice_tail(n = 1) %>%
  ungroup()

p <- percent_control_threshold_plotting_data %>%
  mutate(across(c(ic_50_threshold, `Compound Concentration (nM)`), ~as.character(.x) %>% fct_inseq(ordered = TRUE))) %>%
  ggplot(aes(`Percent Control`, frac_below_threshold, color = `Compound Concentration (nM)`)) +
  geom_line() +
  # geom_point() +
  facet_wrap(~ic_50_threshold) +
  labs(y = "Fraction of measurements below IC50 threshold")

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

```


```{r}
pseudo_ic50_comp_stats <- pseudo_ic50_comp %>%
  drop_na(pseudo_ic50) %>%
  filter(classification != "discordant") %>%
  mutate(
      log10_ic50 = log10(value),
      log10_pseudo_ic50 = log10(pseudo_ic50),
      log_ic50_diff = log10_ic50 - log10_pseudo_ic50,
      diff_greater_than_1 = abs(log_ic50_diff) > 1,
      difference_class = case_when(
          log_ic50_diff > 1 ~ "FP",
          log_ic50_diff < -1 ~ "FN",
          TRUE ~ "concordant"
      )
  )
```


```{r}
stats_per_group <- function(grouped_df) {
    grouped_df %>%
    summarize(
        n = n(),
        n_binding_class = sum(classification == "binding"),
        frac_binding_class = n_binding_class / n,
        n_diff_greater_than_1 = sum(diff_greater_than_1, na.rm = TRUE),
        frac_diff_greater_than_1 = n_diff_greater_than_1 / n,
        mean_log_ic50_diff = mean(log_ic50_diff, na.rm = TRUE),
        sd_log_ic50_diff = sd(log_ic50_diff, na.rm = TRUE),
        n_false_positives = sum(difference_class == "FP", na.rm = TRUE),
        n_false_negatives = sum(difference_class == "FN", na.rm = TRUE),
        frac_false_positives = n_false_positives / n,
        frac_false_negatives = n_false_negatives / n,
        .groups = "drop"
    )
}

pseudo_ic50_comp_stats_per_compound <- pseudo_ic50_comp_stats %>%
    group_by(lspci_id, name) %>%
    filter(n() > 10) %>%
    stats_per_group()
```

```{r}
library(ggrepel)
p <- pseudo_ic50_comp_stats_per_compound %>%
  mutate(
    outside_prediction_interval = {
      mdl <- lm(
        frac_diff_greater_than_1 ~ frac_binding_class,
        data = .
      )
      pred_interval <- predict(mdl, interval = "prediction", level = 0.8)
      frac_diff_greater_than_1 < pred_interval[, "lwr"] | frac_diff_greater_than_1 > pred_interval[, "upr"]
    }
  ) %>%
  ggplot(
    aes(
      frac_binding_class, frac_diff_greater_than_1,
      text = name
    )
  ) +
    geom_point(shape = 16) +
    geom_text_repel(
      aes(label = name),
      data = ~.x %>%
        mutate(name = if_else(!outside_prediction_interval, "", name)),
      # min.segment.length = 0,
      max.overlaps = Inf,
    ) +
    guides(color = "none") +
    labs(
        x = "Fraction of 'binding' dose-response curves",
        y = "Fraction of pseudo IC50s with >10x difference",
        title = "Per compound"
    )


plotly::ggplotly(p)

ggsave(
    "plots/eurofins_ic50_vs_kinomescan_pseudo_ic50_differences_per_compound.pdf", p,
    width = 10, height = 7
)
```

Now do the same per kinase

```{r}

pseudo_ic50_comp_stats_per_kinase <- pseudo_ic50_comp_stats %>%
    group_by(`DiscoveRx Gene Symbol`) %>%
    filter(n() > 10) %>%
    stats_per_group() %>%
    inner_join(
      kinomescan_kinase_info
    )
```

```{r}
lm(
  frac_diff_greater_than_1 ~ frac_binding_class,
  data = pseudo_ic50_comp_stats_per_kinase
) %>%
  summary()

cor.test(
  ~frac_diff_greater_than_1 + frac_binding_class,
  data = pseudo_ic50_comp_stats_per_kinase
)

p <- pseudo_ic50_comp_stats_per_kinase %>%
  mutate(
    outside_prediction_interval = {
      mdl <- lm(
        frac_diff_greater_than_1 ~ frac_binding_class,
        data = .
      )
      pred_interval <- predict(mdl, interval = "prediction", level = 0.8)
      frac_diff_greater_than_1 < pred_interval[, "lwr"] | frac_diff_greater_than_1 > pred_interval[, "upr"]
    }
  ) %>%
  ggplot(
    aes(
      frac_binding_class, frac_diff_greater_than_1,
      color = `Kinase Form` == "Wild Type",
      text = `DiscoveRx Gene Symbol`
    )
  ) +
    geom_point(shape = 16) +
    geom_text_repel(
      aes(label = `DiscoveRx Gene Symbol`),
      data = ~.x %>%
        mutate(`DiscoveRx Gene Symbol` = if_else(!outside_prediction_interval, "", `DiscoveRx Gene Symbol`)),
      # min.segment.length = 0,
      max.overlaps = Inf,
      box.padding = 0.5,
      max.iter = 50000, max.time = 50,
      force = 1.5,
      seed = 42
    ) +
    guides(color = "none") +
    labs(
        x = "Fraction of 'binding' dose-response curves",
        y = "Fraction of pseudo IC50s with >10x difference",
        title = "Per kinase"
    )

plotly::ggplotly(p)

ggsave(
    "plots/eurofins_ic50_vs_kinomescan_pseudo_ic50_differences_per_kinase.pdf", p,
    width = 10, height = 7
)
```

Establish confidence intervals for IC50 values given the Eurofins data

Doing it in bins of %control values as well as cumulatively

```{r}
single_dose_comp_binned <- single_dose_comp %>%
  filter(dataset == "original_repeat_replaced") %>%
  mutate(
    percent_control_bin = cut(
      `Percent Control`,
      breaks = c(-Inf, seq(from = 10, to = 90, by = 10), Inf)
    ),
    kd_bin = cut(
      chembl_kd,
      breaks = c(-Inf, 10**seq(-1, 4), Inf)
    )
  )

percent_control_ic50_ci <- single_dose_comp_binned %>%
  group_by(`Compound Concentration (nM)`) %>%
  filter(n() > 1000) %>%
  ungroup() %>%
  filter(
    classification != "discordant",
    !qpcr_miss
  ) %>%
  group_by(`Compound Concentration (nM)`, percent_control_bin) %>%
  summarize(
    n = n(),
    ci = quantile(value, probs = c(0.05, .5, 0.95)) %>%
      enframe("quantile", "value") %>%
      # pivot_wider(names_from = quantile, values_from = value, names_prefix = "ci_") %>%
      list(),
    .groups = "drop"
  ) %>%
  unnest(ci)

p <- percent_control_ic50_ci %>%
  ggplot(
    aes(percent_control_bin, value, color = quantile, group = quantile)
  ) +
  geom_point() +
  geom_line() +
  facet_wrap(~`Compound Concentration (nM)`) +
  scale_y_log10()

percent_control_ic50_ci_ic50 <- single_dose_comp_binned %>%
  group_by(`Compound Concentration (nM)`) %>%
  filter(n() > 1000) %>%
  ungroup() %>%
  filter(
    classification != "discordant",
    !qpcr_miss
  ) %>%
  group_by(`Compound Concentration (nM)`, ic_50_bin) %>%
  summarize(
    n = n(),
    ci = quantile(`Percent Control`, probs = c(0.05, .25, .5, .75, 0.95)) %>%
      enframe("quantile", "value") %>%
      # pivot_wider(names_from = quantile, values_from = value, names_prefix = "ci_") %>%
      list(),
    .groups = "drop"
  ) %>%
  unnest(ci)

p <- percent_control_ic50_ci_ic50 %>%
  ggplot(
    aes(ic_50_bin, value, color = quantile, group = quantile)
  ) +
  geom_point() +
  geom_line() +
  facet_wrap(~`Compound Concentration (nM)`) +
  scale_y_log10()

# library(ggbeeswarm)
# p <- single_dose_comp_binned %>%
#   ggplot(
#     aes(ic_50_bin, value)
#   ) +
#   ggrastr::rasterize(geom_beeswarm()) +
#   facet_wrap(~`Compound Concentration (nM)`) +
#   scale_y_log10()


p <- single_dose_comp_binned %>%
  group_by(`Compound Concentration (nM)`) %>%
  filter(n() > 1000) %>%
  ungroup() %>%
  mutate(
    across(
      `Compound Concentration (nM)`,
      ~ordered(
        .x,
        levels = sort(unique(.x)),
        labels = paste0(sort(unique(.x)), " nM")
      )
    )
  ) %>%
  filter(
    classification != "discordant",
    !qpcr_miss
  ) %>%
  ggplot(
    aes(
      `Percent Control`, value,
      color = after_stat(quantile) %>% {
        factor(
          .,
          levels = sort(unique(.)),
          labels = scales::label_percent()(sort(unique(.)))
        )
      }
      # color = fct_inseq(as.character(after_stat(quantile)), ordered = TRUE)
    )
  ) +
  geom_quantile(
    method = "rqss",
    quantiles = c(.05, .25, .5, .75, .95),
    lambda = 10
  ) +
  # geom_point() +
  # geom_smooth(method = "lm", se = FALSE) +
  scale_y_log10(
    breaks = 10**(0:4)
  ) +
  scale_color_manual(
    values = c(
      "5%" = "red",
      "25%" = "orange",
      "50%" = "black",
      "75%" = "orange",
      "95%" = "red"
    )
  ) +
  labs(
    x = "% control (Eurofins)",
    y = "IC50 (ChEMBL)",
    color = "Quantile"
  ) +
  facet_wrap(~`Compound Concentration (nM)`)

ggsave(
  "plots/percent_control_ic50_ci_quantile_plot.pdf",
  p, width = 8, height = 6
)

p <- single_dose_comp_binned %>%
  group_by(`Compound Concentration (nM)`) %>%
  filter(n() > 1000) %>%
  ungroup() %>%
  mutate(
    across(
      `Compound Concentration (nM)`,
      ~ordered(
        .x,
        levels = sort(unique(.x)),
        labels = paste0(sort(unique(.x)), " nM")
      )
    )
  ) %>%
  filter(
    classification != "discordant",
    !qpcr_miss
  ) %>%
  ggplot(
    aes(
      value, `Percent Control`,
      color = after_stat(quantile) %>% {
        factor(
          .,
          levels = sort(unique(.)),
          labels = scales::label_percent()(sort(unique(.)))
        )
      }
      # color = fct_inseq(as.character(after_stat(quantile)), ordered = TRUE)
    )
  ) +
  geom_quantile(
    method = "rqss",
    quantiles = c(.05, .25, .5, .75, .95),
    lambda = 10
  ) +
  # geom_point() +
  # geom_smooth(method = "lm", se = FALSE) +
  scale_x_log10(
    breaks = 10**(0:4)
  ) +
  scale_color_manual(
    values = c(
      "5%" = "red",
      "25%" = "orange",
      "50%" = "black",
      "75%" = "orange",
      "95%" = "red"
    )
  ) +
  labs(
    x = "IC50 (ChEMBL)",
    y = "% control (Eurofins)",
    color = "Quantile"
  ) +
  facet_wrap(~`Compound Concentration (nM)`)

ggsave(
  "plots/percent_control_ic50_ci_quantile_plot_reverse.pdf",
  p, width = 8, height = 6
)
```

```{r}
KD_THRESHOLDS <- c(10, 100, 1000, 5000)
library(ROCit)

single_dose_comp_binned_filtered <- single_dose_comp_binned %>%
  group_by(`Compound Concentration (nM)`) %>%
  filter(n() > 1000) %>%
  ungroup() %>%
  mutate(
    across(
      `Compound Concentration (nM)`,
      ~ordered(
        .x,
        levels = sort(unique(.x)),
        labels = paste0(sort(unique(.x)), " nM")
      )
    )
  ) %>%
  filter(
    classification != "discordant",
    !qpcr_miss,
    !exclude_target
  )

single_dose_comp_rocs <- single_dose_comp_binned_filtered %>%
  crossing(
    kd_threshold = KD_THRESHOLDS
  ) %>%
  group_by(`Compound Concentration (nM)`, kd_threshold) %>%
  summarize(
    roc = rocit(
      100 - `Percent Control`,
      # `Percent Control`,
      chembl_kd < cur_group()$kd_threshold,
      method = "empirical"
    ) %>%
      list(),
    .groups = "drop"
  ) %>%
  mutate(
    # roc_ci = map(roc, ciROC),
    roc_measures = map(
      roc,
      ~measureit(.x, measure = c("SENS", "SPEC", "PREC", "PPV")) %>%
        unclass() %>%
        as_tibble()
    )
  )


single_dose_comp_roc_values <- single_dose_comp_rocs %>%
  transmute(
    `Compound Concentration (nM)`,
    kd_threshold = fct_inorder(
      paste(kd_threshold, "nM")
    ),
    roc_values = map(roc, plot)
  )

fn_fp_weights <- c("2xFN" = 2/3, "equal" = .5, "2xFP" = 1/3)
compute_weighted_score <- function(df, fn_fp_weight) {
  df %>%
    mutate(
      score = 2 * (fn_fp_weight * (TP / (TP + FN)) + (1 - fn_fp_weight) * (TN / (TN + FP))) - 1,
      score_naive = fn_fp_weight * FN + (1 - fn_fp_weight) * FP
    )
}

single_dose_comp_scores_weighted <- single_dose_comp_rocs %>%
  cross_join(
    enframe(
      fn_fp_weights,
      name = "weight_name",
      value = "fn_fp_weight"
    )
  ) %>%
  mutate(
    scores = map2(roc_measures, fn_fp_weight, compute_weighted_score)
  ) %>%
  select(-roc, -roc_measures)


single_dose_comp_scores_weighted_diagnostics <- single_dose_comp_scores_weighted %>%
  unnest(scores) %>%
  pivot_longer(
    c(TP, TN, FP, FN, score, score_naive, SENS, SPEC, PREC, PPV),
    names_to = "variable",
    values_to = "value"
  ) %>%
  filter(
    case_when(
      variable %in% c("score", "score_naive") ~ TRUE,
      TRUE ~ weight_name == "equal"
    )
  ) %>%
  mutate(
    Cutoff = 100 - Cutoff
  )

single_dose_comp_scores_weighted_diagnostic_optim <- single_dose_comp_scores_weighted_diagnostics %>%
  filter(str_starts(variable, "score")) %>%
  group_nest(
    variable, `Compound Concentration (nM)`, kd_threshold, weight_name
  ) %>%
  mutate(
    data = map2(
      data, variable,
      \(df, variable) {
        if (variable == "score")
          slice_max(df, value)
        else
          slice_min(df, value)
      }
    )
  ) %>%
  unnest(data)

measure_names <- tribble(
  ~abbreviation, ~name, ~name_long,
  "SENS", "Sensitivity", "Sensitivity TP / (TP + FN)",
  "SPEC", "Specificity", "Specificity TN / (TN + FP)",
  "PREC", "Precision", "Precision TP / (TP + FP)",
  "score", "Weighted Youden", "Weighted Youden",
  "score_naive", "Weighted Score", "Weighted Score",
  "TP", "True Positives", "True Positives TP",
  "TN", "True Negatives", "True Negatives TN",
  "FP", "False Positives", "False Positives FP",
  "FN", "False Negatives", "False Negatives FN"
)

library(ggrepel)
p <- single_dose_comp_scores_weighted_diagnostics %>%
  filter(
    variable %in% c("score", "SENS", "SPEC", "PREC")
  ) %>%
  left_join(
    measure_names,
    by = c("variable" = "abbreviation")
  ) %>%
  ggplot(
    aes(
      Cutoff, value, color = name_long,
      linetype = weight_name
    )
  ) +
  geom_line() +
  geom_point(
    data = single_dose_comp_scores_weighted_diagnostic_optim %>%
      filter(variable %in% c("score")) %>%
      left_join(
        measure_names,
        by = c("variable" = "abbreviation")
      ),
    show.legend = FALSE
  ) +
  geom_text_repel(
    aes(label = signif(Cutoff, 2)),
    data = single_dose_comp_scores_weighted_diagnostic_optim %>%
      filter(variable %in% c("score")) %>%
      left_join(
        measure_names,
        by = c("variable" = "abbreviation")
      ),
    size = 3,
    show.legend = FALSE
  ) +
  scale_linetype_manual(
    values = c(
      "2xFN" = "dotted",
      "equal" = "solid",
      "2xFP" = "dashed"
    )
  ) +
  facet_grid(vars(kd_threshold), vars(`Compound Concentration (nM)`)) +
  labs(
    x = "Cutoff", y = "Value",
    color = "Variable", linetype = "Weights"
  )

ggsave(
  file.path("plots", "kd_roc_weighted_scores.pdf"),
  width = 8, height = 5
)

p <- single_dose_comp_scores_weighted_diagnostics %>%
  filter(
    variable %in% c("score_naive", "FP", "FN")
  ) %>%
  left_join(
    measure_names,
    by = c("variable" = "abbreviation")
  ) %>%
  ggplot(
    aes(
      Cutoff, value, color = name_long,
      linetype = weight_name
    )
  ) +
  geom_line() +
  geom_point(
    data = single_dose_comp_scores_weighted_diagnostic_optim %>%
      filter(variable %in% c("score_naive")) %>%
      left_join(
        measure_names,
        by = c("variable" = "abbreviation")
      ),
    show.legend = FALSE
  ) +
  geom_text_repel(
    aes(label = signif(Cutoff, 2)),
    data = single_dose_comp_scores_weighted_diagnostic_optim %>%
      filter(variable %in% c("score_naive")) %>%
      left_join(
        measure_names,
        by = c("variable" = "abbreviation")
      ),
    size = 3,
    show.legend = FALSE
  ) +
  scale_linetype_manual(
    values = c(
      "2xFN" = "dotted",
      "equal" = "solid",
      "2xFP" = "dashed"
    )
  ) +
  scale_y_sqrt() +
  facet_grid(vars(kd_threshold), vars(`Compound Concentration (nM)`)) +
  labs(
    x = "Cutoff", y = "Value",
    color = "Variable", linetype = "Weights"
  )

ggsave(
  file.path("plots", "kd_roc_weighted_score_naive.pdf"),
  width = 8, height = 5
)

```


```{r}
param_order <- c(
  "Cutoff", "FP", "FN", "SENS", "SPEC", "PREC"
)

optimal_params_table <- single_dose_comp_scores_weighted_diagnostics %>%
  inner_join(
    single_dose_comp_scores_weighted_diagnostic_optim %>%
      filter(variable == "score") %>%
      select(-variable, -value)
  ) %>%
  bind_rows(
    single_dose_comp_scores_weighted_diagnostic_optim %>%
      filter(variable == "score") %>%
      transmute(
        `Compound Concentration (nM)`, kd_threshold,
        weight_name,
        value = Cutoff, variable = "Cutoff"
      ) %>%
      distinct()
  ) %>%
  filter(
    weight_name == "equal",
    variable %in% c("Cutoff", "FP", "FN", "SENS", "SPEC", "PREC")
  ) %>%
  select(
    `Compound Concentration (nM)`, kd_threshold,
    variable, value
  ) %>%
  mutate(
    value = case_when(
      variable %in% c("SENS", "SPEC", "PREC") ~ sprintf("%.2f", value),
      TRUE ~ sprintf("%.0f", value)
    ),
    variable = factor(variable, levels = param_order),
    kd_threshold = fct_inseq(as.character(kd_threshold)) %>%
      fct_relabel(\(x) paste0("Target affinity <", x, " nM"))
  ) %>%
  arrange(
    kd_threshold,
    `Compound Concentration (nM)`,
    variable
  ) %>%
  pivot_wider(
    names_from = `Compound Concentration (nM)`,
    values_from = value
  )

library(gt)
optimal_params_gt <- optimal_params_table %>%
  gt(
    groupname_col = "kd_threshold",
    rowname_col = "variable"
  ) %>%
  cols_align(
    align = "right",
    columns = -variable
  ) %>%
  cols_align(
    align = "left",
    columns = variable
  ) %>%
  tab_stub_indent(rows = everything(), indent = 2) %>%
  data_color(
    columns = everything(),
    rows = variable %in% c("SENS", "SPEC", "PREC"),
    direction = "row",
    fn = \(x) {
      # scales::col_numeric(
      #   palette = c("white", "red"),
      #   domain = NULL
      # )(as.numeric(x))
      x_fct <- fct_inseq(x, ordered = TRUE)
      scales::col_factor(
        palette = c("white", "red"),
        domain = levels(x_fct),
        ordered = TRUE
      )(x_fct)
    }
  ) %>%
  data_color(
    columns = everything(),
    rows = variable %in% c("FP", "FN"),
    direction = "row",
    fn = \(x) {
      # scales::col_numeric(
      #   palette = c("red", "white"),
      #   domain = NULL
      # )(as.numeric(x))
      x_fct <- fct_inseq(x, ordered = TRUE)
      scales::col_factor(
        palette = c("red", "white"),
        domain = levels(x_fct),
        ordered = TRUE
      )(x_fct)
    }
  ) %>%
  tab_options(
    table.font.names = "Helvetica"
  )

gtsave(
  optimal_params_gt,
  file.path("plots", "kd_roc_weighted_equal_optimal_table.html")
)
```


```{r}

pwalk(
  single_dose_comp_rocs,
  function(`Compound Concentration (nM)`, ic50_threshold, roc, ...) {
    withr::with_pdf(
      paste0("plots/percent_control_ic50_roc_threshold_", ic50_threshold, "_", `Compound Concentration (nM)`, ".pdf"),
      plot(roc)
    )
  }
)

pwalk(
  single_dose_comp_rocs,
  function(`Compound Concentration (nM)`, ic50_threshold, roc, ...) {
    withr::with_pdf(
      paste0("plots/percent_control_ic50_ks_threshold_", ic50_threshold, "_", `Compound Concentration (nM)`, ".pdf"),
      ksplot(roc)
    )
  }
)

```

Table with FP, FN, AUC for each compound concentration and threshold
highlighted best and worst screening concentration for each

```{r}
confusion_tables <- single_dose_comp_binned_filtered %>%
  inner_join(
    single_dose_comp_roc_values_optimal %>%
      select(
        `Compound Concentration (nM)`,
        ic50_threshold,
        `Cutoff`
      ),
    by = "Compound Concentration (nM)",
    relationship = "many-to-many"
  ) %>%
  group_by(`Compound Concentration (nM)`, ic50_threshold) %>%
  summarize(
    confusion_table = list(
      table(
        chembl_ic50_below = value < as.double(str_replace(ic50_threshold, " nM", "")),
        eurofins_below = `Percent Control` < `Cutoff`
      )
    ),
    .groups = "drop"
  )

library(cvms)

# p <- plot_confusion_matrix(
#   confusion_tables$confusion_table[[1]] %>%
#   as_tibble(),
#   target_col = "chembl_ic50_below",
#   prediction_col = "eurofins_below",
#   counts_col = "n"
# )
# ggsave(
#   "test.pdf", p, width = 2.5, height = 2.5
# )

dir.create("plots/percent_control_ic50_roc_confusion_tables", showWarnings = FALSE)
pwalk(
  confusion_tables,
  function(`Compound Concentration (nM)`, ic50_threshold, confusion_table, ...) {
    p <- plot_confusion_matrix(
      confusion_table %>%
        as_tibble(),
      target_col = "chembl_ic50_below",
      prediction_col = "eurofins_below",
      counts_col = "n"
    ) +
      labs(
        subtitle = paste("Predicting IC50 <", ic50_threshold, "\nat", `Compound Concentration (nM)`)
      ) +
      theme(
        plot.subtitle = element_text(size = 8)
      )
    ggsave(
      paste0("plots/percent_control_ic50_roc_confusion_tables/confusion_table_", ic50_threshold, "_", `Compound Concentration (nM)`, ".pdf"),
      p, width = 2.5, height = 3
    )
  }
)



```

Plot precision-recall curves

```{r}
single_dose_comp_scores_weighted_diagnostic_optim
library(ggrepel)
single_dose_comp_pr_plots <- single_dose_comp_rocs %>%
  transmute(
    `Compound Concentration (nM)`,
    kd_threshold = fct_inorder(
      paste(kd_threshold, "nM")
    ),
    roc_measures
  ) %>%
  unnest(roc_measures) %>%
  mutate(across(Cutoff, ~100 - .x)) %>%
  ggplot(
    aes(
      SENS, PREC,
      color = fct_inorder(
        as.character(`Compound Concentration (nM)`)
      )
    )
  ) +
  geom_line() +
  geom_point(
    data = ~.x %>%
      filter(
        Cutoff %in% seq(from = 0, to = 100, by = 25)
      ) %>%
      bind_rows(
        inner_join(
          .x,
          single_dose_comp_scores_weighted_diagnostic_optim %>%
            filter(variable == "score", weight_name == "equal") %>%
            mutate(
              kd_threshold = fct_inorder(
                paste(kd_threshold, "nM")
              )
            )
        )
      ),
    show.legend = FALSE
  ) +
  geom_text_repel(
    aes(label = Cutoff),
    data = ~.x %>%
      mutate(
        Cutoff = if_else(
          Cutoff %in% seq(from = 25, to = 75, by = 25),
          as.character(Cutoff),
          ""
        )
      ) %>%
      anti_join(
        single_dose_comp_scores_weighted_diagnostic_optim %>%
          filter(variable == "score", weight_name == "equal") %>%
          mutate(across(Cutoff, as.character)) %>%
          mutate(
            kd_threshold = fct_inorder(
              paste(kd_threshold, "nM")
            )
          )
      ) %>%
      bind_rows(
        inner_join(
          .x,
          single_dose_comp_scores_weighted_diagnostic_optim %>%
            filter(variable == "score", weight_name == "equal") %>%
            mutate(
              kd_threshold = fct_inorder(
                paste(kd_threshold, "nM")
              )
            )
        ) %>%
        mutate(
          across(Cutoff, ~paste0(.x, "*"))
        )
      ),
    box.padding = 0.5,
    show.legend = FALSE,
    max.iter = 50000,
    max.time = 2,
    force = 3,
    seed = 42,
    max.overlaps = Inf
  ) +
  labs(
    x = "Recall TP / (TP + FN)",
    y = "Precision TP / (TP + FP)",
    color = "Assay concentration",
    title = "Precision-recall curves for prediction of Kd < threshold"
  ) +
  facet_wrap(~kd_threshold) +
  coord_cartesian(xlim = c(0, 1.1))

ggsave(
  "plots/percent_control_kd_pr.pdf",
  single_dose_comp_pr_plots,
  width = 10, height = 4
)

single_dose_comp_roc_plot <- single_dose_comp_roc_values %>%
  mutate(
    across(
      roc_values,
      map,
      ~.x[c("FPR/1-Specificity", "TPR/Sensitivity/Recall", "Cutoff")] %>%
        as_tibble()
    )
  ) %>%
  unnest(roc_values) %>%
  mutate(across(Cutoff, ~100 - .x)) %>%
  ggplot(
    aes(
      `FPR/1-Specificity`, `TPR/Sensitivity/Recall`,
      color = fct_inorder(
        as.character(`Compound Concentration (nM)`)
      )
    )
  ) +
  geom_line() +
  geom_point(
    data = ~.x %>%
      filter(
        Cutoff %in% seq(from = 0, to = 100, by = 25)
      ) %>%
      bind_rows(
        single_dose_comp_roc_values_optimal
      ),
    show.legend = FALSE
  ) +
  geom_text_repel(
    aes(label = Cutoff),
    data = ~.x %>%
      mutate(
        Cutoff = if_else(
          Cutoff %in% seq(from = 25, to = 75, by = 25),
          as.character(Cutoff),
          ""
        )
      ) %>%
      # If optimal value is equal to one of the plotted
      # values don't repeat it
      anti_join(
        single_dose_comp_roc_values_optimal %>%
          mutate(across(Cutoff, as.character))
      ) %>%
      bind_rows(
        single_dose_comp_roc_values_optimal %>%
          mutate(across(Cutoff, ~paste0(.x, "*")))
      ),
    # point.padding = unit(20, "pt"),
    # min.segment.length = 0,
    box.padding = 0.5,
    show.legend = FALSE,
    max.iter = 50000,
    max.time = 2,
    force = 3,
    seed = 42,
    max.overlaps = Inf
  ) +
  labs(
    x = "1 - Specificity TN / (TN + FP)",
    y = "Sensitivity TP / (TP + FN)",
    color = "Assay concentration",
    title = "ROC curve for predicting an Kd below threshold"
  ) +
  facet_wrap(~ic50_threshold) +
  coord_cartesian(xlim = c(-0.1, 1), ylim = c(0, 1.1))

ggsave(
  "plots/percent_control_ic50_roc.pdf",
  single_dose_comp_roc_plot,
  width = 12, height = 4
)
```