gene_set_analysis_plotting.Rmd

---
title: "Peripheral neuropathy gene set analysis plotting"
author: "Clemens Hug"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

library(here)
library(tidyverse)
library(data.table)
library(synExtra)
library(qs)
library(powerjoin)
library(plotly)
library(ggrepel)

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

## Reading GSE results

```{r}
topgo_res_pert <- syn("syn64461453") %>%
  qread()

fgsea_res_pert <- syn("syn64461445") %>%
  qread()

fgsea_res_rescue <- syn("syn64461444") %>%
  qread()

de_long_selected <- syn("syn64582734") %>%
  read_csv()

feasible_genes <- syn("syn64582742") %>%
  read_csv()

fgsea_overrep <- syn("syn64582959") %>%
  qread()

fgsea_intersection_overrep <- syn("syn64586672") %>%
  qread()

rag_ion_channel_genes <- syn("syn64745787") %>%
  readxl::read_excel() %>%
  pivot_longer(
    everything(),
    names_to = "gene_set",
    values_to = "gene_symbol_mouse"
  ) %>%
  drop_na()

ptx_rescue_list <- syn("syn64608137") %>%
  read_csv()

ptx_rescue_types <- syn("syn64608138") %>%
  read_csv()
```



## Plot Reactome PTX only and comparison with combination treatment


```{r}
fgsea_res_single_treatments <- fgsea_res_pert %>%
  filter(
    metric == "signed_padj"
  ) %>%
  separate_wider_delim(
    pathway,
    delim = "_",
    names = c("database", "pathway"),
    too_many = "merge"
  ) %>%
  filter(
    database == "REACTOME",
    padj < 0.01
  ) %>%
  mutate(
    sign = if_else(NES > 0, "pos", "neg")
  )

ps <- fgsea_res_single_treatments %>%
  group_nest(
    contrast
  ) %>%
  mutate(
    p = map2(
      contrast, data,
      \(x, y) {
        y %>%
          mutate(
            pathway = fct_reorder(pathway, NES)
          ) %>%
          ggplot(
            aes(
              x = NES,
              y = pathway,
              fill = sign
            )
          ) +
          geom_col() +
          geom_text(
            aes(
              x = if_else(sign == "pos", -.05, .05),
              hjust = if_else(sign == "pos", 1, 0),
              label = pathway
            ),
            size = 2
          ) +
          scale_fill_discrete(guide = "none", direction = -1) +
          theme(
            axis.text.y = element_blank(),
            panel.grid.major.y = element_blank(),
          ) +
          labs(
            x = "Normalized enrichment score",
            y = NULL,
            title = paste0("Reactome pathways enriched in ", x, " treatment"),
            subtitle = "FDR < 0.01"
          ) + {
            if (!any(y$sign == "neg"))
              lims(x = c(-max(abs(y$NES)), NA))
          }
      }
    )
  )

pwalk(
  ps,
  \(contrast, p, ...) {
    ggsave(
      file.path("plots", paste0("reactome_", contrast, "_bars.pdf")),
      p, width = 8, height = 8
    )
  }
)
```

### Comparison with combination treatment

```{r}
all_contrasts <- fgsea_res_pert$contrast %>%
  unique()

contrast_plotting_pairs <- bind_rows(
  tibble(
    contrast_1 = "PTX",
    contrast_2 = all_contrasts[
      str_ends(all_contrasts, fixed("+ PTX")) |
      str_ends(all_contrasts, fixed("+ PTX vs PTX"))
    ]
  ),
  tibble(
    contrast_1 = "GNE-495 + PTX",
    contrast_2 = "JNK-inh + PTX"
  )
)

fgsea_res_ptx_vs_combo <- fgsea_res_pert %>%
  filter(
    metric == "signed_padj"
  ) %>%
  separate_wider_delim(
    pathway,
    delim = "_",
    names = c("database", "pathway"),
    too_many = "merge"
  ) %>%
  mutate(
    signed_padj = -sign(NES) * log10(padj)
  ) %>%
  filter(
    database == "REACTOME"
  ) %>% {
    d <- .
    power_inner_join(
      contrast_plotting_pairs,
      d,
      by = join_by(contrast_1 == contrast),
      suffix = c("", "_1"),
      check = check_specs(
        unmatched_keys_left = "warn"
      )
    ) %>%
    power_inner_join(
      d,
      by = join_by(contrast_2 == contrast, metric, database, pathway),
      suffix = c("_1", "_2"),
      check = check_specs(
        unmatched_keys_left = "warn"
      )
    )
  } %>%
  mutate(
    sig_class = case_when(
      padj_1 < .05 & padj_2 < .05 ~ "both",
      padj_1 < .05 ~ "1",
      padj_2 < .05 ~ "2",
      TRUE ~ "none"
    )
  )

sig_class_colors <- c(
  "both" = "purple",
  "1" = "red",
  "2" = "blue",
  "none" = "grey"
)

ps <- fgsea_res_ptx_vs_combo %>%
  group_nest(contrast_1, contrast_2) %>%
  mutate(
    p = pmap(
      .,
      \(contrast_1, contrast_2, data, ...) {
        ggplot(
          data %>%
            arrange(match(sig_class, names(sig_class_colors))),
          aes(
            x = signed_padj_1,
            y = signed_padj_2
          )
        ) +
        geom_point(
          aes(
            color = sig_class,
            text = pathway
          ),
          alpha = .8
        ) +
        scale_color_manual(
          name = "Significant\nenrichment\nFDR < 0.05",
          values = sig_class_colors,
          breaks = names(sig_class_colors),
          labels = c(
            "Both",
            contrast_1,
            contrast_2,
            "None"
          )
        ) +
        theme_minimal() +
        labs(
          x = paste("Enrichment", contrast_1),
          y = paste("Enrichment", contrast_2)
        )
      }
    )
  )

pwalk(
  ps,
  \(contrast_1, contrast_2, p, ...) {
    ggsave(
      file.path("plots", paste0("reactome_", contrast_1, "_vs_", contrast_2, "_scatter_padj.pdf")),
      p, width = 7, height = 6
    )
    htmlwidgets::saveWidget(
      ggplotly(p),
      file.path("plots", paste0("reactome_", contrast_1, "_vs_", contrast_2, "_scatter_padj.html"))
    )
  }
)


ps <- fgsea_res_ptx_vs_combo %>%
  group_nest(contrast_ptx, contrast_combo) %>%
  mutate(
    p = map2(
      contrast_combo, data,
      \(x, y) {
        ggplot(
          y %>%
            arrange(match(sig_class, names(sig_class_colors))),
          aes(
            x = NES_ptx,
            y = NES_combo
          )
        ) +
        geom_point(
          aes(
            color = sig_class,
            text = pathway
          ),
          alpha = .8
        ) +
        scale_color_manual(
          values = sig_class_colors
        )
      }
    )
  )

pwalk(
  ps,
  \(contrast_ptx, contrast_combo, p, ...) {
    ggsave(
      file.path("plots", paste0("reactome_", contrast_ptx, "_vs_", contrast_combo, "_scatter.pdf")),
      p, width = 6, height = 6
    )
    htmlwidgets::saveWidget(
      ggplotly(p),
      file.path("plots", paste0("reactome_", contrast_ptx, "_vs_", contrast_combo, "_scatter.html"))
    )
  }
)
```

