Ace::Sequence - Examine ACeDB Sequence Objects

    # open database connection and get an Ace::Object sequence
    use Ace::Sequence;

    $db  = Ace->connect(-host => 'stein.cshl.org',-port => 200005);
    $obj = $db->fetch(Predicted_gene => 'ZK154.3');

    # Wrap it in an Ace::Sequence object 
    $seq = Ace::Sequence->new($obj);

    # Find all the exons
    @exons = $seq->features('exon');

    # Find all the exons predicted by various versions of "genefinder"
    @exons = $seq->features('exon:genefinder.*');

    # Iterate through the exons, printing their start, end and DNA
    for my $exon (@exons) {
      print join "\t",$exon->start,$exon->end,$exon->dna,"\n";
    }

    # Find the region 1000 kb upstream of the first exon
    $sub = Ace::Sequence->new(-seq=>$exons[0],
                              -offset=>-1000,-length=>1000);

    # Find all features in that area
    @features = $sub->features;

    # Print its DNA
    print $sub->dna;

    # Create a new Sequence object from the first 500 kb of chromosome 1
    $seq = Ace::Sequence->new(-name=>'CHROMOSOME_I',-db=>$db,
			      -offset=>0,-length=>500_000);

    # Get the GFF dump as a text string
    $gff = $seq->gff;

    # Limit dump to Predicted_genes
    $gff_genes = $seq->gff(-features=>'Predicted_gene');

    # Return a GFF object (using optional GFF.pm module from Sanger)
    $gff_obj = $seq->GFF;

 $seq = Ace::Sequence->new($object);

 $seq = Ace::Sequence->new(-source  => $object,
                           -offset  => $offset,
                           -length  => $length,
			   -refseq  => $reference_sequence);

 $seq = Ace::Sequence->new(-name    => $name,
			   -db      => $db,
                           -offset  => $offset,
                           -length  => $length,
			   -refseq  => $reference_sequence);

  $name = $seq->asString;

  $source = $seq->source_seq;

  $parent = $seq->parent_seq;

  $refseq = $seq->refseq;

  $seq->refseq($new_ref);

  $start = $seq->start;

  $end = $seq->end;

  $offset = $seq->offset;

  $length = $seq->length;

  $length = $seq->abslength;

  $strand = $seq->strand;

  $dna = $seq->dna;

  $name = $seq->name;

  $parent = $seq->parent;

  @children = $seq->get_children();

  @features = $seq->features;
  @features = $seq->features('exon','intron','Predicted_gene');
  @features = $seq->features('exon:GeneFinder','Predicted_gene:hand.*');

  type:subtype

  @features = sort { $a->start <=> $b->start } $seq->features;

  my $list = $seq->feature_list();

  $gff = $seq->gff();
  $gff = $seq->gff(-abs      => 1,
                   -features => ['exon','intron:GeneFinder']);

  $gff_object = $seq->gff;
  $gff_object = $seq->gff(-abs      => 1,
                   -features => ['exon','intron:GeneFinder']);

 $abs = $seq->absolute;
 $abs = $seq->absolute(1);

  $merge = $seq->automerge;
  $seq->automerge(0);

  $db = $seq->db;

package Ace::Sequence;
use strict;

use Carp;
use strict;
use Ace 1.50 qw(:DEFAULT rearrange);
use Ace::Sequence::FeatureList;
use Ace::Sequence::Feature;
use AutoLoader 'AUTOLOAD';
use vars '$VERSION';
my %CACHE;

$VERSION = '1.51';

use constant CACHE => 1;

use overload
  '""'       => 'asString',
  cmp        => 'cmp',
;

# synonym: stop = end
*stop = \&end;
*abs = \&absolute;
*source_seq = \&source;
*source_tag = \&subtype;
*primary_tag = \&type;

my %plusminus = (	 '+' => '-',
		 '-' => '+',
		 '.' => '.');

# internal keys
#    parent    => reference Sequence in "+" strand
#    p_offset  => our start in the parent
#    length    => our length
#    strand    => our strand (+ or -)
#    refseq    => reference Sequence for coordinate system

