kinase_families.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() %>%
  filter(dataset == "original_repeat_replaced")

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

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

single_dose <- syn("syn26486828") %>%
  fread()

partition_index_max <- syn("syn52576129") %>%
  read_csv() %>%
  filter(dataset == "original_repeat_replaced")

structural_similarity_clusters <- syn("syn63660132") %>%
  read_csv()

download.file("https://github.com/sorgerlab/famplex/raw/43732660e42fd1746c2d1e61552bda0a90cbd664/relations.csv", destfile = "~/data/famplex_relations.csv")
famplex_relations <- read_csv("~/data/famplex_relations.csv", col_names = c("namespace", "id", "relation", "family_namespace", "family"))
```

Extract the kinases from the famplex relations.

Use list of preselected family kinases

```{r}
famplex_kinases <- famplex_relations %>%
  filter(
    relation == "isa",
    id %in% kinomescan_kinase_info$hgnc_symbol,
    family %in% c(
      "AKT", "AMPK_alpha", "CDK", "CSNK1",
      "ERK", "FGFR", "GSK3", "PDGFR", "ERBB",
      "JAK", "P90RSK", "MAP3K", "TAOK", "MAP2K",
      "PAK", "ROCK", "p38", "ADRBK", "DDR", "EPHA",
      "EPHB", "GRK", "IKK_family", "JNK", "NTRK",
      "PKC", "SRC", "VEGFR", "PIM", "SIK", "DYRK",
      "RSTK1", "RSTK2", "DAPK", "CLK", "RIPK", "PKN", "LATS"
    )
  ) %>%
  distinct(hgnc_symbol = id, family) %>%
  # manually add HIP kinases
  bind_rows(
    tibble(
      hgnc_symbol = c("HIPK1", "HIPK2", "HIPK3", "HIPK4"),
      family = "HIP"
    )
  )
```

Clean pseudo IC50s by retaining only concordant values and dose-response curves
that are not discordant.

```{r}
kinomescan_pseudo_kds_clean <- kinomescan_pseudo_kds %>%
  filter(
    classification != "discordant"
  ) %>%
  # Remove any mutants
  inner_join(
    kinomescan_kinase_info %>%
      filter(`Kinase Form` == "Wild Type") %>%
      select(`DiscoveRx Gene Symbol`)
  ) %>%
  # Select kinase form (complex or PTM) with the strongest inhibition
  filter(
    `DiscoveRx Gene Symbol` %in% {
      group_by(., `DiscoveRx Gene Symbol`, hgnc_symbol) %>%
        summarize(
          frac_binding = sum(classification == "binding") / n(),
          .groups = "drop"
        ) %>%
        group_by(hgnc_symbol) %>%
        arrange(desc(frac_binding)) %>%
        slice_head(n = 1) %>%
        pull(`DiscoveRx Gene Symbol`)
    }
  ) %>%
  # group_by(lspci_id, name, hgnc_symbol, `Entrez Gene Symbol`) %>%
  # arrange(classification == "binding", pseudo_ic50) %>%
  # slice_head(n = 1) %>%
  # ungroup() %>%
  power_inner_join(
    famplex_kinases %>%
      select(hgnc_symbol, family),
    check = check_specs(
      duplicate_keys_right = "warn"
    )
  )
```

Compute number of compounds inhibiting each kinase

```{r}
kinase_stats <- kinomescan_pseudo_kds_clean %>%
  group_by(hgnc_symbol, family, `DiscoveRx Gene Symbol`) %>%
  summarize(
    n_binding = sum(classification == "binding"),
    .groups = "drop"
  ) %>%
  arrange(desc(n_binding))
```

Group kinases by family and compute some summary stats for each family member,
like the number of compounds inhibiting it

```{r}
kinase_stats_by_family <- kinase_stats %>%
  mutate(
    across(c(hgnc_symbol, `DiscoveRx Gene Symbol`), fct_inorder),
    family = factor(
      family,
      levels = group_by(., family) %>%
        summarize(max_spread = max(n_binding) - min(n_binding), .groups = "drop") %>%
        arrange(desc(max_spread)) %>%
        pull(family)
    )
  )

