package Gene;
require Exporter;

use strict;
use GenomicsUtils;
use Description;


sub new;

# Constructor for gene object. Returns Gene object.  Parameters should
# be specified in a form locusId => x, $start => ss and so on:
# 
#  scaffoldRef - reference to Scaffold object, 
#  begin, end, strand, locusId, version, posId, cogInfoId
#--------------------------------------------------------------
# The following functions return basic info about genes and are
# relatively self-explanatory.

sub start();
sub stop(); 
sub begin();
sub end();
# start/stop returns the position of the first/last bp in the
# start/stop codon for each Locus.
# begin/end returns the first and last bp in the gene (begin<end)
# Note: in the DB, begin and end are stored but they are named start
# and stop (@#%$&!)

sub length(); # in aa not bases
sub strand(); # Returns '+' or '-'
sub locusId();
sub posId();
sub scaffoldId();
sub taxonomyId();
sub taxonomyName();
sub taxonomyShort();
sub scaffoldInfo();
sub toString(); # Returns some information about the gene.
sub gff; #returns a brief description of the gene in GFF format
sub annotation();
sub annotationRef();# Returns list/reference to a list of annotation objects for the gene
sub synonym(); # Returns reference to a list of synonyms for the gene. Alternatively, synonym(x) returns just the synonym of type x
sub cogFunCode();# Returns Riley funtional code for COGs
sub cogDesc(); #Returns the description of COG
sub cog();
sub ec(); #returns a ref to array of EC numbers
sub go(); #returns a ref to array of GO terms
sub ipr(); #returns a ref to array of InterPro domains
sub bestBlastpRef(); #returns the best blastp hit, bestBlastpRef(1) to retrieve best non-paralog hit
sub codonUsageTable();
sub addCodonUsage();
sub isPlasmid();
sub description();
sub descriptionRef(); # Returns list/reference to a list of descriptions [annotations] for the gene.
sub scaffold(); # Returns reference to scaffold object this gene belongs to.

#--------------------------------------------------------------
# Sequence accession methods

sub dna(); # Returns nucleotide sequence of the gene
sub protein(); # Returns amino acid sequence

sub getUpstream; # Returns upstream region for this gene
#parameters:
#downstream - if downstream is true, returns downstream region
#length - how many bps? (default -> 150)
#offset - offset by how many bp upstream(downstream)? (default -> 0)
#         an offset of -3 will return a sequence including the
#         start(stop) codon
#truncate - if truncate is true, remove overlap with adjacent coding regions
#NOTE - this function no longer tries to find the first gene in the operon
#use [$gene->getOperon()]->[0]->getUpstream() instead


#--------------------------------------------------------------
# These functions return info on how this gene is related to others in
# the DB.

sub getOrthologListRef($);
# Only parameter is self

sub getOrthologList($);
# Returns a list of orthologous genes

sub getOrtholog($$);
# Returns the ortholog from specified taxonomyId, or undef

sub getBlastPListRef();
sub getBlastPList();
# Returns list of references to a hash containing field name/field
# value pairs in BlastP table for this gene

sub getBestHitListRef;
sub getBestHitList;

#--------------------------------------------------------------

sub nextGene();
sub previousGene();
# Returns object reference to the next/previous gene along the
# chromosome.

sub upstreamGene();
sub downstreamGene();
# Similar to next and previous, but returns ref to the next gene
# up/downstream on the same strand.  If there are no more adjacent
# genes on the same strand, returns undef.

sub getOperon;
#returns a list of genes known or predicted to be in an operon
#with the query gene (including the query itself)
#if the map parameter is defined (and evaluates to true), then
#a hash of lists is returned, where each list is a single operon
#including the query, and the indexes are tuIds

#--------------------------------------------------------------
# Internal functions used only by the module

sub getDescription($);
# Returns a list of annotation objects based descriptionId
# Parameters:
#   descriptionId   		  

sub setIndex($);
# Sets the index of the gene in the list of genes ordered by begin
# position in the scaffold

sub index();
# Returns index of the gene in a list of genes

sub indexAssigned();
# Returns true if index is assigned to the gene

#--------------------------------------------------------------

