#!/usr/bin/perl -w

=head1 SYNOPSIS
  
  tandemgenes.pl - functions to locate tandem duplicate genes in genomes
    using seed set of predicted exons (GFF input), genome fasta

  d.gilbert, may 2007

=head1 run bunch

#!/bin/tcsh
# tandyrun.sh

source `/bio/argos/ROOT/bin/argos-env`
set em=/bio/bio-grid/mb/EVidenceModeler
set pa=/bio/bio-grid/mb/PASA/

setenv EVM_HOME $em
setenv PASA_HOME $pa
setenv PASAHOME $pa
set path= ( $path $PASA_HOME/bin )

set tandy=/bio/bio-grid/dpulex/prots/tandemgenes.pl
set data=/bio/bio-grid/mb/EVidenceModeler/daphd
cd $data

#  (scaffold_*)
foreach sc (scaffold_9*)
  echo "# tandemfind $sc"
  pushd $sc
  $tandy -debug -act extract -gff dpulex1_predict.gff -genome dpulex1.fa \
     jex//CDS:JGI gex//CDS:NCBI_GNO dex//CDS:DGIL_SNO > dpulex1_exons.fa  
  $tandy -debug -act reduce -fa dpulex1_exons.fa  > dpulex1_exons.nr
  $tandy -debug -act tandy -genome dpulex1.fa -query dpulex1_exons.nr \
     -mineval 1e-30  > dpulex1_exons.nr.gff
  popd
end

#? use exon or CDS feature? some only w/ CDS


=cut

use lib '/bio/argos/common/perl/lib','/bio/argos/common/system-local/perl/lib';#debug
use strict;
use Bio::DB::GFF;
use Bio::DB::Fasta;
use File::Basename;
use Getopt::Long;

#? patch_feature_extend();
use vars qw/%EQlegend/;

## NEED Parent in aggregator to distinguish separate matches by same prot
my @IDGROUPS = qw(Parent Target ID); # with this cant get target()->name !!!
my @aggregators = qw( alignment coding match processed_transcript ); #matchparts 

my $BINSIZE = 5000;
my $NEARDIST= 15000; #?

my $ADDLOC2ID=1;
my $MINEVAL= 1e-10;
my ($MINSIZE, $MAXSIZE,$REDUCE_GROUPONLY);
my ($igequal, $igaltequal, $igunequal, $ignear, $igclass, $iginsideof, $iginnear, $ignoreclass);
$iginsideof=1; # is $eq == 3 ; match is inside other known matches

my ($debug, $force, $fasta, $queryfa, $gff, $location, $action, 
    $fasta_db, $database, $user, $password, ) =('')x30;

my $optok= GetOptions(
  "database=s", \$database, "user=s",  \$user, "password=s", \$password,
  "gff=s", \$gff, "fasta|genome=s",\$fasta, "queryfasta=s", \$queryfa,
  "location=s",  \$location,
  "action=s",  \$action,
  "debug!",  \$debug, "force!", \$force,
  "ignoreclass=s", \$ignoreclass,
  "igequal|equal!",  \$igequal, "igaltequal|altequal!", \$igaltequal, 
    "igunequal|unequal!", \$igunequal, 
    "iginside|insideof!", \$iginsideof, "igclass", \$igclass,
    "ignear", \$ignear, "iginnear|innear!", \$iginnear,
  "MINEVAL=s", \$MINEVAL, "MINSIZE=n", \$MINSIZE, "MAXSIZE=n", \$MAXSIZE,
  "NEARDIST=s", \$NEARDIST, "GROUPONLY!", \$REDUCE_GROUPONLY,
  );

# print STDERR "# EQ_legend: ", join( ", ",  map{ $_."=".$EQlegend{$_}} sort{$a<=>$b}keys %EQlegend),"\n";

die <<USAGE unless($optok && ($database||$gff||$fasta) && $action);
usage: script [-options]   name1//feature-type/type2  name2//type2 name3/type3

   script -act extract -gff my.gff -fa my.fa exons//exon:JGI  > exons.fa
   script -act reduce -fa my.fa  > exons.nr
   script -act find -genome genome.fa -query exons.nr > exon-find.out
       script -act tandy -genome genome.fa -query exons.nr -mineval 1e-30 > exon-find.gff
       
   script -act overlap -gff my.gff jgigenes//mRNA:JGI  est//HSP:EST  dmgene//HSP:modDM
    
   -action = fasta-extract, reduce|nonredundant-fasta, find-tandem|tandy, overlap-count, ...

   -gff  file.gff|-   
   -fasta|genome file.fa [genome.fa] ; -queryfasta query.fa
   -location= scaffold:start-end for gff
   -database= mysql-dsn, lucene /path/to/index; for gff
      -user= , -pass= for mysql
 find act:
   -igequal, -igaltequal, -igunequal, -iginside, -ignear, -iginnear : ignore these exon classes
   -ignoreclass= dont compare to this ID class (^IDpattern)
   -minEval $MINEVAL  : skip low scores in find exons
   -nearDist $NEARDIST : distance to mark as near region unknown match (possible tandem)
 reduce act:   
   -grouponly : dont remove nr, mark group only 
 extract act:   
   -minsize n , -maxsize n  
   
 requirements:  blat (reduce, find),  Bioperl/DB/GFF,Fasta  
                exonerate (tandy final tandem validate?)
USAGE

my @typeset;
my %typeset = map {
  my ($name, $rest) = split("//", $_);
  push @typeset, $name;
  ($name => [ split("/", $rest) ])
} @ARGV;


if($action =~ /fasta|extract/) {
  my $seg = openGFFdb($location); 
  fasta_extract($seg);
  
} elsif($action =~ /reduce|nonredundant/) {
  blat_reduce($fasta);
  
} elsif($action =~ /find|tandy|tandx/) {
  my $flags=""; 
  $flags.="noprint" if($action =~ /tandy|tandx/);
  
  my $exon_table= exon_finder( $fasta, $queryfa, $location, $flags);
  ## find tandems strategy 1: find -2 near dupls, show all in region ?
  
  if($action =~ /tandy/ && @$exon_table > 0) {
    tandem_gff( $exon_table, $fasta, $queryfa, $location);
    
  } elsif($action =~ /tandx/ && @$exon_table > 0) {
    tandem_exonerate( $exon_table, $fasta, $queryfa, $location);
  }
  
} elsif($action =~ /overlap/) {
  my $seg = openGFFdb($location); 
  overlap_count($seg);
  
} elsif($action) {
  warn "# unknown action: $action\n";
}

exit;

#--------------------------------

sub openGFFdb
{
  my($location)=@_;
  #return unless($gff || $database);
  
      ## there is a damn bug in some version of BioPerl GFF (mem db) 
      ## where ID/Parent/disid are missing
      ## is this related to refseq range entry missing?
      
  # globals: database, gff, fasta, user/password, @aggregators, @IDGROUPS
  my $db = "";
  if($gff) {
    $db= Bio::DB::GFF->new( -adaptor => "memory",
		      -gff => $gff, -fasta => $fasta,
		      -aggregators => \@aggregators, -preferred_groups => \@IDGROUPS,  );
		      
  } elsif($database =~ /lucene/) {
    $db= Bio::DB::GFF->new( -adaptor => "lucene",
		      -dsn => $database, 
		      -aggregators => \@aggregators,  -preferred_groups => \@IDGROUPS, );
		      
  } elsif($database) {
    my $dsn = ($database =~ /dbi:/) ? $database : "dbi:mysql:$database";
    $db= Bio::DB::GFF->new( #-adaptor => "mysql",
          -dsn => $dsn,
		      -user => $user, -password => $password,
		      -aggregators => \@aggregators, -preferred_groups => \@IDGROUPS,  );
  } else {
    warn "# No GFF database\n"; return;
  }  
 
  my $seg= $db;
  if($location) {
    my($ref,$start,$stop)= split(/[:,\-\.]+/, $location);
    ($seg)= $db->segment(-ref=>$ref,-start=>$start,-stop=>$stop);
    warn "# Location restriction: $location\n";
  }

  return $seg;
}


sub get_dna {
  my($fasta, $ref, $start, $stop)= @_; #, $fasta_db_ref
  unless(-f $fasta && $ref && $stop>0) {
    warn "need fasta file, ref-ID, start, stop location\n"; return;
    }
  
  ## my $fasta_db= $$fasta_db_ref if(ref $fasta_db_ref); # ??  
  my $havedb= ref $fasta_db;
  if($havedb) {
    $havedb= $fasta_db->index_name() eq $fasta_db->index_name($fasta,-d $fasta);
    }
  unless($havedb) {
    my $db = eval {Bio::DB::Fasta->new($fasta);}
      or warn "$@\nCan't open sequence file(s). \n" and return;
    $fasta_db= $db;  ## $fasta_db_ref= \$fasta_db;
    }
    
  my $seq = $fasta_db->seq($ref, $start => $stop) 
      or warn "cant locate seq of $ref:$start-$stop\n" and return;
  $seq= $seq->seq if(ref $seq); # is this weird bioperl change here or not
  return $seq;
}


