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

3
  Copyright (c) 1999-2013 The European Bioinformatics Institute and
4
  Genome Research Limited.  All rights reserved.
5

6 7 8 9 10 11 12 13
  This software is distributed under a modified Apache license.
  For license details, please see

    http://www.ensembl.org/info/about/code_licence.html

=head1 CONTACT

  Please email comments or questions to the public Ensembl
14
  developers list at <dev@ensembl.org>.
15 16 17 18 19

  Questions may also be sent to the Ensembl help desk at
  <helpdesk@ensembl.org>.

=cut
20 21 22

=head1 NAME

23 24
Bio::EnsEMBL::DBSQL::SliceAdaptor - A database aware adaptor responsible for
the creation of Slice objects.
25 26 27

=head1 SYNOPSIS

Graham McVicker's avatar
Graham McVicker committed
28
  use Bio::EnsEMBL::Utils::Slice qw(split_Slices);
29 30 31
  use Bio::EnsEMBL::Registry;

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

36 37
  $slice_adaptor =
    Bio::EnsEMBL::Registry->get_adaptor( "human", "core", "slice" );
38

Graham McVicker's avatar
Graham McVicker committed
39
  # get a slice on the entire chromosome X
40
  $chr_slice = $slice_adaptor->fetch_by_region( 'chromosome', 'X' );
Graham McVicker's avatar
Graham McVicker committed
41 42

  # get a slice for each clone in the database
43 44
  foreach $cln_slice ( @{ $slice_adaptor->fetch_all('clone') } ) {
    # do something with clone
45 46
  }

Graham McVicker's avatar
Graham McVicker committed
47
  # get a slice which is part of NT_004321
48 49 50
  $spctg_slice =
    $slice_adaptor->fetch_by_region( 'supercontig', 'NT_004321',
    200_000, 600_000 );
51

52 53
  # get all non-redundant slices from the highest possible coordinate
  # systems
Graham McVicker's avatar
Graham McVicker committed
54 55 56
  $slices = $slice_adaptor->fetch_all('toplevel');

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

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

  # split up a list of slices into smaller slices
63
  $overlap    = 1000;
Graham McVicker's avatar
Graham McVicker committed
64
  $max_length = 1e6;
65
  $slices     = split_Slices( $slices, $max_length, $overlap );
Graham McVicker's avatar
Graham McVicker committed
66 67

  # store a list of slice names in a file
68
  open( FILE, ">$filename" ) or die("Could not open file $filename");
Graham McVicker's avatar
Graham McVicker committed
69 70 71 72 73 74
  foreach my $slice (@$slices) {
    print FILE $slice->name(), "\n";
  }
  close FILE;

  # retreive a list of slices from a file
75 76
  open( FILE, $filename ) or die("Could not open file $filename");
  while ( $name = <FILE> ) {
Graham McVicker's avatar
Graham McVicker committed
77 78 79 80 81
    chomp($name);
    $slice = $slice_adaptor->fetch_by_name($name);
    # do something with slice
  }

82 83
=head1 DESCRIPTION

84 85
This module is responsible for fetching Slices representing genomic
regions from a database.  A Details on how slices can be used are in the
86
Bio::EnsEMBL::Slice module.
87

88
=head1 METHODS
89 90 91 92

=cut


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

97

98 99
use Bio::EnsEMBL::DBSQL::BaseAdaptor;
use Bio::EnsEMBL::Slice;
100
use Bio::EnsEMBL::CircularSlice;
101
use Bio::EnsEMBL::Mapper;
Ian Longden's avatar
Ian Longden committed
102
use Bio::EnsEMBL::LRGSlice;
103
use Bio::EnsEMBL::Utils::Exception qw(throw deprecate warning stack_trace_dump);
104
use Bio::EnsEMBL::ProjectionSegment;
Andy Yates's avatar
Andy Yates committed
105
use Scalar::Util qw/looks_like_number/;
106
use Bio::EnsEMBL::Utils::Scalar qw/assert_integer/;
107

108
@ISA = ('Bio::EnsEMBL::DBSQL::BaseAdaptor');
109

