#! /usr/bin/perl
use strict;
use Getopt::Std;
use FindBin;
use lib "$FindBin::Bin/.";
use File::Temp qw/ tempfile tempdir /;

my %opts=(f=>undef, t=>0, r=>0.9);
getopts('f:t:r:', \%opts);

die "usage:repro_gtf.pl  gtf newgtf" if (@ARGV<2); 	
my $thr= $opts{f};
my $ratio=$opts{r};
my $futher_file=$opts{f};


my $gff=shift;
my $newgtf=shift;



open my $hgff, "<$gff";
open my $hnewgtf, ">$newgtf"; ##use this temoporatorily
my $geneid;
my $trid;
my $trseq;
my $chid;
my $fpkm;
my $exons;
my $transcripts;
my %geneid2seq;

my $tr_start;
my $tr_end;
my $fpkm_thr=0.1;

while (<$hgff>) {
	my @s= split /\t/;
	my $gg=$s[8];	
	if ( $s[2] eq "transcript" ) {
		if ( $geneid) {
		    $geneid2seq{ $geneid }->{"transcript"}->{$trid} ->{"pos"} = [$tr_start, $tr_end];
		    $geneid2seq{ $geneid }->{"transcript"}->{$trid} ->{"fpkm"} =$fpkm;
		    $geneid2seq{ $geneid }->{"chr"}=$chid;

		}
		$chid=$s[0];
		$trseq="";
		$tr_start=$s[3];
		$tr_end=$s[4];

		#$exons=1;	
		my @g=split($s[8], /;/);
		($geneid)= $gg =~ /gene_id\s\"(\S+)\";/;
		($trid)= $gg =~ /transcript_id\s\"(\S+)\";/;
		($fpkm)= $gg =~ /FPKM\s\"(\S+)\";/;
	}
	if ( $s[2] eq "exon") {
		my ($fpkm2)= $gg =~ /FPKM\s\"(\S+)\";/;
		next if ($fpkm2) < $fpkm_thr;
		$geneid2seq{ $geneid }-> {"exon"} -> {join("\t", $s[0],$s[3],$s[4])} =1; 

		if (exists $geneid2seq{ $geneid }->{"transcript"}->{$trid} ->{"exon"} )  {
		   push(@{$geneid2seq{ $geneid }->{"transcript"}->{$trid} ->{"exon"}}, $s[3],$s[4] ); 
		} else {
		   $geneid2seq{ $geneid }->{"transcript"}->{$trid} ->{"exon"}=[$s[3],$s[4] ]; 
		}
	}
	
}
if ( $geneid) {
    $geneid2seq{ $geneid }->{"transcript"}->{$trid} ->{"pos"} = [$tr_start, $tr_end];
    $geneid2seq{ $geneid }->{"transcript"}->{$trid} ->{"fpkm"} =$fpkm;
    $geneid2seq{ $geneid }->{"chr"}=$chid;
}

close $hgff;

foreach my $ke (keys %geneid2seq) {
    my %introns;

    for $trid ( keys %{$geneid2seq{ $ke }->{"transcript"}} ) {

	my @tr_pos =@{ $geneid2seq{ $ke }->{"transcript"}->{$trid} ->{"pos"} };
	my @exons;
	@exons  =@{ $geneid2seq{ $ke }->{"transcript"}->{$trid} ->{"exon"} }  if ( exists $geneid2seq{ $ke }->{"transcript"}->{$trid} ->{"exon"}) ;
	
 	if ( $geneid2seq{ $ke }->{"transcript"}->{$trid} ->{"fpkm"} < $fpkm_thr) {
	    delete $geneid2seq{ $ke }->{"transcript"}->{$trid};
	    next;
	}

	if ($tr_pos[0] > $tr_pos[1] ) {
		print STDERR "reverse transcribled gene\n";
		exit;
	}
	for (my $i=1; $i < @exons-1; $i +=2) {
	    if ( $exons[$i]+1 <$exons[$i+1]) {
		$introns{ join("\t", $exons[$i]+1, $exons[$i+1]-1)}=1;
	   }
	} 
    }
    my $ts=scalar keys %{$geneid2seq{ $ke }->{"transcript"}};
    my $es=scalar keys %{$geneid2seq{ $ke }->{"exon"}};
    next if ($ts <1 || $es <1);

    my $num_intron=scalar keys(%introns);


    #print join("\t", $ke."_From_".$geneid2seq{$ke}->{"chr"}, $ts, $es, $num_intron), "\n";
    print join("\t", $ke, $ts, $es, $num_intron), "\n";
}

