#!/usr/bin/perl -w 
use strict;

my $usage = "Usage: $0 <gene models gff> <pol-II fold change gff> <SI-threshold>\n";
my $geneGff = $ARGV[0] or die $usage;
my $polIIFoldChangeGff = $ARGV[1] or die $usage;
my $threshold = $ARGV[2] or die $usage;
my $buffer = 300;
my($geneSet,$source) = readGeneGffFile($geneGff);
my $polIIData = readPolIIData($polIIFoldChangeGff);

foreach my $chrom (keys %{$geneSet}) {
    foreach my $geneInfo (@{$geneSet->{$chrom}}) {
	my($start,$stop,$strand,$id,$name,$note) = @{$geneInfo};
	my @overlappingValues;
	my $maxTss = 0;
	foreach my $tile (@{$polIIData->{$chrom}}) {
	    my($pos,$foldChange) = @{$tile};
	    if(($start - $buffer<= $pos)&&($pos <= $stop + $buffer)) {
		# store values for the median
		push(@overlappingValues,$foldChange);
		if($strand eq "+") {
		    if(($start - $buffer<= $pos)&&
		       ($pos <= $start + $buffer)) {
			# near the TSS
			if($foldChange >= $maxTss) {
			    # store the max value near the TSS
			    $maxTss = $foldChange;
			}
		    }
		} else {
		    if(($stop - $buffer<= $pos)&&
		       ($pos <= $stop + $buffer)) {
			# near the TSS             
			if($foldChange >= $maxTss) {
			    $maxTss = $foldChange;
			}
		    }
		}
	    }
	}
	my $median = getMedian(@overlappingValues);
	my $StallingIndex = $median != 0 ? $maxTss/$median : $maxTss/1.0;
	$StallingIndex = sprintf("%.3f",$StallingIndex);
	if($StallingIndex >= $threshold) {
	    print "$id\n";
	}
    }
}

sub getMedian {
    my @terms = @_;
    if(@terms) {
        my @sortedTerms = sort {$a <=> $b} @terms;
        my $size = scalar(@sortedTerms);
        my $middle = int($size/2);
        return $sortedTerms[$middle];
    } else {
        return 0.0;
    }
}

sub readPolIIData {
    my($polIIFile)=@_;
    my %polIIdata;
    open(PFILE,$polIIFile);
    while(<PFILE>) {
        chomp;
	my($chrom,$source,$type,$start,$stop,$foldChange,$strand,$phase,$info)=split(/\t/,$_);
	my $pos = int(($start+$stop)/2);
	# 1st fundamental snag:
	$chrom = "chr" . $chrom;
        push(@{$polIIdata{$chrom}},[$pos,$foldChange]);
    }
    return \%polIIdata;
}

sub readGeneGffFile {
    my($gffFile) = @_;
    my %dataSet;
    my %added;
    my $thisSource="";
    open(GFILE,$gffFile);
    while(<GFILE>) {
        unless(/\#/) {
            chomp;
            my($chrom,$source,$type,$start,$stop,$score,$strand,$phase,$info)=split(/\t/,$_);
            if($type eq "mRNA") {
                unless($thisSource) {
                    $thisSource = $source;
                }
                my $id;
                if(my($thisId) = $info =~ /ID=(.*?);/) {
                    $id = $thisId;
                } elsif(($thisId) = $info =~ /ID=(.*)/) {
                    $id = $thisId;
                } else {
                    die "could not parse line $_\n";
                }
                my($name,$note);
                my($thisName) = $info =~ /Name=(.*?);/;
                my($thisNote) = $info =~ /Note=(.*)/;
                my($thisFullName) = $info =~ /fullname=(.*?);/;
                $note = $thisNote ? $thisNote : $id;
                $note = ((!$thisNote)&&($thisFullName)) ? $thisFullName : $note;
                $name = $thisName ? $thisName : $id;
		# convert transcript ID to gene ID.
		($id) = $id =~ /(CG.*)\-R/;
		# store only unique TSSs.
		unless($added{$chrom.$start.$stop}) {
                    push(@{$dataSet{$chrom}},[$start,$stop,$strand,$id,$name,$note]);
                    $added{$chrom.$start.$strand}++;
                }
            }
        }
    }
    foreach my $chrom (keys %dataSet) {
        @{$dataSet{$chrom}} = sort {$a->[0] <=> $b->[0]} @{$dataSet{$chrom}};
    }
    return(\%dataSet,$thisSource);
}