sub fasta_extract {
  my($seg)= @_;

  my($nex)=(0)x10;
  for my $reference (@typeset) {
    my %didid=(); my $typeref= $typeset{$reference};
    my @feat = $seg->features(-type => $typeset{$reference}); # == any feature type
    for my $feat (@feat) {
    
      my($ab,$ae)= ($feat->abs_start, $feat->abs_end); ($ab,$ae)=($ae,$ab) if($ab>$ae);
      my $size=$ae-$ab+1;
      next if((defined $MINSIZE and $size < $MINSIZE) 
           || (defined $MAXSIZE and $size > $MAXSIZE));
   
      my $seq= $feat->seq;
      next unless($seq);
      $seq= $seq->seq if(ref $seq); #bioperl version change
      $seq =~ s/(\w{60})/$1\n/g;
      
      my ($id, $ispar);
      ($id)= $feat->attributes("ID") and $ispar=0;
      ($id)= $feat->attributes("Parent") and $ispar=1 unless($id);
      ($id)= $feat->display_id and $ispar=0 unless($id);
      ($id)= $feat->id and $ispar=0 unless($id);
      ## there is a damn bug in some version of BioPerl GFF (mem db) 
      ## where ID/Parent/disid are missing
      #warn "#dbg ft=",$feat->gff3_string(0) if $debug;
 
      $didid{$id}++; 
      $id.=".".$didid{$id} if($ispar || $didid{$id}>1);
      
      #my $loc=  $feat->ref() .':'. $feat->abs_start .'-'. $feat->abs_end .':'. $feat->strand; 
      my $ref= $feat->ref();
      my $ori= $feat->strand();
      my $loc= "$ab-$ae";
      $id.=":$loc" if($ADDLOC2ID);

      my $doc= join("; ","type=".$feat->type(), "loc=$ref:$loc:$ori"); #$id,
      print ">$id $doc\n$seq\n";  $nex++;
    }
    warn "# fasta_extract n=$nex of $reference//",join(",", @$typeref),"\n"; $nex=0;
  }
  
  # warn "# fasta_extract $nex\n";
}



sub dosystem {
  warn "# ",join(" ",@_),"\n" if $debug;
  return system(@_) if(@_);
}


=item blat find alternate/tandem/lost exons

blat -oneOff=1 -minIdentity=80 scaffold_4.fa jgifm4_exons.fa jgifm4_exons.blatf8 -out=blast8

# list non-same exons (maybe known, unknown)
cat gnomon4_exons.blatf8 | perl -ne 'chomp; @v=split"\t"; ($q,$b,$e,$p,$s)=@v[0,8,9,10,11]; ($q,$qb,$qe)=sp
lit(/[:-]/,$q); \
 $eq=($qb==$b and $qe==$e);  \
 print join("\t",$q,$qb,$qe,$b,$e,$p),"\n" unless($eq or $p>1e-50);' | sort -k1,1 | more

# list located exons, same, known-alt and unknown
 cat gnomon4_exons.blatf8 jgifm4_exons.blatf8 | \
  perl -ne '@v=split"\t"; ($b,$e)=@v[8,9]; @v[8,9]=($e,$b) if($b>$e); \
  print join("\t",@v);' | sort -k9,9n -k10,10nr | \
  perl -ne 'BEGIN{print join("\t",qw(GeneExonID______ OrigB OrigE AtB AtE EQ Eval)),"\n";}\
  chomp; @v=split"\t"; ($q,$b,$e,$p,$s)=@v[0,8,9,10,11]; ($q,$qb,$qe)=split(/[:-]/,$q); \
 ($b,$e)=($e,$b) if($b>$e); ($qb,$qe)=($qe,$qb) if($qb>$qe); \
 $eq=($qb==$b and $qe==$e)?1:0;  \
 print join("\t",$q,$qb,$qe,$b,$e,$eq,$p),"\n" unless( $p>1e-50);' | more

# file format blatf8 == blast f8; genome locs at 8,9
NCBI_GNO_1182044.1:3009787-3009874      scaffold_4      100.00  88      0       0       1       88      
  3009787    3009874 1.1e-42 171.0
NCBI_GNO_1182044.2:3009959-3010051      scaffold_4      100.00  93      0       0       1       93      
  3009959    3010051 8.3e-46 181.0


=cut

sub _isoverlap { # ($gb,$ge, $qb,$qe, $qid) 
  if ($_[0] <= $_[3] && $_[1] >= $_[2]) # gb <= qe && ge >= qb
    { return ($ignoreclass && $_[4] =~ /^($ignoreclass)/) ? -1 : 1; }
  else { return 0; }
  }
sub _isnear {  # ($gb,$ge, $qb,$qe, $qid)
  my $gm= int(($_[0]+$_[1])/2); 
  my $qm= int(($_[2]+$_[3])/2); 
  return (abs($gm - $qm) < $NEARDIST) ? 1 : 0; 
  }




