---
title: "TF Binding Analysis"
output: html_notebook
---

```{r}
library(ggpubr)

# Source the base code and functions
if (Sys.info()[['sysname']] == "Windows") {
  source("C:/Users/Kevin/Desktop/MacAlpine_KM/Code/get_midpoints.R")
  gene.df = read.csv("D:/HARDAC_scripts/sacCer3_genes_for_making_schematic.csv",header = T)
} else {
  gene.df = read.csv("/Users/kmoyung/MacAlpine_KM/HARDAC_scripts/sacCer3_genes_for_making_schematic.csv", header = T)
  source("/Users/kmoyung/MacAlpine_KM/Code/get_midpoints.R")
}
```

## Load the master dataframe

```{r}
master.df <- readRDS("/Users/kmoyung/Desktop/master_diff_length_corrected.RDS")
```

## Compare annotated binding data with actual TF changes

```{r}
change <- 0.5

sig_chrom.df <- master.df %>%
  filter(sp::point.in.polygon(Microarray_FC, -log10(gamma_pvals), curve_data$x, curve_data$y) == 0 & (abs(Microarray_FC) >= 0.5 | between(Microarray_FC, -0.5, 0.5))) %>%
  mutate(ActualBinding = ifelse(abs(Subnuc_Norm) > change, "Yes", "No"))

# Subset for interactions with binding data (direct) and label if there's a significant binding change
direct_all.df <- sig_chrom.df %>%
  filter(MacIsaac_Site == 1 | Rossi_Site == 1 | hasMotif == 1)

# Plot as a proportional bar chart
long.df <- melt(direct_all.df[, c("Mutant", "Locus", "hasMotif", "MacIsaac_Site", "Rossi_Site", "ActualBinding")]) %>%
  filter(value == 1) %>%
  group_by(variable, ActualBinding) %>%
  mutate(n = sum(value))

compare_direct <- ggplot(long.df, aes(x = variable, group = ActualBinding, fill = ActualBinding)) + 
  geom_bar(stat = "count", color = "black", width = 0.5) + 
  # scale_fill_manual(values = c("Yes" = "#542788", "No" = "#8073AC")) + 
  scale_fill_manual(values = c("Yes" = "#542788", "No" = "grey70")) + 
  theme_classic() + 
  theme(panel.grid = element_blank(), 
        axis.text = element_text(color = "black"), 
        legend.position = "none", 
        legend.background = element_blank(),
        legend.box.background = element_rect(colour = "black")) + 
  guides(fill = guide_legend(nrow = 1, override.aes = list(size = 1))) + 
  labs(x = "", y = "Count", fill = "w/ Chromatin Evidence") + 
  scale_x_discrete(labels = c("Motif", "MacIsaac", "Rossi"))

compare_direct
```

## Get total direct vs. indirect effects

```{r}
compare_direct_indirect <- ggplot(sig_chrom.df, aes(x = ActualBinding, fill = ActualBinding)) + 
  geom_bar(stat = "count", color = "black", width = 0.5) + 
  geom_text(stat = "count", aes(label = ..count..), vjust = -1) + 
  scale_x_discrete(labels = c("Indirect", "Direct w/\nChromatin Evidence")) + 
  scale_fill_manual(values = c("Yes" = "#542788", "No" = "#FDB863")) + 
  theme_classic() + 
  theme(panel.grid = element_blank(), 
        legend.position = "none", 
        axis.text = element_text(color = "black")) + 
  labs(x = "", y = "Count") + 
  coord_cartesian(clip = "off")

compare_direct_indirect
```

## Identify promoters with differential binding (pioneer factors)

Step 1: Identify genes with multiple binding TFs (either Fimo motif OR Rossi) 
Step 2: Subset to include sites where at least ONE mutant has a strong occupancy change

