Mercurial > repos > peterjc > mira_assembler
comparison tools/sr_assembly/mira.py @ 1:c947750f82fb draft
Uploaded v0.0.5
| author | peterjc |
|---|---|
| date | Wed, 24 Apr 2013 11:51:50 -0400 |
| parents | d0a5acdf1638 |
| children |
comparison
equal
deleted
inserted
replaced
| 0:d0a5acdf1638 | 1:c947750f82fb |
|---|---|
| 5 import sys | 5 import sys |
| 6 import subprocess | 6 import subprocess |
| 7 import shutil | 7 import shutil |
| 8 import time | 8 import time |
| 9 | 9 |
| 10 WRAPPER_VER = "0.0.5" #Keep in sync with the XML file | |
| 11 | |
| 10 def stop_err(msg, err=1): | 12 def stop_err(msg, err=1): |
| 11 sys.stderr.write(msg+"\n") | 13 sys.stderr.write(msg+"\n") |
| 12 sys.exit(err) | 14 sys.exit(err) |
| 13 | 15 |
| 14 def tcs_to_tabular(old, new): | 16 |
| 15 in_handle = open(old, "rU") | 17 def get_version(): |
| 16 out_handle = open(new, "w") | 18 """Run MIRA to find its version number""" |
| 17 assert in_handle.readline() == "#TCS V1.0\n" | 19 # At the commend line I would use: mira -v | head -n 1 |
| 18 assert in_handle.readline() == "#\n" | 20 # however there is some pipe error when doing that here. |
| 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" | 21 try: |
| 20 assert in_handle.readline() == "#\n" | 22 child = subprocess.Popen(["mira", "-v"], |
| 21 out_handle.write("#%s\n" % "\t".join(["contig", "pasPos", "upadPos", "B", "Q", | 23 stdout=subprocess.PIPE, |
| 22 "tcov", "covA", "covC", "covG", "covT", "cov*", | 24 stderr=subprocess.STDOUT) |
| 23 "qA", "qC", "qG", "qT", "q*", "S", "Tags"])) | 25 except Exception, err: |
| 24 for line in in_handle: | 26 sys.stderr.write("Error invoking command:\n%s\n\n%s\n" % (" ".join(cmd), err)) |
| 25 parts = line.rstrip("\n").split(None,22) | 27 sys.exit(1) |
| 26 assert parts[3] == parts[6] == parts[13] == parts[19] == parts[21] == "|" | 28 ver, tmp = child.communicate() |
| 27 wanted = parts[:3] + parts[4:6]+parts[7:13]+parts[14:19]+parts[20:21]+parts[22:] | 29 del child |
| 28 out_handle.write("%s\n" % "\t".join(wanted)) | 30 return ver.split("\n", 1)[0] |
| 29 out_handle.close() | 31 |
| 30 in_handle.close() | 32 |
| 33 mira_ver = get_version() | |
| 34 if "V3.4." not in mira_ver: | |
| 35 stop_err("This wrapper is for MIRA V3.4, not %s" % mira_ver) | |
| 36 if "-v" in sys.argv: | |
| 37 print "MIRA wrapper version %s," % WRAPPER_VER | |
| 38 print mira_ver | |
| 39 sys.exit(0) | |
| 40 | |
| 31 | 41 |
| 32 def collect_output(temp, name): | 42 def collect_output(temp, name): |
| 33 n3 = (temp, name, name, name) | 43 n3 = (temp, name, name, name) |
| 34 tcs_to_tabular("%s/%s_assembly/%s_d_results/%s_out.tcs" % n3, out_tcs) | 44 f = "%s/%s_assembly/%s_d_results" % (temp, name, name) |
| 35 for old, new in [("%s/%s_assembly/%s_d_results/%s_out.unpadded.fasta" % n3, out_fasta), | 45 if not os.path.isdir(f): |
| 36 ("%s/%s_assembly/%s_d_results/%s_out.unpadded.fasta.qual" % n3, out_qual), | 46 stop_err("Missing output folder") |
| 37 ("%s/%s_assembly/%s_d_results/%s_out.wig" % n3, out_wig), | 47 if not os.listdir(f): |
| 38 ("%s/%s_assembly/%s_d_results/%s_out.caf" % n3, out_caf), | 48 stop_err("Empty output folder") |
| 39 ("%s/%s_assembly/%s_d_results/%s_out.ace" % n3, out_ace)]: | 49 missing = [] |
| 50 for old, new in [("%s/%s_out.unpadded.fasta" % (f, name), out_fasta), | |
| 51 ("%s/%s_out.unpadded.fasta.qual" % (f, name), out_qual), | |
| 52 ("%s/%s_out.wig" % (f, name), out_wig), | |
| 53 ("%s/%s_out.caf" % (f, name), out_caf), | |
| 54 ("%s/%s_out.ace" % (f, name), out_ace)]: | |
| 40 if not os.path.isfile(old): | 55 if not os.path.isfile(old): |
| 41 stop_err("Missing %s output file" % os.path.splitext(old)[-1]) | 56 missing.append(os.path.splitext(old)[-1]) |
| 42 else: | 57 else: |
| 43 shutil.move(old, new) | 58 shutil.move(old, new) |
| 59 if missing: | |
| 60 stop_err("Missing output files: %s" % ", ".join(missing)) | |
| 44 | 61 |
| 45 def clean_up(temp, name): | 62 def clean_up(temp, name): |
| 46 folder = "%s/%s_assembly" % (temp, name) | 63 folder = "%s/%s_assembly" % (temp, name) |
| 47 if os.path.isdir(folder): | 64 if os.path.isdir(folder): |
| 48 shutil.rmtree(folder) | 65 shutil.rmtree(folder) |
| 49 | 66 |
| 50 #TODO - Run MIRA in /tmp or a configurable directory? | 67 #TODO - Run MIRA in /tmp or a configurable directory? |
| 68 #Currently Galaxy puts us somewhere safe like: | |
| 69 #/opt/galaxy-dist/database/job_working_directory/846/ | |
| 51 temp = "." | 70 temp = "." |
| 52 name, out_fasta, out_qual, out_tcs, out_ace, out_caf, out_wig, out_log = sys.argv[1:9] | 71 name, out_fasta, out_qual, out_ace, out_caf, out_wig, out_log = sys.argv[1:8] |
| 72 | |
| 53 start_time = time.time() | 73 start_time = time.time() |
| 74 cmd_list =sys.argv[8:] | |
| 75 cmd = " ".join(cmd_list) | |
| 76 | |
| 77 assert os.path.isdir(temp) | |
| 78 d = "%s_assembly" % name | |
| 79 assert not os.path.isdir(d), "Path %s already exists" % d | |
| 54 try: | 80 try: |
| 55 cmd = " ".join(sys.argv[9:]) | 81 #Check path access |
| 56 child = subprocess.Popen(sys.argv[9:], | 82 os.mkdir(d) |
| 57 stdout=subprocess.PIPE, stderr=subprocess.STDOUT) | 83 except Exception, err: |
| 84 sys.stderr.write("Error making directory %s\n%s" % (d, err)) | |
| 85 sys.exit(1) | |
| 86 | |
| 87 #print os.path.abspath(".") | |
| 88 #print cmd | |
| 89 | |
| 90 handle = open(out_log, "w") | |
| 91 try: | |
| 92 #Run MIRA | |
| 93 child = subprocess.Popen(cmd_list, | |
| 94 stdout=handle, | |
| 95 stderr=subprocess.STDOUT) | |
| 58 except Exception, err: | 96 except Exception, err: |
| 59 sys.stderr.write("Error invoking command:\n%s\n\n%s\n" % (cmd, err)) | 97 sys.stderr.write("Error invoking command:\n%s\n\n%s\n" % (cmd, err)) |
| 98 #TODO - call clean up? | |
| 99 handle.write("Error invoking command:\n%s\n\n%s\n" % (cmd, err)) | |
| 100 handle.close() | |
| 60 sys.exit(1) | 101 sys.exit(1) |
| 61 #Use .communicate as can get deadlocks with .wait(), | 102 #Use .communicate as can get deadlocks with .wait(), |
| 62 stdout, stderr = child.communicate() | 103 stdout, stderr = child.communicate() |
| 104 assert not stdout and not stderr #Should be empty as sent to handle | |
| 63 run_time = time.time() - start_time | 105 run_time = time.time() - start_time |
| 64 return_code = child.returncode | 106 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)) | 107 handle.write("\n\nMIRA took %0.2f minutes\n" % (run_time / 60.0)) |
| 68 print "MIRA took %0.2f minutes" % (run_time / 60.0) | 108 print "MIRA took %0.2f minutes" % (run_time / 60.0) |
| 69 if return_code: | 109 if return_code: |
| 70 handle.write("Return error code %i from command:\n" % return_code) | 110 handle.write("Return error code %i from command:\n" % return_code) |
| 71 handle.write(cmd + "\n") | 111 handle.write(cmd + "\n") |
| 73 clean_up(temp, name) | 113 clean_up(temp, name) |
| 74 stop_err("Return error code %i from command:\n%s" % (return_code, cmd), | 114 stop_err("Return error code %i from command:\n%s" % (return_code, cmd), |
| 75 return_code) | 115 return_code) |
| 76 handle.close() | 116 handle.close() |
| 77 | 117 |
| 118 #print "Collecting output..." | |
| 78 collect_output(temp, name) | 119 collect_output(temp, name) |
| 120 | |
| 121 #print "Cleaning up..." | |
| 79 clean_up(temp, name) | 122 clean_up(temp, name) |
| 123 | |
| 80 print "Done" | 124 print "Done" |