changeset 1:faae21ba5c63 draft

Uploaded
author davidvanzessen
date Tue, 25 Oct 2016 07:28:43 -0400
parents c33d93683a09
children e85fec274cde
files aa_histogram.r merge_and_filter.r pattern_plots.r shm_csr.py shm_csr.r shm_csr.xml wrapper.sh
diffstat 7 files changed, 591 insertions(+), 567 deletions(-) [+]
line wrap: on
line diff
--- a/aa_histogram.r	Thu Oct 13 10:52:24 2016 -0400
+++ b/aa_histogram.r	Tue Oct 25 07:28:43 2016 -0400
@@ -37,7 +37,7 @@
 
         print("---------------- plot ----------------")
 
-        m = ggplot(dat_dt, aes(x=i, y=freq)) + theme(axis.text.x = element_text(angle = 90, hjust = 1))
+        m = ggplot(dat_dt, aes(x=i, y=freq)) + theme(axis.text.x = element_text(angle = 90, hjust = 1), text = element_text(size=13, colour="black"))
         m = m + geom_bar(stat="identity", colour = "black", fill = "darkgrey", alpha=0.8) + scale_x_continuous(breaks=dat_dt$i, labels=dat_dt$i)
         m = m + annotate("segment", x = 0.5, y = -0.05, xend=26.5, yend=-0.05, colour="darkgreen", size=1) + annotate("text", x = 13, y = -0.1, label="FR1")
         m = m + annotate("segment", x = 26.5, y = -0.07, xend=38.5, yend=-0.07, colour="darkblue", size=1) + annotate("text", x = 32.5, y = -0.15, label="CDR1")
--- a/merge_and_filter.r	Thu Oct 13 10:52:24 2016 -0400
+++ b/merge_and_filter.r	Tue Oct 25 07:28:43 2016 -0400
@@ -32,7 +32,6 @@
 	gene_identification$chunk_hit_percentage = (gene_identification$length / gene_identification$ref.length) * 100
 	gene_identification = gene_identification[,c("qseqid", "chunk_hit_percentage", "pident", "qstart", "sseqid")]
 	colnames(gene_identification) = c("Sequence.ID", "chunk_hit_percentage", "nt_hit_percentage", "start_locations", "best_match")
-	
 }
 
 input.sequence.count = nrow(summ)
@@ -62,32 +61,9 @@
 	summ = summ[summ$Functionality != "No results" & summ$Functionality != "unknown (see comment)" & summ$Functionality != "unknown",]
 }
 
-print(paste("Number of sequences after productive filter:", nrow(summ)))
-
-filtering.steps = rbind(filtering.steps, c("After productive filter", nrow(summ)))
-
-splt = strsplit(class.filter, "_")[[1]]
-chunk_hit_threshold = as.numeric(splt[1])
-nt_hit_threshold = as.numeric(splt[2])
-
-higher_than=(summ$chunk_hit_percentage >= chunk_hit_threshold & summ$nt_hit_percentage >= nt_hit_threshold)
-
-unmatched=summ[NULL,c("Sequence.ID", "chunk_hit_percentage", "nt_hit_percentage", "start_locations", "best_match")]
+print(paste("Number of sequences after functionality filter:", nrow(summ)))
 
-if(!all(higher_than, na.rm=T)){ #check for 'not all' because that would mean the unmatched set is empty
-	unmatched = summ[!higher_than,]
-	unmatched = unmatched[,c("Sequence.ID", "chunk_hit_percentage", "nt_hit_percentage", "start_locations", "best_match")]
-	unmatched$best_match = paste("unmatched,", unmatched$best_match)
-	summ[!higher_than,"best_match"] = paste("unmatched,", summ[!higher_than,"best_match"])
-}
-
-if(any(higher_than, na.rm=T)){
-	#summ = summ[higher_than,]
-}
-
-if(nrow(summ) == 0){
-	stop("No data remaining after filter")
-}
+filtering.steps = rbind(filtering.steps, c("After functionality filter", nrow(summ)))
 
 result = merge(summ, mutationanalysis[,!(names(mutationanalysis) %in% names(summ)[-1])], by="Sequence.ID")
 
@@ -114,22 +90,16 @@
 result$JGene = gsub("^Homsap ", "", result$J.GENE.and.allele)
 result$JGene = gsub("[*].*", "", result$JGene)
 
-result$past = do.call(paste, c(result[unlist(strsplit(unique.type, ","))], sep = ":"))
-
-result = result[!(duplicated(result$past)), ]
-
-result = result[,!(names(result) %in% c("past"))]
-
-print(paste("Number of sequences in result after", unique.type, "filtering:", nrow(result)))
-
-filtering.steps = rbind(filtering.steps, c("After duplicate filter", nrow(result)))
-
 print(paste("Number of empty CDR1 sequences:", sum(result$CDR1.IMGT.seq == "")))
 print(paste("Number of empty FR2 sequences:", sum(result$FR2.IMGT.seq == "")))
 print(paste("Number of empty CDR2 sequences:", sum(result$CDR2.IMGT.seq == "")))
 print(paste("Number of empty FR3 sequences:", sum(result$FR3.IMGT.seq == "")))
 
