#!/usr/bin/perl -w
# Expects as input a tab-delimited primers table, with columns named "Locus", "Forward Primer",
# "Reverse Primer", and "Expected Length",
# and a blast output file where the query names are rOrf or fOrf (for the two primers)
# and the DB sequence names are VIMSS scaffoldIds, and an Locus position table
#
# EXAMPLE OF CREATING THE INPUT FILES
# (notice this example keeps non-genomic scaffolds)
#mysql -u test -ptest genomics_test -B -e 'select ss.scaffoldId, ss.sequence from from ScaffoldSeq ss,Scaffold s WHERE s.taxonomyId=211586 AND ss.sequence IS NOT NULL AND s.isActive=1;' | perl -F\\t -ane 'print ">$F[0]\n$F[1]" if $F[0] ne "ss.scaffoldId"' > ~/ShewanellaArrays/ScaffoldDB
#/usr2/people/gtl/bin/blast/formatdb -i ~/ShewanellaArrays/ScaffoldDB -o T -p F -l ~/ShewanellaArrays/ScaffoldDB.log
#/usr2/people/gtl/bin/blast/blastall -p blastn -i ~/ShewanellaArrays/BothPrimers.fasta -o ~/ShewanellaArrays/BothPrimers.blastScaffolds -m 8 -F F -d ~/ShewanellaArrays/ScaffoldDB
#### ignore blastall warnings on bad queries for the "OPEN" primers which are not in fact DNA sequences
#mysql -u test -ptest genomics_test -B -e 'select o.locusId locusId, p.scaffoldId scaffoldId, p.start start, p.stop stop, p.strand strand FROM Locus o, Position p, Scaffold s WHERE o.posId=p.posId AND p.scaffoldId=s.scaffoldId AND s.isActive=1 AND s.taxonomyId=211586 AND o.priority=1' > ~/ShewanellaArrays/Locuspos
#
# BUGS
# does not handle the special case of primers that cross the origin

