SliceAdaptor.pm 46 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 97
use Bio::EnsEMBL::Utils::Exception qw(throw deprecate warning);
use Bio::EnsEMBL::Utils::Cache; #CPAN LRU cache
98

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

101
my $SEQ_REGION_CACHE_SIZE = 1000;
102

103 104
sub new {
  my $caller = shift;
105

106
  my $class = ref($caller) || $caller;
107

108
  my $self = $class->SUPER::new(@_);
109

110 111
  my %name_cache;
  my %id_cache;
112

113 114 115 116 117 118 119 120 121 122 123 124
  tie(%name_cache, 'Bio::EnsEMBL::Utils::Cache', $SEQ_REGION_CACHE_SIZE);
  tie(%id_cache, 'Bio::EnsEMBL::Utils::Cache', $SEQ_REGION_CACHE_SIZE);

  $self->{'_name_cache'} = \%name_cache;
  $self->{'_id_cache'} = \%id_cache;

  return $self;
}


=head2 fetch_by_region

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

               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
169
  Caller     : general
170 171 172

=cut

173 174 175 176
sub fetch_by_region {
  my ($self, $coord_system_name, $seq_region_name,
      $start, $end, $strand, $version) = @_;

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

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

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

185 186 187 188 189 190 191 192 193 194 195 196 197
  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;
    }
198
  }
199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223

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

    $key = lc(join(':',$seq_region_name,$cs->name(), $cs->version));
  } 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;
    }
    $constraint .= "ORDER BY cs.rank DESC";
224 225
  }

226 227
  #check the cache so we only go to the db if necessary
  my $length;
228
  my $name_cache = $self->{'_name_cache'};
229

230
  if($key && exists($name_cache->{$key})) {
231 232
    $length = $name_cache->{$key}->[1];
  } else {
233 234
    my $sth = $self->prepare($sql . " WHERE sr.name = ? AND " .
                             $constraint);
Graham McVicker's avatar
Graham McVicker committed
235

236 237 238 239
    # 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);
240

241 242
    if($sth->rows() == 0) {
      $sth->finish();
Graham McVicker's avatar
Graham McVicker committed
243

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

247 248
      $sth = $self->prepare($sql . " WHERE sr.name LIKE ? AND " . $constraint);
      $sth->execute("$seq_region_name.%", @bind_vals);
249 250

      my $prefix_len = length($seq_region_name) + 1;
251 252
      my $high_ver = undef;
      my $high_cs = $cs;
253 254 255 256

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

257 258 259 260 261
      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);
262

263 264
        # cache values for future reference
        $key = lc(join(':',$tmp_name,$tmp_cs->name(),$tmp_cs->version));
265
        $name_cache->{$key}         = [$id,$tmp_length];
266
        $self->{'_id_cache'}->{$id} = [$tmp_name,$tmp_length,$high_cs];
Graham McVicker's avatar
Graham McVicker committed
267

268
        my $tmp_ver = substr($tmp_name, $prefix_len);
269

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

273 274 275 276 277 278 279 280
        # 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
281 282
        }
      }
283 284
      $sth->finish();

285 286
      $cs = $high_cs;

287
      #return if we did not find any appropriate match:
288
      return undef if(!defined($high_ver));
289 290
    } else {

291 292
      my ($id, $cs_id);
      ($seq_region_name, $id, $length, $cs_id) = $sth->fetchrow_array();
293 294
      $sth->finish();

295 296 297 298 299
      if(!$cs) {
        $cs = $csa->fetch_by_dbID($cs_id);
        $key = lc(join(':',$seq_region_name,$cs->name(), $cs->version));
      }

300 301
      #cache results to speed up future queries
      $name_cache->{$key}         = [$id,$length];
302
      $self->{'_id_cache'}->{$id} = [$seq_region_name, $length, $cs];
303
    }
304
  }
305

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

308 309 310 311
  if($end < $start) {
    throw('start [$start] must be less than or equal to end [$end]');
  }

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



323

324
=head2 fetch_by_name
325

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

337 338 339 340
               Returns undef if no seq_region with the provided name exists in
               the database.

  Returntype : Bio::EnsEMBL::Slice or undef
341 342
  Exceptions : throw if incorrent arg provided
  Caller     : Pipeline
343 344 345

=cut

