#!/usr/bin/perl -w
#converts gnc*.genomes files into a list of similar genomes
#Uses simple cutoffs for %still converged and adjacent, or %conserved in both genomes, on
#adjacent converged pairs of genes, to decide if two genomes are similar.
#
#maxAdj and maxBoth divide the #adjacent or #both pairs by the minimum of the two genome sizes,
#and use all adjacent pairs, while ConvAdj and ConvBoth use convergent pairs only,
#and NSameAdj and NSameBoth use both convergent and divergent pairs (and divide by total).
#
#If for some reason there is insufficient information then one genome should probably be
#thrown out but this does NOT do that.

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

{
    my $usage = "Usage: ClusterGenomes.pl -in=gncAllGenomes\n"
	. "   -maxNSame=10 | -maxAdj=30 | -maxConvAdj=5 | -maxBoth=50 | -maxConvBoth=50  | -maxNSameBoth=50 > genomeClusters\n";

    my ($maxConvAdj, $maxConvBoth, $maxAdj, $maxBoth, $maxNSameAdj, $maxNSameBoth, $in);
    (GetOptions('maxAdj=i' => \$maxAdj,
	       'maxBoth=i' => \$maxBoth,
               'maxConvAdj=i' => \$maxConvAdj,
	       'maxConvBoth=i' => \$maxConvBoth,
               'maxNSameAdj=i' => \$maxNSameAdj,
	       'maxNSameBoth=i' => \$maxNSameBoth,
	       'in=s' => \$in) && @ARGV == 0 && defined($in)) || die $usage;
    $maxConvAdj || $maxConvBoth || $maxBoth || $maxAdj || $maxNSameAdj || $maxNSameBoth || die $usage;

    my %scores = (); # genome1 -> genome2 -> list of values indexed as above
    my %objects = (); # make sure we get objects both ways
    my %names = ();

    my $genomes = Vector::readDelimList($in);
    foreach my $name (qw(TaxId1 TaxId2 Genome1 Genome2 nTot1 nTot2 nBoth nAdj)) {
	die "Input file missing column $name" if !exists $genomes->[0]{$name};
    }

    foreach my $row (@$genomes) {
	$scores{$row->{"TaxId1"}} = {} if !exists $scores{$row->{"TaxId1"}};
	die "Duplicate value for Genome1 -> Genome2" if exists $scores{$row->{"TaxId1"}}{$row->{"TaxId2"}};
	$scores{$row->{"TaxId1"}}{$row->{"TaxId2"}} = $row;
	$objects{$row->{"TaxId1"}} = 1;
	$objects{$row->{"TaxId2"}} = 1;
	$names{$row->{"TaxId1"}} = $row->{"Genome1"};
	$names{$row->{"TaxId2"}} = $row->{"Genome2"};
    }
    print STDERR "Read " . Vector::sum( map {scalar keys %{$scores{$_}}} (keys %scores) ) . " scores\n";

    my %similar = map { $_ => {} } keys %scores;
    foreach my $tax1 (keys %scores) {
	foreach my $tax2 (keys %{ $scores{$tax1} }) {
	    my $v = $scores{$tax1}{$tax2};

	    my $convTot = $v->{"nConvB"}+$v->{"nConvO"}+$v->{"nConvN"};
	    my $divTot = $v->{"nDivB"}+$v->{"nDivO"}+$v->{"nDivN"};
	    my $nsameTot = $convTot + $divTot;
	    my $min = min($v->{"nTot1"}, $v->{"nTot2"}); # new way: normalize by genome size
	    if ($maxAdj) {
		$similar{$tax1}{$tax2} = 1 if $v->{"nAdj"} > ($maxAdj/100.0)*$min;
	    } elsif ($maxBoth) {
		$similar{$tax1}{$tax2} = 1 if $v->{"nBoth"} > ($maxBoth/100.0)*$min;
	    } elsif ($maxConvAdj) {
		$similar{$tax1}{$tax2} = 1 if $v->{"nConvA"} >= ($maxConvAdj/100.0)*$convTot;
	    } elsif ($maxConvBoth) {
		$similar{$tax1}{$tax2} = 1 if $v->{"nConvB"} >= ($maxConvBoth/100.0)*$convTot;
	    } elsif ($maxNSameAdj) {
		$similar{$tax1}{$tax2} = 1 if $v->{"nConvA"} + $v->{"nSameA"} - 1 >= ($maxNSameAdj/100.0)*$nsameTot;
	    } elsif ($maxNSameBoth) {
		$similar{$tax1}{$tax2} = 1 if $v->{"nConvB"} + $v->{"nSameB"} - 1 >= ($maxNSameBoth/100.0)*$nsameTot;
	    }
	}
    }

    print STDERR "Similar pairs: " . Vector::sum( map {scalar keys %{$similar{$_}}} (keys %similar) ) . "\n";
    
    my @objects = keys %objects;
    my ($clustersRef, $objsRef) = Cluster::clusterSingly(\@objects, \%similar);

    print STDERR "Clustered " . scalar (keys %{ $objsRef }) . " genomes into "
	. (scalar keys %{ $clustersRef }) . " clusters\n";

    print join("\t","Cluster","Size","TaxId","Genome"),"\n";
    foreach my $cluster (sort keys %{ $clustersRef}) {
	foreach my $tax1 (@{ $clustersRef->{$cluster} }) {
	    print join("\t",$cluster,scalar @{ $clustersRef->{$cluster} }, $tax1,$names{$tax1}),"\n";
	}
    }
}

