|
76
|
1 library(reshape2)
|
|
|
2
|
|
|
3 args <- commandArgs(trailingOnly = TRUE)
|
|
|
4
|
|
100
|
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]
|
|
81
|
10 NToverview.file = paste(outputdir, "ntoverview.txt", sep="/")
|
|
|
11 NTsum.file = paste(outputdir, "ntsum.txt", sep="/")
|
|
76
|
12 main.html = "index.html"
|
|
|
13
|
|
|
14 setwd(outputdir)
|
|
|
15
|
|
100
|
16 before.unique = read.table(before.unique.file, header=T, sep="\t", fill=T, stringsAsFactors=F, quote="")
|
|
|
17 merged = read.table(merged.file, header=T, sep="\t", fill=T, stringsAsFactors=F, quote="")
|
|
82
|
18 hotspot.analysis.sum = read.table(hotspot.analysis.sum.file, header=F, sep=",", fill=T, stringsAsFactors=F, quote="")
|
|
76
|
19
|
|
102
|
20 #before.unique = before.unique[!grepl("unmatched", before.unique$best_match),]
|
|
82
|
21
|
|
100
|
22 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)
|
|
|
23
|
|
|
24 IDs = before.unique[,c("Sequence.ID", "seq_conc", "best_match", "Functionality")]
|
|
76
|
25 IDs$best_match = as.character(IDs$best_match)
|
|
|
26
|
|
|
27 #dat = data.frame(data.table(dat)[, list(freq=.N), by=c("best_match", "seq_conc")])
|
|
|
28
|
|
100
|
29 dat = data.frame(table(before.unique$seq_conc))
|
|
91
|
30 #dat = data.frame(table(merged$seq_conc, merged$Functionality))
|
|
76
|
31
|
|
87
|
32 #dat = dat[dat$Freq > 1,]
|
|
76
|
33
|
|
91
|
34 #names(dat) = c("seq_conc", "Functionality", "Freq")
|
|
|
35 names(dat) = c("seq_conc", "Freq")
|
|
76
|
36
|
|
|
37 dat$seq_conc = factor(dat$seq_conc)
|
|
|
38
|
|
82
|
39 dat = dat[order(as.character(dat$seq_conc)),]
|
|
76
|
40
|
|
|
41 #writing html from R...
|
|
|
42 td = function(val) { paste("<td>", val, "</td>", sep="") }
|
|
|
43 tr = function(val) { capture.output(cat("<tr>", td(val), "</tr>", sep="")) }
|
|
|
44 make.link = function(id, clss, val) { paste("<a href='", clss, "_", id, ".html'>", val, "</a>", sep="") }
|
|
|
45 tbl = function(df) { res = "<table border='1'>"; for(i in 1:nrow(df)){ res = paste(res, tr(df[i,]), sep=""); }; res = paste(res, "</table>"); }
|
|
|
46
|
|
|
47 cat("<table border='1'>", file=main.html, append=F)
|
|
88
|
48 cat("<caption>CDR1+FR2+CDR2+FR3+CDR3 sequences that show up more than once</caption>", file=main.html, append=T)
|
|
102
|
49 cat("<tr><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></tr>", file=main.html, append=T)
|
|
76
|
50
|
|
99
|
51
|
|
|
52
|
|
|
53 single.sequences=0 #sequence only found once, skipped
|
|
|
54 in.multiple=0 #same sequence across multiple subclasses
|
|
|
55 multiple.in.one=0 #same sequence multiple times in one subclass
|
|
102
|
56 unmatched=0 #all of the sequences are unmatched
|
|
|
57 some.unmatched=0 #one or more sequences in a clone are unmatched
|
|
99
|
58 matched=0 #should be the same als matched sequences
|
|
|
59
|
|
102
|
60 sequence.id.page="by_id.html"
|
|
|
61
|
|
76
|
62 for(i in 1:nrow(dat)){
|
|
98
|
63
|
|
91
|
64 ca1 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^ca1", IDs$best_match),]
|
|
|
65 ca2 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^ca2", IDs$best_match),]
|
|
76
|
66
|
|
91
|
67 cg1 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^cg1", IDs$best_match),]
|
|
|
68 cg2 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^cg2", IDs$best_match),]
|
|
|
69 cg3 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^cg3", IDs$best_match),]
|
|
|
70 cg4 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^cg4", IDs$best_match),]
|
|
76
|
71
|
|
91
|
72 cm = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^cm", IDs$best_match),]
|
|
76
|
73
|
|
102
|
74 un = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^unmatched", IDs$best_match),]
|
|
|
75 allc = rbind(ca1, ca2, cg1, cg2, cg3, cg4, cm, un)
|
|
|
76
|
|
|
77 classes = c(nrow(ca1), nrow(ca2), nrow(cg1), nrow(cg2), nrow(cg3), nrow(cg4), nrow(cm), nrow(un))
|
|
82
|
78
|
|
|
79 classes.sum = sum(classes)
|
|
|
80
|
|
87
|
81 if(classes.sum == 1){
|
|
99
|
82 single.sequences = single.sequences + 1
|
|
82
|
83 next
|
|
99
|
84 }
|
|
|
85
|
|
102
|
86 if(nrow(un) == classes.sum){
|
|
|
87 unmatched = unmatched + 1
|
|
|
88 next
|
|
|
89 }
|
|
|
90
|
|
99
|
91 matched = matched + sum(classes > 0) #count in how many subclasses the sequence occurs.
|
|
|
92
|
|
|
93 if(any(classes == classes.sum)){
|
|
102
|
94 multiple.in.one = multiple.in.one + 1
|
|
|
95 } else if (nrow(un) > 0) {
|
|
|
96 some.unmatched = some.unmatched + 1
|
|
98
|
97 } else {
|
|
102
|
98 in.multiple = in.multiple + 1
|
|
82
|
99 }
|
|
|
100
|
|
76
|
101 id = as.numeric(dat[i,"seq_conc"])
|
|
79
|
102
|
|
92
|
103 functionality = paste(unique(allc[,"Functionality"]), collapse=",")
|
|
102
|
104
|
|
|
105 by.id.row = c()
|
|
|
106
|
|
76
|
107 if(nrow(ca1) > 0){
|
|
|
108 cat(tbl(ca1), file=paste("ca1_", id, ".html", sep=""))
|
|
|
109 }
|
|
|
110
|
|
|
111 if(nrow(ca2) > 0){
|
|
|
112 cat(tbl(ca2), file=paste("ca2_", id, ".html", sep=""))
|
|
|
113 }
|
|
|
114
|
|
|
115 if(nrow(cg1) > 0){
|
|
|
116 cat(tbl(cg1), file=paste("cg1_", id, ".html", sep=""))
|
|
|
117 }
|
|
|
118
|
|
|
119 if(nrow(cg2) > 0){
|
|
|
120 cat(tbl(cg2), file=paste("cg2_", id, ".html", sep=""))
|
|
|
121 }
|
|
|
122
|
|
|
123 if(nrow(cg3) > 0){
|
|
|
124 cat(tbl(cg3), file=paste("cg3_", id, ".html", sep=""))
|
|
|
125 }
|
|
|
126
|
|
|
127 if(nrow(cg4) > 0){
|
|
|
128 cat(tbl(cg4), file=paste("cg4_", id, ".html", sep=""))
|
|
|
129 }
|
|
|
130
|
|
|
131 if(nrow(cm) > 0){
|
|
|
132 cat(tbl(cm), file=paste("cm_", id, ".html", sep=""))
|
|
|
133 }
|
|
102
|
134
|
|
|
135 if(nrow(un) > 0){
|
|
|
136 cat(tbl(un), file=paste("un_", id, ".html", sep=""))
|
|
|
137 }
|
|
76
|
138
|
|
|
139 ca1.html = make.link(id, "ca1", nrow(ca1))
|
|
|
140 ca2.html = make.link(id, "ca2", nrow(ca2))
|
|
|
141
|
|
|
142 cg1.html = make.link(id, "cg1", nrow(cg1))
|
|
|
143 cg2.html = make.link(id, "cg2", nrow(cg2))
|
|
|
144 cg3.html = make.link(id, "cg3", nrow(cg3))
|
|
|
145 cg4.html = make.link(id, "cg4", nrow(cg4))
|
|
|
146
|
|
|
147 cm.html = make.link(id, "cm", nrow(cm))
|
|
|
148
|
|
102
|
149 un.html = make.link(id, "un", nrow(un))
|
|
|
150
|
|
|
151 #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)
|
|
|
152 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)
|
|
76
|
153
|
|
|
154 cat(tr(rw), file=main.html, append=T)
|
|
102
|
155
|
|
|
156
|
|
|
157 for(i in 1:nrow(allc)){ #generate html by id
|
|
|
158 html = make.link(id, allc[i,"best_match"], allc[i,"Sequence.ID"])
|
|
|
159 cat(paste(html, "<br />"), file=sequence.id.page, append=T)
|
|
|
160 }
|
|
76
|
161 }
|
|
|
162
|
|
|
163 cat("</table>", file=main.html, append=T)
|
|
81
|
164
|
|
99
|
165 print(paste("Single sequences:", single.sequences))
|
|
|
166 print(paste("Sequences in multiple subclasses:", in.multiple))
|
|
|
167 print(paste("Multiple sequences in one subclass:", multiple.in.one))
|
|
102
|
168 print(paste("Matched with unmatched:", some.unmatched))
|
|
99
|
169 print(paste("Count that should match 'matched' sequences:", matched))
|
|
81
|
170
|
|
|
171 #ACGT overview
|
|
|
172
|
|
103
|
173 NToverview = merged[!grepl("^unmatched", merged$best_match),]
|
|
100
|
174
|
|
90
|
175 NToverview$seq = paste(NToverview$CDR1.IMGT.seq, NToverview$FR2.IMGT.seq, NToverview$CDR2.IMGT.seq, NToverview$FR3.IMGT.seq, sep="_")
|
|
82
|
176
|
|
81
|
177 NToverview$A = nchar(gsub("[^Aa]", "", NToverview$seq))
|
|
|
178 NToverview$C = nchar(gsub("[^Cc]", "", NToverview$seq))
|
|
|
179 NToverview$G = nchar(gsub("[^Gg]", "", NToverview$seq))
|
|
|
180 NToverview$T = nchar(gsub("[^Tt]", "", NToverview$seq))
|
|
|
181
|
|
82
|
182 #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))
|
|
|
183
|
|
|
184 #NToverview = rbind(NToverview, NTsum)
|
|
|
185
|
|
|
186 NTresult = data.frame(nt=c("A", "C", "T", "G"))
|
|
|
187
|
|
|
188 for(clazz in gene.classes){
|
|
91
|
189 NToverview.sub = NToverview[grepl(paste("^", clazz, sep=""), NToverview$best_match),]
|
|
82
|
190 new.col.x = c(sum(NToverview.sub$A), sum(NToverview.sub$C), sum(NToverview.sub$T), sum(NToverview.sub$G))
|
|
|
191 new.col.y = sum(new.col.x)
|
|
|
192 new.col.z = round(new.col.x / new.col.y * 100, 2)
|
|
|
193
|
|
|
194 tmp = names(NTresult)
|
|
|
195 NTresult = cbind(NTresult, data.frame(new.col.x, new.col.y, new.col.z))
|
|
|
196 names(NTresult) = c(tmp, paste(clazz, c("x", "y", "z"), sep=""))
|
|
|
197 }
|
|
81
|
198
|
|
98
|
199 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)
|
|
|
200
|
|
|
201 NToverview = NToverview[!grepl("unmatched", NToverview$best_match),]
|
|
|
202
|
|
82
|
203 new.col.x = c(sum(NToverview$A), sum(NToverview$C), sum(NToverview$T), sum(NToverview$G))
|
|
|
204 new.col.y = sum(new.col.x)
|
|
|
205 new.col.z = round(new.col.x / new.col.y * 100, 2)
|
|
|
206
|
|
|
207 tmp = names(NTresult)
|
|
|
208 NTresult = cbind(NTresult, data.frame(new.col.x, new.col.y, new.col.z))
|
|
|
209 names(NTresult) = c(tmp, paste("all", c("x", "y", "z"), sep=""))
|
|
81
|
210
|
|
82
|
211 names(hotspot.analysis.sum) = names(NTresult)
|
|
|
212
|
|
|
213 hotspot.analysis.sum = rbind(hotspot.analysis.sum, NTresult)
|
|
|
214
|
|
|
215 write.table(hotspot.analysis.sum, hotspot.analysis.sum.file, quote=F, sep=",", row.names=F, col.names=F, na="0")
|
|
|
216
|
|
81
|
217
|
|
|
218
|
|
|
219
|
|
|
220
|
|
|
221
|
|
|
222
|
|
|
223
|
|
|
224
|
|
|
225
|
|
|
226
|
|
|
227
|
|
|
228
|
|
|
229
|
|
|
230
|
|
|
231
|
|
|
232
|
|
|
233
|
|
|
234
|
|
|
235
|
|
|
236
|
|
|
237
|
|
|
238
|
|
|
239
|
|
|
240
|
|
|
241
|
|
|
242
|
|
|
243
|
|
|
244
|
|
|
245
|
