#!/usr/bin/env perl
use strict;
use Getopt::Long;
use File::Basename qw( fileparse );

my $proj = 'EVM';
my $annotation_version = '1A';
my $src = 'EVM_PASA';
# # PASA_UPDATE: evm.model.CM3.4_scaffold00009.45, single gene model update, valid-1, status:[pasa:asmbl_14713,status:16],
#  valid-1
# # PASA_UPDATE: evm.model.CM3.4_scaffold00009.45.1, alt-splice addition, valid-1, status:[pasa:asmbl_14714,status:25], va
# lid-1
# CM3.4_scaffold00009     .       gene    2205345 2213448 .       +       .       ID=evm.TU.CM3.4_scaffold00009.45;Name=EV
# M%20prediction%20CM3.4_scaffold00009.45
# CM3.4_scaffold00009     .       mRNA    2205345 2213448 .       +       .       ID=evm.model.CM3.4_scaffold00009.45;Pare
# nt=evm.TU.CM3.4_scaffold00009.45
# CM3.4_scaffold00009     .       five_prime_utr  2205345 2205472 .       +       .       ID=evm.model.CM3.4_scaffold00009
# .45.utr5p1;Parent=evm.model.CM3.4_scaffold00009.45
# CM3.4_scaffold00009     .       exon    2205345 2205477 .       +       .       ID=evm.model.CM3.4_scaffold00009.45.exon
# 1;Parent=evm.model.CM3.4_scaffold00009.45
# CM3.4_scaffold00009     .       CDS     2205473 2205477 .       +       0       ID=evm.model.CM3.4_scaffold00009.45.cds1
# ;Parent=evm.model.CM3.4_scaffold00009.45
# CM3.4_scaffold00009     .       exon    2206200 2206223 .       +       .       ID=evm.model.CM3.4_scaffold00009.45.exon
# 2;Parent=evm.model.CM3.4_scaffold00009.45
# CM3.4_scaffold00009     .       CDS     2206200 2206223 .       +       1       ID=evm.model.CM3.4_scaffold00009.45.cds2
# ;Parent=evm.model.CM3.4_scaffold00009.45
# CM3.4_scaffold00009     .       exon    2207089 2207194 .       +       .       ID=evm.model.CM3.4_scaffold00009.45.exon
# 3;Parent=evm.model.CM3.4_scaffold00009.45
# CM3.4_scaffold00009     .       CDS     2207089 2207194 .       +       1       ID=evm.model.CM3.4_scaffold00009.45.cds3
# ;Parent=evm.model.CM3.4_scaffold00009.45
# CM3.4_scaffold00009     .       exon    2207297 2207365 .       +       .       ID=evm.model.CM3.4_scaffold00009.45.exon
# 4;Parent=evm.model.CM3.4_scaffold00009.45
# CM3.4_scaffold00009     .       CDS     2207297 2207365 .       +       0       ID=evm.model.CM3.4_scaffold00009.45.cds4
# ;Parent=evm.model.CM3.4_scaffold00009.45
# CM3.4_scaffold00009     .       exon    2210314 2210807 .       +       .       ID=evm.model.CM3.4_scaffold00009.45.exon
# 5;Parent=evm.model.CM3.4_scaffold00009.45
# CM3.4_scaffold00009     .       CDS     2210314 2210807 .       +       0       ID=evm.model.CM3.4_scaffold00009.45.cds5
# ;Parent=evm.model.CM3.4_scaffold00009.45
# CM3.4_scaffold00009     .       exon    2210986 2211247 .       +       .       ID=evm.model.CM3.4_scaffold00009.45.exon
# 6;Parent=evm.model.CM3.4_scaffold00009.45
# CM3.4_scaffold00009     .       CDS     2210986 2211247 .       +       1       ID=evm.model.CM3.4_scaffold00009.45.cds6
# ;Parent=evm.model.CM3.4_scaffold00009.45
# CM3.4_scaffold00009     .       exon    2212711 2212919 .       +       .       ID=evm.model.CM3.4_scaffold00009.45.exon
# 7;Parent=evm.model.CM3.4_scaffold00009.45
# CM3.4_scaffold00009     .       CDS     2212711 2212919 .       +       0       ID=evm.model.CM3.4_scaffold00009.45.cds7
# ;Parent=evm.model.CM3.4_scaffold00009.45
# CM3.4_scaffold00009     .       exon    2213033 2213448 .       +       .       ID=evm.model.CM3.4_scaffold00009.45.exon
# 8;Parent=evm.model.CM3.4_scaffold00009.45
# CM3.4_scaffold00009     .       CDS     2213033 2213099 .       +       1       ID=evm.model.CM3.4_scaffold00009.45.cds8
# ;Parent=evm.model.CM3.4_scaffold00009.45
# CM3.4_scaffold00009     .       three_prime_utr 2213100 2213448 .       +       .       ID=evm.model.CM3.4_scaffold00009
# .45.utr3p1;Parent=evm.model.CM3.4_scaffold00009.45
# CM3.4_scaffold00009     .       mRNA    2205345 2211520 .       +       .       ID=evm.model.CM3.4_scaffold00009.45.1;Pa
# rent=evm.TU.CM3.4_scaffold00009.45
# CM3.4_scaffold00009     .       five_prime_utr  2205345 2205472 .       +       .       ID=evm.model.CM3.4_scaffold00009
# .45.1.utr5p1;Parent=evm.model.CM3.4_scaffold00009.45.1
# CM3.4_scaffold00009     .       exon    2205345 2205477 .       +       .       ID=evm.model.CM3.4_scaffold00009.45.1.ex
# on1;Parent=evm.model.CM3.4_scaffold00009.45.1
# CM3.4_scaffold00009     .       CDS     2205473 2205477 .       +       0       ID=evm.model.CM3.4_scaffold00009.45.1.cd
# s1;Parent=evm.model.CM3.4_scaffold00009.45.1
# CM3.4_scaffold00009     .       exon    2206200 2206223 .       +       .       ID=evm.model.CM3.4_scaffold00009.45.1.ex
# on2;Parent=evm.model.CM3.4_scaffold00009.45.1
# CM3.4_scaffold00009     .       CDS     2206200 2206223 .       +       1       ID=evm.model.CM3.4_scaffold00009.45.1.cd
# s2;Parent=evm.model.CM3.4_scaffold00009.45.1
# CM3.4_scaffold00009     .       exon    2207089 2207194 .       +       .       ID=evm.model.CM3.4_scaffold00009.45.1.ex
# on3;Parent=evm.model.CM3.4_scaffold00009.45.1
# CM3.4_scaffold00009     .       CDS     2207089 2207194 .       +       1       ID=evm.model.CM3.4_scaffold00009.45.1.cd
# s3;Parent=evm.model.CM3.4_scaffold00009.45.1
# CM3.4_scaffold00009     .       exon    2207297 2207365 .       +       .       ID=evm.model.CM3.4_scaffold00009.45.1.ex
# on4;Parent=evm.model.CM3.4_scaffold00009.45.1
# CM3.4_scaffold00009     .       CDS     2207297 2207365 .       +       0       ID=evm.model.CM3.4_scaffold00009.45.1.cd
# s4;Parent=evm.model.CM3.4_scaffold00009.45.1
# CM3.4_scaffold00009     .       exon    2210314 2210807 .       +       .       ID=evm.model.CM3.4_scaffold00009.45.1.ex
# on5;Parent=evm.model.CM3.4_scaffold00009.45.1
# CM3.4_scaffold00009     .       CDS     2210314 2210807 .       +       0       ID=evm.model.CM3.4_scaffold00009.45.1.cd
# s5;Parent=evm.model.CM3.4_scaffold00009.45.1
# CM3.4_scaffold00009     .       exon    2210986 2211520 .       +       .       ID=evm.model.CM3.4_scaffold00009.45.1.ex
# on6;Parent=evm.model.CM3.4_scaffold00009.45.1
# CM3.4_scaffold00009     .       CDS     2210986 2211334 .       +       1       ID=evm.model.CM3.4_scaffold00009.45.1.cd
# s6;Parent=evm.model.CM3.4_scaffold00009.45.1
# CM3.4_scaffold00009     .       three_prime_utr 2211335 2211520 .       +       .       ID=evm.model.CM3.4_scaffold00009
# .45.1.utr3p1;Parent=evm.model.CM3.4_scaffold00009.45.1


