#!/usr/bin/perl
# exon2intron.pl
# assume input is GFF version 3
# add pull_intergenes() 


use strict;

my $exontype= $ENV{exon} || "exon"; # maybe exon,CDS : both? problems

my(%exons, %gstart, %gend, %badcall, %allin);
#?? need to collect all? or assume sorted at least by $r
my $lastr="";
my $allout=0;

print "##gff-version 3\n";
while(<>) {
  next unless(/^\w/); chomp;
  my($r,$s,$t,$rb,$re,$v,$o,$p,$a)=split"\t"; 
  if($lastr and $lastr ne $r) { 
    put_introns(); %exons=(); %gstart= %gend= (); %allin=(); 
    }
  my($gid)=m/(?:ID|Parent)=([^;]+)/; 
  if($t =~ /^$exontype/) { 
    #my $bad=(abs($re-$rb)<3); # also filter huge introns?
    #$badcall{$gid}++ if($bad);
    push( @{$exons{$gid}},[$rb,$re,$o,$p, $r,$s,$t,$v,$a]) ; # unless($bad); 
    $gstart{$gid}= $rb unless(defined $gstart{$gid} and $rb>$gstart{$gid});
    $gend{$gid}= $re unless(defined $gend{$gid} and $re<$gend{$gid});
    } 
  $lastr= $r;
}

put_introns();


sub _sort_geneloc {
  my ($ab,$ae)= ($gstart{$a}, $gend{$a}); 
  my ($bb,$be)= ($gstart{$b}, $gend{$b}); 
  return $ab <=> $bb or $ae <=> $be;
}

sub merge_att {
  my($at,$bt)= @_;
  my %at; 
  map{ 
    my($k,$v)=split "=",$_,2; 
    if($v=~/^\d/ and $at{$k} =~ /^\d/) { $at{$k}+= $v; }
    elsif($at{$k} and $v eq $at{$k}) {}
    elsif($at{$k}) { $at{$k} .= ",$v"; }
    elsif($k && $v) { $at{$k}.= $v; }
    } split ";\s*","$at;$bt";
  
  return join( ";", map{ $_."=".$at{$_} }sort keys %at);
}

sub put_introns {

  #?? output by id or by location?
  my @geneids= sort _sort_geneloc keys %exons; # sort keys %exons
  
  foreach my $ig (0 .. $#geneids) {
    my $id= $geneids[$ig] or next;
    
    my $flag="";
    my @ex = @{$exons{$id}}; 
  
    my $geneo=${$ex[0]}[2];
    my $rev= ($geneo eq "-" || $geneo < 0);
    my $nintron= @ex - 1;
  
      #?? want this rev sorting?
    if($rev) { @ex= sort{ $b->[1] <=> $a->[1] } @ex; }
    else { @ex= sort{ $a->[0] <=> $b->[0] } @ex; }
  
    my($lb,$le,$lo,$lp, $lat);
    my($rb,$re,$o,$p,  $r,$s,$t,$v,$at)= @{$ex[0]};
    $at =~ s/Parent=$id;?//;

    my $geneb= $rev ? $re : $rb;
    my $genee= $rev ? ${$ex[-1]}[0] : ${$ex[-1]}[1];
    # my $genespan= 1 + abs($genee - $geneb); # or rev

    if($ig>0) {
      my $lastid= $geneids[$ig-1];
      my $lend= $gend{$lastid};
      my $gstart= $gstart{$id};
      my $att= "ID=IG.${lastid}.${id}"; # add any gene annots?
      putgff($r, $s, "intergenic_region",$lend+1,$gstart-1,".",".",".",$att);
      }
    
    #$cdslen += 1 + abs($re-$rb) - $p; # is there ever phase on exon0 ?
    foreach my $i (1..$nintron) {
      ($lb,$le,$lo,$lp,$lat)= ($rb,$re,$o,$p,$at);
      ($rb,$re,$o,$p,  $r,$s,$t,$v,$at)= @{$ex[$i]};
      $at =~ s/Parent=$id;?//;
      
      my($inb,$ine)= ($rev) ? ($re+1,$lb-1) : ($le+1,$rb-1);
      # my $insize = 1 + $ine - $inb; # 
      # my $insize = $rev ? abs($lb - $re) : abs($rb - $le); # $insize--; # what?
      my $insize = $rev ? $lb - $re : $rb - $le; #? $insize--; # what?
      my $inlm   = $insize % 3;
  
      my $inatt="Parent=$id;in=$i";
      # $inatt .= ";li=$insize.$inlm"; ## ;ID=$id-in$i"; # what; $a ? dont need both ID and in=
      # note, should output in absolute location order, not by geneid && rev
      ## want to keep SOME of both exon annots: option
      # $inatt .= ";$lat" if($lat); $inatt .= ";$at" if ($at);
      $inatt .= ";" . merge_att($lat,$at) if($lat or $at);
       
      putgff( $r,$s,"intron",$inb,$ine,".",$lo,".",$inatt);
      }
    
    
    }

}

sub putgff {
  my @gff=@_; # ($r,$s,$t,$rb,$re,$v,$o,$p,$a)
  print join("\t",@gff),"\n";
}

__END__

## some odd modulus error here ....
# dgbook% gzcat daphgeno/dpulex_jgi060905_JGI_FM5.gff.gz | \
# egrep '^scaffold_([1-9][0-9]|[0-9])    ' | env exon="CDS" $gin | \
# perl -ne'if(m/li=(\d+).(\d)/){ ($l,$m)=($1,$2);  $m{"p".$m}++ if($l<100);} END{foreach $m (sort keys %m){print "$m $m{$m}\n";}}'
#
# # with $insize-- or 1 + $ine-$inb
# p0 17121   << are these extras, or is p1 below losing 200?
# p1 17523
# p2 17285
# 
# # without $insize--
# p0 17285
# p1 16877   <<
# p2 17523
# 


>> really need all of genome marked, annotated
>> as inverse of gene parts: introns, intergenic, utrs, along with
>> exons(CDS_exons, UTR_exons)


To: Don Gilbert <gilbertd@cricket.bio.indiana.edu>
Subject: Re: two things for Daphnia genome


Hello Don,

Thanks for your rapid reply.  We are almost done with analysing the tiling 
arrays data and we are ready to compare the transcribed fragments with the 
predicted gene models.  The GFF files that you produced earlier contains 
the coordinates of the exons and the CDS.  Can GFF files be created that 
contains the coordinates of the predicted introns (with and without EST 
support), of the predicted UTRs and the predicted intergenic spaces?  in 
this way, I will send you stats on the overlap betwen the tiling data and 
the gene model features.  More soon, John