SliceAdaptor.pm 55.2 KB
Newer Older
1 2

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

# POD documentation - main docs before the code

=head1 NAME

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

=head1 SYNOPSIS

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

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

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



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

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

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


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

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

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

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


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

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

    # do something with slice
  }

75 76
=head1 DESCRIPTION

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

81
=head1 CONTACT
82

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

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

87
=head1 METHODS
88 89 90 91

=cut


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

96

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

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

103

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

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

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

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

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

116 117 118 119
  my $seq_region_cache = $self->db->get_SeqRegionCache();

  $self->{'sr_name_cache'} = $seq_region_cache->{'name_cache'};
  $self->{'sr_id_cache'}   = $seq_region_cache->{'id_cache'};
120 121 122 123 124 125 126

  return $self;
}


=head2 fetch_by_region

127
  Arg [1]    : string $coord_system_name (optional)
128
               The name of the coordinate system of the slice to be created
129
               This may be a name of an actual coordinate system or an alias
130 131 132 133
               to a coordinate system.  Valid aliases are 'seqlevel' or
               'toplevel'.
  Arg [2]    : string $seq_region_name
               The name of the sequence region that the slice will be
Graham McVicker's avatar
Graham McVicker committed
134
               created on.
135 136 137 138 139 140 141 142
  Arg [3]    : int $start (optional, default = 1)
               The start of the slice on the sequence region
  Arg [4]    : int $end (optional, default = seq_region length)
               The end of the slice on the sequence region
  Arg [5]    : int $strand (optional, default = 1)
               The orientation of the slice on the sequence region
  Arg [6]    : string $version (optional, default = default version)
               The version of the coordinate system to use (e.g. NCBI33)
143 144
  Arg [7]    : boolean $no_fuzz (optional, default = undef (false))
               If true (non-zero), do not use "fuzzy matching" (see below).
145 146 147
  Example    : $slice = $slice_adaptor->fetch_by_region('chromosome', 'X');
               $slice = $slice_adaptor->fetch_by_region('clone', 'AC008066.4');
  Description: Retrieves a slice on the requested region.  At a minimum the
148 149 150 151 152 153 154 155 156 157 158
               name the name of the seq_region to fetch must be provided.

               If no coordinate system name is provided than a slice on the
               highest ranked coordinate system with a matching
               seq_region_name will be returned.  If a version but no
               coordinate system name is provided, the same behaviour will
               apply, but only coordinate systems of the appropriate version
               are considered.  The same applies if the 'toplevel' coordinate
               system is specified, however in this case the version is
               ignored.  The coordinate system should always be specified if
               it is known, since this is unambiguous and faster.
159 160 161 162 163 164 165

               Some fuzzy matching is performed if no exact match for
               the provided name is found.  This allows clones to be
               fetched even when their version is not known.  For
               example fetch_by_region('clone', 'AC008066') will
               retrieve the sequence_region with name 'AC008066.4'.

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

169 170 171 172 173 174 175
               If the requested seq_region is not found in the database undef
               is returned.

  Returntype : Bio::EnsEMBL::Slice or undef
  Exceptions : throw if no seq_region_name is provided
               throw if invalid coord_system_name is provided
               throw if start > end is provided
176
  Caller     : general
177
  Status     : Stable
178 179 180

=cut

181 182 183 184

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

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

193
  throw('seq_region_name argument is required') if (!defined($seq_region_name));
194

195 196
  my $cs;
  my $csa = $self->db->get_CoordSystemAdaptor();
197

198
  if ($coord_system_name) {
199 200
    $cs = $csa->fetch_by_name($coord_system_name,$version);

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

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

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

223
  }
224 225 226 227 228 229

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

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

236
	      $constraint = "sr.coord_system_id = ?";
237

238
    $key = "$seq_region_name:".$cs->dbID();
239 240 241 242 243 244 245 246 247 248
  } 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;
    }
249
    $constraint .= "ORDER BY cs.rank ASC";
250 251
  }

252
  # check the cache so we only go to the db if necessary
253
  my $length;
254
  my $arr;
255
  if ($key) {
256
    $arr = $self->{'sr_name_cache'}->{$key};
257
  }
258

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

    $sth->execute();
272

273
    if ($sth->rows() == 0) {
274
      $sth->finish();
Graham McVicker's avatar
Graham McVicker committed
275

276 277
      if ($no_fuzz) { return undef }

278 279
      # do fuzzy matching, assuming that we are just missing a version on
      # the end of the seq_region name
280
      
281
      $sth = $self->prepare($sql . " WHERE sr.name LIKE ? AND " . $constraint);
282

283
      $sth->bind_param(1, "$seq_region_name.%", SQL_VARCHAR);
284 285

      if ($cs){
286
	  $sth->bind_param(2, $cs->dbID, SQL_INTEGER);
287 288
      }
      else{
289
	  $sth->bind_param(2, $version, SQL_VARCHAR) if ($version);	
290
      }
291
      $sth->execute();
292 293

      my $prefix_len = length($seq_region_name) + 1;
294 295
      my $high_ver = undef;
      my $high_cs = $cs;
296 297 298 299

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

300 301 302
      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
303
      my $i = 0;
304 305

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

308
        # cache values for future reference
309
        my $arr = [ $id, $tmp_name, $cs_id, $tmp_length ];
310 311
        $self->{'sr_name_cache'}->{"$tmp_name:$cs_id"} = $arr;
        $self->{'sr_id_cache'}->{"$id"} = $arr;
Graham McVicker's avatar
Graham McVicker committed
312

313
        my $tmp_ver = substr($tmp_name, $prefix_len);
314

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

318 319 320 321 322 323 324 325
        # 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
326
        }
