okl_selectivity.Rmd

---
title: "Analysis of OKL selectivity 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_datasets = "syn52504516",
  kinase_info = "syn51286743",
  klaeger_ic50s = "syn51288931",
  nominal_targets = "syn52947779",
  okl_pseudo_kds = "syn51080578",
  okl_classification = "syn51080579"
)

input_files <- map(
  inputs, syn
)

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

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

klaeger_ic50s <- input_files[["klaeger_ic50s"]] %>%
  fread() %>%
  as_tibble()

nominal_targets <- input_files[["nominal_targets"]] %>%
  fread() %>%
  as_tibble() %>%
  power_inner_join(
    single_dose_datasets %>%
      distinct(lspci_id, hmsl_id),
    by = c("HMSLID" = "hmsl_id"),
    check = check_specs(
      unmatched_keys_left = "warn",
      duplicate_keys_right = "warn"
    )
  )

pseudo_kd <- input_files[["okl_pseudo_kds"]] %>%
  fread() %>%
  as_tibble()

classification <- input_files[["okl_classification"]] %>%
  fread() %>%
  as_tibble()
```

For computing CATDS try to reduce noise by censoring Percent Control values of
compound-target pairs that are classified as non-binding. We frequently have
values of 90% activity remaining for non-binding compounds, which is likely
just noise but affects CATDS calculation significantly.

Setting Percent Control to 100% for measurements classified as non-binding
AND where Percent Control >35%

```{r}
catds <- single_dose_datasets %>%
  filter(
    `Compound Concentration (nM)` %in% c(12.5, 100, 1000, 10000),
    !exclude_target,
    !qpcr_miss
  ) %>%
  mutate(
    non_binding_censored = FALSE
  ) %>%
  bind_rows(
    power_left_join(
      ., select(classification, dataset, lspci_id, `DiscoveRx Gene Symbol`, classification),
      by = c("dataset", "lspci_id", "DiscoveRx Gene Symbol"),
      check = check_specs(
        unmatched_keys_left = "warn",
        duplicate_keys_right = "warn"
      )
    ) %>%
    mutate(
      non_binding_censored = TRUE,
      `Percent Control` = if_else(
        classification == "non-binding" & `Percent Control` > 35,
        100, `Percent Control`
      )
    )
  ) %>%
  group_by(non_binding_censored, dataset, lspci_id, `Compound Concentration (nM)`) %>%
  filter(n() > 50) %>%
  mutate(
    catds = (100 - `Percent Control`) / sum(100 - `Percent Control`)
  ) %>%
  ungroup() %>%
  select(dataset, non_binding_censored, everything())

catds_max <- catds %>%
  group_by(dataset, non_binding_censored, lspci_id, `Compound Concentration (nM)`) %>%
  slice_max(
    catds,
    n = 1,
    na_rm = TRUE,
    with_ties = TRUE
  ) %>%
  ungroup()

catds %>%
  group_nest(non_binding_censored) %>%
  pwalk(
    \(non_binding_censored, data) {
      write_csv(
        data %>%
          select(where(negate(is.list))),
        paste0("classification/okl_eurofins_catds", if (non_binding_censored) "_censored" else "", ".csv")
      )
    }
  )

catds_max %>%
  group_nest(non_binding_censored) %>%
  pwalk(
    \(non_binding_censored, data) {
      write_csv(
        data %>%
          select(where(negate(is.list))),
        paste0("classification/okl_eurofins_catds_index_max", if (non_binding_censored) "_censored" else "", ".csv")
      )
    }
  )

```

CATDSmax here can include ties. This happens especially when a compound has
multiple targets with the same potency at the limit of our resolution (e.g. <12.5nM)


```{r}
catds %>%
  filter(
    `Compound Concentration (nM)` %in% c(12.5, 100, 1000, 10000)
  ) %>%
  ggplot(
    aes(catds, color = fct_inseq(as.character(`Compound Concentration (nM)`)))
  ) +
  geom_density() +
  facet_wrap(~dataset) +
  scale_y_log10()

