#!/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 Information;
use Cluster;
use Vector;

{
    my $usage = "GenePhyloScores.pl -taxon=Id -gnc=genomeSim -cluster=similarGenomes\n" .
	"   [-Minscaffoldk=300] [-debug=10] [-host=host] > scoresFile\n" .
	"GenePhyloScores prints to STDOUT the log-scores for the phylo profile method\n".
	"for pairs of adjacent genes from the given genome\n";
    
    my $taxon = "";
    my $minScaffoldK = 300;
    my $clusterFile = "";
    my $debug = 0;
    my $gnc = "";
    my $host = undef;
    (GetOptions('taxon=s' => \$taxon,
		'minScaffoldK=i' => \$minScaffoldK,
		'cluster=s' => \$clusterFile,
		'debug=i' => \$debug,
		'host=s' => \$host,
		'gnc=s' => \$gnc)
     && @ARGV == 0) || die $usage;
    die $usage if !$taxon || !$clusterFile || !$gnc;
    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']);

    print join("\t","Taxon","Gene1","Gene2","Sep","Code","Strand1","Strand2","MI","rProfile","N1","N2"),"\n";

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

	my ($sep, $sepCode) =  GeneMethods::geneSeparation($gene, $next);
	my %orths1 = (); # genome -> list of orths
	my %orths2 = (); # genome -> list of orths
	map { addToHashList(\%orths1, $_->scaffold()->taxonomyId(), $_); } $gene->getOrthologList();
	map { addToHashList(\%orths2, $_->scaffold()->taxonomyId(), $_); } $next->getOrthologList();

	my %cons = (); # genome -> 1, 2, 12, or absent if neither
	foreach my $genomeOrth (keys %orths1) {
	    $cons{$genomeOrth} = (exists $orths2{$genomeOrth} ? 12 : 1);
	}
	foreach my $genomeOrth (keys %orths2) {
	    $cons{$genomeOrth} = 2 if !exists $cons{$genomeOrth};
	}

	# Now compute the mutual information and the correlation, taking clusters into account
	my ($nc1, $nc2, $nc12, $ncNeither, $nc12Miss) = (0,0,0,0,0); # cluster counts for MI
	my ($sumf12,$sumf1,$sumf2,$sum2f1,$sum2f2) = (0,0,0,0,0);

	foreach my $cluster (keys %{ $clusterToTaxonsRef }) {
	    next if $taxonToClusterRef->{$taxon} eq $cluster; # ignore self clusters!
	    my ($n1, $n2, $n12, $nNeither) = (0,0,0,0); # within this cluster

	    foreach my $genomeOrth (@{ $clusterToTaxonsRef->{$cluster} }) {
		if (exists $cons{$genomeOrth}) {
		    my $c = $cons{$genomeOrth};
		    $n1++ if $c == 1 || $c == 12;
		    $n2++ if $c == 2 || $c == 12;
		    $n12++ if $c == 12;
		} else {
		    $nNeither++;
		}
	    }
	    my $tot = scalar @{ $clusterToTaxonsRef->{$cluster} };
	    my ($f1, $f2, $f12) = ($n1/$tot, $n2/$tot, $n12/$tot);
	    $sumf12 += $f12;
	    $sumf1 += $f1;
	    $sumf2 += $f2;
	    $sum2f1 += $f1*$f1;
	    $sum2f2 += $f2*$f2;
	    $nc1++ if $n1 > 0;
	    $nc2++ if $n2 > 0;
	    $ncNeither++ if $n1==0 && $n2==0;
	    if ($n1 > 0 && $n2 > 0) {
		if ($n12 == 0) {
		    $nc12Miss++;
		} else {
		    $nc12++;
		}
	    }
	}
	my $nc1Only = $nc1 - $nc12 - $nc12Miss;
	my $nc2Only = $nc2 - $nc12 - $nc12Miss;
	my $ncTot = $nc1Only + $nc2Only + $nc12 + $nc12Miss + $ncNeither;
	my $info1 = Information::information($nc1,$ncTot-$nc1);
	my $info2 = Information::information($nc2,$ncTot-$nc2);
	my $info12 = Information::information($nc1Only,$nc2Only,$nc12,$nc12Miss,$ncNeither);
	my $MI = $info1 + $info2 - $info12;
	my $corr = 0;
	if ($sumf1 > 0 && $sumf2 > 0 && abs($sumf1-$ncTot) > 1e-5 && abs($sumf2-$ncTot) > 1e-5) {
	    my $meanf1 = $sumf1/$ncTot;
	    my $meanf2 = $sumf2/$ncTot;
	    $corr = (($sumf12/$ncTot) - $meanf1*$meanf2) / sqrt($meanf1*$meanf2*(1-$meanf1)*(1-$meanf2));
	}

	print join("\t",$taxon, $gene->locusId(), $next->locusId(), $sep, $sepCode,
		   $gene->strand(), $next->strand(),
		   $MI, $corr,$nc1,$nc2),"\n";
    }
}
