#!/usr/bin/perl -w
# details on where the gene neighbor score for given pair(s) comes from

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 $usage = "GeneNeighborScoreDetail.pl -genes Gene1,Gene2,...,GeneN  -gnc genomeSim -cluster=similarGenomes\n" .
	"   [-minScaffoldk=300] [-debug=10] [-maxSep=5000] [-host host] > scoresFile\n" .
	"GeneNeighborScoreDetail prints to STDOUT the log-scores for the gene neighbor method\n".
	"for each pair of genes from each genome with them both present (whether adjacent or not)\n".
	"The cluster is indicated so that you can figure out which score was actually used\n";
    
    my $locuslist = "";
    my $minScaffoldK = 300;
    my $clusterFile = "";
    my $debug = 0;
    my $maxSep = 5000;
    my $host=undef;
    my $gnc = "";
    (GetOptions('genes=s' => \$locuslist,
		'minScaffoldK=i' => \$minScaffoldK,
		'cluster=s' => \$clusterFile,
		'debug=i' => \$debug,
		'maxSep=i' => \$maxSep,
		'gnc=s' => \$gnc,
		'host=s' => \$host)
     && @ARGV == 0) || die $usage;
    die $usage if !$locuslist || !$clusterFile;

    my ($clusterToTaxonsRef, $taxonToClusterRef) = Cluster::readCluster($clusterFile);
    my $genomeSimRef = Vector::readDelim($gnc, "TaxId2");

    GenomicsUtils::setDebug($debug) if $debug > 0;
    connect($host);

    my @locuss = split /,/, $locuslist;
    die $usage if scalar @locuss < 2;
    my @genes = ();
    foreach my $locus (@locuss) {
	my $gene = Gene::new('locusId' => $locus);
	die "Bad locus $locus" unless $gene->locusId() eq $locus;
	push @genes, $gene;
    }
    my @taxons = map { $_->scaffold()->taxonomyId() } @genes;
    my $taxon = $taxons[0];
    die "Locuses not from the same genome" unless scalar Vector::uniq(sort(@taxons)) == 1;
    die "Locuses from taxon $taxon, similarity scores from taxon " . ( values %{ $genomeSimRef->{"TaxId1"} } )[0]
	unless $taxon eq ( values %{ $genomeSimRef->{"TaxId1"} } )[0];

    my %orths = (); # gene -> genome -> list of orths

    foreach my $gene ( @genes ) {
	my $scaffold = $gene->scaffold();
	my $locusId = $gene->locusId();
	$orths{$locusId} = {};
	map { addToHashList($orths{$locusId}, $_->scaffold()->taxonomyId(), $_); } $gene->getOrthologList();
    }

    print  join("\t", "Taxon","Gene1","Gene2","OrthCluster","OrthTaxon","Score","Separation","Code", "OTName"),"\n";
    foreach my $gene1 (@genes) {
	foreach my $gene2 (@genes) {
	    next if $gene1->locusId() >= $gene2->locusId();
	    my $orths1 = $orths{$gene1->locusId()};
	    my $orths2 = $orths{$gene2->locusId()};

	    my %consAdj = (); # genome -> separation code if near in that genome
	    my %score = (); # genome -> score if both present in that genome
	    my %sep = (); # genome -> sep if both present in that genome

	    foreach my $genomeOrth (keys %{$orths1}) {
		next if !defined $orths2->{$genomeOrth};
		foreach my $orth1 (@{ $orths1->{$genomeOrth} }) {
		    foreach my $orth2 (@{ $orths2->{$genomeOrth} }) {
			next if $orth1->scaffoldId() ne $orth2->scaffoldId();
			my ($sepOrth, $codeOrth) = GeneMethods::geneSeparation($orth1, $orth2);
			die $genomeOrth if !exists $genomeSimRef->{"nConvA"}{$genomeOrth};
			my $nConvA = $genomeSimRef->{"nConvA"}{$genomeOrth};
			$nConvA = 0.5 if $nConvA < 1;
			my $nConvTot = $genomeSimRef->{"nConvN"}{$genomeOrth} +
			    $genomeSimRef->{"nConvO"}{$genomeOrth} +
			    $genomeSimRef->{"nConvB"}{$genomeOrth};
			$score{$genomeOrth} = ($sepOrth < $maxSep ? log($nConvTot/$nConvA) : 0 );
			$sep{$genomeOrth} = $sepOrth if ($sepOrth < $maxSep);
			$consAdj{$genomeOrth} = $codeOrth if ($sepOrth < $maxSep);
		    }
		}
	    }
	    my $total = 0;
	    foreach my $cluster (sort {$a <=> $b} keys %{ $clusterToTaxonsRef}) {
		my $clusterBest = 0;
		foreach my $genomeOrth (sort {$a <=> $b} @{ $clusterToTaxonsRef->{$cluster} }) {
		    if (exists $score{$genomeOrth}) {
			my $sep = "";
			$sep = $sep{$genomeOrth} if exists $sep{$genomeOrth};
			my $code = "";
			$code = $consAdj{$genomeOrth} if exists $consAdj{$genomeOrth};
			print join("\t", $taxon, $gene1->locusId(), $gene2->locusId(),
				   $cluster, $genomeOrth, $score{$genomeOrth}, $sep, $code,
				   Genome::new(taxonomyId => $genomeOrth)->name() ), "\n";
			$clusterBest = $score{$genomeOrth} if $score{$genomeOrth} > $clusterBest;
		    }
		}
		$total += $clusterBest;
	    }
	    print STDERR join("\t","GNScore",$taxon,$gene1->locusId(),$gene2->locusId(),$total),"\n";
	}
    }
}
