|
0
|
1 args <- commandArgs(trailingOnly = TRUE)
|
|
|
2
|
|
|
3 input = args[1]
|
|
|
4 summaryinput = args[2]
|
|
|
5 outputdir = args[3]
|
|
|
6 setwd(outputdir)
|
|
|
7
|
|
|
8 #dat = read.table("NWK276_MID6_25NT/8_V-REGION-nt-mutation-statistics_NWK276_MID6_25NT_051113.txt", header=T, sep="\t", fill=T, stringsAsFactors=F)
|
|
|
9 dat = read.table(input, header=T, sep="\t", fill=T, stringsAsFactors=F)
|
|
|
10
|
|
|
11 datSum = read.table(summaryinput, header=T, sep="\t", fill=T, stringsAsFactors=F)
|
|
3
|
12 datSum = datSum[,c("Sequence.ID","J.GENE.and.allele", "AA.JUNCTION")]
|
|
0
|
13
|
|
|
14 dat = merge(dat, datSum, by="Sequence.ID", all.x=T)
|
|
|
15
|
|
|
16 if(length(dat$Sequence.ID) == 0){
|
|
|
17 setwd(outputdir)
|
|
|
18 result = data.frame(x = rep(0, 5), y = rep(0, 5), z = rep(NA, 5))
|
|
|
19 row.names(result) = c("Number of Mutations (%)", "Transition (%)", "Transversions (%)", "Transitions at G C (%)", "Targeting of C G (%)")
|
|
|
20 write.table(x=result, file="mutations.txt", sep=",",quote=F,row.names=T,col.names=F)
|
|
|
21 transitionTable = data.frame(A=rep(0, 4),C=rep(0, 4),G=rep(0, 4),T=rep(0, 4))
|
|
|
22 row.names(transitionTable) = c("A", "C", "G", "T")
|
|
|
23 transitionTable["A","A"] = NA
|
|
|
24 transitionTable["C","C"] = NA
|
|
|
25 transitionTable["G","G"] = NA
|
|
|
26 transitionTable["T","T"] = NA
|
|
|
27 write.table(x=transitionTable, file="transitions.txt", sep=",",quote=F,row.names=T,col.names=NA)
|
|
|
28 cat("0", file="n.txt")
|
|
|
29 stop("No data")
|
|
|
30 }
|
|
|
31
|
|
|
32 #print(paste("After filtering on 'productive' and deduplicating based on V and AA JUNCTION there are", length(dat$Sequence.ID), "sequences left"))
|
|
|
33
|
|
|
34 result = data.frame(x = 1:5, y = 1:5, z = 1:5)
|
|
|
35 row.names(result) = c("Number of Mutations (%)", "Transition (%)", "Transversions (%)", "Transitions at G C (%)", "Targeting of C G (%)")
|
|
|
36
|
|
|
37 cleanup_columns = c("FR1.IMGT.c.a",
|
|
|
38 "FR2.IMGT.g.t",
|
|
|
39 "CDR1.IMGT.Nb.of.nucleotides",
|
|
|
40 "CDR2.IMGT.t.a",
|
|
|
41 "FR1.IMGT.c.g",
|
|
|
42 "CDR1.IMGT.c.t",
|
|
|
43 "FR2.IMGT.a.c",
|
|
|
44 "FR2.IMGT.Nb.of.mutations",
|
|
|
45 "FR2.IMGT.g.c",
|
|
|
46 "FR2.IMGT.a.g",
|
|
|
47 "FR3.IMGT.t.a",
|
|
|
48 "FR3.IMGT.t.c",
|
|
|
49 "FR2.IMGT.g.a",
|
|
|
50 "FR3.IMGT.c.g",
|
|
|
51 "FR1.IMGT.Nb.of.mutations",
|
|
|
52 "CDR1.IMGT.g.a",
|
|
|
53 "CDR1.IMGT.t.g",
|
|
|
54 "CDR1.IMGT.g.c",
|
|
|
55 "CDR2.IMGT.Nb.of.nucleotides",
|
|
|
56 "FR2.IMGT.a.t",
|
|
|
57 "CDR1.IMGT.Nb.of.mutations",
|
|
|
58 "CDR1.IMGT.a.g",
|
|
|
59 "FR3.IMGT.a.c",
|
|
|
60 "FR1.IMGT.g.a",
|
|
|
61 "FR3.IMGT.a.g",
|
|
|
62 "FR1.IMGT.a.t",
|
|
|
63 "CDR2.IMGT.a.g",
|
|
|
64 "CDR2.IMGT.Nb.of.mutations",
|
|
|
65 "CDR2.IMGT.g.t",
|
|
|
66 "CDR2.IMGT.a.c",
|
|
|
67 "CDR1.IMGT.t.c",
|
|
|
68 "FR3.IMGT.g.c",
|
|
|
69 "FR1.IMGT.g.t",
|
|
|
70 "FR3.IMGT.g.t",
|
|
|
71 "CDR1.IMGT.a.t",
|
|
|
72 "FR1.IMGT.a.g",
|
|
|
73 "FR3.IMGT.a.t",
|
|
|
74 "FR3.IMGT.Nb.of.nucleotides",
|
|
|
75 "FR2.IMGT.t.c",
|
|
|
76 "CDR2.IMGT.g.a",
|
|
|
77 "FR2.IMGT.t.a",
|
|
|
78 "CDR1.IMGT.t.a",
|
|
|
79 "FR2.IMGT.t.g",
|
|
|
80 "FR3.IMGT.t.g",
|
|
|
81 "FR2.IMGT.Nb.of.nucleotides",
|
|
|
82 "FR1.IMGT.t.a",
|
|
|
83 "FR1.IMGT.t.g",
|
|
|
84 "FR3.IMGT.c.t",
|
|
|
85 "FR1.IMGT.t.c",
|
|
|
86 "CDR2.IMGT.a.t",
|
|
|
87 "FR2.IMGT.c.t",
|
|
|
88 "CDR1.IMGT.g.t",
|
|
|
89 "CDR2.IMGT.t.g",
|
|
|
90 "FR1.IMGT.Nb.of.nucleotides",
|
|
|
91 "CDR1.IMGT.c.g",
|
|
|
92 "CDR2.IMGT.t.c",
|
|
|
93 "FR3.IMGT.g.a",
|
|
|
94 "CDR1.IMGT.a.c",
|
|
|
95 "FR2.IMGT.c.a",
|
|
|
96 "FR3.IMGT.Nb.of.mutations",
|
|
|
97 "FR2.IMGT.c.g",
|
|
|
98 "CDR2.IMGT.g.c",
|
|
|
99 "FR1.IMGT.g.c",
|
|
|
100 "CDR2.IMGT.c.t",
|
|
|
101 "FR3.IMGT.c.a",
|
|
|
102 "CDR1.IMGT.c.a",
|
|
|
103 "CDR2.IMGT.c.g",
|
|
|
104 "CDR2.IMGT.c.a",
|
|
|
105 "FR1.IMGT.c.t")
|
|
|
106
|
|
|
107 for(col in cleanup_columns){
|
|
|
108 dat[,col] = gsub("\\(.*\\)", "", dat[,col])
|
|
|
109 #dat[dat[,col] == "",] = "0"
|
|
|
110 dat[,col] = as.numeric(dat[,col])
|
|
|
111 dat[is.na(dat[,col]),] = 0
|
|
|
112 }
|
|
|
113
|
|
3
|
114 dat$VGene = gsub("^Homsap ", "", dat$V.GENE.and.allele)
|
|
|
115 dat$VGene = gsub("[*].*", "", dat$VGene)
|
|
|
116 dat$JGene = gsub("^Homsap ", "", dat$J.GENE.and.allele)
|
|
|
117 dat$JGene = gsub("[*].*", "", dat$JGene)
|
|
|
118
|
|
|
119 dat$past = paste(dat$AA.JUNCTION, dat$VGene, dat$JGene, (dat$FR1.IMGT.Nb.of.mutations + dat$CDR1.IMGT.Nb.of.mutations + dat$FR2.IMGT.Nb.of.mutations + dat$CDR2.IMGT.Nb.of.mutations + dat$FR3.IMGT.Nb.of.mutations))
|
|
|
120
|
|
|
121 dat = dat[!duplicated(dat$past), ]
|
|
|
122
|
|
0
|
123 VRegionMutations = sum(dat$FR1.IMGT.Nb.of.mutations +
|
|
|
124 dat$CDR1.IMGT.Nb.of.mutations +
|
|
|
125 dat$FR2.IMGT.Nb.of.mutations +
|
|
|
126 dat$CDR2.IMGT.Nb.of.mutations +
|
|
|
127 dat$FR3.IMGT.Nb.of.mutations)
|
|
|
128
|
|
|
129 VRegionNucleotides = sum( dat$FR1.IMGT.Nb.of.nucleotides +
|
|
|
130 dat$CDR1.IMGT.Nb.of.nucleotides +
|
|
|
131 dat$FR2.IMGT.Nb.of.nucleotides +
|
|
|
132 dat$CDR2.IMGT.Nb.of.nucleotides +
|
|
|
133 dat$FR3.IMGT.Nb.of.nucleotides)
|
|
|
134
|
|
|
135 result[1,"x"] = VRegionMutations
|
|
|
136 result[1,"y"] = VRegionNucleotides
|
|
|
137 result[1,"z"] = round(result[1,"x"] / result[1,"y"] * 100, digits=1)
|
|
|
138
|
|
|
139 transitionMutations = sum(dat$FR1.IMGT.a.g +
|
|
|
140 dat$FR1.IMGT.g.a +
|
|
|
141 dat$FR1.IMGT.c.t +
|
|
|
142 dat$FR1.IMGT.t.c +
|
|
|
143 dat$CDR1.IMGT.a.g +
|
|
|
144 dat$CDR1.IMGT.g.a +
|
|
|
145 dat$CDR1.IMGT.c.t +
|
|
|
146 dat$CDR1.IMGT.t.c +
|
|
|
147 dat$FR2.IMGT.a.g +
|
|
|
148 dat$FR2.IMGT.g.a +
|
|
|
149 dat$FR2.IMGT.c.t +
|
|
|
150 dat$FR2.IMGT.t.c +
|
|
|
151 dat$CDR2.IMGT.a.g +
|
|
|
152 dat$CDR2.IMGT.g.a +
|
|
|
153 dat$CDR2.IMGT.c.t +
|
|
|
154 dat$CDR2.IMGT.t.c +
|
|
|
155 dat$FR3.IMGT.a.g +
|
|
|
156 dat$FR3.IMGT.g.a +
|
|
|
157 dat$FR3.IMGT.c.t +
|
|
|
158 dat$FR3.IMGT.t.c)
|
|
|
159
|
|
|
160 result[2,"x"] = transitionMutations
|
|
|
161 result[2,"y"] = VRegionMutations
|
|
|
162 result[2,"z"] = round(result[2,"x"] / result[2,"y"] * 100, digits=1)
|
|
|
163
|
|
|
164 transversionMutations = sum( dat$FR1.IMGT.a.c +
|
|
|
165 dat$FR1.IMGT.c.a +
|
|
|
166 dat$FR1.IMGT.a.t +
|
|
|
167 dat$FR1.IMGT.t.a +
|
|
|
168 dat$FR1.IMGT.g.c +
|
|
|
169 dat$FR1.IMGT.c.g +
|
|
|
170 dat$FR1.IMGT.g.t +
|
|
|
171 dat$FR1.IMGT.t.g +
|
|
|
172 dat$CDR1.IMGT.a.c +
|
|
|
173 dat$CDR1.IMGT.c.a +
|
|
|
174 dat$CDR1.IMGT.a.t +
|
|
|
175 dat$CDR1.IMGT.t.a +
|
|
|
176 dat$CDR1.IMGT.g.c +
|
|
|
177 dat$CDR1.IMGT.c.g +
|
|
|
178 dat$CDR1.IMGT.g.t +
|
|
|
179 dat$CDR1.IMGT.t.g +
|
|
|
180 dat$FR2.IMGT.a.c +
|
|
|
181 dat$FR2.IMGT.c.a +
|
|
|
182 dat$FR2.IMGT.a.t +
|
|
|
183 dat$FR2.IMGT.t.a +
|
|
|
184 dat$FR2.IMGT.g.c +
|
|
|
185 dat$FR2.IMGT.c.g +
|
|
|
186 dat$FR2.IMGT.g.t +
|
|
|
187 dat$FR2.IMGT.t.g +
|
|
|
188 dat$CDR2.IMGT.a.c +
|
|
|
189 dat$CDR2.IMGT.c.a +
|
|
|
190 dat$CDR2.IMGT.a.t +
|
|
|
191 dat$CDR2.IMGT.t.a +
|
|
|
192 dat$CDR2.IMGT.g.c +
|
|
|
193 dat$CDR2.IMGT.c.g +
|
|
|
194 dat$CDR2.IMGT.g.t +
|
|
|
195 dat$CDR2.IMGT.t.g +
|
|
|
196 dat$FR3.IMGT.a.c +
|
|
|
197 dat$FR3.IMGT.c.a +
|
|
|
198 dat$FR3.IMGT.a.t +
|
|
|
199 dat$FR3.IMGT.t.a +
|
|
|
200 dat$FR3.IMGT.g.c +
|
|
|
201 dat$FR3.IMGT.c.g +
|
|
|
202 dat$FR3.IMGT.g.t +
|
|
|
203 dat$FR3.IMGT.t.g)
|
|
|
204
|
|
|
205 result[3,"x"] = transversionMutations
|
|
|
206 result[3,"y"] = VRegionMutations
|
|
|
207 result[3,"z"] = round(result[3,"x"] / result[3,"y"] * 100, digits=1)
|
|
|
208
|
|
|
209 transitionMutationsAtGC = sum(dat$FR1.IMGT.g.a +
|
|
|
210 dat$FR1.IMGT.c.t +
|
|
|
211 dat$CDR1.IMGT.g.a +
|
|
|
212 dat$CDR1.IMGT.c.t +
|
|
|
213 dat$FR2.IMGT.g.a +
|
|
|
214 dat$FR2.IMGT.c.t +
|
|
|
215 dat$CDR2.IMGT.g.a +
|
|
|
216 dat$CDR2.IMGT.c.t +
|
|
|
217 dat$FR3.IMGT.g.a +
|
|
|
218 dat$FR3.IMGT.c.t)
|
|
|
219
|
|
|
220 totalMutationsAtGC = sum(dat$FR1.IMGT.g.a +
|
|
|
221 dat$FR1.IMGT.c.t +
|
|
|
222 dat$FR1.IMGT.c.a +
|
|
|
223 dat$FR1.IMGT.g.c +
|
|
|
224 dat$FR1.IMGT.c.g +
|
|
|
225 dat$FR1.IMGT.g.t +
|
|
|
226 dat$CDR1.IMGT.g.a +
|
|
|
227 dat$CDR1.IMGT.c.t +
|
|
|
228 dat$CDR1.IMGT.c.a +
|
|
|
229 dat$CDR1.IMGT.g.c +
|
|
|
230 dat$CDR1.IMGT.c.g +
|
|
|
231 dat$CDR1.IMGT.g.t +
|
|
|
232 dat$FR2.IMGT.g.a +
|
|
|
233 dat$FR2.IMGT.c.t +
|
|
|
234 dat$FR2.IMGT.c.a +
|
|
|
235 dat$FR2.IMGT.g.c +
|
|
|
236 dat$FR2.IMGT.c.g +
|
|
|
237 dat$FR2.IMGT.g.t +
|
|
|
238 dat$CDR2.IMGT.g.a +
|
|
|
239 dat$CDR2.IMGT.c.t +
|
|
|
240 dat$CDR2.IMGT.c.a +
|
|
|
241 dat$CDR2.IMGT.g.c +
|
|
|
242 dat$CDR2.IMGT.c.g +
|
|
|
243 dat$CDR2.IMGT.g.t +
|
|
|
244 dat$FR3.IMGT.g.a +
|
|
|
245 dat$FR3.IMGT.c.t +
|
|
|
246 dat$FR3.IMGT.c.a +
|
|
|
247 dat$FR3.IMGT.g.c +
|
|
|
248 dat$FR3.IMGT.c.g +
|
|
|
249 dat$FR3.IMGT.g.t)
|
|
|
250
|
|
|
251 result[4,"x"] = transitionMutationsAtGC
|
|
|
252 result[4,"y"] = totalMutationsAtGC
|
|
|
253 result[4,"z"] = round(result[4,"x"] / result[4,"y"] * 100, digits=1)
|
|
|
254
|
|
|
255 result[5,"x"] = totalMutationsAtGC
|
|
|
256 result[5,"y"] = VRegionMutations
|
|
|
257 result[5,"z"] = round(result[5,"x"] / result[5,"y"] * 100, digits=1)
|
|
|
258
|
|
|
259
|
|
|
260 #transitiontable
|
|
|
261
|
|
|
262 transitionTable = data.frame(A=1:4,C=1:4,G=1:4,T=1:4)
|
|
|
263 row.names(transitionTable) = c("A", "C", "G", "T")
|
|
|
264 transitionTable["A","A"] = NA
|
|
|
265 transitionTable["C","C"] = NA
|
|
|
266 transitionTable["G","G"] = NA
|
|
|
267 transitionTable["T","T"] = NA
|
|
|
268 nts = c("a", "c", "g", "t")
|
|
|
269
|
|
|
270
|
|
|
271 for(nt1 in nts){
|
|
|
272 for(nt2 in nts){
|
|
|
273 if(nt1 == nt2){
|
|
|
274 next
|
|
|
275 }
|
|
|
276 NT1 = LETTERS[letters == nt1]
|
|
|
277 NT2 = LETTERS[letters == nt2]
|
|
|
278 FR1 = paste("FR1.IMGT.", nt1, ".", nt2, sep="")
|
|
|
279 CDR1 = paste("CDR1.IMGT.", nt1, ".", nt2, sep="")
|
|
|
280 FR2 = paste("FR2.IMGT.", nt1, ".", nt2, sep="")
|
|
|
281 CDR2 = paste("CDR2.IMGT.", nt1, ".", nt2, sep="")
|
|
|
282 FR3 = paste("FR3.IMGT.", nt1, ".", nt2, sep="")
|
|
|
283 transitionTable[NT1,NT2] = sum( dat[,FR1] +
|
|
|
284 dat[,CDR1] +
|
|
|
285 dat[,FR2] +
|
|
|
286 dat[,CDR2] +
|
|
|
287 dat[,FR3])
|
|
|
288 }
|
|
|
289 }
|
|
|
290
|
|
|
291 setwd(outputdir)
|
|
|
292 write.table(x=result, file="mutations.txt", sep=",",quote=F,row.names=T,col.names=F)
|
|
|
293 write.table(x=transitionTable, file="transitions.txt", sep=",",quote=F,row.names=T,col.names=NA)
|
|
|
294 cat(length(dat$Sequence.ID), file="n.txt")
|
|
|
295
|
|
|
296
|
|
2
|
297 weblogo = dat[,c("Sequence.ID", "VGene")]
|
|
|
298 weblogo$VGene = gsub("\\*.*", "", weblogo$VGene)
|
|
0
|
299
|
|
2
|
300 rs12 = read.table("HS12RSS.txt", sep="\t", header=TRUE)
|
|
|
301 rs23 = read.table("HS23RSS.txt", sep="\t", header=TRUE)
|
|
0
|
302
|
|
2
|
303 result12 = merge(weblogo, rs12, by.x="VGene", by.y="Gene")
|
|
|
304 result23 = merge(weblogo, rs23, by.x="VGene", by.y="Gene")
|
|
0
|
305
|
|
2
|
306 write.table(x=result12$Sequence, file="weblogo_in_rs12.txt", sep=",",quote=F,row.names=F,col.names=F)
|
|
|
307 write.table(x=result23$Sequence, file="weblogo_in_rs23.txt", sep=",",quote=F,row.names=F,col.names=F)
|
|
0
|
308
|
