SliceAdaptor.pm 55.1 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
  Arg [7]    : boolean $no_fuzz (optional, default = undef (false))
               If true (non-zero), do not use "fuzzy matching" (see below).
140 141 142
  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
143 144 145 146 147 148 149 150 151 152 153
               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.
154 155 156 157 158 159 160

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

161 162 163
               The fuzzy matching can be turned off by setting the
               $no_fuzz argument to a true value.

164 165 166 167 168 169 170
               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
171
  Caller     : general
172
  Status     : Stable
173 174 175

=cut

176 177 178 179

#
# ARNE: This subroutine needs simplification!! 
#
180
sub fetch_by_region {
181 182 183
  my ( $self, $coord_system_name, $seq_region_name, $start, $end,
       $strand, $version, $no_fuzz )
    = @_;
184

185 186
  $start  = 1 if (!defined($start));
  $strand = 1 if (!defined($strand));
187

188
  throw('seq_region_name argument is required') if (!defined($seq_region_name));
189

190 191
  my $cs;
  my $csa = $self->db->get_CoordSystemAdaptor();
192

193
  if ($coord_system_name) {
194 195
    $cs = $csa->fetch_by_name($coord_system_name,$version);

196
    ## REMOVE THESE THREE LINES WHEN STICKLEBACK DB IS FIXED!
197 198 199 200 201 202 203
    ## Anne/ap5 (2007-10-09):
    # The problem was that the stickleback genebuild called the
    # chromosomes 'groups', which meant they weren't being picked out by
    # the karyotype drawing code.  Apparently they are usually called
    # 'groups' in the stickleback community, even though they really are
    # chromosomes!

204 205 206 207
    if( !$cs && $coord_system_name eq 'chromosome' ) {
      $cs = $csa->fetch_by_name('group',$version);
    }

208
    if (!$cs) {
209 210 211 212
      throw("Unknown coordinate system:\n name='$coord_system_name' " .
            "version='$version'\n");
    }
    # fetching by toplevel is same as fetching w/o name or version
213
    if ($cs->is_top_level()) {
214 215 216
      $cs = undef;
      $version = undef;
    }
217

218
  }
219 220 221 222 223 224

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

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

231
	      $constraint = "sr.coord_system_id = ?";
232

233
    $key = "$seq_region_name:".$cs->dbID();
234 235 236 237 238 239 240 241 242 243
  } 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;
    }
244
    $constraint .= "ORDER BY cs.rank ASC";
245 246
  }

247
  # check the cache so we only go to the db if necessary
248
  my $length;
249
  my $arr;
250
  if ($key) {
251
    $arr = $self->{'sr_name_cache'}->{$key};
252
  }
253

254
  if ($arr) {
255
    $length = $arr->[3];
256
  } else {
257
    my $sth = $self->prepare($sql . " WHERE sr.name = ? AND " . $constraint);
258
    $sth->bind_param(1, $seq_region_name, SQL_VARCHAR);
259
    if ($cs){
260
	$sth->bind_param(2, $cs->dbID, SQL_INTEGER);
261 262
    }
    else{
263
	$sth->bind_param(2, $version, SQL_VARCHAR)  if ($version);	
264
    }
265 266

    $sth->execute();
267

268
    if ($sth->rows() == 0) {
269
      $sth->finish();
Graham McVicker's avatar
Graham McVicker committed
270

271 272
      if ($no_fuzz) { return undef }

273 274
      # do fuzzy matching, assuming that we are just missing a version on
      # the end of the seq_region name
275
      
276
      $sth = $self->prepare($sql . " WHERE sr.name LIKE ? AND " . $constraint);
277

278
      $sth->bind_param(1, "$seq_region_name.%", SQL_VARCHAR);
279 280

      if ($cs){
281
	  $sth->bind_param(2, $cs->dbID, SQL_INTEGER);
282 283
      }
      else{
284
	  $sth->bind_param(2, $version, SQL_VARCHAR) if ($version);	
285
      }
286
      $sth->execute();
287 288

      my $prefix_len = length($seq_region_name) + 1;
289 290
      my $high_ver = undef;
      my $high_cs = $cs;
291 292 293 294

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

295 296 297
      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
298
      my $i = 0;
299 300

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

303
        # cache values for future reference
304
        my $arr = [ $id, $tmp_name, $cs_id, $tmp_length ];
305 306
        $self->{'sr_name_cache'}->{"$tmp_name:$cs_id"} = $arr;
        $self->{'sr_id_cache'}->{"$id"} = $arr;
Graham McVicker's avatar
Graham McVicker committed
307

308
        my $tmp_ver = substr($tmp_name, $prefix_len);
309

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

313 314 315 316 317 318 319 320
        # 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
321
        }
322 323

        $i++;
Graham McVicker's avatar
Graham McVicker committed
324
      }
