SliceAdaptor.pm 87.8 KB
Newer Older
1
=head1 LICENSE
2

3
Copyright [1999-2013] Wellcome Trust Sanger Institute and the EMBL-European Bioinformatics Institute
4

5 6 7 8 9 10 11 12 13 14 15 16 17
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

     http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.

=cut
18 19 20 21 22


=head1 CONTACT

  Please email comments or questions to the public Ensembl
Magali Ruffier's avatar
Magali Ruffier committed
23
  developers list at <http://lists.ensembl.org/mailman/listinfo/dev>.
24 25

  Questions may also be sent to the Ensembl help desk at
Magali Ruffier's avatar
Magali Ruffier committed
26
  <http://www.ensembl.org/Help/Contact>.
27 28

=cut
29 30 31

=head1 NAME

32 33
Bio::EnsEMBL::DBSQL::SliceAdaptor - A database aware adaptor responsible for
the creation of Slice objects.
34 35 36

=head1 SYNOPSIS

Graham McVicker's avatar
Graham McVicker committed
37
  use Bio::EnsEMBL::Utils::Slice qw(split_Slices);
38 39 40
  use Bio::EnsEMBL::Registry;

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

45 46
  $slice_adaptor =
    Bio::EnsEMBL::Registry->get_adaptor( "human", "core", "slice" );
47

Graham McVicker's avatar
Graham McVicker committed
48
  # get a slice on the entire chromosome X
49
  $chr_slice = $slice_adaptor->fetch_by_region( 'chromosome', 'X' );
Graham McVicker's avatar
Graham McVicker committed
50 51

  # get a slice for each clone in the database
52 53
  foreach $cln_slice ( @{ $slice_adaptor->fetch_all('clone') } ) {
    # do something with clone
54 55
  }

Graham McVicker's avatar
Graham McVicker committed
56
  # get a slice which is part of NT_004321
57 58 59
  $spctg_slice =
    $slice_adaptor->fetch_by_region( 'supercontig', 'NT_004321',
    200_000, 600_000 );
60

61 62
  # get all non-redundant slices from the highest possible coordinate
  # systems
Graham McVicker's avatar
Graham McVicker committed
63 64 65
  $slices = $slice_adaptor->fetch_all('toplevel');

  # include non-reference regions
66
  $slices = $slice_adaptor->fetch_all( 'toplevel', undef, 1 );
Graham McVicker's avatar
Graham McVicker committed
67 68

  # include non-duplicate regions
69
  $slices = $slice_adaptor->fetch_all( 'toplevel', undef, 0, 1 );
Graham McVicker's avatar
Graham McVicker committed
70 71

  # split up a list of slices into smaller slices
72
  $overlap    = 1000;
Graham McVicker's avatar
Graham McVicker committed
73
  $max_length = 1e6;
74
  $slices     = split_Slices( $slices, $max_length, $overlap );
Graham McVicker's avatar
Graham McVicker committed
75 76

  # store a list of slice names in a file
77
  open( FILE, ">$filename" ) or die("Could not open file $filename");
Graham McVicker's avatar
Graham McVicker committed
78 79 80 81 82 83
  foreach my $slice (@$slices) {
    print FILE $slice->name(), "\n";
  }
  close FILE;

  # retreive a list of slices from a file
84 85
  open( FILE, $filename ) or die("Could not open file $filename");
  while ( $name = <FILE> ) {
Graham McVicker's avatar
Graham McVicker committed
86 87 88 89 90
    chomp($name);
    $slice = $slice_adaptor->fetch_by_name($name);
    # do something with slice
  }

91 92
=head1 DESCRIPTION

93 94
This module is responsible for fetching Slices representing genomic
regions from a database.  A Details on how slices can be used are in the
95
Bio::EnsEMBL::Slice module.
96

97
=head1 METHODS
98 99 100 101

=cut


102
package Bio::EnsEMBL::DBSQL::SliceAdaptor;
103 104 105
use vars qw(@ISA);
use strict;

106

