Mercurial > repos > peterjc > blast2go
diff tools/ncbi_blast_plus/blast2go.py @ 1:0f159cf346c8
Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
author | peterjc |
---|---|
date | Tue, 07 Jun 2011 16:29:28 -0400 |
parents | cd52c931b325 |
children | a5594772282c |
line wrap: on
line diff
--- a/tools/ncbi_blast_plus/blast2go.py Tue Jun 07 16:28:31 2011 -0400 +++ b/tools/ncbi_blast_plus/blast2go.py Tue Jun 07 16:29:28 2011 -0400 @@ -6,8 +6,15 @@ * Blast2GO properties filename (settings file) * Output tabular filename +Sadly b2g4pipe v2.3.5 cannot cope with current style large BLAST XML +files (e.g. from BLAST 2.2.25+), so we have to reformat these to +avoid it crashing with a Java heap space OutOfMemoryError. + +As part of this reformatting, we check for BLASTP or BLASTX output +(otherwise raise an error), and print the query count. + It then calls the Java command line tool, and moves the output file to -the location Galaxy is expecting. +the location Galaxy is expecting, and removes the tempory XML file. """ import sys import os @@ -27,12 +34,66 @@ xml_file, prop_file, tabular_file = sys.argv[1:] +#We should have write access here: +tmp_xml_file = tabular_file + ".tmp.xml" + if not os.path.isfile(xml_file): stop_err("Input BLAST XML file not found: %s" % xml_file) if not os.path.isfile(prop_file): stop_err("Blast2GO configuration file not found: %s" % prop_file) +def prepare_xml(original_xml, mangled_xml): + """Reformat BLAST XML to suit Blast2GO. + + Blast2GO can't cope with 1000s of <Iteration> tags within a + single <BlastResult> tag, so instead split this into one + full XML record per interation (i.e. per query). This gives + a concatenated XML file mimicing old versions of BLAST. + + This also checks for BLASTP or BLASTX output, and outputs + the number of queries. Galaxy will show this as "info". + """ + in_handle = open(original_xml) + footer = " </BlastOutput_iterations>\n</BlastOutput>\n" + header = "" + while True: + line = in_handle.readline() + if not line: + #No hits? + stop_err("Problem with XML file?") + if line.strip() == "<Iteration>": + break + header += line + + if "<BlastOutput_program>blastx</BlastOutput_program>" in header: + print "BLASTX output identified" + elif "<BlastOutput_program>blastp</BlastOutput_program>" in header: + print "BLASTP output identified" + else: + in_handle.close() + stop_err("Expect BLASTP or BLASTX output") + + out_handle = open(mangled_xml, "w") + out_handle.write(header) + out_handle.write(line) + count = 1 + while True: + line = in_handle.readline() + if not line: + break + elif line.strip() == "<Iteration>": + #Insert footer/header + out_handle.write(footer) + out_handle.write(header) + count += 1 + out_handle.write(line) + + out_handle.close() + in_handle.close() + print "Input has %i queries" % count + + def run(cmd): #Avoid using shell=True when we call subprocess to ensure if the Python #script is killed, so too is the child process. @@ -44,10 +105,11 @@ stdout, stderr = child.communicate() return_code = child.returncode if return_code: + cmd_str = " ".join(cmd) if stderr and stdout: - stop_err("Return code %i from command:\n%s\n\n%s\n\n%s" % (return_code, err, stdout, stderr)) + stop_err("Return code %i from command:\n%s\n\n%s\n\n%s" % (return_code, cmd_str, stdout, stderr)) else: - stop_err("Return code %i from command:\n%s\n%s" % (return_code, err, stderr)) + stop_err("Return code %i from command:\n%s\n%s" % (return_code, cmd_str, stderr)) #For early diagnostics, else: print stdout @@ -56,16 +118,25 @@ if not os.path.isfile(blast2go_jar): stop_err("Blast2GO JAR file not found: %s" % blast2go_jar) -#We will have write access whereever the output should be, +prepare_xml(xml_file, tmp_xml_file) +#print "XML file prepared for Blast2GO" + +#We will have write access wherever the output should be, #so we'll ask Blast2GO to use that as the stem for its output #(it will append .annot to the filename) cmd = ["java", "-jar", blast2go_jar, - "-in", xml_file, + "-in", tmp_xml_file, "-prop", prop_file, - "-out", tabular_file, - "-a"] + "-out", tabular_file, #Used as base name for output files + "-a", # Generate *.annot tabular file + #"-img", # Generate images, feature not in v2.3.5 + ] +#print " ".join(cmd) run(cmd) +#Remove the temp XML file +os.remove(tmp_xml_file) + out_file = tabular_file + ".annot" if not os.path.isfile(out_file): stop_err("ERROR - No output annotation file from Blast2GO")