my $DATA_DIR = "/data/genomics/data/genomes/";

sub new
{
    my %params = @_;#%{$_[0]};

    die "No locusId declared in Gene::new()" unless defined $params{'locusId'};
    my $self = {};
    bless $self;
    
    #get the best version if none specified
    if (! defined($params{'version'})){
	my $qRef = GenomicsUtils::query( qq{SELECT version FROM Locus WHERE locusId=$params{'locusId'} AND priority=1} );
	$params{'version'} = ${$qRef->[0]}[0];
    }
    $self->{locusId_} = $params{'locusId'};
    $self->{version_} = $params{'version'};

    if (! ( defined($params{'begin'}) &&
	    defined($params{'end'}) &&
	    defined($params{'strand'}) &&
	    defined($params{'posId'}) ) ){
	my $qRef = GenomicsUtils::query( qq{SELECT p.begin, p.end, p.strand, o.posId, c.cogInfoId FROM Position p, Locus o LEFT JOIN COG c ON (o.locusId=c.locusId AND o.version=c.version) WHERE o.locusId=$params{'locusId'} AND o.posId=p.posId AND o.version = $params{'version'}} ) if $params{'version'};

	( $self->{start_},
	  $self->{stop_},
	  $self->{strand_},
	  $self->{posId_},
	  $self->{cog_} ) = @{$qRef->[0]};

	return $self;
    }

    $self->{scaffold_} = $params{'scaffoldRef'};
    $self->{scaffold_} = Scaffold::new(scaffoldId => $params{'scaffoldId'})
        if (!defined $self->{scaffold_} && exists $params{'scaffoldId'});
    $self->{start_} = $params{'begin'};
    $self->{stop_} = $params{'end'};
    $self->{strand_} = $params{'strand'};
    $self->{posId_} = $params{'posId'};
    $self->{cog_} = defined( $params{'cogInfoId'} ) ? $params{'cogInfoId'} : "";
    $self->{index_} = $params{'index'};
    $self->{synonym_} = $params{'synonym'};
    $self->{orthologListRef_} = $params{'orthologListRef'} if exists $params{'orthologListRef'};
    $self->{ec_} = $params{'ec'};
    $self->{go_} = $params{'go'};
    $self->{ipr_} = $params{'ipr'};
    $self->{cogDesc_} = $params{'cogDesc'};
    $self->{cogFun_} = $params{'cogFun'};
    $self->{bestBlastp_} = $params{'bestBlastp'};

    return $self;
}

sub begin() 
{
    my $self = shift;
    return $self->{start_};
}

sub end()
{
    my $self = shift;
    return $self->{stop_};
}

sub start() 
{
    my $self = shift;
    return $self->{strand_} eq "+" ? $self->{start_} : $self->{stop_};
}

sub stop()
{
    my $self = shift;
    return $self->{strand_} eq "+" ? $self->{stop_} : $self->{start_};
}

sub length(){
    my $self = shift;
    return ($self->{stop_} - $self->{start_} + 1)/3;
    &GenomicsUtils::warn("Gene length not multiple of 3 in Gene::length()",2);
}

sub strand() 
{
    my $self = shift;
    return $self->{strand_};
}

sub locusId()
{
    my $self = shift;
    return $self->{locusId_};
}

sub version()
{
    my $self = shift;
    return $self->{version_};
}

sub posId()
{
    my $self = shift;
    return $self->{posId_};
}

sub toString() 
{
    my $self = shift;
    return "[".$self->{locusId_}."] ".$self->{start_}."..".$self->{stop_}." ".$self->{strand_};
}

sub dna() 
{
    my $self = shift;
    if( !defined($self->{dna_} ) ) {
	my $locusId=$self->{locusId_};
	my $version=$self->{version_}; 
	($self->{dna_}) = GenomicsUtils::dbh()->selectrow_array( qq/SELECT sequence FROM LocusSeq WHERE locusId=$locusId and version=$version/ );
    }

    return $self->{dna_};
}

