---
title: "MNase-seq for MicroC"
author: "Elisa Oberbeckmann"
date: "10/17/2022"
output: html_document
---
```{r import bam data, echo=FALSE}
# generate coverage files by using bam-files
# plot fragment length distribution
setwd("/Users/Elisa/Documents/For Manuscripts/Quililan, Oberbeckmann et al./Plots/bam-files/")
# Coverage file - paired-end #
rm(list=ls())
library(GenomicAlignments)
library(rtracklayer)
library(zoo)
list.files(".",".bam") -> all.files

for (file in all.files) {
  gr <- readGAlignmentPairs(file)
  gr <- granges(gr)
  flt <- gr@ranges@width>124 & gr@ranges@width<206
  gra <- gr[flt,]
  gra <- resize(gra, 50, fix = "center")
  cov <- coverage(gra)
  gd <- as(cov, "GRanges")
  
  gran.data <- gd
  chr.gran <- split(gran.data, seqnames(gran.data))
  res.smooth.lfc <- NULL
  for (chr in seqnames(gran.data)@values) {

  temp.metdat <- elementMetadata(chr.gran[[chr]])[,1]
    
  smooth.lfc <- rollapply(temp.metdat, width = 20, FUN = mean,
    align = "center", partial = T) 

  res.smooth.lfc <- append(res.smooth.lfc, smooth.lfc)
}
  gran.data$score <- res.smooth.lfc
  gran.data -> gd
  
  export.bw(gd, gsub(".bam","_cov.bw",file))
  
    pdf(gsub(".bam",".pdf",file), 6,4) 
  h <- hist(gr@ranges@width[!gr@ranges@width>550], breaks = 100 , xlim =c(50,550), xlab = "fragment length / bp", main = file, cex.lab = 1.5, cex.axis = 1.5)
  i <- h$breaks[which.max(h$counts)]
  text(i+100, y = max(h$counts), labels = paste(i, "bp"))
  abline(v = 125,                      
      col = "red",
      lwd = 3)
  abline(v = 205,                       
       col = "red",
        lwd = 3)
  dev.off() 
} 

```

```{r alignment, echo=TRUE}
rm(list=ls())
library(zoo)
library(rtracklayer)

overview <- read.csv("~/Documents/For Manuscripts/Quililan, Oberbeckmann et al./Plots/bam-files/MCC-Kim.csv")# document with sample conditions
load("Alignments_for_pGPV_X.rds") # table with ChIP-exo peaks
pGPV_X[5:6] -> pGPV_X

overview -> sc.files
coverageWindowsStranded <- function(centers, window.size, coverage) {
  
  centers <- centers[centers$seqnames %in% names(coverage),]
  
  result <- sapply(names(coverage), function(x) {
    my.cov <- coverage[[x]]
    my.peaks <- centers[centers$seqnames==x,]
    mw.views <- Views(my.cov, start=my.peaks$start+4-ceiling(window.size/2), width=window.size+1)
    ## remove out-of bounds views
    flt <- start(mw.views)>0 & end(mw.views) < length(my.cov)
    mw.views <- mw.views[flt,]
    my.peaks <- my.peaks[flt,]
    if (length(mw.views) > 0) {
      mat <- as.matrix(mw.views)
      colnames(mat) <- seq(from=(0-ceiling(window.size/2)), to=0+ceiling(window.size/2))
      rownames(mat) <- rownames(my.peaks)
      return(mat)
    } else {
      return(NULL)
    }
  })
  mat <- Reduce(rbind, result)
  centers <- centers[rownames(centers) %in% rownames(mat),]
  na.omit(match(rownames(centers), rownames(mat)) )-> o
  centers <- centers[o,]
  mat[centers$strand=="-",] <- t(apply(mat[centers$strand=="-",],1,rev))
  colnames(mat) <- seq(-ceiling(win.size/2),ceiling(win.size/2))
  rownames(mat) <- centers$name
  mat
}

norm.square <- function(mat) {t(apply(mat, 1, function(x){x/(sum(x))}))}
win.size <- 3000
all.tidy <- c()
df2 <- data.frame(matrix(NA, nrow = 3001, ncol = 0))
for (x in names(pGPV_X)) {
  pGPV_X[[x]] -> centers
  df1 <- c()
  for (i in 1:nrow(sc.files)) {
    import.bw(sc.files$source.file[i]) -> gr
    coverage(gr, weight = "score") -> cov
    names(cov) <- paste0("chr",sub("chr","", names(cov)))
    
    mat <- coverageWindowsStranded(centers, win.size, cov)
    mat.norm <- norm.square(mat)
    vec <- colMeans(mat.norm, na.rm = T)
    vec/mean(vec) -> vec
    df1 <- cbind(df1,vec)
  }
  as.data.frame(df1) -> df
  colnames(df) <- paste(overview$name, x)
  df2 <- cbind(df2, df)
  df_tidy <- c()
  for (i in 1:ncol(df)){
    df_tidy <- rbind(df_tidy, data.frame(pos=as.numeric(names(vec)), occ = df[,i], 
                                         barrier = rep(overview$barrier[i], nrow(df)), 
                                         remodeler = rep(overview$remodeler[i], nrow(df)), 
                                         replicate = rep(overview$replicate[i], nrow(df)), 
                                         density = rep(overview$density[i], nrow(df)),
                                         alignment = rep(x, nrow(df))))
  }
  all.tidy <- rbind(all.tidy, df_tidy)
}
save(df2, file = "Table_peak_calling.rds") # used for peak calling
save(all.tidy, file = "Tidy_table_composite_Reb1-Abf1.rds") # used for plots
```


