Mercurial > repos > greg > gene_family_classifier

comparison gene_family_classifier.py @ 96:ba924b724a8d draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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)