Mercurial > repos > peterjc > mira4_assembler

comparison tools/mira4/mira4_convert.py @ 21:4abe8d59a438 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded v0.0.4 preview 1; fix getting BAM without MAF
author peterjc
date Tue, 28 Oct 2014 08:29:59 -0400
parents
children
comparison
equal deleted inserted replaced
20:aeb3e35f8236 21:4abe8d59a438
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.")