catds_max %>%
  filter(
    `Compound Concentration (nM)` %in% c(12.5, 100, 1000, 10000)
  ) %>%
  ggplot(
    aes(catds, color = fct_inseq(as.character(`Compound Concentration (nM)`)))
  ) +
  geom_density() +
  facet_wrap(~dataset + non_binding_censored)

catds_max %>%
  ggplot(
    aes(catds, color = fct_inseq(as.character(`Compound Concentration (nM)`)))
  ) +
  stat_ecdf() +
  facet_wrap(~dataset + non_binding_censored) +
  lims(x = c(0, 1))

catds %>%
  filter(name == "TOFACITINIB", str_detect(hgnc_symbol, "JAK")) %>%
  View()

catds %>%
  mutate(target_engagement = 100 - `Percent Control`) %>%
  filter(name == "TOFACITINIB", `Compound Concentration (nM)` == 12.5) %>%
  View()
```


Compute partition index

Analysis of kinase inhibitor selectivity using a thermodynamics-based partition index
Alan C Cheng  1 , John Eksterowicz, Stephanie Geuns-Meyer, Yaxiong Sun


```{r}
partition_index <- pseudo_kd %>%
  # filter(pseudo_ic50_relation == "=") %>%
  filter(
    !exclude_target,
    classification != "discordant"
  ) %>%
  group_by(dataset, lspci_id) %>%
  filter(n() > 50) %>%
  mutate(
    partition_index = (1 / pseudo_kd) / sum(1 / pseudo_kd, na.rm = TRUE)
  ) %>%
  ungroup()

partition_index_max <- partition_index %>%
  group_by(
    dataset, lspci_id, hmsl_id, library, name
  ) %>%
  slice_max(
    partition_index,
    n = 1,
    na_rm = TRUE,
    with_ties = TRUE
  )

write_csv(
  partition_index %>%
    select(where(negate(is.list))),
  "classification/okl_eurofins_partition_index.csv"
)

write_csv(
  partition_index_max %>%
    select(where(negate(is.list))),
  "classification/okl_eurofins_partition_index_max.csv"
)
```

Same as CATDSmax, partition_index_max can include ties.

```{r}
partition_index_max %>%
  ggplot(aes(partition_index_max)) +
    geom_histogram() +
    facet_wrap(~dataset + library)

partition_index_max %>%
  ggplot(aes(partition_index_max, color = dataset)) +
    stat_ecdf() +
    lims(x = c(0, 1)) +
    facet_wrap(~library)
```


```{r}
klaeger_partition_index <- klaeger_ic50s %>%
  # filter(IC50_relation == "=") %>%
  group_by(lspci_id, Compound) %>%
  mutate(
    partition_index = (1 / IC50) / sum(1 / IC50, na.rm = TRUE)
  ) %>%
  ungroup()

klaeger_partition_index_max <- klaeger_partition_index %>%
  group_by(lspci_id, Compound) %>%
  summarize(
    partition_index_max = max(partition_index, na.rm = TRUE),
    .groups = "drop"
  )
```



```{r}
synStoreMany(
  c(
    "classification/okl_eurofins_partition_index.csv",
    "classification/okl_eurofins_partition_index_max.csv",
    "classification/okl_eurofins_catds.csv",
    "classification/okl_eurofins_catds_index_max.csv",
    "classification/okl_eurofins_catds_censored.csv",
    "classification/okl_eurofins_catds_index_max_censored.csv"
  ),
  parentId = "syn18508401",
  used = unname(inputs),
  executed = "https://github.com/labsyspharm/okl-analysis/blob/main/okl_selectivity.Rmd",
  forceVersion = FALSE
)

```

```{r}
klaeger_partition_index_max %>%
  ggplot(aes(partition_index_max)) +
    stat_ecdf() +
    lims(x = c(0, 1))

shared_targets_klaeger_okl <- intersect(
  klaeger_ic50s$Name,
  pseudo_ic50$hgnc_symbol
)

