Mercurial > repos > peterjc > mira4_assembler

comparison tools/mira4/mira4.py @ 0:6a88b42ce6b9 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded v0.0.4, previously only on the TestToolShed
author peterjc
date Fri, 21 Nov 2014 06:42:56 -0500
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:6a88b42ce6b9
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 import tempfile
10 from optparse import OptionParser
11
12 #Do we need any PYTHONPATH magic?
13 from mira4_make_bam import make_bam
14
15 WRAPPER_VER = "0.0.4" #Keep in sync with the XML file
16
17 def stop_err(msg, err=1):
18 sys.stderr.write(msg+"\n")
19 sys.exit(err)
20
21
22 def get_version(mira_binary):
23 """Run MIRA to find its version number"""
24 # At the commend line I would use: mira -v | head -n 1
25 # however there is some pipe error when doing that here.
26 cmd = [mira_binary, "-v"]
27 try:
28 child = subprocess.Popen(cmd,
29 stdout=subprocess.PIPE,
30 stderr=subprocess.STDOUT)
31 except Exception, err:
32 sys.stderr.write("Error invoking command:\n%s\n\n%s\n" % (" ".join(cmd), err))
33 sys.exit(1)
34 ver, tmp = child.communicate()
35 del child
36 return ver.split("\n", 1)[0].strip()
37
38 #Parse Command Line
39 usage = """Galaxy MIRA4 wrapper script v%s - use as follows:
40
41 $ python mira4.py ...
42
43 This will run the MIRA binary and collect its output files as directed.
44 """ % WRAPPER_VER
45 parser = OptionParser(usage=usage)
46 parser.add_option("-m", "--manifest", dest="manifest",
47 default=None, metavar="FILE",
48 help="MIRA manifest filename")
49 parser.add_option("--maf", dest="maf",
50 default="-", metavar="FILE",
51 help="MIRA MAF output filename")
52 parser.add_option("--bam", dest="bam",
53 default="-", metavar="FILE",
54 help="Unpadded BAM output filename")
55 parser.add_option("--fasta", dest="fasta",
56 default="-", metavar="FILE",
57 help="Unpadded FASTA output filename")
58 parser.add_option("--log", dest="log",
59 default="-", metavar="FILE",
60 help="MIRA logging output filename")
61 parser.add_option("-v", "--version", dest="version",
62 default=False, action="store_true",
63 help="Show version and quit")
64 options, args = parser.parse_args()
65 manifest = options.manifest
66 out_maf = options.maf
67 out_bam = options.bam
68 out_fasta = options.fasta
69 out_log = options.log
70
71 try:
72 mira_path = os.environ["MIRA4"]
73 except KeyError:
74 stop_err("Environment variable $MIRA4 not set")
75 mira_binary = os.path.join(mira_path, "mira")
76 if not os.path.isfile(mira_binary):
77 stop_err("Missing mira under $MIRA4, %r\nFolder contained: %s"
78 % (mira_binary, ", ".join(os.listdir(mira_path))))
79 mira_convert = os.path.join(mira_path, "miraconvert")
80 if not os.path.isfile(mira_convert):
81 stop_err("Missing miraconvert under $MIRA4, %r\nFolder contained: %s"
82 % (mira_convert, ", ".join(os.listdir(mira_path))))
83
84 mira_ver = get_version(mira_binary)
85 if not mira_ver.strip().startswith("4.0"):
86 stop_err("This wrapper is for MIRA V4.0, not:\n%s\n%s" % (mira_ver, mira_binary))
87 mira_convert_ver = get_version(mira_convert)
88 if not mira_convert_ver.strip().startswith("4.0"):
89 stop_err("This wrapper is for MIRA V4.0, not:\n%s\n%s" % (mira_ver, mira_convert))
90 if options.version:
91 print "%s, MIRA wrapper version %s" % (mira_ver, WRAPPER_VER)
92 if mira_ver != mira_convert_ver:
93 print "WARNING: miraconvert %s" % mira_convert_ver
94 sys.exit(0)
95
96 if not manifest:
97 stop_err("Manifest is required")
98 elif not os.path.isfile(manifest):
99 stop_err("Missing input MIRA manifest file: %r" % manifest)
100
101
102 try:
103 threads = int(os.environ.get("GALAXY_SLOTS", "1"))
104 except ValueError:
105 threads = 1
106 assert 1 <= threads, threads
107
108
109 def override_temp(manifest):
110 """Override ``-DI:trt=/tmp`` in manifest with environment variable.
111
112 Currently MIRA 4 does not allow envronment variables like ``$TMP``
113 inside the manifest, which is a problem if you need to override
114 the default at run time.
115
116 The tool XML will ``/tmp`` and we replace that here with
117 ``tempfile.gettempdir()`` which will respect $TMPDIR, $TEMP, $TMP
118 as explained in the Python standard library documentation:
119 http://docs.python.org/2/library/tempfile.html#tempfile.tempdir
120
121 By default MIRA 4 would write its temporary files within the output
122 folder, which is a problem if that is a network drive.
123 """
124 handle = open(manifest, "r")
125 text = handle.read()
126 handle.close()
127
128 #At time of writing, this is at the end of a file,
129 #but could be followed by a space in future...
130 text = text.replace("-DI:trt=/tmp", "-DI:trt=" + tempfile.gettempdir())
131
132 #Want to try to ensure this gets written to disk before MIRA attempts
133 #to open it - any networked file system may impose a delay...
134 handle = open(manifest, "w")
135 handle.write(text)
136 handle.flush()
137 os.fsync(handle.fileno())
138 handle.close()
139
140
141 def log_manifest(manifest):
142 """Write the manifest file to stderr."""
143 sys.stderr.write("\n%s\nManifest file\n%s\n" % ("="*60, "="*60))
144 with open(manifest) as h:
145 for line in h:
146 sys.stderr.write(line)
147 sys.stderr.write("\n%s\nEnd of manifest\n%s\n" % ("="*60, "="*60))
148
149
150 def collect_output(temp, name, handle):
151 """Moves files to the output filenames (global variables)."""
152 n3 = (temp, name, name, name)
153 f = "%s/%s_assembly/%s_d_results" % (temp, name, name)
154 if not os.path.isdir(f):
155 log_manifest(manifest)
156 stop_err("Missing output folder")
157 if not os.listdir(f):
158 log_manifest(manifest)
159 stop_err("Empty output folder")
160 missing = []
161
162 old_maf = "%s/%s_out.maf" % (f, name)
163 if not os.path.isfile(old_maf):
164 #Triggered extractLargeContigs.sh?
165 old_maf = "%s/%s_LargeContigs_out.maf" % (f, name)
166
167 #De novo or single strain mapping,
168 old_fasta = "%s/%s_out.unpadded.fasta" % (f, name)
169 ref_fasta = "%s/%s_out.padded.fasta" % (f, name)
170 if not os.path.isfile(old_fasta):
171 #Mapping (StrainX versus reference) or de novo
172 old_fasta = "%s/%s_out_StrainX.unpadded.fasta" % (f, name)
173 ref_fasta = "%s/%s_out_StrainX.padded.fasta" % (f, name)
174 if not os.path.isfile(old_fasta):
175 old_fasta = "%s/%s_out_ReferenceStrain.unpadded.fasta" % (f, name)
176 ref_fasta = "%s/%s_out_ReferenceStrain.padded.fasta" % (f, name)
177
178
179 missing = False
180 for old, new in [(old_maf, out_maf),
181 (old_fasta, out_fasta)]:
182 if not os.path.isfile(old):
183 missing = True
184 elif not new or new == "-":
185 handle.write("Ignoring %s\n" % old)
186 else:
187 handle.write("Capturing %s\n" % old)
188 shutil.move(old, new)
189 if missing:
190 log_manifest(manifest)
191 sys.stderr.write("Contents of %r:\n" % f)
192 for filename in sorted(os.listdir(f)):
193 sys.stderr.write("%s\n" % filename)
194
195 #For mapping mode, probably most people would expect a BAM file
196 #using the reference FASTA file...
197 if out_bam and out_bam != "-":
198 if out_maf and out_maf != "-":
199 msg = make_bam(mira_convert, out_maf, ref_fasta, out_bam, handle)
200 else:
201 #Not collecting the MAF file, use original location
202 msg = make_bam(mira_convert, old_maf, ref_fasta, out_bam, handle)
203 if msg:
204 stop_err(msg)
205
206 def clean_up(temp, name):
207 folder = "%s/%s_assembly" % (temp, name)
208 if os.path.isdir(folder):
209 shutil.rmtree(folder)
210
211 #TODO - Run MIRA in /tmp or a configurable directory?
212 #Currently Galaxy puts us somewhere safe like:
213 #/opt/galaxy-dist/database/job_working_directory/846/
214 temp = "."
215
216 name = "MIRA"
217
218 override_temp(manifest)
219
220 start_time = time.time()
221 cmd_list = [mira_binary, "-t", str(threads), manifest]
222 cmd = " ".join(cmd_list)
223
224 assert os.path.isdir(temp)
225 d = "%s_assembly" % name
226 #This can fail on my development machine if stale folders exist
227 #under Galaxy's .../database/job_working_directory/ tree:
228 assert not os.path.isdir(d), "Path %r already exists:\n%s" % (d, os.path.abspath(d))
229 try:
230 #Check path access
231 os.mkdir(d)
232 except Exception, err:
233 log_manifest(manifest)
234 sys.stderr.write("Error making directory %s\n%s" % (d, err))
235 sys.exit(1)
236
237 #print os.path.abspath(".")
238 #print cmd
239
240 if out_log and out_log != "-":
241 handle = open(out_log, "w")
242 else:
243 handle = open(os.devnull, "w")
244 handle.write("======================== MIRA manifest (instructions) ========================\n")
245 m = open(manifest, "rU")
246 for line in m:
247 handle.write(line)
248 m.close()
249 del m
250 handle.write("\n")
251 handle.write("============================ Starting MIRA now ===============================\n")
252 handle.flush()
253 try:
254 #Run MIRA
255 child = subprocess.Popen(cmd_list,
256 stdout=handle,
257 stderr=subprocess.STDOUT)
258 except Exception, err:
259 log_manifest(manifest)
260 sys.stderr.write("Error invoking command:\n%s\n\n%s\n" % (cmd, err))
261 #TODO - call clean up?
262 handle.write("Error invoking command:\n%s\n\n%s\n" % (cmd, err))
263 handle.close()
264 sys.exit(1)
265 #Use .communicate as can get deadlocks with .wait(),
266 stdout, stderr = child.communicate()
267 assert not stdout and not stderr #Should be empty as sent to handle
268 run_time = time.time() - start_time
269 return_code = child.returncode
270 handle.write("\n")
271 handle.write("============================ MIRA has finished ===============================\n")
272 handle.write("MIRA took %0.2f hours\n" % (run_time / 3600.0))
273 if return_code:
274 print "MIRA took %0.2f hours" % (run_time / 3600.0)
275 handle.write("Return error code %i from command:\n" % return_code)
276 handle.write(cmd + "\n")
277 handle.close()
278 clean_up(temp, name)
279 log_manifest(manifest)
280 stop_err("Return error code %i from command:\n%s" % (return_code, cmd),
281 return_code)
282 handle.flush()
283
284 if os.path.isfile("MIRA_assembly/MIRA_d_results/ec.log"):
285 handle.write("\n")
286 handle.write("====================== Extract Large Contigs failed ==========================\n")
287 e = open("MIRA_assembly/MIRA_d_results/ec.log", "rU")
288 for line in e:
289 handle.write(line)
290 e.close()
291 handle.write("============================ (end of ec.log) =================================\n")
292 handle.flush()
293
294 #print "Collecting output..."
295 start_time = time.time()
296 collect_output(temp, name, handle)
297 collect_time = time.time() - start_time
298 handle.write("MIRA took %0.2f hours; collecting output %0.2f minutes\n" % (run_time / 3600.0, collect_time / 60.0))
299 print("MIRA took %0.2f hours; collecting output %0.2f minutes\n" % (run_time / 3600.0, collect_time / 60.0))
300
301 if os.path.isfile("MIRA_assembly/MIRA_d_results/ec.log"):
302 #Treat as an error, but doing this AFTER collect_output
303 sys.stderr.write("Extract Large Contigs failed\n")
304 handle.write("Extract Large Contigs failed\n")
305 handle.close()
306 sys.exit(1)
307
308 #print "Cleaning up..."
309 clean_up(temp, name)
310
311 handle.write("\nDone\n")
312 handle.close()
313 print("Done")