Updated script

- processed reviewer comments

power_analysis.R

+################################################################################
+###bar-1 eQTL study 
+################################################################################
+
+
+###Set your work directory
+setwd("E:/OneDrive - Wageningen University & Research/Shared_data_Mark/Projects/Biometris/R")
+workwd <- getwd()
+filename <- "4temp"
+
+###Load pre-made functions
+#uses eQTL pipeline functions https://git.wur.nl/mark_sterken/eQTL_pipeline
+#     transcriptomics functions https://git.wur.nl/mark_sterken/Transcriptomics.workbench
+#     genetic map functions https://git.wur.nl/mark_sterken/genetic.map
+git_dir <- "E:/OneDrive - Wageningen University & Research/Shared_data_Mark/Projects_R_Zone/Git"
+source(paste(git_dir,"/NEMA_functions/Loader.R",sep=""))
+
+###Set data locations
+
+
+
+################################################################################
+###Dependencies
+################################################################################
+
+install <- FALSE
+#.libPaths(.libPaths("C:/Program Files/R/R-3.4.2/library"))
+if(install){
+    install.packages("tidyverse")
+    install.packages("colorspace")
+    install.packages("RColorBrewer")
+    install.packages("BiocInstaller",repos="http://bioconductor.org/packages/3.3/bioc")
+    #source("http://www.bioconductor.org/biocLite.R") ; biocLite("limma") ; biocLite("statmod")
+    install.packages("gridExtra")
+    install.packages("VennDiagram")
+    install.packages("readxl")
+    install.packages("writexl")
+    install.packages("rmarkdown")
+}
+
+###load
+library("colorspace")
+library("RColorBrewer")
+library(limma)
+library(gridExtra)
+library("VennDiagram")
+library(readxl)
+library("rmarkdown")
+library(tidyverse)
+library(writexl)
+
+################################################################################
+###Plotting theme, colours
+################################################################################
+
+
+###Set plotting theme
+presentation <- theme(axis.text.x = element_text(size=10, face="bold", color="black"),
+                      axis.text.y = element_text(size=10, face="bold", color="black"),
+                      axis.title.x = element_text(size=20, face="bold", color="black"),
+                      axis.title.y = element_text(size=20, face="bold", color="black"),
+                      strip.text.x = element_text(size=20, face="bold", color="black"),
+                      strip.text.y = element_text(size=20, face="bold", color="black"),
+                      plot.title = element_text(size=24, face="bold"),
+                      panel.background = element_rect(fill = "white",color="black"),
+                      panel.grid.major = element_line(colour = "grey80"),
+                      panel.grid.minor = element_blank(),
+                      legend.position = "right")
+
+blank_theme <- theme(plot.background = element_blank(),
+                     panel.grid.major = element_blank(),
+                     panel.grid.minor = element_blank(),
+                     panel.border = element_blank(),
+                     panel.background = element_blank(),
+                     axis.title.x = element_blank(),
+                     axis.title.y = element_blank(),
+                     axis.text.x = element_blank(),
+                     axis.text.y = element_blank(),
+                     axis.ticks = element_blank())
+
+
+###Here you can set colours for plotting in theme using ggplot2
+#display.brewer.all()
+myColors <- c(brewer.pal(9,"Set1")[c(2,5,5,3)],"black","grey")
+
+
+names(myColors) <- c("CB4856","N2","N2bar1","trans","cis","RIL")
+
+colScale <- scale_colour_manual(name = "Treatment",values = myColors)
+fillScale <- scale_fill_manual(name = "Treatment",values = myColors)
+
+
+
+
+
+
+############################################################################
+###Power analysis (Supplementary table 4)
+############################################################################
+
+###load stuff
+data <- read_excel("Table S1.xlsx",sheet=1)
+
+strain.map <- data.matrix(wur.pop.map[,toupper(colnames(wur.pop.map)) %in% data$Line])
+
+###Power analysis
+biom.simulation <- QTL.map.1.sim(strain.map=strain.map,strain.marker=wur.pop.marker,threshold=3.4,n_per_marker=100)
+
+biom.simulation2 <- QTL.map.1.sim(strain.map=strain.map,strain.marker=wur.pop.marker,threshold=3.0,n_per_marker=100)
+
+
+write_xlsx(biom.simulation,path = "Supplementary_tablex1-power_analysis_pertemp.xlsx")
+write_xlsx(biom.simulation,path = "Supplementary_tablex2-power_analysis_plasticity.xlsx")
+
+###inspect file
+bar.simulation
+bar.simulation[7,]
+
+