SliceAdaptor.pm

#
# Ensembl module for Bio::EnsEMBL::DBSQL::SliceAdaptor
#
#
# Copyright (c) 2004 Ensembl
#
# You may distribute this module under the same terms as perl itself

# POD documentation - main docs before the code

=head1 NAME

Bio::EnsEMBL::DBSQL::SliceAdaptor - A database aware adaptor responsible for
the creation of Slice objects.

=head1 SYNOPSIS

  use Bio::EnsEMBL::Utils::Slice qw(split_Slices);
  use Bio::EnsEMBL::Registry;

  Bio::EnsEMBL::Registry->load_registry_from_db(
              -host => 'ensembldb.ensembl.org',
              -user => 'anonymous'
  );

  $slice_adaptor = Bio::EnsEMBL::Registry->get_adaptor("human", "core", "slice");



  # get a slice on the entire chromosome X
  $chr_slice = $slice_adaptor->fetch_by_region('chromosome','X');

  # get a slice for each clone in the database
  foreach $cln_slice (@{$slice_adaptor->fetch_all('clone')}) {
    #do something with clone
  }

  # get a slice which is part of NT_004321
  $spctg_slice = $slice_adaptor->fetch_by_region('supercontig','NT_004321',
                                                    200_000, 600_000);


  # get all non-redundant slices from the highest possible coordinate systems
  $slices = $slice_adaptor->fetch_all('toplevel');

  # include non-reference regions
  $slices = $slice_adaptor->fetch_all('toplevel',undef,1);

  # include non-duplicate regions
  $slices = $slice_adaptor->fetch_all('toplevel', undef, 0, 1);

  # split up a list of slices into smaller slices
  $overlap = 1000;
  $max_length = 1e6;
  $slices = split_Slices($slices, $max_length, $overlap);


  # store a list of slice names in a file
  open(FILE, ">$filename") or die("Could not open file $filename");
  foreach my $slice (@$slices) {
    print FILE $slice->name(), "\n";
  }
  close FILE;

  # retreive a list of slices from a file
  open(FILE, $filename) or die("Could not open file $filename");
  while($name = <FILE>) {
    chomp($name);
    $slice = $slice_adaptor->fetch_by_name($name);

    # do something with slice
  }

=head1 DESCRIPTION

This module is responsible for fetching Slices representing genomic regions
from a database.  A Details on how slices can be used are in the
Bio::EnsEMBL::Slice module.

=head1 CONTACT

This module is part of the Ensembl project http://www.ensembl.org

For more information email <ensembl-dev@ebi.ac.uk>

=head1 METHODS

=cut


package Bio::EnsEMBL::DBSQL::SliceAdaptor;
use vars qw(@ISA);
use strict;


use Bio::EnsEMBL::DBSQL::BaseAdaptor;
use Bio::EnsEMBL::Slice;
use Bio::EnsEMBL::Mapper;

use Bio::EnsEMBL::Utils::Exception qw(throw deprecate warning stack_trace_dump);


@ISA = ('Bio::EnsEMBL::DBSQL::BaseAdaptor');

sub new {
  my $caller = shift;

  my $class = ref($caller) || $caller;

  my $self = $class->SUPER::new(@_);

  # use a cache which is shared and also used by the assembly
  # mapper adaptor

  my $seq_region_cache = $self->db->get_SeqRegionCache();

  $self->{'sr_name_cache'} = $seq_region_cache->{'name_cache'};
  $self->{'sr_id_cache'}   = $seq_region_cache->{'id_cache'};

  return $self;
}


