SliceAdaptor.pm 49.6 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 20
  use Bio::EnsEMBL::DBSQL::DBAdaptor;
  use Bio::EnsEMBL::Utils::Slice qw(split_Slices);
21

Graham McVicker's avatar
Graham McVicker committed
22
  $db = Bio::EnsEMBL::DBSQL::DBAdaptor->new(...);
23

Graham McVicker's avatar
Graham McVicker committed
24
  $slice_adaptor = $db->get_SliceAdaptor();
25

Graham McVicker's avatar
Graham McVicker committed
26 27 28 29 30
  # 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')}) {
31 32 33
    #do something with clone
  }

Graham McVicker's avatar
Graham McVicker committed
34 35
  # get a slice which is part of NT_004321
  $spctg_slice = $slice_adaptor->fetch_by_region('supercontig','NT_004321',
36
                                                    200_000, 600_000);
37 38


Graham McVicker's avatar
Graham McVicker committed
39 40 41 42 43 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
  # 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
  }

70 71
=head1 DESCRIPTION

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

76
=head1 CONTACT
77

Graham McVicker's avatar
Graham McVicker committed
78 79 80
This module is part of the Ensembl project http://www.ensembl.org

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

82
=head1 METHODS
83 84 85 86

=cut


87
package Bio::EnsEMBL::DBSQL::SliceAdaptor;
88 89 90
use vars qw(@ISA);
use strict;

91

92 93
use Bio::EnsEMBL::DBSQL::BaseAdaptor;
use Bio::EnsEMBL::Slice;
94
use Bio::EnsEMBL::Mapper;
95

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

98

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

101 102
sub new {
  my $caller = shift;
103

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

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

108 109 110 111 112 113
  # use a cache which is shared and also used by the assembly
  # mapper adaptor
  my $seq_region_cache = $self->db->get_SeqRegionCache();

  $self->{'sr_name_cache'} = $seq_region_cache->{'name_cache'};
  $self->{'sr_id_cache'}   = $seq_region_cache->{'id_cache'};
114 115 116 117 118 119 120

  return $self;
}


=head2 fetch_by_region

121
  Arg [1]    : string $coord_system_name (optional)
122
               The name of the coordinate system of the slice to be created
123
               This may be a name of an actual coordinate system or an alias
124 125 126 127
               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
128
               created on.
129 130 131 132 133 134 135 136
  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)
137 138 139
  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
140 141 142 143 144 145 146 147 148 149 150
               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.
151 152 153 154 155 156 157 158 159 160 161 162 163 164

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

               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
165
  Caller     : general
166
  Status     : Stable
167 168 169

=cut

170 171 172 173

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

178 179
  $start  = 1 if (!defined($start));
  $strand = 1 if (!defined($strand));
180

181
  throw('seq_region_name argument is required') if (!defined($seq_region_name));
182

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

186
  if ($coord_system_name) {
187 188
    $cs = $csa->fetch_by_name($coord_system_name,$version);

189
    ## REMOVE THESE THREE LINES WHEN STICKLEBACK DB IS FIXED!
190 191 192 193
    if( !$cs && $coord_system_name eq 'chromosome' ) {
      $cs = $csa->fetch_by_name('group',$version);
    }

194
    if (!$cs) {
195 196 197 198 199
      throw("Unknown coordinate system:\n name='$coord_system_name' " .
            "version='$version'\n");
    }

    # fetching by toplevel is same as fetching w/o name or version
200
    if ($cs->is_top_level()) {
201 202 203
      $cs = undef;
      $version = undef;
    }
204
  }
205 206 207 208 209 210

  my $constraint;
  my $sql;
  my @bind_vals;
  my $key;

