SliceAdaptor.pm 64.5 KB
Newer Older
1
=head1 LICENSE
2

3
  Copyright (c) 1999-2010 The European Bioinformatics Institute and
4
  Genome Research Limited.  All rights reserved.
5

6 7 8 9 10 11 12 13
  This software is distributed under a modified Apache license.
  For license details, please see

    http://www.ensembl.org/info/about/code_licence.html

=head1 CONTACT

  Please email comments or questions to the public Ensembl
14
  developers list at <dev@ensembl.org>.
15 16 17 18 19

  Questions may also be sent to the Ensembl help desk at
  <helpdesk@ensembl.org>.

=cut
20 21 22

=head1 NAME

23 24
Bio::EnsEMBL::DBSQL::SliceAdaptor - A database aware adaptor responsible for
the creation of Slice objects.
25 26 27

=head1 SYNOPSIS

Graham McVicker's avatar
Graham McVicker committed
28
  use Bio::EnsEMBL::Utils::Slice qw(split_Slices);
29 30 31
  use Bio::EnsEMBL::Registry;

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

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

Graham McVicker's avatar
Graham McVicker committed
39
  # get a slice on the entire chromosome X
40
  $chr_slice = $slice_adaptor->fetch_by_region( 'chromosome', 'X' );
Graham McVicker's avatar
Graham McVicker committed
41 42

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

Graham McVicker's avatar
Graham McVicker committed
47
  # get a slice which is part of NT_004321
48 49 50
  $spctg_slice =
    $slice_adaptor->fetch_by_region( 'supercontig', 'NT_004321',
    200_000, 600_000 );
51

52 53
  # get all non-redundant slices from the highest possible coordinate
  # systems
Graham McVicker's avatar
Graham McVicker committed
54 55 56
  $slices = $slice_adaptor->fetch_all('toplevel');

  # include non-reference regions
57
  $slices = $slice_adaptor->fetch_all( 'toplevel', undef, 1 );
Graham McVicker's avatar
Graham McVicker committed
58 59

  # include non-duplicate regions
60
  $slices = $slice_adaptor->fetch_all( 'toplevel', undef, 0, 1 );
Graham McVicker's avatar
Graham McVicker committed
61 62

  # split up a list of slices into smaller slices
63
  $overlap    = 1000;
Graham McVicker's avatar
Graham McVicker committed
64
  $max_length = 1e6;
65
  $slices     = split_Slices( $slices, $max_length, $overlap );
Graham McVicker's avatar
Graham McVicker committed
66 67

  # store a list of slice names in a file
68
  open( FILE, ">$filename" ) or die("Could not open file $filename");
Graham McVicker's avatar
Graham McVicker committed
69 70 71 72 73 74
  foreach my $slice (@$slices) {
    print FILE $slice->name(), "\n";
  }
  close FILE;

  # retreive a list of slices from a file
75 76
  open( FILE, $filename ) or die("Could not open file $filename");
  while ( $name = <FILE> ) {
Graham McVicker's avatar
Graham McVicker committed
77 78 79 80 81
    chomp($name);
    $slice = $slice_adaptor->fetch_by_name($name);
    # do something with slice
  }

82 83
=head1 DESCRIPTION

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

88
=head1 METHODS
89 90 91 92

=cut


93
package Bio::EnsEMBL::DBSQL::SliceAdaptor;
94 95 96
use vars qw(@ISA);
use strict;

97

98 99
use Bio::EnsEMBL::DBSQL::BaseAdaptor;
use Bio::EnsEMBL::Slice;
100
use Bio::EnsEMBL::CircularSlice;
101
use Bio::EnsEMBL::Mapper;
Ian Longden's avatar
Ian Longden committed
102
use Bio::EnsEMBL::LRGSlice;
103
use Bio::EnsEMBL::Utils::Exception qw(throw deprecate warning stack_trace_dump);
104

105

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

108 109
sub new {
  my $caller = shift;
110

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

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

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

118 119 120 121
  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'};
122

Ian Longden's avatar
Ian Longden committed
123 124 125 126

  $self->{'lrg_region_test'} = undef;
  my $meta_container = $self->db->get_MetaContainer();
  my @values = $meta_container->list_value_by_key("LRG");
127
  if(scalar(@values) and $values[0]->[0]){
Ian Longden's avatar
Ian Longden committed
128 129
    $self->{'lrg_region_test'} = $values[0]->[0];
  }
130 131 132 133 134 135
  return $self;
}