-if(empty.region.filter == "FR1"){
+if(empty.region.filter == "leader"){
+	result = result[result$FR1.IMGT.seq != "" & result$CDR1.IMGT.seq != "" & result$FR2.IMGT.seq != "" & result$CDR2.IMGT.seq != "" & result$FR3.IMGT.seq != "", ]
+	print(paste("Number of sequences after empty FR1, CDR1, FR2, CDR2 and FR3 column filter:", nrow(result)))
+	filtering.steps = rbind(filtering.steps, c("After empty FR1, CDR1, FR2, CDR2, FR3 filter", nrow(result)))
+} else if(empty.region.filter == "FR1"){
 	result = result[result$CDR1.IMGT.seq != "" & result$FR2.IMGT.seq != "" & result$CDR2.IMGT.seq != "" & result$FR3.IMGT.seq != "", ]
 	print(paste("Number of sequences after empty CDR1, FR2, CDR2 and FR3 column filter:", nrow(result)))
 	filtering.steps = rbind(filtering.steps, c("After empty CDR1, FR2, CDR2, FR3 filter", nrow(result)))
@@ -143,7 +113,12 @@
 	filtering.steps = rbind(filtering.steps, c("After empty CDR2, FR3 filter", nrow(result)))
 }
 
-if(empty.region.filter == "FR1"){
+print(paste("Number of sequences in result after CDR/FR filtering:", nrow(result)))
+print(paste("Number of matched sequences in result after CDR/FR filtering:", nrow(result[!grepl("unmatched", result$best_match),])))
+
+if(empty.region.filter == "leader"){
+	result = result[!(grepl("n|N", result$FR1.IMGT.seq) | grepl("n|N", result$FR2.IMGT.seq) | grepl("n|N", result$FR3.IMGT.seq) | grepl("n|N", result$CDR1.IMGT.seq) | grepl("n|N", result$CDR2.IMGT.seq) | grepl("n|N", result$CDR3.IMGT.seq)),]
+} else if(empty.region.filter == "FR1"){
 	result = result[!(grepl("n|N", result$FR2.IMGT.seq) | grepl("n|N", result$FR3.IMGT.seq) | grepl("n|N", result$CDR1.IMGT.seq) | grepl("n|N", result$CDR2.IMGT.seq) | grepl("n|N", result$CDR3.IMGT.seq)),]
 } else if(empty.region.filter == "CDR1"){
 	result = result[!(grepl("n|N", result$FR2.IMGT.seq) | grepl("n|N", result$FR3.IMGT.seq) | grepl("n|N", result$CDR2.IMGT.seq) | grepl("n|N", result$CDR3.IMGT.seq)),]
@@ -171,7 +146,9 @@
 if(filter.unique != "no"){
 	clmns = names(result)
 	
-	if(empty.region.filter == "FR1"){
+	if(empty.region.filter == "leader"){
+		result$unique.def = paste(result$FR1.IMGT.seq, result$CDR1.IMGT.seq, result$FR2.IMGT.seq, result$CDR2.IMGT.seq, result$FR3.IMGT.seq, result$CDR3.IMGT.seq)
+	} else if(empty.region.filter == "FR1"){
 		result$unique.def = paste(result$CDR1.IMGT.seq, result$FR2.IMGT.seq, result$CDR2.IMGT.seq, result$FR3.IMGT.seq, result$CDR3.IMGT.seq)
 	} else if(empty.region.filter == "CDR1"){
 		rresult$unique.def = paste(result$FR2.IMGT.seq, result$CDR2.IMGT.seq, result$FR3.IMGT.seq, result$CDR3.IMGT.seq)
@@ -199,10 +176,41 @@
 	#write.table(inputdata.removed, "unique_removed.csv", sep=",",quote=F,row.names=F,col.names=T)
 }
 
-print(paste("Number of sequences in result after CDR/FR filtering:", nrow(result)))
-print(paste("Number of matched sequences in result after CDR/FR filtering:", nrow(result[!grepl("unmatched", result$best_match),])))
+filtering.steps = rbind(filtering.steps, c("After filter unique sequences", nrow(result)))
+
+
+splt = strsplit(class.filter, "_")[[1]]
+chunk_hit_threshold = as.numeric(splt[1])
+nt_hit_threshold = as.numeric(splt[2])
+
+higher_than=(summ$chunk_hit_percentage >= chunk_hit_threshold & summ$nt_hit_percentage >= nt_hit_threshold)
+
+unmatched=summ[NULL,c("Sequence.ID", "chunk_hit_percentage", "nt_hit_percentage", "start_locations", "best_match")]
+
+if(!all(higher_than, na.rm=T)){ #check for 'not all' because that would mean the unmatched set is empty
+	unmatched = summ[!higher_than,]
+	unmatched = unmatched[,c("Sequence.ID", "chunk_hit_percentage", "nt_hit_percentage", "start_locations", "best_match")]
+	unmatched$best_match = paste("unmatched,", unmatched$best_match)
+	summ[!higher_than,"best_match"] = paste("unmatched,", summ[!higher_than,"best_match"])
+}
 
-filtering.steps = rbind(filtering.steps, c("After unique filter", nrow(result)))
+if(any(higher_than, na.rm=T)){
+	#summ = summ[higher_than,]
+}
+
+if(nrow(summ) == 0){
+	stop("No data remaining after filter")
+}
+
+result$past = do.call(paste, c(result[unlist(strsplit(unique.type, ","))], sep = ":"))
+
+result = result[!(duplicated(result$past)), ]
+
+result = result[,!(names(result) %in% c("past"))]
+
+print(paste("Number of sequences in result after", unique.type, "filtering:", nrow(result)))
+
+filtering.steps = rbind(filtering.steps, c("After remove duplicates based on filter", nrow(result)))
 
 print(paste("Number of rows in result:", nrow(result)))
 print(paste("Number of rows in unmatched:", nrow(unmatched)))
--- a/pattern_plots.r	Thu Oct 13 10:52:24 2016 -0400
+++ b/pattern_plots.r	Tue Oct 25 07:28:43 2016 -0400
@@ -45,7 +45,7 @@
 write.table(data1, plot1.txt, quote=F, sep="\t", na="", row.names=F, col.names=T)
 
 p = ggplot(data1, aes(Class, value)) + geom_bar(aes(fill=Type), stat="identity", position="dodge", colour = "black") + ylab("% of mutations") + guides(fill=guide_legend(title=NULL))
-p = p + theme(panel.background = element_rect(fill = "white", colour="black")) + scale_fill_manual(values=c("RGYW.WRCY" = "white", "TW.WA" = "blue4"))
+p = p + theme(panel.background = element_rect(fill = "white", colour="black"),text = element_text(size=13, colour="black")) + scale_fill_manual(values=c("RGYW.WRCY" = "white", "TW.WA" = "blue4"))
 #p = p + scale_colour_manual(values=c("RGYW.WRCY" = "black", "TW.WA" = "blue4"))
 png(filename=plot1.png, width=480, height=300)
 print(p)
@@ -79,7 +79,7 @@
 write.table(data2, plot2.txt, quote=F, sep="\t", na="", row.names=F, col.names=T)
 
 p = ggplot(data2, aes(x=Class, y=value, fill=Type)) + geom_bar(position="fill", stat="identity", colour = "black") + scale_y_continuous(labels=percent_format()) + guides(fill=guide_legend(title=NULL)) + ylab("% of mutations")
-p = p + theme(panel.background = element_rect(fill = "white", colour="black")) + scale_fill_manual(values=c("A/T" = "blue4", "G/C transversions" = "gray74", "G/C transitions" = "white"))
+p = p + theme(panel.background = element_rect(fill = "white", colour="black"), text = element_text(size=13, colour="black")) + scale_fill_manual(values=c("A/T" = "blue4", "G/C transversions" = "gray74", "G/C transitions" = "white"))
 #p = p + scale_colour_manual(values=c("A/T" = "blue4", "G/C transversions" = "gray74", "G/C transitions" = "black"))
 png(filename=plot2.png, width=480, height=300)
 print(p)
@@ -115,7 +115,7 @@
 write.table(data3, plot3.txt, quote=F, sep="\t", na="", row.names=F, col.names=T)
 
 p = ggplot(data3, aes(Class, value)) + geom_bar(aes(fill=Type), stat="identity", position="dodge", colour = "black") + ylab("% of nucleotides") + guides(fill=guide_legend(title=NULL))
-p = p + theme(panel.background = element_rect(fill = "white", colour="black")) + scale_fill_manual(values=c("A/T" = "blue4", "G/C transversions" = "gray74", "G/C transitions" = "white"))
+p = p + theme(panel.background = element_rect(fill = "white", colour="black"), text = element_text(size=13, colour="black")) + scale_fill_manual(values=c("A/T" = "blue4", "G/C transversions" = "gray74", "G/C transitions" = "white"))
 #p = p + scale_colour_manual(values=c("A/T" = "blue4", "G/C transversions" = "gray74", "G/C transitions" = "black"))
 png(filename=plot3.png, width=480, height=300)
 print(p)
--- a/shm_csr.py	Thu Oct 13 10:52:24 2016 -0400
+++ b/shm_csr.py	Tue Oct 25 07:28:43 2016 -0400
@@ -7,15 +7,14 @@
 parser.add_argument("--input",
 					help="The '7_V-REGION-mutation-and-AA-change-table' and '10_V-REGION-mutation-hotspots' merged together, with an added 'best_match' annotation")
 parser.add_argument("--genes", help="The genes available in the 'best_match' column")
-parser.add_argument("--includefr1", help="Should the mutation/nucleotides in the FR1 region be included?")
+parser.add_argument("--empty_region_filter", help="Where does the sequence start?", choices=['leader', 'FR1', 'CDR1', 'FR2'])
 parser.add_argument("--output", help="Output file")
 
 args = parser.parse_args()
 
 infile = args.input
 genes = str(args.genes).split(",")
-print "includefr1 =", args.includefr1
-include_fr1 = True if args.includefr1 == "yes" else False
+empty_region_filter = args.empty_region_filter
 outfile = args.output
 
 genedic = dict()
@@ -59,16 +58,14 @@
 		ID = linesplt[IDIndex]
 		genedic[ID] = linesplt[best_matchIndex]
 		try:
-			if linesplt[fr1Index] != "NA":
-				mutationdic[ID + "_FR1"] = [mutationMatcher.match(x).groups() for x in linesplt[fr1Index].split("|") if x] if include_fr1 else []
-			else:
-				mutationdic[ID + "_FR1"] = []
-			mutationdic[ID + "_CDR1"] = [mutationMatcher.match(x).groups() for x in linesplt[cdr1Index].split("|") if x] if linesplt[cdr1Index] != "NA" else []
-			mutationdic[ID + "_FR2"] = [mutationMatcher.match(x).groups() for x in linesplt[fr2Index].split("|") if x] if linesplt[fr2Index] != "NA" else []
-			mutationdic[ID + "_CDR2"] = [mutationMatcher.match(x).groups() for x in linesplt[cdr2Index].split("|") if x] if linesplt[cdr2Index] != "NA" else []
+			mutationdic[ID + "_FR1"] = [mutationMatcher.match(x).groups() for x in linesplt[fr1Index].split("|") if x] if (linesplt[fr1Index] != "NA" and empty_region_filter == "leader") else []
+			mutationdic[ID + "_CDR1"] = [mutationMatcher.match(x).groups() for x in linesplt[cdr1Index].split("|") if x] if (linesplt[cdr1Index] != "NA" and empty_region_filter in ["leader", "FR1"]) else []
+			mutationdic[ID + "_FR2"] = [mutationMatcher.match(x).groups() for x in linesplt[fr2Index].split("|") if x] if (linesplt[fr2Index] != "NA" and empty_region_filter in ["leader", "FR1", "CDR1"]) else []
+			mutationdic[ID + "_CDR2"] = [mutationMatcher.match(x).groups() for x in linesplt[cdr2Index].split("|") if x] if (linesplt[cdr2Index] != "NA") else []
 			mutationdic[ID + "_FR2-CDR2"] = mutationdic[ID + "_FR2"] + mutationdic[ID + "_CDR2"]
 			mutationdic[ID + "_FR3"] = [mutationMatcher.match(x).groups() for x in linesplt[fr3Index].split("|") if x] if linesplt[fr3Index] != "NA" else []
-		except e:
+		except Exception as e:
+			print "Something went wrong while processing this line:"
 			print linesplt
 			print linecount
 			print e
@@ -189,8 +186,6 @@
 		linesplt = line.split("\t")
 		gene = linesplt[best_matchIndex]
 		ID = linesplt[IDIndex]
-		if ID == "ca2":
-			print linesplt
 		RGYW = [(int(x), int(y), z) for (x, y, z) in
 				[hotspotMatcher.match(x).groups() for x in linesplt[aggctatIndex].split("|") if x]]
 		WRCY = [(int(x), int(y), z) for (x, y, z) in
@@ -201,8 +196,7 @@
 			  [hotspotMatcher.match(x).groups() for x in linesplt[tatIndex].split("|") if x]]
 		RGYWCount[ID], WRCYCount[ID], WACount[ID], TWCount[ID] = 0, 0, 0, 0
 
-		mutationList = (mutationdic[ID + "_FR1"] if include_fr1 else []) + mutationdic[ID + "_CDR1"] + mutationdic[
-			ID + "_FR2"] + mutationdic[ID + "_CDR2"] + mutationdic[ID + "_FR3"]
+		mutationList = mutationdic[ID + "_FR1"] + mutationdic[ID + "_CDR1"] + mutationdic[ID + "_FR2"] + mutationdic[ID + "_CDR2"] + mutationdic[ID + "_FR3"]
 		for mutation in mutationList:
 			frm, where, to, AAfrm, AAwhere, AAto, junk = mutation
 			mutation_in_RGYW = any([(start <= int(where) <= end) for (start, end, region) in RGYW])
@@ -271,7 +265,7 @@
 			o.write(typ + " (%)")
 			curr = dic[typ]
 			for gene in genes:
-				geneMatcher = geneMatchers[gene] #re.compile("^" + gene + ".*") #recompile every loop....
+				geneMatcher = geneMatchers[gene]
 				if valuedic[gene + "_" + fname] is 0:
 					o.write(",0,0,0")
 				else:
--- a/shm_csr.r	Thu Oct 13 10:52:24 2016 -0400
+++ b/shm_csr.r	Tue Oct 25 07:28:43 2016 -0400
@@ -1,493 +1,516 @@
-library(data.table)
-library(ggplot2)
-library(reshape2)
-
-args <- commandArgs(trailingOnly = TRUE)
-
-input = args[1]
-genes = unlist(strsplit(args[2], ","))
-outputdir = args[3]
-include_fr1 = ifelse(args[4] == "yes", T, F)
-setwd(outputdir)
-
-dat = read.table(input, header=T, sep="\t", fill=T, stringsAsFactors=F)
-
-if(length(dat$Sequence.ID) == 0){
-  setwd(outputdir)
-  result = data.frame(x = rep(0, 5), y = rep(0, 5), z = rep(NA, 5))
-  row.names(result) = c("Number of Mutations (%)", "Transition (%)", "Transversions (%)", "Transitions at G C (%)", "Targeting of C G (%)")
-  write.table(x=result, file="mutations.txt", sep=",",quote=F,row.names=T,col.names=F)
-  transitionTable = data.frame(A=rep(0, 4),C=rep(0, 4),G=rep(0, 4),T=rep(0, 4))
-  row.names(transitionTable) = c("A", "C", "G", "T")
-  transitionTable["A","A"] = NA
-  transitionTable["C","C"] = NA
-  transitionTable["G","G"] = NA
-  transitionTable["T","T"] = NA
-  write.table(x=transitionTable, file="transitions.txt", sep=",",quote=F,row.names=T,col.names=NA)
-  cat("0", file="n.txt")
-  stop("No data")
-}
-
-cleanup_columns = c("FR1.IMGT.c.a",
-                    "FR2.IMGT.g.t",
-                    "CDR1.IMGT.Nb.of.nucleotides",
-                    "CDR2.IMGT.t.a",
-                    "FR1.IMGT.c.g",
-                    "CDR1.IMGT.c.t",
-                    "FR2.IMGT.a.c",
-                    "FR2.IMGT.Nb.of.mutations",
-                    "FR2.IMGT.g.c",
-                    "FR2.IMGT.a.g",
-                    "FR3.IMGT.t.a",
-                    "FR3.IMGT.t.c",
-                    "FR2.IMGT.g.a",
-                    "FR3.IMGT.c.g",
-                    "FR1.IMGT.Nb.of.mutations",
-                    "CDR1.IMGT.g.a",
-                    "CDR1.IMGT.t.g",
-                    "CDR1.IMGT.g.c",
-                    "CDR2.IMGT.Nb.of.nucleotides",
-                    "FR2.IMGT.a.t",
-                    "CDR1.IMGT.Nb.of.mutations",
-                    "CDR3.IMGT.Nb.of.nucleotides",
-                    "CDR1.IMGT.a.g",
-                    "FR3.IMGT.a.c",
-                    "FR1.IMGT.g.a",
-                    "FR3.IMGT.a.g",
-                    "FR1.IMGT.a.t",
-                    "CDR2.IMGT.a.g",
-                    "CDR2.IMGT.Nb.of.mutations",
-                    "CDR2.IMGT.g.t",
-                    "CDR2.IMGT.a.c",
-                    "CDR1.IMGT.t.c",
-                    "FR3.IMGT.g.c",
-                    "FR1.IMGT.g.t",
-                    "FR3.IMGT.g.t",
-                    "CDR1.IMGT.a.t",
-                    "FR1.IMGT.a.g",
-                    "FR3.IMGT.a.t",
-                    "FR3.IMGT.Nb.of.nucleotides",
-                    "FR2.IMGT.t.c",
-                    "CDR2.IMGT.g.a",
-                    "FR2.IMGT.t.a",
-                    "CDR1.IMGT.t.a",
-                    "FR2.IMGT.t.g",
-                    "FR3.IMGT.t.g",
-                    "FR2.IMGT.Nb.of.nucleotides",
-                    "FR1.IMGT.t.a",
-                    "FR1.IMGT.t.g",
-                    "FR3.IMGT.c.t",
-                    "FR1.IMGT.t.c",
-                    "CDR2.IMGT.a.t",
-                    "FR2.IMGT.c.t",
-                    "CDR1.IMGT.g.t",
-                    "CDR2.IMGT.t.g",
-                    "FR1.IMGT.Nb.of.nucleotides",
-                    "CDR1.IMGT.c.g",
-                    "CDR2.IMGT.t.c",
-                    "FR3.IMGT.g.a",
-                    "CDR1.IMGT.a.c",
-                    "FR2.IMGT.c.a",
-                    "FR3.IMGT.Nb.of.mutations",
-                    "FR2.IMGT.c.g",
-                    "CDR2.IMGT.g.c",
-                    "FR1.IMGT.g.c",
-                    "CDR2.IMGT.c.t",
-                    "FR3.IMGT.c.a",
-                    "CDR1.IMGT.c.a",
-                    "CDR2.IMGT.c.g",
-                    "CDR2.IMGT.c.a",
-                    "FR1.IMGT.c.t",
-                    "FR1.IMGT.Nb.of.silent.mutations",
-                    "FR2.IMGT.Nb.of.silent.mutations",
-                    "FR3.IMGT.Nb.of.silent.mutations",
-                    "FR1.IMGT.Nb.of.nonsilent.mutations",
-                    "FR2.IMGT.Nb.of.nonsilent.mutations",
-                    "FR3.IMGT.Nb.of.nonsilent.mutations")
-
-
-print("Cleaning up columns")
-for(col in cleanup_columns){
-  dat[,col] = gsub("\\(.*\\)", "", dat[,col])
-  #dat[dat[,col] == "",] = "0"
-  dat[,col] = as.numeric(dat[,col])
-  dat[is.na(dat[,col]),col] = 0
-}
-
-regions = c("FR1", "CDR1", "FR2", "CDR2", "FR3")
-if(!include_fr1){
-	regions = c("CDR1", "FR2", "CDR2", "FR3")
-}
-
-sum_by_row = function(x, columns) { sum(as.numeric(x[columns]), na.rm=T) }
-
-print("aggregating data into new columns")
-
-VRegionMutations_columns = paste(regions, ".IMGT.Nb.of.mutations", sep="")
-dat$VRegionMutations =  apply(dat, FUN=sum_by_row, 1, columns=VRegionMutations_columns)
-
-VRegionNucleotides_columns = paste(regions, ".IMGT.Nb.of.nucleotides", sep="")
-dat$FR3.IMGT.Nb.of.nucleotides = nchar(dat$FR3.IMGT.seq)
-dat$VRegionNucleotides =  apply(dat, FUN=sum_by_row, 1, columns=VRegionNucleotides_columns)
-
-transitionMutations_columns = paste(rep(regions, each=4), c(".IMGT.a.g", ".IMGT.g.a", ".IMGT.c.t", ".IMGT.t.c"), sep="")
-dat$transitionMutations = apply(dat, FUN=sum_by_row, 1, columns=transitionMutations_columns)
-
-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="")
-dat$transversionMutations = apply(dat, FUN=sum_by_row, 1, columns=transversionMutations_columns)
-
-
-transitionMutationsAtGC_columns = paste(rep(regions, each=2), c(".IMGT.g.a",".IMGT.c.t"), sep="")
-dat$transitionMutationsAtGC = apply(dat, FUN=sum_by_row, 1, columns=transitionMutationsAtGC_columns)
-
-
-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="")
-#totalMutationsAtGC_columns = paste(rep(regions, each=6), c(".IMGT.g.a",".IMGT.c.t",".IMGT.c.a",".IMGT.c.g",".IMGT.g.t"), sep="")
-dat$totalMutationsAtGC = apply(dat, FUN=sum_by_row, 1, columns=totalMutationsAtGC_columns)
-
-transitionMutationsAtAT_columns = paste(rep(regions, each=2), c(".IMGT.a.g",".IMGT.t.c"), sep="")
-dat$transitionMutationsAtAT = apply(dat, FUN=sum_by_row, 1, columns=transitionMutationsAtAT_columns)
-
-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="")
-#totalMutationsAtAT_columns = paste(rep(regions, each=5), c(".IMGT.a.g",".IMGT.t.c",".IMGT.a.c",".IMGT.g.c",".IMGT.t.g"), sep="")
-dat$totalMutationsAtAT = apply(dat, FUN=sum_by_row, 1, columns=totalMutationsAtAT_columns)
-
-
-FRRegions = regions[grepl("FR", regions)]
-CDRRegions = regions[grepl("CDR", regions)]
-
-FR_silentMutations_columns = paste(FRRegions, ".IMGT.Nb.of.silent.mutations", sep="")
-dat$silentMutationsFR = apply(dat, FUN=sum_by_row, 1, columns=FR_silentMutations_columns)
-
-CDR_silentMutations_columns = paste(CDRRegions, ".IMGT.Nb.of.silent.mutations", sep="")
-dat$silentMutationsCDR = apply(dat, FUN=sum_by_row, 1, columns=CDR_silentMutations_columns)
-
-FR_nonSilentMutations_columns = paste(FRRegions, ".IMGT.Nb.of.nonsilent.mutations", sep="")
-dat$nonSilentMutationsFR = apply(dat, FUN=sum_by_row, 1, columns=FR_nonSilentMutations_columns)
-
-CDR_nonSilentMutations_columns = paste(CDRRegions, ".IMGT.Nb.of.nonsilent.mutations", sep="")
-dat$nonSilentMutationsCDR = apply(dat, FUN=sum_by_row, 1, columns=CDR_nonSilentMutations_columns)
-
-mutation.sum.columns = c("Sequence.ID", "VRegionMutations", "VRegionNucleotides", "transitionMutations", "transversionMutations", "transitionMutationsAtGC", "transitionMutationsAtAT", "silentMutationsFR", "nonSilentMutationsFR", "silentMutationsCDR", "nonSilentMutationsCDR")
-
-write.table(dat[,mutation.sum.columns], "mutation_by_id.txt", sep="\t",quote=F,row.names=F,col.names=T)
-
-setwd(outputdir)
-
-base.order = data.frame(base=c("A", "T", "C", "G"), order=1:4)
-
-calculate_result = function(i, gene, dat, matrx, f, fname, name){
-	tmp = dat[grepl(paste("^", gene, ".*", sep=""), dat$best_match),]
-
-	j = i - 1
-	x = (j * 3) + 1
-	y = (j * 3) + 2
-	z = (j * 3) + 3
-	 
-	if(nrow(tmp) > 0){
-	  
-		if(fname == "sum"){
-		matrx[1,x] = round(f(tmp$VRegionMutations, na.rm=T), digits=1)
-		matrx[1,y] = round(f(tmp$VRegionNucleotides, na.rm=T), digits=1)
-		matrx[1,z] = round(f(matrx[1,x] / matrx[1,y]) * 100, digits=1)
-		} else {
-		matrx[1,x] = round(f(tmp$VRegionMutations, na.rm=T), digits=1)
-		matrx[1,y] = round(f(tmp$VRegionNucleotides, na.rm=T), digits=1)
-		matrx[1,z] = round(f(tmp$VRegionMutations / tmp$VRegionNucleotides) * 100, digits=1)
-		}
-
-		matrx[2,x] = round(f(tmp$transitionMutations, na.rm=T), digits=1)
-		matrx[2,y] = round(f(tmp$VRegionMutations, na.rm=T), digits=1)
-		matrx[2,z] = round(matrx[2,x] / matrx[2,y] * 100, digits=1)
-
-		matrx[3,x] = round(f(tmp$transversionMutations, na.rm=T), digits=1)
-		matrx[3,y] = round(f(tmp$VRegionMutations, na.rm=T), digits=1)
-		matrx[3,z] = round(matrx[3,x] / matrx[3,y] * 100, digits=1)
-
-		matrx[4,x] = round(f(tmp$transitionMutationsAtGC, na.rm=T), digits=1)
-		matrx[4,y] = round(f(tmp$totalMutationsAtGC, na.rm=T), digits=1)
-		matrx[4,z] = round(matrx[4,x] / matrx[4,y] * 100, digits=1)
-
-		matrx[5,x] = round(f(tmp$totalMutationsAtGC, na.rm=T), digits=1)
-		matrx[5,y] = round(f(tmp$VRegionMutations, na.rm=T), digits=1)
-		matrx[5,z] = round(matrx[5,x] / matrx[5,y] * 100, digits=1)
-
-		matrx[6,x] = round(f(tmp$transitionMutationsAtAT, na.rm=T), digits=1)
-		matrx[6,y] = round(f(tmp$totalMutationsAtAT, na.rm=T), digits=1)
-		matrx[6,z] = round(matrx[6,x] / matrx[6,y] * 100, digits=1)
-
-		matrx[7,x] = round(f(tmp$totalMutationsAtAT, na.rm=T), digits=1)
-		matrx[7,y] = round(f(tmp$VRegionMutations, na.rm=T), digits=1)
-		matrx[7,z] = round(matrx[7,x] / matrx[7,y] * 100, digits=1)
-
-		matrx[8,x] = round(f(tmp$nonSilentMutationsFR, na.rm=T), digits=1)
-		matrx[8,y] = round(f(tmp$silentMutationsFR, na.rm=T), digits=1)
-		matrx[8,z] = round(matrx[8,x] / matrx[8,y], digits=1)
-
-		matrx[9,x] = round(f(tmp$nonSilentMutationsCDR, na.rm=T), digits=1)
-		matrx[9,y] = round(f(tmp$silentMutationsCDR, na.rm=T), digits=1)
-		matrx[9,z] = round(matrx[9,x] / matrx[9,y], digits=1)
-
-		if(fname == "sum"){
-			matrx[10,x] = round(f(rowSums(tmp[,c("FR2.IMGT.Nb.of.nucleotides", "FR3.IMGT.Nb.of.nucleotides")], na.rm=T)), digits=1)
-			matrx[10,y] = round(f(tmp$VRegionNucleotides, na.rm=T), digits=1)
-			matrx[10,z] = round(matrx[10,x] / matrx[10,y] * 100, digits=1)
-
-			matrx[11,x] = round(f(rowSums(tmp[,c("CDR1.IMGT.Nb.of.nucleotides", "CDR2.IMGT.Nb.of.nucleotides")], na.rm=T)), digits=1)
-			matrx[11,y] = round(f(tmp$VRegionNucleotides, na.rm=T), digits=1)
-			matrx[11,z] = round(matrx[11,x] / matrx[11,y] * 100, digits=1)
-		}
-	}
-  
-	transitionTable = data.frame(A=zeros,C=zeros,G=zeros,T=zeros)
-	row.names(transitionTable) = c("A", "C", "G", "T")
-	transitionTable["A","A"] = NA
-	transitionTable["C","C"] = NA
-	transitionTable["G","G"] = NA
-	transitionTable["T","T"] = NA
-
-	if(nrow(tmp) > 0){
-		for(nt1 in nts){
-			for(nt2 in nts){
-				if(nt1 == nt2){
-					next
-				}
-				NT1 = LETTERS[letters == nt1]
-				NT2 = LETTERS[letters == nt2]
-				FR1 = paste("FR1.IMGT.", nt1, ".", nt2, sep="")
-				CDR1 = paste("CDR1.IMGT.", nt1, ".", nt2, sep="")
-				FR2 = paste("FR2.IMGT.", nt1, ".", nt2, sep="")
-				CDR2 = paste("CDR2.IMGT.", nt1, ".", nt2, sep="")
-				FR3 = paste("FR3.IMGT.", nt1, ".", nt2, sep="")
-				if(include_fr1){
-					transitionTable[NT1,NT2] = sum(tmp[,c(FR1, CDR1, FR2, CDR2, FR3)])
-				} else {
-					transitionTable[NT1,NT2] = sum(tmp[,c(CDR1, FR2, CDR2, FR3)])
-				}
-			}
-		}
-		transition = transitionTable
-		transition$id = names(transition)
-		
-		transition2 = melt(transition, id.vars="id")
-		
-		transition2 = merge(transition2, base.order, by.x="id", by.y="base")
-		transition2 = merge(transition2, base.order, by.x="variable", by.y="base")
-
-		transition2[is.na(transition2$value),]$value = 0
-		
-		if(!all(transition2$value == 0)){ #having rows of data but a transition table filled with 0 is bad
-
-			print("Plotting stacked transition")
-
-			png(filename=paste("transitions_stacked_", name, ".png", sep=""))
-			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
-			p = p + xlab("From base") + ylab("To base") + ggtitle("Mutations frequency from base to base") + guides(fill=guide_legend(title=NULL))
-			p = p + theme(panel.background = element_rect(fill = "white", colour="black")) + scale_fill_manual(values=c("A" = "blue4", "G" = "lightblue1", "C" = "olivedrab3", "T" = "olivedrab4"))
-			#p = p + scale_colour_manual(values=c("A" = "black", "G" = "black", "C" = "black", "T" = "black"))
-			print(p)
-			dev.off()
-
-			print("Plotting heatmap transition")
-
-			png(filename=paste("transitions_heatmap_", name, ".png", sep=""))
-			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
-			p = p + xlab("From base") + ylab("To base") + ggtitle("Mutations frequency from base to base")  + theme(panel.background = element_rect(fill = "white", colour="black"))
-			print(p)
-			dev.off()
-		} else {
-			print("No data to plot")
-		}
-	}
-
-	#print(paste("writing value file: ", name, "_", fname, "_value.txt" ,sep=""))
-
-	write.table(x=transitionTable, file=paste("transitions_", name ,"_", fname, ".txt", sep=""), sep=",",quote=F,row.names=T,col.names=NA)
-	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)
-
-	cat(matrx[1,x], file=paste(name, "_", fname, "_value.txt" ,sep=""))
-	cat(nrow(tmp), file=paste(name, "_", fname, "_n.txt" ,sep=""))
-
-	#print(paste(fname, name, nrow(tmp)))
-
-	matrx
-}
-
-nts = c("a", "c", "g", "t")
-zeros=rep(0, 4)
-
-funcs = c(median, sum, mean)
-fnames = c("median", "sum", "mean")
-
-print("Creating result tables")
-
-for(i in 1:length(funcs)){
-	func = funcs[[i]]
-	fname = fnames[[i]]
-	
-	rows = 9
-	if(fname == "sum"){
-		rows = 11
-	}
-	matrx = matrix(data = 0, ncol=((length(genes) + 1) * 3),nrow=rows)
-	
-	for(i in 1:length(genes)){
-		print(paste("Creating table for", fname, genes[i]))
-		matrx = calculate_result(i, genes[i], dat, matrx, func, fname, genes[i])
-	}
-
-	matrx = calculate_result(i + 1, ".*", dat[!grepl("unmatched", dat$best_match),], matrx, func, fname, name="all")
-	
-	result = data.frame(matrx)
-	if(fname == "sum"){
-		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")
-	} else {
-		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)")
-	}
-
-	write.table(x=result, file=paste("mutations_", fname, ".txt", sep=""), sep=",",quote=F,row.names=T,col.names=F)
-}
-
-print("Adding median number of mutations to sum table")
-
-sum.table = read.table("mutations_sum.txt", sep=",", header=F)
-median.table = read.table("mutations_median.txt", sep=",", header=F)
-
-new.table = sum.table[1,]
-new.table[2,] = median.table[1,]
-new.table[3:12,] = sum.table[2:11,]
-new.table[,1] = as.character(new.table[,1])
-new.table[2,1] = "Median of Number of Mutations (%)"
-
-#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"),]
-
-write.table(x=new.table, file="mutations_sum.txt", sep=",",quote=F,row.names=F,col.names=F)
-
-
-print("Plotting IGA piechart")
-
-dat = dat[!grepl("^unmatched", dat$best_match),]
-
-#blegh
-genesForPlot = dat[grepl("IGA", dat$best_match),]$best_match
-if(length(genesForPlot) > 0){
-	genesForPlot = data.frame(table(genesForPlot))
-	colnames(genesForPlot) = c("Gene","Freq")
-	genesForPlot$label = paste(genesForPlot$Gene, "-", genesForPlot$Freq)
-
-	pc = ggplot(genesForPlot, aes(x = factor(1), y=Freq, fill=Gene))
-	pc = pc + geom_bar(width = 1, stat = "identity") + scale_fill_manual(labels=genesForPlot$label, values=c("IGA1" = "lightblue1", "IGA2" = "blue4"))
-	pc = pc + coord_polar(theta="y")
-	pc = pc + theme(panel.background = element_rect(fill = "white", colour="black"))
-	pc = pc + xlab(" ") + ylab(" ") + ggtitle(paste("IGA subclasses", "( n =", sum(genesForPlot$Freq), ")"))
-	write.table(genesForPlot, "IGA.txt", sep="\t",quote=F,row.names=F,col.names=T)
-
-	png(filename="IGA.png")
-	print(pc)
-	dev.off()
-}
-
-print("Plotting IGG piechart")
-
-genesForPlot = dat[grepl("IGG", dat$best_match),]$best_match
-if(length(genesForPlot) > 0){
-	genesForPlot = data.frame(table(genesForPlot))
-	colnames(genesForPlot) = c("Gene","Freq")
-	genesForPlot$label = paste(genesForPlot$Gene, "-", genesForPlot$Freq)
-
-	pc = ggplot(genesForPlot, aes(x = factor(1), y=Freq, fill=Gene))
-	pc = pc + geom_bar(width = 1, stat = "identity") + scale_fill_manual(labels=genesForPlot$label, values=c("IGG1" = "olivedrab3", "IGG2" = "red", "IGG3" = "gold", "IGG4" = "darkred"))
-	pc = pc + coord_polar(theta="y") 
-	pc = pc + theme(panel.background = element_rect(fill = "white", colour="black"))
-	pc = pc + xlab(" ") + ylab(" ") + ggtitle(paste("IGG subclasses", "( n =", sum(genesForPlot$Freq), ")"))
-	write.table(genesForPlot, "IGG.txt", sep="\t",quote=F,row.names=F,col.names=T)
-
-	png(filename="IGG.png")
-	print(pc)
-	dev.off()
-}
-
-
-print("Plotting scatterplot")
-
-dat$percentage_mutations = round(dat$VRegionMutations / dat$VRegionNucleotides * 100, 2)
-
-p = ggplot(dat, aes(best_match, percentage_mutations))
-p = p + geom_point(aes(colour=best_match), position="jitter") + geom_boxplot(aes(middle=mean(percentage_mutations)), alpha=0.1, outlier.shape = NA)
-p = p + xlab("Subclass") + ylab("Frequency") + ggtitle("Frequency scatter plot") + theme(panel.background = element_rect(fill = "white", colour="black"))
-p = p + scale_fill_manual(values=c("IGA1" = "lightblue1", "IGA2" = "blue4", "IGG1" = "olivedrab3", "IGG2" = "red", "IGG3" = "gold", "IGG4" = "darkred", "IGM" = "black"))
-p = p + scale_colour_manual(values=c("IGA1" = "lightblue1", "IGA2" = "blue4", "IGG1" = "olivedrab3", "IGG2" = "red", "IGG3" = "gold", "IGG4" = "darkred", "IGM" = "black"))
-
-png(filename="scatter.png")
-print(p)
-dev.off()
-
-write.table(dat[,c("Sequence.ID", "best_match", "VRegionMutations", "VRegionNucleotides", "percentage_mutations")], "scatter.txt", sep="\t",quote=F,row.names=F,col.names=T)
-
-write.table(dat, input, sep="\t",quote=F,row.names=F,col.names=T)
-
-
-print("Plotting frequency ranges plot")
-
-dat$best_match_class = substr(dat$best_match, 0, 3)
-freq_labels = c("0", "0-2", "2-5", "5-10", "10-15", "15-20", "20")
-dat$frequency_bins = cut(dat$percentage_mutations, breaks=c(-Inf, 0, 2,5,10,15,20, Inf), labels=freq_labels)
-
-frequency_bins_sum = data.frame(data.table(dat)[, list(class_sum=sum(.N)), by=c("best_match_class")])
-
-frequency_bins_data = data.frame(data.table(dat)[, list(frequency_count=.N), by=c("best_match_class", "frequency_bins")])
-
-frequency_bins_data = merge(frequency_bins_data, frequency_bins_sum, by="best_match_class")
-
-frequency_bins_data$frequency = round(frequency_bins_data$frequency_count / frequency_bins_data$class_sum * 100, 2)
-
-p = ggplot(frequency_bins_data, aes(frequency_bins, frequency))
-p = p + geom_bar(aes(fill=best_match_class), stat="identity", position="dodge") + theme(panel.background = element_rect(fill = "white", colour="black"))
-p = p + xlab("Frequency ranges") + ylab("Frequency") + ggtitle("Mutation Frequencies by class") + scale_fill_manual(values=c("IGA" = "blue4", "IGG" = "olivedrab3", "IGM" = "black"))
-
-png(filename="frequency_ranges.png")
-print(p)
-dev.off()
-
-frequency_bins_data_by_class = frequency_bins_data
-
-write.table(frequency_bins_data_by_class, "frequency_ranges_classes.txt", sep="\t",quote=F,row.names=F,col.names=T)
-
-frequency_bins_data = data.frame(data.table(dat)[, list(frequency_count=.N), by=c("best_match", "best_match_class", "frequency_bins")])
-
-frequency_bins_data = merge(frequency_bins_data, frequency_bins_sum, by="best_match_class")
-
-frequency_bins_data$frequency = round(frequency_bins_data$frequency_count / frequency_bins_data$class_sum * 100, 2)
-
-write.table(frequency_bins_data, "frequency_ranges_subclasses.txt", sep="\t",quote=F,row.names=F,col.names=T)
-
-
-#frequency_bins_data_by_class
-#frequency_ranges_subclasses.txt
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+library(data.table)
+library(ggplot2)
+library(reshape2)
+
+args <- commandArgs(trailingOnly = TRUE)
+
+input = args[1]
+genes = unlist(strsplit(args[2], ","))
+outputdir = args[3]
+empty.region.filter = args[4]
+setwd(outputdir)
+
+dat = read.table(input, header=T, sep="\t", fill=T, stringsAsFactors=F)
+
+if(length(dat$Sequence.ID) == 0){
+  setwd(outputdir)
+  result = data.frame(x = rep(0, 5), y = rep(0, 5), z = rep(NA, 5))
+  row.names(result) = c("Number of Mutations (%)", "Transition (%)", "Transversions (%)", "Transitions at G C (%)", "Targeting of C G (%)")
+  write.table(x=result, file="mutations.txt", sep=",",quote=F,row.names=T,col.names=F)
+  transitionTable = data.frame(A=rep(0, 4),C=rep(0, 4),G=rep(0, 4),T=rep(0, 4))
+  row.names(transitionTable) = c("A", "C", "G", "T")
+  transitionTable["A","A"] = NA
+  transitionTable["C","C"] = NA
+  transitionTable["G","G"] = NA
+  transitionTable["T","T"] = NA
+
+  write.table(x=transitionTable, file="transitions.txt", sep=",",quote=F,row.names=T,col.names=NA)
+  cat("0", file="n.txt")
+  stop("No data")
+}
+
+cleanup_columns = c("FR1.IMGT.c.a",
+					"FR2.IMGT.g.t",
+					"CDR1.IMGT.Nb.of.nucleotides",
+					"CDR2.IMGT.t.a",
+					"FR1.IMGT.c.g",
+					"CDR1.IMGT.c.t",
+					"FR2.IMGT.a.c",
+					"FR2.IMGT.Nb.of.mutations",
+					"FR2.IMGT.g.c",
+					"FR2.IMGT.a.g",
+					"FR3.IMGT.t.a",
+					"FR3.IMGT.t.c",
+					"FR2.IMGT.g.a",
+					"FR3.IMGT.c.g",
+					"FR1.IMGT.Nb.of.mutations",
+					"CDR1.IMGT.g.a",
+					"CDR1.IMGT.t.g",
+					"CDR1.IMGT.g.c",
+					"CDR2.IMGT.Nb.of.nucleotides",
+					"FR2.IMGT.a.t",
+					"CDR1.IMGT.Nb.of.mutations",
+					"CDR3.IMGT.Nb.of.nucleotides",
+					"CDR1.IMGT.a.g",
+					"FR3.IMGT.a.c",
+					"FR1.IMGT.g.a",
+					"FR3.IMGT.a.g",
+					"FR1.IMGT.a.t",
+					"CDR2.IMGT.a.g",
+					"CDR2.IMGT.Nb.of.mutations",
+					"CDR2.IMGT.g.t",
+					"CDR2.IMGT.a.c",
+					"CDR1.IMGT.t.c",
+					"FR3.IMGT.g.c",
+					"FR1.IMGT.g.t",
+					"FR3.IMGT.g.t",
+					"CDR1.IMGT.a.t",
+					"FR1.IMGT.a.g",
+					"FR3.IMGT.a.t",
+					"FR3.IMGT.Nb.of.nucleotides",
+					"FR2.IMGT.t.c",
+					"CDR2.IMGT.g.a",
+					"FR2.IMGT.t.a",
+					"CDR1.IMGT.t.a",
+					"FR2.IMGT.t.g",
+					"FR3.IMGT.t.g",
+					"FR2.IMGT.Nb.of.nucleotides",
+					"FR1.IMGT.t.a",
+					"FR1.IMGT.t.g",
+					"FR3.IMGT.c.t",
+					"FR1.IMGT.t.c",
+					"CDR2.IMGT.a.t",
+					"FR2.IMGT.c.t",
+					"CDR1.IMGT.g.t",
+					"CDR2.IMGT.t.g",
+					"FR1.IMGT.Nb.of.nucleotides",
+					"CDR1.IMGT.c.g",
+					"CDR2.IMGT.t.c",
+					"FR3.IMGT.g.a",
+					"CDR1.IMGT.a.c",
+					"FR2.IMGT.c.a",
+					"FR3.IMGT.Nb.of.mutations",
+					"FR2.IMGT.c.g",
+					"CDR2.IMGT.g.c",
+					"FR1.IMGT.g.c",
+					"CDR2.IMGT.c.t",
+					"FR3.IMGT.c.a",
+					"CDR1.IMGT.c.a",
+					"CDR2.IMGT.c.g",
+					"CDR2.IMGT.c.a",
+					"FR1.IMGT.c.t",
+					"FR1.IMGT.Nb.of.silent.mutations",
+					"FR2.IMGT.Nb.of.silent.mutations",
+					"FR3.IMGT.Nb.of.silent.mutations",
+					"FR1.IMGT.Nb.of.nonsilent.mutations",
+					"FR2.IMGT.Nb.of.nonsilent.mutations",
+					"FR3.IMGT.Nb.of.nonsilent.mutations")
+
+print("Cleaning up columns")
+
+for(col in cleanup_columns){
+  dat[,col] = gsub("\\(.*\\)", "", dat[,col])
+  #dat[dat[,col] == "",] = "0"
+  dat[,col] = as.numeric(dat[,col])
+  dat[is.na(dat[,col]),col] = 0
+}
+
+regions = c("FR1", "CDR1", "FR2", "CDR2", "FR3")
+if(empty.region.filter == "FR1") {
+	regions = c("CDR1", "FR2", "CDR2", "FR3")
+} else if (empty.region.filter == "CDR1") {
+	regions = c("FR2", "CDR2", "FR3", "CDR3")
+} else if (empty.region.filter == "FR2") {
+	regions = c("CDR2", "FR3", "CDR3")
+}
+
+sum_by_row = function(x, columns) { sum(as.numeric(x[columns]), na.rm=T) }
+
+print("aggregating data into new columns")
+
+VRegionMutations_columns = paste(regions, ".IMGT.Nb.of.mutations", sep="")
+dat$VRegionMutations =  apply(dat, FUN=sum_by_row, 1, columns=VRegionMutations_columns)
+
+VRegionNucleotides_columns = paste(regions, ".IMGT.Nb.of.nucleotides", sep="")
+dat$FR3.IMGT.Nb.of.nucleotides = nchar(dat$FR3.IMGT.seq)
+dat$VRegionNucleotides =  apply(dat, FUN=sum_by_row, 1, columns=VRegionNucleotides_columns)
+
+transitionMutations_columns = paste(rep(regions, each=4), c(".IMGT.a.g", ".IMGT.g.a", ".IMGT.c.t", ".IMGT.t.c"), sep="")
+dat$transitionMutations = apply(dat, FUN=sum_by_row, 1, columns=transitionMutations_columns)
+
+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="")
+dat$transversionMutations = apply(dat, FUN=sum_by_row, 1, columns=transversionMutations_columns)
+
+transitionMutationsAtGC_columns = paste(rep(regions, each=2), c(".IMGT.g.a",".IMGT.c.t"), sep="")
+dat$transitionMutationsAtGC = apply(dat, FUN=sum_by_row, 1, columns=transitionMutationsAtGC_columns)
+
+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="")
+#totalMutationsAtGC_columns = paste(rep(regions, each=6), c(".IMGT.g.a",".IMGT.c.t",".IMGT.c.a",".IMGT.c.g",".IMGT.g.t"), sep="")
+dat$totalMutationsAtGC = apply(dat, FUN=sum_by_row, 1, columns=totalMutationsAtGC_columns)
+
+transitionMutationsAtAT_columns = paste(rep(regions, each=2), c(".IMGT.a.g",".IMGT.t.c"), sep="")
+dat$transitionMutationsAtAT = apply(dat, FUN=sum_by_row, 1, columns=transitionMutationsAtAT_columns)
+
+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="")
+#totalMutationsAtAT_columns = paste(rep(regions, each=5), c(".IMGT.a.g",".IMGT.t.c",".IMGT.a.c",".IMGT.g.c",".IMGT.t.g"), sep="")
+dat$totalMutationsAtAT = apply(dat, FUN=sum_by_row, 1, columns=totalMutationsAtAT_columns)
+
+FRRegions = regions[grepl("FR", regions)]
+CDRRegions = regions[grepl("CDR", regions)]
+
+FR_silentMutations_columns = paste(FRRegions, ".IMGT.Nb.of.silent.mutations", sep="")
+dat$silentMutationsFR = apply(dat, FUN=sum_by_row, 1, columns=FR_silentMutations_columns)
+
+CDR_silentMutations_columns = paste(CDRRegions, ".IMGT.Nb.of.silent.mutations", sep="")
+dat$silentMutationsCDR = apply(dat, FUN=sum_by_row, 1, columns=CDR_silentMutations_columns)
+
+FR_nonSilentMutations_columns = paste(FRRegions, ".IMGT.Nb.of.nonsilent.mutations", sep="")
+dat$nonSilentMutationsFR = apply(dat, FUN=sum_by_row, 1, columns=FR_nonSilentMutations_columns)
+
+CDR_nonSilentMutations_columns = paste(CDRRegions, ".IMGT.Nb.of.nonsilent.mutations", sep="")
+dat$nonSilentMutationsCDR = apply(dat, FUN=sum_by_row, 1, columns=CDR_nonSilentMutations_columns)
+
+mutation.sum.columns = c("Sequence.ID", "VRegionMutations", "VRegionNucleotides", "transitionMutations", "transversionMutations", "transitionMutationsAtGC", "transitionMutationsAtAT", "silentMutationsFR", "nonSilentMutationsFR", "silentMutationsCDR", "nonSilentMutationsCDR")
+write.table(dat[,mutation.sum.columns], "mutation_by_id.txt", sep="\t",quote=F,row.names=F,col.names=T)
+
+setwd(outputdir)
+
+base.order = data.frame(base=c("A", "T", "C", "G"), order=1:4)
+
+calculate_result = function(i, gene, dat, matrx, f, fname, name){
+	tmp = dat[grepl(paste("^", gene, ".*", sep=""), dat$best_match),]
+
+	j = i - 1
+	x = (j * 3) + 1
+	y = (j * 3) + 2
+	z = (j * 3) + 3
+
+	if(nrow(tmp) > 0){
+		if(fname == "sum"){
+			matrx[1,x] = round(f(tmp$VRegionMutations, na.rm=T), digits=1)
+			matrx[1,y] = round(f(tmp$VRegionNucleotides, na.rm=T), digits=1)
+			matrx[1,z] = round(f(matrx[1,x] / matrx[1,y]) * 100, digits=1)
+		} else {
+			matrx[1,x] = round(f(tmp$VRegionMutations, na.rm=T), digits=1)
+			matrx[1,y] = round(f(tmp$VRegionNucleotides, na.rm=T), digits=1)
+			matrx[1,z] = round(f(tmp$VRegionMutations / tmp$VRegionNucleotides) * 100, digits=1)
+		}
+
+		matrx[2,x] = round(f(tmp$transitionMutations, na.rm=T), digits=1)
+		matrx[2,y] = round(f(tmp$VRegionMutations, na.rm=T), digits=1)
+		matrx[2,z] = round(matrx[2,x] / matrx[2,y] * 100, digits=1)
+
+		matrx[3,x] = round(f(tmp$transversionMutations, na.rm=T), digits=1)
+		matrx[3,y] = round(f(tmp$VRegionMutations, na.rm=T), digits=1)
+		matrx[3,z] = round(matrx[3,x] / matrx[3,y] * 100, digits=1)
+
+		matrx[4,x] = round(f(tmp$transitionMutationsAtGC, na.rm=T), digits=1)
+		matrx[4,y] = round(f(tmp$totalMutationsAtGC, na.rm=T), digits=1)
+		matrx[4,z] = round(matrx[4,x] / matrx[4,y] * 100, digits=1)
+
+		matrx[5,x] = round(f(tmp$totalMutationsAtGC, na.rm=T), digits=1)
+		matrx[5,y] = round(f(tmp$VRegionMutations, na.rm=T), digits=1)
+		matrx[5,z] = round(matrx[5,x] / matrx[5,y] * 100, digits=1)
+
+		matrx[6,x] = round(f(tmp$transitionMutationsAtAT, na.rm=T), digits=1)
+		matrx[6,y] = round(f(tmp$totalMutationsAtAT, na.rm=T), digits=1)
+		matrx[6,z] = round(matrx[6,x] / matrx[6,y] * 100, digits=1)
+
+		matrx[7,x] = round(f(tmp$totalMutationsAtAT, na.rm=T), digits=1)
+		matrx[7,y] = round(f(tmp$VRegionMutations, na.rm=T), digits=1)
+		matrx[7,z] = round(matrx[7,x] / matrx[7,y] * 100, digits=1)
+
+		matrx[8,x] = round(f(tmp$nonSilentMutationsFR, na.rm=T), digits=1)
+		matrx[8,y] = round(f(tmp$silentMutationsFR, na.rm=T), digits=1)
+		matrx[8,z] = round(matrx[8,x] / matrx[8,y], digits=1)
+
+		matrx[9,x] = round(f(tmp$nonSilentMutationsCDR, na.rm=T), digits=1)
+		matrx[9,y] = round(f(tmp$silentMutationsCDR, na.rm=T), digits=1)
+		matrx[9,z] = round(matrx[9,x] / matrx[9,y], digits=1)
+
+		if(fname == "sum"){
+			matrx[10,x] = round(f(rowSums(tmp[,c("FR2.IMGT.Nb.of.nucleotides", "FR3.IMGT.Nb.of.nucleotides")], na.rm=T)), digits=1)
+			matrx[10,y] = round(f(tmp$VRegionNucleotides, na.rm=T), digits=1)
+			matrx[10,z] = round(matrx[10,x] / matrx[10,y] * 100, digits=1)
+
+			matrx[11,x] = round(f(rowSums(tmp[,c("CDR1.IMGT.Nb.of.nucleotides", "CDR2.IMGT.Nb.of.nucleotides")], na.rm=T)), digits=1)
+			matrx[11,y] = round(f(tmp$VRegionNucleotides, na.rm=T), digits=1)
+			matrx[11,z] = round(matrx[11,x] / matrx[11,y] * 100, digits=1)
+		}
+	}
+  
+	transitionTable = data.frame(A=zeros,C=zeros,G=zeros,T=zeros)
+	row.names(transitionTable) = c("A", "C", "G", "T")
+	transitionTable["A","A"] = NA
+	transitionTable["C","C"] = NA
+	transitionTable["G","G"] = NA
+	transitionTable["T","T"] = NA
+
+	if(nrow(tmp) > 0){
+		for(nt1 in nts){
+			for(nt2 in nts){
+				if(nt1 == nt2){
+					next
+				}
+				NT1 = LETTERS[letters == nt1]
+				NT2 = LETTERS[letters == nt2]
+				FR1 = paste("FR1.IMGT.", nt1, ".", nt2, sep="")
+				CDR1 = paste("CDR1.IMGT.", nt1, ".", nt2, sep="")
+				FR2 = paste("FR2.IMGT.", nt1, ".", nt2, sep="")
+				CDR2 = paste("CDR2.IMGT.", nt1, ".", nt2, sep="")
+				FR3 = paste("FR3.IMGT.", nt1, ".", nt2, sep="")
+				if (empty.region.filter == "leader"){
+					transitionTable[NT1,NT2] = sum(tmp[,c(FR1, CDR1, FR2, CDR2, FR3)])
+				} else if (empty.region.filter == "FR1") {
+					transitionTable[NT1,NT2] = sum(tmp[,c(CDR1, FR2, CDR2, FR3)])
+				} else if (empty.region.filter == "CDR1") {
+					transitionTable[NT1,NT2] = sum(tmp[,c(FR2, CDR2, FR3)])
+				} else if (empty.region.filter == "FR2") {
+					transitionTable[NT1,NT2] = sum(tmp[,c(CDR2, FR3)])
+				}
+			}
+		}
+		transition = transitionTable
+		transition$id = names(transition)
+		
+		transition2 = melt(transition, id.vars="id")
+
+		transition2 = merge(transition2, base.order, by.x="id", by.y="base")
+
+		transition2 = merge(transition2, base.order, by.x="variable", by.y="base")
+
+		transition2[is.na(transition2$value),]$value = 0
+
+		if(any(transition2$value == 0)){ #having rows of data but a transition table filled with 0 is bad
+			print("Plotting stacked transition")
+			png(filename=paste("transitions_stacked_", name, ".png", sep=""))
+			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
+			p = p + xlab("From base") + ylab("To base") + ggtitle("Mutations frequency from base to base") + guides(fill=guide_legend(title=NULL))
+			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"))
+			#p = p + scale_colour_manual(values=c("A" = "black", "G" = "black", "C" = "black", "T" = "black"))
+			print(p)
+			dev.off()
+
+			print("Plotting heatmap transition")
+
+			png(filename=paste("transitions_heatmap_", name, ".png", sep=""))
+			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
+			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"))
+			print(p)
+			dev.off()
+		} else {
+			print("No data to plot")
+		}
+	}
+
+	#print(paste("writing value file: ", name, "_", fname, "_value.txt" ,sep=""))
+	write.table(x=transitionTable, file=paste("transitions_", name ,"_", fname, ".txt", sep=""), sep=",",quote=F,row.names=T,col.names=NA)
+	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)
+	cat(matrx[1,x], file=paste(name, "_", fname, "_value.txt" ,sep=""))
+	cat(nrow(tmp), file=paste(name, "_", fname, "_n.txt" ,sep=""))
+	#print(paste(fname, name, nrow(tmp)))
+	matrx
+}
+nts = c("a", "c", "g", "t")
+zeros=rep(0, 4)
+funcs = c(median, sum, mean)
+fnames = c("median", "sum", "mean")
+
+print("Creating result tables")
+
+for(i in 1:length(funcs)){
+	func = funcs[[i]]
+	fname = fnames[[i]]
+	
+	rows = 9
+	if(fname == "sum"){
+		rows = 11
+	}
+	matrx = matrix(data = 0, ncol=((length(genes) + 1) * 3),nrow=rows)
+	for(i in 1:length(genes)){
+		print(paste("Creating table for", fname, genes[i]))
+		matrx = calculate_result(i, genes[i], dat, matrx, func, fname, genes[i])
+	}
+	matrx = calculate_result(i + 1, ".*", dat[!grepl("unmatched", dat$best_match),], matrx, func, fname, name="all")
+
+	result = data.frame(matrx)
+	if(fname == "sum"){
+		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")
+	} else {
+		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)")
+	}
+	write.table(x=result, file=paste("mutations_", fname, ".txt", sep=""), sep=",",quote=F,row.names=T,col.names=F)
+}
+
+print("Adding median number of mutations to sum table")
+sum.table = read.table("mutations_sum.txt", sep=",", header=F)
+median.table = read.table("mutations_median.txt", sep=",", header=F)
+
+new.table = sum.table[1,]
+new.table[2,] = median.table[1,]
+new.table[3:12,] = sum.table[2:11,]
+new.table[,1] = as.character(new.table[,1])
+new.table[2,1] = "Median of Number of Mutations (%)"
+
+#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"),]
+
+write.table(x=new.table, file="mutations_sum.txt", sep=",",quote=F,row.names=F,col.names=F)
+
+print("Plotting IGA piechart")
+
+dat = dat[!grepl("^unmatched", dat$best_match),]
+
+#blegh
+
+genesForPlot = dat[grepl("IGA", dat$best_match),]$best_match
+
+if(length(genesForPlot) > 0){
+	genesForPlot = data.frame(table(genesForPlot))
+	colnames(genesForPlot) = c("Gene","Freq")
+	genesForPlot$label = paste(genesForPlot$Gene, "-", genesForPlot$Freq)
+
+	pc = ggplot(genesForPlot, aes(x = factor(1), y=Freq, fill=Gene))
+	pc = pc + geom_bar(width = 1, stat = "identity") + scale_fill_manual(labels=genesForPlot$label, values=c("IGA1" = "lightblue1", "IGA2" = "blue4"))
+	pc = pc + coord_polar(theta="y") + scale_y_continuous(breaks=NULL)
+	pc = pc + theme(panel.background = element_rect(fill = "white", colour="black"), text = element_text(size=13, colour="black"))
+	pc = pc + xlab(" ") + ylab(" ") + ggtitle(paste("IGA subclasses", "( n =", sum(genesForPlot$Freq), ")"))
+	write.table(genesForPlot, "IGA.txt", sep="\t",quote=F,row.names=F,col.names=T)
+
+	png(filename="IGA.png")
+	print(pc)
+	dev.off()
+}
+
+print("Plotting IGG piechart")
+
+genesForPlot = dat[grepl("IGG", dat$best_match),]$best_match
+
+if(length(genesForPlot) > 0){
+	genesForPlot = data.frame(table(genesForPlot))
+	colnames(genesForPlot) = c("Gene","Freq")
+	genesForPlot$label = paste(genesForPlot$Gene, "-", genesForPlot$Freq)
+
+	pc = ggplot(genesForPlot, aes(x = factor(1), y=Freq, fill=Gene))
+	pc = pc + geom_bar(width = 1, stat = "identity") + scale_fill_manual(labels=genesForPlot$label, values=c("IGG1" = "olivedrab3", "IGG2" = "red", "IGG3" = "gold", "IGG4" = "darkred"))
+	pc = pc + coord_polar(theta="y") + scale_y_continuous(breaks=NULL)
+	pc = pc + theme(panel.background = element_rect(fill = "white", colour="black"), text = element_text(size=13, colour="black"))
+	pc = pc + xlab(" ") + ylab(" ") + ggtitle(paste("IGG subclasses", "( n =", sum(genesForPlot$Freq), ")"))
+	write.table(genesForPlot, "IGG.txt", sep="\t",quote=F,row.names=F,col.names=T)
+
+	png(filename="IGG.png")
+	print(pc)
+	dev.off()
+}
+
+print("Plotting scatterplot")
+
+dat$percentage_mutations = round(dat$VRegionMutations / dat$VRegionNucleotides * 100, 2)
+dat.clss = dat
+
+dat.clss$best_match = substr(dat.clss$best_match, 0, 3)
+
+dat.clss = rbind(dat, dat.clss)
+
+p = ggplot(dat.clss, aes(best_match, percentage_mutations))
+p = p + geom_point(aes(colour=best_match), position="jitter") + geom_boxplot(aes(middle=mean(percentage_mutations)), alpha=0.1, outlier.shape = NA)
+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"))
+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"))
+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"))
+
+png(filename="scatter.png")
+print(p)
+dev.off()
+
+write.table(dat[,c("Sequence.ID", "best_match", "VRegionMutations", "VRegionNucleotides", "percentage_mutations")], "scatter.txt", sep="\t",quote=F,row.names=F,col.names=T)
+
+write.table(dat, input, sep="\t",quote=F,row.names=F,col.names=T)
+
+print("Plotting frequency ranges plot")
+
+dat$best_match_class = substr(dat$best_match, 0, 3)
+freq_labels = c("0", "0-2", "2-5", "5-10", "10-15", "15-20", "20")
+dat$frequency_bins = cut(dat$percentage_mutations, breaks=c(-Inf, 0, 2,5,10,15,20, Inf), labels=freq_labels)
+
+frequency_bins_sum = data.frame(data.table(dat)[, list(class_sum=sum(.N)), by=c("best_match_class")])
+
+frequency_bins_data = data.frame(data.table(dat)[, list(frequency_count=.N), by=c("best_match_class", "frequency_bins")])
+
+frequency_bins_data = merge(frequency_bins_data, frequency_bins_sum, by="best_match_class")
+
+frequency_bins_data$frequency = round(frequency_bins_data$frequency_count / frequency_bins_data$class_sum * 100, 2)
+
+p = ggplot(frequency_bins_data, aes(frequency_bins, frequency))
+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"))
+p = p + xlab("Frequency ranges") + ylab("Frequency") + ggtitle("Mutation Frequencies by class") + scale_fill_manual(values=c("IGA" = "blue4", "IGG" = "olivedrab3", "IGM" = "black"))
+
+png(filename="frequency_ranges.png")
+print(p)
+dev.off()
+
+frequency_bins_data_by_class = frequency_bins_data
+
+write.table(frequency_bins_data_by_class, "frequency_ranges_classes.txt", sep="\t",quote=F,row.names=F,col.names=T)
+
+frequency_bins_data = data.frame(data.table(dat)[, list(frequency_count=.N), by=c("best_match", "best_match_class", "frequency_bins")])
+
+frequency_bins_data = merge(frequency_bins_data, frequency_bins_sum, by="best_match_class")
+
+frequency_bins_data$frequency = round(frequency_bins_data$frequency_count / frequency_bins_data$class_sum * 100, 2)
+
+write.table(frequency_bins_data, "frequency_ranges_subclasses.txt", sep="\t",quote=F,row.names=F,col.names=T)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
--- a/shm_csr.xml	Thu Oct 13 10:52:24 2016 -0400
+++ b/shm_csr.xml	Tue Oct 25 07:28:43 2016 -0400
@@ -1,24 +1,20 @@
 <tool id="shm_csr" name="SHM &amp; CSR pipeline" version="1.0">
 	<description></description>
 	<command interpreter="bash">
