#!/usr/bin/perl -w
# Hapinfo to methylation haplotype load (MHL)
# Run the script to the Hapinfo directory
# Contact: Kun Zhang
# Version 1.3
# Update: 2016-02-29

use strict;
use Cwd;
die &USAGE if @ARGV <1;
my %mch_load_matrix;
my %probe_HMH_samples;
my %hap_count_matrix;
my $hapinfList=shift @ARGV;
open FF,$hapinfList;
chomp(my @hapInfo_files=<FF>);
close FF;

# making sure this region have at least 3 CpGs coverage and at least 3 reads coverage
my $min_hap_count = 3;
my $min_hap_length = 1;

my @sample_list;
foreach my $hapInfo_file(sort @hapInfo_files){
        my @line=split /\//,$hapInfo_file;
	my $sample_name = $line[$#line];
	$sample_name =~ s/.hapInfo.txt//;
	push(@sample_list, $sample_name);
	open(INFILE, "$hapInfo_file") || die("Error in opening $hapInfo_file!");
	while(my $line = <INFILE>){
		chop($line);
		my @fields = split(/\t/, $line);
		next if(scalar(@fields)<4);
		my $probeID = $fields[0];
		my $hapString = $fields[1];
		next if(length($hapString)<$min_hap_length);		
		$hap_count_matrix{$probeID}->{$sample_name}->{$hapString}+=$fields[2];
	}
	close(INFILE);
}

my @unmethylated_haps= ("T", "TT", "TTT", "TTTT", "TTTTT","TTTTTT","TTTTTTT","TTTTTTTT","TTTTTTTTT");
my @methylated_haps  = ("C", "CC", "CCC", "CCCC", "CCCCC","CCCCCC","CCCCCCC","CCCCCCCC","CCCCCCCCC");

print "Probe_id\t", join("\t", sort @sample_list), "\n";
foreach my $probeID (keys(%hap_count_matrix)){
	print "$probeID";
	foreach my $sample_name (sort @sample_list){
		my $mc_total=0;
		my $ct_total=0;
		my $total_hap_counts=0;
		foreach my $hapString (keys(%{$hap_count_matrix{$probeID}->{$sample_name}})){
			$total_hap_counts+=$hap_count_matrix{$probeID}->{$sample_name}->{$hapString};
			for(my $i = 0; $i < length($hapString); $i++){
				my $sub_hapString = substr($hapString,$i,1);
				next if($sub_hapString =~ /[NAG]/i);
				$mc_total+=$hap_count_matrix{$probeID}->{$sample_name}->{$hapString} if($sub_hapString eq "C");
				$ct_total+=$hap_count_matrix{$probeID}->{$sample_name}->{$hapString};
			}
		}
		if($total_hap_counts < $min_hap_count){
			print "\tNA";
		}else{
			print "\t", $mc_total/$ct_total;
		}
	}
	print "\n";
}

sub USAGE{
print "\nperl $0 Hapinfo_File_list > Ouput.txt\n";
print "Just use: ls *hapInfo.txt > Hapinfo_File_list to Get Hapinfo_File_list\n";
}
