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