# object constructor
# usually called like this:
# $seq = Ace::Sequence->new($object);
# but can be called like this:
# $seq = Ace::Sequence->new(-db=>$db,-name=>$name);
# or
# $seq = Ace::Sequence->new(-seq    => $object,
#                           -offset => $offset,
#                           -length => $length,
#                           -ref    => $refseq
#                           );
# $refseq, if provided, will be used to establish the coordinate
# system.  Otherwise the first base pair will be set to 1.
sub new {
  my $pack = shift;
  my ($seq,$start,$end,$offset,$length,$refseq,$db) = 
    rearrange([
	       ['SEQ','SEQUENCE','SOURCE'],
	      'START',
	       ['END','STOP'],
	       ['OFFSET','OFF'],
	       ['LENGTH','LEN'],
	       'REFSEQ',
	       ['DATABASE','DB'],
	      ],@_);

  # Object must have a parent sequence and/or a reference
  # sequence.  In some cases, the parent sequence will be the
  # object itself.  The reference sequence is used to set up
  # the frame of reference for the coordinate system.

  # fetch the sequence object if we don't have it already
  croak "Please provide either a Sequence object or a database and name"
    unless ref($seq) || ($seq && $db);

  # convert start into offset
  $offset = $start - 1 if defined($start) and !defined($offset);

  # convert stop/end into length
  $length = ($end > $start) ? $end - $offset : $end - $offset - 2
    if defined($end) && !defined($length);

  # if just a string is passed, try to fetch a Sequence object
  my $obj = ref($seq) ? $seq : $db->fetch('Sequence'=>$seq);
  unless ($obj) {
    Ace->error("No Sequence named $obj found in database");
    return;
  }

  # get parent coordinates and length of this sequence
  # the parent is an Ace Sequence object in the "+" strand
  my ($parent,$p_offset,$p_length,$strand) = find_parent($obj);
  return unless $parent;

  # handle negative strands
  my $r_strand = $strand;
  my $r_offset = $p_offset;
  $offset ||= 0;
  $offset *= -1 if $strand < 0;

  # handle feature objects
  $offset += $obj->offset if $obj->can('smapped');

  # get source
  my $source = $obj->can('smapped') ? $obj->source : $obj;

  # store the object into our instance variables
  my $self = bless {
		    obj        => $source,
		    offset     => $offset,
		    length     => $length || $p_length,
		    parent     => $parent,
		    p_offset   => $p_offset,
		    refseq     => [$source,$r_offset,$r_strand],
		    strand     => $strand,
		    absolute   => 0,
		    automerge  => 1,
		   },$pack;

  # set the reference sequence
  eval { $self->refseq($refseq) } or return if defined $refseq;

  # wheww!
  return $self;
}

# return the "source" object that the user offset from
sub source {
  $_[0]->{obj};
}

# return the parent object
sub parent { $_[0]->{parent} }

# return the length
#sub length { $_[0]->{length} }
sub length { 
  my $self = shift;
  my ($start,$end) = ($self->start,$self->end);
  return $end - $start + ($end > $start ? 1 : -1);  # for stupid 1-based adjustments
}

sub reversed {  return shift->strand < 0; }

sub automerge {
  my $self = shift;
  my $d = $self->{automerge};
  $self->{automerge} = shift if @_;
  $d;
}