library(ggpattern)
library(ggrastr)
p <- kinase_stats_by_family %>%
  inner_join(
    kinomescan_kinase_info %>%
      select(`DiscoveRx Gene Symbol`, `Expression System`, `Kinase Construct`) %>%
      mutate(across(`DiscoveRx Gene Symbol`, factor, levels = kinase_stats$`DiscoveRx Gene Symbol`))
  ) %>%
  mutate(
    across(c(`Kinase Construct`, `Expression System`), factor),
  ) %>%
  ggplot(
    aes(
      x = `DiscoveRx Gene Symbol`,
      y = n_binding,
      fill = `Expression System`,
      pattern = `Kinase Construct`
    )
  ) +
    coord_flip() +
    # geom_col() +
    geom_col_pattern(color = "black") +
    scale_pattern_manual(values = c("Full Length" = "none", "Partial Length" = "circle")) +
    facet_wrap(~family, scales = "free") +
    theme_minimal() +
    theme(legend.key.size = unit(1, "cm")) +
    labs(x = "Kinase", y = "Number of compounds inhibiting kinase")

ggsave(
  "plots/kinase_families_n_binding.pdf", p,
  width = 16,
  height = 12,
)
```

```{r}
library(seriation)
co_cluster_mats <- function(mats) {
  set.seed(42)
  dms <- map(mats, dist, method = "euclidean")
  dm <- reduce(dms, `+`)
  cr <- hclust(dm, method = "average") %>%
    reorder(dm, method = "olo")
  cr
}

cluster_mat <- function(mat) {
  set.seed(42)
  dr <- dist(mat, method = "euclidean")
  cr <- hclust(dr, method = "average") %>%
    reorder(dr, method = "olo")
  cr
}

library(pheatmap)
plot_mat <- function(mat, cr = NULL, cc = NULL) {
  viridis_map <- viridis::plasma(100)
  pheatmap::pheatmap(
    mat,
    show_rownames = FALSE, show_colnames = TRUE,
    color = rev(viridis_map),
    cluster_rows = if (is.null(cr)) cluster_mat(mat) else cr,
    cluster_cols = if (is.null(cc)) cluster_mat(t(mat)) else cc,
    treeheight_row = 20, treeheight_col = 20,
    silent = TRUE
  )
}

generate_mat <- function(targets) {
  all_doses <- single_dose %>%
    filter(
      `DiscoveRx Gene Symbol` %in% targets
    ) %>%
    select(
      lspci_id, `DiscoveRx Gene Symbol`, `Compound Concentration (nM)`, `Percent Control`
    ) %>%
    arrange(lspci_id, `DiscoveRx Gene Symbol`, `Compound Concentration (nM)`) %>%
    group_by(`Compound Concentration (nM)`) %>%
    filter(n() > 10) %>%
    ungroup()
  measured_in_all <- all_doses %>%
    group_by(lspci_id) %>%
    summarize(
      n = n(),
      .groups = "drop"
    ) %>%
    filter(n == max(n)) %>%
    pull(lspci_id)
  all_doses_mats <- all_doses %>%
    mutate(across(`Compound Concentration (nM)`, ~factor(as.character(.x), levels = as.character(sort(unique(.x)))))) %>%
    mutate(
      across(`Compound Concentration (nM)`, fct_relabel, paste0, " nM")
    ) %>%
    filter(lspci_id %in% measured_in_all) %>%
    group_by(`Compound Concentration (nM)`) %>%
    filter(n() > 10) %>%
    summarize(
      mat = cur_data() %>%
        pivot_wider(names_from = `DiscoveRx Gene Symbol`, values_from = `Percent Control`) %>%
        column_to_rownames("lspci_id") %>%
        as.matrix() %>%
        list(),
      .groups = "drop"
    ) %>% {
      set_names(.$mat, .$`Compound Concentration (nM)`)
    }
  # all_doses_mats
  heatmaps <- map(all_doses_mats, plot_mat, cr = co_cluster_mats(all_doses_mats), cc = co_cluster_mats(map(all_doses_mats, t)))
  # patchwork::wrap_plots(map(heatmaps, "gtable"))
  # egg::ggarrange(
  #   plots = map(heatmaps, "gtable"),
  #   labels = names(heatmaps)
  # )
  cowplot::plot_grid(
    plotlist = map(heatmaps, "gtable"),
    labels = names(heatmaps)
  )
  # do.call(gridExtra::grid.arrange, heatmaps)
  # heatmaps
}
```


```{r}
heatmaps <- generate_mat(
  kinase_stats %>%
    filter(family == "HIP") %>%
    pull(`DiscoveRx Gene Symbol`)
)

