Contents

1 Setup and data

source("../utils/utils.R")
   ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
   ✔ dplyr     1.1.0     ✔ readr     2.1.4
   ✔ forcats   1.0.0     ✔ stringr   1.5.0
   ✔ ggplot2   3.4.1     ✔ tibble    3.2.0
   ✔ 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 ]8;;http://conflicted.r-lib.org/conflicted package]8;; 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
   
   
   
   
   Registered S3 method overwritten by 'gplots':
     method         from 
     reorder.factor gdata
   
   ChIPseeker v1.34.1  For help: https://guangchuangyu.github.io/software/ChIPseeker
   
   If you use ChIPseeker in published research, please cite:
   Qianwen Wang, Ming Li, Tianzhi Wu, Li Zhan, Lin Li, Meijun Chen, Wenqin Xie, Zijing Xie, Erqiang Hu, Shuangbin Xu, Guangchuang Yu. Exploring epigenomic datasets by ChIPseeker. Current Protocols 2022, 2(10): e585
   
   Loading required package: graph
   
   
   Attaching package: 'graph'
   
   
   The following object is masked from 'package:Biostrings':
   
       complement
   
   
   The following object is masked from 'package:stringr':
   
       boundary
   
   
   Loading required package: Biobase
   
   Welcome to Bioconductor
   
       Vignettes contain introductory material; view with
       'browseVignettes()'. To cite Bioconductor, see
       'citation("Biobase")', and for packages 'citation("pkgname")'.
   
   
   Loading required package: GO.db
   
   Loading required package: AnnotationDbi
   
   
   Attaching package: 'AnnotationDbi'
   
   
   The following object is masked from 'package:dplyr':
   
       select
   
   
   Loading required package: SparseM
   
   
   Attaching package: 'SparseM'
   
   
   The following object is masked from 'package:base':
   
       backsolve
   
   
   
   groupGOTerms:    GOBPTerm, GOMFTerm, GOCCTerm environments built.
   
   
   Attaching package: 'topGO'
   
   
   The following object is masked from 'package:grid':
   
       depth
   
   
   The following object is masked from 'package:IRanges':
   
       members
   
   
   Loading required package: GenomicFeatures
   
   
   Attaching package: 'GenomicFeatures'
   
   
   The following object is masked from 'package:topGO':
   
       genes
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 2A Examples of variants CTCF 6-8h

# counts per line
ll_ctcf = get_counts_per_line("ctcf/68") 
   [1] "399_399_1"
   [1] "399_399_2"
   [1] "vgn_28_1"
   [1] "vgn_28_2"
   [1] "vgn_307_1"
   [1] "vgn_307_2"
   [1] "vgn_399_1"
   [1] "vgn_399_2"
   [1] "vgn_57_1"
   [1] "vgn_57_2"
   [1] "vgn_639_1"
   [1] "vgn_639_2"
   [1] "vgn_712_1"
   [1] "vgn_712_2"
   [1] "vgn_714_1"
   [1] "vgn_714_2"
   [1] "vgn_852_1"
   [1] "vgn_852_2"
   [1] "vgn_vgn_1"
   [1] "vgn_vgn_2"
#vid = "chr3R_9194649"
vid = "chr3R_9196035"
l = plot_ai_and_read_depth_for_variant(cht %>% filter(condition == "ctcf/68"), ll_ctcf, vid) 
   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.
outf = file.path(outdir_fig_main, paste0("Fig2A_example_ctcf_tot_counts.pdf"))
ggsave(outf, l[[1]], width = 3, height = 3)

outf = file.path(outdir_fig_main, paste0("Fig2A_example_ctcf_allele_ratios.pdf"))
ggsave(outf, l[[2]], width = 3, height = 3)
tmp = l[[3]]

tmp = cht %>% filter(peak_id == "chr3R_9196039_ctcf/68") %>% mutate(dist = TEST.SNP.POS - peak_summit)
p = ggplot(tmp, aes(x = dist, y = AI_abs, color = signif_strongAI)) + geom_point(size = 1) + 
  theme_bw() +
  scale_color_manual(values = c("darkgrey", cbPalette[2])) +
  xlab("Distance to peak summit") +
  ylab("Allele imbalance")
p

outf = file.path(outdir_fig_main, paste0("Fig2A_example_ctcf_all_variants_for_peak.pdf"))
ggsave(outf, p, width = 4, height = 2)

3 Figure 2B Examples of variants Bin 6-8h

ll_bin = get_counts_per_line("bin/68") 
   [1] "399_399_1"
   [1] "399_399_2"
   [1] "vgn_28_1"
   [1] "vgn_28_2"
   [1] "vgn_307_1"
   [1] "vgn_307_2"
   [1] "vgn_399_1"
   [1] "vgn_399_2"
   [1] "vgn_57_1"
   [1] "vgn_57_2"
   [1] "vgn_639_1"
   [1] "vgn_639_2"
   [1] "vgn_712_1"
   [1] "vgn_712_2"
   [1] "vgn_714_1"
   [1] "vgn_714_2"
   [1] "vgn_852_1"
   [1] "vgn_852_2"
   [1] "vgn_vgn_1"
   [1] "vgn_vgn_2"
vid = "chr3R_26862967"

l = plot_ai_and_read_depth_for_variant(cht %>% filter(condition == "bin/68"), ll_bin, vid) 


outf = file.path(outdir_fig_main, paste0("Fig2B_example_bin68_tot_counts.pdf"))
ggsave(outf, l[[1]], width = 3, height = 3)

outf = file.path(outdir_fig_main, paste0("Fig2B_example_bin68_allele_ratios.pdf"))
ggsave(outf, l[[2]], width = 3, height = 3)

l[[3]] %>% select(sample_id, TEST.SNP.HAPLOTYPE) %>% arrange(TEST.SNP.HAPLOTYPE)
# peak - chr3R_26863200_bin/68
tmp = cht %>% filter(peak_id == "chr3R_26863200_bin/68") %>% mutate(dist = TEST.SNP.POS - peak_summit)
p = ggplot(tmp, aes(x = dist, y = AI_abs, color = signif_strongAI)) + geom_point(size = 1) + 
  theme_bw() +
  scale_color_manual(values = c("darkgrey", cbPalette[2])) +
  xlab("Distance to peak summit") +
  ylab("Allele imbalance")
p

outf = file.path(outdir_fig_main, paste0("Fig2B_example_bin_all_variants_for_peak.pdf"))
ggsave(outf, p, width = 4, height = 2)

Get variants completely linked with chr3R_26862967 (peak chr3R_26863200_bin/68)

variants = cht %>% filter(condition == "bin/68" & peak_id == "chr3R_26863200_bin/68" & signif_strongAI) %>% select(snp_id) %>% unlist()

lines_var = lapply(ll_bin, function(l) l %>% mutate(snp_id = paste(CHROM, TEST.SNP.POS, sep = "_")) %>% filter(snp_id %in% variants)) %>% bind_rows()

lines_var %>% select(snp_id, sample_id, TEST.SNP.GENOTYPE) %>% unique() %>% spread(sample_id, TEST.SNP.GENOTYPE)
lines_var %<>% select(snp_id, sample_id, TEST.SNP.GENOTYPE) %>% unique() %>% 
  group_by(snp_id) %>%
  summarize(genotype_line = paste(TEST.SNP.GENOTYPE, collapse = "."))
  
lines_var %<>% group_by(genotype_line) %>% mutate(n = n())

4 Figure 2C QQplots

min.p = 1e-20

for(ab_tp in ab_tp_list) {
  print(ab_tp)
  
  filenames = c("cht_results.txt", "cht_results_as.txt", "cht_results_bnb.txt", "cht_results_permuted.txt")
  data_labels = c("CHT full", "CHT AS", "CHT BNB", "CHT permuted")
  names(data_labels) = filenames
  
  res_qq = lapply(filenames, function(f) {cht =  load_cht_results(ab_tp, file_name = f);
                                          lab = data_labels[f];
                                          N = cht %>% filter(signif_strongAI) %>% select(snp_id) %>% unique() %>% nrow()
                                          lab = paste(lab, paste0("N=", N), sep = "\n")
                                          cat(lab)
                                          get_qq_values(cht, lab)
                                          }) %>% bind_rows()
  
  res_qq$label = factor(res_qq$label, levels = unique(res_qq$label))
  
  p = ggplot(res_qq, aes(x, y, color = label)) + 
    geom_point_rast(size = 0.2) + 
    geom_abline(intercept = 0, slope = 1) +
    xlab("Null -log10 p-values") + 
    ylab("Observed -log10 p-values") +
    theme_bw() +
    guides(colour = guide_legend(override.aes = list(size=6))) +
    scale_color_manual(name = "", values = cbPalette) +
    theme(axis.text.y = element_text(size=16), axis.text.x = element_text(size=16), 
          axis.title.x = element_text(size=16), axis.title.y = element_text(size=16),
          legend.text=element_text(size=14), legend.title=element_text(size=16),
          legend.key = element_rect(size = 5),
          legend.key.size = unit(3, 'lines'))
  
  print(p)
  outf = file.path(outdir_fig_main, paste0("Fig2C_qqplot_rastr_", gsub("\\/", "_", ab_tp), ".pdf"))
  ggsave(outf, p, width = 6, height = 4)
  
}
   [1] "twi/24"
   CHT full
   N=13542CHT AS
   N=15790CHT BNB
   N=501CHT permuted
   N=91
   Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
   ℹ Please use the `linewidth` argument instead.
   This warning is displayed once every 8 hours.
   Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
   generated.

   [1] "ctcf/68"
   CHT full
   N=12789CHT AS
   N=13139CHT BNB
   N=4788CHT permuted
   N=255

   [1] "mef2/68"
   CHT full
   N=7956CHT AS
   N=9021CHT BNB
   N=1780CHT permuted
   N=40

   [1] "mef2/1012"
   CHT full
   N=9889CHT AS
   N=11548CHT BNB
   N=2346CHT permuted
   N=83

   [1] "bin/68"
   CHT full
   N=5801CHT AS
   N=6812CHT BNB
   N=672CHT permuted
   N=16

   [1] "bin/1012"
   CHT full
   N=4253CHT AS
   N=4719CHT BNB
   N=398CHT permuted
   N=12

