| blob | ce39fd63c0f0f28c37cc26286f8727a1d188f40f |
1 #! /usr/bin/env python
2 import os
3 import glob
4 import copy
5 import math
6 import random
7 import subprocess
8 from subprocess import call
10 def find_accuracy(fn_truth, fn_assemby, name, method = "panda", printnice = True):
11 numseqs = 0
12 numassebled = 0
13 numcorrect = 0
14 ftruth = open(fn_truth)
15 seq_len_map = {}
16 line = ftruth.readline()
17 while line!="":
18 ss = line.strip().split()
19 seq_len_map[ss[0]] = int(ss[1])
20 numseqs = numseqs + 1
21 line = ftruth.readline()
22 ftruth.close()
24 fin = open(fn_assemby)
25 line = fin.readline()
26 while line!="":
27 seqname = line.strip()
28 if method == "shera":
29 seqname = seqname[1:].split("_")[0]
30 seq = fin.readline().strip()
31 as_seq_len = len(seq)
32 true_seq_len = int(seq_len_map[seqname])
33 numassebled = numassebled + 1
34 if true_seq_len == as_seq_len:
35 numcorrect = numcorrect + 1
36 else:
37 if method == "panda":
38 seqname = seqname[1:].split(":")[0]
39 print(seqname)
40 elif method == "flash":
41 seqname = seqname[1:]
42 elif method == "cope":
43 #@cp0 EU861894-140/1_EU861894-140/2 28 172
44 seqname = seqname[1:].split()[1].split("/")[0]
45 else:
46 seqname = seqname[1:].split("/")[0]
47 seq = fin.readline().strip()
48 fin.readline()
49 qual = fin.readline().strip()
50 as_seq_len = len(seq)
51 true_seq_len = int(seq_len_map[seqname])
52 numassebled = numassebled + 1
53 if true_seq_len == as_seq_len:
54 numcorrect = numcorrect + 1
55 line = fin.readline()
56 fin.close()
58 if numseqs == 0:
59 numseqs = 1
60 if numassebled == 0:
61 numassebled = 1
63 if printnice:
64 print("Assebly method:" + method)
65 print("Total seqs:" + repr(numseqs))
66 print("Assembled seqs:" + repr(numassebled))
67 print("Corrected seqs:" + repr(numcorrect))
68 print("Cor / Total:" + repr(float(numcorrect)/float(numseqs)))
69 print("Cor / Assemled:" + repr(float(numcorrect)/float(numassebled)))
70 else:
71 if numcorrect == 0:
72 numcorrect = numseqs - numassebled
73 print(name + " " +repr(numseqs)+" "+repr(numassebled)+" "+repr(numcorrect)+" "+repr(float(numcorrect)/float(numseqs))+" " + repr(float(numcorrect)/float(numseqs)) + " " + repr(1- float(numcorrect)/float(numseqs)))
74 else:
75 print(name + " " +repr(numseqs)+" "+repr(numassebled)+" "+repr(numcorrect)+" "+repr(float(numcorrect)/float(numseqs))+" " + repr(float(numcorrect)/float(numassebled)) + " " + repr(1- float(numcorrect)/float(numassebled)))
77 return numcorrect
79 def test_pear(forward, reverse, output, pvalue, minoverlap, maxlen, minlen, mintrimlen, minquality, maxuncalled, scoremethod, empirical_freqs, truelenfile, testname, testw):
80 if empirical_freqs == "yes":
81 #print "../src/pear" + " " + "-f" + " " + forward + " " + "-r" + " " + reverse + " " + "-o" + " " + output + " " + "-p" + " " + pvalue + " " + "-v" + " " + minoverlap + " " + "-m" + " " + maxlen + " " + "-n" + " " + minlen + " " + "-t" + " " + mintrimlen + " " + "-q" + " " + minquality + " " + "-u" + " " + maxuncalled + " " + "-s" + " " + scoremethod + " " + "-j" + " " + "2" + " " + "-g" + " " + testw + " " + "-y" + " " + "200000000"
82 call(["../src/pear","-f",forward,"-r",reverse,"-o", output,"-p", pvalue, "-v", minoverlap, "-m", maxlen, "-n", minlen, "-t", mintrimlen, "-q", minquality, "-u", maxuncalled, "-s", scoremethod, "-j", "2", "-g", testw, "-y", "200000000" ], stdout=open(os.devnull, "w"), stderr=subprocess.STDOUT)
83 else:
84 #print "../src/pear" + "-f" + forward + "-r" + reverse + "-o" + output + "-p" + pvalue + "-v" + minoverlap + "-m" + maxlen + "-n" + minlen + "-t" + mintrimlen + "-q", minquality + "-u", maxuncalled + "-s" + scoremethod + "-j" + "2" + "-g" + testw + "-y" + "200000000 -e"
85 call(["../src/pear","-f",forward,"-r",reverse,"-o", output,"-p", pvalue, "-v", minoverlap, "-m", maxlen, "-n", minlen, "-t", mintrimlen, "-q", minquality, "-u", maxuncalled, "-s", scoremethod, "-j", "2", "-g", testw, "-y", "200000000", "-e"], stdout=open(os.devnull, "w"), stderr=subprocess.STDOUT)
86 return find_accuracy(fn_truth = truelenfile, fn_assemby = output+".assembled.fastq", name = testname, method = "pear", printnice = False)
88 def publication_data():
89 print("Testname" + " " +"Reads"+" "+"Assembled"+" "+"Correct"+" "+"Correct/Reads"+" " + "Correct/Assembled" + " fp_rate")
90 num101bpt_nt = test_pear(forward = "16S_101bp1.fq.through.fastq",
91 reverse = "16S_101bp2.fq.through.fastq",
92 output = "to1",
93 pvalue = "1.0",
94 minoverlap = "1",
95 maxlen = "500",
96 minlen = "50",
97 mintrimlen = "50",
98 minquality = "0",
99 maxuncalled = "1",
100 scoremethod = "2",
101 empirical_freqs = "yes",
102 truelenfile = "16S_101bp_len.txt",
103 testname = "101bpt_nt",
104 testw = "1")
105 num101bpt_t01 = test_pear(forward = "16S_101bp1.fq.through.fastq",
106 reverse = "16S_101bp2.fq.through.fastq",
107 output = "to1",
108 pvalue = "0.01",
109 minoverlap = "1",
110 maxlen = "500",
111 minlen = "50",
112 mintrimlen = "50",
113 minquality = "0",
114 maxuncalled = "1",
115 scoremethod = "2",
116 empirical_freqs = "yes",
117 truelenfile = "16S_101bp_len.txt",
118 testname = "101bpt_t01",
119 testw = "1")
120 num101bpt_t01 = test_pear(forward = "16S_101bp1.fq.through.fastq",
121 reverse = "16S_101bp2.fq.through.fastq",
122 output = "to1",
123 pvalue = "0.01",
124 minoverlap = "1",
125 maxlen = "500",
126 minlen = "50",
127 mintrimlen = "50",
128 minquality = "0",
129 maxuncalled = "1",
130 scoremethod = "2",
131 empirical_freqs = "yes",
132 truelenfile = "16S_101bp_len.txt",
133 testname = "101bpt_t2_01",
134 testw = "2")
135 num101bp_nt = test_pear(forward = "16S_101bp1.fq",
136 reverse = "16S_101bp2.fq",
137 output = "to1",
138 pvalue = "1.0",
139 minoverlap = "1",
140 maxlen = "500",
141 minlen = "50",
142 mintrimlen = "50",
143 minquality = "0",
144 maxuncalled = "1",
145 scoremethod = "2",
146 empirical_freqs = "yes",
147 truelenfile = "16S_101bp_len.txt",
148 testname = "101bp_nt",
149 testw = "1")
150 num101bp_t01 = test_pear(forward = "16S_101bp1.fq",
151 reverse = "16S_101bp2.fq",
152 output = "to1",
153 pvalue = "0.01",
154 minoverlap = "1",
155 maxlen = "500",
156 minlen = "50",
157 mintrimlen = "50",
158 minquality = "0",
159 maxuncalled = "1",
160 scoremethod = "2",
161 empirical_freqs = "yes",
162 truelenfile = "16S_101bp_len.txt",
163 testname = "101bp_t01",
164 testw = "1")
165 num101bp_t01 = test_pear(forward = "16S_101bp1.fq",
166 reverse = "16S_101bp2.fq",
167 output = "to1",
168 pvalue = "0.01",
169 minoverlap = "1",
170 maxlen = "500",
171 minlen = "50",
172 mintrimlen = "50",
173 minquality = "0",
174 maxuncalled = "1",
175 scoremethod = "2",
176 empirical_freqs = "yes",
177 truelenfile = "16S_101bp_len.txt",
178 testname = "101bp_t2_01",
179 testw = "2")
180 num101bp_nt = test_pear(forward = "16S_150bp1.fq",
181 reverse = "16S_150bp2.fq",
182 output = "to1",
183 pvalue = "1.0",
184 minoverlap = "1",
185 maxlen = "500",
186 minlen = "50",
187 mintrimlen = "50",
188 minquality = "0",
189 maxuncalled = "1",
190 scoremethod = "2",
191 empirical_freqs = "yes",
192 truelenfile = "16S_150bp_len.txt",
193 testname = "150bp_nt",
194 testw = "1")
195 num101bp_t01 = test_pear(forward = "16S_150bp1.fq",
196 reverse = "16S_150bp2.fq",
197 output = "to1",
198 pvalue = "0.01",
199 minoverlap = "1",
200 maxlen = "500",
201 minlen = "50",
202 mintrimlen = "50",
203 minquality = "0",
204 maxuncalled = "1",
205 scoremethod = "2",
206 empirical_freqs = "yes",
207 truelenfile = "16S_150bp_len.txt",
208 testname = "150bp_t01",
209 testw = "1")
210 num101bp_t01 = test_pear(forward = "16S_150bp1.fq",
211 reverse = "16S_150bp2.fq",
212 output = "to1",
213 pvalue = "0.01",
214 minoverlap = "1",
215 maxlen = "500",
216 minlen = "50",
217 mintrimlen = "50",
218 minquality = "0",
219 maxuncalled = "1",
220 scoremethod = "2",
221 empirical_freqs = "yes",
222 truelenfile = "16S_150bp_len.txt",
223 testname = "150bp_t2_01",
224 testw = "2")
225 num101bp_nt = test_pear(forward = "16S_165bp1.fq",
226 reverse = "16S_165bp2.fq",
227 output = "to1",
228 pvalue = "1.0",
229 minoverlap = "1",
230 maxlen = "500",
231 minlen = "50",
232 mintrimlen = "50",
233 minquality = "0",
234 maxuncalled = "1",
235 scoremethod = "2",
236 empirical_freqs = "yes",
237 truelenfile = "16S_165bp_len.txt",
238 testname = "165bp_nt",
239 testw = "1")
240 num101bp_t01 = test_pear(forward = "16S_165bp1.fq",
241 reverse = "16S_165bp2.fq",
242 output = "to1",
243 pvalue = "0.01",
244 minoverlap = "1",
245 maxlen = "500",
246 minlen = "50",
247 mintrimlen = "50",
248 minquality = "0",
249 maxuncalled = "1",
250 scoremethod = "2",
251 empirical_freqs = "yes",
252 truelenfile = "16S_165bp_len.txt",
253 testname = "165bp_t01",
254 testw = "1")
255 num101bp_t01 = test_pear(forward = "16S_165bp1.fq",
256 reverse = "16S_165bp2.fq",
257 output = "to1",
258 pvalue = "0.01",
259 minoverlap = "1",
260 maxlen = "500",
261 minlen = "50",
262 mintrimlen = "50",
263 minquality = "0",
264 maxuncalled = "1",
265 scoremethod = "2",
266 empirical_freqs = "yes",
267 truelenfile = "16S_165bp_len.txt",
268 testname = "165bp_t2_01",
269 testw = "2")
270 num101bp_nt = test_pear(forward = "16S_180bp1.fq",
271 reverse = "16S_180bp2.fq",
272 output = "to1",
273 pvalue = "1.0",
274 minoverlap = "1",
275 maxlen = "500",
276 minlen = "50",
277 mintrimlen = "50",
278 minquality = "0",
279 maxuncalled = "1",
280 scoremethod = "2",
281 empirical_freqs = "yes",
282 truelenfile = "16S_180bp_len.txt",
283 testname = "180bp_nt",
284 testw = "1")
285 num101bp_t01 = test_pear(forward = "16S_180bp1.fq",
286 reverse = "16S_180bp2.fq",
287 output = "to1",
288 pvalue = "0.01",
289 minoverlap = "1",
290 maxlen = "500",
291 minlen = "50",
292 mintrimlen = "50",
293 minquality = "0",
294 maxuncalled = "1",
295 scoremethod = "2",
296 empirical_freqs = "yes",
297 truelenfile = "16S_180bp_len.txt",
298 testname = "180bp_t01",
299 testw = "1")
300 num101bp_t01 = test_pear(forward = "16S_180bp1.fq",
301 reverse = "16S_180bp2.fq",
302 output = "to1",
303 pvalue = "0.01",
304 minoverlap = "1",
305 maxlen = "500",
306 minlen = "50",
307 mintrimlen = "50",
308 minquality = "0",
309 maxuncalled = "1",
310 scoremethod = "2",
311 empirical_freqs = "yes",
312 truelenfile = "16S_180bp_len.txt",
313 testname = "180bp_t2_01",
314 testw = "2")
315 num101bp_nt = test_pear(forward = "16S_190bp1.fq",
316 reverse = "16S_190bp2.fq",
317 output = "to1",
318 pvalue = "1.0",
319 minoverlap = "1",
320 maxlen = "500",
321 minlen = "50",
322 mintrimlen = "50",
323 minquality = "0",
324 maxuncalled = "1",
325 scoremethod = "2",
326 empirical_freqs = "yes",
327 truelenfile = "16S_190bp_len.txt",
328 testname = "190bp_nt",
329 testw = "1")
330 num101bp_t01 = test_pear(forward = "16S_190bp1.fq",
331 reverse = "16S_190bp2.fq",
332 output = "to1",
333 pvalue = "0.01",
334 minoverlap = "1",
335 maxlen = "500",
336 minlen = "50",
337 mintrimlen = "50",
338 minquality = "0",
339 maxuncalled = "1",
340 scoremethod = "2",
341 empirical_freqs = "yes",
342 truelenfile = "16S_190bp_len.txt",
343 testname = "190bp_t01",
344 testw = "1")
345 num101bp_t01 = test_pear(forward = "16S_190bp1.fq",
346 reverse = "16S_190bp2.fq",
347 output = "to1",
348 pvalue = "0.01",
349 minoverlap = "1",
350 maxlen = "500",
351 minlen = "50",
352 mintrimlen = "50",
353 minquality = "0",
354 maxuncalled = "1",
355 scoremethod = "2",
356 empirical_freqs = "yes",
357 truelenfile = "16S_190bp_len.txt",
358 testname = "190bp_t2_01",
359 testw = "2")
360 num101bp_nt = test_pear(forward = "16S_250bp1.fq",
361 reverse = "16S_250bp2.fq",
362 output = "to1",
363 pvalue = "1.0",
364 minoverlap = "1",
365 maxlen = "500",
366 minlen = "50",
367 mintrimlen = "50",
368 minquality = "0",
369 maxuncalled = "1",
370 scoremethod = "2",
371 empirical_freqs = "yes",
372 truelenfile = "16S_250bp_len.txt",
373 testname = "250bp_nt",
374 testw = "1")
375 num101bp_t01 = test_pear(forward = "16S_250bp1.fq",
376 reverse = "16S_250bp2.fq",
377 output = "to1",
378 pvalue = "0.01",
379 minoverlap = "1",
380 maxlen = "500",
381 minlen = "50",
382 mintrimlen = "50",
383 minquality = "0",
384 maxuncalled = "1",
385 scoremethod = "2",
386 empirical_freqs = "yes",
387 truelenfile = "16S_250bp_len.txt",
388 testname = "250bp_t01",
389 testw = "1")
390 num101bp_t01 = test_pear(forward = "16S_250bp1.fq",
391 reverse = "16S_250bp2.fq",
392 output = "to1",
393 pvalue = "0.01",
394 minoverlap = "1",
395 maxlen = "500",
396 minlen = "50",
397 mintrimlen = "50",
398 minquality = "0",
399 maxuncalled = "1",
400 scoremethod = "2",
401 empirical_freqs = "yes",
402 truelenfile = "16S_250bp_len.txt",
403 testname = "250bp_t2_01",
404 testw = "2")
405 os.remove("to1.assembled.fastq")
406 os.remove("to1.discarded.fastq")
407 os.remove("to1.unassembled.forward.fastq")
408 os.remove("to1.unassembled.reverse.fastq")
410 def test_cope():
411 find_accuracy(fn_truth = "16S_101bp_len.txt", fn_assemby = "", name= "101bpt", method = "cope", printnice = True)
413 def main():
414 #101bp through, should have corrected assembled sequences =
415 print("Testname" + " " +"Reads"+" "+"Assembled"+" "+"Correct"+" "+"Correct/Reads"+" " + "Correct/Assembled" + " fp_rate")
416 num101bpt_t1 = test_pear(forward = "16S_101bp1.fq.through.fastq",
417 reverse = "16S_101bp2.fq.through.fastq",
418 output = "to1",
419 pvalue = "0.01",
420 minoverlap = "1",
421 maxlen = "500",
422 minlen = "50",
423 mintrimlen = "50",
424 minquality = "10",
425 maxuncalled = "0.05",
426 scoremethod = "2",
427 empirical_freqs = "yes",
428 truelenfile = "16S_101bp_len.txt",
429 testname = "101bpt_t1",
430 testw = "1")
431 num150bp_t2 = test_pear(forward = "16S_150bp1.fq",
432 reverse = "16S_150bp2.fq",
433 output = "to1",
434 pvalue = "0.01",
435 minoverlap = "1",
436 maxlen = "500",
437 minlen = "50",
438 mintrimlen = "50",
439 minquality = "10",
440 maxuncalled = "0.05",
441 scoremethod = "2",
442 empirical_freqs = "yes",
443 truelenfile = "16S_150bp_len.txt",
444 testname = "150bp_t2",
445 testw = "2")
446 num165bp_t1_o10 = test_pear(forward = "16S_165bp1.fq",
447 reverse = "16S_165bp2.fq",
448 output = "to1",
449 pvalue = "0.01",
450 minoverlap = "10",
451 maxlen = "500",
452 minlen = "50",
453 mintrimlen = "50",
454 minquality = "10",
455 maxuncalled = "0.05",
456 scoremethod = "2",
457 empirical_freqs = "yes",
458 truelenfile = "16S_165bp_len.txt",
459 testname = "165bp_t1_o10",
460 testw = "1")
463 num180bp_t1 = test_pear(forward = "16S_180bp1.fq",
464 reverse = "16S_180bp2.fq",
465 output = "to1",
466 pvalue = "0.01",
467 minoverlap = "1",
468 maxlen = "500",
469 minlen = "50",
470 mintrimlen = "50",
471 minquality = "10",
472 maxuncalled = "0.05",
473 scoremethod = "2",
474 empirical_freqs = "yes",
475 truelenfile = "16S_180bp_len.txt",
476 testname = "180bp_t1",
477 testw = "1")
479 num180bp_t1_001 = test_pear(forward = "16S_180bp1.fq",
480 reverse = "16S_180bp2.fq",
481 output = "to1",
482 pvalue = "0.001",
483 minoverlap = "1",
484 maxlen = "500",
485 minlen = "50",
486 mintrimlen = "50",
487 minquality = "10",
488 maxuncalled = "0.05",
489 scoremethod = "2",
490 empirical_freqs = "yes",
491 truelenfile = "16S_180bp_len.txt",
492 testname = "180bp_t1_001",
493 testw = "1")
495 num180bp_t2_001 = test_pear(forward = "16S_180bp1.fq",
496 reverse = "16S_180bp2.fq",
497 output = "to1",
498 pvalue = "0.001",
499 minoverlap = "1",
500 maxlen = "500",
501 minlen = "50",
502 mintrimlen = "50",
503 minquality = "10",
504 maxuncalled = "0.05",
505 scoremethod = "2",
506 empirical_freqs = "yes",
507 truelenfile = "16S_180bp_len.txt",
508 testname = "180bp_t2_001",
509 testw = "2")
511 num190bp_nt = test_pear(forward = "16S_190bp1.fq",
512 reverse = "16S_190bp2.fq",
513 output = "to1",
514 pvalue = "1.0",
515 minoverlap = "1",
516 maxlen = "500",
517 minlen = "50",
518 mintrimlen = "50",
519 minquality = "10",
520 maxuncalled = "0.05",
521 scoremethod = "2",
522 empirical_freqs = "yes",
523 truelenfile = "16S_190bp_len.txt",
524 testname = "190bp_nt",
525 testw = "2")
526 num250bp_t2 = test_pear(forward = "16S_250bp1.fq",
527 reverse = "16S_250bp2.fq",
528 output = "to1",
529 pvalue = "0.01",
530 minoverlap = "1",
531 maxlen = "500",
532 minlen = "50",
533 mintrimlen = "50",
534 minquality = "10",
535 maxuncalled = "0.05",
536 scoremethod = "2",
537 empirical_freqs = "yes",
538 truelenfile = "16S_250bp_len.txt",
539 testname = "250bp_t2",
540 testw = "2")
542 os.remove("to1.assembled.fastq")
543 os.remove("to1.discarded.fastq")
544 os.remove("to1.unassembled.forward.fastq")
545 os.remove("to1.unassembled.reverse.fastq")
548 if (num101bpt_t1 == 33022 and
549 num150bp_t2 == 28228 and
550 num165bp_t1_o10 == 25817 and
551 num180bp_t1 == 18115 and
552 num180bp_t1_001 == 14679 and
553 num180bp_t2_001 == 15532 and
554 num190bp_nt == 17112 and
555 num250bp_t2 == 23063):
556 print("Tests passed!")
557 else:
558 print("Warnning: tests failed!")
561 main()
562 #publication_data()