#!/usr/bin/perl

#delete peaks with enrichment fold less than "v"

use strict;
use warnings;
use diagnostics;

my $usage = qq{
    filterPeaks5.2.pl

    -----------------------------
    mandatory parameters:
    
    -f input_file               input file 
    -c control_file            	control file    
    
    -----------------------------
    optional parameters:
    
    -w flank	  (delault 0) flanks to add to window length in contol data
    -v enrichment (delault 2)	enrichment fold over control data
    -t threshold (default 0)	threshold on density	
    -n name (defaut "filter_") output name
};

if(scalar(@ARGV) == 0){
    print $usage;
    exit(0);
}

## mandatory arguments

my $input_file = "";
my $control_file = "";
my $window = 1;
my $seuil = 0;
my $name = "filter_";

my $fold = 2;

## optional arguments

## parse command line arguments

while(scalar(@ARGV) > 0){
    my $this_arg = shift @ARGV;
    if ( $this_arg eq '-h') {print "$usage\n"; exit; }

    elsif ( $this_arg eq '-f') {$input_file = shift @ARGV;}
    elsif ( $this_arg eq '-c') {$control_file = shift @ARGV;}

    elsif ( $this_arg eq '-w') {$window = shift @ARGV;}
    elsif ( $this_arg eq '-v') {$fold = shift @ARGV;}
    elsif ( $this_arg eq '-t') {$seuil = shift @ARGV;}
    elsif ( $this_arg eq '-n') {$name = shift @ARGV;}

    elsif ( $this_arg =~ m/^-/ ) { print "unknown flag: $this_arg\n";}
}


if ( $input_file eq ""){
    die "you should specify input file\n";
}
if( $control_file eq ""){
    die "you should specify control file\n";
}

if ($window<=0) {
	$window = 1;
}

## print parameters

print "\n-----------------\n\n";
print "input file: $input_file\n";
print "control file: $control_file\n";
print "window length: $window\n";

print "\n-----------------\n\n";



my %contrlL; 

my %contrlR; 
my %scor; 


open CTRL, "< $control_file" || die "$control_file: $!\n";
<CTRL>;
while(<CTRL>){
	chomp;
    	my @cur_fields = split(/\t/,$_);
    	my $chr = $cur_fields[1];    	
    	unless ($chr=~/chr/) {
    		$chr = "chr".$chr;
    	}
    	my $fp = $cur_fields[2];
    	my $ep = $cur_fields[3];
    	my $score = $cur_fields[5];
    	if (exists ($contrlL{$chr})) {
    		push (@{$contrlL{$chr}},$fp);
    		push (@{$contrlR{$chr}},$ep);
    		push (@{$scor{$chr}},$score);
    	}
    	else {
    		my @arrL = ($fp);
    		my @arrR = ($ep);
    		my @arrS = ($score);
    		$contrlL{$chr} = \@arrL;   
    		$contrlR{$chr} = \@arrR; 
    		$scor{$chr} = \@arrS; 
    	}
    	
}
close CTRL;

open INPUT, "< $input_file" || die "$input_file: $!\n";
#<INPUT>;
$_ = <INPUT>;

$input_file =~ s/.*[\\\/]//;
my $output_fname = $name.$input_file;

if(-e $output_fname){ unlink("$output_fname");}
open OUT, "> $output_fname" || die "$output_fname: $!\n";


print OUT;
while(<INPUT>){
    chomp;
    my @cur_fields = split(/\t/,$_);
    my $corr = 0;
    if ($cur_fields[0] =~/chr/) {
	  $corr = -1;  
    }
    my $chr = $cur_fields[1+$corr];
    unless ($chr=~/chr/) {
        $chr = "chr".$chr;
    }
    my $fp = $cur_fields[2+$corr];
    my $ep = $cur_fields[3+$corr];
    my $score = $cur_fields[5+$corr];
    
    if ($score<=$seuil) {
	next;
    }  

    my @arrL = @{$contrlL{$chr}};
    my @arrR = @{$contrlR{$chr}};
    my @arrS = @{$scor{$chr}};

    my $check = &check($fp,$ep,$score,\@arrL,\@arrR,\@arrS,$window,$fold);
    if ($check) {
    	print OUT;
    	print OUT "\n";
    }
	else {
		print ;	print "\n";
	}
    
}

close INPUT;
close OUT;


sub check{
	my ($x,$y,$score,$arrL,$arrR,$arrS,$window,$fold) = @_;
	
	for my $i (1..@{$arrL}-1) {
		if (abs($x-$arrL->[$i])<100*$window){
			my $dist = &dist($x,$y,$arrL->[$i],$arrR->[$i]);
			my $ratio = $score/($arrS->[$i]);
			if ((abs($dist)<= $window)&&($ratio<=$fold)) {
				#print "$score / $arrS->[$i] = $ratio, dist = $dist";	print "\n";
				return 0 ;
			}
		}	
	}
	1;
}



sub minDist {
	my ($x,$y,$arrL,$arrR) = @_;
	
	my $minD = abs($x-$arrL->[0]); 
	for my $i (1..@{$arrL}-1) {
		#print "HERE: $arrL->[$i],$arrR->[$i]\n" if (($x == 132731207)&&($y==132732533));
		if (abs($x-$arrL->[$i])<100*$window){
		print "$arrL->[$i],$arrR->[$i]\n" if (($x == 132731207)&&($y==132732533));
			$minD = &min ($minD,&dist($x,$y,$arrL->[$i],$arrR->[$i]));
		}	
	}
	$minD;
}

sub min {
	my ($a,$b) = @_;
	if ($a<$b) {
		return $a;
	}
	$b;
}

sub dist {
	my ($x,$y,$a,$b) = @_;
	if ($y<$a) {
		return $a-$y;
	}
	if ($b<$x) {
			return $x-$b;
	}
	0;
}
