---
title: "ChIP-FIRE correlation in variable FIRE region"
output: 
  pdf_document: 
    number_sections: yes
---

Here we try to identify variable bins for individuals or populations using FIRE, and see how variation correlate with other epigenetics marks. I used signals from peak alone, normalized by subtracting input, and quantile normalized at bin level. 

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo=FALSE, warning=FALSE, message=FALSE, fig.align='center', fig.height = 10, fig.width = 10)
```

```{r load_library}
library("data.table")
library("qvalue")
library("ggplot2")
library("pander")
library("reshape2")
library("pheatmap")
library("RColorBrewer")
library("limma")

options(scipen = 999)
```

#Identification of variable bins
```{r}
#here we read previous data instead.
bin_sig <- read.table("FIRE_LCL20_40Kb.sig_bin.txt", head=F)
```

We first filter out bins whose FIRE are not significant in any individuals. There are `r nrow(bin_sig)` significant bins for further test.

##limma test
```{r limma_test}
indv_test <- read.table("indv_test_limma_FIRE_LCL20_40Kb.tsv", head=T)
```

There are `r nrow(indv_test)` bins out of total `r nrow(bin_sig)` bins are significant under adjusted p value 0.1.

#Correlation between FIRE and ChIP-seq signal

We use bin-wise shuffle as control.

```{r corr_load_bin, results="hide", cache=TRUE}
corr <- read.table("FIRE_ChIP_LCL20_40Kb.cor", head=F)
colnames(corr) <- c("bin", "scc", "mark")
corr$status <- "Varaible"

rand_corr <- data.frame(fread("FIRE_ChIP_LCL20_40Kb.rand.cor", stringsAsFactors=TRUE))
colnames(rand_corr) <- c("bin", "status", "scc", "mark")
rand_corr <- rand_corr[c(1, 3, 4, 2)]

corr <- rbind(rand_corr, corr)

corr$rand <- ifelse( grepl("rand", corr$status), "Random", "Observed" )
```

```{r}
corr_median <- read.table("FIRE_ChIP_LCL20_40Kb.rand.median", head=F)
colnames(corr_median) <- c("scc", "absscc", "mark")
```

## Raw Signal
```{r corr_compare_bin_raw, cache=TRUE}
pval <- sapply(
  split(
    corr[corr$rand=="Observed",], corr[corr$rand=="Observed", "mark"]
    ), 
  function(x) paste(
    "P=", format.pval(
      sum(median(x$scc)<corr_median[corr_median$mark==x$mark[1], "scc"]) / nrow(corr_median[corr_median$mark==x$mark[1],]), eps=1/nrow(corr_median[corr_median$mark==x$mark[1],])
      )
    )
  )
pval <-  data.frame(mark=names(pval), pval=pval)

ggplot( corr, aes(scc) ) + 
  geom_line(stat="density", alpha=0.5, aes(color=status, linetype=rand)) + 
  facet_wrap(~mark,scales="free_y") + 
  scale_color_manual(values = ifelse( grepl("rand", levels(corr$status)), "gray", "red"), guide = "none") +
  geom_text(data = pval,aes(x = 0.5, y = 0.1, label = pval)) + 
  theme_bw() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) + 
  theme(aspect.ratio = 1) + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.title=element_blank()) + 
  theme(text = element_text(size=12)) +
  ylab("Density") + 
  xlab("Spearman Correlation")
```

```{r}
save.image("corr_FIRE_ChIP.RData")
```

## Absolute Signal
```{r corr_compare_bin, cache=TRUE}
corr$scc <- abs(corr$scc)

pval <- sapply(
  split(
    corr[corr$rand=="Observed",], corr[corr$rand=="Observed", "mark"]
    ), 
  function(x) paste(
    "P=", format.pval(
      sum(median(x$scc)<corr_median[corr_median$mark==x$mark[1], "absscc"]) / nrow(corr_median[corr_median$mark==x$mark[1],]), eps=1/nrow(corr_median[corr_median$mark==x$mark[1],])
      )
    )
  )
pval <-  data.frame(mark=names(pval), pval=pval)

ggplot( corr, aes(scc) ) + 
  geom_line(stat="density", alpha=0.5, aes(color=status, linetype=rand)) + 
  facet_wrap(~mark,scales="free_y") + 
  scale_color_manual(values = ifelse( grepl("rand", levels(corr$status)), "gray", "red"), guide = "none") +
  geom_text(data = pval,aes(x = 0.5, y = 0.1, label = pval)) + 
  theme_bw() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) + 
  theme(aspect.ratio = 1) + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.title=element_blank()) + 
  theme(text = element_text(size=12)) +
  ylab("Density") + 
  xlab("Absolute Spearman Correlation")
```