sub exon_finder {
  my($genomefa,$queryfa,$location,$flags)=@_;
  
  unless($genomefa && $queryfa && -f $genomefa && -f $queryfa) {
    warn "# Missing: genome $genomefa or query $queryfa\n"; return;
    }
    
  my $outfile= "$queryfa.blatf8";
#   (my $outfile = $queryfa) =~ s/\.\w+$//;
#   $outfile .= "_". basename($genomefa);
#   $outfile =~ s/\.\w+//;
#   $outfile .= ".blatf8";
  
  print "# exon_finder\n";
  my( %ab, @table, @allids, %altids, @groups, %ingroup);
  my( $doprint, $dosave)= (1,1); #??
  $flags ||="";
  $doprint=0 if $flags =~ /noprint/;
  $dosave=0 if $flags =~ /nosave/;
  
  warn "# grep >ids $queryfa\n" if $debug;
  open(GREP, "grep '^>' $queryfa|") 
    or die "grep $queryfa";
  while(<GREP>){ 
    my($id)=m/>(\S+)/; 
    push(@allids,$id) if($id); 
    if(m/altids=(\S+)/) { $altids{$id}= $1; $altids{$id}=~ s/pm,//; } # FIXME: pm, = partial match
    #if(m/altids=(\S+)/){ my @alt= split(/,/, $1); $altids{$id}= \@alt; }
    } 
  close(GREP);

  ## FIXME: option for megablast, parameter choices
  my $cmd="blat -oneOff=1 -minIdentity=80 $genomefa $queryfa -out=blast8 $outfile";
  if (!$force && -f $outfile) {
    warn "# reusing old $outfile;\n# $cmd\n";
  } else {
    my $ok= dosystem($cmd);
  }
  
  unless(-f $outfile) {
    warn "# No blat result found: $outfile\n"; return;
    }
    
  ## find tandems strategy 1: find -2 near dupls, show all in region ?
    
  ## FIXME: some of these eq==0 matches to 1mary exon are == alt exon with error
  ##  i.e. alt exon wasn't perfect match to 1st exon, show here as eq==0 to alt, w/ overlap
  ## fixed: add $ref/scaffold column to sort: -k2,2
  ## fixed: add >bitscore column to sort  : -k12,12nr
  
  my (%hits, %genehits);
  my $fixlocperl=  ##?? can we do w/o this ; losses strand info
  '@v=split"\t"; ($b,$e)=@v[8,9]; if($b>$e){ @v[8,9]=($e,$b); $v[0].="-";} print join("\t",@v);';
  #'@v=split"\t"; ($b,$e)=@v[8,9]; @v[8,9]=($e,$b) if($b>$e); print join("\t",@v);';
  my $sortblast8= 'sort -k2,2 -k9,9n -k10,10nr -k12,12nr'; # genome.ref,b,e ; bitscore
  
  my($locref,$locstart,$locstop)=(0)x10;
  if($location) {
    ($locref,$locstart,$locstop)= split(/[:,\-\.]+/, $location);
    print "# Location restriction: $location\n";
    }

  print "# EQ_legend: ", join( ", ",  map{ $_."=".$EQlegend{$_}} sort{$a<=>$b}keys %EQlegend),"\n"
      . "#".join("\t",qw(Genome___ GeneExonID______ OrigB OrigE AtB AtE Or EQ Eval AltIDs)),"\n"
      if($doprint);
      
  open(F, "cat $outfile | perl -ne '$fixlocperl' | $sortblast8 |")
    or die "processing $outfile";

  while(<F>){
    chomp; my @v=split"\t"; 
    my($qid,$ref,$b,$e,$p,$s)=@v[0,1,8,9,10,11]; 

    my $strand= ($qid =~ s/-$//) ? "-" : "+";
    my $wantloc=1;
    if($locstop>0) {
        ## be careful here; may want some below checks done before skipping location
      next unless($ref eq $locref);  # Ok always?
      $wantloc= _isoverlap($b,$e,$locstart,$locstop,"XXlocXX");  
      # next if($wantloc==0);  # ** defer to later; after _isnear() checks
    }
    
    my $alt= $altids{$qid}; # || ".";  
    my @alt= ($alt) ? split(",",$alt) : (); 
    $alt ||= ".";
    
    my($qb,$qe);
    ($qid,$qb,$qe)=split(/[:-]/,$qid); 
    ($qb,$qe)=($qe,$qb) if($qb>$qe);
    
    #####my $strand= ($b>$e) ? "-" : "+";
    ($b,$e)=($e,$b) if($b>$e);    ### already done in sort ** NEED for strand info ** ??????
    
    my $eq= _isoverlap($b,$e,$qb,$qe,$qid); # ($qb<=$e && $qe>=$b)? 1 : 0; 
    
    # check altids for eq >> known exon locs
    my ($altnear_id,$alteq_id)=('') x 9;
    my $altnear=0;
    foreach my $aid (@alt) {
      my($ab,$ae);
      ($aid,$ab,$ae)=split(/[:-]/,$aid); 
      ($ab,$ae)=($ae,$ab) if($ab>$ae);
      my $aeq= _isoverlap($b,$e,$ab,$ae,$aid);  
      if($aeq!=0) {
        $alteq_id=$aid;  $alt =~ s/$aid/\*$aid/; # mark all matches
        $eq=2 unless($eq>0 || $aeq<=0);
        # last;
        }
      else {
        if( _isnear($b,$e,$ab,$ae,$aid)) { $altnear= 1; $altnear_id=$aid; }
        }            
      }
    
    ## also, keep list of $eq, $alteq_id and screen out ! $eq inside these matches
    ## should be sorted to allow

    my $overlap=0; my $overlap_id;
    my $isnear=0;
       
    if ($eq>0) { 
      my $loc= [$qb,$qe,$qid,$eq]; ## need to add $alt locs from above ??
      my @bins= (int($qb/$BINSIZE) .. int($qe/$BINSIZE));
      foreach my $ib (@bins) { push @{$hits{$ib}}, $loc; } 
      
    } else {  #FIXME; overlap with a following Better exon not captured; sort by >bitscore?
      my $overlap_eq=1;
      my @bins= (int($b/$BINSIZE) .. int($e/$BINSIZE));
      BINS: foreach my $ib (@bins) { 
        if($hits{$ib}) { 
          foreach my $be (@{$hits{$ib}}) {
            my @be=@$be; ## {$hits{$ib}}; # fixme: list of @be
            my $heq= _isoverlap($b,$e,$be[0],$be[1],"XXignoreXX");    
            if($heq>0) {
              my $oid= $overlap_id= $be[2];  
              $overlap_eq= $be[3]; # if $oeq < 1, no overlap ?? or set -2 instead of 4
              $alt =~ s/$oid/\@$oid/; # mark all matches
              $overlap=1; 
              last BINS;
              }
            }
        }
      }
      
      # $eq=3 if $overlap; ## looks like we can drop these by default
      if($overlap) {
        $eq=3;
        unless($overlap_id eq $qid || $overlap_id eq $altnear_id) {
          $isnear= _isnear($b,$e,$qb,$qe,$qid);
          if($isnear||$altnear) { $eq= ($overlap_eq < 1) ? -2 : 4; }
          if($altnear) { $alt =~ s/$altnear_id/\*$altnear_id/;  }
          # FIXME: these inside+near shouldn't be printed separately,
          #   but mark geneid in overlapping main exon
          }
        }
    }

    ## ?? is this right??
    if($eq<=0 && @alt) { ## add to hits so we don't show all alt matches also
      my $loc= [$b,$e,$qid, $eq];  
      my @bins= (int($b/$BINSIZE) .. int($e/$BINSIZE));
      foreach my $ib (@bins) { push @{$hits{$ib}}, $loc; } # FIXME: should push(hits{ib},$loc)
    }
    
    ## if $eq == 0, want to check if near: need to look thru bins?
    ## if $eq == 3 (inside known), also check for _isnear and flag if these are known but tandems?
    
    if($eq<=0) {
      $isnear= ($eq<0) ? 0 : _isnear($b,$e,$qb,$qe,$qid);
      $eq= -2 if($isnear||$altnear); #??

#     } elsif ($eq == 3) { # do above
#       $isnear= _isnear($b,$e,$qb,$qe,$qid);
#       $eq= 4 if($isnear||$altnear); #??
#       # FIXME: these inside+near shouldn't be printed separately,
#       #   but mark geneid in overlapping main exon

    }
      
#   # not what we want; _isnear same as _isoverlap where eq == 1      
#     } elsif ($eq == 1) {
#       $isnear= _isnear($b,$e,$qb,$qe,$qid); 
#       $eq= 5 if($isnear||$altnear); #??
#     } elsif ($eq == 2) {
#       $isnear= _isnear($b,$e,$qb,$qe,$qid);
#       $eq= 6 if($isnear||$altnear); #??
    
    ## for $eq <= 0; mark regions of multi-exon matches (same geneid)
    ## try this for all $eq states? when _isnear
    if($eq<=0 || $eq>=4) { ## skip -1 here
      my $MINEXON_GENE=2; # or 3?
      my($geneid,$exnum)= $qid =~ m/^(.+)\.(\d+)$/;
      $exnum ||= 1; #??
      my @bins= (int($b/$BINSIZE) .. int($e/$BINSIZE));
      foreach my $ib (@bins) { 
        # need to clean this hash out over local regions;  like %hits @bins?
        $genehits{$geneid}{$ib}{$exnum}++; 
        my $nexons= scalar(keys %{$genehits{$geneid}{$ib}}); 
        if($nexons>=$MINEXON_GENE) {  ## really need to pre-mark exons 1,2 if get here .. already printed tho
          my $mark=('+') x $nexons;  
          unless($alt =~ s/\+*$geneid/$mark$geneid/){ $alt =~ s/$/,$mark$geneid/; }
          $alt=~ s/\.,//;
          }
        } 
    }
    
    next if($wantloc==0);  # defered skip from above; need _isnear() check before
  

BEGIN { # define here so we can update easily     
  %EQlegend= (
    0,"unknown", 1,"equal", 2,"altequal", 3,"insideof", 
    4,"near+overlap", # 5,"equal+near", 6,"altequal+near",  ##9,"error",
    -1,"ignoredmatch", -2, "near");
}; 
    
    unless( $p > $MINEVAL
      || ($eq == 0 && $igunequal) 
      || ($eq == 1 && $igequal)
      || ($eq == 2 && $igaltequal)
      || ($eq == 3 && $iginsideof)
      || ($eq == 4 && $iginnear)
      || ($eq == -2 && $ignear)
      || ($eq == -1 && $igclass)
      ) {   
      my $pfmt= sprintf("%8.1e",$p); #? same as before '7.8e-63' '2.8e-106'
      push(@table, [$ref,$qid,$qb,$qe,$b,$e,$strand,$eq,$p,$alt]) if($dosave);
      print join("\t",$ref,$qid,$qb,$qe,$b,$e,$strand,$eq,$pfmt,$alt),"\n" if($doprint);
      }

    
    }
  close(F);
  return \@table;
}


sub printgenes_gff {
  my($genes, $geneloc, @lastid)=@_;
  my %skipid= map { $_,1 } @lastid; # may be partial
  
  foreach my $geneid (sort keys %$genes) {
    next if($skipid{$geneid});
    my($gref,$gb,$ge,$gstrand, $geq, $galt)= @{ $geneloc->{$geneid} };
    
    my $matchID= "td_$geneid"; # FIXME: option
    (my $geqtype= $EQlegend{$geq}) =~ s/\W/_/g;
    my $gat= "ID=$matchID"; 
    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"; 
      }
    print join("\t",$gref,"tandy.".$geqtype,"match", $gb, $ge,".",$gstrand,".",$gat),"\n";

    my @exons= @{ $genes->{$geneid} };
    foreach my $ex (@exons) {
      my($exnum, $xid, $ref, $b,$e, $p, $eq, $xalt)= @$ex;
      (my $xeqtype= $EQlegend{$eq}) =~ s/\W/_/g;
      my $xat= "Parent=$matchID;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";
      }
      
    ## clear hash
    delete $genes->{$geneid}; ## unless( grep {$_ eq $geneid} @lastid );
    delete $geneloc->{$geneid}; # this too
  }
    
}



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

  #my $gffout= "$queryfa.find.gff"; 
  ##warn "# tandem_gff:  $gffout\n"; # if $debug

  print "##gff-version 3\n";
  print "# tandemgenes-version 1\n";
  print "# finding unlocated genes in genome from exon-collection\n";
  print "# query: $queryfa; genome: $genomefa\n";
  print "# Location restriction: $location\n" if($location);
  print "# source-legend: ", join( ", ",  
    map{ my $s= $EQlegend{$_}; $s =~ s/\W/_/g; "tandy.".$s; } 
      sort{$a<=>$b} keys %EQlegend),"\n\n";

  
  my (%genes, %geneloc, $atid, $lastid, $lastid2, $lastid3);
  my $nx= scalar(@$exon_table);
  foreach my $ix (0..$nx-1) {
  
    my($ref,$qid,$qb,$qe,$b,$e,$strand,$eq,$p,$alt)= @{ $exon_table->[$ix] };
    
    next unless($eq == -2 || $eq == 0  || $eq == 4); 

    ####(my $eqtype= $EQlegend{$eq}) =~ s/\W/_/g;
    
    if($alt =~ m/(\+{2,})([^,+]+)/) {
      my ($pluss, $xid)= ($1, $2);
      
      ## xid may have location: ID=Dappu1_FM5_96539.2:169647-169896 ;; drop it
      my($geneid,$exnum)= ($xid,1);
      if($xid =~ m/^(.+)\.(\d+)/) { $geneid=$1; $exnum=$2; }
      
      unless(exists $genes{$geneid}) { $genes{$geneid}=[]; }
      push( @{ $genes{$geneid}}, [$exnum, $xid, $ref,$b,$e, $p, $eq, $alt]);
      unless(exists $geneloc{$geneid}) { $geneloc{$geneid}=[$ref,$b,$e,$strand,$eq, $alt]; }
      else {
        $geneloc{$geneid}[1]= $b if( $b<$geneloc{$geneid}[1]);
        $geneloc{$geneid}[2]= $e if( $e>$geneloc{$geneid}[2]);
        }

      # backtrack in exontab for other exons: $ix-1,-2,-3...
      if($pluss eq '++') {
        my $need=2;
        for (my $jx= $ix-1; $jx>$ix-30; $jx--) {
          #my
          my($jref,$jqid,$jqb,$jqe,$jb,$je,$jstrand,$jeq,$jp,$jalt)= @{$exon_table->[$jx]};
          if ( ($jeq<=0 || $jeq>=4)  
            && ($jqid =~ /^$geneid\.(\d+)/ || $jalt =~ /^$geneid\.(\d+)/ )
             ){
            my($jexnum, $jxid);
            $jexnum= $1; $jxid= $geneid.".".$jexnum;
            push( @{ $genes{$geneid}}, [$jexnum,$jxid,$jref,$jb,$je, $jp, $jeq, $jalt]);
            $geneloc{$geneid}[1]= $jb if( $jb<$geneloc{$geneid}[1]);
            $geneloc{$geneid}[2]= $je if( $je>$geneloc{$geneid}[2]);
            --$need; last if($need<1);
            }
          }
        }
 
 
      $lastid3= $lastid2; $lastid2= $lastid; $lastid= $geneid;
      printgenes_gff(\%genes, \%geneloc, $lastid, $lastid2, $lastid3) 
        if($lastid3 && $lastid ne $lastid2 && $lastid2 ne $lastid3);
      ## think this is bad; storing only by geneid in hash causes mutli-loc matches
      #   to be joined in one geneloc .. bad idea; delete geneloc{geneid} cure?
      
      }
    }

  printgenes_gff(\%genes, \%geneloc);
}

