SliceAdaptor.pm 54.5 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 stack_trace_dump);
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
  # use a cache which is shared and also used by the assembly
  # mapper adaptor
110

111 112 113 114
  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'};
115 116 117 118 119 120 121

  return $self;
}


=head2 fetch_by_region

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

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

=cut

171 172 173 174

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

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

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

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

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

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

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

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

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

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

218
	      $constraint = "sr.coord_system_id = ?";
219

220
    $key = "$seq_region_name:".$cs->dbID();
221 222 223 224 225 226 227 228 229 230
  } 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;
    }
231
    $constraint .= "ORDER BY cs.rank ASC";
232 233
  }

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

241
  if ($arr) {
242
    $length = $arr->[3];
243
  } else {
244
    my $sth = $self->prepare($sql . " WHERE sr.name = ? AND " . $constraint);
245 246 247
    # 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.
248
    $sth->bind_param(1, "$seq_region_name", SQL_VARCHAR);
249
    if ($cs){
250
	$sth->bind_param(2, $cs->dbID, SQL_INTEGER);
251 252
    }
    else{
253
	$sth->bind_param(2, $version, SQL_VARCHAR)  if ($version);	
254
    }
255 256

    $sth->execute();
257

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

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

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

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

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

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

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

Felix Kokocinski's avatar
Felix Kokocinski committed
286
      my $i = 0;
287 288

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

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

296
        my $tmp_ver = substr($tmp_name, $prefix_len);
297

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

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

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

315 316 317 318 319
      # 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");
      }

320 321
      $cs = $high_cs;

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

325 326
    } else {

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

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

      $cs = $csa->fetch_by_dbID( $cs_id );
337
    }
338
  }
339

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

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

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



357

358
=head2 fetch_by_name
359

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

371 372 373 374
               Returns undef if no seq_region with the provided name exists in
               the database.

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

=cut

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

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

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

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

396
  # Rearrange arguments to suit fetch_by_region
397

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



411 412 413 414 415
=head2 fetch_by_seq_region_id

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

=cut

430
sub fetch_by_seq_region_id {
431
  my ($self, $seq_region_id,$start,$end,$strand) = @_;
432

433
  my $arr = $self->{'sr_id_cache'}->{ $seq_region_id };
434 435
  my ($name, $length, $cs);

436 437 438 439
  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);
440 441 442 443 444
  } else {
    my $sth = $self->prepare("SELECT name, length, coord_system_id " .
                             "FROM seq_region " .
                             "WHERE seq_region_id = ?");

445 446
    $sth->bind_param(1,$seq_region_id,SQL_INTEGER);
    $sth->execute();
447

448
    return undef if($sth->rows() == 0);
449 450 451 452 453 454 455 456

    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
457
    my $arr = [ $seq_region_id, $name, $cs_id, $length ];
458 459 460

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

463 464 465
  return Bio::EnsEMBL::Slice->new(-COORD_SYSTEM      => $cs,
                                  -SEQ_REGION_NAME   => $name,
                                  -SEQ_REGION_LENGTH => $length,
466 467 468
                                  -START             => $start || 1,
                                  -END               => $end || $length,
                                  -STRAND            => $strand || 1,
469
                                  -ADAPTOR           => $self);
Arne Stabenau's avatar
Arne Stabenau committed
470 471
}

Graham McVicker's avatar
Graham McVicker committed
472

473 474

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

Arne Stabenau's avatar
Arne Stabenau committed
476
  Arg [1]    : Bio::EnsEMBL::Slice $slice
477 478 479 480 481 482 483 484 485 486
               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
487
  Status     : Stable
Arne Stabenau's avatar
Arne Stabenau committed
488 489

=cut
Graham McVicker's avatar
Graham McVicker committed
490

