Mercurial > repos > bornea > prohits_dotplot_generator

comparison Dotplot_Release/R_dotPlot_hc.R @ 3:bc752a05f16d draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author bornea
date Tue, 15 Mar 2016 15:25:15 -0400
parents
children
comparison
equal deleted inserted replaced
2:cfe2edb1c5d8 3:bc752a05f16d
1 #!/usr/bin/env Rscript
2
3 args <- commandArgs(trailingOnly = TRUE)
4
5 pheatmapj_loc <- paste(args[6],"pheatmap_j.R",sep="/")
6 heatmap2j_loc <- paste(args[6],"heatmap_2j.R",sep="/")
7
8 library('latticeExtra')
9 library('RColorBrewer')
10 library('grid')
11 library(reshape2)
12 library('gplots')
13 library('gtools')
14 source(pheatmapj_loc)
15 source(heatmap2j_loc)
16
17 data.file <- read.table("SC_data.txt", sep="\t", header=TRUE, row.names=1) ### import spectral count data
18 data.file2 <- read.table("FDR_data.txt", sep="\t", header=TRUE, row.names=1) ### import FDR count data
19
20 #setting parameters
21
22 Sfirst=as.numeric(args[1]) #first FDR cutoff
23 Ssecond=as.numeric(args[2]) #second FDR cutoff
24 maxp=as.integer(args[3]) #maximum value for a spectral count
25 methd <- args[4]
26 dist_methd <- args[5]
27
28 #determine bait and prey ordering
29
30 dist_bait <- dist(as.matrix(t(data.file)), method= dist_methd) # "euclidean", "maximum", "manhattan", "canberra", "binary" or "minkowski"
31 dist_prey <- dist(as.matrix(data.file), method= dist_methd)
32
33 if(methd == "ward"){
34 dist_bait <- dist_bait^2 #comment out this line and the next if not using Ward's method of clustering
35 dist_prey <- dist_prey^2
36 }
37
38 hc_bait <- hclust(dist_bait, method = methd) # method = "average", "single", "complete", "ward", "mcquitty", "median" or "centroid"
39 hc_prey <- hclust(dist_prey, method = methd)
40
41 data.file = data.file[hc_prey$order, , drop = FALSE]
42 data.file = data.file[, hc_bait$order, drop = FALSE]
43 data.file2 = data.file2[hc_prey$order, , drop = FALSE]
44 data.file2 = data.file2[, hc_bait$order, drop = FALSE]
45
46 x_ord=factor(row.names(data.file), levels=row.names(data.file))
47 y_ord=factor(names(data.file[1,]), levels=names(data.file[1,]))
48
49 df<-data.frame(y=rep(y_ord, nrow(data.file))
50 ,x=rep(x_ord, each=ncol(data.file))
51 ,z1=as.vector(t(data.file)) # Circle color
52 ,z2=as.vector(t(data.file/apply(data.file,1,max))) # Circle size
53 ,z3=as.vector(t(data.file2)) # FDR
54 )
55
56 df$z1[df$z1>maxp] <- maxp #maximum value for spectral count
57 df$z2[df$z2==0] <- NA
58 df$z3[df$z3>Ssecond] <- 0.05*maxp
59 df$z3[df$z3<=Ssecond & df$z3>Sfirst] <- 0.5*maxp
60 df$z3[df$z3<=Sfirst] <- 1*maxp
61 df$z4 <- df$z1
62 df$z4[df$z4==0] <- 0
63 df$z4[df$z4>0] <- 2.5
64
65 # The labeling for the colorkey
66
67 labelat = c(0, maxp)
68 labeltext = c(0, maxp)
69
70 # color scheme to use
71
72 nmb.colors<-maxp
73 z.colors<-grey(rev(seq(0,0.9,0.9/nmb.colors))) #grayscale color scale
74
75 #plot dotplot
76
77 pl <- levelplot(z1~x*y, data=df
78 ,col.regions =z.colors #terrain.colors(100)
79 ,scales = list(x = list(rot = 90), y=list(cex=0.8), tck=0) # rotates X,Y labels and changes scale
80 ,colorkey = FALSE
81 #,colorkey = list(space="bottom", width=1.5, height=0.3, labels=list(at = labelat, labels = labeltext)) #put colorkey at top with my labeling scheme
82 ,xlab="Prey", ylab="Bait"
83 ,panel=function(x,y,z,...,col.regions){
84 print(x)
85 z.c<-df$z1[ (df$x %in% as.character(x)) & (df$y %in% y)]
86 z.2<-df$z2[ (df$x %in% as.character(x)) & (df$y %in% y)]
87 z.3<-df$z3
88 z.4<-df$z4
89 panel.xyplot(x,y
90 ,as.table=TRUE
91 ,pch=21 # point type to use (circles in this case)
92 ,cex=((z.2-min(z.2,na.rm=TRUE))/(max(z.2,na.rm=TRUE)-min(z.2,na.rm=TRUE)))*3 #circle size
93 ,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
94 ,col=z.colors[1+z.3] # border colors
95 ,lex=z.4 #border thickness
96 )
97 }
98 #,main="Fold change" # graph main title
99 )
100 if(ncol(data.file) > 4) ht=3.5+(0.36*((ncol(data.file)-1)-4)) else ht=3.5
101 if(nrow(data.file) > 20) wd=8.25+(0.29*(nrow(data.file) -20)) else wd=5.7+(0.28*(nrow(data.file) -10))
102 pdf("dotplot.pdf", onefile = FALSE, paper = "special", height = ht, width = wd, pointsize = 2)
103 print(pl)
104 dev.off()
105
106 #plot bait vs prey heatmap
107
108 heat_df <- acast(df, y~x, value.var="z1")
109 heat_df <- apply(heat_df, 2, rev)
110
111 if(ncol(data.file) > 4) ht=3.5+(0.1*((ncol(data.file)-1)-4)) else ht=3.5
112 if(nrow(data.file) > 20) wd=8.25+(0.1*(nrow(data.file)-20)) else wd=5+(0.1*(nrow(data.file)-10))
113 pdf("heatmap_borders.pdf", onefile = FALSE, paper = "special", height = ht, width = wd, pointsize = 2)
114 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))
115 dev.off()
116
117 pdf("heatmap_no_borders.pdf", onefile = FALSE, paper = "special", height = ht, width = wd, pointsize = 2)
118 pheatmap_j(heat_df, scale="none", border_color=NA, cluster_rows=FALSE, cluster_cols=FALSE, col=colorRampPalette(c("#FFFFFF", brewer.pal(9,"Blues")))(100))
119 dev.off()
120
121 #plot bait vs bait heatmap using dist matrix
122 dist_bait <- dist_bait/max(dist_bait)
123 pdf("bait2bait.pdf", onefile = FALSE, paper = "special")
124 heatmap_2j(as.matrix(dist_bait), trace="none", scale="none", density.info="none", col=rev(colorRampPalette(c("#FFFFFF", brewer.pal(9,"Blues")))(100)), xMin=0, xMax=1, margins=c(1.5*max(nchar(rownames(as.matrix(dist_bait)))),1.5*max(nchar(colnames(as.matrix(dist_bait))))))
125 dev.off()