Thank you very much for this quick and thorough answer!
Hi,
first of all, I have to admit I am new to Bioperl and Perl in general, so please be patient if I am asking a stupid question.
I want to calculate how much of my query sequence (in %) are matched by the hits in my BLAST output. A formula would look like this: (length of the BLAST alignment) / (length of the query sequence)
The length of the query is given in the header of each BLAST result as "Length = .." .
Currently I am accessing the data in the BLAST output like this (the e-value and other attributes works fine):
#!/usr/bin/perl
use strict;
use warnings;
use Bio::SearchIO;
my $report_obj = new Bio::SearchIO(-format => 'blast',
-file => 'test.blastx');
while (my $result = $report_obj -> next_result) {
while (my $hit = $result -> next_hit) {
# get the e-value
my $evalue = $hit -> expect;
print "\n" . 'e: ' . $evalue . "\n";
while ( my $hsp = $hit->next_hsp ) {
# get the query length
my $length_q = $hsp -> length('query');
print 'query length: ' . $length_q . "\n";
}
}
}
Does $hsp also have an attribute like query length, if so: how do I access it?
Or is there another way of getting the information I want to extract?
I appreciate every help :)
1 answer
All of the methods for Bio::SearchIO are summarised in the HOW-TO.
The length of the query sequence (as opposed to the length of the query in HSPs) is found by:
$result->query_length
What you are extracting using:
$hsp->length('query')
is the length of query sequence participating in the alignment, minus gaps.
"Length of the BLAST alignment" is a little imprecise. There may be > 1 HSP for each query/hit pair (in different frames, for example). You probably do not want to sum the length of all HSPs. This is why we typically deal with hits on a "per HSP" basis. So the closest thing to what you described: length of the BLAST alignment / length of the query sequence is probably:
$hsp->length('hit') / $result->query_length
assuming you do not want to count gaps as HSP length.
Anyway - have a good read of the HOW-TO and everything will be made clear.
Log in to answer this question.