SliceAdaptor.pm 57.4 KB
Newer Older
1 2

#
Simon Potter's avatar
pod  
Simon Potter committed
3
# Ensembl module for Bio::EnsEMBL::DBSQL::SliceAdaptor
4 5
#
#
Graham McVicker's avatar
Graham McVicker committed
6
# Copyright (c) 2004 Ensembl
7 8 9 10 11 12 13
#
# You may distribute this module under the same terms as perl itself

# POD documentation - main docs before the code

=head1 NAME

14 15
Bio::EnsEMBL::DBSQL::SliceAdaptor - A database aware adaptor responsible for
the creation of Slice objects.
16 17 18

=head1 SYNOPSIS

Graham McVicker's avatar
Graham McVicker committed
19
  use Bio::EnsEMBL::Utils::Slice qw(split_Slices);
20 21 22 23 24 25 26 27
  use Bio::EnsEMBL::Registry;

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

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



Graham McVicker's avatar
Graham McVicker committed
31 32 33 34 35
  # get a slice on the entire chromosome X
  $chr_slice = $slice_adaptor->fetch_by_region('chromosome','X');

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

Graham McVicker's avatar
Graham McVicker committed
39 40
  # get a slice which is part of NT_004321
  $spctg_slice = $slice_adaptor->fetch_by_region('supercontig','NT_004321',
41
                                                    200_000, 600_000);
42 43


Graham McVicker's avatar
Graham McVicker committed
44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74
  # get all non-redundant slices from the highest possible coordinate systems
  $slices = $slice_adaptor->fetch_all('toplevel');

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

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

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


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

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

    # do something with slice
  }

75 76
=head1 DESCRIPTION

77
This module is responsible for fetching Slices representing genomic regions
Graham McVicker's avatar
Graham McVicker committed
78
from a database.  A Details on how slices can be used are in the
79
Bio::EnsEMBL::Slice module.
80

81
=head1 CONTACT
82

Graham McVicker's avatar
Graham McVicker committed
83 84 85
This module is part of the Ensembl project http://www.ensembl.org

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

87
=head1 METHODS
88 89 90 91

=cut


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

96

97 98
use Bio::EnsEMBL::DBSQL::BaseAdaptor;
use Bio::EnsEMBL::Slice;
99
use Bio::EnsEMBL::Mapper;
100

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

103

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

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

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

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

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

116 117 118 119
  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'};
120 121 122 123 124 125 126

  return $self;
}


=head2 fetch_by_region

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

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

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

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

=cut

181 182 183 184

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    $sth->execute();
282

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

286 287
      if ($no_fuzz) { return undef }

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

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

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

302
      $sth->execute();
303 304

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

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

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

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

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

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

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

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

360 361
      $cs = $high_cs;

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

365 366
    } else {

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

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

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

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

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



402

403
=head2 fetch_by_name
404

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

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

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

=cut

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

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

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

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

441
  # Rearrange arguments to suit fetch_by_region
442

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



456 457 458 459 460
=head2 fetch_by_seq_region_id

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

=cut

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

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

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

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

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

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

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

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

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

507 508 509 510 511 512 513 514 515
  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
516

Graham McVicker's avatar
Graham McVicker committed
517

518 519

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

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

=cut
Graham McVicker's avatar
Graham McVicker committed
535

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

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

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

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

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

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

568 569 570 571 572 573
  }

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

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

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

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


Arne Stabenau's avatar
Arne Stabenau committed
583

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

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

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

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

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

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

=cut

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

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

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

  return [] if(!$orig_cs);

  my $sth;
651

652
  my %bad_vals=();
653

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

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

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

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

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

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

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

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

712
  my $cache_count = 0;
713 714 715

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

=cut

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

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

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

839 840 841 842 843 844
=head2 fetch_by_chr_band

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

=cut

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

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

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

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

868 869 870 871 872 873
  my ( $slice_start, $slice_end ) = $sth->fetchrow_array;

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

  throw("Band not recognised in database");
877
} ## end sub fetch_by_chr_band
878 879 880



Web Admin's avatar