110 111
sub new {
  my $caller = shift;
112

113
  my $class = ref($caller) || $caller;
114

115
  my $self = $class->SUPER::new(@_);
116

117 118
  # use a cache which is shared and also used by the assembly
  # mapper adaptor
119

120 121 122 123
  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'};
124

Ian Longden's avatar
Ian Longden committed
125 126 127
  $self->{'lrg_region_test'} = undef;
  my $meta_container = $self->db->get_MetaContainer();
  my @values = $meta_container->list_value_by_key("LRG");
128
  if(scalar(@values) and $values[0]->[0]){
Ian Longden's avatar
Ian Longden committed
129 130
    $self->{'lrg_region_test'} = $values[0]->[0];
  }
131 132 133 134 135 136
  return $self;
}


=head2 fetch_by_region

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

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

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

179 180 181 182 183 184 185
               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
186
  Caller     : general
187
  Status     : Stable
188 189 190

=cut

191 192 193 194

#
# ARNE: This subroutine needs simplification!! 
#
195
sub fetch_by_region {
196 197 198
  my ( $self, $coord_system_name, $seq_region_name, $start, $end,
       $strand, $version, $no_fuzz )
    = @_;
199

200 201 202
  assert_integer($start, 'start') if $start;
  assert_integer($end, 'end') if $end;

203 204
  if ( !defined($start) )  { $start  = 1 }
  if ( !defined($strand) ) { $strand = 1 }
205

206 207 208
  if ( !defined($seq_region_name) ) {
    throw('seq_region_name argument is required');
  }
209

210 211
  my $cs;
  my $csa = $self->db->get_CoordSystemAdaptor();
212

213 214
  if ( defined($coord_system_name) ) {
    $cs = $csa->fetch_by_name( $coord_system_name, $version );
215

216
    ## REMOVE THESE THREE LINES WHEN STICKLEBACK DB IS FIXED!
217 218 219 220 221 222 223
    ## 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!

224 225
    if ( !defined($cs) && $coord_system_name eq 'chromosome' ) {
      $cs = $csa->fetch_by_name( 'group', $version );
226 227
    }

228
    if ( !defined($cs) ) {
229
      throw( sprintf( "Unknown coordinate system:\n"
230
                        . "name='%s' version='%s'\n",
231
                      $coord_system_name, $version ) );
232
    }
233

234
    # fetching by toplevel is same as fetching w/o name or version
235 236
    if ( $cs->is_top_level() ) {
      $cs      = undef;
237 238
      $version = undef;
    }
239

240
  } ## end if ( defined($coord_system_name...))
241 242 243

  my $constraint;
  my $sql;
244
  my @bind_params;
245 246
  my $key;

247 248 249 250
  if ( defined($cs) ) {
    $sql = sprintf( "SELECT sr.name, sr.seq_region_id, sr.length, %d "
                      . "FROM seq_region sr ",
                    $cs->dbID() );
251

252 253
    $constraint = "AND sr.coord_system_id = ?";
    push( @bind_params, [ $cs->dbID(), SQL_INTEGER ] );
254

255
    $key = "$seq_region_name:" . $cs->dbID();
256
  } else {
257 258 259
    $sql =
      "SELECT sr.name, sr.seq_region_id, sr.length, cs.coord_system_id "
      . "FROM seq_region sr, coord_system cs ";
260

261 262 263 264 265
    $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) ) {
266
      $constraint .= "AND cs.version = ? ";
267
      push( @bind_params, [ $version, SQL_VARCHAR ] );
268
    }
269

270
    $constraint .= "ORDER BY cs.rank ASC";
271 272
  }

273
  # check the cache so we only go to the db if necessary
274
  my $length;
275
  my $arr;
276

277 278 279
  if ( defined($key) ) { $arr = $self->{'sr_name_cache'}->{$key} }

  if ( defined($arr) ) {
280
    $length = $arr->[3];
281
  } else {
282 283 284 285 286 287 288 289
    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] );
290
    }
291 292

    $sth->execute();
293

