Mercurial > repos > davidvanzessen > mutation_analysis

comparison sequence_overview.r @ 0:8a5a2abbb870 draft default tip

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