Read Blast Output Directly From Stdout With Bioperl

Question

[Pasta]

Hi there,

I finally decided myself to try Perl for a project and I got a question concerning BioPerl. I would like to read the result (blastxml) of a blast query directly from STDOUT. Is it possible ?

Here is my attempt:

   my $command = "echo '>$SeqName\n$_[0]' | blastall -a 4 -p blastn -d $db  $options";

   open $fh,"$command |" || die("cannot run fasta cmd of $command: $!\n");

   my $seqout = Bio::SeqIO->new(-fh => $fh,  -format => 'blastxml');

Also, I manage to catch STDOUT using IO::CaptureOutput, but I cant feed it into Bio::SeqIO. Maybe there is something than can be done with that.

Any help will be appreciated. Thanks

[Neilfws]

As Fwip says in their answer below: this will not work because you're using Bio::SeqIO. BLAST XML is not a sequence format. Suggest you read the Bioperl documentation and try some simple examples to get comfortable.

Answer 1

[Fwip]

First, I usually capture input like this:

my $command = "echo '>$SeqName\n$_[0]' | blastall -a 4 -p blastn -d $db  $options";
my $result = system($command);

And then $result has whatever came from STDOUT.

I am not aware of a SeqIO class for BlastXML input. Personally, I use the StandAloneBlastPlus module, but it looks like you are using the older interface to BLAST, in which case you may find more success using StandAloneBlast. I am not sure if these write to a temporary file or not.

The Bio::SearchIO module may be what you're looking for instead of SeqIO.

[Pasta]

Perfect ! I simply had to use SearchIO instead. Thanks for your time

Answer 2

[Damian Kao]

Most likely it is because you are feeding an incomplete XML file into the BioPerl parser.

XML format contain start and end tags that informs the parser of the data structure. At the very least, the stdout file will be missing the end tags for the topmost start tags (probably BlastOutput and BlastOutput_interations tags).

Here is an example of one iteration of the blast xml output:

<Iteration>
  <Iteration_iter-num>1</Iteration_iter-num>
  <Iteration_query-ID>Query_1</Iteration_query-ID>
  <Iteration_query-def>queryID</Iteration_query-def>
  <Iteration_query-len>3440</Iteration_query-len>
  <Iteration_hits></Iteration_hits>
  <Iteration_stat>
    <Statistics>
      <Statistics_db-num>407788</Statistics_db-num>
      <Statistics_db-len>154603577</Statistics_db-len>
      <Statistics_hsp-len>127</Statistics_hsp-len>
      <Statistics_eff-space>104767976519</Statistics_eff-space>
      <Statistics_kappa>0.041</Statistics_kappa>
      <Statistics_lambda>0.267</Statistics_lambda>
      <Statistics_entropy>0.14</Statistics_entropy>
    </Statistics>
  </Iteration_stat>
  <Iteration_message>No hits found</Iteration_message>
</Iteration>

If blastall is outputting xml lines one at a time, you will be feeding the parser one of the above lines at a time. Which will make no sense at all to the parser. I suggest you read up on the [?]XML format[?]

I would configure blastall to output a tab delimited format or text format for easier parsing.

If you really want to use XML, you can try to write a data buffer into your perl script. The data buffer will have to:

1. save the initial blast settings/version information
2. check if current multi-line xml iteration of blast is finished from stdout and save it to the buffer
3. append the initial blast settings/version information to the current iteration
4. append the appropriate end tags
5. feed the formatted data to the BioPerl parser

[Pasta]

No worries, the XML output is fine. I just want to feed the XML from STDIN to the BioPerl parser

[Damian Kao]

I understand that it is fine. But the way XML is parsed means you need to give it a complete data structure. STDOUT from blastall will feed the parser partial data.