MCmax25_asinT_groupStats.Rmd

---
title: "Nov3_MCmax25DMR_asinT_Stats"
author: "Shelly Trigg"
date: "11/03/2019"
output: rmarkdown::github_document
---

This script was run with Gannet mounted


load libraries
```{r}
library(gplots)
library(ggplot2)
library(dplyr)
library(broom)
library(RColorBrewer)
```

read in data
```{r}
DMRs <- read.table("../DMR_cov_in_0.75_SamplesPerCategory/DMR250bp_MCmax25_cov5x_rms_results_filtered.tsv", header = TRUE, sep = "\t")


#write out bedfile
write.table(DMRs[,1:4],"DMR250bp_MCmax25_cov5x_0.75xGroup.bed", sep = "\t",quote = FALSE, row.names = FALSE, col.names = FALSE )

```

Make a unique ID column 
```{r}
#for all ambient sample comparison
DMRs$ID <- paste(DMRs$chr,":",DMRs$start,"-",DMRs$end, sep = "")
DMRs$ID <- gsub("__.*__.*:",":",DMRs$ID)
```

reformat data for calculating group effect
```{r}
#reformat all to long format
DMRs_STACKED <- tidyr::gather(DMRs[,7:27], "group_name", "perc.meth",1:20)

#make group name column (e.g. 'CTRL_16C_26psu', '8C_32psu', etc.)
DMRs_STACKED$group_name <- gsub("methylation_level_","",DMRs_STACKED$group_name)

#make temperature column
DMRs_STACKED$temp <- gsub("C_.*","",DMRs_STACKED$group_name)
DMRs_STACKED$temp <- gsub("CTRL_","", DMRs_STACKED$temp)
DMRs_STACKED$temp <- as.factor(DMRs_STACKED$temp)

#make a salinity column
DMRs_STACKED$salinity <- gsub(".*C_","",DMRs_STACKED$group_name)
DMRs_STACKED$salinity <- gsub("psu_.*","",DMRs_STACKED$salinity)
DMRs_STACKED$salinity <- as.factor(DMRs_STACKED$salinity)

#make a treatment column (salinity and temp info)
DMRs_STACKED$treatment <- paste0(DMRs_STACKED$temp,"C_", DMRs_STACKED$salinity, "psu")

#make sample number column
DMRs_STACKED$sample_number <- gsub(".*_","",DMRs_STACKED$group_name)

#remove sample number from group name
DMRs_STACKED$group_name <- gsub("u_.*","u", DMRs_STACKED$group_name)

#create column for infestation status that's called "biotic_stress"
for (i in 1:nrow(DMRs_STACKED)){
  if(substr(DMRs_STACKED$group_name[i], 1,4) == "CTRL"){
    DMRs_STACKED$biotic_stress[i] <- "no sea lice"
  }
  else{DMRs_STACKED$biotic_stress[i] <- "sea lice"}
}
 

```

plot distribution of % methylation in all DMRs in all samples
```{r}
jpeg("DMR_percmeth_hist.jpg", width = 800, height = 800)
ggplot(DMRs_STACKED) + geom_histogram(aes(perc.meth, group = group_name, color = group_name,fill = group_name), bins = 10, position = "identity", alpha = 0.5) + theme_bw()
dev.off()


```


arc sin sqrt transform the data
```{r}
#arcsin sqrt transformation function
asinTransform <- function(p) { asin(sqrt(p))}

#arcsin transform data 
#day 10
DMRs_STACKED_asin <- DMRs_STACKED
DMRs_STACKED_asin$perc.meth <- asinTransform(DMRs_STACKED_asin$perc.meth)

```

plot distribution of TRANSFORMED % methylation in all DMRs in all samples
```{r}
jpeg("DMR_Tpercmeth_hist.jpg", width = 800, height = 800)
ggplot(DMRs_STACKED_asin) + geom_histogram(aes(perc.meth, group = group_name, color = group_name,fill = group_name), bins = 10, position = "identity", alpha = 0.5) + theme_bw()
dev.off()

```

