annotate Dotplot_Release/R_dotPlot_nc.R @ 0:dfa3436beb67 draft

Uploaded
author bornea
date Fri, 29 Jan 2016 09:56:02 -0500
parents
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
1 #!/usr/bin/env Rscript
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
2
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
3 args <- commandArgs(trailingOnly = TRUE)
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
4
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
5 pheatmapj_loc <- paste(args[9],"pheatmap_j.R",sep="/")
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
6
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
7 library('latticeExtra')
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
8 library('RColorBrewer')
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
9 library('grid')
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
10 library(reshape2)
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
11 source(pheatmapj_loc)
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
12
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
13 data.file <- read.table("SC_data.txt", sep="\t", header=TRUE, row.names=1) ### import spectral count data
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
14 data.file2 <- read.table("FDR_data.txt", sep="\t", header=TRUE, row.names=1) ### import FDR count data
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
15 bait_l <- scan(args[4], what="") ### import bait list
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
16 if(args[5] == 0) prey_l <- scan(args[6], what="") ### import prey list
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
17 methd <- args[7]
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
18 dist_methd <- args[8]
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
19
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
20 #setting parameters
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
21
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
22 Sfirst=as.numeric(args[1]) #first FDR cutoff
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
23 Ssecond=as.numeric(args[2]) #second FDR cutoff
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
24 maxp=as.integer(args[3]) #maximum value for a spectral count
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
25
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
26 #extract only needed data
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
27
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
28 if(args[5] == 0){
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
29 remove <- vector()
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
30 remove <- prey_l[prey_l %in% row.names(data.file)]
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
31 prey_l <- prey_l[prey_l %in% remove]
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
32 remove <- bait_l[bait_l %in% names(data.file)]
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
33 bait_l <- bait_l[bait_l %in% remove]
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
34 data.file <- data.file[prey_l, bait_l]
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
35 data.file2 <- data.file2[prey_l, bait_l]
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
36 } else{
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
37 remove <- vector()
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
38 remove <- bait_l[bait_l %in% names(data.file)]
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
39 bait_l <- bait_l[bait_l %in% remove]
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
40 data.file <- data.file[, bait_l]
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
41 data.file2 <- data.file2[, bait_l]
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
42 prey_keep = apply(data.file2, 1, function(x) sum(x<=Sfirst) >= 1)
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
43 data.file <- data.file[prey_keep,]
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
44 data.file2 <- data.file2[prey_keep,]
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
45 }
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
46
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
47 #determine bait and prey ordering
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
48
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
49 y_ord=factor(names(data.file[1,]),levels=bait_l)
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
50
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
51 if(args[5] == 0){
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
52 x_ord=factor(rownames(data.file),levels=prey_l)
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
53 } else {
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
54
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
55 data.file <- data.file[which(rowSums(data.file) > 0),]
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
56 dist_prey <- dist(as.matrix(data.file), method= dist_methd)
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
57
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
58 if(methd == "ward"){
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
59 dist_prey <- dist_prey^2
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
60 }
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
61
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
62 hc_prey <- hclust(dist_prey, method = methd)
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
63
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
64 data.file = data.file[hc_prey$order, , drop = FALSE]
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
65 data.file2 = data.file2[hc_prey$order, , drop = FALSE]
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
66
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
67 x_ord=factor(row.names(data.file), levels=row.names(data.file))
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
68 }
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
69
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
70 df<-data.frame(y=rep(y_ord, nrow(data.file))
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
71 ,x=rep(x_ord, each=ncol(data.file))
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
72 ,z1=as.vector(t(data.file)) # Circle color
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
73 ,z2=as.vector(t(data.file/apply(data.file,1,max))) # Circle size
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
74 ,z3=as.vector(t(data.file2)) # FDR
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
75 )
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
76
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
77 df$z1[df$z1>maxp] <- maxp #maximum value for spectral count
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
78 df$z2[df$z2==0] <- NA
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
79 df$z3[df$z3>Ssecond] <- 0.05*maxp
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
80 df$z3[df$z3<=Ssecond & df$z3>Sfirst] <- 0.5*maxp
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
81 df$z3[df$z3<=Sfirst] <- 1*maxp
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
82 df$z4 <- df$z1
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
83 df$z4[df$z4==0] <- 0
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
84 df$z4[df$z4>0] <- 2.5
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
85
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
86 # The labeling for the colorkey
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
87
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
88 labelat = c(0, maxp)
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
89 labeltext = c(0, maxp)
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
90
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
91 # color scheme to use
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
92
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
93 nmb.colors<-maxp
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
94 z.colors<-grey(rev(seq(0,0.9,0.9/nmb.colors))) #grayscale color scale
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
95
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
96 #plot dotplot
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
97
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
98 pl <- levelplot(z1~x*y, data=df
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
99 ,col.regions =z.colors #terrain.colors(100)
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
100 ,scales = list(x = list(rot = 90), y=list(cex=0.8), tck=0) # rotates X,Y labels and changes scale
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
101 ,colorkey = FALSE
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
102 #,colorkey = list(space="bottom", width=1.5, height=0.3, labels=list(at = labelat, labels = labeltext)) #put colorkey at top with my labeling scheme
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
103 ,xlab="Prey", ylab="Bait"
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
104 ,panel=function(x,y,z,...,col.regions){
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
105 print(x)
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
106 z.c<-df$z1[ (df$x %in% as.character(x)) & (df$y %in% y)]
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
107 z.2<-df$z2[ (df$x %in% as.character(x)) & (df$y %in% y)]
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
108 z.3<-df$z3
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
109 z.4<-df$z4
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
110 panel.xyplot(x,y
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
111 ,as.table=TRUE
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
112 ,pch=21 # point type to use (circles in this case)
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
113 ,cex=((z.2-min(z.2,na.rm=TRUE))/(max(z.2,na.rm=TRUE)-min(z.2,na.rm=TRUE)))*3 #circle size
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
114 ,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
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
115 ,col=z.colors[1+z.3] # border colors
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
116 ,lex=z.4 #border thickness
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
117 )
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
118 }
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
119 #,main="Fold change" # graph main title
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
120 )
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
121 if(ncol(data.file) > 4) ht=3.5+(0.36*((ncol(data.file)-1)-4)) else ht=3.5
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
122 if(nrow(data.file) > 20) wd=8.25+(0.29*(nrow(data.file)-20)) else wd=5.7+(0.28*(nrow(data.file)-10))
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
123 pdf("dotplot.pdf", onefile = FALSE, paper = "special", height = ht, width = wd, pointsize = 2)
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
124 print(pl)
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
125 dev.off()
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
126
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
127 #plot heatmap
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
128
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
129 heat_df <- acast(df, y~x, value.var="z1")
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
130 heat_df <- apply(heat_df, 2, rev)
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
131
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
132 if(ncol(data.file) > 4) ht=3.5+(0.1*((ncol(data.file)-1)-4)) else ht=3.5
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
133 if(nrow(data.file) > 20) wd=8.25+(0.1*(nrow(data.file)-20)) else wd=5+(0.1*(nrow(data.file)-10))
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
134 pdf("heatmap_borders.pdf", onefile = FALSE, paper = "special", height = ht, width = wd, pointsize = 2)
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
135 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))
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
136 dev.off()
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
137
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
138 pdf("heatmap_no_borders.pdf", onefile = FALSE, paper = "special", height = ht, width = wd, pointsize = 2)
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
139 pheatmap_j(heat_df, scale="none", border_color=NA, cluster_rows=FALSE, cluster_cols=FALSE, col=colorRampPalette(c("#FFFFFF", brewer.pal(9,"Blues")))(100))
dfa3436beb67 Uploaded
bornea
parents:
diff changeset
140 dev.off()