#!/usr/bin/perl -w
#list of adjacent genes across 1 genome or all genomes

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

{
    my $usage = "GeneNeighbors.pl [-genome=Genome -minScaffoldK=300 -host=host]\n" .
	"GeneNeighbors prints to STDOUT pairs of adjacent genes from the genome or from all genomes\n" .
	"Output file is tab-delimited with columns\n" .
	"Direction is Convergent, Divergent, or Same\n";
    
    my $genomeName = "";
    my $minScaffoldK = 300;
    my $host="localhost";
    (GetOptions('genome=s' => \$genomeName,
		'minScaffoldK=i'=> \$minScaffoldK,
		'host=s' => \$host)
     && @ARGV == 0) || die $usage;
    connect($host);
    my @genomes = ($genomeName eq "") ? Genome::getGenomeList()
	: ( Genome::new("taxonomyId" => $genomeName) );
    print STDOUT join("\t", qw/Genome Scaffold Locus1 Begin1 End1 Strand1
		      Locus2 Begin2 End2 Strand2 Direction Separation/),"\n";
    foreach my $genome (@genomes) {
	my @genes = $genome->genes();
	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 ($separation, $code) = GeneMethods::geneSeparation($gene,$next);
	    print join("\t", $genome->taxonomyId(), $scaffold->scaffoldId(),
		       $gene->locusId(), $gene->begin(), $gene->end(), $gene->strand(),
		       $next->locusId(), $next->begin(), $next->end(), $next->strand(),
		       $code, $separation ),"\n";
	}
    }
    disconnect();
}