## Run anova on TRANSFORMED data to assess group differences for each DMR
**hypothesis: There is no effect from biotic stress**
```{r}
#run an 1way ANOVA on 8C 26psu infested vs. 8C 26psu contrl
df_inf8C <- data.frame(DMRs_STACKED_asin[which(DMRs_STACKED_asin$temp == 8 & DMRs_STACKED_asin$salinity == 26),])

#remove ID if NA = true for both control samples
inf8C_rm <- rownames(DMR_m[which(is.na(DMR_m[,"methylation_level_CTRL_8C_26psu_1"]) & is.na(DMR_m[,"methylation_level_CTRL_8C_26psu_2"])),])

DMR_1way_aov_8Cinf <- df_inf8C[!(df_inf8C$ID%in% inf8C_rm),] %>% group_by(ID) %>% do(meth_aov_models = aov(perc.meth ~ biotic_stress, data =  . ))

#summarize ANOVA data
DMR_1way_aov_8Cinf_modelsumm <- glance(DMR_1way_aov_8Cinf, meth_aov_models)

#try 2way temp x biotic stress (only 26psu samples)
df_inf816C <- data.frame(DMRs_STACKED_asin[which(DMRs_STACKED_asin$salinity == 26),])

#remove ID if NA = true for both control samples
inf8C_rm <- rownames(DMR_m[which(is.na(DMR_m[,"methylation_level_CTRL_8C_26psu_1"]) & is.na(DMR_m[,"methylation_level_CTRL_8C_26psu_2"])),])
inf16C_rm <- rownames(DMR_m[which(is.na(DMR_m[,"methylation_level_CTRL_16C_26psu_1"]) & is.na(DMR_m[,"methylation_level_CTRL_16C_26psu_2"])),])



DMR_2way_aov_816Cinf <- df_inf816C[!(df_inf816C$ID%in% c(inf8C_rm,inf16C_rm)),] %>% group_by(ID) %>% do(meth_aov_models = aov(perc.meth ~ temp*biotic_stress, data =  . ))

DMR_2way_tuk_816Cinf <- df_inf816C[!(df_inf816C$ID%in% c(inf8C_rm,inf16C_rm)),] %>% group_by(ID) %>% do(meth_tuk_models = TukeyHSD(aov(perc.meth ~ biotic_stress*temp, data =  . )))

#summarize ANOVA data
DMR_2way_aov_816Cinf_modelsumm <- glance(DMR_2way_aov_816Cinf, meth_aov_models)

DMR_2way_tuk_816Cinf_modelsumm <- tidy(DMR_2way_tuk_816Cinf, meth_tuk_models)

#plot
DMR_inf816C_mean_m <- as.matrix(data.frame(cbind(MeanCTRL8C,MeanCTRL16C,Mean8C26psu,Mean16C26psu )))

DMR_inf816C_mean_m <- DMR_inf816C_mean_m[which(rownames(DMR_inf816C_mean_m) %in% pull(DMR_2way_tuk_816Cinf_modelsumm[which(DMR_2way_tuk_816Cinf_modelsumm$adj.p.value < 0.05),],ID)),]

colnames(DMR_inf816C_mean_m) <- gsub("methylation_level_","",colnames(DMR_inf816C_mean_m))


#simplify rownames and reorder columns 
rownames(tukTxS_0.05sal_DMR_mean_m) <- order(rownames(tukTxS_0.05sal_DMR_mean_m))

#define column (group) colors
ColSideColors <- cbind(group = c(rep("lightcyan1",1),rep("thistle1",1),rep("cyan",1),rep("magenta",1)))

#define palate for heatmap
col <- colorRampPalette(brewer.pal(5, "RdBu"))(256)


#heatmap with clustered groups and simplified numbers
heatmap.2(DMR_inf816C_mean_m,margins = c(10,20), cexRow = 1.5, cexCol = 1,ColSideColors = ColSideColors,Colv = NA,col = rev(col), na.color = "darkgray",sepwidth=c(0.025,0.025), sepcolor="darkgray",colsep=1:ncol(DMR_inf816C_mean_m),rowsep=1:nrow(DMR_inf816C_mean_m), density.info = "none", trace = "none", scale = "row")





#run an 1way ANOVA on 16C 26psu infested vs. 16C 26psu control



#find overlapping significant DMRs
```





**hypothesis: There is no effect from biotic stress**

```{r}
#run ANOVA testing for infestion effect pooling samples
DMR_1way_aov_inf <- DMRs_STACKED_asin %>% group_by(ID) %>%
do(meth_aov_models = aov(perc.meth~biotic_stress, data =  . ))
#summarize ANOVA data
DMR_1way_aov_inf_modelsumm <- glance(DMR_1way_aov_inf, meth_aov_models)
write.csv(DMR_1way_aov_inf_modelsumm, "DMR_MCmax25_1wayAOV_infest_modelsumm.csv", row.names = FALSE, quote = FALSE)





```

plot DMRs significant at 0.01
```{r}
jpeg("DMR_MCmax25DMR_Taov0.01InfestPercMeth.jpg", width = 10, height = 4, units = "in", res = 300)
p <- ggplot(data = DMRs_STACKED[which(DMRs_STACKED$ID %in% pull(DMR_1way_aov_inf_modelsumm[which(DMR_1way_aov_inf_modelsumm$p.value < 0.01),],ID)),],aes(x = biotic_stress,y = perc.meth, color = biotic_stress)) + facet_wrap(~ID, scale = "free") + theme_bw() + theme(axis.text.x = element_text(size = 7,angle = 45, hjust =1),axis.title=element_text(size=12,face="bold"))
p + geom_jitter(width = 0.3) +  ggtitle("DMRs that show an infestation effect significant at ANOVA p.value < 0.01")
dev.off()
```

plot DMRs significant at 0.05
```{r}
jpeg("DMR_MCmax25DMR_Taov0.05InfestPercMeth.jpg", width = 17, height = 12, units = "in", res = 300)
p <- ggplot(data = DMRs_STACKED[which(DMRs_STACKED$ID %in% pull(DMR_1way_aov_inf_modelsumm[which(DMR_1way_aov_inf_modelsumm$p.value < 0.05),],ID)),],aes(x = biotic_stress,y = perc.meth, color = biotic_stress)) + facet_wrap(~ID, scale = "free") + theme_bw() + theme(axis.text.x = element_text(size = 7,angle = 45, hjust =1),axis.title=element_text(size=12,face="bold"))
p + geom_jitter(width = 0.3) + ggtitle("DMRs that show an infestation effect significant at ANOVA p.value < 0.05")
dev.off()
```

plot DMRs significant at 0.1
```{r}
jpeg("DMR_MCmax25DMR_Taov0.1InfestPercMeth.jpg", width = 17, height = 12, units = "in", res = 300)
p <- ggplot(data = DMRs_STACKED[which(DMRs_STACKED$ID %in% pull(DMR_1way_aov_inf_modelsumm[which(DMR_1way_aov_inf_modelsumm$p.value < 0.1),],ID)),],aes(x = biotic_stress,y = perc.meth, color = biotic_stress)) + facet_wrap(~ID, scale = "free") + theme_bw() + theme(axis.text.x = element_text(size = 7,angle = 45, hjust =1),axis.title=element_text(size=12,face="bold"))
p + geom_jitter(width = 0.3) + ggtitle("DMRs that show an infestation effect significant at ANOVA p.value < 0.1")
dev.off()
```

