perl/etc/justonce/bstanno.pl

   1 #!/bin/env perl
   2 =pod
   3 Author: Hu Xuesong @ BIOPIC <galaxy001@gmail.com>
   4 Version: 1.0.0 @ 20120330
   5 =cut
   6 use strict;
   7 use warnings;
   8 #use Data::Dump qw(ddx);
   9
  10 die "Usage: $0 <nt header> <blast file>\n" if @ARGV<2;
  11 my $header=shift;
  12 my $blastf=shift;
  13
  14 my %Anno;
  15 open H,'<',$header or die "Error opening $header : $!\n";
  16 while (<H>) {
  17         chomp;
  18         s/^>// or warn "[x][$_]\n";
  19         my @dat=split /\001/;
  20         for (@dat) {
  21                 /^gi(\S*) (.*)/ or die;
  22                 $Anno{$1}=$2;
  23 #print "[@dat] ----- $1,$2\n";
  24         }
  25 }
  26 close H;
  27 warn "[!]Load Hearers done.\n";
  28
  29 $|=1;
  30 open O,'>',$blastf.'.anot' or die "Error opening $blastf.anot : $!\n";
  31 print O "# qseqid sacc length evalue mismatch annot bitscore qstart qend sstart send btop sseqid\n";
  32
  33 sub deal($) {
  34         my $Aref=$_[0];
  35         my $Emin=(sort {$b->[10] <=> $a->[10] || $a->[7] <=> $b->[7]} @$Aref)[0];       # len Desc, E Asc
  36 #ddx $Aref;
  37 #print '-'x5,"\n";
  38 #ddx $Emin;
  39 #print '-'x75,"\n";
  40         my ($qseqid,$sseqid,$sacc,$qstart,$qend,$sstart,$send,$evalue,$bitscore,$score,$length,$pident,$nident,$mismatch,$positive,$gapopen,$gaps,$btop)=@$Emin;
  41         my $annot=$sseqid;
  42         $annot =~ s/^gi//;
  43         if (exists $Anno{$annot}) {
  44                 $annot = '['.$Anno{$annot}.']';
  45                 print O join("\t",$qseqid,$sacc,$length,$evalue,$mismatch,$annot,$bitscore,$qstart,$qend,$sstart,$send,$btop,$sseqid),"\n";
  46         } else {
  47                 print STDERR '^';
  48 warn "$annot\n";
  49         }
  50         return;
  51 }
  52
  53 my $lastID='';
  54 my @lastDat;
  55 open B,'<',$blastf or die "Error opening $blastf : $!\n";
  56 while (<B>) {
  57         chomp;
  58         my @t=split /\t/;
  59         if ($t[0] eq $lastID) {
  60                 push @lastDat,\@t;
  61         } else {
  62                 &deal(\@lastDat) if @lastDat>0;
  63                 @lastDat = (\@t);
  64                 $lastID = $t[0];
  65         }
  66 }
  67
  68 close B;
  69 close O;
  70