Figure_1

Contents

1 Setup and data

source("../utils/utils.R")

   ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
   ✔ dplyr     1.1.2     ✔ readr     2.1.4
   ✔ forcats   1.0.0     ✔ stringr   1.5.0
   ✔ ggplot2   3.4.4     ✔ tibble    3.2.1
   ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
   ✔ purrr     1.0.1     
   ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
   ✖ dplyr::filter() masks stats::filter()
   ✖ dplyr::lag()    masks stats::lag()
   ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
   
   Attaching package: 'magrittr'
   
   
   The following object is masked from 'package:purrr':
   
       set_names
   
   
   The following object is masked from 'package:tidyr':
   
       extract
   
   
   Loading required package: GenomicRanges
   
   Loading required package: stats4
   
   Loading required package: BiocGenerics
   
   
   Attaching package: 'BiocGenerics'
   
   
   The following objects are masked from 'package:lubridate':
   
       intersect, setdiff, union
   
   
   The following objects are masked from 'package:dplyr':
   
       combine, intersect, setdiff, union
   
   
   The following objects are masked from 'package:stats':
   
       IQR, mad, sd, var, xtabs
   
   
   The following objects are masked from 'package:base':
   
       anyDuplicated, aperm, append, as.data.frame, basename, cbind,
       colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
       get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
       match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
       Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort,
       table, tapply, union, unique, unsplit, which.max, which.min
   
   
   Loading required package: S4Vectors
   
   
   Attaching package: 'S4Vectors'
   
   
   The following objects are masked from 'package:lubridate':
   
       second, second<-
   
   
   The following objects are masked from 'package:dplyr':
   
       first, rename
   
   
   The following object is masked from 'package:tidyr':
   
       expand
   
   
   The following objects are masked from 'package:base':
   
       expand.grid, I, unname
   
   
   Loading required package: IRanges
   
   
   Attaching package: 'IRanges'
   
   
   The following object is masked from 'package:lubridate':
   
       %within%
   
   
   The following objects are masked from 'package:dplyr':
   
       collapse, desc, slice
   
   
   The following object is masked from 'package:purrr':
   
       reduce
   
   
   Loading required package: GenomeInfoDb
   
   
   Attaching package: 'GenomicRanges'
   
   
   The following object is masked from 'package:magrittr':
   
       subtract
   
   
   Loading required package: grid
   
   Loading required package: Biostrings
   
   Loading required package: XVector
   
   
   Attaching package: 'XVector'
   
   
   The following object is masked from 'package:purrr':
   
       compact
   
   
   
   Attaching package: 'Biostrings'
   
   
   The following object is masked from 'package:grid':
   
       pattern
   
   
   The following object is masked from 'package:base':
   
       strsplit
   
   
   
   Attaching package: 'gridExtra'
   
   
   The following object is masked from 'package:BiocGenerics':
   
       combine
   
   
   The following object is masked from 'package:dplyr':
   
       combine
   
   
   
   Attaching package: 'data.table'
   
   
   The following object is masked from 'package:GenomicRanges':
   
       shift
   
   
   The following object is masked from 'package:IRanges':
   
       shift
   
   
   The following objects are masked from 'package:S4Vectors':
   
       first, second
   
   
   The following objects are masked from 'package:lubridate':
   
       hour, isoweek, mday, minute, month, quarter, second, wday, week,
       yday, year
   
   
   The following objects are masked from 'package:dplyr':
   
       between, first, last
   
   
   The following object is masked from 'package:purrr':
   
       transpose

config = load_config()

# load CHT results
cht_full = lapply(ab_tp_list, function(ab_tp) load_cht_results(ab_tp, remove_chr = F)) %>% bind_rows()
cht = cht_full %>% filter(!TEST.SNP.CHROM %in% c("chrX", "chrY", "chrM"))
cht_sign = cht %>% filter(signif_strongAI) 

# genes and promoters
genes = load_genes()
promoters = resize(genes, width = 1000, fix = "start")

# combined motif set (all TFs, peaks + alleles)
fimo = get_full_motif_sets(cht, ab_tp_list)
# only alleles
fimo_alleles  = lapply(ab_tp_list, function(ab_tp) parse_motifs_in_two_alleles(ab_tp, cht)) %>% bind_rows() 

2 Figure 1 C motif logos

outf_base = file.path(outdir_fig_main, "/motif_logos/")
f = config$data$motif_databases$denovo_motifs
pfms = read_meme(f)

