diff --git a/tools/metag_tools/megablast_wrapper.py b/tools/metag_tools/megablast_wrapper.py index da1a7ef0890..44e8d5306a5 100644 --- a/tools/metag_tools/megablast_wrapper.py +++ b/tools/metag_tools/megablast_wrapper.py @@ -1,53 +1,15 @@ #! /usr/bin/python """ run megablast for metagenomics data -database builds file -usage: %prog database/reference query output -Wen-Yu Chung """ import sys, os, tempfile DB_LOC = "/depot/data2/galaxy/blastdb.loc" -def print_bed_format(fields, file_handle): - #print " ".join(fields) - print >> file_handle, "\t".join(fields) - return 0 - -def print_maf_format(fields, file_handle): - score = fields.pop(0) - print >> file_handle, "a score=%s" %score - print >> file_handle, "s %s" %" ".join(fields) - print >> file_handle - return 0 - -def parse_megablast_output(filename, output_format, file_handle): - fields = [] - file_to_be_parsed = open(filename, 'r') - for i, line in enumerate(file_to_be_parsed): - line = line.strip('\r\n') - if (not line.startswith("#")): - [query_id, subject_id, iden, align_length, mismatches, gaps, q_start, q_end, s_start, s_end, evalue, bit_score] = line.split() - if (int(s_start) > int(s_end)): - strand = "-" - temp = s_start - s_start = s_end - s_end = temp - else: - strand = "+" - if (output_format == "bed"): - fields = [subject_id, s_start, s_end, query_id, str(0), strand] - print_bed_format(fields, file_handle) - else: - fields = [bit_score, subject_id, s_start, str(int(align_length)-int(gaps)), strand, "srcSize", "text"] - print_maf_format(fields, file_handle) - return 0 - def __main__(): # file I/O db_build = sys.argv[1] query_filename = sys.argv[2] - #output_format = sys.argv[3] output_filename = sys.argv[3] # megablast parameters @@ -63,44 +25,43 @@ def __main__(): db = {} db_file = open(DB_LOC, "r") for i, line in enumerate(db_file): - line = line.strip('\r\n') + line = line.rstrip('\r\n') fields = line.split() db[(fields[0])] = [] for j in xrange(1, len(fields)): db[(fields[0])].append(fields[j]) - #print "\n".join(db[(db_build)]) # prepare to run megablast if os.system('which megablast 2>&1'): print >> sys.stderr, "Cannot locate megablast." for chunk in db[(db_build)]: - #chunk = db[(db_build)][0] # test - #if 1: # test - #if os.path.exists('/usr/bin/megablast'): megablast_output_file = tempfile.NamedTemporaryFile('w') megablast_output_filename = megablast_output_file.name megablast_arguments = ["megablast", "-d", chunk, "-i", query_filename, "-o", megablast_output_filename] - megablast_parameters = ["-m", "8", "-D", "3", "-a", "1"] + megablast_parameters = ["-m", "8", "-D", "3", "-a", "8"] megablast_user_inputs = ["-W", mega_word_size, "-p", mega_iden_cutoff, "-t", mega_disc_word, "-N", mega_disc_type, "-F", mega_filter] megablast_command = " ".join(megablast_arguments) + " " + " ".join(megablast_parameters) + " " + " ".join(megablast_user_inputs) + " 2>&1" - #print megablast_command os.system(megablast_command) - #parse_megablast_output(megablast_output_filename, output_format, output_file) - megablast_output_handle = open(megablast_output_filename, 'r') - for i, line in enumerate(megablast_output_handle): - line = line.strip('\r\n') + for i, line in enumerate( file(megablast_output_filename) ): + line = line.rstrip('\r\n') fields = line.split() if (not line.startswith("#")): # replace subject id with gi number # remove this after re-build blastdb subject_id_fields = fields[1].split('|') - gi = subject_id_fields[1] + if len(subject_id_fields) > 1: gi = subject_id_fields[1] + else: gi = subject_id_fields[0] fields[1] = gi print >> output_file, "\t".join(fields) - - megablast_output_file.close() output_file.close() - + + # megablast generates a file called error.log, if empty, delete it, if not, show the contents + if os.path.exists('./error.log'): + for i, line in enumerate( file('./error.log') ): + line = line.rstrip('\r\n') + print >> sys.stdout, line + os.remove('./error.log') + if __name__ == "__main__" : __main__() \ No newline at end of file diff --git a/tools/metag_tools/megablast_wrapper.xml b/tools/metag_tools/megablast_wrapper.xml index 9aca6b9f490..22a32c618a9 100644 --- a/tools/metag_tools/megablast_wrapper.xml +++ b/tools/metag_tools/megablast_wrapper.xml @@ -1,19 +1,12 @@ - on genomic data + for Metagenomics Projects megablast_wrapper.py $source_select $input_query $output1 $word_size $iden_cutoff $disc_word $disc_type $filter_query - - - - + + @@ -33,12 +26,10 @@ - .. class:: warningmark