|
0
|
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"
|
|
7
|
13 empty.region.filter = args[6]
|
|
|
14
|
|
0
|
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
|
|
7
|
24 if(empty.region.filter == "leader"){
|
|
31
|
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, before.unique$CDR3.IMGT.seq)
|
|
7
|
26 } else if(empty.region.filter == "FR1"){
|
|
31
|
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, before.unique$CDR3.IMGT.seq)
|
|
7
|
28 } else if(empty.region.filter == "CDR1"){
|
|
31
|
29 before.unique$seq_conc = paste(before.unique$FR2.IMGT.seq, before.unique$CDR2.IMGT.seq, before.unique$FR3.IMGT.seq, before.unique$CDR3.IMGT.seq)
|
|
7
|
30 } else if(empty.region.filter == "FR2"){
|
|
31
|
31 before.unique$seq_conc = paste(before.unique$CDR2.IMGT.seq, before.unique$FR3.IMGT.seq, before.unique$CDR3.IMGT.seq)
|
|
7
|
32 }
|
|
0
|
33
|
|
|
34 IDs = before.unique[,c("Sequence.ID", "seq_conc", "best_match", "Functionality")]
|
|
|
35 IDs$best_match = as.character(IDs$best_match)
|
|
|
36
|
|
7
|
37 dat = data.frame(table(before.unique$seq_conc))
|
|
0
|
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)
|
|
7
|
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
|
|
0
|
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
|
|
7
|
257 #NToverview = merged[!grepl("^unmatched", merged$best_match),]
|
|
|
258 NToverview = merged
|
|
0
|
259
|
|
7
|
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 }
|
|
0
|
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){
|
|
8
|
282 print(paste("class:", clazz))
|
|
0
|
283 NToverview.sub = NToverview[grepl(paste("^", clazz, sep=""), NToverview$best_match),]
|
|
8
|
284 print(paste("nrow:", nrow(NToverview.sub)))
|
|
0
|
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
|
