nzs_analysis_meanDist.Rmd

---
title: "Untitled"
author: "David Corujo"
date: "2023-10-03"
output: html_document
editor_options:
  chunk_output_type: console
---

# Load packages
```{r}
library(regioneReloaded)
library(ggplot2)
library(patchwork)
library(reshape2)
library(tidyverse)
```

# Plot nzs and zscore for ALL comparisons

## Function: individual plots
```{r}
funPlot <- function(full_df, Bname, xlabs = xlabs) {
ggplot(full_df[full_df$B == Bname,], aes(x = fraction, y = score, color = type, group = type)) +
  facet_grid(type~B, scale = "free") +
  geom_point() +
  geom_line() +
  scale_y_continuous(limits = function(x){c(min(0, x) * 1.1, max(0, x) * 1.1)}) +
  scale_color_manual(values = c("green4", "darkorange")) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) +
  scale_x_discrete(labels = xlabs)
}
```

## Function: create list of plots, wrap and save
```{r}
wrapPlots <- function(r1, permList_zs, permList_nzs, ran_title = "RaR") {
  
  # Obtain permTest from list
  cwPT <- permList_zs[[r1]]
  cwPTn <- permList_nzs[[r1]]
  
  # Obtain n regions for plotting later
  xlabs <- unlist(lapply(cwPT@multiOverlaps, FUN = function (x) unique(x$n_regionA)))

  # Extract matrices of ZS and nZS
  matX<-t(cwPT@matrix$GMat)
  matX<-as.data.frame(matX[order(rownames(matX)), order(colnames(matX))])
  matX

  matXn<-t(cwPTn@matrix$GMat)
  matXn<-as.data.frame(matXn[order(rownames(matXn)), order(colnames(matXn))])
  matXn

  # Conver to dataframe to prepare for tidyr funcions
  df <- data.frame(matX)
  df$fraction <- rownames(df)
  
  df2 <- data.frame(matXn)
  df2$fraction <- rownames(df2)
  
  # Melt ZS and nZS dfs into long format
  meltmat <- tidyr::pivot_longer(df, cols = !fraction, names_to = "B", values_to = "score")
  meltmat$type <- "zscore"
  
  meltmat2 <- tidyr::pivot_longer(df2, cols = !fraction, names_to = "B", values_to = "score")
  meltmat2$type <- "n_zscore"
  
  # Join into a single df
  full_df <- full_join(meltmat, 
            meltmat2)
  
  # Reorder factors
  full_df$type <- factor(full_df$type, levels = c("zscore", "n_zscore"))
  
  all_B <- unique(full_df$B)
  
  # Generate a list of plots using plotting function
  plotList <- lapply(all_B, FUN = funPlot, full_df = full_df,  xlabs = xlabs)
  return(plotList)
}
```

# Check results from HPC run and load region sets
```{r}
regSets <- readRDS("regionSets/ENCODE_filterScore.RDS")

permRes_RaR <- readRDS("hpcResults/permResList_RaR_meanDist.RDS")
permRes_ReR <- readRDS("hpcResults/permResList_ReR_meanDist.RDS")
permRes_ReG <- readRDS("hpcResults/permResList_ReG_meanDist.RDS")
```

# Generate matrix and plot for all tests

```{r}
permList_zs_RaR <- lapply(permRes_RaR, FUN = makeCrosswiseMatrix, symm_matrix = F, zs.type = "z_score")
permList_nzs_RaR <- lapply(permRes_RaR, FUN = makeCrosswiseMatrix, symm_matrix = F)

permList_zs_ReR <- lapply(permRes_ReR, FUN = makeCrosswiseMatrix, symm_matrix = F, zs.type = "z_score")
permList_nzs_ReR <- lapply(permRes_ReR, FUN = makeCrosswiseMatrix, symm_matrix = F)

permList_zs_ReG <- lapply(permRes_ReG, FUN = makeCrosswiseMatrix, symm_matrix = F, zs.type = "z_score")
permList_nzs_ReG <- lapply(permRes_ReG, FUN = makeCrosswiseMatrix, symm_matrix = F)
```