dir.create("plots/kinase_families_heatmaps", showWarnings = FALSE)
walk(
  kinase_stats$family %>%
    unique(),
    ~ggsave(
      paste0("plots/kinase_families_heatmaps/", .x, ".pdf"),
      generate_mat(
        kinase_stats %>%
          filter(family == .x) %>%
          pull(`DiscoveRx Gene Symbol`)
      )
    )
)
```


## CDK upset

```{r}
cdk_whitelist <- c(
  "CDK4", "CDK6", "CDK7",
  # "PFTAIRE2", "PFTAIRE1",
  "CDK14", "CDK15",
  "CDK16", "CDK17",
  "CDK18", "CDK3",
  "CDK2", "CDK5",
  "CDK13", "CDK8",
  "CDK8", "CDK19"
)
setdiff(cdk_whitelist, kinomescan_kinase_info$hgnc_symbol)

kinomescan_pseudo_kds_cdk <- kinomescan_pseudo_kds %>%
  filter(
    classification != "discordant",
    hgnc_symbol %in% cdk_whitelist
  ) %>%
  mutate(
    symbol = if_else(
      str_detect(`DiscoveRx Gene Symbol`, "CDK"), `DiscoveRx Gene Symbol`, hgnc_symbol
    )
  )

cdk_upset_data <- tibble(
  kd_cutoff = c(100, 1000)
) %>%
  mutate(
    upset_data = map(
      kd_cutoff,
      \(x) kinomescan_pseudo_kds_cdk %>%
        mutate(binding = pseudo_kd <= x) %>%
        group_by(lspci_id) %>%
        filter(sum(binding) > 0) %>%
        ungroup() %>%
        select(hmsl_id, name, lspci_id, symbol, binding) %>%
        pivot_wider(names_from = symbol, values_from = binding) %>%
        mutate(across(-c(name, lspci_id, hmsl_id), \(y) replace_na(y, FALSE)))
    ),
    p = map(
      upset_data,
      \(x) ComplexUpset::upset(
          select(x, -name, -hmsl_id) %>%
            column_to_rownames("lspci_id"),
          colnames(x) %>%
            setdiff(c("name", "lspci_id", "hmsl_id")),
          min_size = 1,
          name = "Compound targets",
          height_ratio = 1.5,
          width_ratio = .25
        )
    )
  )

pwalk(
  cdk_upset_data,
  \(kd_cutoff, p, ...) {
    ggsave(
      paste0("plots/cdk_upset_", kd_cutoff, ".pdf"),
      p,
      width = 12,
      height = 8
    )
  }
)
```

Including only OKL compounds when computing the most inhibited kinases
but compute number of hit kinases for all compounds

```{r}
kinomescan_pseudo_kds_clean2 <- kinomescan_pseudo_kds %>%
  filter(
    hgnc_symbol != "",
    classification != "discordant",
    pseudo_kd_concordant != "discordant"
  )

CUTOFFS <- c(100, 1000)

most_inhibited_raw <- kinomescan_pseudo_kds_clean2 %>%
  filter(library == "OKL") %>%
  group_by(hgnc_symbol, `DiscoveRx Gene Symbol`) %>%
  summarize(
    n_inhibitors = map_int(
      set_names(CUTOFFS),
      ~ sum(pseudo_kd <= as.numeric(.x), na.rm = TRUE)
    ) %>%
      list(),
    .groups = "drop"
  )

most_inhibited <- most_inhibited_raw %>%
  unnest_longer(n_inhibitors, indices_to = "cutoff", values_to = "n_inhibitors") %>%
  mutate(across(cutoff, as.numeric)) %>%
  arrange(desc(n_inhibitors))

most_targets_raw <- kinomescan_pseudo_kds_clean2 %>%
  group_by(lspci_id, hmsl_id, library, name) %>%
  summarize(
    n_targets = map_int(
      set_names(CUTOFFS),
      ~ sum(pseudo_kd <= as.numeric(.x), na.rm = TRUE)
    ) %>%
      list(),
    .groups = "drop"
  )

