#!/usr/bin/env perl
use warnings;
use strict;

my $f_chr_len = $ARGV[0]; # chromosome length file for whole genome (e.g., ./data/rg_human_hg19.len)
my $dir_out = $ARGV[1]; # the files in this folder with suffixs of "fends_GC" and "indexes_mapped_uniquely"

######################## get chr len
my %chr_lens;
open IN, "$f_chr_len";
while (my $line = <IN>) {
  chomp $line; $line =~ s/^\s+//; $line =~ s/\s+$//;
  my @items = split(/\s+/, $line);
  my $chr = $items[0];
  $chr =~ s/chr//;
  $chr_lens{$chr} = $items[1];
}
close IN;

for my $chr (keys %chr_lens) {

  #################### reads mapped uniquely
  my %mapped;
  my $f_mapped = "$dir_out/$chr.indexes_mapped_uniquely";
  if (! -e $f_mapped) { next; }
  open IN, "$f_mapped" or die $!;
  while (my $line = <IN>) {
    chomp $line; $line =~ s/^\s+//; $line =~ s/\s+$//;
    $mapped{$line} = 1;
  }
  close IN;

  my $f_out = "$dir_out/$chr.fends_mappability";
  open OUT, ">$f_out";

  my $chr_len = $chr_lens{$chr};
  my %GCs;
  my $f_GC = "$dir_out/$chr.fends_GC";
  if (! -e $f_GC) { next; }
  open IN, "$f_GC" or die $!;
  while (my $line = <IN>) {
    chomp $line; $line =~ s/^\s+//; $line =~ s/\s+$//;
    my @items = split(/\s+/, $line);
    my $pos = $items[0] + 1; # the coordinate in fends_GC is 0-based.
    
    my ($map1, $map2);  
    my $before_1 = $pos - 500;
    my $before_2 = $pos - 1;
    if ($before_1 < 1) { 
      $map1 = 'NA';
    } else {
      my ($fz, $fm);
      $fz = $fm = 0;
      for (my $i = $before_1; $i <= $before_2; $i++) {
        if ($i % 10 == 1) { $fm++; }
        if (exists $mapped{$i}) { $fz++; }
      }
      $map1 = sprintf "%0.3f", $fz / $fm;
    }
    ###############################
    my $after_1 = $pos + 1;
    my $after_2 = $pos + 500;
    if ($after_2 > $chr_len) {
      $map2 = "NA";
    } else {
      my ($fz, $fm);
      $fz = $fm = 0;
      for (my $i = $after_1; $i <= $after_2; $i++) {
        if ($i % 10 == 1) { $fm++; }
        if (exists $mapped{$i}) { $fz++; }
      }
      $map2 = sprintf "%0.3f", $fz / $fm;
    }
    print OUT "$pos\t$map1\t$map2\n";
  }close IN;
  close OUT;
}

# End