```{r}
change <- -1.4
fillcol <- "#87AEDF"

# Subset for interactions with binding data (direct) and loci with at least one mutant with a strong TF change
direct_master.df <- master.df %>%
  filter(Rossi_Site == 1 | MacIsaac_Site == 1) %>%
  group_by(Locus) %>%
  # filter(n() > 1 & any(Subnuc_Norm < change))
  filter(n() > 1 & any(TF_Change < change))

# Create a column highlighting whether or not it's the key TF
# direct_master.df$KeyTF <- ifelse(direct_master.df$Subnuc_Norm < change, "Pioneer TF", "Other TF")
direct_master.df$KeyTF <- ifelse(direct_master.df$TF_Change < change, "Pioneer TF", "Other TF")

num_tfs <- ggplot(subset(direct_master.df, KeyTF == "Pioneer TF"), aes(x = reorder(stringr::str_to_title(Mutant), Mutant, function(x) + length(x)))) + 
  geom_bar(stat = "count", color = "black", fill = fillcol) + 
  theme_classic() + 
  theme(legend.position = "none",
        plot.title = element_text(face = "bold", color = "black", size = 12),
        axis.text = element_text(color = "black")) + 
  labs(x = "Pioneer TF", y = "# of Promoters") + 
  coord_cartesian(expand = c(0,0)) + 
  coord_flip()

num_tfs
# ggsave(patch, filename = "/Users/kmoyung/Desktop/num_tfs.png", dpi = 800, width = 3, height = 3)
```

## Plot changes in TF occupancy at each locus

```{r}
selected_loci <- c("MET10", "GCV3", "DSE1", "BDF1", "QCR10", "THI74")
# selected_loci <- c("DSE1", "GCV3", "MET10", "QCR10")
# selected_loci2 <- c("DSE1", "BDF1", "QCR10")

# Relabel mutants
direct_master.df$Label <- tolower(paste0(tolower(direct_master.df$Mutant), "∆"))

# occupancy_change <- ggplot(subset(direct_master.df, Locus %in% selected_loci), aes(x = reorder(Label, Subnuc_Norm), y = Subnuc_Norm, fill = KeyTF)) + 
#   theme_bw() + 
#   geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.5) + 
#   geom_segment(aes(y = 0, yend = Subnuc_Norm, xend = reorder(Label, Subnuc_Norm), color = KeyTF), alpha = 1, linewidth = 5) + 
#   geom_point(shape = 21, color = NA, size = 2) + 
#   theme(axis.text.x = element_blank(),
#         axis.ticks.x = element_blank(),
#         axis.text.y = element_text(color = "black", size = 12), 
#         legend.position = "none", 
#         panel.grid = element_blank(), 
#         plot.title = element_text(hjust = 0.5, face = "bold"), 
#         strip.background = element_rect(color = NA, fill = NA), 
#         strip.text = element_text(color = "black", face = "bold.italic")) + 
#   scale_fill_manual(values = c("Pioneer TF" = fillcol, "Other TF" = "grey70")) + 
#   scale_color_manual(values = c("Pioneer TF" = fillcol, "Other TF" = "grey70")) + 
#   labs(x = "", y = "TF Occupancy \u0394") + 
#   geom_text(data = subset(direct_master.df, Locus %in% selected_loci), aes(label = Label, x = reorder(Label, Subnuc_Norm), group = Locus, y = ifelse(Subnuc_Norm > 0, Subnuc_Norm + 0.25, Subnuc_Norm - 0.25)), position = position_dodge2(preserve = "single", width = 0.75), label.size = NA, size = 4, fontface = "italic") +
#   facet_wrap(~ Locus, scales = "free_x")

# Use TF Change, but if it doesn't exist then default to Subnuc_Norm
direct_master.df$TF_LFC <- ifelse(is.na(direct_master.df$TF_Change), direct_master.df$Subnuc_Norm, direct_master.df$TF_Change)

direct_master.df$KeyTF <- ifelse(direct_master.df$TF_LFC < change, "Pioneer TF", "Other TF")

occupancy_change <- ggplot(subset(direct_master.df, Locus %in% selected_loci), aes(x = reorder(Label, TF_LFC), y = TF_LFC, fill = KeyTF)) + 
  theme_bw() + 
  geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.5) + 
  geom_segment(aes(y = 0, yend = TF_LFC, xend = reorder(Label, TF_LFC), color = KeyTF), alpha = 1, linewidth = 5) + 
  geom_point(shape = 21, color = NA, size = 2) + 
  theme(axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.text.y = element_text(color = "black", size = 12), 
        legend.position = "none", 
        panel.grid = element_blank(), 
        plot.title = element_text(hjust = 0.5, face = "bold"), 
        strip.background = element_rect(color = NA, fill = NA), 
        strip.text = element_text(color = "black", face = "bold.italic")) + 
  scale_fill_manual(values = c("Pioneer TF" = fillcol, "Other TF" = "grey70")) + 
  scale_color_manual(values = c("Pioneer TF" = fillcol, "Other TF" = "grey70")) + 
  labs(x = "", y = "TF Occupancy \u0394") + 
  geom_text(data = subset(direct_master.df, Locus %in% selected_loci), aes(label = Label, x = reorder(Label, TF_LFC), group = Locus, y = ifelse(TF_LFC > 0, TF_LFC + 0.25, TF_LFC - 0.25)), position = position_dodge2(preserve = "single", width = 0.75), label.size = NA, size = 4, fontface = "italic") +
  facet_wrap(~ Locus, scales = "free_x")

occupancy_change
# ggsave(patch, filename = "/Users/kmoyung/Desktop/tf_binding_barplot.png", dpi = 800, width = 6, height = 3)
```