327 328

        $i++;
Graham McVicker's avatar
Graham McVicker committed
329
      }
330 331
      $sth->finish();

332 333 334 335 336
      # 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");
      }

337 338
      $cs = $high_cs;

339
      #return if we did not find any appropriate match:
340 341
      return undef if (!defined($high_ver));

342 343
    } else {

344 345
      my ($id, $cs_id);
      ($seq_region_name, $id, $length, $cs_id) = $sth->fetchrow_array();
346
      $sth->finish();
347

348 349
      # cahce to speed up for future queries
      my $arr = [ $id, $seq_region_name, $cs_id, $length ];
350 351
      $self->{'sr_name_cache'}->{"$seq_region_name:$cs_id"} = $arr;
      $self->{'sr_id_cache'}->{"$id"} = $arr;
352
      $cs = $csa->fetch_by_dbID( $cs_id );
353
    }
354
  }
355

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

358 359
  if($end+1 < $start) {
    throw("start [$start] must be less than or equal to end+1 [$end+1]");
360 361
  }

362
  return Bio::EnsEMBL::Slice->new(-COORD_SYSTEM      => $cs,
363 364 365 366 367 368
                                  -SEQ_REGION_NAME   => $seq_region_name,
                                  -SEQ_REGION_LENGTH => $length,
                                  -START             => $start,
                                  -END               => $end,
                                  -STRAND            => $strand,
                                  -ADAPTOR           => $self);
369 370 371 372
}



373

374
=head2 fetch_by_name
375

Graham McVicker's avatar
Graham McVicker committed
376
  Arg [1]    : string $name
377 378 379 380 381 382 383 384 385
  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
386

387 388 389 390
               Returns undef if no seq_region with the provided name exists in
               the database.

  Returntype : Bio::EnsEMBL::Slice or undef
391 392
  Exceptions : throw if incorrent arg provided
  Caller     : Pipeline
393
  Status     : Stable
394 395 396

=cut

397 398 399
sub fetch_by_name {
  my $self = shift;
  my $name = shift;
Graham McVicker's avatar
Graham McVicker committed
400

401 402
  if(!$name) {
    throw("name argument is required");
Graham McVicker's avatar
Graham McVicker committed
403
  }
404 405 406

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

407
  if(scalar(@array) < 3 || scalar(@array) > 6) {
408
    throw("Malformed slice name [$name].  Format is " .
Laura Clarke's avatar
 
Laura Clarke committed
409
        "coord_system:version:name:start:end:strand");
410 411
  }

412
  # Rearrange arguments to suit fetch_by_region
413

414 415 416 417 418 419 420 421 422
  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
423
}
Graham McVicker's avatar
Graham McVicker committed
424 425 426



427 428 429 430 431
=head2 fetch_by_seq_region_id

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

=cut

446
sub fetch_by_seq_region_id {
447
  my ($self, $seq_region_id,$start,$end,$strand) = @_;
448

449
  my $arr = $self->{'sr_id_cache'}->{ $seq_region_id };
450 451
  my ($name, $length, $cs);

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

461 462
    $sth->bind_param(1,$seq_region_id,SQL_INTEGER);
    $sth->execute();
463

464
    return undef if($sth->rows() == 0);
465 466 467 468 469 470 471 472

    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
473
    my $arr = [ $seq_region_id, $name, $cs_id, $length ];
474 475 476

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

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

Graham McVicker's avatar
Graham McVicker committed
488

489 490

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

Arne Stabenau's avatar
Arne Stabenau committed
492
  Arg [1]    : Bio::EnsEMBL::Slice $slice
493 494 495 496 497 498 499 500 501 502
               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
503
  Status     : Stable
Arne Stabenau's avatar
Arne Stabenau committed
504 505

=cut
Graham McVicker's avatar
Graham McVicker committed
506

507 508 509
sub get_seq_region_id {
  my $self = shift;
  my $slice = shift;
Graham McVicker's avatar
Graham McVicker committed
510

511 512 513
  if(!$slice || !ref($slice) || !$slice->isa('Bio::EnsEMBL::Slice')) {
    throw('Slice argument is required');
  }
514
  
515
  my $seq_region_name = $slice->seq_region_name();
516
  my $key = $seq_region_name.":".$slice->coord_system->dbID();
517
  my $arr = $self->{'sr_name_cache'}->{"$key"};
518

519 520
  if( $arr ) {
    return $arr->[0];
521 522
  }

523
  my $cs_id = $slice->coord_system->dbID();
524 525 526 527 528 529

  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
530 531 532
  $sth->bind_param(1,"$seq_region_name",SQL_VARCHAR);
  $sth->bind_param(2,$cs_id,SQL_INTEGER);
  $sth->execute();
533 534 535

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

539 540 541 542 543 544
  }

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

  #cache information for future requests
545
  $arr = [ $seq_region_id, $seq_region_name, $cs_id, $length ];
546 547 548 549

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

550
  return $seq_region_id;
Arne Stabenau's avatar
Arne Stabenau committed
551
}
552 553