-		wrapper.sh $in_file custom $out_file $out_file.files_path ${in_file.name} ${include_fr1} $functionality $unique $naive_output_ca $naive_output_cg $naive_output_cm $filter_uniques $class_filter $empty_region_filter
+		wrapper.sh $in_file custom $out_file $out_file.files_path ${in_file.name} "-" $functionality $unique $naive_output_ca $naive_output_cg $naive_output_cm $filter_uniques $class_filter $empty_region_filter
 	</command>
 	<inputs>
 		<param name="in_file" type="data" label="IMGT zip file to be analysed" />
 		<param name="empty_region_filter" type="select" label="Sequence starts at" help="" >
+			<option value="leader" selected="true">Leader: include FR1, CDR1, FR2, CDR2, FR3 in filters</option>
 			<option value="FR1" selected="true">FR1: include CDR1,FR2,CDR2,FR3 in filters</option>
 			<option value="CDR1">CDR1: include FR2,CDR2,FR3 in filters</option>
 			<option value="FR2">FR2: include CDR2,FR3 in filters</option>
 		</param>
-		<param name="include_fr1" type="select" label="Include mutations in FR1 region" help="" >
-			<option value="yes">yes</option>
-			<option value="no" selected="true">no</option>
-		</param>
 		<param name="functionality" type="select" label="Functionality filter" help="" >
