本地BLAST比对后,如果使用outfmt 5参数的话,会产生一个xml格式的文件,里面的比对信息不像tabular(outfmt 6)那样简显,但是对比信息却很完整。简单列举一些常用的信息。 1. 每条序列的所有对比信息是以标记<Iteration>开始,</Iteration>结束,之间框在一起的则是一条序列的所有比对结果信息 2.
#!/usr/bin/perl -w
use strict;
open my $fh, "blast.xml" or die "Cannot open the blast file\n";
local $/ = "<\/Iteration>\n";
while(<$fh>){
my $query = $1 if ($_ =~ /<Iteration_query-def>(\S+)<\/Iteration_query-def>/);
my $query_len = $1 if ($_ =~ /<Iteration_query-len>(\S+)<\/Iteration_query-len>/);
while ($_ =~ /<Hit>(.*?)<\/Hit>/sg){
chomp(my $hit = $1);
my $id = $1 if ($hit =~ /<Hit_id>(\S+)<\/Hit_id>/);
my $def = $1 if ($hit =~ /<Hit_def>(.*?)<\/Hit_def>/);
my $num = 0;
while ($hit =~ /<Hsp>(.*?)<\/Hsp>/sg){
$num++;
my $hsp = $1;
my $query_from = $1 if ($hsp =~ /<Hsp_query-from>(\S+)<\/Hsp_query-from>/);
my $query_to = $1 if ($hsp =~ /<Hsp_query-to>(\S+)<\/Hsp_query-to>/);
my $hit_frame = $1 if ($hsp =~ /<Hsp_query-frame>(\S+)<\/Hsp_query-frame>/);
my $qseq = $1 if ($hsp =~ /<Hsp_qseq>(\S+)<\/Hsp_qseq>/);
my $align_len = $1 if ($hsp =~ /<Hsp_align-len>(\S+)<\/Hsp_align-len>/);
print ">$query;orf$num len=$align_len frame:$hit_frame start:$query_from end:$query_to $id $def\n";
}
}
}
close $fh;
最后的输出打印出来的结果是临时的,可以随意修改
本文出自于http://www.bioinfo-scrounger.com转载请注明出处