|
2
|
1 #usage $0 minValue maxValue chip.peaks control.peaks outputFile.png output.txt
|
|
|
2 args <- commandArgs(TRUE)
|
|
|
3 minValue <- type.convert(args[2])
|
|
|
4 maxValue <- type.convert(args[3])
|
|
|
5 dataTable <-suppressWarnings(read.table(args[4], header=FALSE));
|
|
|
6 chip<-data.frame(dataTable)
|
|
|
7 dataTable <-suppressWarnings(read.table(args[5], header=FALSE));
|
|
|
8 control<-data.frame(dataTable)
|
|
|
9 x <-c((minValue-0.5):(maxValue+0.5))
|
|
|
10 breaks <- c(0,x,1000)
|
|
|
11 controlHist <-hist(control$V6, breaks=breaks, right=FALSE, plot=FALSE )
|
|
|
12 chipHist <-hist(chip$V6, breaks=breaks, right=FALSE, plot=FALSE )
|
|
|
13
|
|
|
14 ifPDF <- 0
|
|
|
15 if (length(args)>=8) {
|
|
|
16 ifPDF=args[8]
|
|
|
17 }
|
|
|
18 if (ifPDF==1) {
|
|
|
19 pdf(file = args[6], width = 8, height = 7, pointsize = 20, bg = "white")
|
|
|
20 } else {
|
|
|
21 png(filename = args[6], width = 580, height = 580, units = "px", pointsize = 20, bg = "white", res = NA)
|
|
|
22 }
|
|
|
23 suppressWarnings(plot(controlHist$mids,controlHist$counts,xlab = "Peak height",xlim=c(minValue-0.5,maxValue+0.5), ylab = "Peak count",pch=17, col = colors()[328], log = "y"))
|
|
|
24 suppressWarnings(points(chipHist$mids,chipHist$counts,xlab = "peak height",xlim=c(minValue-0.5,maxValue+0.5),ylab = "peak count",pch=15, col = colors()[131], log = "y"))
|
|
|
25 legend(maxValue*0.7,y = max(chipHist$counts)*0.7, bty="n",c("ChIP","Control"), cex=1, col = c(colors()[131],colors()[328]), lty = c(-1, -1), pch = c(15, 17))
|
|
|
26
|
|
|
27 dev.off()
|
|
|
28
|
|
|
29 sink(args[7], append=FALSE, split=FALSE)
|
|
|
30 cat (paste("peak height","# peaks in ChIP","# peaks in Control","#control/chip","Cumulative #control/chip (FDR)","\n",sep='\t'))
|
|
|
31 for (xval in c(minValue:maxValue)) {
|
|
|
32 for (i in (1:length(chipHist$mids))) {
|
|
|
33 if (xval==chipHist$mids[i]) {
|
|
|
34 ychip <- chipHist$counts[i]
|
|
|
35 }
|
|
|
36 }
|
|
|
37 for (i in (1:length(controlHist$mids))) {
|
|
|
38 if (xval==controlHist$mids[i]) {
|
|
|
39 ycontrol <- controlHist$counts[i]
|
|
|
40 }
|
|
|
41 }
|
|
|
42 ychipCum = 1;
|
|
|
43 ycontrolCum = 1 ;
|
|
|
44 for (i in (1:length(chipHist$mids))) {
|
|
|
45 if (xval<=chipHist$mids[i]) {
|
|
|
46 ychipCum <- ychipCum+chipHist$counts[i]
|
|
|
47 }
|
|
|
48 }
|
|
|
49 for (i in (1:length(controlHist$mids))) {
|
|
|
50 if (xval<=controlHist$mids[i]) {
|
|
|
51 ycontrolCum <- ycontrolCum+controlHist$counts[i]
|
|
|
52 }
|
|
|
53 }
|
|
|
54 cumRatio= ycontrolCum/ychipCum
|
|
|
55 if(cumRatio>1){cumRatio = 1}
|
|
|
56 cat (paste(xval,ychip,ycontrol,ycontrol/ychip,cumRatio,"\n",sep='\t'))
|
|
|
57 }
|