346 347 348
sub fetch_by_name {
  my $self = shift;
  my $name = shift;
Graham McVicker's avatar
Graham McVicker committed
349

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

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

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

  my ($cs_name, $cs_version, $seq_region, $start, $end, $strand) = @array;

  return $self->fetch_by_region($cs_name,$seq_region, $start,
364
                                $end, $strand, $cs_version);
Graham McVicker's avatar
Graham McVicker committed
365
}
Graham McVicker's avatar
Graham McVicker committed
366 367 368



369 370 371 372 373 374 375 376
=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.
377 378
               Returns undef if no such slice exists.
  Returntype : Bio::EnsEMBL::Slice or undef
Graham McVicker's avatar
Graham McVicker committed
379 380
  Exceptions : none
  Caller     : general
Graham McVicker's avatar
Graham McVicker committed
381 382 383

=cut

384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399
sub fetch_by_seq_region_id {
  my ($self, $seq_region_id) = @_;

  my $id_cache = $self->{'_id_cache'};

  my ($name, $length, $cs);

  if(exists $id_cache->{$seq_region_id}) {
    ($name, $length, $cs) = @{$id_cache->{$seq_region_id}};
  } else {
    my $sth = $self->prepare("SELECT name, length, coord_system_id " .
                             "FROM seq_region " .
                             "WHERE seq_region_id = ?");

    $sth->execute($seq_region_id);

400
    return undef if($sth->rows() == 0);
401 402 403 404 405 406 407 408 409 410 411 412 413

    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
    $id_cache->{$seq_region_id} = [$name, $length, $cs];
    my $key = lc(join(':', $name, $cs->name, $cs->version));
    $self->{'_name_cache'}->{$key} = [$seq_region_id, $length];
  }

414 415 416 417 418 419 420
  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
421 422
}

Graham McVicker's avatar
Graham McVicker committed
423

424 425

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

Arne Stabenau's avatar
Arne Stabenau committed
427
  Arg [1]    : Bio::EnsEMBL::Slice $slice
428 429 430 431 432 433 434 435 436 437
               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
Arne Stabenau's avatar
Arne Stabenau committed
438 439

=cut
Graham McVicker's avatar
Graham McVicker committed
440

441 442 443
sub get_seq_region_id {
  my $self = shift;
  my $slice = shift;
Graham McVicker's avatar
Graham McVicker committed
444

445 446 447 448 449 450 451 452 453
  if(!$slice || !ref($slice) || !$slice->isa('Bio::EnsEMBL::Slice')) {
    throw('Slice argument is required');
  }

  my $cs_name = $slice->coord_system->name();
  my $cs_version = $slice->coord_system->version();
  my $seq_region_name = $slice->seq_region_name();

  my $key = lc(join(':', $seq_region_name,$cs_name,$cs_version));
454

455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486
  my $name_cache = $self->{'_name_cache'};

  if(exists($name_cache->{$key})) {
    return $name_cache->{$key}->[0];
  }

  my $csa = $self->db()->get_CoordSystemAdaptor();
  my $coord_system = $csa->fetch_by_name($cs_name,$cs_version);

  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
  $sth->execute("$seq_region_name", $coord_system->dbID());

  if($sth->rows() != 1) {
    throw("Non existant or ambigous seq_region:\n" .
          "  coord_system=[$cs_name],\n" .
          "  name=[$seq_region_name],\n" .
          "  version=[$cs_version]");
  }

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

  #cache information for future requests
  $name_cache->{$key} = [$seq_region_id, $length];
  $self->{'_id_cache'}->{$seq_region_id} =
    [$seq_region_name, $length, $coord_system];

  return $seq_region_id;
Arne Stabenau's avatar
Arne Stabenau committed
487
}
488 489


Arne Stabenau's avatar
Arne Stabenau committed
490

491 492 493 494 495 496 497 498 499
=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
500
  Arg [3]    : bool $include_non_reference (optional)
501 502 503
               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
504 505
  Arg [4]    : int $include_duplicates (optional)
               If set duplicate regions will be returned.
506

507 508 509
  Example    : @chromos = @{$slice_adaptor->fetch_all('chromosome','NCBI33')};
               @contigs = @{$slice_adaptor->fetch_all('contig')};

Graham McVicker's avatar
Graham McVicker committed
510 511 512 513 514
               # 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)};
515 516 517 518 519 520

  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.
521 522
               If the coordinate system with the provided name and version
               does not exist an empty list is returned.
523
               If the coordinate system name provided is 'toplevel', all
524
               non-redundant toplevel slices are returned.
525 526 527 528

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

529
  Returntype : listref of Bio::EnsEMBL::Slices