107 108
use Bio::EnsEMBL::DBSQL::BaseAdaptor;
use Bio::EnsEMBL::Slice;
109
use Bio::EnsEMBL::CircularSlice;
110
use Bio::EnsEMBL::Mapper;
Ian Longden's avatar
Ian Longden committed
111
use Bio::EnsEMBL::LRGSlice;
112
use Bio::EnsEMBL::Utils::Exception qw(throw deprecate warning stack_trace_dump);
113
use Bio::EnsEMBL::ProjectionSegment;
Andy Yates's avatar
Andy Yates committed
114
use Scalar::Util qw/looks_like_number/;
115
use Bio::EnsEMBL::Utils::Scalar qw/assert_integer/;
116

117
@ISA = ('Bio::EnsEMBL::DBSQL::BaseAdaptor');
118

119 120
sub new {
  my $caller = shift;
121

122
  my $class = ref($caller) || $caller;
123

124
  my $self = $class->SUPER::new(@_);
125

126 127
  # use a cache which is shared and also used by the assembly
  # mapper adaptor
128

129 130 131 132
  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'};
133

Ian Longden's avatar
Ian Longden committed
134 135 136
  $self->{'lrg_region_test'} = undef;
  my $meta_container = $self->db->get_MetaContainer();
  my @values = $meta_container->list_value_by_key("LRG");
137
  if(scalar(@values) and $values[0]->[0]){
Ian Longden's avatar
Ian Longden committed
138 139
    $self->{'lrg_region_test'} = $values[0]->[0];
  }
140 141 142 143 144 145
  return $self;
}


=head2 fetch_by_region

146
  Arg [1]    : string $coord_system_name (optional)
147
               The name of the coordinate system of the slice to be created
148
               This may be a name of an actual coordinate system or an alias
149 150 151 152
               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
153
               created on.
154 155 156 157 158 159 160 161
  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)
162 163
  Arg [7]    : boolean $no_fuzz (optional, default = undef (false))
               If true (non-zero), do not use "fuzzy matching" (see below).
164 165 166
  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
167 168 169 170 171 172 173 174 175 176 177
               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.
178 179 180 181 182 183 184

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

185 186 187
               The fuzzy matching can be turned off by setting the
               $no_fuzz argument to a true value.

188 189 190 191 192 193 194
               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
195
  Caller     : general
196
  Status     : Stable
197 198 199

=cut

200 201 202 203