=item tandem_exonerate

  for input table like this, create 'cdna' and use exonerate to locate fully on genome?

  #table =
  # Genome___      GeneExonID______        OrigB   OrigE   AtB     AtE     EQ      Eval    AltIDs
  scaffold_4      Dappu1_FM5_234437.8     663878  664019  651718  651855  -2      5.0e-53 .        ** putative rev-tandem; 6 of 8 exon matches
  scaffold_4      Dappu1_FM5_234437.4     662308  662457  653140  653225  -2      3.0e-38 .
  scaffold_4      Dappu1_FM5_234437.3     661775  662013  653520  653745  -2      2.3e-103        .
  scaffold_4      Dappu1_FM5_234437.2     661649  661752  653776  653879  -2      4.1e-51 .
  scaffold_4      Dappu1_FM5_234437.1     661211  661538  653988  654116  -2      1.5e-58 .
  scaffold_4      Dappu1_FM5_234437.1     661211  661538  654118  654319  -2      3.6e-85 .

=cut

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

  my $regionfile= "$queryfa.generegion"; #??
  my $cdnafile= "$queryfa.cdnas";
  warn "# tandem_exonerate: making inputs $regionfile, $cdnafile\n"; # if $debug
  
  my (%genes,%geneloc);
  foreach my $ex (@$exon_table) {
    my($ref,$qid,$qb,$qe,$b,$e,$strand,$eq,$p,$alt)= @$ex;
    
    #next unless($eq == -2 || $eq == 0); #?? screen here or not? do we want  $eq == 0 also?
    next unless($eq <= 0); #?? screen here or not? do we want  $eq == 0 also?
    
    # arg; putative gene may have exons from 2+ geneid's 
    # need to separate by contiguous location .. 
      ## xid may have location: ID=Dappu1_FM5_96539.2:169647-169896 ;; drop it
    my($geneid,$exnum)= $qid =~ m/^(.+)\.(\d+)/;
    $exnum ||= 1; #??
    
    ## arg; need to worry about exon order, etc.  not just random set of exons/gene
    my $exdna= get_dna( $genomefa, $ref, $qb, $qe);
    push(@{ $genes{$geneid}}, [$exnum, $exdna,$qid,$ref,$b,$e]);
    
    unless(exists $geneloc{$geneid}) { $geneloc{$geneid}=[$ref,$b,$e]; }
    else {
      $geneloc{$geneid}[1]= $b if( $b<$geneloc{$geneid}[1]);
      $geneloc{$geneid}[2]= $e if( $e>$geneloc{$geneid}[2]);
      }
  }

  open(GENOME,">$regionfile") or die "> $regionfile";
  open(CDNA,">$cdnafile") or die "> $cdnafile";
  foreach my $geneid (sort keys %genes) {
    my($gref,$gb,$ge)= @{ $geneloc{$geneid} };
    ($gb,$ge)= ($gb-1000, $ge+1000);  # expand some .. $NEARDIST ??
    
    my $generegion= get_dna($genomefa, $gref, $gb, $ge);
    print GENOME ">${geneid}_region loc=$gref:$gb-$ge\n$generegion\n";
    
    my @cdna=(); my $chead=""; my $cdna=""; my $exmiss="nnnnnn";
    my @exons= @{ $genes{$geneid} };
    foreach my $ex (@exons) {
      my($exnum,$exdna,$qid,$ref,$b,$e)= @$ex;
      $chead .= "loc.$exnum=$ref:$b-$e";
      $cdna .= $exdna ."\n";
      
      ## sample exnum is 8,4,3,2,1,1, ... 8
      ## this is wrong ; input order is sorted right, but spacers are unknown
      # while( $cdna[$exnum-1] ) { $exnum++; }
      # $cdna[$exnum-1]= $exdna;  
      # print CDNA ">${geneid}_exon.$exnum loc=$ref:$b-$e\n$exdna\n"; 
        #?? want to try this way? ; exonerate wont stitch these into one alignment tho.
    }
    # do we write separate EST-like .fa instead of one (partial) cDNA ??
    # foreach (@cdna) { $_= $exmiss unless($_); }
    # my $cdna= join("",@cdna);
    print CDNA ">${geneid}_cdna\n$cdna\n";
  }
  close(GENOME); close(CDNA);
  
  warn "# tandem_exonerate: 
  exonerate --bestn 1 --model cdna2genome $cdnafile $regionfile
  \n";
}


=item blat reduce redundant exons

 to turn all exon fasta set into non-redundant set
  blat -fastMap scaffold_4.fa  gnomon4_exons.fa gnomon4_exons.blatf8 -out=psl

 for exon map w/ error:
  blat -oneOff=1 -minIdentity=80 scaffold_4.fa  gnomon4_exons.fa gnomon4_exons.blatf8 
     ?? -out=blast8 or psl

blat -fastMap all_exons.fa all_exons.fa stdout | perl -ne
'@v=split"\t";($a,$b)= sort @v[9,13]; $a{$a}++; $ab{$a}{$b}++; $ab{$b}{$a}++; 
END{ @a= reverse sort keys %ab; foreach $a (@a) { 
  @bb= sort keys %{$ab{$a}}; 
  @b=(); foreach (@bb) { push(@b,$_) if($ab{$a}{$_});} 
  print join(" ",@b),"\n" if(@b); 
  foreach $b (@b) { $ab{$a}{$b}=0; $ab{$b}{$a}=0; } } 
}' > ! common_exons.idlist2

  # blat table cols (psl)
  
  match  mismat  xxx  nnn  qg qg tg tg ori  Qname Qsize Qb Qe  Tname Tsize Tb Te blkn  blkxxx
    0     1       2    3    4  5 6  7   8    9    10    11 12   13    14   15 16  17   18..
    
    
=cut

