#!/usr/bin/perl
# repeatfilter.perl
# apply (dros) identified/predicted TE repeat GFF to tandy results
# -- blat output?  -- tandy.gff

use strict;
use warnings;
use Getopt::Long;

our $BINSIZE  = 1000 ; #was# 5000;
my ($repeats,$repeatlist,$gff,$action,$act,$ok);
my ($ACT_DROP,$ACT_KEEP,$ACT_MARK)=(1..9);
my $debug=1;

my $optok= GetOptions(
  "repeats=s", \$repeats, 
  "gff=s", \$gff,  
  "action=s", \$action, 
  "debug!", \$debug, 
  );

die 
"usage:  perl repeatfilter -act keep|drop|mark -rep repeats.gff -gff tandy.gff > tandyfilt.gff
" unless($optok and $action and -f $repeats);

$act= ($action =~ /keep/) ? $ACT_KEEP : ($action =~ /mark/) ? $ACT_MARK : $ACT_DROP;
$ok = ($repeats =~ /.gz$/) ? open(R,"gunzip -c $repeats |") : open(R,$repeats);
die "bad $repeats" unless($ok);
$repeatlist= collect_repeats(*R); close(R);

my $nr=0;
my $gffh= *STDIN;
$ok = ($gff =~ /.gz$/) ? open($gffh,"gunzip -c $gff |") 
      : ($gff =~ /^stdin|-/) ? $gff= *STDIN
      : open($gffh,$gff);
      
die "bad $gff" unless($ok);
while(<$gffh>){
  unless(/^\w/){ next if(/^(#n |$)/); print and next; }
  my($ref,$src,$typ,$ab,$ae,@gff)= split"\t";
  my $isrep= repeat_overlaps($ref,$ab,$ae);
  if($isrep) {
    $nr++;
    next if($act == $ACT_DROP);
    s/$/;terepeat=1/; # $ACT_MARK
  } else { 
    next if($act == $ACT_KEEP); 
  }
  print;
}  
warn"#repeats found=$nr\n" if $debug;


sub repeat_overlaps
{
  my($ref,$ab,$ae)= @_;    
  return 0 unless($repeatlist->{$ref});
  my @bins= (int($ab/$BINSIZE) .. int($ae/$BINSIZE));
  LBIN: foreach my $ib (@bins) {
    $repeatlist->{$ref}{$ib} or next;
    my @locs= @{$repeatlist->{$ref}{$ib}};
    foreach my $rloc (@locs) {
      my ($lb,$le)= @{$rloc}[0,1];
      if(_isoverlap($ab,$ae,$lb,$le)) { return 1; }
      }
    }
  return 0;
}

sub collect_repeats
{
  my($gff)= @_;
  my %repeats=(); my $nr=0;
  while(<$gff>){
    next unless(/^\w/);
    my($ref,$src,$typ,$ab,$ae)= split"\t";
    my $rloc= [$ab,$ae,$ref]; $nr++;
    my @bins= (int($ab/$BINSIZE) .. int($ae/$BINSIZE));
    foreach my $ib (@bins) {
      push( @{$repeats{$ref}{$ib}}, $rloc);
      }
    }
  warn"#collect_repeats=$nr\n" if $debug;
  return \%repeats;
}


sub _isoverlap {
  my($gb,$ge, $qb,$qe)= @_; 
  return ($gb <= $qe && $ge >= $qb) ? 1 : 0;
}