# **MODIFY THIS CHUNK**
library(here)
kundaje_dir <- trimws(readr::read_lines(here("code/AK_PROJ_DIR.txt")))
doc_id <- "04c"
out <- here( "output/03-chrombpnet/03-syntax/", doc_id); dir.create(out, recursive = TRUE)
figout <- here( "figures/03-chrombpnet/03-syntax", doc_id, "/"); dir.create(figout, recursive = TRUE)
chrombpnet_dir <- here( "output/03-chrombpnet")Here, we load the results of in silico marginalizations to assess motif cooperativity/synergy and syntax constraints.
Firstly, for each composite motif, we focus on the cluster with the most instances, and evaluate synergy in that cell type. We will assign motif pairs as having synergy or not, and having hard or soft syntax (or both or neither).
Secondly, we run the in silico marginalizations for every composite motif in every cell type, to more generally assess specificity of the results and investigate model predictions on motifs for TFs that are not expressed to be co-expressed.
library(dplyr)
library(tidyr)
library(ggplot2)
library(readr)
library(scales)
library(glue)
library(purrr)
library(stringr)
library(ggrepel)
library(cowplot)
library(pheatmap)
library(ggseqlogo)
library(universalmotif)
library(TFBSTools)
# for interactive stuff
library(plotly)
library(reactable)
library(shiny)
script_path <- here( "code/utils/")
source(file.path(script_path, "plotting_config.R"))
source(file.path(script_path, "hdma_palettes.R"))
source(file.path(script_path, "sj_scRNAseq_helpers.R"))
source(file.path(script_path, "chrombpnet_utils.R"))
ggplot2::theme_set(theme_BOR() + theme(strip.background = element_blank()))# load hits
hits_all <- readRDS(file.path(chrombpnet_dir, "03-syntax/01/hits.Rds")) %>% bind_rows()
# load intermediate objects
load(file.path(chrombpnet_dir, "03-syntax/01/intermediate_obj.Rda"))Load merged modisco reports:
cwm_means <- read_tsv(here("output/03-chrombpnet/02-compendium/modisco_compiled/cwm_sums_df.tsv"))## Rows: 6362 Columns: 8
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (4): cluster, pattern_class, motif, component_pattern
## dbl (4): cwm_sum, abs_cwm_sum, cwm_mean, abs_cwm_mean
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.modisco_merged <- read_tsv(here("output/03-chrombpnet/02-compendium/modisco_compiled_anno/modisco_merged_reports.tsv")) %>%
# get short cluster ID
left_join(cluster_meta %>% dplyr::select(Cluster, Cluster_ChromBPNet),
by = c("component_celltype" = "Cluster_ChromBPNet")) %>%
left_join(cwm_means, by = "component_pattern")## Rows: 6362 Columns: 6
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (5): merged_pattern, component_celltype, pattern_class, pattern, compone...
## dbl (1): n_seqlets
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.head(modisco_merged)## # A tibble: 6 × 14
## merged_pattern component_celltype pattern_class.x pattern n_seqlets
## <chr> <chr> <chr> <chr> <dbl>
## 1 neg.Average_12__merged_p… Spleen_c5 neg_patterns patter… 88
## 2 neg.Average_12__merged_p… Spleen_c0 neg_patterns patter… 461
## 3 neg.Average_12__merged_p… Skin_c5 neg_patterns patter… 302
## 4 neg.Average_12__merged_p… Heart_c13 neg_patterns patter… 128
## 5 neg.Average_12__merged_p… Skin_c7 neg_patterns patter… 708
## 6 neg.Average_12__merged_p… Heart_c5 neg_patterns patter… 47
## # ℹ 9 more variables: component_pattern <chr>, Cluster <chr>, cluster <chr>,
## # pattern_class.y <chr>, motif <chr>, cwm_sum <dbl>, abs_cwm_sum <dbl>,
## # cwm_mean <dbl>, abs_cwm_mean <dbl>Path to the full pipeline for clusters of focus:
ism_dir <- file.path(chrombpnet_dir, "03-syntax/04b/in_silico_marginalization_full/")
fs::dir_exists(ism_dir)## /oak/stanford/groups/wjg/skim/projects/HDMA-public/output/03-chrombpnet/03-syntax/04b/in_silico_marginalization_full/
## TRUEList of composites to test:
compo_to_test <- read_tsv(here("code/03-chrombpnet/03-syntax/04b-composites_to_test.tsv")) %>%
filter(test_spacing == "Y")## Rows: 157 Columns: 26
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (15): motif_name, Cluster, category, annotation_broad, component_motifA,...
## dbl (8): n_hits_per_cluster, idx_uniq, variant_idx, total_hits, start_A, en...
## lgl (3): composite_cwm_fwd, componentA_cwm_fwd, componentB_cwm_fwd
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.compo_to_test$motif_name_safe <- gsub("\\/", "\\.", compo_to_test$motif_name)
compo_to_test$done <- ifelse(file.exists(file.path(ism_dir, compo_to_test$motif_name_safe, compo_to_test$Cluster, ".done")),
TRUE,
FALSE)
glue("{sum(compo_to_test$done)/nrow(compo_to_test)*100}%")## 100%compo_to_test <- compo_to_test %>%
dplyr::rename(seqA = test_seqA, seqB = test_seqB)Check for duplicates:
# sort the two motif sequences alphabetically so we can find duplicates even
# when the order is swapped
pairs <- map_dfr(1:nrow(compo_to_test), function(i) {
motifs <- compo_to_test[i, c("seqA", "seqB")] %>% unlist() %>% sort()
data.frame(motif_name = compo_to_test$motif_name[i],
A = motifs[[1]],
B = motifs[[2]])
})
# should be empty if no duplicates.
pairs %>% group_by(A, B) %>% mutate(n = n()) %>% filter(n > 1)## # A tibble: 0 × 4
## # Groups: A, B [0]
## # ℹ 4 variables: motif_name <chr>, A <chr>, B <chr>, n <int>Outputs of the in silico marginalization experiments:
result_outs <- pmap_dfr(list(compo_to_test$motif_name_safe, compo_to_test$Cluster, compo_to_test$motif_name),
~ data.table::fread(file.path(ism_dir, ..1, ..2, "results.tsv"), data.table = FALSE) %>%
mutate(motif_name = ..3)) %>%
left_join(compo_to_test %>% dplyr::select(motif_name, category), by = "motif_name") %>%
dplyr::select(-seqA_palindrome, -seqB_palindrome)Z-scores on the marginal joint effects:
z_scores <- pmap_dfr(list(compo_to_test$motif_name_safe, compo_to_test$Cluster, compo_to_test$motif_name),
function(motif_safe, cluster, motif) {
df <- data.table::fread(file.path(ism_dir, motif_safe, cluster, "Z_scores.tsv"), data.table = FALSE) %>%
mutate(motif_name = motif)
df
}) %>%
left_join(compo_to_test %>% dplyr::select(motif_name, category), by = "motif_name")Summarized effects for solo motifs:
solo_df <- pmap_dfr(list(compo_to_test$motif_name_safe, compo_to_test$Cluster, compo_to_test$motif_name),
~ data.table::fread(file.path(ism_dir, ..1, ..2, "summarized_effects_solo.tsv"),
data.table = FALSE) %>%
mutate(motif_name = ..3)) %>%
left_join(compo_to_test %>% dplyr::select(motif_name, category), by = "motif_name") %>%
dplyr::rename(orientation = key)Summarized effects for spaced motifs:
spacing_df <- pmap_dfr(list(compo_to_test$motif_name_safe, compo_to_test$Cluster, compo_to_test$motif_name),
~ data.table::fread(file.path(ism_dir, ..1, ..2, "summarized_effects_per_spacing.tsv"), data.table = FALSE) %>%
mutate(motif_name = ..3)) %>%
left_join(compo_to_test %>% dplyr::select(motif_name, category, seqA, seqB), by = "motif_name") %>%
dplyr::rename(orientation = key) %>%
rowwise() %>%
mutate(dist_center_to_center = round(nchar(seqA)/2) + spacing + round(nchar(seqB)/2))All experiments in all cell types:
ism_dir_all <- file.path(chrombpnet_dir, "03-syntax/04b/in_silico_marginalization")
fs::dir_exists(ism_dir_all)## /oak/stanford/groups/wjg/skim/projects/HDMA-public/output/03-chrombpnet/03-syntax/04b/in_silico_marginalization
## TRUEList of composites to test:
# get all tests performed (all motif pairs x all clusters)
all_tests <- expand_grid(
motif_name_safe = compo_to_test$motif_name_safe,
cluster = chrombpnet_models_keep$Cluster
)
all_tests$done <- ifelse(fs::file_exists(file.path(ism_dir_all, all_tests$motif_name_safe, all_tests$cluster, ".done")),
TRUE,
FALSE)
sum(all_tests$done)## [1] 26082sum(all_tests$done)/nrow(all_tests)## [1] 1nrow(all_tests)-sum(all_tests$done)## [1] 0all_tests2 <- all_tests %>% filter(done) %>%
left_join(compo_to_test %>% dplyr::select(-Cluster, -done), by = "motif_name_safe")Outputs of the in silico marginalization experiments:
result_outs_allexp <- map2_dfr(all_tests2$motif_name_safe, all_tests2$cluster,
~ data.table::fread(file.path(ism_dir_all, .x, .y, "results.tsv"), data.table = FALSE)
%>% mutate(cluster = .y)) %>%
left_join(compo_to_test %>% dplyr::select(motif_name, category), by = "motif_name") %>%
dplyr::select(-seqA_palindrome, -seqB_palindrome)
nrow(result_outs_allexp)## [1] 26082Outputs of the solo tests:
solo_df_allexp <- pmap_dfr(list(all_tests2$motif_name_safe, all_tests2$cluster, all_tests2$motif_name),
~ data.table::fread(file.path(ism_dir_all, ..1, ..2, "summarized_effects_solo.tsv"),
data.table = FALSE) %>%
mutate(motif_name = ..3, cluster = ..2))
head(solo_df_allexp)## key counts_after_mean counts_after_sd effect_mean effect_sd motif_name
## 1 A 2.405999 0.1453060 0.021411831 0.010908893 1|AP-2_BHLH
## 2 B 2.379369 0.1503053 -0.005217277 0.010656357 1|AP-2_BHLH
## 3 A + B 4.785368 0.2955811 0.016194545 0.020554082 1|AP-2_BHLH
## 4 A 2.436762 0.1218845 -0.001816668 0.001396795 1|AP-2_BHLH
## 5 B 2.429812 0.1221980 -0.008766225 0.003088963 1|AP-2_BHLH
## 6 A + B 4.866573 0.2440730 -0.010582877 0.004265720 1|AP-2_BHLH
## cluster
## 1 Adrenal_c0
## 2 Adrenal_c0
## 3 Adrenal_c0
## 4 Adrenal_c1
## 5 Adrenal_c1
## 6 Adrenal_c1Save intermediates:
save(spacing_df, solo_df, z_scores, result_outs, compo_to_test, result_outs_allexp, solo_df_allexp,
file = glue("{out}/intermediate_obj.Rda"))blues <- RColorBrewer::brewer.pal(6, "Blues")
purples <- RColorBrewer::brewer.pal(6, "Purples")
palette_orientation_hetero <- c('A,B HT' = blues[3],
'B,A HT' = blues[4],
'HH' = blues[5],
'TT' = blues[6])
palette_orientation_homo <- c('HT' = purples[2],
'HH' = purples[4],
'TT' = purples[6])
palette_coop <- c("homo - not cooperative" = "gray70",
"hetero - not cooperative" = "gray70",
"homo - cooperative" = "#7696f5",
"hetero - cooperative" = "#ab7ec2",
"homo - specific" = "royalblue3",
"hetero - specific" = "darkorchid4")
shapemap_result <- c("no synergy" = 1,
"soft" = 24,
"hard & soft" = 23,
"hard" = 22)
cmap_category["not cooperative"] <- "gray80"
palette_motifs <- colorRampPalette(RColorBrewer::brewer.pal(9, "Blues")[3:9])(n=nrow(compo_to_test))We want to combine results into one dataframe:
# join results
all_results_best <- compo_to_test %>%
dplyr::select(motif_name, Cluster, n_hits_per_cluster, idx_uniq, annotation_broad, component_motifA, component_motifB, seqA, seqB) %>%
left_join(result_outs, by = "motif_name") %>%
left_join(z_scores %>% dplyr::select(-category),
by = c("motif_name", "best_orientation" = "orientation", "best_spacing" = "spacing")) %>%
left_join(spacing_df %>%
dplyr::select(-category, -effect_mean, -effect_sd, -seqA, -seqB) %>%
dplyr::rename(best_counts_after_mean = counts_after_mean,
best_counts_after_sd = counts_after_sd),
by = c("motif_name", "best_orientation" = "orientation", "best_spacing" = "spacing")) %>%
left_join(solo_df %>%
dplyr::filter(orientation == "A + B") %>%
dplyr::select(motif_name,
sum_counts_after_mean = counts_after_mean,
sum_counts_after_sd = counts_after_sd,
sum_effect = effect_mean)) %>%
# the maximum difference between joint vs independent effects across all arrangements
mutate(joint_vs_ind_max = joint_effect - sum_effect)## Joining with `by = join_by(motif_name)`# adjust p-value
all_results_best$joint_vs_sum_p_value_adj <- p.adjust(all_results_best$joint_vs_sum_p_value, method = "BH")We want to calculate the difference in joint/independent effects for each arrangement;
all_results_spacing <- spacing_df %>%
dplyr::rename(joint_effect = effect_mean) %>%
left_join(solo_df %>% filter(orientation == "A + B") %>%
dplyr::select(-orientation, -counts_after_mean, -counts_after_sd, -effect_sd, sum_effect = effect_mean),
by = c("motif_name", "category")) %>%
mutate(joint_vs_ind_per_arrangement = joint_effect - sum_effect)Similarly, for the set of tests of all motifs in all cell types, we compute difference between joint and independent effects:
# adjust p-value
result_outs_allexp$joint_vs_sum_p_value_adj <- p.adjust(result_outs_allexp$joint_vs_sum_p_value, method = "BH")
all_results_allexp <- result_outs_allexp %>%
left_join(solo_df_allexp %>% filter(key == "A + B") %>%
mutate(sum_effect = effect_mean) %>%
dplyr::select(-key, -counts_after_mean, -counts_after_sd, -effect_sd, -effect_mean),
by = c("motif_name", "cluster")) %>%
mutate(joint_vs_ind_max = joint_effect - sum_effect) %>%
mutate(log10padj = -log10(joint_vs_sum_p_value_adj))Let’s annotate the best orientation for each motif:
all_results_spacing$is_best_orientation <- FALSE
get_best_orientation <- function(motif) {
all_results_best %>% filter(motif_name == motif) %>% pull(best_orientation)
}
motifs <- unique(all_results_spacing$motif_name)
print(motifs)## [1] "1|AP-2_BHLH"
## [2] "10|BHLH:ATOH/NEUROD_HD:DLX/LHX#1"
## [3] "104|BZIP:ATF/CREB_ETS"
## [4] "105|BZIP:ATF/CREB_ETS:ELF/ETV"
## [5] "106|BZIP:ATF/CREB_MEF2"
## [6] "110|BZIP:ATF/CREB_SP/KLF#2"
## [7] "112|BZIP:ATF/CREB_TEAD#2"
## [8] "12|BHLH:ATOH/NEUROD_MEF2"
## [9] "121|BZIP:CEBP_FOX#2"
## [10] "123|BZIP:CEBP_NFI"
## [11] "124|BZIP:CEBP_NR2/RXR"
## [12] "125|BZIP:CEBP_TEAD"
## [13] "128|BZIP:DBP/HLF_FOX"
## [14] "130|BZIP:DBP/HLF_NKX"
## [15] "139|BZIP:FOSL/JUND_ETS"
## [16] "14|BHLH:MXI1_BHLH:TWIST/ZBTB"
## [17] "140|BZIP:FOSL/JUND_FOX"
## [18] "141|BZIP:FOSL/JUND_IKZF1"
## [19] "142|BZIP:FOSL/JUND_NR:ESRRA/NR5A#1"
## [20] "144|BZIP:FOSL/JUND_TEAD"
## [21] "146|BZIP:MAFG/MAFF_BZIP:MAFG/MAFF"
## [22] "147|BZIP:MAFG/MAFF_TEAD"
## [23] "15|BHLH:NHLH/TFAP_NFY"
## [24] "152|CTCF_TBP"
## [25] "156|EBF_ETS:ELF/ETV"
## [26] "16|BHLH:OLIG/NEUROG_MEF2"
## [27] "162|ETS:ELF/ETV_ETS:ELF/ETV#1"
## [28] "163|ETS:ELF/ETV_FOX"
## [29] "164|ETS:ELF/ETV_NFY"
## [30] "165|ETS:ELF/ETV_SP/KLF#1"
## [31] "17|BHLH:OLIG/NEUROG_SOX"
## [32] "177|ETS:ELF/SPIB_ETS:ELF/SPIB#2"
## [33] "185|ETS:ELF/SPIB_FOX"
## [34] "187|ETS:ELF/SPIB_GATA#2"
## [35] "191|ETS:ELF/SPIB_NFI#4"
## [36] "193|ETS:ELF/SPIB_NFY"
## [37] "199|ETS:ELF/SPIB_NR:ESRRA/NR5A#5"
## [38] "201|ETS:ELF/SPIB_RUNX#1"
## [39] "208|ETS:ELF/SPIB_SOX#6"
## [40] "209|ETS:ELF/SPIB_SP/KLF"
## [41] "215|ETS_ETS#1"
## [42] "217|ETS_FOX#1"
## [43] "219|ETS_IKZF1"
## [44] "221|ETS_NFAT#2"
## [45] "224|ETS_NR:ESRRA/NR5A#2"
## [46] "225|ETS_SOX"
## [47] "226|ETS_SP/KLF#1"
## [48] "240|FOX_FOX#3"
## [49] "242|FOX_GATA#2"
## [50] "245|FOX_NFI"
## [51] "25|BHLH:TFAP/MYOD_BHLH:TFAP/MYOD#2"
## [52] "253|FOX_NKX#7"
## [53] "255|FOX_NR:HNF4A/HNF4G"
## [54] "256|FOX_NR:NR2F/RXRG"
## [55] "257|FOX_RFX"
## [56] "258|FOX_SOX"
## [57] "263|GATA1_TAL1"
## [58] "266|GATA_GATA#3"
## [59] "271|GATA_HD:MEIS/TGIF/TBX#1"
## [60] "273|GATA_MEF2#1"
## [61] "278|GATA_NFI"
## [62] "280|GATA_NR:ESRRA/NR5A#2"
## [63] "281|GATA_NR:NR4A/NR2C"
## [64] "283|GATA_SP/KLF#2"
## [65] "285|GATA_TEAD"
## [66] "288|GRHL_GRHL#1"
## [67] "290|GRHL_TEAD#1"
## [68] "30|BHLH:TFAP/MYOD_BHLH:TWIST/ZBTB"
## [69] "301|HD:DLX/LHX_HD:MEIS/TGIF/TBX"
## [70] "303|HD:DLX/LHX_POU"
## [71] "307|HD:MEIS/TGIF/TBX_HD:MEIS/TGIF/TBX#1"
## [72] "309|HD:MEIS/TGIF/TBX_HD:PAX/VSX"
## [73] "310|HD:MEIS/TGIF/TBX_HD:PITX/OTX"
## [74] "311|HD:MEIS/TGIF/TBX_MEF2"
## [75] "32|BHLH:TFAP/MYOD_BZIP:FOSL/JUND#2"
## [76] "321|HD:PITX/OTX_HD:PITX/OTX#4"
## [77] "322|HD:PITX/OTX_HD:SIX/ZNF"
## [78] "323|HD:PITX/OTX_SP/KLF"
## [79] "326|HD_HD#1"
## [80] "33|BHLH:TFAP/MYOD_ETS"
## [81] "330|HD_NFI#1"
## [82] "338|IKZF1_NR"
## [83] "339|IKZF1_RUNX"
## [84] "34|BHLH:TFAP/MYOD_ETS:ELF/SPIB"
## [85] "340|IRF/STAT_IRF/STAT#1"
## [86] "345|MEF2_NFI"
## [87] "349|MEF2_TBP"
## [88] "35|BHLH:TFAP/MYOD_FOX#1"
## [89] "351|NFAT_NFAT#1"
## [90] "354|NFAT_SOX#1"
## [91] "357|NFI_NR2E1"
## [92] "368|NFY_NFY#1"
## [93] "371|NFY_SP/KLF#2"
## [94] "38|BHLH:TFAP/MYOD_GATA"
## [95] "382|NKX_NKX#5"
## [96] "389|NKX_TEAD#2"
## [97] "39|BHLH:TFAP/MYOD_HD:PBX/PKNOX"
## [98] "393|NR:ESRRA/NR5A_SOX"
## [99] "395|NR:HNF4A/HNF4G_NR:HNF4A/HNF4G"
## [100] "40|BHLH:TFAP/MYOD_ROR"
## [101] "402|NR_NR#2"
## [102] "409|PAX_SOX"
## [103] "41|BHLH:TFAP/MYOD_SP/KLF"
## [104] "413|POU4_POU4"
## [105] "414|POU_POU#1"
## [106] "416|POU_RUNX"
## [107] "422|RFX_RUNX"
## [108] "423|RFX_SOX"
## [109] "424|RFX_SP/KLF"
## [110] "428|RUNX_RUNX#1"
## [111] "430|RUNX_TCF7L/LEF"
## [112] "434|SOX_SOX#1"
## [113] "439|SP/KLF_SP/KLF#3"
## [114] "446|SRF_TEAD"
## [115] "452|TEAD_TEAD"
## [116] "46|BHLH:TWIST/ZBTB_BHLH:TWIST/ZBTB"
## [117] "467|p53_p53#1"
## [118] "47|BHLH:TWIST/ZBTB_HD"
## [119] "48|BHLH:TWIST/ZBTB_HD:PAX/VSX"
## [120] "5|BHLH:ASCL/TCF_ETS:ELF/SPIB"
## [121] "55|BHLH:USF/BHLHE_ETS:ELF/SPIB"
## [122] "56|BHLH:USF/BHLHE_NFY"
## [123] "59|BHLH_BHLH:ATOH/NEUROD"
## [124] "6|BHLH:ASCL/TCF_TCF7L/LEF"
## [125] "60|BHLH_BHLH:MXI1"
## [126] "63|BHLH_BHLH:TWIST/ZBTB#2"
## [127] "65|BHLH_BZIP:ATF/CREB#1"
## [128] "69|BHLH_EBF"
## [129] "71|BHLH_HD#2"
## [130] "72|BHLH_HD:MEIS/TGIF/TBX#1"
## [131] "74|BHLH_HD:PBX1/PDX1"
## [132] "75|BHLH_HD:PITX/OTX#1"
## [133] "79|BHLH_IKZF1"
## [134] "80|BHLH_MEF2#1"
## [135] "84|BHLH_NFI"
## [136] "9|BHLH:ATOH/NEUROD_BHLH:ATOH/NEUROD"
## [137] "91|BHLH_RUNX#3"
## [138] "97|BHLH_p53"for (i in seq_along(motifs)) {
best_ori <- get_best_orientation(motifs[i])
all_results_spacing[all_results_spacing$motif_name == motifs[i] & all_results_spacing$orientation == best_ori, ]$is_best_orientation <- TRUE
}Let’s compute the max difference in joint vs ind. effects at the best orientation, in the 20-150bp range:
tmp <- all_results_spacing %>%
filter(is_best_orientation) %>%
filter(spacing %in% 20:150) %>%
group_by(motif_name) %>%
slice_max(n = 1, order_by = joint_effect) %>%
mutate(joint_vs_ind_med_distance = joint_effect - sum_effect) %>%
dplyr::select(motif_name, joint_vs_ind_med_distance) %>%
ungroup()
all_results_best <- all_results_best %>% left_join(tmp, by = "motif_name")What’s the distribution of counts?
# first, get the mean
motif_order_mean <- spacing_df %>%
group_by(motif_name) %>%
summarize(mean_counts = mean(counts_after_mean)) %>%
arrange(desc(mean_counts)) %>%
pull(motif_name)
spacing_df %>%
mutate(motif_name = factor(motif_name, levels = motif_order_mean)) %>%
mutate(upper = counts_after_mean + counts_after_sd,
lower = counts_after_mean - counts_after_sd) %>%
ggplot(aes(x = motif_name, y = counts_after_mean)) +
geom_point(aes(fill = category), shape = 21, color = "white", size = 2, alpha = 0.5) +
# geom_errorbar(aes(color = category, ymin = lower, ymax = upper), width = 0.1) +
scale_fill_manual(values = cmap_category) +
# scale_color_manual(values = cmap_category) +
facet_grid(. ~ category, space = "free_x", scales = "free_x") +
theme(strip.background = element_blank(),
axis.text.x = element_text(size = 9)) +
rotate_x() +
ylab("Mean predicted log-counts for AB") +
ggtitle("Distribution of predicted counts for all spacings/orientations of each comp. motif") +
no_legend()