# hash of codons to one-letter amino acid codes
sub codonUsageTable() {
    my @list = (TCA => 'S',TCG => 'S',TCC => 'S',TCT => 'S',
		TTT => 'F',TTC => 'F',TTA => 'L',TTG => 'L',
		TAT => 'Y',TAC => 'Y',TAA => '*',TAG => '*',
		TGT => 'C',TGC => 'C',TGA => '*',TGG => 'W',
		CTA => 'L',CTG => 'L',CTC => 'L',CTT => 'L',
		CCA => 'P',CCG => 'P',CCC => 'P',CCT => 'P',
		CAT => 'H',CAC => 'H',CAA => 'Q',CAG => 'Q',
		CGA => 'R',CGG => 'R',CGC => 'R',CGT => 'R',
		ATT => 'I',ATC => 'I',ATA => 'I',ATG => 'M',
		ACA => 'T',ACG => 'T',ACC => 'T',ACT => 'T',
		AAT => 'N',AAC => 'N',AAA => 'K',AAG => 'K',
		AGT => 'S',AGC => 'S',AGA => 'R',AGG => 'R',
		GTA => 'V',GTG => 'V',GTC => 'V',GTT => 'V',
		GCA => 'A',GCG => 'A',GCC => 'A',GCT => 'A',
		GAT => 'D',GAC => 'D',GAA => 'E',GAG => 'E',
		GGA => 'G',GGG => 'G',GGC => 'G',GGT => 'G');
    return @list;
}

# assumes all of the entries in the hash exist
# returns 0 if the gene has a frameshift preventing codon usage determination
sub addCodonUsage() {
    my ($self, $hash) = @_;
    my $length = $self->end() - $self->begin() + 1;
    return 0 if ($length % 3 != 0);
    my $seq = $self->dna();
    for (my $i = 0; $i < CORE::length($seq); $i += 3) {
	$hash->{substr($seq, $i, 3)}++;
    }
    return 1;
}

sub protein()
{
    my $self = shift;
    if( !defined($self->{protein_} ) ) {
	my $locusId=$self->{locusId_};
	my $version=$self->{version_}; 
	($self->{protein_}) = GenomicsUtils::dbh()->selectrow_array( qq/SELECT sequence FROM AASeq WHERE locusId=$locusId and version=$version/ );
    }

    return $self->{protein_};
}

sub annotation() 
{
    return @{&annotationRef(@_)};
}

sub annotationRef()
{
    my $self = shift;
    if( !defined( $self->{annotation_} ) ) {
	$self->{annotation_} = getDescription( $self->{locusId_} );
    }

    return $self->{annotation_};
}

sub description() 
{
    return @{&descriptionRef(@_)};
}

sub descriptionRef()
{
    my $self = shift;
    my $annotationList = $self->annotationRef();
    my @description;
    foreach my $annotation (@$annotationList) {
	push @description, $annotation->description();
    }
    return \@description;
}

sub synonym() 
{
    my $self = shift;
    my @params;
    if (@_ > 0){
	@params = @_;
    }
    my $locusId = $self->{locusId_};
    my $version = $self->{version_};
    if( !defined( $self->{synonym_} ) ) {
	my $synRef = GenomicsUtils::query( qq{SELECT name,type FROM Synonym WHERE locusId='$locusId' and version=$version} );
	my %syn = map {$_->[1],$_->[0]} @{$synRef};
	$self->{synonym_} = \%syn;

    }
    if (@params){
	my @results = ();
	foreach my $type (@params){
	    if (defined($self->{synonym_}{$type})){
		push @results, $self->{synonym_}{$type};
	    }
	}
	return @results;
    }
    return $self->{synonym_};
}