most_targets <- most_targets_raw %>%
  unnest_longer(n_targets, indices_to = "cutoff", values_to = "n_targets") %>%
  mutate(across(cutoff, as.numeric)) %>%
  arrange(desc(n_targets))
```


```{r}
library(ComplexHeatmap)
library(circlize)
library(rlang)
library(seriation)
library(viridisLite)
library(impute)

cluster_fun_eucl <- function(mat, sample_in_col = TRUE) {
  if (!sample_in_col) {
    mat <- t(mat)
  }
  mat_imp <- impute.knn(
    mat, rng.seed = 42
  )[["data"]]
  if (!sample_in_col) {
    mat_imp <- t(mat_imp)
    mat <- t(mat)
  }
  # browser()
  dist_mat <- dist(mat_imp)
  # dist_mat <- as.dist(mat)
  clust <- hclust(dist_mat, method = "average")
  reorder(clust, dist_mat, method = "OLO")
}

createAnnotation <- function(df, colors = list(), which = "column") {
  # Load necessary libraries
  color_maps <- list()

  # Iterate over each column in the dataframe
  for (col_name in names(df)) {
    # Check if the column is numeric
    if (is.numeric(df[[col_name]])) {
      # Create a color mapping function for numeric columns
      if (is.null(colors[[col_name]]))
        colors[[col_name]] <- c("white", "red")
      if (is.function(colors[[col_name]])) {
        col_fun <- colors[[col_name]]
      } else if (is.character(colors[[col_name]])) {
        n <- length(colors[[col_name]])
        if (n == 1) {
          col_fun <- colorRamp2(
            seq(from = min(df[[col_name]]), to = max(df[[col_name]]), length.out = 100),
            viridis(100, option = colors[[col_name]])
          )
        } else {
          col_fun <- colorRamp2(
            seq(from = min(df[[col_name]]), to = max(df[[col_name]]), length.out = n),
            colors[[col_name]]
          )
        }
      } else
        stop("Dont know how to handle colors for column ", col_name)
      color_maps[[col_name]] <- col_fun
    } else {
      if (is.character(colors[[col_name]])) {
        color_maps[[col_name]] <- colors[[col_name]]
        next
      }
      col_fun <- colors[[col_name]] %||% rainbow
      # Create a named vector of colors for categorical columns
      values <- df[[col_name]]
      unique_values <- if (is.factor(values)) levels(values) else unique(df[[col_name]])
      named_colors <- set_names(col_fun(length(unique_values)), unique_values)
      color_maps[[col_name]] <- named_colors
    }
  }

  # Combine all annotations
  combined_annotations <- HeatmapAnnotation(df = df, col = color_maps, which = which)

  return(combined_annotations)
}

```

```{r}

cdk_hms <- tibble(
  kd_cutoff = c(100, 1000)
) %>%
  mutate(
    data = map(
      kd_cutoff,
      \(x) kinomescan_pseudo_kds_cdk %>%
        mutate(
          binding = pseudo_kd <= x
        ) %>%
        group_by(lspci_id) %>%
        filter(sum(binding) > 0) %>%
        ungroup()
    ),
    hm = map2(
      data, kd_cutoff,
      \(x, co) {
        m <- x %>%
          transmute(symbol, lspci_id, pseudo_kd = log10(pseudo_kd)) %>%
          pivot_wider(names_from = lspci_id, values_from = pseudo_kd) %>%
          column_to_rownames("symbol") %>%
          as.matrix()
        col_anno <- most_targets %>%
          filter(cutoff == co) %>%
          slice(
            match(colnames(m), as.character(lspci_id))
          ) %>%
          select(lspci_id, n_targets) %>%
          power_inner_join(
            x %>%
              group_by(lspci_id) %>%
              summarize(n_targets_cdk = sum(binding), .groups = "drop"),
            by = "lspci_id",
            check = check_specs(
              duplicate_keys_right = "warn",
              duplicate_keys_left = "warn",
              unmatched_keys_left = "warn"
            )
          ) %>%
          column_to_rownames("lspci_id")
        row_anno <- most_inhibited %>%
          filter(cutoff == co) %>%
          power_inner_join(
            kinomescan_pseudo_kds_cdk %>%
              select(symbol, `DiscoveRx Gene Symbol`) %>%
              distinct(),
            by = "DiscoveRx Gene Symbol",
            check = check_specs(
              duplicate_keys_right = "warn",
              duplicate_keys_left = "warn",
              unmatched_keys_right = "warn"
            )
          ) %>%
          slice(
            match(rownames(m), symbol)
          ) %>%
          select(symbol, n_inhibitors) %>%
          column_to_rownames("symbol")
        # browser()
        Heatmap(
          m,
          name = "Pseudo Kd",
          col = viridisLite::viridis(100, direction = -1),
          cluster_columns = cluster_fun_eucl,
          cluster_rows = cluster_fun_eucl,
          left_annotation = createAnnotation(
            row_anno,
            which = "row"
          ),
          top_annotation = createAnnotation(
            col_anno,
            which = "column",
            colors = list(n_targets_cdk = "magma")
          ),
          show_column_names = FALSE,
          # cell_fun = \(j, i, x, y, width, height, fill) {
          #   # browser()
          #   if (!is.na(mat[i, j]) && mat[i, j] < log10(cu))
          #     grid.text("*", x, y, gp = gpar(fontsize = 10, col = better_contrast(fill, c("white", "black"))))
          # }
        )
      }
    )
  )