#
# ARNE: This subroutine needs simplification!! 
#
204
sub fetch_by_region {
205 206 207
  my ( $self, $coord_system_name, $seq_region_name, $start, $end,
       $strand, $version, $no_fuzz )
    = @_;
208

209 210 211
  assert_integer($start, 'start') if $start;
  assert_integer($end, 'end') if $end;

212 213
  if ( !defined($start) )  { $start  = 1 }
  if ( !defined($strand) ) { $strand = 1 }
214

215 216 217
  if ( !defined($seq_region_name) ) {
    throw('seq_region_name argument is required');
  }
218

219 220
  my $cs;
  my $csa = $self->db->get_CoordSystemAdaptor();
221

222 223
  if ( defined($coord_system_name) ) {
    $cs = $csa->fetch_by_name( $coord_system_name, $version );
224

225
    ## REMOVE THESE THREE LINES WHEN STICKLEBACK DB IS FIXED!
226 227 228 229 230 231 232
    ## 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!

233 234
    if ( !defined($cs) && $coord_system_name eq 'chromosome' ) {
      $cs = $csa->fetch_by_name( 'group', $version );
235 236
    }

237
    if ( !defined($cs) ) {
238
      throw( sprintf( "Unknown coordinate system:\n"
239
                        . "name='%s' version='%s'\n",
240
                      $coord_system_name, $version ) );
241
    }
242

243
    # fetching by toplevel is same as fetching w/o name or version
244 245
    if ( $cs->is_top_level() ) {
      $cs      = undef;
246 247
      $version = undef;
    }
248

249
  } ## end if ( defined($coord_system_name...))
250 251 252

  my $constraint;
  my $sql;
253
  my @bind_params;
254 255
  my $key;

256 257 258 259
  if ( defined($cs) ) {
    $sql = sprintf( "SELECT sr.name, sr.seq_region_id, sr.length, %d "
                      . "FROM seq_region sr ",
                    $cs->dbID() );
260

261 262
    $constraint = "AND sr.coord_system_id = ?";
    push( @bind_params, [ $cs->dbID(), SQL_INTEGER ] );
263

264
    $key = "$seq_region_name:" . $cs->dbID();
265
  } else {
266 267 268
    $sql =
      "SELECT sr.name, sr.seq_region_id, sr.length, cs.coord_system_id "
      . "FROM seq_region sr, coord_system cs ";
269

270 271 272 273 274
    $constraint = "AND sr.coord_system_id = cs.coord_system_id "
      . "AND cs.species_id = ? ";
    push( @bind_params, [ $self->species_id(), SQL_INTEGER ] );

    if ( defined($version) ) {
275
      $constraint .= "AND cs.version = ? ";
276
      push( @bind_params, [ $version, SQL_VARCHAR ] );
277
    }
278

279
    $constraint .= "ORDER BY cs.rank ASC";
280 281
  }

282
  # check the cache so we only go to the db if necessary
283
  my $length;
284
  my $arr;
285

286 287 288
  if ( defined($key) ) { $arr = $self->{'sr_name_cache'}->{$key} }

  if ( defined($arr) ) {
289
    $length = $arr->[3];
290
  } else {
291 292 293 294 295 296 297 298
    my $sth =
      $self->prepare( $sql . "WHERE sr.name = ? " . $constraint );

    unshift( @bind_params, [ $seq_region_name, SQL_VARCHAR ] );

    my $pos = 0;
    foreach my $param (@bind_params) {
      $sth->bind_param( ++$pos, $param->[0], $param->[1] );
299
    }
300 301

    $sth->execute();
302

303
    if ( $sth->rows() == 0 ) {
304
      $sth->finish();
Graham McVicker's avatar
Graham McVicker committed
305

306 307

      # try synonyms
308
      my $syn_sql_sth = $self->prepare("select s.name, cs.name, cs.version from seq_region s join seq_region_synonym ss using (seq_region_id) join coord_system cs using (coord_system_id) where ss.synonym = ? and cs.species_id =?");
309
      $syn_sql_sth->bind_param(1, $seq_region_name, SQL_VARCHAR);
310
      $syn_sql_sth->bind_param(2, $self->species_id(), SQL_INTEGER);
311
      $syn_sql_sth->execute();
312 313
      my ($new_name, $new_coord_system, $new_version);
      $syn_sql_sth->bind_columns( \$new_name, \$new_coord_system, \$new_version);
314
            
315
      if($syn_sql_sth->fetch){
316
        $syn_sql_sth->finish;
317
        if ((not defined($cs)) || ($cs->name eq $new_coord_system)) {
318
            return $self->fetch_by_region($new_coord_system, $new_name, $start, $end, $strand, $new_version, $no_fuzz);
319
        } elsif ($cs->name ne $new_coord_system) {
320 321 322 323 324
            warning("Searched for a known feature on coordinate system: ".$cs->dbID." but found it on: ".$new_coord_system.
            "\n No result returned, consider searching without coordinate system or use toplevel.");
            return;
        }
        
325 326 327 328
      }
      $syn_sql_sth->finish;


329 330
      if ($no_fuzz) { return undef }

331 332
      # Do fuzzy matching, assuming that we are just missing a version
      # on the end of the seq_region name.
333

334 335
      $sth =
        $self->prepare( $sql . " WHERE sr.name LIKE ? " . $constraint );
336

337 338 339 340 341 342
      $bind_params[0] =
        [ sprintf( '%s.%%', $seq_region_name ), SQL_VARCHAR ];

      $pos = 0;
      foreach my $param (@bind_params) {
        $sth->bind_param( ++$pos, $param->[0], $param->[1] );
343
      }
344

345
      $sth->execute();
346 347

      my $prefix_len = length($seq_region_name) + 1;
348 349
      my $high_ver   = undef;
      my $high_cs    = $cs;
350

351 352
      # Find the fuzzy-matched seq_region with the highest postfix
      # (which ought to be a version).
353

354 355
      my ( $tmp_name, $id, $tmp_length, $cs_id );
      $sth->bind_columns( \( $tmp_name, $id, $tmp_length, $cs_id ) );
356

Felix Kokocinski's avatar
Felix Kokocinski committed
357
      my $i = 0;
358

359 360 361
      while ( $sth->fetch ) {
        my $tmp_cs =
          ( defined($cs) ? $cs : $csa->fetch_by_dbID($cs_id) );
362

363
        # cache values for future reference
364
        my $arr = [ $id, $tmp_name, $cs_id, $tmp_length ];
365
        $self->{'sr_name_cache'}->{"$tmp_name:$cs_id"} = $arr;
366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384
        $self->{'sr_id_cache'}->{"$id"}                = $arr;

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

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

        # 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
385
        }
386 387

        $i++;
388
      } ## end while ( $sth->fetch )
389 390
      $sth->finish();

391
      # warn if fuzzy matching found more than one result
392 393 394 395 396 397 398
      if ( $i > 1 ) {
        warning(
          sprintf(
            "Fuzzy matching of seq_region_name "
              . "returned more than one result.\n"
              . "You might want to check whether the returned seq_region\n"
              . "(%s:%s) is the one you intended to fetch.\n",
399
            $high_cs->name(), $seq_region_name ) );
400 401
      }

402 403
      $cs = $high_cs;

404 405
      # return if we did not find any appropriate match:
      if ( !defined($high_ver) ) { return undef }
406

407 408
    } else {

409 410 411
      my ( $id, $cs_id );
      ( $seq_region_name, $id, $length, $cs_id ) =
        $sth->fetchrow_array();
412
      $sth->finish();
413

414
      # cache to speed up for future queries
415
      my $arr = [ $id, $seq_region_name, $cs_id, $length ];
416
      $self->{'sr_name_cache'}->{"$seq_region_name:$cs_id"} = $arr;
417 418
      $self->{'sr_id_cache'}->{"$id"}                       = $arr;
      $cs = $csa->fetch_by_dbID($cs_id);
419
    }
420
  } ## end else [ if ( defined($arr) ) ]
