annotate shm_csr.r @ 4:477e95b098fd draft

Uploaded

1 library(data.table)
2 library(ggplot2)
3 library(reshape2)
4 
5 args <- commandArgs(trailingOnly = TRUE)
6 
7 input = args[1]
8 genes = unlist(strsplit(args[2], ","))
9 outputdir = args[3]
10 empty.region.filter = args[4]
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 
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 cleanup_columns = c("FR1.IMGT.c.a",
33 					"FR2.IMGT.g.t",
34 					"CDR1.IMGT.Nb.of.nucleotides",
35 					"CDR2.IMGT.t.a",
36 					"FR1.IMGT.c.g",
37 					"CDR1.IMGT.c.t",
38 					"FR2.IMGT.a.c",
39 					"FR2.IMGT.Nb.of.mutations",
40 					"FR2.IMGT.g.c",
41 					"FR2.IMGT.a.g",
42 					"FR3.IMGT.t.a",
43 					"FR3.IMGT.t.c",
44 					"FR2.IMGT.g.a",
45 					"FR3.IMGT.c.g",
46 					"FR1.IMGT.Nb.of.mutations",
47 					"CDR1.IMGT.g.a",
48 					"CDR1.IMGT.t.g",
49 					"CDR1.IMGT.g.c",
50 					"CDR2.IMGT.Nb.of.nucleotides",
51 					"FR2.IMGT.a.t",
52 					"CDR1.IMGT.Nb.of.mutations",
53 					"CDR3.IMGT.Nb.of.nucleotides",
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",
00 					"CDR2.IMGT.c.a",
01 					"FR1.IMGT.c.t",
02 					"FR1.IMGT.Nb.of.silent.mutations",
03 					"FR2.IMGT.Nb.of.silent.mutations",
04 					"FR3.IMGT.Nb.of.silent.mutations",
05 					"FR1.IMGT.Nb.of.nonsilent.mutations",
06 					"FR2.IMGT.Nb.of.nonsilent.mutations",
07 					"FR3.IMGT.Nb.of.nonsilent.mutations")
08 
09 print("Cleaning up columns")
10 
11 for(col in cleanup_columns){
12   dat[,col] = gsub("\\(.*\\)", "", dat[,col])
13   #dat[dat[,col] == "",] = "0"
14   dat[,col] = as.numeric(dat[,col])
15   dat[is.na(dat[,col]),col] = 0
16 }
17 
18 regions = c("FR1", "CDR1", "FR2", "CDR2", "FR3")
19 if(empty.region.filter == "FR1") {
20 	regions = c("CDR1", "FR2", "CDR2", "FR3")
21 } else if (empty.region.filter == "CDR1") {
22 	regions = c("FR2", "CDR2", "FR3", "CDR3")
23 } else if (empty.region.filter == "FR2") {
24 	regions = c("CDR2", "FR3", "CDR3")
25 }
26 
27 sum_by_row = function(x, columns) { sum(as.numeric(x[columns]), na.rm=T) }
28 
29 print("aggregating data into new columns")
30 
31 VRegionMutations_columns = paste(regions, ".IMGT.Nb.of.mutations", sep="")
32 dat$VRegionMutations =  apply(dat, FUN=sum_by_row, 1, columns=VRegionMutations_columns)
33 
34 VRegionNucleotides_columns = paste(regions, ".IMGT.Nb.of.nucleotides", sep="")
35 dat$FR3.IMGT.Nb.of.nucleotides = nchar(dat$FR3.IMGT.seq)
36 dat$VRegionNucleotides =  apply(dat, FUN=sum_by_row, 1, columns=VRegionNucleotides_columns)
37 
38 transitionMutations_columns = paste(rep(regions, each=4), c(".IMGT.a.g", ".IMGT.g.a", ".IMGT.c.t", ".IMGT.t.c"), sep="")
39 dat$transitionMutations = apply(dat, FUN=sum_by_row, 1, columns=transitionMutations_columns)
40 
41 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="")
42 dat$transversionMutations = apply(dat, FUN=sum_by_row, 1, columns=transversionMutations_columns)
43 
44 transitionMutationsAtGC_columns = paste(rep(regions, each=2), c(".IMGT.g.a",".IMGT.c.t"), sep="")
45 dat$transitionMutationsAtGC = apply(dat, FUN=sum_by_row, 1, columns=transitionMutationsAtGC_columns)
46 
47 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="")
48 #totalMutationsAtGC_columns = paste(rep(regions, each=6), c(".IMGT.g.a",".IMGT.c.t",".IMGT.c.a",".IMGT.c.g",".IMGT.g.t"), sep="")
49 dat$totalMutationsAtGC = apply(dat, FUN=sum_by_row, 1, columns=totalMutationsAtGC_columns)
50 
51 transitionMutationsAtAT_columns = paste(rep(regions, each=2), c(".IMGT.a.g",".IMGT.t.c"), sep="")
52 dat$transitionMutationsAtAT = apply(dat, FUN=sum_by_row, 1, columns=transitionMutationsAtAT_columns)
53 
54 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="")
55 #totalMutationsAtAT_columns = paste(rep(regions, each=5), c(".IMGT.a.g",".IMGT.t.c",".IMGT.a.c",".IMGT.g.c",".IMGT.t.g"), sep="")
56 dat$totalMutationsAtAT = apply(dat, FUN=sum_by_row, 1, columns=totalMutationsAtAT_columns)
57 
58 FRRegions = regions[grepl("FR", regions)]
59 CDRRegions = regions[grepl("CDR", regions)]
60 
61 FR_silentMutations_columns = paste(FRRegions, ".IMGT.Nb.of.silent.mutations", sep="")
62 dat$silentMutationsFR = apply(dat, FUN=sum_by_row, 1, columns=FR_silentMutations_columns)
63 
64 CDR_silentMutations_columns = paste(CDRRegions, ".IMGT.Nb.of.silent.mutations", sep="")
65 dat$silentMutationsCDR = apply(dat, FUN=sum_by_row, 1, columns=CDR_silentMutations_columns)
66 
67 FR_nonSilentMutations_columns = paste(FRRegions, ".IMGT.Nb.of.nonsilent.mutations", sep="")
68 dat$nonSilentMutationsFR = apply(dat, FUN=sum_by_row, 1, columns=FR_nonSilentMutations_columns)
69 
70 CDR_nonSilentMutations_columns = paste(CDRRegions, ".IMGT.Nb.of.nonsilent.mutations", sep="")
71 dat$nonSilentMutationsCDR = apply(dat, FUN=sum_by_row, 1, columns=CDR_nonSilentMutations_columns)
72 
73 mutation.sum.columns = c("Sequence.ID", "VRegionMutations", "VRegionNucleotides", "transitionMutations", "transversionMutations", "transitionMutationsAtGC", "transitionMutationsAtAT", "silentMutationsFR", "nonSilentMutationsFR", "silentMutationsCDR", "nonSilentMutationsCDR")
74 write.table(dat[,mutation.sum.columns], "mutation_by_id.txt", sep="\t",quote=F,row.names=F,col.names=T)
75 
76 setwd(outputdir)
77 
78 base.order = data.frame(base=c("A", "T", "C", "G"), order=1:4)
79 
80 calculate_result = function(i, gene, dat, matrx, f, fname, name){
81 	tmp = dat[grepl(paste("^", gene, ".*", sep=""), dat$best_match),]
82 
83 	j = i - 1
84 	x = (j * 3) + 1
85 	y = (j * 3) + 2
86 	z = (j * 3) + 3
87 
88 	if(nrow(tmp) > 0){
89 		if(fname == "sum"){
90 			matrx[1,x] = round(f(tmp$VRegionMutations, na.rm=T), digits=1)
91 			matrx[1,y] = round(f(tmp$VRegionNucleotides, na.rm=T), digits=1)
92 			matrx[1,z] = round(f(matrx[1,x] / matrx[1,y]) * 100, digits=1)
93 		} else {
94 			matrx[1,x] = round(f(tmp$VRegionMutations, na.rm=T), digits=1)
95 			matrx[1,y] = round(f(tmp$VRegionNucleotides, na.rm=T), digits=1)
96 			matrx[1,z] = round(f(tmp$VRegionMutations / tmp$VRegionNucleotides) * 100, digits=1)
97 		}
98 
99 		matrx[2,x] = round(f(tmp$transitionMutations, na.rm=T), digits=1)
00 		matrx[2,y] = round(f(tmp$VRegionMutations, na.rm=T), digits=1)
01 		matrx[2,z] = round(matrx[2,x] / matrx[2,y] * 100, digits=1)
02 
03 		matrx[3,x] = round(f(tmp$transversionMutations, na.rm=T), digits=1)
04 		matrx[3,y] = round(f(tmp$VRegionMutations, na.rm=T), digits=1)
05 		matrx[3,z] = round(matrx[3,x] / matrx[3,y] * 100, digits=1)
06 
07 		matrx[4,x] = round(f(tmp$transitionMutationsAtGC, na.rm=T), digits=1)
08 		matrx[4,y] = round(f(tmp$totalMutationsAtGC, na.rm=T), digits=1)
09 		matrx[4,z] = round(matrx[4,x] / matrx[4,y] * 100, digits=1)
10 
11 		matrx[5,x] = round(f(tmp$totalMutationsAtGC, na.rm=T), digits=1)
12 		matrx[5,y] = round(f(tmp$VRegionMutations, na.rm=T), digits=1)
13 		matrx[5,z] = round(matrx[5,x] / matrx[5,y] * 100, digits=1)
14 
15 		matrx[6,x] = round(f(tmp$transitionMutationsAtAT, na.rm=T), digits=1)
16 		matrx[6,y] = round(f(tmp$totalMutationsAtAT, na.rm=T), digits=1)
17 		matrx[6,z] = round(matrx[6,x] / matrx[6,y] * 100, digits=1)
18 
19 		matrx[7,x] = round(f(tmp$totalMutationsAtAT, na.rm=T), digits=1)
20 		matrx[7,y] = round(f(tmp$VRegionMutations, na.rm=T), digits=1)
21 		matrx[7,z] = round(matrx[7,x] / matrx[7,y] * 100, digits=1)
22 
23 		matrx[8,x] = round(f(tmp$nonSilentMutationsFR, na.rm=T), digits=1)
24 		matrx[8,y] = round(f(tmp$silentMutationsFR, na.rm=T), digits=1)
25 		matrx[8,z] = round(matrx[8,x] / matrx[8,y], digits=1)
26 
27 		matrx[9,x] = round(f(tmp$nonSilentMutationsCDR, na.rm=T), digits=1)
28 		matrx[9,y] = round(f(tmp$silentMutationsCDR, na.rm=T), digits=1)
29 		matrx[9,z] = round(matrx[9,x] / matrx[9,y], digits=1)
30 
31 		if(fname == "sum"){
32 			matrx[10,x] = round(f(rowSums(tmp[,c("FR2.IMGT.Nb.of.nucleotides", "FR3.IMGT.Nb.of.nucleotides")], na.rm=T)), digits=1)
33 			matrx[10,y] = round(f(tmp$VRegionNucleotides, na.rm=T), digits=1)
34 			matrx[10,z] = round(matrx[10,x] / matrx[10,y] * 100, digits=1)
35 
36 			matrx[11,x] = round(f(rowSums(tmp[,c("CDR1.IMGT.Nb.of.nucleotides", "CDR2.IMGT.Nb.of.nucleotides")], na.rm=T)), digits=1)
37 			matrx[11,y] = round(f(tmp$VRegionNucleotides, na.rm=T), digits=1)
38 			matrx[11,z] = round(matrx[11,x] / matrx[11,y] * 100, digits=1)
39 		}
40 	}
41   
42 	transitionTable = data.frame(A=zeros,C=zeros,G=zeros,T=zeros)
43 	row.names(transitionTable) = c("A", "C", "G", "T")
44 	transitionTable["A","A"] = NA
45 	transitionTable["C","C"] = NA
46 	transitionTable["G","G"] = NA
47 	transitionTable["T","T"] = NA
48 
49 	if(nrow(tmp) > 0){
50 		for(nt1 in nts){
51 			for(nt2 in nts){
52 				if(nt1 == nt2){
53 					next
54 				}
55 				NT1 = LETTERS[letters == nt1]
56 				NT2 = LETTERS[letters == nt2]
57 				FR1 = paste("FR1.IMGT.", nt1, ".", nt2, sep="")
58 				CDR1 = paste("CDR1.IMGT.", nt1, ".", nt2, sep="")
59 				FR2 = paste("FR2.IMGT.", nt1, ".", nt2, sep="")
60 				CDR2 = paste("CDR2.IMGT.", nt1, ".", nt2, sep="")
61 				FR3 = paste("FR3.IMGT.", nt1, ".", nt2, sep="")
62 				if (empty.region.filter == "leader"){
63 					transitionTable[NT1,NT2] = sum(tmp[,c(FR1, CDR1, FR2, CDR2, FR3)])
64 				} else if (empty.region.filter == "FR1") {
65 					transitionTable[NT1,NT2] = sum(tmp[,c(CDR1, FR2, CDR2, FR3)])
66 				} else if (empty.region.filter == "CDR1") {
67 					transitionTable[NT1,NT2] = sum(tmp[,c(FR2, CDR2, FR3)])
68 				} else if (empty.region.filter == "FR2") {
69 					transitionTable[NT1,NT2] = sum(tmp[,c(CDR2, FR3)])
70 				}
71 			}
72 		}
73 		transition = transitionTable
74 		transition$id = names(transition)
75 		
76 		transition2 = melt(transition, id.vars="id")
77 
78 		transition2 = merge(transition2, base.order, by.x="id", by.y="base")
79 
80 		transition2 = merge(transition2, base.order, by.x="variable", by.y="base")
81 
82 		transition2[is.na(transition2$value),]$value = 0
83 		
84 		print(transition2)
85 
86 		if(any(transition2$value != 0)){ #having rows of data but a transition table filled with 0 is bad
87 			print("Plotting stacked transition")
88 			png(filename=paste("transitions_stacked_", name, ".png", sep=""))
89 			p = ggplot(transition2, aes(factor(reorder(id, order.x)), y=value, fill=factor(reorder(variable, order.y)))) + geom_bar(position="fill", stat="identity", colour="black") #stacked bar
90 			p = p + xlab("From base") + ylab("To base") + ggtitle("Mutations frequency from base to base") + guides(fill=guide_legend(title=NULL))
91 			p = p + theme(panel.background = element_rect(fill = "white", colour="black"), text = element_text(size=13, colour="black")) + scale_fill_manual(values=c("A" = "blue4", "G" = "lightblue1", "C" = "olivedrab3", "T" = "olivedrab4"))
92 			#p = p + scale_colour_manual(values=c("A" = "black", "G" = "black", "C" = "black", "T" = "black"))
93 			print(p)
94 			dev.off()
95 
96 			print("Plotting heatmap transition")
97 
98 			png(filename=paste("transitions_heatmap_", name, ".png", sep=""))
99 			p = ggplot(transition2, aes(factor(reorder(id, order.x)), factor(reorder(variable, order.y)))) + geom_tile(aes(fill = value)) + scale_fill_gradient(low="white", high="steelblue") #heatmap
00 			p = p + xlab("From base") + ylab("To base") + ggtitle("Mutations frequency from base to base")  + theme(panel.background = element_rect(fill = "white", colour="black"), text = element_text(size=13, colour="black"))
01 			print(p)
02 			dev.off()
03 		} else {
04 			print("No data to plot")
05 		}
06 	}
07 
08 	#print(paste("writing value file: ", name, "_", fname, "_value.txt" ,sep=""))
09 	write.table(x=transitionTable, file=paste("transitions_", name ,"_", fname, ".txt", sep=""), sep=",",quote=F,row.names=T,col.names=NA)
10 	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)
11 	cat(matrx[1,x], file=paste(name, "_", fname, "_value.txt" ,sep=""))
12 	cat(nrow(tmp), file=paste(name, "_", fname, "_n.txt" ,sep=""))
13 	#print(paste(fname, name, nrow(tmp)))
14 	matrx
15 }
16 nts = c("a", "c", "g", "t")
17 zeros=rep(0, 4)
18 funcs = c(median, sum, mean)
19 fnames = c("median", "sum", "mean")
20 
21 print("Creating result tables")
22 
23 for(i in 1:length(funcs)){
24 	func = funcs[[i]]
25 	fname = fnames[[i]]
26 	
27 	rows = 9
28 	if(fname == "sum"){
29 		rows = 11
30 	}
31 	matrx = matrix(data = 0, ncol=((length(genes) + 1) * 3),nrow=rows)
32 	for(i in 1:length(genes)){
33 		print(paste("Creating table for", fname, genes[i]))
34 		matrx = calculate_result(i, genes[i], dat, matrx, func, fname, genes[i])
35 	}
36 	matrx = calculate_result(i + 1, ".*", dat[!grepl("unmatched", dat$best_match),], matrx, func, fname, name="all")
37 
38 	result = data.frame(matrx)
39 	if(fname == "sum"){
40 		row.names(result) = c("Number of Mutations (%)", "Transitions (%)", "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)", "nt in FR", "nt in CDR")
41 	} else {
42 		row.names(result) = c("Number of Mutations (%)", "Transitions (%)", "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)")
43 	}
44 	write.table(x=result, file=paste("mutations_", fname, ".txt", sep=""), sep=",",quote=F,row.names=T,col.names=F)
45 }
46 
47 print("Adding median number of mutations to sum table")
48 sum.table = read.table("mutations_sum.txt", sep=",", header=F)
49 median.table = read.table("mutations_median.txt", sep=",", header=F)
50 
51 new.table = sum.table[1,]
52 new.table[2,] = median.table[1,]
53 new.table[3:12,] = sum.table[2:11,]
54 new.table[,1] = as.character(new.table[,1])
55 new.table[2,1] = "Median of Number of Mutations (%)"
56 
57 #sum.table = sum.table[c("Number of Mutations (%)", "Median of 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)", "nt in FR", "nt in CDR"),]
58 
59 write.table(x=new.table, file="mutations_sum.txt", sep=",",quote=F,row.names=F,col.names=F)
60 
61 print("Plotting IGA piechart")
62 
63 dat = dat[!grepl("^unmatched", dat$best_match),]
64 
65 #blegh
66 
67 genesForPlot = dat[grepl("IGA", dat$best_match),]$best_match
68 
69 if(length(genesForPlot) > 0){
70 	genesForPlot = data.frame(table(genesForPlot))
71 	colnames(genesForPlot) = c("Gene","Freq")
72 	genesForPlot$label = paste(genesForPlot$Gene, "-", genesForPlot$Freq)
73 
74 	pc = ggplot(genesForPlot, aes(x = factor(1), y=Freq, fill=Gene))
75 	pc = pc + geom_bar(width = 1, stat = "identity") + scale_fill_manual(labels=genesForPlot$label, values=c("IGA1" = "lightblue1", "IGA2" = "blue4"))
76 	pc = pc + coord_polar(theta="y") + scale_y_continuous(breaks=NULL)
77 	pc = pc + theme(panel.background = element_rect(fill = "white", colour="black"), text = element_text(size=13, colour="black"))
78 	pc = pc + xlab(" ") + ylab(" ") + ggtitle(paste("IGA subclasses", "( n =", sum(genesForPlot$Freq), ")"))
79 	write.table(genesForPlot, "IGA_pie.txt", sep="\t",quote=F,row.names=F,col.names=T)
80 
81 	png(filename="IGA.png")
82 	print(pc)
83 	dev.off()
84 }
85 
86 print("Plotting IGG piechart")
87 
88 genesForPlot = dat[grepl("IGG", dat$best_match),]$best_match
89 
90 if(length(genesForPlot) > 0){
91 	genesForPlot = data.frame(table(genesForPlot))
92 	colnames(genesForPlot) = c("Gene","Freq")
93 	genesForPlot$label = paste(genesForPlot$Gene, "-", genesForPlot$Freq)
94 
95 	pc = ggplot(genesForPlot, aes(x = factor(1), y=Freq, fill=Gene))
96 	pc = pc + geom_bar(width = 1, stat = "identity") + scale_fill_manual(labels=genesForPlot$label, values=c("IGG1" = "olivedrab3", "IGG2" = "red", "IGG3" = "gold", "IGG4" = "darkred"))
97 	pc = pc + coord_polar(theta="y") + scale_y_continuous(breaks=NULL)
98 	pc = pc + theme(panel.background = element_rect(fill = "white", colour="black"), text = element_text(size=13, colour="black"))
99 	pc = pc + xlab(" ") + ylab(" ") + ggtitle(paste("IGG subclasses", "( n =", sum(genesForPlot$Freq), ")"))
00 	write.table(genesForPlot, "IGG_pie.txt", sep="\t",quote=F,row.names=F,col.names=T)
01 
02 	png(filename="IGG.png")
03 	print(pc)
04 	dev.off()
05 }
06 
07 print("Plotting scatterplot")
08 
09 dat$percentage_mutations = round(dat$VRegionMutations / dat$VRegionNucleotides * 100, 2)
10 dat.clss = dat
11 
12 dat.clss$best_match = substr(dat.clss$best_match, 0, 3)
13 
14 dat.clss = rbind(dat, dat.clss)
15 
16 p = ggplot(dat.clss, aes(best_match, percentage_mutations))
17 p = p + geom_point(aes(colour=best_match), position="jitter") + geom_boxplot(aes(middle=mean(percentage_mutations)), alpha=0.1, outlier.shape = NA)
18 p = p + xlab("Subclass") + ylab("Frequency") + ggtitle("Frequency scatter plot") + theme(panel.background = element_rect(fill = "white", colour="black"), text = element_text(size=13, colour="black"))
19 p = p + scale_fill_manual(values=c("IGA" = "blue4", "IGA1" = "lightblue1", "IGA2" = "blue4", "IGG" = "olivedrab3", "IGG1" = "olivedrab3", "IGG2" = "red", "IGG3" = "gold", "IGG4" = "darkred", "IGM" = "darkviolet", "all" = "blue4"))
20 p = p + scale_colour_manual(values=c("IGA" = "blue4", "IGA1" = "lightblue1", "IGA2" = "blue4", "IGG" = "olivedrab3", "IGG1" = "olivedrab3", "IGG2" = "red", "IGG3" = "gold", "IGG4" = "darkred", "IGM" = "darkviolet", "all" = "blue4"))
21 
22 png(filename="scatter.png")
23 print(p)
24 dev.off()
25 
26 write.table(dat[,c("Sequence.ID", "best_match", "VRegionMutations", "VRegionNucleotides", "percentage_mutations")], "scatter.txt", sep="\t",quote=F,row.names=F,col.names=T)
27 
28 write.table(dat, input, sep="\t",quote=F,row.names=F,col.names=T)
29 
30 print("Plotting frequency ranges plot")
31 
32 dat$best_match_class = substr(dat$best_match, 0, 3)
33 freq_labels = c("0", "0-2", "2-5", "5-10", "10-15", "15-20", "20")
34 dat$frequency_bins = cut(dat$percentage_mutations, breaks=c(-Inf, 0, 2,5,10,15,20, Inf), labels=freq_labels)
35 
36 frequency_bins_sum = data.frame(data.table(dat)[, list(class_sum=sum(.N)), by=c("best_match_class")])
37 
38 frequency_bins_data = data.frame(data.table(dat)[, list(frequency_count=.N), by=c("best_match_class", "frequency_bins")])
39 
40 frequency_bins_data = merge(frequency_bins_data, frequency_bins_sum, by="best_match_class")
41 
42 frequency_bins_data$frequency = round(frequency_bins_data$frequency_count / frequency_bins_data$class_sum * 100, 2)
43 
44 p = ggplot(frequency_bins_data, aes(frequency_bins, frequency))
45 p = p + geom_bar(aes(fill=best_match_class), stat="identity", position="dodge") + theme(panel.background = element_rect(fill = "white", colour="black"), text = element_text(size=13, colour="black"))
46 p = p + xlab("Frequency ranges") + ylab("Frequency") + ggtitle("Mutation Frequencies by class") + scale_fill_manual(values=c("IGA" = "blue4", "IGG" = "olivedrab3", "IGM" = "black", "all" = "blue4"))
47 
48 png(filename="frequency_ranges.png")
49 print(p)
50 dev.off()
51 
52 frequency_bins_data_by_class = frequency_bins_data
53 
54 write.table(frequency_bins_data_by_class, "frequency_ranges_classes.txt", sep="\t",quote=F,row.names=F,col.names=T)
55 
56 frequency_bins_data = data.frame(data.table(dat)[, list(frequency_count=.N), by=c("best_match", "best_match_class", "frequency_bins")])
57 
58 frequency_bins_data = merge(frequency_bins_data, frequency_bins_sum, by="best_match_class")
59 
60 frequency_bins_data$frequency = round(frequency_bins_data$frequency_count / frequency_bins_data$class_sum * 100, 2)
61 
62 write.table(frequency_bins_data, "frequency_ranges_subclasses.txt", sep="\t",quote=F,row.names=F,col.names=T)
63 
64 
65 
66 
67 
68 
69 
70 
71 
72 
73 
74 
75 
76 
77 
78 
79 
80 
81 
82 
83 
84 
85 
86 
87 
88 
89 
90 
91 
92 
93 
94 
95 
96 
97 
98 
99 
00 
01 
02 
03 
04 
05 
06 
07 
08 
09 
10 
11 
12 
13 
14 
15 
16 
17 
18