5 Figure 2D Number of affected peaks (with distance of top variant to summit)

# Number of peaks with variants per condition
cht %<>% 
  group_by(condition) %>%
  mutate(n_peaks = length(unique(peak_id))) %>%
  ungroup()
 
# Number of peaks affected by cis variation
cht_sel = cht %>% 
  filter(signif_strongAI) %>%
  group_by(condition) %>%
  mutate(n_AI_peaks = length(unique(peak_id)),
         share_AI_peaks = n_AI_peaks / mean(n_peaks))   

# Select closest top variant to peak summit (min p-value) 
cht_sel %<>%  
  group_by(condition, peak_id) %>%
  mutate(min_pval = min(P.VALUE)) %>% # select variants with min p-value
  filter(P.VALUE == min_pval) %>%
  mutate(min_dist2summit = min(dist2summit)) %>%
  filter(dist2summit == min_dist2summit) %>%
  mutate(bin_dist = cut(dist2summit, breaks = c(0, 250, 1250, 2500), labels = c("in peak", "<= 1kB", "> 1kB"), include.lowest = T))


df_sum = cht_sel %>%
  group_by(condition, bin_dist, share_AI_peaks) %>%
  summarise(n_bin_dist = n(), share_bin_dist = n_bin_dist / mean(n_AI_peaks), tot_AI = mean(n_AI_peaks))
   `summarise()` has grouped output by 'condition', 'bin_dist'. You can override
   using the `.groups` argument.
df_sum$label = factor(ab_tp_labels[df_sum$condition], levels = ab_tp_labels)
df_sum$bin_dist = factor(df_sum$bin_dist, levels = rev(levels(df_sum$bin_dist)))

label_df = df_sum %>% 
  group_by(label, share_AI_peaks) %>%
  summarize(n_AI_peaks = sum(n_bin_dist), share_AI_peaks = round(share_AI_peaks, 2)) %>%
  unique()
   Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
   dplyr 1.1.0.
   ℹ Please use `reframe()` instead.
   ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
     always returns an ungrouped data frame and adjust accordingly.
   Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
   generated.
   `summarise()` has grouped output by 'label', 'share_AI_peaks'. You can override
   using the `.groups` argument.
p = ggplot(df_sum, aes(x = label, y = n_bin_dist, fill = bin_dist)) +
  geom_bar(stat = "identity", width = 0.6) +
  scale_fill_manual(values = cbPalette, name = "Distance to summit\n(nearest top-sign. variant)") +
  annotate("text", x = 1:6, y = label_df$n_AI_peaks + 30, label = label_df$share_AI_peaks, size = 5) +
  xlab("") +
  ylab("# of peaks") +
  theme_bw() +
  theme(axis.text.x = element_text(size=16, angle = 45, hjust = 1, color = TFcols), axis.text.y = element_text(size=16), 
        axis.title.x = element_text(size=18), axis.title.y = element_text(size=16),
        legend.text=element_text(size=14), legend.title=element_text(size=14))
   Warning: Vectorized input to `element_text()` is not officially supported.
   ℹ Results may be unexpected or may change in future versions of ggplot2.
outf = file.path(outdir_fig_main, paste0("Fig2D_number_aff_peaks.pdf"))
ggsave(outf, p, width = 7, height = 5)

p

df_sum

6 Figure 2 E ATAC-Seq

# CHT results

cht = lapply(ab_tp_list, function(ab_tp) load_cht_results(ab_tp)) %>% bind_rows()
cht_gr = cht %>% mutate(chr = TEST.SNP.CHROM, start = TEST.SNP.POS, end = TEST.SNP.POS) %>% GRanges()

# F1 ATAC-seq results

# coordinates of ATAC peaks
atac_regions_path = "/g/furlong/project/68_F1_cisreg_ichip/data/F1_paper_multiom/ATAC_feature_location_dm6.bed"
atac_regions = import(atac_regions_path, format = "bed")

# imbalance info for ATAC peaks
atac_ai_path = "/g/furlong/project/68_F1_cisreg_ichip/data/F1_paper_multiom/ATAC_all_peaks_atac_X.txt"
atac_ai = read.delim(atac_ai_path)

# additional annotations for allele imbalance
atac_ai_sum = atac_ai %>% 
  mutate(AI = padj < 0.01 & abs(0.5 - meanprop) > 0.1) %>%
  group_by(feature, time) %>% 
  summarize(n_AI = sum(AI), 
            AI_peak = any(AI),
            mean_AI = abs(0.5 - mean(meanprop)), 
            max_AI = max(abs(0.5 - meanprop)),
            n_het_lines =n())
   `summarise()` has grouped output by 'feature'. You can override using the
   `.groups` argument.
# combine AII info with peak coordinates
atac_f1_df = merge(data.frame(atac_regions), atac_ai_sum, by.x = "name", by.y = "feature") %>% filter(!seqnames %in% c("chrX", "chrY"))
atac_ai_gr = GRanges(atac_f1_df)


# Do the analysis per TF and ATAC time-point


df = data.frame()

for (ab_tp in ab_tp_list) {
  
  tp = gsub("-", "", timepoints[ab_tp])
  
  # get TF peaks
  AI_peaks = get_peaks_from_cht(ab_tp, cht)
  names(mcols(AI_peaks)) = paste0(names(mcols(AI_peaks)), ".tf")
  
  # get ATAC peaks for selected time-point
  AI_atac = atac_ai_gr[atac_ai_gr$time == tp]
  names(mcols(AI_atac)) = paste0(names(mcols(AI_atac)), ".atac")
  
  ov = findOverlaps(AI_peaks, AI_atac)
  
  
  res = cbind.data.frame(mcols(AI_peaks[queryHits(ov)]), mcols(AI_atac[subjectHits(ov)]))
  
  print(ab_tp)
  wt = wilcox.test(mean_AI.atac ~ AI_peak.tf, res, alternative = "less")
  print(wt)
  wt = wilcox.test(max_AI.atac ~ AI_peak.tf, res, alternative = "less")
  print(wt)
  
  df = rbind.data.frame(df, res)
  
}
   `summarise()` has grouped output by 'peak_id', 'condition', 'chr', 'start'. You
   can override using the `.groups` argument.
   [1] "twi/24"
   
    Wilcoxon rank sum test with continuity correction
   
   data:  mean_AI.atac by AI_peak.tf
   W = 1408786, p-value <2e-16
   alternative hypothesis: true location shift is less than 0
   
   
    Wilcoxon rank sum test with continuity correction
   
   data:  max_AI.atac by AI_peak.tf
   W = 1317170, p-value <2e-16
   alternative hypothesis: true location shift is less than 0
   `summarise()` has grouped output by 'peak_id', 'condition', 'chr', 'start'. You
   can override using the `.groups` argument.
   [1] "ctcf/68"
   
    Wilcoxon rank sum test with continuity correction
   
   data:  mean_AI.atac by AI_peak.tf
   W = 493598, p-value = 1.2e-11
   alternative hypothesis: true location shift is less than 0
   
   
    Wilcoxon rank sum test with continuity correction
   
   data:  max_AI.atac by AI_peak.tf
   W = 444450, p-value <2e-16
   alternative hypothesis: true location shift is less than 0
   `summarise()` has grouped output by 'peak_id', 'condition', 'chr', 'start'. You
   can override using the `.groups` argument.
   [1] "mef2/68"
   
    Wilcoxon rank sum test with continuity correction
   
   data:  mean_AI.atac by AI_peak.tf
   W = 947570, p-value <2e-16
   alternative hypothesis: true location shift is less than 0
   
   
    Wilcoxon rank sum test with continuity correction
   
   data:  max_AI.atac by AI_peak.tf
   W = 879208, p-value <2e-16
   alternative hypothesis: true location shift is less than 0
   `summarise()` has grouped output by 'peak_id', 'condition', 'chr', 'start'. You
   can override using the `.groups` argument.
   [1] "mef2/1012"
   
    Wilcoxon rank sum test with continuity correction
   
   data:  mean_AI.atac by AI_peak.tf
   W = 1304888, p-value <2e-16
   alternative hypothesis: true location shift is less than 0
   
   
    Wilcoxon rank sum test with continuity correction
   
   data:  max_AI.atac by AI_peak.tf
   W = 1237730, p-value <2e-16
   alternative hypothesis: true location shift is less than 0
   `summarise()` has grouped output by 'peak_id', 'condition', 'chr', 'start'. You
   can override using the `.groups` argument.
   [1] "bin/68"
   
    Wilcoxon rank sum test with continuity correction
   
   data:  mean_AI.atac by AI_peak.tf
   W = 416654, p-value <2e-16
   alternative hypothesis: true location shift is less than 0
   
   
    Wilcoxon rank sum test with continuity correction
   
   data:  max_AI.atac by AI_peak.tf
   W = 417385, p-value <2e-16
   alternative hypothesis: true location shift is less than 0
   `summarise()` has grouped output by 'peak_id', 'condition', 'chr', 'start'. You
   can override using the `.groups` argument.
   [1] "bin/1012"
   
    Wilcoxon rank sum test with continuity correction
   
   data:  mean_AI.atac by AI_peak.tf
   W = 397406, p-value = 3.2e-07
   alternative hypothesis: true location shift is less than 0
   
   
    Wilcoxon rank sum test with continuity correction
   
   data:  max_AI.atac by AI_peak.tf
   W = 369466, p-value = 1.5e-11
   alternative hypothesis: true location shift is less than 0
