---
title: "SimulationsCreate"
author: "Ankur Chakravarthy"
date: "19 November 2017"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## R Markdown

Create datasets

```{r}

library(dplyr)
library(ggplot2)
library(reshape2)

Percentages <- c(0.001/100, 0.01/100, 0.1/100, 1/100, 10/100)
DMRs <- c(1,10,100,1000)
Coverage <- c(10,100,1000,10000)

TabData <- expand.grid(Percentages, DMRs, Coverage)
colnames(TabData) <- c("Abundance","DMRs","Coverage")


load("SimulationOutput_10Iters.RData")  # Created from NewSim_10Iters.R

#The structure of these objects is
#Each row is an iteration
#Each column is a DMR 
#Each cell represents number of times we detected that DMR at the given coverage
#So what we need to do is every time a DMR was detected, we set it to 1, sum the rows to see how many DMRs were detected, set that iteration to 1 if the n(DMR) > threshold, and finally sum and divide by 10 to give p(DMR)


S2 <- SamplingTests.list
for(i in 1:length(S2)) {
  
  S2[[i]][S2[[i]] > 1] <- 1
  
}

S2 <- lapply(S2, function(x) rowSums(x))
S2 <- lapply(S2, function(x) ifelse(x > 0, 1, 0))
S2 <- lapply(S2, function(x) sum(x)/10)

TabData$AtleastOne <- unlist(S2)

#Then 5 DMRs 


S2 <- SamplingTests.list

for(i in 1:length(S2)) {
  
  S2[[i]][S2[[i]] > 1] <- 1
  
}

S2 <- lapply(S2, function(x) rowSums(x))
S2 <- lapply(S2, function(x) ifelse(x > 5, 1, 0))
S2 <- lapply(S2, function(x) sum(x)/10)

TabData$AtleastFive <- unlist(S2)


#Then 10 DMRs 

S2 <- SamplingTests.list

for(i in 1:length(S2)) {
  
  S2[[i]][S2[[i]] > 1] <- 1
  
}

S2 <- lapply(S2, function(x) rowSums(x))
S2 <- lapply(S2, function(x) ifelse(x > 10, 1, 0))
S2 <- lapply(S2, function(x) sum(x)/10)

TabData$AtleastTen <- unlist(S2)


```

## Including Plots

Plot out

```{r}

library(ggplot2)

pdf(width = 10, height = 8, file = "SimulationPDFs.pdf")
qplot(data = TabData, y = AtleastOne, facets= DMRs~ (100 * Abundance), x = Coverage, geom = "blank")+
  geom_line(colour = I("dodgerblue3"), size = I(1))+
  geom_point(size = I(2), colour = I("orange"))+
  scale_x_log10(breaks = c(10,100,1000,10000))+
  ylab("Probability of detecting 1+ Epimutations")+
  xlab("Coverage")+
  theme_bw()
  

qplot(data = TabData, y = AtleastFive, facets= DMRs~ (100 *Abundance), x = Coverage, geom = "blank")+
  geom_line(colour = I("dodgerblue3"), size = I(1))+
  geom_point(size = I(2), colour = I("orange"))+
  scale_x_log10(breaks = c(10,100,1000,10000))+
  ylab("Probability of detecting 5+ Epimutations")+
  xlab("Coverage")+
  theme_bw()



qplot(data = TabData, y = AtleastTen, facets= DMRs~ (100 *Abundance), x = Coverage, geom = "blank")+
  geom_line(colour = I("dodgerblue3"), size = I(1))+
  geom_point(size = I(2), colour = I("orange"))+
  scale_x_log10(breaks = c(10,100,1000,10000))+
  ylab("Probability of detecting 10+ Epimutations")+
  xlab("Coverage")+
  theme_bw()

dev.off()

```