# return reference sequence
sub refseq {
  my $self = shift;
  my $prev = $self->{refseq};
  if (@_) {
    my $refseq = shift;
    my $arrayref;

  BLOCK: {
      last BLOCK unless defined ($refseq);

      if (ref($refseq) && ref($refseq) eq 'ARRAY') {
	$arrayref = $refseq;
	last BLOCK;
      }

      if (ref($refseq) && ($refseq->can('smapped'))) {
	croak "Reference sequence has no common ancestor with sequence"
	  unless $self->parent eq $refseq->parent;
	my ($a,$b,$c) = @{$refseq->{refseq}};
	#	$b += $refseq->offset;
	$b += $refseq->offset;
	$arrayref = [$refseq,$b,$refseq->strand];
	last BLOCK;
      }


      # look up reference sequence in database if we aren't given
      # database object already
      $refseq = $self->db->fetch('Sequence' => $refseq)
	unless $refseq->isa('Ace::Object');
      croak "Invalid reference sequence" unless $refseq;

      # find position of ref sequence in parent strand
      my ($r_parent,$r_offset,$r_length,$r_strand) = find_parent($refseq);
      croak "Reference sequence has no common ancestor with sequence" 
	unless $r_parent eq $self->{parent};

      # set to array reference containing this information
      $arrayref = [$refseq,$r_offset,$r_strand];
    }
    $self->{refseq} = $arrayref;
  }
  return unless $prev;
  return $self->parent if $self->absolute;
  return wantarray ? @{$prev} : $prev->[0];
}

# return strand
sub strand { return $_[0]->{strand} }

# return reference strand
sub r_strand {
  my $self = shift;
  return "+1" if $self->absolute;
  if (my ($ref,$r_offset,$r_strand) = $self->refseq) {
    return $r_strand;
  } else {
    return $self->{strand}
  }
}

sub offset { $_[0]->{offset} }
sub p_offset { $_[0]->{p_offset} }

sub smapped { 1; }
sub type    { 'Sequence' }
sub subtype { }

sub debug {
  my $self = shift;
  my $d = $self->{_debug};
  $self->{_debug} = shift if @_;
  $d;
}

# return the database this sequence is associated with
sub db {
  return Ace->name2db($_[0]->{db} ||= $_[0]->source->db);
}

sub start {
  my ($self,$abs) = @_;
  $abs = $self->absolute unless defined $abs;
  return $self->{p_offset} + $self->{offset} + 1 if $abs;

  if ($self->refseq) {
    my ($ref,$r_offset,$r_strand) = $self->refseq;
    return $r_strand < 0 ? 1 + $r_offset - ($self->{p_offset} + $self->{offset})
                         : 1 + $self->{p_offset} + $self->{offset} - $r_offset;
  }

  else {
    return $self->{offset} +1;
  }

}

sub end { 
  my ($self,$abs) = @_;
  my $start = $self->start($abs);
  my $f = $self->{length} > 0 ? 1 : -1;  # for stupid 1-based adjustments
  if ($abs && $self->refseq ne $self->parent) {
    my $r_strand = $self->r_strand;
    return $start - $self->{length} + $f 
      if $r_strand < 0 or $self->{strand} < 0 or $self->{length} < 0;
    return  $start + $self->{length} - $f
  }
  return  $start + $self->{length} - $f if $self->r_strand eq $self->{strand};
  return  $start - $self->{length} + $f;
}

# turn on absolute coordinates (relative to reference sequence)
sub absolute {
  my $self = shift;
  my $prev = $self->{absolute};
  $self->{absolute} = $_[0] if defined $_[0];
  return $prev;
}

# human readable string (for debugging)
sub asString {
  my $self = shift;
  if ($self->absolute) {
    return join '',$self->parent,'/',$self->start,',',$self->end;
  } elsif (my $ref = $self->refseq){
    my $label = $ref->isa('Ace::Sequence::Feature') ? $ref->info : "$ref";
    return join '',$label,'/',$self->start,',',$self->end;

  } else {
    join '',$self->source,'/',$self->start,',',$self->end;
  }
}

sub cmp {
  my ($self,$arg,$reversed) = @_;
  if (ref($arg) and $arg->isa('Ace::Sequence')) {
    my $cmp = $self->parent cmp $arg->parent 
      || $self->start <=> $arg->start;
    return $reversed ? -$cmp : $cmp;
  }
  my $name = $self->asString;
  return $reversed ? $arg cmp $name : $name cmp $arg;
}

