#!/usr/common/usg/languages/perl/5.16.0/bin/perl -w
#
# kmerHistAnalyzer.pl by Jarrod Chapman <jchapman@lbl.gov> Tue Sep 20 04:01:31 PDT 2011
# Copyright 2011 Jarrod Chapman. All rights reserved.
#

use Getopt::Std;
my %opts = ();
my $validLine = getopts('h:SP:', \%opts);
my @required = ("h",);
my $nRequired = 0;
map {$nRequired += exists($opts{$_})} @required;
$validLine &= ($nRequired == @required);
if ($validLine != 1) {
    print "Usage: ./kmerHistAnalyzer.pl <-h histFile (2 column format expected)> <<-S(econdaryPeak?)>> <<-P setPeak>>\n";
    exit;
}

my $histFile = $opts{"h"};

my $setPeak = 0;
if (exists($opts{"P"})) {
    $setPeak = $opts{"P"};
}

my $useSecondary = 0;
if (exists($opts{"S"})) {
    $useSecondary = 1;
}

open (H,$histFile) || die "Couldn't open $histFile\n";

my %hist = ();
while (my $line = <H>) {
    chomp $line;
    my @cols = split(/\s+/,$line);
    $hist{$cols[0]} = $cols[1];
}
close H;

my @freqs = sort {$a <=> $b} (keys(%hist));
my $nFreqs = scalar(@freqs);

my @localMaxima = ();
my $globalMaximum = 0;
my $peakFreq = 0;
for (my $i = 1; $i < $nFreqs-1; $i++) {
    if (($hist{$freqs[$i]} > $hist{$freqs[$i+1]}) && ($hist{$freqs[$i]} > $hist{$freqs[$i-1]})) {
	push(@localMaxima,$freqs[$i]);
#	print "$freqs[$i]\t$hist{$freqs[$i]}\n";
	if ($hist{$freqs[$i]} > $globalMaximum) {
	    $peakFreq = $freqs[$i];
	    $globalMaximum = $hist{$peakFreq};
	}
    }
}

if ($setPeak) {
    $peakFreq = $setPeak;
    $globalMaximum = $hist{$peakFreq};
} elsif ($useSecondary) {
    $peakFreq = $localMaxima[1];
    $globalMaximum = $hist{$peakFreq};
}

print "# Peak depth:  $peakFreq ($globalMaximum)\n";

my $minFreq = 0.25*$peakFreq;

my $genomeEstimate = 0;

while (my ($f,$c) = each(%hist)) {
    unless ($f > $minFreq) {
	next;
    }

    $genomeEstimate += ($c/2)*($f/$peakFreq);
}

my $genomeMB = sprintf("%.2f",$genomeEstimate/1e6);
print "# Estimated genome size: $genomeMB MB\n";

my $cumul = 0;
my @markers = (0.5,1,2,5,10,20,50,100,200,500,1000,2000,5000,10000,20000,50000,100000);
my $currentMarker = 0;
my $nMarkers = scalar(@markers);
my $keepMarking = 1;

foreach my $f (@freqs) {
    unless ($f > $minFreq) {
	next;
    }

    my $c = $hist{$f};

    my $normF = $f/$peakFreq;

    if ($keepMarking && ($normF > $markers[$currentMarker])) {
	my $lastpG = sprintf("%.2f",100*$cumul/$genomeEstimate);
	print "# $lastpG% <= $markers[$currentMarker]X\n";
	$currentMarker++;
	if ($currentMarker >= $nMarkers) {
	    $keepMarking = 0;
	}
    }

    $cumul += ($c/2)*($f/$peakFreq);
    my $pG = $cumul/$genomeEstimate;

    print "$normF\t$pG\n";


}
