SliceAdaptor.pm 48.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 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 180 181
  $start  = 1 if(!defined($start));
  $strand = 1 if(!defined($strand));

  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 187 188 189 190 191 192 193 194 195 196 197 198
  if($coord_system_name) {
    $cs = $csa->fetch_by_name($coord_system_name,$version);

    if(!$cs) {
      throw("Unknown coordinate system:\n name='$coord_system_name' " .
            "version='$version'\n");
    }

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

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

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

    $constraint = "sr.coord_system_id = ?";

213
    $key = "$seq_region_name:".$cs->dbID();
214 215 216 217 218 219 220 221 222 223
  } 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;
    }
224
    $constraint .= "ORDER BY cs.rank ASC";
225 226
  }

227 228
  #check the cache so we only go to the db if necessary
  my $length;
229 230
  my $arr;
  if( $key ) {
231
    $arr = $self->{'sr_name_cache'}->{$key};
232
  }
233

234 235
  if( $arr ) {
    $length = $arr->[3];
236
  } else {
237 238
    my $sth = $self->prepare($sql . " WHERE sr.name = ? AND " .
                             $constraint);
Graham McVicker's avatar
Graham McVicker committed
239

240 241 242 243
    # 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.
    $sth->execute("$seq_region_name", @bind_vals);
244

245 246
    if($sth->rows() == 0) {
      $sth->finish();
Graham McVicker's avatar
Graham McVicker committed
247

248
      #do fuzzy matching, assuming that we are just missing a version on
249
      #the end of the seq_region name
Graham McVicker's avatar
Graham McVicker committed
250

251 252
      $sth = $self->prepare($sql . " WHERE sr.name LIKE ? AND " . $constraint);
      $sth->execute("$seq_region_name.%", @bind_vals);
253 254

      my $prefix_len = length($seq_region_name) + 1;
255 256
      my $high_ver = undef;
      my $high_cs = $cs;
257 258 259 260

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

261 262 263 264 265
      my ($tmp_name, $id, $tmp_length, $cs_id);
      $sth->bind_columns(\$tmp_name, \$id, \$tmp_length, \$cs_id);

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

267
        # cache values for future reference
268
        my $arr = [ $id, $tmp_name, $cs_id, $tmp_length ];
269 270
        $self->{'sr_name_cache'}->{"$tmp_name:$cs_id"} = $arr;
        $self->{'sr_id_cache'}->{"$id"} = $arr;
Graham McVicker's avatar
Graham McVicker committed
271

272
        my $tmp_ver = substr($tmp_name, $prefix_len);
273

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

277 278 279 280 281 282 283 284
        # 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
285 286
        }
      }
287 288
      $sth->finish();

289 290
      $cs = $high_cs;

291
      #return if we did not find any appropriate match:
292
      return undef if(!defined($high_ver));
293 294
    } else {

295 296
      my ($id, $cs_id);
      ($seq_region_name, $id, $length, $cs_id) = $sth->fetchrow_array();
297
      $sth->finish();
298

299 300
      # cahce to speed up for future queries
      my $arr = [ $id, $seq_region_name, $cs_id, $length ];
301 302
      $self->{'sr_name_cache'}->{"$seq_region_name:$cs_id"} = $arr;
      $self->{'sr_id_cache'}->{"$id"} = $arr;
303 304

      $cs = $csa->fetch_by_dbID( $cs_id );
305
    }
306
  }
307

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

310 311
  if($end+1 < $start) {
    throw("start [$start] must be less than or equal to end+1 [$end+1]");
312 313
  }

314
  return Bio::EnsEMBL::Slice->new(-COORD_SYSTEM      => $cs,
315 316 317 318 319 320
                                  -SEQ_REGION_NAME   => $seq_region_name,
                                  -SEQ_REGION_LENGTH => $length,
                                  -START             => $start,
                                  -END               => $end,
                                  -STRAND            => $strand,
                                  -ADAPTOR           => $self);