kd_klaeger_okl <- bind_rows(
  pseudo_kd %>%
    filter(
      hgnc_symbol %in% shared_targets_klaeger_okl,
      !exclude_target,
      classification != "discordant"
    ) %>%
    transmute(
      lspci_id, hgnc_symbol, name, dataset,
      pseudo_ic50, pseudo_ic50_relation, pseudo_ic50_concordant,
      source = "okl"
    ),
  klaeger_ic50s %>%
    filter(
      Name %in% shared_targets_klaeger_okl
    ) %>%
    transmute(
      lspci_id, hgnc_symbol = Name, name = Compound,
      dataset = "klaeger",
      pseudo_ic50 = IC50, pseudo_ic50_relation = IC50_relation,
      pseudo_ic50_concordant = "concordant",
      source = "klaeger"
    )
) %>%
  mutate(
    source_dataset = paste(source, dataset, sep = "_")
  )

partition_index_klaeger_okl <- ic50_klaeger_okl %>%
  # filter(pseudo_ic50_relation == "=") %>%
  group_by(source, dataset, source_dataset, lspci_id) %>%
  # filter(n() > 50) %>%
  mutate(
    partition_index = (1 / pseudo_ic50) / sum(1 / pseudo_ic50, na.rm = TRUE)
  ) %>%
  ungroup()

partition_index_max_klaeger_okl <- partition_index_klaeger_okl %>%
  group_by(source, dataset, source_dataset, lspci_id) %>%
  summarize(
    partition_index_max = max(partition_index, na.rm = TRUE),
    .groups = "drop"
  )

partition_index_max_klaeger_okl %>%
  ggplot(aes(partition_index_max, color = source_dataset)) +
    stat_ecdf() +
    lims(x = c(0, 1))

p <- partition_index_max_klaeger_okl %>%
  filter(
    lspci_id %in% {
      group_by(., lspci_id) %>%
      filter(uniqueN(source) == 2) %>%
      pull(lspci_id)
    }
  ) %>%
  ggplot(aes(partition_index_max, color = source_dataset)) +
    stat_ecdf() +
    lims(x = c(0, 1)) +
    labs(
      x = "Partition index",
      y = "Cumulative fraction of compounds",
      title = "Compounds and targets shared between Klaeger and OKL",
      color = "Dataset"
    )

ggsave(
  file.path("plots", "partition_index_klaeger_okl_shared_compounds_and_targets_ecdf.pdf"),
  p, width = 6, height = 4
)

p <- partition_index_max_klaeger_okl %>%
  filter(
    lspci_id %in% {
      group_by(., lspci_id) %>%
      filter(uniqueN(source) == 2) %>%
      pull(lspci_id)
    }
  ) %>%
  select(-source, -dataset) %>%
  pivot_wider(names_from = source_dataset, values_from = partition_index_max) %>%
  ggplot(aes(klaeger_klaeger, okl_original_repeat_replaced)) +
    geom_point(shape = 16, alpha = 0.5) +
    geom_smooth(method = "lm") +
    lims(x = c(0, 1), y = c(0, 1)) +
    labs(x = "Partition index Klaeger", y = "Partition index OKL (repeat)")


ggsave(
  file.path("plots", "partition_index_klaeger_vs_okl_shared_compounds_and_targets.pdf"),
  p, width = 5, height = 5
)

partition_index_max_klaeger_okl %>%
  filter(
    lspci_id %in% {
      group_by(., lspci_id) %>%
      filter(uniqueN(source) == 2) %>%
      pull(lspci_id)
    }
  ) %>%
  select(-source, -dataset) %>%
  pivot_wider(names_from = source_dataset, values_from = partition_index_max) %>%
  {cor.test(.$klaeger_klaeger, .$okl_original_repeat_replaced)}


partition_index_max_klaeger_okl %>%
  filter(
    lspci_id %in% {
      group_by(., lspci_id) %>%
      filter(uniqueN(source) == 2) %>%
      pull(lspci_id)
    }
  ) %>%
  ggplot(aes(partition_index_max, fill = source_dataset)) +
    geom_bar(stat = "bin", color = "black", alpha = 0.5, position = "identity") +
    lims(x = c(0, 1))


```