1
$\begingroup$

I have to use the old version of Blast (not Blast+) to compare DNA sequences and I have to know the number of all hits. They are not given in the header and I have to know them for many many of blast outputs and its nearly impossible to count them (or to export them in excel). I think I probably have to write a 'Perl Script' and I did my best but I am not familiar with programming and it is not working. Furthermore I had a look on the internet and they are always comparing it to BIO Search. I have e code, which for me looks pretty good,but it also done a comparison with BIO Search, what I don't need/want.

use strict;
use warnings;
use Bio::SearchIO;

Usage information

die "Usage: $0 <BLAST-report-file> <number-of-top-hits> <output-file>\n", if (@ARGV != 3);

my ($infile,$numHits,$outfile) = @ARGV;
    print "Parsing the BLAST result ...";
    my $in = Bio::SearchIO->new(-format => 'blast', -file => $infile);
    open (OUT,">$outfile") or die "Cannot open $outfile: $!";

# print the header info for tab-deliminated columns
print OUT "query_name\tquery_length\taccession_number\tlength\tdescription\tE value\tbit score\tframe\tquery_start\t";
print OUT "query_end\thit_start\thit_end\tpositives\tidentical\n";

# extraction of information for each result recursively
while ( my $result = $in->next_result ) {
    # the name of the query sequence
    print OUT $result->query_name . "\t";

        # the length of the query sequence
        print OUT $result->query_length;

        # output "no hits found" if there is no hits
        if ( $result->num_hits == 0 ) {
    		print OUT "\tNo hits found\n";
        	} else {
    		my $count = 0;

                # process each hit recursively
        while (my $hit = $result->next_hit) {
            print OUT "\t" if ($count > 0);
                            # get the accession numbers of the hits
    			print OUT "\t" . $hit->accession . "\t";
                        # get the lengths of the hit sequences
                        print OUT $hit->length . "\t";
                            # get the description of the hit sequences
    			print OUT $hit->description . "\t";
                        # get the E value of the hit
            print OUT $hit->significance . "\t";
                            #get the bit score of the hit
    			print OUT $hit->bits . "\t";

                        my $hspcount = 0;

                        # process the top HSP for the top hit
            while (my $hsp = $hit->next_hsp) {
                            print OUT "\t\t\t\t\t\t\t", if ($hspcount > 0);
                            	# get the frame of the query sequence
    				print OUT $hsp->query->frame . "\t";
                                # get the start and the end of the query sequence in the alignment
                print OUT $hsp->start('query') . "\t" . $hsp->end('query'). "\t";
                                # get the start and the end of the hit sequence in the alignment
                print OUT $hsp->start('hit') . "\t" . $hsp->end('hit') . "\t";
                                # get the similarity value
                printf OUT "%.1f" , ($hsp->frac_conserved * 100);
    				print OUT "%\t";
                                    # get the identity value
    				printf OUT "%.1f" , ($hsp->frac_identical * 100);
                    print OUT "%\n";
                                $hspcount++;
                            }
    			$count++;

                        # flow control for the number of hits needed
            last if ($count == $numHits);
        }
        }
}
close OUT;
print " DONE!!!\n";

May someone help me, please. I need a program, which first tells me the header and then the number of hits. I am a little bit desperate and would be so happy if someone can help me. That would be so nice. Greetings Tina

$\endgroup$
8
  • $\begingroup$ Could you provide more info on how it's not working? I mean does it run but prints wrong output, or it won't run at all? $\endgroup$ Mar 20 '15 at 10:58
  • $\begingroup$ If your issue is related to your perl script then this is not on topic in biology-se. If you are interested in general procedure to get total number of hits you can see the above link. $\endgroup$
    – WYSIWYG
    Mar 20 '15 at 11:40
  • $\begingroup$ Can you paste the output format? $\endgroup$
    – WYSIWYG
    Mar 20 '15 at 11:46
  • 1
    $\begingroup$ Thanks for the fast answer. The script i didn't write by myself, I just found it on the internet. I don't get any output, because it works with BioSearch (you have to download). But I think with little changes it could work...but i am not sure (4th line??). Maybe there is a easier way. I just wanted to ask here, if someone has an idea how it will work or how i can reach my aim, counting the hits. In the link you posted is written, that the hits are in the header, but there are not, in my case. $\endgroup$
    – Nina
    Mar 20 '15 at 12:12
  • 1
    $\begingroup$ The scrip won't run because it has a missing package that it wants to use (use Bio::SearchIO;) so you either download it , or use a different script $\endgroup$ Mar 20 '15 at 12:19