df %<>% 
  mutate(label = factor(ab_tp_labels[condition.tf], levels = ab_tp_labels))

df %>% group_by(condition.tf, AI_peak.tf) %>%
  summarize(mean(mean_AI.atac, mean(max_AI.atac)))
   `summarise()` has grouped output by 'condition.tf'. You can override using the
   `.groups` argument.
p = ggplot(df, aes(x = AI_peak.tf,  y = max_AI.atac)) +
  geom_violin(fill = "darkblue", alpha = 0.3) +
  geom_boxplot(width = 0.4, outlier.size = 0.1, fill = "darkblue", alpha = 0.7) +
  theme_bw() +
  stat_compare_means() +
  xlab("AI TF peaks (F1)") +
  ylab("Allele imbalance of ATAC peaks (F1)") +
  theme(axis.text.y = element_text(size=12), axis.text.x = element_text(size=12), 
        axis.title.x = element_text(size=14), axis.title.y = element_text(size=14))

p

ggsave(file.path(outdir_fig_main, "Fig2E_AI_ATAC_combined.pdf"), p, width = 3, height = 5)

7 Figure 2F Fisher’s Test - TSS distal vs proximal

# only quantified peaks for CHT
peaks = lapply(ab_tp_list, function(ab_tp) get_peaks_from_cht(ab_tp, cht, as_granges = T))
   `summarise()` has grouped output by 'peak_id', 'condition', 'chr', 'start'. You
   can override using the `.groups` argument.
   `summarise()` has grouped output by 'peak_id', 'condition', 'chr', 'start'. You
   can override using the `.groups` argument.
   `summarise()` has grouped output by 'peak_id', 'condition', 'chr', 'start'. You
   can override using the `.groups` argument.
   `summarise()` has grouped output by 'peak_id', 'condition', 'chr', 'start'. You
   can override using the `.groups` argument.
   `summarise()` has grouped output by 'peak_id', 'condition', 'chr', 'start'. You
   can override using the `.groups` argument.
   `summarise()` has grouped output by 'peak_id', 'condition', 'chr', 'start'. You
   can override using the `.groups` argument.
names(peaks) = ab_tp_list

peaks = lapply(peaks, function(x) {ov = findOverlaps(x, promoters)
                                   x$TSS = "TSS-distal"
                                   x$TSS[unique(queryHits(ov))] = "TSS-proximal"
                                   x})

peaks_df = lapply(peaks, function(x) as.data.frame(x)) %>% bind_rows()


fts = lapply(ab_tp_list, function(ab_tp) {fisher_test_two_groups(peaks_df, ab_tp,
                                                                 group1 = "AI_peak", group1_val = c(TRUE, FALSE), 
                                                                 group2 = "TSS", group2_val = c("TSS-distal", "TSS-proximal"))}) %>%
  bind_rows()


fts$label = ab_tp_labels
fts$label = factor(fts$label, levels = ab_tp_labels)

fts %<>% mutate(ratio_label = paste(round(r1, 2), round(r2, 2), sep = ":"))

print(fts)
     condition                 comp1                          comp2 n11 n12  n21
   1    twi/24 AI_peak_TRUE_vs_FALSE TSS_TSS-distal_vs_TSS-proximal 563 280 2031
   2   ctcf/68 AI_peak_TRUE_vs_FALSE TSS_TSS-distal_vs_TSS-proximal 344 237  854
   3   mef2/68 AI_peak_TRUE_vs_FALSE TSS_TSS-distal_vs_TSS-proximal 342 220 1689
   4 mef2/1012 AI_peak_TRUE_vs_FALSE TSS_TSS-distal_vs_TSS-proximal 466 219 2079
   5    bin/68 AI_peak_TRUE_vs_FALSE TSS_TSS-distal_vs_TSS-proximal 229 155 1001
   6  bin/1012 AI_peak_TRUE_vs_FALSE TSS_TSS-distal_vs_TSS-proximal 189 124 1164
      n22      r1      r2       pval odds_ratio        label ratio_label
   1 2480 0.66785 0.45023 1.7999e-31     2.4548    Twi, 2-4h   0.67:0.45
   2 1722 0.59208 0.33152 1.1049e-30     2.9257   CTCF, 6-8h   0.59:0.33
   3 2497 0.60854 0.40349 4.9357e-20     2.2978   Mef2, 6-8h    0.61:0.4
   4 2880 0.68029 0.41924 4.6753e-38     2.9471 Mef2, 10-12h   0.68:0.42
   5 1999 0.59635 0.33367 9.5351e-23     2.9494    Bin, 6-8h    0.6:0.33
   6 1884 0.60383 0.38189 5.0732e-14     2.4663  Bin, 10-12h    0.6:0.38
p = ggplot(fts, aes(x = label, y = odds_ratio)) +
  geom_hline(yintercept = 1, color = "darkred") +
  geom_bar(aes(fill = -log10(pval)), color = "darkblue", stat = "identity", position = "dodge", width = 0.5) +
  #scale_fill_manual(name = "", values = c("darkblue", "darkgrey"), labels = c("AI peaks", "non-AI peaks")) +
  geom_text(aes(label = round(r1, 2), x = label, y = odds_ratio + 0.1), data = fts, size = 6) +
  ylab("Fisher's Test Odds Ratio") +
  theme_bw() +
  theme(axis.text.x = element_text(size = 16, angle = 45, hjust = 1, colour = TFcols),
        axis.text.y = element_text(size = 14),
        axis.title.x = element_blank(),
        axis.title.y = element_text(size = 16),
        legend.text = element_text(size=14),
        legend.title = element_text(size=14))
   Warning: Vectorized input to `element_text()` is not officially supported.
   ℹ Results may be unexpected or may change in future versions of ggplot2.
p

outf = file.path(outdir_fig_main, paste0("Fig2F_tss_dist_prox_fisher.pdf"))
ggsave(outf, p, width = 5, height = 5)
res = lapply(c("TSS-distal", "TSS-proximal"), function(TSS_type) {
                                                      ps = lapply(peaks, function(x) {x[x$TSS == TSS_type]})
                                                      lapply(ab_tp_list, function(ab_tp) fisher_test_AIpeaks_cobinding(ab_tp, ps) %>%
                                                               mutate(TSS = TSS_type)) %>% bind_rows()}) %>% 
  bind_rows()  

res$label = ab_tp_labels
res$label = factor(res$label, levels = ab_tp_labels)