use strict;
use lib "../stat";
use Vector;
use lib  "../Perl/modules";
use GenomicsUtils;
use Getopt::Long;
use POSIX; # for floor and ceil
{
    my $primers = "";
    my $blast = "";
    my $posFile = "";
    my $maxOff = 4;
    my $maxPCRk = 25;
    my $lengthErr = 25; # warn if length error is >= this value
    my $minOverlap = 100; # warn if overlap is under this value

    (GetOptions('primers=s' => \$primers,
		'blast=s' => \$blast,
		'pos=s' => \$posFile,
		'maxOff=i' => \$maxOff,
		'maxPCRk=i' => \$maxPCRk,
		'minOverlap=i' => \$minOverlap,
		'lengthTolerance=i' => \$lengthErr) && @ARGV == 0 && $primers && $blast && $posFile)
	|| die "Usage: ParsePrimers.pl -primers Primers.CVS -blast Primers.blast -pos Locuspos\n"
	       . " [-maxOff 4 -maxPCRk 25 -lengthTolerance 25 -minOverlap 100]\n";
    
    my $ptable = readDelim($primers,"Locus");
    foreach my $col ("Locus", "Expected Length", "Forward Primer", "Reverse Primer") {
	die "Primers must have a $col field" if !exists $ptable->{$col};
    }
    my $pos = readDelim($posFile, "locusId");
    foreach my $col ("locusId", "scaffoldId", "start", "stop", "strand") { # assume strand = +-, start < stop
	die "Positions must have a $col field" if !exists $pos->{$col};
    }
    my $posDiv = 10*1000;
    my %posToORFs = ();
    foreach my $locus (keys %{ $pos->{"locusId"} }) {
	my $index1 = floor($pos->{"start"}{$locus}/$posDiv);
	my $index2 = floor($pos->{"stop"}{$locus}/$posDiv);
	&addToHashList(\%posToORFs, $index1, $locus);
	&addToHashList(\%posToORFs, $index2, $locus) if $index1 != $index2;
    }

    my %locusToVIMSS = ();
    my %primerHits = (); # r or f . locusId -> list of scaffold, begin, end (implying strand), nOff
    my ($SCAFFOLD,$BEGIN,$END,$NOFF) = (0,1,2,3);
    
    open (BLAST, "<", $blast) || die "Cannot read blast file $blast";
    while(<BLAST>) {
	s/[\r\n]+$//;
	my ($query, $subject, $identity, $alignment, $mismatch, $gap, $qBegin, $qEnd, $sBegin, $sEnd, $evalue, $score) = split /\t/;
	my $forward;
	my $locus;
	if ($query =~ m/^f(.*)$/) {
	    $forward = 1; $locus = $1;
	} elsif ($query =~ m/^r(.*)$/) {
	    $forward = 0; $locus = $1;
	} else {
	    die "Cannot parse query named $query";
	}
	my $querySeq = $ptable->{$forward ? "Forward Primer" : "Reverse Primer"}{$locus};
	die "Unknown locus $locus in query named $query" if !defined($querySeq) || $querySeq eq "";
	my $nOff = $mismatch + length($querySeq)-($qEnd-$qBegin+1);
	addToHashList(\%primerHits, $query, [ $subject, $sBegin, $sEnd, $nOff ]) if ($nOff <= $maxOff);
    }
    close(BLAST) || die "Error closing blast file $blast";

    print join("\t", "Locus","VIMSS","Pos","LenDiff","overlap","nPos","nVIMSS","nForward","nReverse","nMismatch"),"\n";

    foreach my $locus (sort keys %{ $ptable->{"Locus"} }) {
	my @fHits = exists $primerHits{"f".$locus} ? @{ $primerHits{"f" . $locus} } : ();
	my @rHits = exists $primerHits{"r".$locus} ? @{ $primerHits{"r" . $locus} } : ();
	# we build these up as we go - note VIMSS is empty if no hits
	my @pos = ();
	my @VIMSS = (); # not a 1:1 mapping to positions
	my @lenDiff = ();
	my @lengthWarn = (); # do not print length discrepancies if they are from multiple hits
	my @strandWarn = ();
	my @overlap = ();
	my $mismatch = 0;
	foreach my $fHit (@fHits) {
	    foreach my $rHit (@rHits) {
		# the pointing together case, from the plus strand perspective, is begin1 end1 ... end2 begin2,
		# where 1=forward, 2=reverse, or vice versa
		if (abs($fHit->[$END] - $rHit->[$END]) < 1000*$maxPCRk
		    && $fHit->[$SCAFFOLD] eq $rHit->[$SCAFFOLD]
		    && ( $fHit->[$END] > $fHit->[$BEGIN] xor $rHit->[$END] > $rHit->[$BEGIN] ) # different strands
		    && ( $fHit->[$END] > $fHit->[$BEGIN] xor $fHit->[$END] > $rHit->[$END])) {
		    if ($fHit->[$NOFF] == 0 && $rHit->[$NOFF] == 0) {
			push @pos, join(":", $fHit->[$SCAFFOLD], $fHit->[$BEGIN], $rHit->[$BEGIN]);
			my $vimssThis = "none";
			my $minRange = Vector::min($fHit->[$BEGIN], $fHit->[$END], $rHit->[$BEGIN], $rHit->[$END]);
			my $maxRange = Vector::max($fHit->[$BEGIN], $fHit->[$END], $rHit->[$BEGIN], $rHit->[$END]);
			my $firstIndex = floor($minRange/$posDiv);
			my $lastIndex = ceil($maxRange/$posDiv);
			my %locusSeen = ();
			foreach my $index ($firstIndex..$lastIndex) {
			    if (exists $posToORFs{$index}) {
				foreach my $vimssId (@{ $posToORFs{$index} }) {
				    next if $locusSeen{$vimssId};
				    $locusSeen{$vimssId} = 1;
				    my $vimssStart = $pos->{"start"}{$vimssId};
				    my $vimssEnd = $pos->{"stop"}{$vimssId};
				    if ($vimssEnd >= $minRange && $vimssStart <= $maxRange
					&& $pos->{"scaffoldId"}{$vimssId} eq $fHit->[$SCAFFOLD]) {
					$vimssThis = $vimssId;
					push @VIMSS, $vimssId;
					if ($pos->{"strand"}{$vimssId} eq "+" xor $fHit->[$BEGIN] < $fHit->[$END]) {
					    push @strandWarn, "Warning: locus $locus VIMSS and $vimssId are on"
						. " opposite strands: "
						. join(" ","Locus hit", 
						       "s".$fHit->[$SCAFFOLD], $fHit->[$BEGIN], $fHit->[$END],
						       $rHit->[$BEGIN], $rHit->[$END], "VIMSS",
						       $vimssStart,$vimssEnd,$pos->{"strand"}{$vimssId}) . "\n";
					}
					push @overlap, min( abs($minRange-max($vimssStart,$vimssEnd)),
							    abs($maxRange-min($vimssStart,$vimssEnd)) );
				    }
				}
			    }
			} # end foreach index
			my $len = abs($fHit->[$BEGIN] - $rHit->[$BEGIN]);
			push @lenDiff,  $len - $ptable->{"Expected Length"}{$locus};
			push @lengthWarn, "Warning: locus $locus VIMSS $vimssThis"
			    . " hitting " . join(":","s".$fHit->[$SCAFFOLD],$fHit->[$BEGIN],$fHit->[$END],
						 $rHit->[$BEGIN],$rHit->[$END])
			    .  " length discrepancy: $len vs expected "
			    . $ptable->{"Expected Length"}{$locus} . "\n"
			    if (abs($len - $ptable->{"Expected Length"}{$locus}) > $lengthErr);
		    } else { # else NOFF != 0
			$mismatch++;
		    }
		} # end if hits point towards each other
	    } # end foreach rHit
	} # end foreach fHit
	my $nVIMSS = scalar @VIMSS;
	print join("\t", $locus, join(",", @VIMSS), join(",", @pos), join(",", @lenDiff), join(",", @overlap),
		   scalar @pos, $nVIMSS,
		   scalar @fHits, scalar @rHits, $mismatch), "\n";
	print STDERR "Warning: locus $locus matches $nVIMSS VIMSS ids; primers match " . scalar (@pos) . " locations"
	    . " and $mismatch additional non-perfect hits\n"
	    if scalar @pos != 1 || scalar @VIMSS != 1;
	if (scalar @pos == 1) {
	    map { print STDERR $_ } @strandWarn;
	    map { print STDERR $_ } @lengthWarn;
	    print STDERR "locus $locus matches VIMSS $VIMSS[0] uniquely but with overlap only " . $overlap[0] . "\n"
		if $nVIMSS == 1 && $overlap[0] < $minOverlap;
	}
    }
}
