A human iPSC-derived sensory neuron platform for high-throughput discovery of neuroprotectants against chemotherapy-induced peripheral neuropathy.

Code ↔ Paper

9 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 9 matches

[1] § Results › Molecular profiling of paclitaxel-induced neurotoxicity ↔ rnaseq/plot_dge.Rmd, lines 49–77 · score 0.88 · TUBA1A, TUBB2B, TUBB4B, ACAT2, ACSS2, ASS1
[2] § Results › Molecular profiling of paclitaxel-induced neurotoxicity ↔ gene_set_analysis_plotting.Rmd, lines 393–516 · score 0.82 · fatty acid metabolism, gap junctions, Rho GTPase, Golgi, trafficking, vesicle
[3] § Results › Molecular profiling of paclitaxel-induced neurotoxicity ↔ gene_set_analysis_plotting.Rmd, lines 393–516 · score 0.64 · extracellular matrix, gene expression, JNK Inh, interactions, metabolism, Lipid
[4] § STAR★Methods › Method details › Gene set enrichment analysis from RNA sequencing data ↔ gene_set_analysis_plotting.Rmd, lines 3148–3207 · score 0.57 · KEGG medicus, log2 fold, Hallmark, Reactome, enrichment, gene
[5] § STAR★Methods › Method details › siRNA transfection and qPCR validation ↔ marker_expression.Rmd, lines 74–134 · score 0.56 · MAP4K6, MAP4K4, sensory neurons, TNIK
[6] § Results › Molecular profiling of paclitaxel-induced neurotoxicity ↔ gene_set_analysis_plotting.Rmd, lines 759–883 · score 0.54 · fatty acid, cholesterol, apoptosis, metabolism, phospholipids, activation
[7] § Results › Molecular profiling of paclitaxel-induced neurotoxicity ↔ gene_set_analysis_plotting.Rmd, lines 759–883 · score 0.54 · junction assembly, modification, GTPase, regulation, phosphoprotein, mediated
[8] § STAR★Methods › Method details › Identification of differentially expressed genes ↔ rnaseq/deseq.Rmd, lines 868–939 · score 0.51 · sequencing batch, DESeq2, Wald, Salmon, gene
[9] § Results › High-throughput screen of kinase inhibitors identified 19 neuroprotective compounds and three STE20 kinases as drivers of axon degeneration ↔ marker_expression.Rmd, lines 74–134 · score 0.51 · MAP4K6, MAP4K4, sensory neuron, TNIK

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

R Markdown · 3,506 lines · 86 KB · no license · 5 matches