pwalk(
  cdk_hms,
  function(hm, kd_cutoff, ...) {
    withr::with_pdf(
      file.path("plots", paste0("cdk_heatmap_", kd_cutoff, ".pdf")),
      print(hm), width = 8, height = 5
    )
  }
)

```

```{r}
kd_log_mat <- kinomescan_pseudo_kds_clean %>%
  filter(lspci_id %in% structural_similarity_clusters$lspci_id) %>%
  select(lspci_id, `DiscoveRx Gene Symbol`, pseudo_kd) %>%
  pivot_wider(names_from = `DiscoveRx Gene Symbol`, values_from = pseudo_kd) %>%
  column_to_rownames("lspci_id") %>%
  as.matrix() %>%
  log10()

row_anno_data <- partition_index_max %>%
  distinct(
    lspci_id,
    name,
    partition_index
  ) %>%
  power_left_join(
    most_targets %>%
      filter(cutoff == 100) %>%
      select(lspci_id, n_targets),
    by = "lspci_id",
    check = check_specs(
      duplicate_keys_right = "warn",
      duplicate_keys_left = "warn",
      unmatched_keys_left = "warn"
    )
  ) %>%
  power_inner_join(
    structural_similarity_clusters %>%
      transmute(
        lspci_id,
        cluster = factor(cluster)
      ),
    by = "lspci_id",
    check = check_specs(
      duplicate_keys_right = "warn",
      duplicate_keys_left = "warn",
      unmatched_keys_right = "warn"
    )
  ) %>%
  slice(match(rownames(kd_log_mat), lspci_id))

row_anno <- createAnnotation(
  row_anno_data %>%
    select(-name) %>%
    column_to_rownames("lspci_id"),
  which = "row",
  colors = list(
    cluster = \(x) ggokabeito::palette_okabe_ito(seq_len(x))
    # cluster = \(x) RColorBrewer::brewer.pal(x, "Set2")
  )
)

col_anno_data <- most_inhibited %>%
  filter(cutoff == 100) %>%
  select(`DiscoveRx Gene Symbol`, n_inhibitors) %>%
  power_inner_join(
    kinomescan_kinase_info %>%
      select(`DiscoveRx Gene Symbol`, Group),
    by = "DiscoveRx Gene Symbol",
    check = check_specs(
      unmatched_keys_left = "warn",
      duplicate_keys_right = "warn",
      duplicate_keys_left = "warn",
    )
  ) %>%
  slice(match(colnames(kd_log_mat), `DiscoveRx Gene Symbol`))

col_anno <- createAnnotation(
  col_anno_data %>%
    column_to_rownames("DiscoveRx Gene Symbol"),
  which = "column"
)

hm <- Heatmap(
  kd_log_mat,
  name = "log10 Pseudo Kd",
  col = viridisLite::viridis(100, direction = -1),
  cluster_columns = \(x) cluster_fun_eucl(x, sample_in_col = FALSE),
  cluster_rows = \(x) cluster_fun_eucl(x, sample_in_col = FALSE),
  left_annotation = row_anno,
  top_annotation = col_anno,
  show_column_names = FALSE,
  show_row_names = FALSE
)
```