#!/usr/bin/perl -w

use lib "/usr2/people/sminovit/Perl";
use GenomicsUtils;
use Genome;
use Scaffold;
use Gene;
use Getopt::Long;
use strict;


my $helpStr=qq{

  Instruction:

    getOperons.060303.pl <scaffoldId> <outFile> <output format> (options)

	<scaffoldId>
	<outFile>
	<output format>:
	        1. for default
		2. for dimerFinder

  Options:

    --gnm or -g: set minimum gnm score, default = 2.5
    --interGenicL or -l: set minimum intergenic length, default = 55

  Description:

	getOperons.060303.pl takes a scaffoldId and returns predicted operons. 
	Currently only gene neighbor method is implemented.

	Fixing bugs from getOperon.pretty

};

unless (@ARGV >= 3){die "$helpStr";}
elsif ($ARGV[0] =~/-/){die "$helpStr";}
elsif ($ARGV[0] !~ /\d+/){print "\n\nScaffold ID: $ARGV[0] not found!!!!\n\n"; die "$helpStr";}

my $scaffoldId = $ARGV[0];
my $outFile = $ARGV[1];
my $outFormat = $ARGV[2];
GenomicsUtils::connect();

#Determine the two parameters: minimum GNM score and minimum Intergenic length
my $minGnm = 0;
my $minIL = 0;
GetOptions( 'gnm|g=i' => \$minGnm,
	    'intergenicLength|l=i' => \$minIL 
);
$minGnm = 2.5 unless $minGnm;
$minIL = 55 unless $minIL;

my $scaffoldR = Scaffold::new($scaffoldId);
my $NC_num = $scaffoldR->file();
my $fileDir = "/usr2/people/alm/data/operons/method/gnm";

#number of operons expected = final count of $directions *2;

my %gnmPairs = ();
my $directions = 0; #keeps track on how many times direction changes

open INFILE, "$fileDir/$scaffoldId.gnm" or die "Can't open $fileDir/$scaffoldId.gnm: $!\n";
while(<INFILE>){
    my ($locusA, $locusB, $consecutive, $interGenicL, $strandA, $strandB, $gnmScore) = split ' ';
    next if (!$consecutive); # only look at a gene pair if the two genes are consecutive
    if (abs($interGenicL)>100000){
	$interGenicL = 100; #arbitrary intergenic length to deal  with those wrapped around genes
    }
    if ( ($strandA ne $strandB) ){
	$gnmPairs{$locusA."=".$locusB} = 0;
	$directions++;
    }
    elsif ( ($gnmScore < $minGnm) && (abs($interGenicL) > $minIL) ){
	$gnmPairs{$locusA."=".$locusB} = 0;
    }
    else{
	$gnmPairs{$locusA."=".$locusB} = 1;
    }
}

close INFILE;

 # The larger $gnmScore is, the more likely the two genes are on the same operon

my %operons = (); #initiate a hash of operons
my %checkedOrf = (); 
my $operonId = 1;

my $genesR = $scaffoldR->genes(); #assuming that genes() return a ref to array of genes that are the best set (highest priority) of the scaffold and ordered by their beginning positions.

for (my $i = 0; $i < @$genesR; $i++){
    if ($i == ( scalar(@$genesR)-1 ) ){
	push @{ $operons{"$operonId"} }, $genesR->[$i]->locusId();
	last;
    }	
   
    push @{ $operons{"$operonId"} }, $genesR->[$i]->locusId();
    my $dsGene = $genesR->[$i+1]->locusId();
    $operonId++ if $gnmPairs{$genesR->[$i]->locusId()."=".$dsGene}==0;
   
}

open OUTFILE, ">$outFile" or die "Can't write to $outFile: $!\n";

if ($outFormat == 1){
    print OUTFILE "\#Arkin Group Operon Prediction\n";
    print OUTFILE "\#gnm score maximum: $minGnm\n";
    print OUTFILE "\#InterGenic length maximum: $minIL\n";
    print OUTFILE "\#$NC_num\n";
    print OUTFILE join "\t", "#Operon", "Begin", "End", "Strand  Info\n";

}
foreach (sort {$a<=>$b} keys %operons){
    my $firstOrfId = $operons{$_}[0];
    my $lastOrfId = $operons{$_}[-1];
    my $firstGene = $scaffoldR->gene($firstOrfId);
    my $lastGene = $scaffoldR->gene($lastOrfId);
    my $locusStrand = $firstGene->strand();
    my $locusBegin = $firstGene->begin();
    my $locusEnd = $lastGene->end();
    print OUTFILE join "\t", ( join "\|", @{ $operons{$_} } ), $locusBegin, $locusEnd, $locusStrand, "Locus".$_."\n";


}

#map{print OUTFILE join "\t", $_, @{$operons{$_}},"\n";}sort {$a<=>$b}keys %operons;
close OUTFILE;
GenomicsUtils::disconnect();
exit;
    