### Highlighting certain pathways

```{r}
highlighted_pathways <- tribble(
  ~contrast_1, ~contrast_2, ~pathways,
  "GNE-495 + PTX", "JNK-inh + PTX", c(
    "METABOLISM_OF_LIPIDS",
    "CHOLESTEROL_BIOSYNTHESIS",
    "ACTIVATION_OF_GENE_EXPRESSION_BY_SREBF_SREBP",
    "METABOLISM_OF_STEROIDS",
    "REGULATION_OF_CHOLESTEROL_BIOSYNTHESIS_BY_SREBP_SREBF",
    "FATTY_ACID_METABOLISM",
    "FATTY_ACYL_COA_BIOSYNTHESIS",
    "PHOSPHOLIPID_METABOLISM",
    "NEUTROPHIL_DEGRANULATION",
    "VESICLE_MEDIATED_TRANSPORT",
    "ER_TO_GOLGI_ANTEROGRADE_TRANSPORT",
    "CELL_CYCLE_MITOTIC",
    "RHO_GTPASE_EFFECTORS",
    "CYTOKINE_SIGNALING_IN_IMMUNE_SYSTEM",
    "EXTRACELLULAR_MATRIX_ORGANIZATION",
    "COPI_MEDIATED_ANTEROGRADE_TRANSPORT",
    "SYNDECAN_INTERACTIONS",
    ##
    "ACTIVATION_OF_GENE_EXPRESSION_BY_SREBF_SREBP",
    "FATTY_ACID_METABOLISM",
    "METABOLISM_OF_LIPIDS",
    "TRNA_AMINOACYLATION",
    "CELL_CYCLE",
    "REGULATION_OF_EXPRESSION_OF_SLITS_AND_ROBOS",
    "GAP_JUNCTION_TRAFFICKING_AND_REGULATION",
    "ER_TO_GOLGI_ANTEROGRADE_TRANSPORT",
    "COPI_MEDIATED_ANTEROGRADE_TRANSPORT"
  ) %>%
    unique()
)


ps <- fgsea_res_ptx_vs_combo %>%
  group_nest(contrast_1, contrast_2) %>%
  inner_join(
    highlighted_pathways,
    by = c("contrast_1", "contrast_2")
  ) %>%
  mutate(
    p = pmap(
      .,
      \(contrast_1, contrast_2, data, pathways, ...) {
        ggplot(
          data %>%
            arrange(match(sig_class, names(sig_class_colors))),
          aes(
            x = signed_padj_1,
            y = signed_padj_2
          )
        ) +
        geom_point(
          aes(
            color = sig_class,
            text = pathway
          ),
          alpha = .8
        ) +
        geom_text_repel(
          data = \(x) mutate(
              x,
              pathway_short = if_else(
                pathway %in% pathways,
                pathway,
                ""
              ) %>%
                str_replace_all(fixed("_"), " ") %>%
                str_wrap(20)
            ),
          aes(
            label = pathway_short
          ),
          size = 2,
          # force = 3,
          max.iter = 1e6,
          seed = 42,
          max.overlaps = Inf,
          min.segment.length = 0.3
        ) +
        scale_color_manual(
          name = "Significant\nenrichment\nFDR < 0.05",
          values = sig_class_colors,
          breaks = names(sig_class_colors),
          labels = c(
            "Both",
            contrast_1,
            contrast_2,
            "None"
          )
        ) +
        scale_x_continuous(
          expand = expansion(mult = c(.18, .05))
        ) +
        scale_y_continuous(
          expand = expansion(mult = c(.1, .05))
        ) +
        theme_minimal() +
        labs(
          x = paste("Enrichment", contrast_1),
          y = paste("Enrichment", contrast_2)
        )
      }
    )
  )

pwalk(
  ps,
  \(contrast_1, contrast_2, p, ...) {
    ggsave(
      file.path("plots", paste0("reactome_", contrast_1, "_vs_", contrast_2, "_scatter_padj_highlighted.pdf")),
      p, width = 7.5, height = 6
    )
    # htmlwidgets::saveWidget(
    #   ggplotly(p),
    #   file.path("plots", paste0("reactome_", contrast_ptx, "_vs_", contrast_combo, "_scatter_padj.html"))
    # )
  }
)
```


## Heatmap of Reactome PTX and combo


```{r}

library(seriation)
cluster_df <- function(df, row_var, col_var, value_var, values_fill = 0) {
  # browser()
  mat <- df %>%
    select({{row_var}}, {{col_var}}, {{value_var}}) %>%
    pivot_wider(names_from = {{col_var}}, values_from = {{value_var}}, values_fill = values_fill) %>%
    column_to_rownames(rlang::as_name(rlang::enquo(row_var)))
  # browser()
  if (rlang::is_bare_numeric(pull(df, {{value_var}}))) {
    dist_rows <- dist(mat, method = "euclidian")
    dist_cols <- dist(t(mat), method = "euclidian")
  } else {
    # browser()
    dist_rows <- cluster::daisy(mat, metric = "gower")
    dist_cols <- t(mat) %>%
      as.data.frame() %>%
      mutate(across(everything(), \(x) factor(x, levels = levels(pull(df, {{value_var}}))))) %>%
      cluster::daisy(metric = "gower")
  }
  clust_rows <- hclust(dist_rows, method = "average") %>%
      reorder(dist_rows, method = "olo")
  clust_cols <- hclust(dist_cols, method = "average") %>%
      reorder(dist_cols, method = "olo")
  df %>%
    mutate(
      "{{row_var}}" := factor({{row_var}}, levels = clust_rows$labels[clust_rows$order]),
      "{{col_var}}" := factor({{col_var}}, levels = clust_cols$labels[clust_cols$order])
    )
}


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)
  # dist_mat <- as.dist(mat)
  clust <- hclust(dist_mat, method = "average")
  reorder(clust, dist_mat, method = "OLO")
}

```

