SliceAdaptor.pm 58.3 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 14 15 16 17 18 19
  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
  developers list at <ensembl-dev@ebi.ac.uk>.

  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::Mapper;
Ian Longden's avatar
Ian Longden committed
101
use Bio::EnsEMBL::LRGSlice;
102
use Bio::EnsEMBL::Utils::Exception qw(throw deprecate warning stack_trace_dump);
103

104

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

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

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

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

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

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

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

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


=head2 fetch_by_region

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

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

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

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

=cut

189 190 191 192

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    $sth->execute();
290

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

294 295
      if ($no_fuzz) { return undef }

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

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

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

310
      $sth->execute();
311 312

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

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

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

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

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

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

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

356
      # warn if fuzzy matching found more than one result
357 358 359 360 361 362 363 364 365
      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
          ) );
366 367
      }

368 369
      $cs = $high_cs;

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

373 374
    } else {

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

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

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

390 391 392 393 394 395
  if ( $end + 1 < $start ) {
    throw(
           sprintf(
                  "start [%d] must be less than or equal to end+1 [%d]",
                  $start, $end + 1
           ) );
396 397
  }

Ian Longden's avatar
Ian Longden committed
398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418
  if(defined($self->{'lrg_region_test'}) and substr($cs->name,0,3) eq $self->{'lrg_region_test'}){
    return
      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{
    return
      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} );
  }
419
} ## end sub fetch_by_region
420 421 422



423

424
=head2 fetch_by_name
425

Graham McVicker's avatar
Graham McVicker committed
426
  Arg [1]    : string $name
427 428 429 430 431 432 433 434 435
  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
436

437 438 439 440
               Returns undef if no seq_region with the provided name exists in
               the database.

  Returntype : Bio::EnsEMBL::Slice or undef
441 442
  Exceptions : throw if incorrent arg provided
  Caller     : Pipeline
443
  Status     : Stable
444 445 446

=cut

447 448 449
sub fetch_by_name {
  my $self = shift;
  my $name = shift;
Graham McVicker's avatar
Graham McVicker committed
450

451 452
  if(!$name) {
    throw("name argument is required");
Graham McVicker's avatar
Graham McVicker committed
453
  }
454 455 456

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

457
  if(scalar(@array) < 3 || scalar(@array) > 6) {
458
    throw("Malformed slice name [$name].  Format is " .
Laura Clarke's avatar
 
Laura Clarke committed
459
        "coord_system:version:name:start:end:strand");
460 461
  }

462
  # Rearrange arguments to suit fetch_by_region
463

464 465 466 467 468 469 470 471 472
  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
473
}
Graham McVicker's avatar
Graham McVicker committed
474 475 476



477 478 479 480 481
=head2 fetch_by_seq_region_id

  Arg [1]    : string $seq_region_id
               The internal identifier of the seq_region to create this slice
               on
482 483 484
  Arg [2]    : optional start
  Arg [3]    : optional end
  Arg [4]    : optional strand
485 486 487
  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.
488 489
               Returns undef if no such slice exists.
  Returntype : Bio::EnsEMBL::Slice or undef
Graham McVicker's avatar
Graham McVicker committed
490 491
  Exceptions : none
  Caller     : general
492
  Status     : Stable
Graham McVicker's avatar
Graham McVicker committed
493 494 495

=cut

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

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

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

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

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

516
    ( $name, $cs_id, $length ) = $sth->fetchrow_array();
517 518 519 520 521
    $sth->finish();

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

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

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

528
  return
529 530 531 532 533 534 535 536
    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} );
537
} ## end sub fetch_by_seq_region_id
Arne Stabenau's avatar
Arne Stabenau committed
538

Graham McVicker's avatar
Graham McVicker committed
539

540 541

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

Arne Stabenau's avatar
Arne Stabenau committed
543
  Arg [1]    : Bio::EnsEMBL::Slice $slice
544 545 546 547 548 549 550 551 552 553
               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
554
  Status     : Stable
Arne Stabenau's avatar
Arne Stabenau committed
555 556

=cut
Graham McVicker's avatar
Graham McVicker committed
557

558 559 560
sub get_seq_region_id {
  my $self = shift;
  my $slice = shift;
Graham McVicker's avatar
Graham McVicker committed
561

562
  if(!$slice || !ref($slice) || (!$slice->isa('Bio::EnsEMBL::Slice') && !$slice->isa('Bio::EnsEMBL::LRGSlice'))) {
563 564
    throw('Slice argument is required');
  }
565
  
566
  my $seq_region_name = $slice->seq_region_name();
567
  my $key = $seq_region_name.":".$slice->coord_system->dbID();
568
  my $arr = $self->{'sr_name_cache'}->{"$key"};
569

570 571
  if( $arr ) {
    return $arr->[0];
572 573
  }

574
  my $cs_id = $slice->coord_system->dbID();
575 576 577 578 579 580

  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
581 582 583
  $sth->bind_param(1,"$seq_region_name",SQL_VARCHAR);
  $sth->bind_param(2,$cs_id,SQL_INTEGER);
  $sth->execute();
584 585 586

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

590 591 592 593 594 595
  }

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

  #cache information for future requests