## ggplot figures

```{r ggplot individual replicates, echo = T}
rm(list=ls())
library(tidyverse)
load("Tidy_table_composite_Reb1-Abf1.rds")
col <- c("darkblue",  "red3") 


all.tidy[all.tidy$density == "80",] -> data
data$remodeler = factor(data$remodeler, levels=c("in vivo",  "INO80", "Chd1", "ISW2",  "RSC","no remodeler"))
data$alignment = factor(data$alignment, levels=c( "Abf1-ChIPexo", "Reb1-ChIPexo"))
data$barrier = factor(data$barrier, levels=c( "none", "Abf1-Reb1", "in vivo"))
data$replicate = factor(data$replicate, levels=c( "1", "2")) 
na.omit(data) -> data


p <- ggplot(data, aes(x=pos, y=occ, color = replicate, group =replicate )) +

  stat_summary(fun = 'mean', geom = 'line',  size = 0.8, alpha = 0.9) +
  facet_grid(alignment ~remodeler) +
  scale_color_manual("replicate", values = col) +
  scale_y_continuous("nucleosome density") +
  scale_x_continuous("distance to indicated site / bp", limits = c(-1500,1500)) +
  theme_bw()+
  theme(strip.background = element_blank(),legend.position = "top", 
        strip.text.x = element_text(size=8, face = "bold"),
        strip.text.y = element_text(size=8, face = "bold"))

ggsave("Extended-data-Fig7.pdf", plot = p, width = 29, height = 12, units = "cm")
```

```{r ggplot merges Abf1, echo = T}
rm(list=ls())
library(tidyverse)
library(viridis)
load("Tidy_table_composite_Reb1-Abf1.rds") 
col <- c("darkblue", "dodgerblue", "red2", "black")


all.tidy[all.tidy$density == "80",] -> data
data$replicate = factor(data$replicate, levels=c( "1", "2")) 
data$remodeler = factor(data$remodeler, levels=c("in vivo", "INO80", "ISW2", "Chd1", "RSC", "no remodeler"))
data$alignment = factor(data$alignment, levels=c( "Abf1-ChIPexo"))
na.omit(data) -> data


p <- ggplot(data, aes(x=pos, y=occ)) +
  stat_summary(fun = 'mean', geom = 'line',  size = 0.7, alpha = 0.8) +
  facet_grid(.~ remodeler) +
  scale_color_manual("replicate", values = col) +
  scale_y_continuous("nucleosome density") + 
  scale_x_continuous("distance to Reb1 site / bp", limits = c(-1500,1500)) +
  theme_bw()+
  theme(strip.background = element_blank(),legend.position = "top", 
        strip.text.x = element_text(size=8, face = "bold"),
        strip.text.y = element_text(size=8, face = "bold"))


ggsave("Fig3_Abf1.pdf", plot = p, width = 28, height = 6.15, units = "cm")
```