325 326
      $sth->finish();

327 328 329 330 331
      # 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");
      }

332 333
      $cs = $high_cs;

334
      #return if we did not find any appropriate match:
335 336
      return undef if (!defined($high_ver));

337 338
    } else {

339 340
      my ($id, $cs_id);
      ($seq_region_name, $id, $length, $cs_id) = $sth->fetchrow_array();
341
      $sth->finish();
342

343 344
      # cahce to speed up for future queries
      my $arr = [ $id, $seq_region_name, $cs_id, $length ];
345 346
      $self->{'sr_name_cache'}->{"$seq_region_name:$cs_id"} = $arr;
      $self->{'sr_id_cache'}->{"$id"} = $arr;
347
      $cs = $csa->fetch_by_dbID( $cs_id );
348
    }
349
  }
350

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

353 354
  if($end+1 < $start) {
    throw("start [$start] must be less than or equal to end+1 [$end+1]");
355 356
  }

357
  return Bio::EnsEMBL::Slice->new(-COORD_SYSTEM      => $cs,
358 359 360 361 362 363
                                  -SEQ_REGION_NAME   => $seq_region_name,
                                  -SEQ_REGION_LENGTH => $length,
                                  -START             => $start,
                                  -END               => $end,
                                  -STRAND            => $strand,
                                  -ADAPTOR           => $self);
364 365 366 367
}



368

369
=head2 fetch_by_name
370

Graham McVicker's avatar
Graham McVicker committed
371
  Arg [1]    : string $name
372 373 374 375 376 377 378 379 380
  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
381

382 383 384 385
               Returns undef if no seq_region with the provided name exists in
               the database.

  Returntype : Bio::EnsEMBL::Slice or undef
386 387
  Exceptions : throw if incorrent arg provided
  Caller     : Pipeline
388
  Status     : Stable
389 390 391

=cut

392 393 394
sub fetch_by_name {
  my $self = shift;
  my $name = shift;
Graham McVicker's avatar
Graham McVicker committed
395

396 397
  if(!$name) {
    throw("name argument is required");
Graham McVicker's avatar
Graham McVicker committed
398
  }
399 400 401

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

402
  if(scalar(@array) < 3 || scalar(@array) > 6) {
403
    throw("Malformed slice name [$name].  Format is " .
Laura Clarke's avatar
 
Laura Clarke committed
404
        "coord_system:version:name:start:end:strand");
405 406
  }

407
  # Rearrange arguments to suit fetch_by_region
408

409 410 411 412 413 414 415 416 417
  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
418
}
Graham McVicker's avatar
Graham McVicker committed
419 420 421



422 423 424 425 426
=head2 fetch_by_seq_region_id

  Arg [1]    : string $seq_region_id
               The internal identifier of the seq_region to create this slice
               on
427 428 429
  Arg [2]    : optional start
  Arg [3]    : optional end
  Arg [4]    : optional strand
430 431 432
  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.
433 434
               Returns undef if no such slice exists.
  Returntype : Bio::EnsEMBL::Slice or undef
Graham McVicker's avatar
Graham McVicker committed
435 436
  Exceptions : none
  Caller     : general
437
  Status     : Stable
Graham McVicker's avatar
Graham McVicker committed
438 439 440

=cut

