#!/usr/bin/perl
# tanpatt.pl : tandem duplicate gene run patterns
# cat dpul_gnoannot_pdup.gff | sort -k1,1 -k4,4n | perl tanpatt.pl | more

sub putc {
 print shift; $np++; print "\n" if ($np % 100 == 0);
}

sub binc { 
  my $b=shift; 
  if($b>1000){$b=1000*int($b/1000);}
  elsif($b>100){$b=100*int($b/100);} 
  elsif($b>10){$b=10*int($b/10);}
  return $b;
}
  
sub countc {
  my $c=shift; # or $po?
  if($c =~ /[A-Z]/) { $sa++; $thisl= $c; }
  elsif($c =~ /[o\.]/) { $so++;  }
  if($c eq $lastc) { # 1let run
    if($c =~ /[A-Z]/) { $sa1++; $ra++;  }
    elsif($c =~ /[o\.]/) { $so1++; $ro++; }
  } else {
    if($lastc =~ /[A-Z]/) { $lastl= $lastc; if($ra) { $lra{ binc($ra)}++; $nra++; $sra+=$ra; } $ra=0; }
    elsif($lastc =~ /[o\.]/) { if($ro) { $lro{ binc($ro)}++; $nro++; $sro+=$ro; } $ro=0; }
    
    #    if( $rh{$c} && $rh{$c1} ) { $r2++; }
    if( ($lastc =~ /[A-Z]/ || $c =~ /[A-Z]/) && $lastl && $thisl && $lastl ne $thisl ) {
      if( grep({$lastl eq $_} @lastc) && grep({$thisl eq $_} @lastc) ) { $sa2++; $ra2++; }
      shift(@lastc); push(@lastc,$c); # prefilled w/ WINDOW x .; change only at transition??
      }
  }

  if($c) {
    # shift(@lastc); push(@lastc,$c);
    $sc++;
    }
  $lastc= $c;
}

sub putmeans {
  $pa= $sa/$sc; $po= $so/$sc;  # useful
  $pa1= $sa1/$sc; $po1= $so1/$sc; #? not useful
  $pa2= $sa2/$sc;
  
  $mra= ($nra>0) ? $sra/$nra : 0; # useful
  $mro= ($nro>0) ? $sro/$nro : 0;
  
  ($rin,$ris)=(0,0);
  @runs= sort{ $a <=> $b } keys %lra;
  foreach $r (@runs) {
    $rin++ if ($lra{$r}); # $rin += $lra{$r};
    $ris += $r * $lra{$r};
  }
  $mria= ($rin)? $ris/$rin : 0;
  $runla= join(",", map{ "$_:$lra{$_}" } @runs[0..15]);
  
 ($rin,$ris)=(0,0);
  @runs= sort{ $a <=> $b } keys %lro;
  foreach $r (@runs) {
    $rin++ if ($lro{$r}); # $rin += $lra{$r};
    $ris += $r * $lro{$r};
  }
  $mrio= ($rin)? $ris/$rin : 0;
  $runlo= join(",", map{ "$_:$lro{$_}" } @runs[0..15]);
  
  %v=( 
    # meanA => $ma, meanO => $mo,  
    # sdA => $da, sdO => $do,
    meanrun_AA => $mra, #? redundant w/ pct_A
    meanrun_oo => $mro, 
    run_AA => $mria,
    run_oo => $mrio,
    # runAB => $mr2,
    pct_A => $pa, 
    pct_O => $po,
    pct_AA => $pa1, 
    pct_AB => $pa2,
    pct_OO => $po1,
    nGene => $sc,  
    #nline => $nl,
    );
  print "# statistics of runs\n# ";
  foreach $v (sort keys %v) { $fm=($v =~/pct/)?".3f":".2f"; printf "$v: %$fm, ", $v{$v}; } 
  print"\n";
  print "# runl_AA: $runla\n";
  print "# runl_oo: $runlo\n";
  
}

sub putmark {
  my $po= shift;
  if($po eq ".") { $c=$po; } 
  elsif($c{$po}) { $c=$c{$po}; }
  elsif($pc{$po}>1) {  $c= shift @cfree; $c||="x"; $c{$po}=$c; } 
  else { $c= "o"; }
  unless( $n<=$WINDOW ) { putc($c); countc($c) }
}

# fixme: reset at scaffold
sub resetruns {
  while( $po= shift @pw ) { putmark($po); }
  $didreset= $n > $WINDOW;
  countc(0) if $didreset;
  %pc=(); $n=0; %c=(); 
  $lastl = $thisl = 0;
  @lastc=(); push(@lastc, (".") x 10); 
  @pw=();    push(@pw, (".") x $WINDOW); 
  @cfree=(); foreach my $i (0..25) { push(@cfree, chr($i+65)); }
  return $didreset;
}

BEGIN{ 
  $WINDOW=20;
  resetruns();
  #push(@pw, (".") x $WINDOW); 
  #foreach $i (0..25) { push(@cfree, chr($i+65)); }
  print <<"EOT";
# Pattern of Tandem duplicate gene runs. 
# 
# This list shows duplicate gene order in schematic,
# for genes classed as the same paralog in runs (A..Z Letter), 
# or a paralog not in a run (o), or singleton gene (.)
# Scaffold breaks are shown as (^).
# 
# A sliding window of $WINDOW genes is used, and Letters
# @cfree
# are recycled, so "AA" following "AA" by $WINDOW or more
# denote different paralogs.

EOT
}

## cele has alt-tr; need to filter overlapped genes ****

while (<>) {
  next unless(/^\w/);
  my($ref,$src,$tp,$rb,$re,$rp,$ro)=split"\t";
  next if($ref eq $lastr && $rb < $lre && $re > $lrb && $ro eq $lro); # alt-tr
  if($ref ne $lastr) { if(resetruns()) { putc('|');} }
  ($lastr, $lrb, $lre, $lro)= ($ref, $rb, $re, $ro);
  
  $p= (m/paralog=([^;\s]+)/) ? $1 : "."; # ditto for tandy=
  #$p= (m/tandy=([^;\s]+)/) ? $1 : "."; # ditto for tandy=

  $pc{$p}++; push(@pw,$p);
  $n++; 
    
  $po= shift @pw; putmark($po);

  unless(grep { $_ eq $po } @pw) {
    delete $pc{$po}; 
    $cf= delete $c{$po}; 
    push(@cfree, $cf) if($cf and $cf ne ".");
  }
  
}
print "\n";
putmeans();