596
  $arr = [ $seq_region_id, $seq_region_name, $cs_id, $length ];
597 598 599 600

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

601
  return $seq_region_id;
Arne Stabenau's avatar
Arne Stabenau committed
602
}
603 604


Arne Stabenau's avatar
Arne Stabenau committed
605

606 607 608 609 610 611 612 613 614
=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
615
  Arg [3]    : bool $include_non_reference (optional)
616 617 618
               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
619 620
  Arg [4]    : int $include_duplicates (optional)
               If set duplicate regions will be returned.
621
               
622
               NOTE: if you do not use this option and you have a PAR
623 624 625
               (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.
626

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

Graham McVicker's avatar
Graham McVicker committed
630 631 632 633 634
               # 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)};
635 636 637 638 639 640

  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.
641 642
               If the coordinate system with the provided name and version
               does not exist an empty list is returned.
643
               If the coordinate system name provided is 'toplevel', all
644 645
               non-redundant toplevel slices are returned (note that any
               coord_system_version argument is ignored in that case).
646 647 648 649

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

650
  Returntype : listref of Bio::EnsEMBL::Slices
651
  Exceptions : none
652
  Caller     : general
653
  Status     : Stable
Web Admin's avatar
Web Admin committed
654 655 656

=cut

657 658 659 660 661
sub fetch_all {
  my $self = shift;
  my $cs_name = shift;
  my $cs_version = shift || '';

662
  my ($include_non_reference, $include_duplicates) = @_;
Web Admin's avatar
Web Admin committed
663

664 665 666
  #
  # verify existance of requested coord system and get its id
  #
667 668 669 670 671 672
  my $csa       = $self->db->get_CoordSystemAdaptor();
  my $orig_cs   = $csa->fetch_by_name($cs_name, $cs_version);

  return [] if(!$orig_cs);

  my $sth;
673

674
  my %bad_vals=();
675

676 677 678
  #
  # Get a hash of non reference seq regions
  #
679 680 681 682 683 684 685 686 687 688 689
  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 = ?' );

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

693
    my ($seq_region_id);
694 695 696
    $sth2->bind_columns( \$seq_region_id );

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

700 701
    $sth2->finish();
  }
702

703
  #
704
  # Retrieve the seq_regions from the database
705
  #
706 707 708 709 710 711 712 713 714 715 716 717 718
  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 );
719 720
    $sth->execute();
  } else {
721
    $sth =
722 723 724 725
      $self->prepare(   'SELECT sr.seq_region_id, sr.name, '
                      . 'sr.length, sr.coord_system_id '
                      . 'FROM seq_region sr '
                      . 'WHERE sr.coord_system_id = ?' );
726 727 728

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

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

734
  my $cache_count = 0;
735 736 737

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

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

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

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

753
        $cache_count++;
754
      }
755

756 757 758 759 760 761 762 763
      my $slice = Bio::EnsEMBL::Slice->new_fast({
	  'start'           => 1,
          'end'             => $length,
          'strand'          => 1,
         'seq_region_name'  => $name,
         'seq_region_length'=> $length,
         'coord_system'     => $cs,
         'adaptor'          => $self});
764 765

      if(!defined($include_duplicates) or !$include_duplicates){
766 767 768 769 770 771 772 773 774 775 776 777
        # 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];
            }
778
          }
779 780 781
        } else {
          # no duplicate regions
          push @out, $slice;
782
        }
783 784
      } else {
        # we want duplicates anyway so do not do any checks
785
        push @out, $slice;
786
      }
787
    }
788
  }
789
  $sth->finish();
790

791
  return \@out;
Web Admin's avatar
Web Admin committed
792 793
}

794 795 796 797 798 799 800 801 802 803
=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

804
sub is_toplevel {
805
  my $self = shift;
806
  my $id   = shift;
Web Admin's avatar
Web Admin committed
807

808 809 810 811 812 813 814
  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 );
815
  $sth->execute();
816 817 818 819 820

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

  while ( $sth->fetch ) {
Ian Longden's avatar
Ian Longden committed
821
    $sth->finish;
822 823
    return 1;
  }
824

Ian Longden's avatar
Ian Longden committed
825
  $sth->finish;
826 827
  return 0;
}
828

Web Admin's avatar
Web Admin committed
829 830 831 832
=head2 fetch_by_band

 Title   : fetch_by_band
 Usage   :
833
 Function: Does not work please use fetch_by_chr_band
Web Admin's avatar
Web Admin committed
834
 Example :
835
 Returns : Bio::EnsEMBL::Slice
Web Admin's avatar
Web Admin committed
836
 Args    : the band name
837
 Status     : AT RISK
Web Admin's avatar
Web Admin committed
838 839 840 841 842 843

=cut

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

844
  my $sth = $self->dbc->prepare
Web Admin's avatar
Web Admin committed
845 846 847 848
        ("select s.name,max(k.seq_region_id)-min(k.seq_region_id, min(k.seq_region_start), max(k.seq_region_id) " .
         "