plot heatmap of DMRs significant at 0.1
```{r}
#create matrix for all ambient samples

DMR_m <- as.matrix(DMRs[,7:26])
rownames(DMR_m) <- DMRs$ID


#calculate group means
#group means
Mean16C32psu <- rowMeans(DMR_m[,grep("16C_32psu", colnames(DMR_m))], na.rm = TRUE)

Mean16C26psu <- rowMeans(DMR_m[,grep("16C_26psu", colnames(DMR_m))], na.rm = TRUE)

Mean8C32psu <- rowMeans(DMR_m[,grep("8C_32psu", colnames(DMR_m))], na.rm = TRUE)

Mean8C26psu <- rowMeans(DMR_m[,grep("8C_26psu", colnames(DMR_m))], na.rm = TRUE)

MeanCTRL8C <- rowMeans(DMR_m[,grep("CTRL_8C_26psu", colnames(DMR_m))], na.rm = TRUE)

MeanCTRL16C <- rowMeans(DMR_m[,grep("CTRL_16C_26psu", colnames(DMR_m))], na.rm = TRUE)

DMR_mean_m <- as.matrix(data.frame(cbind(Mean8C26psu,Mean8C32psu,Mean16C26psu,Mean16C32psu,MeanCTRL8C, MeanCTRL16C)))


##ANOVA data
aov_0.1_DMR_m <- DMR_m[which(rownames(DMR_m) %in% pull(DMR_1way_aov_inf_modelsumm[which(DMR_1way_aov_inf_modelsumm$p.value < 0.1),],ID)),]

colnames(aov_0.1_DMR_m) <- gsub("methylation_level_","",colnames(aov_0.1_DMR_m))

ColSideColors <- cbind(group = c(rep("plum2",4),rep("plum4",4),rep("green1",4),rep("green3",4),rep("magenta",2), rep("cyan",2)))

jpeg("DMR_MCmax25DMR_Taov0.1_infest_heatmap.jpg", width = 800, height = 1000)
heatmap.2(aov_0.1_DMR_m,margins = c(10,20), cexRow = 1.2, cexCol = 1,ColSideColors = ColSideColors, Colv=NA, col = bluered, na.color = "black", density.info = "none", trace = "none", scale = "row")
dev.off()


#plot heatmap of mean methylation
uninf_mean <- rowMeans(aov_0.1_DMR_m[,grep("CTRL", colnames(aov_0.1_DMR_m))], na.rm = TRUE)
inf_mean <- rowMeans(aov_0.1_DMR_m[,-grep("CTRL", colnames(aov_0.1_DMR_m))], na.rm = TRUE)

Mean16C32psu <- rowMeans(aov_0.1_DMR_m[,grep("16C_32psu", colnames(aov_0.1_DMR_m))], na.rm = TRUE)

Mean16C26psu <- rowMeans(aov_0.1_DMR_m[,grep("16C_26psu", colnames(aov_0.1_DMR_m))], na.rm = TRUE)

Mean8C32psu <- rowMeans(aov_0.1_DMR_m[,grep("8C_32psu", colnames(aov_0.1_DMR_m))], na.rm = TRUE)

Mean8C26psu <- rowMeans(aov_0.1_DMR_m[,grep("8C_26psu", colnames(aov_0.1_DMR_m))], na.rm = TRUE)


aov_0.1_DMR_mean_m <- as.matrix(data.frame(cbind(Mean8C26psu,Mean8C32psu,Mean16C26psu,Mean16C32psu,uninf_mean)))

jpeg("DMR_MCmax25DMR_Taov0.1_infest_mean_heatmap.jpg", width = 8, height = 10, units = "in", res = 300)
ColSideColors <-  c("cyan", "darkcyan", "magenta", "darkmagenta","darkgray")
heatmap.2(aov_0.1_DMR_mean_m,margins = c(10,20), cexRow = 1.2, cexCol = 1, Colv=NA, ColSideColors = ColSideColors,col = bluered, na.color = "black", density.info = "none", trace = "none", scale = "row")
dev.off()
```

For the 24 DMRs significant at anova p.value of 0.1, do any show significant effect of temperature, salinity, or their interaction?



**Hypotehsis: there is no effect from temperature, salinity, or their interaction**
```{r}
DMR_2way_aov_tempxsal_nC <- DMRs_STACKED_asin[,-grep("CTRL",DMRs_STACKED_asin$group_name)] %>% group_by(ID) %>%
do(meth_aov_models = aov(perc.meth~temp*salinity, data =  . ))
#summarize ANOVA data
DMR_2way_aov_tempxsal_nC_modelsumm <- glance(DMR_2way_aov_tempxsal_nC, meth_aov_models)
```
plot significant regions
```{r}
jpeg("DMR_MCmax25DMR_Taov0.01TxSPercMeth.jpg", width = 17, height = 12, units = "in", res = 300)
p <- ggplot(data = DMRs_STACKED[which(DMRs_STACKED$ID %in% pull(DMR_2way_aov_tempxsal_nC_modelsumm[which(DMR_2way_aov_tempxsal_nC_modelsumm$p.value < 0.01),],ID) & substr(DMRs_STACKED$group_name, 1,4)!="CTRL"),],aes(x = group_name,y = perc.meth, color = group_name)) + facet_wrap(~ID, scale = "free") + theme_bw() + theme(axis.text.x = element_text(size = 7,angle = 45, hjust =1),axis.title=element_text(size=12,face="bold"))
p + geom_jitter(width = 0.3) +  ggtitle("DMRs that show an overall TxS model effect significant at ANOVA p.value < 0.01")
dev.off()

jpeg("DMR_MCmax25DMR_Taov0.05TxSPercMeth.jpg", width = 17, height = 12, units = "in", res = 300)
p <- ggplot(data = DMRs_STACKED[which(DMRs_STACKED$ID %in% pull(DMR_2way_aov_tempxsal_nC_modelsumm[which(DMR_2way_aov_tempxsal_nC_modelsumm$p.value < 0.05),],ID) & substr(DMRs_STACKED$group_name, 1,4)!="CTRL"),],aes(x = group_name,y = perc.meth, color = group_name)) 

p + geom_boxplot()+ geom_jitter(width = 0.3) + ggtitle("DMRs that show an overall TxS model effect significant at ANOVA p.value < 0.05") + facet_wrap(~ID, scale = "free") +theme_bw() + theme(axis.text.x = element_text(size = 7,angle = 45, hjust =1),axis.title=element_text(size=12,face="bold"))
dev.off()
```