sub blat_reduce {
  my($infasta)=@_;
  my($outfasta)= "$infasta.nr";
  my(%ab, %abpart, @allids, @pgroups, %inpgroup, @groups, %ingroup);
  
  warn "# grep >ids $infasta\n" if $debug;
  open(GREP, "grep '^>' $infasta|") 
    or die "grep $infasta";
  while(<GREP>){ my($id)=m/>(\S+)/; push(@allids,$id) if($id); } 
  close(GREP);

my $RMIS=2;

  warn "# blat -fastMap $infasta  $infasta -out=psl stdout\n" if $debug;
  open(BLAT, "blat -fastMap $infasta  $infasta -out=psl stdout |")
    or die "blat $infasta";
  # psLayout version 3
  while(<BLAT>){
    my @v=split"\t"; next unless(/^\d+/ && @v > 13);
    my($a,$b)= sort @v[9,13]; # ids
    my($mat, $mis, $at,$ab,$ae, $bt,$bb,$be)= @v[0, 1, 10,11,12, 14,15,16]; # sizes
    
    #$abmat{$a}{$b}=$mat; $ablen{$a}=$at; $ablen{$b}= $bt;
    ## need to handle case where $at <<< $bt or vs; short vs long
    if($at >= $bt-$RMIS && $bt >= $at-$RMIS && $mat >= $at-$RMIS) { #? should use some %align, like 90%
       $ab{$a}{$b}= $mat;
       $ab{$b}{$a}= $mat;
       
    } else { ## keep all others ; if ($at >= $bt-$RMIS && $mat >= $bt-$RMIS) { 
      # here $b is ~full in $a but may be <<shorter
       $abpart{$a}{$b}= $mat;
    }
    
    ## FIXME: really should prefer longest exon if a!=b diff length
    ## also want reciprocal score test?
    # $ab{$a}{$b}++;  $ab{$b}{$a}++; #??
    }
  close(BLAT);
  my $nab= scalar(keys %ab);
  warn "# done blat -fastMap $infasta; nab= $nab\n" if $debug;

## FIXME: modify this; keep sizes, align-size and reduce only when truly same exon
## but store altid for all close matches

  my $ngroup=0;
  foreach my $a (sort keys %ab) { 
    my @b=(); 
    foreach my $b (sort keys %{$ab{$a}}) { 
      #push(@b,$b) if( defined $ab{$a}{$b} && defined $ab{$b}{$a}); 
      #delete $ab{$a}{$b}; delete $ab{$b}{$a}; # = 0; 
      push(@b,$b) if( $ab{$a}{$b}>0 && $ab{$b}{$a}>0 ); 
      $ab{$a}{$b}= $ab{$b}{$a} = 0;
      } 
    if(@b>1) {
      push(@groups, \@b);
      $ngroup= @groups; 
      @ingroup{@b}= ($ngroup) x scalar(@b);
      }
    }
  my $npgroup=0;
  foreach my $a (sort keys %abpart) { 
    my @b=(); 
    foreach my $b (sort keys %{$abpart{$a}}) { 
      push(@b,$b) if( $abpart{$a}{$b} >0 ); 
      $abpart{$a}{$b}= 0; 
      } 
    if(@b>1) {
      push(@pgroups, \@b);
      $npgroup= @pgroups; 
      @inpgroup{@b}= ($npgroup) x scalar(@b);
      }
    }
  
  # my %singles= map{ $_ => 1 } grep {!exists $ingroup{$_}} @allids;
  my %singles= map{ $_ => 1 } grep {! $ingroup{$_}} @allids;
  my $ninput   = scalar(@allids);
  my $ngrouped = scalar (keys %ingroup);
  my $nsingle  = scalar (keys %singles);
  my ($nprint, $print)=(0)x10;
  my $pgrouped = scalar (keys %inpgroup);

  warn "# reduce $infasta n=$ninput to $ngroup [w/$ngrouped, part=$pgrouped] + $nsingle [x1] \n" if $debug;
  open(FA,"$infasta"); 
  # open(OUT,">$outfasta") or die ">$outfasta"; # ? write to stdout instead
  while(<FA>){
    if(/^>(\S+)/){ 
      my $id=$1;
      if(my $ig= $ingroup{$id}) { 
        my $h= join ",", @{$groups[$ig-1]};
        $print=0; 
        if($h) { 
          s/$/; altids=$h/ ; $print=1; 
          }
        unless($REDUCE_GROUPONLY) { $groups[$ig-1]= []; } # done group
      } else {
        $print= exists $singles{$id}; # should be always true here
      }

        ## add in $inpgroup ; $pgroups[] here for partial/subset matches
      if(my $ig= $inpgroup{$id}) { 
        my $h= join ",", @{$pgroups[$ig-1]};
        if($h) {
          unless(s/(altids=[^\s;]+)/$1,pm,$h/) { s/$/; altids=pm,$h/ ; }  
        }
      }
      
      $nprint++ if $print;
      }
    print $_ if($print);
    }
  # close(OUT); 
  close(FA);
  
  warn "# blat_reduced $infasta\n";
  warn "# in ids=$ninput; out ids=$nprint; grouped=$ngrouped [$ngroup]; partial=$pgrouped; single=$nsingle\n";
  
  # return ($outfasta);
}



sub overlap_count {
  my($seg)= @_;

  print join("\t", qw(ref type gene_id comp overlaps equals distinct)),"\n";
  
  for my $reference (@typeset) {
    my @feat = $seg->features(-type => $typeset{$reference});
    for my $feat (@feat) {
      for my $typeset (@typeset) {
        ##next if ($typeset eq $reference); # dgg
        my @overlaps = $feat->overlapping_features(-type => $typeset{$typeset});
        my @equals; # dgg
        my @distinct;
        if (@overlaps) {
          @distinct = $overlaps[0];
          for (my $i = 1 ; $i < @overlaps ; $i++) {
            my $b = $overlaps[$i];
            ## print "b=",$b->display_id,",",$b->type,"\n"; #DEBUG
            my $overlapped = 0;
            for my $a (@distinct) { #?? should this be exons of $feat
            ### for my $a (map{$_->sub_SeqFeature()}@distinct) { #?? exons of $feat
            ### print "a=",$a->display_id,",",$a->type,"\n"; #DEBUG
              if ($a->equals($b)) {
                push @equals, $b;
                $overlapped = 1;
                last;
              } elsif ($a->overlaps($b)) {
                $a->extend($b); # dgg what is this? sub union()? overlap_extent ?
                $overlapped = 1;
                last;
              }
            }
            unless ($overlapped) {
              push @distinct, $b;
            }
          }
        }
        
        my ($gname)= ($feat->attributes("Name"),$feat->display_id); 
        my ($gtype)= $feat->type(); 
        print join("\t", $reference, $gtype, $gname, #was: id, 
          $typeset, scalar(@overlaps),  scalar(@equals), scalar(@distinct)), "\n";
      }
    }
  }

}


sub patch_feature_extend {
  local $^W = 0;
  eval <<'END';
use Bio::DB::GFF::RelSegment;
sub Bio::DB::GFF::RelSegment::extend {
  my $self     = shift;
  my($other,$so) = @_;
  return unless $other->isa('Bio::DB::GFF::RelSegment');
  return if $self->abs_ref ne $other->abs_ref;
  my ($low,$high);
  foreach ($self,$other) {
  $low  = $_->abs_start  if !defined($low)  or $low  > $_->abs_start;
  $high = $_->abs_end if !defined($high) or $high < $_->abs_end;
  }
  $self->{start}= $self->_abs2rel($low);
  $self->{stop}= $self->_abs2rel($high);
}
END
  warn $@ if $@;
}

__END__

=item   Daphnia Hemoglobin region, using JGI FM5 exons only  

../tandemgenes.pl -debug -act find -genome scaffold_4.fa -query jgifm4_exons.fa -MINEVAL 1e-30
# grep >ids jgifm4_exons.fa
# echo reusing old jgifm4_exons.fa.blatf8;
echo blat -oneOff=1 -minIdentity=80 scaffold_4.fa jgifm4_exons.fa -out=blast8 jgifm4_exons.fa.blatf8;
reusing old jgifm4_exons.fa.blatf8
blat -oneOff=1 -minIdentity=80 scaffold_4.fa jgifm4_exons.fa -out=blast8 jgifm4_exons.fa.blatf8

# exon_finder
# EQ_legend: -2=near, -1=ignoredmatch, 0=unequal, 1=equal, 2=altequal, 3=insideof
#Genome___      GeneExonID______        OrigB   OrigE   AtB     AtE     EQ      Eval    AltIDs
scaffold_4      Dappu1_FM5_221075.1     765     884     765     884     1       1.6e-60 .
scaffold_4      Dappu1_FM5_221075.2     941     1027    941     1027    1       2.2e-42 .
scaffold_4      Dappu1_FM5_221075.7     1585    1667    1585    1667    1       1.8e-39 .
...

Hb1>scaffold_4      Dappu1_FM5_96311.1      2366641 2366746 2366641 2366746 1       1.9e-53 .   
Hb1>scaffold_4      Dappu1_FM5_221319.2     2370276 2370436 2366863 2367023 -2      1.6e-74 .
Hb1>scaffold_4      Dappu1_FM5_96311.2      2366863 2367023 2366863 2367023 1       6.1e-87 .
Hb1>scaffold_4      Dappu1_FM5_96311.3      2367103 2367328 2367103 2367328 1       4.2e-125        .
Hb1>scaffold_4      Dappu1_FM5_207029.4     2376738 2376862 2367424 2367548 -2      7.9e-36 .
Hb1>scaffold_4      Dappu1_FM5_234835.4     2373419 2373543 2367424 2367548 -2      5.1e-35 .
Hb1>scaffold_4      Dappu1_FM5_96311.4      2367424 2367548 2367424 2367548 1       2.5e-65 .
Hb1>scaffold_4      Dappu1_FM5_96311.5      2367630 2367752 2367630 2367752 1       4.2e-64 .
Hb1>scaffold_4      Dappu1_FM5_221319.6     2371238 2371466 2367831 2368059 -2      2.2e-105        .
Hb1>scaffold_4      Dappu1_FM5_96311.6      2367831 2368059 2367831 2368059 1       9.8e-127        .
Hb1>scaffold_4      Dappu1_FM5_96311.7      2368133 2368355 2368133 2368355 1       6.4e-122        .

