#!/usr/bin/perl
#########################################################################################
#  This script is used to combine all of the normalized_read_count file together, and  	#
#  then do the methylation enrichment analysis for each TE across all the samples using #
#  the following equation:								#
#	DNA methylation enrichment = log2(obs/exp) = log2(read_density/baseline)	#
#  Usage:										#
#	perl MethylationEnrichmentAnalysis.pl <sample type file> <s1.normalized.rp ...  #
#		sn.normalized.rp>							# 
#  Format of sample type file:								#
#	sampleName	sampleType							#
#	for example:									#
#	Medip_H1EsB1_2lanes.all.normalized.rp	H1ES					#
#	Medip_H1EsB2_2lanes.all.normalized.rp	H1ES					#
#	Medip_HuAC_2lanes.all.normalized.rp	Cortex					#
#	Medip_HuAG_2lanes.all.normalized.rp	GE					#
#	Medip_HuAT_2lanes.all.normalized.rp	FetalBrain				#
#											#
#########################################################################################

use warnings;
use strict;

if(@ARGV < 2)
{
	die "Usage: perl $0 <sample_type> <s1.normalized.rp ... sn.normalized.rp>";
}


my $sampletype = shift@ARGV;
my @files = @ARGV;
my %weight;
my %numtype;
my %type;
my $totaltype;

open(ST, $sampletype) ||die "Couldn't open the sample type file $sampletype";
while(<ST>)
{
	chomp;
	my @info = split/\t/, $_;
	$type{$info[0]} = $info[1];
	$numtype{$info[1]} ++;
}

foreach my $sample(keys%type)
{
	my $tmp = $type{$sample};
	$weight{$sample} = $numtype{$tmp};
} 	
foreach my $eachtype(keys %numtype)
{
	$totaltype++;
}


my %methyl;
foreach my $file(@files)
{
	chomp $file;
	open(IN, $file) ||die "Couldn't open $file\n";
	while(my $line = <IN>)
	{
		chomp $line;
		my @tmp = split/\t/, $line;
		my $rp = $tmp[0];
		my $value = $tmp[2];
		push @{$methyl{$rp}}, $value;
	}
	close IN;
}

print "repeat_family @files\n";
foreach my $key (sort{$a cmp $b} keys %methyl)
{
	print "$key";
	my $mean;
	for(my $i = 0; $i<= $#files; $i++)
	{
		$mean += $methyl{$key}[$i]/$weight{$files[$i]};
	}
	$mean = $mean/$totaltype;

	for(my $j = 0; $j <= $#files; $j++)
	{
		my $ratio = $methyl{$key}[$j]/$mean;
		my $lg =log($ratio)/log(2);
		print " $lg";
	}
	print "\n";
}