## peak calling

```{r peak calling and distance calculation, echo = T}
rm(list=ls())
library(plyr)
library(tidyverse)
library(ggplot2)
library(RColorBrewer)
library(Hmisc)

find_peaks <- function (x, m = 3){
  shape <- diff(sign(diff(x, na.pad = FALSE)))
  pks <- sapply(which(shape < 0), FUN = function(i){
    z <- i - m + 1
    z <- ifelse(z > 0, z, 1)
    w <- i + m + 1
    w <- ifelse(w < length(x), w, length(x))
    if(all(x[c(z : i, (i + 2) : w)] <= x[i + 1])) return(i + 1) else return(numeric(0))
  })
  pks <- unlist(pks)
  pks
}

overview <- read.csv("~/Documents/For Manuscripts/Quililan, Oberbeckmann et al./Plots/bam-files/MCC-Kim.csv")
overview[!str_detect(overview$name, pattern="no remodeler"),] -> overview
overview[!str_detect(overview$name, pattern="none"),] -> overview

load("Table_peak_calling.rds")
df2[,!str_detect(colnames(df2), pattern="no remodeler")] -> df2
df2[,str_detect(colnames(df2), pattern=" Abf1-ChIP")] -> df1
df2[,str_detect(colnames(df2), pattern=" Reb1-ChIP")] -> df2

all <- list(df1, df2)
names(all) <- c("Abf1", "Reb1")

ticks <- c("-1500", "-750", "0", "750", "1500")
peak.table <- c(1:8)

pdf("All_peaks_MNase-seq_exo.pdf")  # Plot aligned MNase-seq data and call peaks (red dots)
par(mfrow = c(2, 2))

for (barrier in names(all)){
  all[[barrier]] -> df
 for (i in 1:ncol(df)) {
  df[,i] -> w
  names(w) <- rownames(df)
  p <- find_peaks(w, m = 50)
  p[!between(as.numeric(names(p)), -100, 100)] -> p
  p[!between(as.numeric(names(p)), -1500, -779)]-> p
  p[!between(as.numeric(names(p)), 747, 1500)] -> p
  peak.table <- rbind(peak.table, p)
  #ind <- rep(1, length(w))
  #ind[p] <- 2
  colnames(df2)[i] -> titel
  plot(names(w), w, type = 'l', main = titel, cex.main = 0.7, ylab = "dyad density", xlab = paste0("distance to ", barrier,  " / bp"), xaxt = "n", cex.axis = 0.8, xlim=c(-1500,1500), ylim=c(0.2,2))
       
  points(names(p), w[p], col = 'red', pch = 20)
  axis(1, at=ticks,labels=ticks, cex.axis = 0.8)
  text(as.numeric(names(p)) + 10, w[p]*1.1, labels = names(p), cex = 0.7)

  } 
}
dev.off()

for(i in 1:nrow(peak.table)){
    peak.table[i,][duplicated(peak.table[i,])] <- NA 
} # remove duplicates and replace with NA

peak.table <- peak.table[-1,]  - 1500 
rownames(peak.table) <- c(colnames(df1), colnames(df2))
#apply(peak.table, 2, function(x) ifelse(between(x,-20, 20), NA, x)) -> peak.table

all.peaks.pos <- c(1:4)
all.peaks.neg <- c(1:4)
for (n in 1:nrow(peak.table)) {
  pos <- na.omit(peak.table[n,][na.omit(peak.table[n,]) > 0])
  all.peaks.pos <- rbind(all.peaks.pos, pos)
  neg <- na.omit(peak.table[n,][peak.table[n,] < 0])
  all.peaks.neg <- rbind(all.peaks.neg, neg)
}

all.peaks.n <- c(1:4)
for (i in 1:nrow(all.peaks.neg)){
  all.peaks.neg[i,][duplicated(all.peaks.neg[i,])] <- NA 
  all <- sort(abs(all.peaks.neg[i,]))
  all.peaks.n <- rbind(all.peaks.n, all)
}
for(i in 1:nrow(all.peaks.n)){
  all.peaks.n[i,][duplicated(all.peaks.n[i,])] <- NA 
}
all.peaks.n <- all.peaks.n[-1:-2,]

for(i in 1:nrow(all.peaks.pos)){
  all.peaks.pos[i,][duplicated(all.peaks.pos[i,])] <- NA 
}
all.peaks.pos <- all.peaks.pos[-1,]

table <- list(all.peaks.n, all.peaks.pos)
names(table) <- c("negative", "positive")
NRL <- list()

for (a in 1:length(table)) {
  werte <- table[[a]]
  rownames(werte) <- c(colnames(df1),colnames(df2))
  werte <- werte[,1:4]
  colnames(werte) <- 1:4
  dist.all <- c()
  DtB <- c()
  for (w in 1:nrow(werte)){ 
      x <- werte[w,]
      x[1] -> value
      DtB <- rbind(DtB,value)
      dist <- c()
    for(i in 1:length(x)){
        er <- i+1
        sr <- er-1
        dist <-  append(dist, x[er] - x[sr])
    }
      as.vector(dist) -> dist
      dist.all <- rbind(dist.all, dist)
  }
  rownames(dist.all) <- c(colnames(df1), colnames(df2))
  dist.all <- dist.all  - 147
  DtB <- DtB - 73
  #apply(dist.all, 2, function(x) ifelse(x < 8, NA, x)) -> dist.all
  apply(dist.all, 2, function(x) ifelse(x > 70, NA, x)) -> dist.all
  dist.all <- cbind(dist.all, DtB)
  NRL[[names(table)[a]]] <- dist.all
}
colnames(NRL$negative) <- c(-1:-4, "-DtB")
colnames(NRL$positive) <- c(1:4, "DtB")

overview[overview$name %in% gsub(" Abf1-ChIPexo", "", colnames(df1)),] -> overview
overview <- rbind(overview, overview)
overview$alignment <- c(rep("Abf1", nrow(overview)/2), rep("Reb1", nrow(overview)/2)) 
overview.all <- cbind(overview, NRL[[1]][,c(1:3,5)], NRL[[2]][,c(1:3,5)])

save(overview.all, file = "Overview_with_distances.rds")
overview.all -> assembly

### all linker ####
assembly$`3`[assembly$remodeler == 'RSC'] <- NA
assembly$`-3`[assembly$remodeler == 'RSC'] <- NA
assembly$`2`[assembly$remodeler == 'RSC'] <- NA
assembly$`-2`[assembly$remodeler == 'RSC'] <- NA
assembly$`1`[assembly$remodeler == 'RSC'] <- NA
assembly$`-1`[assembly$remodeler == 'RSC'] <- NA
assembly$`3`[assembly$`3`< 2] <- NA
assembly$`-3`[assembly$`-3`< 2] <- NA
assembly$`DtB`[assembly$remodeler == 'RSC'& assembly$barrier == 'Abf1'] <- NA

len.all <- length(assembly$`1`)
linker1 <- rep("linker 1", len.all)
linker2 <- rep("linker 2", len.all)
linker3 <- rep("linker 3", len.all)
dtb <- ifelse(assembly$barrier == 'Reb1', 'dist. to Reb1', 'dist. to Abf1')
dir <- c(rep("upstream", len.all), rep("downstream", len.all))
rem <- as.character(assembly$remodeler)
b <- as.character(assembly$barrier)
den<- assembly$density
rep <- assembly$replicate
al <- assembly$alignment


linker <- c(assembly$`1`, assembly$`-1`, assembly$`2`, assembly$`-2`, assembly$`3`, assembly$`-3`, assembly$`-DtB`, assembly$DtB )
direction <- rep(dir, 4)
linker.type <- c(linker1,linker1, linker2, linker2, linker3, linker3, dtb, dtb)
remodeler <- rep(rem, 8)
barrier <- rep(b, 8)
density<- rep(den, 8)
replicate <- rep(rep, 8)
alignment <- rep(al, 8)

tidy.table <- data.frame(remodeler, barrier, replicate, linker.type, direction, linker, density, alignment)
save(tidy.table,file="Tidy_table_all_linker.rds")
```