=head2 fetch_by_region

  Arg [1]    : string $coord_system_name (optional)
               The name of the coordinate system of the slice to be created
               This may be a name of an actual coordinate system or an alias
               to a coordinate system.  Valid aliases are 'seqlevel' or
               'toplevel'.
  Arg [2]    : string $seq_region_name
               The name of the sequence region that the slice will be
               created on.
  Arg [3]    : int $start (optional, default = 1)
               The start of the slice on the sequence region
  Arg [4]    : int $end (optional, default = seq_region length)
               The end of the slice on the sequence region
  Arg [5]    : int $strand (optional, default = 1)
               The orientation of the slice on the sequence region
  Arg [6]    : string $version (optional, default = default version)
               The version of the coordinate system to use (e.g. NCBI33)
  Arg [7]    : boolean $no_fuzz (optional, default = undef (false))
               If true (non-zero), do not use "fuzzy matching" (see below).
  Example    : $slice = $slice_adaptor->fetch_by_region('chromosome', 'X');
               $slice = $slice_adaptor->fetch_by_region('clone', 'AC008066.4');
  Description: Retrieves a slice on the requested region.  At a minimum the
               name the name of the seq_region to fetch must be provided.

               If no coordinate system name is provided than a slice on the
               highest ranked coordinate system with a matching
               seq_region_name will be returned.  If a version but no
               coordinate system name is provided, the same behaviour will
               apply, but only coordinate systems of the appropriate version
               are considered.  The same applies if the 'toplevel' coordinate
               system is specified, however in this case the version is
               ignored.  The coordinate system should always be specified if
               it is known, since this is unambiguous and faster.

               Some fuzzy matching is performed if no exact match for
               the provided name is found.  This allows clones to be
               fetched even when their version is not known.  For
               example fetch_by_region('clone', 'AC008066') will
               retrieve the sequence_region with name 'AC008066.4'.

               The fuzzy matching can be turned off by setting the
               $no_fuzz argument to a true value.

               If the requested seq_region is not found in the database undef
               is returned.

  Returntype : Bio::EnsEMBL::Slice or undef
  Exceptions : throw if no seq_region_name is provided
               throw if invalid coord_system_name is provided
               throw if start > end is provided
  Caller     : general
  Status     : Stable

=cut


