annotate consol_fit.py @ 8:b108f7fe9493 draft default tip

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