294
    if ( $sth->rows() == 0 ) {
295
      $sth->finish();
Graham McVicker's avatar
Graham McVicker committed
296

297 298

      # try synonyms
299
      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 =?");
300
      $syn_sql_sth->bind_param(1, $seq_region_name, SQL_VARCHAR);
301
      $syn_sql_sth->bind_param(2, $self->species_id(), SQL_INTEGER);
302
      $syn_sql_sth->execute();
303 304
      my ($new_name, $new_coord_system, $new_version);
      $syn_sql_sth->bind_columns( \$new_name, \$new_coord_system, \$new_version);
305
            
306
      if($syn_sql_sth->fetch){
307
        $syn_sql_sth->finish;
308
        if ((not defined($cs)) || ($cs->name eq $new_coord_system)) {
309
            return $self->fetch_by_region($new_coord_system, $new_name, $start, $end, $strand, $new_version, $no_fuzz);
310
        } elsif ($cs->name ne $new_coord_system) {
311 312 313 314 315
            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;
        }
        
316 317 318 319
      }
      $syn_sql_sth->finish;


320 321
      if ($no_fuzz) { return undef }

322 323
      # Do fuzzy matching, assuming that we are just missing a version
      # on the end of the seq_region name.
324

325 326
      $sth =
        $self->prepare( $sql . " WHERE sr.name LIKE ? " . $constraint );
327

328 329 330 331 332 333
      $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] );
334
      }
335

336
      $sth->execute();
337 338

      my $prefix_len = length($seq_region_name) + 1;
339 340
      my $high_ver   = undef;
      my $high_cs    = $cs;
341

342 343
      # Find the fuzzy-matched seq_region with the highest postfix
      # (which ought to be a version).
344

345 346
      my ( $tmp_name, $id, $tmp_length, $cs_id );
      $sth->bind_columns( \( $tmp_name, $id, $tmp_length, $cs_id ) );
347

Felix Kokocinski's avatar
Felix Kokocinski committed
348
      my $i = 0;
349

350 351 352
      while ( $sth->fetch ) {
        my $tmp_cs =
          ( defined($cs) ? $cs : $csa->fetch_by_dbID($cs_id) );
353

354
        # cache values for future reference
355
        my $arr = [ $id, $tmp_name, $cs_id, $tmp_length ];
356
        $self->{'sr_name_cache'}->{"$tmp_name:$cs_id"} = $arr;
357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375
        $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
376
        }
377 378

        $i++;
379
      } ## end while ( $sth->fetch )
380 381
      $sth->finish();

382
      # warn if fuzzy matching found more than one result
383 384 385 386 387 388 389
      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",
390
            $high_cs->name(), $seq_region_name ) );
391 392
      }

393 394
      $cs = $high_cs;

395 396
      # return if we did not find any appropriate match:
      if ( !defined($high_ver) ) { return undef }
397

398 399
    } else {

400 401 402
      my ( $id, $cs_id );
      ( $seq_region_name, $id, $length, $cs_id ) =
        $sth->fetchrow_array();
403
      $sth->finish();
404

405
      # cache to speed up for future queries
406
      my $arr = [ $id, $seq_region_name, $cs_id, $length ];
407
      $self->{'sr_name_cache'}->{"$seq_region_name:$cs_id"} = $arr;
408 409
      $self->{'sr_id_cache'}->{"$id"}                       = $arr;
      $cs = $csa->fetch_by_dbID($cs_id);
410
    }
411
  } ## end else [ if ( defined($arr) ) ]
412

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

415 416
  #If this was given then check if we've got a circular seq region otherwise
  #let it fall through to the normal Slice method
417
  if ( $end + 1 < $start ) {
418
    my $cs_id = $cs->dbID();
419
    my $seq_region_id = $self->{'sr_name_cache'}->{"$seq_region_name:$cs_id"}->[0];
420 421 422 423 424 425 426 427 428 429 430 431 432
    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;
    }
433 434
  }

435 436 437
  if ( defined( $self->{'lrg_region_test'} )
       and substr( $cs->name, 0, 3 ) eq $self->{'lrg_region_test'} )
  {
Ian Longden's avatar
Ian Longden committed
438
    return
439 440 441 442 443 444 445 446
      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
447
    return
448 449 450 451 452 453 454 455
      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
456
  }
457
} ## end sub fetch_by_region
458

459 460 461
=head2 fetch_by_toplevel_location

  Arg [1]     : string $location
462
                Ensembl formatted location. Can be a format like 
Andy Yates's avatar
Andy Yates committed
463 464 465 466 467 468
                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<->.
469 470
  Arg[2]      : boolean $no_warnings
                Suppress warnings from this method
471 472
  Arg[3]      : boolean $no_fuzz
                Stop fuzzy matching of sequence regions from occuring