321 322 323 324
}



325

326
=head2 fetch_by_name
327

Graham McVicker's avatar
Graham McVicker committed
328
  Arg [1]    : string $name
329 330 331 332 333 334 335 336 337
  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
338

339 340 341 342
               Returns undef if no seq_region with the provided name exists in
               the database.

  Returntype : Bio::EnsEMBL::Slice or undef
343 344
  Exceptions : throw if incorrent arg provided
  Caller     : Pipeline
345
  Status     : Stable
346 347 348

=cut

349 350 351
sub fetch_by_name {
  my $self = shift;
  my $name = shift;
Graham McVicker's avatar
Graham McVicker committed
352

353 354
  if(!$name) {
    throw("name argument is required");
Graham McVicker's avatar
Graham McVicker committed
355
  }
356 357 358

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

359
  if(scalar(@array) < 3 || scalar(@array) > 6) {
360
    throw("Malformed slice name [$name].  Format is " .
Laura Clarke's avatar
 
Laura Clarke committed
361
        "coord_system:version:name:start:end:strand");
362 363
  }

364
  # Rearrange arguments to suit fetch_by_region
365

366 367 368 369 370 371 372 373 374
  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
375
}
Graham McVicker's avatar
Graham McVicker committed
376 377 378



379 380 381 382 383 384 385 386
=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.
387 388
               Returns undef if no such slice exists.
  Returntype : Bio::EnsEMBL::Slice or undef
Graham McVicker's avatar
Graham McVicker committed
389 390
  Exceptions : none
  Caller     : general
391
  Status     : Stable
Graham McVicker's avatar
Graham McVicker committed
392 393 394

=cut

395 396 397
sub fetch_by_seq_region_id {
  my ($self, $seq_region_id) = @_;

398
  my $arr = $self->{'sr_id_cache'}->{ $seq_region_id };
399 400
  my ($name, $length, $cs);

401 402 403 404
  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);
405 406 407 408 409 410 411
  } else {
    my $sth = $self->prepare("SELECT name, length, coord_system_id " .
                             "FROM seq_region " .
                             "WHERE seq_region_id = ?");

    $sth->execute($seq_region_id);

412
    return undef if($sth->rows() == 0);
413 414 415 416 417 418 419 420

    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
421
    my $arr = [ $seq_region_id, $name, $cs_id, $length ];
422 423 424

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

427 428 429 430 431 432 433
  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
434 435
}

Graham McVicker's avatar
Graham McVicker committed
436

437 438

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

Arne Stabenau's avatar
Arne Stabenau committed
440
  Arg [1]    : Bio::EnsEMBL::Slice $slice
441 442 443 444 445 446 447 448 449 450
               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
451
  Status     : Stable
Arne Stabenau's avatar
Arne Stabenau committed
452 453

=cut
Graham McVicker's avatar
Graham McVicker committed
454

455 456 457
sub get_seq_region_id {
  my $self = shift;
  my $slice = shift;
Graham McVicker's avatar
Graham McVicker committed
458

459 460 461
  if(!$slice || !ref($slice) || !$slice->isa('Bio::EnsEMBL::Slice')) {
    throw('Slice argument is required');
  }
462
  
463
  my $seq_region_name = $slice->seq_region_name();
464
  my $key = $seq_region_name.":".$slice->coord_system->dbID();
465
  my $arr = $self->{'sr_name_cache'}->{"$key"};
466

467 468
  if( $arr ) {
    return $arr->[0];
469 470
  }

471
  my $cs_id = $slice->coord_system->dbID();
472 473 474 475 476 477

  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
478
  $sth->execute("$seq_region_name", $cs_id );
479 480 481

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

485 486 487 488 489 490
  }

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

  #cache information for future requests
491
  $arr = [ $seq_region_id, $seq_region_name, $cs_id, $length ];
492 493 494 495

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

496
  return $seq_region_id;
