Mercurial > repos > kaymccoy > calculate_fitnesses

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/calc_fitness.py	Thu Aug 11 18:34:33 2016 -0400
@@ -0,0 +1,553 @@
+# A translation of calc_fitness.pl into python! For analysis of Tn-Seq.
+# This script requires BioPython, which in turn has a good number of dependencies (some optional but very helpful).
+# How to install BioPython and a list of its dependencies can be found here: http://biopython.org/DIST/docs/install/Installation.html
+# To see what future edits / tests are planned for this script, search for the phrase "in the future".
+
+
+
+
+
+
+
+
+
+
+##### ARGUMENTS #####
+
+def print_usage():
+	print "\n" + "You are missing one or more required flags. A complete list of flags accepted by calc_fitness is as follows:" + "\n\n"
+	print "\033[1m" + "Required" + "\033[0m" + "\n"
+	print "-ref" + "\t\t" + "The name of the reference genome file, in GenBank format." + "\n"
+	print "-t1" + "\t\t" + "The name of the bowtie mapfile from time 1." + "\n"
+	print "-t2" + "\t\t" + "The name of the bowtie mapfile from time 2." + "\n"
+	print "-out" + "\t\t" + "Name of a file to enter the .csv output." + "\n"
+	print "\n"
+	print "\033[1m" + "Optional" + "\033[0m" + "\n"
+	print "-expansion" + "\t\t" + "Expansion factor (default: 250)" + "\n"
+	print "-d" + "\t\t" + "All reads being analyzed are downstream of the transposon" + "\n"
+	print "-reads1" + "\t\t" + "The number of reads to be used to calculate the correction factor for time 0." + "\n\t\t" + "(default counted from bowtie output)" + "\n"
+	print "-reads2" + "\t\t" + "The number of reads to be used to calculate the correction factor for time 6." + "\n\t\t" + "(default counted from bowtie output)" + "\n"
+	print "-cutoff" + "\t\t" + "Discard any positions where the average of counted transcripts at time 0 and time 1 is below this number (default 0)" + "\n"
+	print "-cutoff2" + "\t\t" + "Discard any positions within the normalization genes where the average of counted transcripts at time 0 and time 1 is below this number (default 0)" + "\n"
+	print "-strand" + "\t\t" + "Use only the specified strand (+ or -) when counting transcripts (default: both)" + "\n"
+	print "-normalize" + "\t" + "A file that contains a list of genes that should have a fitness of 1" + "\n"
+	print "-maxweight" + "\t" + "The maximum weight a transposon gene can have in normalization calculations" + "\n"
+	print "-multiply" + "\t" + "Multiply all fitness scores by a certain value (e.g., the fitness of a knockout). You should normalize the data." + "\n"
+	print "-ef" + "\t\t" + "Exclude insertions that occur in the first N amount (%) of gene--becuase may not affect gene function." + "\n"
+	print "-el" + "\t\t" + "Exclude insertions in the last N amount (%) of the gene--considering truncation may not affect gene function." + "\n"
+	print "-wig" + "\t\t" + "Create a wiggle file for viewing in a genome browser. Provide a filename." + "\n"
+	print "\n"
+
+import argparse
+parser = argparse.ArgumentParser()
+parser.add_argument("-ref", action="store", dest="ref_genome")
+parser.add_argument("-t1", action="store", dest="mapfile1")
+parser.add_argument("-t2", action="store", dest="mapfile2")
+parser.add_argument("-out", action="store", dest="outfile")
+parser.add_argument("-out2", action="store", dest="outfile2")
+parser.add_argument("-expansion", action="store", dest="expansion_factor")
+parser.add_argument("-d", action="store", dest="downstream")
+parser.add_argument("-reads1", action="store", dest="reads1")
+parser.add_argument("-reads2", action="store", dest="reads2")
+parser.add_argument("-cutoff", action="store", dest="cutoff")
+parser.add_argument("-cutoff2", action="store", dest="cutoff2")
+parser.add_argument("-strand", action="store", dest="usestrand")
+parser.add_argument("-normalize", action="store", dest="normalize")
+parser.add_argument("-maxweight", action="store", dest="max_weight")
+parser.add_argument("-multiply", action="store", dest="multiply")
+parser.add_argument("-ef", action="store", dest="exclude_first")
+parser.add_argument("-el", action="store", dest="exclude_last")
+parser.add_argument("-wig", action="store", dest="wig")
+arguments = parser.parse_args()
+
+# Checks that all the required arguments have actually been entered
+# The reference genome is required to find what gene each insertion falls within, the mapfiles are required because the relative number of insertions between time 2 and time 1 are used to calculate each insertion location's fitness, and the outfile is required so you... get an outfile.
+
+if (not arguments.ref_genome or not arguments.mapfile1 or not arguments.mapfile2 or not arguments.outfile):
+	print_usage()
+	quit()
+
+# Sets the default value of the expansion factor to 250, which is not a particularly meaningful number; the expansion factor just needs to have some value so that rough fitness calculations can be run.
+# Measuring and inputting your own expansion factor is highly recommended, as without it the fitness calculations will not be nearly as accurate. See the equations used to calculate fitness if you're confused about that.
+
+if (not arguments.expansion_factor):
+	arguments.expansion_factor = 250
+
+# Sets the default maximum weight of a transposon gene; 75 is similarly not a particularly meaningful number.
+
+if (not arguments.max_weight):
+	arguments.max_weight = 75
+
+# Sets the default value of cutoff to 0; cutoff exists to discard positions with a low number of counted transcripts, because fitnesses calculated from them may not be very accurate, by the same reasoning that studies with low sample sizes are innacurate.
+
+if (not arguments.cutoff):
+	arguments.cutoff = 0
+
+# Sets the default value of cutoff2 to 10; cutoff2 exists to discard positions within normalization genes with a low number of counted transcripts, because fitnesses calculated from them similarly may not be very accurate.
+# This only has an effect if it's larger than cutoff, since the normalization step references a list of insertions already affected by cutoff.
+
+if (not arguments.cutoff2):
+	arguments.cutoff2 = 10
+
+#Sets the default value of strand to "both", indicating you'll use reads from both Watson & Crick strands.
+
+if (not arguments.usestrand):
+	arguments.usestrand = "both"
+
+
+
+
+
+
+##### PARSING THE REFERENCE GENOME #####
+
+def get_time():
+	import datetime
+	return datetime.datetime.now().time()
+print "\n" + "Starting: " + str(get_time()) + "\n"
+
+# Uses Biopython's SeqIO to parse the reference genome file for features; these come from its feature table, which lists all known features and their relevance - if they're a gene, a recognizable repeat, etc.
+# This is the reason that all reference genome files must be in standard GenBank format!
+
+from Bio import SeqIO
+import os.path
+handle = open(arguments.ref_genome, "rU")
+for record in SeqIO.parse(handle, "genbank"):
+    refname = record.id
+    features = record.features
+handle.close()
+
+# Makes a dictionary out of each feature that's a gene - with its gene name, start location, end location, and strand as keys to their values. Then makes a list out of all those dictionaries for ease of accessing later on.
+
+feature_list = []
+for feature in features:
+	if feature.type == "gene":
+		gene = feature.qualifiers["locus_tag"]
+		strand = feature.location.strand
+		start = float(feature.location.start)
+		end = float(feature.location.end)
+
+# Exclude_first and exclude_last are used here to exclude whatever percentage of the genes you like from calculations; e.g. a value of 0.1 for exclude_last would exclude the last 10% of all genes!
+# This can be useful because insertions at the very start or end of genes often don't actually break its function.
+
+		if (arguments.exclude_first):
+			start += (end - start) * float(arguments.exclude_first)
+		if (arguments.exclude_last):
+			end -= (end - start) * float(arguments.exclude_last)
+		feature_dictionary = {"gene": gene, "start": start, "end": end, "strand": strand}
+		feature_list.append(feature_dictionary)
+
+print "Done generating feature lookup: " + str(get_time()) + "\n"
+
+
+
+
+
+
+
+
+
+
+##### PARSING THE MAPFILES #####
+
+with open(arguments.mapfile1) as file:
+	r1 = file.readlines()
+with open(arguments.mapfile2) as file:
+	r2 = file.readlines()
+
+# When called, goes through each line of the mapfile (each line is a compressed read) to find the strand, count, and position of the read. It may be helpful to look at how the mapfiles are formatted to understand how this code finds them.
+# The strand of a read is + (Watson) or - (Crick).
+# The count is the number of times a particular read was actually sequenced (counted by fastx_collapser before calc_fitness.py is even called).
+# The position is the position of the sequence when alligned to the reference genome (found by bowtie, also before calc_fitness.py is called). This is the position of the FIRST nucleotide in the sequence, when mapped.
+# Also records the total numbers of reads for the plus and minus strands, and the number of plus and minus sites from the mapfiles - not for fitness calculations, just for reference, as they're printed later.
+
+def read_mapfile(reads):
+	plus_total = 0
+	minus_total = 0
+	plus_counts = {"total": 0, "sites": 0}
+	minus_counts = {"total": 0, "sites": 0}
+	for read in reads:
+		if "-" in read.split()[0]:
+			strand = read.split()[1]
+			count = float(read.split()[0].split("-")[1])
+			position = float(read.split()[3])
+		else:
+			continue
+
+# If for some reason you want to skip all reads from one of the strands - for example, if you wanted to compare the two strands - that's done here.
+
+		if arguments.usestrand != "both" and strand != arguments.usestrand:
+			continue
+
+# Makes a dictionary for the + strand, with each insert position as a key and the number of insertions there as its corresponding value.
+# Here we run into a tricky problem, in that the actual location of an insertion would be right before the position of the first nucleotide in a flanking sequence "following" the insertion, but right after the last nucleotide in a flanking sequence "preceding" the insertion.
+# Thus the location of an insertion from a following flanking sequence would be the same as the mapped position of the sequence, but the location of an insertion from a preceding flanking sequence would be its mapped position + the length of the flanking sequence - 2.
+# The -2 comes from a fake "TA" in the read, caused by MmeI cutting unevenly at an "AT" site so that an extra "TA" is created, artifically adding +2 to the flanking sequence length.
+# However, it is currently impossible to distinguish between following and preceding flanking sequences, and so we split the difference by adding (the length of the sequence - 2) to only the + strand reads.
+# In the future we will try not adding anything to both, or adding 1/2(the sequence length -2) to both + an - strands, to see if that improves calculations.
+# The 'if arguments.downstream' option is there for a different method than the standard, in which reads are only taken from downstream of the transposon.
+
+		if (strand == "+"):
+			sequence_length = len(read.split()[4])
+			if arguments.downstream:
+				position += 0
+			else:
+				position += (sequence_length - 2)
+			plus_counts["total"] += count
+			plus_counts["sites"] += 1
+			if position in plus_counts:
+				plus_counts[position] += count
+			else:
+				plus_counts[position] = count
+
+# Makes a dictionary for the - strand, with each insert position as a key and the number of insertions there as its corresponding value.
+
+		else:
+			minus_counts["total"] += count
+			minus_counts["sites"] += 1
+			if position in minus_counts:
+				minus_counts[position] += count
+			else:
+				minus_counts[position] = count
+	return (plus_counts, minus_counts)
+
+# Calls read_mapfile(reads) to parse arguments.reads1 and arguments.reads2 (your reads from t1 and t2).
+
+(plus_ref_1, minus_ref_1) = read_mapfile(r1)
+print "Read first file: " + str(get_time()) + "\n"
+(plus_ref_2, minus_ref_2) = read_mapfile(r2)
+print "Read second file: " + str(get_time()) + "\n"
+
+# Prints the number of + and - reads from reads1 and reads2, as well as the number of + and - sites found from the mapfiles, for reference when debugging etc.
+# The number of sites is the length of a given dictionary of sites - 1 because its last key, "total", isn't actually a site.
+
+print "Reads:" + "\n"
+print "1: + " + str(plus_ref_1["total"]) + " - " + str(minus_ref_1["total"]) + "\n"
+print "2: + " + str(plus_ref_2["total"]) + " - " + str(minus_ref_2["total"]) + "\n"
+print "Sites:" + "\n"
+print "1: + " + str(plus_ref_1["sites"]) + " - " + str(minus_ref_1["sites"]) + "\n"
+print "2: + " + str(plus_ref_2["sites"]) + " - " + str(minus_ref_2["sites"]) + "\n"
+
+
+
+
+
+
+
+
+
+
+##### FITNESS CALCULATIONS #####
+
+# If reads1 and reads2 weren't specified in the command line, sets them as the total number of reads (found in read_mapfile())
+
+if not arguments.reads1:
+	arguments.reads1 = plus_ref_1["total"] + minus_ref_1["total"]
+if not arguments.reads2:
+	arguments.reads2 = plus_ref_2["total"] + minus_ref_2["total"]
+
+# Calculates the correction factors for reads from t1 and t2; cfactor1 and cfactor2 are the number of reads from t1 and t2 respectively divided by total, which is the average number of reads between the two.
+# For example, if there were more reads from t2 than t1 - in a 6:4 ratio - cfactor2 would be 1.2 and cfactor1 would be 0.8
+# This is used later on to correct for pipetting errors, or any other error that would cause unequal amounts of DNA from t1 and t2 to be sequenced so that an unequal amount of reads is produced
+# However, because slightly more or less DNA shouldn't change the % frequency of a mutation (it would only change the absolute number of counts) this may not be necessary; in the future we'll try removing this.
+
+total = (float(arguments.reads1) + float(arguments.reads2))/2
+cfactor1 = float(arguments.reads1)/total
+cfactor2 = float(arguments.reads2)/total
+print "Cfactor 1: " + str(cfactor1) + "\n"
+print "Cfactor 2: " + str(cfactor2) + "\n"
+import math
+import csv
+results = [["position", "strand", "count_1", "count_2", "ratio", "mt_freq_t1", "mt_freq_t2", "pop_freq_t1", "pop_freq_t2", "gene", "D", "W", "nW"]]
+genic = 0
+total_inserts = 0
+with open(arguments.ref_genome, "r") as file:
+	firstline = file.readline()
+genomelength = firstline.split()[2]
+i = 0
+while i < float(genomelength):
+
+# At each possible location for an insertion in the genome, counts the number of actual insertions at t1 and which strand(s) the corresponding reads came from.
+
+	c1 = 0
+	if i in plus_ref_1:
+		c1 = float(plus_ref_1[i])
+		strand = "+/"
+		if i in minus_ref_1:
+			c1 += float(minus_ref_1[i])
+			strand = "b/"
+	elif i in minus_ref_1:
+		c1 = float(minus_ref_1[i])
+		strand = "-/"
+
+# If there were no insertions at a certain location at t1 just continues to the next location; there can't be any comparison to make between t1 and t2 if there are no t1 insertions!
+
+	else:
+		i += 1
+		continue
+
+# At each location where there was an insertion at t1, counts the number of insertions at t2 and which strand(s) the corresponding reads came from.
+
+	c2 = 0
+	if i in plus_ref_2:
+		c2 = float(plus_ref_2[i])
+		if i in minus_ref_2:
+			c2 += float(minus_ref_2[i])
+			strand += "b"
+		else:
+			strand += "+"
+	elif i in minus_ref_2:
+		c2 = float(minus_ref_2[i])
+		strand += "-"
+
+# Corrects with cfactor1 and cfactor2; as noted above, this may not be necessary and we'll see what removing it does in the future. Possible alternate / non-corrected code is shown in the comments just below.
+
+ 	c1 /= cfactor1
+ 	if c2 != 0:
+ 		c2 /= cfactor2
+ 		ratio = c2/c1
+ 	else:
+ 		c2 = 0
+ 		ratio = 0
+
+#	if c2 != 0:
+#		ratio = c2/c1
+#	else:
+#		ratio = 0
+
+# Passes by all insertions with a number of reads smaller than the cutoff, as they may lead to inaccurate fitness calculations.
+
+	if (c1 + c2)/2 < float(arguments.cutoff):
+		i+= 1
+		continue
+
+# Calculates each insertion's frequency within the populations at t1 and t2. If cfactor correction were removed, "total" would have to be changed to "float(arguments.reads1)" and "float(arguments.reads2)" respectively, as shown in the comments just below.
+
+ 	mt_freq_t1 = c1/total
+ 	mt_freq_t2 = c2/total
+#	mt_freq_t1 = c1/float(arguments.reads1)
+#	mt_freq_t2 = c2/float(arguments.reads2)
+	pop_freq_t1 = 1 - mt_freq_t1
+	pop_freq_t2 = 1 - mt_freq_t2
+
+# Calculates each insertion's fitness! This is from the fitness equation log((frequency of mutation @ time 2 / frequency of mutation @ time 1)*expansion factor)/log((frequency of population without the mutation @ time 2 / frequency of population without the mutation @ time 1)*expansion factor)
+
+	w = 0
+	if mt_freq_t2 != 0:
+		top_w = math.log(mt_freq_t2*(float(arguments.expansion_factor)/mt_freq_t1))
+		bot_w = math.log(pop_freq_t2*(float(arguments.expansion_factor)/pop_freq_t1))
+		w = top_w/bot_w
+
+# Checks which gene locus the insertion falls within, and records that.
+
+	gene = ''
+	for feature_dictionary in feature_list:
+		if feature_dictionary["start"] <= i and feature_dictionary["end"] >= i:
+			gene = "".join(feature_dictionary["gene"])
+			genic += 1
+			break
+	total_inserts += 1
+
+# Writes all relevant information on each insertion and its fitness to a cvs file: the location of the insertion, its strand, c1, c2, etc. (the variable names are self-explanatiory)
+# w is written twice, because the second w will be normalized if normalization is called for, thus becoming nW.
+
+	row = [i, strand, c1, c2, ratio, mt_freq_t1, mt_freq_t2, pop_freq_t1, pop_freq_t2, gene, arguments.expansion_factor, w, w]
+	results.append(row)
+	i += 1
+with open(arguments.outfile, "wb") as csvfile:
+    writer = csv.writer(csvfile)
+    writer.writerows(results)
+
+# Prints the time when done comparing mapfiles - that is, when done with all fitness calculations but normalization - as well as the number of genic inserts and number of total inserts.
+
+print "Done comparing mapfiles " + str(get_time()) + "\n"
+print "Genic: " + str(genic) + "\n"
+print "Total: " + str(total_inserts) + "\n"
+
+
+
+
+
+
+
+
+
+
+##### NORMALIZATION #####
+
+# If making a WIG file is requested in the arguments, starts a string to be added to and then written to the WIG file with a typical WIG file header.
+# The header is just in a typical WIG file format; if you'd like to look into this more UCSC has notes on formatting WIG files on their site.
+
+if (arguments.wig):
+	wigstring = "track type=wiggle_0 name=" + arguments.wig + "\n" + "variableStep chrom=" + refname + "\n"
+
+# If a file's given for normalization, starts normalization; this corrects for anything that would cause all the fitness values to be too high or too low.
+
+if (arguments.normalize):
+
+# Makes a list of the genes in the normalization file, which should all be transposon genes (these naturally ought to have a fitness value of exactly 1, because transposons are generally non-coding DNA)
+
+	with open(arguments.normalize) as file:
+		transposon_genes = file.read().splitlines()
+	print "Normalize genes loaded" + "\n"
+	blank_ws = 0
+	sum = 0
+	count = 0
+	weights = []
+	scores = []
+	for list in results:
+
+# Finds all insertions within one of the normalization genes that also have a w value; gets their c1 and c2 values (the number of insertions at t1 and t2) and takes the average of that!
+# The average is later used as the "weight" of an insertion location's fitness - if it's had more insertions, it should weigh proportionally more towards the average fitness of insertions within the normalization genes.
+
+		if list[9] != '' and list[9] in transposon_genes and list[11]:
+			c1 = list[2]
+			c2 = list[3]
+			score = list[11]
+			avg = (c1 + c2)/2
+
+# Skips over those insertion locations with too few insertions - their fitness values are less accurate because they're based on such small insertion numbers.
+
+			if float(c1) >= float(arguments.cutoff2):
+			#if avg >= float(arguments.cutoff2):
+
+# Sets a max weight, to prevent insertion location scores with huge weights from unbalancing the normalization.
+
+				if (avg >= float(arguments.max_weight)):
+					avg = float(arguments.max_weight)
+
+# Tallies how many w values are 0 within the blank_ws value; you might get many transposon genes with a w value of 0 if a bottleneck occurs, which is especially common with in vivo experiments.
+# For example, when studying a nasal infection in a mouse model, what bacteria "sticks" and is able to survive and what bacteria is swallowed and killed or otherwise flushed out tends to be a matter
+# of chance not fitness; all mutants with an insertion in a specific transposon gene could be flushed out by chance!
+
+				if score == 0:
+					blank_ws += 1
+
+# Adds the fitness values of the insertions within normalization genes together and increments count so their average fitness (sum/count) can be calculated later on
+
+				sum += score
+				count += 1
+
+# Records the weights of the fitness values of the insertion locations in corresponding lists - for example, weights[2] would be the weight of the fitness value at score[2]
+
+				weights.append(avg)
+				scores.append(score)
+
+				print str(list[9]) + " " + str(score) + " " + str(c1)
+
+# Counts and removes all "blank" fitness values of normalization genes - those that = 0 - because they most likely don't really have a fitness value of 0, and you just happened to not get any reads from that location at t2.
+
+	blank_count = 0
+	original_count = len(scores)
+	i = 0
+	while i < original_count:
+		w_value = scores[i]
+		if w_value == 0:
+			blank_count += 1
+			weights.pop[i]
+			scores.pop[i]
+			i-=1
+		i += 1
+
+# If no normalization genes can pass the cutoff, normalization cannot occur, so this ends the script advises the user to try again and lower cutoff and/or cutoff2.
+
+	if len(scores) == 0:
+		print 'ERROR: The normalization genes do not have enough reads to pass cutoff and/or cutoff2; please lower one or both of those arguments.' + "\n"
+		quit()
+
+# Prints the number of of blank fitness values found and removed for reference. Writes the percentage to a file so it can be referenced for aggregate analysis.
+
+	pc_blank_normals = float(blank_count) / float(original_count)
+	print "# blank out of " + str(original_count) + ": " + str(pc_blank_normals) + "\n"
+	with open(arguments.outfile2, "w") as f:
+		f.write("blanks: " + str(pc_blank_normals) + "\n" + "total: " + str(total) + "\n" + "refname: " + refname)
+
+
+# Finds "average" - the average fitness value for an insertion within the transposon genes - and "weighted_average" - the average fitness value for an insertion within the transposon genes weighted by how many insertions each had.
+
+	average = sum / count
+	i = 0
+	weighted_sum = 0
+	weight_sum = 0
+	while i < len(weights):
+		weighted_sum += weights[i]*scores[i]
+		weight_sum += weights[i]
+		i += 1
+	weighted_average = weighted_sum/weight_sum
+
+# Prints the regular average, weighted average, and total insertions for reference
+
+	print "Normalization step:" + "\n"
+	print "Regular average: " + str(average) + "\n"
+	print "Weighted Average: " + str(weighted_average) + "\n"
+	print "Total Insertions: " + str(count) + "\n"
+
+# The actual normalization happens here; every fitness score is divided by the average fitness found for genes that should have a value of 1.
+# For example, if the average fitness for genes was too low overall - let's say 0.97 within the normalization geness - every fitness would be proportionally raised.
+
+	old_ws = 0
+	new_ws = 0
+	wcount = 0
+	for list in results:
+		if list[11] == 'W':
+			continue
+		new_w = float(list[11])/weighted_average
+
+# Sometimes you want to multiply all the fitness values by a constant; this does that.
+# For example you might multiply all the values by a constant for a genetic interaction screen - where Tn-Seq is performed as usual except there's one background knockout all the mutants share. This is
+# because independent mutations should have a fitness value that's equal to their individual fitness values multipled, but related mutations will deviate from that; to find those deviations you'd multiply
+# all the fitness values from mutants from a normal library by the fitness of the background knockout and compare that to the fitness values found from the knockout library!
+
+		if arguments.multiply:
+			new_w *= float(arguments.multiply)
+
+# Records the old w score for reference, and adds it to a total sum of all w scores (so that the old w mean and new w mean can be printed later).
+
+		if float(list[11]) > 0:
+			old_ws += float(list[11])
+			new_ws += new_w
+			wcount += 1
+
+# Writes the new w score into the results list of lists.
+
+		list[12] = new_w
+
+# Adds a line to wiglist for each insertion position, with the insertion position and its new w value.
+
+		if (arguments.wig):
+			wigstring += str(list[0]) + " " + str(new_w) + "\n"
+
+# Prints the old w mean and new w mean for reference.
+
+	old_w_mean = old_ws / wcount
+	new_w_mean = new_ws / wcount
+	print "Old W Average: " + str(old_w_mean) + "\n"
+	print "New W Average: " + str(new_w_mean) + "\n"
+
+# Overwrites the old file with the normalized file.
+
+with open(arguments.outfile, "wb") as csvfile:
+    writer = csv.writer(csvfile)
+    writer.writerows(results)
+
+# If a WIG file was requested, actually creates the WIG file and writes wiglist to it
+# So what's written here is the WIG header plus each insertion position and it's new w value if normalization were called for, and each insertion position and its unnormalized w value if normalization were not called for.
+
+if (arguments.wig):
+	if (arguments.normalize):
+		with open(arguments.wig, "wb") as wigfile:
+			wigfile.write(wigstring)
+	else:
+		for list in results:
+			wigstring += str(list[0]) + " " + str(list[11]) + "\n"
+		with open(arguments.wig, "wb") as wigfile:
+				wigfile.write(wigstring)
+
+
+
+
+
+#FITOPSpy="-ef .0 -el .10 -cutoff 0"
+#
+#python ../script/calc_fitness.py $FITOPSpy -wig gview/py_L1_2394eVI_Gluc.wig -t1 alignments/L1_2394eVI_Input.map -t2 alignments/L1_2394eVI_Gluc_T2.map -ref=NC_003028b2.gbk -out results/py_L1_2394eVI_Gluc.csv -out2 results/py_2_L1_2394eVI_Gluc.txt -expansion 675 -normalize tigr4_normal.txt
+#python ../script/calc_fitness.py $FITOPSpy -wig gview/py_L3_2394eVI_Gluc.wig -t1 alignments/L3_2394eVI_Input.map -t2 alignments/L3_2394eVI_Gluc_T2.map -ref=NC_003028b2.gbk -out results/py_L3_2394eVI_Gluc.csv -out2 results/py_2_L3_2394eVI_Gluc.txt -expansion 244 -normalize tigr4_normal.txt
\ No newline at end of file