sub ec(){
    my $self = shift;
    my $locusId = $self->{locusId_}; 
    my $version = $self->{version_};
    if (!defined( $self->{ec_} ) ){
	my $ecRef = GenomicsUtils::query( qq{SELECT Locus2Ec.ecNum, ECInfo.name, Locus2Ec.evidence FROM Locus2Ec LEFT JOIN ECInfo ON Locus2Ec.ecNum=ECInfo.ecNum WHERE locusId='$locusId' AND version=$version} );
	my @ec = ();
	if ($ecRef){
	    foreach (@$ecRef){
		my $evidence = defined($_->[2])?$_->[2]:"KEGG";

		push @ec, $_->[0]."\t".$_->[1]."\t"."(".$evidence.")" if $_->[0];
	    }

	}
	$self->{ec_} = \@ec;
    }
    return $self->{ec_};
	
	
}
sub ipr(){
    my $self = shift;
    my $locusId = $self->{locusId_}; 
    my $version = $self->{version_};
    if (!defined( $self->{ipr_} ) ){
	my $iprRef = GenomicsUtils::query( qq{SELECT distinct(iprId), iprName FROM InterPro WHERE locusId='$locusId' AND iprId <> "NULL"} );
	my @ipr = ();
	if ($iprRef){
	    foreach (@$iprRef){
		push @ipr, $_->[0].": ".$_->[1];
	    }

	}
	$self->{ipr_} = \@ipr;
    }
    return $self->{ipr_};
	
	
}
sub go(){
    my $self = shift;
    my $locusId = $self->{locusId_}; 
    my $version = $self->{version_};
    if (!defined( $self->{go_} ) ){
	my $goRef = GenomicsUtils::query( qq{SELECT t.term_type, t.acc, t.name FROM Locus2Go o, term t WHERE o.locusId='$locusId' AND o.goId=t.id ORDER BY term_type} );
	my @go = ();
	if ($goRef){
	    foreach (@$goRef){
		$_->[0] =~ s/molecular_function/[M]/; 
		$_->[0] =~ s/cellular_component/[C]/;
		$_->[0] =~ s/biological_process/[B]/;
		push @go, $_->[0]." ".$_->[1]." ".$_->[2];
	    }

	}
	$self->{go_} = \@go;
    }
    return $self->{go_};
	
	
}
sub cogFunCode()
{
    my $self = shift; 
    $self->{cogFun_} = "" unless $self->{cog_};
    if( !defined($self->{cogFun_}) ) {
	my $cogInfoId = $self->cog();
	my $funRef = GenomicsUtils::query( qq/SELECT funCode, description from COGInfo WHERE cogInfoId=$cogInfoId/);
	$self->{cogFun_} = defined($funRef->[0][0])?$funRef->[0][0]:"";
	$self->{cogDesc_} = defined($funRef->[0][1])?$funRef->[0][1]:"";
    }
    return $self->{cogFun_};
}
sub cogDesc()
{
    my $self = shift;
    $self->{cogDesc_} = "" unless $self->{cog_};
    
    if( !defined($self->{cogDesc_}) ) {
	my $cogInfoId = $self->cog();
	my $funRef = GenomicsUtils::query( qq/SELECT funCode, description from COGInfo WHERE cogInfoId=$cogInfoId/);
	$self->{cogFun_} = defined($funRef->[0][0])?$funRef->[0][0]:"";
	$self->{cogDesc_} = defined($funRef->[0][1])?$funRef->[0][1]:"";
    }
    return $self->{cogDesc_};
}
sub cog()
{
    my $self = shift;
    if( !defined($self->{cog_}) ) {
	my $locusId=$self->{locusId_}; 
	my $cogRef = GenomicsUtils::query(qq{SELECT c.cogInfoId from Locus o left join COG c on (o.locusId=c.locusId and o.version=c.version) WHERE o.locusId=$locusId });
	$self->{cog_} = defined($cogRef->[0][0])?$cogRef->[0][0]:"";
    }

    return $self->{cog_};
}
sub isPlasmid()
{
    my $self = shift;
    my $isPlasmid = 0;
    my $sId = $self->scaffoldId();
    my $dataR = GenomicsUtils::query(qq{SELECT isGenomic FROM Scaffold WHERE scaffoldId=$sId});
    $isPlasmid = 1 if $dataR->[0][0] == 0;
    return $isPlasmid;

}
sub bestBlastpRef()
{
    my $self = shift;
    my $excludeParalog = @_?$_[0]:0;
    if (!defined($self->{bestBlastp_})){
	my $locusId = $self->{locusId_};
	my $taxId = $self->taxonomyId();
	my $blastRef;
	if ($excludeParalog){
	    $blastRef = GenomicsUtils::query(qq{SELECT b.subject, b.score, b.identity, b.evalue, b.alignLength from BLASTp b, Locus o, Scaffold s where b.locusId=$locusId AND b.subject<>$locusId AND o.locusId=b.subject AND o.priority=1 AND s.scaffoldId=o.scaffoldId AND s.isActive=1 AND s.taxonomyId<>$taxId order by score desc limit 1;});
	    
	}
	else{
	    $blastRef = GenomicsUtils::query(qq{SELECT b.subject, b.score, b.identity, b.evalue, b.alignLength from BLASTp b, Locus o, Scaffold s where b.locusId=$locusId AND b.subject<>$locusId AND o.locusId=b.subject AND o.priority=1 AND s.scaffoldId=o.scaffoldId AND s.isActive=1 order by score desc limit 1;})
	}
	if (defined($blastRef->[0][0])){
	    my $tempObj = Gene::new(locusId=>$blastRef->[0][0]);
	    my %hitInfo = (
			   subject => $tempObj,
			   score => $blastRef->[0][1],
			   identity => $blastRef->[0][2],
			   evalue => $blastRef->[0][3],
			   align => $blastRef->[0][4]
			   );
	    $self->{bestBlastp_} = \%hitInfo;

	}

    }
    return $self->{bestBlastp_};
}

sub getDescription($)
{
    my $locus = shift;
    my $annotations = GenomicsUtils::query( qq{SELECT descriptionId, description, created, source, locusId, version from Description WHERE locusId='$locus' } );

    my @annotations;
    foreach my $refRow (@$annotations) {
	my ($descriptionId, $description, $created, $source, $locusId, $version ) = @$refRow;
	push @annotations, Description::new( $descriptionId, $description, $created, $source, $locusId, $version );
    }

    return \@annotations;
}

sub getOrthologList($){
    return @{&getOrthologListRef(@_)};
}

sub getOrthologListRef($){
    my $self = shift;
    if (!exists $self->{orthologListRef_}) {
	my @orthologs = GenomicsUtils::queryList( qq/ SELECT distinct(locus2) FROM Ortholog WHERE
						  locus1=$self->{locusId_} /  );
	if (scalar @orthologs) {
	    $self->{orthologListRef_} = Scaffold::fetchGenesListRef( locusId => \@orthologs);
	} else {
	    $self->{orthologListRef_} = [];
	}
    }
    return $self->{orthologListRef_};
}

sub getOrtholog($$){
    my ($self, $taxonomyId) = @_;
    my $orthList = $self->getOrthologListRef();

    foreach my $orth (@$orthList){
	return $orth if $orth->scaffold()->taxonomyId() == $taxonomyId;
    }

    return undef;
}


sub scaffold()
{
    my $self = shift;
    if( !defined($self->{scaffold_}) ) {
	my $locusId = $self->locusId();
	my $version = $self->{version_};
	my $scaffoldId = GenomicsUtils::dbh()->selectrow_array( qq/SELECT scaffoldId 
from Locus where locusId=$locusId and version=$version/ );
	$self->{scaffold_} = Scaffold::new( scaffoldId => $scaffoldId );
	print STDERR "Gene.pm(line 396): get scaffoldId for a single locus\n"
	    if GenomicsUtils::getDebug() >= 10;
    }
    return $self->{scaffold_};
}

sub scaffoldId() {
    my $self = shift;
    return $self->scaffold()->scaffoldId();
}
sub scaffoldInfo() {
    my $self = shift;
    return $self->scaffold()->comment();
}
sub taxonomyId(){
    my $self = shift;
    return $self->scaffold()->taxonomyId();
}
sub taxonomyName(){
    my $self = shift;
    return $self->scaffold()->taxonomyName();
}
sub taxonomyShort(){
    my $self = shift;
    return $self->scaffold()->taxonomyShort();
}
sub index()
{
    my $self = shift;
    if( !defined($self->{index_}) ) {
	my $locusId = $self->locusId();
	my $g = $self->scaffold()->genesMapRef()->{$locusId};
	die "Logic error in Gene::index()\n" unless defined( $g->{index_} );
# 	my $genes = $self->scaffold()->genesRef();

# 	my $index = 0;
# 	foreach my $g (@$genes) {
# 	    last if ( $g->locusId()==$locusId );
# 	    ++$index;
# 	}
	$self->{index_} = $g->{index_};
    }
    return $self->{index_};
}