Arne Stabenau's avatar
Arne Stabenau committed
554

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

576 577 578
  Example    : @chromos = @{$slice_adaptor->fetch_all('chromosome','NCBI33')};
               @contigs = @{$slice_adaptor->fetch_all('contig')};

Graham McVicker's avatar
Graham McVicker committed
579 580 581 582 583
               # 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)};
584 585 586 587 588 589

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

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

599
  Returntype : listref of Bio::EnsEMBL::Slices
600
  Exceptions : none
601
  Caller     : general
602
  Status     : Stable
Web Admin's avatar
Web Admin committed
603 604 605

=cut

606 607 608 609 610
sub fetch_all {
  my $self = shift;
  my $cs_name = shift;
  my $cs_version = shift || '';

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

613 614 615
  #
  # verify existance of requested coord system and get its id
  #
616 617 618 619 620 621
  my $csa       = $self->db->get_CoordSystemAdaptor();
  my $orig_cs   = $csa->fetch_by_name($cs_name, $cs_version);

  return [] if(!$orig_cs);

  my $sth;
622

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

664 665
  my ($seq_region_id, $name, $length, $cs_id);
  $sth->bind_columns(\$seq_region_id, \$name, \$length, \$cs_id);
666

667
  my $cache_count = 0;
668 669 670

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

674
      if(!$cs) {
675
        throw("seq_region $name references non-existent coord_system $cs_id.");
676
      }
677

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

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

686
        $cache_count++;
687
      }
688 689 690 691 692 693 694 695 696 697 698

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

724
  return \@out;
Web Admin's avatar
Web Admin committed
725 726
}

727 728 729 730 731 732 733 734 735 736 737 738 739
=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
740

741 742 743 744 745 746 747
  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
748
    $sth->finish;
749 750
    return 1;
  }
Ian Longden's avatar
Ian Longden committed
751
  $sth->finish;
752 753
  return 0;
}
754

Web Admin's avatar
Web Admin committed
755 756 757 758
=head2 fetch_by_band

 Title   : fetch_by_band
 Usage   :
759
 Function: Does not work please use fetch_by_chr_band
Web Admin's avatar
Web Admin committed
760
 Example :
761
 Returns : Bio::EnsEMBL::Slice
Web Admin's avatar
Web Admin committed
762
 Args    : the band name
763
 Status     : AT RISK
Web Admin's avatar
Web Admin committed
764 765 766 767 768 769

=cut

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

770
  my $sth = $self->dbc->prepare
Web Admin's avatar
Web Admin committed
771 772 773 774
        ("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");

775 776
  $sth->bind_param(1,"$band%",SQL_VARCHAR);
  $sth->execute();
Web Admin's avatar
Web Admin committed
777 778 779 780 781 782 783 784 785 786
  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");
}

787 788 789 790 791 792
=head2 fetch_by_chr_band

 Title   : fetch_by_chr_band
 Usage   :
 Function: create a Slice representing a series of bands
 Example :
793
 Returns : Bio::EnsEMBL::Slice
794
 Args    : the band name
795
 Status     : Stable
796 797 798 799 800 801

=cut

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

Web Admin's avatar
Web Admin committed
802
  my $chr_slice = $self->fetch_by_region('toplevel', $chr);
803 804
  my $seq_region_id = $self->get_seq_region_id($chr_slice);

805
  my $sth = $self->dbc->prepare
Web Admin's avatar
Web Admin committed
806
        ("select min(k.seq_region_start), max(k.seq_region_end) " .
807 808 809
         "from karyotype as k " .
         "where k.seq_region_id = ? and k.band like ?");

810 811 812
  $sth->bind_param(1,$seq_region_id,SQL_INTEGER);
  $sth->bind_param(2,"$band%",SQL_VARCHAR);
  $sth->execute();
813 814 815
  my ( $slice_start, $slice_end) = $sth->fetchrow_array;

  if(defined $slice_start) {
Web Admin's avatar
Web Admin committed
816
    return $self->fetch_by_region('toplevel',$chr,$slice_start,$slice_end);
817 818 819
  }

  throw("Band not recognised in database");
820 821 822 823
}



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

=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
857
=head2 fetch_by_transcript_stable_id
858

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

=cut

Graham McVicker's avatar
Graham McVicker committed
877
sub fetch_by_transcript_stable_id{
878 879
  my ($self,$transcriptid,$size) = @_;

880 881
  throw('Transcript argument is required.') if(!$transcriptid);

882
  my $ta = $self->db->get_TranscriptAdaptor;
883
  my $transcript = $ta->fetch_by_stable_id($transcriptid);
884

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

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

Graham McVicker's avatar
Graham McVicker committed
890 891