530
  Exceptions : none
531
  Caller     : general
Web Admin's avatar
Web Admin committed
532 533 534

=cut

535 536 537 538 539
sub fetch_all {
  my $self = shift;
  my $cs_name = shift;
  my $cs_version = shift || '';

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

542 543 544
  #
  # verify existance of requested coord system and get its id
  #
545 546 547 548 549 550
  my $csa       = $self->db->get_CoordSystemAdaptor();
  my $orig_cs   = $csa->fetch_by_name($cs_name, $cs_version);

  return [] if(!$orig_cs);

  my $sth;
551

552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570
  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();
  }
571
  #
572
  # Retrieve the seq_regions from the database
573
  #
574
  if($orig_cs->is_top_level()) {
575 576 577 578 579 580 581 582
   $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");
583 584
    $sth->execute();
  } else {
585 586 587 588 589
     $sth =
       $self->prepare('SELECT seq_region_id, name, length, coord_system_id ' .
                      'FROM   seq_region ' .
                      'WHERE  coord_system_id =?');
     $sth->execute($orig_cs->dbID);
590
  }
591

592 593
  my ($seq_region_id, $name, $length, $cs_id);
  $sth->bind_columns(\$seq_region_id, \$name, \$length, \$cs_id);
594 595 596

  my $name_cache = $self->{'_name_cache'};
  my $id_cache   = $self->{'_id_cache'};
597
  my $cache_count = 0;
598 599 600

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

604
      if(!$cs) {
605
        throw("seq_region $name references non-existent coord_system $cs_id.");
606
      }
607

608
      my $cs_key = lc($cs->name().':'.$cs_version);
609

610 611 612
      #cache values for future reference, but stop adding to the cache once we
      #we know we have filled it up
      if($cache_count < $SEQ_REGION_CACHE_SIZE) {
613 614 615 616
        my $key = lc($name) . ':'. $cs_key;
        $name_cache->{$key} = [$seq_region_id, $length];
        $id_cache->{$seq_region_id} = [$name, $length, $cs];
        $cache_count++;
617
      }
618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635

      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){
        #do not include duplicates
        my @dup = @{$self->fetch_normalized_slice_projection($slice)};
        foreach my $dup_test( @dup){
          if($dup_test->[2]->get_seq_region_id == $slice->get_seq_region_id){
            push @out, $dup_test->[2];
          }
        }
636
      }
637 638
      else {
        push @out, $slice;
639
      }
640
    }
641
  }
642
  $sth->finish();
643

644
  return \@out;
Web Admin's avatar
Web Admin committed
645 646 647
}


648 649 650 651
sub deleteObj {
  my $self = shift;

  $self->SUPER::deleteObj;
652

653 654 655 656 657 658 659
  $self->{'_id_cache'} = {};
  $self->{'_name_cache'} = {};
  $self->{'_exc_cache'} = {};
}



Web Admin's avatar
Web Admin committed
660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689
=head2 fetch_by_band

 Title   : fetch_by_band
 Usage   :
 Function: create a Slice representing a series of bands
 Example :
 Returns :
 Args    : the band name

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

690 691 692 693 694 695 696 697 698 699 700 701 702 703
=head2 fetch_by_chr_band

 Title   : fetch_by_chr_band
 Usage   :
 Function: create a Slice representing a series of bands
 Example :
 Returns :
 Args    : the band name

=cut

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

Web Admin's avatar
Web Admin committed
704
  my $chr_slice = $self->fetch_by_region('toplevel', $chr);
705 706 707 708

  my $seq_region_id = $self->get_seq_region_id($chr_slice);

  my $sth = $self->db->prepare
Web Admin's avatar
Web Admin committed
709
        ("select min(k.seq_region_start), max(k.seq_region_end) " .
710 711 712 713 714 715 716
         "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
717
    return $self->fetch_by_region('toplevel',$chr,$slice_start,$slice_end);
718 719 720
  }

  throw("Band not recognised in database");
721 722 723 724
}



Web Admin's avatar
Web Admin committed
725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756
=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

=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
757
=head2 fetch_by_transcript_stable_id
758

Graham McVicker's avatar
Graham McVicker committed
759 760 761 762 763 764 765 766 767
  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 
768 769
               number of basepairs on either side of the
               transcript.
Graham McVicker's avatar
Graham McVicker committed
770
  Returntype : Bio::EnsEMBL::Slice
771
  Exceptions : Thrown if the transcript is not in the database.
