|
25
|
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.5" #Keep in sync with the XML file
|
|
|
11
|
|
|
12 def sys_exit(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 sys_exit("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 KeyError:
|
|
|
43 sys_exit("Environment variable $MIRA4 not set")
|
|
|
44 mira_binary = os.path.join(mira_path, "mirabait")
|
|
|
45 if not os.path.isfile(mira_binary):
|
|
|
46 sys_exit("Missing mirabait under $MIRA4, %r\nFolder contained: %s"
|
|
|
47 % (mira_binary, ", ".join(os.listdir(mira_path))))
|
|
|
48 mira_ver = get_version(mira_binary)
|
|
|
49 if not mira_ver.strip().startswith("4.0"):
|
|
|
50 sys_exit("This wrapper is for MIRA V4.0, not:\n%s" % mira_ver)
|
|
|
51 if "-v" in sys.argv or "--version" in sys.argv:
|
|
|
52 print "%s, MIRA wrapper version %s" % (mira_ver, WRAPPER_VER)
|
|
|
53 sys.exit(0)
|
|
|
54
|
|
|
55
|
|
|
56 format, output_choice, strand_choice, kmer_length, min_occurance, bait_file, in_file, out_file = sys.argv[1:]
|
|
|
57
|
|
|
58 if format.startswith("fastq"):
|
|
|
59 format = "fastq"
|
|
|
60 elif format == "mira":
|
|
|
61 format = "maf"
|
|
|
62 elif format != "fasta":
|
|
|
63 sys_exit("Was not expected format %r" % format)
|
|
|
64
|
|
|
65 assert out_file.endswith(".dat")
|
|
|
66 out_file_stem = out_file[:-4]
|
|
|
67
|
|
|
68 cmd_list = [mira_binary, "-f", format, "-t", format,
|
|
|
69 "-k", kmer_length, "-n", min_occurance,
|
|
|
70 bait_file, in_file, out_file_stem]
|
|
|
71 if output_choice == "pos":
|
|
|
72 pass
|
|
|
73 elif output_choice == "neg":
|
|
|
74 #Invert the selection...
|
|
|
75 cmd_list.insert(1, "-i")
|
|
|
76 else:
|
|
|
77 sys_exit("Output choice should be 'pos' or 'neg', not %r" % output_choice)
|
|
|
78 if strand_choice == "both":
|
|
|
79 pass
|
|
|
80 elif strand_choice == "fwd":
|
|
|
81 #Ingore reverse strand...
|
|
|
82 cmd_list.insert(1, "-r")
|
|
|
83 else:
|
|
|
84 sys_exit("Strand choice should be 'both' or 'fwd', not %r" % strand_choice)
|
|
|
85
|
|
|
86 cmd = " ".join(cmd_list)
|
|
|
87 #print cmd
|
|
|
88 start_time = time.time()
|
|
|
89 try:
|
|
|
90 #Run MIRA
|
|
|
91 child = subprocess.Popen(cmd_list,
|
|
|
92 stdout=subprocess.PIPE,
|
|
|
93 stderr=subprocess.STDOUT)
|
|
|
94 except Exception, err:
|
|
|
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 sys_exit("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 sys_exit("Missing output file from mirabait: %s" % out_tmp)
|
|
|
114 shutil.move(out_tmp, out_file)
|
