SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Blast HSP tiling BioPerl Module Query swttalyan Bioinformatics 0 08-24-2015 10:27 AM
bioperl blast parsing problem Mark Bioinformatics 9 04-27-2013 05:22 PM
Multiple blast results in the same object with BIOPERL angeloulivieri Bioinformatics 0 04-18-2013 02:46 AM
BLAST+ creating custom blast database and using blast+ filtering features deniz Bioinformatics 2 10-26-2012 12:04 PM
Parsing BLAST results using BioPerl Ben Saville Bioinformatics 8 08-24-2010 08:43 AM

Reply
 
Thread Tools
Old 10-26-2015, 10:29 AM   #1
uloeber
Member
 
Location: Germany

Join Date: Mar 2013
Posts: 35
Default Filtering BLAST results using BioPerl

Hi all,
I try to fix the following issue:
I have a BLAST output from the command line, using PacBio data as queries and a custom database. I have some filtering criteria for my hits, e.g. >db_1 should at least cover 80bp, while >db_2 should cover at least 1000bp.

In the end I'd like to have a graphical representation of my queries and the corresponding (filtered) hits.
A refinement would be, if the alignment of db_2 is in total longer than the alignment of db_3, db_3 should be discarded...

I tried to use Bio::Search::Result::ResultI first. I ended up with a very ugly script, which works and processes my custom blast tab, but I even failed to visualize my results.
Code:
while( my $result = $in->next_result ) {
  ## $result is a Bio::Search::Result::ResultI compliant object
  while( my $hit = $result->next_hit ) {
    ## $hit is a Bio::Search::Hit::HitI compliant object
    while( my $hsp = $hit->next_hsp ) {
      ## $hsp is a Bio::Search::HSP::HSPI compliant object
      if( $hit->name =~ m/someexpression/ ) {
       if($hit->name =~ m/db_2/ && $hsp->length('hit')>=90){
                my $aln = $hsp->get_aln;
                my $alnIO = Bio::AlignIO->new(-format =>"msf", -file => ">hsp.msf");
                $alnIO->write_aln($aln);
                print   $result->query_name, "\t",
                  $hit->name, "\t",
                  $hsp->length('total'), "\t",
                  $hsp->percent_identity, "\t",
                  $hsp->start('query'), "\t",
                  $hsp->end('query'), "\t",
                  $hsp->start('hit'), "\t",
                  $hsp->end('hit'), "\t",
                  $alnIO, "\n";
        }
       elsif($hit->name =~ m/db_2/ && $hsp->length('hit')>=80){
                my $aln = $hsp->get_aln;
                my $alnIO = Bio::AlignIO->new(-format =>"msf", -file => ">hsp.msf");
                $alnIO->write_aln($aln);
                print   $result->query_name, "\t",
                  $hit->name, "\t",
                  $hsp->length('total'), "\t",
                  $hsp->percent_identity, "\t",
...
Next thing I gave a try was Bio::SearchIO::Writer::HTMLResultWriter;
Code:
my $in = new Bio::SearchIO(     -format => 'blast', 
                                -file   => $ARGV[0]);


sub hsp_filter {
        my $hsp = shift;
        return 1 if ($hit->name =~ m/db_1/ && $hsp->length('hit')>=70);
        return 1 if ($hit->name =~ m/db_2/ && $hsp->length('hit')>=65);
 ...

}

my $writer = Bio::SearchIO::Writer::HTMLResultWriter ->new (-filters => {'HSP' => \&hsp_filter});
my $out = Bio::SearchIO ->new (-writer =>$writer);
$out ->write_result($in->next_result);
The problem is, the syntax doesn't work. Of cause, because I mixed up hits and hsp. Is it possible to combine hit and hsp filters?

I know, that this question is very specific. But maybe anybody struggles on similar issues parsing blast outputs, trying to filter them somehow and visualize the results. I hope I could make my problem clear and maybe anybody could give me a suggestion how to solve it.

Thanks in advance!

ps: it doesn't matter whether I parse xml, pure blast or tables, just want to get it to run
uloeber is offline   Reply With Quote
Old 10-29-2015, 09:18 AM   #2
uloeber
Member
 
Location: Germany

Join Date: Mar 2013
Posts: 35
Default

Okay, I fixed most of the issues... is anybody aware of how to write multiple panels in an output file using
while( my $result = $searchio->next_result )
with
Bio::Graphics::Panel->new
??

I will post my complete script, when I found a solution and am done
uloeber is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 04:48 PM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2018, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO