#!/usr/bin/perl -w
# compute pseudogenes (good hits to known ORFs).
# In the future, should check to see if the hit is an RNA gene (and label those differently)...

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

{
    my $usage = "MakePseudogenes.pl -taxon=nnnn [-host=host] > FastaFile\n".
	"MakePseudogenes.pl prints to STDOUT the intergenic regions over 50 bases in length\n";
    
    my $taxon = "";
    my $host = undef;
    my $debug = 0;
    (GetOptions('taxon=s' => \$taxon,
		'host=s' => \$host,
		'debug=i' => \$debug)
     && @ARGV ==0
     && $taxon) || die $usage;

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

    my $genome = Genome::new("taxonomyId" => $taxon);
    my @genes = $genome->genes();


    foreach my $gene (@genes) {
	my $scaffold = $gene->scaffold();
	my $scaffoldSeqRef = $scaffold->sequenceRef();
	my $next = $gene->nextGene();
	if($next) {
	    my ($sep, $sepCode) = GeneMethods::geneSeparation($gene,$next);
	    my $left = $gene->end()+1;
	    my $right = $next->begin()-1;
	    my $seq = ($left < $right) ? substr($$scaffoldSeqRef, $left, $right-$left+1 ) : "";
	    $seq = &GenomicsUtils::reverseComplement($seq) if $gene->strand() eq "-";
	    my $title = join(":", $taxon, $scaffold->scaffoldId(), $gene->locusId(), $next->locusId(),
			     $sep, $sepCode,length($seq));
	    print STDOUT GenomicsUtils::formatFASTA($title,$seq) if length($seq) > 50;
	}
    }
}

	    
