#!/usr/bin/perl

=item tandem algo3

  1. find interesting ++ marked exon/gene
  
  2. collect all exons of ++marked in region i - 100 .. i + 200 or such
      save as geneid->@exons; this should include exons of tandem dupls.
      
  3. add in same/almost same exon locations for more geneids to add common id
      redo step 2. to collect new geneid exons in region
      
  4. segregate exons to gene models in tandem region (how??)
      - using exon.nums from ids, break at change in exon.num sequence?  1,2,3,4,5 | 2,3,4,5
      - using max intron size?  max/ave of step b/n base exons?
  
  5. mark as done all gene/exons in region and move on (step 1)

=cut

sub print_genegroup_gff {
  my( $genes, $geneloc, $commonid )=@_;
  
  foreach my $geneid ($commonid, sort keys %$genes) {
    
    next unless exists $geneloc->{$geneid};
    my($gref,$gb,$ge,$gstrand, $geq, $galt, $gnexons, $gnaltexon)= @{ $geneloc->{$geneid} };
    
    my $matchID= "td_$geneid"; # FIXME: option
    my $gat= "ID=$matchID";
    
    (my $geqtype= $EQlegend{$geq}) =~ s/\W/_/g;
    if($commonid eq $geneid) {
      $geqtype= "common"; print "\n";
      }
    $gat .= ";tdclass=$geqtype"; 
    if($galt) { 
      $galt=~s/[\+\*\@]//g; # really should compress this to single hashed id and match w/ others
      my @galt= split ",",$galt;
      foreach (@galt) { s/\:[\d\-]+$//; s/\.[\d]+$//; } # drop .exnum:location
      my %galt= map{$_,1}@galt; $galt= join".",sort keys %galt; # put in std order
      $gat.= ";altid=$galt"; 
      }
    $gat.=";nexons=$gnexons" if($gnexons);
    $gat.=";altexons=$gnaltexon" if($gnaltexon);
    
    print join("\t",$gref,"tandy.".$geqtype,"match", $gb, $ge,".",$gstrand,".",$gat),"\n";
    delete $geneloc->{$geneid};
    
    next unless exists $genes->{$geneid};
    my($l_ref, $l_b, $l_e)=(0)x3;
    my @exons= @{ $genes->{$geneid} };
    foreach my $ex (sort _exlocsort @exons) {
      my($exnum, $xid, $ref, $b,$e, $p, $eq, $xalt)= @$ex;
      
      next if($ref eq $l_ref && $b eq $l_b && $e eq $l_e);
      ($l_ref, $l_b, $l_e)= ($ref,$b,$e);
      
      (my $xeqtype= $EQlegend{$eq}) =~ s/\W/_/g;
      my $xat= "Parent=$matchID;tdclass=$xeqtype;xid=$xid;ix=$exnum"; ##;eq=$eq
      if($galt) { #? or xalt 
        $xat.= ";altpar=$galt";  
        }
      print join("\t",$gref,"tandy.".$xeqtype,"match_part",$b,$e,$p,$gstrand,".",$xat),"\n";
      }   
    delete $genes->{$geneid};
  }    
}


sub print_genegroup {
  my($ggroup, $commonid, $exon_table, $atx)= @_;
  
  ## all of ggroup is presumptive common-gene-set
  ## need to 4. separate out gene models, add a common-id/tag, adjust exons?
  
  my %genes=();
  my %geneloc=();
  my ($nexons, $naltexons, @altids);
  foreach my $gg (@$ggroup) {
    my($ix, $geneid, $exons)= @ $gg;
    
    foreach my $exon (@$exons) {
      my($ref,$qid,$qb,$qe,$tb,$te,$strand,$eq,$p,$alt, @altexons)= @{ $exon };
      ## qid may have location: ID=Dappu1_FM5_96539.2:169647-169896 ; drop it
      my($trid,$exnum)= ($qid,1);
      if($qid =~ m/^(.+)\.(\d+)/) { $trid=$1; $exnum=$2; }
      
      $nexons++;
      $naltexons += scalar(@altexons) if(@altexons);
      push @altids, map { s/[:].*$//; $_; } split ",",$alt if($alt);
      
      foreach my $id ($commonid,$trid) {
      unless(exists $geneloc{$id}) {
        $geneloc{$id}=[$ref,$tb,$te,$strand, $eq, $alt]; 
        }
      else {
        $geneloc{$id}[1]= $tb if( $tb<$geneloc{$id}[1]);
        $geneloc{$id}[2]= $te if( $te>$geneloc{$id}[2]);
        }
      }

      #?? handle @altexons 
      foreach my $id ($trid) { # ?? , $trid
        unless(exists $genes{$id}) { $genes{$id}=[]; }
        push( @{ $genes{$id}}, [$exnum, $qid, $ref, $tb, $te, $p, $eq, $alt]);
        ## FIXME: ^ with common ids, can have many exons at same location ; test 
      }
      
    }
  }
  
  # my %altids= map {$_,1} @altids;
  $geneloc{$commonid}->[5]= join(",", @altids); ## (sort keys %altids)); #[5]== $alt 
  $geneloc{$commonid}->[6]= $nexons;
  $geneloc{$commonid}->[7]= $naltexons;
  
  print_genegroup_gff( \%genes, \%geneloc, $commonid);
}

sub _sameloc1 { # ref1,b1,e1 vs ref2,b2,e2; equal location with some +/- 10b slop
  my($ra,$ba,$ea, $rb, $bb, $eb)= @_;
  return ($ra eq $rb 
    && int($ba/SAMEBASE) == int($bb/SAMEBASE) 
    && int($ea/SAMEBASE) == int($eb/SAMEBASE));
}


sub tandem3 {
  my($exon_table, $genomefa, $queryfa, $location, $idmap)=@_;

  tandy_gffhead($exon_table, $genomefa, $queryfa, $location, $idmap);
  
  my (%genes, %geneloc, @genegroup, %bin_save, $atid);

  my $nx= scalar(@$exon_table);
  foreach my $ix (0..$nx-1) {
  
    my $thisexon=  $exon_table->[$ix];
    my($ref,$qid,$qb,$qe,$tb,$te,$strand,$eq,$p,$alt)= @{ $thisexon };
    # $thisexon->[7] == $eq;
    
    next unless($eq == -2 || $eq == 0  || $eq == 4); # ok? 
    next if($ref =~/_$/); #mark did this one
    next unless($alt =~ m/[^+]?(\+\+)([^,+]+)/);  # only stop at first mark
    #my ($pluss, $xgid)= ($1, $2);
    my $xgid= $1.$2;
    
    my $commonid="";
    @genegroup=(); # only collect over this marked gene? then print @genegroup
    my @galt=  split ",",$alt ;
    unshift(@galt, $xgid); # do 1st ??
    my %gidmark;
    my %gids= map{ 
      my $mk=(m/[\+]/)?1:0; s/[\+\*\@]//g;  s/\:[\d\-]+$//; 
      $gidmark{$_}++ if($mk); $_,1; 
      } @galt;
    
    #while($alt =~ m/[^+]?(\+\+)([^,+]+)/g) # only stop at first mark
    while (my $xid= shift(@galt))
    {
      
      ## ? dont do all ids here, just ++marked ones ?, but we probably want eq == 1,2 gene types 
      
      #my ($pluss, $xid)= ($1, $2);
      my $pluss = ($xid =~ m/[\+]/) ? 1 : 0; # $xid =~ tr/\+/\+/;
      $xid =~ s/[\+\*\@]//g;  
      $xid =~ s/\:[\d\-]+$//; # drop :location; keep .exnum
      my($geneid,$exnum)= ($xid,1);
      if($xid =~ m/^(.+)\.(\d+)/) { $geneid=$1; $exnum=$2; }
      
#       unless($pluss) { $pluss= $gidmark{$geneid}; }
#       unless($pluss) {
#         my @balt= keys %{$idmap{$geneid}};
#         foreach my $bg (@balt) { if($gidmark{$bg}) { $pluss=1; last; }}
#       }
#      next unless($pluss); #?
     
      unless($commonid) {
        $commonid= $idmap->{$geneid}{common} || $geneid; 
        $commonid= "c.".$commonid;
        }
        
      my @exons=(); my @altexons= ();
      
      # 2. forestep collect all geneid exons
      ## my($ix0, $ix1)= ($ix-150, $ix+150);
      ## my($ix0, $ix1)= ($ix-1500, $ix+1500);
      my($ix0, $ix1)= ($ix-500, $ix+500);
      $ix0= 0 if($ix0<0);  $ix1= $nx if($ix1>$nx);
      
      for (my $jx= $ix0; $jx<$ix1; $jx++) {
        my $exon = $exon_table->[$jx];
        
        my($jref,$jqid,$jqb,$jqe,$jb,$je,$jstrand,$jeq,$jp,$jalt)= @$exon;
        next if($jref =~/_$/); #mark did this one
        
        if ( ($jx == $ix) 
          || ($jqid =~ /$geneid\.(\d+)/ || $jalt =~ /$geneid\.(\d+)/ ) ){
          my $jexnum= ($jx == $ix) ? $exnum : $1;
          my $jxid= $geneid.".".$jexnum;
          my @jalt=split ",",$jalt;
          foreach (@jalt) { 
            my $mk=0; $mk=1 if(s/[\+]//g); s/[\*\@]//g;  s/\:[\d\-]+$//;
            $gidmark{$_}++ if($mk); 
            push(@galt, $_) unless($gids{$_}++); 
            } 
          $exon->[0].="_"; # mark done
          push( @exons, $exon);
          }
        }

      # 3. exons at same/almost same locs
if(0) {      
      for (my $jx= $ix0; $jx<$ix1; $jx++) {
        my $altexon = $exon_table->[$jx];
        
        my($jref,$jqid,$jqb,$jqe,$jb,$je,$jstrand,$jeq,$jp,$jalt)= @$altexon;
        next if($jref =~/_$/); #mark did this one
   
        foreach my $ex (@exons,@altexons) {
          (my $eref= $ex->[0]) =~ s/_$//;
          my $issame= _sameloc1($jref,$jb,$je, $eref, $ex->[4], $ex->[5]);
          if($issame) {
            push @altexons, $altexon;
            $altexon->[0].="_"; # mark done
            push @{$ex}, $altexon;
            last;
            }
          }
        }
}

      # $genex{$geneid}= \@exons;  # realy should combine all at this $ix marker
      push @genegroup, [$ix, $geneid, \@exons ] if(@exons);
      }
    
    print_genegroup( \@genegroup, $commonid, $exon_table, $ix) if @genegroup;
    }
    
}


1;