# HG changeset patch # User greg # Date 1488210480 18000 # Node ID ba924b724a8da80ddb00b998490c300f2a13db1a # Parent 431a63fced688ecfd7c7929defec6a61130614bb Uploaded diff -r 431a63fced68 -r ba924b724a8d gene_family_classifier.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/gene_family_classifier.py Mon Feb 27 10:48:00 2017 -0500 @@ -0,0 +1,156 @@ +#!/usr/bin/env python +import argparse +import os +import shutil +import subprocess +import sys +import tempfile + +BUFF_SIZE = 1048576 +OUTPUT_DIR = 'geneFamilyClassification_dir' + +parser = argparse.ArgumentParser() +parser.add_argument('--input', dest='input', help="Input dataset") +parser.add_argument('--scaffold', dest='scaffold', default='mode', help='Orthogroups or gene families proteins scaffold') +parser.add_argument('--method', dest='method', help='Protein clustering method') +parser.add_argument('--classifier', dest='classifier', help='Protein classification method') +parser.add_argument('--config_dir', dest='config_dir', help='Directory containing default configuration files') +parser.add_argument('--num_threads', dest='num_threads', type=int, help='Number of threads to use for execution') +parser.add_argument('--super_orthogroups', dest='super_orthogroups', default=None, help='Super orthogroups clustering specification') +parser.add_argument('--single_copy_custom', dest='single_copy_custom', default=None, help='Custom single copy orthogroup configuration') +parser.add_argument('--single_copy_taxa', dest='single_copy_taxa', type=int, default=None, help='"Minimum single copy taxa required in orthogroup') +parser.add_argument('--taxa_present', dest='taxa_present', type=int, default=None, help='Minimum taxa required in single copy orthogroup') +parser.add_argument('--orthogroup_fasta', dest='orthogroup_fasta', default=None, help='Flag to create orthogroup sequences') +parser.add_argument('--coding_sequences', dest='coding_sequences', default=None, help='Flag to create orthogroup coding sequences') +parser.add_argument('--save_hmmscan_log', dest='save_hmmscan_log', default=None, help='Flag to save the hmmscan log') +parser.add_argument('--hmmscan_log', dest='hmmscan_log', default=None, help='hmmscan log file') +parser.add_argument('--output_ptortho', dest='output_ptortho', default=None, help='Output for orthogroups') +parser.add_argument('--output_ptorthocs', dest='output_ptorthocs', default=None, help='Output for orthogroups with corresponding coding sequences') +parser.add_argument('--output_ptsco', dest='output_ptsco', default=None, help='Output for single copy orthogroups') + +args = parser.parse_args() + + +def get_stderr_exception(tmp_err, tmp_stderr, tmp_out, tmp_stdout, include_stdout=False): + tmp_stderr.close() + """ + Return a stderr string of reasonable size. + """ + # Get stderr, allowing for case where it's very large. + tmp_stderr = open(tmp_err, 'rb') + stderr_str = '' + buffsize = BUFF_SIZE + try: + while True: + stderr_str += tmp_stderr.read(buffsize) + if not stderr_str or len(stderr_str) % buffsize != 0: + break + except OverflowError: + pass + tmp_stderr.close() + if include_stdout: + tmp_stdout = open(tmp_out, 'rb') + stdout_str = '' + buffsize = BUFF_SIZE + try: + while True: + stdout_str += tmp_stdout.read(buffsize) + if not stdout_str or len(stdout_str) % buffsize != 0: + break + except OverflowError: + pass + tmp_stdout.close() + if include_stdout: + return 'STDOUT\n%s\n\nSTDERR\n%s\n' % (stdout_str, stderr_str) + return stderr_str + + +def move_directory_files(source_dir, destination_dir): + source_directory = os.path.abspath(source_dir) + destination_directory = os.path.abspath(destination_dir) + if not os.path.isdir(destination_directory): + os.makedirs(destination_directory) + for dir_entry in os.listdir(source_directory): + source_entry = os.path.join(source_directory, dir_entry) + shutil.move(source_entry, destination_directory) + + +def stop_err(msg): + sys.stderr.write(msg) + sys.exit(1) + + +def write_html_output(output, title, dir): + fh = open(output, 'wb') + fh.write('\n%s\n\n\n

\n

\n\n\n') + fh.close() + + +# Define command response buffers. +tmp_out = tempfile.NamedTemporaryFile().name +tmp_stdout = open(tmp_out, 'wb') +tmp_err = tempfile.NamedTemporaryFile().name +tmp_stderr = open(tmp_err, 'wb') +# Build the command line. +cmd = 'GeneFamilyClassifier' +cmd += ' --proteins %s' % args.input +cmd += ' --scaffold %s' % args.scaffold +cmd += ' --method %s' % args.method +cmd += ' --classifier %s' % args.classifier +cmd += ' --config_dir %s' % args.config_dir +cmd += ' --num_threads %d' % args.num_threads +if args.super_orthogroups is not None: + cmd += ' --super_orthogroups %s' % args.super_orthogroups +if args.single_copy_custom is not None: + cmd += ' --single_copy_custom %s' % args.single_copy_custom +if args.single_copy_taxa is not None: + cmd += ' --single_copy_taxa %d' % args.single_copy_taxa +if args.taxa_present is not None: + cmd += ' --taxa_present %d' % args.taxa_present +if args.orthogroup_fasta is None: + create_ortho_sequences = False +else: + create_ortho_sequences = True + cmd += ' --orthogroup_fasta' +if args.coding_sequences is None: + create_corresponding_coding_sequences = False +else: + create_corresponding_coding_sequences = True + cmd += ' --coding_sequences %s' % args.coding_sequences +# Run the command. +proc = subprocess.Popen(args=cmd, stderr=tmp_stderr, stdout=tmp_stdout, shell=True) +rc = proc.wait() +# Handle execution errors. +if rc != 0: + error_message = get_stderr_exception(tmp_err, tmp_stderr, tmp_out, tmp_stdout) + stop_err( error_message ) +# Handle hmmscan.log output. +if args.classifier in ['hmmscan', 'both']: + src_hmmscan_log = os.path.join(OUTPUT_DIR, 'hmmscan.log') + if os.path.exists(src_hmmscan_log): + if args.save_hmmscan_log is None: + os.remove(args.hmmscan_log) + else: + shutil.move(src_hmmscan_log, args.hmmscan_log) +# Handle orthogroups outputs. +if create_ortho_sequences: + if create_corresponding_coding_sequences: + out_file = args.output_ptorthocs + title = 'Orthogroups with corresponding coding sequences' + else: + out_file = args.output_ptortho + title = 'Orthogroups' + orthogroups_fasta_src_dir = os.path.join(OUTPUT_DIR, 'orthogroups_fasta') + orthogroups_fasta_dest_dir = out_file.files_path + move_directory_files(orthogroups_fasta_src_dir, orthogroups_fasta_src_dir) + write_html_output(out_file, title, orthogroups_fasta_dest_dir) +# Handle single copy orthogroup outputs. +if args.output_ptsco is not None: + single_copy_fasta_src_dir = os.path.join(OUTPUT_DIR, 'single_copy_fasta') + single_copy_fasta_dest_dir = args.output_ptsco.files_path + title = 'Single copy orthogroups' + move_directory_files(single_copy_fasta_src_dir, single_copy_fasta_dest_dir) + write_html_output(args.output_ptsco, title, single_copy_fasta_dest_dir)