plot heatmap of TxS 0.05 sig DMRs
```{r}
aov_0.05TxS_DMR_m <- DMR_m[which(rownames(DMR_m) %in% pull(DMR_2way_aov_tempxsal_nC_modelsumm[which(DMR_2way_aov_tempxsal_nC_modelsumm$p.value < 0.05),],ID)),]

colnames(aov_0.05TxS_DMR_m) <- gsub("methylation_level_","",colnames(aov_0.05TxS_DMR_m))

ColSideColors <- cbind(group = c(rep("magenta",4),rep("darkmagenta",4),rep("cyan",4),rep("darkcyan",4),rep("thistle1",2), rep("lightcyan1",2)))

jpeg("DMR_MCmax25DMR_Taov0.05_TxS_heatmap.jpg", width = 800, height = 1000)
heatmap.2(aov_0.05TxS_DMR_m,margins = c(10,20), cexRow = 1.2, cexCol = 1,ColSideColors = ColSideColors, Colv=NA, col = bluered, na.color = "black", density.info = "none", trace = "none", scale = "row")
dev.off()

#group means
Mean16C32psu <- rowMeans(aov_0.05TxS_DMR_m[,grep("16C_32psu", colnames(aov_0.05TxS_DMR_m))], na.rm = TRUE)

Mean16C26psu <- rowMeans(aov_0.05TxS_DMR_m[,grep("16C_26psu", colnames(aov_0.05TxS_DMR_m))], na.rm = TRUE)

Mean8C32psu <- rowMeans(aov_0.05TxS_DMR_m[,grep("8C_32psu", colnames(aov_0.05TxS_DMR_m))], na.rm = TRUE)

Mean8C26psu <- rowMeans(aov_0.05TxS_DMR_m[,grep("8C_26psu", colnames(aov_0.05TxS_DMR_m))], na.rm = TRUE)

MeanCTRL8C <- rowMeans(aov_0.05TxS_DMR_m[,grep("CTRL_8C_26psu", colnames(aov_0.05TxS_DMR_m))], na.rm = TRUE)

MeanCTRL16C <- rowMeans(aov_0.05TxS_DMR_m[,grep("CTRL_16C_26psu", colnames(aov_0.05TxS_DMR_m))], na.rm = TRUE)


aov_0.05TxS_DMR_mean_m <- as.matrix(data.frame(cbind(Mean8C26psu,Mean8C32psu,Mean16C26psu,Mean16C32psu,MeanCTRL8C, MeanCTRL16C)))

ColSideColors <- cbind(group = c(rep("cyan",1),rep("darkcyan",1),rep("magenta",1),rep("darkmagenta",1), rep("lightcyan1",1),rep("thistle1",1)))
col <- colorRampPalette(brewer.pal(10, "RdYlBu"))(256)
jpeg("DMR_MCmax25DMR_Taov0.05_TxSmean_heatmap.jpg", width = 800, height = 1000)
heatmap.2(aov_0.05TxS_DMR_mean_m,margins = c(10,20), cexRow = 1.2, cexCol = 1,ColSideColors = ColSideColors, Colv=NA, col = rev(col), na.color = "black", density.info = "none", trace = "none", scale = "row")
dev.off()



###try a looser cutoff of 0.1
aov_0.1TxS_DMR_m <- DMR_m[which(rownames(DMR_m) %in% pull(DMR_2way_aov_tempxsal_nC_modelsumm[which(DMR_2way_aov_tempxsal_nC_modelsumm$p.value < 0.1),],ID)),]

colnames(aov_0.1TxS_DMR_m) <- gsub("methylation_level_","",colnames(aov_0.1TxS_DMR_m))


#group means
Mean16C32psu <- rowMeans(aov_0.1TxS_DMR_m[,grep("16C_32psu", colnames(aov_0.1TxS_DMR_m))], na.rm = TRUE)

Mean16C26psu <- rowMeans(aov_0.1TxS_DMR_m[,grep("16C_26psu", colnames(aov_0.1TxS_DMR_m))], na.rm = TRUE)

Mean8C32psu <- rowMeans(aov_0.1TxS_DMR_m[,grep("8C_32psu", colnames(aov_0.1TxS_DMR_m))], na.rm = TRUE)

Mean8C26psu <- rowMeans(aov_0.1TxS_DMR_m[,grep("8C_26psu", colnames(aov_0.1TxS_DMR_m))], na.rm = TRUE)

MeanCTRL8C <- rowMeans(aov_0.1TxS_DMR_m[,grep("CTRL_8C_26psu", colnames(aov_0.1TxS_DMR_m))], na.rm = TRUE)

MeanCTRL16C <- rowMeans(aov_0.1TxS_DMR_m[,grep("CTRL_16C_26psu", colnames(aov_0.1TxS_DMR_m))], na.rm = TRUE)


aov_0.1TxS_DMR_mean_m <- as.matrix(data.frame(cbind(MeanCTRL8C, Mean8C26psu,Mean8C32psu,MeanCTRL16C,Mean16C26psu,Mean16C32psu)))
ColSideColors <- cbind(group = c(rep("lightcyan1",1),rep("cyan",1),rep("darkcyan",1),rep("thistle1",1),rep("magenta",1),rep("darkmagenta",1)))
col <- colorRampPalette(brewer.pal(5, "RdBu"))(256)
jpeg("DMR_MCmax25DMR_Taov0.1_TxSmean_heatmap.jpg", width = 800, height = 1000)
heatmap.2(aov_0.1TxS_DMR_mean_m,margins = c(10,20), cexRow = 1.2, cexCol = 1,ColSideColors = ColSideColors, Colv=NA, col = rev(col), na.color = "black",sepwidth=c(0.025,0.025), sepcolor="darkgray",colsep=1:ncol(aov_0.1TxS_DMR_mean_m),rowsep=1:nrow(aov_0.1TxS_DMR_mean_m), density.info = "none", trace = "none", scale = "row")
dev.off()
```