211
  if ($cs) {
212
    push @bind_vals, $cs->dbID();
213
    $sql = "SELECT sr.name, sr.seq_region_id, sr.length, " .
214 215 216 217
           $cs->dbID() ." FROM seq_region sr ";

    $constraint = "sr.coord_system_id = ?";

218
    $key = "$seq_region_name:".$cs->dbID();
219 220 221 222 223 224 225 226 227 228
  } else {
    $sql = "SELECT sr.name, sr.seq_region_id, sr.length, " .
           "       cs.coord_system_id " .
           "FROM   seq_region sr, coord_system cs ";

    $constraint = "sr.coord_system_id = cs.coord_system_id ";
    if($version) {
      $constraint .= "AND cs.version = ? ";
      push @bind_vals, $version;
    }
229
    $constraint .= "ORDER BY cs.rank ASC";
230 231
  }

232
  # check the cache so we only go to the db if necessary
233
  my $length;
234
  my $arr;
235
  if ($key) {
236
    $arr = $self->{'sr_name_cache'}->{$key};
237
  }
238

239
  if ($arr) {
240
    $length = $arr->[3];
241
  } else {
242
    my $sth = $self->prepare($sql . " WHERE sr.name = ? AND " . $constraint);
243 244 245
    # Quotes around "$seq_region_name" are needed so that mysql does not
    # treat chromosomes like '6' as an int.  This was doing horrible
    # inexact matches like '6DR51', '6_UN', etc.
246
    $sth->bind_param(1, "$seq_region_name", SQL_VARCHAR);
247
    if ($cs){
248
	$sth->bind_param(2, $cs->dbID, SQL_INTEGER);
249 250
    }
    else{
251
	$sth->bind_param(2, $version, SQL_VARCHAR)  if ($version);	
252
    }
253 254

    $sth->execute();
255

256
    if ($sth->rows() == 0) {
257
      $sth->finish();
Graham McVicker's avatar
Graham McVicker committed
258

259 260
      # do fuzzy matching, assuming that we are just missing a version on
      # the end of the seq_region name
261
      
262
      $sth = $self->prepare($sql . " WHERE sr.name LIKE ? AND " . $constraint);
263

264
      $sth->bind_param(1, "$seq_region_name.%", SQL_VARCHAR);
265 266

      if ($cs){
267
	  $sth->bind_param(2, $cs->dbID, SQL_INTEGER);
268 269
      }
      else{
270
	  $sth->bind_param(2, $version, SQL_VARCHAR) if ($version);	
271
      }
272
      $sth->execute();
273 274

      my $prefix_len = length($seq_region_name) + 1;
275 276
      my $high_ver = undef;
      my $high_cs = $cs;
277 278 279 280

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

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

284 285 286
      my $i;

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

289
        # cache values for future reference
290
        my $arr = [ $id, $tmp_name, $cs_id, $tmp_length ];
291 292
        $self->{'sr_name_cache'}->{"$tmp_name:$cs_id"} = $arr;
        $self->{'sr_id_cache'}->{"$id"} = $arr;
Graham McVicker's avatar
Graham McVicker committed
293

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

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

299 300 301 302 303 304 305 306
        # 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
307
        }
308 309

        $i++;
Graham McVicker's avatar
Graham McVicker committed
310
      }
311 312
      $sth->finish();

313 314 315 316 317
      # warn if fuzzy matching found more than one result
      if ($i > 1) {
        warning("Fuzzy matching of seq_region_name returned more than one result.\nYou might want to check whether the returned seq_region\n(".$high_cs->name.":$seq_region_name) is the one you intended to fetch.\n");
      }

318 319
      $cs = $high_cs;

320
      #return if we did not find any appropriate match:
321 322
      return undef if (!defined($high_ver));

323 324
    } else {

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

329 330
      # cahce to speed up for future queries
      my $arr = [ $id, $seq_region_name, $cs_id, $length ];
331 332
      $self->{'sr_name_cache'}->{"$seq_region_name:$cs_id"} = $arr;
      $self->{'sr_id_cache'}->{"$id"} = $arr;
333 334

      $cs = $csa->fetch_by_dbID( $cs_id );
335
    }
336
  }
337

338
  $end = $length if(!defined($end));
Graham McVicker's avatar
Graham McVicker committed
339