Hb2>scaffold_4      Dappu1_FM5_221319.1     2370075 2370184 2370075 2370184 1       5.8e-56 .
Hb2>scaffold_4      Dappu1_FM5_221319.2     2370276 2370436 2370276 2370436 1       7.0e-87 .
Hb2>scaffold_4      Dappu1_FM5_221319.3     2370514 2370739 2370514 2370739 1       9.0e-125        .
Hb2>scaffold_4      Dappu1_FM5_221319.4     2370819 2370943 2370819 2370943 1       3.2e-65 .
Hb2>scaffold_4      Dappu1_FM5_221319.5     2371045 2371167 2371045 2371167 1       1.2e-63 .
Hb2>scaffold_4      Dappu1_FM5_221319.6     2371238 2371466 2371238 2371466 1       1.1e-126        .
Hb2>scaffold_4      Dappu1_FM5_221319.7     2371535 2371701 2371535 2371701 1       1.8e-89 .

Hb3>scaffold_4      Dappu1_FM5_234835.2     2372901 2373058 2372901 2373058 1       5.6e-85 .
Hb3>scaffold_4      Dappu1_FM5_207029.3     2376433 2376658 2373121 2373346 -2      3.8e-122        .
Hb3>scaffold_4      Dappu1_FM5_234835.3     2373121 2373346 2373121 2373346 1       2.8e-125        .
Hb3>scaffold_4      Dappu1_FM5_207029.4     2376738 2376862 2373419 2373543 -2      1.1e-64 .
Hb3>scaffold_4      Dappu1_FM5_234835.4     2373419 2373543 2373419 2373543 1       1.9e-65 .
Hb3>scaffold_4      Dappu1_FM5_234835.5     2373673 2373795 2373673 2373795 1       1.3e-63 .
Hb3>scaffold_4      Dappu1_FM5_234835.6     2373868 2374093 2373868 2374093 1       3.2e-125        .
Hb3>scaffold_4      Dappu1_FM5_207029.7     2382097 2382308 2374171 2374288 -2      3.7e-41 .
Hb3>scaffold_4      Dappu1_FM5_234835.7     2374171 2374287 2374171 2374287 1       6.8e-60 .

Hb4,Hb5>scaffold_4      Dappu1_FM5_207029.1     2376038 2376140 2376038 2376140 1       5.6e-52 .
        scaffold_4      Dappu1_FM5_234835.2     2372901 2373058 2376211 2376368 -2      8.7e-79 .
        scaffold_4      Dappu1_FM5_42627.2      2386245 2386402 2376211 2376340 0       1.6e-68 .
Hb4,Hb5>scaffold_4      Dappu1_FM5_207029.2     2376214 2376368 2376214 2376368 1       3.1e-83 .
Hb4,Hb5>scaffold_4      Dappu1_FM5_207029.3     2376433 2376658 2376433 2376658 1       2.8e-125        .
Hb4,Hb5>scaffold_4      Dappu1_FM5_207029.4     2376738 2376862 2376738 2376862 1       1.7e-65 .
Hb4,Hb5>scaffold_4      Dappu1_FM5_207029.5     2376975 2377097 2376975 2377097 1       1.2e-63 .
        scaffold_4      Dappu1_FM5_234835.6     2373868 2374093 2377165 2377390 -2      2.8e-100        .
        scaffold_4      Dappu1_FM5_42627.6      2387211 2387436 2377165 2377390 0       4.1e-94 .
Hb4,Hb5>scaffold_4      Dappu1_FM5_207029.6     2381813 2382038 2377165 2377389 -2      8.4e-117        .
        scaffold_4      Dappu1_FM5_96311.6      2367831 2368059 2377165 2377357 -2      5.5e-76 .
        scaffold_4      Dappu1_FM5_221319.6     2371238 2371466 2377165 2377281 -2      1.6e-47 .
        scaffold_4      Dappu1_FM5_230334.6     2389863 2390088 2377259 2377345 0       7.9e-31 .
        scaffold_4      Dappu1_FM5_234835.7     2374171 2374287 2377451 2377558 -2      4.9e-35 .
Hb4,Hb5>scaffold_4      Dappu1_FM5_207029.7     2382097 2382308 2377451 2377552 -2      1.2e-50 .
Hb4,Hb5>scaffold_4      Dappu1_FM5_207029.1     2376038 2376140 2380722 2380824 -2      3.2e-51 .
        scaffold_4      Dappu1_FM5_42627.2      2386245 2386402 2380895 2381052 -2      1.8e-83 .
        scaffold_4      Dappu1_FM5_234835.2     2372901 2373058 2380895 2381046 -2      1.9e-71 .
Hb4,Hb5>scaffold_4      Dappu1_FM5_207029.2     2376214 2376368 2380898 2381046 -2      6.9e-76 .
        scaffold_4      Dappu1_FM5_230334.2     2388940 2389097 2380927 2381052 -2      2.8e-43 .
Hb4,Hb5>scaffold_4      Dappu1_FM5_207029.3     2376433 2376658 2381117 2381323 -2      4.5e-102        .
        scaffold_4      Dappu1_FM5_234835.3     2373121 2373346 2381117 2381323 -2      6.0e-99 .
        scaffold_4      Dappu1_FM5_42627.3      2386466 2386691 2381117 2381314 -2      4.5e-100        .
        scaffold_4      Dappu1_FM5_230334.3     2389174 2389399 2381127 2381299 -2      2.6e-72 .
        scaffold_4      Dappu1_FM5_96311.3      2367103 2367328 2381131 2381342 0       3.4e-90 .
        scaffold_4      Dappu1_FM5_221319.3     2370514 2370739 2381167 2381321 0       1.3e-69 .
Hb4,Hb5>scaffold_4      Dappu1_FM5_207029.4     2376738 2376862 2381416 2381540 -2      1.3e-62 .
        scaffold_4      Dappu1_FM5_234835.4     2373419 2373543 2381416 2381540 -2      8.3e-62 .
        scaffold_4      Dappu1_FM5_42627.4      2386765 2386889 2381425 2381540 -2      1.8e-59 .
        scaffold_4      Dappu1_FM5_42627.5      2387018 2387140 2381630 2381745 -2      9.8e-40 .
        scaffold_4      Dappu1_FM5_234835.5     2373673 2373795 2381630 2381721 -2      1.2e-41 .
Hb4,Hb5>scaffold_4      Dappu1_FM5_207029.5     2376975 2377097 2381633 2381745 -2      1.0e-52 .
Hb4,Hb5>scaffold_4      Dappu1_FM5_207029.6     2381813 2382038 2381813 2382038 1       3.7e-125        .
Hb4,Hb5>scaffold_4      Dappu1_FM5_207029.7     2382097 2382308 2382097 2382308 1       4.4e-115        .

Hb6,Hb7>scaffold_4      Dappu1_FM5_42627.1      2383418 2383492 2383418 2383492 1       1.2e-35 .
        scaffold_4      Dappu1_FM5_230334.3     2389174 2389399 2383802 2383994 -2      9.8e-87 .
        scaffold_4      Dappu1_FM5_234835.3     2373121 2373346 2383811 2383973 0       7.0e-64 .
        scaffold_4      Dappu1_FM5_207029.3     2376433 2376658 2383813 2383973 -2      1.6e-63 .
Hb6,Hb7>scaffold_4      Dappu1_FM5_42627.3      2386466 2386691 2383813 2383973 -2      9.8e-62 .
        scaffold_4      Dappu1_FM5_221319.3     2370514 2370739 2383871 2383976 0       3.4e-44 .
        scaffold_4      Dappu1_FM5_96311.3      2367103 2367328 2383871 2383952 0       2.6e-34 .
        scaffold_4      Dappu1_FM5_207029.1     2376038 2376140 2386072 2386174 0       2.0e-50 .
Hb6,Hb7>scaffold_4      Dappu1_FM5_42627.2      2386245 2386402 2386245 2386402 1       4.4e-85 .
        scaffold_4      Dappu1_FM5_207029.3     2376433 2376658 2386466 2386691 0       4.0e-121        .
        scaffold_4      Dappu1_FM5_234835.3     2373121 2373346 2386466 2386691 0       5.3e-118        .
Hb6,Hb7>scaffold_4      Dappu1_FM5_42627.3      2386466 2386691 2386466 2386691 1       2.5e-125        .
Hb6,Hb7>scaffold_4      Dappu1_FM5_42627.4      2386765 2386889 2386765 2386889 1       1.5e-65 .
Hb6,Hb7>scaffold_4      Dappu1_FM5_42627.5      2387018 2387140 2387018 2387140 1       6.2e-64 .
        scaffold_4      Dappu1_FM5_234835.6     2373868 2374093 2387211 2387436 0       5.5e-107        .
Hb6,Hb7>scaffold_4      Dappu1_FM5_42627.6      2387211 2387436 2387211 2387436 1       4.2e-125        .
        scaffold_4      Dappu1_FM5_42627.7      2387508 2387624 2387508 2387624 1       4.1e-60 .

Hb8>scaffold_4      Dappu1_FM5_230334.1     2388721 2388828 2388721 2388828 1       4.6e-55 .
Hb8>scaffold_4      Dappu1_FM5_230334.2     2388940 2389097 2388940 2389097 1       8.3e-85 .
Hb8>scaffold_4      Dappu1_FM5_230334.3     2389174 2389399 2389174 2389399 1       1.0e-124        .
Hb8>scaffold_4      Dappu1_FM5_230334.4     2389476 2389600 2389476 2389600 1       5.3e-65 .
Hb8>scaffold_4      Dappu1_FM5_230334.5     2389670 2389792 2389670 2389792 1       7.1e-64 .
Hb8>scaffold_4      Dappu1_FM5_230334.6     2389863 2390088 2389863 2390088 1       1.3e-124        .
Hb8>scaffold_4      Dappu1_FM5_230334.7     2390156 2390359 2390156 2390359 1       1.4e-110        .



