Mercurial > repos > davidvanzessen > mutation_analysis

changeset 65:ae8b721a2964 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author davidvanzessen
date Mon, 04 Apr 2016 04:25:46 -0400
parents 0fdd90f7c654
children 88e0e7665086
files merge_and_filter.r mutation_analysis.py.bak mutation_analysis.r.bak
diffstat 3 files changed, 3 insertions(+), 721 deletions(-) [+]
line wrap: on
line diff
--- a/merge_and_filter.r	Fri Apr 01 08:54:24 2016 -0400
+++ b/merge_and_filter.r	Mon Apr 04 04:25:46 2016 -0400
@@ -103,14 +103,13 @@
 	} else {
 		result$unique.def = paste(result$CDR1.IMGT, result$CDR2.IMGT, result$CDR3.IMGT, result$FR2.IMGT, result$FR3.IMGT)
 	}
-	result.filtered = result[duplicated(result$unique.def) | duplicated(result$unique.def, fromLast=T),]
-	fltr = result$Sequence.ID %in% result.filtered$Sequence.ID
+
 	#fltr = result$unique.def %in% result.filtered$unique.def
 		
 	if(grepl("keep", filter_unique)){
-		result = result[!fltr,]
+		result = result[!duplicated(result$unique.def),]
 	} else {
-		result = result[fltr,]
+		result = result[duplicated(result$unique.def) | duplicated(result$unique.def, fromLast=T),]
 		result = result[!duplicated(result$unique.def),]
 	}
 	