491 492 493
sub get_seq_region_id {
  my $self = shift;
  my $slice = shift;
Graham McVicker's avatar
Graham McVicker committed
494

495 496 497
  if(!$slice || !ref($slice) || !$slice->isa('Bio::EnsEMBL::Slice')) {
    throw('Slice argument is required');
  }
498
  
499
  my $seq_region_name = $slice->seq_region_name();
500
  my $key = $seq_region_name.":".$slice->coord_system->dbID();
501
  my $arr = $self->{'sr_name_cache'}->{"$key"};
502

503 504
  if( $arr ) {
    return $arr->[0];
505 506
  }

507
  my $cs_id = $slice->coord_system->dbID();
508 509 510 511 512 513

  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
514 515 516
  $sth->bind_param(1,"$seq_region_name",SQL_VARCHAR);
  $sth->bind_param(2,$cs_id,SQL_INTEGER);
  $sth->execute();
517 518 519

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

523 524 525 526 527 528
  }

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

  #cache information for future requests
529
  $arr = [ $seq_region_id, $seq_region_name, $cs_id, $length ];
530 531 532 533

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

534
  return $seq_region_id;
Arne Stabenau's avatar
Arne Stabenau committed
535
}
536 537


Arne Stabenau's avatar
Arne Stabenau committed
538

539 540 541 542 543 544 545 546 547
=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
548
  Arg [3]    : bool $include_non_reference (optional)
549 550 551
               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
552 553
  Arg [4]    : int $include_duplicates (optional)
               If set duplicate regions will be returned.
554
               
555
               NOTE: if you do not use this option and you have a PAR
556 557 558
               (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.
559

560 561 562
  Example    : @chromos = @{$slice_adaptor->fetch_all('chromosome','NCBI33')};
               @contigs = @{$slice_adaptor->fetch_all('contig')};

Graham McVicker's avatar
Graham McVicker committed
563 564 565 566 567
               # 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)};
568 569 570 571 572 573

  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.
574 575
               If the coordinate system with the provided name and version
               does not exist an empty list is returned.
576
               If the coordinate system name provided is 'toplevel', all
577
               non-redundant toplevel slices are returned.
578 579 580 581

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

582
  Returntype : listref of Bio::EnsEMBL::Slices
583
  Exceptions : none
584
  Caller     : general
585
  Status     : Stable
Web Admin's avatar
Web Admin committed
586 587 588

=cut

589 590 591 592 593
sub fetch_all {
  my $self = shift;
  my $cs_name = shift;
  my $cs_version = shift || '';

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

596 597 598
  #
  # verify existance of requested coord system and get its id
  #
599 600 601 602 603 604
  my $csa       = $self->db->get_CoordSystemAdaptor();
  my $orig_cs   = $csa->fetch_by_name($cs_name, $cs_version);

  return [] if(!$orig_cs);

  my $sth;
605

606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624
  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();
  }
625
  #
626
  # Retrieve the seq_regions from the database
627
  #
628
  if($orig_cs->is_top_level()) {
629 630 631 632 633 634 635 636
   $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");
637 638
    $sth->execute();
  } else {
639 640 641 642
     $sth =
       $self->prepare('SELECT seq_region_id, name, length, coord_system_id ' .
                      'FROM   seq_region ' .
                      'WHERE  coord_system_id =?');
643 644
     $sth->bind_param(1,$orig_cs->dbID,SQL_INTEGER);
     $sth->execute();
645
  }
646

647 648
  my ($seq_region_id, $name, $length, $cs_id);
  $sth->bind_columns(\$seq_region_id, \$name, \$length, \$cs_id);
649

650
  my $cache_count = 0;
651 652 653

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

657
      if(!$cs) {
658
        throw("seq_region $name references non-existent coord_system $cs_id.");
659
      }
660

661 662
      #cache values for future reference, but stop adding to the cache once we
      #we know we have filled it up
663
      if($cache_count < $Bio::EnsEMBL::Utils::SeqRegionCache::SEQ_REGION_CACHE_SIZE) {
664
        my $arr = [ $seq_region_id, $name, $cs_id, $length ];
665 666 667

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

669
        $cache_count++;
670
      }
671 672 673 674 675 676 677 678 679 680 681

      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){
682 683 684 685 686 687 688 689 690 691 692 693
        # 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];
            }
694
          }
695 696 697
        } else {
          # no duplicate regions
          push @out, $slice;
698
        }
699 700
      } else {
        # we want duplicates anyway so do not do any checks
701
        push @out, $slice;
702
      }
703
    }
704
  }