441
sub fetch_by_seq_region_id {
442
  my ($self, $seq_region_id,$start,$end,$strand) = @_;
443

444
  my $arr = $self->{'sr_id_cache'}->{ $seq_region_id };
445 446
  my ($name, $length, $cs);

447
  if( $arr &&  defined($arr->[2])) {
448 449 450
    my $cs_id;
    ($name, $cs_id, $length ) = ( $arr->[1], $arr->[2], $arr->[3] );
    $cs = $self->db->get_CoordSystemAdaptor->fetch_by_dbID($cs_id);
451 452 453 454 455
  } else {
    my $sth = $self->prepare("SELECT name, length, coord_system_id " .
                             "FROM seq_region " .
                             "WHERE seq_region_id = ?");

456 457
    $sth->bind_param(1,$seq_region_id,SQL_INTEGER);
    $sth->execute();
458

459
    return undef if($sth->rows() == 0);
460 461 462 463 464 465 466 467

    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
468
    my $arr = [ $seq_region_id, $name, $cs_id, $length ];
469 470 471

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

474 475 476
  return Bio::EnsEMBL::Slice->new(-COORD_SYSTEM      => $cs,
                                  -SEQ_REGION_NAME   => $name,
                                  -SEQ_REGION_LENGTH => $length,
477 478 479
                                  -START             => $start || 1,
                                  -END               => $end || $length,
                                  -STRAND            => $strand || 1,
480
                                  -ADAPTOR           => $self);
Arne Stabenau's avatar
Arne Stabenau committed
481 482
}

Graham McVicker's avatar
Graham McVicker committed
483

484 485

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

Arne Stabenau's avatar
Arne Stabenau committed
487
  Arg [1]    : Bio::EnsEMBL::Slice $slice
488 489 490 491 492 493 494 495 496 497
               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
498
  Status     : Stable
Arne Stabenau's avatar
Arne Stabenau committed
499 500

=cut
Graham McVicker's avatar
Graham McVicker committed
501

502 503 504
sub get_seq_region_id {
  my $self = shift;
  my $slice = shift;
Graham McVicker's avatar
Graham McVicker committed
505

506 507 508
  if(!$slice || !ref($slice) || !$slice->isa('Bio::EnsEMBL::Slice')) {
    throw('Slice argument is required');
  }
509
  
510
  my $seq_region_name = $slice->seq_region_name();
511
  my $key = $seq_region_name.":".$slice->coord_system->dbID();
512
  my $arr = $self->{'sr_name_cache'}->{"$key"};
513

514 515
  if( $arr ) {
    return $arr->[0];
516 517
  }

518
  my $cs_id = $slice->coord_system->dbID();
519 520 521 522 523 524

  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
525 526 527
  $sth->bind_param(1,"$seq_region_name",SQL_VARCHAR);
  $sth->bind_param(2,$cs_id,SQL_INTEGER);
  $sth->execute();
528 529 530

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

534 535 536 537 538 539
  }

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

  #cache information for future requests
540
  $arr = [ $seq_region_id, $seq_region_name, $cs_id, $length ];
541 542 543 544

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

545
  return $seq_region_id;
Arne Stabenau's avatar
Arne Stabenau committed
546
}
547 548


Arne Stabenau's avatar
Arne Stabenau committed
549

550 551 552 553 554 555 556 557 558
=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
559
  Arg [3]    : bool $include_non_reference (optional)
560 561 562
               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
563 564
  Arg [4]    : int $include_duplicates (optional)
               If set duplicate regions will be returned.
565
               
566
               NOTE: if you do not use this option and you have a PAR
567 568 569
               (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.
570

571 572 573
  Example    : @chromos = @{$slice_adaptor->fetch_all('chromosome','NCBI33')};
               @contigs = @{$slice_adaptor->fetch_all('contig')};

Graham McVicker's avatar
Graham McVicker committed
574 575 576 577 578
               # 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)};
579 580 581 582 583 584

  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.
585 586
               If the coordinate system with the provided name and version
               does not exist an empty list is returned.
587
               If the coordinate system name provided is 'toplevel', all
588 589
               non-redundant toplevel slices are returned (note that any
               coord_system_version argument is ignored in that case).
590 591 592 593

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

594
  Returntype : listref of Bio::EnsEMBL::Slices
595
  Exceptions : none
596
  Caller     : general
597
  Status     : Stable
Web Admin's avatar
Web Admin committed
598 599 600

=cut

601 602 603 604 605
sub fetch_all {
  my $self = shift;
  my $cs_name = shift;
  my $cs_version = shift || '';

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

608 609 610
  #
  # verify existance of requested coord system and get its id
  #
611 612 613 614 615 616
  my $csa       = $self->db->get_CoordSystemAdaptor();
  my $orig_cs   = $csa->fetch_by_name($cs_name, $cs_version);

  return [] if(!$orig_cs);

  my $sth;
617

618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636
  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();
  }
637
  #
638
  # Retrieve the seq_regions from the database
639
  #
640
  if($orig_cs->is_top_level()) {
641 642 643 644 645 646 647 648
   $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");
649 650
    $sth->execute();
  } else {
651 652 653 654
     $sth =
       $self->prepare('SELECT seq_region_id, name, length, coord_system_id ' .
                      'FROM   seq_region ' .
                      'WHERE  coord_system_id =?');
655 656
     $sth->bind_param(1,$orig_cs->dbID,SQL_INTEGER);
     $sth->execute();
657
  }
658

659 660
  my ($seq_region_id, $name, $length, $cs_id);
  $sth->bind_columns(\$seq_region_id, \$name, \$length, \$cs_id);
661

662
  my $cache_count = 0;
663 664 665

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

669
      if(!$cs) {
670
        throw("seq_region $name references non-existent coord_system $cs_id.");
671
      }
672

673 674
      #cache values for future reference, but stop adding to the cache once we
      #we know we have filled it up
675
      if($cache_count < $Bio::EnsEMBL::Utils::SeqRegionCache::SEQ_REGION_CACHE_SIZE) {
676
        my $arr = [ $seq_region_id, $name, $cs_id, $length ];
677 678 679

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

681
        $cache_count++;
682
      }
683 684 685 686 687 688 689 690 691 692 693

      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){
694 695 696 697 698 699 700 701 702 703 704 705
        # 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];
            }