340 341
  if($end+1 < $start) {
    throw("start [$start] must be less than or equal to end+1 [$end+1]");
342 343
  }

344
  return Bio::EnsEMBL::Slice->new(-COORD_SYSTEM      => $cs,
345 346 347 348 349 350
                                  -SEQ_REGION_NAME   => $seq_region_name,
                                  -SEQ_REGION_LENGTH => $length,
                                  -START             => $start,
                                  -END               => $end,
                                  -STRAND            => $strand,
                                  -ADAPTOR           => $self);
351 352 353 354
}



355

356
=head2 fetch_by_name
357

Graham McVicker's avatar
Graham McVicker committed
358
  Arg [1]    : string $name
359 360 361 362 363 364 365 366 367
  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
368

369 370 371 372
               Returns undef if no seq_region with the provided name exists in
               the database.

  Returntype : Bio::EnsEMBL::Slice or undef
373 374
  Exceptions : throw if incorrent arg provided
  Caller     : Pipeline
375
  Status     : Stable
376 377 378

=cut

379 380 381
sub fetch_by_name {
  my $self = shift;
  my $name = shift;
Graham McVicker's avatar
Graham McVicker committed
382

383 384
  if(!$name) {
    throw("name argument is required");
Graham McVicker's avatar
Graham McVicker committed
385
  }
386 387 388

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

389
  if(scalar(@array) < 3 || scalar(@array) > 6) {
390
    throw("Malformed slice name [$name].  Format is " .
Laura Clarke's avatar
 
Laura Clarke committed
391
        "coord_system:version:name:start:end:strand");
392 393
  }

394
  # Rearrange arguments to suit fetch_by_region
395

396 397 398 399 400 401 402 403 404
  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
405
}
Graham McVicker's avatar
Graham McVicker committed
406 407 408



409 410 411 412 413 414 415 416
=head2 fetch_by_seq_region_id

  Arg [1]    : string $seq_region_id
               The internal identifier of the seq_region to create this slice
               on
  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.
417 418
               Returns undef if no such slice exists.
  Returntype : Bio::EnsEMBL::Slice or undef
Graham McVicker's avatar
Graham McVicker committed
419 420
  Exceptions : none
  Caller     : general
421
  Status     : Stable
Graham McVicker's avatar
Graham McVicker committed
422 423 424

=cut

425 426 427
sub fetch_by_seq_region_id {
  my ($self, $seq_region_id) = @_;

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

431 432 433 434
  if( $arr ) {
    my $cs_id;
    ($name, $cs_id, $length ) = ( $arr->[1], $arr->[2], $arr->[3] );
    $cs = $self->db->get_CoordSystemAdaptor->fetch_by_dbID($cs_id);
435 436 437 438 439
  } else {
    my $sth = $self->prepare("SELECT name, length, coord_system_id " .
                             "FROM seq_region " .
                             "WHERE seq_region_id = ?");

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

443
    return undef if($sth->rows() == 0);
444 445 446 447 448 449 450 451

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

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

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

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

458 459 460 461 462 463 464
  return Bio::EnsEMBL::Slice->new(-COORD_SYSTEM      => $cs,
                                  -SEQ_REGION_NAME   => $name,
                                  -SEQ_REGION_LENGTH => $length,
                                  -START             => 1,
                                  -END               => $length,
                                  -STRAND            => 1,
                                  -ADAPTOR           => $self);
Arne Stabenau's avatar
Arne Stabenau committed
465 466
}

Graham McVicker's avatar
Graham McVicker committed
467

468 469

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

Arne Stabenau's avatar
Arne Stabenau committed
471
  Arg [1]    : Bio::EnsEMBL::Slice $slice
472 473 474 475 476 477 478 479 480 481
               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
482
  Status     : Stable
Arne Stabenau's avatar
Arne Stabenau committed
483 484

=cut
Graham McVicker's avatar
Graham McVicker committed
485

486 487 488
sub get_seq_region_id {
  my $self = shift;
  my $slice = shift;
Graham McVicker's avatar
Graham McVicker committed
489

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

498 499
  if( $arr ) {
    return $arr->[0];
500 501
  }

502
  my $cs_id = $slice->coord_system->dbID();
503 504 505 506 507 508

  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
509 510 511
  $sth->bind_param(1,"$seq_region_name",SQL_VARCHAR);
  $sth->bind_param(2,$cs_id,SQL_INTEGER);
  $sth->execute();
512 513 514

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

518 519 520 521 522 523
  }

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

  #cache information for future requests
524
  $arr = [ $seq_region_id, $seq_region_name, $cs_id, $length ];
525 526 527 528

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

529
  return $seq_region_id;
Arne Stabenau's avatar
Arne Stabenau committed
530
}
531 532