421

422
  if ( !defined($end) ) { $end = $length }
Graham McVicker's avatar
Graham McVicker committed
423

424 425
  #If this was given then check if we've got a circular seq region otherwise
  #let it fall through to the normal Slice method
426
  if ( $end + 1 < $start ) {
427
    my $cs_id = $cs->dbID();
428
    my $seq_region_id = $self->{'sr_name_cache'}->{"$seq_region_name:$cs_id"}->[0];
429 430 431 432 433 434 435 436 437 438 439 440 441
    if($self->is_circular($seq_region_id)) {
      my $new_sl =
        Bio::EnsEMBL::CircularSlice->new(
                                     -COORD_SYSTEM    => $cs,
                                     -SEQ_REGION_NAME => $seq_region_name,
                                     -SEQ_REGION_LENGTH => $length,
                                     -START             => $start,
                                     -END               => $end,
                                     -STRAND            => 1,
                                     -ADAPTOR           => $self );
  
      return $new_sl;
    }
442 443
  }

444 445 446
  if ( defined( $self->{'lrg_region_test'} )
       and substr( $cs->name, 0, 3 ) eq $self->{'lrg_region_test'} )
  {
Ian Longden's avatar
Ian Longden committed
447
    return
448 449 450 451 452 453 454 455
      Bio::EnsEMBL::LRGSlice->new( -COORD_SYSTEM    => $cs,
                                   -SEQ_REGION_NAME => $seq_region_name,
                                   -SEQ_REGION_LENGTH => $length,
                                   -START             => $start,
                                   -END               => $end,
                                   -STRAND            => $strand,
                                   -ADAPTOR           => $self );
  } else {
Ian Longden's avatar
Ian Longden committed
456
    return
457 458 459 460 461 462 463 464
      Bio::EnsEMBL::Slice->new_fast( {
                                  'coord_system'    => $cs,
                                  'seq_region_name' => $seq_region_name,
                                  'seq_region_length' => $length,
                                  'start'             => $start,
                                  'end'               => $end,
                                  'strand'            => $strand,
                                  'adaptor'           => $self } );
Ian Longden's avatar
Ian Longden committed
465
  }
