annotate Dotplot_Release/R_dotPlot.R @ 15:a6c081fcc0a9 draft

Uploaded
author bornea
date Wed, 16 Mar 2016 12:10:35 -0400
parents bc752a05f16d
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
3
bc752a05f16d Uploaded
bornea
parents:
diff changeset
1 #!/usr/bin/env Rscript
bc752a05f16d Uploaded
bornea
parents:
diff changeset
2
bc752a05f16d Uploaded
bornea
parents:
diff changeset
3 args <- commandArgs(trailingOnly = TRUE)
bc752a05f16d Uploaded
bornea
parents:
diff changeset
4
bc752a05f16d Uploaded
bornea
parents:
diff changeset
5 library('latticeExtra')
bc752a05f16d Uploaded
bornea
parents:
diff changeset
6 library('colorRamps')
bc752a05f16d Uploaded
bornea
parents:
diff changeset
7
bc752a05f16d Uploaded
bornea
parents:
diff changeset
8 data.file <- read.table("SC_data.txt", sep="\t", header=TRUE, row.names=1) ### import spectral count data
bc752a05f16d Uploaded
bornea
parents:
diff changeset
9 data.file2 <- read.table("FDR_data.txt", sep="\t", header=TRUE, row.names=1) ### import FDR count data
bc752a05f16d Uploaded
bornea
parents:
diff changeset
10 data.file3 <- read.table("clustered_matrix.txt", sep="\t", header=TRUE, row.names=1) ### import clustered matrix
bc752a05f16d Uploaded
bornea
parents:
diff changeset
11 data.file4 <- scan("singletons.txt", what="", sep="\n", strip.white=T) ### import singleton data
bc752a05f16d Uploaded
bornea
parents:
diff changeset
12
bc752a05f16d Uploaded
bornea
parents:
diff changeset
13 #setting parameters
bc752a05f16d Uploaded
bornea
parents:
diff changeset
14
bc752a05f16d Uploaded
bornea
parents:
diff changeset
15 Sfirst=as.numeric(args[1]) #first FDR cutoff
bc752a05f16d Uploaded
bornea
parents:
diff changeset
16 Ssecond=as.numeric(args[2]) #second FDR cutoff
bc752a05f16d Uploaded
bornea
parents:
diff changeset
17 maxp=as.integer(args[3]) #maximum value for a spectral count
bc752a05f16d Uploaded
bornea
parents:
diff changeset
18
bc752a05f16d Uploaded
bornea
parents:
diff changeset
19 #calculate column and row lengths
bc752a05f16d Uploaded
bornea
parents:
diff changeset
20
bc752a05f16d Uploaded
bornea
parents:
diff changeset
21 #determine bait and prey ordering
bc752a05f16d Uploaded
bornea
parents:
diff changeset
22
bc752a05f16d Uploaded
bornea
parents:
diff changeset
23 bait_levels=names(data.file3)
bc752a05f16d Uploaded
bornea
parents:
diff changeset
24 prey_levels=c(rownames(data.file3),data.file4)
bc752a05f16d Uploaded
bornea
parents:
diff changeset
25
bc752a05f16d Uploaded
bornea
parents:
diff changeset
26 x_ord=factor(row.names(data.file),levels=prey_levels)
bc752a05f16d Uploaded
bornea
parents:
diff changeset
27 y_ord=factor(names(data.file),levels=bait_levels)
bc752a05f16d Uploaded
bornea
parents:
diff changeset
28
bc752a05f16d Uploaded
bornea
parents:
diff changeset
29 df<-data.frame(y=rep(y_ord,nrow(data.file))
bc752a05f16d Uploaded
bornea
parents:
diff changeset
30 ,x=rep(x_ord, each=ncol(data.file))
bc752a05f16d Uploaded
bornea
parents:
diff changeset
31 ,z1=as.vector(t(data.file)) # Circle color
bc752a05f16d Uploaded
bornea
parents:
diff changeset
32 ,z2=as.vector(t(data.file/apply(data.file,1,max))) # Circle size
bc752a05f16d Uploaded
bornea
parents:
diff changeset
33 ,z3=as.vector(t(data.file2)) # FDR
bc752a05f16d Uploaded
bornea
parents:
diff changeset
34 )
bc752a05f16d Uploaded
bornea
parents:
diff changeset
35
bc752a05f16d Uploaded
bornea
parents:
diff changeset
36 df$z1[df$z1>maxp] <- maxp #maximum value for spectral count
bc752a05f16d Uploaded
bornea
parents:
diff changeset
37 df$z2[df$z2==0] <- NA
bc752a05f16d Uploaded
bornea
parents:
diff changeset
38 df$z3[df$z3>Ssecond] <- 0.05*maxp
bc752a05f16d Uploaded
bornea
parents:
diff changeset
39 df$z3[df$z3<=Ssecond & df$z3>Sfirst] <- 0.5*maxp
bc752a05f16d Uploaded
bornea
parents:
diff changeset
40 df$z3[df$z3<=Sfirst] <- 1*maxp
bc752a05f16d Uploaded
bornea
parents:
diff changeset
41 df$z4 <- df$z1
bc752a05f16d Uploaded
bornea
parents:
diff changeset
42 df$z4[df$z4==0] <- 0
bc752a05f16d Uploaded
bornea
parents:
diff changeset
43 df$z4[df$z4>0] <- 2.5
bc752a05f16d Uploaded
bornea
parents:
diff changeset
44
bc752a05f16d Uploaded
bornea
parents:
diff changeset
45 # The labeling for the colorkey
bc752a05f16d Uploaded
bornea
parents:
diff changeset
46
bc752a05f16d Uploaded
bornea
parents:
diff changeset
47 labelat = c(0, maxp)
bc752a05f16d Uploaded
bornea
parents:
diff changeset
48 labeltext = c(0, maxp)
bc752a05f16d Uploaded
bornea
parents:
diff changeset
49
bc752a05f16d Uploaded
bornea
parents:
diff changeset
50 # color scheme to use
bc752a05f16d Uploaded
bornea
parents:
diff changeset
51
bc752a05f16d Uploaded
bornea
parents:
diff changeset
52 nmb.colors<-maxp
bc752a05f16d Uploaded
bornea
parents:
diff changeset
53 z.colors<-grey(rev(seq(0,0.9,0.9/nmb.colors))) #grayscale color scale
bc752a05f16d Uploaded
bornea
parents:
diff changeset
54
bc752a05f16d Uploaded
bornea
parents:
diff changeset
55 #plot
bc752a05f16d Uploaded
bornea
parents:
diff changeset
56
bc752a05f16d Uploaded
bornea
parents:
diff changeset
57 pl <- levelplot(z1~x*y, data=df
bc752a05f16d Uploaded
bornea
parents:
diff changeset
58 ,col.regions =z.colors #terrain.colors(100)
bc752a05f16d Uploaded
bornea
parents:
diff changeset
59 ,scales = list(x = list(rot = 90), y=list(cex=0.8), tck=0) # rotates X,Y labels and changes scale
bc752a05f16d Uploaded
bornea
parents:
diff changeset
60 ,colorkey = FALSE
bc752a05f16d Uploaded
bornea
parents:
diff changeset
61 ,xlab="Prey", ylab="Bait"
bc752a05f16d Uploaded
bornea
parents:
diff changeset
62 ,panel=function(x,y,z,...,col.regions){
bc752a05f16d Uploaded
bornea
parents:
diff changeset
63 print(x)
bc752a05f16d Uploaded
bornea
parents:
diff changeset
64 z.c<-df$z1[ (df$x %in% as.character(x)) & (df$y %in% y)]
bc752a05f16d Uploaded
bornea
parents:
diff changeset
65 z.2<-df$z2[ (df$x %in% as.character(x)) & (df$y %in% y)]
bc752a05f16d Uploaded
bornea
parents:
diff changeset
66 z.3<-df$z3
bc752a05f16d Uploaded
bornea
parents:
diff changeset
67 z.4<-df$z4
bc752a05f16d Uploaded
bornea
parents:
diff changeset
68 panel.xyplot(x,y
bc752a05f16d Uploaded
bornea
parents:
diff changeset
69 ,as.table=TRUE
bc752a05f16d Uploaded
bornea
parents:
diff changeset
70 ,pch=21 # point type to use (circles in this case)
bc752a05f16d Uploaded
bornea
parents:
diff changeset
71 ,cex=((z.2-min(z.2,na.rm=TRUE))/(max(z.2,na.rm=TRUE)-min(z.2,na.rm=TRUE)))*3 #circle size
bc752a05f16d Uploaded
bornea
parents:
diff changeset
72 ,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
bc752a05f16d Uploaded
bornea
parents:
diff changeset
73 ,col=z.colors[1+z.3] # border colors
bc752a05f16d Uploaded
bornea
parents:
diff changeset
74 ,lex=z.4 #border thickness
bc752a05f16d Uploaded
bornea
parents:
diff changeset
75 )
bc752a05f16d Uploaded
bornea
parents:
diff changeset
76 }
bc752a05f16d Uploaded
bornea
parents:
diff changeset
77 #,main="Fold change" # graph main title
bc752a05f16d Uploaded
bornea
parents:
diff changeset
78 )
bc752a05f16d Uploaded
bornea
parents:
diff changeset
79 if(ncol(data.file) > 4) ht=3.5+(0.36*((ncol(data.file)-1)-4)) else ht=3.5
bc752a05f16d Uploaded
bornea
parents:
diff changeset
80 if(nrow(data.file) > 20) wd=8.25+(0.29*(nrow(data.file)-20)) else wd=5+(0.28*(nrow(data.file)-10))
bc752a05f16d Uploaded
bornea
parents:
diff changeset
81 pdf("dotplot.pdf", onefile = FALSE, paper = "special", height = ht, width = wd, pointsize = 2)
bc752a05f16d Uploaded
bornea
parents:
diff changeset
82 print(pl)
bc752a05f16d Uploaded
bornea
parents:
diff changeset
83 dev.off()