---
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("syn70748503") %>%
read_csv()
fgsea_overrep <- syn("syn64582959") %>%
qread()
fgsea_overrep_single_treatments <- syn("syn71710650") %>%
qread()
fgsea_overrep_both_doses <- syn("syn71703313") %>%
read_csv()
fgsea_overrep_clusters <- syn("syn71761437") %>%
qread()
fgsea_intersection_overrep <- syn("syn64586672") %>%
qread()
indra_res <- syn("syn64908970") %>%
read_csv()
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()
```
```{r}
database_abbreviations <- c(
"REACTOME" = "R",
"HALLMARK" = "H",
"PID" = "P",
"KEGG_MEDICUS" = "K",
"GOBP" = "BP",
"GOMF" = "MF",
"GOCC" = "CC"
)
abbreviate_prefix <- function(strings, abbr_map = database_abbreviations) {
# Create regex pattern from abbreviation names
pattern <- paste0("^(", paste(names(abbr_map), collapse = "|"), ")_")
# Replace matching prefixes with abbreviations
str_replace(strings, pattern, function(x) {
prefix <- str_remove(x, "_$")
paste0(abbr_map[prefix], "_")
})
}
recode_batch <- \(x) recode(x, batch_1 = "(3uM)", batch_2 = "(1uM)", batch_1_riki = "(3uM, Riki)")
```
## Plot Reactome PTX only and comparison with combination treatment
```{r}
fgsea_res_single_treatments <- fgsea_res_pert %>%
filter(
metric == "signed_padj"
) %>%
mutate(
contrast_batch = paste(contrast, recode_batch(batch)),
database = str_extract(pathway, "^(REACTOME|HALLMARK|PID|KEGG_MEDICUS|GOBP|GOMF|GOCC)"),
pathway = abbreviate_prefix(pathway)
) %>%
filter(
padj < 0.01
) %>%
mutate(
sign = if_else(NES > 0, "pos", "neg")
)
ps <- fgsea_res_single_treatments %>%
group_nest(
batch, contrast_batch, contrast
) %>%
mutate(
p = map2(
contrast_batch, 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("Pathways enriched in ", x, " treatment"),
subtitle = "FDR < 0.01"
) + {
if (!any(y$sign == "neg"))
lims(x = c(-max(abs(y$NES)), NA))
}
}
)
)
dir.create(
here("plots", "gene_set_analysis"),
showWarnings = FALSE
)
pwalk(
ps,
\(contrast_batch, p, ...) {
ggsave(
here("plots", "gene_set_analysis", paste0("fgsea_single_perturbation_", contrast_batch, "_bars.pdf")),
p, width = 8, height = 8
)
}
)
ps <- fgsea_res_single_treatments %>%
filter(!str_starts(database, coll("GO"))) %>%
group_nest(
batch, contrast_batch, contrast
) %>%
mutate(
p = map2(
contrast_batch, 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("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_batch, p, ...) {
ggsave(
here("plots", "gene_set_analysis", paste0("fgsea_single_perturbation_", contrast_batch, "_no_go_bars.pdf")),
p, width = 8, height = 8
)
}
)
```
### Compare batch 1 vs batch 2
```{r}
fgsea_res_pert_vs <- fgsea_res_pert %>%
mutate(
database = str_extract(pathway, "^(REACTOME|HALLMARK|PID|KEGG_MEDICUS)"),
pathway = str_remove(pathway, "^(REACTOME|HALLMARK|PID|KEGG_MEDICUS)_"),
sign = if_else(NES > 0, "pos", "neg"),
signed_padj = -sign(NES) * log10(padj)
) %>%
pivot_wider(
names_from = batch,
values_from = c(pval, padj, log2err, ES, NES, leadingEdge, sign, size, signed_padj)
)
p <- fgsea_res_pert_vs %>%
ggplot(
aes(
x = signed_padj_batch1,
y = signed_padj_batch2
)
) +
geom_point(
shape = 16, alpha = .7
) +
facet_grid(
contrast ~ metric
)
```
### 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)
) %>% {
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, batch),
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(batch, 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,
\(batch, contrast_1, contrast_2, p, ...) {
ggsave(
here("plots", "gene_set_analysis", paste0("fgsea_single_perturbation_", contrast_1, "_vs_", contrast_2, "_", batch, "_scatter_padj.pdf")),
p, width = 7, height = 6
)
htmlwidgets::saveWidget(
ggplotly(p),
here("plots", "gene_set_analysis", paste0("fgsea_single_perturbation_", contrast_1, "_vs_", contrast_2, "_", batch, "_scatter_padj.html")),
selfcontained = FALSE
)
}
)
```
### 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"))
# )
}
)
```
### Dotplot of single comparison FGSEA
```{r}
fgsea_single_dot_data <- fgsea_res_pert %>%
filter(
metric == "signed_padj",
batch != "batch_1_riki"
) %>%
mutate(
contrast_batch = paste(contrast, recode_batch(batch)),
database = str_extract(pathway, "^(REACTOME|HALLMARK|PID|KEGG_MEDICUS|GOBP|GOMF|GOCC)"),
pathway = abbreviate_prefix(pathway)
) %>%
semi_join(
filter(
.,
padj < .01,
contrast %in% c("PTX", "GNE-495 + PTX vs PTX", "JNK-inh + PTX vs PTX")
),
by = "pathway"
) %>%
mutate(
padj_bin = cut(
padj, breaks = c(-Inf, 0.001, 0.01, 0.05, Inf), labels = c("<0.001", "<0.01", "<0.05", ">=0.05"),
ordered_result = TRUE
)
)
fgsea_single_dot_data_filtered <- fgsea_single_dot_data %>%
filter(
contrast == "PTX" |
str_detect(contrast, coll("vs PTX")),
!str_starts(database, coll("GO"))
) %>%
mutate(
contrast_batch = paste(contrast, recode_batch(batch))
) %>%
cluster_df(
row_var = pathway,
col_var = contrast_batch,
value_var = NES
)
est_perc <- quantile(fgsea_single_dot_data_filtered$NES, c(.05, .95), na.rm = TRUE)
abs_max_est_perc <- max(abs(est_perc))
fgsea_single_dot_plot <- ggplot(
fgsea_single_dot_data_filtered,
aes(
x = contrast,
y = pathway,
size = padj_bin,
color = NES
)
) +
geom_point() +
paletteer::scale_color_paletteer_c(
"ggthemes::Orange-Blue-White Diverging",
limits = c(-abs_max_est_perc, abs_max_est_perc),
oob = scales::squish
) +
scale_size_manual(
values = c(
"<0.001" = 5,
"<0.01" = 4,
"<0.05" = 3,
">=0.05" = 1
)
) +
facet_wrap(
~batch,
axes = "all",
axis.labels = "margins"
) +
labs(
x = "Contrast",
y = "Pathway",
color = "Normalized\nenrichment\nscore"
) +
theme_minimal() +
theme(
axis.text.x = element_text(angle = 45, hjust = 1)
)
fgsea_single_dot_plot
ggsave(
file.path("plots", "fgsea_single_perturbation_dotplot.pdf"),
fgsea_single_dot_plot,
width = 10,
height = 12,
device = cairo_pdf
)
```
## 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")
}
cluster_fun_binary <- function(mat, sample_in_col = TRUE) {
# Convert logical/character to binary if needed
if (is.logical(mat) || all(mat %in% c("TRUE", "FALSE", TRUE, FALSE))) {
mat <- matrix(as.logical(mat), nrow = nrow(mat), ncol = ncol(mat))
mat <- mat * 1 # Convert to 0/1
}
# Use binary distance (Jaccard is good for binary data)
# method = "binary" uses Jaccard distance for binary data
dist_mat <- dist(mat, method = "binary")
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,
!(contrast == "PTX" & batch == "batch_2"),
batch != "batch_1_riki"
) %>%
mutate(
database = str_extract(pathway, "^(REACTOME|HALLMARK|PID|KEGG_MEDICUS|GOBP|GOMF|GOCC)"),
pathway = abbreviate_prefix(pathway),
contrast_batch = paste(contrast, recode_batch(batch)),
signed_padj = -sign(NES) * log10(padj)
) %>%
filter(
!str_starts(database, coll("GO"))
) %>%
cluster_df(
pathway, contrast_batch, signed_padj, values_fill = 0
)
fgsea_res_pert_mat <- fgsea_res_pert_hm_data %>%
arrange(pathway, contrast_batch) %>%
select(pathway, contrast_batch, signed_padj) %>%
pivot_wider(names_from = contrast_batch, 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), ]
)
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_contrast_batch_order <- fgsea_res_pert_hm_data %>%
distinct(batch, contrast = factor(contrast, levels = fgsea_hm_conditions), contrast_batch) %>%
arrange(contrast, batch)
hm <- Heatmap(
fgsea_res_pert_mat[as.character(fgsea_res_pert_hm_pathways), hm_contrast_batch_order$contrast_batch],
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"
)
)
fgsea_picked_to_highlight <- c(
"METABOLISM_OF_LIPIDS",
"METABOLISM_OF_STEROIDS",
"CHOLESTEROL_BIOSYNTHESIS",
"ACTIVATION_OF_GENE_EXPRESSION_BY_SREBF_SREBP",
"REGULATION_OF_CHOLESTEROL_BIOSYNTHESIS_BY_SREBP_SREBF",
"REGULATION_OF_EXPRESSION_OF_SLITS_AND_ROBOS",
"SIGNALING_BY_ROBO_RECEPTORS",
"FATTY_ACID_METABOLISM",
"FATTY_ACYL_COA_BIOSYNTHESIS",
"PHOSPHOLIPID_METABOLISM",
"GLYCEROPHOSPHOLIPID_BIOSYNTHESIS",
"GAP_JUNCTION_ASSEMBLY",
"KINESINS",
"GAP_JUNCTION_TRAFFICKING_AND_REGULATION",
"ORGANELLE_BIOGENESIS_AND_MAINTENANCE",
"RHO_GTPASES_ACTIVATE_IQGAPS",
"COPI_INDEPENDENT_GOLGI_TO_ER_RETROGRADE_TRAFFIC",
"RHO_GTPASES_ACTIVATE_FORMINS",
"HDACS_DEACETYLATE_HISTONES",
"COPI_DEPENDENT_GOLGI_TO_ER_RETROGRADE_TRAFFIC",
"GOLGI_TO_ER_RETROGRADE_TRANSPORT",
"INTRA_GOLGI_AND_RETROGRADE_GOLGI_TO_ER_TRAFFIC",
"SIGNALING_BY_RHO_GTPASES_MIRO_GTPASES_AND_RHOBTB3",
"TRANSPORT_TO_THE_GOLGI_AND_SUBSEQUENT_MODIFICATION",
"RHO_GTPASE_EFFECTORS",
"APOPTOSIS",
"ER_TO_GOLGI_ANTEROGRADE_TRANSPORT",
"VESICLE_MEDIATED_TRANSPORT",
"MEMBRANE_TRAFFICKING"
)
mat <- fgsea_res_pert_mat[as.character(fgsea_res_pert_hm_pathways), hm_contrast_batch_order$contrast_batch]
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 = fgsea_picked_to_highlight %>%
{
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", "gene_set_analysis", "fgsea_single_perturbation_heatmap_nes_k8.pdf"),
draw(hm),
width = 6, height = 8
)
InteractiveComplexHeatmap::htShiny(draw(hm), save = here("plots/reactome_pert_heatmap_k8"))
hm_with_names <- 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 = TRUE,
row_names_gp = gpar(fontsize = 8),
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"
)
)
withr::with_pdf(
here("plots", "reactome_pert_heatmap_k8_with_names.pdf"),
draw(hm_with_names),
width = 9, height = 16
)
```
## 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_collection %in% c("H") |
gs_subcollection %in% c(
"CP:PID",
"CP:REACTOME",
"CP:KEGG_MEDICUS",
"GO:BP",
"GO:MF",
"GO:CC"
),
!str_detect(gs_name, coll("MEDICUS_PATHOGEN")),
!str_detect(gs_name, coll("MEDICUS_VARIANT"))
) %>%
mutate(
database = str_extract(gs_name, "^(REACTOME|HALLMARK|PID|KEGG_MEDICUS|GOBP|GOMF|GOCC)"),
pathway = abbreviate_prefix(gs_name)
)
msigdbr_of_interest_mat <- msigdbr_of_interest %>%
transmute(gs_name, gene_symbol, dummy = TRUE) %>%
distinct() %>%
pivot_wider(
names_from = gs_name,
values_from = dummy,
values_fill = FALSE
) %>%
column_to_rownames("gene_symbol") %>%
as.matrix()
msigdbr_of_interest_sim <- msigdbr_of_interest_mat %>%
proxy::dist(method = "Jaccard", by_rows = FALSE)
qsave(
msigdbr_of_interest_sim,
here("results", "msigdbr_of_interest_jaccard_sim.qs")
)
```
```{r}
fgsea_hm_conditions <- c(
"PTX", "GNE-495", "JNK-inh",
"GNE-495 + PTX", "JNK-inh + PTX",
"GNE-495 + PTX vs PTX", "JNK-inh + PTX vs PTX"
)
de_pert_hm_data <- de_long_selected %>%
filter(
contrast %in% fgsea_hm_conditions
) %>%
mutate(
signed_padj = -sign(logFC) * log10(FDR_min_zero),
contrast_batch = paste(contrast, recode_batch(batch))
)
fgsea_hm_order <- de_pert_hm_data %>%
distinct(
batch,
contrast = factor(contrast, levels = fgsea_hm_conditions),
contrast_batch
) %>%
arrange(contrast, batch)
gene_id_gene_symbol_map <- de_pert_hm_data %>%
distinct(gene_id, gene_name) %>%
{with(., set_names(gene_name, gene_id))}
plot_pathway_heatmap <- function(pathway, value_col = signed_padj) {
gene_sets <- msigdbr_of_interest %>%
filter(.data$pathway == .env$pathway) %>%
pull(ensembl_gene) %>%
unique()
de_mat <- de_pert_hm_data %>%
filter(gene_id %in% gene_sets) %>%
select(gene_id, contrast_batch, {{value_col}}) %>%
pivot_wider(names_from = contrast_batch, values_from = {{value_col}}, values_fill = 0) %>%
column_to_rownames("gene_id") %>%
as.matrix() %>% {
.[, fgsea_hm_order$contrast_batch]
}
# row_meta <- fgsea_res_pert %>%
# filter(
# pathway == !!pathway,
# metric == "signed_padj"
# ) %>%
# select(compound, gene_id = leadingEdge) %>%
# unchop(gene_id) %>%
# mutate(
# gene_id = factor(gene_id, levels = rownames(de_mat)),
# decoy = 1L
# ) %>%
# pivot_wider(
# names_from = compound, values_from = decoy,
# names_prefix = "leading_edge_",
# id_expand = TRUE, values_fill = 0L
# ) %>%
# column_to_rownames("gene_id") %>% {
# .[rownames(de_mat), ]
# }
mat_max_abs <- max(abs(quantile(de_mat, c(.025, .975), na.rm = TRUE)))
hm <- Heatmap(
de_mat,
name = "Signed -log10 FDR",
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,
column_split = fgsea_hm_order$contrast,
row_labels = gene_id_gene_symbol_map[rownames(de_mat)],
# right_annotation = HeatmapAnnotation(
# df = row_meta,
# which = "row"
# )
)
return(hm)
}
plot_pathway_heatmap("R_P38MAPK_EVENTS")
library(ComplexHeatmap)
de_pert_rescue_hm_data <- tibble(
pathway = c(
"METABOLISM_OF_LIPIDS",
"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) {
tibble(
gene_id = rownames(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}
bar_plot_colors <- c(
up = "#d73027",
down = "#4575b4"
)
fgsea_overrep_data <- fgsea_overrep %>%
mutate(
compound_batch = paste(compound, recode_batch(batch)),
database = str_extract(pathway, "^(REACTOME|HALLMARK|PID|KEGG_MEDICUS|GOBP|GOMF|GOCC)"),
pathway = abbreviate_prefix(pathway)
) %>%
mutate(
# sign = if_else(NES > 0, "pos", "neg"),
signed_padj = -if_else(perturbed == "up", 1, -1) * log10(padj),
pathway_direction = paste(perturbed, pathway, sep = "_")
)
ps <- fgsea_overrep_data %>%
group_nest(
perturbed_threshold, compound_batch, rescue_definition
) %>%
mutate(
p = pmap(
list(compound_batch, data, rescue_definition),
\(x, y, z) {
y %>%
filter(
padj < .01
) %>%
mutate(
across(pathway_direction, \(x) fct_reorder(x, signed_padj))
) %>%
ggplot(
aes(
x = signed_padj,
y = pathway_direction,
fill = perturbed
)
) +
geom_col() +
geom_text(
aes(
x = if_else(perturbed == "up", -.05, .05),
hjust = if_else(perturbed == "up", 1, 0),
label = pathway
),
size = 2
) +
scale_fill_manual(
values = alpha(bar_plot_colors, .7)
) +
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_batch, p, rescue_definition, perturbed_threshold, ...) {
n <- nrow(layer_data(p))
h <- pmin(49, 2 + .1 * n)
ggsave(
here("plots", "gene_set_analysis", paste0("fgsea_single_rescue_", rescue_definition, "_", perturbed_threshold, "_", compound_batch, "_bars.pdf")),
p, width = 8, height = h
)
}
)
```
## vs PTX
```{r}
fgsea_overrep_vs_ptx_pert <- fgsea_overrep %>%
filter(
rescue_definition != "perturbed"
) %>%
mutate(
compound_batch = paste(compound, recode_batch(batch)),
database = str_extract(pathway, "^(REACTOME|HALLMARK|PID|KEGG_MEDICUS|GOBP|GOMF|GOCC)"),
pathway = abbreviate_prefix(pathway)
) %>%
mutate(
# sign = if_else(NES > 0, "pos", "neg"),
signed_padj = -if_else(perturbed == "up", 1, -1) * log10(padj),
pathway_direction = paste(perturbed, pathway, sep = "_")
) %>%
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", "batch"),
suffix = c("_overrep", "_PTX"),
check = check_specs(
unmatched_keys_left = "info",
unmatched_keys_right = "warn",
duplicate_keys_right = "warn"
)
)
ps <- fgsea_overrep_vs_ptx_pert %>%
group_nest(
rescue_definition, perturbed_threshold
) %>%
mutate(
p = map(
data,
\(d) {
ggplot(
d,
aes(
x = signed_padj_PTX,
y = signed_padj_overrep
)
) +
geom_point(
aes(
text = pathway
)
) +
facet_grid(
vars(compound), vars(batch)
)
}
)
)
pwalk(
ps,
\(rescue_definition, perturbed_threshold, p, ...) {
ggsave(
here("plots", "gene_set_analysis", paste0("fgsea_rescue_overrep_vs_ptx_fgsea_", rescue_definition, "_", perturbed_threshold, ".pdf")),
p, width = 10, height = 8
)
}
)
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 %>%
mutate(
database = str_extract(pathway, "^(REACTOME|HALLMARK|PID|KEGG_MEDICUS|GOBP|GOMF|GOCC)"),
pathway = abbreviate_prefix(pathway)
) %>%
mutate(
# sign = if_else(NES > 0, "pos", "neg"),
signed_padj = -if_else(perturbed == "up", 1, -1) * log10(padj),
pathway_direction = paste(perturbed, pathway, sep = "_")
)
ps <- fgsea_intersection_overrep_data %>%
filter(perturbed %in% c("up", "down")) %>%
group_nest(
perturbed_threshold, batch, rescue_type, rescue_definition
) %>%
mutate(
p = pmap(
list(rescue_type, rescue_definition, data),
\(x, y, z) {
d <- z %>%
filter(
padj < .01
) %>%
mutate(
across(pathway_direction, \(x) fct_reorder(x, signed_padj))
)
ggplot(
d,
aes(
x = signed_padj,
y = pathway_direction,
fill = perturbed
)
) +
geom_col() +
geom_text(
aes(
x = if_else(perturbed == "up", -.05, .05),
hjust = if_else(perturbed == "up", 1, 0),
label = pathway
),
size = 2
) +
scale_fill_manual(
values = alpha(bar_plot_colors, .7)
) +
theme(
axis.text.y = element_blank(),
panel.grid.major.y = element_blank(),
) + {
if (!any(d$perturbed == "down"))
lims(x = c(-max(abs(d$signed_padj)), NA))
} +
labs(
x = "Signed adjusted p-value",
y = NULL,
title = paste0("Pathways enriched in genes rescued by ", x),
subtitle = paste("FDR < 0.01", y)
)
}
)
)
pwalk(
ps,
\(perturbed_threshold, batch, rescue_definition, rescue_type, p, ...) {
ggsave(
here("plots", "gene_set_analysis", paste0("fgsea_rescue_overrep_bars_", rescue_definition, "_", rescue_type, "_", batch, "_", perturbed_threshold, ".pdf")),
p, width = 10, height = 8
)
}
)
```
### Dot plot
```{r}
fgsea_intersection_dotplot_data <- fgsea_intersection_overrep_data %>%
mutate(
padj_bin = cut(
padj,
breaks = c(0, .01, .05, 1),
labels = c("FDR < 0.01", "0.01 < FDR < 0.05", "ns")
),
estimate = log2(foldEnrichment)
) %>%
filter(
rescue_type %in% c("GNE-495 only", "JNK-inh only"),
perturbed %in% c("up", "down")
) %>%
# inner_join(
# fgsea_intersection_overrep_data_single,
# by = c("rescue_definition", "rescue_type", "pathway_db")
# ) %>%
# mutate(
# directed_estimate = if_else(direction == "up", estimate, -estimate)
# ) %>%
group_nest(
perturbed_threshold, batch, rescue_definition
) %>%
mutate(
p = map(
data,
\(x) {
d <- group_by(
x,
pathway
) %>%
filter(
any(padj < .01)
) %>%
ungroup()
dmats <- d %>%
select(rescue_type, pathway, perturbed, estimate) %>%
arrange(pathway) %>%
split(.$perturbed) %>%
map(
\(y) select(y, -perturbed) %>%
pivot_wider(
names_from = rescue_type,
values_from = estimate
) %>%
column_to_rownames("pathway") %>%
as.matrix() %>%
dist()
)
dmat <- reduce(dmats, `+`)
browser()
clust <- hclust(dmat, method = "average") %>%
reorder(dmat, method = "OLO")
# browser()
ggplot(
d %>%
mutate(
pathway = factor(pathway, levels = clust$labels[clust$order])
),
aes(
x = rescue_type,
y = pathway,
color = estimate,
size = padj_bin
)
) +
geom_point(shape = 16) +
scale_color_viridis_c(trans = "pseudo_log") +
scale_size_manual(
values = c("FDR < 0.01" = 4, "0.01 < FDR < 0.05" = 2, "ns" = 1)
) +
# scale_color_distiller(palette = "RdBu") +
facet_wrap(~perturbed)
}
)
)
```
Check what genes are in JNK-inh rescued mitosis related pathways
```{r}
x <- fgsea_intersection_overrep_data %>%
arrange(p.value_fgsea) %>%
filter(
rescue_definition == "significant_opposite",
rescue_type == "JNK-inh only"
)
x %>%
filter(
pathway == "METALLOPROTEASE_DUBS"
) %>%
chuck("overlapGenes", 1) %>%
{gene_id_gene_symbol_map[.]}
gene_id_gene_symbol_map[x$overlapGenes[[1]]]
```
```{r}
flag_similar_adjacent_leaves <- function(clust, sim_matrix, threshold) {
# Get the dendrogram representation
dend <- as.dendrogram(clust)
# Function to find adjacent leaves and check their similarity
find_adjacent_pairs <- function(node, parent_leaves = NULL) {
if (is.leaf(node)) {
# For leaf nodes, return the label
return(list(leaves = labels(node), pairs = data.frame()))
} else {
# For internal nodes, process left and right children
left_result <- find_adjacent_pairs(node[[1]], parent_leaves)
right_result <- find_adjacent_pairs(node[[2]], left_result$leaves)
# Combine leaves from both children
all_leaves <- c(left_result$leaves, right_result$leaves)
# Identify pairs of leaves that are adjacent across the left-right boundary
# (the rightmost leaf of left subtree and leftmost leaf of right subtree)
boundary_pairs <- data.frame(
leaf1 = tail(left_result$leaves, 1),
leaf2 = head(right_result$leaves, 1),
stringsAsFactors = FALSE
)
# Calculate similarity for these boundary pairs
if (nrow(boundary_pairs) > 0) {
boundary_pairs$similarity <- sim_matrix[boundary_pairs$leaf1, boundary_pairs$leaf2]
# Flag pairs below threshold
boundary_pairs <- boundary_pairs[boundary_pairs$similarity <= threshold, ]
}
# Combine all pairs found
all_pairs <- rbind(left_result$pairs, right_result$pairs, boundary_pairs)
return(list(leaves = all_leaves, pairs = all_pairs))
}
}
# Start the recursive search
result <- find_adjacent_pairs(dend)
return(result$pairs)
}
```
```{r}
msigdbr_of_interest_sim_trans <- msigdbr_of_interest_sim
attr(msigdbr_of_interest_sim_trans, "Labels") <- distinct(
msigdbr_of_interest, gs_name, database, pathway
) %>%
transmute(
pathway_db = paste(str_sub(database, 1, 1), pathway, sep = "_"),
gs_name
) %>%
slice(match(labels(msigdbr_of_interest_sim_trans), gs_name)) %>%
pull(pathway_db)
msigdbr_of_interest_sim_trans <- as.matrix(msigdbr_of_interest_sim_trans)
fgsea_intersection_overrep_data_single <- fgsea_intersection_overrep_data %>%
filter(
perturbed %in% c("up", "down")
) %>%
mutate(
# Replacing -Inf with -1 for log2 fold enrichment
estimate = log2(if_else(foldEnrichment == 0, .5, foldEnrichment))
) %>%
group_by(
perturbed_threshold, batch, rescue_definition, rescue_type, database, pathway
) %>%
summarize(
direction = case_when(
all(padj >= .01) ~ "none",
padj[1] < .2 * padj[2] | padj[2] < .2 * padj[1] ~ perturbed[which.min(padj)],
TRUE ~ "both"
),
estimate = estimate[order(padj)[1]],
padj = min(padj),
# estimate = estimate[which.min(padj)],
.groups = "drop"
)
fgsea_blacklist <- c(
"HCMV EARLY EVENTS",
"HCMV INFECTION",
"INFECTIOUS DISEASE"
)
fgsea_intersection_dotplot_data_single <- fgsea_intersection_overrep_data_single %>%
mutate(
padj_bin = cut(
padj,
breaks = c(0, .001, .01, 1),
labels = c("FDR < 0.001", "FDR < 0.01", "ns")
),
# across(
# pathway_db,
# \(x) str_remove_all(
# x, coll("_")
# ) %>%
# str_sub(start = 3)
# )
) %>%
filter(
!database %in% c("KEGG", "GOBP", "GOMF", "GOCC"),
rescue_type %in% c("GNE-495 only", "JNK-inh only"),
) %>%
mutate(
directed_estimate = if_else(direction == "up", estimate, -estimate),
rescute_type_batch = paste(rescue_type, recode_batch(batch), sep = "_")
) %>%
group_nest(
perturbed_threshold, rescue_definition
) %>%
crossing(
p_threshold = c(0.005, 0.001, .01)
) %>%
mutate(
res = map2(
data, p_threshold,
\(x, y) {
# browser()
d <- group_by(
x,
pathway
) %>%
filter(
any(padj < y)
) %>%
ungroup()
dwide <- d %>%
select(rescute_type_batch, pathway, directed_estimate) %>%
arrange(pathway) %>%
pivot_wider(
names_from = rescute_type_batch,
values_from = directed_estimate
)
dmat <- dwide %>%
column_to_rownames("pathway") %>%
as.matrix() %>%
dist()
# browser()
clust <- hclust(dmat, method = "average") %>%
reorder(dmat, method = "OLO")
row_labels <- clust$labels[clust$order] %>% {
set_names(
str_replace_all(., coll("_"), " ") %>%
str_sub(start = 2),
.
)
}
# sim_hm <- msigdbr_of_interest_sim_trans[
# rev(clust$labels[clust$order]), rev(clust$labels[clust$order])
# ] %>%
# pheatmap(
# cluster_rows = FALSE,
# cluster_cols = FALSE
# )
est_perc <- quantile(d$directed_estimate, c(.05, .95))
abs_max_est_perc <- max(abs(est_perc))
p <- ggplot(
d %>%
mutate(
pathway = factor(pathway, levels = clust$labels[clust$order])
),
aes(
x = rescute_type_batch,
y = pathway,
color = directed_estimate,
size = padj_bin
)
) +
geom_point(shape = 16) +
# scale_color_viridis_c(trans = "pseudo_log") +
# scale_color_distiller(
# limits = c(-abs_max_est_perc, abs_max_est_perc),
# palette = "RdYlBu", trans = "pseudo_log",
# oob = scales::squish
# ) +
paletteer::scale_color_paletteer_c(
palette = "pals::ocean.balance",
limits = c(-abs_max_est_perc, abs_max_est_perc),
oob = scales::squish
) +
scale_size_manual(
values = rev(c("FDR < 0.001" = 6, "FDR < 0.01" = 4, "ns" = 2))
) +
scale_x_discrete(position = "top") +
scale_y_discrete(
breaks = names(row_labels),
labels = unname(row_labels)
) +
theme(
axis.text.x = element_text(angle = 45, hjust = 0)
) +
labs(
x = "Rescue by",
y = NULL,
color = "Rescue of\ngenes",
size = NULL
)
list(
p = p
# sim_hm = sim_hm
) %>%
map(list)
}
)
) %>%
select(-data) %>%
unnest_wider(
res
)
flag_similar_adjacent_leaves(clust, msigdbr_of_interest_sim_trans, .6)
pwalk(
fgsea_intersection_dotplot_data_single,
\(p, rescue_definition, p_threshold, perturbed_threshold, ...) {
ggsave(
here("plots", "gene_set_analysis", paste0("fgsea_rescue_overrep_dotplot_", rescue_definition, "_", perturbed_threshold, "_", p_threshold, ".pdf")),
p[[1]], width = 9.8, height = 9
)
# withr::with_pdf(
# file.path("plots", "overrep", paste0("overrep_intersection_dotplot_", rescue_definition, "_", p_threshold, "_similarity.pdf")),
# draw(sim_hm[[1]]),
# width = 9.8, height = 9
# )
}
)
fgsea_intersection_dotplot_data_single <- fgsea_intersection_overrep_data_single %>%
mutate(
padj_bin = cut(
padj,
breaks = c(0, .01, .05, 1),
labels = c("FDR < 0.01", "FDR < 0.05", "ns")
),
# across(
# pathway_db,
# \(x) str_remove_all(
# x, coll("_")
# ) %>%
# str_sub(start = 3)
# )
) %>%
filter(
!database %in% c("KEGG", "GOBP", "GOMF", "GOCC"),
rescue_type %in% c("GNE-495 only", "JNK-inh only"),,
batch != "batch_1_riki"
) %>%
mutate(
directed_estimate = if_else(direction == "up", estimate, -estimate),
rescute_type_batch = paste(rescue_type, recode_batch(batch), sep = "_")
) %>%
group_nest(
perturbed_threshold, rescue_definition
) %>%
crossing(
p_threshold = c(0.005, 0.001, .01)
) %>%
mutate(
res = map2(
data, p_threshold,
possibly(\(x, y) {
# browser()
d <- group_by(
x,
pathway
) %>%
filter(
any(padj < y)
) %>%
ungroup()
dwide <- d %>%
select(rescute_type_batch, pathway, directed_estimate) %>%
arrange(pathway) %>%
pivot_wider(
names_from = rescute_type_batch,
values_from = directed_estimate
)
dmat <- dwide %>%
column_to_rownames("pathway") %>%
as.matrix() %>%
dist()
# browser()
clust <- hclust(dmat, method = "average") %>%
reorder(dmat, method = "OLO")
row_labels <- clust$labels[clust$order] %>% {
set_names(
str_replace_all(., coll("_"), " ") %>%
str_sub(start = 2),
.
)
}
# sim_hm <- msigdbr_of_interest_sim_trans[
# rev(clust$labels[clust$order]), rev(clust$labels[clust$order])
# ] %>%
# pheatmap(
# cluster_rows = FALSE,
# cluster_cols = FALSE
# )
est_perc <- quantile(d$directed_estimate, c(.05, .95))
abs_max_est_perc <- max(abs(est_perc))
p <- ggplot(
d %>%
mutate(
pathway = factor(pathway, levels = clust$labels[clust$order])
),
aes(
x = rescute_type_batch,
y = pathway,
color = directed_estimate,
size = padj_bin
)
) +
geom_point(shape = 16) +
# scale_color_viridis_c(trans = "pseudo_log") +
# scale_color_distiller(
# limits = c(-abs_max_est_perc, abs_max_est_perc),
# palette = "RdYlBu", trans = "pseudo_log",
# oob = scales::squish
# ) +
paletteer::scale_color_paletteer_c(
palette = "pals::ocean.balance",
limits = c(-abs_max_est_perc, abs_max_est_perc),
oob = scales::squish
) +
scale_size_manual(
values = rev(c("FDR < 0.01" = 6, "FDR < 0.05" = 4, "ns" = 2))
) +
scale_x_discrete(position = "top") +
scale_y_discrete(
breaks = names(row_labels),
labels = unname(row_labels)
) +
theme(
axis.text.x = element_text(angle = 45, hjust = 0)
) +
labs(
x = "Rescue by",
y = NULL,
color = "Rescue of\ngenes",
size = NULL
)
list(
p = p
# sim_hm = sim_hm
) %>%
map(list)
})
)
) %>%
select(-data) %>%
unnest_wider(
res
)
pwalk(
fgsea_intersection_dotplot_data_single %>%
filter(map_lgl(p, \(x) !is.null(x))),
\(p, rescue_definition, p_threshold, perturbed_threshold, ...) {
ggsave(
here("plots", "gene_set_analysis", paste0("fgsea_rescue_overrep_dotplot_batch1_2_", rescue_definition, "_", perturbed_threshold, "_", p_threshold, ".pdf")),
p[[1]], width = 9.8, height = 9
)
# withr::with_pdf(
# file.path("plots", "overrep", paste0("overrep_intersection_dotplot_", rescue_definition, "_", p_threshold, "_similarity.pdf")),
# draw(sim_hm[[1]]),
# width = 9.8, height = 9
# )
}
)
```
### 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"))
)
}
)
fgsea_intersection_overrep_gne_vs_jnk <- fgsea_intersection_overrep_data %>%
filter(
rescue_definition == "significant_opposite",
perturbed %in% c("up", "down"),
rescue_type %in% c("GNE-495 only", "JNK-inh only")
) %>%
select(
perturbed, rescue_type, database, pathway, pathway_db, pathway_db_direction,
signed_padj, estimate
) %>%
pivot_wider(
names_from = rescue_type,
values_from = c(signed_padj, estimate)
)
p <- fgsea_intersection_overrep_gne_vs_jnk %>%
ggplot(
aes(
x = `signed_padj_JNK-inh only`,
y = `signed_padj_GNE-495 only`,
text = pathway_db
)
) +
geom_point()
p
plotly::ggplotly(p)
```
### 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
) %>%
distinct(gene_symbol, ensembl_gene) %>%
semi_join(
de_pert_hm_data,
by = c("ensembl_gene" = "gene_id")
)
),
mat = map(
genes,
\(x) de_pert_hm_data %>%
semi_join(
x,
by = c("gene_id" = "ensembl_gene")
) %>%
select(gene_id, contrast, logFC) %>%
pivot_wider(names_from = contrast, values_from = logFC, values_fill = 0) %>%
column_to_rownames("gene_id") %>%
as.matrix() %>% {
.[x$ensembl_gene, fgsea_hm_conditions]
}
),
row_meta = map2(
pathway, genes,
\(x, y) {
power_left_join(
y,
fgsea_res_rescue_data %>%
filter(pathway == x) %>%
select(compound, gene_id = leadingEdge) %>%
unchop(gene_id) %>%
mutate(decoy = 1L) %>%
pivot_wider(
names_from = compound, values_from = decoy,
names_prefix = "leading_edge_",
values_fill = 0L
),
by = c("ensembl_gene" = "gene_id"),
check = check_specs(
unmatched_keys_right = "warn",
duplicate_keys_right = "warn",
duplicate_keys_left = "warn"
)
) %>%
mutate(
across(
starts_with("leading_edge_"),
\(x) if_else(is.na(x), 0L, x)
)
)
}
),
hm = pmap(
list(pathway, mat, row_meta, genes),
\(x, y, z, g) {
# browser()
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,
row_labels = g$gene_symbol,
right_annotation = HeatmapAnnotation(
df = z %>%
select(ensembl_gene, starts_with("leading_edge_")) %>%
column_to_rownames("ensembl_gene") %>%
as.matrix(),
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_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}
```
### INDRA plots
```{r}
idnra_res_plot_data <- indra_res %>%
filter(
rescue_type %in% c(
"GNE-495 only", "JNK-inh only"
),
str_starts(query_type, fixed("indra"))
) %>%
mutate(
signed_q = if_else(
perturbed == "down",
1, -1
) * log10(q)
) %>%
separate_wider_delim(
curie,
delim = ":",
names = c("database", "id")
) %>%
filter(
database %in% c("fplx", "hgnc"),
q < .01
) %>%
mutate(
sign = if_else(perturbed == "up", "pos", "neg")
)
ps <- idnra_res_plot_data %>%
group_nest(
rescue_type, query_type
) %>%
mutate(
p = map2(
rescue_type, data,
\(x, y) {
y %>%
mutate(
name = fct_reorder(name, signed_q)
) %>%
ggplot(
aes(
x = signed_q,
y = name,
fill = sign
)
) +
geom_col() +
geom_text(
aes(
x = if_else(sign == "pos", -.05, .05),
hjust = if_else(sign == "pos", 1, 0),
label = name
),
size = 2
) +
scale_fill_discrete(guide = "none", direction = -1) +
theme(
axis.text.y = element_blank(),
panel.grid.major.y = element_blank(),
) +
labs(
x = "Signed log10 p-value",
y = NULL,
title = paste0(query_type, " analysis of genes rescued by ", x),
subtitle = "FDR < 0.01"
)
}
)
)
pwalk(
ps,
\(rescue_type, query_type, p, ...) {
ggsave(
file.path("plots", paste0(query_type, "_", rescue_type, "_bars.pdf")),
p, width = 8, height = 8
)
}
)
```
## Check if p38 target TFs are differentially expressed
```{r}
atf6_targets <- clipr::read_clip()
de_pert_hm_data_atf6 <- de_pert_hm_data %>%
filter(
gene_name %in% atf6_targets,
gene_id %in% {
de_pert_hm_data %>%
filter(
FDR < .05
) %>%
pull(gene_id)
}
)
de_pert_hm_data_atf6_mat <- de_pert_hm_data_atf6 %>%
dplyr::select(
gene_id, contrast, logFC
) %>%
pivot_wider(names_from = contrast, values_from = logFC, values_fill = 0) %>%
column_to_rownames("gene_id") %>%
as.matrix() %>%
t() %>%
scale() %>%
t() %>% {
.[, fgsea_hm_conditions]
}
hm <- Heatmap(
de_pert_hm_data_atf6_mat,
cluster_rows = cluster_fun_eucl,
cluster_columns = FALSE,
show_row_names = TRUE
)
hm
pca_fun <- function(mat, meta) {
pca <- prcomp(t(mat), center = TRUE, scale = FALSE)
pca_x <- as_tibble(pca$x) %>%
bind_cols(
meta
)
list(
x = pca_x,
sd = broom::tidy(pca, matrix = "eigenvalues"),
rotation = broom::tidy(pca, matrix = "rotation")
)
}
de_pert_hm_data_atf6_pca <- pca_fun(
de_pert_hm_data_atf6_mat,
tibble(contrast = fgsea_hm_conditions)
)
p <- ggplot(
de_pert_hm_data_atf6_pca$x,
aes(
x = PC1, y = PC3, color = contrast
)
) +
geom_point() +
labs(
x = paste0("PC1 (", round(de_pert_hm_data_atf6_pca$sd$percent[1] * 100, 2), "%)"),
y = paste0("PC3 (", round(de_pert_hm_data_atf6_pca$sd$percent[3] * 100, 2), "%)")
) +
theme_minimal()
hm <- Heatmap(
de_pert_hm_data_atf6_pca$x %>%
column_to_rownames("contrast") %>%
# scale() %>%
t(),
cluster_rows = cluster_fun_eucl,
cluster_columns = FALSE,
show_row_names = TRUE
)
hm
```
VX-702
DORAMAPIMOD
LOSMAPIMOD
## Visualize overlap between rescue type lists
```{r}
ptx_rescue_overlap <- ptx_rescue_types %>%
filter(
rescue_type %in% c("GNE-495 only", "JNK-inh only", "both")
) %>%
mutate(
# rescue_type_direction = paste(
# rescue_type, perturbed
# ),
dummy = TRUE,
across(batch, recode_batch)
) %>%
select(-c(`GNE-495`, `JNK-inh`, perturbed)) %>%
pivot_wider(
names_from = c(rescue_type, batch),
values_from = dummy,
values_fill = list(dummy = FALSE)
)
ptx_rescue_overlap_hm <- ptx_rescue_overlap %>%
group_nest(
perturbed_threshold, rescue_definition
) %>%
mutate(
hm = map(
data,
\(x) {
mat <- x %>%
# mutate(across(-gene_id, \(y) if_else(y, "TRUE", "FALSE"))) %>%
column_to_rownames("gene_id") %>%
as.matrix()
# browser()
Heatmap(
mat * 1L,
name = "Rescue",
col = c(`0` = "white", `1` = "blue"),
# cluster_rows = cluster_fun_binary,
# cluster_columns = cluster_fun_binary,
clustering_distance_rows = "binary",
clustering_distance_columns = "binary",
show_row_names = FALSE
)
}
)
)
pwalk(
ptx_rescue_overlap_hm,
\(perturbed_threshold, rescue_definition, hm, ...) {
withr::with_pdf(
file.path("plots", paste0("ptx_rescue_gene_set_overlap_hm_", rescue_definition, "_", perturbed_threshold, ".pdf")),
draw(hm),
width = 6, height = 10
)
}
)
```
## Dot plot of overrepresentation both doses
```{r}
fgsea_overrep_both_doses %>%
filter(padj < .05) %>%
count(name, direction)
fgsea_overrep_single_treatments %>%
filter(padj < .05) %>%
count(contrast, direction, significant_threshold, batch) %>%
print(n = Inf)
```
name direction n
<chr> <chr> <int>
1 GNE-495 + PTX vs PTX both 1350
2 GNE-495 + PTX vs PTX down 747
3 GNE-495 + PTX vs PTX up 762
4 JNK-inh + PTX vs PTX both 2955
5 JNK-inh + PTX vs PTX down 1923
6 JNK-inh + PTX vs PTX up 1404
```{r}
combine_updown <- function(data, direction_col, fold_change_col, p_col, group_cols, p_threshold = .05) {
group_syms <- syms(group_cols)
combine_impl <- function(d, g) {
direction <- if (all(d[[p_col]] < p_threshold)) {
"both"
} else if (all(d[[p_col]] >= p_threshold)) {
"none"
} else {
d[[direction_col]][which.min(d[[p_col]])]
}
min_p_idx <- which.min(d[[p_col]])
fold_change <- if (direction == "both") {
NA_real_
} else {
d[[fold_change_col]][min_p_idx]
}
p_value <- d[[p_col]][min_p_idx]
tibble(
direction = direction,
fold_change = fold_change,
p_value = p_value
)
}
data %>%
group_by(!!!group_syms) %>%
group_modify(combine_impl) %>%
ungroup()
}
fgsea_overrep_ptx_updown_combined <- fgsea_overrep_single_treatments %>%
filter(
contrast == "PTX",
direction %in% c("up", "down"),
significant_threshold == .05,
batch %in% c("batch_1", "batch_2")
) %>%
mutate(
log2_foldEnrichment = log2(if_else(foldEnrichment == 0, .5, foldEnrichment))
) %>%
combine_updown(
direction_col = "direction",
fold_change_col = "log2_foldEnrichment",
p_col = "padj",
group_cols = c("contrast", "pathway", "batch"),
p_threshold = .01
)
fgsea_overrep_both_doses_updown_combined <- fgsea_overrep_both_doses %>%
filter(
direction %in% c("up", "down")
) %>%
mutate(
# Replacing -Inf with -1 for log2 fold enrichment
estimate = log2(if_else(foldEnrichment == 0, .5, foldEnrichment))
) %>%
combine_updown(
direction_col = "direction",
fold_change_col = "log2_foldEnrichment",
p_col = "padj",
group_cols = c("name", "pathway"),
p_threshold = .01
)
fgsea_overrep_both_doses_updown_combined %>%
filter(is.na(p_value))
```
Can't use combination strategy so well for PTX because many gene sets are enriched
in both up and down directions. Using "both"
```{r}
fgsea_overrep_ptx <- fgsea_overrep_single_treatments %>%
filter(
contrast == "PTX",
direction == "both",
significant_threshold == .05,
batch == "batch_1"
)
fgsea_overrep_both_doses_filtered <- fgsea_overrep_both_doses %>%
filter(
direction == "both"
)
fgsea_overrep_both_doses_and_ptx <- bind_rows(
fgsea_overrep_ptx,
fgsea_overrep_both_doses_filtered %>%
rename(contrast = name)
) %>%
mutate(
# Replacing -Inf with -1 for log2 fold enrichment
log2_foldEnrichment = log2(if_else(foldEnrichment == 0, .5, foldEnrichment))
)
fgsea_overrep_both_doses_updown_combined_dot_plots <- fgsea_overrep_both_doses_and_ptx %>%
mutate(
padj_bin = cut(
padj,
breaks = c(0, .001, .01, .05, 1),
labels = c("<0.001", "<0.01", "<0.05", "ns")
),
database = str_extract(pathway, "^(REACTOME|HALLMARK|PID|KEGG_MEDICUS|GOBP|GOMF|GOCC)"),
pathway = abbreviate_prefix(pathway)
) %>%
filter(
database %in% c("REACTOME", "HALLMARK", "PID", "KEGG_MEDICUS")
) %>%
crossing(
p_threshold = c(0.001, .01, .05)
) %>%
group_nest(p_threshold) %>%
mutate(
res = map2(
data, p_threshold,
\(x, y) {
# browser()
# d <- semi_join(
# x,
# group_by(x, contrast) %>%
# filter(padj < y),
# by = "pathway"
# )
# d <- semi_join(
# x,
# group_by(x, contrast) %>%
# arrange(padj) %>%
# slice_head(n = 20),
# by = "pathway"
# )
d <- semi_join(
x,
group_by(x, pathway) %>%
filter(padj < y) %>%
summarize(groups = paste(unique(contrast), collapse = "_"), p_mean = 1 / sum(1 / padj, na.rm = TRUE)) %>%
group_by(groups) %>%
arrange(p_mean) %>%
slice_head(n = 20),
by = "pathway"
)
dwide <- d %>%
arrange(pathway) %>%
pivot_wider(
id_cols = pathway,
names_from = contrast,
values_from = log2_foldEnrichment
)
dmat <- dwide %>%
column_to_rownames("pathway") %>%
as.matrix() %>%
dist()
# browser()
clust <- hclust(dmat, method = "average") %>%
reorder(dmat, method = "OLO")
row_labels <- clust$labels[clust$order] %>% {
set_names(
str_replace_all(., coll("_"), " ") %>%
str_sub(start = 2),
.
)
}
# sim_hm <- msigdbr_of_interest_sim_trans[
# rev(clust$labels[clust$order]), rev(clust$labels[clust$order])
# ] %>%
# pheatmap(
# cluster_rows = FALSE,
# cluster_cols = FALSE
# )
est_perc <- quantile(d$log2_foldEnrichment, c(.05, .95))
abs_max_est_perc <- max(abs(est_perc))
p <- ggplot(
d %>%
mutate(
pathway = factor(pathway, levels = clust$labels[clust$order])
),
aes(
x = contrast,
y = pathway,
color = log2_foldEnrichment,
size = padj_bin
)
) +
geom_point(shape = 16) +
# scale_color_viridis_c(trans = "pseudo_log") +
# scale_color_distiller(
# limits = c(-abs_max_est_perc, abs_max_est_perc),
# palette = "RdYlBu", trans = "pseudo_log",
# oob = scales::squish
# ) +
paletteer::scale_color_paletteer_c(
palette = "pals::ocean.balance",
limits = c(-abs_max_est_perc, abs_max_est_perc),
oob = scales::squish
) +
scale_size_manual(
values = rev(c("<0.001" = 6, "<0.01" = 4, "<0.05" = 3, "ns" = 2))
) +
scale_x_discrete(position = "top") +
scale_y_discrete(
breaks = names(row_labels),
labels = unname(row_labels)
) +
theme(
axis.text.x = element_text(angle = 45, hjust = 0)
) +
labs(
x = NULL,
y = NULL,
size = NULL
)
list(
p = p,
clust_data = d
# sim_hm = sim_hm
) %>%
map(list)
}
)
) %>%
select(-data) %>%
unnest_wider(
res
)
flag_similar_adjacent_leaves(clust, msigdbr_of_interest_sim_trans, .6)
pwalk(
fgsea_intersection_dotplot_data_single,
\(p, rescue_definition, p_threshold, perturbed_threshold, ...) {
ggsave(
here("plots", "gene_set_analysis", paste0("fgsea_rescue_overrep_dotplot_", rescue_definition, "_", perturbed_threshold, "_", p_threshold, ".pdf")),
p[[1]], width = 9.8, height = 9
)
# withr::with_pdf(
# file.path("plots", "overrep", paste0("overrep_intersection_dotplot_", rescue_definition, "_", p_threshold, "_similarity.pdf")),
# draw(sim_hm[[1]]),
# width = 9.8, height = 9
# )
}
)
```
### Up down separately
```{r}
fgsea_overrep_ptx <- fgsea_overrep_single_treatments %>%
filter(
contrast == "PTX",
direction != "both",
significant_threshold == .05,
batch == "batch_1"
)
fgsea_overrep_both_doses_filtered <- fgsea_overrep_both_doses %>%
filter(
direction != "both"
)
fgsea_overrep_both_doses_and_ptx <- bind_rows(
fgsea_overrep_ptx,
fgsea_overrep_both_doses_filtered %>%
rename(contrast = name)
) %>%
mutate(
# Replacing -Inf with -1 for log2 fold enrichment
log2_foldEnrichment = log2(if_else(foldEnrichment == 0, .5, foldEnrichment)),
directed_enrichment = pmax(
log2_foldEnrichment,
0
) * if_else(direction == "up", 1, -1)
)
fgsea_overrep_both_doses_updown_combined_dot_plots <- fgsea_overrep_both_doses_and_ptx %>%
mutate(
padj_bin = cut(
padj,
breaks = c(0, .001, .01, .05, 1),
labels = c("<0.001", "<0.01", "<0.05", "ns")
),
database = str_extract(pathway, "^(REACTOME|HALLMARK|PID|KEGG_MEDICUS|GOBP|GOMF|GOCC)"),
pathway = abbreviate_prefix(pathway)
) %>%
filter(
database %in% c("REACTOME", "HALLMARK", "PID", "KEGG_MEDICUS")
) %>%
crossing(
p_threshold = c(0.001, .01, .05)
) %>%
group_nest(p_threshold) %>%
mutate(
res = map2(
data, p_threshold,
\(x, y) {
# browser()
# d <- semi_join(
# x,
# group_by(x, contrast) %>%
# filter(padj < y),
# by = "pathway"
# )
# d <- semi_join(
# x,
# group_by(x, contrast) %>%
# arrange(padj) %>%
# slice_head(n = 20),
# by = "pathway"
# )
d <- semi_join(
x,
group_by(x, pathway) %>%
filter(padj < y, is_main) %>%
summarize(groups = paste(sort(unique(contrast)), collapse = "_"), p_mean = 1 / sum(1 / padj, na.rm = TRUE)) %>%
group_by(groups) %>%
arrange(p_mean) %>%
slice_head(n = 10) %>%
bind_rows(tibble(pathway = "P_P38_MK2_PATHWAY")),
by = "pathway"
)
dwide <- d %>%
arrange(pathway) %>%
pivot_wider(
id_cols = pathway,
names_from = c(direction, contrast),
values_from = directed_enrichment
)
dmat <- dwide %>%
column_to_rownames("pathway") %>%
as.matrix() %>%
dist()
# browser()
clust <- hclust(dmat, method = "average") %>%
reorder(dmat, method = "OLO")
row_labels <- clust$labels[clust$order] %>% {
set_names(
str_replace_all(., coll("_"), " ") %>%
str_sub(start = 2),
.
)
}
# sim_hm <- msigdbr_of_interest_sim_trans[
# rev(clust$labels[clust$order]), rev(clust$labels[clust$order])
# ] %>%
# pheatmap(
# cluster_rows = FALSE,
# cluster_cols = FALSE
# )
est_perc <- quantile(d$directed_enrichment, c(.05, .95))
abs_max_est_perc <- max(abs(est_perc))
p <- ggplot(
d %>%
mutate(
pathway = factor(pathway, levels = clust$labels[clust$order])
),
aes(
x = contrast,
y = pathway,
color = directed_enrichment,
size = padj_bin,
shape = direction
)
) +
geom_text(
aes(
label = if_else(direction == "up", "◖", "◗"),
hjust = if_else(direction == "up", .5, .5)
),
vjust = .5
) +
guides(
size = guide_legend(override.aes = list(shape = "◖"))
) +
paletteer::scale_color_paletteer_c(
palette = "ggthemes::Orange-Blue-White Diverging",
limits = c(-abs_max_est_perc, abs_max_est_perc),
oob = scales::squish,
direction = -1
) +
scale_size_manual(
values = rev(c("<0.001" = 8, "<0.01" = 6, "<0.05" = 4, "ns" = 2))
) +
scale_x_discrete(position = "top") +
scale_y_discrete(
breaks = names(row_labels),
labels = unname(row_labels)
) +
theme(
axis.text.x = element_text(angle = 45, hjust = 0),
panel.grid.major.x = element_blank()
) +
labs(
x = NULL,
y = NULL,
size = "FDR",
color = "log2 fold\nenrichment"
)
list(
p = p,
clust_data = d
# sim_hm = sim_hm
) %>%
map(list)
}
)
) %>%
select(-data) %>%
unnest_wider(
res
)
pwalk(
fgsea_overrep_both_doses_updown_combined_dot_plots,
\(p_threshold, p, ...) {
ggsave(
file.path(
"plots", "gene_set_analysis", paste0("fgsea_overrep_both_doses_updown_dotplot_separate_", p_threshold, ".pdf")
), p[[1]], width = 6, height = 12, dev = cairo_pdf
)
}
)
```
```{r}
fgsea_overrep_both_doses_and_ptx %>%
filter(
str_detect(pathway, coll("p38", ignore_case = TRUE))
) %>%
arrange(pval) %>%
View()
fgsea_res_pert %>%
filter(
str_detect(pathway, coll("p38", ignore_case = TRUE)),
metric == "signed_padj"
) %>%
View()
```
## Enrichment bar graphs on DE gene clusters
```{r}
fgsea_overrep_clusters
```
TODO:
1. Heatmap of logFC all DE genes, including JNK+GNE double/triple treatments
2. Dot plot of PTX, GNE vs DMSO, JNK vs DMOS, double treatments vs PTX (FGSEA enrichments)
3. Scatter plots of FGSEA enrichments, highlighting interesting ones, labelling quadrants
4. PCAs and volcano plo

gene_set_analysis_plotting.Rmd at commit 3329343, no license · at the source

OverviewAuthors: Veselina Petrova1,2, Caitlin E Mills3,4, Clemens Hug3, Aysel Cetinkaya-Fisgin5, Jennifer Splaine6, Sepideh Fouladzadeh7, Sara Hakim1,2, Rasheen Powell1,2, Shannon Zhen1, Mirra Chung3,4, Gary A Bradshaw3, Tao Deng8, Ilyas Singec8, Qing Wang9,10, Riki Kawaguchi9,10, Harathi Jonnagaddala11, Lee B Barrett1,2, Jennifer A Smith6, Marian Kalocsay11, Benjamin M Gyori7,12, Ahmet Hoke5, Peter K Sorger3,4, Clifford J Woolf1,2
ORCID iDs: Rasheen Powell, Shannon Zhen, Ilyas Singec, Marian Kalocsay, Clifford J Woolf
F.M. Kirby Neurobiology Center, Program in Neurobiology, Boston Children’s Hospital, Boston, MA, USA
Department of Neurobiology, Harvard Medical School, Boston, MA, USA
Laboratory of Systems Pharmacology, Harvard Program in Therapeutic Science, Harvard Medical School, Boston, MA, USA
Department of Systems Biology, Harvard Medical School, Boston, MA, USA
Department of Neurology, Neuromuscular Division, Johns Hopkins School of Medicine, Baltimore, MD, USA
ICCB-Longwood Screening Facility, Harvard Medical School, 250 Longwood Avenue, Boston, MA, USA
Department of Bioengineering, Northeastern University, Boston, MA, USA
National Center for Advancing Translational Sciences (NCATS), Division of Preclinical Innovation, Stem Cell Translation Laboratory (SCTL), National Institutes of Health (NIH), Rockville, MD, USA
Department of Neurology, David Geffen School of Medicine, University of California, Los Angeles, Los Angeles, CA, USA
Center for Neurobehavioral Genetics, Semel Institute for Neuroscience and Human Behavior, University of California, Los Angeles, Los Angeles, CA, USA
Department of Experimental Radiation Oncology, The University of Texas MD Anderson Cancer Center, Houston, TX, USA
Khoury College of Computer Sciences, Northeastern University, Boston, MA, USA
Institutions: Boston Children's Hospital (United States); Harvard University (United States); Johns Hopkins University (United States); Johns Hopkins Medicine (United States); Northeastern University (United States); National Institutes of Health (United States); National Center for Advancing Translational Sciences (United States); University of California, Los Angeles (United States); The University of Texas MD Anderson Cancer Center (United States)
Journal: Cell reports. Medicine, volume 7, issue 5, article 102787
Dates: received 11 December 2025; accepted 8 April 2026; published online 6 May 2026; in print May 2026
Type: Research article · Language: English
License: CC BY-NC-ND
Identifiers: DOI 10.1016/j.xcrm.2026.102787 · PMID 42097147 · PMCID PMC13198234 · OpenAlex W7160418959
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: human (organism), mouse (organism), other condition (population)
Methods: Statistics, Machine learning, Smoothing, state filtering, decompositions
Keywords: chemotherapy-induced peripheral neuropathy, iPSC-derived sensory neurons, high-throughput screening, neuroprotective small molecules, axon degeneration, STE20 kinases
MeSH: Antineoplastic Agents*, High-Throughput Screening Assays*, Induced Pluripotent Stem Cells*, Neuroprotective Agents*, Peripheral Nervous System Diseases*, Sensory Receptor Cells*, Animals, Axons, Disease Models, Animal, Drug Discovery, Ganglia, Spinal, Humans, Mice, Paclitaxel, Protein Kinase Inhibitors (* major topic)
Topic: Cancer Treatment and Pharmacology (Oncology, Medicine), according to OpenAlex
Funding: Defense Advanced Research Projects Agency Defense Sciences Office (HR0011-19-2-0022); Defense Advanced Research Projects Agency; National Institute of Neurological Disorders and Stroke (R35NS105076, 5P50HD105351); Dr. Miriam and Sheldon G. Adelson Medical Research Foundation; Advanced Research Projects Agency for Health; Cancer Prevention and Research Institute of Texas (RR220032); Pharmaceutical Research and Manufacturers of America Foundation; National Institutes of Health
Citations: cited by 1 paper (Europe PMC); 78 references in the paper
Research resources: Chicken anti-Neurofilament H RRID:AB_11212161, Donkey anti-mouse Alexa Fluor 488 RRID:AB_141607, Rabbit anti-TAC1 RRID:AB_1623286, Rabbit anti-phospho-cJun (Ser73) RRID:AB_2129575, Donkey anti-rabbit Alex Fluor 568 RRID:AB_2534017, Donkey anti-rabbit Alexa Fluor 647 RRID:AB_2536183, Rabbit anti-PGP9.5 RRID:AB_2862594, Donkey anti-chicken Alexa Fluor 647 RRID:AB_2921074, Hoechst 33342 RRID:AB_3675235, Mouse anti-βIII-tubulin RRID:AB_477590, Rabbit anti-CGRP RRID:AB_572217, Anti-β-tubulin Alexa Fluor 555 RRID:AB_823665, Rabbit anti-Peripherin RRID:AB_90725, Rabbit anti-Brn3a RRID:AB_92154, R v4.4.2 RRID:SCR_001905, GraphPad Prism RRID:SCR_002798, limma v3.62.2 RRID:SCR_010943, DESeq2 v1.46.0 RRID:SCR_015687, MetaXpress Custom Module Editor RRID:SCR_016654, tximport v1.34.0 RRID:SCR_016752, Salmon v1.10.3 RRID:SCR_017036, Seaborn (Python) RRID:SCR_018132, the IncuCyte RRID:SCR_019874, RRID:SCR_020294, fgsea v1.32.4 RRID:SCR_020938, msigdbr v24.1.0 RRID:SCR_022870, IncuCyte RRID:SCR_025411
AbstractThe abstract is not reproduced here: the paper's license (CC BY-NC-ND) does not allow it. Read it in the paper, at the publisher or on Europe PMC.

Repository

Its files are read in the Code ↔ Paper reader above, with 9 matches between paragraphs and lines of code.

labsyspharm/peripheral-neuropathy

License: none: the authors keep all their rights

State: the link answers, verified on 28 September 2026

Evidence: files inventoried

Commit: 33293433fa769ac9d5a43bded03ad322ed2771c1, 2 December 2025

Languages: R (7)

Size: 8 files, 7 scripts

Software Heritage: not archived

Found in: “Data and code availability”

Holds: 7 notebooks

Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation

Tools: data.table (7 files), tidyverse (7 files), circlize (5 files), ComplexHeatmap (5 files), broom (3 files), Plotly (3 files), DESeq2 (1 file), edgeR (1 file), emmeans (1 file), ggplot2 (1 file), patchwork (1 file), psych (1 file)

Availability: 1 check, the latest on 28 September 2026: the link answers

28 September 2026: the link answers

7 files

de_clustering.Rmd, R, 747 lines
de_comparison.Rmd, R, 206 lines
gene_set_analysis.Rmd, R, 3,522 lines
gene_set_analysis_plotting.Rmd, R, 3,506 lines, 5 matches
marker_expression.Rmd, R, 408 lines, 2 matches
rnaseq/deseq.Rmd, R, 1,590 lines, 1 match
rnaseq/plot_dge.Rmd, R, 204 lines, 1 match

At the source: github.com/labsyspharm/peripheral-neuropathy

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
7 scripts, each with its path and the digest of its content;
9 matches between paragraphs of the paper and lines of the code (method lexical-v1);
neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

geo:GSE298389, at NCBI GEO; found in “Data and code availability”

Code and data availability statement

The paper has a code and data availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

it points to a dataset: NCBI GEO GSE298389
it points to the authors' code: labsyspharm/peripheral-neuropathy
it says that the data are available on request
it says that the code is available on request

Read it in the paper: doi.org/10.1016/j.xcrm.2026.102787.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 2, 28 September 2026

Authors: added Clifford J Woolf (0000-0002-6636-3897); removed Clifford J Woolf

Version 1, 28 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 23 authors, 6 keywords, 15 MeSH terms, 8 funders, 77 references, 27 RRIDs.

Cite

This paper

Petrova, V., Mills, C. E., Hug, C., Cetinkaya-Fisgin, A., Splaine, J., Fouladzadeh, S., Hakim, S., Powell, R., Zhen, S., Chung, M., Bradshaw, G. A., Deng, T., Singec, I., Wang, Q., Kawaguchi, R., Jonnagaddala, H., Barrett, L. B., Smith, J. A., Kalocsay, M., . . . Woolf, C. J. (2026). A human iPSC-derived sensory neuron platform for high-throughput discovery of neuroprotectants against chemotherapy-induced peripheral neuropathy. Cell reports. Medicine, 7(5), 102787. https://doi.org/10.1016/j.xcrm.2026.102787

BibTeX

@article{petrova2026human, author = {Petrova, Veselina and Mills, Caitlin E and Hug, Clemens and Cetinkaya-Fisgin, Aysel and Splaine, Jennifer and Fouladzadeh, Sepideh and Hakim, Sara and Powell, Rasheen and Zhen, Shannon and Chung, Mirra and Bradshaw, Gary A and Deng, Tao and Singec, Ilyas and Wang, Qing and Kawaguchi, Riki and Jonnagaddala, Harathi and Barrett, Lee B and Smith, Jennifer A and Kalocsay, Marian and Gyori, Benjamin M and Hoke, Ahmet and Sorger, Peter K and Woolf, Clifford J}, title = {{A human iPSC-derived sensory neuron platform for high-throughput discovery of neuroprotectants against chemotherapy-induced peripheral neuropathy}}, journal = {Cell reports. Medicine}, year = {2026}, month = may, volume = {7}, number = {5}, pages = {102787}, publisher = {Elsevier}, issn = {2666-3791}, doi = {10.1016/j.xcrm.2026.102787}, url = {https://doi.org/10.1016/j.xcrm.2026.102787}, pmid = {42097147}, pmcid = {PMC13198234} }

RIS

TY - JOUR AU - Petrova, Veselina AU - Mills, Caitlin E AU - Hug, Clemens AU - Cetinkaya-Fisgin, Aysel AU - Splaine, Jennifer AU - Fouladzadeh, Sepideh AU - Hakim, Sara AU - Powell, Rasheen AU - Zhen, Shannon AU - Chung, Mirra AU - Bradshaw, Gary A AU - Deng, Tao AU - Singec, Ilyas AU - Wang, Qing AU - Kawaguchi, Riki AU - Jonnagaddala, Harathi AU - Barrett, Lee B AU - Smith, Jennifer A AU - Kalocsay, Marian AU - Gyori, Benjamin M AU - Hoke, Ahmet AU - Sorger, Peter K AU - Woolf, Clifford J TI - A human iPSC-derived sensory neuron platform for high-throughput discovery of neuroprotectants against chemotherapy-induced peripheral neuropathy T2 - Cell reports. Medicine J2 - Cell Rep Med PY - 2026 DA - 2026/05/06 VL - 7 IS - 5 SP - 102787 SN - 2666-3791 PB - Elsevier DO - 10.1016/j.xcrm.2026.102787 UR - https://doi.org/10.1016/j.xcrm.2026.102787 LA - en ER -

CSL-JSON

{ "id": "10.1016/j.xcrm.2026.102787", "type": "article-journal", "title": "A human iPSC-derived sensory neuron platform for high-throughput discovery of neuroprotectants against chemotherapy-induced peripheral neuropathy", "container-title": "Cell reports. Medicine", "author": [ { "family": "Petrova", "given": "Veselina" }, { "family": "Mills", "given": "Caitlin E" }, { "family": "Hug", "given": "Clemens" }, { "family": "Cetinkaya-Fisgin", "given": "Aysel" }, { "family": "Splaine", "given": "Jennifer" }, { "family": "Fouladzadeh", "given": "Sepideh" }, { "family": "Hakim", "given": "Sara" }, { "family": "Powell", "given": "Rasheen" }, { "family": "Zhen", "given": "Shannon" }, { "family": "Chung", "given": "Mirra" }, { "family": "Bradshaw", "given": "Gary A" }, { "family": "Deng", "given": "Tao" }, { "family": "Singec", "given": "Ilyas" }, { "family": "Wang", "given": "Qing" }, { "family": "Kawaguchi", "given": "Riki" }, { "family": "Jonnagaddala", "given": "Harathi" }, { "family": "Barrett", "given": "Lee B" }, { "family": "Smith", "given": "Jennifer A" }, { "family": "Kalocsay", "given": "Marian" }, { "family": "Gyori", "given": "Benjamin M" }, { "family": "Hoke", "given": "Ahmet" }, { "family": "Sorger", "given": "Peter K" }, { "family": "Woolf", "given": "Clifford J" } ], "container-title-short": "Cell Rep Med", "volume": "7", "issue": "5", "page": "102787", "DOI": "10.1016/j.xcrm.2026.102787", "PMID": "42097147", "PMCID": "PMC13198234", "ISSN": "2666-3791", "publisher": "Elsevier", "URL": "https://doi.org/10.1016/j.xcrm.2026.102787", "language": "en", "issued": { "date-parts": [ [ 2026, 5, 6 ] ] } }

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1093/bioinformatics/btag592 [code]: Network-based stratification of allele-specific expression reveals patient subgroups in Huntington's disease.
Journal: Bioinformatics (Oxford, England)
In common: edgeR, broom, circlize, 8 other tools, other condition, 1 reference
[2] doi:10.1038/s42003-026-10180-5 [code]: Modeling Friedreich's ataxia with Bergmann glia-enriched human cerebellar organoids.
Journal: Communications biology
In common: circlize, DESeq2, ComplexHeatmap, 3 other tools, other condition, 1 reference, author Ilyas Singeç
[3] doi:10.1016/j.isci.2026.115657 [code]: Integration of machine learning to develop a disulfidptosis model for predicting glioma prognosis, immunotherapy response, and drug.
Journal: iScience
In common: edgeR, psych, broom, 7 other tools, other condition
[4] doi:10.1016/j.xcrm.2026.102766 [code]: A longitudinal single-cell and spatial multiomic atlas of pediatric high-grade glioma.
Journal: Cell reports. Medicine
In common: edgeR, circlize, DESeq2, 6 other tools, other condition, 2 references
[5] doi:10.1038/s42003-026-10957-8 [code]: Brain defence by the extracellular matrix protein Cochlin.
Journal: Communications biology
In common: edgeR, broom, circlize, 7 other tools, mouse
[6] doi:10.1038/s41467-026-70243-3 [code]: TYK2 mediates neuroinflammation in Alzheimer's disease brains with TDP-43 pathology.
Journal: Nature communications
In common: broom, circlize, ComplexHeatmap, 3 other tools, 1 reference, author Marian Kalocsay
[7] doi:10.1016/j.xcrm.2026.102682 [code]: TET CpG sequence-context-specific DNA demethylation shapes progression of IDH-mutant gliomas.
Journal: Cell reports. Medicine
In common: edgeR, broom, circlize, 6 other tools, other condition, 1 reference
[8] doi:10.1038/s41467-026-74753-y [code]: A human-specific microRNA controls the timing of excitatory synaptogenesis.
Journal: Nature communications
In common: edgeR, circlize, emmeans, 6 other tools, 1 reference
[9] doi:10.1038/s41386-026-02406-1 [code]: Functional genomic profiling of schizophrenia-associated genes reveals key microglial regulators.
Journal: Neuropsychopharmacology : official publication of the American College of Neuropsychopharmacology
In common: broom, circlize, DESeq2, 6 other tools, 1 reference
[10] doi:10.1016/j.isci.2026.115573 [code]: Female cortical cellular mosaicism underlies shared MeCP2 and PCB impacted gene pathways.
Journal: iScience
In common: edgeR, broom, circlize, 6 other tools, mouse

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.