Arne Stabenau's avatar
Arne Stabenau committed
497
}
498 499


Arne Stabenau's avatar
Arne Stabenau committed
500

501 502 503 504 505 506 507 508 509
=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
510
  Arg [3]    : bool $include_non_reference (optional)
511 512 513
               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
514 515
  Arg [4]    : int $include_duplicates (optional)
               If set duplicate regions will be returned.
516

517 518 519
  Example    : @chromos = @{$slice_adaptor->fetch_all('chromosome','NCBI33')};
               @contigs = @{$slice_adaptor->fetch_all('contig')};

Graham McVicker's avatar
Graham McVicker committed
520 521 522 523 524
               # 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)};
525 526 527 528 529 530

  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.
531 532
               If the coordinate system with the provided name and version
               does not exist an empty list is returned.
533
               If the coordinate system name provided is 'toplevel', all
534
               non-redundant toplevel slices are returned.
535 536 537 538

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

539
  Returntype : listref of Bio::EnsEMBL::Slices
540
  Exceptions : none
541
  Caller     : general
542
  Status     : Stable
Web Admin's avatar
Web Admin committed
543 544 545

=cut

546 547 548 549 550
sub fetch_all {
  my $self = shift;
  my $cs_name = shift;
  my $cs_version = shift || '';

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

553 554 555
  #
  # verify existance of requested coord system and get its id
  #
556 557 558 559 560 561
  my $csa       = $self->db->get_CoordSystemAdaptor();
  my $orig_cs   = $csa->fetch_by_name($cs_name, $cs_version);

  return [] if(!$orig_cs);

  my $sth;
562

563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581
  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();
  }
582
  #
583
  # Retrieve the seq_regions from the database
584
  #
585
  if($orig_cs->is_top_level()) {
586 587 588 589 590 591 592 593
   $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");
594 595
    $sth->execute();
  } else {
596 597 598 599 600
     $sth =
       $self->prepare('SELECT seq_region_id, name, length, coord_system_id ' .
                      'FROM   seq_region ' .
                      'WHERE  coord_system_id =?');
     $sth->execute($orig_cs->dbID);
601
  }
602

603 604
  my ($seq_region_id, $name, $length, $cs_id);
  $sth->bind_columns(\$seq_region_id, \$name, \$length, \$cs_id);
605

606
  my $cache_count = 0;
607 608 609

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

613
      if(!$cs) {
614
        throw("seq_region $name references non-existent coord_system $cs_id.");
615
      }
616

617 618
      #cache values for future reference, but stop adding to the cache once we
      #we know we have filled it up
619
      if($cache_count < $Bio::EnsEMBL::Utils::SeqRegionCache::SEQ_REGION_CACHE_SIZE) {
620
        my $arr = [ $seq_region_id, $name, $cs_id, $length ];
621 622 623

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

625
        $cache_count++;
626
      }
627 628 629 630 631 632 633 634 635 636 637

      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){
638 639 640 641 642 643 644 645 646 647 648 649
        # 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];
            }
650
          }
651 652 653
        } else {
          # no duplicate regions
          push @out, $slice;
654
        }
655 656
      } else {
        # we want duplicates anyway so do not do any checks
657
        push @out, $slice;
658
      }
659
    }
660
  }
661
  $sth->finish();
662

663
  return \@out;
Web Admin's avatar
Web Admin committed
664 665 666
}


667

Web Admin's avatar
Web Admin committed
668 669 670 671 672 673
=head2 fetch_by_band

 Title   : fetch_by_band
 Usage   :
 Function: create a Slice representing a series of bands
 Example :
674
 Returns : Bio::EnsEMBL::Slice
Web Admin's avatar
Web Admin committed
675
 Args    : the band name
676
 Status     : Stable
Web Admin's avatar
Web Admin committed
677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698

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

  $sth->execute( "$band%" );
  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");
}

