4
$\begingroup$

I have run a blastx of metagenomic databases (raw illumina reads) using the nr database. Unfortunately, I forgot to add the --outfmt 6 argument to the code and got the traditional output.

Could I parse this output into --outfmt 6?

Improve this question
$\endgroup$
4
  • $\begingroup$ By "metagenomic databases" do you mean collections of reads or assembled contigs? BLAST searches of metagenomic reads vs nr is a very bad idea. Much better tools exist for classifying, etc. $\endgroup$
    – Daniel Standage
    Sep 6, 2017 at 18:32
  • $\begingroup$ @ Daniel Standage I refer to a collection of metagenomic reads. Could you suggest any other classification method for shotgun metagenomic reads? $\endgroup$
    – andresito
    Sep 6, 2017 at 20:24
  • $\begingroup$ Sure! Although you should take my recommendations with a grain of salt, since this isn't really my area of expertise. But tools like Kraken and Kaiju seem to be pretty popular. $\endgroup$
    – Daniel Standage
    Sep 6, 2017 at 23:41
  • $\begingroup$ I think the standard BLAST output is a moving target so the exact parsing will depend on the version of BLAST you use; anyway I would recommend re-running with the --outfmt 6 flag unless you don't have the computational resources. If it's too slow you can use DIAMOND as a replacement for BLAST with very similar commands $\endgroup$
    – Chris_Rands
    Sep 7, 2017 at 9:53

2 Answers 2

Reset to default
3
$\begingroup$

Biopython has a parser for Blast's default text format (which I assume is the one you have). Just use the blast-text string to specify this format in SearchIO's parse method.

from Bio import SearchIO
results = SearchIO.parse('myfile.txt', 'blast-text')
for entry in results:
    # extract the information you want
    query = entry.id
    for hit in entry.hits:
        subject = hit.id
        for hsp in hit.hsps:
            identity = hsp.ident_num *100 / entry.seq_len
            aln_len = hsp.aln.get_alignment_length()
            mismatches = aln_len - hsp.ident_num # not sure about this
            gaps = hsp.gap_num
            q_start = hsp.query_start + 1
            q_end = hsp.query_end
            s_start = hsp.hit_start + 1
            s_end = hsp.hit_end
            e_value = hsp.evalue
            bit_score = int(hsp.bitscore)
            # print in outfmt 6 format
            print('\t'.join([str(x) for x in [query,
                                      subject,
                                      identity,
                                      aln_len,
                                      mismatches,
                                      gaps,
                                      q_start,
                                      q_end,
                                      s_start,
                                      s_end,
                                      e_value,
                                      bit_score]]))

Have a look at Biopython's docs for more information on how to parse Blast results.

Improve this answer
$\endgroup$
2
  • $\begingroup$ This is a very versatile approach, thanks. $\endgroup$
    – bedeabc
    Sep 8, 2017 at 12:12
  • $\begingroup$ @andresito, please remember to accept the answer that helped you the most $\endgroup$
    – mgalardini
    Sep 13, 2017 at 18:53
3
$\begingroup$

The blast_formatter command allows you to produce output in any of the formats supported by BLAST. However, it requires that you run blastx with --outfmt 11 (ASN.1 format). I'm not sure there are any tools that will do the conversion from the full default BLAST report. This could potentially be scripted, but would be complex enough that you're probably better off re-running BLAST.

Improve this answer
$\endgroup$
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.