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

Andy Yates's avatar
Andy Yates committed
3
  Copyright (c) 1999-2012 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);
Andy Yates's avatar
Andy Yates committed
104
use Scalar::Util qw/looks_like_number/;
105

106
@ISA = ('Bio::EnsEMBL::DBSQL::BaseAdaptor');
107

108 109
sub new {
  my $caller = shift;
110

111
  my $class = ref($caller) || $caller;
112

113
  my $self = $class->SUPER::new(@_);
114

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

118 119 120 121
  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'};
122

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


=head2 fetch_by_region

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

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

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

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

=cut

189 190 191 192

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

198 199
  if ( !defined($start) )  { $start  = 1 }
  if ( !defined($strand) ) { $strand = 1 }
200

201 202 203
  if ( !defined($seq_region_name) ) {
    throw('seq_region_name argument is required');
  }
204

205 206
  my $cs;
  my $csa = $self->db->get_CoordSystemAdaptor();
207

208 209
  if ( defined($coord_system_name) ) {
    $cs = $csa->fetch_by_name( $coord_system_name, $version );
210

211
    ## REMOVE THESE THREE LINES WHEN STICKLEBACK DB IS FIXED!
212 213 214 215 216 217 218
    ## 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!

219 220
    if ( !defined($cs) && $coord_system_name eq 'chromosome' ) {
      $cs = $csa->fetch_by_name( 'group', $version );
221 222
    }

223
    if ( !defined($cs) ) {
224
      throw( sprintf( "Unknown coordinate system:\n"
225
                        . "name='%s' version='%s'\n",
226
                      $coord_system_name, $version ) );
227
    }
228

229
    # fetching by toplevel is same as fetching w/o name or version
230 231
    if ( $cs->is_top_level() ) {
      $cs      = undef;
232 233
      $version = undef;
    }
234

235
  } ## end if ( defined($coord_system_name...))
236 237 238

  my $constraint;
  my $sql;
239
  my @bind_params;
240 241
  my $key;

242 243 244 245
  if ( defined($cs) ) {
    $sql = sprintf( "SELECT sr.name, sr.seq_region_id, sr.length, %d "
                      . "FROM seq_region sr ",
                    $cs->dbID() );
246

247 248
    $constraint = "AND sr.coord_system_id = ?";
    push( @bind_params, [ $cs->dbID(), SQL_INTEGER ] );
249

250
    $key = "$seq_region_name:" . $cs->dbID();
251
  } else {
252 253 254
    $sql =
      "SELECT sr.name, sr.seq_region_id, sr.length, cs.coord_system_id "
      . "FROM seq_region sr, coord_system cs ";
255

256 257 258 259 260
    $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) ) {
261
      $constraint .= "AND cs.version = ? ";
262
      push( @bind_params, [ $version, SQL_VARCHAR ] );
263
    }
264

265
    $constraint .= "ORDER BY cs.rank ASC";
266 267
  }

268
  # check the cache so we only go to the db if necessary
269
  my $length;
270
  my $arr;
271

272 273 274
  if ( defined($key) ) { $arr = $self->{'sr_name_cache'}->{$key} }

  if ( defined($arr) ) {
275
    $length = $arr->[3];
276
  } else {
277 278 279 280 281 282 283 284
    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] );
285
    }
286 287

    $sth->execute();
288

289
    if ( $sth->rows() == 0 ) {
290
      $sth->finish();
Graham McVicker's avatar
Graham McVicker committed
291

292 293 294 295 296 297 298 299 300 301 302 303 304 305

      # try synonyms
      my $syn_sql_sth = $self->prepare("select s.name from seq_region s, seq_region_synonym ss where s.seq_region_id = ss.seq_region_id and ss.synonym = ?");
      $syn_sql_sth->bind_param(1, $seq_region_name, SQL_VARCHAR);
      $syn_sql_sth->execute();
      my $new_name;
      $syn_sql_sth->bind_columns( \$new_name);
      if($syn_sql_sth->fetch){
	$syn_sql_sth->finish;
	return $self->fetch_by_region($coord_system_name, $new_name, $start, $end, $strand, $version, $no_fuzz);
      }
      $syn_sql_sth->finish;


306 307
      if ($no_fuzz) { return undef }

308 309
      # Do fuzzy matching, assuming that we are just missing a version
      # on the end of the seq_region name.
310

311 312
      $sth =
        $self->prepare( $sql . " WHERE sr.name LIKE ? " . $constraint );
313

314 315 316 317 318 319
      $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] );