run tukey
```{r}
#subset data for DMRs with TxS 2way ANOVA p.value < 0.1 and run tukey HSD for each DMR
DMR_tuk_tempxsal_nC <- DMRs_STACKED_asin[,-grep("CTRL",DMRs_STACKED_asin$group_name)] %>% group_by(ID) %>%
do(meth_tuk_models = TukeyHSD(aov(perc.meth~temp*salinity, data =  . )))
#summarize ANOVA data
DMR_2way_aov_tempxsal_nC_modelsumm <- glance(DMR_2way_aov_tempxsal_nC, meth_tuk_models)
  
#summarize TukeyHSD results
DMR_tuk_tempxsal_nC_summ <- tidy(DMR_tuk_tempxsal_nC, meth_tuk_models)

#write out DMRs sig @ p <= 0.05
write.csv(DMR_tuk_tempxsal_nC_summ[which(DMR_tuk_tempxsal_nC_summ$adj.p.value <= 0.05),],"DMR_tuk_tempxsal_nC_summ.csv", row.names = FALSE, quote = FALSE )
```


plot chromosomes with sig Tuk DMRs
```{r}

chr_data <- read.table("../../data/chr_name_conversion.txt", sep = "\t", header = FALSE)
colnames(chr_data) <- c("chr.num","ssa.num","chr")

sig_chr <- rownames(DMR_mean_m[which(rownames(DMR_mean_m) %in% pull(DMR_tuk_tempxsal_nC_summ[which(DMR_tuk_tempxsal_nC_summ$adj.p.value < 0.05),],ID)),])
sig_chr <- gsub(":.*","",sig_chr)
sig_chr <- data.frame(table(sig_chr))
colnames(sig_chr)[1] <- "chr"
sig_chr <- merge(chr_data, sig_chr, by = "chr", all = TRUE)
sig_chr[is.na(sig_chr$Freq)] <- 0
sig_chr$perc <- sig_chr$Freq/52
sig_chr$type <- "sig"
ggplot(sig_chr, aes(x= chr)) + geom_bar(aes(y = perc),stat = "identity")



all_chr_cov <- read.table("../../data/chr_5x0.75grp_cov.txt", sep = " ", header = FALSE)
colnames(all_chr_cov) <- c("Freq", "chr")
all_chr_cov <- merge(chr_data,all_chr_cov, by = "chr", all = TRUE)
all_chr_cov$type <- "bkgd"
all_chr_cov$perc <- all_chr_cov$Freq/sum(all_chr_cov$Freq)

chr_sig_bkgd <- rbind(sig_chr,all_chr_cov)

#plot chromosome coverage by DMRs vs. background
ggplot(chr_sig_bkgd, aes(x= ssa.num, group = type, fill = type)) + geom_bar(aes(y = perc),stat = "identity", position = "dodge") + theme(axis.text.x = element_text(size = 10, hjust= 1,angle = 90)) 


  
```


try plotting heatmap of tukey