sub setIndex($){
    my $self = shift;
    $self->{index_} = $_[0];
}

sub indexAssigned()
{
    my $self = shift;
    return defined($self->{index_});
}

sub upstreamGene(){
    my $self = shift;
    my $g = ( $self->{strand_} eq '+') ? $self->previousGene() :
	$self->nextGene();

    return undef unless defined $g;
    return ($self->{strand_} eq $g->strand()) ? $g : undef;
}

sub downstreamGene(){
    my $self = shift;
    my $g = ( $self->{strand_} eq '+') ? $self->nextGene() :
	$self->previousGene();

    return undef unless defined $g;
    return ($self->{strand_} eq $g->strand()) ? $g : undef;
}

sub nextGene(){
    my $self = shift;
    my $index = $self->index();
    my $genes = $self->scaffold()->genesRef();
    if( $index<@$genes-1 ) {
	return $genes->[$index+1];
    } 
    if ( $self->scaffold()->isCircular() &&
	 !$self->scaffold()->isPartial() ){
	return $genes->[0];
    }
    return undef;
}

sub previousGene(){
    my $self = shift;
    my $index = $self->index();
    my $genes = $self->scaffold()->genesRef();
    if( $index>0 ) {
	my $gene = $$genes[$index - 1];
	return $gene;
    } 
    if ( $self->scaffold()->isCircular() &&
	 !$self->scaffold()->isPartial() ){
	return $genes->[@$genes-1];
    }	
    return undef;
}

sub getUpstream{
    my $self = shift;
    my %params = @_;
    my $downstream = defined($params{'downstream'}) &&                        
	$params{'downstream'};

    #how many bps?
    my $length = defined($params{'length'}) ?
	$params{'length'} : 150;
    #how far upstream? 0 stops before the first base of the start codon
    my $offset = defined($params{'offset'}) ?
	$params{'offset'} : 0;

    my $plusStrand = ($self->{'strand_'} eq '+');
    #if "downstream" is passed in, treat this like the upstream of the
    #antisense strand (but we won't reverseComplement later)
    $plusStrand = !$plusStrand if $downstream;

    my ($begin,$end);
    if ($plusStrand){
	$begin = $self->begin() - $length - $offset;
    }else{
	$begin = $self->end() + $offset + 1;
    }
    $end = $begin + $length - 1;

    #truncate regions overlapped by other Locuses

    #OK - this ignores the possibility of a gene further
    #(up/down)stream that also overlaps our gene start
    #I don't know if this actually happens, but imagine it
    #could esp. for lousy gene models
    if (defined $params{'truncate'} && $params{'truncate'}){
	if ($plusStrand){
	    my $upGene = $self->previousGene();
	    $begin = ($begin < $upGene->end() + 1) ?
		$upGene->end() + 1 : $begin;
	}else{
	    my $upGene = $self->nextGene();
	    $end = ($end > $upGene->begin() - 1) ?
		$upGene->begin() - 1 : $end;
	}
    }

    return "" if $end < $begin;

    my $seq = $self->scaffold()->sequenceRef();
    #DB indexing starts from 1 not 0
    my $up = substr( $$seq, $begin - 1, $end - $begin + 1 );
    return &GenomicsUtils::reverseComplement( $up ) if ( (!$plusStrand && !$downstream) ||
							 ($plusStrand && $downstream) );
    return $up;
}

sub setSynonym($) 
{
    my $self = shift;
    my $synRef = shift;

    $self->{synonym_} = $synRef;
}

#REMOVED - setScaffold not needed b/c scaffold is preloaded (see Scaffold:new())
#sub setScaffold($) 
#{
#    my $self = shift;
#    my $scaffoldRef = shift;
#
#    $self->{scaffold_} = $scaffoldRef;
#}

sub getBlastPList(){
    return @{&getBlastPListRef(@_)};
}

