#!/usr/bin/perl -w
#lists of conserved neighboring genes in 1 genome, and similarities to other genomes

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

{
    my $usage = "Usage: GeneNeighborsConserved.pl -taxon=TaxonId\n" .
	"                [-minScaffoldK=300 -maxSep=5000 -out=gncTaxon -debug=10 -host=localhost]\n" .
	"GeneNeighborsConserved prints to STDOUT the similarities of pairs of genomes\n" .
	"According to #pairs of adjacent genes with orthologous (and nearby) pairs\n" .
	"Output files are gncTaxon.locuss and gncTaxon.genomes, tab-delimited with columns\n";

    my $taxon = ""; # actually the starting taxonomy id
    my $minScaffoldK = 300;
    my $maxSep = 5000;
    my $out = "";
    my $debug = 0;
    my $host = undef;
    (GetOptions('taxon=s' => \$taxon,
		'minScaffoldK=i' => \$minScaffoldK,
		'maxSep=i' => \$maxSep,
		'out=s' => \$out,
		'host=s' => \$host,
		'debug=i' => \$debug)
     && @ARGV == 0) || die $usage;
    die $usage if !$taxon;
    $out = "gnc".$taxon if !$out;
    GenomicsUtils::setDebug($debug) if $debug > 0;
    connect($host);
    
    my $genome = Genome::new("taxonomyId" => $taxon);

    open LocusS, '>', "$out.locuss" || die "Cannot write to $out.locuss";
    print LocusS join("\t", qw/Genome1 Genome2 Gene Next PairConserved PairNear PairOrientation/),"\n";

    # per pair (indexed by the main gene's locus ID):
    my %nexts = (); # locusId -> gene, next
    my %sepCodes = ();
    my %separations = ();

    my %genomesOrth = (); # list of genomes with orthologs for any gene (to produce neither)

    # per pair, per other genome (1 and 2 is shorthand for main gene and next gene)
    # gene -> hash of genomes (all values are 1)
    my %orthStatus = (); # gene -> genome -> 1, 2, or 12 (neither not included)
    my %orth12sep = (); # gene -> genome -> (separation, codes) # only for the 12 (both) genomes
    my ($SEP,$CODE) = (0,1);

    my %scScore = ( "Convergent" => 0, "Divergent" => 0, "Same" => 1 );

    my @genes = $genome->genes(preload => ['ortholog']);
    foreach my $gene ( @genes ){
	# nextGene is per-scaffold and handles circular and partial scaffolds correctly
	next unless my $next = $gene->nextGene();
	my $scaffold = $gene->scaffold();
	next unless $scaffold->length() >= $minScaffoldK*1000;
	my $locusId = $gene->locusId();
	$nexts{$locusId} = [$gene,$next];
	($separations{$locusId}, $sepCodes{$locusId}) = GeneMethods::geneSeparation($gene, $next);

	$orthStatus{$locusId} = {};
	$orth12sep{$locusId} = {};
	    
	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();
	
	foreach my $genomeOrth (keys %orths1) {
	    $genomesOrth{$genomeOrth} = 1;
	    if (defined $orths2{$genomeOrth}) {
		$orthStatus{$locusId}{$genomeOrth} = 12;
		foreach my $orth1 (@{ $orths1{$genomeOrth} }) {
		    foreach my $orth2 (@{ $orths2{$genomeOrth} }) {
			next if $orth1->scaffoldId() ne $orth2->scaffoldId();
			my ($sep, $code) = GeneMethods::geneSeparation($orth1, $orth2);
			# Ignore orthologous sepCode because of issues with strand bias...
			    #|| $scScore{$code} > $scScore{ $orth12sep{$locusId}{$genomeOrth}[$CODE] }
			    #|| ($scScore{$orth12sep{$locusId}{$genomeOrth}[$CODE]} == $scScore{$code}
			if (!defined $orth12sep{$locusId}{$genomeOrth}
			    || $sep < $orth12sep{$locusId}{$genomeOrth}[$SEP]) {
			    $orth12sep{$locusId}{$genomeOrth} = [$sep, $code];
			}
		    }
		}
	    } else {
		$orthStatus{$locusId}{$genomeOrth} = 1;
	    }
	}
	foreach my $genomeOrth (keys %orths2) {
	    $genomesOrth{$genomeOrth} = 1;
	    next if $genomeOrth eq $genome || defined $orths1{$genomeOrth};
	    $orthStatus{$locusId}{$genomeOrth} = 2;
	}
	foreach my $genomeOrth (sort keys %{ $orthStatus{$locusId}}) {
	    print LocusS join("\t", $genome->taxonomyId(), $genomeOrth, $locusId, $next->locusId(),
			    ($orthStatus{$locusId}{$genomeOrth}==12 ? 1 : 0),
			    (exists $orth12sep{$locusId}{$genomeOrth} &&
			     $orth12sep{$locusId}{$genomeOrth}[$SEP] <= $maxSep ? 1 : 0),
			    (exists $orth12sep{$locusId}{$genomeOrth} ?
			     $orth12sep{$locusId}{$genomeOrth}[$CODE] : "None")
			    ),"\n";
	}
    }
    close(LocusS) || die "Cannot close $out.locuss";

    # simNeighbor records pairs with the given orientation in the source genome
    # which are conserved (any orientation) or adjacent (original orientation AND nearby)...
    # simProfile is the same summed over all orientations
    my %simNeighbor = (); # taxon2 -> sepCode -> nBoth, nOne, nNeither, nAdj
    my %simProfile = (); # taxon2 -> nBoth, nOne, nNeither, nAdj (ignoring orientation)
    my ($BOTH, $ONE, $NEITHER, $ADJ) = 0..3;

    foreach my $locusId (keys %nexts) {
	my ($gene, $next) = @{ $nexts{$locusId} };
	my $sepCode = $sepCodes{$locusId};
	foreach my $taxon2 (keys %genomesOrth) {
	    $simProfile{$taxon2} = [0,0,0,0]
		if !defined $simProfile{$taxon2};
	    if (!defined $simNeighbor{$taxon2}) {
		$simNeighbor{$taxon2} = { "Same" => [0,0,0,0],
					  "Convergent" => [0,0,0,0],
					  "Divergent" => [0,0,0,0] };
	    }
	}
	foreach my $taxon2 (keys %{ $orthStatus{$locusId}}) {
	    my $shared = ( $orthStatus{$locusId}{$taxon2} == 12 ? $BOTH : $ONE );
	    $simProfile{$taxon2}[$shared]++;
	    $simNeighbor{$taxon2}{$sepCode}[$shared]++;
	    if (defined $orth12sep{$locusId}{$taxon2}
		&& $orth12sep{$locusId}{$taxon2}[$SEP] <= $maxSep) {
		$simProfile{$taxon2}[$ADJ]++;
		$simNeighbor{$taxon2}{$sepCode}[$ADJ]++;
		    # ignore sepCode because it makes things tricky with strand bias...
		    # if $orth12sep{$locusId}{$taxon2}[$CODE] eq $sepCode;
	    }
	}
	foreach my $taxon2 (keys %genomesOrth) {
	    if (!defined $orthStatus{$locusId}{$taxon2}) {
		$simProfile{$taxon2}[$NEITHER]++;
		$simNeighbor{$taxon2}{$sepCode}[$NEITHER]++;
	    }
	}
    }

    open(GENOMES, '>', "$out.genomes") || die "Cannot write to $out.genomes";
    print GENOMES join("\t", qw/TaxId1 Genome1 TaxId2 Genome2 nTot1 nTot2 nBoth nOne nNeither nAdj
		       nSameB nSameO nSameN nSameA
		       nConvB nConvO nConvN nConvA
		       nDivB  nDivO  nDivN  nDivA/), "\n";
    foreach my $taxon2 (sort keys %simProfile) {
	if ($taxon2 ne $taxon) {
	    print GENOMES join("\t", $taxon, $genome->name(),
			       $taxon2, Genome::new("taxonomyId"=>$taxon2)->name(),
			       scalar (keys %nexts), # genes considered in this genome
			       # this should give all genes considered by the ortholog table for taxon2
			       &queryScalar("SELECT COUNT(*) FROM Locus o,Position p,Scaffold s"
					    . " WHERE o.posId=p.posId AND o.priority=1 AND p.scaffoldId=s.scaffoldId"
					    . " AND s.isGenomic=1 AND s.isActive=1 AND s.taxonomyId=$taxon2"),
			       @{ $simProfile{$taxon2} },
			       @{ $simNeighbor{$taxon2}{"Same"} },
			       @{ $simNeighbor{$taxon2}{"Convergent"} },
			       @{ $simNeighbor{$taxon2}{"Divergent"} } ), "\n";
	}
    }
    close(GENOMES) || die "Cannot close $out.genomes";
}