--- a/mutation_analysis.py.bak	Fri Apr 01 08:54:24 2016 -0400
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,248 +0,0 @@
-from collections import defaultdict
-import re
-import argparse
-
-parser = argparse.ArgumentParser()
-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="The genes available in the 'best_match' column")
-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
-outfile = args.output
-
-genedic = dict()
-
-mutationdic = dict()
-mutationMatcher = re.compile("^(.)(\d+).(.),?(.)?(\d+)?.?(.)?(.?.?.?.?.?)?")
-linecount = 0
-
-IDIndex = 0
-best_matchIndex = 0
-fr1Index = 0
-cdr1Index = 0
-fr2Index = 0
-cdr2Index = 0
-fr3Index = 0
-first = True
-IDlist = []
-mutationList = []
-mutationListByID = {}
-cdr1LengthDic = {}
-cdr2LengthDic = {}
-
-with open(infile, 'r') as i:
-	for line in i:
-		if first:
-			linesplt = line.split("\t")
-			IDIndex = linesplt.index("Sequence.ID")
-			best_matchIndex = linesplt.index("best_match")
-			fr1Index = linesplt.index("FR1.IMGT")
-			cdr1Index = linesplt.index("CDR1.IMGT")
-			fr2Index = linesplt.index("FR2.IMGT")
-			cdr2Index = linesplt.index("CDR2.IMGT")
-			fr3Index = linesplt.index("FR3.IMGT")
-			cdr1LengthIndex = linesplt.index("CDR1.IMGT.Nb.of.nucleotides")
-			cdr2LengthIndex = linesplt.index("CDR2.IMGT.Nb.of.nucleotides")
-			first = False
-			continue
-		linecount += 1
-		linesplt = line.split("\t")
-		ID = linesplt[IDIndex]
-		genedic[ID] = linesplt[best_matchIndex]
-		mutationdic[ID + "_FR1"] = [mutationMatcher.match(x).groups() for x in linesplt[fr1Index].split("|") if x] if include_fr1 else []
-		mutationdic[ID + "_CDR1"] = [mutationMatcher.match(x).groups() for x in linesplt[cdr1Index].split("|") if x]
-		mutationdic[ID + "_FR2"] = [mutationMatcher.match(x).groups() for x in linesplt[fr2Index].split("|") if x]
-		mutationdic[ID + "_CDR2"] = [mutationMatcher.match(x).groups() for x in linesplt[cdr2Index].split("|") if x]
-		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]
-
-		mutationList += mutationdic[ID + "_FR1"] + mutationdic[ID + "_CDR1"] + mutationdic[ID + "_FR2"] + mutationdic[ID + "_CDR2"] + mutationdic[ID + "_FR3"]
-		mutationListByID[ID] = mutationdic[ID + "_FR1"] + mutationdic[ID + "_CDR1"] + mutationdic[ID + "_FR2"] + mutationdic[ID + "_CDR2"] + mutationdic[ID + "_FR3"]
-
-		cdr1Length = linesplt[cdr1LengthIndex]
-		cdr2Length = linesplt[cdr2LengthIndex]
-
-		cdr1LengthDic[ID] = int(cdr1Length) / 3
-		cdr2LengthDic[ID] = int(cdr2Length) / 3
-
-		IDlist += [ID]
-
-AALength = (int(max(mutationList, key=lambda i: int(i[4]) if i[4] else 0)[4]) + 1)  # [4] is the position of the AA mutation, None if silent
-
-AA_mutation = [0] * AALength
-AA_mutation_empty = AA_mutation[:]
-
-aa_mutations_by_id_file = outfile[:outfile.rindex("/")] + "/aa_id_mutations.txt"
-with open(aa_mutations_by_id_file, 'w') as o:
-	for ID in mutationListByID.keys():
-		AA_mutation_for_ID = AA_mutation_empty[:]
-		for mutation in mutationListByID[ID]:
-			if mutation[4]:
-				AA_mutation[int(mutation[4])] += 1
-				AA_mutation_for_ID[int(mutation[4])] += 1
-		o.write(ID + "\t" + "\t".join([str(x) for x in AA_mutation_for_ID[1:]]) + "\n")
-
-
-
-#absent AA stuff
-absentAACDR1Dic = defaultdict(list)
-absentAACDR1Dic[5] = range(29,36)
-absentAACDR1Dic[6] = range(29,35)
-absentAACDR1Dic[7] = range(30,35)
-absentAACDR1Dic[8] = range(30,34)
-absentAACDR1Dic[9] = range(31,34)
-absentAACDR1Dic[10] = range(31,33)
-absentAACDR1Dic[11] = [32]
-
-absentAACDR2Dic = defaultdict(list)
-absentAACDR2Dic[0] = range(55,65)
-absentAACDR2Dic[1] = range(56,65)
-absentAACDR2Dic[2] = range(56,64)
-absentAACDR2Dic[3] = range(57,64)
-absentAACDR2Dic[4] = range(57,63)
-absentAACDR2Dic[5] = range(58,63)
-absentAACDR2Dic[6] = range(58,62)
-absentAACDR2Dic[7] = range(59,62)
-absentAACDR2Dic[8] = range(59,61)
-absentAACDR2Dic[9] = [60]
-
-absentAA = [len(IDlist)] * (AALength-1)
-for k, cdr1Length in cdr1LengthDic.iteritems():
-	for c in absentAACDR1Dic[cdr1Length]:
-		absentAA[c] -= 1
-
-for k, cdr2Length in cdr2LengthDic.iteritems():
-	for c in absentAACDR2Dic[cdr2Length]:
-		absentAA[c] -= 1
-
-
-aa_mutations_by_id_file = outfile[:outfile.rindex("/")] + "/absent_aa_id.txt"
-with open(aa_mutations_by_id_file, 'w') as o:
-	o.write("ID\tcdr1length\tcdr2length\t" + "\t".join([str(x) for x in range(1,AALength-1)]) + "\n")
-	for ID in IDlist:
-		absentAAbyID = [1] * (AALength-1)
-		cdr1Length = cdr1LengthDic[ID]
-		for c in absentAACDR1Dic[cdr1Length]:
-			absentAAbyID[c] -= 1
-
-		cdr2Length = cdr2LengthDic[ID]
-		for c in absentAACDR2Dic[cdr2Length]:
-			absentAAbyID[c] -= 1
-		o.write(ID + "\t" + str(cdr1Length) + "\t" + str(cdr2Length) + "\t" + "\t".join([str(x) for x in absentAAbyID]) + "\n")
-
-
-
-aa_mutations_file = outfile[:outfile.rindex("/")] + "/aa_mutations.txt"
-with open(aa_mutations_file, 'w') as o:
-	o.write("row.name\t" + "\t".join([str(x) for x in range(1, AALength-1)]) + "\n")
-	o.write("mutations.at.position\t" + "\t".join([str(x) for x in AA_mutation[1:]]) + "\n")
-	o.write("AA.at.position\t" + "\t".join([str(x) for x in absentAA]) + "\n")
-
-if linecount == 0:
-	print "No data, exiting"
-	with open(outfile, 'w') as o:
-		o.write("RGYW (%)," + ("0,0,0\n" * len(genes)))
-		o.write("WRCY (%)," + ("0,0,0\n" * len(genes)))
-		o.write("WA (%)," + ("0,0,0\n" * len(genes)))
-		o.write("TW (%)," + ("0,0,0\n" * len(genes)))
-	import sys
-
-	sys.exit()
-
-hotspotMatcher = re.compile("[actg]+,(\d+)-(\d+)\((.*)\)")
-RGYWCount = {g: 0 for g in genes}
-WRCYCount = {g: 0 for g in genes}
-WACount = {g: 0 for g in genes}
-TWCount = {g: 0 for g in genes}
-
-IDIndex = 0
-ataIndex = 0
-tatIndex = 0
-aggctatIndex = 0
-atagcctIndex = 0
-first = True
-with open(infile, 'r') as i:
-	for line in i:
-		if first:
-			linesplt = line.split("\t")
-			ataIndex = linesplt.index("X.a.t.a")
-			tatIndex = linesplt.index("t.a.t.")
-			aggctatIndex = linesplt.index("X.a.g.g.c.t..a.t.")
-			atagcctIndex = linesplt.index("X.a.t..a.g.c.c.t.")
-			first = False
-			continue
-		linesplt = line.split("\t")
-		gene = linesplt[best_matchIndex]
-		ID = linesplt[IDIndex]
-		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
-				[hotspotMatcher.match(x).groups() for x in linesplt[atagcctIndex].split("|") if x]]
-		WA = [(int(x), int(y), z) for (x, y, z) in
-			  [hotspotMatcher.match(x).groups() for x in linesplt[ataIndex].split("|") if x]]
-		TW = [(int(x), int(y), z) for (x, y, z) in
-			  [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"]
-		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])
-			mutation_in_WRCY = any([(start <= int(where) <= end) for (start, end, region) in WRCY])
-			mutation_in_WA = any([(start <= int(where) <= end) for (start, end, region) in WA])
-			mutation_in_TW = any([(start <= int(where) <= end) for (start, end, region) in TW])
-
-			in_how_many_motifs = sum([mutation_in_RGYW, mutation_in_WRCY, mutation_in_WA, mutation_in_TW])
-
-			if in_how_many_motifs > 0:
-				RGYWCount[ID] += (1.0 * int(mutation_in_RGYW)) / in_how_many_motifs
-				WRCYCount[ID] += (1.0 * int(mutation_in_WRCY)) / in_how_many_motifs
-				WACount[ID] += (1.0 * int(mutation_in_WA)) / in_how_many_motifs
-				TWCount[ID] += (1.0 * int(mutation_in_TW)) / in_how_many_motifs
-
-directory = outfile[:outfile.rfind("/") + 1]
-value = 0
-valuedic = dict()
-for gene in genes:
-	with open(directory + gene + "_value.txt", 'r') as v:
-		valuedic[gene] = int(v.readlines()[0].rstrip())
-with open(directory + "total_value.txt", 'r') as v:
-	valuedic["total"] = int(v.readlines()[0].rstrip())
-
-dic = {"RGYW": RGYWCount, "WRCY": WRCYCount, "WA": WACount, "TW": TWCount}
-arr = ["RGYW", "WRCY", "WA", "TW"]
-with open(outfile, 'w') as o:
-	for typ in arr:
-		o.write(typ + " (%)")
-		curr = dic[typ]
-		for gene in genes:
-			geneMatcher = re.compile(".*" + gene + ".*")
-			if valuedic[gene] is 0:
-				o.write(",0,0,0")
-			else:
-				x = int(round(sum([curr[x] for x in [y for y, z in genedic.iteritems() if geneMatcher.match(z)]])))
-				y = valuedic[gene]
-				z = str(round(x / float(valuedic[gene]) * 100, 1))
-				o.write("," + str(x) + "," + str(y) + "," + z)
-		# for total
-		x = int(round(sum([y for x, y in curr.iteritems()])))
-		y = valuedic["total"]
-		z = str(round(x / float(valuedic["total"]) * 100, 1))
-		o.write("," + str(x) + "," + str(y) + "," + z + "\n")
-
-
-# for testing
-seq_motif_file = outfile[:outfile.rindex("/")] + "/motif_per_seq.txt"
-with open(seq_motif_file, 'w') as o:
-	o.write("ID\tRGYWC\tWRCY\tWA\tTW\n")
-	for ID in IDlist:
-		o.write(ID + "\t" + str(round(RGYWCount[ID], 2)) + "\t" + str(round(WRCYCount[ID], 2)) + "\t" + str(
-			round(WACount[ID], 2)) + "\t" + str(round(TWCount[ID], 2)) + "\n")
--- a/mutation_analysis.r.bak	Fri Apr 01 08:54:24 2016 -0400
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,469 +0,0 @@
-library(data.table)
-library(ggplot2)
-
-args <- commandArgs(trailingOnly = TRUE)
-
-input = args[1]
-genes = unlist(strsplit(args[2], ","))
-outputdir = args[3]
-print(args[4])
-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",
-                    "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")
-
-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]),] = 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) }
-
-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$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)
-
-nts = c("a", "c", "g", "t")
-zeros=rep(0, 4)
-matrx = matrix(data = 0, ncol=((length(genes) + 1) * 3),nrow=9)
-for(i in 1:length(genes)){
-  gene = genes[i]
-  tmp = dat[grepl(paste(".*", gene, ".*", sep=""), dat$best_match),]
-  if(gene == "."){
-    tmp = dat
-  }
-  j = i - 1
-  x = (j * 3) + 1
-  y = (j * 3) + 2
-  z = (j * 3) + 3
-  matrx[1,x] = sum(tmp$VRegionMutations)
-  matrx[1,y] = sum(tmp$VRegionNucleotides)
-  matrx[1,z] = round(matrx[1,x] / matrx[1,y] * 100, digits=1)
-
-  matrx[2,x] = sum(tmp$transitionMutations)
-  matrx[2,y] = sum(tmp$VRegionMutations)
-  matrx[2,z] = round(matrx[2,x] / matrx[2,y] * 100, digits=1)
-
-  matrx[3,x] = sum(tmp$transversionMutations)
-  matrx[3,y] = sum(tmp$VRegionMutations)
-  matrx[3,z] = round(matrx[3,x] / matrx[3,y] * 100, digits=1)
-
-  matrx[4,x] = sum(tmp$transitionMutationsAtGC)
-  matrx[4,y] = sum(tmp$totalMutationsAtGC)
-  matrx[4,z] = round(matrx[4,x] / matrx[4,y] * 100, digits=1)
-
-  matrx[5,x] = sum(tmp$totalMutationsAtGC)
-  matrx[5,y] = sum(tmp$VRegionMutations)
-  matrx[5,z] = round(matrx[5,x] / matrx[5,y] * 100, digits=1)
-
-  matrx[6,x] = sum(tmp$transitionMutationsAtAT)
-  matrx[6,y] = sum(tmp$totalMutationsAtAT)
-  matrx[6,z] = round(matrx[6,x] / matrx[6,y] * 100, digits=1)
-
-  matrx[7,x] = sum(tmp$totalMutationsAtAT)
-  matrx[7,y] = sum(tmp$VRegionMutations)
-  matrx[7,z] = round(matrx[7,x] / matrx[7,y] * 100, digits=1)
-
-  matrx[8,x] = sum(tmp$nonSilentMutationsFR)
-  matrx[8,y] = sum(tmp$silentMutationsFR)
-  matrx[8,z] = round(matrx[8,x] / matrx[8,y], digits=1)
-
-  matrx[9,x] = sum(tmp$nonSilentMutationsCDR)
-  matrx[9,y] = sum(tmp$silentMutationsCDR)
-  matrx[9,z] = round(matrx[9,x] / matrx[9,y], 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)])
-				}
-			}
-		}
-  }
-  
-  
-  write.table(x=transitionTable, file=paste("transitions_", gene ,".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_", gene ,".txt", sep=""), sep="\t",quote=F,row.names=F,col.names=T)
-  
-  cat(matrx[1,x], file=paste(gene, "_value.txt" ,sep=""))
-  cat(length(tmp$Sequence.ID), file=paste(gene, "_n.txt" ,sep=""))
-}
-
-#again for all of the data
-tmp = dat
-j = i
-x = (j * 3) + 1
-y = (j * 3) + 2
-z = (j * 3) + 3
-matrx[1,x] = sum(tmp$VRegionMutations)
-matrx[1,y] = sum(tmp$VRegionNucleotides)
-matrx[1,z] = round(matrx[1,x] / matrx[1,y] * 100, digits=1)
-
-matrx[2,x] = sum(tmp$transitionMutations)
-matrx[2,y] = sum(tmp$VRegionMutations)
-matrx[2,z] = round(matrx[2,x] / matrx[2,y] * 100, digits=1)
-
-matrx[3,x] = sum(tmp$transversionMutations)
-matrx[3,y] = sum(tmp$VRegionMutations)
-matrx[3,z] = round(matrx[3,x] / matrx[3,y] * 100, digits=1)
-
-matrx[4,x] = sum(tmp$transitionMutationsAtGC)
-matrx[4,y] = sum(tmp$totalMutationsAtGC)
-matrx[4,z] = round(matrx[4,x] / matrx[4,y] * 100, digits=1)
-
-matrx[5,x] = sum(tmp$totalMutationsAtGC)
-matrx[5,y] = sum(tmp$VRegionMutations)
-matrx[5,z] = round(matrx[5,x] / matrx[5,y] * 100, digits=1)
-
-matrx[6,x] = sum(tmp$transitionMutationsAtAT)
-matrx[6,y] = sum(tmp$totalMutationsAtAT)
-matrx[6,z] = round(matrx[6,x] / matrx[6,y] * 100, digits=1)
-
-matrx[7,x] = sum(tmp$totalMutationsAtAT)
-matrx[7,y] = sum(tmp$VRegionMutations)
-matrx[7,z] = round(matrx[7,x] / matrx[7,y] * 100, digits=1)
-
-matrx[8,x] = sum(tmp$nonSilentMutationsFR)
-matrx[8,y] = sum(tmp$silentMutationsFR)
-matrx[8,z] = round(matrx[8,x] / matrx[8,y], digits=1)
-
-matrx[9,x] = sum(tmp$nonSilentMutationsCDR)
-matrx[9,y] = sum(tmp$silentMutationsCDR)
-matrx[9,z] = round(matrx[9,x] / matrx[9,y], digits=1)
-
-transitionTable = data.frame(A=1:4,C=1:4,G=1:4,T=1: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
-
-
-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)])
-		}
-	}
-}
-write.table(x=transitionTable, file="transitions.txt", 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="matched_all.txt", sep="\t",quote=F,row.names=F,col.names=T)
-cat(matrx[1,x], file="total_value.txt")
-cat(length(tmp$Sequence.ID), file="total_n.txt")
-
-
-
-result = data.frame(matrx)
-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)")
-
-write.table(x=result, file="mutations.txt", sep=",",quote=F,row.names=T,col.names=F)
-
-
-if (!("ggplot2" %in% rownames(installed.packages()))) {
-	install.packages("ggplot2", repos="http://cran.xl-mirror.nl/") 
-}
-
-
-genesForPlot = gsub("[0-9]", "", dat$best_match)
-genesForPlot = data.frame(table(genesForPlot))
-colnames(genesForPlot) = c("Gene","Freq")
-genesForPlot$label = paste(genesForPlot$Gene, "-", genesForPlot$Freq)
-write.table(genesForPlot, "all.txt", sep="\t",quote=F,row.names=F,col.names=T)
-
-
-pc = ggplot(genesForPlot, aes(x = factor(1), y=Freq, fill=label))
-pc = pc + geom_bar(width = 1, stat = "identity")
-pc = pc + coord_polar(theta="y")
-pc = pc + xlab(" ") + ylab(" ") + ggtitle(paste("Classes", "( n =", sum(genesForPlot$Freq), ")"))
-
-png(filename="all.png")
-pc
-dev.off()
-
-
-#blegh
-genesForPlot = dat[grepl("ca", 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=label))
-	pc = pc + geom_bar(width = 1, stat = "identity")
-	pc = pc + coord_polar(theta="y")
-	pc = pc + xlab(" ") + ylab(" ") + ggtitle(paste("IgA subclasses", "( n =", sum(genesForPlot$Freq), ")"))
-	write.table(genesForPlot, "ca.txt", sep="\t",quote=F,row.names=F,col.names=T)
-
-	png(filename="ca.png")
-	print(pc)
-	dev.off()
-}
-
-genesForPlot = dat[grepl("cg", 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=label))
-	pc = pc + geom_bar(width = 1, stat = "identity")
-	pc = pc + coord_polar(theta="y")
-	pc = pc + xlab(" ") + ylab(" ") + ggtitle(paste("IgG subclasses", "( n =", sum(genesForPlot$Freq), ")"))
-	write.table(genesForPlot, "cg.txt", sep="\t",quote=F,row.names=F,col.names=T)
-
-	png(filename="cg.png")
-	print(pc)
-	dev.off()
-}
-
-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")
-
-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)
-
-
-
-
-
-
-dat$best_match_class = substr(dat$best_match, 0, 2)
-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_data = data.frame(data.table(dat)[, list(frequency_count=.N), by=c("best_match_class", "frequency_bins")])
-
-p = ggplot(frequency_bins_data, aes(frequency_bins, frequency_count))
-p = p + geom_bar(aes(fill=best_match_class), stat="identity", position="dodge")
-p = p + xlab("Frequency ranges") + ylab("Frequency") + ggtitle("Mutation Frequencies by class")
-
-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", "frequency_bins")])
-
-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
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-