Commit 90715bd3 authored by Sterken, Mark's avatar Sterken, Mark
Browse files

Analysis Karengera et al., indirect acting xenobiotics.

parents
.Rproj.user
.Rhistory
.RData
.Ruserdata
################################################################################
###Wetsus toxicology experiment: indirect acting xenobiotics
################################################################################
###Set your work directory
setwd("C:/Users/sterk009/OneDrive - Wageningen University & Research/Shared_data_Mark/Projects/wetsus/20190905_experiment_2_setup/R")
workwd <- getwd()
filename <- "indirect_acting"
###Load pre-made functions
#git_dir <- "C:/Users/sterk009/OneDrive - Wageningen University & Research/Shared_data_Mark/Projects_R_Zone/Git"
#source(paste(git_dir,"/NEMA_functions/Loader.R",sep=""))
#Load custom functions
for(i in 1:length(dir("./Function_scripts/"))){
source(paste("./Function_scripts/", dir("./Function_scripts/")[i], sep = ""))
}
################################################################################
###Dependencies
################################################################################
install <- FALSE
#.libPaths(.libPaths("C:/Program Files/R/R-3.3.3/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("openxlsx")
install.packages("rmarkdown")
}
###load
library("colorspace")
library("RColorBrewer")
library(limma)
library(gridExtra)
library("VennDiagram")
library(openxlsx)
library("rmarkdown")
library(tidyverse)
################################################################################
###Plotting theme, colours
################################################################################
###Set plotting theme
presentation <- theme(axis.text.x = element_text(size=16, face="bold", color="black"),
axis.text.y = element_text(size=16, 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.grid.minor = element_blank(),
panel.grid.major = 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("black",brewer.pal(9,"Set1")[c(3,9)])
names(myColors) <- c("Cis","Trans","none")
colScale <- scale_colour_manual(name = "Condition",values = myColors)
fillScale <- scale_fill_manual(name = "Condition",values = myColors)
################################################################################################################################################################
###load all the data ################################################################################
################################################################################################################################################################
###normalization code, kept here for legacy, need to download archive to run it
# ###load targets and read-in data
# targets <- read.delim("./Data/Target_Wetsus_tox_2.txt"); head(targets)
#
# ###Normalize data
# ###
# rg.norm <- transcriptomics.agilent.norm(targets = targets,array_dir = paste(workwd,"/../../../Projects_data_reagents/MA_elegans",sep=""), save_dir = paste(workwd,"/Normalized/",sep=""))
# trans.int <- transcriptomics.transform.norm(rg.norm,filename=filename,save_dir=paste(workwd,"/Normalized/",sep=""))
#
# ###Checks
# correlsums <- transcriptomics.check.cor(trans.int,filename=filename,save_dir=paste(workwd,"/Normalized/",sep=""))
# transcriptomics.check.genes(trans.int,spot.id=agi.id$gene_public_name,filename=filename,save_dir=paste(workwd,"/Normalized/",sep=""))
#
# ###Make a list
# colnames.sep <- ":"
# colnames.names <- c("number","sample2","sample","treatment","concentration","repeat","batch")
#
# list.data <- transcriptomics.list.to.df(trans.int = trans.int, spot.id=Agilent.Ce.V2$SpotID,colnames.sep = colnames.sep, colnames.names = colnames.names) %>%
# dplyr::select(-sample2) # save(list.data, file="obj_list.data.Rdata")
# save(agi.id,file="obj_agi.id.Rdata")
load("./Data/obj_list.data.Rdata")
load("./Data/obj_agi.id.Rdata")
################################################################################################################################################################
###Batch correction ################################################################################
################################################################################################################################################################
###Batch correction
for(i in 1:45220){
data <- filter(list.data,SpotID==paste("AGIWUR",i,sep=""))%>%
mutate(dye=ifelse(number<17,"cy3","cy5"))
tmp <- summary(lm(log2_intensities~as.factor(batch)+dye+treatment,data=data))[[4]]
list.data <- mutate(list.data,log2_intensities=ifelse(SpotID==paste("AGIWUR",i,sep="") & batch==2,log2_intensities-tmp[2,1],log2_intensities)) %>%
mutate(log2_intensities=ifelse(SpotID==paste("AGIWUR",i,sep="") & batch==3,log2_intensities-tmp[3,1],log2_intensities)) %>%
mutate(log2_intensities=ifelse(SpotID==paste("AGIWUR",i,sep="") & number>12,log2_intensities-tmp[4,1],log2_intensities))
cat(i," ")
}
list.data <- select(list.data, -c(log2_ratio_mean,z_score)) %>%
group_by(SpotID) %>%
mutate(log2_ratio_mean=log2((2^log2_intensities)/(mean(2^log2_intensities)))) %>%
data.frame()
################################################################################
###Summarize data for array express
################################################################################
# ###Array express
# tmp <- separate(data.frame(tmp=c(as.character(unlist(targets$Cy3)),as.character(unlist(targets$Cy5)))),tmp,into=c("sample2","sample","treatment","concentration","repeat","batch"),sep=":")
# tmp <- cbind(FileName=rep(targets$FileName,times=2),label=rep(c("Cy3","Cy5"),each=nrow(targets)),tmp)
#
# Samples <- mutate(tmp,Organism="Caenorhabditis elegans",
# Age=55,
# treatment = ifelse(treatment=="B(a)P","benzo(a)pyrene",ifelse(treatment=="PCB1254","PCB mixture Aroclor 1254",ifelse(treatment=="TCDD","2,3,7,8-tetrachlorodibenzodioxin",ifelse(treatment=="AFB1","aflatoxin B1",ifelse(treatment=="Ctrl_DMSO","control",treatment))))),
# Organism_part="whole organism",
# Environmental_history=paste("grown for 31 hours at 20C after L1 arrest, subsequently underwent ",treatment," ",ifelse(treatment !="control",paste("(",concentration," μM) ",sep=""),""),"treatment for 24 hours.",sep=""),
# Developmental_stage="L4 larva",
# Sex="hermaphrodite",
# Diet="Escherichia coli OP50") %>%
# separate(FileName,into=c("folder","file"),sep="Raw/")
#
#
# write.xlsx(Samples,file="ArrayExpress_experiment.xlsx")
#
#
# data.nu <- dplyr::select(list.data,SpotID,number,log2_intensities) %>%
# mutate(number=paste("Sample",number,"normalised"),SpotID=as.numeric(gsub("AGIWUR","",SpotID))) %>%
# spread(key=number,value=log2_intensities) %>%
# arrange(SpotID)
#
# write.table(data.nu,file="log2_normalised.txt",sep="\t",quote=FALSE)
#
#
################################################################################################################################################################
###Principal component analysis (check only) ################################################################################
################################################################################################################################################################
###PCO
data.test <- filter(list.data) %>%
transcriptomics.df.to.list(variable.test=c("sample","treatment","concentration"),transformation.test="log2_ratio_mean"); #lapply(data.test,head)
tmp <- prcomp(data.test[[1]],scale. = TRUE)
summary(tmp)
plot.data <- as.data.frame(cbind(SampleID=names(data.test$sample),strain=names(data.test$strain),treatment=names(data.test$treatment),concentration_mM=names(data.test$concentration),tmp$rotation))
for(i in 4:ncol(plot.data)){
plot.data[,i] <- as.numeric(as.character(unlist(plot.data[,i])))
}
####Get legend function (not used now)
get_legend<-function(myggplot){
tmp <- ggplot_gtable(ggplot_build(myggplot))
leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
legend <- tmp$grobs[[leg]]
return(legend)
}
###colors
myColors <- c(brewer.pal(9,"Set1")[c(2,1,4,5)])
names(myColors) <- c("control","ENU","HCHO","MMS")
colScale <- scale_colour_manual(name = "Toxin",values = myColors)
fillScale <- scale_fill_manual(name = "Toxin",values = myColors)
ax.pl <- paste(colnames(summary(tmp)$importance)[1:6]," (",round(summary(tmp)$importance[2,1:6]*100,digits=1),"%)",sep="")
ax.plx <- ax.pl[c(1,3,5)]
ax.ply <- ax.pl[c(2,4,6)]
plot.nu <- NULL; plot.nu <- as.list(plot.nu)
plot.nu[[1]] <- ggplot(plot.data, aes(x=PC1,y=PC2,colour=treatment))
plot.nu[[2]] <- ggplot(plot.data, aes(x=PC3,y=PC4,colour=treatment))
plot.nu[[3]] <- ggplot(plot.data, aes(x=PC5,y=PC6,colour=treatment))
LgPl <- get_legend(plot.nu[[1]]+geom_point() + presentation)
for(i in 1:3){
plot.nu[[i]] <- plot.nu[[i]] + geom_point(aes(alpha=0.5)) + presentation + theme(legend.position = "none") + labs(x=ax.plx[i],y=ax.ply[i])
}
###Blank placeholder
BlPl <- ggplot()+geom_blank(aes(1,1))+ blank_theme
pdf(file=paste(workwd,"/20191202_PCO_analysis.pdf",sep=""),width=14,height=3)#,units="in",compression="lzw+p",res=600)
grid.arrange(plot.nu[[1]],plot.nu[[2]],plot.nu[[3]],LgPl,ncol=4)
dev.off()
################################################################################################################################################################
###Linear ################################################################################
################################################################################################################################################################
################################################################################
###R only step
################################################################################
###LM control vs B(a)P
data.test <- filter(list.data,treatment %in% c("Ctrl_DMSO","B(a)P")) %>% #,concentration_mM %in% c(0,5)) %>%
transcriptomics.df.to.list(variable.test=c("treatment","concentration"),transformation.test="log2_intensities"); lapply(data.test,head)
lm_BaP <- transcriptomics.lm(lm_data = data.test[[1]],lm_factors = as.numeric(names(data.test$concentration)))
sum(p.adjust(lm_BaP$model[,1],method="fdr") < 0.05); -log10(max(lm_BaP$model[p.adjust(lm_BaP$model[,1],method="fdr") < 0.05,1]))
###LM control vs PCB1254
data.test <- filter(list.data,treatment %in% c("Ctrl_DMSO","PCB1254")) %>% #,concentration_mM %in% c(0,5)) %>%
transcriptomics.df.to.list(variable.test=c("treatment","concentration"),transformation.test="log2_intensities"); lapply(data.test,head)
lm_PCB1254 <- transcriptomics.lm(lm_data = data.test[[1]],lm_factors = as.numeric(names(data.test$concentration)))
sum(p.adjust(lm_PCB1254$model[,1],method="fdr") < 0.05); -log10(max(lm_PCB1254$model[p.adjust(lm_PCB1254$model[,1],method="fdr") < 0.05,1]))
###LM control vs TCDD
data.test <- filter(list.data,treatment %in% c("Ctrl_DMSO","TCDD")) %>% #,concentration_mM %in% c(0,5)) %>%
transcriptomics.df.to.list(variable.test=c("treatment","concentration"),transformation.test="log2_intensities"); lapply(data.test,head)
lm_TCDD <- transcriptomics.lm(lm_data = data.test[[1]],lm_factors = as.numeric(names(data.test$concentration)))
sum(p.adjust(lm_TCDD$model[,1],method="fdr") < 0.05); -log10(max(lm_TCDD$model[p.adjust(lm_TCDD$model[,1],method="fdr") < 0.05,1]))
###LM control vs AFB1
data.test <- filter(list.data,treatment %in% c("Ctrl_DMSO","AFB1")) %>% #,concentration_mM %in% c(0,5)) %>%
transcriptomics.df.to.list(variable.test=c("treatment","concentration"),transformation.test="log2_intensities"); lapply(data.test,head)
lm_AFB1 <- transcriptomics.lm(lm_data = data.test[[1]],lm_factors = as.numeric(names(data.test$concentration)))
sum(p.adjust(lm_AFB1$model[,1],method="fdr") < 0.05); -log10(max(lm_AFB1$model[p.adjust(lm_AFB1$model[,1],method="fdr") < 0.05,1]))
model.data <- data.frame(cbind(select(agi.id,SpotID,Sequence:gene_sequence_name),
treatment=rep(c("BaP","PCB1254","TCDD","AFB1"),each=45220),
rbind(cbind(significance=lm_BaP$model[,1],effect=lm_BaP$model[,2]),
cbind(significance=lm_PCB1254$model[,1],effect=lm_PCB1254$model[,2]),
cbind(significance=lm_TCDD$model[,1],effect=lm_TCDD$model[,2]),
cbind(significance=lm_AFB1$model[,1],effect=lm_AFB1$model[,2]))))
for(i in c(1,2,3,6,7,8,9,10)){
model.data[,i] <- as.character(unlist(model.data[,i]))
}
for(i in c(4,5,11,12)){
model.data[,i] <- as.numeric(as.character(unlist(model.data[,i])))
}
model.data <- group_by(model.data,treatment) %>%
mutate(significance_FDR=p.adjust(significance,method="BH")) %>%
data.frame() %>%
mutate(sigcolour=ifelse(significance<0.00001,"-log10(p) > 5",ifelse(significance<0.0001,"-log10(p) > 4","-log10(p) < 4")),significance=-log10(significance)) %>%
arrange(treatment,SpotID)
###count
group_by(model.data,treatment) %>%
filter(significance > 5.0) %>%
filter(!duplicated(gene_WBID), !is.na(gene_WBID)) %>%
summarise(n= length(gene_WBID),thr=max(significance_FDR))
save(model.data, file=paste(workwd,"/obj_model.data.Rdata",sep=""))
################################################################################
###visualize: volcano plot & gene lists (not in paper)
################################################################################
load(file="./Data/obj_model.data.Rdata")
myColors <- c(brewer.pal(9,"Set1")[c(1,5,9)])
names(myColors) <- c("-log10(p) > 5","-log10(p) > 4","-log10(p) < 4")
colScale <- scale_colour_manual(name = "Significance",values = myColors)
###Volcano plot
plot.data <- filter(model.data,treatment=="PCB1254")
p1 <- ggplot(plot.data, aes(x=effect,y=significance,colour=sigcolour)) +
geom_point(alpha=0.4) + presentation + theme(legend.position = "none") + colScale + ylim(0,14) +
labs(x="Effect PCB1254",y="Significance (-log10(p))") + scale_x_continuous(breaks = seq(-1,1,by=0.1),labels=seq(-1,1,by=0.1))
plot.data <- filter(model.data,treatment=="TCDD")
p2 <- ggplot(plot.data, aes(x=effect,y=significance,colour=sigcolour)) + colScale + ylim(0,14) +
geom_point(alpha=0.4) + presentation + theme(legend.position = "none") +
labs(x="Effect TCDD",y="Significance (-log10(p))") + scale_x_continuous(breaks = seq(-1,1,by=0.1),labels=seq(-1,1,by=0.1))
plot.data <- filter(model.data,treatment=="AFB1")
p3 <- ggplot(plot.data, aes(x=effect,y=significance,colour=sigcolour)) + colScale + ylim(0,14) +
geom_point(alpha=0.4) + presentation + theme(legend.position = "none") +
labs(x="Effect AFB1",y="Significance (-log10(p))") + scale_x_continuous(breaks = seq(-1,1,by=0.1),labels=seq(-1,1,by=0.1))
plot.data <- filter(model.data,treatment=="BaP")
p4 <- ggplot(plot.data, aes(x=effect,y=significance,colour=sigcolour)) + colScale + ylim(0,14) +
geom_point(alpha=0.4) + presentation +
labs(x="Effect BaP",y="Significance (-log10(p))") + scale_x_continuous(breaks = seq(-1,1,by=0.1),labels=seq(-1,1,by=0.1))
###Plot
annotation.grobA <- title.grob <- textGrob(label = "A",x = unit(0, "lines"),y = unit(0, "lines"),hjust = 0, vjust = 0,gp = gpar(fontsize = 20,fontface="bold"))
annotation.grobB <- title.grob <- textGrob(label = "B",x = unit(0, "lines"),y = unit(0, "lines"),hjust = 0, vjust = 0,gp = gpar(fontsize = 20,fontface="bold"))
annotation.grobC <- title.grob <- textGrob(label = "C",x = unit(0, "lines"),y = unit(0, "lines"),hjust = 0, vjust = 0,gp = gpar(fontsize = 20,fontface="bold"))
annotation.grobD <- title.grob <- textGrob(label = "D",x = unit(0, "lines"),y = unit(0, "lines"),hjust = 0, vjust = 0,gp = gpar(fontsize = 20,fontface="bold"))
p1 <- arrangeGrob(p1,top=annotation.grobA)
p2 <- arrangeGrob(p2,top=annotation.grobB)
p3 <- arrangeGrob(p3,top=annotation.grobC)
p4 <- arrangeGrob(p4,top=annotation.grobD)
pdf(file=paste(workwd,"/20191203-Treatment_volcano.pdf",sep=""),width=12,height=4)#,units="in",compression="lzw+p",res=600)
grid.arrange(p1,p2,p3,p4,ncol=4,widths=c(1,1,1,1.5))
dev.off()
write.xlsx(model.data,file=paste(workwd,"/supplementary_table1-Linear_model_total.xlsx",sep=""))
write.xlsx(filter(model.data,significance>5),file=paste(workwd,"/Highly_significant.xlsx",sep=""))
################################################################################
###Data exploration (figures not in paper)
################################################################################
load(file="./Data/obj_model.data.Rdata")
###for text
filter(model.data,significance >5,gene_WBID != "",!is.na(gene_WBID)) %>%
group_by(treatment) %>%
summarise(n=length(unique(gene_WBID)),nspot=length(unique(SpotID)))
genes <- group_by(model.data,treatment,sign(effect)) %>%
filter(gene_WBID != "",!is.na(gene_WBID)) %>%
filter(significance==max(significance)) %>%
data.frame() %>%
filter(effect != 0)
data.plot <- merge(list.data,genes,by.x=2,by.y=1)
ggplot(data.plot,aes(x=concentration,y=log2_intensities)) +
geom_boxplot() + facet_grid(gene_public_name~treatment.x,scale="free",space="free_x")
genes <- group_by(model.data,treatment,sign(effect)) %>%
filter(gene_WBID != "",!is.na(gene_WBID)) %>%
filter(significance>5) %>%
data.frame() %>%
filter(effect != 0)
data.plot <- filter(list.data,SpotID %in% genes$SpotID) %>%
merge(dplyr::select(genes,c(SpotID,Sequence:gene_sequence_name))) %>%
group_by(gene_public_name,treatment,concentration) %>%
summarise(plotval=median(log2_ratio_mean)) %>%
group_by(gene_public_name) %>%
mutate(sort1=plotval[treatment=="Ctrl_DMSO"],sort2=plotval[treatment=="PCB1254" & concentration == 30],
sort3=plotval[treatment=="AFB1" & concentration == 30],sort4=plotval[treatment=="B(a)P" & concentration == 30]) %>%
arrange(sort1,sort2,sort3,sort4) %>%
data.frame() %>%
mutate(gene_public_name=factor(gene_public_name,levels = unique(gene_public_name))) %>%
mutate(treatment= factor(treatment,levels = c("Ctrl_DMSO","TCDD","AFB1","PCB1254","B(a)P")))
pdf(file=paste(workwd,"/20191203-concentration_plot.pdf",sep=""),width=6,height=12)#,units="in",compression="lzw+p",res=600)
ggplot(data.plot,aes(x=concentration,y=gene_public_name,fill=plotval)) +
geom_raster() + facet_grid(.~treatment,scale="free_x",space="free_x") +
scale_fill_gradientn(colours=brewer.pal(9,"RdYlBu"))
dev.off()
genes <- filter(model.data,grepl("cyp",gene_public_name))
data.plot <- filter(list.data,SpotID %in% genes$SpotID) %>%
merge(dplyr::select(genes,c(SpotID,Sequence:gene_sequence_name))) %>%
group_by(gene_public_name,treatment,concentration) %>%
summarise(plotval=median(log2_ratio_mean)) %>%
group_by(gene_public_name) %>%
mutate(sort1=plotval[treatment=="Ctrl_DMSO"],sort2=plotval[treatment=="PCB1254" & concentration == 30],
sort3=plotval[treatment=="AFB1" & concentration == 30],sort4=plotval[treatment=="B(a)P" & concentration == 30]) %>%
arrange(sort1,sort2,sort3,sort4) %>%
data.frame() %>%
mutate(gene_public_name=factor(gene_public_name,levels = unique(gene_public_name))) %>%
mutate(treatment= factor(treatment,levels = c("Ctrl_DMSO","TCDD","AFB1","PCB1254","B(a)P")))
pdf(file=paste(workwd,"/20191203-cyp_plot.pdf",sep=""),width=6,height=12)#,units="in",compression="lzw+p",res=600)
ggplot(data.plot,aes(x=concentration,y=gene_public_name,fill=plotval)) +
geom_raster() + facet_grid(.~treatment,scale="free_x",space="free_x") +
scale_fill_gradientn(colours=brewer.pal(9,"RdYlBu"))
dev.off()
################################################################################
###visualize: overlap treatments (figure in paper)
################################################################################
load(file="./Data/obj_model.data.Rdata")
###check overlap
data.plot <- filter(model.data,significance>4,treatment!="TCDD") %>%
group_by(treatment) %>%
filter(!duplicated(gene_WBID),!is.na(gene_WBID),gene_WBID != "") %>%
select(gene_WBID,treatment) %>%
data.frame()
data.plot <- split(data.plot[,1],data.plot[,2])
venn.tox <- venn.diagram(data.plot,filename = NULL, fill=brewer.pal(3,"Set1"),sub.fontfamily = "sans", main.fontfamily = "sans",cat.fontfamily="sans",fontfamily="sans")
pdf(file="20191203_figure6a-Venn_all.pdf",width=3,height=3)
grid.newpage()
grid.draw(venn.tox)
dev.off()
###check up and down (not in paper)
data.plot <- filter(model.data,significance>4,sign(effect)==-1) %>%
group_by(treatment) %>%
filter(!duplicated(gene_WBID),!is.na(gene_WBID),gene_WBID != "") %>%
select(gene_WBID,treatment) %>%
data.frame()
data.plot <- split(data.plot[,1],data.plot[,2])
venn.tox <- venn.diagram(data.plot,filename = NULL, fill=brewer.pal(2,"Set1"),sub.fontfamily = "sans", main.fontfamily = "sans",cat.fontfamily="sans",fontfamily="sans")
pdf(file="20191203_figure6b-Venn_down.pdf",width=3,height=3)
grid.newpage()
grid.draw(venn.tox)
dev.off()
data.plot <- filter(model.data,significance>4,sign(effect)==1) %>%
group_by(treatment) %>%
filter(!duplicated(gene_WBID),!is.na(gene_WBID),gene_WBID != "") %>%
select(gene_WBID,treatment) %>%
data.frame()
data.plot <- split(data.plot[,1],data.plot[,2])
venn.tox <- venn.diagram(data.plot,filename = NULL, fill=brewer.pal(4,"Set1"),sub.fontfamily = "sans", main.fontfamily = "sans",cat.fontfamily="sans",fontfamily="sans")
pdf(file="20191203_figure6c-Venn_up.pdf",width=3,height=3)
grid.newpage()
grid.draw(venn.tox)
dev.off()
overlap <- filter(model.data,significance>4) %>%
group_by(treatment) %>%
filter(!duplicated(gene_WBID),!is.na(gene_WBID),gene_WBID != "") %>%
group_by(gene_WBID) %>%
filter(length(gene_WBID)==3) %>%
data.frame()
write.xlsx(overlap,file="genes_in_3_treatments.xlsx")
################################################################################
###Wetsus toxicology experiment
################################################################################
###Set your work directory
setwd("E:/Nemawork/Projects/wetsus/20190905_experiment_2_setup/R/")
workwd <- getwd()
filename <- "wetsus_experiment2"
###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:/Nemawork/Projects_R_zone/Git"
source(paste(git_dir,"/Loader_git.R",sep=""))
###Set data locations
support_git_dir <- paste(git_dir,"/Kammenga_lab_experiments/Wetsus/supporting_files/",sep="")
################################################################################
###Dependencies
################################################################################
install <- FALSE
#.libPaths(.libPaths("C:/Program Files/R/R-3.3.3/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("openxlsx")
install.packages("rmarkdown")
}
###load
library("colorspace")
library("RColorBrewer")
library(limma)
library(gridExtra)
library("VennDiagram")
library(openxlsx)
library("rmarkdown")
library(tidyverse)
################################################################################
###Plotting theme, colours
################################################################################