#
# ARNE: This subroutine needs simplification!! 
#
sub fetch_by_region {
  my ( $self, $coord_system_name, $seq_region_name, $start, $end,
       $strand, $version, $no_fuzz )
    = @_;

  if ( !defined($start) )  { $start  = 1 }
  if ( !defined($strand) ) { $strand = 1 }

  if ( !defined($seq_region_name) ) {
    throw('seq_region_name argument is required');
  }

  my $cs;
  my $csa = $self->db->get_CoordSystemAdaptor();

  if ( defined($coord_system_name) ) {
    $cs = $csa->fetch_by_name( $coord_system_name, $version );

    ## REMOVE THESE THREE LINES WHEN STICKLEBACK DB IS FIXED!
    ## Anne/ap5 (2007-10-09):
    # The problem was that the stickleback genebuild called the
    # chromosomes 'groups', which meant they weren't being picked out by
    # the karyotype drawing code.  Apparently they are usually called
    # 'groups' in the stickleback community, even though they really are
    # chromosomes!

    if ( !defined($cs) && $coord_system_name eq 'chromosome' ) {
      $cs = $csa->fetch_by_name( 'group', $version );
    }

    if ( !defined($cs) ) {
      throw(
             sprintf( "Unknown coordinate system:\n"
                        . "name='%s' version='%s'\n",
                      $coord_system_name, $version
             ) );
    }

    # fetching by toplevel is same as fetching w/o name or version
    if ( $cs->is_top_level() ) {
      $cs      = undef;
      $version = undef;
    }

  } ## end if ( defined($coord_system_name...

  my $constraint;
  my $sql;
  my @bind_params;
  my $key;

  if ( defined($cs) ) {
    $sql = sprintf( "SELECT sr.name, sr.seq_region_id, sr.length, %d "
                      . "FROM seq_region sr ",
                    $cs->dbID() );

    $constraint = "AND sr.coord_system_id = ?";
    push( @bind_params, [ $cs->dbID(), SQL_INTEGER ] );

    $key = "$seq_region_name:" . $cs->dbID();
  } else {
    $sql =
      "SELECT sr.name, sr.seq_region_id, sr.length, cs.coord_system_id "
      . "FROM seq_region sr, coord_system cs ";

    $constraint = "AND sr.coord_system_id = cs.coord_system_id "
      . "AND cs.species_id = ? ";
    push( @bind_params, [ $self->species_id(), SQL_INTEGER ] );

    if ( defined($version) ) {
      $constraint .= "AND cs.version = ? ";
      push( @bind_params, [ $version, SQL_VARCHAR ] );
    }

    $constraint .= "ORDER BY cs.rank ASC";
  }

  # check the cache so we only go to the db if necessary
  my $length;
  my $arr;

  if ( defined($key) ) { $arr = $self->{'sr_name_cache'}->{$key} }

  if ( defined($arr) ) {
    $length = $arr->[3];
  } else {
    my $sth =
      $self->prepare( $sql . "WHERE sr.name = ? " . $constraint );

    unshift( @bind_params, [ $seq_region_name, SQL_VARCHAR ] );

    my $pos = 0;
    foreach my $param (@bind_params) {
      $sth->bind_param( ++$pos, $param->[0], $param->[1] );
    }

    $sth->execute();

    if ( $sth->rows() == 0 ) {
      $sth->finish();

      if ($no_fuzz) { return undef }

      # Do fuzzy matching, assuming that we are just missing a version
      # on the end of the seq_region name.

      $sth =
        $self->prepare( $sql . " WHERE sr.name LIKE ? " . $constraint );

      $bind_params[0] =
        [ sprintf( '%s.%%', $seq_region_name ), SQL_VARCHAR ];

      $pos = 0;
      foreach my $param (@bind_params) {
        $sth->bind_param( ++$pos, $param->[0], $param->[1] );
      }

      $sth->execute();

      my $prefix_len = length($seq_region_name) + 1;
      my $high_ver   = undef;
      my $high_cs    = $cs;

      # Find the fuzzy-matched seq_region with the highest postfix
      # (which ought to be a version).

      my ( $tmp_name, $id, $tmp_length, $cs_id );
      $sth->bind_columns( \( $tmp_name, $id, $tmp_length, $cs_id ) );

      my $i = 0;

      while ( $sth->fetch ) {
        my $tmp_cs =
          ( defined($cs) ? $cs : $csa->fetch_by_dbID($cs_id) );

        # cache values for future reference
        my $arr = [ $id, $tmp_name, $cs_id, $tmp_length ];
        $self->{'sr_name_cache'}->{"$tmp_name:$cs_id"} = $arr;
        $self->{'sr_id_cache'}->{"$id"}                = $arr;

        my $tmp_ver = substr( $tmp_name, $prefix_len );

        # skip versions which are non-numeric and apparently not
        # versions
        if ( $tmp_ver !~ /^\d+$/ ) { next }

        # take version with highest num, if two versions match take one
        # with highest ranked coord system (lowest num)
        if ( !defined($high_ver)
          || $tmp_ver > $high_ver
          || ( $tmp_ver == $high_ver && $tmp_cs->rank < $high_cs->rank )
          )
        {
          $seq_region_name = $tmp_name;
          $length          = $tmp_length;
          $high_ver        = $tmp_ver;
          $high_cs         = $tmp_cs;
        }

        $i++;
      } ## end while ( $sth->fetch )
      $sth->finish();

      # warn if fuzzy matching found more than one result
      if ( $i > 1 ) {
        warning(
          sprintf(
            "Fuzzy matching of seq_region_name "
              . "returned more than one result.\n"
              . "You might want to check whether the returned seq_region\n"
              . "(%s:%s) is the one you intended to fetch.\n",
            $high_cs->name(), $seq_region_name
          ) );
      }

      $cs = $high_cs;

      # return if we did not find any appropriate match:
      if ( !defined($high_ver) ) { return undef }

    } else {

      my ( $id, $cs_id );
      ( $seq_region_name, $id, $length, $cs_id ) =
        $sth->fetchrow_array();
      $sth->finish();

      # cache to speed up for future queries
      my $arr = [ $id, $seq_region_name, $cs_id, $length ];
      $self->{'sr_name_cache'}->{"$seq_region_name:$cs_id"} = $arr;
      $self->{'sr_id_cache'}->{"$id"}                       = $arr;
      $cs = $csa->fetch_by_dbID($cs_id);
    }
  } ## end else [ if ( defined($arr) )

  if ( !defined($end) ) { $end = $length }

  if ( $end + 1 < $start ) {
    throw(
           sprintf(
                  "start [%d] must be less than or equal to end+1 [%d]",
                  $start, $end + 1
           ) );
  }

  return
    Bio::EnsEMBL::Slice->new( -COORD_SYSTEM      => $cs,
                              -SEQ_REGION_NAME   => $seq_region_name,
                              -SEQ_REGION_LENGTH => $length,
                              -START             => $start,
                              -END               => $end,
                              -STRAND            => $strand,
                              -ADAPTOR           => $self );
} ## end sub fetch_by_region




