Mercurial > repos > peterjc > mira_assembler

comparison tools/sr_assembly/mira.py @ 0:03b240624b5a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Migrated tool version 0.0.1 from old tool shed archive to new tool shed repository
author peterjc
date Tue, 07 Jun 2011 16:23:51 -0400
parents
children 298f5c1d9521
comparison
equal deleted inserted replaced
-1:000000000000 0:03b240624b5a
1 #!/usr/bin/env python
2 """A simple wrapper script to call MIRA and collect its output.
3 """
4 import os
5 import sys
6 import subprocess
7 import shutil
8 import time
9
10 def stop_err(msg, err=1):
11 sys.stderr.write(msg+"\n")
12 sys.exit(err)
13
14 def tcs_to_tabular(old, new):
15 in_handle = open(old, "rU")
16 out_handle = open(new, "w")
17 assert in_handle.readline() == "#TCS V1.0\n"
18 assert in_handle.readline() == "#\n"
19 assert in_handle.readline() == "# contig name padPos upadPos | B Q | tcov covA covC covG covT cov* | qA qC qG qT q* | S | Tags\n"
20 assert in_handle.readline() == "#\n"
21 out_handle.write("#%s\n" % "\t".join(["contig", "pasPos", "upadPos", "B", "Q",
22 "tcov", "covA", "covC", "covG", "covT", "cov*",
23 "qA", "qC", "qG", "qT", "q*", "S", "Tags"]))
24 for line in in_handle:
25 parts = line.rstrip("\n").split(None,22)
26 assert parts[3] == parts[6] == parts[13] == parts[19] == parts[21] == "|"
27 wanted = parts[:3] + parts[4:6]+parts[7:13]+parts[14:19]+parts[20:21]+parts[22:]
28 out_handle.write("%s\n" % "\t".join(wanted))
29 out_handle.close()
30 in_handle.close()
31
32 def collect_output(temp, name):
33 n3 = (temp, name, name, name)
34 tcs_to_tabular("%s/%s_assembly/%s_d_results/%s_out.tcs" % n3, out_tcs)
35 for old, new in [("%s/%s_assembly/%s_d_results/%s_out.unpadded.fasta" % n3, out_fasta),
36 ("%s/%s_assembly/%s_d_results/%s_out.unpadded.fasta.qual" % n3, out_qual),
37 ("%s/%s_assembly/%s_d_results/%s_out.wig" % n3, out_wig),
38 ("%s/%s_assembly/%s_d_results/%s_out.caf" % n3, out_caf),
39 ("%s/%s_assembly/%s_d_results/%s_out.ace" % n3, out_ace)]:
40 if not os.path.isfile(old):
41 stop_err("Missing %s output file" % os.path.splitext(old)[-1])
42 else:
43 shutil.move(old, new)
44
45 def clean_up(temp, name):
46 folder = "%s/%s_assembly" % (temp, name)
47 if os.path.isdir(folder):
48 shutil.rmtree(folder)
49
50 #TODO - Run MIRA in /tmp or a configurable directory?
51 temp = "."
52 name, out_fasta, out_qual, out_tcs, out_ace, out_caf, out_wig, out_log = sys.argv[1:9]
53 start_time = time.time()
54 try:
55 cmd = " ".join(sys.argv[9:])
56 child = subprocess.Popen(sys.argv[9:],
57 stdout=subprocess.PIPE, stderr=subprocess.STDOUT)
58 except Exception, err:
59 sys.stderr.write("Error invoking command:\n%s\n\n%s\n" % (cmd, err))
60 sys.exit(1)
61 #Use .communicate as can get deadlocks with .wait(),
62 stdout, stderr = child.communicate()
63 run_time = time.time() - start_time
64 return_code = child.returncode
65 handle = open(out_log, "w")
66 handle.write(stdout)
67 handle.write("\n\nMIRA took %0.2f minutes\n" % (run_time / 60.0))
68 print "MIRA took %0.2f minutes" % (run_time / 60.0)
69 if return_code:
70 handle.write("Return error code %i from command:\n" % return_code)
71 handle.write(cmd + "\n")
72 handle.close()
73 clean_up(temp, name)
74 stop_err("Return error code %i from command:\n%s" % (return_code, cmd),
75 return_code)
76 handle.close()
77
78 collect_output(temp, name)
79 clean_up(temp, name)
80 print "Done"