-			<option value="productive" selected="true">Productive: Keep "productive" and "productive (see comment)"</option>
-			<option value="unproductive">Unproductive: Keep "unproductive" and "unproductive (see comment)"</option>
-			<option value="remove_unknown">Remove "unknown" and "unknown (see comment)"</option>
-			<option value="dont_filter">Don't filter</option>
+			<option value="productive" selected="true">Productive (Productive and Productive see comment)</option>
+			<option value="unproductive">Unproductive (Unproductive and Unproductive see comment)</option>
+			<option value="remove_unknown">Productive and Unproductive (Productive, Productive see comment, Unproductive, Unproductive and Unproductive see comment)</option>
 		</param>
 		<param name="filter_uniques" type="select" label="Filter unique sequences" help="See below for an example.">
 			<option value="remove">Remove uniques (Based on nucleotide sequence + C)</option>
@@ -42,7 +38,7 @@
 			<option value="60_55">>60% class and >55% subclass</option>
 			<option value="70_0">>70% class</option>
 			<option value="60_0">>60% class</option>
-			<option value="101_101">No class</option>
+			<option value="101_101">Do not assign (sub)class</option>
 		</param>
 		<conditional name="naive_output_cond">
 			<param name="naive_output" type="select" label="Output new IMGT archives per class into your history?">