```{r}
fgsea_hm_conditions <- c(
  "PTX", "GNE-495", "JNK-inh",
  "GNE-495 + PTX", "JNK-inh + PTX"
)

fgsea_res_pert_hm_data <- fgsea_res_pert %>%
  filter(
    metric == "signed_padj",
    contrast %in% fgsea_hm_conditions
  ) %>%
  separate_wider_delim(
    pathway,
    delim = "_",
    names = c("database", "pathway"),
    too_many = "merge"
  ) %>%
  mutate(
    signed_padj = -sign(NES) * log10(padj)
  ) %>%
  filter(
    database == "REACTOME"
  ) %>%
  cluster_df(
    pathway, contrast, signed_padj, values_fill = 0
  )

fgsea_res_pert_mat <- fgsea_res_pert_hm_data %>%
  arrange(pathway, contrast) %>%
  select(pathway, contrast, signed_padj) %>%
  pivot_wider(names_from = contrast, values_from = signed_padj, values_fill = 0) %>%
  column_to_rownames("pathway") %>%
  as.matrix()

fgsea_res_pert_hm_pathways <- fgsea_res_pert_hm_data %>%
  filter(
    padj < 0.01
  ) %>%
  pull(pathway) %>%
  unique()
```


```{r}
row_clust <- cluster_fun_eucl(
  fgsea_res_pert_mat[as.character(fgsea_res_pert_hm_pathways), fgsea_hm_conditions]
)

row_clust_df_long <- tibble(k = 4:10) %>%
  mutate(
    res = map(k, \(x) enframe(cutree(row_clust, k = x), name = "pathway", value = "cluster"))
  ) %>%
  unnest(res)

row_clust_df <- row_clust_df_long %>%
  mutate(
    cluster = fct_inseq(as.character(cluster))
  ) %>%
  pivot_wider(
    names_from = k,
    values_from = cluster,
    names_prefix = "k_"
  )
```

k = 8 seems to work well


```{r}
dir.create(here("results"), showWarnings = FALSE)
write_csv(
  row_clust_df_long %>%
    filter(k == 8),
  here("results", "reactome_pert_clusters_k8.csv")
)

```

```{r}
library(ComplexHeatmap)
hm <- Heatmap(
  fgsea_res_pert_mat[as.character(fgsea_res_pert_hm_pathways), fgsea_hm_conditions],
  cluster_rows = cluster_fun_eucl,
  cluster_columns = FALSE,
  show_row_names = FALSE,
  left_annotation = HeatmapAnnotation(
    df = column_to_rownames(row_clust_df, "pathway")[as.character(fgsea_res_pert_hm_pathways), ],
    which = "row"
  )
)

mat <- fgsea_res_pert_mat[as.character(fgsea_res_pert_hm_pathways), fgsea_hm_conditions]
mat_max_abs <- max(abs(quantile(mat, c(.025, .975), na.rm = TRUE)))
hm <- Heatmap(
  mat,
  col = circlize::colorRamp2(seq(from = -mat_max_abs, to = mat_max_abs, length.out = 51), paletteer::paletteer_c("pals::ocean.balance", n = 51)),
  name = "Normalized enrichment score",
  cluster_rows = cluster_fun_eucl,
  cluster_columns = FALSE,
  show_row_names = FALSE,
  left_annotation = HeatmapAnnotation(
    df = row_clust_df %>%
        select(pathway, cluster = k_8) %>%
        slice(match(as.character(fgsea_res_pert_hm_pathways), pathway)) %>%
        column_to_rownames("pathway"),
    col = list(
      cluster = set_names(
        ggokabeito::palette_okabe_ito(1:8), 1:8
      )
    ),
    which = "row"
  ),
  right_annotation = HeatmapAnnotation(
    pathways = c(
      "ACTIVATION_OF_GENE_EXPRESSION_BY_SREBF_SREBP",
      "FATTY_ACID_METABOLISM",
      "METABOLISM_OF_LIPIDS",
      "TRNA_AMINOACYLATION",
      "CELL_CYCLE",
      "REGULATION_OF_EXPRESSION_OF_SLITS_AND_ROBOS",
      "GAP_JUNCTION_TRAFFICKING_AND_REGULATION",
      "ER_TO_GOLGI_ANTEROGRADE_TRANSPORT",
      "COPI_MEDIATED_ANTEROGRADE_TRANSPORT"
    ) %>%
      {
        anno_mark(
          at = match(., fgsea_res_pert_hm_pathways),
          labels = str_wrap(
            str_replace_all(., "_", " "),
            20
          ),
          labels_gp = gpar(fontsize = 6),
          padding = unit(2, "mm")
        )
      },
    which = "row"
  )
)

withr::with_pdf(
  here("plots", "reactome_pert_heatmap_k8.pdf"),
  draw(hm),
  width = 6, height = 8
)

InteractiveComplexHeatmap::htShiny(hm, save = here("plots/reactome_pert_heatmap_k8"))
```


## Plotting rescue vector FGSEA results

```{r}
fgsea_res_rescue_data <- fgsea_res_rescue %>%
  filter(metric == "signed_padj_rescue_clamped") %>%
  separate_wider_delim(
    pathway,
    delim = "_",
    names = c("database", "pathway"),
    too_many = "merge"
  ) %>%
  mutate(
    sign = if_else(NES > 0, "pos", "neg"),
    pathway_db = paste(str_sub(database, 1, 1), pathway, sep = "_")
  )


ps <- fgsea_res_rescue_data %>%
  group_nest(
    compound
  ) %>%
  mutate(
    p = map2(
      compound, data,
      \(x, y) {
        y %>%
          filter(
            padj < .05
          ) %>%
          mutate(
            across(pathway_db, \(x) fct_reorder(x, NES))
          ) %>%
          ggplot(
            aes(
              x = NES,
              y = pathway_db,
              fill = sign
            )
          ) +
          geom_col() +
          geom_text(
            aes(
              x = if_else(sign == "pos", -.05, .05),
              hjust = if_else(sign == "pos", 1, 0),
              label = pathway_db
            ),
            size = 2
          ) +
          scale_fill_discrete(guide = "none", direction = -1) +
          theme(
            axis.text.y = element_blank(),
            panel.grid.major.y = element_blank(),
          ) +
          labs(
            x = "Normalized enrichment score",
            y = NULL,
            title = paste0("Pathways enriched in genes rescued by ", x),
            subtitle = "FDR < 0.05"
          ) + {
            if (!any(y$sign == "neg"))
              lims(x = c(-max(abs(y$NES)), NA))
          }
      }
    )
  )

pwalk(
  ps,
  \(compound, p, ...) {
    ggsave(
      file.path("plots", paste0("reactome_", compound, "_bars_rescue.pdf")),
      p, width = 8, height = 8
    )
  }
)
```


### Comparing rescue vector enrichment between combos