```{r}

tukTxS_0.05sal_DMR_mean_m <- DMR_mean_m[which(rownames(DMR_mean_m) %in% pull(DMR_tuk_tempxsal_nC_summ[which(DMR_tuk_tempxsal_nC_summ$adj.p.value < 0.05 & DMR_tuk_tempxsal_nC_summ$term == "salinity"),],ID)),]

colnames(tukTxS_0.05sal_DMR_mean_m) <- gsub("methylation_level_","",colnames(tukTxS_0.05sal_DMR_mean_m))

#reorder columns so samples with same salinities are next to each other
tukTxS_0.05sal_DMR_mean_m <- tukTxS_0.05sal_DMR_mean_m[,c(2,4,1,3,5,6)]

#simplify rownames and reorder columns 
rownames(tukTxS_0.05sal_DMR_mean_m) <- order(rownames(tukTxS_0.05sal_DMR_mean_m))

#define column (group) colors
ColSideColors <- cbind(group = c(rep("darkcyan",1),rep("darkmagenta",1),rep("cyan",1),rep("magenta",1),rep("lightcyan1",1),rep("thistle1",1)))

#define palate for heatmap
col <- colorRampPalette(brewer.pal(5, "RdBu"))(256)


#heatmap with clustered groups and simplified numbers
jpeg("DMR_MCmax25DMR_TtukTxS0.05sal_mean_heatmap.jpg", width = 800, height = 1000)
heatmap.2(tukTxS_0.05sal_DMR_mean_m,margins = c(10,20), cexRow = 1.5, cexCol = 1,ColSideColors = ColSideColors, Colv= NA,col = rev(col), na.color = "darkgray",sepwidth=c(0.025,0.025), sepcolor="darkgray",colsep=1:ncol(tukTxS_0.05sal_DMR_mean_m),rowsep=1:nrow(tukTxS_0.05sal_DMR_mean_m), density.info = "none", trace = "none", scale = "row")
dev.off()


#heatmap.2(tukTxS_0.05sal_DMR_mean_m,margins = c(10,20), cexRow = 1.2, cexCol = 1,ColSideColors = ColSideColors, col = rev(col), na.color = "black",sepwidth=c(0.025,0.025), sepcolor="darkgray",colsep=1:ncol(tukTxS_0.05sal_DMR_mean_m),rowsep=1:nrow(tukTxS_0.05sal_DMR_mean_m), density.info = "none", trace = "none", scale = "row")

#temp < 0.05
tukTxS_0.05temp_DMR_mean_m <- DMR_mean_m[which(rownames(DMR_mean_m) %in% pull(DMR_tuk_tempxsal_nC_summ[which(DMR_tuk_tempxsal_nC_summ$adj.p.value < 0.05 & DMR_tuk_tempxsal_nC_summ$term == "temp"),],ID)),]

colnames(tukTxS_0.05temp_DMR_mean_m) <- gsub("methylation_level_","",colnames(tukTxS_0.05temp_DMR_mean_m))

#simplify row names
rownames(tukTxS_0.05temp_DMR_mean_m) <- order(rownames(tukTxS_0.05temp_DMR_mean_m))

#define column (group) colors
ColSideColors <- cbind(group = c(rep("cyan",1),rep("darkcyan",1),rep("magenta",1),rep("darkmagenta",1),rep("lightcyan1",1),rep("thistle1",1)))

col <- colorRampPalette(brewer.pal(5, "RdBu"))(256)
jpeg("DMR_MCmax25DMR_TtukTxS0.05temp_mean_heatmap.jpg", width = 800, height = 1000)
heatmap.2(tukTxS_0.05temp_DMR_mean_m,margins = c(10,20), cexRow = 2, cexCol = 1,ColSideColors = ColSideColors, Colv=NA, col = rev(col), na.color = "black",sepwidth=c(0.025,0.025), sepcolor="darkgray",colsep=1:ncol(tukTxS_0.05temp_DMR_mean_m),rowsep=1:nrow(tukTxS_0.05temp_DMR_mean_m), density.info = "none", trace = "none", scale = "row")
dev.off()

```

```{r}
DMR_2wayAOVtuk0.05_bed <- DMRs[which(DMRs$ID %in% pull(DMR_tuk_tempxsal_nC_summ[which(DMR_tuk_tempxsal_nC_summ$adj.p.value < 0.05),],ID)),1:4]

write.table(DMR_2wayAOVtuk0.05_bed, "DMR_2wayAOVTxStuk0.05.bed", sep = "\t", quote = FALSE, row.names = FALSE, col.names = FALSE)


```




**Hypotehsis: there is no effect from temperature, salinity, or their interaction**

```{r}
#subset data for DMRs that show infestation effect significant at ANOVA p.val < 0.1 and excluding control samples
DMR_STACKED_asin_inf_sig_0.1 <- DMRs_STACKED_asin[which(DMRs_STACKED_asin$ID %in% pull(DMR_1way_aov_inf_modelsumm[which(DMR_1way_aov_inf_modelsumm$p.value < 0.1),],ID) & substr(DMRs_STACKED_asin$group_name,1,4)!="CTRL"),]

#run 2way temp x salinity ANOVA on each DMR
DMR_2way_aov_TxS <- DMR_STACKED_asin_inf_sig_0.1 %>% group_by(ID) %>%
do(meth_aov_models = aov(perc.meth~temp*salinity, data =  . ))

#create ANOVA summary tables
DMR_2way_aov_TxS_modelsumm <- glance(DMR_2way_aov_TxS, meth_aov_models)
DMR_2way_aov_TxS_coeff <- tidy(DMR_2way_aov_TxS, meth_aov_models)

#subset data for DMRs with TxS 2way ANOVA p.value < 0.1 and run tukey HSD for each DMR
DMR_tuk_TxS <- DMR_STACKED_asin_inf_sig_0.1[which(DMR_STACKED_asin_inf_sig_0.1$ID %in% pull(DMR_2way_aov_TxS_modelsumm[which(DMR_2way_aov_TxS_modelsumm$p.value < 0.1),], ID)),] %>% group_by(ID) %>% do(meth_tuk_models = TukeyHSD(aov(perc.meth~temp*salinity, data =  . )))

#summarize TukeyHSD results
DMR_tuk_TxS_summ <- tidy(DMR_tuk_TxS, meth_tuk_models)
```

Write out 2way ANOVA and Tukey data summaries 
```{r}

write.csv(DMR_2way_aov_TxS_modelsumm, "DMR_MCmax25_2wayAOV_TxS_modelsumm.csv", row.names = FALSE, quote = FALSE)

write.csv(DMR_tuk_TxS_summ, "DMR_MCmax25_DMR_TukHSD_TxS_modelsumm.csv", row.names = FALSE, quote = FALSE)

```


Create bed file for 24 DMRs with infestation effect significant at ANOVA p.value < 0.1
```{r}
DMR_1way_aov_inf_modelsumm_4bed <- merge(DMRs[,c("ID","chr", "start","end")], DMR_1way_aov_inf_modelsumm[which(DMR_1way_aov_inf_modelsumm$p.value < 0.1),], by = "ID")

DMR_1way_aov_inf_modelsumm_4bed <- DMR_1way_aov_inf_modelsumm_4bed[,2:4]

write.table(DMR_1way_aov_inf_modelsumm_4bed, "DMR_2way_aov_infest_modelsumm_4bed.txt", sep = "\t", quote = FALSE, row.names = FALSE)

write.table(DMR_1way_aov_inf_modelsumm_4bed, "DMR_1way_aov0.1_infest_modelsumm_4bed.txt", sep = "\t", quote = FALSE, row.names = FALSE)
```