sub getBlastPListRef(){
    my $self = shift;
    my $locusId = $self->{locusId_};
    my $version = $self->{version_};

    if (!defined $self->{blastp_}){
	my $statement = qq(select b.locusId,b.version,subject,sVersion,identity,alignLength,mismatch,gap,qBegin,qEnd,sBegin,sEnd,evalue,score,s.scaffoldId,taxonomyId from BLASTp b,Scaffold s,Locus o where b.locusId=$locusId and b.version=$version and subject=o.locusId and sVersion=o.version and o.scaffoldId=s.scaffoldId);
	my $sth = GenomicsUtils::dbh()->prepare($statement);
	$sth->execute();
	my @resultSet;
	while( my $row_hash = $sth->fetchrow_hashref() ) {
	    push @resultSet, $row_hash; 
	}
	$sth->finish;
	$self->{blastp_} = \@resultSet;
    }
    return $self->{blastp_};
}

sub getSubjectBlastPList(){
    return @{&getSubjectBlastPListRef(@_)};
}

sub getSubjectBlastPListRef(){
    my $self = shift;
    my $locusId = $self->{locusId_};
    my $version = $self->{version_};

    #add taxId here in query    
    my $statement = qq(select * from BlastP where subject=$locusId and version=$version);
    my $sth = GenomicsUtils::dbh()->prepare($statement);
    $sth->execute();
    my @resultSet;
    while( my $row_hash = $sth->fetchrow_hashref() ) {
	push @resultSet, $row_hash; 
    }
    $sth->finish;
    return \@resultSet;
}

sub getBestHitList{
    return @{&getBestHitListRef(@_)};
}

sub getBestHitListRef{
    my $self = shift;
    my @orthologs = GenomicsUtils::queryList( qq/ SELECT distinct(subjectOrfId) FROM BestHit
					      WHERE queryOrfId=$self->{locusId_} /  );

    return Scaffold::fetchGenesListRef( locusId => \@orthologs, @_ ) if @orthologs;
    &warn("Gene: No BestHits for gene $self->{locusId_}",2);
    return [];
}

sub getOperon{
    my $self = shift;
    my %params = @_;
    my %locusMap;
    my $locusId = $self->locusId();

    #make a hash of all genes in an operon w/ query
    my @locusIds = &queryList( qq{ SELECT o2.locusId FROM Locus2Operon o1, Locus2Operon o2 WHERE o1.locusId=$locusId AND o1.tuId=o2.tuId GROUP BY o2.locusId } );
    my %genes = &Scaffold::fetchGenesMap( locusId => [@locusIds] );

    #return a hash of operons indexed by tuId
    if ((defined $params{map}) && $params{map}){

	#get key,value pairs
	my @TUs = @{&query( qq{ SELECT o2.tuId, o2.locusId FROM Locus2Operon o1, Locus2Operon o2, Position p, Locus o WHERE o1.locusId=$locusId AND o1.tuId=o2.tuId AND o2.locusId=o.locusId AND o.posId=p.posId ORDER BY p.start} )};
	#sort by most upstream
	@TUs = reverse @TUs if $self->strand() eq '-';
	my @locusMap = map {@$_} @TUs;
	my $key;
	my $isKey = 0;

	#make hash using genes already retrieved from DB
	foreach my $val (@locusMap){
	    $isKey = !$isKey;
	    if ($isKey){
		$key = $val;
		next;
	    }
	    push @{$locusMap{$key}}, $genes{$val};
	}
	return %locusMap;
    }

    #sort by most upstream
    my @genes = sort {$a->start() <=> $b->start()} values %genes;
    if ($self->strand() eq '-'){
	return reverse @genes;
    }
    return @genes;
}

sub gff{
    my $self = shift;
    my @features = ( "VIMSS".$self->locusId(), #seqname
		     "VIMSS_DB", #source
		     "gene model", #feature
		     $self->begin(), #start (must be <= end)
		     $self->end(), #end
		     ".", #score
		     $self->strand(), #strand
		     ".", #frame
		     join ";", map { "Synonym$_ "."\"".${$self->synonym()}{$_}."\"" }
		     sort keys %{$self->synonym()} ); #attribute


   return join "\t", @features;
}

sub DESTROY()
{
    
}

1;
