#!/usr/bin/perl -w
#gene neighbor scores for list of adjacent genes across a genome

use strict;
use lib "../Perl/modules";
use Gene;
use GeneMethods;
use Scaffold;
use Genome;
use GenomicsUtils;
use Getopt::Long;
use lib "../stat";
use Vector;
use Cluster;

my $control = "nsame"; # the default

# empirical likelihood value for one genome
# returns a positive value for the suprise value of presence if $on,
# and a negative value for the suprise value of absence if !$on
sub llscore($$$) {
    my ($genomeSimRef, $genomeOrth, $on) = @_;
    die $genomeOrth if !exists $genomeSimRef->{"nTot2"}{$genomeOrth};
    my ($bottom, $top);
    if ($control eq "both") {
	$bottom = $genomeSimRef->{"nBoth"}{$genomeOrth};
	$bottom = .5 if $bottom == 0;
	$top = $genomeSimRef->{"nAdj"}{$genomeOrth};
    } elsif ($control eq "min") {
	my $nTot1 = $genomeSimRef->{"nTot1"}{$genomeOrth};
	my $nTot2 = $genomeSimRef->{"nTot2"}{$genomeOrth};
	$bottom = ($nTot1 < $nTot2 ? $nTot1 : $nTot2);
	$top = $genomeSimRef->{"nAdj"}{$genomeOrth};
    } elsif ($control eq "nsame") {
	$bottom = $genomeSimRef->{"nConvB"}{$genomeOrth} + $genomeSimRef->{"nDivB"}{$genomeOrth};
	$top = $genomeSimRef->{"nConvA"}{$genomeOrth} +$genomeSimRef->{"nDivA"}{$genomeOrth};
    }
    $top = 0.5 if $top < 1;
    $bottom = 1 if $bottom < 1;
    $top = $bottom - 0.5 if $bottom eq $top;
    my $value = log( $on ? $bottom/$top : 1-($top/$bottom) );
#    print STDERR join("\t","on",$on,"bottom",$bottom,"top",$top,
#		      "nConvB",$genomeSimRef->{"nConvB"}{$genomeOrth},
#		      "nDivB",$genomeSimRef->{"nDivB"}{$genomeOrth},
#		      "nConvA",$genomeSimRef->{"nConvA"}{$genomeOrth},
#		      "nDivA",$genomeSimRef->{"nDivA"}{$genomeOrth},
#		      "result",$value),"\n";
    return($value);
}

