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