#!/usr/bin/perl -w
# Load predicted operons into the database.
# Actually produces a file of SQL commands to execute (including references to the files
# produced, to upload them). Those SQL commands must then be run on the database server.

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

{
    my $operonObject = 2; # code in the ObjectDescription table
    my $usage = "LoadOperons.pl -dir=~/operons/out [-use bOp] [-taxon=Name] [-out=dir] [-prefix=gnc]\n"
	. "   [-host=host] [-debug=10] > sqlCommands\n";

    my $taxonSpec = "";
    my $dir = "";
    my $out = "";
    my $prefix = "gnc";
    my $host=undef;
    my $debug=0;
    my $use="bOp";
    (GetOptions('taxon=s' => \$taxonSpec,
		'dir=s' => \$dir,
		'out=s' => \$out,
		'prefix=s' => \$prefix,
		'debug=i' => \$debug,
		'use=s' => \$use,
		'host=s' => \$host)
     && @ARGV == 0) || die $usage;
    die $usage if !$dir;
    die "Not a directory: $dir" if !(-d $dir);
    $out = $dir if !$out;
    die "Not a directory: $out" if !(-d $out);

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

    my @genomes = ();
    if ($taxonSpec) {
	@genomes = ( Genome::new("taxonomyId" => $taxonSpec) );
    } else {
	@genomes = Genome::getGenomeList();
    }
    
    my $Locus2OperonFile = "$out/" . $prefix . ($taxonSpec ? $taxonSpec : "AllGenomes") . ".Locus2Operon";
    my $OperonFile = "$out/" . $prefix . ($taxonSpec ? $taxonSpec : "AllGenomes") . ".Operon";
    my $OperonPosFile = "$out/" . $prefix . ($taxonSpec ? $taxonSpec : "AllGenomes") . ".OperonPosition";

    print "# Files to be read by SQL being created in:\n# $Locus2OperonFile\n# $OperonFile\n# $OperonPosFile\n";
    open(O2O, '>', $Locus2OperonFile) || die "Cannot open $Locus2OperonFile";
    open(OPERON, '>', $OperonFile) || die "Cannot open $OperonFile";
    open(OPOS, '>', $OperonPosFile) || die "Cannot open $OperonPosFile";

    my $tuId = GenomicsUtils::queryScalar("SELECT MAX(tuId) FROM Operon");
    my $posId = GenomicsUtils::queryScalar("SELECT MAX(posId) FROM Position");
    $tuId = 0 unless $tuId > 0;
    $posId = 0 unless $posId > 0;
    print "# Maximum known Operon.tuId: $tuId\n";
    print "# Maximum known Position.posId: $posId\n";

    foreach my $genome (@genomes) {
	my @genes = $genome->genes();
	my %genes = map { $_->locusId() => $_ } @genes;

	my $taxon = $genome->taxonomyId();
	my $caiFile = "$dir/$prefix$taxon.cai";
	my $predsFile = "$dir/$prefix$taxon.pred";

	if( !(-e $caiFile) || !(-e $predsFile)) {
	    print STDERR "# Warning: no operons predictions for $taxon\n";
	    next;
	}
	print STDERR "# Setting up operons tables for $taxon\n";

	my $caisRef = Vector::readDelimList($caiFile);
	my $predsRef = Vector::readDelimList($predsFile);
	die "Invalid use $use" if !exists $predsRef->[0]->{$use};

	my %sameOp = (); # stores locus1->locus2->1 if they are in the same operon
	my %locuss = (); # list of locuss

	foreach my $caisLine (@{ $caisRef }) {
	    my $locusId1 = $caisLine->{"Gene1"};
	    $locusId1 =~ s/ //g;
	    my $locusId2 = $caisLine->{"Gene2"};
	    $locusId2 =~ s/ //g;
	    $locuss{ $locusId1 } = 1;
	    $locuss{ $locusId2 } = 1;
	}
	foreach my $predsLine (@{ $predsRef }) {
	    my $locusId1 = $predsLine->{"Gene1"};
	    my $locusId2 = $predsLine->{"Gene2"};
	    $locusId1 =~ s/ //g;
	    $locusId2 =~ s/ //g;
	    die $locusId1 if !exists $locuss{$locusId1};
	    die $locusId2 if !exists $locuss{$locusId2};

	    my $val = $predsLine->{$use};

	    # allow true-false or p from 0 to 1 threshold .5
	    $val = 1 if ($val =~ m/TRUE/);
	    $val = 0 if ($val =~ m/FALSE/);

	    if ($val >= 0.5) {
		$sameOp{$locusId1} = {} if !exists $sameOp{$locusId1};
		$sameOp{$locusId1}{$locusId2} = 1;
	    }
	}
	my @locuss = sort {$a <=> $b} keys %locuss;
	my ($clustersRef, $objsRef) = Cluster::clusterSingly(\@locuss, \%sameOp);

	foreach my $cluster (sort {$a <=> $b} keys %{ $clustersRef }) {
	    $tuId++;
	    $posId++;
	    my @opOrfs = @{ $clustersRef->{$cluster} };
	    my $gene1 = $genes{$opOrfs[0]};
	    my $strand = $gene1->strand();
	    my $scaffoldId = $gene1->scaffoldId();
	    my $begin = Vector::min(map { $genes{$_}->begin() } @opOrfs );
	    my $end = Vector::max(map { $genes{$_}->end() } @opOrfs );
	    print OPOS join("\t",$posId,$scaffoldId,$operonObject,$strand,$begin,$end),"\n";
	    foreach my $locus (@opOrfs) {
		print O2O "$locus\t$tuId\n";
	    }
	    print OPERON "$tuId\t$posId\tVIMSS Predicted\n";
	}
    }

    close(O2O) || die "Cannot close $Locus2OperonFile";
    close(OPERON) || die "Cannot close $OperonFile";
    close(OPOS) || die "Cannot close $OperonPosFile";

    print "use genomics_test;\n";
    foreach my $genome (@genomes) {
	my $taxon = $genome->taxonomyId();
	my $where = "Locus2Operon.tuId=Operon.tuId AND Operon.posId=Position.posId"
	    . " AND Position.scaffoldId=Scaffold.scaffoldId AND Scaffold.taxonomyId=$taxon"
	    . " AND Operon.evidence='VIMSS Predicted'";
	my $nToDel = queryScalar("SELECT COUNT(*) FROM Locus2Operon, Operon, Position, Scaffold WHERE $where");
	print STDERR "# For taxon $taxon, have $nToDel to delete\n";
	print "DELETE FROM Locus2Operon, Operon, Position USING Locus2Operon,Operon,Position,Scaffold WHERE $where;\n";
    }
    print "LOAD DATA INFILE '" . $OperonPosFile . "' INTO TABLE Position (posId,scaffoldId,objectId,strand,begin,end);\n";
    print "LOAD DATA INFILE '" . $OperonFile . "' INTO TABLE Operon (tuId,posId,evidence);\n";
    print "LOAD DATA INFILE '" . $Locus2OperonFile . "' INTO TABLE Locus2Operon (locusId,tuId);\n";
}