plot DMRs with salinity effect significant at adj.p.value < 0.1 TukHSD
```{r}
jpeg("DMR_MCmax25DMR_Taov0.1SalPercMeth.jpg", width = 10, height = 4, units = "in", res = 300)
p <- ggplot(data = DMRs_STACKED[which(DMRs_STACKED$ID %in% pull(DMR_tuk_TxS_summ[which(DMR_tuk_TxS_summ$adj.p.value < 0.1 & DMR_tuk_TxS_summ$term == "salinity"),],ID) & substr(DMRs_STACKED$group_name,1,4)!="CTRL"),],aes(x = group_name,y = perc.meth, color = group_name)) + facet_wrap(~ID, scale = "free") + theme_bw() + theme(axis.text.x = element_text(size = 7,angle = 45, hjust =1),axis.title=element_text(size=12,face="bold"))
p + geom_jitter(width = 0.3) + ggtitle("DMRs that show a salinity effect significant at TukeyHSD p.value < 0.1")
dev.off()
```

plot DMRs with temperature effect significant at adj.p.value < 0.1 TukHSD
```{r}
jpeg("DMR_MCmax25DMR_Taov0.1TempPercMeth.jpg", width = 10, height = 4, units = "in", res = 300)
p <- ggplot(data = DMRs_STACKED[which(DMRs_STACKED$ID %in% pull(DMR_tuk_TxS_summ[which(DMR_tuk_TxS_summ$adj.p.value < 0.1 & DMR_tuk_TxS_summ$term == "temp"),],ID) & substr(DMRs_STACKED$group_name,1,4)!="CTRL"),],aes(x = group_name,y = perc.meth, color = group_name)) + facet_wrap(~ID, scale = "free") + theme_bw() + theme(axis.text.x = element_text(size = 7,angle = 45, hjust =1),axis.title=element_text(size=12,face="bold"))
p + geom_jitter(width = 0.3) + ggtitle("DMRs that show a temperature effect significant at TukeyHSD p.value < 0.1")
dev.off()
```

plot DMRs with Temp x Salinity interaction effect significant at 0.1 TukHSD
```{r}
jpeg("DMR_MCmax25DMR_Taov0.1TxSPercMeth.jpg", width = 10, height = 8, units = "in", res = 300)
p <- ggplot(data = DMRs_STACKED[which(DMRs_STACKED$ID %in% pull(DMR_tuk_TxS_summ[which(DMR_tuk_TxS_summ$adj.p.value < 0.1 & DMR_tuk_TxS_summ$term == "temp:salinity"),],ID) & substr(DMRs_STACKED$group_name,1,4)!="CTRL"),],aes(x = group_name,y = perc.meth, color = group_name)) + facet_wrap(~ID, scale = "free") + theme_bw() + theme(axis.text.x = element_text(size = 7,angle = 45, hjust =1),axis.title=element_text(size=12,face="bold"))
p + geom_jitter(width = 0.3) + ggtitle("DMRs that show a temp x salinity effect significant at TukeyHSD p.value < 0.1")
dev.off()
```

**hypothesis: There is no group difference between the following samples:**

- 16C\_26psu_\infested 
- 16C\_32psu\_infested
- 8C\_26psu\_infested 
- 8C\_32psu\_infested 
- 16C\_26psu\_NOTinfested
- 8C\_26psu\_NOTinfested 

```{r}
#run 1 way ANOVA on group name (e.g. 'CTRL_16C_26psu', '8C_32psu', etc.)
DMR_aov = DMRs_STACKED_asin %>% group_by(ID) %>%
do(meth_aov_models = aov(perc.meth ~ group_name, data =  . ))
#summarize ANOVA data
DMR_aov_modelsumm <- glance(DMR_aov, meth_aov_models)
```

**hypothesis: There is no effect from salinity, temperature, biotic stress, or their interactions**
```{r}
#run 3 way ANOVA
DMR_3way_aov <- 
DMRs_STACKED_asin %>% group_by(ID) %>%
do(meth_aov_models = aov(perc.meth~temp*salinity*biotic_stress, data =  . ))
#summarize 3way ANOVA data
DMR_3way_aov_modelsumm <- tidy(DMR_3way_aov, meth_aov_models)

#run tukey test on 3way ANOVA
DMR_3way_tuk <- 
DMRs_STACKED_asin %>% group_by(ID) %>%
do(tuk_aov_models = TukeyHSD(aov(perc.meth~temp*salinity*biotic_stress, data =  . )))
#summarize tukey results
DMR_3way_tuk_modelsumm <- tidy(DMR_3way_tuk, tuk_aov_models)
#subset tukey results with adj.p.value < 0.05
DMR_3way_tuk_modelsumm_0.05 <- DMR_3way_tuk_modelsumm[which(DMR_3way_tuk_modelsumm$adj.p.value < 0.05),]

#for DMRs with p.val < 0.05 for overall model, get term p.val; I think this step is redundant with the Tukey test
DMR_3way_aov_modelsumm_0.05 <- DMR_3way_aov_modelsumm[which(DMR_3way_aov_modelsumm$ID %in% pull(DMR_aov_modelsumm[which(DMR_aov_modelsumm$p.value < 0.05),],ID)),]
```

Subset out DMRs that showed an infestation effect significant at p < 0.05 in 3way ANOVA
```{r}
DMR_3way_aov_inf_modelsumm_0.05 <- DMR_3way_aov_modelsumm_0.05[which(DMR_3way_aov_modelsumm_0.05$ID %in% pull(DMR_3way_aov_modelsumm_0.05[which(DMR_3way_aov_modelsumm_0.05$term == "biotic_stress" & DMR_3way_aov_modelsumm_0.05$p.value < 0.05 ),],ID)),]
```