## Plot expression fold changes between pioneer and non-pioneer factors

```{r}
library(ggpubr)

pioneer_factors <- unique(subset(direct_master.df, KeyTF == "Pioneer TF")$Mutant)

direct_master.df$KeyTF <- ifelse(direct_master.df$Subnuc_Norm < -0.5, "Pioneer TF", "Other TF")

compare_expression <- ggplot(subset(direct_master.df, !is.na(KeyTF)), aes(x = KeyTF, y = Microarray_FC, color = KeyTF)) + 
  geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.5) + 
  geom_boxplot(outlier.shape = NA, width = 0.3, staplewidth = 0.5, position= position_nudge(x=0.5)) +
  geom_jitter(shape = 16, position = position_jitter(0.1), aes(color = KeyTF), alpha = 0.6, size = 1) +
  theme_classic() + 
  stat_compare_means(method = "t.test", label.y = -3.5, aes(label = paste0("p < ", after_stat(p.format)))) + 
  labs(x = "", y = expression(Expression~log[2]~FC)) + 
  theme(legend.position = "none", 
        axis.text.y = element_text(color = "black", size = 12), 
        axis.text.x = element_text(color = "black")) + 
  scale_fill_manual(values = c("Pioneer TF" = fillcol, "Other TF" = "grey70")) + 
  scale_color_manual(values = c("Pioneer TF" = fillcol, "Other TF" = "grey70"))
  
compare_expression

compare_tf_occ <- ggplot(direct_master.df, aes(x = KeyTF, y = TF_Change, color = KeyTF)) + 
  geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.5) + 
  geom_boxplot(outlier.shape = NA, width = 0.5, position= position_nudge(x=-.5)) +
  geom_jitter(shape = 16, position = position_jitter(0.2), aes(color = KeyTF), alpha = 0.8, size = 2) +
  theme_classic() + 
  stat_compare_means(method = "wilcox.test", label.y = -3) + 
  labs(x = "", y = "TF Occupancy (log2FC)") + 
  theme(legend.position = "none", 
        axis.text = element_text(color = "black")) + 
  scale_fill_manual(values = c("Pioneer TF" = fillcol, "Other TF" = "grey70")) + 
  scale_color_manual(values = c("Pioneer TF" = fillcol, "Other TF" = "grey70"))

# compare_expression + compare_tf_occ
# ggsave(compare_expression, filename = "/Users/kmoyung/Desktop/compare_expression.png", dpi = 800, width = 3, height = 3)
```

## Save direct/indirect TF binding data as supplemental figure

```{r}
supp_tf_binding <- compare_direct_indirect + compare_direct + plot_annotation(tag_levels = "A") & theme(plot.tag = element_text(face = "bold"))
supp_tf_binding

ggsave(supp_tf_binding, filename = "/Users/kmoyung/supp_figures/direct_indirect.png", dpi = 800)
```

## Patch figures

```{r}
# design <- "ACDE \n BFFF \n #FFF"
design <- "AEEE \n AEEE \n BEEE"
# design <- "ACCC \n BCCC"

tf_figure <-  num_tfs + compare_expression + occupancy_change + plot_layout(design = design)
tf_figure

# ggsave(tf_figure, filename = "/Users/kmoyung/Desktop/binding_pioneer_tf_figure.png", dpi = 800, width = 13, height = 10)
ggsave(tf_figure, filename = "/Users/kmoyung/Desktop/binding_pioneer_tf_figure_new.png", dpi = 800, width = 8, height = 6)
```