=head2 fetch_by_name

  Arg [1]    : string $name
  Example    : $name  = 'chromosome:NCBI34:X:1000000:2000000:1';
               $slice = $slice_adaptor->fetch_by_name($name);
               $slice2 = $slice_adaptor->fetch_by_name($slice3->name());
  Description: Fetches a slice using a slice name (i.e. the value returned by
               the Slice::name method).  This is useful if you wish to 
               store a unique identifier for a slice in a file or database or
               pass a slice over a network.
               Slice::name allows you to serialise/marshall a slice and this
               method allows you to deserialise/unmarshal it.

               Returns undef if no seq_region with the provided name exists in
               the database.

  Returntype : Bio::EnsEMBL::Slice or undef
  Exceptions : throw if incorrent arg provided
  Caller     : Pipeline
  Status     : Stable

=cut

sub fetch_by_name {
  my $self = shift;
  my $name = shift;

  if(!$name) {
    throw("name argument is required");
  }

  my @array = split(/:/,$name);

  if(scalar(@array) < 3 || scalar(@array) > 6) {
    throw("Malformed slice name [$name].  Format is " .
        "coord_system:version:name:start:end:strand");
  }

  # Rearrange arguments to suit fetch_by_region

  my @targetarray;

  $targetarray[0]=$array[0];
  $targetarray[5]=(($array[1]&&$array[1] ne "")?$array[1]:undef);
  $targetarray[1]=(($array[2]&&$array[2] ne "")?$array[2]:undef);
  $targetarray[2]=(($array[3]&&$array[3] ne "")?$array[3]:undef);
  $targetarray[3]=(($array[4]&&$array[4] ne "")?$array[4]:undef);
  $targetarray[4]=(($array[5]&&$array[5] ne "")?$array[5]:undef);
  return $self->fetch_by_region(@targetarray);
}



=head2 fetch_by_seq_region_id

  Arg [1]    : string $seq_region_id
               The internal identifier of the seq_region to create this slice
               on
  Arg [2]    : optional start
  Arg [3]    : optional end
  Arg [4]    : optional strand
  Example    : $slice = $slice_adaptor->fetch_by_seq_region_id(34413);
  Description: Creates a slice object of an entire seq_region using the
               seq_region internal identifier to resolve the seq_region.
               Returns undef if no such slice exists.
  Returntype : Bio::EnsEMBL::Slice or undef
  Exceptions : none
  Caller     : general
  Status     : Stable

=cut

sub fetch_by_seq_region_id {
  my ( $self, $seq_region_id, $start, $end, $strand ) = @_;

  my $arr = $self->{'sr_id_cache'}->{$seq_region_id};
  my ( $name, $length, $cs, $cs_id );

  if ( $arr && defined( $arr->[2] ) ) {
    ( $name, $cs_id, $length ) = ( $arr->[1], $arr->[2], $arr->[3] );
    $cs = $self->db->get_CoordSystemAdaptor->fetch_by_dbID($cs_id);
  } else {
    my $sth =
      $self->prepare(   "SELECT sr.name, sr.coord_system_id, sr.length "
                      . "FROM seq_region sr "
                      . "WHERE sr.seq_region_id = ? " );

    $sth->bind_param( 1, $seq_region_id, SQL_INTEGER );
    $sth->execute();

    if ( $sth->rows() == 0 ) { return undef }

    ( $name, $cs_id, $length ) = $sth->fetchrow_array();
    $sth->finish();

    $cs = $self->db->get_CoordSystemAdaptor->fetch_by_dbID($cs_id);

    #cache results to speed up repeated queries
    my $arr = [ $seq_region_id, $name, $cs_id, $length ];

    $self->{'sr_name_cache'}->{"$name:$cs_id"} = $arr;
    $self->{'sr_id_cache'}->{"$seq_region_id"} = $arr;
  }

  return
    Bio::EnsEMBL::Slice->new( -COORD_SYSTEM      => $cs,
                              -SEQ_REGION_NAME   => $name,
                              -SEQ_REGION_LENGTH => $length,
                              -START             => $start || 1,
                              -END               => $end || $length,
                              -STRAND            => $strand || 1,
                              -ADAPTOR           => $self );
} ## end sub fetch_by_seq_region_id