=cut

=item output partway into this

../tandemgenes.pl -debug -act find -genome scaffold_4.fa -query jgifm4_exons.fa \
  -MINEVAL 1e-30  -uneq -neardist 20000   \
  -loc scaffold_4:651000-689000

# grep >ids jgifm4_exons.fa
# exon_finder
# echo reusing old jgifm4_exons.fa.blatf8;
echo blat -oneOff=1 -minIdentity=80 scaffold_4.fa jgifm4_exons.fa -out=blast8 jgifm4_exons.fa.blatf8;
reusing old jgifm4_exons.fa.blatf8
blat -oneOff=1 -minIdentity=80 scaffold_4.fa jgifm4_exons.fa -out=blast8 jgifm4_exons.fa.blatf8

# Location restriction: scaffold_4:651000-689000
# EQ_legend: -2=near, -1=ignoredmatch, 0=unequal, 1=equal, 2=altequal, 3=insideof
#Genome___      GeneExonID______        OrigB   OrigE   AtB     AtE     EQ      Eval    AltIDs

scaffold_4      Dappu1_FM5_234437.8     663878  664019  651718  651855  -2      5.0e-53 .        ** putative rev-tandem; 6 of 8 exon matches
scaffold_4      Dappu1_FM5_234437.4     662308  662457  653140  653225  -2      3.0e-38 .
scaffold_4      Dappu1_FM5_234437.3     661775  662013  653520  653745  -2      2.3e-103        .
scaffold_4      Dappu1_FM5_234437.2     661649  661752  653776  653879  -2      4.1e-51 .
scaffold_4      Dappu1_FM5_234437.1     661211  661538  653988  654116  -2      1.5e-58 .
scaffold_4      Dappu1_FM5_234437.1     661211  661538  654118  654319  -2      3.6e-85 .
scaffold_4      Dappu1_FM5_234437.1     661211  661538  661211  661538  1       1.9e-184        .
scaffold_4      Dappu1_FM5_234437.2     661649  661752  661649  661752  1       2.3e-52 .
scaffold_4      Dappu1_FM5_234437.3     661775  662013  661775  662013  1       5.8e-132        .
scaffold_4      Dappu1_FM5_234437.4     662308  662457  662308  662457  1       1.4e-79 .
scaffold_4      Dappu1_FM5_234437.5     662502  662572  662502  662572  1       3.1e-33 .
scaffold_4      Dappu1_FM5_234437.6     662643  662747  662643  662747  1       6.2e-53 .
scaffold_4      Dappu1_FM5_234437.8     663878  664019  663878  664019  1       1.7e-74 .
scaffold_4      Dappu1_FM5_234438.2     668790  668880  668790  668880  1       1.1e-44 .
scaffold_4      Dappu1_FM5_234438.3     670034  670107  670034  670107  1       8.3e-35 .
scaffold_4      Dappu1_FM5_234439.1     670211  670297  670211  670297  1       7.0e-42 .
scaffold_4      Dappu1_FM5_234439.2     670361  670513  670361  670513  1       1.5e-80 .
scaffold_4      Dappu1_FM5_234439.3     670588  670664  670588  670664  1       1.7e-36 .
scaffold_4      Dappu1_FM5_234437.8     663878  664019  671502  671633  -2      7.8e-63 .
scaffold_4      Dappu1_FM5_234441.1     674594  676312  674594  676312  1       0.0e+00 .
scaffold_4      Dappu1_FM5_234441.2     676431  676591  676431  676591  1       2.2e-86 .
scaffold_4      Dappu1_FM5_234441.3     676738  676999  676738  676999  1       9.1e-147        .
scaffold_4      Dappu1_FM5_234442.2     677101  678567  677101  678567  1       0.0e+00 .
scaffold_4      Dappu1_FM5_234443.1     678581  680583  678581  680583  1       0.0e+00 .
scaffold_4      Dappu1_FM5_234443.3     681114  681296  681114  681296  1       2.1e-100        .
scaffold_4      Dappu1_FM5_234443.4     682577  682663  682577  682663  1       3.7e-43 .
scaffold_4      Dappu1_FM5_234443.5     683610  683732  683610  683732  1       2.2e-64 .
scaffold_4      Dappu1_FM5_234443.6     684772  684856  684772  684856  1       9.3e-42 .
scaffold_4      Dappu1_FM5_234443.7     685363  685721  685363  685721  1       4.0e-204        .


=cut

=item tandem gene look

  exon blat/megablast compare
  
  ## modified to pull exons; do this elsewhere ??
  $bg/mb/glean-gene/bin/glean-overlap \
  -gff $em/daphc/scaffold_4/dpulex1_predict.gff -fa $em/daphc/scaffold_4/dpulex1.fa \
  -act fasta cds//CDS:JGI > ! scaf4_exons.fa
  
  blat scaffold_4.fa scaf4_exons.fa scaf4_exons.blat
  
  blat -oneOff=1 -minIdentity=80 scaffold_4.fa scaf4_exons.fa scaf4_exons.blat2
   # -oneOff=1 adds 300; -minId=80 adds 1100; both add 2300 ; -minId=60 adds 100
  
  blat -oneOff=1 -minIdentity=60 scaffold_4.fa scaf4_exons.fa -out=blast8 scaf4_exons.blatf8
   
  $nb/megablast -g T -W 12 -t 21 -D 3 \
    -i scaf4_exons.fa -d scaffold_4.fa | sort -k1,1 -k9,9n > scaf4_exons.mblast
      .. disc-megablast this way finds 2x blat
   
  cat scaf4_exons.blat | perl -ne'@v=split"\t"; print"$v[9]\n";' | sort | uniq -c |\
  sort -k1,1nr | grep -v ' 1 ' | sort -k2,2 > scaf4_exons.duplist
    .. ^ revise this part for full compare of dup exons with dpulex1_predict.gff
    .. look for dups of most gene exons
    .. find which gene models overlap dups: none, same as primary (error), other (full|partial) model
    .. ignore, count or treat same the distant dupls and near dupls?

=item blat for exons of tandem dupl (7 or 8 Car in this region)

dupl region: scaffold_4:2423251..2448250

genewise model mixes 4 tandem Car 5-exon genes
Dappu1_FM5_42005  	scaffold_4:2431108..2440394 (+ strand)

:: scaf4_exons.blat2
psLayout version 3

match   mis-    rep.    N's     Q gap   Q gap   T gap   T gap   strand  Q               Q       Q       Q       T               T       T       T       block   blockSizes      qStarts  tStarts
        match   match           count   bases   count   bases           name            size    start   end     name            size    start   end     count
--------
223     0       0       0       0       0       0       0       +       Dappu1_FM5_42005.1:2431108..2431330     223     0       223     scaffold_4      3075709 2431107 2431330 1       223,    0,      2431107,
194     17      0       0       1       12      1       16      +       Dappu1_FM5_42005.1:2431108..2431330     223     0       223     scaffold_4      3075709 2428250 2428477 2       35,176, 0,47,   2428250,2428301,
142     0       0       0       0       0       0       0       +       Dappu1_FM5_42005.2:2431417..2431558     142     0       142     scaffold_4      3075709 2431416 2431558 1       142,    0,      2431416,
128     7       0       0       0       0       0       0       +       Dappu1_FM5_42005.2:2431417..2431558     142     0       135     scaffold_4      3075709 2428563 2428698 1       135,    0,      2428563,
39      7       0       0       0       0       0       0       +       Dappu1_FM5_42005.2:2431417..2431558     142     89      135     scaffold_4      3075709 2466663 2466709 1       46,     89,     2466663,
106     0       0       0       0       0       0       0       +       Dappu1_FM5_42005.3:2437524..2437629     106     0       106     scaffold_4      3075709 2437523 2437629 1       106,    0,      2437523,
231     0       0       0       0       0       0       0       +       Dappu1_FM5_42005.4:2437709..2437939     231     0       231     scaffold_4      3075709 2437708 2437939 1       231,    0,      2437708,
126     0       0       0       0       0       0       0       +       Dappu1_FM5_42005.5:2440269..2440394     126     0       126     scaffold_4      3075709 2440268 2440394 1       126,    0,      2440268,
102     7       0       0       0       0       0       0       +       Dappu1_FM5_42005.5:2440269..2440394     126     0       109     scaffold_4      3075709 2435445 2435554 1       109,    0,      2435445,

melon.%
 $nb/megablast -g T -W 12 -t 21 -i scaf4_exons.fa -d scaffold_4.fa -D 3 | grep Dappu1_FM5_42005 | sort -k1,1 -k9,9n
