comparison mutation_analysis.r.bak @ 53:7290a88ea202 draft

Uploaded
author davidvanzessen
date Mon, 29 Feb 2016 10:49:39 -0500
parents
children
comparison
equal deleted inserted replaced
52:d3542f87a304 53:7290a88ea202
1 library(data.table)
2 library(ggplot2)
3
4 args <- commandArgs(trailingOnly = TRUE)
5
6 input = args[1]
7 genes = unlist(strsplit(args[2], ","))
8 outputdir = args[3]
9 print(args[4])
10 include_fr1 = ifelse(args[4] == "yes", T, F)
11 setwd(outputdir)
12
13 dat = read.table(input, header=T, sep="\t", fill=T, stringsAsFactors=F)
14
15 if(length(dat$Sequence.ID) == 0){
16 setwd(outputdir)
17 result = data.frame(x = rep(0, 5), y = rep(0, 5), z = rep(NA, 5))
18 row.names(result) = c("Number of Mutations (%)", "Transition (%)", "Transversions (%)", "Transitions at G C (%)", "Targeting of C G (%)")
19 write.table(x=result, file="mutations.txt", sep=",",quote=F,row.names=T,col.names=F)
20 transitionTable = data.frame(A=rep(0, 4),C=rep(0, 4),G=rep(0, 4),T=rep(0, 4))
21 row.names(transitionTable) = c("A", "C", "G", "T")
22 transitionTable["A","A"] = NA
23 transitionTable["C","C"] = NA
24 transitionTable["G","G"] = NA
25 transitionTable["T","T"] = NA
26 write.table(x=transitionTable, file="transitions.txt", sep=",",quote=F,row.names=T,col.names=NA)
27 cat("0", file="n.txt")
28 stop("No data")
29 }
30
31
32
33 cleanup_columns = c("FR1.IMGT.c.a",
34 "FR2.IMGT.g.t",
35 "CDR1.IMGT.Nb.of.nucleotides",
36 "CDR2.IMGT.t.a",
37 "FR1.IMGT.c.g",
38 "CDR1.IMGT.c.t",
39 "FR2.IMGT.a.c",
40 "FR2.IMGT.Nb.of.mutations",
41 "FR2.IMGT.g.c",
42 "FR2.IMGT.a.g",
43 "FR3.IMGT.t.a",
44 "FR3.IMGT.t.c",
45 "FR2.IMGT.g.a",
46 "FR3.IMGT.c.g",
47 "FR1.IMGT.Nb.of.mutations",
48 "CDR1.IMGT.g.a",
49 "CDR1.IMGT.t.g",
50 "CDR1.IMGT.g.c",
51 "CDR2.IMGT.Nb.of.nucleotides",
52 "FR2.IMGT.a.t",
53 "CDR1.IMGT.Nb.of.mutations",
54 "CDR1.IMGT.a.g",
55 "FR3.IMGT.a.c",
56 "FR1.IMGT.g.a",
57 "FR3.IMGT.a.g",
58 "FR1.IMGT.a.t",
59 "CDR2.IMGT.a.g",
60 "CDR2.IMGT.Nb.of.mutations",
61 "CDR2.IMGT.g.t",
62 "CDR2.IMGT.a.c",
63 "CDR1.IMGT.t.c",
64 "FR3.IMGT.g.c",
65 "FR1.IMGT.g.t",
66 "FR3.IMGT.g.t",
67 "CDR1.IMGT.a.t",
68 "FR1.IMGT.a.g",
69 "FR3.IMGT.a.t",
70 "FR3.IMGT.Nb.of.nucleotides",
71 "FR2.IMGT.t.c",
72 "CDR2.IMGT.g.a",
73 "FR2.IMGT.t.a",
74 "CDR1.IMGT.t.a",
75 "FR2.IMGT.t.g",
76 "FR3.IMGT.t.g",
77 "FR2.IMGT.Nb.of.nucleotides",
78 "FR1.IMGT.t.a",
79 "FR1.IMGT.t.g",
80 "FR3.IMGT.c.t",
81 "FR1.IMGT.t.c",
82 "CDR2.IMGT.a.t",
83 "FR2.IMGT.c.t",
84 "CDR1.IMGT.g.t",
85 "CDR2.IMGT.t.g",
86 "FR1.IMGT.Nb.of.nucleotides",
87 "CDR1.IMGT.c.g",
88 "CDR2.IMGT.t.c",
89 "FR3.IMGT.g.a",
90 "CDR1.IMGT.a.c",
91 "FR2.IMGT.c.a",
92 "FR3.IMGT.Nb.of.mutations",
93 "FR2.IMGT.c.g",
94 "CDR2.IMGT.g.c",
95 "FR1.IMGT.g.c",
96 "CDR2.IMGT.c.t",
97 "FR3.IMGT.c.a",
98 "CDR1.IMGT.c.a",
99 "CDR2.IMGT.c.g",
100 "CDR2.IMGT.c.a",
101 "FR1.IMGT.c.t",
102 "FR1.IMGT.Nb.of.silent.mutations",
103 "FR2.IMGT.Nb.of.silent.mutations",
104 "FR3.IMGT.Nb.of.silent.mutations",
105 "FR1.IMGT.Nb.of.nonsilent.mutations",
106 "FR2.IMGT.Nb.of.nonsilent.mutations",
107 "FR3.IMGT.Nb.of.nonsilent.mutations")
108
109 for(col in cleanup_columns){
110 dat[,col] = gsub("\\(.*\\)", "", dat[,col])
111 #dat[dat[,col] == "",] = "0"
112 dat[,col] = as.numeric(dat[,col])
113 dat[is.na(dat[,col]),] = 0
114 }
115
116 regions = c("FR1", "CDR1", "FR2", "CDR2", "FR3")
117 if(!include_fr1){
118 regions = c("CDR1", "FR2", "CDR2", "FR3")
119 }
120
121 sum_by_row = function(x, columns) { sum(as.numeric(x[columns]), na.rm=T) }
122
123 VRegionMutations_columns = paste(regions, ".IMGT.Nb.of.mutations", sep="")
124 dat$VRegionMutations = apply(dat, FUN=sum_by_row, 1, columns=VRegionMutations_columns)
125
126 VRegionNucleotides_columns = paste(regions, ".IMGT.Nb.of.nucleotides", sep="")
127 dat$VRegionNucleotides = apply(dat, FUN=sum_by_row, 1, columns=VRegionNucleotides_columns)
128
129 transitionMutations_columns = paste(rep(regions, each=4), c(".IMGT.a.g", ".IMGT.g.a", ".IMGT.c.t", ".IMGT.t.c"), sep="")
130 dat$transitionMutations = apply(dat, FUN=sum_by_row, 1, columns=transitionMutations_columns)
131
132 transversionMutations_columns = paste(rep(regions, each=8), c(".IMGT.a.c",".IMGT.c.a",".IMGT.a.t",".IMGT.t.a",".IMGT.g.c",".IMGT.c.g",".IMGT.g.t",".IMGT.t.g"), sep="")
133 dat$transversionMutations = apply(dat, FUN=sum_by_row, 1, columns=transversionMutations_columns)
134
135
136 transitionMutationsAtGC_columns = paste(rep(regions, each=2), c(".IMGT.g.a",".IMGT.c.t"), sep="")
137 dat$transitionMutationsAtGC = apply(dat, FUN=sum_by_row, 1, columns=transitionMutationsAtGC_columns)
138
139
140 totalMutationsAtGC_columns = paste(rep(regions, each=6), c(".IMGT.c.g",".IMGT.c.t",".IMGT.c.a",".IMGT.g.c",".IMGT.g.a",".IMGT.g.t"), sep="")
141 #totalMutationsAtGC_columns = paste(rep(regions, each=6), c(".IMGT.g.a",".IMGT.c.t",".IMGT.c.a",".IMGT.c.g",".IMGT.g.t"), sep="")
142 dat$totalMutationsAtGC = apply(dat, FUN=sum_by_row, 1, columns=totalMutationsAtGC_columns)
143
144 transitionMutationsAtAT_columns = paste(rep(regions, each=2), c(".IMGT.a.g",".IMGT.t.c"), sep="")
145 dat$transitionMutationsAtAT = apply(dat, FUN=sum_by_row, 1, columns=transitionMutationsAtAT_columns)
146
147 totalMutationsAtAT_columns = paste(rep(regions, each=6), c(".IMGT.a.g",".IMGT.a.c",".IMGT.a.t",".IMGT.t.g",".IMGT.t.c",".IMGT.t.a"), sep="")
148 #totalMutationsAtAT_columns = paste(rep(regions, each=5), c(".IMGT.a.g",".IMGT.t.c",".IMGT.a.c",".IMGT.g.c",".IMGT.t.g"), sep="")
149 dat$totalMutationsAtAT = apply(dat, FUN=sum_by_row, 1, columns=totalMutationsAtAT_columns)
150
151
152 FRRegions = regions[grepl("FR", regions)]
153 CDRRegions = regions[grepl("CDR", regions)]
154
155 FR_silentMutations_columns = paste(FRRegions, ".IMGT.Nb.of.silent.mutations", sep="")
156 dat$silentMutationsFR = apply(dat, FUN=sum_by_row, 1, columns=FR_silentMutations_columns)
157
158 CDR_silentMutations_columns = paste(CDRRegions, ".IMGT.Nb.of.silent.mutations", sep="")
159 dat$silentMutationsCDR = apply(dat, FUN=sum_by_row, 1, columns=CDR_silentMutations_columns)
160
161 FR_nonSilentMutations_columns = paste(FRRegions, ".IMGT.Nb.of.nonsilent.mutations", sep="")
162 dat$nonSilentMutationsFR = apply(dat, FUN=sum_by_row, 1, columns=FR_nonSilentMutations_columns)
163
164 CDR_nonSilentMutations_columns = paste(CDRRegions, ".IMGT.Nb.of.nonsilent.mutations", sep="")
165 dat$nonSilentMutationsCDR = apply(dat, FUN=sum_by_row, 1, columns=CDR_nonSilentMutations_columns)
166
167 mutation.sum.columns = c("Sequence.ID", "VRegionMutations", "VRegionNucleotides", "transitionMutations", "transversionMutations", "transitionMutationsAtGC", "transitionMutationsAtAT", "silentMutationsFR", "nonSilentMutationsFR", "silentMutationsCDR", "nonSilentMutationsCDR")
168
169 write.table(dat[,mutation.sum.columns], "mutation_by_id.txt", sep="\t",quote=F,row.names=F,col.names=T)
170
171 setwd(outputdir)
172
173 nts = c("a", "c", "g", "t")
174 zeros=rep(0, 4)
175 matrx = matrix(data = 0, ncol=((length(genes) + 1) * 3),nrow=9)
176 for(i in 1:length(genes)){
177 gene = genes[i]
178 tmp = dat[grepl(paste(".*", gene, ".*", sep=""), dat$best_match),]
179 if(gene == "."){
180 tmp = dat
181 }
182 j = i - 1
183 x = (j * 3) + 1
184 y = (j * 3) + 2
185 z = (j * 3) + 3
186 matrx[1,x] = sum(tmp$VRegionMutations)
187 matrx[1,y] = sum(tmp$VRegionNucleotides)
188 matrx[1,z] = round(matrx[1,x] / matrx[1,y] * 100, digits=1)
189
190 matrx[2,x] = sum(tmp$transitionMutations)
191 matrx[2,y] = sum(tmp$VRegionMutations)
192 matrx[2,z] = round(matrx[2,x] / matrx[2,y] * 100, digits=1)
193
194 matrx[3,x] = sum(tmp$transversionMutations)
195 matrx[3,y] = sum(tmp$VRegionMutations)
196 matrx[3,z] = round(matrx[3,x] / matrx[3,y] * 100, digits=1)
197
198 matrx[4,x] = sum(tmp$transitionMutationsAtGC)
199 matrx[4,y] = sum(tmp$totalMutationsAtGC)
200 matrx[4,z] = round(matrx[4,x] / matrx[4,y] * 100, digits=1)
201
202 matrx[5,x] = sum(tmp$totalMutationsAtGC)
203 matrx[5,y] = sum(tmp$VRegionMutations)
204 matrx[5,z] = round(matrx[5,x] / matrx[5,y] * 100, digits=1)
205
206 matrx[6,x] = sum(tmp$transitionMutationsAtAT)
207 matrx[6,y] = sum(tmp$totalMutationsAtAT)
208 matrx[6,z] = round(matrx[6,x] / matrx[6,y] * 100, digits=1)
209
210 matrx[7,x] = sum(tmp$totalMutationsAtAT)
211 matrx[7,y] = sum(tmp$VRegionMutations)
212 matrx[7,z] = round(matrx[7,x] / matrx[7,y] * 100, digits=1)
213
214 matrx[8,x] = sum(tmp$nonSilentMutationsFR)
215 matrx[8,y] = sum(tmp$silentMutationsFR)
216 matrx[8,z] = round(matrx[8,x] / matrx[8,y], digits=1)
217
218 matrx[9,x] = sum(tmp$nonSilentMutationsCDR)
219 matrx[9,y] = sum(tmp$silentMutationsCDR)
220 matrx[9,z] = round(matrx[9,x] / matrx[9,y], digits=1)
221
222
223 transitionTable = data.frame(A=zeros,C=zeros,G=zeros,T=zeros)
224 row.names(transitionTable) = c("A", "C", "G", "T")
225 transitionTable["A","A"] = NA
226 transitionTable["C","C"] = NA
227 transitionTable["G","G"] = NA
228 transitionTable["T","T"] = NA
229
230 if(nrow(tmp) > 0){
231 for(nt1 in nts){
232 for(nt2 in nts){
233 if(nt1 == nt2){
234 next
235 }
236 NT1 = LETTERS[letters == nt1]
237 NT2 = LETTERS[letters == nt2]
238 FR1 = paste("FR1.IMGT.", nt1, ".", nt2, sep="")
239 CDR1 = paste("CDR1.IMGT.", nt1, ".", nt2, sep="")
240 FR2 = paste("FR2.IMGT.", nt1, ".", nt2, sep="")
241 CDR2 = paste("CDR2.IMGT.", nt1, ".", nt2, sep="")
242 FR3 = paste("FR3.IMGT.", nt1, ".", nt2, sep="")
243 if(include_fr1){
244 transitionTable[NT1,NT2] = sum(tmp[,c(FR1, CDR1, FR2, CDR2, FR3)])
245 } else {
246 transitionTable[NT1,NT2] = sum(tmp[,c(CDR1, FR2, CDR2, FR3)])
247 }
248 }
249 }
250 }
251
252
253 write.table(x=transitionTable, file=paste("transitions_", gene ,".txt", sep=""), sep=",",quote=F,row.names=T,col.names=NA)
254 write.table(x=tmp[,c("Sequence.ID", "best_match", "chunk_hit_percentage", "nt_hit_percentage", "start_locations")], file=paste("matched_", gene ,".txt", sep=""), sep="\t",quote=F,row.names=F,col.names=T)
255
256 cat(matrx[1,x], file=paste(gene, "_value.txt" ,sep=""))
257 cat(length(tmp$Sequence.ID), file=paste(gene, "_n.txt" ,sep=""))
258 }
259
260 #again for all of the data
261 tmp = dat
262 j = i
263 x = (j * 3) + 1
264 y = (j * 3) + 2
265 z = (j * 3) + 3
266 matrx[1,x] = sum(tmp$VRegionMutations)
267 matrx[1,y] = sum(tmp$VRegionNucleotides)
268 matrx[1,z] = round(matrx[1,x] / matrx[1,y] * 100, digits=1)
269
270 matrx[2,x] = sum(tmp$transitionMutations)
271 matrx[2,y] = sum(tmp$VRegionMutations)
272 matrx[2,z] = round(matrx[2,x] / matrx[2,y] * 100, digits=1)
273
274 matrx[3,x] = sum(tmp$transversionMutations)
275 matrx[3,y] = sum(tmp$VRegionMutations)
276 matrx[3,z] = round(matrx[3,x] / matrx[3,y] * 100, digits=1)
277
278 matrx[4,x] = sum(tmp$transitionMutationsAtGC)
279 matrx[4,y] = sum(tmp$totalMutationsAtGC)
280 matrx[4,z] = round(matrx[4,x] / matrx[4,y] * 100, digits=1)
281
282 matrx[5,x] = sum(tmp$totalMutationsAtGC)
283 matrx[5,y] = sum(tmp$VRegionMutations)
284 matrx[5,z] = round(matrx[5,x] / matrx[5,y] * 100, digits=1)
285
286 matrx[6,x] = sum(tmp$transitionMutationsAtAT)
287 matrx[6,y] = sum(tmp$totalMutationsAtAT)
288 matrx[6,z] = round(matrx[6,x] / matrx[6,y] * 100, digits=1)
289
290 matrx[7,x] = sum(tmp$totalMutationsAtAT)
291 matrx[7,y] = sum(tmp$VRegionMutations)
292 matrx[7,z] = round(matrx[7,x] / matrx[7,y] * 100, digits=1)
293
294 matrx[8,x] = sum(tmp$nonSilentMutationsFR)
295 matrx[8,y] = sum(tmp$silentMutationsFR)
296 matrx[8,z] = round(matrx[8,x] / matrx[8,y], digits=1)
297
298 matrx[9,x] = sum(tmp$nonSilentMutationsCDR)
299 matrx[9,y] = sum(tmp$silentMutationsCDR)
300 matrx[9,z] = round(matrx[9,x] / matrx[9,y], digits=1)
301
302 transitionTable = data.frame(A=1:4,C=1:4,G=1:4,T=1:4)
303 row.names(transitionTable) = c("A", "C", "G", "T")
304 transitionTable["A","A"] = NA
305 transitionTable["C","C"] = NA
306 transitionTable["G","G"] = NA
307 transitionTable["T","T"] = NA
308
309
310 for(nt1 in nts){
311 for(nt2 in nts){
312 if(nt1 == nt2){
313 next
314 }
315 NT1 = LETTERS[letters == nt1]
316 NT2 = LETTERS[letters == nt2]
317 FR1 = paste("FR1.IMGT.", nt1, ".", nt2, sep="")
318 CDR1 = paste("CDR1.IMGT.", nt1, ".", nt2, sep="")
319 FR2 = paste("FR2.IMGT.", nt1, ".", nt2, sep="")
320 CDR2 = paste("CDR2.IMGT.", nt1, ".", nt2, sep="")
321 FR3 = paste("FR3.IMGT.", nt1, ".", nt2, sep="")
322 if(include_fr1){
323 transitionTable[NT1,NT2] = sum(tmp[,c(FR1, CDR1, FR2, CDR2, FR3)])
324 } else {
325 transitionTable[NT1,NT2] = sum(tmp[,c(CDR1, FR2, CDR2, FR3)])
326 }
327 }
328 }
329 write.table(x=transitionTable, file="transitions.txt", sep=",",quote=F,row.names=T,col.names=NA)
330 write.table(x=tmp[,c("Sequence.ID", "best_match", "chunk_hit_percentage", "nt_hit_percentage", "start_locations")], file="matched_all.txt", sep="\t",quote=F,row.names=F,col.names=T)
331 cat(matrx[1,x], file="total_value.txt")
332 cat(length(tmp$Sequence.ID), file="total_n.txt")
333
334
335
336 result = data.frame(matrx)
337 row.names(result) = c("Number of Mutations (%)", "Transition (%)", "Transversions (%)", "Transitions at G C (%)", "Targeting of C G (%)", "Transitions at A T (%)", "Targeting of A T (%)", "FR R/S (ratio)", "CDR R/S (ratio)")
338
339 write.table(x=result, file="mutations.txt", sep=",",quote=F,row.names=T,col.names=F)
340
341
342 if (!("ggplot2" %in% rownames(installed.packages()))) {
343 install.packages("ggplot2", repos="http://cran.xl-mirror.nl/")
344 }
345
346
347 genesForPlot = gsub("[0-9]", "", dat$best_match)
348 genesForPlot = data.frame(table(genesForPlot))
349 colnames(genesForPlot) = c("Gene","Freq")
350 genesForPlot$label = paste(genesForPlot$Gene, "-", genesForPlot$Freq)
351 write.table(genesForPlot, "all.txt", sep="\t",quote=F,row.names=F,col.names=T)
352
353
354 pc = ggplot(genesForPlot, aes(x = factor(1), y=Freq, fill=label))
355 pc = pc + geom_bar(width = 1, stat = "identity")
356 pc = pc + coord_polar(theta="y")
357 pc = pc + xlab(" ") + ylab(" ") + ggtitle(paste("Classes", "( n =", sum(genesForPlot$Freq), ")"))
358
359 png(filename="all.png")
360 pc
361 dev.off()
362
363
364 #blegh
365 genesForPlot = dat[grepl("ca", dat$best_match),]$best_match
366 if(length(genesForPlot) > 0){
367 genesForPlot = data.frame(table(genesForPlot))
368 colnames(genesForPlot) = c("Gene","Freq")
369 genesForPlot$label = paste(genesForPlot$Gene, "-", genesForPlot$Freq)
370
371 pc = ggplot(genesForPlot, aes(x = factor(1), y=Freq, fill=label))
372 pc = pc + geom_bar(width = 1, stat = "identity")
373 pc = pc + coord_polar(theta="y")
374 pc = pc + xlab(" ") + ylab(" ") + ggtitle(paste("IgA subclasses", "( n =", sum(genesForPlot$Freq), ")"))
375 write.table(genesForPlot, "ca.txt", sep="\t",quote=F,row.names=F,col.names=T)
376
377 png(filename="ca.png")
378 print(pc)
379 dev.off()
380 }
381
382 genesForPlot = dat[grepl("cg", dat$best_match),]$best_match
383 if(length(genesForPlot) > 0){
384 genesForPlot = data.frame(table(genesForPlot))
385 colnames(genesForPlot) = c("Gene","Freq")
386 genesForPlot$label = paste(genesForPlot$Gene, "-", genesForPlot$Freq)
387
388 pc = ggplot(genesForPlot, aes(x = factor(1), y=Freq, fill=label))
389 pc = pc + geom_bar(width = 1, stat = "identity")
390 pc = pc + coord_polar(theta="y")
391 pc = pc + xlab(" ") + ylab(" ") + ggtitle(paste("IgG subclasses", "( n =", sum(genesForPlot$Freq), ")"))
392 write.table(genesForPlot, "cg.txt", sep="\t",quote=F,row.names=F,col.names=T)
393
394 png(filename="cg.png")
395 print(pc)
396 dev.off()
397 }
398
399 dat$percentage_mutations = round(dat$VRegionMutations / dat$VRegionNucleotides * 100, 2)
400
401 p = ggplot(dat, aes(best_match, percentage_mutations))
402 p = p + geom_point(aes(colour=best_match), position="jitter") + geom_boxplot(aes(middle=mean(percentage_mutations)), alpha=0.1, outlier.shape = NA)
403 p = p + xlab("Subclass") + ylab("Frequency") + ggtitle("Frequency scatter plot")
404
405 png(filename="scatter.png")
406 print(p)
407 dev.off()
408
409 write.table(dat[,c("Sequence.ID", "best_match", "VRegionMutations", "VRegionNucleotides", "percentage_mutations")], "scatter.txt", sep="\t",quote=F,row.names=F,col.names=T)
410
411 write.table(dat, input, sep="\t",quote=F,row.names=F,col.names=T)
412
413
414
415
416
417
418 dat$best_match_class = substr(dat$best_match, 0, 2)
419 freq_labels = c("0", "0-2", "2-5", "5-10", "10-15", "15-20", "20")
420 dat$frequency_bins = cut(dat$percentage_mutations, breaks=c(-Inf, 0, 2,5,10,15,20, Inf), labels=freq_labels)
421
422 frequency_bins_data = data.frame(data.table(dat)[, list(frequency_count=.N), by=c("best_match_class", "frequency_bins")])
423
424 p = ggplot(frequency_bins_data, aes(frequency_bins, frequency_count))
425 p = p + geom_bar(aes(fill=best_match_class), stat="identity", position="dodge")
426 p = p + xlab("Frequency ranges") + ylab("Frequency") + ggtitle("Mutation Frequencies by class")
427
428 png(filename="frequency_ranges.png")
429 print(p)
430 dev.off()
431
432 frequency_bins_data_by_class = frequency_bins_data
433
434 write.table(frequency_bins_data_by_class, "frequency_ranges_classes.txt", sep="\t",quote=F,row.names=F,col.names=T)
435
436 frequency_bins_data = data.frame(data.table(dat)[, list(frequency_count=.N), by=c("best_match", "frequency_bins")])
437
438 write.table(frequency_bins_data, "frequency_ranges_subclasses.txt", sep="\t",quote=F,row.names=F,col.names=T)
439
440
441 #frequency_bins_data_by_class
442 #frequency_ranges_subclasses.txt
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469