466
} ## end sub fetch_by_region
467

468 469 470
=head2 fetch_by_toplevel_location

  Arg [1]     : string $location
471
                Ensembl formatted location. Can be a format like 
Andy Yates's avatar
Andy Yates committed
472 473 474 475 476 477
                C<name:start-end>, C<name:start..end>, C<name:start:end>, 
                C<name:start>, C<name>. We can also support strand 
                specification as a +/- or 1/-1. 
                
                Location names must be separated by a C<:>. All others can be
                separated by C<..>, C<:> or C<->.
478 479
  Arg[2]      : boolean $no_warnings
                Suppress warnings from this method
480 481
  Arg[3]      : boolean $no_fuzz
                Stop fuzzy matching of sequence regions from occuring
482 483 484 485
  Arg[4]      : boolean $ucsc
                If we are unsuccessful at retriving a location retry taking any 
                possible chr prefix into account e.g. chrX and X are treated as
                equivalents
486
  Example     : my $slice = $sa->fetch_by_toplevel_location('X:1-10000')
Andy Yates's avatar
Andy Yates committed
487
                my $slice = $sa->fetch_by_toplevel_location('X:1-10000:-1')
488 489 490
  Description : Converts an Ensembl location/region into the sequence region
                name, start and end and passes them onto C<fetch_by_region()>. 
                The code assumes that the required slice is on the top level
491 492
                coordinate system. The code assumes that location formatting
                is not perfect and will perform basic cleanup before parsing.
493
  Returntype  : Bio::EnsEMBL::Slice
494 495
  Exceptions  : If $location is false otherwise see C<fetch_by_location()>
                or C<fetch_by_region()>
496 497 498 499 500 501
  Caller      : General
  Status      : Beta

=cut

sub fetch_by_toplevel_location {
502 503
  my ($self, $location, $no_warnings, $no_fuzz, $ucsc) = @_;
  return $self->fetch_by_location($location, 'toplevel', undef, $no_warnings, $no_fuzz, $ucsc);
504 505 506 507 508 509 510 511 512 513 514
}

=head2 fetch_by_location

  Arg [1]     : string $location
                Ensembl formatted location. Can be a format like 
                C<name:start-end>, C<name:start..end>, C<name:start:end>, 
                C<name:start>, C<name>. We can also support strand 
                specification as a +/- or 1/-1. 
                
                Location names must be separated by a C<:>. All others can be
515
                separated by C<..>, C<:>, C<_> or C<->.
516 517 518 519 520 521 522 523
  Arg[2]      : String $coord_system_name
                The coordinate system to retrieve
  Arg[3]      : String $coord_system_version
                Optional parameter. Version of the coordinate system to fetch
  Arg[4]      : boolean $no_warnings
                Suppress warnings from this method
  Arg[5]      : boolean $no_fuzz
                Stop fuzzy matching of sequence regions from occuring
524 525 526 527
  Arg[6]      : boolean $ucsc
                If we are unsuccessful at retriving a location retry taking any 
                possible chr prefix into account e.g. chrX and X are treated as
                equivalents
528 529
  Example     : my $slice = $sa->fetch_by_location('X:1-10000','chromosome')
                my $slice = $sa->fetch_by_location('X:1-10000:-1','toplevel')
530 531
  Description : Converts an Ensembl location/region into the sequence region
                name, start and end and passes them onto C<fetch_by_region()>. 