Arne Stabenau's avatar
Arne Stabenau committed
533

534 535 536 537 538 539 540 541 542
=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
543
  Arg [3]    : bool $include_non_reference (optional)
544 545 546
               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
547 548
  Arg [4]    : int $include_duplicates (optional)
               If set duplicate regions will be returned.
549

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

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

  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.
564 565
               If the coordinate system with the provided name and version
               does not exist an empty list is returned.
566
               If the coordinate system name provided is 'toplevel', all
567
               non-redundant toplevel slices are returned.
568 569 570 571

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

572
  Returntype : listref of Bio::EnsEMBL::Slices
573
  Exceptions : none
574
  Caller     : general
575
  Status     : Stable
Web Admin's avatar
Web Admin committed
576 577 578

=cut

579 580 581 582 583
sub fetch_all {
  my $self = shift;
  my $cs_name = shift;
  my $cs_version = shift || '';

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

586 587 588
  #
  # verify existance of requested coord system and get its id
  #
589 590 591 592 593 594
  my $csa       = $self->db->get_CoordSystemAdaptor();
  my $orig_cs   = $csa->fetch_by_name($cs_name, $cs_version);

  return [] if(!$orig_cs);

  my $sth;
595

596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614
  my %bad_vals=();
  #
  # Get a hash of non reference seq regions
  #
  if(!$include_non_reference){
    my $sth2 = $self->prepare(
			      "SELECT sr.seq_region_id ".
			      "FROM seq_region sr, seq_region_attrib sra, attrib_type at ".
			      " WHERE at.code='non_ref'".
			      "  AND sra.seq_region_id=sr.seq_region_id ".
			      "  AND at.attrib_type_id=sra.attrib_type_id " );
    $sth2->execute();
    my ($seq_region_id);
    $sth2->bind_columns(\$seq_region_id);
    while($sth2->fetch()) {
      $bad_vals{$seq_region_id} = 1;
    }
    $sth2->finish();
  }
615
  #
616
  # Retrieve the seq_regions from the database
617
  #
618
  if($orig_cs->is_top_level()) {
619 620 621 622 623 624 625 626
   $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 " .
                      "WHERE at.code='toplevel' " .
                      "AND   at.attrib_type_id=sra.attrib_type_id " .
                      "AND   sra.seq_region_id=sr.seq_region_id");
627 628
    $sth->execute();
  } else {
629 630 631 632
     $sth =
       $self->prepare('SELECT seq_region_id, name, length, coord_system_id ' .
                      'FROM   seq_region ' .
                      'WHERE  coord_system_id =?');
633 634
     $sth->bind_param(1,$orig_cs->dbID,SQL_INTEGER);
     $sth->execute();
635
  }
636

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

640
  my $cache_count = 0;
641 642 643

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

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

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

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

659
        $cache_count++;
660
      }
661 662 663 664 665 666 667 668 669 670 671

      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){
672 673 674 675 676 677 678 679 680 681 682 683
        # 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];
            }
684
          }
685 686 687
        } else {
          # no duplicate regions
          push @out, $slice;
688
        }
689 690
      } else {
        # we want duplicates anyway so do not do any checks
691
        push @out, $slice;
692
      }
693
    }
694
  }
695
  $sth->finish();
696

697
  return \@out;
Web Admin's avatar
Web Admin committed
698 699 700
}


