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