#!/usr/bin/perl -w
use strict;

# PCR error rate - pure guesswork.  At the very least set this to 1e-15 in
# order to prevent underflows in the floating point code as it's not
# particularly rigorous.
#
# The affect of reducing this is to increase Q1 scores and consequentially
# also give more variability to U2 (SAM Q2/Q3) scores.

my $pPCR = 2e-5;  # 1 in 50,000

#-----------------------------------------------------------------------------
# SAM line generation
sub output_sam {
    my ($end, $name,$rname,$rg,$pos,$strand,$s1,$q1,$s2,$q2) = @_;
    my $flag = 1; #paired
    my $cigar = length($s1) . "M";
    my $aux="RG:Z:$rg\tE2:Z:$s2\tU2:Z:$q2";

    $flag += 0x10 if (($strand eq "-") != ($end == 1));
    $flag += $end ? 0x80 : 0x40;

    print "$name\t$flag\t$rname\t$pos\t255\t$cigar\t=\t0\t0\t$s1\t$q1\t$aux\n";
}

#-----------------------------------------------------------------------------
# Arg handling
if (scalar(@ARGV) < 1) {
    print STDERR "Usage: like2sam library_file.gz [a.like ...]\n";
    exit 1
}
my $rg = shift(@ARGV); #Library file
open(my $fh, "gzip -cd < $rg |") || die "$rg";

my $line = 0;

print "\@SQ\tSN:NA19240_chr20_50000_10049999_hap1\tLN:10000000\n";

#-----------------------------------------------------------------------------
# Main processing loop

while (<>) {
    my @F = split("\t", $_);
    my $lib_name = <$fh>;
    my $lib_seq  = <$fh>;

    $line++;

    # Skip filtered data
    next if scalar(@F) <= 4;

    # Parse library name line
    $lib_name =~ />(\S+).*Strand (.) .* (\d+)--(\d+)/;
    my ($rname,$strand,$start,$end) = ($1,$2,$3,$4);

    if ($strand eq "-") {
	my ($rst,$ren) = ($rname =~ /_(\d+)_(\d+)_/);
	my $len = $end - $start;

	$start = ($ren-$rst) - $end + 2;
	$end   = $start + $len;
    } else {
	$start++;
    }

    # This is the correct name from illumina format data, but we seem to have
    # collisions. I expect the names have been generated randomly in
    # the simulated test set.
    #my $name = "$F[0]_$F[1]_$F[2]_$F[3]";
    my $name = "s_$line";

    my $seq1  = ""; my $seq2  = "";
    my $qual1 = ""; my $qual2 = "";
    for (my $i = 4; $i <= $#F; $i++) {
	my $s1 = "";   my $s2 = "";
	my $q1 = -100; my $q2 = -100; my $q3 = -100;
	my @q = split(" ", $F[$i]);
	my $p_tot = 0;

	for (my $j = 0; $j < 4; $j++) {
	    my $prob = exp(-$q[$j]); 
	    $p_tot += $prob;

	    if ($prob > $q1) {
		$q3 = $q2;
		$q2 = $q1;   $s2 = $s1;
		$q1 = $prob; $s1 = $j;
	    } elsif ($prob > $q2) {
		$q3 = $q2;
		$q2 = $prob; $s2 = $j;
	    } elsif ($prob > $q3) {
		$q3 = $prob;
	    }
	}

	# Renormalise and add in a chance of an external event (eg PCR)
	# causing a mutation prior to sequencing.
	$q1 = (1-$pPCR)*$q1/$p_tot + $pPCR/3;
	$q2 = (1-$pPCR)*$q2/$p_tot + $pPCR/3;
	$q3 = (1-$pPCR)*$q3/$p_tot + $pPCR/3;

	# Q2 defined as a conditional probability
	#$q2 = int(-10*log(1 - $q2/(1-$q1))/log(10) + .5);

	# Q2 defines as per-SAM, the ratio of 2nd to 3rd likelihoods.
	# Given we cannot work out Q3, we estimate it as the remainder.
	# This works out as Q2=10*($p[1]-$p[0])/ln(10)) if @p was sorted in
	# order of best[0] to poorest[3] calls. However we can't add in the
	# fixed PCR error via this route.
	#
	$q2 = int(10*log($q2/$q3)/log(10)+.5);

	$q1 = int(-10*log(1-$q1)/log(10) + .5);

	$q1 = 126-33 if ($q1 > 126-33);
	$q2 = 126-33 if ($q2 > 126-33);
#	$q1 = 255-33 if ($q1 > 255-33);
#	$q2 = 255-33 if ($q2 > 255-33);
	$qual1 .= chr(33+$q1);
	$qual2 .= chr(33+$q2);

	$seq1 .= $s1;
	$seq2 .= $s2;
    }
    
    $seq1 =~ tr/0123/ACGT/;
    $seq2 =~ tr/0123/ACGT/;

    # Paired end, so split.
    my $len = length($seq1) / 2;

    if ($strand eq "-") {
	output_sam(1, $name, $rname, $rg, $start, $strand,
		   substr($seq1, $len, $len), substr($qual1, $len, $len),
		   substr($seq2, $len, $len), substr($qual2, $len, $len));

	$seq1  = reverse(substr($seq1,  0, $len));
	$seq1  =~ tr/acgtACGT/tgcaTGCA/;
	$qual1 = reverse(substr($qual1, 0, $len));
	$seq2  = reverse(substr($seq2,  0, $len));
	$seq2  =~ tr/acgtACGT/tgcaTGCA/;
	$qual2 = reverse(substr($qual2, 0, $len));
	output_sam(0, $name, $rname, $rg, $end-$len, $strand,
		   $seq1, $qual1, $seq2, $qual2);
    } else {
	output_sam(0, $name, $rname, $rg, $start, $strand,
		   substr($seq1, 0,    $len), substr($qual1, 0,    $len),
		   substr($seq2, 0,    $len), substr($qual2, 0,    $len));

	$seq1  = reverse(substr($seq1,  $len, $len));
	$seq1  =~ tr/acgtACGT/tgcaTGCA/;
	$qual1 = reverse(substr($qual1, $len, $len));
	$seq2  = reverse(substr($seq2,  $len, $len));
	$seq2  =~ tr/acgtACGT/tgcaTGCA/;
	$qual2 = reverse(substr($qual2, $len, $len));
	output_sam(1, $name, $rname, $rg, $end-$len+1, $strand,
		   $seq1, $qual1, $seq2, $qual2);
    }
}

close($fh);

