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

3 4
  Copyright (c) 1999-2009 The European Bioinformatics Institute and
  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;
101

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 122 123 124 125 126 127

  return $self;
}


=head2 fetch_by_region

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

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

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

170 171 172 173 174 175 176
               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
177
  Caller     : general
178
  Status     : Stable
179 180 181

=cut

182 183 184 185

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

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

194 195 196
  if ( !defined($seq_region_name) ) {
    throw('seq_region_name argument is required');
  }
197

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

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

204
    ## REMOVE THESE THREE LINES WHEN STICKLEBACK DB IS FIXED!
205 206 207 208 209 210 211
    ## 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!

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

216 217 218 219 220 221
    if ( !defined($cs) ) {
      throw(
             sprintf( "Unknown coordinate system:\n"
                        . "name='%s' version='%s'\n",
                      $coord_system_name, $version
             ) );
222
    }
223

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

230
  } ## end if ( defined($coord_system_name...
231 232 233

  my $constraint;
  my $sql;
234
  my @bind_params;
235 236
  my $key;

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

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

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

251 252 253 254 255
    $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) ) {
256
      $constraint .= "AND cs.version = ? ";
257
      push( @bind_params, [ $version, SQL_VARCHAR ] );
258
    }
259

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

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

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

  if ( defined($arr) ) {
270
    $length = $arr->[3];
271
  } else {
272 273 274 275 276 277 278 279
    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] );
280
    }
281 282

    $sth->execute();
283

284
    if ( $sth->rows() == 0 ) {
285
      $sth->finish();
Graham McVicker's avatar
Graham McVicker committed
286

287 288
      if ($no_fuzz) { return undef }

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

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

295 296 297 298 299 300
      $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] );
301
      }
302

303
      $sth->execute();
304 305

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

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

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

Felix Kokocinski's avatar
Felix Kokocinski committed
315
      my $i = 0;
316

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

321
        # cache values for future reference
322
        my $arr = [ $id, $tmp_name, $cs_id, $tmp_length ];
323
        $self->{'sr_name_cache'}->{"$tmp_name:$cs_id"} = $arr;
324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342
        $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
343
        }
344 345

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

349
      # warn if fuzzy matching found more than one result
350 351 352 353 354 355 356 357 358
      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
          ) );
359 360
      }

361 362
      $cs = $high_cs;

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

366 367
    } else {

368 369 370
      my ( $id, $cs_id );
      ( $seq_region_name, $id, $length, $cs_id ) =
        $sth->fetchrow_array();
371
      $sth->finish();
372

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

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

383 384 385 386 387 388
  if ( $end + 1 < $start ) {
    throw(
           sprintf(
                  "start [%d] must be less than or equal to end+1 [%d]",
                  $start, $end + 1
           ) );
389 390
  }

391 392 393 394 395 396 397 398 399
  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
400 401 402



403

404
=head2 fetch_by_name
405

Graham McVicker's avatar
Graham McVicker committed
406
  Arg [1]    : string $name
407 408 409 410 411 412 413 414 415
  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
416

417 418 419 420
               Returns undef if no seq_region with the provided name exists in
               the database.

  Returntype : Bio::EnsEMBL::Slice or undef
421 422
  Exceptions : throw if incorrent arg provided
  Caller     : Pipeline
423
  Status     : Stable
424 425 426

=cut

427 428 429
sub fetch_by_name {
  my $self = shift;
  my $name = shift;
Graham McVicker's avatar
Graham McVicker committed
430

431 432
  if(!$name) {
    throw("name argument is required");
Graham McVicker's avatar
Graham McVicker committed
433
  }
434 435 436

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

437
  if(scalar(@array) < 3 || scalar(@array) > 6) {
438
    throw("Malformed slice name [$name].  Format is " .
Laura Clarke's avatar
 
Laura Clarke committed
439
        "coord_system:version:name:start:end:strand");
440 441
  }

442
  # Rearrange arguments to suit fetch_by_region
443

444 445 446 447 448 449 450 451 452
  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
453
}
Graham McVicker's avatar
Graham McVicker committed
454 455 456



457 458 459 460 461
=head2 fetch_by_seq_region_id

  Arg [1]    : string $seq_region_id
               The internal identifier of the seq_region to create this slice
               on
462 463 464
  Arg [2]    : optional start
  Arg [3]    : optional end
  Arg [4]    : optional strand
465 466 467
  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.
468 469
               Returns undef if no such slice exists.
  Returntype : Bio::EnsEMBL::Slice or undef
Graham McVicker's avatar
Graham McVicker committed
470 471
  Exceptions : none
  Caller     : general
472
  Status     : Stable
Graham McVicker's avatar
Graham McVicker committed
473 474 475

=cut

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

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

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

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

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

496
    ( $name, $cs_id, $length ) = $sth->fetchrow_array();
497 498 499 500 501
    $sth->finish();

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

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

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

508 509 510 511 512 513 514 515 516
  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
Arne Stabenau's avatar
Arne Stabenau committed
517

Graham McVicker's avatar
Graham McVicker committed
518

519 520

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

Arne Stabenau's avatar
Arne Stabenau committed
522
  Arg [1]    : Bio::EnsEMBL::Slice $slice
523 524 525 526 527 528 529 530 531 532
               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
533
  Status     : Stable
Arne Stabenau's avatar
Arne Stabenau committed
534 535

=cut
Graham McVicker's avatar
Graham McVicker committed
536

537 538 539
sub get_seq_region_id {
  my $self = shift;
  my $slice = shift;
Graham McVicker's avatar
Graham McVicker committed
540

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

549 550
  if( $arr ) {
    return $arr->[0];
551 552
  }

553
  my $cs_id = $slice->coord_system->dbID();
554 555 556 557 558 559

  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
560 561 562
  $sth->bind_param(1,"$seq_region_name",SQL_VARCHAR);
  $sth->bind_param(2,$cs_id,SQL_INTEGER);
  $sth->execute();
563 564 565

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

569 570 571 572 573 574
  }

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

  #cache information for future requests
575
  $arr = [ $seq_region_id, $seq_region_name, $cs_id, $length ];
576 577 578 579

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

580
  return $seq_region_id;
Arne Stabenau's avatar
Arne Stabenau committed
581
}
582 583


