perl/etc/justonce/bstannoA.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 "# Count qseqid sacc length evalue mismatch annot bitscore qstart qend sstart send btop sseqid\n";
  32
  33 sub deal($) {
  34         my $Aref=$_[0];
  35         my (%Dat,@Arr);
  36         for (@$Aref) {
  37                 push @{$Dat{$_->[2]}},$_;
  38         }
  39         for my $k (keys %Dat) {
  40                 my $Emin=(sort {$b->[10] <=> $a->[10] || $a->[7] <=> $b->[7]} @{$Dat{$k}})[0];  # len Desc, E Asc
  41                 push @Arr,[$Emin,scalar @{$Dat{$k}}];
  42         }
  43         for (@Arr) {
  44                 my $Emin=$_->[0];
  45                 my ($qseqid,$sseqid,$sacc,$qstart,$qend,$sstart,$send,$evalue,$bitscore,$score,$length,$pident,$nident,$mismatch,$positive,$gapopen,$gaps,$btop)=@$Emin;
  46                 my $annot=$sseqid;
  47                 $annot =~ s/^gi//;
  48                 if (exists $Anno{$annot}) {
  49                         $annot = '['.$Anno{$annot}.']';
  50                         print O join("\t",$_->[1],$qseqid,$sacc,$length,$evalue,$mismatch,$annot,$bitscore,$qstart,$qend,$sstart,$send,$btop,$sseqid),"\n";
  51                 } else {
  52                         print STDERR '^';
  53 warn "$annot\n";
  54                 }
  55         }
  56         return;
  57 }
  58
  59 my $lastID='';
  60 my @lastDat;
  61 open B,'<',$blastf or die "Error opening $blastf : $!\n";
  62 while (<B>) {
  63         chomp;
  64         my @t=split /\t/;
  65         if ($t[0] eq $lastID) {
  66                 push @lastDat,\@t;
  67         } else {
  68                 &deal(\@lastDat) if @lastDat>0;
  69                 @lastDat = (\@t);
  70                 $lastID = $t[0];
  71         }
  72 }
  73
  74 close B;
  75 close O;
  76