print(res)
      condition              comp1                       comp2 n11 n12  n21  n22
   1     twi/24 isAI_TRUE_vs_FALSE is_co_binding_TRUE_vs_FALSE 216 347  927 1104
   2    ctcf/68 isAI_TRUE_vs_FALSE is_co_binding_TRUE_vs_FALSE 133 211  479  375
   3    mef2/68 isAI_TRUE_vs_FALSE is_co_binding_TRUE_vs_FALSE 190 152 1122  567
   4  mef2/1012 isAI_TRUE_vs_FALSE is_co_binding_TRUE_vs_FALSE 220 246 1123  956
   5     bin/68 isAI_TRUE_vs_FALSE is_co_binding_TRUE_vs_FALSE 137  92  720  281
   6   bin/1012 isAI_TRUE_vs_FALSE is_co_binding_TRUE_vs_FALSE 104  85  748  416
   7     twi/24 isAI_TRUE_vs_FALSE is_co_binding_TRUE_vs_FALSE 214  66 2137  343
   8    ctcf/68 isAI_TRUE_vs_FALSE is_co_binding_TRUE_vs_FALSE 184  53 1616  106
   9    mef2/68 isAI_TRUE_vs_FALSE is_co_binding_TRUE_vs_FALSE 173  47 2196  301
   10 mef2/1012 isAI_TRUE_vs_FALSE is_co_binding_TRUE_vs_FALSE 152  67 2386  494
   11    bin/68 isAI_TRUE_vs_FALSE is_co_binding_TRUE_vs_FALSE 134  21 1821  178
   12  bin/1012 isAI_TRUE_vs_FALSE is_co_binding_TRUE_vs_FALSE 101  23 1697  187
           r1      r2       pval odds_ratio          TSS        label
   1  0.38366 0.45643 2.1353e-03    0.74142   TSS-distal    Twi, 2-4h
   2  0.38663 0.56089 5.1600e-08    0.49376   TSS-distal   CTCF, 6-8h
   3  0.55556 0.66430 1.5295e-04    0.63181   TSS-distal   Mef2, 6-8h
   4  0.47210 0.54016 8.8023e-03    0.76141   TSS-distal Mef2, 10-12h
   5  0.59825 0.71928 4.4116e-04    0.58142   TSS-distal    Bin, 6-8h
   6  0.55026 0.64261 1.8332e-02    0.68069   TSS-distal  Bin, 10-12h
   7  0.76429 0.86169 3.9172e-05    0.52056 TSS-proximal    Twi, 2-4h
   8  0.77637 0.93844 1.0712e-13    0.22798 TSS-proximal   CTCF, 6-8h
   9  0.78636 0.87946 2.0268e-04    0.50467 TSS-proximal   Mef2, 6-8h
   10 0.69406 0.82847 2.7561e-06    0.46984 TSS-proximal Mef2, 10-12h
   11 0.86452 0.91096 6.0911e-02    0.62389 TSS-proximal    Bin, 6-8h
   12 0.81452 0.90074 5.6340e-03    0.48412 TSS-proximal  Bin, 10-12h
p = ggplot(res, aes(x = label, y = odds_ratio)) +
  facet_wrap(~TSS) +
  geom_hline(yintercept = 1, color = "darkred") +
  geom_bar(aes(fill = -log10(pval)), color = "darkblue", stat = "identity", position = "dodge", width = 0.5) +
  #scale_fill_manual(name = "", values = c("darkblue", "darkgrey"), labels = c("AI peaks", "non-AI peaks")) +
  geom_text(aes(label = round(r1, 2), x = label, y = odds_ratio + 0.05), data = res, size = 6) +
  ylab("Fisher's Test Odds Ratio") +
  ylim(0,3) +
  theme_bw() +
  theme(axis.text.x = element_text(size = 16, angle = 45, hjust = 1, colour = TFcols),
        axis.text.y = element_text(size = 14),
        axis.title.x = element_blank(),
        axis.title.y = element_text(size = 16),
        legend.text = element_text(size=14),
        legend.title = element_text(size=14))
   Warning: Vectorized input to `element_text()` is not officially supported.
   ℹ Results may be unexpected or may change in future versions of ggplot2.
p