699 700 701 702 703 704
=head2 fetch_by_chr_band

 Title   : fetch_by_chr_band
 Usage   :
 Function: create a Slice representing a series of bands
 Example :
705
 Returns : Bio::EnsEMBL::Slice
706
 Args    : the band name
707
 Status     : Stable
708 709 710 711 712 713

=cut

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

Web Admin's avatar
Web Admin committed
714
  my $chr_slice = $self->fetch_by_region('toplevel', $chr);
715 716 717
  my $seq_region_id = $self->get_seq_region_id($chr_slice);

  my $sth = $self->db->prepare
Web Admin's avatar
Web Admin committed
718
        ("select min(k.seq_region_start), max(k.seq_region_end) " .
719 720 721 722 723 724 725
         "from karyotype as k " .
         "where k.seq_region_id = ? and k.band like ?");

  $sth->execute( $seq_region_id, "$band%" );
  my ( $slice_start, $slice_end) = $sth->fetchrow_array;

  if(defined $slice_start) {
Web Admin's avatar
Web Admin committed
726
    return $self->fetch_by_region('toplevel',$chr,$slice_start,$slice_end);
727 728 729
  }

  throw("Band not recognised in database");
730 731 732 733
}



Web Admin's avatar
Web Admin committed
734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749
=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
750
  Status     : Stable
Web Admin's avatar
Web Admin committed
751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766

=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
767
=head2 fetch_by_transcript_stable_id
768

Graham McVicker's avatar
Graham McVicker committed
769 770 771 772 773 774 775 776 777
  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 
778 779
               number of basepairs on either side of the
               transcript.
Graham McVicker's avatar
Graham McVicker committed
780
  Returntype : Bio::EnsEMBL::Slice
781
  Exceptions : Thrown if the transcript is not in the database.
Graham McVicker's avatar
Graham McVicker committed
782
  Caller     : general
783
  Status     : Stable
784 785 786

=cut

Graham McVicker's avatar
Graham McVicker committed
787
sub fetch_by_transcript_stable_id{
788 789
  my ($self,$transcriptid,$size) = @_;

790 791
  throw('Transcript argument is required.') if(!$transcriptid);

792
  my $ta = $self->db->get_TranscriptAdaptor;
793
  my $transcript = $ta->fetch_by_stable_id($transcriptid);
794

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

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

Graham McVicker's avatar
Graham McVicker committed
800 801


Graham McVicker's avatar
Graham McVicker committed
802 803
=head2 fetch_by_transcript_id

Graham McVicker's avatar
Graham McVicker committed
804 805 806 807 808 809 810
  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);
811 812 813 814
  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
815
  Returntype : Bio::EnsEMBL::Slice
816 817
  Exceptions : throw on incorrect args
               throw if transcript is not in database
Graham McVicker's avatar
Graham McVicker committed
818
  Caller     : general
819
  Status     : Stable
820 821 822

=cut

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

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

828 829 830 831 832 833
  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);
834 835
}

836 837


Graham McVicker's avatar
Graham McVicker committed
838
=head2 fetch_by_gene_stable_id
839

Graham McVicker's avatar
Graham McVicker committed
840
  Arg [1]    : string $geneid
841
               The stable id of the gene around which the slice is
Graham McVicker's avatar
Graham McVicker committed
842 843
               desired
  Arg [2]    : (optional) int $size
Graham McVicker's avatar
Graham McVicker committed
844
               The length of the flanking regions the slice should encompass
Graham McVicker's avatar
Graham McVicker committed
845 846
               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
847 848 849
  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.
850
               The slice will be created in the genes native coordinate system.
Graham McVicker's avatar
Graham McVicker committed
851
  Returntype : Bio::EnsEMBL::Slice
852 853
  Exceptions : throw on incorrect args
               throw if transcript does not exist
Graham McVicker's avatar
Graham McVicker committed
854
  Caller     : general
855
  Status     : Stable
856 857 858

=cut