705
  $sth->finish();
706

707
  return \@out;
Web Admin's avatar
Web Admin committed
708 709
}

710 711 712 713 714 715 716 717 718 719 720 721 722
=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

sub is_toplevel{
  my $self = shift;
  my $id = shift;
Web Admin's avatar
Web Admin committed
723

724 725 726 727 728 729 730
  my $sth = $self->prepare("SELECT at.code from seq_region_attrib sra, attrib_type at WHERE sra.seq_region_id = $id AND  at.attrib_type_id = sra.attrib_type_id AND at.code = 'toplevel'");
  
  $sth->execute();
  
  my $code = undef;
  $sth->bind_columns(\$code);
  while($sth->fetch){
Ian Longden's avatar
Ian Longden committed
731
    $sth->finish;
732 733
    return 1;
  }
Ian Longden's avatar
Ian Longden committed
734
  $sth->finish;
735 736
  return 0;
}
737

Web Admin's avatar
Web Admin committed
738 739 740 741 742 743
=head2 fetch_by_band

 Title   : fetch_by_band
 Usage   :
 Function: create a Slice representing a series of bands
 Example :
744
 Returns : Bio::EnsEMBL::Slice
Web Admin's avatar
Web Admin committed
745
 Args    : the band name
746
 Status     : Stable
Web Admin's avatar
Web Admin committed
747 748 749 750 751 752 753 754 755 756 757

=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");

758 759
  $sth->bind_param(1,"$band%",SQL_VARCHAR);
  $sth->execute();
Web Admin's avatar
Web Admin committed
760 761 762 763 764 765 766 767 768 769
  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");
}

770 771 772 773 774 775
=head2 fetch_by_chr_band

 Title   : fetch_by_chr_band
 Usage   :
 Function: create a Slice representing a series of bands
 Example :
776
 Returns : Bio::EnsEMBL::Slice
777
 Args    : the band name
778
 Status     : Stable
779 780 781 782 783 784

=cut

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

Web Admin's avatar
Web Admin committed
785
  my $chr_slice = $self->fetch_by_region('toplevel', $chr);
786 787 788
  my $seq_region_id = $self->get_seq_region_id($chr_slice);

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

793 794 795
  $sth->bind_param(1,$seq_region_id,SQL_INTEGER);
  $sth->bind_param(2,"$band%",SQL_VARCHAR);
  $sth->execute();
796 797 798
  my ( $slice_start, $slice_end) = $sth->fetchrow_array;

  if(defined $slice_start) {
Web Admin's avatar
Web Admin committed
799
    return $self->fetch_by_region('toplevel',$chr,$slice_start,$slice_end);
800 801 802
  }

  throw("Band not recognised in database");
803 804 805 806
}



Web Admin's avatar
Web Admin committed
807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822
=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
823
  Status     : Stable
Web Admin's avatar
Web Admin committed
824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839

=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
840
=head2 fetch_by_transcript_stable_id
841

Graham McVicker's avatar
Graham McVicker committed
842 843 844 845 846 847 848 849 850
  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 
851 852
               number of basepairs on either side of the
               transcript.
Graham McVicker's avatar
Graham McVicker committed
853
  Returntype : Bio::EnsEMBL::Slice
854
  Exceptions : Thrown if the transcript is not in the 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_stable_id{
861 862
  my ($self,$transcriptid,$size) = @_;

863 864
  throw('Transcript argument is required.') if(!$transcriptid);

865
  my $ta = $self->db->get_TranscriptAdaptor;
866
  my $transcript = $ta->fetch_by_stable_id($transcriptid);
867

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

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

Graham McVicker's avatar
Graham McVicker committed
873 874


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

Graham McVicker's avatar
Graham McVicker committed
877 878 879 880 881 882 883
  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);
884 885 886 887
  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
888
  Returntype : Bio::EnsEMBL::Slice
889 890
  Exceptions : throw on incorrect args
               throw if transcript is not in database
Graham McVicker's avatar
Graham McVicker committed
891
  Caller     : general
892
  Status     : Stable
893 894 895

=cut

Graham McVicker's avatar
Graham McVicker committed
896
sub fetch_by_transcript_id {
897