annotate [APliBio]Nebula tools suite/Nebula/AnnotateGenes/geneDist.R @ 0:2ec3ba0e9e70 draft

Uploaded
author alermine
date Thu, 25 Oct 2012 08:18:25 -0400
parents
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
1 args <- commandArgs()
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
2 input <- args[4]
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
3 pngFile <- args[5]
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
4 dataTable <-read.table(file=input, header=TRUE);
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
5 chip.data<-data.frame(dataTable)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
6 ifReg <- 0
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
7 if (length(unique(chip.data$Reg))>1) {
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
8 ifReg <- 1
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
9 }
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
10
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
11 ifPDF <- 0
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
12 bootstrap <- 1
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
13
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
14 if (length(args)>=8) {
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
15 ifPDF=args[8]
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
16 bootstrap=args[9]
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
17 }
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
18 if (length(args)==7 & args[7]==1) {
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
19 ifPDF=1
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
20 }
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
21
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
22 ifControl <- 0
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
23 if (length(args)>=7 & args[7]!=1 & args[7]!=0) {
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
24 dataTable <-read.table(file=args[7], header=TRUE);
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
25 control.data<-data.frame(dataTable)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
26 ifControl <- 1
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
27
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
28
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
29 }
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
30
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
31 error.bar <- function(x, y, upper, lower=upper, length=0.1,...){
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
32 if(length(x) != length(y) | length(y) !=length(lower) | length(lower) != length(upper))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
33 stop("vectors must be same length")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
34 arrows(x,y+upper, x, y-lower, angle=90, code=3, length=length, ...)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
35 }
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
36
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
37
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
38
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
39 logFile <- args[6]
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
40 sink(logFile, append=FALSE, split=FALSE)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
41
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
42 if (ifReg & ifControl) {
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
43
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
44 if (ifPDF==1) {
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
45 pdf(file = pngFile, width = 14, height = 13, pointsize = 20, bg = "white")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
46 } else {
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
47 png(filename = pngFile, width = 1140, height = 840, units = "px", pointsize = 20, bg = "white", res = NA)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
48 plot(1:10)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
49 }
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
50 op <- par(mfrow = c(3,1))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
51 } else {
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
52 if (ifReg | ifControl) {
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
53
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
54 if (ifPDF==1) {
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
55 pdf(file = pngFile, width = 20, height = 17, pointsize = 20, bg = "white")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
56 } else {
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
57 png(filename = pngFile, width = 1580, height = 1230, units = "px", pointsize = 20, bg = "white", res = NA)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
58 plot(1:10)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
59 }
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
60 op <- par(mfrow = c(2,1))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
61 } else {
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
62 if (ifPDF==1) {
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
63 pdf(file = pngFile, width = 22, height = 8, pointsize = 20, bg = "white")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
64 } else {
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
65 png(filename = pngFile, width = 1580, height = 530, units = "px", pointsize = 20, bg = "white", res = NA)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
66 plot(1:10)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
67 }
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
68 op <- par(mfrow = c(1,1))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
69 }
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
70 }
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
71 myColor <- 1
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
72 myColor[1] <- colors()[131]
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
73 myColor[2] <- colors()[59]
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
74 myColor[3] <- colors()[76]
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
75 myColor[4] <- colors()[88]
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
76 myColor[5] <- colors()[17]
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
77 myColor[6] <- colors()[565]
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
78 myColor[7] <- colors()[454]
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
79 myColor[8] <- colors()[401]
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
80 myColor[9] <- colors()[99]
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
81 myColorControl <- 1
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
82 myColorControl[1] <- colors()[24]
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
83 myColorControl[2] <- colors()[278]
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
84 myColorControl[3] <- colors()[305]
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
85 myColorControl[4] <- colors()[219]
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
86 myColorControl[5] <- colors()[343]
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
87 myColorControl[6] <- colors()[245]
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
88 myLevels <- 0
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
89 nn <- colnames(chip.data)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
90 cc <- c(1:41)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
91 colnames(chip.data) <- cc
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
92 countTotal <- length(unique(chip.data$"1"))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
93 tt <- which(chip.data$"16">0)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
94 countTotalEnh <- length(unique(chip.data$"1"[tt]))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
95 tt <- which(chip.data$"10">0)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
96 countTotalProm <- length(unique(chip.data$"1"[tt]))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
97 tt <- which(chip.data$"12">0)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
98 countTotalImDown <- length(unique(chip.data$"1"[tt]))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
99 tt <- which(chip.data$"18">0)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
100 countTotalIntra <- length(unique(chip.data$"1"[tt]))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
101 tt <- which(chip.data$"20">0)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
102 countTotal5kbD <- length(unique(chip.data$"1"[tt]))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
103 tt <- which(chip.data$"28">0)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
104 countTotal1IntronI <- length(unique(chip.data$"1"[tt]))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
105 tt <- which(chip.data$"32">0)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
106 countTotalExonsI <- length(unique(chip.data$"1"[tt]))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
107 tt <- which(chip.data$"36">0)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
108 countTotalIntronsI <- length(unique(chip.data$"1"[tt]))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
109 tt <- which(chip.data$"40">0)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
110 countTotalJonctionsI <- length(unique(chip.data$"1"[tt]))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
111
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
112 names <- c("Enh.","Prom.","Imm.Down.","Intrag.","GeneDown.","F.Intron","Exons","2,3,etc.Introns","E.I.Junctions")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
113 yChIPTotal <- c (countTotalEnh/countTotal,countTotalProm/countTotal, countTotalImDown/countTotal,countTotalIntra/countTotal,countTotal5kbD/countTotal,countTotal1IntronI/countTotal,countTotalExonsI/countTotal,countTotalIntronsI/countTotal,countTotalJonctionsI/countTotal)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
114 if(!ifReg&!ifControl) {
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
115 par(mar=c(5.1, 7.1, 4.1, 2.1))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
116 barplot(yChIPTotal,xlab="",beside=TRUE, col=c(myColor), names.arg=c(names),ylab="Proportion of genes with a peak")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
117 cat ("Proportion of genes with a peak in a given genomic region:\n")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
118 cat (paste(c(names),sep='\t'))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
119 cat("\n")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
120 cat (paste(c(yChIPTotal),sep='\t'))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
121 cat("\n\n")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
122 }
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
123 if (ifControl) {
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
124 if (bootstrap>1) {
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
125 yControlTotalMatrix <- NULL
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
126 }
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
127 for (fileNumber in 1:bootstrap) {
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
128
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
129 if (fileNumber>=2) {
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
130 dataTable <-read.table(file=paste(args[7],fileNumber,sep=""), header=TRUE);
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
131 control.data<-data.frame(dataTable)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
132 }
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
133
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
134 colnames(control.data) <- cc
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
135 countTotalCntr <- length(unique(control.data$"1"))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
136 tt <- which(control.data$"16">0)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
137 countTotalEnhCntr <- length(unique(control.data$"1"[tt]))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
138 tt <- which(control.data$"10">0)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
139 countTotalPromCntr <- length(unique(control.data$"1"[tt]))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
140 tt <- which(control.data$"12">0)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
141 countTotalImDownCntr <- length(unique(control.data$"1"[tt]))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
142 tt <- which(control.data$"18">0)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
143 countTotalIntraCntr <- length(unique(control.data$"1"[tt]))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
144 tt <- which(control.data$"20">0)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
145 countTotal5kbDCntr <- length(unique(control.data$"1"[tt]))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
146 tt <- which(control.data$"28">0)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
147 countTotal1IntronICntr <- length(unique(control.data$"1"[tt]))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
148 tt <- which(control.data$"32">0)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
149 countTotalExonsICntr <- length(unique(control.data$"1"[tt]))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
150 tt <- which(control.data$"36">0)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
151 countTotalIntronsICntr <- length(unique(control.data$"1"[tt]))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
152 tt <- which(control.data$"40">0)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
153 countTotalJonctionsICntr <- length(unique(control.data$"1"[tt]))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
154 yControlTotal <- c (countTotalEnhCntr/countTotalCntr,countTotalPromCntr/countTotalCntr, countTotalImDownCntr/countTotalCntr,countTotalIntraCntr/countTotalCntr,countTotal5kbDCntr/countTotalCntr,countTotal1IntronICntr/countTotalCntr,countTotalExonsICntr/countTotalCntr,countTotalIntronsICntr/countTotalCntr,countTotalJonctionsICntr/countTotalCntr)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
155 if (bootstrap>1) {
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
156 yControlTotalMatrix <- rbind(yControlTotalMatrix,yControlTotal)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
157 }
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
158 }
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
159 if (bootstrap>1) {
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
160 yControlTotal <- colMeans(yControlTotalMatrix)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
161 Nrows <- nrow(yControlTotalMatrix)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
162 colVars <- Nrows/(Nrows-1) * (colMeans(yControlTotalMatrix*yControlTotalMatrix)-colMeans(yControlTotalMatrix)^2)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
163 }
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
164 if (!ifReg) {
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
165 cum = matrix( 0, nrow=2, ncol=length(names), byrow = TRUE)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
166 for (i in c(1:length(names))) {
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
167 cum[1,i] <- yChIPTotal[i]
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
168 cum[2,i] <- yControlTotal[i]
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
169 }
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
170 if (bootstrap>1) {
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
171 wiskers <- matrix(c(colVars-colVars,sqrt(colVars)),2,length(names),byrow=TRUE)*1.96
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
172 }
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
173 par(mar=c(5.1, 7.1, 4.1, 2.1))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
174 barx <- barplot(cum,xlab="",beside=TRUE, col=c(myColor[6],myColor[5]), legend = c("ChIP","Control"), names.arg=c(names),ylab="Proportion of genes with a peak")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
175 cat ("Proportion of genes with a peak from the ChIP dataset in a given genomic region:\n")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
176 cat (paste(c(names),sep='\t'))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
177 cat("\n")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
178 cat (paste(c(yChIPTotal),sep='\t'))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
179 cat("\n\n")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
180 cat ("Proportion of genes with a peak from the Control dataset in a given genomic region:\n")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
181 cat (paste(c(names),sep='\t'))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
182 cat("\n")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
183 cat (paste(c(yControlTotal),sep='\t'))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
184 cat("\n\n")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
185 if (bootstrap>1) {
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
186 error.bar(barx,cum,wiskers)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
187 cat ("Standard deviation for the Control dataset in a given genomic region:\n")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
188 cat (paste(c(names),sep='\t'))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
189 cat("\n")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
190 cat (paste(c(sqrt(colVars)),sep='\t'))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
191 cat("\n\n")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
192 }
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
193 enrich <- NULL
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
194 for (i in c(1:length(names))) {
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
195 enrich[i] <- yChIPTotal[i]/yControlTotal[i];
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
196 }
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
197 barplot(enrich-1,xlab="",beside=TRUE, col=c(myColor), names.arg=c(names),ylab="Enrichment in comparison\nwith the control dataset", yaxt="n")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
198 minX <- min(enrich-1)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
199 maxX <- max(enrich-1)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
200 x = seq(length=11, from=round(minX*10)/10, by=round((maxX-minX)*10)/100)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
201 axis(2, at=x,labels=x+1, las=2)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
202
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
203 cat ("Enrichment of genomic regions, ChIP peaks vs Control Peaks:\n")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
204 cat (paste(c(names),sep='\t'))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
205 cat("\n")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
206 cat (paste(c(yChIPTotal/yControlTotal),sep='\t'))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
207 cat("\n\n")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
208 if (bootstrap>1) {
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
209 z <- (yChIPTotal-yControlTotal)/sqrt(colVars)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
210 pvalues <- 2*pnorm(-abs(z))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
211 cat ("Two-side P-values for each genomic category:\n")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
212 cat (paste(c(names),sep='\t'))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
213 cat("\n")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
214 cat (paste(c(yChIPTotal/yControlTotal),sep='\t'))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
215 cat("\n\n")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
216 }
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
217 }
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
218 }
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
219 if (ifReg) {
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
220 n.types <- length(levels(chip.data$"6"))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
221 myLevels <- levels(chip.data$"6")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
222 nlev <- length(names)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
223
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
224 if (ifControl) {
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
225 cum = matrix( 0, nrow=length(myLevels)+1, ncol=nlev, byrow = TRUE)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
226 cumEnrichTotal = matrix( 0, nrow=length(myLevels), ncol=nlev, byrow = TRUE)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
227 cumEnrichControl = matrix( 0, nrow=length(myLevels), ncol=nlev, byrow = TRUE)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
228 }else {
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
229 cum = matrix( 0, nrow=length(myLevels), ncol=nlev, byrow = TRUE)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
230 cumEnrichTotal = matrix( 0, nrow=length(myLevels), ncol=nlev, byrow = TRUE)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
231 }
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
232 colReg <-NULL
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
233 for (r in c(1:length(myLevels))) {
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
234 tt <- which(chip.data$"6"==myLevels[r])
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
235 subset.data <- (chip.data[tt,])
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
236 countTotalSubset <- length(unique(subset.data$"1"))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
237 tt <- which(subset.data$"16">0)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
238 countTotalEnhSubset <- length(unique(subset.data$"1"[tt]))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
239 tt <- which(subset.data$"10">0)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
240 countTotalPromSubset <- length(unique(subset.data$"1"[tt]))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
241 tt <- which(subset.data$"12">0)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
242 countTotalImDownSubset <- length(unique(subset.data$"1"[tt]))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
243 tt <- which(subset.data$"18">0)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
244 countTotalIntraSubset <- length(unique(subset.data$"1"[tt]))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
245 tt <- which(subset.data$"20">0)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
246 countTotal5kbDSubset <- length(unique(subset.data$"1"[tt]))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
247 tt <- which(subset.data$"28">0)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
248 countTotal1IntronISubset <- length(unique(subset.data$"1"[tt]))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
249 tt <- which(subset.data$"32">0)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
250 countTotalExonsISubset <- length(unique(subset.data$"1"[tt]))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
251 tt <- which(subset.data$"36">0)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
252 countTotalIntronsISubset <- length(unique(subset.data$"1"[tt]))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
253 tt <- which(subset.data$"40">0)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
254 countTotalJonctionsISubset <- length(unique(subset.data$"1"[tt]))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
255 ySubsetTotal <- c (countTotalEnhSubset/countTotalSubset,countTotalPromSubset/countTotalSubset, countTotalImDownSubset/countTotalSubset,countTotalIntraSubset/countTotalSubset,countTotal5kbDSubset/countTotalSubset,countTotal1IntronISubset/countTotalSubset,countTotalExonsISubset/countTotalSubset,countTotalIntronsISubset/countTotalSubset,countTotalJonctionsISubset/countTotalSubset)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
256 for (i in c(1:nlev)) {
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
257 cum[r,i] <- ySubsetTotal[i]
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
258 cumEnrichTotal[r,i] <- ySubsetTotal[i]/yChIPTotal[i]
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
259 if (ifControl) {
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
260 cumEnrichControl[r,i] <- ySubsetTotal[i]/yControlTotal[i]
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
261 }
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
262 }
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
263 colReg[r]<-myColor[r+3]
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
264 }
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
265 if (ifControl) {
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
266 for (i in c(1:nlev)) {
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
267 cum[4,i] <- yControlTotal[i]
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
268 }
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
269 }
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
270
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
271 cat ("Proportion of genes with a peak from the ChIP dataset in a given genomic region:\n")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
272 for (r in c(1:length(myLevels))) {
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
273 cat (paste(myLevels[r],":\n",sep=""))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
274 cat (paste(c(names),sep='\t'))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
275 cat("\n")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
276 cat (paste(c(cum[r,] ),sep='\t'))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
277 cat("\n")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
278 }
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
279 cat("\n")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
280 par(mar=c(5.1, 7.1, 4.1, 2.1))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
281 if (ifControl) {
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
282 barx <- barplot(cum,xlab="",beside=TRUE, col=c(colReg, colors()[334]), legend = c(myLevels,"Control"), names.arg=c(names),ylab="Proportion of genes with a peak")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
283 cat ("Proportion of genes with a peak from the Control dataset in a given genomic region:\n")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
284 cat (paste(c(names),sep='\t'))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
285 cat("\n")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
286 cat (paste(c(yControlTotal),sep='\t'))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
287 cat("\n\n")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
288 if (bootstrap>1) {
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
289 wiskers <- cum-cum
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
290 wiskers[nrow(wiskers),] <- sqrt(colVars)*1.96
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
291 error.bar(barx,cum,wiskers)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
292 cat ("Standard deviation for the Control dataset in a given genomic region:\n")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
293 cat (paste(c(names),sep='\t'))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
294 cat("\n")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
295 cat (paste(c(sqrt(colVars)),sep='\t'))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
296 cat("\n\n")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
297 }
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
298 } else {
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
299 barplot(cum,xlab="",beside=TRUE, col=c(colReg), legend = c(myLevels), names.arg=c(names),ylab="Proportion of genes with a peak")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
300 }
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
301 barplot(cumEnrichTotal-1,xlab="",beside=TRUE, col=c(colReg), names.arg=c(names),ylab="Enrichment in comparison\nwith the total gene set", yaxt="n")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
302 minX <- min(cumEnrichTotal-1)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
303 maxX <- max(cumEnrichTotal-1)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
304 x = seq(length=11, from=round(minX*10)/10, by=round((maxX-minX)*10)/100)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
305 axis(2, at=x,labels=x+1, las=2)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
306
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
307 cat ("Enrichment of genomic regions, Transcriptional categories vs All Genes:\n")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
308 for (r in c(1:length(myLevels))) {
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
309 cat (paste(myLevels[r],":\n",sep=""))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
310 cat (paste(c(names),sep='\t'))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
311 cat("\n")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
312 cat (paste(c(cumEnrichTotal[r,]),sep='\t'))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
313 cat("\n")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
314 }
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
315 cat("\n")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
316
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
317 if (ifControl) {
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
318 barplot(cumEnrichControl-1,xlab="",beside=TRUE, col=c(colReg), names.arg=c(names),ylab="Enrichment in comparison\nwith control", yaxt="n")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
319 minX <- min(cumEnrichControl-1)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
320 maxX <- max(cumEnrichControl-1)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
321 x = seq(length=11, from=round(minX*10)/10, by=round((maxX-minX)*10)/100)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
322 axis(2, at=x,labels=x+1, las=2)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
323 cat ("Enrichment of genomic regions, ChIP peaks vs Control Peaks:\n")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
324 for (r in c(1:length(myLevels))) {
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
325 cat (paste(myLevels[r],":\n",sep=""))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
326 cat (paste(c(names),sep='\t'))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
327 cat("\n")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
328 cat (paste(c(cumEnrichControl[r,]),sep='\t'))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
329 cat("\n")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
330 }
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
331 cat("\n")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
332 if (bootstrap>1) {
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
333 cat ("Two-side P-values for each genomic category:\n")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
334 for (r in c(1:length(myLevels))) {
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
335 z <- (cum[r,]-yControlTotal)/sqrt(colVars)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
336 pvalues <- 2*pnorm(-abs(z))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
337 cat (paste(myLevels[r],":\n",sep=""))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
338 cat (paste(c(names),sep='\t'))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
339 cat("\n")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
340 cat (paste(c(pvalues),sep='\t'))
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
341 cat("\n")
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
342 }
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
343 }
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
344 }
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
345 }
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
346 sink() #stop sinking :)
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
347 dev.off()
2ec3ba0e9e70 Uploaded
alermine
parents:
diff changeset
348