859 860
sub fetch_by_gene_stable_id {
  my ($self,$geneid,$size) = @_;
861

862
  throw('Gene argument is required.') if(!$geneid);
863

864 865
  my $gene_adaptor = $self->db->get_GeneAdaptor();
  my $gene = $gene_adaptor->fetch_by_stable_id($geneid);
866

867
  throw("Gene [$geneid] does not exist in DB.") if(!$gene);
868

869
  return $self->fetch_by_Feature($gene, $size);
870 871 872
}


873

874
=head2 fetch_by_Feature
Graham McVicker's avatar
Graham McVicker committed
875

876 877 878 879
  Arg [1]    : Bio::EnsEMBL::Feature $feat
               The feature to fetch the slice around
  Arg [2]    : int size (optional)
               The desired number of flanking basepairs around the feature.
880 881
               The size may also be provided as a percentage of the feature 
               size such as 200% or 80.5%.
882 883 884 885 886 887
  Example    : $slice = $slice_adaptor->fetch_by_Feature($feat, 100);
  Description: Retrieves a slice around a specific feature.  All this really
               does is return a resized version of the slice that the feature
               is already on. Note that slices returned from this method
               are always on the forward strand of the seq_region regardless of
               the strandedness of the feature passed in.
Graham McVicker's avatar
Graham McVicker committed
888
  Returntype : Bio::EnsEMBL::Slice
889 890 891 892
  Exceptions : throw if the feature does not have an attached slice
               throw if feature argument is not provided
  Caller     : fetch_by_gene_stable_id, fetch_by_transcript_stable_id,
               fetch_by_gene_id, fetch_by_transcript_id
893
  Status     : Stable
Graham McVicker's avatar
Graham McVicker committed
894 895 896

=cut

897 898
sub fetch_by_Feature{
  my ($self, $feature, $size) = @_;
Graham McVicker's avatar
Graham McVicker committed
899

900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931
  $size ||= 0;

  if(!ref($feature) || !$feature->isa('Bio::EnsEMBL::Feature')) {
    throw('Feature argument expected.');
  }

  my $slice = $feature->slice();
  if(!$slice || !$slice->isa('Bio::EnsEMBL::Slice')) {
    throw('Feature must be attached to a valid slice.');
  }

  my $fstart = $feature->start();
  my $fend   = $feature->end();
  if(!defined($fstart) || !defined($fend)) {
    throw('Feature must have defined start and end.');
  }

  #convert the feature slice coordinates to seq_region coordinates
  my $slice_start  = $slice->start();
  my $slice_end    = $slice->end();
  my $slice_strand = $slice->strand();
  if($slice_start != 1 || $slice_strand != 1) {
    if($slice_strand == 1) {
      $fstart = $fstart + $slice_start - 1;
      $fend   = $fend   + $slice_start - 1;
    } else {
      my $tmp_start = $fstart;
      $fstart = $slice_end - $fend      + 1;
      $fend   = $slice_end - $tmp_start + 1;
    }
  }

932 933 934 935
  ## Size may be stored as a %age of the length of the feature
  ## Size = 100% gives no context
  ## Size = 200% gives context - 50% the size of the feature either side of 
  ## feature
Web Admin's avatar
Web Admin committed
936 937 938

  $size = int( ($1-100)/200 * ($fend-$fstart+1) ) if( $size =~/([\d+\.]+)%/ );

939
  #return a new slice covering the region of the feature
940
  my $S = Bio::EnsEMBL::Slice->new
941 942 943 944 945 946 947
    (-seq_region_name   => $slice->seq_region_name,
     -seq_region_length => $slice->seq_region_length,
     -coord_system      => $slice->coord_system,
     -start             => $fstart - $size,
     -end               => $fend + $size,
     -strand            => 1,
     -adaptor           => $self);
948 949
  $S->{'_raw_feature_strand'}  = $feature->strand * $slice_strand if $feature->can('strand');
  return $S;
950 951