view release on metacpan or  search on metacpan

Ace/Sequence.pm  view on Meta::CPAN

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;

=head1 NAME

Ace::Sequence - Examine ACeDB Sequence Objects

=head1 SYNOPSIS

    # 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;

=head1 DESCRIPTION

I<Ace::Sequence>, and its allied classes L<Ace::Sequence::Feature> and
L<Ace::Sequence::FeatureList>, provide a convenient interface to the
ACeDB Sequence classes and the GFF sequence feature file format.

Using this class, you can define a region of the genome by using a
landmark (sequenced clone, link, superlink, predicted gene), an offset
from that landmark, and a distance.  Offsets and distances can be
positive or negative.  This will return an I<Ace::Sequence> object.
Once a region is defined, you may retrieve its DNA sequence, or query
the database for any features that may be contained within this
region.  Features can be returned as objects (using the
I<Ace::Sequence::Feature> class), as GFF text-only dumps, or in the
form of the GFF class defined by the Sanger Centre's GFF.pm module.

This class builds on top of L<Ace> and L<Ace::Object>.  Please see
their manual pages before consulting this one.

=head1 Creating New Ace::Sequence Objects, the new() Method

 $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);

In order to create an I<Ace::Sequence> you will need an active I<Ace>
database accessor.  Sequence regions are defined using a "source"
sequence, an offset, and a length.  Optionally, you may also provide a
"reference sequence" to establish the coordinate system for all
inquiries.  Sequences may be generated from existing I<Ace::Object>