|
4
|
1 # Consol_fit! It's a script & it'll consolidate your fitness values if you got them from a looping trimming pipeline instead of the standard split-by-transposon pipeline. That's it!
|
|
|
2
|
|
|
3 # Test: python ../script/consol_fit.py -calctxt results/py_2_L3_2394eVI_Gluc.txt -wig gview/consol_L3_2394eVI_Gluc.wig -i results/py_L3_2394eVI_Gluc.csv -out results/consol_L3_2394eVI_Gluc.csv -out2 results/py_2_L3_2394eVI_Gluc.csv -normalize tigr4_normal.txt
|
|
|
4
|
|
|
5 # Test: python ../script/consol_fit.py -calctxt results/py_2_L3_2394eVI_Gluc.txt -wig gview/consol_L3_2394eVI_Gluc.wig -i results/galaxy_test.csv -out results/consol_L3_2394eVI_Gluc.csv -out2 results/py_2_L3_2394eVI_Gluc.csv -normalize tigr4_normal.txt
|
|
|
6
|
|
|
7 import math
|
|
|
8 import csv
|
|
|
9
|
|
|
10
|
|
|
11
|
|
|
12
|
|
|
13
|
|
|
14
|
|
|
15
|
|
|
16
|
|
|
17
|
|
|
18
|
|
|
19 ##### ARGUMENTS #####
|
|
|
20
|
|
|
21 def print_usage():
|
|
|
22 print "\n" + "You are missing one or more required flags. A complete list of flags accepted by calc_fitness is as follows:" + "\n\n"
|
|
|
23 print "\033[1m" + "Required" + "\033[0m" + "\n"
|
|
|
24 print "-i" + "\t\t" + "The calc_fit file to be consolidated" + "\n"
|
|
|
25 print "-out" + "\t\t" + "Name of a file to enter the .csv output." + "\n"
|
|
|
26 print "-out2" + "\t\t" + "Name of a file to put the percent blank score in (used in aggregate)." + "\n"
|
|
|
27 print "-calctxt" + "\t\t" + "The txt file output from calc_fit" + "\n"
|
|
|
28 print "-normalize" + "\t" + "A file that contains a list of genes that should have a fitness of 1" + "\n"
|
|
|
29 print "\n"
|
|
|
30 print "\033[1m" + "Optional" + "\033[0m" + "\n"
|
|
|
31 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"
|
|
|
32 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"
|
|
|
33 print "-wig" + "\t\t" + "Create a wiggle file for viewing in a genome browser. Provide a filename." + "\n"
|
|
|
34 print "-maxweight" + "\t" + "The maximum weight a transposon gene can have in normalization calculations" + "\n"
|
|
|
35 print "-multiply" + "\t" + "Multiply all fitness scores by a certain value (e.g., the fitness of a knockout). You should normalize the data." + "\n"
|
|
|
36 print "\n"
|
|
|
37
|
|
|
38 import argparse
|
|
|
39 parser = argparse.ArgumentParser()
|
|
|
40 parser.add_argument("-calctxt", action="store", dest="calctxt")
|
|
|
41 parser.add_argument("-normalize", action="store", dest="normalize")
|
|
|
42 parser.add_argument("-i", action="store", dest="input")
|
|
|
43 parser.add_argument("-out", action="store", dest="outfile")
|
|
|
44 parser.add_argument("-out2", action="store", dest="outfile2")
|
|
|
45 parser.add_argument("-cutoff", action="store", dest="cutoff")
|
|
|
46 parser.add_argument("-cutoff2", action="store", dest="cutoff2")
|
|
|
47 parser.add_argument("-wig", action="store", dest="wig")
|
|
|
48 parser.add_argument("-maxweight", action="store", dest="max_weight")
|
|
|
49 parser.add_argument("-multiply", action="store", dest="multiply")
|
|
|
50 arguments = parser.parse_args()
|
|
|
51
|
|
|
52 # Checks that all the required arguments have actually been entered
|
|
|
53
|
|
|
54 if (not arguments.input or not arguments.outfile or not arguments.calctxt):
|
|
|
55 print_usage()
|
|
|
56 quit()
|
|
|
57
|
|
|
58 #
|
|
|
59
|
|
|
60 if (not arguments.max_weight):
|
|
|
61 arguments.max_weight = 75
|
|
|
62
|
|
|
63 #
|
|
|
64
|
|
|
65 if (not arguments.cutoff):
|
|
|
66 arguments.cutoff = 0
|
|
|
67
|
|
|
68 # 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.
|
|
|
69 # This only has an effect if it's larger than cutoff, since the normalization step references a list of insertions already affected by cutoff.
|
|
|
70
|
|
|
71 if (not arguments.cutoff2):
|
|
|
72 arguments.cutoff2 = 10
|
|
|
73
|
|
|
74 #Gets total & refname from calc_fit outfile2
|
|
|
75
|
|
|
76 with open(arguments.calctxt) as file:
|
|
|
77 calctxt = file.readlines()
|
|
|
78 total = float(calctxt[1].split()[1])
|
|
|
79 refname = calctxt[2].split()[1]
|
|
|
80
|
|
|
81
|
|
|
82
|
|
|
83
|
|
|
84
|
|
|
85
|
|
|
86
|
|
|
87
|
|
|
88
|
|
|
89
|
|
|
90 ##### CONSOLIDATING THE CALC_FIT FILE #####
|
|
|
91
|
|
|
92 with open(arguments.input) as file:
|
|
|
93 input = file.readlines()
|
|
|
94 results = [["position", "strand", "count_1", "count_2", "ratio", "mt_freq_t1", "mt_freq_t2", "pop_freq_t1", "pop_freq_t2", "gene", "D", "W", "nW"]]
|
|
|
95 i = 1
|
|
|
96 d = float(input[i].split(",")[10])
|
|
|
97 while i < len(input):
|
|
|
98 position = float(input[i].split(",")[0])
|
|
|
99 strands = input[i].split(",")[1]
|
|
|
100 c1 = float(input[i].split(",")[2])
|
|
|
101 c2 = float(input[i].split(",")[3])
|
|
|
102 gene = input[i].split(",")[9]
|
|
|
103 while i + 1 < len(input) and float(input[i+1].split(",")[0]) - position <= 4:
|
|
|
104 if i + 1 < len(input):
|
|
|
105 i += 1
|
|
|
106 c1 += float(input[i].split(",")[2])
|
|
|
107 c2 += float(input[i].split(",")[3])
|
|
|
108 strands = input[i].split(",")[1]
|
|
|
109 if strands[0] == 'b':
|
|
|
110 new_strands = 'b/'
|
|
|
111 elif strands[0] == '+':
|
|
|
112 if input[i].split(",")[1][0] == 'b':
|
|
|
113 new_strands = 'b/'
|
|
|
114 elif input[i].split(",")[1][0] == '+':
|
|
|
115 new_strands = '+/'
|
|
|
116 elif input[i].split(",")[1][0] == '-':
|
|
|
117 new_strands = 'b/'
|
|
|
118 elif strands[0] == '-':
|
|
|
119 if input[i].split(",")[1][0] == 'b':
|
|
|
120 new_strands = 'b/'
|
|
|
121 elif input[i].split(",")[1][0] == '+':
|
|
|
122 new_strands = 'b/'
|
|
|
123 elif input[i].split(",")[1][0] == '-':
|
|
|
124 new_strands = '-/'
|
|
|
125 if len(strands) == 3:
|
|
|
126 if len(input[i].split(",")[1]) < 3:
|
|
|
127 new_strands += strands[2]
|
|
|
128 elif strands[0] == 'b':
|
|
|
129 new_strands += 'b'
|
|
|
130 elif strands[0] == '+':
|
|
|
131 if input[i].split(",")[1][2] == 'b':
|
|
|
132 new_strands += 'b'
|
|
|
133 elif input[i].split(",")[1][2] == '+':
|
|
|
134 new_strands += '+'
|
|
|
135 elif input[i].split(",")[1][2] == '-':
|
|
|
136 new_strands += 'b'
|
|
|
137 elif strands[0] == '-':
|
|
|
138 if input[i].split(",")[1][2] == 'b':
|
|
|
139 new_strands += 'b'
|
|
|
140 elif input[i].split(",")[1][2] == '+':
|
|
|
141 new_strands += 'b'
|
|
|
142 elif input[i].split(",")[1][2] == '-':
|
|
|
143 new_strands += '-'
|
|
|
144 else:
|
|
|
145 if len(input[i].split(",")[1]) == 3:
|
|
|
146 new_strands += input[i].split(",")[1][2]
|
|
|
147 strands = new_strands
|
|
|
148 i +=1
|
|
|
149 if c2 != 0:
|
|
|
150 ratio = c2/c1
|
|
|
151 else:
|
|
|
152 ratio = 0
|
|
|
153 mt_freq_t1 = c1/total
|
|
|
154 mt_freq_t2 = c2/total
|
|
|
155 pop_freq_t1 = 1 - mt_freq_t1
|
|
|
156 pop_freq_t2 = 1 - mt_freq_t2
|
|
|
157 w = 0
|
|
|
158 if mt_freq_t2 != 0:
|
|
|
159 top_w = math.log(mt_freq_t2*(d/mt_freq_t1))
|
|
|
160 bot_w = math.log(pop_freq_t2*(d/pop_freq_t1))
|
|
|
161 w = top_w/bot_w
|
|
|
162 row = [position, strands, c1, c2, ratio, mt_freq_t1, mt_freq_t2, pop_freq_t1, pop_freq_t2, gene, d, w, w]
|
|
|
163 results.append(row)
|
|
|
164 with open(arguments.outfile, "wb") as csvfile:
|
|
|
165 writer = csv.writer(csvfile)
|
|
|
166 writer.writerows(results)
|
|
|
167
|
|
|
168
|
|
|
169
|
|
|
170
|
|
|
171
|
|
|
172
|
|
|
173
|
|
|
174
|
|
|
175
|
|
|
176
|
|
|
177 ##### REDOING NORMALIZATION #####
|
|
|
178
|
|
|
179 # 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.
|
|
|
180 # 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.
|
|
|
181
|
|
|
182 if (arguments.wig):
|
|
|
183 wigstring = "track type=wiggle_0 name=" + arguments.wig + "\n" + "variableStep chrom=" + refname + "\n"
|
|
|
184
|
|
|
185 # 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.
|
|
|
186
|
|
|
187 if (arguments.normalize):
|
|
|
188
|
|
|
189 # 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)
|
|
|
190
|
|
|
191 with open(arguments.normalize) as file:
|
|
|
192 transposon_genes = file.read().splitlines()
|
|
|
193 print "Normalize genes loaded" + "\n"
|
|
|
194 blank_ws = 0
|
|
|
195 sum = 0
|
|
|
196 count = 0
|
|
|
197 weights = []
|
|
|
198 scores = []
|
|
|
199 for list in results:
|
|
|
200
|
|
|
201 # 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!
|
|
|
202 # 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.
|
|
|
203
|
|
|
204 if list[9] != '' and list[9] in transposon_genes and list[11]:
|
|
|
205 c1 = list[2]
|
|
|
206 c2 = list[3]
|
|
|
207 score = list[11]
|
|
|
208 avg = (c1 + c2)/2
|
|
|
209
|
|
|
210 # Skips over those insertion locations with too few insertions - their fitness values are less accurate because they're based on such small insertion numbers.
|
|
|
211
|
|
|
212 if float(c1) >= float(arguments.cutoff2):
|
|
|
213
|
|
|
214 # Sets a max weight, to prevent insertion location scores with huge weights from unbalancing the normalization.
|
|
|
215
|
|
|
216 if (avg >= float(arguments.max_weight)):
|
|
|
217 avg = float(arguments.max_weight)
|
|
|
218
|
|
|
219 # 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.
|
|
|
220 # 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
|
|
|
221 # of chance not fitness; all mutants with an insertion in a specific transposon gene could be flushed out by chance!
|
|
|
222
|
|
|
223 if score == 0:
|
|
|
224 blank_ws += 1
|
|
|
225
|
|
|
226 # 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
|
|
|
227
|
|
|
228 sum += score
|
|
|
229 count += 1
|
|
|
230
|
|
|
231 # 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]
|
|
|
232
|
|
|
233 weights.append(avg)
|
|
|
234 scores.append(score)
|
|
|
235
|
|
|
236 print str(list[9]) + " " + str(score) + " " + str(c1)
|
|
|
237
|
|
|
238 # 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.
|
|
|
239
|
|
|
240 blank_count = 0
|
|
|
241 original_count = len(scores)
|
|
|
242 i = 0
|
|
|
243 while i < original_count:
|
|
|
244 w_value = scores[i]
|
|
|
245 if w_value == 0:
|
|
|
246 blank_count += 1
|
|
|
247 weights.pop[i]
|
|
|
248 scores.pop[i]
|
|
|
249 i-=1
|
|
|
250 i += 1
|
|
|
251
|
|
|
252 # 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.
|
|
|
253
|
|
|
254 if len(scores) == 0:
|
|
|
255 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"
|
|
|
256 quit()
|
|
|
257
|
|
|
258 # 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.
|
|
|
259
|
|
|
260 pc_blank_normals = float(blank_count) / float(original_count)
|
|
|
261 print "# blank out of " + str(original_count) + ": " + str(pc_blank_normals) + "\n"
|
|
|
262 with open(arguments.outfile2, "w") as f:
|
|
|
263 f.write("blanks: " + str(pc_blank_normals) + "\n" + "total: " + str(total) + "\n" + "refname: " + refname)
|
|
|
264
|
|
|
265
|
|
|
266 # 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.
|
|
|
267
|
|
|
268 average = sum / count
|
|
|
269 i = 0
|
|
|
270 weighted_sum = 0
|
|
|
271 weight_sum = 0
|
|
|
272 while i < len(weights):
|
|
|
273 weighted_sum += weights[i]*scores[i]
|
|
|
274 weight_sum += weights[i]
|
|
|
275 i += 1
|
|
|
276 weighted_average = weighted_sum/weight_sum
|
|
|
277
|
|
|
278 # Prints the regular average, weighted average, and total insertions for reference
|
|
|
279
|
|
|
280 print "Normalization step:" + "\n"
|
|
|
281 print "Regular average: " + str(average) + "\n"
|
|
|
282 print "Weighted Average: " + str(weighted_average) + "\n"
|
|
|
283 print "Total Insertions: " + str(count) + "\n"
|
|
|
284
|
|
|
285 # The actual normalization happens here; every fitness score is divided by the average fitness found for genes that should have a value of 1.
|
|
|
286 # 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.
|
|
|
287
|
|
|
288 old_ws = 0
|
|
|
289 new_ws = 0
|
|
|
290 wcount = 0
|
|
|
291 for list in results:
|
|
|
292 if list[11] == 'W':
|
|
|
293 continue
|
|
|
294 new_w = float(list[11])/weighted_average
|
|
|
295
|
|
|
296 # Sometimes you want to multiply all the fitness values by a constant; this does that.
|
|
|
297 # 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
|
|
|
298 # 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
|
|
|
299 # 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!
|
|
|
300
|
|
|
301 if arguments.multiply:
|
|
|
302 new_w *= float(arguments.multiply)
|
|
|
303
|
|
|
304 # 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).
|
|
|
305
|
|
|
306 if float(list[11]) > 0:
|
|
|
307 old_ws += float(list[11])
|
|
|
308 new_ws += new_w
|
|
|
309 wcount += 1
|
|
|
310
|
|
|
311 # Writes the new w score into the results list of lists.
|
|
|
312
|
|
|
313 list[12] = new_w
|
|
|
314
|
|
|
315 # Adds a line to wiglist for each insertion position, with the insertion position and its new w value.
|
|
|
316
|
|
|
317 if (arguments.wig):
|
|
|
318 wigstring += str(list[0]) + " " + str(new_w) + "\n"
|
|
|
319
|
|
|
320 # Prints the old w mean and new w mean for reference.
|
|
|
321
|
|
|
322 old_w_mean = old_ws / wcount
|
|
|
323 new_w_mean = new_ws / wcount
|
|
|
324 print "Old W Average: " + str(old_w_mean) + "\n"
|
|
|
325 print "New W Average: " + str(new_w_mean) + "\n"
|
|
|
326
|
|
|
327 # Overwrites the old file with the normalized file.
|
|
|
328
|
|
|
329 with open(arguments.outfile, "wb") as csvfile:
|
|
|
330 writer = csv.writer(csvfile)
|
|
|
331 writer.writerows(results)
|
|
|
332
|
|
|
333 # If a WIG file was requested, actually creates the WIG file and writes wiglist to it
|
|
|
334 # 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.
|
|
|
335
|
|
|
336 if (arguments.wig):
|
|
|
337 if (arguments.normalize):
|
|
|
338 with open(arguments.wig, "wb") as wigfile:
|
|
|
339 wigfile.write(wigstring)
|
|
|
340 else:
|
|
|
341 for list in results:
|
|
|
342 wigstring += str(list[0]) + " " + str(list[11]) + "\n"
|
|
|
343 with open(arguments.wig, "wb") as wigfile:
|
|
|
344 wigfile.write(wigstring) |