```{r}
fgsea_rescue_combo_vs_combo <- fgsea_res_rescue_data %>%
  mutate(
    signed_padj = -sign(NES) * log10(padj)
  ) %>%
  # filter(
  #   database == "REACTOME"
  # ) %>%
  pivot_wider(
    names_from = compound,
    values_from = -c(metric, compound, pathway, database)
  ) %>%
  mutate(
    sig_class = case_when(
      `padj_GNE-495` < .05 & `padj_JNK-inh` < .05 ~ "both",
      `padj_GNE-495` < .05 ~ "GNE-495",
      `padj_JNK-inh` < .05 ~ "JNK-inh",
      TRUE ~ "none"
    )
  )

sig_class_colors <- c(
  "both" = "purple",
  "GNE-495" = "red",
  "JNK-inh" = "blue",
  "none" = "grey"
)

p <- fgsea_rescue_combo_vs_combo %>%
  arrange(match(sig_class, names(sig_class_colors))) %>%
  ggplot(
    aes(
      x = `signed_padj_GNE-495`,
      y = `signed_padj_JNK-inh`
    )
  ) +
  geom_point(
    aes(
      color = sig_class,
      text = pathway
    ),
    alpha = .8
  ) +
  scale_color_manual(
    values = sig_class_colors
  ) +
  labs(
    x = paste("Enrichment GNE-495 rescue"),
    y = paste("Enrichment JNK-inh rescue")
  )

ggsave(
  file.path("plots", paste0("reactome_rescue_gne_vs_jnk_scatter_padj.pdf")),
  p, width = 8, height = 8
)
htmlwidgets::saveWidget(
  ggplotly(p),
  file.path("plots", paste0("reactome_rescue_gne_vs_jnk_scatter_padj.html"))
)


```

### Comparing rescue vector enrichment against PTX alone

```{r}
fgsea_rescue_combo_vs_tx <- fgsea_res_rescue_data %>%
  mutate(
    signed_padj = -sign(NES) * log10(padj)
  ) %>%
  power_inner_join(
    fgsea_res_pert %>%
      filter(
        contrast == "PTX",
        metric == "signed_padj"
      ) %>%
      mutate(
        signed_padj = -sign(NES) * log10(padj)
      ) %>%
      separate_wider_delim(
        pathway,
        delim = "_",
        names = c("database", "pathway"),
        too_many = "merge"
      ),
    by = c("database", "pathway"),
    suffix = c("_rescue", "_PTX"),
    check = check_specs(
      unmatched_keys_left = "warn",
      unmatched_keys_right = "warn",
      duplicate_keys_right = "warn"
    )
  )

p <- fgsea_rescue_combo_vs_tx %>%
  ggplot(
    aes(
      x = signed_padj_PTX,
      y = signed_padj_rescue
    )
  ) +
  geom_point(
    aes(
      text = pathway
    )
  ) +
  facet_wrap(~compound)


htmlwidgets::saveWidget(
  ggplotly(p),
  file.path("plots", "reactome_rescue_vs_ptx_scatter_padj.html")
)
```


## Plot LFC heatmap of genes in enriched pathways PTX vs combo
```{r}
library(msigdbr)
all_msigdbr <- msigdbr()

msigdbr_of_interest <- all_msigdbr %>%
  filter(
    gs_cat == "H" |
    gs_subcat %in% c(
      "CP:PID",
      "CP:KEGG",
      "CP:REACTOME"
      # "GO:BP",
      # "GO:MF"
    ),
    ensembl_gene %in% feasible_genes$gene_id
  ) %>%
  separate_wider_delim(
    gs_name,
    delim = "_",
    names = c("database", "pathway"),
    too_many = "merge"
  )
```


```{r}
fgsea_hm_conditions <- c(
  "PTX", "GNE-495", "JNK-inh",
  "GNE-495 + PTX", "JNK-inh + PTX"
)

de_pert_hm_data <- de_long_selected %>%
  filter(
    contrast %in% fgsea_hm_conditions
  ) %>%
  mutate(
    signed_padj = -sign(logFC) * log10(FDR_min_zero)
  )

library(ComplexHeatmap)
de_pert_rescue_hm_data <- tibble(
  pathway = c(
    "RECYCLING_PATHWAY_OF_L1",
    "REGULATION_OF_EXPRESSION_OF_SLITS_AND_ROBOS",
    "SRP_DEPENDENT_COTRANSLATIONAL_PROTEIN_TARGETING_TO_MEMBRANE",
    # "EPH_EPHRIN_SIGNALIN",
    "EPHB_MEDIATED_FORWARD_SIGNALING",
    "P38MAPK_EVENTS",
    "RHO_GTPASE_EFFECTORS"
  )
) %>%
  mutate(
    genes = map(
      pathway,
      \(x) msigdbr_of_interest %>%
        filter(
          pathway == x
        ) %>%
        pull(ensembl_gene) %>%
        unique()
    ),
    mat = map(
      genes,
      \(x) de_pert_hm_data %>%
        filter(
          gene_id %in% x
        ) %>%
        select(gene_id, contrast, logFC) %>%
        pivot_wider(names_from = contrast, values_from = logFC, values_fill = 0) %>%
        column_to_rownames("gene_id") %>%
        as.matrix() %>% {
          .[, fgsea_hm_conditions]
        }
    ),
    row_meta = map2(
      pathway, mat,
      \(x, y) {
        fgsea_res_rescue_data %>%
          filter(pathway == x) %>%
          select(compound, gene_id = leadingEdge) %>%
          unchop(gene_id) %>%
          mutate(
            gene_id = factor(gene_id, levels = rownames(y)),
            decoy = 1L
          ) %>%
          pivot_wider(
            names_from = compound, values_from = decoy,
            names_prefix = "leading_edgge_",
            id_expand = TRUE, values_fill = 0L
          ) %>%
          column_to_rownames("gene_id") %>% {
            .[rownames(y), ]
          }
      }
    ),
    hm = pmap(
      list(pathway, mat, row_meta),
      \(x, y, z) {
        mat_max_abs <- max(abs(quantile(y, c(.025, .975), na.rm = TRUE)))
        Heatmap(
          y,
          name = "log2 fold change",
          col = circlize::colorRamp2(seq(from = -mat_max_abs, to = mat_max_abs, length.out = 51), paletteer::paletteer_c("pals::ocean.balance", n = 51)),
          cluster_rows = cluster_fun_eucl,
          cluster_columns = FALSE,
          right_annotation = HeatmapAnnotation(
            df = z,
            which = "row"
          )
        )
      }
    )
  )

dir.create(
  here("plots", "de_heatmaps"),
  showWarnings = FALSE
)

pwalk(
  de_pert_rescue_hm_data,
  \(pathway, hm, ...) {
    withr::with_pdf(
      file.path("plots", "de_heatmaps", paste0("de_heatmap_", pathway, ".pdf")),
      draw(hm),
      width = 6, height = 10
    )
  }
)
```

## Overrepresentation analysis of rescue gene sets

### Plotting single conditions as bars