320
      }
321

322
      $sth->execute();
323 324

      my $prefix_len = length($seq_region_name) + 1;
325 326
      my $high_ver   = undef;
      my $high_cs    = $cs;
327

328 329
      # Find the fuzzy-matched seq_region with the highest postfix
      # (which ought to be a version).
330

331 332
      my ( $tmp_name, $id, $tmp_length, $cs_id );
      $sth->bind_columns( \( $tmp_name, $id, $tmp_length, $cs_id ) );
333

Felix Kokocinski's avatar
Felix Kokocinski committed
334
      my $i = 0;
335

336 337 338
      while ( $sth->fetch ) {
        my $tmp_cs =
          ( defined($cs) ? $cs : $csa->fetch_by_dbID($cs_id) );
339

340
        # cache values for future reference
341
        my $arr = [ $id, $tmp_name, $cs_id, $tmp_length ];
342
        $self->{'sr_name_cache'}->{"$tmp_name:$cs_id"} = $arr;
343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361
        $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
362
        }
363 364

        $i++;
365
      } ## end while ( $sth->fetch )
366 367
      $sth->finish();

368
      # warn if fuzzy matching found more than one result
369 370 371 372 373 374 375
      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",
376
            $high_cs->name(), $seq_region_name ) );
377 378
      }

379 380
      $cs = $high_cs;

381 382
      # return if we did not find any appropriate match:
      if ( !defined($high_ver) ) { return undef }
383

384 385
    } else {

386 387 388
      my ( $id, $cs_id );
      ( $seq_region_name, $id, $length, $cs_id ) =
        $sth->fetchrow_array();
389
      $sth->finish();
390

391
      # cache to speed up for future queries
392
      my $arr = [ $id, $seq_region_name, $cs_id, $length ];
393
      $self->{'sr_name_cache'}->{"$seq_region_name:$cs_id"} = $arr;
394 395
      $self->{'sr_id_cache'}->{"$id"}                       = $arr;
      $cs = $csa->fetch_by_dbID($cs_id);
396
    }
397
  } ## end else [ if ( defined($arr) ) ]
398

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

401 402
  #If this was given then check if we've got a circular seq region otherwise
  #let it fall through to the normal Slice method
403
  if ( $end + 1 < $start ) {
404
    my $cs_id = $cs->dbID();
405
    my $seq_region_id = $self->{'sr_name_cache'}->{"$seq_region_name:$cs_id"}->[0];
406 407 408 409 410 411 412 413 414 415 416 417 418
    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;
    }
419 420
  }

421 422 423
  if ( defined( $self->{'lrg_region_test'} )
       and substr( $cs->name, 0, 3 ) eq $self->{'lrg_region_test'} )
  {
Ian Longden's avatar
Ian Longden committed
424
    return
425 426 427 428 429 430 431 432
      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
433
    return
434 435 436 437 438 439 440 441
      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
442
  }
443
} ## end sub fetch_by_region
444

445 446 447
=head2 fetch_by_toplevel_location

  Arg [1]     : string $location
448
                Ensembl formatted location. Can be a format like 
Andy Yates's avatar
Andy Yates committed
449 450 451 452 453 454
                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<->.
455 456
  Arg[2]      : boolean $no_warnings
                Suppress warnings from this method
457
  Example     : my $slice = $sa->fetch_by_toplevel_location('X:1-10000')
Andy Yates's avatar
Andy Yates committed
458
                my $slice = $sa->fetch_by_toplevel_location('X:1-10000:-1')
459 460 461
  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
462 463
                coordinate system. The code assumes that location formatting
                is not perfect and will perform basic cleanup before parsing.
464 465 466 467 468 469 470 471
  Returntype  : Bio::EnsEMBL::Slice
  Exceptions  : If $location is false otherwise see C<fetch_by_region()>
  Caller      : General
  Status      : Beta

=cut

sub fetch_by_toplevel_location {
472
  my ($self, $location, $no_warnings) = @_;
473

474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512
  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";
  }
  
  my $coord_system_name = 'toplevel';
  my $slice = $self->fetch_by_region($coord_system_name, $seq_region_name, $start, $end, $strand, undef, 0);
  return unless $slice;
  
  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
                separated by C<..>, C<:> or C<->.
  Arg[2]      : boolean $no_warnings
                Suppress warnings from this method
513 514
  Arg[3]      : boolean $no_errors
                Supress errors being thrown from this method