## Check for TF binding sites with a GAIN of TF as a result of deletion

```{r}
change <- 1
fillcol <- "#DA95CC"

# Subset for interactions with binding data (direct)
direct_master.df <- master.df %>%
  filter(Rossi_Site == 1) %>%
  # filter(hasMotif == 1) %>%
  group_by(Locus) %>%
  filter(n() > 1 & any(Subnuc_Norm > change))

# Create a column highlighting whether or not it's the key TF
direct_master.df$KeyTF <- ifelse(direct_master.df$Subnuc_Norm > change, "Repressor TF", "Other TF")

# num_sites <- ggplot(direct_master.df, aes(x = reorder(Locus, Locus, function(x) - length(x)))) + 
#   geom_bar(stat = "count", color = "black") + 
#   theme_classic() + 
#   theme(legend.position = "none", 
#         axis.text.x = element_text(angle = 70, vjust = 1, hjust = 1)) + 
#   labs(x = "Gene", y = "Count", title = "# of TFs bound to a promoter")

num_tfs <- ggplot(subset(direct_master.df, KeyTF == "Repressor TF"), aes(x = reorder(stringr::str_to_title(Mutant), Mutant, function(x) + length(x)))) + 
  geom_bar(stat = "count", color = "black", fill = fillcol) + 
  theme_classic() + 
  theme(legend.position = "none",
        plot.title = element_text(face = "bold", color = "black", size = 12),
        axis.text = element_text(color = "black")) + 
  labs(x = "TF", y = "# of Binding Sites", title = "Inhibitory TF Activity") + 
  coord_cartesian(expand = c(0,0)) + 
  coord_flip()

patch <- num_tfs
patch
ggsave(patch, filename = "/Users/kmoyung/Desktop/num_tfs_gain.png", dpi = 800, width = 3, height = 2)
```
```{r}
selected_loci <- c("ANB1", "CIT2", "HTA1", "MRS4", "SER3", "YLR339C")
patch <- ggplot(subset(direct_master.df, Locus %in% selected_loci), aes(x = Locus, y = Subnuc_Norm, fill = KeyTF)) + 
  theme_bw() + 
  geom_bar(stat = "identity", position = position_dodge2(preserve = "single"), color = "black", width = 0.75) + 
  facet_wrap(~ Locus, scales = "free_x", ncol = 3, nrow = 2) + 
  theme(axis.text.x = element_blank(), 
        axis.ticks.x = element_blank(),
        axis.text.y = element_text(color = "black"), 
        legend.position = "none") + 
  scale_fill_manual(values = c("Repressor TF" = fillcol, "Other TF" = "grey70")) + 
  labs(x = "", y = "TF Occupancy \u0394") + 
  geom_label_repel(data = subset(direct_master.df, Locus %in% selected_loci & Subnuc_Norm > change), aes(label = paste0(tolower(Mutant), "\u0394")), label.padding = 0.1, 
                   show.legend = F, 
                   label.size = NA, 
                   fill = NA, 
                   force_pull = 1, 
                   fontface = "italic", 
                   size = 4, 
                   nudge_x = -0.5, 
                   arrow = arrow(length = unit(0.15, 'cm')), 
                   arrow.fill = "black")

patch
ggsave(patch, filename = "/Users/kmoyung/Desktop/tf_binding_barplot_gain.png", dpi = 800, width = 6, height = 3)
```

## Plot expression fold changes between gain and non-gain TFs

```{r}
compare_expression <- ggplot(direct_master.df, aes(x = KeyTF, y = Microarray_FC, color = KeyTF)) + 
  geom_boxplot(outlier.shape = NA, width = 0.5) +
  geom_jitter(shape = 16, position = position_jitter(0.2), aes(color = KeyTF), alpha = 0.8, size = 2) +
  theme_classic() + 
  stat_compare_means(method = "wilcox.test", label.y = -3) + 
  labs(x = "", y = "Expression (log2FC)") + 
  theme(legend.position = "none") + 
  scale_fill_manual(values = c("Repressor TF" = fillcol, "Other TF" = "grey70")) + 
  scale_color_manual(values = c("Repressor TF" = fillcol, "Other TF" = "grey70"))
                       
compare_expression
ggsave(compare_expression, filename = "/Users/kmoyung/Desktop/compare_expression_gain.png", dpi = 800, width = 3, height = 3)
```