```{r}
fgsea_overrep_data <- fgsea_overrep %>%
  separate_wider_delim(
    pathway,
    delim = "_",
    names = c("database", "pathway"),
    too_many = "merge"
  ) %>%
  mutate(
    # sign = if_else(NES > 0, "pos", "neg"),
    signed_padj = -if_else(perturbed == "up", 1, -1) * log10(padj_fgsea),
    pathway_db = paste(str_sub(database, 1, 1), pathway, sep = "_"),
    pathway_db_direction = paste(perturbed, pathway_db, sep = "_")
  )

ps <- fgsea_overrep_data %>%
  group_nest(
    compound, rescue_definition
  ) %>%
  mutate(
    p = pmap(
      list(compound, data, rescue_definition),
      \(x, y, z) {
        y %>%
          filter(
            padj_fgsea < .01
          ) %>%
          mutate(
            across(pathway_db_direction, \(x) fct_reorder(x, signed_padj))
          ) %>%
          ggplot(
            aes(
              x = signed_padj,
              y = pathway_db_direction,
              fill = perturbed
            )
          ) +
          geom_col() +
          geom_text(
            aes(
              x = if_else(perturbed == "up", -.05, .05),
              hjust = if_else(perturbed == "up", 1, 0),
              label = pathway_db
            ),
            size = 2
          ) +
          scale_fill_discrete(guide = "none", direction = -1) +
          theme(
            axis.text.y = element_blank(),
            panel.grid.major.y = element_blank(),
          ) +
          labs(
            x = "Signed adjusted p-value",
            y = NULL,
            title = paste0("Pathways enriched in genes ", z, " with ", x),
            subtitle = "FDR < 0.01"
          )
      }
    )
  )

dir.create(
  here("plots", "overrep"),
  showWarnings = FALSE
)

pwalk(
  ps,
  \(compound, p, rescue_definition, ...) {
    n <- nrow(layer_data(p))
    h <- 2 + .1 * n
    ggsave(
      file.path("plots", "overrep", paste0("overrep_", rescue_definition, "_", compound, "_bars.pdf")),
      p, width = 8, height = h
    )
  }
)

```

## vs PTX

```{r}
fgsea_overrep_vs_ptx_pert <- fgsea_overrep %>%
  filter(
    rescue_definition != "perturbed"
  ) %>%
  separate_wider_delim(
    pathway,
    delim = "_",
    names = c("database", "pathway"),
    too_many = "merge"
  ) %>%
  mutate(
    signed_padj = -if_else(perturbed == "up", 1, -1) * log10(padj)
  ) %>%
  power_inner_join(
    fgsea_res_single_treatments %>%
      filter(contrast == "PTX", metric == "signed_padj") %>%
      mutate(
        signed_padj = -if_else(NES > 0, 1, -1) * log10(padj)
      ),
    by = c("database", "pathway"),
    suffix = c("_overrep", "_PTX"),
    check = check_specs(
      unmatched_keys_left = "warn",
      unmatched_keys_right = "warn",
      duplicate_keys_right = "warn"
    )
  )

p <- fgsea_overrep_vs_ptx_pert %>%
  ggplot(
    aes(
      x = signed_padj_PTX,
      y = signed_padj_overrep
    )
  ) +
  geom_point(
    aes(
      text = pathway
    )
  ) +
  facet_grid(
   vars(compound), vars(rescue_definition)
  )

htmlwidgets::saveWidget(
  ggplotly(p),
  file.path("plots", "overrep", "overrep_vs_fgsea_ptx_scatter_padj.html")
)
```

## Plotting intersections of rescue gene sets

Have to be careful about perturbed `both_up_down` class, messes up signed_padj

### intersections bar charts

```{r}
fgsea_intersection_overrep_data <- fgsea_intersection_overrep %>%
  separate_wider_delim(
    pathway,
    delim = "_",
    names = c("database", "pathway"),
    too_many = "merge"
  ) %>%
  mutate(
    # sign = if_else(NES > 0, "pos", "neg"),
    signed_padj = -if_else(perturbed %in% c("both_up_down", "up"), 1, -1) * log10(padj_fgsea),
    pathway_db = paste(str_sub(database, 1, 1), pathway, sep = "_"),
    pathway_db_direction = paste(perturbed, pathway_db, sep = "_")
  )

fgsea_intersection_overrep_data %>%
  filter(
    rescue_definition == "significant_opposite",
    pathway == "SRP_DEPENDENT_COTRANSLATIONAL_PROTEIN_TARGETING_TO_MEMBRANE"
  ) %>%
  View()

ps <- fgsea_intersection_overrep_data %>%
  filter(perturbed %in% c("up", "down")) %>%
  group_nest(
    rescue_type, rescue_definition
  ) %>%
  mutate(
    p = pmap(
      list(rescue_type, rescue_definition, data),
      \(x, y, z) {
        z %>%
          filter(
            padj_fgsea < .01
          ) %>%
          mutate(
            across(pathway_db_direction, \(x) fct_reorder(x, signed_padj))
          ) %>%
          ggplot(
            aes(
              x = estimate,
              y = pathway_db_direction,
              fill = perturbed
            )
          ) +
          geom_col() +
          geom_text(
            aes(
              x = if_else(perturbed == "up", -.05, .05),
              hjust = if_else(perturbed == "up", 1, 0),
              label = pathway_db
            ),
            size = 2
          ) +
          scale_fill_discrete(guide = "none", direction = -1) +
          theme(
            axis.text.y = element_blank(),
            panel.grid.major.y = element_blank(),
          ) +
          labs(
            x = "Signed adjusted p-value",
            y = NULL,
            title = paste0("Pathways enriched in genes rescued by ", x),
            subtitle = paste("FDR < 0.01", y)
          )
      }
    )
  )

dir.create(
  here("plots", "overrep"),
  showWarnings = FALSE
)

pwalk(
  ps,
  \(rescue_definition, rescue_type, p, ...) {
    ggsave(
      file.path("plots", "overrep", paste0("overrep_intersection_", rescue_definition, "_", rescue_type, "_bars.pdf")),
      p, width = 8, height = 8
    )
  }
)

```

### intersections versus each other and vs PTX


