#!/usr/bin/perl
=head2 SYNOPSIS

	Description: Considering that bismark reports the methylation in both the strands,
				 the methylation info for the symmetrical CHG sites were merged.
             
=head2 AUTHOR

	Kang Zhang  Email:  demut@foxmail.com
	
=cut

open FA,"<$ARGV[0]" or die;
open IN,"<$ARGV[1]" or die;

while(<FA>){
	chomp;
	if(/^>(\S+)/){
		$id=$1;
	}else{
		$seq{$id}.=$_;
	}
}


while(<IN>){
	if(/^</){
		chomp;
		print "$_\tTotal:Max-[(Pos1-Depth):(Pos2-Depth)]\tMethy:Max-[(Pos1-Depth):(Pos2-Depth)]\n";
		next;}
	chomp;
	my @arr=split;
	my $pos=$arr[1]-1;
	my $site=substr($seq{$arr[0]},$pos,1);
	if($site eq "G" or $site eq "g"){
		my $p=$arr[1]-2;
		$countAll{$arr[0]}{$p}+=$arr[2];
		$count{$arr[0]}{$p}+=$arr[3];
		$countAllmax{$arr[0]}{$p}{2}=$arr[2];
		$countmax{$arr[0]}{$p}{2}=$arr[3];
	}elsif($site eq "C" or $site eq "c"){
		$countAll{$arr[0]}{$arr[1]}+=$arr[2];
		$count{$arr[0]}{$arr[1]}+=$arr[3];
		$countAllmax{$arr[0]}{$arr[1]}{1}=$arr[2];
		$countmax{$arr[0]}{$arr[1]}{1}=$arr[3];
	}else{
		print STDERR "Unknown site! $site...\n";
		next;
	}
}

foreach my $chr (sort keys %countAll){
	foreach my $pos (sort {$a <=> $b} keys %{$countAll{$chr}}){
		my ($a , $b)=($countAllmax{$chr}{$pos}{1},$countAllmax{$chr}{$pos}{2});
		my $max="-";
		if($a >=$b){
			$max=$a;
		}else{$max=$b;}
		if(!defined $countmax{$chr}{$pos}{1}){$countmax{$chr}{$pos}{1}=0;}
		if(!defined $countmax{$chr}{$pos}{2}){$countmax{$chr}{$pos}{2}=0;}
		my ($x , $y)=($countmax{$chr}{$pos}{1},$countmax{$chr}{$pos}{2});
		my $maxx="-";
		if($x >=$y){
			$maxx=$x;
		}else{$maxx=$y;}
		if(!defined $count{$chr}{$pos}){
			$count{$chr}{$pos}=0;
		}
		print "$chr\t$pos\t$countAll{$chr}{$pos}\t$count{$chr}{$pos}\t$max-($a:$b)\t$maxx-($x:$y)\n";
	}
}
