Mercurial > repos > alermine > nebula

comparison [APliBio]Nebula tools suite/Nebula/MakeTSSdist/makeTSSdist.R @ 0:2ec3ba0e9e70 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author alermine
date Thu, 25 Oct 2012 08:18:25 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:2ec3ba0e9e70
1 #usage $0 STEP RIGHT chipPeaks outputFile.png output.txt [controlPeaks] [1 for pdf]
2 args <- commandArgs()
3 print (args)
4 myStep <- type.convert(args[4])
5 maxValue <- type.convert(args[5])
6
7 dataTable <-read.table(file=paste(args[6],".genes.ClosestPeakDist", sep=""), header=TRUE);
8 chip.genes.ClosestPeakDist<-data.frame(dataTable)
9 ifReg <- 0
10 if (length(unique(chip.genes.ClosestPeakDist$Reg))>1) {
11 ifReg <- 1
12 }
13 ifControl <- 0
14
15
16 ifPDF <- 0
17 if (length(args)>=10) {
18 ifPDF=args[10]
19 }
20 if (length(args)==9 & args[9]==1) {
21 ifPDF=1
22 }
23
24 library(Hmisc)
25
26 if (length(args)>=9 & args[9]!=1 & args[9]!=0) {
27 dataTable <-read.table(file=paste(args[9],".genes.ClosestPeakDist", sep=""), header=TRUE);
28 control.genes.ClosestPeakDist<-data.frame(dataTable)
29 ifControl <- 1
30 }
31 if (ifReg & ifControl) {
32 if (ifPDF==1) {
33 pdf(file = args[7], width = 19, height = 8, pointsize = 20, bg = "white")
34 } else {
35 png(filename = args[7], width = 1440 , height = 680, units = "px", pointsize = 20, bg = "white", res = NA)
36 plot(1:10)
37 }
38 op <- par(mfrow = c(2,3))
39 } else {
40 if (ifPDF==1) {
41 pdf(file = args[7], width = 10, height = 13, pointsize = 20, bg = "white")
42 } else {
43 png(filename = args[7], width = 680, height = 880, units = "px", pointsize = 20, bg = "white", res = NA)
44 plot(1:10)
45 }
46 # plot(1:10)
47 op <- par(mfrow = c(2,1))
48 }
49 myColor <- 1
50 myColor[1] <- colors()[131]
51 myColor[2] <- "darkolivegreen3"
52 myColor[3] <- "azure4"
53 myColor[4] <- "royalblue3"
54 myColor[5] <- colors()[17]
55
56 myColorControl <- 1
57
58 myColorControl[1] <- colors()[24]
59 myColorControl[2] <- colors()[278]
60 myColorControl[3] <- colors()[305]
61 myColorControl[4] <- colors()[219]
62 myColorControl[5] <- colors()[343]
63
64
65
66 #for cumulative:
67 dist_real_f <- chip.genes.ClosestPeakDist
68 if (ifControl) {
69 dist_control_f <- control.genes.ClosestPeakDist
70 }
71 step <- myStep
72 lim <- maxValue
73 x <- 0
74 count <- 1
75 countL <-1
76 n.types <- 1
77 myLevels <- 0
78 countTotalCont <- 0
79 countTotal <-0
80 countLCont <- 0
81 cumTotalCont <- 0
82 if (ifReg) {
83 n.types <- length(levels(chip.genes.ClosestPeakDist$Reg))
84 myLevels <- levels(chip.genes.ClosestPeakDist$Reg)
85 cum = matrix( 0, nrow=lim/step +1, ncol=n.types, byrow = TRUE)
86 for (i in c(1:n.types)) {
87 t <- which ((dist_real_f$Reg==myLevels[i]))
88 countL[i] <- length(t)
89 }
90 count <-1
91 for (i in seq(length=lim/step +1, from=0, by=step)) {
92 for (t in c(1:n.types)) {
93 tt <- which ((dist_real_f$Reg==myLevels[t])&(dist_real_f$Dist<=i)&(dist_real_f$Dist>=-i))
94 cum[count,t] <- length(tt)
95 }
96 x[count] <- i
97 count <- count + 1
98 }
99 ymax <- max(cum[,1]/countL[1], na.rm=TRUE)
100 for (i in c(2:n.types)) {
101 ymax <- max(ymax,max(cum[,i]/countL[i], na.rm=TRUE))
102 }
103 myLocCol <- myColor[2]
104
105 par(mar=c(5.1, 7.1, 4.1, 2.1))
106
107 plot (x,cum[,1]/countL[1] ,col = myColor[2],type="l", main="",xlab="",ylab="", lwd = 2, xlim = c(0, lim),xaxt="n" , ylim=c(0,ymax))
108 for (i in c(2:n.types)) {
109 colorr <- i+1
110 myLocCol <- c(myLocCol,myColor[colorr])
111 lines (x,cum[,i]/countL[i] ,col = myColor[colorr],type="l", lwd = 2)
112 print (myColor[colorr])
113 }
114
115 gradi <- 1000
116 if (lim>10000) {
117 gradi <- 10000
118 }
119 if (lim<3000) {
120 gradi <- 500
121 }
122 axisx <- seq(length=lim/gradi+1, from=0, by=gradi)
123 axisxlab <- paste(axisx/1000,"", sep = "")
124 axis(1, at=axisx,labels=axisxlab , las=1, cex.axis=1)
125 ymax <- max(cum[,i]/countL[i], na.rm=TRUE)
126
127 minor.tick(nx=5,tick.ratio=0.5)
128
129 legend(lim*0.45, ymax*0.45, myLevels, cex=1, lwd = 2, bty = "n", col = myLocCol, lty = c(1), pt.bg= c(myLocCol) , merge = TRUE)
130
131 title( main="",xlab="Distance from TSS (Kb)",ylab="Proportion of genes with a peak\nat a given distance (cumulative)")
132
133 if (ifControl) {
134 count <-1
135 n.types <- length(levels(control.genes.ClosestPeakDist$Reg))
136 myLevels <- levels(control.genes.ClosestPeakDist$Reg)
137 cumCont = matrix( 0, nrow=lim/step +1, ncol=n.types, byrow = TRUE)
138 for (i in c(1:n.types)) {
139 t <- which ((dist_control_f$Reg==myLevels[i]))
140 countLCont[i] <- length(t)
141 }
142 for (i in seq(length=lim/step +1, from=0, by=step)) {
143 for (t in c(1:n.types)) {
144 tt <- which ((dist_control_f$Reg==myLevels[t])&(dist_control_f$Dist<=i)&(dist_control_f$Dist>=-i))
145 cumCont[count,t] <- length(tt)
146 }
147 x[count] <- i
148 count <- count + 1
149 }
150 ymax <- max(cumCont[,1]/countLCont[1], na.rm=TRUE)
151 for (i in c(2:n.types)) {
152 ymax <- max(ymax,max(cumCont[,i]/countLCont[i], na.rm=TRUE))
153 }
154 myLocColCntr <- myColorControl[2]
155 plot (x,cumCont[,1]/countLCont[1] ,col = myLocColCntr[1],type="l", main="",xlab="",ylab="", lwd = 2, xlim = c(0, lim),xaxt="n" , ylim=c(0,ymax))
156 for (i in c(2:n.types)) {
157 colorr <- i+1
158 myLocColCntr <- c(myLocColCntr,myColorControl[colorr])
159 lines (x,cumCont[,i]/countLCont[i] ,col = myColorControl[colorr],type="l", lwd = 2)
160 }
161 if (lim>10000) {
162 gradi <- 10000
163 }
164 if (lim<3000) {
165 gradi <- 500
166 }
167 axisx <- seq(length=lim/gradi+1, from=0, by=gradi)
168 axisxlab <- paste(axisx/1000, sep = "")
169 axis(1, at=axisx,labels=axisxlab , las=1, cex.axis=1)
170 minor.tick(nx=5,tick.ratio=0.5)
171 legend(lim*0.45, ymax*0.45, myLevels, cex=1 , lwd = 2, bty = "n", col = myLocColCntr, lty = c(1), pt.bg= c(myLocCol) , merge = TRUE)
172 title( main="",xlab="Distance from TSS (Kb)",ylab="Proportion of genes with a peak\nat a given distance (cumulative)")
173 #real_vs_control_cumulative:
174 count <-1
175 countTotal <- length(dist_real_f$Reg)
176 cumTotal <- 0
177 for (i in seq(length=lim/step +1, from=0, by=step)) {
178 t <- which ((dist_real_f$Dist<=i)&(dist_real_f$Dist>=-i))
179 cumTotal[count] <- length(t)
180 x[count] <- i
181 count <- count + 1
182 }
183 plot (x,cumTotal/countTotal ,col = myColor[1],type="l", main="",xlab="",ylab="", lwd = 2, xlim = c(0, lim),xaxt="n" )
184 gradi <- 1000
185 if (lim>10000) {
186 gradi <- 10000
187 }
188 if (lim<3000) {
189 gradi <- 500
190 }
191 axisx <- seq(length=lim/gradi+1, from=0, by=gradi)
192 axisxlab <- paste(axisx/1000, sep = "")
193 axis(1, at=axisx,labels=axisxlab , las=1, cex.axis=1)
194 ymax <- max(cumTotal/countTotal, na.rm=TRUE)
195 minor.tick(nx=5,tick.ratio=0.5)
196 countTotalCont <- length(dist_control_f$Reg)
197 cumTotalCont <- 0
198 count <- 1
199 for (i in seq(length=lim/step +1, from=0, by=step)) {
200 t <- which ((dist_control_f$Dist<=i)&(dist_control_f$Dist>=-i))
201 cumTotalCont[count] <- length(t)
202 x[count] <- i
203 count <- count + 1
204 }
205 lines (x,cumTotalCont/countTotalCont ,col = colors()[328],type="l", lwd = 2)
206 legend(lim*0.45, ymax*0.45, c("ChIP","Control"), cex=1 , lwd = 2, bty = "n", col = c(myColor[1], colors()[328]), lty = c(1), pt.bg= c(myColor[1], colors()[328]) , merge = TRUE)
207 title( main="",xlab="Distance from TSS (Kb)",ylab="Proportion of genes with a peak\nat a given distance (cumulative)")
208 }
209 } else {
210 countTotal <- length(dist_real_f$Reg)
211 cumTotal <- 0
212 count <-1
213
214 gradi <- 1000
215 if (lim>10000) {
216 gradi <- 10000
217 }
218 if (lim<3000) {
219 gradi <- 500
220 }
221
222 for (i in seq(length=lim/step +1, from=0, by=step)) {
223 t <- which ((dist_real_f$Dist<=i)&(dist_real_f$Dist>=-i))
224 cumTotal[count] <- length(t)
225 x[count] <- i
226 count <- count + 1
227 }
228 par(mar=c(5.1, 7.1, 4.1, 2.1))
229
230 plot (x,cumTotal/countTotal ,col = myColor[1],type="l", main="",xlab="",ylab="", lwd = 2, xlim = c(0, lim),xaxt="n" )
231 axisx <- seq(length=lim/gradi+1, from=0, by=gradi)
232 axisxlab <- paste(axisx/1000, sep = "")
233 axis(1, at=axisx,labels=axisxlab , las=1, cex.axis=1)
234 title( main="",xlab="Distance from TSS (Kb)",ylab="Proportion of genes with a peak\nat a given distance (cumulative)")
235 ymax <- max(cumTotal/countTotal, na.rm=TRUE)
236 if (ifControl) {
237 countTotalCont <- length(dist_control_f$Reg)
238 cumTotalCont <- 0
239 count <- 1
240 for (i in seq(length=lim/step +1, from=0, by=step)) {
241 t <- which ((dist_control_f$Dist<=i)&(dist_control_f$Dist>=-i))
242 cumTotalCont[count] <- length(t)
243 x[count] <- i
244 count <- count + 1
245 }
246 lines (x,cumTotalCont/countTotalCont ,col = colors()[328],type="l", lwd = 2)
247 legend(lim*0.45, ymax*0.45, c("ChIP","Control"), cex=1 , lwd = 2, bty = "n", col = c(myColor[1], colors()[328]), lty = c(1), pt.bg= c(myColor[1], colors()[328]) , merge = TRUE)
248 } else {
249 legend(lim*0.45, ymax*0.45, c("ChIP"), cex=1 , lwd = 2, bty = "n", col = c(myColor[1]), lty = c(1), pt.bg= c(myColor[1]) , merge = TRUE)
250 }
251 }
252
253 sink(args[8], append=FALSE, split=FALSE)
254 if (ifReg) {
255 if (ifControl) {
256 cat (paste("Dist_TSS","% genes w/ a peak in ChIP","% genes w/ a peak in control",sep='\t'))
257 cat("\t")
258 for (i in c(1:n.types)) {
259 cat(paste("% ", myLevels[i]," genes w/ a peak in ChIP", sep=""))
260 cat("\t")
261 }
262
263 for (i in c(1:n.types)) {
264 cat(paste("% ", myLevels[i]," genes w/ a peak in Control", sep=""))
265 cat("\t")
266 }
267 cat("\n")
268 for (i in c(1:length(x))) {
269 cat(paste(x[i],cumTotal[i]/countTotal,cumTotalCont[i]/countTotalCont,sep="\t"))
270 cat("\t")
271 for (t in c(1:n.types)) {
272 cat(paste(cum[i,t]/countL[t],"\t", sep=""))
273 }
274 for (t in c(1:n.types)) {
275 cat(paste(cumCont[i,t]/countLCont[t],"\t", sep=""))
276 }
277 cat("\n")
278 }
279 }else {
280 cat (paste("Dist_TSS","\t",sep=''))
281 for (i in c(1:n.types)) {
282 cat(paste("% ", myLevels[i]," genes w/ a peak in ChIP", "\t", sep=""))
283 }
284 cat("\n")
285 for (i in c(1:length(x))) {
286 cat(paste(x[i],"\t",sep=""))
287 for (t in c(1:n.types)) {
288 cat(paste(cum[i,t]/countL[t],"\t", sep=""))
289 }
290 cat("\n")
291 }
292 }
293 } else {
294 if (ifControl) {
295 cat (paste("Dist_TSS","% genes w/ a peak in ChIP","% genes w/ a peak in control",sep='\t'))
296 cat("\n")
297 for (i in c(1:length(x))) {
298 cat(paste(x[i],cumTotal[i]/countTotal,cumTotalCont[i]/countTotalCont,sep="\t"))
299 cat("\n")
300 }
301 }else {
302 cat (paste("Dist_TSS","% genes w/ a peak in ChIP",sep='\t'))
303 cat("\n")
304 for (i in c(1:length(x))) {
305 cat(paste(x[i],cumTotal[i]/countTotal,sep="\t"))
306 cat("\n")
307 }
308
309 }
310 }
311
312
313 #stop sinking:
314 sink()
315
316 #around TSS:
317 lim <- maxValue
318 step <- myStep
319 my_breaks <- seq(length=lim/step*2 +1, from=-lim, by=step)
320 chip.genes <- read.table(file=paste(args[6],".genes", sep=""), header=TRUE) ;
321 dist_real_f <- chip.genes
322 if (ifControl) {
323 control.genes <- read.table(file=paste(args[9],".genes", sep=""), header=TRUE) ;
324 dist_control_f<-data.frame(control.genes)
325 }
326 if (ifReg) {
327 #n.types <- length(levels(chip.genes.ClosestPeakDist$Reg))
328 #myLevels <- levels(dist_real_f$Reg)
329
330 t<- which (dist_real_f$Reg==myLevels[1])
331 values_real <-dist_real_f$Dist[t]
332 hTSSreal = hist(values_real,plot=FALSE,breaks = c(min(values_real),my_breaks,max(values_real)) )
333 ymax <- max(hTSSreal$density, na.rm=TRUE)
334 for (i in c(2:n.types)) {
335 t<- which (dist_real_f$Reg==myLevels[i])
336 values_real <-dist_real_f$Dist[t]
337 hTSSreal = hist(values_real,plot=FALSE,breaks = c(min(values_real),my_breaks,max(values_real)) )
338 ymax <- max(ymax,max(hTSSreal$density, na.rm=TRUE))
339 }
340
341
342 t<- which (dist_real_f$Reg==myLevels[1])
343 values_real <-dist_real_f$Dist[t]
344 hTSSreal = hist(values_real,plot=FALSE,breaks = c(min(values_real),my_breaks,max(values_real)) )
345 plot (hTSSreal$mids,hTSSreal$density,col = myLocCol[1],type="l", main="",xlab="",ylab="", lwd = 2, xlim = c(-lim, lim),ylim = c(0, ymax), xaxt="n" )
346
347 title( main="",xlab="Distance from TSS (Kb)",ylab="Proportion of genes with a peak\nat a given distance (density)")
348
349 for (i in c(2:n.types)) {
350 t<- which (dist_real_f$Reg==myLevels[i])
351 values_real <-dist_real_f$Dist[t]
352 hTSSreal = hist(values_real,plot=FALSE,breaks = c(min(values_real),my_breaks,max(values_real)) )
353 lines (hTSSreal$mids,hTSSreal$density,col = myLocCol[i],type="l", lwd = 2)
354 }
355 legend(lim*0.1, ymax*0.9, myLevels, cex=1 , lwd = 2, bty = "n", col = myLocCol, lty = c(1), pt.bg= c(myLocCol) , merge = TRUE)
356
357 gradi <- 1000
358 if (lim>10000) {
359 gradi <- 10000
360 }
361 if (lim<3000) {
362 gradi <- 500
363 }
364
365 axisx <- seq(length=2*lim/gradi+1, from=-lim, by=gradi)
366 axisxlab <- paste(axisx/1000, sep = "")
367 axis(1, at=axisx,labels=axisxlab , las=1, cex.axis=1)
368
369
370 #minor.tick(nx=10,tick.ratio=0.5)
371 if (ifControl) {
372 t<- which (dist_control_f$Reg==myLevels[1])
373 values_control <-dist_control_f$Dist[t]
374 hTSScontrol= hist(values_control,plot=FALSE,breaks = c(min(values_control),my_breaks,max(values_control)) )
375 ymax <- max(hTSScontrol$density, na.rm=TRUE)
376 for (i in c(2:n.types)) {
377 t<- which (dist_control_f$Reg==myLevels[i])
378 values_control <-dist_control_f$Dist[t]
379 hTSScontrol = hist(values_control,plot=FALSE,breaks = c(min(values_control),my_breaks,max(values_control)) )
380 ymax <- max(ymax,max(hTSScontrol$density, na.rm=TRUE))
381 }
382 t<- which (dist_control_f$Reg==myLevels[1])
383 values_control <-dist_control_f$Dist[t]
384 hTSScontrol= hist(values_control,plot=FALSE,breaks = c(min(values_control),my_breaks,max(values_control)) )
385 plot (hTSScontrol$mids,hTSScontrol$density,col = myLocColCntr[1],type="l", main="",xlab="",ylab="", lwd = 2, xlim = c(-lim, lim),ylim = c(0, ymax),xaxt="n" )
386 title( main="",xlab="Distance from TSS (Kb)",ylab="Proportion of genes with a peak\nat a given distance (density)")
387 for (i in c(2:n.types)) {
388 t<- which (dist_control_f$Reg==myLevels[i])
389 values_control <-dist_control_f$Dist[t]
390 hTSScontrol = hist(values_control,plot=FALSE,breaks = c(min(values_control),my_breaks,max(values_control)) )
391 lines (hTSScontrol$mids,hTSScontrol$density,col = myLocColCntr[i],type="l", lwd = 2)
392 }
393
394 gradi <- 1000
395 if (lim>10000) {
396 gradi <- 10000
397 }
398 if (lim<3000) {
399 gradi <- 500
400 }
401
402 axisx <- seq(length=2*lim/gradi+1, from=-lim, by=gradi)
403 axisxlab <- paste(axisx/1000, sep = "")
404 axis(1, at=axisx,labels=axisxlab , las=1, cex.axis=1)
405
406 legend(lim*0.1, ymax*0.9, myLevels, cex=1 , lwd = 2, bty = "n", col = myLocColCntr, lty = c(1), pt.bg= c(myLocCol) , merge = TRUE)
407
408 # minor.tick(nx=10,tick.ratio=0.5)
409 #control vs real
410 values_real <-dist_real_f$Dist
411 hTSSreal = hist(values_real, plot=FALSE, breaks = c(min(values_real),my_breaks,max(values_real)) )
412 plot (hTSSreal$mids,hTSSreal$density,col = myColor[1],type="l", main="",xlab="",ylab="", lwd = 2, xlim = c(-lim, lim),xaxt="n")
413 title( main="",xlab="Distance from TSS (Kb)",ylab="Proportion of genes with a peak\nat a given distance (density)")
414 ymax <- max(hTSSreal$density, na.rm=TRUE)
415 values_control <-dist_control_f$Dist
416 hTSScontrol = hist(values_control, plot=FALSE, breaks = c(min(values_control),my_breaks,max(values_control)) )
417 lines (hTSScontrol$mids,hTSScontrol$density,col = colors()[328],type="l", lwd = 2)
418 legend(lim*0.2, ymax*0.9, c("ChIP","Control"), cex=1 , lwd = 2, bty = "n", col = c(myColor[1], colors()[328]), lty = c(1), pt.bg= c(myColor[1], colors()[328]) , merge = TRUE)
419
420 gradi <- 1000
421 if (lim>10000) {
422 gradi <- 10000
423 }
424 if (lim<3000) {
425 gradi <- 500
426 }
427
428 axisx <- seq(length=2*lim/gradi+1, from=-lim, by=gradi)
429 axisxlab <- paste(axisx/1000, sep = "")
430 axis(1, at=axisx,labels=axisxlab , las=1, cex.axis=1)
431
432
433
434 # minor.tick(nx=10,tick.ratio=0.5)
435 }
436 } else {
437 values_real <-dist_real_f$Dist
438 hTSSreal = hist(values_real, plot=FALSE, breaks = c(min(values_real),my_breaks,max(values_real)) )
439 plot (hTSSreal$mids,hTSSreal$density,col = myColor[1],type="l", main="",xlab="",ylab="", lwd = 2, xlim = c(-lim, lim),xaxt="n")
440 title( main="",xlab="Distance from TSS (Kb)",ylab="Proportion of genes with a peak\nat a given distance (density)")
441 ymax <- max(hTSSreal$density, na.rm=TRUE)
442 if (ifControl) {
443 values_control <-dist_control_f$Dist
444 hTSScontrol = hist(values_control, plot=FALSE, breaks = c(min(values_control),my_breaks,max(values_control)) )
445 lines (hTSScontrol$mids,hTSScontrol$density,col = colors()[328],type="l", lwd = 2)
446 legend(lim*0.2, ymax*0.9, c("ChIP","Control"), cex=1 , lwd = 2, bty = "n", col = c(myColor[1], colors()[328]), lty = c(1), pt.bg= c(myColor[1], colors()[328]) , merge = TRUE)
447 } else {
448 legend(lim*0.2, ymax*0.9, c("ChIP"), cex=1 , lwd = 2, bty = "n", col = c(myColor[1]), lty = c(1), pt.bg= c(myColor[1]) , merge = TRUE)
449 }
450
451 gradi <- 1000
452 if (lim>10000) {
453 gradi <- 10000
454 }
455 if (lim<3000) {
456 gradi <- 500
457 }
458
459 axisx <- seq(length=2*lim/gradi+1, from=-lim, by=gradi)
460 axisxlab <- paste(axisx/1000, sep = "")
461 axis(1, at=axisx,labels=axisxlab , las=1, cex.axis=1)
462
463
464 # minor.tick(nx=10,tick.ratio=0.5)
465 }
466 dev.off()
467 q();
468 cat (paste("peak height","# peaks in ChIP","# peaks in Control","#control/chip","\n",sep='\t'))
469 for (xval in c(minValue:maxValue)) {
470 for (i in (1:length(chipHist$mids))) {
471 if (xval==chipHist$mids[i]) {
472 ychip <- chipHist$counts[i]
473 }
474 }
475 for (i in (1:length(controlHist$mids))) {
476 if (xval==controlHist$mids[i]) {
477 ycontrol <- controlHist$counts[i]
478 }
479 }
480 cat (paste(xval,ychip,ycontrol,ycontrol/ychip,"\n",sep='\t'))
481 }