```{r}
fgsea_intersection_vs_ptx <- fgsea_intersection_overrep_data %>%
  power_inner_join(
    fgsea_res_pert %>%
      separate_wider_delim(
        pathway,
        delim = "_",
        names = c("database", "pathway"),
        too_many = "merge"
      ) %>%
      mutate(
        signed_padj = -sign(NES) * log10(padj)
      ) %>%
      filter(contrast == "PTX", metric == "signed_padj"),
    by = c("database", "pathway"),
    suffix = c("_rescue", "_ptx"),
    check = check_specs(
      unmatched_keys_left = "warn",
      unmatched_keys_right = "warn",
      duplicate_keys_right = "warn"
    )
  )

ps <- fgsea_intersection_vs_ptx %>%
  filter(perturbed %in% c("up", "down")) %>%
  group_nest(
    rescue_definition
  ) %>%
  rowwise() %>%
  mutate(
    p = list(
    ggplot(
      data,
      aes(
        signed_padj_ptx,
        signed_padj_rescue,
        color = padj_fgsea < .05
      )
    ) +
    geom_point(
      aes(
        text = paste(
          pathway,
          database,
          "\n",
          paste(
            "padj rescue:", signif(padj_fgsea, 2),
            "overlap:", overlap,
            "gene set size:", size_ptx,
            "odds ratio rescue:", signif(estimate, 2),
            "padj ptx:", signif(padj, 2),
            "NES:", signif(NES, 2)
          )
        )
      ),
      alpha = .5,
      shape = 16
    ) +
    scale_color_manual(
      values = c(
        `TRUE` = "red",
        `FALSE` = "black"
        # up = "red",
        # down = "blue",
        # both_up_down = "black"
      )
    ) +
    ggh4x::facet_wrap2(~rescue_type, scales = "free", axes = "y")
    )
  ) %>%
  ungroup()

pwalk(
  ps,
  \(rescue_definition, p, ...) {
    ggsave(
      file.path("plots", "overrep", paste0("overrep_intersection_vs_ptx_", rescue_definition, ".pdf")),
      p, width = 8, height = 8
    )
    htmlwidgets::saveWidget(
      ggplotly(p),
      file.path("plots", "overrep", paste0("overrep_intersection_vs_ptx_", rescue_definition, ".html"))
    )
  }
)
```


### Heatmap of intersection enrichments

```{r}
intersection_overrep_hm_rescure_type_order <- c(
  "both", "GNE-495 only", "JNK-inh only", "no rescue"
)

de_pert_hm_data <- de_long_selected %>%
  filter(
    contrast %in% fgsea_hm_conditions
  ) %>%
  mutate(
    signed_padj = -sign(logFC) * log10(FDR_min_zero)
  )

library(ComplexHeatmap)
de_pert_rescue_hm_data <- tibble(
  pathway = c(
    "RECYCLING_PATHWAY_OF_L1",
    "REGULATION_OF_EXPRESSION_OF_SLITS_AND_ROBOS",
    "SRP_DEPENDENT_COTRANSLATIONAL_PROTEIN_TARGETING_TO_MEMBRANE",
    # "EPH_EPHRIN_SIGNALIN",
    "EPHB_MEDIATED_FORWARD_SIGNALING",
    "P38MAPK_EVENTS",
    "RHO_GTPASE_EFFECTORS",
    "RAC1_PATHWAY",
    "RHO_GTPASES_ACTIVATE_ROCKS"
  )
) %>%
  mutate(
    genes = map(
      pathway,
      \(x) msigdbr_of_interest %>%
        filter(
          pathway == x
        ) %>%
        pull(ensembl_gene) %>%
        unique()
    ),
    mat = map(
      genes,
      \(x) de_pert_hm_data %>%
        filter(
          gene_id %in% x
        ) %>%
        select(gene_id, con trast, logFC) %>%
        pivot_wider(names_from = contrast, values_from = logFC, values_fill = 0) %>%
        column_to_rownames("gene_id") %>%
        as.matrix() %>% {
          .[, fgsea_hm_conditions]
        }
    ),
    row_meta = map2(
      pathway, mat,
      \(x, y) {
        fgsea_res_rescue_data %>%
          filter(pathway == x) %>%
          select(compound, gene_id = leadingEdge) %>%
          unchop(gene_id) %>%
          mutate(
            gene_id = factor(gene_id, levels = rownames(y)),
            decoy = 1L
          ) %>%
          pivot_wider(
            names_from = compound, values_from = decoy,
            names_prefix = "leading_edgge_",
            id_expand = TRUE, values_fill = 0L
          ) %>%
          column_to_rownames("gene_id") %>% {
            .[rownames(y), ]
          }
      }
    ),
    hm = pmap(
      list(pathway, mat, row_meta),
      \(x, y, z) {
        mat_max_abs <- max(abs(quantile(y, c(.025, .975), na.rm = TRUE)))
        Heatmap(
          y,
          name = "log2 fold change",
          col = circlize::colorRamp2(seq(from = -mat_max_abs, to = mat_max_abs, length.out = 51), paletteer::paletteer_c("pals::ocean.balance", n = 51)),
          cluster_rows = cluster_fun_eucl,
          cluster_columns = FALSE,
          right_annotation = HeatmapAnnotation(
            df = z,
            which = "row"
          )
        )
      }
    )
  )

dir.create(
  here("plots", "de_heatmaps"),
  showWarnings = FALSE
)

pwalk(
  de_pert_rescue_hm_data,
  \(pathway, hm, ...) {
    withr::with_pdf(
      file.path("plots", "de_heatmaps", paste0("de_heatmap_", pathway, ".pdf")),
      draw(hm),
      width = 6, height = 10
    )
  }
)

```

## TopGO

check if sometimes the same pathway is enriched in both up and downregulated genes

```{r}
topgo_res_pert %>%
  pivot_wider(
    names_from = direction,
    values_from = -c(contrast, direction, ontology, GO.ID, Term, term_unique, Annotated)
  ) %>%
  filter(
    fisher_double_down < .05,
    fisher_double_up < .05
  ) %>%
  View()
```

Make directed topgo results object by taking the log10 p-value
and the sign corresponds to the direction of the enrichment

If the same pathway is enriched in both directions, the pathway
occurs twice with different signs

```{r}
topgo_res_pert_directed <- topgo_res_pert %>%
  mutate(
    log_odds_ratio = {Significant / Expected} %>%
      if_else(
        !is.finite(.) | . == 0,
        sign(.) * min(.[. > 0]) * .5,
        .
      ) %>%
      log10()
  ) %>%
  group_by(
    contrast, ontology, GO.ID, Term, term_unique, Annotated
  ) %>%
  slice(
    if (sum(fisher_double < .05) < 2) {
      order(fisher_double)[1]
    } else {
      order(fisher_double)
    }
  ) %>%
  mutate(
    picked_direction = if (n() < 2) direction else "both",
    signed_padj = -sign(log_odds_ratio) * log10(fisher_double)
  ) %>%
  ungroup()

```


```{r}
all_contrasts <- topgo_res_pert_directed$contrast %>%
  unique()

contrast_plotting_pairs <- bind_rows(
  tibble(
    contrast_1 = "PTX",
    contrast_2 = all_contrasts[
      str_ends(all_contrasts, fixed("+ PTX")) |
      str_ends(all_contrasts, fixed("+ PTX vs PTX"))
    ]
  ),
  tibble(
    contrast_1 = "GNE-495 + PTX",
    contrast_2 = "JNK-inh + PTX"
  )
)
```