473 474 475 476
  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
477
  Example     : my $slice = $sa->fetch_by_toplevel_location('X:1-10000')
Andy Yates's avatar
Andy Yates committed
478
                my $slice = $sa->fetch_by_toplevel_location('X:1-10000:-1')
479 480 481
  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
482 483
                coordinate system. The code assumes that location formatting
                is not perfect and will perform basic cleanup before parsing.
484
  Returntype  : Bio::EnsEMBL::Slice
485 486
  Exceptions  : If $location is false otherwise see C<fetch_by_location()>
                or C<fetch_by_region()>
487 488 489 490 491 492
  Caller      : General
  Status      : Beta

=cut

sub fetch_by_toplevel_location {
493 494
  my ($self, $location, $no_warnings, $no_fuzz, $ucsc) = @_;
  return $self->fetch_by_location($location, 'toplevel', undef, $no_warnings, $no_fuzz, $ucsc);
495 496 497 498 499 500 501 502 503 504 505
}

=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
506
                separated by C<..>, C<:>, C<_> or C<->.
507 508 509 510 511 512 513 514
  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
515 516 517 518
  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
519 520 521 522 523 524 525 526 527 528 529 530 531 532
  Example     : my $slice = $sa->fetch_by_toplevel_location('X:1-10000')
                my $slice = $sa->fetch_by_toplevel_location('X:1-10000:-1')
  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
                coordinate system. The code assumes that location formatting
                is not perfect and will perform basic cleanup before parsing.
  Returntype  : Bio::EnsEMBL::Slice
  Exceptions  : If $location or coordinate system is false otherwise 
                see C<fetch_by_region()>
  Caller      : General
  Status      : Beta

=cut
533

534
sub fetch_by_location {
535
  my ($self, $location, $coord_system_name, $coord_system_version, $no_warnings, $no_fuzz, $ucsc) = @_;
536 537 538
  
  throw "No coordinate system name specified" unless $coord_system_name;
  
539 540 541 542 543 544 545 546 547 548
  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";
  }
  
549
  my $slice = $self->fetch_by_region($coord_system_name, $seq_region_name, $start, $end, $strand, $coord_system_version, $no_fuzz);
550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565
  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
    }
  }
566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588
  
  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
589
                separated by C<..>, C<:> C<_>, or C<->.
590 591
  Arg[2]      : boolean $no_warnings
                Suppress warnings from this method
592 593
  Arg[3]      : boolean $no_errors
                Supress errors being thrown from this method
594 595 596 597 598 599 600 601 602
  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 {
603
  my ($self, $location, $no_warnings, $no_errors) = @_;
604 605 606
  
  throw 'You must specify a location' if ! $location;
  
607 608 609 610
  #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;
611
  my $strand_regex = qr/[+-1]|-1/xms;
612
  
613
  my $regex = qr/^((?:\w|\.|_|-)+) \s* :? \s* ($number_regex)? $separator_regex ($number_regex)? $separator_regex ($strand_regex)? $/xms;
614 615
  my ($seq_region_name, $start, $end, $strand);
  if(($seq_region_name, $start, $end, $strand) = $location =~ $regex) {
Andy Yates's avatar
Andy Yates committed
616 617 618 619 620 621 622
    
    if(defined $strand) {
      if(!looks_like_number($strand)) {
        $strand = ($strand eq '+') ? 1 : -1;
      }
    }
    
623
    if(defined $start) {
624
      $start =~ s/$number_seps_regex//g; 
625 626 627 628
      if($start < 1) {
        warning "Start was less than 1 (${start}) which is not allowed. Resetting to 1"  if ! $no_warnings;
        $start = 1;
      }
629
    }
630
    if(defined $end) {
631
      $end =~ s/$number_seps_regex//g;
632
      if($end < 1) {
633
        throw "Cannot request negative or 0 end indexes through this interface. Given $end but expected something greater than 0" unless $no_errors;
634
      }
635 636 637
    }
    
    if(defined $start && defined $end && $start > $end) {
638
      throw "Cannot request a slice whose start is greater than its end. Start: $start. End: $end" unless $no_errors;
639 640
    }
  }
641 642
  
  return ($seq_region_name, $start, $end, $strand);
643 644
}

645
=head2 fetch_by_region_unique
646

647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668
  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.
669

