0%

Blast+ xml格式解读

本地BLAST比对后,如果使用outfmt 5参数的话,会产生一个xml格式的文件,里面的比对信息不像tabular(outfmt 6)那样简显,但是对比信息却很完整。简单列举一些常用的信息。 1. 每条序列的所有对比信息是以标记<Iteration>开始,</Iteration>结束,之间框在一起的则是一条序列的所有比对结果信息 2. 后面接的数字表示序列的编号,1代表第一条序列 3. 序列的header 4. 序列的长度 2. 每条序列的比对结果是以标记<Iteration_hits>开始,</Iteration_hits>结束,框架下面还有<Hit>表示序列的每个比对结果(是指一条序列比对上对应的一条满足阈值的序列),而<Hsp>则表示每个比对结果中的某一块的比对结果(比如一条序列上有好几处跟目标序列比对上了) 1. <Hsp_num>比对上的序列的编号 2. 比对结果的打分 3. <Hsp_evalue>比对结果的evalue值 4. 目标序列比对上的开始位置 5. 目标序列比对上的结束位置 6. 比对上序列的开始位置 7. 比对上序列的结束位置 8. 目标序列翻译的起始位置 9. 比对上序列翻译的起始位置 10. 比对上序列那部分的长度 11. <Hsp_qseq>目标序列比对上的那部分序列信息 12. <Hsp_hseq>比对上序列的那部分序列信息 3. 简单的归纳就这样了,具体可翻看文档 4. 知道这些标记的含义,我们就可以用perl将自己需要的信息提取出来,然后输出成表格。

    #!/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转载请注明出处