```{r}
topgo_res_ptx_vs_combo <- topgo_res_pert_directed %>%
  filter(
    ontology == "BP"
  ) %>% {
    d <- .
    power_inner_join(
      contrast_plotting_pairs,
      d,
      by = join_by(contrast_1 == contrast),
      suffix = c("", "_1"),
      check = check_specs(
        unmatched_keys_left = "warn"
      )
    ) %>%
    power_inner_join(
      d,
      by = join_by(contrast_2 == contrast, ontology, GO.ID, Term, term_unique, Annotated),
      suffix = c("_1", "_2"),
      check = check_specs(
        unmatched_keys_left = "warn"
      )
    )
  } %>%
  mutate(
    sig_class = case_when(
      fisher_double_1 < .05 & fisher_double_2 < .05 ~ "both",
      fisher_double_1 < .05 ~ "1",
      fisher_double_2 < .05 ~ "2",
      TRUE ~ "none"
    )
  )

sig_class_colors <- c(
  "both" = "purple",
  "1" = "red",
  "2" = "blue",
  "none" = "grey"
)

ps <- topgo_res_ptx_vs_combo %>%
  group_nest(contrast_1, contrast_2) %>%
  mutate(
    p = pmap(
      .,
      \(contrast_1, contrast_2, data, ...) {
        ggplot(
          data %>%
            arrange(match(sig_class, names(sig_class_colors))),
          aes(
            x = signed_padj_1,
            y = signed_padj_2
          )
        ) +
        geom_point(
          aes(
            color = sig_class,
            text = Term
          ),
          alpha = .8
        ) +
        scale_color_manual(
          values = sig_class_colors
        ) +
        labs(
          x = paste("Enrichment", contrast_1),
          y = paste("Enrichment", contrast_2)
        )
      }
    )
  )

pwalk(
  ps,
  \(contrast_1, contrast_2, p, ...) {
    ggsave(
      file.path("plots", paste0("topgo_", contrast_1, "_vs_", contrast_2, "_scatter_padj.pdf")),
      p, width = 8, height = 8
    )
    htmlwidgets::saveWidget(
      ggplotly(p),
      file.path("plots", paste0("topgo_", contrast_1, "_vs_", contrast_2, "_scatter_padj.html"))
    )
  }
)


```

## Create Excel table

1. Differential expression
2. FGSEA on raw perturbations
3. FGSEA on rescue vectors
4. Overrepresentation on rescue gene sets
5. Overrepresentation on rescue gene set intersections

```{r}
gene_id_symbol_map <- de_long_selected %>%
  distinct(gene_id, gene_name)

gene_id_symbol_vec <- with(
  gene_id_symbol_map,
  set_names(gene_name, gene_id)
)

excel_de <- de_long_selected %>%
  select(-FDR_min_zero)

excel_fgsea_raw <- fgsea_res_pert %>%
  filter(metric == "signed_padj") %>%
  separate_wider_delim(
    pathway,
    delim = "_",
    names = c("database", "pathway"),
    too_many = "merge"
  ) %>%
  select(-c(metric, log2err, ES)) %>%
  mutate(
    contrast = if_else(str_ends(contrast, fixed("vs PTX")), contrast, paste(contrast, "vs DMSO")),
    leadingEdge = map_chr(
      leadingEdge,
      \(x) str_flatten(gene_id_symbol_vec[x], collapse = ",")
    )
  )

excel_fgsea_rescue <- fgsea_res_rescue %>%
  filter(metric == "signed_padj_rescue_clamped") %>%
  separate_wider_delim(
    pathway,
    delim = "_",
    names = c("database", "pathway"),
    too_many = "merge"
  ) %>%
  select(-c(metric, log2err, ES)) %>%
  mutate(
    leadingEdge = map_chr(
      leadingEdge,
      \(x) str_flatten(gene_id_symbol_vec[x], collapse = ",")
    )
  )

excel_overrep_rescue <- fgsea_overrep %>%
  filter(rescue_definition == "significant_opposite") %>%
  separate_wider_delim(
    pathway,
    delim = "_",
    names = c("database", "pathway"),
    too_many = "merge"
  ) %>%
  transmute(
    compound, ptx_direction = perturbed,
    database, pathway,
    pval = p.value_fisher, padj = padj_fisher,
    odds_ratio = estimate,
    size, overlap,
    leadingEdge = map_chr(
      overlapGenes,
      \(x) str_flatten(gene_id_symbol_vec[x], collapse = ",")
    )
  )

excel_overrep_intersections <- fgsea_intersection_overrep %>%
  filter(rescue_definition == "significant_opposite") %>%
  separate_wider_delim(
    pathway,
    delim = "_",
    names = c("database", "pathway"),
    too_many = "merge"
  ) %>%
  transmute(
    rescue_type, ptx_direction = perturbed,
    database, pathway,
    pval = p.value_fisher, padj = padj_fisher,
    odds_ratio = estimate,
    size, overlap,
    leadingEdge = map_chr(
      overlapGenes,
      \(x) str_flatten(gene_id_symbol_vec[x], collapse = ",")
    )
  )
```

```{r}
library(openxlsx)

openxlsx::write.xlsx(
  list(
    "Differential expression" = excel_de,
    "FGSEA perturbations" = excel_fgsea_raw,
    "FGSEA rescue vectors" = excel_fgsea_rescue,
    "Overrep rescue gene sets" = excel_overrep_rescue,
    "Overrep rescue intersections" = excel_overrep_intersections
  ),
  file = here("results", "gene_enrichment_analysis_results.xlsx"),
  asTable = TRUE
)
```


## RAGs and ion channels