532 533
                The code assumes that location formatting is not perfect and 
                will perform basic cleanup before parsing.
534 535 536 537 538 539 540
  Returntype  : Bio::EnsEMBL::Slice
  Exceptions  : If $location or coordinate system is false otherwise 
                see C<fetch_by_region()>
  Caller      : General
  Status      : Beta

=cut
541

542
sub fetch_by_location {
543
  my ($self, $location, $coord_system_name, $coord_system_version, $no_warnings, $no_fuzz, $ucsc) = @_;
544 545 546
  
  throw "No coordinate system name specified" unless $coord_system_name;
  
547 548 549 550 551 552 553 554 555 556
  my ($seq_region_name, $start, $end, $strand) = $self->parse_location_to_values($location, $no_warnings);

  if(! $seq_region_name) {
    return;
  }
    
  if(defined $start && defined $end && $start > $end) {
    throw "Cannot request a slice whose start is greater than its end. Start: $start. End: $end";
  }
  
557
  my $slice = $self->fetch_by_region($coord_system_name, $seq_region_name, $start, $end, $strand, $coord_system_version, $no_fuzz);
558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573
  if(! defined $slice) {
    if($ucsc) {
      my $ucsc_seq_region_name = $seq_region_name;
      $ucsc_seq_region_name =~ s/^chr//;
      if($ucsc_seq_region_name ne $seq_region_name) {
        $slice = $self->fetch_by_region($coord_system_name, $ucsc_seq_region_name,  $start, $end, $strand, $coord_system_version, $no_fuzz);
        return if ! defined $slice; #if we had no slice still then bail
      }
      else {
        return; #If it was not different then we didn't have the prefix so just return (same bail as before)
      }
    }
    else {
      return; #We didn't have a slice and no UCSC specifics are being triggered
    }
  }
574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596
  
  my $srl = $slice->seq_region_length();
  my $name = $slice->seq_region_name();
  if(defined $start && $start > $srl) {
    throw "Cannot request a slice whose start ($start) is greater than $srl for $name.";
  }
  if(defined $end && $end > $srl) {
    warning "Requested end ($end) is greater than $srl for $name. Resetting to $srl" if ! $no_warnings;
    $slice->{end} = $srl;
  }
  
  return $slice;
}

=head2 parse_location_to_values

  Arg [1]     : string $location
                Ensembl formatted location. Can be a format like 
                C<name:start-end>, C<name:start..end>, C<name:start:end>, 
                C<name:start>, C<name>. We can also support strand 
                specification as a +/- or 1/-1. 
                
                Location names must be separated by a C<:>. All others can be
597
                separated by C<..>, C<:> C<_>, or C<->.
598 599
  Arg[2]      : boolean $no_warnings
                Suppress warnings from this method
600 601
  Arg[3]      : boolean $no_errors
                Supress errors being thrown from this method
