|
7
|
1 #!/usr/bin/env python
|
|
|
2 """Wrapper for Promoter 2.0 for use in Galaxy.
|
|
|
3
|
|
|
4 This script takes exactly three command line arguments:
|
|
|
5 * number of threads
|
|
|
6 * an input DNA FASTA filename
|
|
|
7 * output tabular filename.
|
|
|
8
|
|
|
9 It calls the Promoter 2.0 binary (e.g. .../promoter-2.0/bin/promoter_Linux,
|
|
|
10 bypassing the Perl wrapper script 'promoter' which imposes a significant
|
|
|
11 performace overhead for no benefit here (we don't need HTML output for
|
|
|
12 example).
|
|
|
13
|
|
|
14 The main feature is this Python wrapper script parsers the bespoke
|
|
|
15 tabular output from Promoter 2.0 and reformats it into a Galaxy friendly
|
|
|
16 tab separated table.
|
|
|
17
|
|
|
18 Additionally, in order to take advantage of multiple cores the input FASTA
|
|
|
19 file is broken into chunks and multiple copies of promoter run at once.
|
|
|
20 This can be used in combination with the job-splitting available in Galaxy.
|
|
|
21
|
|
|
22 Note that rewriting the FASTA input file allows us to avoid a bug in
|
|
|
23 promoter 2 with long descriptions in the FASTA header line (over 200
|
|
|
24 characters) which produces stray fragements of the description in the
|
|
|
25 output file, making parsing non-trivial.
|
|
|
26
|
|
|
27 TODO - Automatically extract the sequence containing a promoter prediction?
|
|
|
28 """
|
|
|
29 import sys
|
|
|
30 import os
|
|
|
31 import commands
|
|
|
32 import tempfile
|
|
29
|
33 from seq_analysis_utils import split_fasta, run_jobs, thread_count
|
|
7
|
34
|
|
|
35 FASTA_CHUNK = 500
|
|
|
36
|
|
28
|
37 if "-v" in sys.argv or "--version" in sys.argv:
|
|
|
38 sys.exit(os.system("promoter -V"))
|
|
|
39
|
|
7
|
40 if len(sys.argv) != 4:
|
|
29
|
41 sys.exit("Require three arguments, number of threads (int), input DNA FASTA file & output tabular file. "
|
|
|
42 "Got %i arguments." % (len(sys.argv) - 1))
|
|
7
|
43
|
|
29
|
44 num_threads = thread_count(sys.argv[3], default=4)
|
|
7
|
45 fasta_file = os.path.abspath(sys.argv[2])
|
|
|
46 tabular_file = os.path.abspath(sys.argv[3])
|
|
|
47
|
|
|
48 tmp_dir = tempfile.mkdtemp()
|
|
|
49
|
|
29
|
50
|
|
7
|
51 def get_path_and_binary():
|
|
29
|
52 platform = commands.getoutput("uname") # e.g. Linux
|
|
7
|
53 shell_script = commands.getoutput("which promoter")
|
|
|
54 if not os.path.isfile(shell_script):
|
|
29
|
55 sys.exit("ERROR: Missing promoter executable shell script")
|
|
7
|
56 path = None
|
|
|
57 for line in open(shell_script):
|
|
29
|
58 if line.startswith("setenv"): # could then be tab or space!
|
|
7
|
59 parts = line.rstrip().split(None, 2)
|
|
|
60 if parts[0] == "setenv" and parts[1] == "PROM":
|
|
|
61 path = parts[2]
|
|
|
62 if not path:
|
|
29
|
63 sys.exit("ERROR: Could not find promoter path (PROM) in %r" % shell_script)
|
|
7
|
64 if not os.path.isdir(path):
|
|
29
|
65 sys.exit("ERROR: %r is not a directory" % path)
|
|
7
|
66 bin = "%s/bin/promoter_%s" % (path, platform)
|
|
|
67 if not os.path.isfile(bin):
|
|
29
|
68 sys.exit("ERROR: Missing promoter binary %r" % bin)
|
|
7
|
69 return path, bin
|
|
|
70
|
|
29
|
71
|
|
7
|
72 def make_tabular(raw_handle, out_handle):
|
|
|
73 """Parse text output into tabular, return query count."""
|
|
|
74 identifier = None
|
|
|
75 queries = 0
|
|
|
76 for line in raw_handle:
|
|
29
|
77 # print repr(line)
|
|
7
|
78 if not line.strip() or line == "Promoter prediction:\n":
|
|
|
79 pass
|
|
|
80 elif line[0] != " ":
|
|
29
|
81 identifier = line.strip().replace("\t", " ").split(None, 1)[0]
|
|
7
|
82 queries += 1
|
|
|
83 elif line == " No promoter predicted\n":
|
|
29
|
84 # End of a record
|
|
7
|
85 identifier = None
|
|
|
86 elif line == " Position Score Likelihood\n":
|
|
|
87 assert identifier
|
|
|
88 else:
|
|
|
89 try:
|
|
29
|
90 position, score, likelihood = line.strip().split(None, 2)
|
|
7
|
91 except ValueError:
|
|
|
92 print "WARNING: Problem with line: %r" % line
|
|
|
93 continue
|
|
29
|
94 # sys.exit("ERROR: Problem with line: %r" % line)
|
|
7
|
95 if likelihood not in ["ignored",
|
|
|
96 "Marginal prediction",
|
|
|
97 "Medium likely prediction",
|
|
|
98 "Highly likely prediction"]:
|
|
29
|
99 sys.exit("ERROR: Problem with line: %r" % line)
|
|
7
|
100 out_handle.write("%s\t%s\t%s\t%s\n" % (identifier, position, score, likelihood))
|
|
|
101 return queries
|
|
29
|
102
|
|
7
|
103 working_dir, bin = get_path_and_binary()
|
|
|
104
|
|
|
105 if not os.path.isfile(fasta_file):
|
|
29
|
106 sys.exit("ERROR: Missing input FASTA file %r" % fasta_file)
|
|
7
|
107
|
|
29
|
108 # Note that if the input FASTA file contains no sequences,
|
|
|
109 # split_fasta returns an empty list (i.e. zero temp files).
|
|
|
110 # We deliberately omit the FASTA descriptions to avoid a
|
|
|
111 # bug in promoter2 with descriptions over 200 characters.
|
|
7
|
112 fasta_files = split_fasta(fasta_file, os.path.join(tmp_dir, "promoter"), FASTA_CHUNK, keep_descr=False)
|
|
29
|
113 temp_files = [f + ".out" for f in fasta_files]
|
|
7
|
114 jobs = ["%s %s > %s" % (bin, fasta, temp)
|
|
|
115 for fasta, temp in zip(fasta_files, temp_files)]
|
|
|
116
|
|
29
|
117
|
|
7
|
118 def clean_up(file_list):
|
|
|
119 for f in file_list:
|
|
|
120 if os.path.isfile(f):
|
|
|
121 os.remove(f)
|
|
|
122 try:
|
|
|
123 os.rmdir(tmp_dir)
|
|
29
|
124 except Exception:
|
|
7
|
125 pass
|
|
|
126
|
|
|
127 if len(jobs) > 1 and num_threads > 1:
|
|
29
|
128 # A small "info" message for Galaxy to show the user.
|
|
7
|
129 print "Using %i threads for %i tasks" % (min(num_threads, len(jobs)), len(jobs))
|
|
|
130 cur_dir = os.path.abspath(os.curdir)
|
|
|
131 os.chdir(working_dir)
|
|
|
132 results = run_jobs(jobs, num_threads)
|
|
|
133 os.chdir(cur_dir)
|
|
|
134 for fasta, temp, cmd in zip(fasta_files, temp_files, jobs):
|
|
|
135 error_level = results[cmd]
|
|
|
136 if error_level:
|
|
|
137 try:
|
|
|
138 output = open(temp).readline()
|
|
|
139 except IOError:
|
|
|
140 output = ""
|
|
|
141 clean_up(fasta_files + temp_files)
|
|
29
|
142 sys.exit("One or more tasks failed, e.g. %i from %r gave:\n%s" % (error_level, cmd, output),
|
|
7
|
143 error_level)
|
|
|
144
|
|
|
145 del results
|
|
|
146 del jobs
|
|
|
147
|
|
|
148 out_handle = open(tabular_file, "w")
|
|
|
149 out_handle.write("#Identifier\tPosition\tScore\tLikelihood\n")
|
|
|
150 queries = 0
|
|
|
151 for temp in temp_files:
|
|
|
152 data_handle = open(temp)
|
|
|
153 count = make_tabular(data_handle, out_handle)
|
|
|
154 data_handle.close()
|
|
|
155 if not count:
|
|
|
156 clean_up(fasta_files + temp_files)
|
|
29
|
157 sys.exit("No output from promoter2")
|
|
7
|
158 queries += count
|
|
|
159 out_handle.close()
|
|
|
160
|
|
|
161 clean_up(fasta_files + temp_files)
|
|
|
162 print "Results for %i queries" % queries
|