Graham McVicker's avatar
Graham McVicker committed
772
  Caller     : general
773 774 775

=cut

Graham McVicker's avatar
Graham McVicker committed
776
sub fetch_by_transcript_stable_id{
777 778
  my ($self,$transcriptid,$size) = @_;

779 780
  throw('Transcript argument is required.') if(!$transcriptid);

781
  my $ta = $self->db->get_TranscriptAdaptor;
782
  my $transcript = $ta->fetch_by_stable_id($transcriptid);
783

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

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

Graham McVicker's avatar
Graham McVicker committed
789 790


Graham McVicker's avatar
Graham McVicker committed
791 792
=head2 fetch_by_transcript_id

Graham McVicker's avatar
Graham McVicker committed
793 794 795 796 797 798 799
  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);
800 801 802 803
  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
804
  Returntype : Bio::EnsEMBL::Slice
805 806
  Exceptions : throw on incorrect args
               throw if transcript is not in database
Graham McVicker's avatar
Graham McVicker committed
807
  Caller     : general
808 809 810

=cut

Graham McVicker's avatar
Graham McVicker committed
811
sub fetch_by_transcript_id {
812
  my ($self,$transcriptid,$size) = @_;
Graham McVicker's avatar
Graham McVicker committed
813

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

816 817 818 819 820 821
  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);
822 823
}

824 825


Graham McVicker's avatar
Graham McVicker committed
826
=head2 fetch_by_gene_stable_id
827

Graham McVicker's avatar
Graham McVicker committed
828
  Arg [1]    : string $geneid
829
               The stable id of the gene around which the slice is
Graham McVicker's avatar
Graham McVicker committed
830 831
               desired
  Arg [2]    : (optional) int $size
Graham McVicker's avatar
Graham McVicker committed
832
               The length of the flanking regions the slice should encompass
Graham McVicker's avatar
Graham McVicker committed
833 834
               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
835 836 837
  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.
838
               The slice will be created in the genes native coordinate system.
Graham McVicker's avatar
Graham McVicker committed
839
  Returntype : Bio::EnsEMBL::Slice
840 841
  Exceptions : throw on incorrect args
               throw if transcript does not exist
Graham McVicker's avatar
Graham McVicker committed
842
  Caller     : general
843 844 845

=cut

846 847
sub fetch_by_gene_stable_id {
  my ($self,$geneid,$size) = @_;
848

849
  throw('Gene argument is required.') if(!$geneid);
850

851 852
  my $gene_adaptor = $self->db->get_GeneAdaptor();
  my $gene = $gene_adaptor->fetch_by_stable_id($geneid);
853

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

856
  return $self->fetch_by_Feature($gene, $size);
857 858 859
}


860

861
=head2 fetch_by_Feature
Graham McVicker's avatar
Graham McVicker committed
862

863 864 865 866
  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.
867 868
               The size may also be provided as a percentage of the feature 
               size such as 200% or 80.5%.
869 870 871 872 873 874
  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
875
  Returntype : Bio::EnsEMBL::Slice
876 877 878 879
  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
Graham McVicker's avatar
Graham McVicker committed
880 881 882

=cut

883 884
sub fetch_by_Feature{
  my ($self, $feature, $size) = @_;
Graham McVicker's avatar
Graham McVicker committed
885

886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917
  $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;
    }
  }

918 919 920 921
  ## 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
922 923 924

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

925
  #return a new slice covering the region of the feature
926 927 928 929 930 931 932 933
  return Bio::EnsEMBL::Slice->new
    (-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);
934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965
}



=head2 fetch_by_misc_feature_attribute

  Arg [1]    : string $attribute_type
               The code of the attribute type
  Arg [2]    : (optional) string $attribute_value
               The value of the attribute to fetch by
  Arg [3]    : (optional) int $size
               The amount of flanking region around the misc feature desired.
  Example    : $slice = $sa->fetch_by_misc_feature_attribute('superctg',
                                                             'NT_030871');
               $slice = $sa->fetch_by_misc_feature_attribute('synonym',
                                                             'AL00012311',
                                                             $flanking);
  Description: Fetches a slice around a MiscFeature with a particular
               attribute type and value. If no value is specified then
               the feature with the particular attribute is used.
               If no size is specified then 0 is used.
  Returntype : Bio::EnsEMBL::Slice
  Exceptions : Throw if no feature with the specified attribute type and value
               exists in the database
               Warning if multiple features with the specified attribute type
               and value exist in the database.
  Caller     : webcode

=cut

