|
0
|
1 args <- commandArgs(trailingOnly = TRUE)
|
|
|
2
|
|
|
3 input = args[1]
|
|
4
|
4 genes = unlist(strsplit(args[2], ","))
|
|
0
|
5 outputdir = args[3]
|
|
|
6 setwd(outputdir)
|
|
|
7
|
|
|
8 dat = read.table(input, header=T, sep="\t", fill=T, stringsAsFactors=F)
|
|
|
9
|
|
|
10 if(length(dat$Sequence.ID) == 0){
|
|
4
|
11 setwd(outputdir)
|
|
|
12 result = data.frame(x = rep(0, 5), y = rep(0, 5), z = rep(NA, 5))
|
|
|
13 row.names(result) = c("Number of Mutations (%)", "Transition (%)", "Transversions (%)", "Transitions at G C (%)", "Targeting of C G (%)")
|
|
|
14 write.table(x=result, file="mutations.txt", sep=",",quote=F,row.names=T,col.names=F)
|
|
|
15 transitionTable = data.frame(A=rep(0, 4),C=rep(0, 4),G=rep(0, 4),T=rep(0, 4))
|
|
|
16 row.names(transitionTable) = c("A", "C", "G", "T")
|
|
|
17 transitionTable["A","A"] = NA
|
|
|
18 transitionTable["C","C"] = NA
|
|
|
19 transitionTable["G","G"] = NA
|
|
|
20 transitionTable["T","T"] = NA
|
|
|
21 write.table(x=transitionTable, file="transitions.txt", sep=",",quote=F,row.names=T,col.names=NA)
|
|
|
22 cat("0", file="n.txt")
|
|
|
23 stop("No data")
|
|
0
|
24 }
|
|
|
25
|
|
|
26
|
|
|
27
|
|
|
28 cleanup_columns = c("FR1.IMGT.c.a",
|
|
|
29 "FR2.IMGT.g.t",
|
|
|
30 "CDR1.IMGT.Nb.of.nucleotides",
|
|
|
31 "CDR2.IMGT.t.a",
|
|
|
32 "FR1.IMGT.c.g",
|
|
|
33 "CDR1.IMGT.c.t",
|
|
|
34 "FR2.IMGT.a.c",
|
|
|
35 "FR2.IMGT.Nb.of.mutations",
|
|
|
36 "FR2.IMGT.g.c",
|
|
|
37 "FR2.IMGT.a.g",
|
|
|
38 "FR3.IMGT.t.a",
|
|
|
39 "FR3.IMGT.t.c",
|
|
|
40 "FR2.IMGT.g.a",
|
|
|
41 "FR3.IMGT.c.g",
|
|
|
42 "FR1.IMGT.Nb.of.mutations",
|
|
|
43 "CDR1.IMGT.g.a",
|
|
|
44 "CDR1.IMGT.t.g",
|
|
|
45 "CDR1.IMGT.g.c",
|
|
|
46 "CDR2.IMGT.Nb.of.nucleotides",
|
|
|
47 "FR2.IMGT.a.t",
|
|
|
48 "CDR1.IMGT.Nb.of.mutations",
|
|
|
49 "CDR1.IMGT.a.g",
|
|
|
50 "FR3.IMGT.a.c",
|
|
|
51 "FR1.IMGT.g.a",
|
|
|
52 "FR3.IMGT.a.g",
|
|
|
53 "FR1.IMGT.a.t",
|
|
|
54 "CDR2.IMGT.a.g",
|
|
|
55 "CDR2.IMGT.Nb.of.mutations",
|
|
|
56 "CDR2.IMGT.g.t",
|
|
|
57 "CDR2.IMGT.a.c",
|
|
|
58 "CDR1.IMGT.t.c",
|
|
|
59 "FR3.IMGT.g.c",
|
|
|
60 "FR1.IMGT.g.t",
|
|
|
61 "FR3.IMGT.g.t",
|
|
|
62 "CDR1.IMGT.a.t",
|
|
|
63 "FR1.IMGT.a.g",
|
|
|
64 "FR3.IMGT.a.t",
|
|
|
65 "FR3.IMGT.Nb.of.nucleotides",
|
|
|
66 "FR2.IMGT.t.c",
|
|
|
67 "CDR2.IMGT.g.a",
|
|
|
68 "FR2.IMGT.t.a",
|
|
|
69 "CDR1.IMGT.t.a",
|
|
|
70 "FR2.IMGT.t.g",
|
|
|
71 "FR3.IMGT.t.g",
|
|
|
72 "FR2.IMGT.Nb.of.nucleotides",
|
|
|
73 "FR1.IMGT.t.a",
|
|
|
74 "FR1.IMGT.t.g",
|
|
|
75 "FR3.IMGT.c.t",
|
|
|
76 "FR1.IMGT.t.c",
|
|
|
77 "CDR2.IMGT.a.t",
|
|
|
78 "FR2.IMGT.c.t",
|
|
|
79 "CDR1.IMGT.g.t",
|
|
|
80 "CDR2.IMGT.t.g",
|
|
|
81 "FR1.IMGT.Nb.of.nucleotides",
|
|
|
82 "CDR1.IMGT.c.g",
|
|
|
83 "CDR2.IMGT.t.c",
|
|
|
84 "FR3.IMGT.g.a",
|
|
|
85 "CDR1.IMGT.a.c",
|
|
|
86 "FR2.IMGT.c.a",
|
|
|
87 "FR3.IMGT.Nb.of.mutations",
|
|
|
88 "FR2.IMGT.c.g",
|
|
|
89 "CDR2.IMGT.g.c",
|
|
|
90 "FR1.IMGT.g.c",
|
|
|
91 "CDR2.IMGT.c.t",
|
|
|
92 "FR3.IMGT.c.a",
|
|
|
93 "CDR1.IMGT.c.a",
|
|
|
94 "CDR2.IMGT.c.g",
|
|
|
95 "CDR2.IMGT.c.a",
|
|
|
96 "FR1.IMGT.c.t")
|
|
|
97
|
|
|
98 for(col in cleanup_columns){
|
|
|
99 dat[,col] = gsub("\\(.*\\)", "", dat[,col])
|
|
|
100 #dat[dat[,col] == "",] = "0"
|
|
|
101 dat[,col] = as.numeric(dat[,col])
|
|
|
102 dat[is.na(dat[,col]),] = 0
|
|
|
103 }
|
|
|
104
|
|
7
|
105 dat$VRegionMutations = dat$CDR1.IMGT.Nb.of.mutations +
|
|
0
|
106 dat$FR2.IMGT.Nb.of.mutations +
|
|
|
107 dat$CDR2.IMGT.Nb.of.mutations +
|
|
4
|
108 dat$FR3.IMGT.Nb.of.mutations
|
|
0
|
109
|
|
7
|
110 dat$VRegionNucleotides = dat$CDR1.IMGT.Nb.of.nucleotides +
|
|
0
|
111 dat$FR2.IMGT.Nb.of.nucleotides +
|
|
|
112 dat$CDR2.IMGT.Nb.of.nucleotides +
|
|
4
|
113 dat$FR3.IMGT.Nb.of.nucleotides
|
|
0
|
114
|
|
7
|
115 dat$transitionMutations = dat$CDR1.IMGT.a.g +
|
|
0
|
116 dat$CDR1.IMGT.g.a +
|
|
|
117 dat$CDR1.IMGT.c.t +
|
|
|
118 dat$CDR1.IMGT.t.c +
|
|
|
119 dat$FR2.IMGT.a.g +
|
|
|
120 dat$FR2.IMGT.g.a +
|
|
|
121 dat$FR2.IMGT.c.t +
|
|
|
122 dat$FR2.IMGT.t.c +
|
|
|
123 dat$CDR2.IMGT.a.g +
|
|
|
124 dat$CDR2.IMGT.g.a +
|
|
|
125 dat$CDR2.IMGT.c.t +
|
|
|
126 dat$CDR2.IMGT.t.c +
|
|
|
127 dat$FR3.IMGT.a.g +
|
|
|
128 dat$FR3.IMGT.g.a +
|
|
|
129 dat$FR3.IMGT.c.t +
|
|
4
|
130 dat$FR3.IMGT.t.c
|
|
0
|
131
|
|
7
|
132 dat$transversionMutations = dat$CDR1.IMGT.a.c +
|
|
4
|
133 dat$CDR1.IMGT.c.a +
|
|
|
134 dat$CDR1.IMGT.a.t +
|
|
|
135 dat$CDR1.IMGT.t.a +
|
|
|
136 dat$CDR1.IMGT.g.c +
|
|
|
137 dat$CDR1.IMGT.c.g +
|
|
|
138 dat$CDR1.IMGT.g.t +
|
|
|
139 dat$CDR1.IMGT.t.g +
|
|
|
140 dat$FR2.IMGT.a.c +
|
|
|
141 dat$FR2.IMGT.c.a +
|
|
|
142 dat$FR2.IMGT.a.t +
|
|
|
143 dat$FR2.IMGT.t.a +
|
|
|
144 dat$FR2.IMGT.g.c +
|
|
|
145 dat$FR2.IMGT.c.g +
|
|
|
146 dat$FR2.IMGT.g.t +
|
|
|
147 dat$FR2.IMGT.t.g +
|
|
|
148 dat$CDR2.IMGT.a.c +
|
|
|
149 dat$CDR2.IMGT.c.a +
|
|
|
150 dat$CDR2.IMGT.a.t +
|
|
|
151 dat$CDR2.IMGT.t.a +
|
|
|
152 dat$CDR2.IMGT.g.c +
|
|
|
153 dat$CDR2.IMGT.c.g +
|
|
|
154 dat$CDR2.IMGT.g.t +
|
|
|
155 dat$CDR2.IMGT.t.g +
|
|
|
156 dat$FR3.IMGT.a.c +
|
|
|
157 dat$FR3.IMGT.c.a +
|
|
|
158 dat$FR3.IMGT.a.t +
|
|
|
159 dat$FR3.IMGT.t.a +
|
|
|
160 dat$FR3.IMGT.g.c +
|
|
|
161 dat$FR3.IMGT.c.g +
|
|
|
162 dat$FR3.IMGT.g.t +
|
|
|
163 dat$FR3.IMGT.t.g
|
|
0
|
164
|
|
|
165
|
|
7
|
166 dat$transitionMutationsAtGC = dat$CDR1.IMGT.g.a +
|
|
0
|
167 dat$CDR1.IMGT.c.t +
|
|
|
168 dat$FR2.IMGT.g.a +
|
|
|
169 dat$FR2.IMGT.c.t +
|
|
|
170 dat$CDR2.IMGT.g.a +
|
|
|
171 dat$CDR2.IMGT.c.t +
|
|
|
172 dat$FR3.IMGT.g.a +
|
|
4
|
173 dat$FR3.IMGT.c.t
|
|
0
|
174
|
|
7
|
175 dat$totalMutationsAtGC = dat$CDR1.IMGT.g.a +
|
|
0
|
176 dat$CDR1.IMGT.c.t +
|
|
|
177 dat$CDR1.IMGT.c.a +
|
|
|
178 dat$CDR1.IMGT.g.c +
|
|
|
179 dat$CDR1.IMGT.c.g +
|
|
|
180 dat$CDR1.IMGT.g.t +
|
|
|
181 dat$FR2.IMGT.g.a +
|
|
|
182 dat$FR2.IMGT.c.t +
|
|
|
183 dat$FR2.IMGT.c.a +
|
|
|
184 dat$FR2.IMGT.g.c +
|
|
|
185 dat$FR2.IMGT.c.g +
|
|
|
186 dat$FR2.IMGT.g.t +
|
|
|
187 dat$CDR2.IMGT.g.a +
|
|
|
188 dat$CDR2.IMGT.c.t +
|
|
|
189 dat$CDR2.IMGT.c.a +
|
|
|
190 dat$CDR2.IMGT.g.c +
|
|
|
191 dat$CDR2.IMGT.c.g +
|
|
|
192 dat$CDR2.IMGT.g.t +
|
|
|
193 dat$FR3.IMGT.g.a +
|
|
|
194 dat$FR3.IMGT.c.t +
|
|
|
195 dat$FR3.IMGT.c.a +
|
|
|
196 dat$FR3.IMGT.g.c +
|
|
|
197 dat$FR3.IMGT.c.g +
|
|
4
|
198 dat$FR3.IMGT.g.t
|
|
0
|
199
|
|
|
200
|
|
|
201
|
|
4
|
202 setwd(outputdir)
|
|
|
203
|
|
|
204 matrx = matrix(data = 0, ncol=((length(genes) + 1) * 3),nrow=5)
|
|
|
205 for(i in 1:length(genes)){
|
|
|
206 gene = genes[i]
|
|
|
207 tmp = dat[grepl(paste(".*", gene, ".*", sep=""), dat$best_match),]
|
|
|
208 if(gene == "."){
|
|
|
209 tmp = dat
|
|
|
210 }
|
|
|
211 if(length(tmp) == 0){
|
|
|
212 cat("0", file=paste(gene, "_value.txt" ,sep=""))
|
|
|
213 next
|
|
|
214 }
|
|
|
215 j = i - 1
|
|
|
216 x = (j * 3) + 1
|
|
|
217 y = (j * 3) + 2
|
|
|
218 z = (j * 3) + 3
|
|
|
219 matrx[1,x] = sum(tmp$VRegionMutations)
|
|
|
220 matrx[1,y] = sum(tmp$VRegionNucleotides)
|
|
|
221 matrx[1,z] = round(matrx[1,x] / matrx[1,y] * 100, digits=1)
|
|
|
222 matrx[2,x] = sum(tmp$transitionMutations)
|
|
|
223 matrx[2,y] = sum(tmp$VRegionMutations)
|
|
|
224 matrx[2,z] = round(matrx[2,x] / matrx[2,y] * 100, digits=1)
|
|
|
225 matrx[3,x] = sum(tmp$transversionMutations)
|
|
|
226 matrx[3,y] = sum(tmp$VRegionMutations)
|
|
|
227 matrx[3,z] = round(matrx[3,x] / matrx[3,y] * 100, digits=1)
|
|
|
228 matrx[4,x] = sum(tmp$transitionMutationsAtGC)
|
|
|
229 matrx[4,y] = sum(tmp$totalMutationsAtGC)
|
|
|
230 matrx[4,z] = round(matrx[4,x] / matrx[4,y] * 100, digits=1)
|
|
|
231 matrx[5,x] = sum(tmp$totalMutationsAtGC)
|
|
|
232 matrx[5,y] = sum(tmp$VRegionMutations)
|
|
|
233 matrx[5,z] = round(matrx[5,x] / matrx[5,y] * 100, digits=1)
|
|
|
234
|
|
|
235 transitionTable = data.frame(A=1:4,C=1:4,G=1:4,T=1:4)
|
|
|
236 row.names(transitionTable) = c("A", "C", "G", "T")
|
|
|
237 transitionTable["A","A"] = NA
|
|
|
238 transitionTable["C","C"] = NA
|
|
|
239 transitionTable["G","G"] = NA
|
|
|
240 transitionTable["T","T"] = NA
|
|
|
241 nts = c("a", "c", "g", "t")
|
|
|
242
|
|
|
243
|
|
|
244 for(nt1 in nts){
|
|
|
245 for(nt2 in nts){
|
|
|
246 if(nt1 == nt2){
|
|
|
247 next
|
|
|
248 }
|
|
|
249 NT1 = LETTERS[letters == nt1]
|
|
|
250 NT2 = LETTERS[letters == nt2]
|
|
7
|
251 FR1 = 0 #paste("FR1.IMGT.", nt1, ".", nt2, sep="")
|
|
4
|
252 CDR1 = paste("CDR1.IMGT.", nt1, ".", nt2, sep="")
|
|
|
253 FR2 = paste("FR2.IMGT.", nt1, ".", nt2, sep="")
|
|
|
254 CDR2 = paste("CDR2.IMGT.", nt1, ".", nt2, sep="")
|
|
|
255 FR3 = paste("FR3.IMGT.", nt1, ".", nt2, sep="")
|
|
7
|
256 transitionTable[NT1,NT2] = sum( tmp[,CDR1] +
|
|
|
257 tmp[,FR2] +
|
|
|
258 tmp[,CDR2] +
|
|
|
259 tmp[,FR3])
|
|
4
|
260 }
|
|
|
261 }
|
|
|
262 write.table(x=transitionTable, file=paste("transitions_", gene ,".txt", sep=""), sep=",",quote=F,row.names=T,col.names=NA)
|
|
|
263 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)
|
|
|
264
|
|
|
265 cat(matrx[1,x], file=paste(gene, "_value.txt" ,sep=""))
|
|
|
266 cat(length(tmp$Sequence.ID), file=paste(gene, "_n.txt" ,sep=""))
|
|
|
267 }
|
|
|
268
|
|
|
269 #again for all of the data
|
|
|
270 tmp = dat
|
|
|
271 j = i
|
|
|
272 x = (j * 3) + 1
|
|
|
273 y = (j * 3) + 2
|
|
|
274 z = (j * 3) + 3
|
|
|
275 matrx[1,x] = sum(tmp$VRegionMutations)
|
|
|
276 matrx[1,y] = sum(tmp$VRegionNucleotides)
|
|
|
277 matrx[1,z] = round(matrx[1,x] / matrx[1,y] * 100, digits=1)
|
|
|
278 matrx[2,x] = sum(tmp$transitionMutations)
|
|
|
279 matrx[2,y] = sum(tmp$VRegionMutations)
|
|
|
280 matrx[2,z] = round(matrx[2,x] / matrx[2,y] * 100, digits=1)
|
|
|
281 matrx[3,x] = sum(tmp$transversionMutations)
|
|
|
282 matrx[3,y] = sum(tmp$VRegionMutations)
|
|
|
283 matrx[3,z] = round(matrx[3,x] / matrx[3,y] * 100, digits=1)
|
|
|
284 matrx[4,x] = sum(tmp$transitionMutationsAtGC)
|
|
|
285 matrx[4,y] = sum(tmp$totalMutationsAtGC)
|
|
|
286 matrx[4,z] = round(matrx[4,x] / matrx[4,y] * 100, digits=1)
|
|
|
287 matrx[5,x] = sum(tmp$totalMutationsAtGC)
|
|
|
288 matrx[5,y] = sum(tmp$VRegionMutations)
|
|
|
289 matrx[5,z] = round(matrx[5,x] / matrx[5,y] * 100, digits=1)
|
|
0
|
290
|
|
|
291 transitionTable = data.frame(A=1:4,C=1:4,G=1:4,T=1:4)
|
|
|
292 row.names(transitionTable) = c("A", "C", "G", "T")
|
|
|
293 transitionTable["A","A"] = NA
|
|
|
294 transitionTable["C","C"] = NA
|
|
|
295 transitionTable["G","G"] = NA
|
|
|
296 transitionTable["T","T"] = NA
|
|
|
297 nts = c("a", "c", "g", "t")
|
|
|
298
|
|
|
299
|
|
|
300 for(nt1 in nts){
|
|
4
|
301 for(nt2 in nts){
|
|
|
302 if(nt1 == nt2){
|
|
|
303 next
|
|
|
304 }
|
|
|
305 NT1 = LETTERS[letters == nt1]
|
|
|
306 NT2 = LETTERS[letters == nt2]
|
|
|
307 FR1 = paste("FR1.IMGT.", nt1, ".", nt2, sep="")
|
|
|
308 CDR1 = paste("CDR1.IMGT.", nt1, ".", nt2, sep="")
|
|
|
309 FR2 = paste("FR2.IMGT.", nt1, ".", nt2, sep="")
|
|
|
310 CDR2 = paste("CDR2.IMGT.", nt1, ".", nt2, sep="")
|
|
|
311 FR3 = paste("FR3.IMGT.", nt1, ".", nt2, sep="")
|
|
7
|
312 transitionTable[NT1,NT2] = sum( tmp[,CDR1] +
|
|
|
313 tmp[,FR2] +
|
|
|
314 tmp[,CDR2] +
|
|
|
315 tmp[,FR3])
|
|
4
|
316 }
|
|
|
317 }
|
|
|
318 write.table(x=transitionTable, file="transitions.txt", sep=",",quote=F,row.names=T,col.names=NA)
|
|
|
319 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)
|
|
|
320 cat(matrx[1,x], file="total_value.txt")
|
|
|
321 cat(length(tmp$Sequence.ID), file="total_n.txt")
|
|
|
322
|
|
|
323
|
|
|
324
|
|
|
325 result = data.frame(matrx)
|
|
|
326 row.names(result) = c("Number of Mutations (%)", "Transition (%)", "Transversions (%)", "Transitions at G C (%)", "Targeting of C.G (%)")
|
|
|
327
|
|
|
328 write.table(x=result, file="mutations.txt", sep=",",quote=F,row.names=T,col.names=F)
|
|
|
329
|
|
|
330
|
|
|
331 if (!("ggplot2" %in% rownames(installed.packages()))) {
|
|
|
332 install.packages("ggplot2", repos="http://cran.xl-mirror.nl/")
|
|
|
333 }
|
|
|
334 library(ggplot2)
|
|
|
335
|
|
|
336 genesForPlot = gsub("[0-9]", "", dat$best_match)
|
|
|
337 genesForPlot = data.frame(table(genesForPlot))
|
|
|
338 colnames(genesForPlot) = c("Gene","Freq")
|
|
|
339 genesForPlot$label = paste(genesForPlot$Gene, "-", genesForPlot$Freq)
|
|
|
340
|
|
|
341 pc = ggplot(genesForPlot, aes(x = factor(1), y=Freq, fill=label))
|
|
|
342 pc = pc + geom_bar(width = 1, stat = "identity")
|
|
|
343 pc = pc + coord_polar(theta="y")
|
|
|
344 pc = pc + xlab(" ") + ylab(" ") + ggtitle(paste("IgA", "( n =", sum(genesForPlot$Freq), ")"))
|
|
|
345
|
|
|
346 png(filename="all.png")
|
|
|
347 pc
|
|
|
348 dev.off()
|
|
|
349
|
|
|
350
|
|
|
351 #blegh
|
|
|
352 genesForPlot = dat[grepl("ca", dat$best_match),]$best_match
|
|
|
353 if(length(genesForPlot) > 0){
|
|
|
354 genesForPlot = data.frame(table(genesForPlot))
|
|
|
355 colnames(genesForPlot) = c("Gene","Freq")
|
|
|
356 genesForPlot$label = paste(genesForPlot$Gene, "-", genesForPlot$Freq)
|
|
|
357
|
|
|
358 pc = ggplot(genesForPlot, aes(x = factor(1), y=Freq, fill=label))
|
|
|
359 pc = pc + geom_bar(width = 1, stat = "identity")
|
|
|
360 pc = pc + coord_polar(theta="y")
|
|
|
361 pc = pc + xlab(" ") + ylab(" ") + ggtitle(paste("IgA", "( n =", sum(genesForPlot$Freq), ")"))
|
|
|
362
|
|
|
363
|
|
|
364 png(filename="ca.png")
|
|
|
365 print(pc)
|
|
|
366 dev.off()
|
|
0
|
367 }
|
|
|
368
|
|
4
|
369 genesForPlot = dat[grepl("cg", dat$best_match),]$best_match
|
|
|
370 if(length(genesForPlot) > 0){
|
|
|
371 genesForPlot = data.frame(table(genesForPlot))
|
|
|
372 colnames(genesForPlot) = c("Gene","Freq")
|
|
|
373 genesForPlot$label = paste(genesForPlot$Gene, "-", genesForPlot$Freq)
|
|
|
374
|
|
|
375 pc = ggplot(genesForPlot, aes(x = factor(1), y=Freq, fill=label))
|
|
|
376 pc = pc + geom_bar(width = 1, stat = "identity")
|
|
|
377 pc = pc + coord_polar(theta="y")
|
|
|
378 pc = pc + xlab(" ") + ylab(" ") + ggtitle(paste("IgG", "( n =", sum(genesForPlot$Freq), ")"))
|
|
0
|
379
|
|
|
380
|
|
4
|
381 png(filename="cg.png")
|
|
|
382 print(pc)
|
|
|
383 dev.off()
|
|
|
384 }
|
