Mercurial > repos > peterjc > mira4_assembler

comparison tools/mira4_0/mira4.py @ 2:4eb32a3d67d1 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
v0.0.8 - renamed folder, added note about mirabait
author peterjc
date Wed, 02 Sep 2015 07:46:29 -0400
parents
children 1713289d9908
comparison
equal deleted inserted replaced
1:70248e6e3efc 2:4eb32a3d67d1
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 sys_exit(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 sys_exit("Environment variable $MIRA4 not set")
75 mira_binary = os.path.join(mira_path, "mira")
76 if not os.path.isfile(mira_binary):
77 sys_exit("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 sys_exit("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 sys_exit("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 sys_exit("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 sys_exit("Manifest is required")
98 elif not os.path.isfile(manifest):
99 sys_exit("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 sys_exit("Missing output folder")
157 if not os.listdir(f):
158 log_manifest(manifest)
159 sys_exit("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 sys_exit(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 sys_exit("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")