## Run lapply

```{r}
setNames <- c("MAFF_ENCFF005YUC", "FOXA1_ENCFF011QFM", "RAD21_ENCFF155CEQ", "POLR2A_ENCFF159PYD", 
              "CTCF_ENCFF199YFA", "H3K9me3_ENCFF372HCL", "H3K27ac_ENCFF392KDI", "H3K4me3_ENCFF982DUT")

lapply(setNames, function(x) {wrapPlots(r1 = x, permList_zs_RaR, permList_nzs_RaR, ran_title = "RaR_meanDist")})
lapply(setNames, function(x) {wrapPlots(r1 = x, permList_zs_ReR, permList_nzs_ReR, ran_title = "ReR_meanDist")})
lapply(setNames, function(x) {wrapPlots(r1 = x, permList_zs_ReG, permList_nzs_ReG, ran_title = "ReG_meanDist")})
```

# Plot nzs and zscore for a specific comparison

```{r}
r1 <- "H3K4me3_ENCFF982DUT"

# Obtain permTest from list
cwPT <- permList_zs_RaR[[r1]]
cwPTn <- permList_nzs_RaR[[r1]]

# Obtain n regions for plotting later
xlabs <- unlist(lapply(cwPT@multiOverlaps, FUN = function (x) unique(x$n_regionA)))

# Extract matrices of ZS and nZS
matX<-t(cwPT@matrix$GMat)
matX<-as.data.frame(matX[order(rownames(matX)), order(colnames(matX))])
matX

matXn<-t(cwPTn@matrix$GMat)
matXn<-as.data.frame(matXn[order(rownames(matXn)), order(colnames(matXn))])
matXn

# Conver to dataframe to prepare for tidyr funcions
df <- data.frame(matX)
df$fraction <- rownames(df)

df2 <- data.frame(matXn)
df2$fraction <- rownames(df2)

# Melt ZS and nZS dfs into long format
meltmat <- tidyr::pivot_longer(df, cols = !fraction, names_to = "B", values_to = "score")
meltmat$type <- "zscore"

meltmat2 <- tidyr::pivot_longer(df2, cols = !fraction, names_to = "B", values_to = "score")
meltmat2$type <- "n_zscore"

# Join into a single df
full_df <- full_join(meltmat, 
          meltmat2)

# Reorder factors
full_df$type <- factor(full_df$type, levels = c("zscore", "n_zscore"))

all_B <- unique(full_df$B)
```

## Function to plot selected comparisons
```{r}
# Function to plot selected comparisons

paperPlot <- function(r2, full_df) {
  p1 <- full_df %>%
  filter(B == r2) %>%
  ggplot(aes(x = fraction, y = score, color = type, group = type)) +
  facet_grid(type~B, scale = "free") +
  geom_point() +
  geom_line() +
  scale_y_continuous(limits = function(x){c(min(0, x) * 1.1, max(0, x) * 1.5)}) +
  scale_color_manual(values = c("green4", "darkorange")) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
        legend.position = "none") +
  scale_x_discrete(labels = xlabs) +
  labs(subtitle = r1)
  return(p1)
}
```
# Generate plotList of selected comparisons and save
```{r}
# Generate plotlist of selected comparisons
# Check names
unique(full_df$B)

setNames <- c("POLR2A_ENCFF159PYD", "POLR2A_ENCFF354VWZ", 
              "POLR2G_ENCFF485SII", "POLR2G_ENCFF551IJP",
              "H3K27ac_ENCFF392KDI", "H3K9me3_ENCFF372HCL")

plotList <- lapply(setNames, FUN = function(x) paperPlot(r2 = x, full_df = full_df))

wrap_plots(plotList, ncol = 2)

saveRDS(plotList, file = "Test_zscore/plots_nzs_H3K4me3_RaR_meandist.RDS")
```