=head2 fetch_by_region

136
  Arg [1]    : string $coord_system_name (optional)
137
               The name of the coordinate system of the slice to be created
138
               This may be a name of an actual coordinate system or an alias
139 140 141 142
               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
Graham McVicker's avatar
Graham McVicker committed
143
               created on.
144 145 146 147 148 149 150 151
  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)
152 153
  Arg [7]    : boolean $no_fuzz (optional, default = undef (false))
               If true (non-zero), do not use "fuzzy matching" (see below).
154 155 156
  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
157 158 159 160 161 162 163 164 165 166 167
               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.
168 169 170 171 172 173 174

               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'.

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

178 179 180 181 182 183 184
               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
185
  Caller     : general
186
  Status     : Stable
187 188 189

=cut

190 191 192 193

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

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

202 203 204
  if ( !defined($seq_region_name) ) {
    throw('seq_region_name argument is required');
  }
205

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

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

212
    ## REMOVE THESE THREE LINES WHEN STICKLEBACK DB IS FIXED!
213 214 215 216 217 218 219
    ## 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!

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

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

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

236
  } ## end if ( defined($coord_system_name...))
237 238 239

  my $constraint;
  my $sql;
240
  my @bind_params;
241 242
  my $key;

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

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

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

257 258 259 260 261
    $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) ) {
262
      $constraint .= "AND cs.version = ? ";
263
      push( @bind_params, [ $version, SQL_VARCHAR ] );
264
    }
265

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

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

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

  if ( defined($arr) ) {
276
    $length = $arr->[3];
277
  } else {
278 279 280 281 282 283 284 285
    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] );
286
    }
287 288

    $sth->execute();
289

290
    if ( $sth->rows() == 0 ) {
291
      $sth->finish();
Graham McVicker's avatar
Graham McVicker committed
292

293 294
      if ($no_fuzz) { return undef }

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

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

301 302 303 304 305 306
      $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] );
307
      }
308

309
      $sth->execute();
310 311

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

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

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

Felix Kokocinski's avatar
Felix Kokocinski committed
321
      my $i = 0;
322

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

327
        # cache values for future reference
328
        my $arr = [ $id, $tmp_name, $cs_id, $tmp_length ];
329
        $self->{'sr_name_cache'}->{"$tmp_name:$cs_id"} = $arr;
330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348
        $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;
Graham McVicker's avatar
Graham McVicker committed
349
        }
350 351

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

355
      # warn if fuzzy matching found more than one result
356 357 358 359 360 361 362
      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",
363
            $high_cs->name(), $seq_region_name ) );
364 365
      }

366 367
      $cs = $high_cs;

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

371 372
    } else {

373 374 375
      my ( $id, $cs_id );
      ( $seq_region_name, $id, $length, $cs_id ) =
        $sth->fetchrow_array();
376
      $sth->finish();
377

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

386
  if ( !defined($end) ) { $end = $length }
Graham McVicker's avatar
Graham McVicker committed
387

388
  if ( $end + 1 < $start ) {
389 390 391 392 393 394 395 396 397 398 399
    my $new_sl =
      Bio::EnsEMBL::CircularSlice->new(
                                   -COORD_SYSTEM    => $cs,
                                   -SEQ_REGION_NAME => $seq_region_name,
                                   -SEQ_REGION_LENGTH => $length,
                                   -START             => $start,
                                   -END               => $end,
                                   -STRAND            => 1,
                                   -ADAPTOR           => $self );

    return $new_sl;
400 401
  }

402 403 404
  if ( defined( $self->{'lrg_region_test'} )
       and substr( $cs->name, 0, 3 ) eq $self->{'lrg_region_test'} )
  {
Ian Longden's avatar
Ian Longden committed
405
    return
406 407 408 409 410 411 412 413
      Bio::EnsEMBL::LRGSlice->new( -COORD_SYSTEM    => $cs,
                                   -SEQ_REGION_NAME => $seq_region_name,
                                   -SEQ_REGION_LENGTH => $length,
                                   -START             => $start,
                                   -END               => $end,
                                   -STRAND            => $strand,
                                   -ADAPTOR           => $self );
  } else {
Ian Longden's avatar
Ian Longden committed
414
    return
415 416 417 418 419 420 421 422
      Bio::EnsEMBL::Slice->new_fast( {
                                  'coord_system'    => $cs,
                                  'seq_region_name' => $seq_region_name,
                                  'seq_region_length' => $length,
                                  'start'             => $start,
                                  'end'               => $end,
                                  'strand'            => $strand,
                                  'adaptor'           => $self } );
Ian Longden's avatar
Ian Longden committed
423
  }
424
} ## end sub fetch_by_region
425

