#!/usr/bin/perl -w
use strict;

#  A perl script for extracting the base compositions on CpG sites from the mapped sequencing reads
# USAGE: ./toMethylFreq.pl <sample_mapped.txt > sample_methylFreq.txt
#   sample_mapped.txt: the mapped sequencing reads generated by MapBisReadsTargetsSpeedy.pl
#
# Written by Kun Zhang (kzhang@bioeng.ucsd.edu), Aug. 2008;
# Last modified: 06/30/2009
#


my %methylFreqTable;
my $offset_2;
while(my $line = <STDIN>){
	chop($line);
	my @fields = split(/\t/, $line);
	for(my $i=4; $i<scalar(@fields); $i+=2){
		$methylFreqTable{$fields[1]}->{$fields[$i]}->{$fields[2]}->{$fields[$i+1]}++;
		if($fields[2] eq "C"){
			$offset_2 = $fields[3]-$fields[$i]-1;
		}else{
			$offset_2 = $fields[$i]-$fields[3]+1;
		}
		#$methylFreqTable{$fields[1]}->{$fields[$i]}->{$fields[2]}->{'2'}->{substr($fields[0],$offset_2,2)}++;
	}
}

foreach my $index (sort keys(%methylFreqTable)){
	foreach my $pos (sort keys(%{$methylFreqTable{$index}})){
		foreach my $strand (sort keys(%{$methylFreqTable{$index}->{$pos}})){
			print "$index\t$pos\t$strand";
			my $total =0;
			foreach my $allele (sort keys(%{$methylFreqTable{$index}->{$pos}->{$strand}})){
				$total += $methylFreqTable{$index}->{$pos}->{$strand}->{$allele};
			}
			print "\t$total";
			foreach my $allele (sort keys(%{$methylFreqTable{$index}->{$pos}->{$strand}})){
				print "\t$allele\t", $methylFreqTable{$index}->{$pos}->{$strand}->{$allele}/$total;
			}
			print "\n";
		}
	}
}
