#!/usr/bin/perl -w
use strict;
# A perl script for converting a methylFreq file to a UCSC BED file. 
# USAGE: ./methylFreq2BED.pl sample_methylFreq.txt > sample_methylFreq.BED.txt
#   sample_methylFreq.txt: the methylFreq file generated by toMethylFreq.pl
# Format of the BED file:
# Column 1: chromosome
# Column 2: chromosome_start
# Column 3: chromosome_end
# Column 4: methylation level [0-1]
# column 5: sequencing depth
# column 6: strand, please ignored since the reads from both strands were combined
# Column 7: chromosome_start
# Column 8: chromosome_end
# Column 9: pesudo color, red is methylated, and green is unmethylated
#
# Written by Kun Zhang (kzhang@bioeng.ucsd.edu), Aug. 2008
#

my $minDepth = 10;
my $methylFreqFile = $ARGV[0];

my @palette=("0,240,0", "30,210,0","60,180,0","90,150,0","120,120,0","150,90,0","180,60,0","210,0,0");
my %methylTable;

my @words = split(/[_\.]/, $methylFreqFile);
print "track name=\"", $words[0], "\" description=\"Methylation level\" visibility=2 useScore=1 itemRgb=\"On\"\n";
open(INFILE,"$methylFreqFile") || die("Error in opening file $methylFreqFile!\n");
while(my $line = <INFILE>){
	chop($line);
	my @fields = split(/\t/, $line);
	#next if($fields[3]<$minDepth);
	my ($chr,$chr_start, $chr_end) = split(/[:\-]/, $fields[0]);
	my $target_offset = $fields[1];
	my $strand = $fields[2] eq 'W' ? '+' : '-';
	my %alleleFreq;
	my $Fct;
	for(my $i=4; $i<scalar(@fields); $i+=2){
		$alleleFreq{$fields[$i]}=$fields[$i+1];			
		$Fct += $fields[$i+1] if($fields[$i]=~ /[CT]/);
	}
	next if (!$Fct || $Fct < 0.9);
	my $chr_pos = $chr_start+$target_offset;
	my $index="$chr:$chr_pos";
	$alleleFreq{'C'} =0 if(!$alleleFreq{'C'});
	$methylTable{$index}->{'C'}	=  $methylTable{$index}->{'C'} ? 
				$methylTable{$index}->{'C'} + $fields[3]*$alleleFreq{'C'} : 
				$fields[3]*$alleleFreq{'C'};
	$methylTable{$index}->{'CT'} += $methylTable{$index}->{'CT'} ? 
				$methylTable{$index}->{'CT'} + $fields[3]*$Fct : 
				$fields[3]*$Fct;
}
close(INFILE);

foreach my $index(keys(%methylTable)){
	next if($methylTable{$index}->{'CT'}<$minDepth);
	my ($chr,$chr_pos) = split(/:/, $index);	
	my $methylLevel = sprintf("%4.3f", $methylTable{$index}->{'C'} ? $methylTable{$index}->{'C'}/$methylTable{$index}->{'CT'} : 0);
	print "$chr\t", 
			$chr_pos-1, "\t", 
			$chr_pos, "\t",
			"$methylLevel\t",
			int($methylTable{$index}->{'CT'}), "\t", 
			"+\t",
			$chr_pos-1, "\t",
            $chr_pos, "\t",
			$palette[int($methylLevel*8-0.0001)],"\n";
}
