use strict;


my $chr=$ARGV[0];

sub rc {
    my $str = $_[0];
    # reverse complement
    $str =~ tr/ACGTacgt/TGCAtgca/;
    return scalar reverse $str;
}


sub rvc{
    my $str = $_[0];
    my @a=split/\_/,$str;
    my $rv=&rc($a[0]);
    my $rd=&rc($a[1]);
    my $out="$rv".'_'."$rd";
}



####collect the data; Table from Seplyarskiy et al 2021
open (IN,"/net/home/vseplyarsky/TOPMEDs/tracks_to_correlate/signatures_last5/components_intensities_10kb.txt.mod.epi")||die;
my %rt; my %rec; my %H3k4me3; my %H3k4me1; my $all; my %all;
<IN>;
while (<IN>){
    chomp; my @a=split/\s/;
    $all++;
    ${$all{'rt'}}{$a[-4]}++; ### replication time
    ${$all{'rec'}}{$a[-5]}++; ### Recombination
    ${$all{'me3'}}{$a[-10]}++; ### H3k4me3
    ${$all{'me1'}}{$a[-11]}++; ### H3k4me1
}

#### epigenetic bins boundares
my %bin; my $bi=5;
foreach my $tr (sort keys %all ){
    my $lb=1; my $ls;
    foreach my $v (sort {$a <=> $b} keys %{$all{$tr}}){
	$ls+=${$all{$tr}}{$v};
#	print "$tr $v $lb $ls\n";
	if ($lb*($all/$bi)<=$ls && $lb<5){${$bin{$tr}}{$lb}=$v;$lb++;}}}

my %tr;

my %conv;

$conv{'rt'}=-4;
$conv{'rec'}=-5;
$conv{'me3'}=-10;
$conv{'me1'}=-11;


####atribute every site to epegenetic bin data; Table from Seplyarskiy et al 2021 
open (IN,"/net/home/vseplyarsky/TOPMEDs/tracks_to_correlate/signatures_last5/components_intensities_10kb.txt.mod.epi")||die;
<IN>;
while (<IN>){
    chomp; my @a=split/\s/;
    next unless $a[0]==$chr;
    foreach my $tr (sort keys %conv){
	my $b=0;
	my $p=$conv{$tr};
	for my $i (1..$bi-1){
	    if ($a[$p]>=${$bin{$tr}}{$i}){$b=$i}
	}
#	print "$tr $a[1] $b =${$bin{$tr}}{'1'} ${$bin{$tr}}{'2'}\n";
	${$bin{$tr}}{$a[1]}=$b;

    }
}




my %pol;
my %exep;
### exclude high frequency SNVs
open (pol,"/net/home/vseplyarsky/GNOMAD_Tufts/results/$chr\_common.pol")||die;
while (<pol>){
    chomp;
    my @a=split/\s/;
    ${$exep{$a[0]}}{$a[1]}=1;
}

### GNOMAD rare SNVs
open (pol,"/net/home/vseplyarsky/GNOMAD_Tufts/results/$chr.pol")||die;
while (<pol>){
    chomp;
    my @a=split/\s/;
    ${$pol{$a[0]}}{$a[1]}=1;
}


my %obs; my %exp; my %site;
#### synonimouse sites
    open (IN,"gunzip -c /net/home/vseplyarsky/GNOMAD_Tufts/script/updates/data/$chr\_syn_full_table.gz|");
    while (<IN>){
	chomp; my @a=split/\s/;
#	print "$a[0] $exp{$a[0]}\n";
#	next unless $exp{$a[0]}>=$b4;

	my $mu=$a[2];

	substr ($mu,0,1)='';
        substr ($mu,1,1)='';
	next if ${$exep{$a[1]}}{$mu}; 
	my $pol=0;
	if ( ${$pol{$a[1]}}{$mu}){$pol=1}
	my $v=1-exp(-$a[3]);
#	print "$v $a[3]\n";
        my $mu=$a[2];
	my $an=substr($a[2],1,1);
	if ($an ne 'C' && $an ne 'T'){$mu=&rvc($mu); $a[2]=&rvc($a[2])}
#	if ($a[2]=~/CG/){$mu=substr($a[2],1)}
	$a[2]=$mu;
	my $c=int($a[1]/10000)*10000;
	foreach my $tr (sort keys %conv){
	    my $b=   ${$bin{$tr}}{$c};
#	    print "$b $tr $c $a[2] $v $pol\n";
	    ${${$obs{$tr}}{$a[2]}}{$b}+=$pol;
	    ${${$exp{$tr}}{$a[2]}}{$b}+=$v;
	    ${${$site{$tr}}{$a[2]}}{$b}++;
	}
}



#### print table with rates
foreach my $tr (sort keys %site){
    open (OUT,">/net/home/vseplyarsky/GNOMAD_Tufts/script/updates/data/$chr\_$tr\_binning");
    foreach my $mu (sort keys %{$site{$tr}}){
	for my $b (0..$bi-1){
	    my $r=${${$obs{$tr}}{$mu}}{$b}/${${$site{$tr}}{$mu}}{$b};
            my $r1=${${$exp{$tr}}{$mu}}{$b}/${${$site{$tr}}{$mu}}{$b};
         

print OUT "$mu $b $r $r1 ${${$obs{$tr}}{$mu}}{$b} ${${$exp{$tr}}{$mu}}{$b} ${${$site{$tr}}{$mu}}{$b}\n"
	}}}
