Mercurial > repos > bornea > dotplot_runner
diff Dotplot_Release/R_dotPlot_nc.R @ 0:dfa3436beb67 draft
Uploaded
author | bornea |
---|---|
date | Fri, 29 Jan 2016 09:56:02 -0500 |
parents | |
children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Dotplot_Release/R_dotPlot_nc.R Fri Jan 29 09:56:02 2016 -0500 @@ -0,0 +1,140 @@ +#!/usr/bin/env Rscript + +args <- commandArgs(trailingOnly = TRUE) + +pheatmapj_loc <- paste(args[9],"pheatmap_j.R",sep="/") + +library('latticeExtra') +library('RColorBrewer') +library('grid') +library(reshape2) +source(pheatmapj_loc) + +data.file <- read.table("SC_data.txt", sep="\t", header=TRUE, row.names=1) ### import spectral count data +data.file2 <- read.table("FDR_data.txt", sep="\t", header=TRUE, row.names=1) ### import FDR count data +bait_l <- scan(args[4], what="") ### import bait list +if(args[5] == 0) prey_l <- scan(args[6], what="") ### import prey list +methd <- args[7] +dist_methd <- args[8] + +#setting parameters + +Sfirst=as.numeric(args[1]) #first FDR cutoff +Ssecond=as.numeric(args[2]) #second FDR cutoff +maxp=as.integer(args[3]) #maximum value for a spectral count + +#extract only needed data + +if(args[5] == 0){ + remove <- vector() + remove <- prey_l[prey_l %in% row.names(data.file)] + prey_l <- prey_l[prey_l %in% remove] + remove <- bait_l[bait_l %in% names(data.file)] + bait_l <- bait_l[bait_l %in% remove] + data.file <- data.file[prey_l, bait_l] + data.file2 <- data.file2[prey_l, bait_l] +} else{ + remove <- vector() + remove <- bait_l[bait_l %in% names(data.file)] + bait_l <- bait_l[bait_l %in% remove] + data.file <- data.file[, bait_l] + data.file2 <- data.file2[, bait_l] + prey_keep = apply(data.file2, 1, function(x) sum(x<=Sfirst) >= 1) + data.file <- data.file[prey_keep,] + data.file2 <- data.file2[prey_keep,] +} + +#determine bait and prey ordering + +y_ord=factor(names(data.file[1,]),levels=bait_l) + +if(args[5] == 0){ + x_ord=factor(rownames(data.file),levels=prey_l) +} else { + + data.file <- data.file[which(rowSums(data.file) > 0),] + dist_prey <- dist(as.matrix(data.file), method= dist_methd) + + if(methd == "ward"){ + dist_prey <- dist_prey^2 + } + + hc_prey <- hclust(dist_prey, method = methd) + + data.file = data.file[hc_prey$order, , drop = FALSE] + data.file2 = data.file2[hc_prey$order, , drop = FALSE] + + x_ord=factor(row.names(data.file), levels=row.names(data.file)) +} + +df<-data.frame(y=rep(y_ord, nrow(data.file)) + ,x=rep(x_ord, each=ncol(data.file)) + ,z1=as.vector(t(data.file)) # Circle color + ,z2=as.vector(t(data.file/apply(data.file,1,max))) # Circle size + ,z3=as.vector(t(data.file2)) # FDR +) + +df$z1[df$z1>maxp] <- maxp #maximum value for spectral count +df$z2[df$z2==0] <- NA +df$z3[df$z3>Ssecond] <- 0.05*maxp +df$z3[df$z3<=Ssecond & df$z3>Sfirst] <- 0.5*maxp +df$z3[df$z3<=Sfirst] <- 1*maxp +df$z4 <- df$z1 +df$z4[df$z4==0] <- 0 +df$z4[df$z4>0] <- 2.5 + +# The labeling for the colorkey + +labelat = c(0, maxp) +labeltext = c(0, maxp) + +# color scheme to use + +nmb.colors<-maxp +z.colors<-grey(rev(seq(0,0.9,0.9/nmb.colors))) #grayscale color scale + +#plot dotplot + +pl <- levelplot(z1~x*y, data=df + ,col.regions =z.colors #terrain.colors(100) + ,scales = list(x = list(rot = 90), y=list(cex=0.8), tck=0) # rotates X,Y labels and changes scale + ,colorkey = FALSE + #,colorkey = list(space="bottom", width=1.5, height=0.3, labels=list(at = labelat, labels = labeltext)) #put colorkey at top with my labeling scheme + ,xlab="Prey", ylab="Bait" + ,panel=function(x,y,z,...,col.regions){ + print(x) + z.c<-df$z1[ (df$x %in% as.character(x)) & (df$y %in% y)] + z.2<-df$z2[ (df$x %in% as.character(x)) & (df$y %in% y)] + z.3<-df$z3 + z.4<-df$z4 + panel.xyplot(x,y + ,as.table=TRUE + ,pch=21 # point type to use (circles in this case) + ,cex=((z.2-min(z.2,na.rm=TRUE))/(max(z.2,na.rm=TRUE)-min(z.2,na.rm=TRUE)))*3 #circle size + ,fill=z.colors[floor((z.c-min(z.c,na.rm=TRUE))*nmb.colors/(max(z.c,na.rm=TRUE)-min(z.c,na.rm=TRUE)))+1] # circle colors + ,col=z.colors[1+z.3] # border colors + ,lex=z.4 #border thickness + ) + } + #,main="Fold change" # graph main title + ) +if(ncol(data.file) > 4) ht=3.5+(0.36*((ncol(data.file)-1)-4)) else ht=3.5 +if(nrow(data.file) > 20) wd=8.25+(0.29*(nrow(data.file)-20)) else wd=5.7+(0.28*(nrow(data.file)-10)) +pdf("dotplot.pdf", onefile = FALSE, paper = "special", height = ht, width = wd, pointsize = 2) +print(pl) +dev.off() + +#plot heatmap + +heat_df <- acast(df, y~x, value.var="z1") +heat_df <- apply(heat_df, 2, rev) + +if(ncol(data.file) > 4) ht=3.5+(0.1*((ncol(data.file)-1)-4)) else ht=3.5 +if(nrow(data.file) > 20) wd=8.25+(0.1*(nrow(data.file)-20)) else wd=5+(0.1*(nrow(data.file)-10)) +pdf("heatmap_borders.pdf", onefile = FALSE, paper = "special", height = ht, width = wd, pointsize = 2) +pheatmap_j(heat_df, scale="none", border_color="black", border_width = 0.1, cluster_rows=FALSE, cluster_cols=FALSE, col=colorRampPalette(c("#FFFFFF", brewer.pal(9,"Blues")))(100)) +dev.off() + +pdf("heatmap_no_borders.pdf", onefile = FALSE, paper = "special", height = ht, width = wd, pointsize = 2) +pheatmap_j(heat_df, scale="none", border_color=NA, cluster_rows=FALSE, cluster_cols=FALSE, col=colorRampPalette(c("#FFFFFF", brewer.pal(9,"Blues")))(100)) +dev.off()