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

my $num_chrs = $ARGV[0]; # Human 23 and Mouse 20
my $dir_rg_fa = $ARGV[1]; # directory of reference genome(e.g., chr1.fa) 
my $dir_out = $ARGV[2]; # output directory

my $cutting_code = "GATC"; # for dpnii enzyme
my $fends_len = 200; #length of fragment length for GC
$cutting_code = uc $cutting_code;
my $cutting_code_len = length($cutting_code);


if (! -d $dir_out) { `mkdir $dir_out`; }

for (my $c = 1; $c <= $num_chrs; $c++) {
  my $chr = $c;
  if ($c == $num_chrs) { $chr = 'X'; }
  my $f_fa = "$dir_rg_fa/chr$chr.fa";

  open IN, "$f_fa" or die $!;
  my $n = 0;
  my $seq;
  while (my $line = <IN>) {
    $n++; 
    if ($n == 1) { next; }
    chomp $line; $line =~ s/^\s+//; $line =~ s/\s+$//;
    $seq .= $line;
  }
  close IN;

  my $len = length($seq);
  my $f_out = "$dir_out/$chr.fends_GC";
  open OUT, ">$f_out";

  for (my $i = 0; $i < ($len - $cutting_code_len -1); $i++) {
    my $target = uc substr($seq, $i, $cutting_code_len);
    if ($target eq 'GATC') {
      # [i-200, i-1]
      my $before_1 = $i - $fends_len;
      my $before_2 = $i - 1;
      my $gc1 = 0;
      if ($before_1 < 0) {
        $gc1 = 'NA';
      } else {
        my $seq1 = substr($seq, $before_1, $fends_len);
        $gc1 = sprintf "%0.3f", get_GC($seq1);
      }
      # [i+1, i+200]
      my $after_1 = $i + 1;
      my $after_2 = $i + $fends_len;
      my $gc2 = 0;
      if ($after_2 >= $len) {
        $gc2 = 'NA';
      } else {
        my $seq2 = substr($seq, $after_1, $fends_len);
        $gc2 = sprintf "%0.3f", get_GC($seq2);
      }
      print OUT "$i\t$gc1\t$gc2\n";
    }
  }
  close OUT;
}

sub get_GC{
  my $seq = shift;
  $seq = uc $seq;
  my $seq_len = length($seq);
  my $num_GC = 0;
  for my $c (split //, $seq) {
    if ($c eq 'G' || $c eq 'C') { $num_GC++; }
  }
  $num_GC /= $seq_len;

  return $num_GC;
}


# End

