#! /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 fasta gtf newgtf" if (@ARGV<2); 	
my $thr= $opts{f};
my $ratio=$opts{r};
my $futher_file=$opts{f};

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

my %allgenes;
open my $hfa, "<$fasta";
while (<$hfa>) {
	if ($_ =~/^>(\S+)/ ) {
		$allgenes{$1}=1;	
	}
}


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;

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 }->{"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+)\";/;
	}
	if ( $s[2] eq "exon") {
		($exons)= $_ =~ /exon_number\s\"(\S+)\";/;
		$geneid2seq{ $geneid }-> {"exon"} -> {join("\t", $s[0],$s[3],$s[4])} =1; 

		if (exists $geneid2seq{ $geneid }->{"transcript"}->{$trid} ->{'exon'} )  {
		  #$geneid2seq{ $geneid }->{"transcript"}->{$trid} ->{'exon'}=[ @{$geneid2seq{ $geneid }->{"transcript"}->{$trid} ->{'exon'}}, $s[3],$s[4] ]; 
		  push(@{$geneid2seq{ $geneid }->{"transcript"}->{$trid} ->{'exon'}}, $s[3],$s[4] ); 
		} else {
		  $geneid2seq{ $geneid }->{"transcript"}->{$trid} ->{'exon'}=[$s[3],$s[4] ]; 
		}
	}

	if ( exists $allgenes{$geneid."_From_".$chid} ) {
		my $nename=$geneid."_From_".$chid;
		my $newline=$_;
		$newline=~ s/$geneid/$nename/g;
		print {$hnewgtf} $newline;
	}
	
}
if ( $geneid) {
$geneid2seq{ $geneid }->{"transcript"}->{$trid} ->{'pos'} = [$tr_start, $tr_end];
$geneid2seq{ $geneid }->{"chr"}=$chid;
}

close $hgff;

foreach my $ke (keys %geneid2seq) {
    my $ts=scalar keys %{$geneid2seq{ $ke }->{"transcript"}};
    my $es=scalar keys %{$geneid2seq{ $ke }->{"exon"}};
    my %introns;
    for $trid ( keys %{$geneid2seq{ $ke }->{"transcript"}} ) {

	my @tr_pos =@{ $geneid2seq{ $ke }->{"transcript"}->{$trid} ->{'pos'} };
	my @exons  =@{ $geneid2seq{ $ke }->{"transcript"}->{$trid} ->{'exon'} };

	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 $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";
}
close $hfa;