```{r}
library(biomaRt)

# mouse = useMart("ensembl", dataset = "mmusculus_gene_ensembl")
human = useMart("ensembl", dataset = "hsapiens_gene_ensembl")

# Install packages if you haven't already
# install.packages("httr")
# install.packages("jsonlite")

library(httr)
library(jsonlite)

# Define a function to query mygene.info for a given gene symbol (or alias)
query_mygene <- function(query_term, scopes = "symbol,alias", species = "mouse", fields = "ensembl.gene") {
  base_url <- "http://mygene.info/v3/query"

  # Build the query parameters
  params <- list(
    q = query_term,
    scopes = scopes,
    species = species,
    fields = fields
  )

  # Make the GET request
  res <- POST(url = base_url, body = params, encode = "json")

  # Check for a successful request (status code 200)
  if (status_code(res) == 200) {
    # Parse the JSON content
    data <- fromJSON(content(res, as = "text", encoding = "UTF-8"))
    return(data)
  } else {
    stop("Query failed with status: ", status_code(res))
  }
}

# # mouse_gene_id_mapping1 <- getBM(
# #   attributes = c("ensembl_gene_id", "external_gene_name"),
# #   filters = c("external_gene_name"),
# #   values = unique(rag_ion_channel_genes$gene_symbol_mouse),
# #   mart = mouse
# # )

# # mouse_gene_id_mapping2 <- getBM(
# #   attributes = c("ensembl_gene_id", "mgi_symbol"),
# #   filters = c("mgi_symbol"),
# #   values = unique(rag_ion_channel_genes$gene_symbol_mouse),
# #   mart = mouse
# # )

# mouse_gene_id_mapping <- bind_rows(
#   synonym = dplyr::select(
#     mouse_gene_id_mapping1,
#     ensembl_gene_id_mouse = ensembl_gene_id,
#     gene_symbol_mouse = external_gene_name
#   ),
#   mgi_symbol = dplyr::select(
#     mouse_gene_id_mapping2,
#     ensembl_gene_id_mouse = ensembl_gene_id,
#     gene_symbol_mouse = mgi_symbol
#   )
# )

mouse_gene_id_mapping <- query_mygene(
  unique(rag_ion_channel_genes$gene_symbol_mouse)
) %>%
  transmute(
    gene_symbol_mouse = query,
    ensembl_gene_id_mouse = map(
      ensembl,
      \(x) if (class(x) == "list")
        tibble(gene = x[[1]])
      else
        x
    )
  ) %>%
  unnest(ensembl_gene_id_mouse)

mouse_homologue_mapping <- getHomologs(
  unique(mouse_gene_id_mapping$gene),
  "mus_musculus",
  "homo_sapiens"
)

human_gene_id_symbol_map <- getBM(
  attributes = c("ensembl_gene_id", "hgnc_symbol"),
  filters = "ensembl_gene_id",
  values = unique(na.omit(mouse_homologue_mapping$hsapiens_homolog_ensembl_gene)),
  mart = human
)

rag_ion_channel_genes_mapped <- rag_ion_channel_genes %>%
  left_join(
    mouse_gene_id_mapping %>%
      dplyr::distinct(gene_symbol_mouse, ensembl_gene_id_mouse = gene),
    by = "gene_symbol_mouse"
  ) %>%
  left_join(
    mouse_homologue_mapping %>%
      dplyr::select(
        ensembl_gene_id_mouse = ensembl_gene_id,
        ensembl_gene_id_human = hsapiens_homolog_ensembl_gene
      ),
    by = "ensembl_gene_id_mouse"
  ) %>%
  left_join(
    human_gene_id_symbol_map %>%
      dplyr::select(
        ensembl_gene_id_human = ensembl_gene_id,
        hgnc_symbol
      ),
    by = "ensembl_gene_id_human"
  ) %>%
  mutate(
    gene_display_name = coalesce(
      hgnc_symbol,
      gene_symbol_mouse,
      ensembl_gene_id_human,
      ensembl_gene_id_mouse
    )
  ) %>%
  group_by(gene_display_name) %>%
  mutate(
    unique_gene_id = if (n() > 1)
      paste0(gene_display_name, "-", seq_len(n()))
    else
      gene_display_name
  ) %>%
  ungroup()

```


```{r}
rag_hm_conditions <- c(
  "PTX", "GNE-495", "JNK-inh",
  "GNE-495 + PTX", "JNK-inh + PTX"
)

rag_ion_channel_genes_valid <- rag_ion_channel_genes_mapped %>%
  drop_na(
    ensembl_gene_id_human
  ) %>%
  filter(ensembl_gene_id_human %in% de_long_selected$gene_id)

de_rag_hm_data <- de_long_selected %>%
  filter(
    contrast %in% rag_hm_conditions
  ) %>%
  inner_join(
    rag_ion_channel_genes_valid,
    by = c("gene_id" = "ensembl_gene_id_human")
  ) %>%
  mutate(
    signed_padj = -sign(logFC) * log10(FDR_min_zero)
  )

de_rag_hm_mat <- de_rag_hm_data %>%
  dplyr::select(
    unique_gene_id, contrast, logFC
  ) %>%
  pivot_wider(names_from = contrast, values_from = logFC) %>%
  column_to_rownames("unique_gene_id") %>%
  as.matrix() %>% {
    .[
      rag_ion_channel_genes_valid$unique_gene_id,
      rag_hm_conditions
    ]
  }

library(ComplexHeatmap)
hm <- Heatmap(
  de_rag_hm_mat,
    # t() %>%
    # scale() %>%
    # t(),
  cluster_rows = cluster_fun_eucl,
  cluster_columns = FALSE,
  show_row_names = TRUE,
  row_split = rag_ion_channel_genes_valid$gene_set
)

withr::with_pdf(
  file.path("plots", "de_heatmaps", "de_heatmap_rag_ion_channels_lfc.pdf"),
  draw(hm),
  width = 6, height = 10
)

de_rag_hm_mat <- de_rag_hm_data %>%
  dplyr::select(
    unique_gene_id, contrast, signed_padj
  ) %>%
  pivot_wider(names_from = contrast, values_from = signed_padj) %>%
  column_to_rownames("unique_gene_id") %>%
  as.matrix() %>% {
    .[
      rag_ion_channel_genes_valid$unique_gene_id,
      rag_hm_conditions
    ]
  }

library(ComplexHeatmap)
hm <- Heatmap(
  # de_rag_hm_mat,
  # Clamp to -10 to 10
  de_rag_hm_mat %>% {
    pmax(pmin(., 10), -10)
  },
  name = "-sign(log2FC) * log10(padj)",

    # t() %>%
    # scale() %>%
    # t(),
  col = circlize::colorRamp2(seq(-10, 10, length.out = 51), paletteer::paletteer_c("pals::ocean.balance", n = 51)),
  cluster_rows = cluster_fun_eucl,
  cluster_columns = FALSE,
  show_row_names = FALSE,
  row_split = rag_ion_channel_genes_valid$gene_set
)

withr::with_pdf(
  file.path("plots", "de_heatmaps", "de_heatmap_rag_ion_channels_padj.pdf"),
  draw(hm),
  width = 6, height = 10
)

hm <- Heatmap(
  de_rag_hm_mat,
  # Clamp to -10 to 10
  # de_rag_hm_mat %>% {
  #   pmax(pmin(., 10), -10)
  # },
  name = "-sign(log2FC) * log10(padj)",

    # t() %>%
    # scale() %>%
    # t(),
  col = circlize::colorRamp2(seq(-10, 10, length.out = 51), paletteer::paletteer_c("pals::ocean.balance", n = 51)),
  cluster_rows = cluster_fun_eucl,
  cluster_columns = FALSE,
  show_row_names = FALSE,
  row_split = rag_ion_channel_genes_valid$gene_set
)

withr::with_pdf(
  file.path("plots", "de_heatmaps", "de_heatmap_rag_ion_channels_padj2.pdf"),
  draw(hm),
  width = 6, height = 10
)
```

Overlap of RAGs and ion channels with GNE-495-only rescue genes


```{r}



```