---
title: "Identification of varible FIRE region in 20 individuals"
output:
  pdf_document: default
---

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

```{r}
library(limma)
library(ggplot2)
library(reshape2)
library(preprocessCore)

set.seed(0)
```

Here we identified var FIRE bins at 40Kb resolution removing all bins within 200Kb of large SVs (>10K) and contain any SVs in itself.

#Load data and cleaning
```{r load_data,cache=TRUE}
dat <- read.table("raw_dat/FIRE_LCL20_Reps_40Kb_no_sv.mat", head=T)
row.names(dat) <- dat$bin
dat <- dat[, -1]
dat_z <- apply(dat, 2, function(x) (x-mean(x))/sd(x)) #z score transformation
dat_z <- dat_z[apply(dat_z>1.645, 1, sum) > 1,] # one-side p-value < 0.05
dat <- dat[rownames(dat) %in% rownames(dat_z),]
indv <- gsub("_Rep.", "", colnames(dat)) 
```

#Identify varibale bins using limma

```{r, fig.height=8, fig.width=8, eval=FALSE}
#Check if FIRE for each bin is normal distributed
dat_m <- melt(as.matrix(dat))
colnames(dat_m) <- c("bin", "replicate", "DI")
dat_m$indv <- gsub("_Rep.","",dat_m$replicate)
ggplot(dat_m[dat_m$bin %in% sample(dat_m$bin, 9),], aes(DI)) + geom_line(stat="density") + facet_wrap(~bin, scales = "free")
```

```{r limma_test,cache=TRUE}
design <- model.matrix(~indv)
fit <- lmFit(dat,design)
fit <- eBayes(fit)
write.table(topTable(fit, number = Inf, p.value = 1), "indv_test_limma_FIRE_LCL20_40Kb.tsv", quote=F, sep="\t", row.names=T, col.names=T)
```

There are `r nrow(topTable(fit, number = Inf, p.value = 0.1))` bins out of total `r nrow(dat)` bins are significant under adjusted p value 0.1.