515 516 517 518 519 520 521 522 523
  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 {
524
  my ($self, $location, $no_warnings, $no_errors) = @_;
525 526 527
  
  throw 'You must specify a location' if ! $location;
  
528
  #cleanup any nomenclature like 1_000 or 1 000 or 1,000
529
  my $number_seps_regex = qr/\s+|,|_/;
530 531 532
  my $separator_regex = qr/(?:-|[.]{2}|\:)?/;
  my $number_regex = qr/[0-9,_ E]+/xms;
  my $strand_regex = qr/[+-1]|-1/xms;
533
  
534
  my $regex = qr/^((?:\w|\.|_|-)+) \s* :? \s* ($number_regex)? $separator_regex ($number_regex)? $separator_regex ($strand_regex)? $/xms;
535 536
  my ($seq_region_name, $start, $end, $strand);
  if(($seq_region_name, $start, $end, $strand) = $location =~ $regex) {
Andy Yates's avatar
Andy Yates committed
537 538 539 540 541 542 543
    
    if(defined $strand) {
      if(!looks_like_number($strand)) {
        $strand = ($strand eq '+') ? 1 : -1;
      }
    }
    
544
    if(defined $start) {
545
      $start =~ s/$number_seps_regex//g; 
546 547 548 549
      if($start < 1) {
        warning "Start was less than 1 (${start}) which is not allowed. Resetting to 1"  if ! $no_warnings;
        $start = 1;
      }
550
    }
551
    if(defined $end) {
552
      $end =~ s/$number_seps_regex//g;
553
      if($end < 1) {
554
        throw "Cannot request negative or 0 end indexes through this interface. Given $end but expected something greater than 0" unless $no_errors;
555
      }
556 557 558
    }
    
    if(defined $start && defined $end && $start > $end) {
559
      throw "Cannot request a slice whose start is greater than its end. Start: $start. End: $end" unless $no_errors;
560 561
    }
  }
562 563
  
  return ($seq_region_name, $start, $end, $strand);
564 565
}

566
=head2 fetch_by_region_unique
567

568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589
  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.
590

591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624
               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;

625
  my @out   = ();
626
  my $slice = $self->fetch_by_region(@_);
627

628

629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646
  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] );
647 648
      }
    }
649 650
  } else {
    @out = ($slice);
651
  }
652

653
  return \@out;
654
} ## end sub fetch_by_region_unique
655

656
=head2 fetch_by_name
657

Graham McVicker's avatar
Graham McVicker committed
658
  Arg [1]    : string $name
659 660 661 662 663 664 665 666 667
  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
668

669 670 671 672
               Returns undef if no seq_region with the provided name exists in
               the database.

  Returntype : Bio::EnsEMBL::Slice or undef
673 674
  Exceptions : throw if incorrent arg provided
  Caller     : Pipeline
675
  Status     : Stable
676 677 678

=cut

679 680 681
sub fetch_by_name {
  my $self = shift;
  my $name = shift;
Graham McVicker's avatar
Graham McVicker committed
682

683 684
  if(!$name) {
    throw("name argument is required");
Graham McVicker's avatar
Graham McVicker committed
685
  }
686 687 688

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

689
  if(scalar(@array) < 3 || scalar(@array) > 6) {
690
    throw("Malformed slice name [$name].  Format is " .
Laura Clarke's avatar
 
Laura Clarke committed
691
        "coord_system:version:name:start:end:strand");
692 693
  }

694
  # Rearrange arguments to suit fetch_by_region
695

696 697 698 699 700 701 702 703 704
  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
705
}
Graham McVicker's avatar
Graham McVicker committed
706 707 708



709 710 711 712 713
=head2 fetch_by_seq_region_id

  Arg [1]    : string $seq_region_id
               The internal identifier of the seq_region to create this slice
               on
714 715 716
  Arg [2]    : optional start
  Arg [3]    : optional end
  Arg [4]    : optional strand
717 718 719
  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.
720 721
               Returns undef if no such slice exists.
  Returntype : Bio::EnsEMBL::Slice or undef
Graham McVicker's avatar
Graham McVicker committed
722 723
  Exceptions : none
  Caller     : general
724
  Status     : Stable
Graham McVicker's avatar
Graham McVicker committed
725 726 727

=cut

