#	Chromosome	Gene	Feature	Start	END	Strand	CG_count	CHG_Count	CHH_count	Total_meth_count	C_count	Sequence
#Gm10	Glyma13g27990.1	CDS_1	31097691	31097752	+	1	1	6	8	9	TTAAGGAAATTATAAATCAAAAAATCAAATTTGACTCGTCACTGGAGTACTTGTCCATTTGG

##</code><code>##

assembly	position	strand	class	mc	h
10	8691	+	CG	12	12
10	8692	-	CG	3	3
10	8705	+	CHG	7	18
10	8707	-	CHG	2	4
10	8717	+	CG	25	25
10	8718	-	CG	6	6
10	8984	+	CG	7	7
10	8991	+	CG	6	9

##</code><code>##

 
%methhash = ();
use Benchmark;
    $t0 = Benchmark->new;
use lib qw(.);
use Tie::Autotie 'Tie::IxHash';
tie %methhash, "Tie::IxHash";
open(IN1,"Methylation.gtf");
@genome = grep{/^Gm10/} <IN1>;
open(INPUT, "mc_11268_10.txt");
while($line = <INPUT>){
chomp($line);
$count++;
if($count == 1){
next;
}
else{
@meth = split("\t",$line);
 $meth[0] = "Gm".$meth[0];
foreach $genes(@genome){
	chomp($genes);
	@gene_line = split("\t",$genes);
	if($gene_line[0] eq $meth[0]){
	$start = $gene_line[3] ;
	$end = $gene_line[4];
		if(($meth[1]>=$start)&&($meth[1]<=$end)){
				if($meth[3] eq "CHH"){
			${$methhash{$genes}}[0]+=1;
						}
			elsif($meth[3] eq "CHG"){
				${$methhash{$genes}}[1]+=1;
			}
			elsif($meth[3] eq "CG"){
				${$methhash{$genes}}[2]+=1;
					}
		}
		
		else{next;}
		}
		else{next;}
		@gene_line = ();}
		@meth=();
		}
		}
		open(FILEIT,">res_11268_10.txt");
print FILEIT join("\t",qw(Chromosome Gene Feature Start END CG_count_by_C CHG_count_by_C CHH_count_by_C)),"\n";
foreach $ge (sort keys %methhash) {
	@final = split("\t",$ge);	
	if($final[10]==0){
	$chh = $chg = $cg = "INF";
	}
	else{
	$chh = sprintf("%.3f",(${$methhash{$ge}}[0])/$final[10]);
	$chg = sprintf("%.3f",(${$methhash{$ge}}[1])/$final[10]);
	$cg = sprintf("%.3f",(${$methhash{$ge}}[2])/$final[10]);
    $total_count = sprintf("%.3f",((${$methhash{$ge}}[2])+(${$methhash{$ge}}[1])+(${$methhash{$ge}}[0]))/$final[10]);
	}
	print FILEIT join("\t",(@final[0..4],$cg,$chg,$chh)),"\n";
	$chh = $chg=$cg=$total_count = '';
}
	
	$t1 = Benchmark->new;
$td = timediff($t1, $t0);
print "the code took:",timestr($td),"\n";