#!/usr/common/usg/languages/perl/5.16.0/bin/perl -w
#
# prefixer.pl by Jarrod Chapman <jchapman@lbl.gov> Fri Aug 19 10:59:53 PDT 2011
# Copyright 2011 Jarrod Chapman. All rights reserved.
#

use Getopt::Std;
my %opts = ();
my $validLine = getopts('s:p:m:FGQ:', \%opts);
my @required = ("s","p","m");
my $nRequired = 0;
map {$nRequired += exists($opts{$_})} @required;
$validLine &= ($nRequired == @required);
if ($validLine != 1) {
    print "Usage: ./prefixer2.pl <-p minPrefixLen:maxPrefixLen> <-s seqFile> <-m merSize> <<-F(astQ?)>> <<-G(zipped?)>> <<-Q qualOffset>>\n";
    exit;
}

my ($minPrefixLen,$maxPrefixLen) = $opts{"p"} =~ /^(\d+)\:(\d+)$/;
my $seqFile = $opts{"s"};
my $merSize = $opts{"m"};

my $fastQ = 0;
if (exists($opts{"F"})) {
    $fastQ = 1;
}
my $gzipped = 0;
if (exists($opts{"G"})) {
    $gzipped = 1;
}

my $qualOffset = 64;
if (exists($opts{"Q"})) {
    $qualOffset = $opts{"Q"};
}

my $q2sym = chr(2+$qualOffset);

my $maxSamples = 100000;
my $sampleRate = 1000;
my $seq = "";
my $qualBuffer = "";

if ($gzipped) {
    open (S,"gunzip -c $seqFile |") || die "Couldn't open pipe to $seqFile\n";
} else {
    open (S,$seqFile) || die "Couldn't open $seqFile\n";
}

my $nSamples = 0;
while (my $line = <S>) {

    my $nts = undef;
    my $quals = undef;
    if ($fastQ) {
	$line = <S>;
	chomp $line;
	$nts = $line;
	$line = <S>;
	$line = <S>;
	chomp $line;
	$quals = $line;
    } else {
	chomp $line;
	my @cols = split(/\:/,$line);
	$nts = $cols[5];
	$quals = $cols[6];
    }
    unless (rand($sampleRate) <= 1) {
	next;
    }
    $seq .= "X$nts";
    $qualBuffer .= "X$quals";
    $nSamples++;
    if ($nSamples >= $maxSamples) {
	last;
    }
}
close S;

my %prefixDistinct = ();
my %prefixCounts = ();
my %prefixPTotal = ();    

my $matchString = "[ACGT]{$merSize}";
my %distinctMers = ();

while ($seq =~ /($matchString)/g) {
    my $mer = $1;

    pos($seq) -= ($merSize-1);
    my $sCoord = pos($seq)-1;
    my $merQual = substr($qualBuffer,$sCoord,$merSize);

    if ($merQual =~ /$q2sym/) {
	next;
    }

    my $newMer = 1;
    if (exists($distinctMers{$mer})) {
	$newMer = 0;
    } else {
	$distinctMers{$mer} = 1;
    }

    my @q = split(/ */,$merQual);
    my @pCorrect = map {  1.0 - 10**( -(ord($_)-$qualOffset)/10.0) } @q;
    my $pTrue = 1.0;
    foreach my $p (@pCorrect) {
	$pTrue *= $p;
    }

    for (my $pLen = $minPrefixLen; $pLen <= $maxPrefixLen; $pLen++) {
	my ($prefix) = $mer =~ /^(.{$pLen})/;
	$prefixCounts{$prefix}++;
	if ($newMer) {
	    $prefixDistinct{$prefix}++;
	}
	$prefixPTotal{$prefix} += $pTrue;
    }
}

while (my ($p,$c) = each(%prefixCounts)) {
    my $pT = $prefixPTotal{$p};
    my $pMean = $pT/$c;
    my $eMean = 1-$pMean;
    my $nDistinct = $prefixDistinct{$p};

    print "$p\t$c\t$nDistinct\t$eMean\n";
}



    
    