426
=head2 fetch_by_region_unique
427

428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449
  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_unique('chromosome', 'HSCHR6_MHC_COX');
  Description: Retrieves a slice on the requested region but returns only the unique
               parts of the slice.  At a minimum the
               name the name of the seq_region to fetch must be provided.
450

451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484
               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 : listref Bio::EnsEMBL::Slice
  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

sub fetch_by_region_unique {
  my $self = shift;

485
  my @out   = ();
486
  my $slice = $self->fetch_by_region(@_);
487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505

  if ( !exists( $self->{'asm_exc_cache'} ) ) {
    $self->_build_exception_cache();
  }

  if ( exists(
          $self->{'asm_exc_cache'}->{ $self->get_seq_region_id($slice) }
       ) )
  {
    # 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] );
506 507 508
      }
    }
  }
509

510
  return \@out;
511
} ## end sub fetch_by_region_unique
512

513
=head2 fetch_by_name
514

Graham McVicker's avatar
Graham McVicker committed
515
  Arg [1]    : string $name
516 517 518 519 520 521 522 523 524
  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.
Graham McVicker's avatar
Graham McVicker committed
525

526 527 528 529
               Returns undef if no seq_region with the provided name exists in
               the database.

  Returntype : Bio::EnsEMBL::Slice or undef
530 531
  Exceptions : throw if incorrent arg provided
  Caller     : Pipeline
532
  Status     : Stable
533 534 535

=cut

536 537 538
sub fetch_by_name {
  my $self = shift;
  my $name = shift;
Graham McVicker's avatar
Graham McVicker committed
539

540 541
  if(!$name) {
    throw("name argument is required");
Graham McVicker's avatar
Graham McVicker committed
542
  }
543 544 545

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

546
  if(scalar(@array) < 3 || scalar(@array) > 6) {
547
    throw("Malformed slice name [$name].  Format is " .
Laura Clarke's avatar
 
Laura Clarke committed
548
        "coord_system:version:name:start:end:strand");
549 550
  }

551
  # Rearrange arguments to suit fetch_by_region
552

553 554 555 556 557 558 559 560 561
  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);
Graham McVicker's avatar
Graham McVicker committed
562
}
Graham McVicker's avatar
Graham McVicker committed
563 564 565



566 567 568 569 570
=head2 fetch_by_seq_region_id

  Arg [1]    : string $seq_region_id
               The internal identifier of the seq_region to create this slice
               on
571 572 573
  Arg [2]    : optional start
  Arg [3]    : optional end
  Arg [4]    : optional strand
574 575 576
  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.
577 578
               Returns undef if no such slice exists.
  Returntype : Bio::EnsEMBL::Slice or undef
Graham McVicker's avatar
Graham McVicker committed
579 580
  Exceptions : none
  Caller     : general
581
  Status     : Stable
Graham McVicker's avatar
Graham McVicker committed
582 583 584

=cut

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

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

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

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

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

605
    ( $name, $cs_id, $length ) = $sth->fetchrow_array();
606 607 608 609 610
    $sth->finish();

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

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

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

617
  return
618 619 620 621 622 623 624 625
    Bio::EnsEMBL::Slice->new_fast({ 
	                      'coord_system'     => $cs,
                              'seq_region_name'  => $name,
                              'seq_region_length'=> $length,
                              'start'            => $start || 1,
                              'end'              => $end || $length,
                              'strand'           => $strand || 1,
                              'adaptor'           => $self} );
626
} ## end sub fetch_by_seq_region_id
Arne Stabenau's avatar
Arne Stabenau committed
627

Graham McVicker's avatar
Graham McVicker committed
628

629 630

=head2 get_seq_region_id
Graham McVicker's avatar
Graham McVicker committed
631

Arne Stabenau's avatar
Arne Stabenau committed
632
  Arg [1]    : Bio::EnsEMBL::Slice $slice
633 634 635 636 637 638 639 640 641 642
               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
643
  Status     : Stable
Arne Stabenau's avatar
Arne Stabenau committed
644 645

=cut
Graham McVicker's avatar
Graham McVicker committed
646