Dappu1_FM5_42005.1:2431108..2431330     scaffold_4      90.48   21      2       0       108     128     1707613 1707633 8.0     26.3
Dappu1_FM5_42005.1:2431108..2431330     scaffold_4      90.54   222     17      4       4       221     2428254 2428475 3e-73    270
Dappu1_FM5_42005.1:2431108..2431330     scaffold_4      100.00  223     0       0       1       223     2431108 2431330 4e-125   442
Dappu1_FM5_42005.1:2431108..2431330     scaffold_4      95.24   21      1       0       103     123     2437055 2437075 0.033   34.2
Dappu1_FM5_42005.1:2431108..2431330     scaffold_4      88.46   26      3       0       103     128     2441232 2441257 2.0     28.2
Dappu1_FM5_42005.1:2431108..2431330     scaffold_4      90.00   30      3       0       100     129     2463760 2463789 0.008   36.2
Dappu1_FM5_42005.1:2431108..2431330     scaffold_4      93.10   29      2       0       94      122     2466347 2466375 1e-04   42.1
Dappu1_FM5_42005.1:2431108..2431330     scaffold_4      93.10   29      2       0       94      122     2478249 2478277 1e-04   42.1
Dappu1_FM5_42005.1:2431108..2431330     scaffold_4      100.00  20      0       0       106     125     2482618 2482637 5e-04   40.1
Dappu1_FM5_42005.1:2431108..2431330     scaffold_4      100.00  17      0       0       109     125     2487164 2487180 0.033   34.2
Dappu1_FM5_42005.1:2431108..2431330     scaffold_4      95.00   20      1       0       202     221     2487257 2487276 0.13    32.2

Dappu1_FM5_42005.2:2431417..2431558     scaffold_4      94.81   135     7       0       1       135     2428564 2428698 4e-56    212
Dappu1_FM5_42005.2:2431417..2431558     scaffold_4      100.00  142     0       0       1       142     2431417 2431558 5e-77    281
Dappu1_FM5_42005.2:2431417..2431558     scaffold_4      90.91   22      2       0       40      61      2478497 2478518 1.2     28.2
Dappu1_FM5_42005.2:2431417..2431558     scaffold_4      94.74   19      1       0       117     135     2478574 2478592 0.31    30.2

Dappu1_FM5_42005.3:2437524..2437629     scaffold_4      100.00  14      0       0       6       19      1708011 1708024 0.89    28.2
Dappu1_FM5_42005.3:2437524..2437629     scaffold_4      87.18   39      5       0       15      53      2428792 2428830 0.001   38.2 ..
Dappu1_FM5_42005.3:2437524..2437629     scaffold_4      100.00  15      0       0       73      87      2428850 2428864 0.23    30.2
Dappu1_FM5_42005.3:2437524..2437629     scaffold_4      87.18   39      5       0       15      53      2431641 2431679 0.001   38.2
Dappu1_FM5_42005.3:2437524..2437629     scaffold_4      100.00  15      0       0       73      87      2431699 2431713 0.23    30.2
Dappu1_FM5_42005.3:2437524..2437629     scaffold_4      100.00  106     0       0       1       106     2437524 2437629 1e-55    210
Dappu1_FM5_42005.3:2437524..2437629     scaffold_4      95.00   20      1       0       5       24      2441641 2441660 0.057   32.2 ..
Dappu1_FM5_42005.3:2437524..2437629     scaffold_4      89.19   37      2       2       11      46      2466784 2466819 2e-04   40.1
Dappu1_FM5_42005.3:2437524..2437629     scaffold_4      85.71   42      6       0       5       46      2483024 2483065 0.004   36.2

Dappu1_FM5_42005.4:2437709..2437939     scaffold_4      100.00  231     0       0       1       231     2437709 2437939 7e-130   458
Dappu1_FM5_42005.4:2437709..2437939     scaffold_4      87.50   32      4       0       190     221     2467150 2467181 0.13    32.2
Dappu1_FM5_42005.4:2437709..2437939     scaffold_4      88.10   42      5       0       190     231     2475317 2475358 4e-05   44.1
Dappu1_FM5_42005.4:2437709..2437939     scaffold_4      87.50   32      4       0       190     221     2479022 2479053 0.13    32.2

Dappu1_FM5_42005.5:2440269..2440394     scaffold_4      93.64   110     5       2       1       109     2435446 2435554 1e-40    161
Dappu1_FM5_42005.5:2440269..2440394     scaffold_4      100.00  126     0       0       1       126     2440269 2440394 2e-67    250


  $nb/megablast -g F -W 12 -t 21 -i scaf4_exons.fa -d scaffold_4.fa -D 3 | grep Dappu1_FM5_42005 | sort -k1,1 -k9,9n
Dappu1_FM5_42005.1:2431108..2431330     scaffold_4      90.54   222     17      4       4       221     2428254 2428475 3e-73    270
Dappu1_FM5_42005.1:2431108..2431330     scaffold_4      100.00  223     0       0       1       223     2431108 2431330 4e-125   442
Dappu1_FM5_42005.1:2431108..2431330     scaffold_4      95.24   21      1       0       103     123     2437055 2437075 0.033   34.2
Dappu1_FM5_42005.1:2431108..2431330     scaffold_4      88.46   26      3       0       103     128     2441232 2441257 2.0     28.2
  >>> outside of this dupl. region
Dappu1_FM5_42005.1:2431108..2431330     scaffold_4      90.00   30      3       0       100     129     2463760 2463789 0.008   36.2
Dappu1_FM5_42005.1:2431108..2431330     scaffold_4      93.10   29      2       0       94      122     2466347 2466375 1e-04   42.1
Dappu1_FM5_42005.1:2431108..2431330     scaffold_4      93.10   29      2       0       94      122     2478249 2478277 1e-04   42.1
Dappu1_FM5_42005.1:2431108..2431330     scaffold_4      100.00  20      0       0       106     125     2482618 2482637 5e-04   40.1
   <<<
   
Dappu1_FM5_42005.2:2431417..2431558     scaffold_4      94.81   135     7       0       1       135     2428564 2428698 4e-56    212
Dappu1_FM5_42005.2:2431417..2431558     scaffold_4      100.00  142     0       0       1       142     2431417 2431558 5e-77    281
Dappu1_FM5_42005.2:2431417..2431558     scaffold_4      84.78   46      7       0       90      135     2466664 2466709 0.005   36.2
Dappu1_FM5_42005.2:2431417..2431558     scaffold_4      90.91   22      2       0       40      61      2478497 2478518 1.2     28.2

Dappu1_FM5_42005.3:2437524..2437629     scaffold_4      100.00  15      0       0       73      87      2431699 2431713 0.23    30.2
Dappu1_FM5_42005.3:2437524..2437629     scaffold_4      100.00  106     0       0       1       106     2437524 2437629 1e-55    210
Dappu1_FM5_42005.3:2437524..2437629     scaffold_4      95.00   20      1       0       5       24      2441641 2441660 0.057   32.2
Dappu1_FM5_42005.3:2437524..2437629     scaffold_4      89.19   37      2       2       11      46      2466784 2466819 2e-04   40.1
Dappu1_FM5_42005.3:2437524..2437629     scaffold_4      85.71   42      6       0       5       46      2483024 2483065 0.004   36.2

Dappu1_FM5_42005.4:2437709..2437939     scaffold_4      100.00  231     0       0       1       231     2437709 2437939 7e-130   458
Dappu1_FM5_42005.4:2437709..2437939     scaffold_4      83.33   42      7       0       190     231     2464546 2464587 2.1     28.2
Dappu1_FM5_42005.4:2437709..2437939     scaffold_4      88.10   42      5       0       190     231     2475317 2475358 4e-05   44.1

Dappu1_FM5_42005.5:2440269..2440394     scaffold_4      93.64   110     5       2       1       109     2435446 2435554 1e-40    161
Dappu1_FM5_42005.5:2440269..2440394     scaffold_4      100.00  126     0       0       1       126     2440269 2440394 2e-67    250

=cut

=item dgg notes 

newscript.pl -act fasta \
-gff $em/daphc/scaffold_4/dpulex1_predict.gff -fa $em/daphc/scaffold_4/dpulex1.fa \
cds//CDS:JGI > ! scaffold_4_exons.fa


cat dpulex1_*.gff | $go --gff - tscan//mRNA:twinscan est//HSP:pasa dmel//HSP:modDM mus//HSP:modMM

cat dpulex1_*.gff | $go --gff - genes//mRNA est//HSP:pasa dmel//HSP:modDM mus//HSP:modMM

cat dpulex1_*.gff | $go --gff - cds//CDS est//HSP:pasa dmel//HSP:modDM mus//HSP:modMM

dgbook% wc *gff
    3618   43406  309019 dpulex1_estmatch.gff
   11721  105482  936840 dpulex1_predict.gff
    5245   54170  538259 dpulex1_prot4hsp.gff

perl -pi.old -e'if(/twinscan/){ s/(ID|Parent)=scaffold_/$1=tscan.s/;s/TRscaffold_/tscan.s/g; }' dpulex1_predict.gff


cat dpulex1_*.gff | perl -ne'($r,$s,$t,@v)=split"\t";print"$t:$s\n";' | sort | uniq -c
3015 CDS:DGIL_SNO
1977 CDS:JGI
3001 CDS:twinscan
1122 HSP:modCE
1758 HSP:modDM
2141 HSP:modMM
 223 HSP:modSC
3617 HSP:pasa
2018 exon:JGI
 449 gene:twinscan
 472 mRNA:DGIL_SNO
 339 mRNA:JGI
 449 mRNA:twinscan



=cut