# input shoud be paired ends soap alignment merged W and C strands#

if(@ARGV<3)
{
	warn "usage: $0 <hg19.W.CG.bed> <hg19.C.CG.bed> <*.rmdup.bsalign>\n";
	exit 1;
}


my $W_CG_bed=shift;
my $C_CG_bed=shift;
open FL,"$W_CG_bed";
my %h1;
while(<FL>)
{
	chomp;
	my @sp=split;
	$h1{$sp[0]."\t".$sp[2]}="";
}
close FL;

open FL,"$C_CG_bed";
my %h2;
while(<FL>)
{       
	chomp;
	my @sp=split;
	$h2{$sp[0]."\t".$sp[2]}="";
}
close FL;

my %ll;
open IN,"$soap_bsalign";
my ($meth_nu_r1,$meth_nu_r2,$unmeth_nu_r1, $unmeth_nu_r2)=(0,0,0,0);
while(<>)
{
	my $line1=$_;
	chomp $line1;
	my $line2=<>;
	chomp $line2;
	my @sp1=split /\t/,$line1;
	my @sp2=split /\t/,$line2;
	my $start=$sp1[8];
	my $end=$sp2[8]+length($sp2[1])-1;
	my $length=$end-$start+1;
	my $len1=30;
	my $len2=30;
	my @base1=split //,$sp1[1];
	my @base2=split //,$sp2[1];
	my $count1=0;
	my $count2=length($sp2[1])-30;
	if($sp1[$#sp1] eq "W")
	{
		for my $site($sp1[8]..$sp1[8]+29)
		{
			if(exists $h1{$sp1[7]."\t".$site})
			{
				my $bas=$base1[$count1];
				if($bas eq "C" || $bas eq "c")
				{
					$meth_nu_r1++;
				}
				elsif($bas eq "T"||$bas eq "t")
				{
					$unmeth_nu_r1++;
				}
			}
			$count1++;
		}

		for my $site2($sp2[8]+length($sp2[1])-30..$sp2[8]+length($sp2[1])-1)
		{       
			if(exists $h1{$sp2[7]."\t".$site2})
			{       
				my $bas=$base2[$count2];
				if($bas eq "C" || $bas eq "c")
				{       
					$meth_nu_r2++;
				}
				elsif($bas eq "T"||$bas eq "t")
				{       
					$unmeth_nu_r2++;
				}
			}
			$count2++;
		}

	}
	elsif($sp1[$#sp1] eq "C")
	{
		for my $site($sp1[8]..$sp1[8]+29)
		{       
			if(exists $h2{$sp1[7]."\t".$site})
			{       
				my $bas=$base1[$count1];
				if($bas eq "C" || $bas eq "c")
				{
					$meth_nu_r1++;
				}
				elsif($bas eq "T"||$bas eq "t")
				{
					$unmeth_nu_r1++;
				}
			}
			$count1++;
		}

		for my $site2($sp2[8]+length($sp2[1])-30..$sp2[8]+length($sp2[1])-1)
		{
			if(exists $h2{$sp2[7]."\t".$site2})
			{
				my $bas=$base2[$count2];
				if($bas eq "C" || $bas eq "c")
				{
					$meth_nu_r2++;
				}
				elsif($bas eq "T"||$bas eq "t")
				{
					$unmeth_nu_r2++;
				}
			}
			$count2++;
		}

	}
}

my $meth1=$meth_nu_r1/($meth_nu_r1+$unmeth_nu_r1);
my $meth2=$meth_nu_r2/($meth_nu_r2+$unmeth_nu_r2);
my $jag_index=($meth1-$meth2)*100/($meth1);
print "$jag_index\n";