602 603 604 605 606 607 608 609 610
  Example			: my ($name, $start, $end, $strand) = $sa->parse_location_to_values('X:1..100:1);
  Description	: Takes in an Ensembl location String and returns the parsed
                values
  Returntype 	: List. Contains name, start, end and strand 

=cut


sub parse_location_to_values {
611
  my ($self, $location, $no_warnings, $no_errors) = @_;
612 613 614
  
  throw 'You must specify a location' if ! $location;
  
615 616 617 618
  #cleanup any nomenclature like 1 000 or 1,000
  my $number_seps_regex = qr/\s+|,/;
  my $separator_regex = qr/(?:-|[.]{2}|\:|_)?/;
  my $number_regex = qr/[0-9, E]+/xms;
619
  my $strand_regex = qr/[+-1]|-1/xms;
620
  
621
  my $regex = qr/^((?:\w|\.|_|-)+) \s* :? \s* ($number_regex)? $separator_regex ($number_regex)? $separator_regex ($strand_regex)? $/xms;
622 623
  my ($seq_region_name, $start, $end, $strand);
  if(($seq_region_name, $start, $end, $strand) = $location =~ $regex) {
Andy Yates's avatar
Andy Yates committed
624 625 626 627 628 629 630
    
    if(defined $strand) {
      if(!looks_like_number($strand)) {
        $strand = ($strand eq '+') ? 1 : -1;
      }
    }
    
631
    if(defined $start) {
632
      $start =~ s/$number_seps_regex//g; 
633 634 635 636
      if($start < 1) {
        warning "Start was less than 1 (${start}) which is not allowed. Resetting to 1"  if ! $no_warnings;
        $start = 1;
      }
637
    }
638
    if(defined $end) {
639
      $end =~ s/$number_seps_regex//g;
640
      if($end < 1) {
641
        throw "Cannot request negative or 0 end indexes through this interface. Given $end but expected something greater than 0" unless $no_errors;
642
      }
643 644 645
    }
    
    if(defined $start && defined $end && $start > $end) {
646
      throw "Cannot request a slice whose start is greater than its end. Start: $start. End: $end" unless $no_errors;
647 648
    }
  }
649 650
  
  return ($seq_region_name, $start, $end, $strand);
651 652
}

653
=head2 fetch_by_region_unique
654

655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676
  Arg [1]    : string $coord_system_name (optional)
               The name of the coordinate system of the slice to be created
               This may be a name of an actual coordinate system or an alias
               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
               created on.
  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)
  Arg [7]    : boolean $no_fuzz (optional, default = undef (false))
               If true (non-zero), do not use "fuzzy matching" (see below).
  Example    : $slice = $slice_adaptor->fetch_by_region_unique('chromosome', 'HSCHR6_MHC_COX');
  Description: Retrieves a slice on the requested region but returns only the unique
               parts of the slice.  At a minimum the
               name the name of the seq_region to fetch must be provided.
677

678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711
               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.

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

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

               If the requested seq_region is not found in the database undef
               is returned.

  Returntype : listref Bio::EnsEMBL::Slice
  Exceptions : throw if no seq_region_name is provided
               throw if invalid coord_system_name is provided
               throw if start > end is provided
  Caller     : general
  Status     : Stable

=cut

sub fetch_by_region_unique {
  my $self = shift;

712
  my @out   = ();
713
  my $slice = $self->fetch_by_region(@_);
714

715

716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733
  if ( !exists( $self->{'asm_exc_cache'} ) ) {
    $self->_build_exception_cache();
  }

  if ( exists(
          $self->{'asm_exc_cache'}->{ $self->get_seq_region_id($slice) }
       ) )
  {
    # 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] );
734 735
      }
    }
736 737
  } else {
    @out = ($slice);
738
  }
739

740
  return \@out;
741
} ## end sub fetch_by_region_unique
742

743
=head2 fetch_by_name
744

Graham McVicker's avatar
Graham McVicker committed
745
  Arg [1]    : string $name
746 747 748 749 750 751 752 753 754
  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
755

756 757 758 759
               Returns undef if no seq_region with the provided name exists in
               the database.

  Returntype : Bio::EnsEMBL::Slice or undef
760 761
  Exceptions : throw if incorrent arg provided
  Caller     : Pipeline
762
  Status     : Stable
763 764 765

=cut

766 767 768
sub fetch_by_name {
  my $self = shift;
  my $name = shift;
Graham McVicker's avatar
Graham McVicker committed
769

770 771
  if(!$name) {
    throw("name argument is required");
Graham McVicker's avatar
Graham McVicker committed
772
  }
773 774 775

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

776
  if(scalar(@array) < 3 || scalar(@array) > 6) {
777
    throw("Malformed slice name [$name].  Format is " .
Laura Clarke's avatar
 
Laura Clarke committed
778
        "coord_system:version:name:start:end:strand");
779 780
  }

781
  # Rearrange arguments to suit fetch_by_region
782

783 784 785 786 787 788 789 790 791
  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
792
}
Graham McVicker's avatar
Graham McVicker committed
793 794 795



