--- a/mutation_analysis.py	Tue Aug 04 08:55:14 2015 -0400
+++ b/mutation_analysis.py	Mon Nov 02 04:52:26 2015 -0500
@@ -2,7 +2,8 @@
 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("--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")
@@ -13,7 +14,7 @@
 genes = str(args.genes).split(",")
 print "includefr1 =", args.includefr1
 include_fr1 = True if args.includefr1 == "yes" else False
-outfile = args.output 
+outfile = args.output
 
 genedic = dict()
 
@@ -28,58 +29,61 @@
 fr2Index = 0
 cdr2Index = 0
 fr3Index = 0
-first=True
+first = True
 IDlist = []
 mutationList = []
 
 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")
-			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"]
-		
-		IDlist += [ID]
+    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")
+            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"]
 
-AA_mutation = [0] * (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
+        IDlist += [ID]
+
+AA_mutation = [0] * (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
 
 for mutation in mutationList:
-	if mutation[4]: #if non silent mutation
-		AA_mutation[int(mutation[4])] += 1
+    if mutation[4]:  # if non silent mutation
+        AA_mutation[int(mutation[4])] += 1
 
 aa_mutations_file = outfile[:outfile.rindex("/")] + "/aa_mutations.txt"
 with open(aa_mutations_file, 'w') as o:
-	o.write(",".join([str(x) for x in AA_mutation]) + "\n")
-    
+    o.write(",".join([str(x) for x in AA_mutation]) + "\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()
+    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}
@@ -94,78 +98,80 @@
 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
+    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 + 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())
+    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")
+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
+# for testing
 seq_motif_file = outfile[:outfile.rindex("/")] + "/motif_per_seq.txt"
-first = True
 with open(seq_motif_file, 'w') as o:
-	for ID in IDlist:
-		if first:
-			o.write("ID\tRGYWC\tWRCY\tWA\tTW\n")
-			first = False
-			continue
-		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")
+    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")