Mercurial > repos > davidvanzessen > mutation_analysis

comparison mutation_analysis.r @ 53:7290a88ea202 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author davidvanzessen
date Mon, 29 Feb 2016 10:49:39 -0500
parents 5c6b9e99d576
children 3636d5aaa127
comparison
equal deleted inserted replaced
52:d3542f87a304 53:7290a88ea202
168 168
169 write.table(dat[,mutation.sum.columns], "mutation_by_id.txt", sep="\t",quote=F,row.names=F,col.names=T) 169 write.table(dat[,mutation.sum.columns], "mutation_by_id.txt", sep="\t",quote=F,row.names=F,col.names=T)
170 170
171 setwd(outputdir) 171 setwd(outputdir)
172 172
173 nts = c("a", "c", "g", "t") 173
174 zeros=rep(0, 4) 174 calculate_result = function(i, gene, dat, matrx, f, fname, name){
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),] 175 tmp = dat[grepl(paste(".*", gene, ".*", sep=""), dat$best_match),]
179 if(gene == "."){ 176
180 tmp = dat
181 }
182 j = i - 1 177 j = i - 1
183 x = (j * 3) + 1 178 x = (j * 3) + 1
184 y = (j * 3) + 2 179 y = (j * 3) + 2
185 z = (j * 3) + 3 180 z = (j * 3) + 3
186 matrx[1,x] = sum(tmp$VRegionMutations) 181
187 matrx[1,y] = sum(tmp$VRegionNucleotides) 182 if(nrow(tmp) > 0){
188 matrx[1,z] = round(matrx[1,x] / matrx[1,y] * 100, digits=1) 183
189 184 matrx[1,x] = round(f(tmp$VRegionMutations, na.rm=T), digits=1)
190 matrx[2,x] = sum(tmp$transitionMutations) 185 matrx[1,y] = round(f(tmp$VRegionNucleotides, na.rm=T), digits=1)
191 matrx[2,y] = sum(tmp$VRegionMutations) 186 matrx[1,z] = round(matrx[1,x] / matrx[1,y] * 100, digits=1)
192 matrx[2,z] = round(matrx[2,x] / matrx[2,y] * 100, digits=1) 187
193 188 matrx[2,x] = round(f(tmp$transitionMutations, na.rm=T), digits=1)
194 matrx[3,x] = sum(tmp$transversionMutations) 189 matrx[2,y] = round(f(tmp$VRegionMutations, na.rm=T), digits=1)
195 matrx[3,y] = sum(tmp$VRegionMutations) 190 matrx[2,z] = round(matrx[2,x] / matrx[2,y] * 100, digits=1)
196 matrx[3,z] = round(matrx[3,x] / matrx[3,y] * 100, digits=1) 191
197 192 matrx[3,x] = round(f(tmp$transversionMutations, na.rm=T), digits=1)
198 matrx[4,x] = sum(tmp$transitionMutationsAtGC) 193 matrx[3,y] = round(f(tmp$VRegionMutations, na.rm=T), digits=1)
199 matrx[4,y] = sum(tmp$totalMutationsAtGC) 194 matrx[3,z] = round(matrx[3,x] / matrx[3,y] * 100, digits=1)
200 matrx[4,z] = round(matrx[4,x] / matrx[4,y] * 100, digits=1) 195
201 196 matrx[4,x] = round(f(tmp$transitionMutationsAtGC, na.rm=T), digits=1)
202 matrx[5,x] = sum(tmp$totalMutationsAtGC) 197 matrx[4,y] = round(f(tmp$totalMutationsAtGC, na.rm=T), digits=1)
203 matrx[5,y] = sum(tmp$VRegionMutations) 198 matrx[4,z] = round(matrx[4,x] / matrx[4,y] * 100, digits=1)
204 matrx[5,z] = round(matrx[5,x] / matrx[5,y] * 100, digits=1) 199
205 200 matrx[5,x] = round(f(tmp$totalMutationsAtGC, na.rm=T), digits=1)
206 matrx[6,x] = sum(tmp$transitionMutationsAtAT) 201 matrx[5,y] = round(f(tmp$VRegionMutations, na.rm=T), digits=1)
207 matrx[6,y] = sum(tmp$totalMutationsAtAT) 202 matrx[5,z] = round(matrx[5,x] / matrx[5,y] * 100, digits=1)
208 matrx[6,z] = round(matrx[6,x] / matrx[6,y] * 100, digits=1) 203
209 204 matrx[6,x] = round(f(tmp$transitionMutationsAtAT, na.rm=T), digits=1)
210 matrx[7,x] = sum(tmp$totalMutationsAtAT) 205 matrx[6,y] = round(f(tmp$totalMutationsAtAT, na.rm=T), digits=1)
211 matrx[7,y] = sum(tmp$VRegionMutations) 206 matrx[6,z] = round(matrx[6,x] / matrx[6,y] * 100, digits=1)
212 matrx[7,z] = round(matrx[7,x] / matrx[7,y] * 100, digits=1) 207
213 208 matrx[7,x] = round(f(tmp$totalMutationsAtAT, na.rm=T), digits=1)
214 matrx[8,x] = sum(tmp$nonSilentMutationsFR) 209 matrx[7,y] = round(f(tmp$VRegionMutations, na.rm=T), digits=1)
215 matrx[8,y] = sum(tmp$silentMutationsFR) 210 matrx[7,z] = round(matrx[7,x] / matrx[7,y] * 100, digits=1)
216 matrx[8,z] = round(matrx[8,x] / matrx[8,y], digits=1) 211
217 212 matrx[8,x] = round(f(tmp$nonSilentMutationsFR, na.rm=T), digits=1)
218 matrx[9,x] = sum(tmp$nonSilentMutationsCDR) 213 matrx[8,y] = round(f(tmp$silentMutationsFR, na.rm=T), digits=1)
219 matrx[9,y] = sum(tmp$silentMutationsCDR) 214 matrx[8,z] = round(matrx[8,x] / matrx[8,y], digits=1)
220 matrx[9,z] = round(matrx[9,x] / matrx[9,y], digits=1) 215
221 216 matrx[9,x] = round(f(tmp$nonSilentMutationsCDR, na.rm=T), digits=1)
217 matrx[9,y] = round(f(tmp$silentMutationsCDR, na.rm=T), digits=1)
218 matrx[9,z] = round(matrx[9,x] / matrx[9,y], digits=1)
219 }
222 220
223 transitionTable = data.frame(A=zeros,C=zeros,G=zeros,T=zeros) 221 transitionTable = data.frame(A=zeros,C=zeros,G=zeros,T=zeros)
224 row.names(transitionTable) = c("A", "C", "G", "T") 222 row.names(transitionTable) = c("A", "C", "G", "T")
225 transitionTable["A","A"] = NA 223 transitionTable["A","A"] = NA
226 transitionTable["C","C"] = NA 224 transitionTable["C","C"] = NA
248 } 246 }
249 } 247 }
250 } 248 }
251 249
252 250
253 write.table(x=transitionTable, file=paste("transitions_", gene ,".txt", sep=""), sep=",",quote=F,row.names=T,col.names=NA) 251 print(paste("writing value file: ", name, "_", fname, "_value.txt" ,sep=""))
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) 252
255 253 write.table(x=transitionTable, file=paste("transitions_", name ,"_", fname, ".txt", sep=""), sep=",",quote=F,row.names=T,col.names=NA)
256 cat(matrx[1,x], file=paste(gene, "_value.txt" ,sep="")) 254 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)
257 cat(length(tmp$Sequence.ID), file=paste(gene, "_n.txt" ,sep="")) 255
258 } 256 cat(matrx[1,x], file=paste(name, "_", fname, "_value.txt" ,sep=""))
259 257 cat(length(tmp$Sequence.ID), file=paste(name, "_", fname, "_n.txt" ,sep=""))
260 #again for all of the data 258
261 tmp = dat 259 matrx
262 j = i 260 }
263 x = (j * 3) + 1 261
264 y = (j * 3) + 2 262 nts = c("a", "c", "g", "t")
265 z = (j * 3) + 3 263 zeros=rep(0, 4)
266 matrx[1,x] = sum(tmp$VRegionMutations) 264
267 matrx[1,y] = sum(tmp$VRegionNucleotides) 265 funcs = c(median, sum, mean)
268 matrx[1,z] = round(matrx[1,x] / matrx[1,y] * 100, digits=1) 266 fnames = c("median", "sum", "mean")
269 267
270 matrx[2,x] = sum(tmp$transitionMutations) 268 for(i in 1:length(funcs)){
271 matrx[2,y] = sum(tmp$VRegionMutations) 269 func = funcs[[i]]
272 matrx[2,z] = round(matrx[2,x] / matrx[2,y] * 100, digits=1) 270 fname = fnames[[i]]
273 271
274 matrx[3,x] = sum(tmp$transversionMutations) 272 matrx = matrix(data = 0, ncol=((length(genes) + 1) * 3),nrow=9)
275 matrx[3,y] = sum(tmp$VRegionMutations) 273
276 matrx[3,z] = round(matrx[3,x] / matrx[3,y] * 100, digits=1) 274 for(i in 1:length(genes)){
277 275 matrx = calculate_result(i, genes[i], dat, matrx, func, fname, genes[i])
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 } 276 }
328 } 277
329 write.table(x=transitionTable, file="transitions.txt", sep=",",quote=F,row.names=T,col.names=NA) 278 matrx = calculate_result(i + 1, ".*", dat, matrx, func, fname, name="all")
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) 279
331 cat(matrx[1,x], file="total_value.txt") 280 result = data.frame(matrx)
332 cat(length(tmp$Sequence.ID), file="total_n.txt") 281 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)")
333 282
334 283 write.table(x=result, file=paste("mutations_", fname, ".txt", sep=""), sep=",",quote=F,row.names=T,col.names=F)
335 284 }
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 285
341 286
342 if (!("ggplot2" %in% rownames(installed.packages()))) { 287 if (!("ggplot2" %in% rownames(installed.packages()))) {
343 install.packages("ggplot2", repos="http://cran.xl-mirror.nl/") 288 install.packages("ggplot2", repos="http://cran.xl-mirror.nl/")
344 } 289 }