701

Web Admin's avatar
Web Admin committed
702 703 704 705 706 707
=head2 fetch_by_band

 Title   : fetch_by_band
 Usage   :
 Function: create a Slice representing a series of bands
 Example :
708
 Returns : Bio::EnsEMBL::Slice
Web Admin's avatar
Web Admin committed
709
 Args    : the band name
710
 Status     : Stable
Web Admin's avatar
Web Admin committed
711 712 713 714 715 716 717 718 719 720 721

=cut

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

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

722 723
  $sth->bind_param(1,"$band%",SQL_VARCHAR);
  $sth->execute();
Web Admin's avatar
Web Admin committed
724 725 726 727 728 729 730 731 732 733
  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");
}

734 735 736 737 738 739
=head2 fetch_by_chr_band

 Title   : fetch_by_chr_band
 Usage   :
 Function: create a Slice representing a series of bands
 Example :
740
 Returns : Bio::EnsEMBL::Slice
741
 Args    : the band name
742
 Status     : Stable
743 744 745 746 747 748

=cut

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

Web Admin's avatar
Web Admin committed
749
  my $chr_slice = $self->fetch_by_region('toplevel', $chr);
750 751 752
  my $seq_region_id = $self->get_seq_region_id($chr_slice);

  my $sth = $self->db->prepare
Web Admin's avatar
Web Admin committed
753
        ("select min(k.seq_region_start), max(k.seq_region_end) " .
754 755 756
         "from karyotype as k " .
         "where k.seq_region_id = ? and k.band like ?");

757 758 759
  $sth->bind_param(1,$seq_region_id,SQL_INTEGER);
  $sth->bind_param(2,"$band%",SQL_VARCHAR);
  $sth->execute();
760 761 762
  my ( $slice_start, $slice_end) = $sth->fetchrow_array;

  if(defined $slice_start) {
Web Admin's avatar
Web Admin committed
763
    return $self->fetch_by_region('toplevel',$chr,$slice_start,$slice_end);
764 765 766
  }

  throw("Band not recognised in database");
767 768 769 770
}



Web Admin's avatar
Web Admin committed
771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786
=head2 fetch_by_exon_stable_id

  Arg [1]    : string $exonid
               The stable id of the exon around which the slice is 
               desired
  Arg [2]    : (optional) int $size
               The length of the flanking regions the slice should encompass 
               on either side of the exon (0 by default)
  Example    : $slc = $sa->fetch_by_exon_stable_id('ENSE00000302930',10);
  Description: Creates a slice around the region of the specified exon. 
               If a context size is given, the slice is extended by that 
               number of basepairs on either side of the
               transcript.
  Returntype : Bio::EnsEMBL::Slice
  Exceptions : Thrown if the exon is not in the database.
  Caller     : general
787
  Status     : Stable
Web Admin's avatar
Web Admin committed
788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803

=cut

sub fetch_by_exon_stable_id{
  my ($self,$exonid,$size) = @_;

  throw('Exon argument is required.') if(!$exonid);

  my $ea = $self->db->get_ExonAdaptor;
  my $exon = $ea->fetch_by_stable_id($exonid);

  throw("Exon [$exonid] does not exist in DB.") if(!$exon);

  return $self->fetch_by_Feature($exon, $size);
}

Graham McVicker's avatar
Graham McVicker committed
804
=head2 fetch_by_transcript_stable_id
805

Graham McVicker's avatar
Graham McVicker committed
806 807 808 809 810 811 812 813 814
  Arg [1]    : string $transcriptid
               The stable id of the transcript around which the slice is 
               desired
  Arg [2]    : (optional) int $size
               The length of the flanking regions the slice should encompass 
               on either side of the transcript (0 by default)
  Example    : $slc = $sa->fetch_by_transcript_stable_id('ENST00000302930',10);
  Description: Creates a slice around the region of the specified transcript. 
               If a context size is given, the slice is extended by that 
815 816
               number of basepairs on either side of the
               transcript.
