Mercurial > repos > peterjc > mira4_assembler

comparison tools/mira4/mira4_bait.py @ 9:302d13490b23 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded v0.0.2 preview 1, BAM output
author peterjc
date Thu, 28 Nov 2013 05:07:59 -0500
parents
children a6a56440567c
comparison
equal deleted inserted replaced
8:2ab1d6f6a8ec 9:302d13490b23
1 #!/usr/bin/env python
2 """A simple wrapper script to call MIRA4's mirabait and collect its output.
3 """
4 import os
5 import sys
6 import subprocess
7 import shutil
8 import time
9
10 WRAPPER_VER = "0.0.1" #Keep in sync with the XML file
11
12 def stop_err(msg, err=1):
13 sys.stderr.write(msg+"\n")
14 sys.exit(err)
15
16
17 def get_version(mira_binary):
18 """Run MIRA to find its version number"""
19 # At the commend line I would use: mira -v | head -n 1
20 # however there is some pipe error when doing that here.
21 cmd = [mira_binary, "-v"]
22 try:
23 child = subprocess.Popen(cmd,
24 stdout=subprocess.PIPE,
25 stderr=subprocess.STDOUT)
26 except Exception, err:
27 sys.stderr.write("Error invoking command:\n%s\n\n%s\n" % (" ".join(cmd), err))
28 sys.exit(1)
29 ver, tmp = child.communicate()
30 del child
31 #Workaround for -v not working in mirabait 4.0RC4
32 if "invalid option" in ver.split("\n", 1)[0]:
33 for line in ver.split("\n", 1):
34 if " version " in line:
35 line = line.split()
36 return line[line.index("version")+1].rstrip(")")
37 stop_err("Could not determine MIRA version:\n%s" % ver)
38 return ver.split("\n", 1)[0]
39
40 try:
41 mira_path = os.environ["MIRA4"]
42 except ImportError:
43 stop_err("Environment variable $MIRA4 not set")
44 mira_binary = os.path.join(mira_path, "mirabait")
45 if not os.path.isfile(mira_binary):
46 stop_err("Missing mirabait under $MIRA4, %r" % mira_binary)
47 mira_ver = get_version(mira_binary)
48 if not mira_ver.strip().startswith("4.0"):
49 stop_err("This wrapper is for MIRA V4.0, not:\n%s" % mira_ver)
50 if "-v" in sys.argv or "--version" in sys.argv:
51 print "%s, MIRA wrapper version %s" % (mira_ver, WRAPPER_VER)
52 sys.exit(0)
53
54
55 format, output_choice, strand_choice, kmer_length, min_occurance, bait_file, in_file, out_file = sys.argv[1:]
56
57 if format.startswith("fastq"):
58 format = "fastq"
59 elif format == "mira":
60 format = "maf"
61 elif format != "fasta":
62 stop_err("Was not expected format %r" % format)
63
64 assert out_file.endswith(".dat")
65 out_file_stem = out_file[:-4]
66
67 cmd_list = [mira_binary, "-f", format, "-t", format,
68 "-k", kmer_length, "-n", min_occurance,
69 bait_file, in_file, out_file_stem]
70 if output_choice == "pos":
71 pass
72 elif output_choice == "neg":
73 #Invert the selection...
74 cmd_list.insert(1, "-i")
75 else:
76 stop_err("Output choice should be 'pos' or 'neg', not %r" % output_choice)
77 if strand_choice == "both":
78 pass
79 elif strand_choice == "fwd":
80 #Ingore reverse strand...
81 cmd_list.insert(1, "-r")
82 else:
83 stop_err("Strand choice should be 'both' or 'fwd', not %r" % strand_choice)
84
85 cmd = " ".join(cmd_list)
86 #print cmd
87 start_time = time.time()
88 try:
89 #Run MIRA
90 child = subprocess.Popen(cmd_list,
91 stdout=subprocess.PIPE,
92 stderr=subprocess.STDOUT)
93 except Exception, err:
94 log_manifest(manifest)
95 sys.stderr.write("Error invoking command:\n%s\n\n%s\n" % (cmd, err))
96 sys.exit(1)
97 #Use .communicate as can get deadlocks with .wait(),
98 stdout, stderr = child.communicate()
99 assert stderr is None # Due to way we ran with subprocess
100 run_time = time.time() - start_time
101 return_code = child.returncode
102 print "mirabait took %0.2f minutes" % (run_time / 60.0)
103
104 if return_code:
105 sys.stderr.write(stdout)
106 stop_err("Return error code %i from command:\n%s" % (return_code, cmd),
107 return_code)
108
109 #Capture output
110 out_tmp = out_file_stem + "." + format
111 if not os.path.isfile(out_tmp):
112 sys.stderr.write(stdout)
113 stop_err("Missing output file from mirabait: %s" % out_tmp)
114 shutil.move(out_tmp, out_file)