# Return the DNA
sub dna {
  my $self = shift;
  return $self->{dna} if $self->{dna};
  my $raw = $self->_query('seqdna');
  $raw=~s/^>.*\n//m;
  $raw=~s/^\/\/.*//mg;
  $raw=~s/\n//g;
  $raw =~ s/\0+\Z//; # blasted nulls!
  my $effective_strand = $self->end >= $self->start ? '+1' : '-1';
  _complement(\$raw) if $self->r_strand ne $effective_strand;
  return $self->{dna} = $raw;
}

# return a gff file
sub gff {
  my $self = shift;
  my ($abs,$features) = rearrange([['ABS','ABSOLUTE'],'FEATURES'],@_);
  $abs = $self->absolute unless defined $abs;

  # can provide list of feature names, such as 'similarity', or 'all' to get 'em all
  #  !THIS IS BROKEN; IT SHOULD LOOK LIKE FEATURE()!
  my $opt = $self->_feature_filter($features);

  my $gff = $self->_gff($opt);
  warn $gff if $self->debug;

  $self->transformGFF(\$gff) unless $abs;
  return $gff;
}

# return a GFF object using the optional GFF.pm module
sub GFF {
  my $self = shift;
  my ($filter,$converter) = @_;  # anonymous subs
  croak "GFF module not installed" unless require GFF;
  require GFF::Filehandle;

  my @lines = grep !/^\/\//,split "\n",$self->gff(@_);
  local *IN;
  local ($^W) = 0;  # prevent complaint by GFF module
  tie *IN,'GFF::Filehandle',\@lines;
  my $gff = GFF::GeneFeatureSet->new;
  $gff->read(\*IN,$filter,$converter) if $gff;
  return $gff;
}

# Get the features table.  Can filter by type/subtype this way:
# features('similarity:EST','annotation:assembly_tag')
sub features {
  my $self = shift;
  my ($filter,$opt) = $self->_make_filter(@_);

  # get raw gff file
  my $gff = $self->gff(-features=>$opt);

  # turn it into a list of features
  my @features = $self->_make_features($gff,$filter);

  if ($self->automerge) {  # automatic merging
    # fetch out constructed transcripts and clones
    my %types = map {lc($_)=>1} (@$opt,@_);
    if ($types{'transcript'}) {
      push @features,$self->_make_transcripts(\@features);
      @features = grep {$_->type !~ /^(intron|exon)$/ } @features;
    }
    push @features,$self->_make_clones(\@features)      if $types{'clone'};
    if ($types{'similarity'}) {
      my @f = $self->_make_alignments(\@features);
      @features = grep {$_->type ne 'similarity'} @features;
      push @features,@f;
    }
  }

  return wantarray ? @features : \@features;
}

# A little bit more complex - assemble a list of "transcripts"
# consisting of Ace::Sequence::Transcript objects.  These objects
# contain a list of exons and introns.
sub transcripts {
  my $self    = shift;
  my $curated = shift;
  my $ef       = $curated ? "exon:curated"   : "exon";
  my $if       = $curated ? "intron:curated" : "intron";
  my $sf       = $curated ? "Sequence:curated" : "Sequence";
  my @features = $self->features($ef,$if,$sf);
  return unless @features;
  return $self->_make_transcripts(\@features);
}

sub _make_transcripts {
  my $self = shift;
  my $features = shift;

  require Ace::Sequence::Transcript;
  my %transcripts;

  for my $feature (@$features) {
    my $transcript = $feature->info;
    next unless $transcript;
    if ($feature->type =~ /^(exon|intron|cds)$/) {
      my $type = $1;
      push @{$transcripts{$transcript}{$type}},$feature;
    } elsif ($feature->type eq 'Sequence') {
      $transcripts{$transcript}{base} ||= $feature;
    }
  }

  # get rid of transcripts without exons
  foreach (keys %transcripts) {
    delete $transcripts{$_} unless exists $transcripts{$_}{exon}
  }

  # map the rest onto Ace::Sequence::Transcript objects
  return map {Ace::Sequence::Transcript->new($transcripts{$_})} keys %transcripts;
}

