Mercurial > repos > davidvanzessen > shm_csr

comparison shm_csr.r @ 43:77a7ac76c7b9 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author davidvanzessen
date Tue, 11 Apr 2017 08:02:17 -0400
parents 1cf60ae234b4
children aa8d37bd1930
comparison
equal deleted inserted replaced
42:1cf60ae234b4 43:77a7ac76c7b9
122 regions = c("FR2", "CDR2", "FR3") 122 regions = c("FR2", "CDR2", "FR3")
123 } else if (empty.region.filter == "FR2") { 123 } else if (empty.region.filter == "FR2") {
124 regions = c("CDR2", "FR3") 124 regions = c("CDR2", "FR3")
125 } 125 }
126 126
127 pdfplots = list() #save() this later to create the pdf plots in another script (maybe avoids the "address (nil), cause memory not mapped")
128
127 sum_by_row = function(x, columns) { sum(as.numeric(x[columns]), na.rm=T) } 129 sum_by_row = function(x, columns) { sum(as.numeric(x[columns]), na.rm=T) }
128 130
129 print("aggregating data into new columns") 131 print("aggregating data into new columns")
130 132
131 VRegionMutations_columns = paste(regions, ".IMGT.Nb.of.mutations", sep="") 133 VRegionMutations_columns = paste(regions, ".IMGT.Nb.of.mutations", sep="")
306 p = p + theme(panel.background = element_rect(fill = "white", colour="black"), text = element_text(size=16, colour="black")) + scale_fill_manual(values=c("A" = "blue4", "G" = "lightblue1", "C" = "olivedrab3", "T" = "olivedrab4")) 308 p = p + theme(panel.background = element_rect(fill = "white", colour="black"), text = element_text(size=16, colour="black")) + scale_fill_manual(values=c("A" = "blue4", "G" = "lightblue1", "C" = "olivedrab3", "T" = "olivedrab4"))
307 #p = p + scale_colour_manual(values=c("A" = "black", "G" = "black", "C" = "black", "T" = "black")) 309 #p = p + scale_colour_manual(values=c("A" = "black", "G" = "black", "C" = "black", "T" = "black"))
308 print(p) 310 print(p)
309 dev.off() 311 dev.off()
310 312
311 ggsave(paste("transitions_stacked_", name, ".pdf", sep="")) 313 pdfplots[[paste("transitions_stacked_", name, ".pdf", sep="")]] <<- p
312 314
313 png(filename=paste("transitions_heatmap_", name, ".png", sep="")) 315 png(filename=paste("transitions_heatmap_", name, ".png", sep=""))
314 p = ggplot(transition2, aes(factor(reorder(variable, -order.y)), factor(reorder(id, -order.x)))) + geom_tile(aes(fill = value)) + scale_fill_gradient(low="white", high="steelblue") #heatmap 316 p = ggplot(transition2, aes(factor(reorder(variable, -order.y)), factor(reorder(id, -order.x)))) + geom_tile(aes(fill = value)) + scale_fill_gradient(low="white", high="steelblue") #heatmap
315 p = p + xlab("To base") + ylab("From Base") + ggtitle("Heatmap transition information") + theme(panel.background = element_rect(fill = "white", colour="black"), text = element_text(size=16, colour="black")) 317 p = p + xlab("To base") + ylab("From Base") + ggtitle("Heatmap transition information") + theme(panel.background = element_rect(fill = "white", colour="black"), text = element_text(size=16, colour="black"))
316 print(p) 318 print(p)
317 dev.off() 319 dev.off()
318 320
319 ggsave(paste("transitions_heatmap_", name, ".pdf", sep="")) 321 pdfplots[[paste("transitions_heatmap_", name, ".pdf", sep="")]] <<- p
320 } else { 322 } else {
321 #print("No data to plot") 323 #print("No data to plot")
322 } 324 }
323 } 325 }
324 326
398 400
399 png(filename="IGA.png") 401 png(filename="IGA.png")
400 print(pc) 402 print(pc)
401 dev.off() 403 dev.off()
402 404
403 ggsave("IGA.pdf", pc) 405 pdfplots[["IGA.pdf"]] <- pc
404 } 406 }
405 407
406 print("Plotting IGG piechart") 408 print("Plotting IGG piechart")
407 409
408 genesForPlot = dat[grepl("IGG", dat$best_match),]$best_match 410 genesForPlot = dat[grepl("IGG", dat$best_match),]$best_match
421 423
422 png(filename="IGG.png") 424 png(filename="IGG.png")
423 print(pc) 425 print(pc)
424 dev.off() 426 dev.off()
425 427
426 ggsave("IGG.pdf", pc) 428 pdfplots[["IGG.pdf"]] <- pc
427 } 429 }
428 430
429 print("Plotting scatterplot") 431 print("Plotting scatterplot")
430 432
431 dat$percentage_mutations = round(dat$VRegionMutations / dat$VRegionNucleotides * 100, 2) 433 dat$percentage_mutations = round(dat$VRegionMutations / dat$VRegionNucleotides * 100, 2)
443 445
444 png(filename="scatter.png") 446 png(filename="scatter.png")
445 print(p) 447 print(p)
446 dev.off() 448 dev.off()
447 449
448 ggsave("scatter.pdf", p) 450 pdfplots[["scatter.pdf"]] <- p
449 451
450 write.table(dat[,c("Sequence.ID", "best_match", "VRegionMutations", "VRegionNucleotides", "percentage_mutations")], "scatter.txt", sep="\t",quote=F,row.names=F,col.names=T) 452 write.table(dat[,c("Sequence.ID", "best_match", "VRegionMutations", "VRegionNucleotides", "percentage_mutations")], "scatter.txt", sep="\t",quote=F,row.names=F,col.names=T)
451 453
452 print("Plotting frequency ranges plot") 454 print("Plotting frequency ranges plot")
453 455
469 471
470 png(filename="frequency_ranges.png") 472 png(filename="frequency_ranges.png")
471 print(p) 473 print(p)
472 dev.off() 474 dev.off()
473 475
474 ggsave("frequency_ranges.pdf", p) 476 pdfplots[["frequency_ranges.pdf"]] <- p
477
478 save(pdfplots, file="pdfplots.RData")
475 479
476 frequency_bins_data_by_class = frequency_bins_data 480 frequency_bins_data_by_class = frequency_bins_data
477 481
478 frequency_bins_data_by_class = frequency_bins_data_by_class[order(frequency_bins_data_by_class$best_match_class, frequency_bins_data_by_class$frequency_bins),] 482 frequency_bins_data_by_class = frequency_bins_data_by_class[order(frequency_bins_data_by_class$best_match_class, frequency_bins_data_by_class$frequency_bins),]
479 483