|
0
|
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")
|