# Reassemble clones from clone left and right ends
sub clones {
  my $self = shift;
  my @clones = $self->features('Clone_left_end','Clone_right_end','Sequence');
  my %clones;
  return unless @clones;
  return $self->_make_clones(\@clones);
}

sub _make_clones {
  my $self = shift;
  my $features = shift;

  my (%clones,@canonical_clones);
  my $start_label = $self->strand < 0 ? 'end' : 'start';
  my $end_label   = $self->strand < 0 ? 'start' : 'end';
  for my $feature (@$features) {
    $clones{$feature->info}{$start_label} = $feature->start if $feature->type eq 'Clone_left_end';
    $clones{$feature->info}{$end_label}   = $feature->start if $feature->type eq 'Clone_right_end';

    if ($feature->type eq 'Sequence') {
      my $info = $feature->info;
      next if $info =~ /LINK|CHROMOSOME|\.\w+$/;
      if ($info->Genomic_canonical(0)) {
	push (@canonical_clones,$info->Clone) if $info->Clone;
      }
    }
  }

  foreach (@canonical_clones) {
    $clones{$_} ||= {};
  }

  my @features;
  my ($r,$r_offset,$r_strand) = $self->refseq;
  my $parent = $self->parent;
  my $abs = $self->absolute;
  if ($abs) {
    $r_offset  = 0;
    $r = $parent;
    $r_strand = '+1';
  }

  # BAD HACK ALERT.  WE DON'T KNOW WHERE THE LEFT END OF THE CLONE IS SO WE USE
  # THE MAGIC NUMBER -99_999_999 to mean "off left end" and
  # +99_999_999 to mean "off right end"
  for my $clone (keys %clones) {
    my $start = $clones{$clone}{start} || -99_999_999;
    my $end   = $clones{$clone}{end}   || +99_999_999;
    my $phony_gff = join "\t",($parent,'Clone','structural',$start,$end,'.','.','.',qq(Clone "$clone"));
    push @features,Ace::Sequence::Feature->new($parent,$r,$r_offset,$r_strand,$abs,$phony_gff);
  }
  return @features;
}

# Assemble a list of "GappedAlignment" objects. These objects
# contain a list of aligned segments.
sub alignments {
  my $self    = shift;
  my @subtypes = @_;
  my @types = map { "similarity:\^$_\$" } @subtypes;
  push @types,'similarity' unless @types;
  return $self->features(@types);
}

sub segments {
  my $self = shift;
  return;
}

sub _make_alignments {
  my $self = shift;
  my $features = shift;
  require Ace::Sequence::GappedAlignment;

  my %homol;

  for my $feature (@$features) {
    next unless $feature->type eq 'similarity';
    my $target = $feature->info;
    my $subtype = $feature->subtype;
    push @{$homol{$target,$subtype}},$feature;
  }

  # map onto Ace::Sequence::GappedAlignment objects
  return map {Ace::Sequence::GappedAlignment->new($homol{$_})} keys %homol;
}

# return list of features quickly
sub feature_list {
  my $self = shift;
  return $self->{'feature_list'} if $self->{'feature_list'};
  return unless my $raw = $self->_query('seqfeatures -version 2 -list');
  return $self->{'feature_list'} = Ace::Sequence::FeatureList->new($raw);
}

