#!/usr/bin/perl

use strict;
use Benchmark;
use Data::Dumper 'Dumper';

my %methrange;
my $t0 = Benchmark->new;

open(IN1,"Methylation.gtf");
while (<IN1>) {
    chomp;
    if ( /^\s*Gm10/ ) {
        my ( $bgn, $end, $methcount ) = (split)[3,4,10];
        $methrange{$bgn}{$end}{$_} = $methcount;
    }
}
my @sorted_meth_starts = sort {$a<=>$b} keys %methrange;
my %methhash;

open(INPUT, "mc_11268_10.txt");
while(<INPUT>) {
    next unless /^(\d+)\t(\d+)\t\S\t(\w+)/;  # skip header line and any blank lines
    chomp;
    my ( $meth, $position, $class ) = ( "Gm$1", $2, $3 );

    my $indx = $#sorted_meth_starts;
    while ( $indx >= 0 ) {
        if ( $position >= $sorted_meth_starts[$indx] ) {
            my $bgn = $sorted_meth_starts[$indx];
            for my $end ( keys %{$methrange{$bgn}} ) {
                if ( $position <= $end ) {
                    for my $match ( keys %{$methrange{$bgn}{$end}} ) {
                        $methhash{$match}{methcount} = $methrange{$bgn}{$end}{$match};
                        $methhash{$match}{$class}++;
                    }
                }
            }
            last;
        }
        $indx--;
    }
}

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";

for my $ge (sort keys %methhash) {
    my ( $chh, $chg, $cg, $total_count );
    if ( $methhash{$ge}{methcount} == 0 ) {
        $chh = $chg = $cg = "INF";
    }
    else {
        my $denom = $methhash{$ge}{methcount};
        $chh = sprintf( "%.3f", $methhash{$ge}{CHH} / $denom );
        $chg = sprintf( "%.3f", $methhash{$ge}{CHG} / $denom );
        $cg  = sprintf( "%.3f", $methhash{$ge}{CG} / $denom );
        $total_count = sprintf( "%.3f", ( $methhash{$ge}{CHH}+$methhash{$ge}{CHG}+$methhash{$ge}{CG}) / $denom );
    }
    print FILEIT join( "\t", ( split /\t/, $ge )[0..4], $cg, $chg, $chh ),"\n";
}
        
my $t1 = Benchmark->new;
my $td = timediff($t1, $t0);
print "the code took:",timestr($td),"\n";