670 671 672 673 674 675 676 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
               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;

704
  my @out   = ();
705
  my $slice = $self->fetch_by_region(@_);
706

707

708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725
  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] );
726 727
      }
    }
728 729
  } else {
    @out = ($slice);
730
  }
731

732
  return \@out;
733
} ## end sub fetch_by_region_unique
734

735
=head2 fetch_by_name
736

Graham McVicker's avatar
Graham McVicker committed
737
  Arg [1]    : string $name
738 739 740 741 742 743 744 745 746
  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
747

748 749 750 751
               Returns undef if no seq_region with the provided name exists in
               the database.

  Returntype : Bio::EnsEMBL::Slice or undef
752 753
  Exceptions : throw if incorrent arg provided
  Caller     : Pipeline
754
  Status     : Stable
755 756 757

=cut

758 759 760
sub fetch_by_name {
  my $self = shift;
  my $name = shift;
Graham McVicker's avatar
Graham McVicker committed
761

762 763
  if(!$name) {
    throw("name argument is required");
Graham McVicker's avatar
Graham McVicker committed
764
  }
765 766 767

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

768
  if(scalar(@array) < 3 || scalar(@array) > 6) {
769
    throw("Malformed slice name [$name].  Format is " .
Laura Clarke's avatar
 
Laura Clarke committed
770
        "coord_system:version:name:start:end:strand");
771 772
  }

773
  # Rearrange arguments to suit fetch_by_region
774

775 776 777 778 779 780 781 782 783
  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
784
}
Graham McVicker's avatar
Graham McVicker committed
785 786 787



788 789 790 791 792
=head2 fetch_by_seq_region_id

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

=cut

807
sub fetch_by_seq_region_id {
808
  my ( $self, $seq_region_id, $start, $end, $strand ) = @_;
809

810 811
  my $arr = $self->{'sr_id_cache'}->{$seq_region_id};
  my ( $name, $length, $cs, $cs_id );
812

813

814 815
  if ( $arr && defined( $arr->[2] ) ) {
    ( $name, $cs_id, $length ) = ( $arr->[1], $arr->[2], $arr->[3] );
816
    $cs = $self->db->get_CoordSystemAdaptor->fetch_by_dbID($cs_id);
817
  } else {
818 819
    my $sth =
      $self->prepare(   "SELECT sr.name, sr.coord_system_id, sr.length "
820
                      . "FROM seq_region sr "
821
                      . "WHERE sr.seq_region_id = ? " );
822

823
    $sth->bind_param( 1, $seq_region_id, SQL_INTEGER );
824
    $sth->execute();
825

826
    if ( $sth->rows() == 0 ) { return undef }
827

828
    ( $name, $cs_id, $length ) = $sth->fetchrow_array();
829 830 831 832 833
    $sth->finish();

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

    #cache results to speed up repeated queries
834
    my $arr = [ $seq_region_id, $name, $cs_id, $length ];
835 836 837

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

840
  return
841 842 843 844 845 846 847 848
    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} );
849
} ## end sub fetch_by_seq_region_id
Arne Stabenau's avatar
Arne Stabenau committed
850

Graham McVicker's avatar
Graham McVicker committed
851

852 853

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

Arne Stabenau's avatar
Arne Stabenau committed
855
  Arg [1]    : Bio::EnsEMBL::Slice $slice
856 857 858 859 860 861 862 863 864 865
               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
866
  Status     : Stable
Arne Stabenau's avatar
Arne Stabenau committed
867 868

=cut
Graham McVicker's avatar
Graham McVicker committed
869

870 871 872
sub get_seq_region_id {
  my $self = shift;
  my $slice = shift;
Graham McVicker's avatar
Graham McVicker committed
873

Ian Longden's avatar
Ian Longden committed
874
  if(!$slice || !ref($slice) || !($slice->isa('Bio::EnsEMBL::Slice') or $slice->isa('Bio::EnsEMBL::LRGSlice'))) {
875 876
    throw('Slice argument is required');
  }
877
  
878
  my $seq_region_name = $slice->seq_region_name();
879
  my $key = $seq_region_name.":".$slice->coord_system->dbID();
880
  my $arr = $self->{'sr_name_cache'}->{"$key"};
881

882 883
  if( $arr ) {
    return $arr->[0];
884 885
  }