=head2 get_seq_region_id

  Arg [1]    : Bio::EnsEMBL::Slice $slice
               The slice to fetch a seq_region_id for
  Example    : $srid = $slice_adaptor->get_seq_region_id($slice);
  Description: Retrieves the seq_region id (in this database) given a slice
               Seq region ids are not stored on the slices themselves
               because they are intended to be somewhat database independant
               and seq_region_ids vary accross databases.
  Returntype : int
  Exceptions : throw if the seq_region of the slice is not in the db
               throw if incorrect arg provided
  Caller     : BaseFeatureAdaptor
  Status     : Stable

=cut

sub get_seq_region_id {
  my $self = shift;
  my $slice = shift;

  if(!$slice || !ref($slice) || !$slice->isa('Bio::EnsEMBL::Slice')) {
    throw('Slice argument is required');
  }
  
  my $seq_region_name = $slice->seq_region_name();
  my $key = $seq_region_name.":".$slice->coord_system->dbID();
  my $arr = $self->{'sr_name_cache'}->{"$key"};

  if( $arr ) {
    return $arr->[0];
  }

  my $cs_id = $slice->coord_system->dbID();

  my $sth = $self->prepare("SELECT seq_region_id, length " .
                           "FROM seq_region " .
                           "WHERE name = ? AND coord_system_id = ?");

  #force seq_region_name cast to string so mysql cannot treat as int
  $sth->bind_param(1,"$seq_region_name",SQL_VARCHAR);
  $sth->bind_param(2,$cs_id,SQL_INTEGER);
  $sth->execute();

  if($sth->rows() != 1) {
    throw("Non existant or ambigous seq_region:\n" .
          "  coord_system=[$cs_id],\n" .
          "  name=[$seq_region_name],\n");

  }

  my($seq_region_id, $length) = $sth->fetchrow_array();
  $sth->finish();

  #cache information for future requests
  $arr = [ $seq_region_id, $seq_region_name, $cs_id, $length ];

  $self->{'sr_name_cache'}->{"$seq_region_name:$cs_id"} = $arr;
  $self->{'sr_id_cache'}->{"$seq_region_id"} = $arr;

  return $seq_region_id;
}



=head2 fetch_all

  Arg [1]    : string $coord_system_name
               The name of the coordinate system to retrieve slices of.
               This may be a name of an acutal coordinate system or an alias
               to a coordinate system.  Valid aliases are 'seqlevel' or
               'toplevel'.
  Arg [2]    : string $coord_system_version (optional)
               The version of the coordinate system to retrieve slices of
  Arg [3]    : bool $include_non_reference (optional)
               If this argument is not provided then only reference slices
               will be returned. If set, both reference and non refeference
               slices will be rerurned.
  Arg [4]    : int $include_duplicates (optional)
               If set duplicate regions will be returned.
               
               NOTE: if you do not use this option and you have a PAR
               (pseudo-autosomal region) at the beginning of your seq_region
               then your slice will not start at position 1, so coordinates
               retrieved from this slice might not be what you expected.

  Example    : @chromos = @{$slice_adaptor->fetch_all('chromosome','NCBI33')};
               @contigs = @{$slice_adaptor->fetch_all('contig')};

               # get even non-reference regions
               @slices = @{$slice_adaptor->fetch_all('toplevel',undef,1)};

               # include duplicate regions (such as pseudo autosomal regions)
               @slices = @{$slice_adaptor->fetch_all('toplevel', undef,0,1)};

  Description: Retrieves slices of all seq_regions for a given coordinate
               system.  This is analagous to the methods fetch_all which were
               formerly on the ChromosomeAdaptor, RawContigAdaptor and
               CloneAdaptor classes.  Slices fetched span the entire
               seq_regions and are on the forward strand.
               If the coordinate system with the provided name and version
               does not exist an empty list is returned.
               If the coordinate system name provided is 'toplevel', all
               non-redundant toplevel slices are returned (note that any
               coord_system_version argument is ignored in that case).

               Retrieved slices can be broken into smaller slices using the
               Bio::EnsEMBL::Utils::Slice module.

  Returntype : listref of Bio::EnsEMBL::Slices
  Exceptions : none
  Caller     : general
  Status     : Stable

