Mercurial > repos > davidvanzessen > argalaxy_tools
comparison sequence_overview.r @ 4:5ffd52fc35c4 draft
Uploaded
author | davidvanzessen |
---|---|
date | Mon, 12 Dec 2016 05:22:37 -0500 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
3:beaa487ecf43 | 4:5ffd52fc35c4 |
---|---|
1 library(reshape2) | |
2 | |
3 args <- commandArgs(trailingOnly = TRUE) | |
4 | |
5 before.unique.file = args[1] | |
6 merged.file = args[2] | |
7 outputdir = args[3] | |
8 gene.classes = unlist(strsplit(args[4], ",")) | |
9 hotspot.analysis.sum.file = args[5] | |
10 NToverview.file = paste(outputdir, "ntoverview.txt", sep="/") | |
11 NTsum.file = paste(outputdir, "ntsum.txt", sep="/") | |
12 main.html = "index.html" | |
13 empty.region.filter = args[6] | |
14 | |
15 | |
16 setwd(outputdir) | |
17 | |
18 before.unique = read.table(before.unique.file, header=T, sep="\t", fill=T, stringsAsFactors=F, quote="") | |
19 merged = read.table(merged.file, header=T, sep="\t", fill=T, stringsAsFactors=F, quote="") | |
20 hotspot.analysis.sum = read.table(hotspot.analysis.sum.file, header=F, sep=",", fill=T, stringsAsFactors=F, quote="") | |
21 | |
22 #before.unique = before.unique[!grepl("unmatched", before.unique$best_match),] | |
23 | |
24 if(empty.region.filter == "leader"){ | |
25 before.unique$seq_conc = paste(before.unique$FR1.IMGT.seq, before.unique$CDR1.IMGT.seq, before.unique$FR2.IMGT.seq, before.unique$CDR2.IMGT.seq, before.unique$FR3.IMGT.seq) | |
26 } else if(empty.region.filter == "FR1"){ | |
27 before.unique$seq_conc = paste(before.unique$CDR1.IMGT.seq, before.unique$FR2.IMGT.seq, before.unique$CDR2.IMGT.seq, before.unique$FR3.IMGT.seq) | |
28 } else if(empty.region.filter == "CDR1"){ | |
29 before.unique$seq_conc = paste(before.unique$FR2.IMGT.seq, before.unique$CDR2.IMGT.seq, before.unique$FR3.IMGT.seq) | |
30 } else if(empty.region.filter == "FR2"){ | |
31 before.unique$seq_conc = paste(before.unique$CDR2.IMGT.seq, before.unique$FR3.IMGT.seq) | |
32 } | |
33 | |
34 IDs = before.unique[,c("Sequence.ID", "seq_conc", "best_match", "Functionality")] | |
35 IDs$best_match = as.character(IDs$best_match) | |
36 | |
37 dat = data.frame(table(before.unique$seq_conc)) | |
38 | |
39 names(dat) = c("seq_conc", "Freq") | |
40 | |
41 dat$seq_conc = factor(dat$seq_conc) | |
42 | |
43 dat = dat[order(as.character(dat$seq_conc)),] | |
44 | |
45 #writing html from R... | |
46 get.bg.color = function(val){ | |
47 if(val %in% c("TRUE", "FALSE", "T", "F")){ #if its a logical value, give the background a green/red color | |
48 return(ifelse(val,"#eafaf1","#f9ebea")) | |
49 } else if (!is.na(as.numeric(val))) { #if its a numerical value, give it a grey tint if its >0 | |
50 return(ifelse(val > 0,"#eaecee","white")) | |
51 } else { | |
52 return("white") | |
53 } | |
54 } | |
55 td = function(val) { | |
56 return(paste("<td bgcolor='", get.bg.color(val), "'>", val, "</td>", sep="")) | |
57 } | |
58 tr = function(val) { | |
59 return(paste(c("<tr>", sapply(val, td), "</tr>"), collapse="")) | |
60 } | |
61 | |
62 make.link = function(id, clss, val) { | |
63 paste("<a href='", clss, "_", id, ".html'>", val, "</a>", sep="") | |
64 } | |
65 tbl = function(df) { | |
66 res = "<table border='1'>" | |
67 for(i in 1:nrow(df)){ | |
68 res = paste(res, tr(df[i,]), sep="") | |
69 } | |
70 res = paste(res, "</table>") | |
71 } | |
72 | |
73 cat("<table border='1' class='pure-table pure-table-striped'>", file=main.html, append=F) | |
74 | |
75 if(empty.region.filter == "leader"){ | |
76 cat("<caption>FR1+CDR1+FR2+CDR2+FR3+CDR3 sequences that show up more than once</caption>", file=main.html, append=T) | |
77 } else if(empty.region.filter == "FR1"){ | |
78 cat("<caption>CDR1+FR2+CDR2+FR3+CDR3 sequences that show up more than once</caption>", file=main.html, append=T) | |
79 } else if(empty.region.filter == "CDR1"){ | |
80 cat("<caption>FR2+CDR2+FR3+CDR3 sequences that show up more than once</caption>", file=main.html, append=T) | |
81 } else if(empty.region.filter == "FR2"){ | |
82 cat("<caption>CDR2+FR3+CDR3 sequences that show up more than once</caption>", file=main.html, append=T) | |
83 } | |
84 | |
85 cat("<tr>", file=main.html, append=T) | |
86 cat("<th>Sequence</th><th>Functionality</th><th>ca1</th><th>ca2</th><th>cg1</th><th>cg2</th><th>cg3</th><th>cg4</th><th>cm</th><th>un</th>", file=main.html, append=T) | |
87 cat("<th>total CA</th><th>total CG</th><th>number of subclasses</th><th>present in both Ca and Cg</th><th>Ca1+Ca2</th>", file=main.html, append=T) | |
88 cat("<th>Cg1+Cg2</th><th>Cg1+Cg3</th><th>Cg1+Cg4</th><th>Cg2+Cg3</th><th>Cg2+Cg4</th><th>Cg3+Cg4</th>", file=main.html, append=T) | |
89 cat("<th>Cg1+Cg2+Cg3</th><th>Cg2+Cg3+Cg4</th><th>Cg1+Cg2+Cg4</th><th>Cg1+Cg3+Cg4</th><th>Cg1+Cg2+Cg3+Cg4</th>", file=main.html, append=T) | |
90 cat("</tr>", file=main.html, append=T) | |
91 | |
92 | |
93 | |
94 single.sequences=0 #sequence only found once, skipped | |
95 in.multiple=0 #same sequence across multiple subclasses | |
96 multiple.in.one=0 #same sequence multiple times in one subclass | |
97 unmatched=0 #all of the sequences are unmatched | |
98 some.unmatched=0 #one or more sequences in a clone are unmatched | |
99 matched=0 #should be the same als matched sequences | |
100 | |
101 sequence.id.page="by_id.html" | |
102 | |
103 for(i in 1:nrow(dat)){ | |
104 | |
105 ca1 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^IGA1", IDs$best_match),] | |
106 ca2 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^IGA2", IDs$best_match),] | |
107 | |
108 cg1 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^IGG1", IDs$best_match),] | |
109 cg2 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^IGG2", IDs$best_match),] | |
110 cg3 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^IGG3", IDs$best_match),] | |
111 cg4 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^IGG4", IDs$best_match),] | |
112 | |
113 cm = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^IGM", IDs$best_match),] | |
114 | |
115 un = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^unmatched", IDs$best_match),] | |
116 allc = rbind(ca1, ca2, cg1, cg2, cg3, cg4, cm, un) | |
117 | |
118 ca1.n = nrow(ca1) | |
119 ca2.n = nrow(ca2) | |
120 | |
121 cg1.n = nrow(cg1) | |
122 cg2.n = nrow(cg2) | |
123 cg3.n = nrow(cg3) | |
124 cg4.n = nrow(cg4) | |
125 | |
126 cm.n = nrow(cm) | |
127 | |
128 un.n = nrow(un) | |
129 | |
130 classes = c(ca1.n, ca2.n, cg1.n, cg2.n, cg3.n, cg4.n, cm.n, un.n) | |
131 | |
132 classes.sum = sum(classes) | |
133 | |
134 if(classes.sum == 1){ | |
135 single.sequences = single.sequences + 1 | |
136 next | |
137 } | |
138 | |
139 if(un.n == classes.sum){ | |
140 unmatched = unmatched + 1 | |
141 next | |
142 } | |
143 | |
144 in.classes = sum(classes > 0) | |
145 | |
146 matched = matched + in.classes #count in how many subclasses the sequence occurs. | |
147 | |
148 if(any(classes == classes.sum)){ | |
149 multiple.in.one = multiple.in.one + 1 | |
150 } else if (un.n > 0) { | |
151 some.unmatched = some.unmatched + 1 | |
152 } else { | |
153 in.multiple = in.multiple + 1 | |
154 } | |
155 | |
156 id = as.numeric(dat[i,"seq_conc"]) | |
157 | |
158 functionality = paste(unique(allc[,"Functionality"]), collapse=",") | |
159 | |
160 by.id.row = c() | |
161 | |
162 if(ca1.n > 0){ | |
163 cat(tbl(ca1), file=paste("IGA1_", id, ".html", sep="")) | |
164 } | |
165 | |
166 if(ca2.n > 0){ | |
167 cat(tbl(ca2), file=paste("IGA2_", id, ".html", sep="")) | |
168 } | |
169 | |
170 if(cg1.n > 0){ | |
171 cat(tbl(cg1), file=paste("IGG1_", id, ".html", sep="")) | |
172 } | |
173 | |
174 if(cg2.n > 0){ | |
175 cat(tbl(cg2), file=paste("IGG2_", id, ".html", sep="")) | |
176 } | |
177 | |
178 if(cg3.n > 0){ | |
179 cat(tbl(cg3), file=paste("IGG3_", id, ".html", sep="")) | |
180 } | |
181 | |
182 if(cg4.n > 0){ | |
183 cat(tbl(cg4), file=paste("IGG4_", id, ".html", sep="")) | |
184 } | |
185 | |
186 if(cm.n > 0){ | |
187 cat(tbl(cm), file=paste("IGM_", id, ".html", sep="")) | |
188 } | |
189 | |
190 if(un.n > 0){ | |
191 cat(tbl(un), file=paste("un_", id, ".html", sep="")) | |
192 } | |
193 | |
194 ca1.html = make.link(id, "IGA1", ca1.n) | |
195 ca2.html = make.link(id, "IGA2", ca2.n) | |
196 | |
197 cg1.html = make.link(id, "IGG1", cg1.n) | |
198 cg2.html = make.link(id, "IGG2", cg2.n) | |
199 cg3.html = make.link(id, "IGG3", cg3.n) | |
200 cg4.html = make.link(id, "IGG4", cg4.n) | |
201 | |
202 cm.html = make.link(id, "IGM", cm.n) | |
203 | |
204 un.html = make.link(id, "un", un.n) | |
205 | |
206 #extra columns | |
207 ca.n = ca1.n + ca2.n | |
208 | |
209 cg.n = cg1.n + cg2.n + cg3.n + cg4.n | |
210 | |
211 #in.classes | |
212 | |
213 in.ca.cg = (ca.n > 0 & cg.n > 0) | |
214 | |
215 in.ca1.ca2 = (ca1.n > 0 & ca2.n > 0) | |
216 | |
217 in.cg1.cg2 = (cg1.n > 0 & cg2.n > 0) | |
218 in.cg1.cg3 = (cg1.n > 0 & cg3.n > 0) | |
219 in.cg1.cg4 = (cg1.n > 0 & cg4.n > 0) | |
220 in.cg2.cg3 = (cg2.n > 0 & cg3.n > 0) | |
221 in.cg2.cg4 = (cg2.n > 0 & cg4.n > 0) | |
222 in.cg3.cg4 = (cg3.n > 0 & cg4.n > 0) | |
223 | |
224 in.cg1.cg2.cg3 = (cg1.n > 0 & cg2.n > 0 & cg3.n > 0) | |
225 in.cg2.cg3.cg4 = (cg2.n > 0 & cg3.n > 0 & cg4.n > 0) | |
226 in.cg1.cg2.cg4 = (cg1.n > 0 & cg2.n > 0 & cg4.n > 0) | |
227 in.cg1.cg3.cg4 = (cg1.n > 0 & cg3.n > 0 & cg4.n > 0) | |
228 | |
229 in.cg.all = (cg1.n > 0 & cg2.n > 0 & cg3.n > 0 & cg4.n > 0) | |
230 | |
231 | |
232 | |
233 | |
234 #rw = c(as.character(dat[i,"seq_conc"]), functionality, ca1.html, ca2.html, cg1.html, cg2.html, cg3.html, cg4.html, cm.html, un.html) | |
235 rw = c(as.character(dat[i,"seq_conc"]), functionality, ca1.html, ca2.html, cg1.html, cg2.html, cg3.html, cg4.html, cm.html, un.html) | |
236 rw = c(rw, ca.n, cg.n, in.classes, in.ca.cg, in.ca1.ca2, in.cg1.cg2, in.cg1.cg3, in.cg1.cg4, in.cg2.cg3, in.cg2.cg4, in.cg3.cg4, in.cg1.cg2.cg3, in.cg2.cg3.cg4, in.cg1.cg2.cg4, in.cg1.cg3.cg4, in.cg.all) | |
237 | |
238 cat(tr(rw), file=main.html, append=T) | |
239 | |
240 | |
241 for(i in 1:nrow(allc)){ #generate html by id | |
242 html = make.link(id, allc[i,"best_match"], allc[i,"Sequence.ID"]) | |
243 cat(paste(html, "<br />"), file=sequence.id.page, append=T) | |
244 } | |
245 } | |
246 | |
247 cat("</table>", file=main.html, append=T) | |
248 | |
249 print(paste("Single sequences:", single.sequences)) | |
250 print(paste("Sequences in multiple subclasses:", in.multiple)) | |
251 print(paste("Multiple sequences in one subclass:", multiple.in.one)) | |
252 print(paste("Matched with unmatched:", some.unmatched)) | |
253 print(paste("Count that should match 'matched' sequences:", matched)) | |
254 | |
255 #ACGT overview | |
256 | |
257 #NToverview = merged[!grepl("^unmatched", merged$best_match),] | |
258 NToverview = merged | |
259 | |
260 if(empty.region.filter == "leader"){ | |
261 NToverview$seq = paste(NToverview$FR1.IMGT.seq, NToverview$CDR1.IMGT.seq, NToverview$FR2.IMGT.seq, NToverview$CDR2.IMGT.seq, NToverview$FR3.IMGT.seq) | |
262 } else if(empty.region.filter == "FR1"){ | |
263 NToverview$seq = paste(NToverview$CDR1.IMGT.seq, NToverview$FR2.IMGT.seq, NToverview$CDR2.IMGT.seq, NToverview$FR3.IMGT.seq) | |
264 } else if(empty.region.filter == "CDR1"){ | |
265 NToverview$seq = paste(NToverview$FR2.IMGT.seq, NToverview$CDR2.IMGT.seq, NToverview$FR3.IMGT.seq) | |
266 } else if(empty.region.filter == "FR2"){ | |
267 NToverview$seq = paste(NToverview$CDR2.IMGT.seq, NToverview$FR3.IMGT.seq) | |
268 } | |
269 | |
270 NToverview$A = nchar(gsub("[^Aa]", "", NToverview$seq)) | |
271 NToverview$C = nchar(gsub("[^Cc]", "", NToverview$seq)) | |
272 NToverview$G = nchar(gsub("[^Gg]", "", NToverview$seq)) | |
273 NToverview$T = nchar(gsub("[^Tt]", "", NToverview$seq)) | |
274 | |
275 #Nsum = data.frame(Sequence.ID="-", best_match="Sum", seq="-", A = sum(NToverview$A), C = sum(NToverview$C), G = sum(NToverview$G), T = sum(NToverview$T)) | |
276 | |
277 #NToverview = rbind(NToverview, NTsum) | |
278 | |
279 NTresult = data.frame(nt=c("A", "C", "T", "G")) | |
280 | |
281 for(clazz in gene.classes){ | |
282 print(paste("class:", clazz)) | |
283 NToverview.sub = NToverview[grepl(paste("^", clazz, sep=""), NToverview$best_match),] | |
284 print(paste("nrow:", nrow(NToverview.sub))) | |
285 new.col.x = c(sum(NToverview.sub$A), sum(NToverview.sub$C), sum(NToverview.sub$T), sum(NToverview.sub$G)) | |
286 new.col.y = sum(new.col.x) | |
287 new.col.z = round(new.col.x / new.col.y * 100, 2) | |
288 | |
289 tmp = names(NTresult) | |
290 NTresult = cbind(NTresult, data.frame(new.col.x, new.col.y, new.col.z)) | |
291 names(NTresult) = c(tmp, paste(clazz, c("x", "y", "z"), sep="")) | |
292 } | |
293 | |
294 write.table(NToverview[,c("Sequence.ID", "best_match", "seq", "A", "C", "G", "T")], NToverview.file, quote=F, sep="\t", row.names=F, col.names=T) | |
295 | |
296 NToverview = NToverview[!grepl("unmatched", NToverview$best_match),] | |
297 | |
298 new.col.x = c(sum(NToverview$A), sum(NToverview$C), sum(NToverview$T), sum(NToverview$G)) | |
299 new.col.y = sum(new.col.x) | |
300 new.col.z = round(new.col.x / new.col.y * 100, 2) | |
301 | |
302 tmp = names(NTresult) | |
303 NTresult = cbind(NTresult, data.frame(new.col.x, new.col.y, new.col.z)) | |
304 names(NTresult) = c(tmp, paste("all", c("x", "y", "z"), sep="")) | |
305 | |
306 names(hotspot.analysis.sum) = names(NTresult) | |
307 | |
308 hotspot.analysis.sum = rbind(hotspot.analysis.sum, NTresult) | |
309 | |
310 write.table(hotspot.analysis.sum, hotspot.analysis.sum.file, quote=F, sep=",", row.names=F, col.names=F, na="0") | |
311 | |
312 | |
313 | |
314 | |
315 | |
316 | |
317 | |
318 | |
319 | |
320 | |
321 | |
322 | |
323 | |
324 | |
325 | |
326 | |
327 | |
328 | |
329 | |
330 | |
331 | |
332 | |
333 | |
334 | |
335 | |
336 | |
337 | |
338 | |
339 | |
340 |