Mercurial > repos > peterjc > mira4_assembler

comparison tools/mira4/mira4_convert.py @ 0:6a88b42ce6b9 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded v0.0.4, previously only on the TestToolShed
author peterjc
date Fri, 21 Nov 2014 06:42:56 -0500
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:6a88b42ce6b9
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.")