#!/usr/bin/perl
# idchains.perl version 2
# input: blat or blast f8 format

=item fixme: subset chains

WBGene00000003|WBGene00000004|WBGene00000006|WBGene00000007|WBGene00000008|WBGene00000009|WBGene00000010
WBGene00000004|WBGene00000006|WBGene00000007|WBGene00000008|WBGene00000010
WBGene00000005|WBGene00000006|WBGene00000007|WBGene00000008|WBGene00000009|WBGene00000010
            ^^^^ the following are essentially subsets of 1st;
            fix to remove these from .idchains
=cut

use Getopt::Long;
$isblat=1; $isblast=0; $removesubsets=1; $usebest=0; 
$eval= 1e-10; $bitscore=0; 
$align=0.1; # align=0.5 is too restrictive: Dmoj/Est6 tandems fail
$skipids="";
$reciprocal=1;

my $optok= GetOptions("blast!",\$isblast, "blat!",\$isblat, 
           "best!", \$usebest, "removesubsets!",\$removesubsets, "skipids=s",\$skipids,
           "reciprocal!",\$reciprocal,
           "eval=s",\$eval, "bitscore=s",\$bitscore, "align=s",\$align); 
die 
"usage: idchains.perl -blast|-blat [opts] < result.blast > result.idchains
opts: -[no]removesubsets -skipids=file -eval=$eval OR -bitscore=$bitscore (blast) OR -align=$align (blat)
      -[no]reciprocal
 " unless($optok and ($isblast || $isblat));

$noreciprocal= ! $reciprocal;
%skipids=();
if($skipids && open(F,$skipids)) { while(<F>){chomp; $skipids{$_}=1;} close(F); }

while(<>){
  chomp; @v=split"\t"; 
  unless($isblast) {    # .. blat format
    next unless(/^\d/);
    ($abmat,$aid,$as,$bid,$bs)= @v[0,9,10,13,14]; # blat psl format query,subject
    next if($skipids{$aid} or $skipids{$bid});
    $nself{$aid}++ and next if($aid eq $bid);
    next if( $abmat < $as*$align || $abmat < $bs*$align || abs($as-$bs)>20 );
    $bbits= $abmat; # for high score
  } else {  # .. blast format  
    next unless(/^\w/);
    ($aid,$bid,$beval,$bbits)= @v[0,1,10,11]; # blast format query,subject
    # using prot22_modXX.blastp out got dingbats on aid,bid : clean
    $aid =~ s/\W+$//; $bid =~ s/\W+$//;
    next if($skipids{$aid} or $skipids{$bid});
    $nself{$aid}++ and next if($aid eq $bid);
    if($bitscore) { next if($bbits < $bitscore); } else { next if ($beval > $eval); }
  }
  $pair{$aid}{$bid}++;  # want reciprocal align 
#   if($usebest) {
#     $pair{$aid}{$bid}= $bbits if ($pair{$aid}{$bid} < $bbits);
#   } else {
#    $pair{$aid}{$bid}++;  # want reciprocal align 
#  }
}
END{ 
  my @aid= sort keys %pair; 
  my $ng=scalar(@aid); 
  my @chains;
  
  foreach my $aid (@aid) { 
    my @bid;
    @bid= sort keys %{$pair{$aid}}; 
#     if($usebest) {
#       @bid= sort{ $pair{$aid}{$b} <=> $pair{$aid}{$a} } keys %{$pair{$aid}}; 
#     } else {
#       @bid= sort keys %{$pair{$aid}}; 
#     }
  
   my @cid=(); 
   BLIST: foreach my $bid (@bid) { 
    next unless($removesubsets or !$didab{$aid}{$bid});
    if($noreciprocal or $pair{$bid}{$aid} ) { #&& !$didab{$aid}{$bid}
      push(@cid, $bid); 
      # last BLIST if($usebest); # skip didab filter? no effect
      # usebest has strange result: very few (50, 200, 800) gene pairs this way
      # versus 1000 .. 6000 w/o; something wrong about its use here.
      $didab{$aid}{$bid}=$didab{$bid}{$aid}=1;  # check recip match
      }
     } 
    
   my %bid= map{ $_=>1; } $aid,@cid; 
   @bid= sort keys %bid;
   
   if ($removesubsets and @bid>1) {
    my $nid= @bid;
    foreach my $chain (@chains) {
      my $did=0; map{ $did++ if $bid{$_};  } @$chain;
      @bid=() and last if($did >= $nid);
    }
    push(@chains, \@bid) if(@bid>1);
   }
   
   if (@bid>1) { print join("|",@bid),"\n"; $na++; $nb += @bid; @nc{@bid}=1; } 
   # count mean, median, .. distrib/ quartiles of chain sizes? e.g. n clusters of 2..9,10,20..90,100..500 
   } 
 $nc=scalar(keys %nc);  $nself= scalar(keys %nself);
 $nsing= $nself-$ng;
 print "# nself=$nself; npaired=$ng; nsingle=$nsing; duplicates: ngroup=$na; nuniqgene=$nc; ngroupedgene=$nb; rmsubsets=$removesubsets \n"; 
}