|
0
|
1 #!/usr/bin/env python
|
|
|
2 """A simple wrapper script to call MIRA and collect its output.
|
|
|
3
|
|
|
4 This focuses on the miraconvert binary.
|
|
|
5 """
|
|
|
6 import os
|
|
|
7 import sys
|
|
|
8 import subprocess
|
|
|
9 import shutil
|
|
|
10 import time
|
|
|
11 import tempfile
|
|
|
12 from optparse import OptionParser
|
|
|
13 try:
|
|
|
14 from io import BytesIO
|
|
|
15 except ImportError:
|
|
|
16 #Should we worry about Python 2.5 or older?
|
|
|
17 from StringIO import StringIO as BytesIO
|
|
|
18
|
|
|
19 #Do we need any PYTHONPATH magic?
|
|
|
20 from mira4_make_bam import depad
|
|
|
21
|
|
|
22 WRAPPER_VER = "0.0.5" #Keep in sync with the XML file
|
|
|
23
|
|
|
24 def stop_err(msg, err=1):
|
|
|
25 sys.stderr.write(msg+"\n")
|
|
|
26 sys.exit(err)
|
|
|
27
|
|
|
28 def run(cmd):
|
|
|
29 #Avoid using shell=True when we call subprocess to ensure if the Python
|
|
|
30 #script is killed, so too is the child process.
|
|
|
31 try:
|
|
|
32 child = subprocess.Popen(cmd, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
|
|
|
33 except Exception, err:
|
|
|
34 stop_err("Error invoking command:\n%s\n\n%s\n" % (" ".join(cmd), err))
|
|
|
35 #Use .communicate as can get deadlocks with .wait(),
|
|
|
36 stdout, stderr = child.communicate()
|
|
|
37 return_code = child.returncode
|
|
|
38 if return_code:
|
|
|
39 if stderr and stdout:
|
|
|
40 stop_err("Return code %i from command:\n%s\n\n%s\n\n%s" % (return_code, err, stdout, stderr))
|
|
|
41 else:
|
|
|
42 stop_err("Return code %i from command:\n%s\n%s" % (return_code, err, stderr))
|
|
|
43
|
|
|
44 def get_version(mira_binary):
|
|
|
45 """Run MIRA to find its version number"""
|
|
|
46 # At the commend line I would use: mira -v | head -n 1
|
|
|
47 # however there is some pipe error when doing that here.
|
|
|
48 cmd = [mira_binary, "-v"]
|
|
|
49 try:
|
|
|
50 child = subprocess.Popen(cmd,
|
|
|
51 stdout=subprocess.PIPE,
|
|
|
52 stderr=subprocess.STDOUT)
|
|
|
53 except Exception, err:
|
|
|
54 sys.stderr.write("Error invoking command:\n%s\n\n%s\n" % (" ".join(cmd), err))
|
|
|
55 sys.exit(1)
|
|
|
56 ver, tmp = child.communicate()
|
|
|
57 del child
|
|
|
58 return ver.split("\n", 1)[0].strip()
|
|
|
59
|
|
|
60 #Parse Command Line
|
|
|
61 usage = """Galaxy MIRA4 wrapper script v%s - use as follows:
|
|
|
62
|
|
|
63 $ python mira4_convert.py ...
|
|
|
64
|
|
|
65 This will run the MIRA miraconvert binary and collect its output files as directed.
|
|
|
66 """ % WRAPPER_VER
|
|
|
67 parser = OptionParser(usage=usage)
|
|
|
68 parser.add_option("--input", dest="input",
|
|
|
69 default=None, metavar="FILE",
|
|
|
70 help="MIRA input filename")
|
|
|
71 parser.add_option("-x", "--min_length", dest="min_length",
|
|
|
72 default="0",
|
|
|
73 help="Minimum contig length")
|
|
|
74 parser.add_option("-y", "--min_cover", dest="min_cover",
|
|
|
75 default="0",
|
|
|
76 help="Minimum average contig coverage")
|
|
|
77 parser.add_option("-z", "--min_reads", dest="min_reads",
|
|
|
78 default="0",
|
|
|
79 help="Minimum reads per contig")
|
|
|
80 parser.add_option("--maf", dest="maf",
|
|
|
81 default="", metavar="FILE",
|
|
|
82 help="MIRA MAF output filename")
|
|
|
83 parser.add_option("--ace", dest="ace",
|
|
|
84 default="", metavar="FILE",
|
|
|
85 help="ACE output filename")
|
|
|
86 parser.add_option("--bam", dest="bam",
|
|
|
87 default="", metavar="FILE",
|
|
|
88 help="Unpadded BAM output filename")
|
|
|
89 parser.add_option("--fasta", dest="fasta",
|
|
|
90 default="", metavar="FILE",
|
|
|
91 help="Unpadded FASTA output filename")
|
|
|
92 parser.add_option("--cstats", dest="cstats",
|
|
|
93 default="", metavar="FILE",
|
|
|
94 help="Contig statistics filename")
|
|
|
95 parser.add_option("-v", "--version", dest="version",
|
|
|
96 default=False, action="store_true",
|
|
|
97 help="Show version and quit")
|
|
|
98 options, args = parser.parse_args()
|
|
|
99 if args:
|
|
|
100 stop_err("Expected options (e.g. --input example.maf), not arguments")
|
|
|
101
|
|
|
102 input_maf = options.input
|
|
|
103 out_maf = options.maf
|
|
|
104 out_bam = options.bam
|
|
|
105 out_fasta = options.fasta
|
|
|
106 out_ace = options.ace
|
|
|
107 out_cstats = options.cstats
|
|
|
108
|
|
|
109 try:
|
|
|
110 mira_path = os.environ["MIRA4"]
|
|
|
111 except KeyError:
|
|
|
112 stop_err("Environment variable $MIRA4 not set")
|
|
|
113 mira_convert = os.path.join(mira_path, "miraconvert")
|
|
|
114 if not os.path.isfile(mira_convert):
|
|
|
115 stop_err("Missing miraconvert under $MIRA4, %r\nFolder contained: %s"
|
|
|
116 % (mira_convert, ", ".join(os.listdir(mira_path))))
|
|
|
117
|
|
|
118 mira_convert_ver = get_version(mira_convert)
|
|
|
119 if not mira_convert_ver.strip().startswith("4.0"):
|
|
|
120 stop_err("This wrapper is for MIRA V4.0, not:\n%s\n%s" % (mira_ver, mira_convert))
|
|
|
121 if options.version:
|
|
|
122 print "%s, MIRA wrapper version %s" % (mira_convert_ver, WRAPPER_VER)
|
|
|
123 sys.exit(0)
|
|
|
124
|
|
|
125 if not input_maf:
|
|
|
126 stop_err("Input MIRA file is required")
|
|
|
127 elif not os.path.isfile(input_maf):
|
|
|
128 stop_err("Missing input MIRA file: %r" % input_maf)
|
|
|
129
|
|
|
130 if not (out_maf or out_bam or out_fasta or out_ace or out_cstats):
|
|
|
131 stop_err("No output requested")
|
|
|
132
|
|
|
133
|
|
|
134 def check_min_int(value, name):
|
|
|
135 try:
|
|
|
136 i = int(value)
|
|
|
137 except:
|
|
|
138 stop_err("Bad %s setting, %r" % (name, value))
|
|
|
139 if i < 0:
|
|
|
140 stop_err("Negative %s setting, %r" % (name, value))
|
|
|
141 return i
|
|
|
142
|
|
|
143 min_length = check_min_int(options.min_length, "minimum length")
|
|
|
144 min_cover = check_min_int(options.min_cover, "minimum cover")
|
|
|
145 min_reads = check_min_int(options.min_reads, "minimum reads")
|
|
|
146
|
|
|
147 #TODO - Run MIRA in /tmp or a configurable directory?
|
|
|
148 #Currently Galaxy puts us somewhere safe like:
|
|
|
149 #/opt/galaxy-dist/database/job_working_directory/846/
|
|
|
150 temp = "."
|
|
|
151
|
|
|
152
|
|
|
153 cmd_list = [mira_convert]
|
|
|
154 if min_length:
|
|
|
155 cmd_list.extend(["-x", str(min_length)])
|
|
|
156 if min_cover:
|
|
|
157 cmd_list.extend(["-y", str(min_cover)])
|
|
|
158 if min_reads:
|
|
|
159 cmd_list.extend(["-z", str(min_reads)])
|
|
|
160 cmd_list.extend(["-f", "maf", input_maf, os.path.join(temp, "converted")])
|
|
|
161 if out_maf:
|
|
|
162 cmd_list.append("maf")
|
|
|
163 if out_bam:
|
|
|
164 cmd_list.append("samnbb")
|
|
|
165 if not out_fasta:
|
|
|
166 #Need this for samtools depad
|
|
|
167 out_fasta = os.path.join(temp, "depadded.fasta")
|
|
|
168 if out_fasta:
|
|
|
169 cmd_list.append("fasta")
|
|
|
170 if out_ace:
|
|
|
171 cmd_list.append("ace")
|
|
|
172 if out_cstats:
|
|
|
173 cmd_list.append("cstats")
|
|
|
174 run(cmd_list)
|
|
|
175
|
|
|
176 def collect(old, new):
|
|
|
177 if not os.path.isfile(old):
|
|
|
178 stop_err("Missing expected output file %s" % old)
|
|
|
179 shutil.move(old, new)
|
|
|
180
|
|
|
181 if out_maf:
|
|
|
182 collect(os.path.join(temp, "converted.maf"), out_maf)
|
|
|
183 if out_fasta:
|
|
|
184 #Can we look at the MAF file to see if there are multiple strains?
|
|
|
185 old = os.path.join(temp, "converted_AllStrains.unpadded.fasta")
|
|
|
186 if os.path.isfile(old):
|
|
|
187 collect(old, out_fasta)
|
|
|
188 else:
|
|
|
189 #Might the output be filtered down to zero contigs?
|
|
|
190 old = os.path.join(temp, "converted.fasta")
|
|
|
191 if not os.path.isfile(old):
|
|
|
192 stop_err("Missing expected output FASTA file")
|
|
|
193 elif os.path.getsize(old) == 0:
|
|
|
194 print("Warning - no contigs (harsh filters?)")
|
|
|
195 collect(old, out_fasta)
|
|
|
196 else:
|
|
|
197 stop_err("Missing expected output FASTA file (only generic file present)")
|
|
|
198 if out_ace:
|
|
|
199 collect(os.path.join(temp, "converted.maf"), out_ace)
|
|
|
200 if out_cstats:
|
|
|
201 collect(os.path.join(temp, "converted_info_contigstats.txt"), out_cstats)
|
|
|
202
|
|
|
203 if out_bam:
|
|
|
204 assert os.path.isfile(out_fasta)
|
|
|
205 old = os.path.join(temp, "converted.samnbb")
|
|
|
206 if not os.path.isfile(old):
|
|
|
207 old = os.path.join(temp, "converted.sam")
|
|
|
208 if not os.path.isfile(old):
|
|
|
209 stop_err("Missing expected intermediate file %s" % old)
|
|
|
210 h = BytesIO()
|
|
|
211 msg = depad(out_fasta, old, out_bam, h)
|
|
|
212 if msg:
|
|
|
213 print(msg)
|
|
|
214 print(h.getvalue())
|
|
|
215 h.close()
|
|
|
216 sys.exit(1)
|
|
|
217 h.close()
|
|
|
218 if out_fasta == os.path.join(temp, "depadded.fasta"):
|
|
|
219 #Not asked for by Galaxy, no longer needed
|
|
|
220 os.remove(out_fasta)
|
|
|
221
|
|
|
222 if min_length or min_cover or min_reads:
|
|
|
223 print("Filtered.")
|
|
|
224 else:
|
|
|
225 print("Converted.")
|
