Mercurial > repos > jbrayet > annotatepeaks

comparison catDist.R @ 8:375c3d21e8fc draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author jbrayet
date Thu, 13 Aug 2015 11:17:15 -0400
parents
children
comparison
equal deleted inserted replaced
7:f7123cad20d8 8:375c3d21e8fc
1 #usage $0 STEP RIGHT chipPeaks outputFile.png output.txt [controlPeaks]
2 args <- commandArgs(TRUE)
3
4 input <- args[2]
5 pngFile <- args[3]
6 dataTable <-read.table(file=input, header=TRUE);
7 chip.data<-data.frame(dataTable)
8 ifReg <- 0
9 if (length(unique(chip.data$Reg))>1) {
10 ifReg <- 1
11 }
12
13 ifPDF <- 0
14 if (length(args)>=5) {
15 ifPDF=args[5]
16 }
17 if (length(args)==4 & args[4]==1) {
18 ifPDF=1
19 }
20
21 ifControl <- 0
22 if (length(args)>=4 & args[4]!=1 & args[4]!=0) {
23 dataTable <-read.table(file=args[4], header=TRUE);
24 control.data<-data.frame(dataTable)
25 ifControl <- 1
26 }
27 if (ifReg & ifControl) {
28
29 if (ifPDF==1) {
30 pdf(file = pngFile, width = 14, height = 13, pointsize = 20, bg = "white")
31 } else {
32 png(filename = pngFile, type="cairo" , width = 1440, height = 1040, units = "px", pointsize = 20, bg = "white", res = NA)
33 plot(1:10)
34 }
35 op <- par(mfrow = c(3,2))
36 } else {
37
38 if (ifPDF==1) {
39 pdf(file = pngFile, width = 22, height = 8, pointsize = 20, bg = "white")
40 } else {
41 png(filename = pngFile, type="cairo" , width = 1580, height = 530, units = "px", pointsize = 20, bg = "white", res = NA)
42 plot(1:10)
43 }
44 op <- par(mfrow = c(1,2))
45 }
46 myColor <- 1
47 myColor[1] <- colors()[131]
48 myColor[2] <- colors()[59]
49 myColor[3] <- colors()[76]
50 myColor[4] <- colors()[88]
51 myColor[5] <- colors()[17]
52 myColor[6] <- colors()[565]
53 myColorControl <- 1
54 myColorControl[1] <- colors()[24]
55 myColorControl[2] <- colors()[278]
56 myColorControl[3] <- colors()[305]
57 myColorControl[4] <- colors()[219]
58 myColorControl[5] <- colors()[343]
59 myColorControl[6] <- colors()[245]
60 myLevels <- 0
61 if (ifReg) {
62 if (ifControl) {
63 #control vs real:
64
65 countTotal <- length(chip.data$Reg)
66 totalChIP <- summary(chip.data$Type)/countTotal
67 tt <- which(chip.data$Type=="intragenic")
68 subset.chip <- chip.data[tt,]
69 countIntra <- length(subset.chip$Reg)
70 intraChip<- summary(subset.chip$TypeIntra)/countTotal
71 nlev <- length(levels(chip.data$Type))
72 countTotalCont <- length(control.data$Reg)
73 totalContr <- summary(control.data$Type)/countTotalCont
74 tt <- which(control.data$Type=="intragenic")
75 subset.control <- control.data[tt,]
76 countIntraCont <- length(subset.control$Reg)
77 intraControl<- summary(subset.control$TypeIntra)/countTotalCont
78 cum = matrix( 0, nrow=2, ncol=nlev, byrow = TRUE)
79 for (i in c(1:nlev)) {
80 cum[1,i] <- totalChIP[i]
81 cum[2,i] <- totalContr[i]
82 }
83
84 labels<-c("GeneDown.", "Enh.", "Imm.Down.", "Interg.", "Intrag.", "Prom.")
85 if (length(labels)==length(levels(chip.data$Type))) {
86 barplot(cum,xlab="",beside=TRUE, col=c(myColor[1],colors()[328]), names.arg=labels,ylab="Proportion of peaks")
87 } else {
88 barplot(cum,xlab="",beside=TRUE, col=c(myColor[1],colors()[328]), names.arg=c(levels(chip.data$Type)),ylab="Proportion of peaks")
89 }
90
91 position <- 'topleft'
92 inset <- c(0.1, 0)
93 legend(position, c("ChIP","Control"), bty="n",fill=c(myColor[1],colors()[328]), inset=inset)
94
95 nlev <- length(levels(subset.chip$TypeIntra))
96 cum = matrix( 0, nrow=2, ncol=nlev, byrow = TRUE)
97 for (i in c(1:nlev)) {
98 cum[1,i] <- intraChip[i]
99 cum[2,i] <- intraControl[i]
100 }
101 barplot(cum,xlab="",beside=TRUE, col=c(myColor[1],colors()[328]), names.arg=c(levels(subset.chip$TypeIntra)),ylab="Proportion of peaks")
102
103 position <- 'topleft'
104 inset <- c(0.1, 0)
105 legend(position, c("ChIP","Control"), bty="n",fill=c(myColor[1],colors()[328]), inset=inset)
106 }
107 n.types <- length(levels(chip.data$Reg))
108 myLevels <- levels(chip.data$Reg)
109 nlev <- length(levels(chip.data$Type))
110 cum = matrix( 0, nrow=length(myLevels), ncol=nlev, byrow = TRUE)
111 countTotal <- length(chip.data$Reg)
112 colReg <-NULL
113 for (r in c(1:length(myLevels))) {
114 tt <- which(chip.data$Reg==myLevels[r])
115 totalChIP <- summary(chip.data$Type[tt])/countTotal
116 for (i in c(1:nlev)) {
117 cum[r,i] <- totalChIP[i]
118 }
119 colReg[r]<-myColor[r+3]
120 }
121
122 labels<-c("GeneDown.", "Enh.", "Imm.Down.", "Interg.", "Intrag.", "Prom.")
123 if (length(labels)==length(levels(chip.data$Type))) {
124 #barplot(cum,xlab="",beside=TRUE, col=c(myColor[1],myColor[5]), names.arg=labels,ylab="Proportion of peaks")
125 barplot(cum,xlab="",beside=TRUE, col=c(colReg), names.arg=labels,ylab="Proportion of peaks")
126 } else {
127 barplot(cum,xlab="",beside=TRUE, col=c(colReg), names.arg=c(levels(chip.data$Type)),ylab="Proportion of peaks")
128 }
129
130 position <- 'topleft'
131 inset <- c(0.1, 0)
132 legend(position, c(myLevels), bty="n",fill=c(colReg), inset=inset)
133
134
135 nlev <- length(levels(chip.data$TypeIntra))
136 cum = matrix( 0, nrow=length(myLevels), ncol=nlev, byrow = TRUE)
137 for (r in c(1:length(myLevels))) {
138 tt <- which(chip.data$Reg==myLevels[r]&chip.data$Type=="intragenic")
139 totalChIP <- summary(chip.data$TypeIntra[tt])/countTotal
140 for (i in c(1:nlev)) {
141 cum[r,i] <- totalChIP[i]
142 }
143 }
144 barplot(cum,xlab="",beside=TRUE, col=c(colReg), names.arg=c(levels(chip.data$TypeIntra)),ylab="Proportion of peaks")
145 position <- 'topleft'
146 inset <- c(0.1, 0)
147 legend(position, c(myLevels), bty="n",fill=c(colReg), inset=inset)
148
149
150 if (ifControl) {
151 nlev <- length(levels(control.data$Type))
152 cum = matrix( 0, nrow=length(myLevels), ncol=nlev, byrow = TRUE)
153 countTotal <- length(control.data$Reg)
154 colReg <-NULL
155 for (r in c(1:length(myLevels))) {
156 tt <- which(control.data$Reg==myLevels[r])
157 totalcontrol <- summary(control.data$Type[tt])/countTotal
158 for (i in c(1:nlev)) {
159 cum[r,i] <- totalcontrol[i]
160 }
161 colReg[r]<-myColorControl[r+3]
162 }
163 labels<-c("GeneDown.", "Enh.", "Imm.Down.", "Interg.", "Intrag.", "Prom.")
164 if (length(labels)==length(levels(chip.data$Type))) {
165 barplot(cum,xlab="",beside=TRUE, col=c(colReg), names.arg=labels,ylab="Proportion of peaks")
166 } else {
167 barplot(cum,xlab="",beside=TRUE, col=c(colReg), names.arg=c(levels(control.data$Type)),ylab="Proportion of peaks")
168 }
169 position <- 'topleft'
170 inset <- c(0.1, 0)
171 legend(position, c(myLevels), bty="n",fill=c(colReg), inset=inset)
172
173 nlev <- length(levels(control.data$TypeIntra))
174 cum = matrix( 0, nrow=length(myLevels), ncol=nlev, byrow = TRUE)
175 for (r in c(1:length(myLevels))) {
176 tt <- which(control.data$Reg==myLevels[r]&control.data$Type=="intragenic")
177 totalcontrol <- summary(control.data$TypeIntra[tt])/countTotal
178 for (i in c(1:nlev)) {
179 cum[r,i] <- totalcontrol[i]
180 }
181 }
182 barplot(cum,xlab="",beside=TRUE, col=c(colReg), names.arg=c(levels(control.data$TypeIntra)),ylab="Proportion of peaks")
183 position <- 'topleft'
184 inset <- c(0.1, 0)
185 legend(position, c(myLevels), bty="n",fill=c(colReg), inset=inset)
186 }
187 } else {
188 countTotal <- length(chip.data$Reg)
189 totalChIP <- summary(chip.data$Type)/countTotal
190 tt <- which(chip.data$Type=="intragenic")
191 subset.chip <- chip.data[tt,]
192 countIntra <- length(subset.chip$Reg)
193 intraChip<- summary(subset.chip$TypeIntra)/countTotal
194 nlev <- length(levels(chip.data$Type))
195 if (ifControl) {
196 countTotalCont <- length(control.data$Reg)
197 totalContr <- summary(control.data$Type)/countTotalCont
198 tt <- which(control.data$Type=="intragenic")
199 subset.control <- control.data[tt,]
200 countIntraCont <- length(subset.control$Reg)
201 intraControl<- summary(subset.control$TypeIntra)/countTotalCont
202 cum = matrix( 0, nrow=2, ncol=nlev, byrow = TRUE)
203 for (i in c(1:nlev)) {
204 cum[1,i] <- totalChIP[i]
205 cum[2,i] <- totalContr[i]
206 }
207
208 labels<-c("GeneDown.", "Enh.", "Imm.Down.", "Interg.", "Intrag.", "Prom.")
209 if (length(labels)==length(levels(chip.data$Type))) {
210 #barplot(cum,xlab="",beside=TRUE, col=c(myColor[1],myColor[5]), names.arg=labels,ylab="Proportion of peaks")
211 barplot(cum,xlab="",beside=TRUE, col=c(myColor[1],colors()[328]), names.arg=labels,ylab="Proportion of peaks")
212 } else {
213 barplot(cum,xlab="",beside=TRUE, col=c(myColor[1],colors()[328]), names.arg=c(levels(chip.data$Type)),ylab="Proportion of peaks")
214 }
215
216 position <- 'topleft'
217 inset <- c(0.1, 0)
218 legend(position,c("ChIP","Control"), bty="n",fill=c(myColor[1],colors()[328]), inset=inset)
219
220 nlev <- length(levels(subset.chip$TypeIntra))
221 cum = matrix( 0, nrow=2, ncol=nlev, byrow = TRUE)
222 for (i in c(1:nlev)) {
223 cum[1,i] <- intraChip[i]
224 cum[2,i] <- intraControl[i]
225 }
226 barplot(cum,xlab="",beside=TRUE, col=c(myColor[1],colors()[328]), names.arg=c(levels(subset.chip$TypeIntra)),ylab="Proportion of peaks")
227 position <- 'topleft'
228 inset <- c(0.1, 0)
229 legend(position,c("ChIP","Control"), bty="n",fill=c(myColor[1],colors()[328]), inset=inset)
230
231 } else {
232 labels<-c("GeneDown.", "Enh.", "Imm.Down.", "Interg.", "Intrag.", "Prom.")
233 if (length(labels)==length(levels(chip.data$Type))) {
234 barplot(totalChIP,xlab="", col=myColor, names.arg=labels,ylab="Proportion of peaks")
235 } else {
236 barplot(totalChIP,xlab="", col=myColor,ylab="Proportion of peaks")
237 }
238 barplot(intraChip,xlab="", col=myColor,ylab="Proportion of peaks")
239 }
240 }
241 dev.off()