```{bash}
mkdir 3wayAOV
```

create abundance plots of DMRs with significant term effects for each term
```{r}
jpeg("3wayAOV/DMR_MCmax25DMR_Taov0.05_salinity_boxplots.jpg", width = 17, height = 7, units = "in", res = 300)
p <- ggplot(data = DMRs_STACKED[which(DMRs_STACKED$ID %in% pull(DMR_3way_aov_inf_modelsumm_0.05[which(DMR_3way_aov_inf_modelsumm_0.05$p.value < 0.05 & DMR_3way_aov_inf_modelsumm_0.05$term == "salinity"),],ID)),],aes(x = group_name,y = perc.meth)) + geom_violin(aes(fill = group_name), trim = FALSE)+ facet_wrap(~ID, scale = "free") + theme_bw() + theme(axis.text.x = element_text(size = 7,angle = 45, hjust =1),axis.title=element_text(size=12,face="bold"))
p + geom_jitter(width = 0.2)
dev.off()

jpeg("3wayAOV/DMR_MCmax25DMR_Taov0.05_tempSal_boxplots.jpg", width = 17, height = 7, units = "in", res = 300)
p <- ggplot(data = DMRs_STACKED[which(DMRs_STACKED$ID %in% pull(DMR_3way_aov_inf_modelsumm_0.05[which(DMR_3way_aov_inf_modelsumm_0.05$p.value < 0.05 & DMR_3way_aov_inf_modelsumm_0.05$term == "temp:salinity"),],ID)),],aes(x = group_name,y = perc.meth)) + geom_violin(aes(fill = group_name), trim = FALSE)+ facet_wrap(~ID, scale = "free") + theme_bw() + theme(axis.text.x = element_text(size = 7,angle = 45, hjust =1),axis.title=element_text(size=12,face="bold"))
p + geom_jitter(width = 0.2)
dev.off()

jpeg("3wayAOV/DMR_MCmax25DMR_Taov0.05_temp_boxplots.jpg", width = 17, height = 7, units = "in", res = 300)
p <- ggplot(data = DMRs_STACKED[which(DMRs_STACKED$ID %in% pull(DMR_3way_aov_inf_modelsumm_0.05[which(DMR_3way_aov_inf_modelsumm_0.05$p.value < 0.05 & DMR_3way_aov_inf_modelsumm_0.05$term == "temp"),],ID)),],aes(x = group_name,y = perc.meth)) + geom_violin(aes(fill = group_name), trim = FALSE)+ facet_wrap(~ID, scale = "free") + theme_bw() + theme(axis.text.x = element_text(size = 7,angle = 45, hjust =1),axis.title=element_text(size=12,face="bold"))
p + geom_jitter(width = 0.2)
dev.off()

jpeg("3wayAOV/DMR_MCmax25DMR_Taov0.05_temp_biostress_boxplots.jpg", width = 17, height = 15, units = "in", res = 300)
p <- ggplot(data = DMRs_STACKED[which(DMRs_STACKED$ID %in% pull(DMR_3way_aov_inf_modelsumm_0.05[which(DMR_3way_aov_inf_modelsumm_0.05$p.value < 0.05 & DMR_3way_aov_inf_modelsumm_0.05$term == "temp:biotic_stress"),],ID)),],aes(x = group_name,y = perc.meth)) + geom_violin(aes(fill = group_name), trim = FALSE)+ facet_wrap(~ID, scale = "free") + theme_bw() + theme(axis.text.x = element_text(size = 7,angle = 45, hjust =1),axis.title=element_text(size=12,face="bold"))
p + geom_jitter(width = 0.2)
dev.off()

#look at DMRs with tukey < 0.05

#first create data frame of unique IDs and their frequency in the results table with Tukey HSD adj.p.value < 0.05
sig_DMR_freq <- data.frame(table(DMR_3way_tuk_modelsumm_0.05$ID))

#subset data with the dataframe above
tuk_0.05_DMR_m <- DMR_m[which(rownames(DMR_m) %in% unlist(unique(DMR_3way_tuk_modelsumm_0.05[which(DMR_3way_tuk_modelsumm_0.05$ID %in% sig_DMR_freq[which(sig_DMR_freq$Freq > 1),"Var1"]),"ID"]))),]

#plot
jpeg("3wayAOV/DMR_MCmax25DMR_tuk0.05_heatmap.jpg", width = 800, height = 1000)
heatmap.2(tuk_0.05_DMR_m,margins = c(5,20), cexRow = 1.2, cexCol = 1,ColSideColors = ColSideColors, Colv=NA, col = bluered, na.color = "black", density.info = "none", trace = "none", scale = "row")
dev.off()

```


Plot DMRs with significant group effect at 1way ANOVA p < 0.05 

first create a directory for saving the images
```{bash}
mkdir 1wayAOVgroup_name
```

```{r}

##ANOVA data
aov_0.05_DMRs_m <- DMR_m[which(rownames(DMR_m) %in% pull(DMR_aov_modelsumm[which(DMR_aov_modelsumm$p.value < 0.05),],ID)),]

ColSideColors <- cbind(pH = c(rep("plum2",4),rep("plum4",4),rep("green1",4),rep("green3",4),rep("magenta",2), rep("cyan",2)))

jpeg("1wayAOVgroup_name/DMR_MCmax25DMR_Taov0.05_heatmap.jpg", width = 800, height = 1000)
heatmap.2(aov_0.05_DMRs_m,margins = c(5,20), cexRow = 1.2, cexCol = 1,ColSideColors = ColSideColors, Colv=NA, col = bluered, na.color = "black", density.info = "none", trace = "none", scale = "row")
dev.off()

```