{
    my $usage = "GeneNeighborScores.pl -taxon=Id -gnc=genomeSim -cluster=similarGenomes\n"
	. "   [-control=nsame|both|min] [-minScaffoldk=300] [-geneWindow=1] [-debug=10] [-maxSep=5000]\n"
	. "   [-host host] > scoresFile\n"
	. "GeneNeighborScores prints to STDOUT the log-scores for the gene neighbor method\n"
	. "for pairs of adjacent genes from the given genome\n";
    
    my $taxon = "";
    my $minScaffoldK = 300;
    my $clusterFile = "";
    my $debug = 0;
    my $maxSep = 5000;
    my $gnc = "";
    my $geneWindow = 1;
    my $host=undef;
    (GetOptions('taxon=s' => \$taxon,
		'minScaffoldK=i' => \$minScaffoldK,
		'cluster=s' => \$clusterFile,
		'debug=i' => \$debug,
		'maxSep=i' => \$maxSep,
		'host=s' => \$host,
		'gnc=s' => \$gnc,
		'control=s' => \$control,
		'geneWindow=i' => \$geneWindow)
     && @ARGV == 0) || die $usage;
    die $usage if !$taxon || !$clusterFile || !$gnc;
    die $usage if $control ne "nsame" && $control ne "both" && $control ne "min";
    GenomicsUtils::setDebug($debug) if $debug > 0;
    connect($host);

    my ($clusterToTaxonsRef, $taxonToClusterRef) = Cluster::readCluster($clusterFile);
    print STDERR join("\t","Read Clusters:", scalar keys %{ $clusterToTaxonsRef },
		      "Genomes:", scalar keys %{ $taxonToClusterRef} ),"\n";
    my $genomeSimRef = Vector::readDelim($gnc, "TaxId2");
    print STDERR join("\t", "Read Genome Similarites:",
		      scalar keys %{ $genomeSimRef->{"TaxId2"} }),"\n";

    my $genome = Genome::new("taxonomyId" => $taxon);
    my @genes = $genome->genes(preload => ['ortholog']);

    # GNWithin is the within-cluster score (within negatives for genomes with the genes present but non-adjacent)
    # GNAll is the score if clustering is ignored
    print join("\t","Taxon","Scaffold","Gene1","Gene2","Consecutive",
	       "Sep","Code","Strand1","Strand2","GNScore","GNWithin","GNAll", "GNMinus"),"\n";

    my %orths = (); # gene -> genome -> list of orths
    
    
    foreach my $gene ( @genes ) {
	my $scaffold = $gene->scaffold();
	next unless $scaffold->length() >= $minScaffoldK*1000;
	my $locusId = $gene->locusId();
	$orths{$locusId} = {};
	map { addToHashList($orths{$locusId}, $_->scaffold()->taxonomyId(), $_); } $gene->getOrthologList();
    }

    foreach my $gene ( @genes ) {
	my $scaffold = $gene->scaffold();
	next unless $scaffold->length() >= $minScaffoldK*1000;
	my $next = $gene;
	my $orths1 = $orths{$gene->locusId()};

	foreach my $count (1..$geneWindow) {
	    last unless $next = $next->nextGene();
	    my ($sep, $sepCode) =  GeneMethods::geneSeparation($gene, $next);
	    my $orths2 = $orths{$next->locusId()};

	    my %consAdj = (); # genome -> 1 if near in that genome
	    my %consBoth = (); # genome -> 1 if both in that genome
	    foreach my $genomeOrth (keys %{$orths1}) {
		if (defined $orths2->{$genomeOrth}) {
		    foreach my $orth1 (@{ $orths1->{$genomeOrth} }) {
			foreach my $orth2 (@{ $orths2->{$genomeOrth} }) {
			    $consBoth{$genomeOrth}=1;
			    if ($orth1->scaffoldId() eq $orth2->scaffoldId()) {
				my ($sepOrth, $codeOrth) = GeneMethods::geneSeparation($orth1, $orth2);
				# ignore codeOrth because it makes things tricky to do controls
				# correctly if there is strand bias
				$consAdj{$genomeOrth}=1 if ($sepOrth < $maxSep);
			    }
			}
		    }
		}
	    }

	    # Now compute the logl value by summing over the best value in each cluster
	    my $llClustered = 0;
	    my $llAll = 0; # does not take clustering into account
	    my $llWithin = 0;
	    my $llMinus = 0; # all with penalties for suprising misses

	    foreach my $cluster (keys %{ $clusterToTaxonsRef }) {
		my $bestll = 0;
		foreach my $genomeOrth (@{ $clusterToTaxonsRef->{$cluster} }) {
		    next if $genomeOrth eq $taxon || !exists $genomeSimRef->{"nTot2"}{$genomeOrth};
		    if (exists $consAdj{$genomeOrth}) {
			my $ll = &llscore($genomeSimRef,$genomeOrth,1);
			$llAll += $ll;
			$bestll = $ll if $ll > $bestll;
		    }
		    if (exists $consBoth{$genomeOrth}) {
			my $ll2 = &llscore($genomeSimRef,$genomeOrth,exists $consAdj{$genomeOrth});
			$llMinus += $ll2;
			$llWithin += $ll2 if $taxonToClusterRef->{$taxon} eq $cluster;
		    }
		}
		# ignore self clusters
		$llClustered += $bestll if $taxonToClusterRef->{$taxon} ne $cluster;

	    }
	    print join("\t",$taxon, $scaffold->scaffoldId(),$gene->locusId(), $next->locusId(),
		       ($count==1?1:0),
		       $sep,$sepCode,
		       $gene->strand(), $next->strand(),$llClustered, $llWithin, $llAll, $llMinus),"\n";
	}
    }
}
