view [APliBio]Nebula tools suite/Nebula/MakeStatsCHiPSeq/makeStatsChIPSeq.R @ 3:1c699789d6d3 draft

Uploaded

#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'))
}