796 797 798 799 800
=head2 fetch_by_seq_region_id

  Arg [1]    : string $seq_region_id
               The internal identifier of the seq_region to create this slice
               on
801 802 803
  Arg [2]    : optional start
  Arg [3]    : optional end
  Arg [4]    : optional strand
804 805 806
  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.
807 808
               Returns undef if no such slice exists.
  Returntype : Bio::EnsEMBL::Slice or undef
Graham McVicker's avatar
Graham McVicker committed
809 810
  Exceptions : none
  Caller     : general
811
  Status     : Stable
Graham McVicker's avatar
Graham McVicker committed
812 813 814

=cut

815
sub fetch_by_seq_region_id {
816
  my ( $self, $seq_region_id, $start, $end, $strand, $check_prior_ids ) = @_;
817

Magali Ruffier's avatar
Magali Ruffier committed
818
  my $csa = $self->db->get_CoordSystemAdaptor();
819 820
  my $arr = $self->{'sr_id_cache'}->{$seq_region_id};
  my ( $name, $length, $cs, $cs_id );
821

822

823 824
  if ( $arr && defined( $arr->[2] ) ) {
    ( $name, $cs_id, $length ) = ( $arr->[1], $arr->[2], $arr->[3] );
825
    $cs = $self->db->get_CoordSystemAdaptor->fetch_by_dbID($cs_id);
826
  } else {
827 828
    my $sth =
      $self->prepare(   "SELECT sr.name, sr.coord_system_id, sr.length "
829
                      . "FROM seq_region sr "
830
                      . "WHERE sr.seq_region_id = ? " );
831

832
    $sth->bind_param( 1, $seq_region_id, SQL_INTEGER );
833
    $sth->execute();
834

835 836 837 838 839 840 841 842 843 844 845 846
    if ( $sth->rows() == 0 ) { 
      # This could have been an old seq region id so see if we can
      # translate it into a more recent version.
      if($check_prior_ids) {
        if(exists $csa->{_external_seq_region_mapping}->{$seq_region_id}) {
          my $new_seq_region_id = $csa->{_external_seq_region_mapping}->{$seq_region_id};
          # No need to pass check prior ids flag because it's a 1 step relationship
          return $self->fetch_by_seq_region_id($new_seq_region_id, $start, $end, $strand);
        }
      }
      return undef;
    }
847

848
    ( $name, $cs_id, $length ) = $sth->fetchrow_array();
849 850 851 852 853
    $sth->finish();

    $cs = $self->db->get_CoordSystemAdaptor->fetch_by_dbID($cs_id);

    #cache results to speed up repeated queries
854
    my $arr = [ $seq_region_id, $name, $cs_id, $length ];
855 856 857

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

860
  return
861 862 863 864 865 866 867 868
    Bio::EnsEMBL::Slice->new_fast({ 
	                      'coord_system'     => $cs,
                              'seq_region_name'  => $name,
                              'seq_region_length'=> $length,
                              'start'            => $start || 1,
                              'end'              => $end || $length,
                              'strand'           => $strand || 1,
                              'adaptor'           => $self} );
869
} ## end sub fetch_by_seq_region_id
Arne Stabenau's avatar
Arne Stabenau committed
870

Graham McVicker's avatar
Graham McVicker committed
871

872 873

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

Arne Stabenau's avatar
Arne Stabenau committed
875
  Arg [1]    : Bio::EnsEMBL::Slice $slice
876 877 878 879 880 881 882 883 884 885
               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
886
  Status     : Stable
Arne Stabenau's avatar
Arne Stabenau committed
887 888

=cut
Graham McVicker's avatar
Graham McVicker committed
889

890 891 892
sub get_seq_region_id {
  my $self = shift;
  my $slice = shift;
Graham McVicker's avatar
Graham McVicker committed
893

Ian Longden's avatar
Ian Longden committed
894
  if(!$slice || !ref($slice) || !($slice->isa('Bio::EnsEMBL::Slice') or $slice->isa('Bio::EnsEMBL::LRGSlice'))) {