# #PROT evm.model.CM3.4_scaffold00009.45 evm.TU.CM3.4_scaffold00009.45    MSSKEKPTLGGTRIKTRKRNITIPLDPASFSDAVVQIYLDNAGDLELI
# ARNIESSELDFSRYGDTFFEVIFTGGRTQPGTTKPDEGERHPYSIIDCEPKREAILPSVIYIQKILRRRPFLIKYLENVMRRFLQSLEFFEENERKKLAIFTALAFSQKLSGLPPETVFQ
# PLLKDNLVAKGLVLSFITDFFKEYLMDNSLDDLISILKRGKVEDNLLEFFPSAKRSAEGFSEHFTKEGLVPLVEYNAKKMFDVKLSEMKSALTTQIAEETDIAEVIETVKHRMRDANLPD
# IEVVRILWDVIMDAVQWSGKNQQQNANSALRQVKTWAELLNTFCSTGKLELELIYKVQMQCYEDAKLMKLFPEIVRSLYDQDVLAEDTILHWFRKGTNPKGRQSFVKALEPFVNWLEEAE
# EEE*

# #PROT evm.model.CM3.4_scaffold00009.45.1 evm.TU.CM3.4_scaffold00009.45  MSSKEKPTLGGTRIKTRKRNITIPLDPASFSDAVVQIYLDNAGDLELI
# ARNIESSELDFSRYGDTFFEVIFTGGRTQPGTTKPDEGERHPYSIIDCEPKREAILPSVIYIQKILRRRPFLIKYLENVMRRFLQSLEFFEENERKKLAIFTALAFSQKLSGLPPETVFQ
# PLLKDNLVAKGLVLSFITDFFKEYLMDNSLDDLISILKRGKVEDNLLEFFPSAKRSAEGFSEHFTKEGLVPLVEYNAKKMFDVKLSEMKSALTTQIAEETDIAEVIETVKHRMRDANLPD
# IEVVRILWDVIMDAVQWSGKNQQQNANSALRQVILHTLAYVPLIGTAIGLPTFSPLLCVC*

