Mercurial > repos > davidvanzessen > shm_csr

comparison shm_csr.r @ 0:c33d93683a09 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author davidvanzessen
date Thu, 13 Oct 2016 10:52:24 -0400
parents
children faae21ba5c63
comparison
equal deleted inserted replaced
-1:000000000000 0:c33d93683a09
1 library(data.table)
2 library(ggplot2)
3 library(reshape2)
4
5 args <- commandArgs(trailingOnly = TRUE)
6
7 input = args[1]
8 genes = unlist(strsplit(args[2], ","))
9 outputdir = args[3]
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 cleanup_columns = c("FR1.IMGT.c.a",
32 "FR2.IMGT.g.t",
33 "CDR1.IMGT.Nb.of.nucleotides",
34 "CDR2.IMGT.t.a",
35 "FR1.IMGT.c.g",
36 "CDR1.IMGT.c.t",
37 "FR2.IMGT.a.c",
38 "FR2.IMGT.Nb.of.mutations",
39 "FR2.IMGT.g.c",
40 "FR2.IMGT.a.g",
41 "FR3.IMGT.t.a",
42 "FR3.IMGT.t.c",
43 "FR2.IMGT.g.a",
44 "FR3.IMGT.c.g",
45 "FR1.IMGT.Nb.of.mutations",
46 "CDR1.IMGT.g.a",
47 "CDR1.IMGT.t.g",
48 "CDR1.IMGT.g.c",
49 "CDR2.IMGT.Nb.of.nucleotides",
50 "FR2.IMGT.a.t",
51 "CDR1.IMGT.Nb.of.mutations",
52 "CDR3.IMGT.Nb.of.nucleotides",
53 "CDR1.IMGT.a.g",
54 "FR3.IMGT.a.c",
55 "FR1.IMGT.g.a",
56 "FR3.IMGT.a.g",
57 "FR1.IMGT.a.t",
58 "CDR2.IMGT.a.g",
59 "CDR2.IMGT.Nb.of.mutations",
60 "CDR2.IMGT.g.t",
61 "CDR2.IMGT.a.c",
62 "CDR1.IMGT.t.c",
63 "FR3.IMGT.g.c",
64 "FR1.IMGT.g.t",
65 "FR3.IMGT.g.t",
66 "CDR1.IMGT.a.t",
67 "FR1.IMGT.a.g",
68 "FR3.IMGT.a.t",
69 "FR3.IMGT.Nb.of.nucleotides",
70 "FR2.IMGT.t.c",
71 "CDR2.IMGT.g.a",
72 "FR2.IMGT.t.a",
73 "CDR1.IMGT.t.a",
74 "FR2.IMGT.t.g",
75 "FR3.IMGT.t.g",
76 "FR2.IMGT.Nb.of.nucleotides",
77 "FR1.IMGT.t.a",
78 "FR1.IMGT.t.g",
79 "FR3.IMGT.c.t",
80 "FR1.IMGT.t.c",
81 "CDR2.IMGT.a.t",
82 "FR2.IMGT.c.t",
83 "CDR1.IMGT.g.t",
84 "CDR2.IMGT.t.g",
85 "FR1.IMGT.Nb.of.nucleotides",
86 "CDR1.IMGT.c.g",
87 "CDR2.IMGT.t.c",
88 "FR3.IMGT.g.a",
89 "CDR1.IMGT.a.c",
90 "FR2.IMGT.c.a",
91 "FR3.IMGT.Nb.of.mutations",
92 "FR2.IMGT.c.g",
93 "CDR2.IMGT.g.c",
94 "FR1.IMGT.g.c",
95 "CDR2.IMGT.c.t",
96 "FR3.IMGT.c.a",
97 "CDR1.IMGT.c.a",
98 "CDR2.IMGT.c.g",
99 "CDR2.IMGT.c.a",
100 "FR1.IMGT.c.t",
101 "FR1.IMGT.Nb.of.silent.mutations",
102 "FR2.IMGT.Nb.of.silent.mutations",
103 "FR3.IMGT.Nb.of.silent.mutations",
104 "FR1.IMGT.Nb.of.nonsilent.mutations",
105 "FR2.IMGT.Nb.of.nonsilent.mutations",
106 "FR3.IMGT.Nb.of.nonsilent.mutations")
107
108
109 print("Cleaning up columns")
110 for(col in cleanup_columns){
111 dat[,col] = gsub("\\(.*\\)", "", dat[,col])
112 #dat[dat[,col] == "",] = "0"
113 dat[,col] = as.numeric(dat[,col])
114 dat[is.na(dat[,col]),col] = 0
115 }
116
117 regions = c("FR1", "CDR1", "FR2", "CDR2", "FR3")
118 if(!include_fr1){
119 regions = c("CDR1", "FR2", "CDR2", "FR3")
120 }
121
122 sum_by_row = function(x, columns) { sum(as.numeric(x[columns]), na.rm=T) }
123
124 print("aggregating data into new columns")
125
126 VRegionMutations_columns = paste(regions, ".IMGT.Nb.of.mutations", sep="")
127 dat$VRegionMutations = apply(dat, FUN=sum_by_row, 1, columns=VRegionMutations_columns)
128
129 VRegionNucleotides_columns = paste(regions, ".IMGT.Nb.of.nucleotides", sep="")
130 dat$FR3.IMGT.Nb.of.nucleotides = nchar(dat$FR3.IMGT.seq)
131 dat$VRegionNucleotides = apply(dat, FUN=sum_by_row, 1, columns=VRegionNucleotides_columns)
132
133 transitionMutations_columns = paste(rep(regions, each=4), c(".IMGT.a.g", ".IMGT.g.a", ".IMGT.c.t", ".IMGT.t.c"), sep="")
134 dat$transitionMutations = apply(dat, FUN=sum_by_row, 1, columns=transitionMutations_columns)
135
136 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="")
137 dat$transversionMutations = apply(dat, FUN=sum_by_row, 1, columns=transversionMutations_columns)
138
139
140 transitionMutationsAtGC_columns = paste(rep(regions, each=2), c(".IMGT.g.a",".IMGT.c.t"), sep="")
141 dat$transitionMutationsAtGC = apply(dat, FUN=sum_by_row, 1, columns=transitionMutationsAtGC_columns)
142
143
144 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="")
145 #totalMutationsAtGC_columns = paste(rep(regions, each=6), c(".IMGT.g.a",".IMGT.c.t",".IMGT.c.a",".IMGT.c.g",".IMGT.g.t"), sep="")
146 dat$totalMutationsAtGC = apply(dat, FUN=sum_by_row, 1, columns=totalMutationsAtGC_columns)
147
148 transitionMutationsAtAT_columns = paste(rep(regions, each=2), c(".IMGT.a.g",".IMGT.t.c"), sep="")
149 dat$transitionMutationsAtAT = apply(dat, FUN=sum_by_row, 1, columns=transitionMutationsAtAT_columns)
150
151 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="")
152 #totalMutationsAtAT_columns = paste(rep(regions, each=5), c(".IMGT.a.g",".IMGT.t.c",".IMGT.a.c",".IMGT.g.c",".IMGT.t.g"), sep="")
153 dat$totalMutationsAtAT = apply(dat, FUN=sum_by_row, 1, columns=totalMutationsAtAT_columns)
154
155
156 FRRegions = regions[grepl("FR", regions)]
157 CDRRegions = regions[grepl("CDR", regions)]
158
159 FR_silentMutations_columns = paste(FRRegions, ".IMGT.Nb.of.silent.mutations", sep="")
160 dat$silentMutationsFR = apply(dat, FUN=sum_by_row, 1, columns=FR_silentMutations_columns)
161
162 CDR_silentMutations_columns = paste(CDRRegions, ".IMGT.Nb.of.silent.mutations", sep="")
163 dat$silentMutationsCDR = apply(dat, FUN=sum_by_row, 1, columns=CDR_silentMutations_columns)
164
165 FR_nonSilentMutations_columns = paste(FRRegions, ".IMGT.Nb.of.nonsilent.mutations", sep="")
166 dat$nonSilentMutationsFR = apply(dat, FUN=sum_by_row, 1, columns=FR_nonSilentMutations_columns)
167
168 CDR_nonSilentMutations_columns = paste(CDRRegions, ".IMGT.Nb.of.nonsilent.mutations", sep="")
169 dat$nonSilentMutationsCDR = apply(dat, FUN=sum_by_row, 1, columns=CDR_nonSilentMutations_columns)
170
171 mutation.sum.columns = c("Sequence.ID", "VRegionMutations", "VRegionNucleotides", "transitionMutations", "transversionMutations", "transitionMutationsAtGC", "transitionMutationsAtAT", "silentMutationsFR", "nonSilentMutationsFR", "silentMutationsCDR", "nonSilentMutationsCDR")
172
173 write.table(dat[,mutation.sum.columns], "mutation_by_id.txt", sep="\t",quote=F,row.names=F,col.names=T)
174
175 setwd(outputdir)
176
177 base.order = data.frame(base=c("A", "T", "C", "G"), order=1:4)
178
179 calculate_result = function(i, gene, dat, matrx, f, fname, name){
180 tmp = dat[grepl(paste("^", gene, ".*", sep=""), dat$best_match),]
181
182 j = i - 1
183 x = (j * 3) + 1
184 y = (j * 3) + 2
185 z = (j * 3) + 3
186
187 if(nrow(tmp) > 0){
188
189 if(fname == "sum"){
190 matrx[1,x] = round(f(tmp$VRegionMutations, na.rm=T), digits=1)
191 matrx[1,y] = round(f(tmp$VRegionNucleotides, na.rm=T), digits=1)
192 matrx[1,z] = round(f(matrx[1,x] / matrx[1,y]) * 100, digits=1)
193 } else {
194 matrx[1,x] = round(f(tmp$VRegionMutations, na.rm=T), digits=1)
195 matrx[1,y] = round(f(tmp$VRegionNucleotides, na.rm=T), digits=1)
196 matrx[1,z] = round(f(tmp$VRegionMutations / tmp$VRegionNucleotides) * 100, digits=1)
197 }
198
199 matrx[2,x] = round(f(tmp$transitionMutations, na.rm=T), digits=1)
200 matrx[2,y] = round(f(tmp$VRegionMutations, na.rm=T), digits=1)
201 matrx[2,z] = round(matrx[2,x] / matrx[2,y] * 100, digits=1)
202
203 matrx[3,x] = round(f(tmp$transversionMutations, na.rm=T), digits=1)
204 matrx[3,y] = round(f(tmp$VRegionMutations, na.rm=T), digits=1)
205 matrx[3,z] = round(matrx[3,x] / matrx[3,y] * 100, digits=1)
206
207 matrx[4,x] = round(f(tmp$transitionMutationsAtGC, na.rm=T), digits=1)
208 matrx[4,y] = round(f(tmp$totalMutationsAtGC, na.rm=T), digits=1)
209 matrx[4,z] = round(matrx[4,x] / matrx[4,y] * 100, digits=1)
210
211 matrx[5,x] = round(f(tmp$totalMutationsAtGC, na.rm=T), digits=1)
212 matrx[5,y] = round(f(tmp$VRegionMutations, na.rm=T), digits=1)
213 matrx[5,z] = round(matrx[5,x] / matrx[5,y] * 100, digits=1)
214
215 matrx[6,x] = round(f(tmp$transitionMutationsAtAT, na.rm=T), digits=1)
216 matrx[6,y] = round(f(tmp$totalMutationsAtAT, na.rm=T), digits=1)
217 matrx[6,z] = round(matrx[6,x] / matrx[6,y] * 100, digits=1)
218
219 matrx[7,x] = round(f(tmp$totalMutationsAtAT, na.rm=T), digits=1)
220 matrx[7,y] = round(f(tmp$VRegionMutations, na.rm=T), digits=1)
221 matrx[7,z] = round(matrx[7,x] / matrx[7,y] * 100, digits=1)
222
223 matrx[8,x] = round(f(tmp$nonSilentMutationsFR, na.rm=T), digits=1)
224 matrx[8,y] = round(f(tmp$silentMutationsFR, na.rm=T), digits=1)
225 matrx[8,z] = round(matrx[8,x] / matrx[8,y], digits=1)
226
227 matrx[9,x] = round(f(tmp$nonSilentMutationsCDR, na.rm=T), digits=1)
228 matrx[9,y] = round(f(tmp$silentMutationsCDR, na.rm=T), digits=1)
229 matrx[9,z] = round(matrx[9,x] / matrx[9,y], digits=1)
230
231 if(fname == "sum"){
232 matrx[10,x] = round(f(rowSums(tmp[,c("FR2.IMGT.Nb.of.nucleotides", "FR3.IMGT.Nb.of.nucleotides")], na.rm=T)), digits=1)
233 matrx[10,y] = round(f(tmp$VRegionNucleotides, na.rm=T), digits=1)
234 matrx[10,z] = round(matrx[10,x] / matrx[10,y] * 100, digits=1)
235
236 matrx[11,x] = round(f(rowSums(tmp[,c("CDR1.IMGT.Nb.of.nucleotides", "CDR2.IMGT.Nb.of.nucleotides")], na.rm=T)), digits=1)
237 matrx[11,y] = round(f(tmp$VRegionNucleotides, na.rm=T), digits=1)
238 matrx[11,z] = round(matrx[11,x] / matrx[11,y] * 100, digits=1)
239 }
240 }
241
242 transitionTable = data.frame(A=zeros,C=zeros,G=zeros,T=zeros)
243 row.names(transitionTable) = c("A", "C", "G", "T")
244 transitionTable["A","A"] = NA
245 transitionTable["C","C"] = NA
246 transitionTable["G","G"] = NA
247 transitionTable["T","T"] = NA
248
249 if(nrow(tmp) > 0){
250 for(nt1 in nts){
251 for(nt2 in nts){
252 if(nt1 == nt2){
253 next
254 }
255 NT1 = LETTERS[letters == nt1]
256 NT2 = LETTERS[letters == nt2]
257 FR1 = paste("FR1.IMGT.", nt1, ".", nt2, sep="")
258 CDR1 = paste("CDR1.IMGT.", nt1, ".", nt2, sep="")
259 FR2 = paste("FR2.IMGT.", nt1, ".", nt2, sep="")
260 CDR2 = paste("CDR2.IMGT.", nt1, ".", nt2, sep="")
261 FR3 = paste("FR3.IMGT.", nt1, ".", nt2, sep="")
262 if(include_fr1){
263 transitionTable[NT1,NT2] = sum(tmp[,c(FR1, CDR1, FR2, CDR2, FR3)])
264 } else {
265 transitionTable[NT1,NT2] = sum(tmp[,c(CDR1, FR2, CDR2, FR3)])
266 }
267 }
268 }
269 transition = transitionTable
270 transition$id = names(transition)
271
272 transition2 = melt(transition, id.vars="id")
273
274 transition2 = merge(transition2, base.order, by.x="id", by.y="base")
275 transition2 = merge(transition2, base.order, by.x="variable", by.y="base")
276
277 transition2[is.na(transition2$value),]$value = 0
278
279 if(!all(transition2$value == 0)){ #having rows of data but a transition table filled with 0 is bad
280
281 print("Plotting stacked transition")
282
283 png(filename=paste("transitions_stacked_", name, ".png", sep=""))
284 p = ggplot(transition2, aes(factor(reorder(id, order.x)), y=value, fill=factor(reorder(variable, order.y)))) + geom_bar(position="fill", stat="identity", colour="black") #stacked bar
285 p = p + xlab("From base") + ylab("To base") + ggtitle("Mutations frequency from base to base") + guides(fill=guide_legend(title=NULL))
286 p = p + theme(panel.background = element_rect(fill = "white", colour="black")) + scale_fill_manual(values=c("A" = "blue4", "G" = "lightblue1", "C" = "olivedrab3", "T" = "olivedrab4"))
287 #p = p + scale_colour_manual(values=c("A" = "black", "G" = "black", "C" = "black", "T" = "black"))
288 print(p)
289 dev.off()
290
291 print("Plotting heatmap transition")
292
293 png(filename=paste("transitions_heatmap_", name, ".png", sep=""))
294 p = ggplot(transition2, aes(factor(reorder(id, order.x)), factor(reorder(variable, order.y)))) + geom_tile(aes(fill = value)) + scale_fill_gradient(low="white", high="steelblue") #heatmap
295 p = p + xlab("From base") + ylab("To base") + ggtitle("Mutations frequency from base to base") + theme(panel.background = element_rect(fill = "white", colour="black"))
296 print(p)
297 dev.off()
298 } else {
299 print("No data to plot")
300 }
301 }
302
303 #print(paste("writing value file: ", name, "_", fname, "_value.txt" ,sep=""))
304
305 write.table(x=transitionTable, file=paste("transitions_", name ,"_", fname, ".txt", sep=""), sep=",",quote=F,row.names=T,col.names=NA)
306 write.table(x=tmp[,c("Sequence.ID", "best_match", "chunk_hit_percentage", "nt_hit_percentage", "start_locations")], file=paste("matched_", name , "_", fname, ".txt", sep=""), sep="\t",quote=F,row.names=F,col.names=T)
307
308 cat(matrx[1,x], file=paste(name, "_", fname, "_value.txt" ,sep=""))
309 cat(nrow(tmp), file=paste(name, "_", fname, "_n.txt" ,sep=""))
310
311 #print(paste(fname, name, nrow(tmp)))
312
313 matrx
314 }
315
316 nts = c("a", "c", "g", "t")
317 zeros=rep(0, 4)
318
319 funcs = c(median, sum, mean)
320 fnames = c("median", "sum", "mean")
321
322 print("Creating result tables")
323
324 for(i in 1:length(funcs)){
325 func = funcs[[i]]
326 fname = fnames[[i]]
327
328 rows = 9
329 if(fname == "sum"){
330 rows = 11
331 }
332 matrx = matrix(data = 0, ncol=((length(genes) + 1) * 3),nrow=rows)
333
334 for(i in 1:length(genes)){
335 print(paste("Creating table for", fname, genes[i]))
336 matrx = calculate_result(i, genes[i], dat, matrx, func, fname, genes[i])
337 }
338
339 matrx = calculate_result(i + 1, ".*", dat[!grepl("unmatched", dat$best_match),], matrx, func, fname, name="all")
340
341 result = data.frame(matrx)
342 if(fname == "sum"){
343 row.names(result) = c("Number of Mutations (%)", "Transitions (%)", "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)", "nt in FR", "nt in CDR")
344 } else {
345 row.names(result) = c("Number of Mutations (%)", "Transitions (%)", "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)")
346 }
347
348 write.table(x=result, file=paste("mutations_", fname, ".txt", sep=""), sep=",",quote=F,row.names=T,col.names=F)
349 }
350
351 print("Adding median number of mutations to sum table")
352
353 sum.table = read.table("mutations_sum.txt", sep=",", header=F)
354 median.table = read.table("mutations_median.txt", sep=",", header=F)
355
356 new.table = sum.table[1,]
357 new.table[2,] = median.table[1,]
358 new.table[3:12,] = sum.table[2:11,]
359 new.table[,1] = as.character(new.table[,1])
360 new.table[2,1] = "Median of Number of Mutations (%)"
361
362 #sum.table = sum.table[c("Number of Mutations (%)", "Median of 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)", "nt in FR", "nt in CDR"),]
363
364 write.table(x=new.table, file="mutations_sum.txt", sep=",",quote=F,row.names=F,col.names=F)
365
366
367 print("Plotting IGA piechart")
368
369 dat = dat[!grepl("^unmatched", dat$best_match),]
370
371 #blegh
372 genesForPlot = dat[grepl("IGA", dat$best_match),]$best_match
373 if(length(genesForPlot) > 0){
374 genesForPlot = data.frame(table(genesForPlot))
375 colnames(genesForPlot) = c("Gene","Freq")
376 genesForPlot$label = paste(genesForPlot$Gene, "-", genesForPlot$Freq)
377
378 pc = ggplot(genesForPlot, aes(x = factor(1), y=Freq, fill=Gene))
379 pc = pc + geom_bar(width = 1, stat = "identity") + scale_fill_manual(labels=genesForPlot$label, values=c("IGA1" = "lightblue1", "IGA2" = "blue4"))
380 pc = pc + coord_polar(theta="y")
381 pc = pc + theme(panel.background = element_rect(fill = "white", colour="black"))
382 pc = pc + xlab(" ") + ylab(" ") + ggtitle(paste("IGA subclasses", "( n =", sum(genesForPlot$Freq), ")"))
383 write.table(genesForPlot, "IGA.txt", sep="\t",quote=F,row.names=F,col.names=T)
384
385 png(filename="IGA.png")
386 print(pc)
387 dev.off()
388 }
389
390 print("Plotting IGG piechart")
391
392 genesForPlot = dat[grepl("IGG", dat$best_match),]$best_match
393 if(length(genesForPlot) > 0){
394 genesForPlot = data.frame(table(genesForPlot))
395 colnames(genesForPlot) = c("Gene","Freq")
396 genesForPlot$label = paste(genesForPlot$Gene, "-", genesForPlot$Freq)
397
398 pc = ggplot(genesForPlot, aes(x = factor(1), y=Freq, fill=Gene))
399 pc = pc + geom_bar(width = 1, stat = "identity") + scale_fill_manual(labels=genesForPlot$label, values=c("IGG1" = "olivedrab3", "IGG2" = "red", "IGG3" = "gold", "IGG4" = "darkred"))
400 pc = pc + coord_polar(theta="y")
401 pc = pc + theme(panel.background = element_rect(fill = "white", colour="black"))
402 pc = pc + xlab(" ") + ylab(" ") + ggtitle(paste("IGG subclasses", "( n =", sum(genesForPlot$Freq), ")"))
403 write.table(genesForPlot, "IGG.txt", sep="\t",quote=F,row.names=F,col.names=T)
404
405 png(filename="IGG.png")
406 print(pc)
407 dev.off()
408 }
409
410
411 print("Plotting scatterplot")
412
413 dat$percentage_mutations = round(dat$VRegionMutations / dat$VRegionNucleotides * 100, 2)
414
415 p = ggplot(dat, aes(best_match, percentage_mutations))
416 p = p + geom_point(aes(colour=best_match), position="jitter") + geom_boxplot(aes(middle=mean(percentage_mutations)), alpha=0.1, outlier.shape = NA)
417 p = p + xlab("Subclass") + ylab("Frequency") + ggtitle("Frequency scatter plot") + theme(panel.background = element_rect(fill = "white", colour="black"))
418 p = p + scale_fill_manual(values=c("IGA1" = "lightblue1", "IGA2" = "blue4", "IGG1" = "olivedrab3", "IGG2" = "red", "IGG3" = "gold", "IGG4" = "darkred", "IGM" = "black"))
419 p = p + scale_colour_manual(values=c("IGA1" = "lightblue1", "IGA2" = "blue4", "IGG1" = "olivedrab3", "IGG2" = "red", "IGG3" = "gold", "IGG4" = "darkred", "IGM" = "black"))
420
421 png(filename="scatter.png")
422 print(p)
423 dev.off()
424
425 write.table(dat[,c("Sequence.ID", "best_match", "VRegionMutations", "VRegionNucleotides", "percentage_mutations")], "scatter.txt", sep="\t",quote=F,row.names=F,col.names=T)
426
427 write.table(dat, input, sep="\t",quote=F,row.names=F,col.names=T)
428
429
430 print("Plotting frequency ranges plot")
431
432 dat$best_match_class = substr(dat$best_match, 0, 3)
433 freq_labels = c("0", "0-2", "2-5", "5-10", "10-15", "15-20", "20")
434 dat$frequency_bins = cut(dat$percentage_mutations, breaks=c(-Inf, 0, 2,5,10,15,20, Inf), labels=freq_labels)
435
436 frequency_bins_sum = data.frame(data.table(dat)[, list(class_sum=sum(.N)), by=c("best_match_class")])
437
438 frequency_bins_data = data.frame(data.table(dat)[, list(frequency_count=.N), by=c("best_match_class", "frequency_bins")])
439
440 frequency_bins_data = merge(frequency_bins_data, frequency_bins_sum, by="best_match_class")
441
442 frequency_bins_data$frequency = round(frequency_bins_data$frequency_count / frequency_bins_data$class_sum * 100, 2)
443
444 p = ggplot(frequency_bins_data, aes(frequency_bins, frequency))
445 p = p + geom_bar(aes(fill=best_match_class), stat="identity", position="dodge") + theme(panel.background = element_rect(fill = "white", colour="black"))
446 p = p + xlab("Frequency ranges") + ylab("Frequency") + ggtitle("Mutation Frequencies by class") + scale_fill_manual(values=c("IGA" = "blue4", "IGG" = "olivedrab3", "IGM" = "black"))
447
448 png(filename="frequency_ranges.png")
449 print(p)
450 dev.off()
451
452 frequency_bins_data_by_class = frequency_bins_data
453
454 write.table(frequency_bins_data_by_class, "frequency_ranges_classes.txt", sep="\t",quote=F,row.names=F,col.names=T)
455
456 frequency_bins_data = data.frame(data.table(dat)[, list(frequency_count=.N), by=c("best_match", "best_match_class", "frequency_bins")])
457
458 frequency_bins_data = merge(frequency_bins_data, frequency_bins_sum, by="best_match_class")
459
460 frequency_bins_data$frequency = round(frequency_bins_data$frequency_count / frequency_bins_data$class_sum * 100, 2)
461
462 write.table(frequency_bins_data, "frequency_ranges_subclasses.txt", sep="\t",quote=F,row.names=F,col.names=T)
463
464
465 #frequency_bins_data_by_class
466 #frequency_ranges_subclasses.txt
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493