sub fetch_by_misc_feature_attribute {
  my ($self, $attrib_type_code, $attrib_value, $size) = @_;
Graham McVicker's avatar
Graham McVicker committed
966

967
  my $mfa = $self->db()->get_MiscFeatureAdaptor();
Web Admin's avatar
Web Admin committed
968

969 970
  my $feats = $mfa->fetch_all_by_attribute_type_value($attrib_type_code,
                                                   $attrib_value);
Web Admin's avatar
Web Admin committed
971

972 973 974 975
  if(@$feats == 0) {
    throw("MiscFeature with $attrib_type_code=$attrib_value does " .
          "not exist in DB.");
  }
Graham McVicker's avatar
Graham McVicker committed
976

977 978 979 980
  if(@$feats > 1) {
    warning("MiscFeature with $attrib_type_code=$attrib_value is " .
            "ambiguous - using first one found.");
  }
981

982
  my ($feat) = @$feats;
983

984
  return $self->fetch_by_Feature($feat, $size);
Graham McVicker's avatar
Graham McVicker committed
985 986
}

987
=head2 fetch_normalized_slice_projection
Graham McVicker's avatar
Graham McVicker committed
988

989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026 1027 1028 1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045
  Arg [1]    : Bio::EnsEMBL::Slice $slice
  Example    :  ( optional )
  Description: gives back a project style result. The returned slices 
               represent the areas to which there are symlinks for the 
               given slice. start, end show which area on given slice is 
               symlinked
  Returntype : [[start,end,$slice][]]
  Exceptions : none
  Caller     : BaseFeatureAdaptor

=cut


sub fetch_normalized_slice_projection {
  my $self = shift;
  my $slice = shift;

  my $slice_seq_region_id = $self->get_seq_region_id( $slice );

  my $result = $self->{'asm_exc_cache'}->{$slice_seq_region_id};

  if(!defined($result)) {
    my $sql = "
      SELECT seq_region_id, seq_region_start, seq_region_end,
             exc_type, exc_seq_region_id, exc_seq_region_start,
             exc_seq_region_end
        FROM assembly_exception
       WHERE seq_region_id = ?";

    my $sth = $self->prepare( $sql );
    $sth->execute( $slice_seq_region_id );
    $result = $sth->fetchall_arrayref();
    $self->{'asm_exc_cache'}->{$slice_seq_region_id} = $result;
  }

  my (@haps, @pars);

  foreach my $row (@$result) {
    my ( $seq_region_id, $seq_region_start, $seq_region_end,
         $exc_type, $exc_seq_region_id, $exc_seq_region_start,
         $exc_seq_region_end ) = @$row;

    # need overlapping PAR and all HAPs if any
    if( $exc_type eq "PAR" ) {
      if( $seq_region_start <= $slice->end() && 
          $seq_region_end >= $slice->start() ) {
        push( @pars, [ $seq_region_start, $seq_region_end, $exc_seq_region_id,
                       $exc_seq_region_start, $exc_seq_region_end ] );
      }
    } else {
      push( @haps, [ $seq_region_start, $seq_region_end, $exc_seq_region_id,
                     $exc_seq_region_start, $exc_seq_region_end ] );
    }
  }

  if(!@pars && !@haps) {
    #just return this slice, there were no haps or pars
1046
    return  [bless ( [1,$slice->length, $slice], "Bio::EnsEMBL::ProjectionSegment")];
1047 1048 1049 1050 1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087
  }

  my @syms;
  if( @haps > 1 ) {
    my @sort_haps = sort { $a->[1] <=> $b->[1] } @haps;
    throw( "More than one HAP region not supported yet" );
  } elsif( @haps == 1 ) {
    my $hap = $haps[0];

    my $seq_reg_slice = $self->fetch_by_seq_region_id($slice_seq_region_id);
    my $exc_slice = $self->fetch_by_seq_region_id( $hap->[2] );

    #
    # lengths of haplotype and reference in db may be different
    # we want to use the maximum possible length for the mapping
    # between the two systems
    #
    my $len1 = $seq_reg_slice->length();
    my $len2 = $exc_slice->length();
    my $max_len = ($len1 > $len2) ? $len1 : $len2;

    #the inserted region can differ in length, but mapped sections
    #need to be same lengths
    my $diff = $hap->[4] - $hap->[1];

    # we want the region of the haplotype INVERTED
    push( @syms, [ 1, $hap->[0]-1, $hap->[2], 1, $hap->[3] - 1