# Need to add GFF2 FITS (geneid-like) output
my $version = 3;
my $source = 0;
my $removestop = 0;
my $addstop = 0;
my $txtag = 'transcript';
my $mrna = 0;
my $ov = 3;
my $fits = 0;
my $fl = 0;
my $out = $$;
my $fname = '';
GetOptions(
	   'f|i|a:s'        => \$fname,
	   'v:s' => \$annotation_version,
	   'o:s' =>\$out,
	   'base|project|species|prefix:s' => \$proj
	  );
my ($base,$path,$ext) = fileparse($fname,qw(\.gff3 \.gff));
die "$fname not a gff3 file" if !$ext;
my %transcripts;
my %tx_gene;
my %genes;
my %tx;
my %proteins;
my $cluster = 0;
my $GENE = '';
my %gene_cluster;
my %starts;
my %ends;
my %cluster_start = undef;
my %cluster_end = undef;
my %clustergff;
my $check_start_stop_codon = 0;
#file is well sorted with gene records before transcripts before CDS and exon
#protein sequences follow 
open IN, "<$fname" or die "Couldn't open $fname: $!\n";
my $original = 1;
while (<IN>) {
  chomp;
  if (m/^#PROT/){
    my ($dummy,$tr,$ge,$pseq) = split;
    $proteins{$tr}=$pseq;
    #print STDERR "$tr\t$pseq\n";
  }
  if ( m/^[# ORIGINAL]/){
    $original = 1;
  }elsif( m/^[# PASA_UPDATE]/){
    $original = 0;
  }elsif( m/^[# ]/ || $_!~/\S/){
    next;
  }
  
  my @gtf = split "\t",$_;
  next if $gtf[2] !~/CDS|exon|gene|mRNA/i;
    
  #print STDERR $gtf[8],"\n";
  if ($gtf[2] eq "gene") {
    if ($gtf[8]=~m/ID=([^;]+)/) {
      #print STDERR "$1\n";
      $genes{$1}=\@gtf;
      if ($1 ne $GENE){
	$cluster++;
	$GENE=$1;	
	$cluster_start{$cluster}=$gtf[3];
	$cluster_end{$cluster}=$gtf[4];
	$clustergff{$cluster}=\@gtf;
      }
    }
  } elsif ($gtf[2] =~/mRNA/) {
    $gtf[2] =$txtag;
    my $gene =0;
    my $id =0;
    if ($gtf[8]=~m/Parent=([^;]+)/) {
      $gene = $1;
    }
    if ($gtf[8]=~m/ID=([^;]+)/) {
      $id = $1;
    }
    if ($id && $gene) {
      $tx_gene{$id}=$gene;
      #$gene_tx{$gene}=$id;
      #print STDERR $gtf[7],"\n";
      $tx{$id}->{gff}=\@gtf;
      $tx{$id}->{id}=$id;
      $tx{$id}->{gene}=$gene;
      $tx{$id}->{original}=$original;
      $gene_cluster{$gene}=$cluster;
      $tx{$id}->{cluster}=$cluster;
      #print STDERR "ID:$id\tgene:$gene\n";
    }
  } elsif ($gtf[2] =~/CDS/i) {
    if ($gtf[8]=~m/Parent=([^;]+)/) {
      push @{$tx{$1}->{cds}},\@gtf;
      if (exists($starts{$gtf[0]."-".$gtf[3]})){
	# if ($cluster ne $gene_cluster{$starts{$gtf[0]."-".$gtf[3]}}){
	#   print STDERR $starts{$gtf[0]."-".$gtf[3]}," with start: ",$gtf[0]."-".$gtf[3]," was cluster ",$gene_cluster{$starts{$gtf[0]."-".$gtf[3]}}," and is now cluster $cluster\n";
	# }
	my $previous_cluster = $gene_cluster{$starts{$gtf[0]."-".$gtf[3]}};
	foreach my $g (keys %gene_cluster){
	  if ($gene_cluster{$g} eq $previous_cluster){
	    $gene_cluster{$g}=$cluster;
	    $cluster_start{$cluster}=$genes{$g}->[3] if $genes{$g}->[3] < $cluster_start{$cluster};
	    $cluster_end{$cluster}=$genes{$g}->[4] if $genes{$g}->[4] > $cluster_end{$cluster};
	  }
	}
      }else{
	$starts{$gtf[0]."-".$gtf[3]}=$GENE;
      }
      if (exists($ends{$gtf[0]."-".$gtf[4]})){
	my $previous_cluster = $gene_cluster{$ends{$gtf[0]."-".$gtf[4]}};
	foreach my $g (keys %gene_cluster){
	  if ($gene_cluster{$g} eq $previous_cluster){
	    $gene_cluster{$g}=$cluster;
	    $cluster_start{$cluster}=$genes{$g}->[3] if $genes{$g}->[3] < $cluster_start{$cluster};
	    $cluster_end{$cluster}=$genes{$g}->[4] if $genes{$g}->[4] > $cluster_end{$cluster};
	  }
	}
      }else{
	$ends{$gtf[0]."-".$gtf[4]}=$GENE;
      }
    }
  } elsif ($gtf[2] =~/exon/i) {
    if ($gtf[8]=~m/Parent=([^;]+)/) {
      push @{$tx{$1}->{exons}},\@gtf;
    }
  }
}
close IN;
my @unsorted;
foreach my $t (keys %tx) {
  #my @CDS = sort {$a->start <=> $b->start} @{$t->{CDS}};
  #$t->{CDS}=\@CDS;
  push @unsorted, $tx{$t};
}
my @sorted_transcripts = sort tsort @unsorted;

my %seen_genes;
my %seen_cds;
my %seen_structure;
my $gene_count = 1;
my %tcount = 1;
my %ccount = 1;
my %pcount = 1;
my $newgenename = '';
open PFASTA,">$$.$proj.$annotation_version.protein.fa" or die "Couldn't open $proj.$annotation_version.protein.fa:$!";
#open CPFASTA,">$$.$proj.$annotation_version.ref_protein.fa" or die "Couldn't open $proj.$annotation_version.ref_protein.fa:$!";
open GFF,">$$.$proj.$annotation_version.gff3" or die "Couldn't open $proj.$annotation_version.gff3: $!";
print GFF "##gff-version 3\n";
my ($sec,$min,$hour,$mday,$mon,$year,$wday,$yday,$isdst)=localtime(time);
printf GFF "# date: %4d-%02d-%02d %02d:%02d:%02d\n",$year+1900,$mon+1,$mday,$hour,$min,$sec;
foreach my $transcript (@sorted_transcripts) {
  my $canon = 0;
  my @tgff = @{$transcript->{gff}};
  my $strand = $tgff[6];
  my $tid = $transcript->{id};
  my $prot = $proteins{$tid};
  #print STDERR "$prot\n";
  #print STDERR "$tid\n" if !$prot;
  my $plen = length($prot);
  my $gene_id = $gene_cluster{$transcript->{gene}};
  $newgenename = sprintf("%s%s%06d",$proj,$annotation_version,$gene_id);
  #my $dots = 0;
  #while($tid=~m/\./g){$dots++;};
  #if ($dots == 4){$canon = 1;}
  if (!exists $seen_genes{$gene_id}) {
    
    $clustergff{$gene_id}->[1]=$src;
    $clustergff{$gene_id}->[3]=$cluster_start{$gene_id};
    $clustergff{$gene_id}->[4]=$cluster_end{$gene_id};
    $clustergff{$gene_id}->[8]=~s/ID=[^;]+/ID=$newgenename/;
    $clustergff{$gene_id}->[8]=~s/;Name=[^;]+//;
    print GFF (join("\t",@{$clustergff{$gene_id}}),"\n");
    $seen_genes{$gene_id}++;
    $tcount{$newgenename} = 1;
    $ccount{$newgenename} = 0;
    $pcount{$newgenename} = 0;
  }
  # check for duplicate structure and choose the PASA_UPDATE over the ORIGINAL
  my $uniq_str = '';
  foreach my $gr (sort gffsort @{$transcript->{cds}}) {
    $uniq_str .=join(':',('c',$gr->[0],$gr->[3],$gr->[4],$gr->[6]));
  }
  $uniq_str .=":";
  foreach my $gr (sort gffsort @{$transcript->{exons}}) {
    $uniq_str .=join(':',('e',$gr->[0],$gr->[3],$gr->[4],$gr->[6]));
  }
  ### transcript sorted by PASA_UPDATE's followed by ORIGINAL
  if (! exists $seen_structure{$uniq_str}){
    #$seen_structure{$uniq_str} = $transcript;
    $seen_structure{$uniq_str}++;
  }else{
    print STDERR "$uniq_str seen already -- skipping\n";
    next;
  }
  my $newcdsid = $newgenename.'C'.$tcount{$newgenename};
  my $newtid = $newgenename.'T'.$tcount{$newgenename}++;
  $transcript->{gff}->[1]=$src;
  $transcript->{gff}->[8]=~s/ID=[^;]+/ID=$newtid/;
  $transcript->{gff}->[8]=~s/Parent=[^;]+/Parent=$newgenename/;
  print GFF (join("\t",@{$transcript->{gff}}),"\n");
 
  my @cds;
  my $count = 1;
  my $ecount = 1;
  my $numexon = 0;
  my $numcds = 0;
    

  my $string = '';
  my $cdslen = 0;
  foreach my $gffrec (sort gffsort @{$transcript->{cds}}) {
    $cdslen+=$gffrec->[4]-$gffrec->[3]+1;
    $numcds++;
    $string .=join(':',($gffrec->[0],$gffrec->[3],$gffrec->[6],$gffrec->[4],$gffrec->[6]));
  }
  if ($strand eq '+') {
    foreach my $gffrec (sort gffsort @{$transcript->{exons}}) {
      my $newexonid = $newtid.'.exon'.$ecount;
      $gffrec->[8]=~s/ID=[^;]+/ID=$newexonid/;
      $gffrec->[8]=~s/Parent=[^;]+/Parent=$newtid/;
      $gffrec->[1]=$src;
      print GFF (join("\t",@{$gffrec}),"\n");
      $ecount++;
    }
  } else {
    my $ecount = scalar @{$transcript->{exons}};
    foreach my $gffrec (sort gffsort @{$transcript->{exons}}) {
      my $newexonid = $newtid.'.exon'.$ecount;
      $gffrec->[8]=~s/ID=[^;]+/ID=$newexonid/;
      $gffrec->[8]=~s/Parent=[^;]+/Parent=$newtid/;
      $gffrec->[1]=$src;
      print GFF (join("\t",@{$gffrec}),"\n");
      $ecount--;
    }
  }
  my $newprotid = $newgenename.'P'.$pcount{$newgenename};

  if (!exists $seen_cds{$string}) {
    $pcount{$newgenename}++;
    $newprotid = $newgenename.'P'.$pcount{$newgenename};
    $seen_cds{$string}=$newprotid;
    my $wrap_length = 60;
    my $charcount = 0;
    print PFASTA ">$newprotid\n";
     #print separate file of canonical protein seqs (EVM consensus models not modified by PASA)
    #print CPFASTA ">$newprotid\n" if $canon;
    
    my @chars = split '',$prot;
    for (my $i=0;$i< $plen-1; $i++) {
      print PFASTA $chars[$i];
      #print CPFASTA $chars[$i] if $canon;
      if (++$charcount >= $wrap_length) {
	print PFASTA "\n" ;
	#print CPFASTA "\n" if $canon;
	$charcount = 0;
      }
    }
    print PFASTA $chars[-1],"\n";
    #print CPFASTA $chars[-1],"\n" if $canon;
  }else{
    $newprotid = $seen_cds{$string};
  }
  #print all CDS records
  $canon = 0;
  
  my $pstart = 1;
  my $cdsend = 3*($plen- ($prot=~/m\*$/?1:0));
  my $pend = $plen;
  my $sum = 0;
  my $first = 1;
  my $numcds = scalar  @{$transcript->{cds}};
  my $cdscount = 0;
  my @sortedcds = sort gffsort @{$transcript->{cds}};
    
# need to match up protein length with cds length. adjust first exon length by phase or remainder then keep summing exon lengths.
  if ($strand eq '+') {
    my $phase = $sortedcds[0]->[7];
    my $rmd = ($cdslen - $phase)%3;
    $first = 1;
    foreach my $gffrec (sort gffsort @{$transcript->{cds}}) {
      my $exonlen = $gffrec->[4]-$gffrec->[3]+1;
      if ($first){
	$pstart=1;
	$exonlen = $gffrec->[4]-$gffrec->[3]+1-$phase;
#	$sum = $gffrec->[7];
	 #print STDERR (join("\t",@{$gffrec}),"\n");
	$first = 0;
      }else{
	$exonlen = $gffrec->[4]-$gffrec->[3]+1;
	$pstart = int(($sum + 1 + 2)/3);
      }
      $sum+=$exonlen;
      $cdscount++;
      if ( $cdscount == $numcds){
	#$sum = $gffrec->[7];
	# print STDERR (join("\t",@{$gffrec}),"\n");
	$pend = $plen;
      }else{
	$pend = int(($sum + 2)/3);
      }
      $pend = int(($sum + 2)/3);
      $gffrec->[8]=~s/ID=[^;]+/ID=$newcdsid/;
      $gffrec->[8]=~s/Parent=[^;]+/Parent=$newtid/;
      $gffrec->[8].=";Target=$newprotid $pstart $pend";
      $gffrec->[1]=$src;
      print GFF (join("\t",@{$gffrec}),"\n");

    }
  } else {
    my $phase = $sortedcds[-1]->[7];
    my $rmd = ($cdslen - $phase)%3;
    $sum = $cdslen+1;
    $cdsend  - ($cdsend%3);
    #print STDERR "plen: $plen\n";
    $pstart = $plen;
    $pend = $plen;
    my $numcds = scalar  @{$transcript->{cds}};
    my $cdscount = 0;
    my $first = 1;
    foreach my $gffrec (sort gffsort @{$transcript->{cds}}) {
      my $exonlen = $gffrec->[4]-$gffrec->[3]+1;
      if ($first){
	$pend=$plen;
	$exonlen = $gffrec->[4]-$gffrec->[3]+1-$rmd;
      # 	$sum = $gffrec->[7];
      # 	 print STDERR (join("\t",@{$gffrec}),"\n");
       	$first = 0;
      }else{
	$exonlen = $gffrec->[4]-$gffrec->[3]+1;
	$pend = int(($sum - 1 +2)/3);
      }
      
      $sum-=$exonlen;
      $cdscount++;
      if ( $cdscount == $numcds){
	#$sum = $gffrec->[7];
	 #print STDERR (join("\t",@{$gffrec}),"\n");
	$pstart = 1;
      }else{
	$pstart = int(($sum+2)/3);
      }
      $gffrec->[8]=~s/ID=[^;]+/ID=$newcdsid/;
      $gffrec->[8]=~s/Parent=[^;]+/Parent=$newtid/;
      $gffrec->[8].=";Target=$newprotid $pstart $pend";
      $gffrec->[1]=$src;
      print GFF (join("\t",@{$gffrec}),"\n");
    }
  }
  $ccount{$newgenename}++;
  
  print GFF "###\n" if ($ov!=2);
}
close GFF;
close PFASTA;
#close CPFASTA;
print STDERR "done!\n";
sub tsort
  {
    $a->{gff}->[0] cmp $b->{gff}->[0]
      ||
	$a->{gff}->[3] <=> $b->{gff}->[3]
	  ||
	    $a->{gff}->[4] <=> $b->{gff}->[4]
	      ||
		$a->{gff}->[6] <=> $b->{gff}->[6]
		  ||
		    $a->{original} <=> $b->{original}
		  
	      }
sub gffsort
  {
    $a->[0] cmp $b->[0]
      ||
	$a->[3] <=> $b->[3]
	  ||
	    $a->[4] <=> $b->[4]
	      ||
		$a->[6] <=> $b->[6]
	      }

sub reverse_gffsort
  {
    $b->[0] cmp $a->[0]
      ||
	$b->[3] <=> $a->[3]
	  ||
	    $b->[4] <=> $a->[4]
	      ||
		$b->[6] <=> $a->[6]
	      }
