annotate gene_family_classifier.py @ 109:06f121a1974e draft

Uploaded
author greg
date Mon, 27 Feb 2017 14:34:41 -0500
parents 82144bc3626a
children 6cb3520313b9
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
96
ba924b724a8d Uploaded
greg
parents:
diff changeset
1 #!/usr/bin/env python
ba924b724a8d Uploaded
greg
parents:
diff changeset
2 import argparse
ba924b724a8d Uploaded
greg
parents:
diff changeset
3 import os
ba924b724a8d Uploaded
greg
parents:
diff changeset
4 import shutil
ba924b724a8d Uploaded
greg
parents:
diff changeset
5 import subprocess
ba924b724a8d Uploaded
greg
parents:
diff changeset
6 import sys
ba924b724a8d Uploaded
greg
parents:
diff changeset
7 import tempfile
ba924b724a8d Uploaded
greg
parents:
diff changeset
8
ba924b724a8d Uploaded
greg
parents:
diff changeset
9 BUFF_SIZE = 1048576
ba924b724a8d Uploaded
greg
parents:
diff changeset
10 OUTPUT_DIR = 'geneFamilyClassification_dir'
ba924b724a8d Uploaded
greg
parents:
diff changeset
11
ba924b724a8d Uploaded
greg
parents:
diff changeset
12 parser = argparse.ArgumentParser()
ba924b724a8d Uploaded
greg
parents:
diff changeset
13 parser.add_argument('--input', dest='input', help="Input dataset")
ba924b724a8d Uploaded
greg
parents:
diff changeset
14 parser.add_argument('--scaffold', dest='scaffold', default='mode', help='Orthogroups or gene families proteins scaffold')
ba924b724a8d Uploaded
greg
parents:
diff changeset
15 parser.add_argument('--method', dest='method', help='Protein clustering method')
ba924b724a8d Uploaded
greg
parents:
diff changeset
16 parser.add_argument('--classifier', dest='classifier', help='Protein classification method')
ba924b724a8d Uploaded
greg
parents:
diff changeset
17 parser.add_argument('--config_dir', dest='config_dir', help='Directory containing default configuration files')
ba924b724a8d Uploaded
greg
parents:
diff changeset
18 parser.add_argument('--num_threads', dest='num_threads', type=int, help='Number of threads to use for execution')
ba924b724a8d Uploaded
greg
parents:
diff changeset
19 parser.add_argument('--super_orthogroups', dest='super_orthogroups', default=None, help='Super orthogroups clustering specification')
ba924b724a8d Uploaded
greg
parents:
diff changeset
20 parser.add_argument('--single_copy_custom', dest='single_copy_custom', default=None, help='Custom single copy orthogroup configuration')
ba924b724a8d Uploaded
greg
parents:
diff changeset
21 parser.add_argument('--single_copy_taxa', dest='single_copy_taxa', type=int, default=None, help='"Minimum single copy taxa required in orthogroup')
ba924b724a8d Uploaded
greg
parents:
diff changeset
22 parser.add_argument('--taxa_present', dest='taxa_present', type=int, default=None, help='Minimum taxa required in single copy orthogroup')
ba924b724a8d Uploaded
greg
parents:
diff changeset
23 parser.add_argument('--orthogroup_fasta', dest='orthogroup_fasta', default=None, help='Flag to create orthogroup sequences')
ba924b724a8d Uploaded
greg
parents:
diff changeset
24 parser.add_argument('--coding_sequences', dest='coding_sequences', default=None, help='Flag to create orthogroup coding sequences')
ba924b724a8d Uploaded
greg
parents:
diff changeset
25 parser.add_argument('--save_hmmscan_log', dest='save_hmmscan_log', default=None, help='Flag to save the hmmscan log')
ba924b724a8d Uploaded
greg
parents:
diff changeset
26 parser.add_argument('--hmmscan_log', dest='hmmscan_log', default=None, help='hmmscan log file')
ba924b724a8d Uploaded
greg
parents:
diff changeset
27 parser.add_argument('--output_ptortho', dest='output_ptortho', default=None, help='Output for orthogroups')
101
cd23ab8b617f Uploaded
greg
parents: 96
diff changeset
28 parser.add_argument('--output_ptortho_dir', dest='output_ptortho_dir', default=None, help='output_ptortho.files_path')
96
ba924b724a8d Uploaded
greg
parents:
diff changeset
29 parser.add_argument('--output_ptorthocs', dest='output_ptorthocs', default=None, help='Output for orthogroups with corresponding coding sequences')
101
cd23ab8b617f Uploaded
greg
parents: 96
diff changeset
30 parser.add_argument('--output_ptorthocs_dir', dest='output_ptorthocs_dir', default=None, help='output_ptorthocs.files_path')
96
ba924b724a8d Uploaded
greg
parents:
diff changeset
31 parser.add_argument('--output_ptsco', dest='output_ptsco', default=None, help='Output for single copy orthogroups')
101
cd23ab8b617f Uploaded
greg
parents: 96
diff changeset
32 parser.add_argument('--output_ptsco_dir', dest='output_ptsco_dir', default=None, help='output_ptsco.files_path')
96
ba924b724a8d Uploaded
greg
parents:
diff changeset
33
ba924b724a8d Uploaded
greg
parents:
diff changeset
34 args = parser.parse_args()
ba924b724a8d Uploaded
greg
parents:
diff changeset
35
ba924b724a8d Uploaded
greg
parents:
diff changeset
36
ba924b724a8d Uploaded
greg
parents:
diff changeset
37 def get_stderr_exception(tmp_err, tmp_stderr, tmp_out, tmp_stdout, include_stdout=False):
ba924b724a8d Uploaded
greg
parents:
diff changeset
38 tmp_stderr.close()
ba924b724a8d Uploaded
greg
parents:
diff changeset
39 """
ba924b724a8d Uploaded
greg
parents:
diff changeset
40 Return a stderr string of reasonable size.
ba924b724a8d Uploaded
greg
parents:
diff changeset
41 """
ba924b724a8d Uploaded
greg
parents:
diff changeset
42 # Get stderr, allowing for case where it's very large.
ba924b724a8d Uploaded
greg
parents:
diff changeset
43 tmp_stderr = open(tmp_err, 'rb')
ba924b724a8d Uploaded
greg
parents:
diff changeset
44 stderr_str = ''
ba924b724a8d Uploaded
greg
parents:
diff changeset
45 buffsize = BUFF_SIZE
ba924b724a8d Uploaded
greg
parents:
diff changeset
46 try:
ba924b724a8d Uploaded
greg
parents:
diff changeset
47 while True:
ba924b724a8d Uploaded
greg
parents:
diff changeset
48 stderr_str += tmp_stderr.read(buffsize)
ba924b724a8d Uploaded
greg
parents:
diff changeset
49 if not stderr_str or len(stderr_str) % buffsize != 0:
ba924b724a8d Uploaded
greg
parents:
diff changeset
50 break
ba924b724a8d Uploaded
greg
parents:
diff changeset
51 except OverflowError:
ba924b724a8d Uploaded
greg
parents:
diff changeset
52 pass
ba924b724a8d Uploaded
greg
parents:
diff changeset
53 tmp_stderr.close()
ba924b724a8d Uploaded
greg
parents:
diff changeset
54 if include_stdout:
ba924b724a8d Uploaded
greg
parents:
diff changeset
55 tmp_stdout = open(tmp_out, 'rb')
ba924b724a8d Uploaded
greg
parents:
diff changeset
56 stdout_str = ''
ba924b724a8d Uploaded
greg
parents:
diff changeset
57 buffsize = BUFF_SIZE
ba924b724a8d Uploaded
greg
parents:
diff changeset
58 try:
ba924b724a8d Uploaded
greg
parents:
diff changeset
59 while True:
ba924b724a8d Uploaded
greg
parents:
diff changeset
60 stdout_str += tmp_stdout.read(buffsize)
ba924b724a8d Uploaded
greg
parents:
diff changeset
61 if not stdout_str or len(stdout_str) % buffsize != 0:
ba924b724a8d Uploaded
greg
parents:
diff changeset
62 break
ba924b724a8d Uploaded
greg
parents:
diff changeset
63 except OverflowError:
ba924b724a8d Uploaded
greg
parents:
diff changeset
64 pass
ba924b724a8d Uploaded
greg
parents:
diff changeset
65 tmp_stdout.close()
ba924b724a8d Uploaded
greg
parents:
diff changeset
66 if include_stdout:
ba924b724a8d Uploaded
greg
parents:
diff changeset
67 return 'STDOUT\n%s\n\nSTDERR\n%s\n' % (stdout_str, stderr_str)
ba924b724a8d Uploaded
greg
parents:
diff changeset
68 return stderr_str
ba924b724a8d Uploaded
greg
parents:
diff changeset
69
ba924b724a8d Uploaded
greg
parents:
diff changeset
70
ba924b724a8d Uploaded
greg
parents:
diff changeset
71 def move_directory_files(source_dir, destination_dir):
ba924b724a8d Uploaded
greg
parents:
diff changeset
72 source_directory = os.path.abspath(source_dir)
ba924b724a8d Uploaded
greg
parents:
diff changeset
73 destination_directory = os.path.abspath(destination_dir)
ba924b724a8d Uploaded
greg
parents:
diff changeset
74 if not os.path.isdir(destination_directory):
ba924b724a8d Uploaded
greg
parents:
diff changeset
75 os.makedirs(destination_directory)
ba924b724a8d Uploaded
greg
parents:
diff changeset
76 for dir_entry in os.listdir(source_directory):
ba924b724a8d Uploaded
greg
parents:
diff changeset
77 source_entry = os.path.join(source_directory, dir_entry)
ba924b724a8d Uploaded
greg
parents:
diff changeset
78 shutil.move(source_entry, destination_directory)
ba924b724a8d Uploaded
greg
parents:
diff changeset
79
ba924b724a8d Uploaded
greg
parents:
diff changeset
80
ba924b724a8d Uploaded
greg
parents:
diff changeset
81 def stop_err(msg):
ba924b724a8d Uploaded
greg
parents:
diff changeset
82 sys.stderr.write(msg)
ba924b724a8d Uploaded
greg
parents:
diff changeset
83 sys.exit(1)
ba924b724a8d Uploaded
greg
parents:
diff changeset
84
ba924b724a8d Uploaded
greg
parents:
diff changeset
85
ba924b724a8d Uploaded
greg
parents:
diff changeset
86 def write_html_output(output, title, dir):
107
bdb3929a4d09 Uploaded
greg
parents: 106
diff changeset
87 with open(output, 'w') as fh:
109
06f121a1974e Uploaded
greg
parents: 108
diff changeset
88 fh.write('<html><head><title>%s</title></head>\n' % title)
06f121a1974e Uploaded
greg
parents: 108
diff changeset
89 fh.write('<body><p/><table>\n')
06f121a1974e Uploaded
greg
parents: 108
diff changeset
90 fh.write('<tr><th>Size</th><th>Name</th></tr>\n')
106
a2a3593a42d4 Uploaded
greg
parents: 105
diff changeset
91 for fname in os.listdir(dir):
109
06f121a1974e Uploaded
greg
parents: 108
diff changeset
92 size = str(os.path.getsixe(fname))
06f121a1974e Uploaded
greg
parents: 108
diff changeset
93 link = '<a href="%s" type="text/plain">%s</a>\n' % (fname, fname)
06f121a1974e Uploaded
greg
parents: 108
diff changeset
94 fh.write('<tr><td>%s</td><td>%s</td></tr>\n' % (size, link))
06f121a1974e Uploaded
greg
parents: 108
diff changeset
95 fh.write('</table></body></html>\n')
96
ba924b724a8d Uploaded
greg
parents:
diff changeset
96
ba924b724a8d Uploaded
greg
parents:
diff changeset
97
ba924b724a8d Uploaded
greg
parents:
diff changeset
98 # Define command response buffers.
ba924b724a8d Uploaded
greg
parents:
diff changeset
99 tmp_out = tempfile.NamedTemporaryFile().name
ba924b724a8d Uploaded
greg
parents:
diff changeset
100 tmp_stdout = open(tmp_out, 'wb')
ba924b724a8d Uploaded
greg
parents:
diff changeset
101 tmp_err = tempfile.NamedTemporaryFile().name
ba924b724a8d Uploaded
greg
parents:
diff changeset
102 tmp_stderr = open(tmp_err, 'wb')
ba924b724a8d Uploaded
greg
parents:
diff changeset
103 # Build the command line.
ba924b724a8d Uploaded
greg
parents:
diff changeset
104 cmd = 'GeneFamilyClassifier'
ba924b724a8d Uploaded
greg
parents:
diff changeset
105 cmd += ' --proteins %s' % args.input
ba924b724a8d Uploaded
greg
parents:
diff changeset
106 cmd += ' --scaffold %s' % args.scaffold
ba924b724a8d Uploaded
greg
parents:
diff changeset
107 cmd += ' --method %s' % args.method
ba924b724a8d Uploaded
greg
parents:
diff changeset
108 cmd += ' --classifier %s' % args.classifier
ba924b724a8d Uploaded
greg
parents:
diff changeset
109 cmd += ' --config_dir %s' % args.config_dir
ba924b724a8d Uploaded
greg
parents:
diff changeset
110 cmd += ' --num_threads %d' % args.num_threads
ba924b724a8d Uploaded
greg
parents:
diff changeset
111 if args.super_orthogroups is not None:
ba924b724a8d Uploaded
greg
parents:
diff changeset
112 cmd += ' --super_orthogroups %s' % args.super_orthogroups
ba924b724a8d Uploaded
greg
parents:
diff changeset
113 if args.single_copy_custom is not None:
ba924b724a8d Uploaded
greg
parents:
diff changeset
114 cmd += ' --single_copy_custom %s' % args.single_copy_custom
ba924b724a8d Uploaded
greg
parents:
diff changeset
115 if args.single_copy_taxa is not None:
ba924b724a8d Uploaded
greg
parents:
diff changeset
116 cmd += ' --single_copy_taxa %d' % args.single_copy_taxa
ba924b724a8d Uploaded
greg
parents:
diff changeset
117 if args.taxa_present is not None:
ba924b724a8d Uploaded
greg
parents:
diff changeset
118 cmd += ' --taxa_present %d' % args.taxa_present
ba924b724a8d Uploaded
greg
parents:
diff changeset
119 if args.orthogroup_fasta is None:
ba924b724a8d Uploaded
greg
parents:
diff changeset
120 create_ortho_sequences = False
ba924b724a8d Uploaded
greg
parents:
diff changeset
121 else:
ba924b724a8d Uploaded
greg
parents:
diff changeset
122 create_ortho_sequences = True
ba924b724a8d Uploaded
greg
parents:
diff changeset
123 cmd += ' --orthogroup_fasta'
ba924b724a8d Uploaded
greg
parents:
diff changeset
124 if args.coding_sequences is None:
ba924b724a8d Uploaded
greg
parents:
diff changeset
125 create_corresponding_coding_sequences = False
ba924b724a8d Uploaded
greg
parents:
diff changeset
126 else:
ba924b724a8d Uploaded
greg
parents:
diff changeset
127 create_corresponding_coding_sequences = True
ba924b724a8d Uploaded
greg
parents:
diff changeset
128 cmd += ' --coding_sequences %s' % args.coding_sequences
ba924b724a8d Uploaded
greg
parents:
diff changeset
129 # Run the command.
ba924b724a8d Uploaded
greg
parents:
diff changeset
130 proc = subprocess.Popen(args=cmd, stderr=tmp_stderr, stdout=tmp_stdout, shell=True)
ba924b724a8d Uploaded
greg
parents:
diff changeset
131 rc = proc.wait()
ba924b724a8d Uploaded
greg
parents:
diff changeset
132 # Handle execution errors.
ba924b724a8d Uploaded
greg
parents:
diff changeset
133 if rc != 0:
ba924b724a8d Uploaded
greg
parents:
diff changeset
134 error_message = get_stderr_exception(tmp_err, tmp_stderr, tmp_out, tmp_stdout)
ba924b724a8d Uploaded
greg
parents:
diff changeset
135 stop_err( error_message )
ba924b724a8d Uploaded
greg
parents:
diff changeset
136 # Handle hmmscan.log output.
ba924b724a8d Uploaded
greg
parents:
diff changeset
137 if args.classifier in ['hmmscan', 'both']:
ba924b724a8d Uploaded
greg
parents:
diff changeset
138 src_hmmscan_log = os.path.join(OUTPUT_DIR, 'hmmscan.log')
ba924b724a8d Uploaded
greg
parents:
diff changeset
139 if os.path.exists(src_hmmscan_log):
ba924b724a8d Uploaded
greg
parents:
diff changeset
140 if args.save_hmmscan_log is None:
ba924b724a8d Uploaded
greg
parents:
diff changeset
141 os.remove(args.hmmscan_log)
ba924b724a8d Uploaded
greg
parents:
diff changeset
142 else:
ba924b724a8d Uploaded
greg
parents:
diff changeset
143 shutil.move(src_hmmscan_log, args.hmmscan_log)
ba924b724a8d Uploaded
greg
parents:
diff changeset
144 # Handle orthogroups outputs.
ba924b724a8d Uploaded
greg
parents:
diff changeset
145 if create_ortho_sequences:
ba924b724a8d Uploaded
greg
parents:
diff changeset
146 if create_corresponding_coding_sequences:
ba924b724a8d Uploaded
greg
parents:
diff changeset
147 out_file = args.output_ptorthocs
101
cd23ab8b617f Uploaded
greg
parents: 96
diff changeset
148 orthogroups_fasta_dest_dir = args.output_ptorthocs_dir
109
06f121a1974e Uploaded
greg
parents: 108
diff changeset
149 title = 'Orthogroups and corresponding coding sequences files'
96
ba924b724a8d Uploaded
greg
parents:
diff changeset
150 else:
ba924b724a8d Uploaded
greg
parents:
diff changeset
151 out_file = args.output_ptortho
101
cd23ab8b617f Uploaded
greg
parents: 96
diff changeset
152 orthogroups_fasta_dest_dir = args.output_ptortho_dir
109
06f121a1974e Uploaded
greg
parents: 108
diff changeset
153 title = 'Orthogroups files'
96
ba924b724a8d Uploaded
greg
parents:
diff changeset
154 orthogroups_fasta_src_dir = os.path.join(OUTPUT_DIR, 'orthogroups_fasta')
103
55052911cac3 Uploaded
greg
parents: 101
diff changeset
155 move_directory_files(orthogroups_fasta_src_dir, orthogroups_fasta_dest_dir)
96
ba924b724a8d Uploaded
greg
parents:
diff changeset
156 write_html_output(out_file, title, orthogroups_fasta_dest_dir)
ba924b724a8d Uploaded
greg
parents:
diff changeset
157 # Handle single copy orthogroup outputs.
ba924b724a8d Uploaded
greg
parents:
diff changeset
158 if args.output_ptsco is not None:
ba924b724a8d Uploaded
greg
parents:
diff changeset
159 single_copy_fasta_src_dir = os.path.join(OUTPUT_DIR, 'single_copy_fasta')
101
cd23ab8b617f Uploaded
greg
parents: 96
diff changeset
160 single_copy_fasta_dest_dir = args.output_ptsco_dir
109
06f121a1974e Uploaded
greg
parents: 108
diff changeset
161 title = 'Single copy orthogroups files'
96
ba924b724a8d Uploaded
greg
parents:
diff changeset
162 move_directory_files(single_copy_fasta_src_dir, single_copy_fasta_dest_dir)
ba924b724a8d Uploaded
greg
parents:
diff changeset
163 write_html_output(args.output_ptsco, title, single_copy_fasta_dest_dir)