# transform a GFF file into the coordinate system of the sequence
sub transformGFF {
  my $self = shift;
  my $gff = shift;
  my $parent  = $self->parent;
  my $strand  = $self->{strand};
  my $source  = $self->source;
  my ($ref_source,$ref_offset,$ref_strand)  = $self->refseq;
  $ref_source ||= $source;
  $ref_strand ||= $strand;

  if ($ref_strand > 0) {
    my $o = defined($ref_offset) ? $ref_offset : ($self->p_offset + $self->offset);
    # find anything that looks like a numeric field and subtract offset from it
    $$gff =~ s/(?<!\")\s+(-?\d+)\s+(-?\d+)/"\t" . ($1 - $o) . "\t" . ($2 - $o)/eg;
    $$gff =~ s/^$parent/$source/mg;
    $$gff =~ s/\#\#sequence-region\s+\S+/##sequence-region $ref_source/m;
    $$gff =~ s/FMAP_FEATURES\s+"\S+"/FMAP_FEATURES "$ref_source"/m;
    return;
  } else {  # strand eq '-'
    my $o = defined($ref_offset) ? (2 + $ref_offset) : (2 + $self->p_offset - $self->offset);
    $$gff =~ s/(?<!\")\s+(-?\d+)\s+(-?\d+)\s+([.\d]+)\s+(\S)/join "\t",'',$o-$2,$o-$1,$3,$plusminus{$4}/eg;
    $$gff =~ s/(Target \"[^\"]+\" )(-?\d+) (-?\d+)/$1 $3 $2/g;
    $$gff =~ s/^$parent/$source/mg;
    $$gff =~ s/\#\#sequence-region\s+\S+\s+(-?\d+)\s+(-?\d+)/"##sequence-region $ref_source " . ($o - $2) . ' ' . ($o - $1) . ' (reversed)'/em;
    $$gff =~ s/FMAP_FEATURES\s+"\S+"\s+(-?\d+)\s+(-?\d+)/"FMAP_FEATURES \"$ref_source\" " . ($o - $2) . ' ' . ($o - $1) . ' (reversed)'/em;
  }

}

# return a name for the object
sub name {
  return shift->source_seq->name;
}

# for compatibility with Ace::Sequence::Feature
sub info {
  return shift->source_seq;
}

###################### internal functions #################
# not necessarily object-oriented!!

# return parent, parent offset and strand
sub find_parent {
  my $obj = shift;

  # first, if we are passed an Ace::Sequence, then we can inherit
  # these settings directly
  return (@{$obj}{qw(parent p_offset length)},$obj->r_strand)
    if $obj->isa('Ace::Sequence');

  # otherwise, if we are passed an Ace::Object, then we must
  # traverse upwards until we find a suitable parent
  return _traverse($obj) if $obj->isa('Ace::Object');

  # otherwise, we don't know what to do...
  croak "Source sequence not an Ace::Object or an Ace::Sequence";
}

sub _get_parent {
  my $obj = shift;
  # ** DANGER DANGER WILL ROBINSON! **
  # This is an experiment in caching parents to speed lookups.  Probably eats memory voraciously.
  return $CACHE{$obj} if CACHE && exists $CACHE{$obj};
  my $p = $obj->get(S_Parent=>2)|| $obj->get(Source=>1);
  return unless $p;
  return CACHE ? $CACHE{$obj} = $p->fetch 
               : $p->fetch;
}

sub _get_children {
  my $obj = shift;
  my @pieces = $obj->get(S_Child=>2);
  return @pieces if @pieces;
  return @pieces = $obj->get('Subsequence');
}

# get sequence, offset and strand of topmost container
sub _traverse {
  my $obj = shift;
  my ($offset,$length);

  # invoke seqget to find the top-level container for this sequence
  my ($tl,$tl_start,$tl_end) = _get_toplevel($obj);
  $tl_start ||= 0;
  $tl_end ||= 0;

  # make it an object
  $tl = ref($obj)->new(-name=>$tl,-class=>'Sequence',-db=>$obj->db);

  $offset += $tl_start - 1;  # offset to beginning of toplevel
  $length ||= abs($tl_end - $tl_start) + 1;
  my $strand = $tl_start < $tl_end ? +1 : -1;

  return ($tl,$offset,$strand < 0 ? ($length,'-1') : ($length,'+1') ) if $length;
}

sub _get_toplevel {
  my $obj = shift;
  my $class = $obj->class;
  my $name  = $obj->name;

  my $smap = $obj->db->raw_query("gif smap -from $class:$name");
  my ($parent,$pstart,$pstop,$tstart,$tstop,$map_type) = 
    $smap =~ /^SMAP\s+(\S+)\s+(\d+)\s+(\d+)\s+(\d+)\s+(\d+)\s+(.+)/;

  $parent ||= '';
  $parent =~ s/^Sequence://;  # remove this in next version of Acedb
  return ($parent,$pstart,$pstop);
}

# create subroutine that filters GFF files for certain feature types
sub _make_filter {
  my $self = shift;
  my $automerge = $self->automerge;

  # parse out the filter
  my %filter;
  foreach (@_) {
    my ($type,$filter) = split(':',$_,2);
    if ($automerge && lc($type) eq 'transcript') {
      @filter{'exon','intron','Sequence','cds'} = ([undef],[undef],[undef],[undef]);
    } elsif ($automerge && lc($type) eq 'clone') {
      @filter{'Clone_left_end','Clone_right_end','Sequence'} = ([undef],[undef],[undef]);
    } else {
      push @{$filter{$type}},$filter;
    }
  }

  # create pattern-match sub
  my $sub;
  my $promiscuous;  # indicates that there is a subtype without a type

  if (%filter) {
    my $s = "sub { my \@d = split(\"\\t\",\$_[0]);\n";
    for my $type (keys %filter) {
      my $expr;
      my $subtypes = $filter{$type};
      if ($type ne '') {
	for my $st (@$subtypes) {
	  $expr .= defined $st ? "return 1 if \$d[2]=~/$type/i && \$d[1]=~/$st/i;\n"
	                       : "return 1 if \$d[2]=~/$type/i;\n"
	}
      } else {  # no type, only subtypes
	$promiscuous++;
	for my $st (@$subtypes) {
	  next unless defined $st;
	  $expr .= "return 1 if \$d[1]=~/$st/i;\n";
	}
      }
      $s .= $expr;
    }
    $s .= "return;\n }";

    $sub = eval $s;
    croak $@ if $@;
  } else {
    $sub = sub { 1; }
  }
  return ($sub,$promiscuous ? [] : [keys %filter]);
}

# turn a GFF file and a filter into a list of Ace::Sequence::Feature objects
sub _make_features {
  my $self = shift;
  my ($gff,$filter) = @_;

  my ($r,$r_offset,$r_strand) = $self->refseq;
  my $parent = $self->parent;
  my $abs    = $self->absolute;
  if ($abs) {
    $r_offset  = 0;
    $r = $parent;
    $r_strand = '+1';
  }
  my @features = map {Ace::Sequence::Feature->new($parent,$r,$r_offset,$r_strand,$abs,$_)}
                 grep !m@^(?:\#|//)@ && $filter->($_),split("\n",$gff);
}


# low level GFF call, no changing absolute to relative coordinates
sub _gff {
  my $self = shift;
  my ($opt,$db) = @_;
  my $data = $self->_query("seqfeatures -version 2 $opt",$db);
  $data =~ s/\0+\Z//;
  return $data; #blasted nulls!
}

# shortcut for running a gif query
sub _query {
  my $self = shift;
  my $command = shift;
  my $db      = shift || $self->db;

  my $parent = $self->parent;
  my $start = $self->start(1);
  my $end   = $self->end(1);
  ($start,$end) = ($end,$start) if $start > $end;  #flippity floppity

  my $coord   = "-coords $start $end";

  # BAD BAD HACK ALERT - CHECKS THE QUERY THAT IS PASSED DOWN
  # ALSO MAKES THINGS INCOMPATIBLE WITH PRIOR 4.9 servers.
#  my $opt     = $command =~ /seqfeatures/ ? '-nodna' : '';
  my $opt = '-noclip';

  my $query = "gif seqget $parent $opt $coord ; $command";
  warn $query if $self->debug;

  return $db->raw_query("gif seqget $parent $opt $coord ; $command");
}

# utility function -- reverse complement
sub _complement {
  my $dna = shift;
  $$dna =~ tr/GATCgatc/CTAGctag/;
  $$dna = scalar reverse $$dna;
}

sub _feature_filter {
  my $self = shift;
  my $features = shift;
  return '' unless $features;
  my $opt = '';
  $opt = '+feature ' . join('|',@$features) if ref($features) eq 'ARRAY' && @$features;
  $opt = "+feature $features" unless ref $features;
  $opt;
}

1;

__END__