706
          }
707 708 709
        } else {
          # no duplicate regions
          push @out, $slice;
710
        }
711 712
      } else {
        # we want duplicates anyway so do not do any checks
713
        push @out, $slice;
714
      }
715
    }
716
  }
717
  $sth->finish();
718

719
  return \@out;
Web Admin's avatar
Web Admin committed
720 721
}

722 723 724 725 726 727 728 729 730 731 732 733 734
=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
735

736 737 738 739 740 741 742
  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
743
    $sth->finish;
744 745
    return 1;
  }
Ian Longden's avatar
Ian Longden committed
746
  $sth->finish;
747 748
  return 0;
}
749

Web Admin's avatar
Web Admin committed
750 751 752 753
=head2 fetch_by_band

 Title   : fetch_by_band
 Usage   :
754
 Function: Does not work please use fetch_by_chr_band
Web Admin's avatar
Web Admin committed
755
 Example :
756
 Returns : Bio::EnsEMBL::Slice
Web Admin's avatar
Web Admin committed
757
 Args    : the band name
758
 Status     : AT RISK
Web Admin's avatar
Web Admin committed
759 760 761 762 763 764

=cut

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

765
  my $sth = $self->dbc->prepare
Web Admin's avatar
Web Admin committed
766 767 768 769
        ("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");

770 771
  $sth->bind_param(1,"$band%",SQL_VARCHAR);
  $sth->execute();
Web Admin's avatar
Web Admin committed
772 773 774 775 776 777 778 779 780 781
  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");
}

782 783 784 785 786 787
=head2 fetch_by_chr_band

 Title   : fetch_by_chr_band
 Usage   :
 Function: create a Slice representing a series of bands
 Example :
788
 Returns : Bio::EnsEMBL::Slice
789
 Args    : the band name
790
 Status     : Stable
791 792 793 794 795 796

=cut

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

Web Admin's avatar
Web Admin committed
797
  my $chr_slice = $self->fetch_by_region('toplevel', $chr);
798 799
  my $seq_region_id = $self->get_seq_region_id($chr_slice);

800
  my $sth = $self->dbc->prepare
Web Admin's avatar
Web Admin committed
801
        ("select min(k.seq_region_start), max(k.seq_region_end) " .
802 803 804
         "from karyotype as k " .
         "where k.seq_region_id = ? and k.band like ?");

805 806 807
  $sth->bind_param(1,$seq_region_id,SQL_INTEGER);
  $sth->bind_param(2,"$band%",SQL_VARCHAR);
  $sth->execute();
808 809 810
  my ( $slice_start, $slice_end) = $sth->fetchrow_array;

  if(defined $slice_start) {
Web Admin's avatar
Web Admin committed
811
    return $self->fetch_by_region('toplevel',$chr,$slice_start,$slice_end);
812 813 814
  }

  throw("Band not recognised in database");
815 816 817 818
}



Web Admin's avatar
Web Admin committed
819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834
=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
835
  Status     : Stable
Web Admin's avatar
Web Admin committed
836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851

=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
852
=head2 fetch_by_transcript_stable_id
853

Graham McVicker's avatar
Graham McVicker committed
854 855 856 857 858 859 860 861 862
  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 
863 864
               number of basepairs on either side of the
               transcript.
Graham McVicker's avatar
Graham McVicker committed
865
  Returntype : Bio::EnsEMBL::Slice
866
  Exceptions : Thrown if the transcript is not in the database.
Graham McVicker's avatar
Graham McVicker committed
867
  Caller     : general
868
  Status     : Stable
869 870 871

=cut

Graham McVicker's avatar
Graham McVicker committed
872
sub fetch_by_transcript_stable_id{
873 874
  my ($self,$transcriptid,$size) = @_;

875 876
  throw('Transcript argument is required.') if(!$transcriptid);

877
  my $ta = $self->db->get_TranscriptAdaptor;
878
  my $transcript = $ta->fetch_by_stable_id($transcriptid);
879

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

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