Mercurial > repos > peterjc > mira4_assembler

comparison tools/mira4/mira4.py @ 20:aeb3e35f8236 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded v0.0.4 preview, made MAF and BAM output optional
author peterjc
date Tue, 10 Jun 2014 10:11:58 -0400
parents 8487d70e82aa
children 4abe8d59a438
comparison
equal deleted inserted replaced
19:8487d70e82aa 20:aeb3e35f8236
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 import tempfile 9 import tempfile
10 from optparse import OptionParser
10 11
11 #Do we need any PYTHONPATH magic? 12 #Do we need any PYTHONPATH magic?
12 from mira4_make_bam import make_bam 13 from mira4_make_bam import make_bam
13 14
14 WRAPPER_VER = "0.0.1" #Keep in sync with the XML file 15 WRAPPER_VER = "0.0.4" #Keep in sync with the XML file
15 16
16 def stop_err(msg, err=1): 17 def stop_err(msg, err=1):
17 sys.stderr.write(msg+"\n") 18 sys.stderr.write(msg+"\n")
18 sys.exit(err) 19 sys.exit(err)
19 20
32 sys.exit(1) 33 sys.exit(1)
33 ver, tmp = child.communicate() 34 ver, tmp = child.communicate()
34 del child 35 del child
35 return ver.split("\n", 1)[0].strip() 36 return ver.split("\n", 1)[0].strip()
36 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 options, args = parser.parse_args()
62 manifest = options.manifest
63 out_maf = options.maf
64 out_bam = options.bam
65 out_fasta = options.fasta
66 out_log = options.log
37 67
38 try: 68 try:
39 mira_path = os.environ["MIRA4"] 69 mira_path = os.environ["MIRA4"]
40 except KeyError: 70 except KeyError:
41 stop_err("Environment variable $MIRA4 not set") 71 stop_err("Environment variable $MIRA4 not set")
58 print "%s, MIRA wrapper version %s" % (mira_ver, WRAPPER_VER) 88 print "%s, MIRA wrapper version %s" % (mira_ver, WRAPPER_VER)
59 if mira_ver != mira_convert_ver: 89 if mira_ver != mira_convert_ver:
60 print "WARNING: miraconvert %s" % mira_convert_ver 90 print "WARNING: miraconvert %s" % mira_convert_ver
61 sys.exit(0) 91 sys.exit(0)
62 92
93 if not manifest:
94 stop_err("Manifest is required")
95 elif not os.path.isfile(manifest):
96 stop_err("Missing input MIRA manifest file: %r" % manifest)
97
63 98
64 try: 99 try:
65 threads = int(os.environ.get("GALAXY_SLOTS", "1")) 100 threads = int(os.environ.get("GALAXY_SLOTS", "1"))
66 except ValueError: 101 except ValueError:
67 threads = 1 102 threads = 1
141 missing = False 176 missing = False
142 for old, new in [(old_maf, out_maf), 177 for old, new in [(old_maf, out_maf),
143 (old_fasta, out_fasta)]: 178 (old_fasta, out_fasta)]:
144 if not os.path.isfile(old): 179 if not os.path.isfile(old):
145 missing = True 180 missing = True
181 elif not new or new == "-":
182 handle.write("Ignoring %s\n" % old)
146 else: 183 else:
147 handle.write("Capturing %s\n" % old) 184 handle.write("Capturing %s\n" % old)
148 shutil.move(old, new) 185 shutil.move(old, new)
149 if missing: 186 if missing:
150 log_manifest(manifest) 187 log_manifest(manifest)
152 for filename in sorted(os.listdir(f)): 189 for filename in sorted(os.listdir(f)):
153 sys.stderr.write("%s\n" % filename) 190 sys.stderr.write("%s\n" % filename)
154 191
155 #For mapping mode, probably most people would expect a BAM file 192 #For mapping mode, probably most people would expect a BAM file
156 #using the reference FASTA file... 193 #using the reference FASTA file...
157 msg = make_bam(mira_convert, out_maf, ref_fasta, out_bam, handle) 194 if out_bam and out_bam != "-":
158 if msg: 195 msg = make_bam(mira_convert, out_maf, ref_fasta, out_bam, handle)
159 stop_err(msg) 196 if msg:
197 stop_err(msg)
160 198
161 def clean_up(temp, name): 199 def clean_up(temp, name):
162 folder = "%s/%s_assembly" % (temp, name) 200 folder = "%s/%s_assembly" % (temp, name)
163 if os.path.isdir(folder): 201 if os.path.isdir(folder):
164 shutil.rmtree(folder) 202 shutil.rmtree(folder)
165 203
166 #TODO - Run MIRA in /tmp or a configurable directory? 204 #TODO - Run MIRA in /tmp or a configurable directory?
167 #Currently Galaxy puts us somewhere safe like: 205 #Currently Galaxy puts us somewhere safe like:
168 #/opt/galaxy-dist/database/job_working_directory/846/ 206 #/opt/galaxy-dist/database/job_working_directory/846/
169 temp = "." 207 temp = "."
170 #name, out_fasta, out_qual, out_ace, out_caf, out_wig, out_log = sys.argv[1:8] 208
171 name = "MIRA" 209 name = "MIRA"
172 manifest, out_maf, out_bam, out_fasta, out_log = sys.argv[1:]
173 210
174 override_temp(manifest) 211 override_temp(manifest)
175 212
176 start_time = time.time() 213 start_time = time.time()
177 #cmd_list =sys.argv[8:]
178 cmd_list = [mira_binary, "-t", str(threads), manifest] 214 cmd_list = [mira_binary, "-t", str(threads), manifest]
179 cmd = " ".join(cmd_list) 215 cmd = " ".join(cmd_list)
180 216
181 assert os.path.isdir(temp) 217 assert os.path.isdir(temp)
182 d = "%s_assembly" % name 218 d = "%s_assembly" % name
190 sys.exit(1) 226 sys.exit(1)
191 227
192 #print os.path.abspath(".") 228 #print os.path.abspath(".")
193 #print cmd 229 #print cmd
194 230
195 handle = open(out_log, "w") 231 if out_log and out_log != "-":
232 handle = open(out_log, "w")
233 else:
234 handle = open(os.devnull, "w")
196 handle.write("======================== MIRA manifest (instructions) ========================\n") 235 handle.write("======================== MIRA manifest (instructions) ========================\n")
197 m = open(manifest, "rU") 236 m = open(manifest, "rU")
198 for line in m: 237 for line in m:
199 handle.write(line) 238 handle.write(line)
200 m.close() 239 m.close()