diff [APliBio]Nebula tools suite/Nebula/MakeStatsCHiPSeq/makeStatsChIPSeq.R @ 0:2ec3ba0e9e70 draft

Uploaded

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/[APliBio]Nebula tools suite/Nebula/MakeStatsCHiPSeq/makeStatsChIPSeq.R	Thu Oct 25 08:18:25 2012 -0400
@@ -0,0 +1,46 @@
+#usage $0 minValue maxValue chip.peaks control.peaks outputFile.png output.txt
+args <- commandArgs()
+minValue <- type.convert(args[4])
+maxValue <- type.convert(args[5])
+dataTable <-read.table(args[6], header=FALSE);
+chip<-data.frame(dataTable)
+dataTable <-read.table(args[7], header=FALSE);
+control<-data.frame(dataTable)
+x <-c((minValue-0.5):(maxValue+0.5))
+breaks <- c(0,x,1000)
+controlHist <-hist(control$V6, breaks=breaks, right=FALSE, plot=FALSE )
+chipHist <-hist(chip$V6, breaks=breaks, right=FALSE, plot=FALSE )
+
+ifPDF <- 0
+if (length(args)>=10) {
+	ifPDF=args[10]
+}
+if (ifPDF==1) {
+	pdf(file = args[8], width = 8, height = 7, pointsize = 20, bg = "white")
+} else {
+	png(filename = args[8], width = 580, height = 580, units = "px", pointsize = 20, bg = "white", res = NA)
+}
+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")
+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")
+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))
+
+#library(Hmisc)
+#minor.tick(nx=10, tick.ratio=0.5)
+
+dev.off()
+
+sink(args[9], append=FALSE, split=FALSE)
+cat (paste("peak height","# peaks in ChIP","# peaks in Control","#control/chip","\n",sep='\t'))
+for (xval in c(minValue:maxValue)) {
+  for (i in (1:length(chipHist$mids))) {     
+     if (xval==chipHist$mids[i]) {
+      ychip <- chipHist$counts[i]
+     }
+  }
+  for (i in (1:length(controlHist$mids))) {     
+     if (xval==controlHist$mids[i]) {
+      ycontrol <- controlHist$counts[i]
+     }
+  }
+  cat (paste(xval,ychip,ycontrol,ycontrol/ychip,"\n",sep='\t'))
+}