--- a/wrapper.sh	Thu Oct 13 10:52:24 2016 -0400
+++ b/wrapper.sh	Tue Oct 25 07:28:43 2016 -0400
@@ -158,13 +158,13 @@
 
 classes="IGA,IGA1,IGA2,IGG,IGG1,IGG2,IGG3,IGG4,IGM,unmatched"
 echo "R mutation analysis"
-Rscript $dir/shm_csr.r $outdir/merged.txt $classes $outdir ${include_fr1} 2>&1
+Rscript $dir/shm_csr.r $outdir/merged.txt $classes $outdir ${empty_region_filter} 2>&1
 
 
 echo "---------------- shm_csr.py ----------------"
 echo "---------------- shm_csr.py ----------------<br />" >> $log
 
-python $dir/shm_csr.py --input $outdir/merged.txt --genes $classes --includefr1 "${include_fr1}" --output $outdir/hotspot_analysis.txt
+python $dir/shm_csr.py --input $outdir/merged.txt --genes $classes --empty_region_filter "${empty_region_filter}" --output $outdir/hotspot_analysis.txt
 
 echo "---------------- aa_histogram.r ----------------"
 echo "---------------- aa_histogram.r ----------------<br />" >> $log
@@ -575,6 +575,9 @@
 
 cd $tmp
 
+echo "Cleaning up *.RData files"
+find $outdir/baseline -name "*.RData" -type f -delete
+
 echo "---------------- naive_output.r ----------------"
 echo "---------------- naive_output.r ----------------<br />" >> $log