Arne Stabenau's avatar
Arne Stabenau committed
584

585 586 587 588 589 590 591 592 593
=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
594
  Arg [3]    : bool $include_non_reference (optional)
595 596 597
               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
598 599
  Arg [4]    : int $include_duplicates (optional)
               If set duplicate regions will be returned.
600
               
601
               NOTE: if you do not use this option and you have a PAR
602 603 604
               (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.
605

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

Graham McVicker's avatar
Graham McVicker committed
609 610 611 612 613
               # 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)};
614 615 616 617 618 619

  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.
620 621
               If the coordinate system with the provided name and version
               does not exist an empty list is returned.
622
               If the coordinate system name provided is 'toplevel', all
623 624
               non-redundant toplevel slices are returned (note that any
               coord_system_version argument is ignored in that case).
625 626 627 628

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

629
  Returntype : listref of Bio::EnsEMBL::Slices
630
  Exceptions : none
631
  Caller     : general
632
  Status     : Stable
Web Admin's avatar
Web Admin committed
633 634 635

=cut

636 637 638 639 640
sub fetch_all {
  my $self = shift;
  my $cs_name = shift;
  my $cs_version = shift || '';

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

643 644 645
  #
  # verify existance of requested coord system and get its id
  #
646 647 648 649 650 651
  my $csa       = $self->db->get_CoordSystemAdaptor();
  my $orig_cs   = $csa->fetch_by_name($cs_name, $cs_version);

  return [] if(!$orig_cs);

  my $sth;
652

653
  my %bad_vals=();
654

655 656 657
  #
  # Get a hash of non reference seq regions
  #
658 659 660 661 662 663 664 665 666 667 668
  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 = ?' );

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

672
    my ($seq_region_id);
673 674 675
    $sth2->bind_columns( \$seq_region_id );

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

679 680
    $sth2->finish();
  }
681

682
  #
683
  # Retrieve the seq_regions from the database
684
  #
685 686 687 688 689 690 691 692 693 694 695 696 697
  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 );
698 699
    $sth->execute();
  } else {
700
    $sth =
701 702 703 704
      $self->prepare(   'SELECT sr.seq_region_id, sr.name, '
                      . 'sr.length, sr.coord_system_id '
                      . 'FROM seq_region sr '
                      . 'WHERE sr.coord_system_id = ?' );
705 706 707

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

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

713
  my $cache_count = 0;
714 715 716

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

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

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

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

732
        $cache_count++;
733
      }
734 735 736 737 738 739 740 741 742 743 744

      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){
745 746 747 748 749 750 751 752 753 754 755 756
        # 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];
            }
757
          }
758 759 760
        } else {
          # no duplicate regions
          push @out, $slice;
761
        }
762 763
      } else {
        # we want duplicates anyway so do not do any checks
764
        push @out, $slice;
765
      }
766
    }
767
  }
768
  $sth->finish();
769

770
  return \@out;
Web Admin's avatar
Web Admin committed
771 772
}

773 774 775 776 777 778 779 780 781 782
=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

783
sub is_toplevel {
784
  my $self = shift;
785
  my $id   = shift;
Web Admin's avatar
Web Admin committed
786

787 788 789 790 791 792 793
  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 );
794
  $sth->execute();
795 796 797 798 799

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

  while ( $sth->fetch ) {
Ian Longden's avatar
Ian Longden committed
800
    $sth->finish;
801 802
    return 1;
  }
803

Ian Longden's avatar
Ian Longden committed
804
  $sth->finish;
805 806
  return 0;
}
807

Web Admin's avatar
Web Admin committed
808 809 810 811
=head2 fetch_by_band

 Title   : fetch_by_band
 Usage   :
812
 Function: Does not work please use fetch_by_chr_band
Web Admin's avatar
Web Admin committed
813
 Example :
814
 Returns : Bio::EnsEMBL::Slice
Web Admin's avatar
Web Admin committed
815
 Args    : the band name
816
 Status     : AT RISK
Web Admin's avatar
Web Admin committed
817 818 819 820 821 822

=cut

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

823
  my $sth = $self->dbc->prepare
Web Admin's avatar
Web Admin committed
824 825 826 827
        ("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");

828 829
  $sth->bind_param(1,"$band%",SQL_VARCHAR);
  $sth->execute();
Web Admin's avatar
Web Admin committed
830 831 832 833 834 835 836 837 838 839
  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");
}

840 841 842 843 844 845
=head2 fetch_by_chr_band

 Title   : fetch_by_chr_band
 Usage   :
 Function: create a Slice representing a series of bands
 Example :
846
 Returns : Bio::EnsEMBL::Slice
847
 Args    : the band name
848
 Status     : Stable
849 850 851 852

=cut

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

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