=cut

sub fetch_all {
  my $self = shift;
  my $cs_name = shift;
  my $cs_version = shift || '';

  my ($include_non_reference, $include_duplicates) = @_;

  #
  # verify existance of requested coord system and get its id
  #
  my $csa       = $self->db->get_CoordSystemAdaptor();
  my $orig_cs   = $csa->fetch_by_name($cs_name, $cs_version);

  return [] if(!$orig_cs);

  my $sth;

  my %bad_vals=();

  #
  # Get a hash of non reference seq regions
  #
  if ( !$include_non_reference ) {
    my $sth2 =
      $self->prepare(   'SELECT sr.seq_region_id '
                      . 'FROM seq_region sr, seq_region_attrib sra, '
                      . 'attrib_type at, coord_system cs '
                      . 'WHERE at.code = "non_ref" '
                      . 'AND sra.seq_region_id = sr.seq_region_id '
                      . 'AND at.attrib_type_id = sra.attrib_type_id '
                      . 'AND sr.coord_system_id = cs.coord_system_id '
                      . 'AND cs.species_id = ?' );

    $sth2->bind_param( 1, $self->species_id(), SQL_INTEGER );
    $sth2->execute();

    my ($seq_region_id);
    $sth2->bind_columns( \$seq_region_id );

    while ( $sth2->fetch() ) {
      $bad_vals{$seq_region_id} = 1;
    }

    $sth2->finish();
  }

  #
  # Retrieve the seq_regions from the database
  #
  if ( $orig_cs->is_top_level() ) {
    $sth =
      $self->prepare(   'SELECT sr.seq_region_id, sr.name, '
                      . 'sr.length, sr.coord_system_id '
                      . 'FROM seq_region sr, seq_region_attrib sra, '
                      . 'attrib_type at, coord_system cs '
                      . 'WHERE at.code = "toplevel" '
                      . 'AND at.attrib_type_id = sra.attrib_type_id '
                      . 'AND sra.seq_region_id = sr.seq_region_id '
                      . 'AND sr.coord_system_id = cs.coord_system_id '
                      . 'AND cs.species_id = ?' );

    $sth->bind_param( 1, $self->species_id(), SQL_INTEGER );
    $sth->execute();
  } else {
    $sth =
      $self->prepare(   'SELECT sr.seq_region_id, sr.name, '
                      . 'sr.length, sr.coord_system_id '
                      . 'FROM seq_region sr '
                      . 'WHERE sr.coord_system_id = ?' );

    $sth->bind_param( 1, $orig_cs->dbID, SQL_INTEGER );
    $sth->execute();
  }

  my ( $seq_region_id, $name, $length, $cs_id );
  $sth->bind_columns( \( $seq_region_id, $name, $length, $cs_id ) );

  my $cache_count = 0;

  my @out;
  while($sth->fetch()) {
    if(!defined($bad_vals{$seq_region_id})){
      my $cs = $csa->fetch_by_dbID($cs_id);

      if(!$cs) {
        throw("seq_region $name references non-existent coord_system $cs_id.");
      }

      #cache values for future reference, but stop adding to the cache once we
      #we know we have filled it up
      if($cache_count < $Bio::EnsEMBL::Utils::SeqRegionCache::SEQ_REGION_CACHE_SIZE) {
        my $arr = [ $seq_region_id, $name, $cs_id, $length ];

        $self->{'sr_name_cache'}->{"$name:$cs_id"} = $arr;
        $self->{'sr_id_cache'}->{"$seq_region_id"} = $arr;

        $cache_count++;
      }

      my $slice = Bio::EnsEMBL::Slice->new
        (-START             => 1,
         -END               => $length,
         -STRAND            => 1,
         -SEQ_REGION_NAME   => $name,
         -SEQ_REGION_LENGTH => $length,
         -COORD_SYSTEM      => $cs,
         -ADAPTOR           => $self);

      if(!defined($include_duplicates) or !$include_duplicates){
        # test if this slice *could* have a duplicate (exception) region
        $self->_build_exception_cache() if(!exists $self->{'asm_exc_cache'});
        if(exists $self->{asm_exc_cache}->{$seq_region_id}) {

          # Dereference symlinked assembly regions.  Take out
          # any regions which are symlinked because these are duplicates
          my @projection = @{$self->fetch_normalized_slice_projection($slice)};
          foreach my $segment ( @projection) {
            if($segment->[2]->seq_region_name() eq $slice->seq_region_name() &&
               $segment->[2]->coord_system->equals($slice->coord_system)) {
              push @out, $segment->[2];
            }
          }
        } else {
          # no duplicate regions
          push @out, $slice;
        }
      } else {
        # we want duplicates anyway so do not do any checks
        push @out, $slice;
      }
    }
  }
  $sth->finish();

  return \@out;
}

