diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000000000000000000000000000000000000..5b6a0652566d10360493952aec6d4a4febc77083
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,4 @@
+.Rproj.user
+.Rhistory
+.RData
+.Ruserdata
diff --git a/Analysis_paper_wastewater_fingerprinting.R b/Analysis_paper_wastewater_fingerprinting.R
new file mode 100644
index 0000000000000000000000000000000000000000..633cb203daa8cb2efbd150219785991085aa2971
--- /dev/null
+++ b/Analysis_paper_wastewater_fingerprinting.R
@@ -0,0 +1,503 @@
+################################################################################
+###Wetsus toxicology experiment: indirect acting xenobiotics
+################################################################################
+
+
+ ###Set your work directory
+ setwd("C:/Users/sterk009/OneDrive - Wageningen University & Research/Shared_data_Mark/Projects/wetsus/20200604_experiment_3/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(brewer.pal(8,"Dark2")[c(8,4,3,7,5,1)])
+
+
+ names(myColors) <- c("Community","Hospital","Nursing home","WTP-Influent","WTP-Effluent","Surface water")
+
+ colScale <- scale_colour_manual(name = "Source",values = myColors)
+ fillScale <- scale_fill_manual(name = "Source",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_3.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","Sample_no","batch","sample","description")
+ #
+ # list.data <- transcriptomics.list.to.df(trans.int = trans.int, spot.id=Agilent.Ce.V2$SpotID,colnames.sep = colnames.sep, colnames.names = colnames.names) %>%
+ #
+ # 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() %>%
+ mutate(sample=gsub("SW.","Surface water",sample),
+ sample=gsub("C","Community",sample),
+ sample=gsub("E","WTP-Effluent",sample),
+ sample=gsub("H1","Hospital",sample),
+ sample=gsub("I","WTP-Influent",sample),
+ sample=gsub("N1","Nursing home",sample))
+
+
+
+################################################################################
+###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("Sample_no","batch","sample","description"),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 = description,
+ # Organism_part="whole organism",
+ # Environmental_history=paste("Synchronized L1 stage nematodes grown for 31 hours at 20C till L4 and were exposted to water samples from ",treatment," for 24 hours without food.",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")
+ #
+ # ###Annotare
+ # 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 (Figure 3) ################################################################################
+################################################################################################################################################################
+
+ ###PCO
+ data.test <- filter(list.data) %>%
+ transcriptomics.df.to.list(variable.test=c("Sample_no","sample","batch"),transformation.test="log2_ratio_mean"); #lapply(data.test,head)
+
+ tmp <- prcomp(data.test[[1]],scale. = TRUE)
+ summary(tmp)
+
+ tool.extrnames <- function(x,y){a <- unlist(strsplit(x,split=":"));return(a[seq(y,length(a),by=length(a)/length(x))])}
+
+ plot.data <- as.data.frame(cbind(SampleID=tool.extrnames(rownames(tmp$rotation),1),
+ sample=tool.extrnames(rownames(tmp$rotation),3),
+ batch=tool.extrnames(rownames(tmp$rotation),2),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)
+ }
+
+ #annotation <- filter(plot.data,treatment=="gld-1 RNAi",PC2 < -0.05, PC1>0.05)
+
+
+ p1 <- ggplot(plot.data, aes(x=PC1,y=PC2,colour=sample)) +
+ geom_point(size=3) + presentation + xlab("PCO 1 (42.4%)") + ylab("PCO 2 (14.2%)") +
+ theme(legend.position = "none") + colScale
+
+ p2 <- ggplot(plot.data, aes(x=PC3,y=PC4,colour=sample)) +
+ geom_point(size=3) + presentation + xlab("PCO 3 (9.0%)") + ylab("PCO 4 (8.2%)") +
+ theme(legend.position = "none") + colScale
+
+ p3 <- ggplot(plot.data, aes(x=PC5,y=PC6,colour=sample)) +
+ geom_point(size=3) + presentation + xlab("PCO 5 (5.4%)") + ylab("PCO 6 (3.7%)") +
+ colScale
+
+
+ pdf(file=paste(workwd,"/Figure3-PCA_new1-6.pdf",sep=""),width=14,height=3.5)#,units="in",compression="lzw+p",res=600)
+ grid.arrange(p1,p2,p3,widths=c(1,1,1.3))
+ dev.off()
+
+ p1 <- ggplot(plot.data, aes(x=PC1,y=PC2,colour=sample)) +
+ geom_point(size=3) + presentation + xlab("PCO 1 (42.4%)") + ylab("PCO 2 (14.2%)") +
+ colScale
+
+ pdf(file=paste(workwd,"/Figure3-PCA_new1-2.pdf",sep=""),width=5.5,height=3.5)#,units="in",compression="lzw+p",res=600)
+ grid.arrange(p1)
+ dev.off()
+
+
+################################################################################################################################################################
+###Linear model (Figure 2) ################################################################################
+################################################################################################################################################################
+
+
+ ################################################################################################################################################################
+ ###Run model (R only step) ##########################################################################
+ ################################################################################################################################################################
+
+
+ ###LM control vs H1
+ data.test <- filter(list.data,sample %in% c("Surface water","Hospital")) %>%
+ transcriptomics.df.to.list(variable.test=c("sample"),transformation.test="log2_intensities"); lapply(data.test,head)
+
+ lm_H1 <- transcriptomics.lm(lm_data = data.test[[1]],lm_factors = data.test$sample)
+ sum(p.adjust(lm_H1$model[,1],method="fdr") < 0.05); -log10(max(lm_H1$model[p.adjust(lm_H1$model[,1],method="fdr") < 0.05,1]))
+
+
+ ###LM control vs N1
+ data.test <- filter(list.data,sample %in% c("Surface water","Nursing home")) %>%
+ transcriptomics.df.to.list(variable.test=c("sample","concentration"),transformation.test="log2_intensities"); lapply(data.test,head)
+
+ lm_N1 <- transcriptomics.lm(lm_data = data.test[[1]],lm_factors = data.test$sample)
+ sum(p.adjust(lm_N1$model[,1],method="fdr") < 0.05); -log10(max(lm_N1$model[p.adjust(lm_N1$model[,1],method="fdr") < 0.05,1]))
+
+
+ ###LM control vs C
+ data.test <- filter(list.data,sample %in% c("Surface water","Community")) %>%
+ transcriptomics.df.to.list(variable.test=c("sample","concentration"),transformation.test="log2_intensities"); lapply(data.test,head)
+
+ lm_Cx <- transcriptomics.lm(lm_data = data.test[[1]],lm_factors = data.test$sample)
+ sum(p.adjust(lm_Cx$model[,1],method="fdr") < 0.05); -log10(max(lm_Cx$model[p.adjust(lm_Cx$model[,1],method="fdr") < 0.05,1]))
+
+
+ ###LM control vs I
+ data.test <- filter(list.data,sample %in% c("Surface water","WTP-Influent")) %>%
+ transcriptomics.df.to.list(variable.test=c("sample","concentration"),transformation.test="log2_intensities"); lapply(data.test,head)
+
+ lm_Ix <- transcriptomics.lm(lm_data = data.test[[1]],lm_factors = data.test$sample)
+ sum(p.adjust(lm_Ix$model[,1],method="fdr") < 0.05); -log10(max(lm_Ix$model[p.adjust(lm_Ix$model[,1],method="fdr") < 0.05,1]))
+
+ ###LM control vs E
+ data.test <- filter(list.data,sample %in% c("Surface water","WTP-Effluent")) %>%
+ transcriptomics.df.to.list(variable.test=c("sample","concentration"),transformation.test="log2_intensities"); lapply(data.test,head)
+
+ lm_Ex <- transcriptomics.lm(lm_data = data.test[[1]],lm_factors = data.test$sample)
+ sum(p.adjust(lm_Ex$model[,1],method="fdr") < 0.05); -log10(max(lm_Ex$model[p.adjust(lm_Ex$model[,1],method="fdr") < 0.05,1]))
+
+
+
+ model.data <- data.frame(cbind(select(agi.id,SpotID,FeatureNum,Sequence:gene_sequence_name),
+ sample=rep(c("Hospital","Nursing home","Community","WTP-Influent","WTP-Effluent"),each=45220),
+ rbind(cbind(significance=lm_H1$model[,1],effect=lm_H1$model[,2]),
+ cbind(significance=lm_N1$model[,1],effect=lm_N1$model[,2]),
+ cbind(significance=lm_Cx$model[,1],effect=lm_Cx$model[,2]),
+ cbind(significance=lm_Ix$model[,1],effect=lm_Ix$model[,2]),
+ cbind(significance=lm_Ex$model[,1],effect=lm_Ex$model[,2]))))
+
+
+ for(i in c(1,3,4,7,8,9,10,11)){
+ model.data[,i] <- as.character(unlist(model.data[,i]))
+ }
+ for(i in c(2,5,6,12,13)){
+ model.data[,i] <- as.numeric(as.character(unlist(model.data[,i])))
+ }
+
+ model.data <- group_by(model.data,sample) %>%
+ mutate(significance_FDR=p.adjust(significance,method="BH")) %>%
+ data.frame() %>%
+ mutate(effect = effect * -1) %>%
+ mutate(sigcolour=ifelse(significance<0.00001,"-log10(p) > 5",ifelse(significance<0.0001,"-log10(p) > 4","-log10(p) < 4")),significance=-log10(significance)) %>%
+ arrange(sample,FeatureNum)
+
+ ###for text
+ group_by(model.data,sample) %>%
+ filter(significance > 4.0) %>%
+ filter(!duplicated(gene_WBID), !is.na(gene_WBID)) %>%
+ summarise(n= length(gene_WBID),thr=max(significance_FDR))
+
+ save(model.data, file="./Data/obj_model.data.Rdata")
+ #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=""))
+
+
+ ################################################################################################################################################################
+ ###Volcano plot (figure 1) ##########################################################################
+ ################################################################################################################################################################
+
+ myColors2 <- c(brewer.pal(9,"Set1")[c(1,5,9)])
+ names(myColors2) <- c("-log10(p) > 5","-log10(p) > 4","-log10(p) < 4")
+
+ colScale2 <- scale_colour_manual(name = "Significance",values = myColors)
+
+
+ ###Volcano plot
+ plot.data <- filter(model.data,sample=="Hospital")
+ p1 <- ggplot(plot.data, aes(x=effect,y=significance,colour=sigcolour)) +
+ geom_point(alpha=0.4) + presentation + theme(legend.position = "none") + colScale2 + ylim(0,14) +
+ labs(x="Effect hospital",y="Significance (-log10(p))")
+
+ plot.data <- filter(model.data,sample=="Nursing home")
+ p2 <- ggplot(plot.data, aes(x=effect,y=significance,colour=sigcolour)) + colScale2 + ylim(0,14) +
+ geom_point(alpha=0.4) + presentation + theme(legend.position = "none") +
+ labs(x="Effect nursing home",y="Significance (-log10(p))")
+
+ plot.data <- filter(model.data,sample=="Community")
+ p3 <- ggplot(plot.data, aes(x=effect,y=significance,colour=sigcolour)) + colScale2 + ylim(0,14) +
+ geom_point(alpha=0.4) + presentation + theme(legend.position = "none") +
+ labs(x="Effect community",y="Significance (-log10(p))")
+
+ plot.data <- filter(model.data,sample=="WTP-Influent")
+ p4 <- ggplot(plot.data, aes(x=effect,y=significance,colour=sigcolour)) + colScale2 + ylim(0,14) +
+ geom_point(alpha=0.4) + presentation + theme(legend.position = "none") +
+ labs(x="Effect influent",y="Significance (-log10(p))")
+
+ plot.data <- filter(model.data,sample=="WTP-Effluent")
+ p5 <- ggplot(plot.data, aes(x=effect,y=significance,colour=sigcolour)) + colScale2 + ylim(0,14) +
+ geom_point(alpha=0.4) + presentation +
+ labs(x="Effect effluent",y="Significance (-log10(p))")
+
+
+
+ ###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"))
+ annotation.grobE <- title.grob <- textGrob(label = "E",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)
+ p5 <- arrangeGrob(p5,top=annotation.grobE)
+
+ pdf(file=paste(workwd,"/Figure1-Treatment_volcano_new.pdf",sep=""),width=16,height=4)#,units="in",compression="lzw+p",res=600)
+ grid.arrange(p1,p2,p3,p4,p5,ncol=5,widths=c(1,1,1,1,1.5))
+ dev.off()
+
+
+ ################################################################################################################################################################
+ ###Venn diagram (figure 4) ##########################################################################
+ ################################################################################################################################################################
+
+ ###check overlap
+ data.plot <- filter(model.data,significance>4,sample != "WTP-Effluent") %>%
+ group_by(sample) %>%
+ filter(!duplicated(gene_WBID),!is.na(gene_WBID),gene_WBID != "") %>%
+ select(gene_WBID,sample) %>%
+ data.frame()
+
+ data.plot <- split(data.plot[,1],data.plot[,2])
+
+
+ venn.tox <- venn.diagram(data.plot,filename = NULL, fill=myColors[1:4],sub.fontfamily = "sans", main.fontfamily = "sans",cat.fontfamily="sans",fontfamily="sans")
+
+ pdf(file="Figure4-Venn.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(sample) %>%
+ # filter(!duplicated(gene_WBID),!is.na(gene_WBID),gene_WBID != "") %>%
+ # select(gene_WBID,sample) %>%
+ # 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="20200621_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(sample) %>%
+ # filter(!duplicated(gene_WBID),!is.na(gene_WBID),gene_WBID != "") %>%
+ # select(gene_WBID,sample) %>%
+ # data.frame()
+ #
+ # data.plot <- split(data.plot[,1],data.plot[,2])
+ #
+ # venn.tox <- venn.diagram(data.plot,filename = NULL, fill=brewer.pal(5,"Set1"),sub.fontfamily = "sans", main.fontfamily = "sans",cat.fontfamily="sans",fontfamily="sans")
+ #
+ # pdf(file="20200621_figure6c-Venn_up.pdf",width=3,height=3)
+ # grid.newpage()
+ # grid.draw(venn.tox)
+ # dev.off()
+ #
+ # overlap <- filter(model.data,significance>4) %>%
+ # group_by(sample) %>%
+ # filter(!duplicated(gene_WBID),!is.na(gene_WBID),gene_WBID != "") %>%
+ # group_by(gene_WBID) %>%
+ # filter(length(gene_WBID)==4) %>%
+ # data.frame()
+ #
+ # write.xlsx(overlap,file="genes_in_4_treatments.xlsx")
+ #
+ #
+
+
+ ################################################################################################################################################################
+ ###DEGs (figure 2) ##########################################################################
+ ################################################################################################################################################################
+
+ overlap <- filter(model.data,significance>4) %>%
+ group_by(sample) %>%
+ filter(!duplicated(gene_WBID),!is.na(gene_WBID),gene_WBID != "") %>%
+ group_by(gene_WBID) %>%
+ filter(length(gene_WBID)==4) %>%
+ data.frame()
+
+ data.plot <- filter(list.data,SpotID %in% overlap$SpotID) %>%
+ merge(dplyr::select(overlap,c(SpotID,Sequence:gene_sequence_name))) %>%
+ group_by(gene_public_name,sample) %>%
+ summarise(Expression=median(log2_ratio_mean)) %>%
+ group_by(gene_public_name) %>%
+ mutate(sort1=Expression[sample=="Surface water"],sort2=Expression[sample=="Nursing home"],
+ sort3=Expression[sample=="Community"],sort4=Expression[sample=="Hospital"],sort5=Expression[sample=="WTP-Influent"]) %>%
+ arrange(sort1,sort2,sort3,sort4,sort5) %>%
+ data.frame() %>%
+ mutate(gene_public_name=factor(gene_public_name,levels = unique(gene_public_name))) %>%
+ mutate(sample= factor(sample,levels = c("Community","Hospital","Nursing home","WTP-Influent","WTP-Effluent","Surface water")))
+
+ pdf(file="Figure2B-Expression_overview.pdf",width=5,height=5)#,units="in",compression="lzw+p",res=600)
+ ggplot(data.plot,aes(x=sample,y=gene_public_name,fill=Expression)) +
+ geom_raster() + labs(x="Sample",y="Genes") + presentation +
+ scale_fill_gradientn(colours=brewer.pal(9,"RdYlBu")) +
+ theme(axis.text.y = element_blank(),axis.text.x = element_text(angle=45,hjust=1),axis.ticks.y = element_blank())
+ dev.off()
+
+ genes <- filter(model.data,grepl("cyp",gene_public_name))
+
+ data.plot <- filter(list.data,SpotID %in% genes$SpotID) %>%
+ merge(dplyr::select(overlap,c(SpotID,Sequence:gene_sequence_name))) %>%
+ group_by(gene_public_name,sample) %>%
+ summarise(Expression=median(log2_ratio_mean)) %>%
+ group_by(gene_public_name) %>%
+ mutate(sort1=Expression[sample=="Surface water"],sort2=Expression[sample=="Nursing home"],
+ sort3=Expression[sample=="Community"],sort4=Expression[sample=="Hospital"],sort5=Expression[sample=="WTP-Influent"]) %>%
+ arrange(sort1,sort2,sort3,sort4,sort5) %>%
+ data.frame() %>%
+ mutate(gene_public_name=factor(gene_public_name,levels = unique(gene_public_name))) %>%
+ mutate(sample= factor(sample,levels = c("Community","Hospital","Nursing home","WTP-Influent","WTP-Effluent","Surface water")))
+
+ pdf(file="Figure-Expression_cypgenes.pdf",width=5,height=5)#,units="in",compression="lzw+p",res=600)
+ ggplot(data.plot,aes(x=sample,y=gene_public_name,fill=Expression)) +
+ geom_raster() + labs(x="Sample",y="Genes") + presentation +
+ scale_fill_gradientn(colours=brewer.pal(9,"RdYlBu")) +
+ theme(axis.text.x = element_text(angle=45,hjust=1))
+ dev.off()
+
+
+
+
\ No newline at end of file
diff --git a/Archive/ArrayExpress_experiment.xlsx b/Archive/ArrayExpress_experiment.xlsx
new file mode 100644
index 0000000000000000000000000000000000000000..2accb3737149a16eee2b64477a8e06699e7cc549
Binary files /dev/null and b/Archive/ArrayExpress_experiment.xlsx differ
diff --git a/Archive/Wetsus_water_samples.R b/Archive/Wetsus_water_samples.R
new file mode 100644
index 0000000000000000000000000000000000000000..9f17530e7eb38192bc67ade7b2e1bd164f404ce3
--- /dev/null
+++ b/Archive/Wetsus_water_samples.R
@@ -0,0 +1,611 @@
+################################################################################
+###Wetsus toxicology experiment
+################################################################################
+
+
+###Set your work directory
+setwd("E:/Nemawork/Projects/wetsus/20200604_experiment_3/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,"/NEMA_functions/Loader.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
+################################################################################
+
+
+###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 ################################################################################
+
+ ###load targets and read-in data
+ #targets <- read.delim("W:/PROJECTS/NEMmasstorage/C.elegans/Data/MA/MA_elegans/Target_Wetsus_tox_2.txt"); head(targets)
+ targets <- read.delim("D:/MA/MA_elegans/Target_Wetsus_tox_3.txt"); head(targets)
+
+
+
+ ###Normalize data
+ ###
+ bad.arrays <- c(5,17)
+ rg.norm <- transcriptomics.agilent.norm(targets = targets,array_dir = "D:/MA/MA_elegans/", save_dir = paste(workwd,"/Normalized/",sep=""), bad.arrays = bad.arrays)
+ 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=""))
+
+ ###Batch correction
+ batch <- as.factor(unlist(strsplit(colnames(trans.int$log2_intensities),split=":"))[seq(3,5*ncol(trans.int$log2_intensities),by=5)])
+ dye <- as.factor(as.numeric(as.character(unlist(strsplit(colnames(trans.int$log2_intensities),split=":"))[seq(1,5*ncol(trans.int$log2_intensities),by=5)]))<13)
+ treatment <- as.factor(unlist(strsplit(colnames(trans.int$log2_intensities),split=":"))[seq(4,5*ncol(trans.int$log2_intensities),by=5)])
+ effects <- fast.summ.transcriptomics(lm(t(trans.int$log2_intensities)~batch+treatment))
+ effects <- effects[,-c(1:(ncol(effects)/2))]
+
+ batch <- paste("batch",batch,sep="")
+ output <- trans.int$log2_intensities
+ for(i in 1:length(unique(batch))){
+ if(sum(unique(batch)[i] %in% colnames(effects))==1){
+ output[,batch == unique(batch)[i]] <- output[,batch == unique(batch)[i]] - effects[,colnames(effects)==unique(batch)[i]]
+ }
+ }
+ output <- 2^output
+
+ trans.int$log2_intensities <- log2(output)
+ trans.int$log2_ratio_mean <- log2(output/apply(output,1,mean))
+ trans.int$z_score <- (output-apply(output,1,mean))/apply(output,1,sd)
+
+ ###Make a list
+ colnames.sep <- ":"
+ colnames.names <- c("number","Sample_no","batch","sample","description")
+
+ list.data <- transcriptomics.list.to.df(trans.int = trans.int, spot.id=Agilent.Ce.V2$SpotID,colnames.sep = colnames.sep, colnames.names = colnames.names)
+
+ save(list.data, file=paste(workwd,"/Normalized/list.data.Rdata",sep=""))
+
+
+
+################################################################################################################################################################
+###Principal component analysis ################################################################################
+
+ load(file=paste(workwd,"/Normalized/list.data.Rdata",sep=""))
+
+ ###PCO
+ data.test <- filter(list.data) %>%
+ transcriptomics.df.to.list(variable.test=c("Sample_no","sample","batch"),transformation.test="log2_ratio_mean"); #lapply(data.test,head)
+
+ tmp <- prcomp(data.test[[1]],scale. = TRUE)
+ summary(tmp)
+
+ tool.extrnames <- function(x,y){a <- unlist(strsplit(x,split=":"));return(a[seq(y,length(a),by=length(a)/length(x))])}
+
+ plot.data <- as.data.frame(cbind(SampleID=tool.extrnames(rownames(tmp$rotation),1),
+ sample=tool.extrnames(rownames(tmp$rotation),3),
+ batch=tool.extrnames(rownames(tmp$rotation),2),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)
+ }
+
+ #annotation <- filter(plot.data,treatment=="gld-1 RNAi",PC2 < -0.05, PC1>0.05)
+
+ plot.data <- mutate(plot.data,sample=gsub("SW.","SW",sample))
+
+ p1 <- ggplot(plot.data, aes(x=PC1,y=PC2,colour=sample)) +
+ geom_point() + presentation + xlab("PCO 1 (42.4%)") + ylab("PCO 2 (14.2%)") +
+ theme(legend.position = "none")
+
+ p2 <- ggplot(plot.data, aes(x=PC3,y=PC4,colour=sample)) +
+ geom_point() + presentation + xlab("PCO 3 (9.0%)") + ylab("PCO 4 (8.2%)") +
+ theme(legend.position = "none")
+
+ p3 <- ggplot(plot.data, aes(x=PC5,y=PC6,colour=sample)) +
+ geom_point() + presentation + xlab("PCO 5 (5.4%)") + ylab("PCO 6 (3.2%)")
+
+
+ pdf(file=paste(workwd,"/Figure2-PCA.pdf",sep=""),width=12,height=4)#,units="in",compression="lzw+p",res=600)
+ grid.arrange(p1,p2,p3,widths=c(1,1,1.3))
+ dev.off()
+
+################################################################################################################################################################
+###Linear model ################################################################################
+
+ load(file=paste(workwd,"/Normalized/list.data.Rdata",sep=""))
+
+
+ ###LM control vs H1
+ data.test <- filter(list.data,sample %in% c("SW1","SW2","SW3","H1")) %>%
+ mutate(sample=gsub("SW.","SW",sample)) %>%
+ transcriptomics.df.to.list(variable.test=c("sample"),transformation.test="log2_intensities"); lapply(data.test,head)
+
+ lm_H1 <- transcriptomics.lm(lm_data = data.test[[1]],lm_factors = data.test$sample)
+ sum(p.adjust(lm_H1$model[,1],method="fdr") < 0.05); -log10(max(lm_H1$model[p.adjust(lm_H1$model[,1],method="fdr") < 0.05,1]))
+
+
+ ###LM control vs N1
+ data.test <- filter(list.data,sample %in% c("SW1","SW2","SW3","N1")) %>%
+ mutate(sample=gsub("SW.","SW",sample)) %>%
+ transcriptomics.df.to.list(variable.test=c("sample","concentration"),transformation.test="log2_intensities"); lapply(data.test,head)
+
+ lm_N1 <- transcriptomics.lm(lm_data = data.test[[1]],lm_factors = data.test$sample)
+ sum(p.adjust(lm_N1$model[,1],method="fdr") < 0.05); -log10(max(lm_N1$model[p.adjust(lm_N1$model[,1],method="fdr") < 0.05,1]))
+
+
+ ###LM control vs C
+ data.test <- filter(list.data,sample %in% c("SW1","SW2","SW3","C")) %>%
+ mutate(sample=gsub("SW.","SW",sample)) %>%
+ transcriptomics.df.to.list(variable.test=c("sample","concentration"),transformation.test="log2_intensities"); lapply(data.test,head)
+
+ lm_Cx <- transcriptomics.lm(lm_data = data.test[[1]],lm_factors = data.test$sample)
+ sum(p.adjust(lm_Cx$model[,1],method="fdr") < 0.05); -log10(max(lm_Cx$model[p.adjust(lm_Cx$model[,1],method="fdr") < 0.05,1]))
+
+
+ ###LM control vs I
+ data.test <- filter(list.data,sample %in% c("SW1","SW2","SW3","I")) %>%
+ mutate(sample=gsub("SW.","SW",sample)) %>%
+ transcriptomics.df.to.list(variable.test=c("sample","concentration"),transformation.test="log2_intensities"); lapply(data.test,head)
+
+ lm_Ix <- transcriptomics.lm(lm_data = data.test[[1]],lm_factors = data.test$sample)
+ sum(p.adjust(lm_Ix$model[,1],method="fdr") < 0.05); -log10(max(lm_Ix$model[p.adjust(lm_Ix$model[,1],method="fdr") < 0.05,1]))
+
+ ###LM control vs E
+ data.test <- filter(list.data,sample %in% c("SW1","SW2","SW3","E")) %>%
+ mutate(sample=gsub("SW.","SW",sample)) %>%
+ transcriptomics.df.to.list(variable.test=c("sample","concentration"),transformation.test="log2_intensities"); lapply(data.test,head)
+
+ lm_Ex <- transcriptomics.lm(lm_data = data.test[[1]],lm_factors = data.test$sample)
+ sum(p.adjust(lm_Ex$model[,1],method="fdr") < 0.05); -log10(max(lm_Ex$model[p.adjust(lm_Ex$model[,1],method="fdr") < 0.05,1]))
+
+ model.data <- data.frame(cbind(select(agi.id,SpotID,FeatureNum,Sequence:gene_sequence_name),
+ sample=rep(c("H1","N1","C","I","E"),each=45220),
+ rbind(cbind(significance=lm_H1$model[,1],effect=lm_H1$model[,2]),
+ cbind(significance=lm_N1$model[,1],effect=lm_N1$model[,2]),
+ cbind(significance=lm_Cx$model[,1],effect=lm_Cx$model[,2]),
+ cbind(significance=lm_Ix$model[,1],effect=lm_Ix$model[,2]),
+ cbind(significance=lm_Ex$model[,1],effect=lm_Ex$model[,2]))))
+
+
+ for(i in c(1,3,4,7,8,9,10,11)){
+ model.data[,i] <- as.character(unlist(model.data[,i]))
+ }
+ for(i in c(2,5,6,12,13)){
+ model.data[,i] <- as.numeric(as.character(unlist(model.data[,i])))
+ }
+
+ model.data <- group_by(model.data,sample) %>%
+ mutate(significance_FDR=p.adjust(significance,method="BH")) %>%
+ data.frame() %>%
+ mutate(effect = effect * -1) %>%
+ mutate(sigcolour=ifelse(significance<0.00001,"-log10(p) > 5",ifelse(significance<0.0001,"-log10(p) > 4","-log10(p) < 4")),significance=-log10(significance)) %>%
+ arrange(sample,FeatureNum)
+
+
+ group_by(model.data,sample) %>%
+ filter(significance > 4.0) %>%
+ filter(!duplicated(gene_WBID), !is.na(gene_WBID)) %>%
+ summarise(n= length(gene_WBID),thr=max(significance_FDR))
+
+
+
+ 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,sample=="H1")
+ 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 hospital",y="Significance (-log10(p))")
+
+ plot.data <- filter(model.data,sample=="N1")
+ 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 nursing home",y="Significance (-log10(p))")
+
+ plot.data <- filter(model.data,sample=="C")
+ 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 community",y="Significance (-log10(p))")
+
+ plot.data <- filter(model.data,sample=="E")
+ p4 <- 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 effluent",y="Significance (-log10(p))")
+
+ plot.data <- filter(model.data,sample=="I")
+ p5 <- ggplot(plot.data, aes(x=effect,y=significance,colour=sigcolour)) + colScale + ylim(0,14) +
+ geom_point(alpha=0.4) + presentation +
+ labs(x="Effect Influent",y="Significance (-log10(p))")
+
+
+
+ ###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"))
+ annotation.grobE <- title.grob <- textGrob(label = "E",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)
+ p5 <- arrangeGrob(p5,top=annotation.grobE)
+
+ pdf(file=paste(workwd,"/20200621-Treatment_volcano.pdf",sep=""),width=16,height=4)#,units="in",compression="lzw+p",res=600)
+ grid.arrange(p1,p2,p3,p4,p5,ncol=5,widths=c(1,1,1,1,1.5))
+ dev.off()
+
+
+ save(model.data, file=paste(workwd,"/obj_model.data.Rdata",sep=""))
+ 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=""))
+
+
+ ###for text
+ filter(model.data,significance >4,gene_WBID != "",!is.na(gene_WBID)) %>%
+ group_by(sample) %>%
+ summarise(n=length(unique(gene_WBID)),nspot=length(unique(SpotID)))
+
+
+
+ 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))) %>%
+ mutate(sample=gsub("SW.","SW",sample)) %>%
+ group_by(gene_public_name,sample) %>%
+ summarise(plotval=median(log2_ratio_mean)) %>%
+ group_by(gene_public_name) %>%
+ mutate(sort1=plotval[sample=="SW"],sort2=plotval[sample=="N1"],
+ sort3=plotval[sample=="C"],sort4=plotval[sample=="H1"],sort5=plotval[sample=="I"]) %>%
+ arrange(sort1,sort2,sort3,sort4,sort5) %>%
+ data.frame() %>%
+ mutate(gene_public_name=factor(gene_public_name,levels = unique(gene_public_name))) %>%
+ mutate(sample= factor(sample,levels = c("SW","E","I","N1","H1","C")))
+
+ pdf(file=paste(workwd,"/20200621-cyp_plot.pdf",sep=""),width=6,height=12)#,units="in",compression="lzw+p",res=600)
+ ggplot(data.plot,aes(x=sample,y=gene_public_name,fill=plotval)) +
+ geom_raster() +
+ scale_fill_gradientn(colours=brewer.pal(9,"RdYlBu"))
+ dev.off()
+
+
+ ###check overlap
+ data.plot <- filter(model.data,significance>4) %>%
+ group_by(sample) %>%
+ filter(!duplicated(gene_WBID),!is.na(gene_WBID),gene_WBID != "") %>%
+ select(gene_WBID,sample) %>%
+ data.frame()
+
+ data.plot <- split(data.plot[,1],data.plot[,2])
+
+ venn.tox <- venn.diagram(data.plot,filename = NULL, fill=brewer.pal(5,"Set1"),sub.fontfamily = "sans", main.fontfamily = "sans",cat.fontfamily="sans",fontfamily="sans")
+
+ pdf(file="20200621_figure6a-Venn_all.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(sample) %>%
+ filter(!duplicated(gene_WBID),!is.na(gene_WBID),gene_WBID != "") %>%
+ select(gene_WBID,sample) %>%
+ 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="20200621_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(sample) %>%
+ filter(!duplicated(gene_WBID),!is.na(gene_WBID),gene_WBID != "") %>%
+ select(gene_WBID,sample) %>%
+ data.frame()
+
+ data.plot <- split(data.plot[,1],data.plot[,2])
+
+ venn.tox <- venn.diagram(data.plot,filename = NULL, fill=brewer.pal(5,"Set1"),sub.fontfamily = "sans", main.fontfamily = "sans",cat.fontfamily="sans",fontfamily="sans")
+
+ pdf(file="20200621_figure6c-Venn_up.pdf",width=3,height=3)
+ grid.newpage()
+ grid.draw(venn.tox)
+ dev.off()
+
+ overlap <- filter(model.data,significance>4) %>%
+ group_by(sample) %>%
+ filter(!duplicated(gene_WBID),!is.na(gene_WBID),gene_WBID != "") %>%
+ group_by(gene_WBID) %>%
+ filter(length(gene_WBID)==4) %>%
+ data.frame()
+
+ write.xlsx(overlap,file="genes_in_4_treatments.xlsx")
+
+
+ data.plot <- filter(list.data,SpotID %in% overlap$SpotID) %>%
+ merge(dplyr::select(overlap,c(SpotID,Sequence:gene_sequence_name))) %>%
+ mutate(sample=gsub("SW.","SW",sample)) %>%
+ group_by(gene_public_name,sample) %>%
+ summarise(plotval=median(log2_ratio_mean)) %>%
+ group_by(gene_public_name) %>%
+ mutate(sort1=plotval[sample=="SW"],sort2=plotval[sample=="N1"],
+ sort3=plotval[sample=="C"],sort4=plotval[sample=="H1"],sort5=plotval[sample=="I"]) %>%
+ arrange(sort1,sort2,sort3,sort4,sort5) %>%
+ data.frame() %>%
+ mutate(gene_public_name=factor(gene_public_name,levels = unique(gene_public_name))) %>%
+ mutate(sample= factor(sample,levels = c("SW","E","I","N1","H1","C")))
+
+ pdf(file=paste(workwd,"/20200621-overlap_plot.pdf",sep=""),width=6,height=24)#,units="in",compression="lzw+p",res=600)
+ ggplot(data.plot,aes(x=sample,y=gene_public_name,fill=plotval)) +
+ geom_raster() +
+ scale_fill_gradientn(colours=brewer.pal(9,"RdYlBu"))
+ dev.off()
+
+################################################################################################################################################################
+###enrichment ################################################################################
+ load(file="./obj_model.data.Rdata")
+
+ input <- filter(model.data,significance >4,!is.na(gene_WBID), gene_WBID != "")
+
+ enrichment <- enrich.easy(gene_list = input$gene_WBID,
+ gene_group = input$sample,
+ gene_total = Agilent.Ce.V2$gene_WBID,
+ database_dir = paste(git_dir,"/../R_Resources/enrich_Celegans",sep="")) %>%
+ mutate(FDR=p.adjust(Significance,method="fdr")) %>%
+ filter(Overlap > 2,FDR < 0.05)
+
+
+ write.xlsx(enrichment,file="Supplementary_file2-Enrichment_lod4.xlsx")
+
+
+
+
+################################################################################################################################################################
+###Associate with substances ################################################################################
+
+ load(file=paste(workwd,"/Normalized/list.data.Rdata",sep=""))
+ load(file="./obj_model.data.Rdata")
+
+ ICP <- read.xlsx(xlsxFile = "../data_samples/ICP_sample_Antoine.xlsx",sheet = 2) %>%
+ tidyr::gather(key=contaminant,value=value,-sample)
+ LCMS <- read.xlsx(xlsxFile = "../data_samples/LC-MS_sample_Antoine.xlsx",sheet = 3) %>%
+ tidyr::gather(key=contaminant,value=value,-sample)
+
+ contaminants <- rbind(ICP,LCMS) %>%
+ dplyr::mutate(value=ifelse(is.na(value),0,value)) %>%
+ dplyr::filter(!contaminant %in% c("Cadmium","Cefotaxime","Sulfadoxin")) %>%
+ tidyr::spread(key=contaminant,value=value)
+
+ list.test <- merge(list.data,contaminants) %>%
+ dplyr::mutate(numberspot=as.numeric(as.character(gsub("AGIWUR","",SpotID)))) %>%
+ dplyr::arrange(number,numberspot) %>%
+ dplyr::select(-numberspot)
+
+ output <- NULL
+ for(i in 2:ncol(contaminants)){
+
+ ###LM control vs H1
+ data.test <- filter(list.test) %>%
+ transcriptomics.df.to.list(variable.test=colnames(contaminants)[i],transformation.test="log2_intensities"); lapply(data.test,head)
+
+
+ vals <- as.numeric(as.character(unlist(names(data.test[[2]]))))
+
+ lm_nu <- transcriptomics.lm(lm_data = data.test[[1]],lm_factors = vals)$model
+
+ sum(p.adjust(lm_nu[,1],method="fdr") < 0.05); min(lm_nu[,1])
+
+ output <- rbind(output,
+ cbind(substance=colnames(contaminants)[i],spotnames=rownames(lm_nu),lm_nu))
+ }
+
+ output <- data.frame(output)
+ output[,1] <- as.character(unlist(output[,1]))
+ output[,2] <- as.character(unlist(output[,2]))
+ output[,3] <- as.numeric(as.character(unlist(output[,3])))
+ output[,4] <- as.numeric(as.character(unlist(output[,4])))
+
+
+ output <- mutate(output,fdr=p.adjust(significance_x,method="fdr")) %>%
+ merge(dplyr::select(agi.id,SpotID,FeatureNum,Sequence:gene_sequence_name),by.x = 2,by.y=1) %>%
+ rename(significance = significance_x, effect = effect_x) %>%
+ arrange(substance,FeatureNum) %>%
+ mutate(significance=-log10(significance))
+
+ table_chems <- group_by(output,substance) %>%
+ summarise(n_spots=sum(significance>4),n_genes=length(unique(gene_public_name[significance>4]))) %>%
+ data.frame()
+
+ write.xlsx(table_chems,file = "table_LCMS_ICP.xlsx")
+ write.xlsx(contaminants,file = "modelled_contaminants.xlsx")
+ write.xlsx(output,file="linear_model_LCMS_ICP.xlsx")
+
+ model.data.chems <- output
+
+ save(model.data.chems,file="obj_model.data.chems.Rdata")
+
+
+
+
+
+ data.plot <- filter(model.data,significance>4,sample !="E") %>%
+ group_by(sample) %>%
+ filter(!duplicated(gene_WBID),!is.na(gene_WBID),gene_WBID != "") %>%
+ select(gene_WBID,sample) %>%
+ data.frame()
+
+ for(i in 1:nrow(table_chems)){
+
+ data.plot2 <- filter(model.data.chems,substance==table_chems[i,1],significance>4) %>%
+ filter(!duplicated(gene_WBID),!is.na(gene_WBID),gene_WBID != "") %>%
+ dplyr::select(gene_WBID,substance) %>%
+ rename(sample=substance)
+
+ data.plot2 <- rbind(data.plot,data.plot2)
+ datavenn <- mutate(data.plot2,value=1) %>%
+ spread(key=sample,value=value) %>%
+ merge(dplyr::filter(dplyr::select(agi.id,chromosome:gene_sequence_name),!duplicated(gene_WBID)),by.x = 1,by.y=5)
+
+ data.plot2 <- split(data.plot2[,1],data.plot2[,2])
+
+ venn.tox <- venn.diagram(data.plot2,filename = NULL, fill=brewer.pal(5,"Set1"),sub.fontfamily = "sans", main.fontfamily = "sans",cat.fontfamily="sans",fontfamily="sans")
+
+ pdf(file=paste("20200712_figure7-Venn_all_",table_chems[i,1],".pdf",sep=""),width=4,height=4)
+ grid.newpage()
+ grid.draw(venn.tox)
+ dev.off()
+
+ write.xlsx(datavenn,file=paste("20201023_Venn_data_",table_chems[i,1],".xlsx",sep=""))
+
+ }
+
+
+ data.plot <- filter(model.data.chems,significance>4,substance %in% c("Benzalkonium.Chloride.C12")) %>%
+ merge(list.test,by.x=1,by.y=3) %>%
+ arrange(effect) %>%
+ mutate(Benzalkonium.Chloride.C12=round(Benzalkonium.Chloride.C12,digits=2)) %>%
+ mutate(gene_public_name = factor(gene_public_name,levels = unique(gene_public_name)),
+ Benzalkonium.Chloride.C12=factor(Benzalkonium.Chloride.C12,levels = unique(sort(Benzalkonium.Chloride.C12)))) %>%
+ filter(!is.na(gene_WBID),gene_WBID != "") %>%
+ group_by(Benzalkonium.Chloride.C12,gene_public_name) %>%
+ summarise(log2_ratio_mean=mean(log2_ratio_mean))
+
+ pdf("Exampe_Benzalkonium.Chloride.C12.pdf",width=8,height=16)
+ ggplot(data.plot,aes(x=Benzalkonium.Chloride.C12,y=gene_public_name,fill=log2_ratio_mean)) +
+ geom_raster() + scale_fill_gradientn(colours=brewer.pal(11,"RdYlBu")) + presentation
+ dev.off()
+
+
+################################################################################################################################################################
+###enrichment ################################################################################
+ load(file="./obj_model.data.chems.Rdata")
+
+ input <- filter(model.data.chems,significance >4,!is.na(gene_WBID), gene_WBID != "")
+
+ enrichment <- enrich.easy(gene_list = input$gene_WBID,
+ gene_group = input$substance,
+ gene_total = Agilent.Ce.V2$gene_WBID,
+ database_dir = paste(git_dir,"/../R_Resources/enrich_Celegans",sep="")) %>%
+ mutate(FDR=p.adjust(Significance,method="fdr")) %>%
+ filter(Overlap > 2,FDR < 0.05)
+
+
+ write.xlsx(enrichment,file="Supplementary_file3-Enrichment_lod4_substances.xlsx")
+
+
+################################################################################################################################################################
+###correlation analysis ################################################################################
+
+
+ data.plot <- data.matrix(contaminants[,-1])
+ rownames(data.plot) <- contaminants[,1]
+
+ data.cor <- cor(data.matrix(data.plot[,-1]))
+ colscale <- color.scale(colors.use=brewer.pal(11,"PiYG"),y.val=as.numeric(data.cor))
+
+ ###Save the figure
+ pdf(file="Correlation_contaminants.pdf",width=16,height=8)#,units="in",compression="lzw+p",res=600)
+ heatmap(data.cor,scale="none",col=colscale,margins=c(10,8))
+ par(fig=c(0,0.3,0.5,1),new=T)
+ scale.legend(input=as.numeric(data.cor),start.zero=F,col.scale=colscale,line.col="black")
+ dev.off()
+
+
+
+################################################################################################################################################################
+###correlation analysis ################################################################################
+
+
+
+ data.plot <- filter(model.data.chems,significance>4,substance %in% c("Benzalkonium.Chloride.C12")) %>%
+ merge(list.test,by.x=1,by.y=3) %>