Bioperl:解析blast结果


下面是一个bioperl小脚本, 用于解析blast的结果.

#!/usr/bin/perl -w
use strict;
use Bio::SearchIO;

my $filename = 'MET30-vs-swissprot.3.WUBLASTP';
my $input = Bio::SearchIO->new(-format => 'blast',
          -file   => $filename);
#GOAL:
# Print out the query, hit names and percent identity, percent query aligned
while( my $r = $input->next_result ) {
  my $qname = $r->query_name;
  my $qlen  = $r->query_length;
  while( my $hit = $r->next_hit ) {
    my $hname = $hit->name;
    while( my $hsp = $hit->next_hsp ) {
      my $percent_id = $hsp->frac_identical * 100;
      my $query_aln_len = $hsp->length('query');
      # another way to do this
      # my $query_aln_len = $hsp->query->length;
      print join("\t", $qname, $hname, sprintf("%.1f",$percent_id),
               sprintf("%.2f",100 * ( $query_aln_len / $qlen))), "\n";
    }
  }
}

一条回应:“Bioperl:解析blast结果”

  1. 你好,我现在想通过blast结果判断基因的拷贝数多少。
    请问在score和e-value上应该怎样去cutoff值作为单拷贝和多拷贝的界限。谢谢~