647 648 649
sub get_seq_region_id {
  my $self = shift;
  my $slice = shift;
Graham McVicker's avatar
Graham McVicker committed
650

Ian Longden's avatar
Ian Longden committed
651
  if(!$slice || !ref($slice) || !($slice->isa('Bio::EnsEMBL::Slice') or $slice->isa('Bio::EnsEMBL::LRGSlice'))) {
652 653
    throw('Slice argument is required');
  }
654
  
655
  my $seq_region_name = $slice->seq_region_name();
656
  my $key = $seq_region_name.":".$slice->coord_system->dbID();
657
  my $arr = $self->{'sr_name_cache'}->{"$key"};
658

659 660
  if( $arr ) {
    return $arr->[0];
661 662
  }

663
  my $cs_id = $slice->coord_system->dbID();
664 665 666 667 668 669

  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
670 671 672
  $sth->bind_param(1,"$seq_region_name",SQL_VARCHAR);
  $sth->bind_param(2,$cs_id,SQL_INTEGER);
  $sth->execute();
673 674 675

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

679 680 681 682 683 684
  }

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

  #cache information for future requests
685
  $arr = [ $seq_region_id, $seq_region_name, $cs_id, $length ];
686 687 688 689

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

690
  return $seq_region_id;
Arne Stabenau's avatar
Arne Stabenau committed
691
}
692 693


Arne Stabenau's avatar
Arne Stabenau committed
694

695 696 697 698 699 700 701 702 703
=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
704
  Arg [3]    : bool $include_non_reference (optional)
705 706 707
               If this argument is not provided then only reference slices
               will be returned. If set, both reference and non refeference
               slices will be rerurned.
Graham McVicker's avatar
Graham McVicker committed
708 709
  Arg [4]    : int $include_duplicates (optional)
               If set duplicate regions will be returned.
710
               
711
               NOTE: if you do not use this option and you have a PAR
712 713 714
               (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.
715

716 717 718 719
  Arg[5]     : bool $include_lrg (optional)  (default 0)
               If set lrg regions will be returned aswell.


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

Graham McVicker's avatar
Graham McVicker committed
723 724 725 726 727
               # 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)};
728 729 730 731 732 733

  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.
734 735
               If the coordinate system with the provided name and version
               does not exist an empty list is returned.
736
               If the coordinate system name provided is 'toplevel', all
737 738
               non-redundant toplevel slices are returned (note that any
               coord_system_version argument is ignored in that case).
739 740 741 742

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

743
  Returntype : listref of Bio::EnsEMBL::Slices
744
  Exceptions : none
745
  Caller     : general
746
  Status     : Stable
Web Admin's avatar
Web Admin committed
747 748 749

=cut

750 751 752 753 754
sub fetch_all {
  my $self = shift;
  my $cs_name = shift;
  my $cs_version = shift || '';

755
  my ($include_non_reference, $include_duplicates, $include_lrg) = @_;
Web Admin's avatar
Web Admin committed
756

757 758 759
  #
  # verify existance of requested coord system and get its id
  #
760 761 762 763 764 765
  my $csa       = $self->db->get_CoordSystemAdaptor();
  my $orig_cs   = $csa->fetch_by_name($cs_name, $cs_version);

  return [] if(!$orig_cs);

  my $sth;
766

767
  my %bad_vals=();
768

769 770 771
  #
  # Get a hash of non reference seq regions
  #
772 773 774 775 776 777 778 779 780 781 782
  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 = ?' );

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

786
    my ($seq_region_id);
787 788 789
    $sth2->bind_columns( \$seq_region_id );

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

793 794
    $sth2->finish();
  }
795

796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822
  #
  # if we do not want lrg's then add them to the bad list;
  #
  if (!$include_lrg) {
    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 = "LRG" '
                      . '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();
  }

823
  #
824
  # Retrieve the seq_regions from the database
825
  #
826 827 828 829 830 831 832 833 834 835 836 837 838
  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 );
839 840
    $sth->execute();
  } else {
841
    $sth =
842 843 844 845
      $self->prepare(   'SELECT sr.seq_region_id, sr.name, '
                      . 'sr.length, sr.coord_system_id '
                      . 'FROM seq_region sr '
                      . 'WHERE sr.coord_system_id = ?' );
846 847 848

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

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

854
  my $cache_count = 0;
855 856 857

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

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

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

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

873
        $cache_count++;
874
      }
875

876 877 878 879 880 881 882 883
      my $slice = Bio::EnsEMBL::Sli