0
$\begingroup$

This script, named blast2table.pl is from the code accompanying the Blast book by Korf & Yandell (there is a 3rd author whom I haven't met).

#!/usr/bin/perl -w
use strict;
use Getopt::Std;
use vars qw($opt_p $opt_b $opt_e $opt_m $opt_n);
    getopts('p:b:e:m:n:');
    my $PERCENT = $opt_p ? $opt_p : 0;
my $BITS    = $opt_b ? $opt_b : 0;
    my $EXPECT  = $opt_e ? $opt_e : 1e30;
my $START   = $opt_m ? $opt_m : 0;
    my $END     = $opt_n ? $opt_n : 1e30;

my ($Query, $Sbjct);
my $HSP = "";
    while (<>) {
        if    (/^Query=\s+(\S+)/) {outputHSP(); $Query = $1}
        elsif (/^>(\S+)/)         {outputHSP(); $Sbjct = $1}
        elsif (/^ Score = /) {
            outputHSP();
            my @stat = ($_);
            while (<>) {
                last unless /\S/;
                push @stat, $_
            }
            my $stats = join("", @stat);
            my ($bits) = $stats =~ /(\d\S+) bits/;
            my ($expect) = $stats =~ /Expect\S* = ([\d\.\+\-e]+)/;
	     	$expect = "1$expect" if $expect =~ /^e/;
            my ($match, $total, $percent) = $stats =~ /Identities = (\d+)\/(\d+) \((\d+)%\)/;
            my $mismatch = $total - $match;

        $HSP = {bits => $bits, 
                    expect => $expect, 
                    mismatch => $mismatch, 
                    percent => $percent, 
                    q_begin => 0, 
                    q_end => 0, 
                    q_align => 0, 
                    s_begin => 0, 
                    s_end => 0, 
                    s_align => ""};
        }
        elsif (/^Query:\s+(\d+)\s+(\S+)\s+(\d+)/) {
            $HSP->{q_begin}  = $1 unless $HSP->{q_begin};
            $HSP->{q_end}    = $3;
            $HSP->{q_align} .= $2;
        }
        elsif (/^Sbjct:\s+(\d+)\s+(\S+)\s+(\d+)/) {
            $HSP->{s_begin}  = $1 unless $HSP->{s_begin};
            $HSP->{s_end}    = $3;
            $HSP->{s_align} .= $2;
    }
}
outputHSP();

sub outputHSP {
    return unless $HSP;
        return if $HSP->{percent}  < $PERCENT;
        return if $HSP->{bits}     < $BITS;
        return if $HSP->{expect}   > $EXPECT;
        return if ($HSP->{q_begin} < $START or $HSP->{q_end} < $START);
        return if ($HSP->{q_begin} > $END   or $HSP->{q_end} > $END);
        foreach my $field ('percent', 'q_align', 'mismatch', 's_align', 'q_begin', 'q_end', 's_begin', 's_end', 'expect', 'bits'){
        print "$field not defined\n" if not defined $HSP->{$field};
        }
        print join("\t", $Query, $Sbjct, $HSP->{percent}, length($HSP->{q_align}), $HSP->{mismatch}, countGaps($HSP->{q_align}) + countGaps($HSP->{s_align}), $HSP->{q_begin}, $HSP->{q_end}, $HSP->{s_begin}, $HSP->{s_end}, $HSP->{expect}, $HSP->{bits}), "\n";
    $HSP = "";
}

sub countGaps {
    my ($string) = @_;
        my $count = 0;
    while ($string =~ /\-+/g) {$count++}
        return $count;
}
$\endgroup$

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