#!/usr/bin/perl

#delete peaks with enrichment fold less than "v"
#outputs cumulative numbers of peaks for each chromosome

use strict;
use warnings;
use diagnostics;

my $usage = qq{
    createLog.pl

    -----------------------------
    mandatory parameters:
    
    -f input_file               input file 
    -c control_file            	control file    
    
    -----------------------------
    optional parameters:
    -r ratio		ratio between ChIP and control sizes
     -n name (defaut "log_") output name
};

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

## mandatory arguments

my $input_file = "";
my $control_file = "";

## optional arguments
my $name = "summary_";
my $ratio = 1;
## 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 '-r') {$ratio = 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";
}

my $output_fname = $name.$input_file;

$output_fname =~ s/.peaks//;

## print parameters

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

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

my %scorControl; 

open CTRL, "< $control_file" || die "$control_file: $!\n";
#<CTRL>;
while(<CTRL>){
	chomp;
	next if (m/id/i);
    	my @cur_fields = split(/\t/,$_);
    	my $corr = 0;
	if ($cur_fields[0] =~/chr/) {
	  $corr = -1;  
	}

    	my $chr = $cur_fields[1+$corr]; 
    	$chr=~s/chr//;    	
    	my $score = $cur_fields[5+$corr];
    	
    	unless (exists ($scorControl{$chr})) {
		my @arrS;    		
		$scorControl{$chr} = \@arrS; 
    	}
    	my $index = &ind($score);   
    	
    	#$a = $ok ? $b : $c;  # get a scalar

unless (exists($scorControl{$chr}->[$index])) {
		$scorControl{$chr}->[$index]= 0;
	}
    	$scorControl{$chr}->[$index]++;    	
    	    	
}
close CTRL;

my %scorChip; 
my $maxScore = 0;

open INPUT, "< $input_file" || die "$input_file: $!\n";
#<INPUT>;
#$_ = <INPUT>;
#print OUT;
while(<INPUT>){
    chomp;
    next if (m/id/i);

    my @cur_fields = split(/\t/,$_);
    my $corr = 0;
    if ($cur_fields[0] =~/chr/) {
	  $corr = -1;  
    }
    my $chr = $cur_fields[1+$corr];
    $chr=~s/chr//;    
    my $score = $cur_fields[5+$corr];    
    unless (exists ($scorChip{$chr})) {
	my @arrS;    		
	$scorChip{$chr} = \@arrS; 
    }
    my $index = &ind($score);    
    $maxScore = $maxScore<$index ?$index : $maxScore;

	unless (exists($scorChip{$chr}->[$index])) {
		$scorChip{$chr}->[$index]= 0;
	}
   # $a = $ok ? $b : $c;  # get a scalar
   $scorChip{$chr}->[$index] ++;    	
}

close INPUT;

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

for my $chr(sort byName keys %scorChip) {
	print OUT "Info:    Current chromosome : $chr\n";

	my @valueToAddChIP;
	my @valueToAddControl;

	my $class = $maxScore;
	$scorChip{$chr}->[$class ] = 0 unless ($scorChip{$chr}->[$class]);			
	$scorControl{$chr}->[$class ] = 0 unless ($scorControl{$chr}->[$class]);			

	#print OUT  "Info:    ",$class,"	",$class,"	",$scorChip{$chr}->[$class],"	",$scorControl{$chr}->[$class],"	",$scorControl{$chr}->[$class]*$ratio,"	",$scorChip{$chr}->[$class]-$scorControl{$chr}->[$class]*$ratio,"\n";
	$valueToAddChIP[$class] = $scorChip{$chr}->[$class];
	$valueToAddControl[$class] = $scorControl{$chr}->[$class];
	
	
	
	for ($class = $maxScore-1; $class>=0; $class--) {
		$scorChip{$chr}->[$class ] = 0 unless ($scorChip{$chr}->[$class]);			
		$scorControl{$chr}->[$class ] = 0 unless ($scorControl{$chr}->[$class]);	
		$valueToAddChIP[$class] = $valueToAddChIP[$class+1];
		$valueToAddControl[$class]=$valueToAddControl[$class+1];
		$valueToAddChIP[$class] += $scorChip{$chr}->[$class];
		$valueToAddControl[$class] += $scorControl{$chr}->[$class];
	}
	
	for $class ( 0..$maxScore) {
		print OUT  "Info:    ",$class,"	",$class,"	",$valueToAddChIP[$class] ,"	",$valueToAddControl[$class] ,"	",$valueToAddControl[$class]*$ratio,"	",$valueToAddChIP[$class]-$valueToAddControl[$class]*$ratio,"\n";
	}	
	
}

close OUT;

sub ind {
	my $score = $_[0];
	$score =~ m/(^\d+)/;
	my $int = $1;
	
	if ($score - $int > 0.5) {
		return $int+1;
	}
	return $int;
}

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

sub byName {
	my $left = $a;
	my $right = $b;
	$left =~ s/.fa.*//;

	if ($left =~ m/chr(.*)/) {
		$left = $1; 
	}

	$right =~ s/.fa.*//;
	if ($right =~ m/chr(.*)/) {
		$right = $1; 
	}
 	return ($left <=> $right) if ($left =~ /^\d/)&&($right =~ /^\d/) ;
	$left cmp $right;
}

