Mercurial > repos > davidvanzessen > mutation_analysis
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 |