Mercurial > repos > peterjc > effectivet3

comparison tools/protein_analysis/effectiveT3.py @ 0:43436379876f

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Migrated tool version 0.0.7 from old tool shed archive to new tool shed repository
author peterjc
date Tue, 07 Jun 2011 17:21:30 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:43436379876f
1 #!/usr/bin/env python
2 """Wrapper for EffectiveT3 v1.0.1 for use in Galaxy.
3
4 This script takes exactly five command line arguments:
5 * model name (e.g. TTSS_STD-1.0.1.jar)
6 * threshold (selective or sensitive)
7 * an input protein FASTA filename
8 * output tabular filename
9
10 It then calls the standalone Effective T3 v1.0.1 program (not the
11 webservice), and reformats the semi-colon separated output into
12 tab separated output for use in Galaxy.
13 """
14 import sys
15 import os
16 import subprocess
17
18 #You may need to edit this to match your local setup,
19 effectiveT3_jar = "/opt/EffectiveT3/TTSS_GUI-1.0.1.jar"
20
21
22 def stop_err(msg, error_level=1):
23 """Print error message to stdout and quit with given error level."""
24 sys.stderr.write("%s\n" % msg)
25 sys.exit(error_level)
26
27 if len(sys.argv) != 5:
28 stop_err("Require four arguments: model, threshold, input protein FASTA file & output tabular file")
29
30 model, threshold, fasta_file, tabular_file = sys.argv[1:]
31
32 if not os.path.isfile(fasta_file):
33 stop_err("Input FASTA file not found: %s" % fasta_file)
34
35 if threshold not in ["selective", "sensitive"] \
36 and not threshold.startswith("cutoff="):
37 stop_err("Threshold should be selective, sensitive, or cutoff=..., not %r" % threshold)
38
39 def clean_tabular(raw_handle, out_handle):
40 """Clean up Effective T3 output to make it tabular."""
41 count = 0
42 positive = 0
43 errors = 0
44 for line in raw_handle:
45 if not line or line.startswith("#") \
46 or line.startswith("Id; Description; Score;"):
47 continue
48 assert line.count(";") >= 3, repr(line)
49 #Normally there will just be three semi-colons, however the
50 #original FASTA file's ID or description might have had
51 #semi-colons in it as well, hence the following hackery:
52 try:
53 id_descr, score, effective = line.rstrip("\r\n").rsplit(";",2)
54 #Cope when there was no FASTA description
55 if "; " not in id_descr and id_descr.endswith(";"):
56 id = id_descr[:-1]
57 descr = ""
58 else:
59 id, descr = id_descr.split("; ",1)
60 except ValueError:
61 stop_err("Problem parsing line:\n%s\n" % line)
62 parts = [s.strip() for s in [id, descr, score, effective]]
63 out_handle.write("\t".join(parts) + "\n")
64 count += 1
65 if float(score) < 0:
66 errors += 1
67 if effective.lower() == "true":
68 positive += 1
69 return count, positive, errors
70
71 def run(cmd):
72 #Avoid using shell=True when we call subprocess to ensure if the Python
73 #script is killed, so too is the child process.
74 try:
75 child = subprocess.Popen(cmd, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
76 except Exception, err:
77 stop_err("Error invoking command:\n%s\n\n%s\n" % (" ".join(cmd), err))
78 #Use .communicate as can get deadlocks with .wait(),
79 stdout, stderr = child.communicate()
80 return_code = child.returncode
81 if return_code:
82 if stderr and stdout:
83 stop_err("Return code %i from command:\n%s\n\n%s\n\n%s" % (return_code, err, stdout, stderr))
84 else:
85 stop_err("Return code %i from command:\n%s\n%s" % (return_code, err, stderr))
86
87 if not os.path.isfile(effectiveT3_jar):
88 stop_err("Effective T3 JAR file not found: %s" % effectiveT3_jar)
89
90 effectiveT3_dir = os.path.dirname(effectiveT3_jar)
91 if not os.path.isdir(effectiveT3_dir):
92 stop_err("Effective T3 folder not found: %s" % effectiveT3_dir)
93
94 effectiveT3_model = os.path.join(effectiveT3_dir, "module", model)
95 if not os.path.isfile(effectiveT3_model):
96 stop_err("Effective T3 model JAR file not found: %s" % effectiveT3_model)
97
98 #We will have write access whereever the output should be,
99 temp_file = os.path.abspath(tabular_file + ".tmp")
100
101 #Use absolute paths since will change current directory...
102 tabular_file = os.path.abspath(tabular_file)
103 fasta_file = os.path.abspath(fasta_file)
104
105 cmd = ["java", "-jar", effectiveT3_jar,
106 "-f", fasta_file,
107 "-m", model,
108 "-t", threshold,
109 "-o", temp_file,
110 "-q"]
111
112 try:
113 #Must run from directory above the module subfolder:
114 os.chdir(effectiveT3_dir)
115 except:
116 stop_err("Could not change to Effective T3 folder: %s" % effectiveT3_dir)
117
118 run(cmd)
119
120 if not os.path.isfile(temp_file):
121 stop_err("ERROR - No output file from Effective T3")
122
123 out_handle = open(tabular_file, "w")
124 out_handle.write("#ID\tDescription\tScore\tEffective\n")
125 data_handle = open(temp_file)
126 count, positive, errors = clean_tabular(data_handle, out_handle)
127 data_handle.close()
128 out_handle.close()
129
130 os.remove(temp_file)
131
132 if errors:
133 print "%i sequences, %i positive, %i errors" \
134 % (count, positive, errors)
135 else:
136 print "%i/%i sequences positive" % (positive, count)
137
138 if count and count==errors:
139 #Galaxy will still allow them to see the output file
140 stop_err("All your sequences gave an error code")