Graham McVicker's avatar
Graham McVicker committed
817
  Returntype : Bio::EnsEMBL::Slice
818
  Exceptions : Thrown if the transcript is not in the database.
Graham McVicker's avatar
Graham McVicker committed
819
  Caller     : general
820
  Status     : Stable
821 822 823

=cut

Graham McVicker's avatar
Graham McVicker committed
824
sub fetch_by_transcript_stable_id{
825 826
  my ($self,$transcriptid,$size) = @_;

827 828
  throw('Transcript argument is required.') if(!$transcriptid);

829
  my $ta = $self->db->get_TranscriptAdaptor;
830
  my $transcript = $ta->fetch_by_stable_id($transcriptid);
831

832 833 834 835
  throw("Transcript [$transcriptid] does not exist in DB.") if(!$transcript);

  return $self->fetch_by_Feature($transcript, $size);
}
836

Graham McVicker's avatar
Graham McVicker committed
837 838


Graham McVicker's avatar
Graham McVicker committed
839 840
=head2 fetch_by_transcript_id

Graham McVicker's avatar
Graham McVicker committed
841 842 843 844 845 846 847
  Arg [1]    : int $transcriptid
               The unique database identifier of the transcript around which 
               the slice is desired
  Arg [2]    : (optional) int $size
               The length of the flanking regions the slice should encompass 
               on either side of the transcript (0 by default)
  Example    : $slc = $sa->fetch_by_transcript_id(24, 1000);
848 849 850 851
  Description: Creates a slice around the region of the specified transcript.
               If a context size is given, the slice is extended by that
               number of basepairs on either side of the
               transcript.
Graham McVicker's avatar
Graham McVicker committed
852
  Returntype : Bio::EnsEMBL::Slice
853 854
  Exceptions : throw on incorrect args
               throw if transcript is not in database
Graham McVicker's avatar
Graham McVicker committed
855
  Caller     : general
856
  Status     : Stable
857 858 859

=cut

Graham McVicker's avatar
Graham McVicker committed
860
sub fetch_by_transcript_id {
861
  my ($self,$transcriptid,$size) = @_;
Graham McVicker's avatar
Graham McVicker committed
862

863
  throw('Transcript id argument is required.') if(!$transcriptid);
Graham McVicker's avatar
Graham McVicker committed
864

865 866 867 868 869 870
  my $transcript_adaptor = $self->db()->get_TranscriptAdaptor();
  my $transcript = $transcript_adaptor->fetch_by_dbID($transcriptid);

  throw("Transcript [$transcriptid] does not exist in DB.") if(!$transcript);

  return $self->fetch_by_Feature($transcript, $size);
871 872
}

873 874


Graham McVicker's avatar
Graham McVicker committed
875
=head2 fetch_by_gene_stable_id
876

Graham McVicker's avatar
Graham McVicker committed
877
  Arg [1]    : string $geneid
878
               The stable id of the gene around which the slice is
Graham McVicker's avatar
Graham McVicker committed
879 880
               desired
  Arg [2]    : (optional) int $size
Graham McVicker's avatar
Graham McVicker committed
881
               The length of the flanking regions the slice should encompass
Graham McVicker's avatar
Graham McVicker committed
882 883
               on either side of the gene (0 by default)
  Example    : $slc = $sa->fetch_by_transcript_stable_id('ENSG00000012123',10);
Graham McVicker's avatar
Graham McVicker committed
884 885 886
  Description: Creates a slice around the region of the specified gene.
               If a context size is given, the slice is extended by that
               number of basepairs on either side of the gene.
887
               The slice will be created in the genes native coordinate system.
Graham McVicker's avatar
Graham McVicker committed
888
  Returntype : Bio::EnsEMBL::Slice
889 890
  Exceptions : throw on incorrect args
               throw if transcript does not exist
Graham McVicker's avatar
Graham McVicker committed
891
  Caller     : general
892
  Status     : Stable
893 894 895

=cut

896 897
sub fetch_by_gene_stable_id {
  my ($self,$geneid,$size) = @_;