outf = file.path(outdir_fig_main, paste0("Fig2G_cobinding_fisher.pdf"))
ggsave(outf, p, width = 8, height = 5)
LS0tCnRpdGxlOiAiRmlndXJlXzIiCm91dHB1dDoKICAgQmlvY1N0eWxlOjpodG1sX2RvY3VtZW50OgogICAgICB0b2M6IHRydWUKICAgICAgZGZfcHJpbnQ6IHBhZ2VkCiAgICAgIHNlbGZfY29udGFpbmVkOiB0cnVlCiAgICAgIGNvZGVfZG93bmxvYWQ6IHRydWUKICAgICAgaGlnaGxpZ2h0OiB0YW5nbwojYmlibGlvZ3JhcGh5OiBrbm5fbWxfaW50cm8uYmliCmVkaXRvcl9vcHRpb25zOiAKICBjaHVua19vdXRwdXRfdHlwZTogaW5saW5lCi0tLQoKYGBge3Igc3R5bGUsIGVjaG89RkFMU0UsIHJlc3VsdHM9ImFzaXMifQpsaWJyYXJ5KCJrbml0ciIpCm9wdGlvbnMoZGlnaXRzID0gMiwgd2lkdGggPSA4MCkKb3B0aW9ucyhiaXRtYXBUeXBlID0gJ2NhaXJvJykKZ29sZGVuX3JhdGlvIDwtICgxICsgc3FydCg1KSkgLyAyCm9wdHNfY2h1bmskc2V0KGVjaG8gPSBUUlVFLCB0aWR5ID0gRkFMU0UsIGluY2x1ZGUgPSBUUlVFLCBjYWNoZSA9IEZBTFNFLAogICAgICAgICAgICAgICBkZXY9YygncG5nJywgJ3BkZicpLCBjb21tZW50ID0gJyAgJywgZHBpID0gMzAwKQoKb3B0aW9ucyhzdHJpbmdzQXNGYWN0b3JzID0gRkFMU0UpCmtuaXRyOjpvcHRzX2NodW5rJHNldChjYWNoZT1GQUxTRSkKb3B0aW9ucyhkaWdpdHMgPSA1KSAgICAgICAgIApgYGAKCiMgU2V0dXAgYW5kIGRhdGEKCmBgYHtyfQpzb3VyY2UoIi4uL3V0aWxzL3V0aWxzLlIiKQpjb25maWcgPSBsb2FkX2NvbmZpZygpCgojIGxvYWQgQ0hUIHJlc3VsdHMKY2h0X2Z1bGwgPSBsYXBwbHkoYWJfdHBfbGlzdCwgZnVuY3Rpb24oYWJfdHApIGxvYWRfY2h0X3Jlc3VsdHMoYWJfdHAsIHJlbW92ZV9jaHIgPSBGKSkgJT4lIGJpbmRfcm93cygpCmNodCA9IGNodF9mdWxsICU+JSBmaWx0ZXIoIVRFU1QuU05QLkNIUk9NICVpbiUgYygiY2hyWCIsICJjaHJZIiwgImNock0iKSkKY2h0X3NpZ24gPSBjaHQgJT4lIGZpbHRlcihzaWduaWZfc3Ryb25nQUkpIAoKIyBnZW5lcyBhbmQgcHJvbW90ZXJzCmdlbmVzID0gbG9hZF9nZW5lcygpCnByb21vdGVycyA9IHJlc2l6ZShnZW5lcywgd2lkdGggPSAxMDAwLCBmaXggPSAic3RhcnQiKQoKIyBjb21iaW5lZCBtb3RpZiBzZXQgKGFsbCBURnMsIHBlYWtzICsgYWxsZWxlcykKZmltbyA9IGdldF9mdWxsX21vdGlmX3NldHMoY2h0LCBhYl90cF9saXN0KQojIG9ubHkgYWxsZWxlcwpmaW1vX2FsbGVsZXMgID0gbGFwcGx5KGFiX3RwX2xpc3QsIGZ1bmN0aW9uKGFiX3RwKSBwYXJzZV9tb3RpZnNfaW5fdHdvX2FsbGVsZXMoYWJfdHAsIGNodCkpICU+JSBiaW5kX3Jvd3MoKSAKCmBgYAoKCgojIEZpZ3VyZSAyQSBFeGFtcGxlcyBvZiB2YXJpYW50cyBDVENGIDYtOGgKCgoKYGBge3J9CiMgY291bnRzIHBlciBsaW5lCmxsX2N0Y2YgPSBnZXRfY291bnRzX3Blcl9saW5lKCJjdGNmLzY4IikgCgojdmlkID0gImNocjNSXzkxOTQ2NDkiCnZpZCA9ICJjaHIzUl85MTk2MDM1IgpsID0gcGxvdF9haV9hbmRfcmVhZF9kZXB0aF9mb3JfdmFyaWFudChjaHQgJT4lIGZpbHRlcihjb25kaXRpb24gPT0gImN0Y2YvNjgiKSwgbGxfY3RjZiwgdmlkKSAKCm91dGYgPSBmaWxlLnBhdGgob3V0ZGlyX2ZpZ19tYWluLCBwYXN0ZTAoIkZpZzJBX2V4YW1wbGVfY3RjZl90b3RfY291bnRzLnBkZiIpKQpnZ3NhdmUob3V0ZiwgbFtbMV1dLCB3aWR0aCA9IDMsIGhlaWdodCA9IDMpCgpvdXRmID0gZmlsZS5wYXRoKG91dGRpcl9maWdfbWFpbiwgcGFzdGUwKCJGaWcyQV9leGFtcGxlX2N0Y2ZfYWxsZWxlX3JhdGlvcy5wZGYiKSkKZ2dzYXZlKG91dGYsIGxbWzJdXSwgd2lkdGggPSAzLCBoZWlnaHQgPSAzKQoKYGBgCgpgYGB7cn0KdG1wID0gbFtbM11dCgp0bXAgPSBjaHQgJT4lIGZpbHRlcihwZWFrX2lkID09ICJjaHIzUl85MTk2MDM5X2N0Y2YvNjgiKSAlPiUgbXV0YXRlKGRpc3QgPSBURVNULlNOUC5QT1MgLSBwZWFrX3N1bW1pdCkKcCA9IGdncGxvdCh0bXAsIGFlcyh4ID0gZGlzdCwgeSA9IEFJX2FicywgY29sb3IgPSBzaWduaWZfc3Ryb25nQUkpKSArIGdlb21fcG9pbnQoc2l6ZSA9IDEpICsgCiAgdGhlbWVfYncoKSArCiAgc2NhbGVfY29sb3JfbWFudWFsKHZhbHVlcyA9IGMoImRhcmtncmV5IiwgY2JQYWxldHRlWzJdKSkgKwogIHhsYWIoIkRpc3RhbmNlIHRvIHBlYWsgc3VtbWl0IikgKwogIHlsYWIoIkFsbGVsZSBpbWJhbGFuY2UiKQpwCm91dGYgPSBmaWxlLnBhdGgob3V0ZGlyX2ZpZ19tYWluLCBwYXN0ZTAoIkZpZzJBX2V4YW1wbGVfY3RjZl9hbGxfdmFyaWFudHNfZm9yX3BlYWsucGRmIikpCmdnc2F2ZShvdXRmLCBwLCB3aWR0aCA9IDQsIGhlaWdodCA9IDIpCmBgYAoKCiMgRmlndXJlIDJCIEV4YW1wbGVzIG9mIHZhcmlhbnRzIEJpbiA2LThoCgpgYGB7cn0KbGxfYmluID0gZ2V0X2NvdW50c19wZXJfbGluZSgiYmluLzY4IikgCgp2aWQgPSAiY2hyM1JfMjY4NjI5NjciCgpsID0gcGxvdF9haV9hbmRfcmVhZF9kZXB0aF9mb3JfdmFyaWFudChjaHQgJT4lIGZpbHRlcihjb25kaXRpb24gPT0gImJpbi82OCIpLCBsbF9iaW4sIHZpZCkgCgoKb3V0ZiA9IGZpbGUucGF0aChvdXRkaXJfZmlnX21haW4sIHBhc3RlMCgiRmlnMkJfZXhhbXBsZV9iaW42OF90b3RfY291bnRzLnBkZiIpKQpnZ3NhdmUob3V0ZiwgbFtbMV1dLCB3aWR0aCA9IDMsIGhlaWdodCA9IDMpCgpvdXRmID0gZmlsZS5wYXRoKG91dGRpcl9maWdfbWFpbiwgcGFzdGUwKCJGaWcyQl9leGFtcGxlX2JpbjY4X2FsbGVsZV9yYXRpb3MucGRmIikpCmdnc2F2ZShvdXRmLCBsW1syXV0sIHdpZHRoID0gMywgaGVpZ2h0ID0gMykKCmxbWzNdXSAlPiUgc2VsZWN0KHNhbXBsZV9pZCwgVEVTVC5TTlAuSEFQTE9UWVBFKSAlPiUgYXJyYW5nZShURVNULlNOUC5IQVBMT1RZUEUpCmBgYAoKYGBge3J9CiMgcGVhayAtIGNocjNSXzI2ODYzMjAwX2Jpbi82OAp0bXAgPSBjaHQgJT4lIGZpbHRlcihwZWFrX2lkID09ICJjaHIzUl8yNjg2MzIwMF9iaW4vNjgiKSAlPiUgbXV0YXRlKGRpc3QgPSBURVNULlNOUC5QT1MgLSBwZWFrX3N1bW1pdCkKcCA9IGdncGxvdCh0bXAsIGFlcyh4ID0gZGlzdCwgeSA9IEFJX2FicywgY29sb3IgPSBzaWduaWZfc3Ryb25nQUkpKSArIGdlb21fcG9pbnQoc2l6ZSA9IDEpICsgCiAgdGhlbWVfYncoKSArCiAgc2NhbGVfY29sb3JfbWFudWFsKHZhbHVlcyA9IGMoImRhcmtncmV5IiwgY2JQYWxldHRlWzJdKSkgKwogIHhsYWIoIkRpc3RhbmNlIHRvIHBlYWsgc3VtbWl0IikgKwogIHlsYWIoIkFsbGVsZSBpbWJhbGFuY2UiKQpwCgpvdXRmID0gZmlsZS5wYXRoKG91dGRpcl9maWdfbWFpbiwgcGFzdGUwKCJGaWcyQl9leGFtcGxlX2Jpbl9hbGxfdmFyaWFudHNfZm9yX3BlYWsucGRmIikpCmdnc2F2ZShvdXRmLCBwLCB3aWR0aCA9IDQsIGhlaWdodCA9IDIpCmBgYApHZXQgdmFyaWFudHMgY29tcGxldGVseSBsaW5rZWQgd2l0aCBjaHIzUl8yNjg2Mjk2NyAocGVhayBjaHIzUl8yNjg2MzIwMF9iaW4vNjgpCgpgYGB7cn0KdmFyaWFudHMgPSBjaHQgJT4lIGZpbHRlcihjb25kaXRpb24gPT0gImJpbi82OCIgJiBwZWFrX2lkID09ICJjaHIzUl8yNjg2MzIwMF9iaW4vNjgiICYgc2lnbmlmX3N0cm9uZ0FJKSAlPiUgc2VsZWN0KHNucF9pZCkgJT4lIHVubGlzdCgpCgpsaW5lc192YXIgPSBsYXBwbHkobGxfYmluLCBmdW5jdGlvbihsKSBsICU+JSBtdXRhdGUoc25wX2lkID0gcGFzdGUoQ0hST00sIFRFU1QuU05QLlBPUywgc2VwID0gIl8iKSkgJT4lIGZpbHRlcihzbnBfaWQgJWluJSB2YXJpYW50cykpICU+JSBiaW5kX3Jvd3MoKQoKbGluZXNfdmFyICU+JSBzZWxlY3Qoc25wX2lkLCBzYW1wbGVfaWQsIFRFU1QuU05QLkdFTk9UWVBFKSAlPiUgdW5pcXVlKCkgJT4lIHNwcmVhZChzYW1wbGVfaWQsIFRFU1QuU05QLkdFTk9UWVBFKQoKbGluZXNfdmFyICU8PiUgc2VsZWN0KHNucF9pZCwgc2FtcGxlX2lkLCBURVNULlNOUC5HRU5PVFlQRSkgJT4lIHVuaXF1ZSgpICU+JSAKICBncm91cF9ieShzbnBfaWQpICU+JQogIHN1bW1hcml6ZShnZW5vdHlwZV9saW5lID0gcGFzdGUoVEVTVC5TTlAuR0VOT1RZUEUsIGNvbGxhcHNlID0gIi4iKSkKICAKbGluZXNfdmFyICU8PiUgZ3JvdXBfYnkoZ2Vub3R5cGVfbGluZSkgJT4lIG11dGF0ZShuID0gbigpKQpgYGAKCgojIEZpZ3VyZSAyQyBRUXBsb3RzCgpgYGB7ciBxcXBsb3QsIGNhY2hlPVRSVUV9Cm1pbi5wID0gMWUtMjAKCmZvcihhYl90cCBpbiBhYl90cF9saXN0KSB7CiAgcHJpbnQoYWJfdHApCiAgCiAgZmlsZW5hbWVzID0gYygiY2h0X3Jlc3VsdHMudHh0IiwgImNodF9yZXN1bHRzX2FzLnR4dCIsICJjaHRfcmVzdWx0c19ibmIudHh0IiwgImNodF9yZXN1bHRzX3Blcm11dGVkLnR4dCIpCiAgZGF0YV9sYWJlbHMgPSBjKCJDSFQgZnVsbCIsICJDSFQgQVMiLCAiQ0hUIEJOQiIsICJDSFQgcGVybXV0ZWQiKQogIG5hbWVzKGRhdGFfbGFiZWxzKSA9IGZpbGVuYW1lcwogIAogIHJlc19xcSA9IGxhcHBseShmaWxlbmFtZXMsIGZ1bmN0aW9uKGYpIHtjaHQgPSAgbG9hZF9jaHRfcmVzdWx0cyhhYl90cCwgZmlsZV9uYW1lID0gZik7CiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGxhYiA9IGRhdGFfbGFiZWxzW2ZdOwogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBOID0gY2h0ICU+JSBmaWx0ZXIoc2lnbmlmX3N0cm9uZ0FJKSAlPiUgc2VsZWN0KHNucF9pZCkgJT4lIHVuaXF1ZSgpICU+JSBucm93KCkKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgbGFiID0gcGFzdGUobGFiLCBwYXN0ZTAoIk49IiwgTiksIHNlcCA9ICJcbiIpCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGNhdChsYWIpCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGdldF9xcV92YWx1ZXMoY2h0LCBsYWIpCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIH0pICU+JSBiaW5kX3Jvd3MoKQogIAogIHJlc19xcSRsYWJlbCA9IGZhY3RvcihyZXNfcXEkbGFiZWwsIGxldmVscyA9IHVuaXF1ZShyZXNfcXEkbGFiZWwpKQogIAogIHAgPSBnZ3Bsb3QocmVzX3FxLCBhZXMoeCwgeSwgY29sb3IgPSBsYWJlbCkpICsgCiAgICBnZW9tX3BvaW50X3Jhc3Qoc2l6ZSA9IDAuMikgKyAKICAgIGdlb21fYWJsaW5lKGludGVyY2VwdCA9IDAsIHNsb3BlID0gMSkgKwogICAgeGxhYigiTnVsbCAtbG9nMTAgcC12YWx1ZXMiKSArIAogICAgeWxhYigiT2JzZXJ2ZWQgLWxvZzEwIHAtdmFsdWVzIikgKwogICAgdGhlbWVfYncoKSArCiAgICBndWlkZXMoY29sb3VyID0gZ3VpZGVfbGVnZW5kKG92ZXJyaWRlLmFlcyA9IGxpc3Qoc2l6ZT02KSkpICsKICAgIHNjYWxlX2NvbG9yX21hbnVhbChuYW1lID0gIiIsIHZhbHVlcyA9IGNiUGFsZXR0ZSkgKwogICAgdGhlbWUoYXhpcy50ZXh0LnkgPSBlbGVtZW50X3RleHQoc2l6ZT0xNiksIGF4aXMudGV4dC54ID0gZWxlbWVudF90ZXh0KHNpemU9MTYpLCAKICAgICAgICAgIGF4aXMudGl0bGUueCA9IGVsZW1lbnRfdGV4dChzaXplPTE2KSwgYXhpcy50aXRsZS55ID0gZWxlbWVudF90ZXh0KHNpemU9MTYpLAogICAgICAgICAgbGVnZW5kLnRleHQ9ZWxlbWVudF90ZXh0KHNpemU9MTQpLCBsZWdlbmQudGl0bGU9ZWxlbWVudF90ZXh0KHNpemU9MTYpLAogICAgICAgICAgbGVnZW5kLmtleSA9IGVsZW1lbnRfcmVjdChzaXplID0gNSksCiAgICAgICAgICBsZWdlbmQua2V5LnNpemUgPSB1bml0KDMsICdsaW5lcycpKQogIAogIHByaW50KHApCiAgb3V0ZiA9IGZpbGUucGF0aChvdXRkaXJfZmlnX21haW4sIHBhc3RlMCgiRmlnMkNfcXFwbG90X3Jhc3RyXyIsIGdzdWIoIlxcLyIsICJfIiwgYWJfdHApLCAiLnBkZiIpKQogIGdnc2F2ZShvdXRmLCBwLCB3aWR0aCA9IDYsIGhlaWdodCA9IDQpCiAgCn0KCmBgYAoKCiMgRmlndXJlIDJEIE51bWJlciBvZiBhZmZlY3RlZCBwZWFrcyAod2l0aCBkaXN0YW5jZSBvZiB0b3AgdmFyaWFudCB0byBzdW1taXQpCgpgYGB7cn0KCiMgTnVtYmVyIG9mIHBlYWtzIHdpdGggdmFyaWFudHMgcGVyIGNvbmRpdGlvbgpjaHQgJTw+JSAKICBncm91cF9ieShjb25kaXRpb24pICU+JQogIG11dGF0ZShuX3BlYWtzID0gbGVuZ3RoKHVuaXF1ZShwZWFrX2lkKSkpICU+JQogIHVuZ3JvdXAoKQogCiMgTnVtYmVyIG9mIHBlYWtzIGFmZmVjdGVkIGJ5IGNpcyB2YXJpYXRpb24KY2h0X3NlbCA9IGNodCAlPiUgCiAgZmlsdGVyKHNpZ25pZl9zdHJvbmdBSSkgJT4lCiAgZ3JvdXBfYnkoY29uZGl0aW9uKSAlPiUKICBtdXRhdGUobl9BSV9wZWFrcyA9IGxlbmd0aCh1bmlxdWUocGVha19pZCkpLAogICAgICAgICBzaGFyZV9BSV9wZWFrcyA9IG5fQUlfcGVha3MgLyBtZWFuKG5fcGVha3MpKSAgIAoKIyBTZWxlY3QgY2xvc2VzdCB0b3AgdmFyaWFudCB0byBwZWFrIHN1bW1pdCAobWluIHAtdmFsdWUpIApjaHRfc2VsICU8PiUgIAogIGdyb3VwX2J5KGNvbmRpdGlvbiwgcGVha19pZCkgJT4lCiAgbXV0YXRlKG1pbl9wdmFsID0gbWluKFAuVkFMVUUpKSAlPiUgIyBzZWxlY3QgdmFyaWFudHMgd2l0aCBtaW4gcC12YWx1ZQogIGZpbHRlcihQLlZBTFVFID09IG1pbl9wdmFsKSAlPiUKICBtdXRhdGUobWluX2Rpc3Qyc3VtbWl0ID0gbWluKGRpc3Qyc3VtbWl0KSkgJT4lCiAgZmlsdGVyKGRpc3Qyc3VtbWl0ID09IG1pbl9kaXN0MnN1bW1pdCkgJT4lCiAgbXV0YXRlKGJpbl9kaXN0ID0gY3V0KGRpc3Qyc3VtbWl0LCBicmVha3MgPSBjKDAsIDI1MCwgMTI1MCwgMjUwMCksIGxhYmVscyA9IGMoImluIHBlYWsiLCAiPD0gMWtCIiwgIj4gMWtCIiksIGluY2x1ZGUubG93ZXN0ID0gVCkpCgoKZGZfc3VtID0gY2h0X3NlbCAlPiUKICBncm91cF9ieShjb25kaXRpb24sIGJpbl9kaXN0LCBzaGFyZV9BSV9wZWFrcykgJT4lCiAgc3VtbWFyaXNlKG5fYmluX2Rpc3QgPSBuKCksIHNoYXJlX2Jpbl9kaXN0ID0gbl9iaW5fZGlzdCAvIG1lYW4obl9BSV9wZWFrcyksIHRvdF9BSSA9IG1lYW4obl9BSV9wZWFrcykpCgpkZl9zdW0kbGFiZWwgPSBmYWN0b3IoYWJfdHBfbGFiZWxzW2RmX3N1bSRjb25kaXRpb25dLCBsZXZlbHMgPSBhYl90cF9sYWJlbHMpCmRmX3N1bSRiaW5fZGlzdCA9IGZhY3RvcihkZl9zdW0kYmluX2Rpc3QsIGxldmVscyA9IHJldihsZXZlbHMoZGZfc3VtJGJpbl9kaXN0KSkpCgpsYWJlbF9kZiA9IGRmX3N1bSAlPiUgCiAgZ3JvdXBfYnkobGFiZWwsIHNoYXJlX0FJX3BlYWtzKSAlPiUKICBzdW1tYXJpemUobl9BSV9wZWFrcyA9IHN1bShuX2Jpbl9kaXN0KSwgc2hhcmVfQUlfcGVha3MgPSByb3VuZChzaGFyZV9BSV9wZWFrcywgMikpICU+JQogIHVuaXF1ZSgpCgpwID0gZ2dwbG90KGRmX3N1bSwgYWVzKHggPSBsYWJlbCwgeSA9IG5fYmluX2Rpc3QsIGZpbGwgPSBiaW5fZGlzdCkpICsKICBnZW9tX2JhcihzdGF0ID0gImlkZW50aXR5Iiwgd2lkdGggPSAwLjYpICsKICBzY2FsZV9maWxsX21hbnVhbCh2YWx1ZXMgPSBjYlBhbGV0dGUsIG5hbWUgPSAiRGlzdGFuY2UgdG8gc3VtbWl0XG4obmVhcmVzdCB0b3Atc2lnbi4gdmFyaWFudCkiKSArCiAgYW5ub3RhdGUoInRleHQiLCB4ID0gMTo2LCB5ID0gbGFiZWxfZGYkbl9BSV9wZWFrcyArIDMwLCBsYWJlbCA9IGxhYmVsX2RmJHNoYXJlX0FJX3BlYWtzLCBzaXplID0gNSkgKwogIHhsYWIoIiIpICsKICB5bGFiKCIjIG9mIHBlYWtzIikgKwogIHRoZW1lX2J3KCkgKwogIHRoZW1lKGF4aXMudGV4dC54ID0gZWxlbWVudF90ZXh0KHNpemU9MTYsIGFuZ2xlID0gNDUsIGhqdXN0ID0gMSwgY29sb3IgPSBURmNvbHMpLCBheGlzLnRleHQueSA9IGVsZW1lbnRfdGV4dChzaXplPTE2KSwgCiAgICAgICAgYXhpcy50aXRsZS54ID0gZWxlbWVudF90ZXh0KHNpemU9MTgpLCBheGlzLnRpdGxlLnkgPSBlbGVtZW50X3RleHQoc2l6ZT0xNiksCiAgICAgICAgbGVnZW5kLnRleHQ9ZWxlbWVudF90ZXh0KHNpemU9MTQpLCBsZWdlbmQudGl0bGU9ZWxlbWVudF90ZXh0KHNpemU9MTQpKQoKCm91dGYgPSBmaWxlLnBhdGgob3V0ZGlyX2ZpZ19tYWluLCBwYXN0ZTAoIkZpZzJEX251bWJlcl9hZmZfcGVha3MucGRmIikpCmdnc2F2ZShvdXRmLCBwLCB3aWR0aCA9IDcsIGhlaWdodCA9IDUpCgpwCmBgYApgYGB7cn0KZGZfc3VtCmBgYAoKIyBGaWd1cmUgMiBFIEFUQUMtU2VxCgoKYGBge3J9CgojIENIVCByZXN1bHRzCgpjaHQgPSBsYXBwbHkoYWJfdHBfbGlzdCwgZnVuY3Rpb24oYWJfdHApIGxvYWRfY2h0X3Jlc3VsdHMoYWJfdHApKSAlPiUgYmluZF9yb3dzKCkKY2h0X2dyID0gY2h0ICU+JSBtdXRhdGUoY2hyID0gVEVTVC5TTlAuQ0hST00sIHN0YXJ0ID0gVEVTVC5TTlAuUE9TLCBlbmQgPSBURVNULlNOUC5QT1MpICU+JSBHUmFuZ2VzKCkKCiMgRjEgQVRBQy1zZXEgcmVzdWx0cwoKIyBjb29yZGluYXRlcyBvZiBBVEFDIHBlYWtzCmF0YWNfcmVnaW9uc19wYXRoID0gIi9nL2Z1cmxvbmcvcHJvamVjdC82OF9GMV9jaXNyZWdfaWNoaXAvZGF0YS9GMV9wYXBlcl9tdWx0aW9tL0FUQUNfZmVhdHVyZV9sb2NhdGlvbl9kbTYuYmVkIgphdGFjX3JlZ2lvbnMgPSBpbXBvcnQoYXRhY19yZWdpb25zX3BhdGgsIGZvcm1hdCA9ICJiZWQiKQoKIyBpbWJhbGFuY2UgaW5mbyBmb3IgQVRBQyBwZWFrcwphdGFjX2FpX3BhdGggPSAiL2cvZnVybG9uZy9wcm9qZWN0LzY4X0YxX2Npc3JlZ19pY2hpcC9kYXRhL0YxX3BhcGVyX211bHRpb20vQVRBQ19hbGxfcGVha3NfYXRhY19YLnR4dCIKYXRhY19haSA9IHJlYWQuZGVsaW0oYXRhY19haV9wYXRoKQoKIyBhZGRpdGlvbmFsIGFubm90YXRpb25zIGZvciBhbGxlbGUgaW1iYWxhbmNlCmF0YWNfYWlfc3VtID0gYXRhY19haSAlPiUgCiAgbXV0YXRlKEFJID0gcGFkaiA8IDAuMDEgJiBhYnMoMC41IC0gbWVhbnByb3ApID4gMC4xKSAlPiUKICBncm91cF9ieShmZWF0dXJlLCB0aW1lKSAlPiUgCiAgc3VtbWFyaXplKG5fQUkgPSBzdW0oQUkpLCAKICAgICAgICAgICAgQUlfcGVhayA9IGFueShBSSksCiAgICAgICAgICAgIG1lYW5fQUkgPSBhYnMoMC41IC0gbWVhbihtZWFucHJvcCkpLCAKICAgICAgICAgICAgbWF4X0FJID0gbWF4KGFicygwLjUgLSBtZWFucHJvcCkpLAogICAgICAgICAgICBuX2hldF9saW5lcyA9bigpKQoKIyBjb21iaW5lIEFJSSBpbmZvIHdpdGggcGVhayBjb29yZGluYXRlcwphdGFjX2YxX2RmID0gbWVyZ2UoZGF0YS5mcmFtZShhdGFjX3JlZ2lvbnMpLCBhdGFjX2FpX3N1bSwgYnkueCA9ICJuYW1lIiwgYnkueSA9ICJmZWF0dXJlIikgJT4lIGZpbHRlcighc2VxbmFtZXMgJWluJSBjKCJjaHJYIiwgImNoclkiKSkKYXRhY19haV9nciA9IEdSYW5nZXMoYXRhY19mMV9kZikKCgojIERvIHRoZSBhbmFseXNpcyBwZXIgVEYgYW5kIEFUQUMgdGltZS1wb2ludAoKCmRmID0gZGF0YS5mcmFtZSgpCgpmb3IgKGFiX3RwIGluIGFiX3RwX2xpc3QpIHsKICAKICB0cCA9IGdzdWIoIi0iLCAiIiwgdGltZXBvaW50c1thYl90cF0pCiAgCiAgIyBnZXQgVEYgcGVha3MKICBBSV9wZWFrcyA9IGdldF9wZWFrc19mcm9tX2NodChhYl90cCwgY2h0KQogIG5hbWVzKG1jb2xzKEFJX3BlYWtzKSkgPSBwYXN0ZTAobmFtZXMobWNvbHMoQUlfcGVha3MpKSwgIi50ZiIpCiAgCiAgIyBnZXQgQVRBQyBwZWFrcyBmb3Igc2VsZWN0ZWQgdGltZS1wb2ludAogIEFJX2F0YWMgPSBhdGFjX2FpX2dyW2F0YWNfYWlfZ3IkdGltZSA9PSB0cF0KICBuYW1lcyhtY29scyhBSV9hdGFjKSkgPSBwYXN0ZTAobmFtZXMobWNvbHMoQUlfYXRhYykpLCAiLmF0YWMiKQogIAogIG92ID0gZmluZE92ZXJsYXBzKEFJX3BlYWtzLCBBSV9hdGFjKQogIAogIAogIHJlcyA9IGNiaW5kLmRhdGEuZnJhbWUobWNvbHMoQUlfcGVha3NbcXVlcnlIaXRzKG92KV0pLCBtY29scyhBSV9hdGFjW3N1YmplY3RIaXRzKG92KV0pKQogIAogIHByaW50KGFiX3RwKQogIHd0ID0gd2lsY294LnRlc3QobWVhbl9BSS5hdGFjIH4gQUlfcGVhay50ZiwgcmVzLCBhbHRlcm5hdGl2ZSA9ICJsZXNzIikKICBwcmludCh3dCkKICB3dCA9IHdpbGNveC50ZXN0KG1heF9BSS5hdGFjIH4gQUlfcGVhay50ZiwgcmVzLCBhbHRlcm5hdGl2ZSA9ICJsZXNzIikKICBwcmludCh3dCkKICAKICBkZiA9IHJiaW5kLmRhdGEuZnJhbWUoZGYsIHJlcykKICAKfQoKZGYgJTw+JSAKICBtdXRhdGUobGFiZWwgPSBmYWN0b3IoYWJfdHBfbGFiZWxzW2NvbmRpdGlvbi50Zl0sIGxldmVscyA9IGFiX3RwX2xhYmVscykpCgpkZiAlPiUgZ3JvdXBfYnkoY29uZGl0aW9uLnRmLCBBSV9wZWFrLnRmKSAlPiUKICBzdW1tYXJpemUobWVhbihtZWFuX0FJLmF0YWMsIG1lYW4obWF4X0FJLmF0YWMpKSkKCnAgPSBnZ3Bsb3QoZGYsIGFlcyh4ID0gQUlfcGVhay50ZiwgIHkgPSBtYXhfQUkuYXRhYykpICsKICBnZW9tX3Zpb2xpbihmaWxsID0gImRhcmtibHVlIiwgYWxwaGEgPSAwLjMpICsKICBnZW9tX2JveHBsb3Qod2lkdGggPSAwLjQsIG91dGxpZXIuc2l6ZSA9IDAuMSwgZmlsbCA9ICJkYXJrYmx1ZSIsIGFscGhhID0gMC43KSArCiAgdGhlbWVfYncoKSArCiAgc3RhdF9jb21wYXJlX21lYW5zKCkgKwogIHhsYWIoIkFJIFRGIHBlYWtzIChGMSkiKSArCiAgeWxhYigiQWxsZWxlIGltYmFsYW5jZSBvZiBBVEFDIHBlYWtzIChGMSkiKSArCiAgdGhlbWUoYXhpcy50ZXh0LnkgPSBlbGVtZW50X3RleHQoc2l6ZT0xMiksIGF4aXMudGV4dC54ID0gZWxlbWVudF90ZXh0KHNpemU9MTIpLCAKICAgICAgICBheGlzLnRpdGxlLnggPSBlbGVtZW50X3RleHQoc2l6ZT0xNCksIGF4aXMudGl0bGUueSA9IGVsZW1lbnRfdGV4dChzaXplPTE0KSkKCnAKZ2dzYXZlKGZpbGUucGF0aChvdXRkaXJfZmlnX21haW4sICJGaWcyRV9BSV9BVEFDX2NvbWJpbmVkLnBkZiIpLCBwLCB3aWR0aCA9IDMsIGhlaWdodCA9IDUpCgoKYGBgCgojIEZpZ3VyZSAyRiBGaXNoZXIncyBUZXN0IC0gVFNTIGRpc3RhbCB2cyBwcm94aW1hbAoKYGBge3J9CgojIG9ubHkgcXVhbnRpZmllZCBwZWFrcyBmb3IgQ0hUCnBlYWtzID0gbGFwcGx5KGFiX3RwX2xpc3QsIGZ1bmN0aW9uKGFiX3RwKSBnZXRfcGVha3NfZnJvbV9jaHQoYWJfdHAsIGNodCwgYXNfZ3JhbmdlcyA9IFQpKQpuYW1lcyhwZWFrcykgPSBhYl90cF9saXN0CgpwZWFrcyA9IGxhcHBseShwZWFrcywgZnVuY3Rpb24oeCkge292ID0gZmluZE92ZXJsYXBzKHgsIHByb21vdGVycykKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICB4JFRTUyA9ICJUU1MtZGlzdGFsIgogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHgkVFNTW3VuaXF1ZShxdWVyeUhpdHMob3YpKV0gPSAiVFNTLXByb3hpbWFsIgogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHh9KQoKcGVha3NfZGYgPSBsYXBwbHkocGVha3MsIGZ1bmN0aW9uKHgpIGFzLmRhdGEuZnJhbWUoeCkpICU+JSBiaW5kX3Jvd3MoKQoKCmZ0cyA9IGxhcHBseShhYl90cF9saXN0LCBmdW5jdGlvbihhYl90cCkge2Zpc2hlcl90ZXN0X3R3b19ncm91cHMocGVha3NfZGYsIGFiX3RwLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGdyb3VwMSA9ICJBSV9wZWFrIiwgZ3JvdXAxX3ZhbCA9IGMoVFJVRSwgRkFMU0UpLCAKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBncm91cDIgPSAiVFNTIiwgZ3JvdXAyX3ZhbCA9IGMoIlRTUy1kaXN0YWwiLCAiVFNTLXByb3hpbWFsIikpfSkgJT4lCiAgYmluZF9yb3dzKCkKCgpmdHMkbGFiZWwgPSBhYl90cF9sYWJlbHMKZnRzJGxhYmVsID0gZmFjdG9yKGZ0cyRsYWJlbCwgbGV2ZWxzID0gYWJfdHBfbGFiZWxzKQoKZnRzICU8PiUgbXV0YXRlKHJhdGlvX2xhYmVsID0gcGFzdGUocm91bmQocjEsIDIpLCByb3VuZChyMiwgMiksIHNlcCA9ICI6IikpCgpwcmludChmdHMpCgpwID0gZ2dwbG90KGZ0cywgYWVzKHggPSBsYWJlbCwgeSA9IG9kZHNfcmF0aW8pKSArCiAgZ2VvbV9obGluZSh5aW50ZXJjZXB0ID0gMSwgY29sb3IgPSAiZGFya3JlZCIpICsKICBnZW9tX2JhcihhZXMoZmlsbCA9IC1sb2cxMChwdmFsKSksIGNvbG9yID0gImRhcmtibHVlIiwgc3RhdCA9ICJpZGVudGl0eSIsIHBvc2l0aW9uID0gImRvZGdlIiwgd2lkdGggPSAwLjUpICsKICAjc2NhbGVfZmlsbF9tYW51YWwobmFtZSA9ICIiLCB2YWx1ZXMgPSBjKCJkYXJrYmx1ZSIsICJkYXJrZ3JleSIpLCBsYWJlbHMgPSBjKCJBSSBwZWFrcyIsICJub24tQUkgcGVha3MiKSkgKwogIGdlb21fdGV4dChhZXMobGFiZWwgPSByb3VuZChyMSwgMiksIHggPSBsYWJlbCwgeSA9IG9kZHNfcmF0aW8gKyAwLjEpLCBkYXRhID0gZnRzLCBzaXplID0gNikgKwogIHlsYWIoIkZpc2hlcidzIFRlc3QgT2RkcyBSYXRpbyIpICsKICB0aGVtZV9idygpICsKICB0aGVtZShheGlzLnRleHQueCA9IGVsZW1lbnRfdGV4dChzaXplID0gMTYsIGFuZ2xlID0gNDUsIGhqdXN0ID0gMSwgY29sb3VyID0gVEZjb2xzKSwKICAgICAgICBheGlzLnRleHQueSA9IGVsZW1lbnRfdGV4dChzaXplID0gMTQpLAogICAgICAgIGF4aXMudGl0bGUueCA9IGVsZW1lbnRfYmxhbmsoKSwKICAgICAgICBheGlzLnRpdGxlLnkgPSBlbGVtZW50X3RleHQoc2l6ZSA9IDE2KSwKICAgICAgICBsZWdlbmQudGV4dCA9IGVsZW1lbnRfdGV4dChzaXplPTE0KSwKICAgICAgICBsZWdlbmQudGl0bGUgPSBlbGVtZW50X3RleHQoc2l6ZT0xNCkpCgoKcApvdXRmID0gZmlsZS5wYXRoKG91dGRpcl9maWdfbWFpbiwgcGFzdGUwKCJGaWcyRl90c3NfZGlzdF9wcm94X2Zpc2hlci5wZGYiKSkKZ2dzYXZlKG91dGYsIHAsIHdpZHRoID0gNSwgaGVpZ2h0ID0gNSkKYGBgCgoKYGBge3J9CnJlcyA9IGxhcHBseShjKCJUU1MtZGlzdGFsIiwgIlRTUy1wcm94aW1hbCIpLCBmdW5jdGlvbihUU1NfdHlwZSkgewogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBwcyA9IGxhcHBseShwZWFrcywgZnVuY3Rpb24oeCkge3hbeCRUU1MgPT0gVFNTX3R5cGVdfSkKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgbGFwcGx5KGFiX3RwX2xpc3QsIGZ1bmN0aW9uKGFiX3RwKSBmaXNoZXJfdGVzdF9BSXBlYWtzX2NvYmluZGluZyhhYl90cCwgcHMpICU+JQogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBtdXRhdGUoVFNTID0gVFNTX3R5cGUpKSAlPiUgYmluZF9yb3dzKCl9KSAlPiUgCiAgYmluZF9yb3dzKCkgIAoKcmVzJGxhYmVsID0gYWJfdHBfbGFiZWxzCnJlcyRsYWJlbCA9IGZhY3RvcihyZXMkbGFiZWwsIGxldmVscyA9IGFiX3RwX2xhYmVscykKCnByaW50KHJlcykKCnAgPSBnZ3Bsb3QocmVzLCBhZXMoeCA9IGxhYmVsLCB5ID0gb2Rkc19yYXRpbykpICsKICBmYWNldF93cmFwKH5UU1MpICsKICBnZW9tX2hsaW5lKHlpbnRlcmNlcHQgPSAxLCBjb2xvciA9ICJkYXJrcmVkIikgKwogIGdlb21fYmFyKGFlcyhmaWxsID0gLWxvZzEwKHB2YWwpKSwgY29sb3IgPSAiZGFya2JsdWUiLCBzdGF0ID0gImlkZW50aXR5IiwgcG9zaXRpb24gPSAiZG9kZ2UiLCB3aWR0aCA9IDAuNSkgKwogICNzY2FsZV9maWxsX21hbnVhbChuYW1lID0gIiIsIHZhbHVlcyA9IGMoImRhcmtibHVlIiwgImRhcmtncmV5IiksIGxhYmVscyA9IGMoIkFJIHBlYWtzIiwgIm5vbi1BSSBwZWFrcyIpKSArCiAgZ2VvbV90ZXh0KGFlcyhsYWJlbCA9IHJvdW5kKHIxLCAyKSwgeCA9IGxhYmVsLCB5ID0gb2Rkc19yYXRpbyArIDAuMDUpLCBkYXRhID0gcmVzLCBzaXplID0gNikgKwogIHlsYWIoIkZpc2hlcidzIFRlc3QgT2RkcyBSYXRpbyIpICsKICB5bGltKDAsMykgKwogIHRoZW1lX2J3KCkgKwogIHRoZW1lKGF4aXMudGV4dC54ID0gZWxlbWVudF90ZXh0KHNpemUgPSAxNiwgYW5nbGUgPSA0NSwgaGp1c3QgPSAxLCBjb2xvdXIgPSBURmNvbHMpLAogICAgICAgIGF4aXMudGV4dC55ID0gZWxlbWVudF90ZXh0KHNpemUgPSAxNCksCiAgICAgICAgYXhpcy50aXRsZS54ID0gZWxlbWVudF9ibGFuaygpLAogICAgICAgIGF4aXMudGl0bGUueSA9IGVsZW1lbnRfdGV4dChzaXplID0gMTYpLAogICAgICAgIGxlZ2VuZC50ZXh0ID0gZWxlbWVudF90ZXh0KHNpemU9MTQpLAogICAgICAgIGxlZ2VuZC50aXRsZSA9IGVsZW1lbnRfdGV4dChzaXplPTE0KSkKCgpwCgpvdXRmID0gZmlsZS5wYXRoKG91dGRpcl9maWdfbWFpbiwgcGFzdGUwKCJGaWcyR19jb2JpbmRpbmdfZmlzaGVyLnBkZiIpKQpnZ3NhdmUob3V0ZiwgcCwgd2lkdGggPSA4LCBoZWlnaHQgPSA1KQpgYGAKCg==