728
sub fetch_by_seq_region_id {
729
  my ( $self, $seq_region_id, $start, $end, $strand ) = @_;
730

731 732
  my $arr = $self->{'sr_id_cache'}->{$seq_region_id};
  my ( $name, $length, $cs, $cs_id );
733

734

735 736
  if ( $arr && defined( $arr->[2] ) ) {
    ( $name, $cs_id, $length ) = ( $arr->[1], $arr->[2], $arr->[3] );
737
    $cs = $self->db->get_CoordSystemAdaptor->fetch_by_dbID($cs_id);
738
  } else {
739 740
    my $sth =
      $self->prepare(   "SELECT sr.name, sr.coord_system_id, sr.length "
741
                      . "FROM seq_region sr "
742
                      . "WHERE sr.seq_region_id = ? " );
743

744
    $sth->bind_param( 1, $seq_region_id, SQL_INTEGER );
745
    $sth->execute();
746

747
    if ( $sth->rows() == 0 ) { return undef }
748

749
    ( $name, $cs_id, $length ) = $sth->fetchrow_array();
750 751 752 753 754
    $sth->finish();

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

    #cache results to speed up repeated queries
755
    my $arr = [ $seq_region_id, $name, $cs_id, $length ];
756 757 758

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

761
  return
762 763 764 765 766 767 768 769
    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} );
770
} ## end sub fetch_by_seq_region_id
Arne Stabenau's avatar
Arne Stabenau committed
771

Graham McVicker's avatar
Graham McVicker committed
772

773 774

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

Arne Stabenau's avatar
Arne Stabenau committed
776
  Arg [1]    : Bio::EnsEMBL::Slice $slice
777 778 779 780 781 782 783 784 785 786
               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
787
  Status     : Stable
Arne Stabenau's avatar
Arne Stabenau committed
788 789

=cut
Graham McVicker's avatar
Graham McVicker committed
790

791 792 793
sub get_seq_region_id {
  my $self = shift;
  my $slice = shift;
Graham McVicker's avatar
Graham McVicker committed
794

Ian Longden's avatar
Ian Longden committed
795
  if(!$slice || !ref($slice) || !($slice->isa('Bio::EnsEMBL::Slice') or $slice->isa('Bio::EnsEMBL::LRGSlice'))) {
796 797
    throw('Slice argument is required');
  }
798
  
799
  my $seq_region_name = $slice->seq_region_name();
800
  my $key = $seq_region_name.":".$slice->coord_system->dbID();
801
  my $arr = $self->{'sr_name_cache'}->{"$key"};
802

803 804
  if( $arr ) {
    return $arr->[0];
805 806
  }

807
  my $cs_id = $slice->coord_system->dbID();
808 809 810 811 812 813

  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
814 815 816
  $sth->bind_param(1,"$seq_region_name",SQL_VARCHAR);
  $sth->bind_param(2,$cs_id,SQL_INTEGER);
  $sth->execute();
817 818 819

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

823 824 825 826 827 828
  }

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

  #cache information for future requests
829
  $arr = [ $seq_region_id, $seq_region_name, $cs_id, $length ];
830 831 832 833

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

834
  return $seq_region_id;
Arne Stabenau's avatar
Arne Stabenau committed
835
}
836 837


Arne Stabenau's avatar
Arne Stabenau committed
838

839 840 841 842 843 844 845 846 847
=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
848
  Arg [3]    : bool $include_non_reference (optional)
849 850 851
               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
852 853
  Arg [4]    : int $include_duplicates (optional)
               If set duplicate regions will be returned.
854
               
855
               NOTE: if you do not use this option and you have a PAR
856 857 858
               (pseudo-autosomal region) at the beginning of your seq_region
               then your slice will not start at position 1, so coordinates
               retrieved from this slice might not be what you expected.
859

860 861 862 863
  Arg[5]     : bool $include_lrg (optional)  (default 0)
               If set lrg regions will be returned aswell.


864 865 866
  Example    : @chromos = @{$slice_adaptor->fetch_all('chromosome','NCBI33')};
               @contigs = @{$slice_adaptor->fetch_all('contig')};

Graham McVicker's avatar
Graham McVicker committed
867 868 869 870 871
               # 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)};
872 873 874 875 876 877

  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.
878 879
               If the coordinate system with the provided name and version
               does not exist an empty list is returned.
880
               If the coordinate system name provided is 'toplevel', all
881 882
               non-redundant toplevel slices are returned (note that any
               coord_system_version argument is ignored in that case).
883 884 885 886

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

887
  Returntype : listref of Bio::EnsEMBL::Slices
888
  Exceptions : none
889
  Caller     : general
890
  Status     : Stable
Web Admin's avatar
Web Admin committed
891 892 893