for(pfm in pfms) {
  
  #motif_name = gsub("\\/", "_", pfm@name)
  motif_name = pfm@name
  print(motif_name)
  outf = file.path(outf_base, paste0(motif_name, ".pdf"))
  
  p = view_motifs(pfm, use.type = "ICM", colour.scheme = letter_colors, sort.positions = T) +
    theme_classic() +
    scale_x_continuous(breaks = 1:ncol(pfm), labels = 1:ncol(pfm), name = "Position") +
    ylab("Bits") +
    theme(axis.text.y = element_text(size=12), axis.text.x = element_text(size=12),
          axis.title.y = element_text(size=12), axis.title.x = element_text(size=12),
          legend.position="none",
          plot.title = element_text(color="black", face="bold", size=16, hjust=0.5))
  
  print(p)
  ggsave(outf, p,  width = 4, height = 2)
  
  
}

   [1] "Twi_denovo"

   Scale for x is already present.
   Adding another scale for x, which will replace the existing scale.

   [1] "CTCF_denovo"

   Scale for x is already present.
   Adding another scale for x, which will replace the existing scale.

   [1] "Mef2_denovo.68h"

   Scale for x is already present.
   Adding another scale for x, which will replace the existing scale.

   [1] "Mef2_denovo.1012h"

   Scale for x is already present.
   Adding another scale for x, which will replace the existing scale.

   [1] "Bin_denovo.68h"

   Scale for x is already present.
   Adding another scale for x, which will replace the existing scale.

   [1] "Bin_denovo.1012h"

   Scale for x is already present.
   Adding another scale for x, which will replace the existing scale.

3 Figure 1 C Motif central enrichment

pp = lapply(ab_tp_list, function(ab_tp) {print(ab_tp); plot_motif_central_enrichment(ab_tp, pwm = motifs_denovo_list[ab_tp])})

   [1] "twi/24"

   Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
   ℹ Please use `linewidth` instead.
   This warning is displayed once every 8 hours.
   Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
   generated.

   [1] "ctcf/68"
   [1] "mef2/68"
   [1] "mef2/1012"
   [1] "bin/68"
   [1] "bin/1012"

p = do.call("grid.arrange", c(pp, ncol= 6))

ggsave(file.path(outdir_fig_main, "Fig1C_motif_central_enrichment.pdf"), p, width = 15, height = 2)

4 Figure 1 D Upset plot

# paths to consensus peak sets for all TFs
ff = lapply(ab_tp_list, function(ab_tp) get_path_consensus_peakset(ab_tp))

# read peak sets
peaks = lapply(ff, function(f) load_consensus_peak_set(f, filter_chr = F)) %>% GRangesList()
names(peaks) = ab_tp_labels

combined_peakset = construct_combined_consensus_set(peaks)

df = combined_peakset[ , 6:ncol(combined_peakset)]
names(df) = ab_tp_labels


#{pdf(file.path(outdir_fig_main, "Fig1D_UpsetPlot_full.pdf"), width = 7, height = 5)

p = upset(df,  order.by = "freq", nintersects = 15, keep.order = T, sets = rev(names(df)),
      sets.bar.color = rev(TFcols), text.scale = 1.5)

print(p)

pdf(file.path(outdir_fig_main, "Fig1D_UpsetPlot_full.pdf"), width = 7, height = 5)
p
dev.off()

   pdf 
     2

5 Figure 1 E Number of variants per peak

# paths to consensus peak sets for all TFs
ff = lapply(ab_tp_list, function(ab_tp) get_path_consensus_peakset(ab_tp))
peaks = lapply(ff, function(f) load_consensus_peak_set(f, filter_chr = F) %>% as.data.frame()) %>% bind_rows() %>% select(peak_id, condition)

# annotate by number of variants per peak
peaks_with_variants = get_number_of_variants_per_peak(cht_full)

   `summarise()` has grouped output by 'peak_id'. You can override using the
   `.groups` argument.

df = merge(peaks, peaks_with_variants, by = c("peak_id", "condition"), all.x = T)
df$n_var[is.na(df$n_var)] = 0

df %<>% mutate(bin_n_var = cut(n_var, breaks = c(-1, 0, 20, 50, 80, Inf), labels = c("0", "1-20", "21-50", "51-80", ">80")))

df_sum = df %>% group_by(condition) %>% 
        mutate(n_tot = n()) %>%
        group_by(condition, bin_n_var) %>%
        summarise(n_bin_var = n(), share_bin_var = n_bin_var / mean(n_tot)) %>%
        mutate(bin_n_var = factor(bin_n_var, levels = rev(levels(bin_n_var))))

   `summarise()` has grouped output by 'condition'. You can override using the
   `.groups` argument.

df_sum$label = factor(ab_tp_labels[df_sum$condition], levels = rev(ab_tp_labels))

p = ggplot(df_sum, aes(x = label, y = share_bin_var, fill = bin_n_var)) +
  geom_bar(stat = "identity", width = 0.6) +
  scale_fill_manual(values = cbPalette, name = "# variants per peak") +
  xlab("") +
  ylab("Share of peaks") +
  theme_bw() +
  ylim(c(0, 1)) +
  coord_flip() +
  theme(axis.text.y = element_text(size=18, color = rev(TFcols)), axis.text.x = element_text(size=16), 
        axis.title.x = element_text(size=18), axis.title.y = element_text(size=16),
        legend.text=element_text(size=16), legend.title=element_text(size=16))

   Warning: Vectorized input to `element_text()` is not officially supported.
   ℹ Results may be unexpected or may change in future versions of ggplot2.

p

ggsave(file.path(outdir_fig_main, "Fig1E_num_vriants_per_peak.pdf"), p, width = 8, height = 3)