##</code><code>##


 Gm10   Glyma10g00200.1 CDS_11  8569    8705    -       2       2       6       10      14      GAAGAACCAAAATCGCTGATTTTGAAAAATAGGGGGACTGAATGTCTTGATTTGGAAATAGGGGGACCAAAATTGCGGATTGAGCATTATAGGGGACCAATAGTGTATTTTAGCCTTTATTTTATTTTGAATGTGTG
 Gm10   Glyma10g0test.1 CDS_11  8569    8705    -       2       2       6       10      14      GAAGAACCAAAATCGCTGATTTTGAAAAATAGGGGGACTGAATGTCTTGATTTGGAAATAGGGGGACCAAAATTGCGGATTGAGCATTATAGGGGACCAATAGTGTATTTTAGCCTTTATTTTATTTTGAATGTGTG
 Gm10   Glyma10g00200.1 CDS_10  8822    8822    -       0       0       0       0       0       G
 Gm10   Glyma10g00200.1 CDS_9   8885    9079    -       2       13      23      38      50      GTATTTGTGTATTGCTCCTGGATCTAGACCTTCCCCTGTTGCTGCTGCTGCTGCAAGTCACAAATTAATCACATGTATTTCATGCTTCGATGAGCGTTCACTTGCATATCATGCTGTTGGATATGGAAGAGGATCTCACATTCCAGCAGTAGCCATTACATCATCAGGCACTGCAGTTTCCAACCTTCTTCCTGC
 Gm10   Glyma10g00200.1 CDS_8   9151    9255    -       2       2       11      15      19      GTTGGATCATGCCAATGAATTATCAAACTCCCTGAAAGAAAGTGCCAATATTAACACTGTTCGGGCATCGCTTATTGTTGAAGAATGCACAAGACTTGGTTTGAT
 Gm10   Glyma10g00200.1 CDS_7   11534   11698   -       1       7       16      24      33      GGAATATACTGATCCTCACAAAGGTTCTCTTAACCTTGCTGTGAAGTTGCCTAATGGTTCTCTAAAGCCAGACCTGGGGCCAAAAACATATATTGCTTATGGATTTCTTCAGGAGCTTGGACGTGGTGATTCAGTGACTAAGCTCCATTGTGATATGTCTGATGC
 Gm10   Glyma10g00200.1 CDS_6   12777   12823   -       0       0       8       8       11      GCTTTTTACACAAATTTCATTCTTCACCTTCATTCTCATGAAGATTA
 Gm10   Glyma10g00200.1 CDS_5   12958   13022   -       1       4       4       9       10      GCTGTAGATCTTCAGTACAAGTATTTAAGGCATTTTCAGTGGCATTGGGAAAAGGGGGAACCGCT
 Gm10   Glyma10g00200.1 CDS_4   13348   13502   -       0       4       20      24      33      GAGGATCCTTATCTTCTTCTTGGGTTCTGAAAAAATAAAAACATCAACTCCCTCCTGTATTCATGCTTGGTGTTGGTTATGTGTCATGTCTGGTCTACTATTGACACATGTAAATCCAAATGCACTAGAGCTGCCCTCAAAAAACTCAATTTGGT
 Gm10   Glyma10g00200.1 CDS_3   13713   13829   -       0       2       13      15      25      GATTGAGTTAGATGAGCTGGAAAGTGTTTCCATTCTGTCAATGACCCTTGCATGGGATGAATTTTCCTTCTCTACTTTTCAAGAAGCCCATTATTCACTTCAAGATTCTCTAGACCA
 Gm10   Glyma10g00200.1 CDS_2   13915   14209   -       8       11      40      59      88      GCTGCTCTTCTCCCTCAGCTGCTCTGTCTCCGGCGTTGACATTGGAGGAGGGACTTGAGAAACTGAAGGAGGCTCTCCAAATCTTGAATTCTCCTTCCCCTTCTTCCCCTACTGGATTCCTTAGGTTTCAGGTGGCGCTCCCTCCCAGTCCTAAGACCTTCACTTTGTTTTGCTCCCAACCCCACTCCTCCTCGGTCTTTCCTCTCATTTATGTTTCCAAGAACGACGCCGACTCTAAATCACTCTATGTCAATTTGGATGATAATGTGAGTCACCGGAAAGAAAGGTTCTTTTT
 Gm10   Glyma10g00200.1 CDS_1   14242   14286   -       2       1       9       12      21      GATCTCGTTCCCCCTCACCACATAACACCCTGTCGTTCCCTTCAC
 Gm10   Glyma10g00210.1 CDS_7   15480   15722   -       4       7       26      37      47      GGTACGATCAAACACATTGGGAGCTGTTGGATTCCTTGGGGATAGCAGAAGAATCAATGTTGCCATCACAAGAGCACGCAAACATTTAGCTTTGGTCTGTGACAGCTCGACTATATGCCACAATACCTTCTTAGCAAGGCTTCTGCGTCATATTAGACACTTTGGTAGGGTGAAGCATGCAGAACCAGGTAGTTTTGGAGGATATGGACTTGGGATGAATCCAATATTACCTTCCATTAATTA
 Gm10   Glyma10g00210.1 CDS_6   16625   16782   -       1       9       15      25      32      GGTGTTAGCCCAACAGCTATTGCAGTGCAATCCCCTTATGTTGCTCAAGTACAACTTTTGAGGGACAAGCTTGATGAATTTCCAGAAGCAGCAGGTACTGAGGTTGCAACCATTGACAGTTTTCAAGGTCGGGAAGCTGATGCAGTAATTTTATCCAT
 Gm10   Glyma10g00210.1 CDS_5   17595   17763   -       2       9       16      27      34      GCCTACTTGGATAACACAATGCCCGCTGCTATTGCTAGATACTAGAATGCCATATGGAAGTCTGTCAGTTGGTTGTGAAGAGCATCTAGACCCGGCTGGAACAGGCTCACTTTATAATGAAGGAGAAGCTGAGATAGTTTTGCAGCATGTATTTTCCTTAATCTATGCC
 Gm10   Glyma10g00210.1 CDS_4   18046   19077   -       12      42      91      145     173     GCGCAACTGTGAAGCTTTAATGCTGCTTCAGAAGAATGGTTTACGAAAGAAGAATCCTTCAATTTCTGTTGTTGCTACACTGTTTGGAGATGGGGAAGATGTTGCATGGCTTGAGAAAAATCATTTGGCTGACTGGGCAGAAGAAAAATTGGATGGAAGATTAGGAAATGAAACCTTTGATGATTCTCAGTGGAGAGCAATTGCAATGGGTTTGAATAAAAAGAGGCCTGTATTGGTTATCCAAGGCCCTCCTGGTACAGGCAAGACTGGTTTGCTCAAGCAACTTATAGCATGTGCTGTTCAGCAAGGTGAAAGGGTTCTTGTTACAGCACCTACTAATGCAGCTGTTGATAACATGGTAGAAAAGCTTTCAAATGTTGGATTAAATATAGTGCGGGTTGGAAATCCAGCTCGTATATCAAAAACAGTGGGATCAAAGTCTTTGGAAGAAATTGTAAATGCTAAGCTTGCAAGTTTTCGAGAAGAGTATGAGAGGAAGAAGTCAGATCTAAGAAAAGATCTAAGACATTGTTTAAGGGATGATTCACTAGCTTCAGGCATACGCCAACTTCTGAAGCAACTGGGAAGGTCACTGAAGAAAAAGGAAAAGCAGACCGTAATTGAAGTTCTGTCTAGTGCTCAAGTTGTGGTTGCCACTAATACTGGAGCAGCTGACCCTTTGGTTCGAAGGCTAGATACCTTTGATTTGGTTGTCATAGATGAAGCGGGACAGGCAATTGAACCCTCTTGCTGGATTCCTATATTGCAGGGAAAGCGCTGCATTCTTGCTGGTGATCAATGCCAACTTGCTCCTGTCATATTATCTAGAAAGGCCTTAGAAGTTGGTCTAGGAATATCTCTACTGGAGAGAGCTGCAACTTTGCATGAAGGGATTCTCACCACTAGGTTAACAACACAATACCGTATGAATGATGCAATTGCTAGTTGGGCTTCAAAGGAGATGTACGGAGGATTATTGAAGTCCTCTGAGACTGTCTTTTCTCATCTTCTAGTAGACTCTCCTTTTGTTAA
 Gm10   Glyma10g00210.1 CDS_3   21109   21392   -       4       8       30      42      52      GGACTTGGGGGAATGCATTTGGTGTTATTCAAGGTTGAAGGGAACCACCGGTTACCACCAACCACTCTTTCACCTGGAGACATGGTTTGTGTGAGAACATATGACAGCATGGGTGCAATTACAACTTCTTGCATACAAGGATTTGTGAACAGTTTTGGGGATGATGGCTATAGTATTACAGTTGCTTTAGAGTCACGCCATGGTGATCCTACATTCTCTAAACTATTTGGGAAGAGTGTGCGCATTGACCGTATTCAAGGACTGGCTGATACACTTACTTATGA
 Gm10   Glyma10g00210.1 CDS_2   24354   24662   -       9       13      32      54      66      GCTCAGCATAGGGCGGTTGTGAGAAAGATAACGCAACCCAAGAGTGTCCAAGGCGTTTTAGGAATGGACTTCGAAAAGGTCAAAGCATTACAGCACAGGATTGACGAATTCACCACCCATATGTCAGAACTACTCCGTATTGAAAGGGATGCTGAGTTGGAGTTTACTCAGGAGGAATTGGATGCTGTTCCTAAACCAGATGATACTTCTGATTCTTCAAAAACGATTGATTTCTTGGTTAGCCATAGCCAGCCTCAACAAGAACTCTGCGACACCATTTGTAATTTAAACGCTATCAGTACCTCTACA

##</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>##

Chromosome      Gene    Feature Start   END     CG_count_by_C   CHG_count_by_C  CHH_count_by_C
 Gm10   Glyma10g00200.1 CDS_11  8569    8705    0.143   0.071   0.000
 Gm10   Glyma10g00200.1 CDS_9   8885    9079    0.040   0.000   0.000
 Gm10   Glyma10g0test.1 CDS_11  8569    8705    0.143   0.071   0.000