--- a/mutation_analysis.py	Fri Jan 29 08:11:31 2016 -0500
+++ b/mutation_analysis.py	Mon Feb 29 10:49:39 2016 -0500
@@ -1,3 +1,4 @@
+from __future__ import division
 from collections import defaultdict
 import re
 import argparse
@@ -56,8 +57,7 @@
 		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 + "_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]
@@ -209,35 +209,68 @@
 				WACount[ID] += (1.0 * int(mutation_in_WA)) / in_how_many_motifs
 				TWCount[ID] += (1.0 * int(mutation_in_TW)) / in_how_many_motifs
 
+
+def mean(lst):
+	return (float(sum(lst)) / len(lst)) if len(lst) > 0 else 0.0
+
+
+def median(lst):
+	lst = sorted(lst)
+	l = len(lst)
+	if l == 0:
+		return 0
+	if l == 1:
+		return lst[0]
+		
+	l = int(l / 2)
+	
+	if len(lst) % 2 == 0:
+		print "list length is", l
+		return float(lst[l] + lst[(l - 1)]) / 2.0
+	else:
+		return lst[l]
+
+funcs = {"mean": mean, "median": median, "sum": sum}
+
 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())
+
+for fname in funcs.keys():
+	for gene in genes:
+		with open(directory + gene + "_" + fname + "_value.txt", 'r') as v:
+			valuedic[gene + "_" + fname] = float(v.readlines()[0].rstrip())
+	with open(directory + "all_" + fname + "_value.txt", 'r') as v:
+		valuedic["total_" + fname] = float(v.readlines()[0].rstrip())
+	
+print valuedic
+	
+def get_xyz(lst, gene, f, fname):
+	x = int(round(f(lst)))
+	y = valuedic[gene + "_" + fname]
+	z = str(round(x / float(valuedic[gene + "_" + fname]) * 100, 1)) if valuedic[gene + "_" + fname] != 0 else "0"
+	return (str(x), str(y), z)
 
 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 fname in funcs.keys():
+	func = funcs[fname]
+	foutfile = outfile[:outfile.rindex("/")] + "/hotspot_analysis_" + fname + ".txt"
+	with open(foutfile, 'w') as o:
+		for typ in arr:
+			o.write(typ + " (%)")
+			curr = dic[typ]
+			for gene in genes:
+				geneMatcher = re.compile(".*" + gene + ".*")
+				if valuedic[gene + "_" + fname] is 0:
+					o.write(",0,0,0")
+				else:
+					x, y, z = get_xyz([curr[x] for x in [y for y, z in genedic.iteritems() if geneMatcher.match(z)]], gene, func, fname)
+					o.write("," + x + "," + y + "," + z)
+			# for total
+			x, y, z = get_xyz([y for x, y in curr.iteritems()], "total", func, fname)
+			o.write("," + x + "," + y + "," + z + "\n")
 
 
 # for testing
@@ -245,5 +278,4 @@
 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")
+		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")