=head2 is_toplevel
  Arg        : int seq_region_id 
  Example    : my $top = $slice_adptor->is_toplevel($seq_region_id)
  Description: Returns 1 if slice is a toplevel slice else 0
  Returntype : int
  Caller     : Slice method is_toplevel
  Status     : At Risk

=cut

sub is_toplevel {
  my $self = shift;
  my $id   = shift;

  my $sth = $self->prepare(
            "SELECT at.code from seq_region_attrib sra, attrib_type at "
              . "WHERE sra.seq_region_id = ? "
              . "AND at.attrib_type_id = sra.attrib_type_id "
              . "AND at.code = 'toplevel'" );

  $sth->bind_param( 1, $id, SQL_INTEGER );
  $sth->execute();

  my $code;
  $sth->bind_columns( \$code );

  while ( $sth->fetch ) {
    $sth->finish;
    return 1;
  }

  $sth->finish;
  return 0;
}

=head2 fetch_by_band

 Title   : fetch_by_band
 Usage   :
 Function: Does not work please use fetch_by_chr_band
 Example :
 Returns : Bio::EnsEMBL::Slice
 Args    : the band name
 Status     : AT RISK

=cut

sub fetch_by_band {
  my ($self,$band) = @_;

  my $sth = $self->dbc->prepare
        ("select s.name,max(k.seq_region_id)-min(k.seq_region_id, min(k.seq_region_start), max(k.seq_region_id) " .
         "from karyotype as k " .
         "where k.band like ? and k.seq_region_id = s.seq_region_id");

  $sth->bind_param(1,"$band%",SQL_VARCHAR);
  $sth->execute();
  my ( $seq_region_name, $discrepancy, $seq_region_start, $seq_region_end) = $sth->fetchrow_array;

  if($seq_region_name && $discrepancy>0) {
    throw("Band maps to multiple seq_regions");
  } else {
    return $self->fetch_by_region('toplevel',$seq_region_name,$seq_region_start,$seq_region_end);
  }
  throw("Band not recognised in database");
}

=head2 fetch_by_chr_band

 Title   : fetch_by_chr_band
 Usage   :
 Function: create a Slice representing a series of bands
 Example :
 Returns : Bio::EnsEMBL::Slice
 Args    : the band name
 Status     : Stable

=cut

sub fetch_by_chr_band {
  my ( $self, $chr, $band ) = @_;

  my $chr_slice = $self->fetch_by_region( 'toplevel', $chr );
  my $seq_region_id = $self->get_seq_region_id($chr_slice);

  my $sth =
    $self->prepare(   'SELECT MIN(k.seq_region_start), '
                    . 'MAX(k.seq_region_end) '
                    . 'FROM karyotype k '
                    . 'WHERE k.seq_region_id = ? '
                    . 'AND k.band LIKE ?' );

  $sth->bind_param( 1, $seq_region_id, SQL_INTEGER );
  $sth->bind_param( 2, "$band%",       SQL_VARCHAR );
  $sth->execute();

  my ( $slice_start, $slice_end ) = $sth->fetchrow_array;

  if ( defined $slice_start ) {
    return
      $self->fetch_by_region( 'toplevel',   $chr,
                              $slice_start, $slice_end );
  }

  throw("Band not recognised in database");