```{r ggplot}
rm(list=ls())
library(tidyverse)
library(ggplot2)
library(RColorBrewer)
library(Hmisc)
library(ggpubr)
library(viridis)

load("Tidy_table_all_linker_Quililan_exo.rds")
tidy.table[tidy.table$density == 80,] -> tidy.table
tidy.table[tidy.table$replicate == 1 | tidy.table$replicate == 2,] -> tidy.table
#tidy.table[!(tidy.table$linker.type == "dist. to Reb1" & tidy.table$linker.type == "dist. to Abf1"),] -> td
#write.csv(td, file = "Source_data_Fig4.csv")


###some variables###
col <- c('#ca0020','#000071')

dodge.posn <- position_jitterdodge(jitter.width = 0.3, dodge.width = 0.6)
n_fun <- function(x){
  return(data.frame(y = mean(x), label = paste0(as.integer(mean(x)))))
}


tidy.table -> data
data$linker.type[data$linker.type == 'dist. to Reb1'] <- 'dist. to GRF'
data$linker.type[data$linker.type == 'dist. to Abf1'] <- 'dist. to GRF'
data$remodeler = factor(data$remodeler, levels=c("in vivo","RSC",'INO80','ISW2','Chd1'))
data$replicate = factor(data$replicate, levels=c('1','2'))


p <- ggplot(data, aes(x = factor(remodeler), y = linker, fill = barrier, color = barrier)) +
  geom_dotplot(binaxis = "y", binwidth = 1.2, dotsize = 8, stackdir = "center", alpha = 0.7) +
  scale_y_continuous("length / bp") + 
  scale_x_discrete("remodeler") +
  scale_color_manual("GRF", values = col) +
  scale_fill_manual("GRF",values = col) + 
  theme_bw() +
  facet_grid(.~linker.type ) +
  theme(strip.background = element_blank(), plot.title = element_text(hjust = 0.5), legend.position = "top", 
        strip.text.x = element_text(size=10),
        axis.text.x = element_text(angle = 45, hjust=0.9), 
        strip.text.y = element_text(size=10))
ggsave("Length-distribution_all_remodeler.pdf", plot = p, width = 15, height = 7, units = "cm")

##########################################
tidy.table -> data
data[-which(data$remodeler=='Chd1' & data$barrier=="Reb1" & data$linker.type=="linker 2" & data$replicate == 1 & data$direction == "downstream"),] -> data
data$linker.type[data$linker.type == 'dist. to Reb1'] <- 'dist. to GRF'
data$linker.type[data$linker.type == 'dist. to Abf1'] <- 'dist. to GRF'
data$remodeler = factor(data$remodeler, levels=c("in vivo","RSC",'INO80','ISW2','Chd1')) 
data$replicate = factor(data$replicate, levels=c('1','2')) 
data$density = factor(data$density, levels=c('80')) 
na.omit(data) -> data

dodge.posn <- position_jitterdodge(jitter.width = 0.3, dodge.width = 0.6)
n_fun <- function(x){
  return(data.frame(y = mean(x), label = paste0(as.integer(mean(x)))))
}

linker.distribution <- ggplot(data, aes(x = remodeler, y = linker)) +
  geom_dotplot(binaxis="y", binwidth = 1.2, stackdir = "center", dotsize = 3) +
  scale_y_continuous("length / bp") +
  scale_x_discrete("parameter") +
  theme_bw() +
  stat_summary(fun = mean, geom = "point", size = 5, shape = 15, show.legend = F,  color = "gray") +
  stat_summary(fun.data = n_fun, geom = "text", size = 3) +
  facet_grid( .~linker.type ) +
  theme(strip.background = element_blank(), plot.title = element_text(hjust = 0.5), 
        axis.text.x = element_text( hjust = 1, angle = 45), 
        strip.text.y = element_text(size=10, angle = -90),
        strip.text.x = element_text(size=10, angle = 0)) #, face = "bold"

ggsave("Length-values_all_remodeler.pdf", plot = linker.distribution, width = 15, height = 7, units = "cm")
```
