Mercurial > repos > mahtabm > ensembl

comparison variant_effect_predictor/Bio/EnsEMBL/DBSQL/SliceAdaptor.pm @ 0:1f6dce3d34e0

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author mahtabm
date Thu, 11 Apr 2013 02:01:53 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:1f6dce3d34e0
1 =head1 LICENSE
2
3 Copyright (c) 1999-2012 The European Bioinformatics Institute and
4 Genome Research Limited. All rights reserved.
5
6 This software is distributed under a modified Apache license.
7 For license details, please see
8
9 http://www.ensembl.org/info/about/code_licence.html
10
11 =head1 CONTACT
12
13 Please email comments or questions to the public Ensembl
14 developers list at <dev@ensembl.org>.
15
16 Questions may also be sent to the Ensembl help desk at
17 <helpdesk@ensembl.org>.
18
19 =cut
20
21 =head1 NAME
22
23 Bio::EnsEMBL::DBSQL::SliceAdaptor - A database aware adaptor responsible for
24 the creation of Slice objects.
25
26 =head1 SYNOPSIS
27
28 use Bio::EnsEMBL::Utils::Slice qw(split_Slices);
29 use Bio::EnsEMBL::Registry;
30
31 Bio::EnsEMBL::Registry->load_registry_from_db(
32 -host => 'ensembldb.ensembl.org',
33 -user => 'anonymous'
34 );
35
36 $slice_adaptor =
37 Bio::EnsEMBL::Registry->get_adaptor( "human", "core", "slice" );
38
39 # get a slice on the entire chromosome X
40 $chr_slice = $slice_adaptor->fetch_by_region( 'chromosome', 'X' );
41
42 # get a slice for each clone in the database
43 foreach $cln_slice ( @{ $slice_adaptor->fetch_all('clone') } ) {
44 # do something with clone
45 }
46
47 # get a slice which is part of NT_004321
48 $spctg_slice =
49 $slice_adaptor->fetch_by_region( 'supercontig', 'NT_004321',
50 200_000, 600_000 );
51
52 # get all non-redundant slices from the highest possible coordinate
53 # systems
54 $slices = $slice_adaptor->fetch_all('toplevel');
55
56 # include non-reference regions
57 $slices = $slice_adaptor->fetch_all( 'toplevel', undef, 1 );
58
59 # include non-duplicate regions
60 $slices = $slice_adaptor->fetch_all( 'toplevel', undef, 0, 1 );
61
62 # split up a list of slices into smaller slices
63 $overlap = 1000;
64 $max_length = 1e6;
65 $slices = split_Slices( $slices, $max_length, $overlap );
66
67 # store a list of slice names in a file
68 open( FILE, ">$filename" ) or die("Could not open file $filename");
69 foreach my $slice (@$slices) {
70 print FILE $slice->name(), "\n";
71 }
72 close FILE;
73
74 # retreive a list of slices from a file
75 open( FILE, $filename ) or die("Could not open file $filename");
76 while ( $name = <FILE> ) {
77 chomp($name);
78 $slice = $slice_adaptor->fetch_by_name($name);
79 # do something with slice
80 }
81
82 =head1 DESCRIPTION
83
84 This module is responsible for fetching Slices representing genomic
85 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 =cut
91
92
93 package Bio::EnsEMBL::DBSQL::SliceAdaptor;
94 use vars qw(@ISA);
95 use strict;
96
97
98 use Bio::EnsEMBL::DBSQL::BaseAdaptor;
99 use Bio::EnsEMBL::Slice;
100 use Bio::EnsEMBL::CircularSlice;
101 use Bio::EnsEMBL::Mapper;
102 use Bio::EnsEMBL::LRGSlice;
103 use Bio::EnsEMBL::Utils::Exception qw(throw deprecate warning stack_trace_dump);
104 use Scalar::Util qw/looks_like_number/;
105
106 @ISA = ('Bio::EnsEMBL::DBSQL::BaseAdaptor');
107
108 sub new {
109 my $caller = shift;
110
111 my $class = ref($caller) || $caller;
112
113 my $self = $class->SUPER::new(@_);
114
115 # use a cache which is shared and also used by the assembly
116 # mapper adaptor
117
118 my $seq_region_cache = $self->db->get_SeqRegionCache();
119
120 $self->{'sr_name_cache'} = $seq_region_cache->{'name_cache'};
121 $self->{'sr_id_cache'} = $seq_region_cache->{'id_cache'};
122
123 $self->{'lrg_region_test'} = undef;
124 my $meta_container = $self->db->get_MetaContainer();
125 my @values = $meta_container->list_value_by_key("LRG");
126 if(scalar(@values) and $values[0]->[0]){
127 $self->{'lrg_region_test'} = $values[0]->[0];
128 }
129 return $self;
130 }
131
132
133 =head2 fetch_by_region
134
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 to a coordinate system. Valid aliases are 'seqlevel' or
139 'toplevel'.
140 Arg [2] : string $seq_region_name
141 The name of the sequence region that the slice will be
142 created on.
143 Arg [3] : int $start (optional, default = 1)
144 The start of the slice on the sequence region
145 Arg [4] : int $end (optional, default = seq_region length)
146 The end of the slice on the sequence region
147 Arg [5] : int $strand (optional, default = 1)
148 The orientation of the slice on the sequence region
149 Arg [6] : string $version (optional, default = default version)
150 The version of the coordinate system to use (e.g. NCBI33)
151 Arg [7] : boolean $no_fuzz (optional, default = undef (false))
152 If true (non-zero), do not use "fuzzy matching" (see below).
153 Example : $slice = $slice_adaptor->fetch_by_region('chromosome', 'X');
154 $slice = $slice_adaptor->fetch_by_region('clone', 'AC008066.4');
155 Description: Retrieves a slice on the requested region. At a minimum the
156 name the name of the seq_region to fetch must be provided.
157
158 If no coordinate system name is provided than a slice on the
159 highest ranked coordinate system with a matching
160 seq_region_name will be returned. If a version but no
161 coordinate system name is provided, the same behaviour will
162 apply, but only coordinate systems of the appropriate version
163 are considered. The same applies if the 'toplevel' coordinate
164 system is specified, however in this case the version is
165 ignored. The coordinate system should always be specified if
166 it is known, since this is unambiguous and faster.
167
168 Some fuzzy matching is performed if no exact match for
169 the provided name is found. This allows clones to be
170 fetched even when their version is not known. For
171 example fetch_by_region('clone', 'AC008066') will
172 retrieve the sequence_region with name 'AC008066.4'.
173
174 The fuzzy matching can be turned off by setting the
175 $no_fuzz argument to a true value.
176
177 If the requested seq_region is not found in the database undef
178 is returned.
179
180 Returntype : Bio::EnsEMBL::Slice or undef
181 Exceptions : throw if no seq_region_name is provided
182 throw if invalid coord_system_name is provided
183 throw if start > end is provided
184 Caller : general
185 Status : Stable
186
187 =cut
188
189
190 #
191 # ARNE: This subroutine needs simplification!!
192 #
193 sub fetch_by_region {
194 my ( $self, $coord_system_name, $seq_region_name, $start, $end,
195 $strand, $version, $no_fuzz )
196 = @_;
197
198 if ( !defined($start) ) { $start = 1 }
199 if ( !defined($strand) ) { $strand = 1 }
200
201 if ( !defined($seq_region_name) ) {
202 throw('seq_region_name argument is required');
203 }
204
205 my $cs;
206 my $csa = $self->db->get_CoordSystemAdaptor();
207
208 if ( defined($coord_system_name) ) {
209 $cs = $csa->fetch_by_name( $coord_system_name, $version );
210
211 ## REMOVE THESE THREE LINES WHEN STICKLEBACK DB IS FIXED!
212 ## Anne/ap5 (2007-10-09):
213 # The problem was that the stickleback genebuild called the
214 # chromosomes 'groups', which meant they weren't being picked out by
215 # the karyotype drawing code. Apparently they are usually called
216 # 'groups' in the stickleback community, even though they really are
217 # chromosomes!
218
219 if ( !defined($cs) && $coord_system_name eq 'chromosome' ) {
220 $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 if ( $cs->is_top_level() ) {
231 $cs = undef;
232 $version = undef;
233 }
234
235 } ## end if ( defined($coord_system_name...))
236
237 my $constraint;
238 my $sql;
239 my @bind_params;
240 my $key;
241
242 if ( defined($cs) ) {
243 $sql = sprintf( "SELECT sr.name, sr.seq_region_id, sr.length, %d "
244 . "FROM seq_region sr ",
245 $cs->dbID() );
246
247 $constraint = "AND sr.coord_system_id = ?";
248 push( @bind_params, [ $cs->dbID(), SQL_INTEGER ] );
249
250 $key = "$seq_region_name:" . $cs->dbID();
251 } else {
252 $sql =
253 "SELECT sr.name, sr.seq_region_id, sr.length, cs.coord_system_id "
254 . "FROM seq_region sr, coord_system cs ";
255
256 $constraint = "AND sr.coord_system_id = cs.coord_system_id "
257 . "AND cs.species_id = ? ";
258 push( @bind_params, [ $self->species_id(), SQL_INTEGER ] );
259
260 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 if ( defined($key) ) { $arr = $self->{'sr_name_cache'}->{$key} }
273
274 if ( defined($arr) ) {
275 $length = $arr->[3];
276 } else {
277 my $sth =
278 $self->prepare( $sql . "WHERE sr.name = ? " . $constraint );
279
280 unshift( @bind_params, [ $seq_region_name, SQL_VARCHAR ] );
281
282 my $pos = 0;
283 foreach my $param (@bind_params) {
284 $sth->bind_param( ++$pos, $param->[0], $param->[1] );
285 }
286
287 $sth->execute();
288
289 if ( $sth->rows() == 0 ) {
290 $sth->finish();
291
292
293 # try synonyms
294 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 = ?");
295 $syn_sql_sth->bind_param(1, $seq_region_name, SQL_VARCHAR);
296 $syn_sql_sth->execute();
297 my $new_name;
298 $syn_sql_sth->bind_columns( \$new_name);
299 if($syn_sql_sth->fetch){
300 $syn_sql_sth->finish;
301 return $self->fetch_by_region($coord_system_name, $new_name, $start, $end, $strand, $version, $no_fuzz);
302 }
303 $syn_sql_sth->finish;
304
305
306 if ($no_fuzz) { return undef }
307
308 # Do fuzzy matching, assuming that we are just missing a version
309 # on the end of the seq_region name.
310
311 $sth =
312 $self->prepare( $sql . " WHERE sr.name LIKE ? " . $constraint );
313
314 $bind_params[0] =
315 [ sprintf( '%s.%%', $seq_region_name ), SQL_VARCHAR ];
316
317 $pos = 0;
318 foreach my $param (@bind_params) {
319 $sth->bind_param( ++$pos, $param->[0], $param->[1] );
320 }
321
322 $sth->execute();
323
324 my $prefix_len = length($seq_region_name) + 1;
325 my $high_ver = undef;
326 my $high_cs = $cs;
327
328 # Find the fuzzy-matched seq_region with the highest postfix
329 # (which ought to be a version).
330
331 my ( $tmp_name, $id, $tmp_length, $cs_id );
332 $sth->bind_columns( \( $tmp_name, $id, $tmp_length, $cs_id ) );
333
334 my $i = 0;
335
336 while ( $sth->fetch ) {
337 my $tmp_cs =
338 ( 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 $self->{'sr_id_cache'}->{"$id"} = $arr;
344
345 my $tmp_ver = substr( $tmp_name, $prefix_len );
346
347 # skip versions which are non-numeric and apparently not
348 # versions
349 if ( $tmp_ver !~ /^\d+$/ ) { next }
350
351 # take version with highest num, if two versions match take one
352 # with highest ranked coord system (lowest num)
353 if ( !defined($high_ver)
354 || $tmp_ver > $high_ver
355 || ( $tmp_ver == $high_ver && $tmp_cs->rank < $high_cs->rank )
356 )
357 {
358 $seq_region_name = $tmp_name;
359 $length = $tmp_length;
360 $high_ver = $tmp_ver;
361 $high_cs = $tmp_cs;
362 }
363
364 $i++;
365 } ## end while ( $sth->fetch )
366 $sth->finish();
367
368 # warn if fuzzy matching found more than one result
369 if ( $i > 1 ) {
370 warning(
371 sprintf(
372 "Fuzzy matching of seq_region_name "
373 . "returned more than one result.\n"
374 . "You might want to check whether the returned seq_region\n"
375 . "(%s:%s) is the one you intended to fetch.\n",
376 $high_cs->name(), $seq_region_name ) );
377 }
378
379 $cs = $high_cs;
380
381 # return if we did not find any appropriate match:
382 if ( !defined($high_ver) ) { return undef }
383
384 } else {
385
386 my ( $id, $cs_id );
387 ( $seq_region_name, $id, $length, $cs_id ) =
388 $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 $self->{'sr_id_cache'}->{"$id"} = $arr;
395 $cs = $csa->fetch_by_dbID($cs_id);
396 }
397 } ## end else [ if ( defined($arr) ) ]
398
399 if ( !defined($end) ) { $end = $length }
400
401 #If this was given then check if we've got a circular seq region otherwise
402 #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 if($self->is_circular($seq_region_id)) {
407 my $new_sl =
408 Bio::EnsEMBL::CircularSlice->new(
409 -COORD_SYSTEM => $cs,
410 -SEQ_REGION_NAME => $seq_region_name,
411 -SEQ_REGION_LENGTH => $length,
412 -START => $start,
413 -END => $end,
414 -STRAND => 1,
415 -ADAPTOR => $self );
416
417 return $new_sl;
418 }
419 }
420
421 if ( defined( $self->{'lrg_region_test'} )
422 and substr( $cs->name, 0, 3 ) eq $self->{'lrg_region_test'} )
423 {
424 return
425 Bio::EnsEMBL::LRGSlice->new( -COORD_SYSTEM => $cs,
426 -SEQ_REGION_NAME => $seq_region_name,
427 -SEQ_REGION_LENGTH => $length,
428 -START => $start,
429 -END => $end,
430 -STRAND => $strand,
431 -ADAPTOR => $self );
432 } else {
433 return
434 Bio::EnsEMBL::Slice->new_fast( {
435 'coord_system' => $cs,
436 'seq_region_name' => $seq_region_name,
437 'seq_region_length' => $length,
438 'start' => $start,
439 'end' => $end,
440 'strand' => $strand,
441 'adaptor' => $self } );
442 }
443 } ## end sub fetch_by_region
444
445 =head2 fetch_by_toplevel_location
446
447 Arg [1] : string $location
448 Ensembl formatted location. Can be a format like
449 C<name:start-end>, C<name:start..end>, C<name:start:end>,
450 C<name:start>, C<name>. We can also support strand
451 specification as a +/- or 1/-1.
452
453 Location names must be separated by a C<:>. All others can be
454 separated by C<..>, C<:> or C<->.
455 Arg[2] : boolean $no_warnings
456 Suppress warnings from this method
457 Example : my $slice = $sa->fetch_by_toplevel_location('X:1-10000')
458 my $slice = $sa->fetch_by_toplevel_location('X:1-10000:-1')
459 Description : Converts an Ensembl location/region into the sequence region
460 name, start and end and passes them onto C<fetch_by_region()>.
461 The code assumes that the required slice is on the top level
462 coordinate system. The code assumes that location formatting
463 is not perfect and will perform basic cleanup before parsing.
464 Returntype : Bio::EnsEMBL::Slice
465 Exceptions : If $location is false otherwise see C<fetch_by_region()>
466 Caller : General
467 Status : Beta
468
469 =cut
470
471 sub fetch_by_toplevel_location {
472 my ($self, $location, $no_warnings) = @_;
473
474 my ($seq_region_name, $start, $end, $strand) = $self->parse_location_to_values($location, $no_warnings);
475
476 if(! $seq_region_name) {
477 return;
478 }
479
480 if(defined $start && defined $end && $start > $end) {
481 throw "Cannot request a slice whose start is greater than its end. Start: $start. End: $end";
482 }
483
484 my $coord_system_name = 'toplevel';
485 my $slice = $self->fetch_by_region($coord_system_name, $seq_region_name, $start, $end, $strand, undef, 0);
486 return unless $slice;
487
488 my $srl = $slice->seq_region_length();
489 my $name = $slice->seq_region_name();
490 if(defined $start && $start > $srl) {
491 throw "Cannot request a slice whose start ($start) is greater than $srl for $name.";
492 }
493 if(defined $end && $end > $srl) {
494 warning "Requested end ($end) is greater than $srl for $name. Resetting to $srl" if ! $no_warnings;
495 $slice->{end} = $srl;
496 }
497
498 return $slice;
499 }
500
501 =head2 parse_location_to_values
502
503 Arg [1] : string $location
504 Ensembl formatted location. Can be a format like
505 C<name:start-end>, C<name:start..end>, C<name:start:end>,
506 C<name:start>, C<name>. We can also support strand
507 specification as a +/- or 1/-1.
508
509 Location names must be separated by a C<:>. All others can be
510 separated by C<..>, C<:> or C<->.
511 Arg[2] : boolean $no_warnings
512 Suppress warnings from this method
513 Example : my ($name, $start, $end, $strand) = $sa->parse_location_to_values('X:1..100:1);
514 Description : Takes in an Ensembl location String and returns the parsed
515 values
516 Returntype : List. Contains name, start, end and strand
517
518 =cut
519
520
521 sub parse_location_to_values {
522 my ($self, $location, $no_warnings, $no_errors) = @_;
523
524 throw 'You must specify a location' if ! $location;
525
526 #cleanup any nomenclature like 1_000 or 1 000 or 1,000
527 my $number_seps_regex = qr/\s+|,|_/;
528 my $separator_regex = qr/(?:-|[.]{2}|\:)?/;
529 my $number_regex = qr/[0-9,_ E]+/xms;
530 my $strand_regex = qr/[+-1]|-1/xms;
531
532 my $regex = qr/^(\w+) \s* :? \s* ($number_regex)? $separator_regex ($number_regex)? $separator_regex ($strand_regex)? $/xms;
533 my ($seq_region_name, $start, $end, $strand);
534 if(($seq_region_name, $start, $end, $strand) = $location =~ $regex) {
535
536 if(defined $strand) {
537 if(!looks_like_number($strand)) {
538 $strand = ($strand eq '+') ? 1 : -1;
539 }
540 }
541
542 if(defined $start) {
543 $start =~ s/$number_seps_regex//g;
544 if($start < 1) {
545 warning "Start was less than 1 (${start}) which is not allowed. Resetting to 1" if ! $no_warnings;
546 $start = 1;
547 }
548 }
549 if(defined $end) {
550 $end =~ s/$number_seps_regex//g;
551 if($end < 1) {
552 throw "Cannot request negative or 0 end indexes through this interface. Given $end but expected something greater than 0" unless $no_errors;
553 }
554 }
555
556 if(defined $start && defined $end && $start > $end) {
557 throw "Cannot request a slice whose start is greater than its end. Start: $start. End: $end" unless $no_errors;
558 }
559 }
560
561 return ($seq_region_name, $start, $end, $strand);
562 }
563
564 =head2 fetch_by_region_unique
565
566 Arg [1] : string $coord_system_name (optional)
567 The name of the coordinate system of the slice to be created
568 This may be a name of an actual coordinate system or an alias
569 to a coordinate system. Valid aliases are 'seqlevel' or
570 'toplevel'.
571 Arg [2] : string $seq_region_name
572 The name of the sequence region that the slice will be
573 created on.
574 Arg [3] : int $start (optional, default = 1)
575 The start of the slice on the sequence region
576 Arg [4] : int $end (optional, default = seq_region length)
577 The end of the slice on the sequence region
578 Arg [5] : int $strand (optional, default = 1)
579 The orientation of the slice on the sequence region
580 Arg [6] : string $version (optional, default = default version)
581 The version of the coordinate system to use (e.g. NCBI33)
582 Arg [7] : boolean $no_fuzz (optional, default = undef (false))
583 If true (non-zero), do not use "fuzzy matching" (see below).
584 Example : $slice = $slice_adaptor->fetch_by_region_unique('chromosome', 'HSCHR6_MHC_COX');
585 Description: Retrieves a slice on the requested region but returns only the unique
586 parts of the slice. At a minimum the
587 name the name of the seq_region to fetch must be provided.
588
589 If no coordinate system name is provided than a slice on the
590 highest ranked coordinate system with a matching
591 seq_region_name will be returned. If a version but no
592 coordinate system name is provided, the same behaviour will
593 apply, but only coordinate systems of the appropriate version
594 are considered. The same applies if the 'toplevel' coordinate
595 system is specified, however in this case the version is
596 ignored. The coordinate system should always be specified if
597 it is known, since this is unambiguous and faster.
598
599 Some fuzzy matching is performed if no exact match for
600 the provided name is found. This allows clones to be
601 fetched even when their version is not known. For
602 example fetch_by_region('clone', 'AC008066') will
603 retrieve the sequence_region with name 'AC008066.4'.
604
605 The fuzzy matching can be turned off by setting the
606 $no_fuzz argument to a true value.
607
608 If the requested seq_region is not found in the database undef
609 is returned.
610
611 Returntype : listref Bio::EnsEMBL::Slice
612 Exceptions : throw if no seq_region_name is provided
613 throw if invalid coord_system_name is provided
614 throw if start > end is provided
615 Caller : general
616 Status : Stable
617
618 =cut
619
620 sub fetch_by_region_unique {
621 my $self = shift;
622
623 my @out = ();
624 my $slice = $self->fetch_by_region(@_);
625
626
627 if ( !exists( $self->{'asm_exc_cache'} ) ) {
628 $self->_build_exception_cache();
629 }
630
631 if ( exists(
632 $self->{'asm_exc_cache'}->{ $self->get_seq_region_id($slice) }
633 ) )
634 {
635 # Dereference symlinked assembly regions. Take out any regions
636 # which are symlinked because these are duplicates.
637 my @projection =
638 @{ $self->fetch_normalized_slice_projection($slice) };
639
640 foreach my $segment (@projection) {
641 if ( $segment->[2]->seq_region_name() eq $slice->seq_region_name()
642 && $segment->[2]->coord_system->equals( $slice->coord_system ) )
643 {
644 push( @out, $segment->[2] );
645 }
646 }
647 }
648
649 return \@out;
650 } ## end sub fetch_by_region_unique
651
652 =head2 fetch_by_name
653
654 Arg [1] : string $name
655 Example : $name = 'chromosome:NCBI34:X:1000000:2000000:1';
656 $slice = $slice_adaptor->fetch_by_name($name);
657 $slice2 = $slice_adaptor->fetch_by_name($slice3->name());
658 Description: Fetches a slice using a slice name (i.e. the value returned by
659 the Slice::name method). This is useful if you wish to
660 store a unique identifier for a slice in a file or database or
661 pass a slice over a network.
662 Slice::name allows you to serialise/marshall a slice and this
663 method allows you to deserialise/unmarshal it.
664
665 Returns undef if no seq_region with the provided name exists in
666 the database.
667
668 Returntype : Bio::EnsEMBL::Slice or undef
669 Exceptions : throw if incorrent arg provided
670 Caller : Pipeline
671 Status : Stable
672
673 =cut
674
675 sub fetch_by_name {
676 my $self = shift;
677 my $name = shift;
678
679 if(!$name) {
680 throw("name argument is required");
681 }
682
683 my @array = split(/:/,$name);
684
685 if(scalar(@array) < 3 || scalar(@array) > 6) {
686 throw("Malformed slice name [$name]. Format is " .
687 "coord_system:version:name:start:end:strand");
688 }
689
690 # Rearrange arguments to suit fetch_by_region
691
692 my @targetarray;
693
694 $targetarray[0]=$array[0];
695 $targetarray[5]=(($array[1]&&$array[1] ne "")?$array[1]:undef);
696 $targetarray[1]=(($array[2]&&$array[2] ne "")?$array[2]:undef);
697 $targetarray[2]=(($array[3]&&$array[3] ne "")?$array[3]:undef);
698 $targetarray[3]=(($array[4]&&$array[4] ne "")?$array[4]:undef);
699 $targetarray[4]=(($array[5]&&$array[5] ne "")?$array[5]:undef);
700 return $self->fetch_by_region(@targetarray);
701 }
702
703
704
705 =head2 fetch_by_seq_region_id
706
707 Arg [1] : string $seq_region_id
708 The internal identifier of the seq_region to create this slice
709 on
710 Arg [2] : optional start
711 Arg [3] : optional end
712 Arg [4] : optional strand
713 Example : $slice = $slice_adaptor->fetch_by_seq_region_id(34413);
714 Description: Creates a slice object of an entire seq_region using the
715 seq_region internal identifier to resolve the seq_region.
716 Returns undef if no such slice exists.
717 Returntype : Bio::EnsEMBL::Slice or undef
718 Exceptions : none
719 Caller : general
720 Status : Stable
721
722 =cut
723
724 sub fetch_by_seq_region_id {
725 my ( $self, $seq_region_id, $start, $end, $strand ) = @_;
726
727 my $arr = $self->{'sr_id_cache'}->{$seq_region_id};
728 my ( $name, $length, $cs, $cs_id );
729
730
731 if ( $arr && defined( $arr->[2] ) ) {
732 ( $name, $cs_id, $length ) = ( $arr->[1], $arr->[2], $arr->[3] );
733 $cs = $self->db->get_CoordSystemAdaptor->fetch_by_dbID($cs_id);
734 } else {
735 my $sth =
736 $self->prepare( "SELECT sr.name, sr.coord_system_id, sr.length "
737 . "FROM seq_region sr "
738 . "WHERE sr.seq_region_id = ? " );
739
740 $sth->bind_param( 1, $seq_region_id, SQL_INTEGER );
741 $sth->execute();
742
743 if ( $sth->rows() == 0 ) { return undef }
744
745 ( $name, $cs_id, $length ) = $sth->fetchrow_array();
746 $sth->finish();
747
748 $cs = $self->db->get_CoordSystemAdaptor->fetch_by_dbID($cs_id);
749
750 #cache results to speed up repeated queries
751 my $arr = [ $seq_region_id, $name, $cs_id, $length ];
752
753 $self->{'sr_name_cache'}->{"$name:$cs_id"} = $arr;
754 $self->{'sr_id_cache'}->{"$seq_region_id"} = $arr;
755 }
756
757 return
758 Bio::EnsEMBL::Slice->new_fast({
759 'coord_system' => $cs,
760 'seq_region_name' => $name,
761 'seq_region_length'=> $length,
762 'start' => $start || 1,
763 'end' => $end || $length,
764 'strand' => $strand || 1,
765 'adaptor' => $self} );
766 } ## end sub fetch_by_seq_region_id
767
768
769
770 =head2 get_seq_region_id
771
772 Arg [1] : Bio::EnsEMBL::Slice $slice
773 The slice to fetch a seq_region_id for
774 Example : $srid = $slice_adaptor->get_seq_region_id($slice);
775 Description: Retrieves the seq_region id (in this database) given a slice
776 Seq region ids are not stored on the slices themselves
777 because they are intended to be somewhat database independant
778 and seq_region_ids vary accross databases.
779 Returntype : int
780 Exceptions : throw if the seq_region of the slice is not in the db
781 throw if incorrect arg provided
782 Caller : BaseFeatureAdaptor
783 Status : Stable
784
785 =cut
786
787 sub get_seq_region_id {
788 my $self = shift;
789 my $slice = shift;
790
791 if(!$slice || !ref($slice) || !($slice->isa('Bio::EnsEMBL::Slice') or $slice->isa('Bio::EnsEMBL::LRGSlice'))) {
792 throw('Slice argument is required');
793 }
794
795 my $seq_region_name = $slice->seq_region_name();
796 my $key = $seq_region_name.":".$slice->coord_system->dbID();
797 my $arr = $self->{'sr_name_cache'}->{"$key"};
798
799 if( $arr ) {
800 return $arr->[0];
801 }
802
803 my $cs_id = $slice->coord_system->dbID();
804
805 my $sth = $self->prepare("SELECT seq_region_id, length " .
806 "FROM seq_region " .
807 "WHERE name = ? AND coord_system_id = ?");
808
809 #force seq_region_name cast to string so mysql cannot treat as int
810 $sth->bind_param(1,"$seq_region_name",SQL_VARCHAR);
811 $sth->bind_param(2,$cs_id,SQL_INTEGER);
812 $sth->execute();
813
814 if($sth->rows() != 1) {
815 throw("Non existant or ambigous seq_region:\n" .
816 " coord_system=[$cs_id],\n" .
817 " name=[$seq_region_name],\n");
818
819 }
820
821 my($seq_region_id, $length) = $sth->fetchrow_array();
822 $sth->finish();
823
824 #cache information for future requests
825 $arr = [ $seq_region_id, $seq_region_name, $cs_id, $length ];
826
827 $self->{'sr_name_cache'}->{"$seq_region_name:$cs_id"} = $arr;
828 $self->{'sr_id_cache'}->{"$seq_region_id"} = $arr;
829
830 return $seq_region_id;
831 }
832
833
834
835 =head2 fetch_all
836
837 Arg [1] : string $coord_system_name
838 The name of the coordinate system to retrieve slices of.
839 This may be a name of an acutal coordinate system or an alias
840 to a coordinate system. Valid aliases are 'seqlevel' or
841 'toplevel'.
842 Arg [2] : string $coord_system_version (optional)
843 The version of the coordinate system to retrieve slices of
844 Arg [3] : bool $include_non_reference (optional)
845 If this argument is not provided then only reference slices
846 will be returned. If set, both reference and non refeference
847 slices will be rerurned.
848 Arg [4] : int $include_duplicates (optional)
849 If set duplicate regions will be returned.
850
851 NOTE: if you do not use this option and you have a PAR
852 (pseudo-autosomal region) at the beginning of your seq_region
853 then your slice will not start at position 1, so coordinates
854 retrieved from this slice might not be what you expected.
855
856 Arg[5] : bool $include_lrg (optional) (default 0)
857 If set lrg regions will be returned aswell.
858
859
860 Example : @chromos = @{$slice_adaptor->fetch_all('chromosome','NCBI33')};
861 @contigs = @{$slice_adaptor->fetch_all('contig')};
862
863 # get even non-reference regions
864 @slices = @{$slice_adaptor->fetch_all('toplevel',undef,1)};
865
866 # include duplicate regions (such as pseudo autosomal regions)
867 @slices = @{$slice_adaptor->fetch_all('toplevel', undef,0,1)};
868
869 Description: Retrieves slices of all seq_regions for a given coordinate
870 system. This is analagous to the methods fetch_all which were
871 formerly on the ChromosomeAdaptor, RawContigAdaptor and
872 CloneAdaptor classes. Slices fetched span the entire
873 seq_regions and are on the forward strand.
874 If the coordinate system with the provided name and version
875 does not exist an empty list is returned.
876 If the coordinate system name provided is 'toplevel', all
877 non-redundant toplevel slices are returned (note that any
878 coord_system_version argument is ignored in that case).
879
880 Retrieved slices can be broken into smaller slices using the
881 Bio::EnsEMBL::Utils::Slice module.
882
883 Returntype : listref of Bio::EnsEMBL::Slices
884 Exceptions : none
885 Caller : general
886 Status : Stable
887
888 =cut
889
890 sub fetch_all {
891 my $self = shift;
892 my $cs_name = shift;
893 my $cs_version = shift || '';
894
895 my ($include_non_reference, $include_duplicates, $include_lrg) = @_;
896
897 #
898 # verify existance of requested coord system and get its id
899 #
900 my $csa = $self->db->get_CoordSystemAdaptor();
901 my $orig_cs = $csa->fetch_by_name($cs_name, $cs_version);
902
903 return [] if ( !$orig_cs );
904
905 my %bad_vals=();
906
907
908 #
909 # Get a hash of non reference seq regions
910 #
911 if ( !$include_non_reference ) {
912 my $sth =
913 $self->prepare( 'SELECT sr.seq_region_id '
914 . 'FROM seq_region sr, seq_region_attrib sra, '
915 . 'attrib_type at, coord_system cs '
916 . 'WHERE at.code = "non_ref" '
917 . 'AND sra.seq_region_id = sr.seq_region_id '
918 . 'AND at.attrib_type_id = sra.attrib_type_id '
919 . 'AND sr.coord_system_id = cs.coord_system_id '
920 . 'AND cs.species_id = ?' );
921
922 $sth->bind_param( 1, $self->species_id(), SQL_INTEGER );
923 $sth->execute();
924
925 my ($seq_region_id);
926 $sth->bind_columns( \$seq_region_id );
927
928 while ( $sth->fetch() ) {
929 $bad_vals{$seq_region_id} = 1;
930 }
931 }
932
933 #
934 # if we do not want lrg's then add them to the bad list;
935 #
936 if ( !$include_lrg ) {
937 my $sth =
938 $self->prepare( 'SELECT sr.seq_region_id '
939 . 'FROM seq_region sr, seq_region_attrib sra, '
940 . 'attrib_type at, coord_system cs '
941 . 'WHERE at.code = "LRG" '
942 . 'AND sra.seq_region_id = sr.seq_region_id '
943 . 'AND at.attrib_type_id = sra.attrib_type_id '
944 . 'AND sr.coord_system_id = cs.coord_system_id '
945 . 'AND cs.species_id = ?' );
946
947 $sth->bind_param( 1, $self->species_id(), SQL_INTEGER );
948 $sth->execute();
949
950 my ($seq_region_id);
951 $sth->bind_columns( \$seq_region_id );
952
953 while ( $sth->fetch() ) {
954 $bad_vals{$seq_region_id} = 1;
955 }
956 }
957
958 #
959 # Retrieve the seq_regions from the database
960 #
961
962 my $sth;
963 if ( $orig_cs->is_top_level() ) {
964 $sth =
965 $self->prepare( 'SELECT sr.seq_region_id, sr.name, '
966 . 'sr.length, sr.coord_system_id '
967 . 'FROM seq_region sr, seq_region_attrib sra, '
968 . 'attrib_type at, coord_system cs '
969 . 'WHERE at.code = "toplevel" '
970 . 'AND at.attrib_type_id = sra.attrib_type_id '
971 . 'AND sra.seq_region_id = sr.seq_region_id '
972 . 'AND sr.coord_system_id = cs.coord_system_id '
973 . 'AND cs.species_id = ?' );
974
975 $sth->bind_param( 1, $self->species_id(), SQL_INTEGER );
976 $sth->execute();
977 } else {
978 $sth =
979 $self->prepare( 'SELECT sr.seq_region_id, sr.name, '
980 . 'sr.length, sr.coord_system_id '
981 . 'FROM seq_region sr '
982 . 'WHERE sr.coord_system_id = ?' );
983
984 $sth->bind_param( 1, $orig_cs->dbID, SQL_INTEGER );
985 $sth->execute();
986 }
987
988 my ( $seq_region_id, $name, $length, $cs_id );
989 $sth->bind_columns( \( $seq_region_id, $name, $length, $cs_id ) );
990
991 my $cache_count = 0;
992
993 my @out;
994 while($sth->fetch()) {
995 if(!defined($bad_vals{$seq_region_id})){
996 my $cs = $csa->fetch_by_dbID($cs_id);
997
998 if(!$cs) {
999 throw("seq_region $name references non-existent coord_system $cs_id.");
1000 }
1001
1002 #cache values for future reference, but stop adding to the cache once we
1003 #we know we have filled it up
1004 if($cache_count < $Bio::EnsEMBL::Utils::SeqRegionCache::SEQ_REGION_CACHE_SIZE) {
1005 my $arr = [ $seq_region_id, $name, $cs_id, $length ];
1006
1007 $self->{'sr_name_cache'}->{"$name:$cs_id"} = $arr;
1008 $self->{'sr_id_cache'}->{"$seq_region_id"} = $arr;
1009
1010 $cache_count++;
1011 }
1012
1013 my $slice = Bio::EnsEMBL::Slice->new_fast({
1014 'start' => 1,
1015 'end' => $length,
1016 'strand' => 1,
1017 'seq_region_name' => $name,
1018 'seq_region_length'=> $length,
1019 'coord_system' => $cs,
1020 'adaptor' => $self});
1021
1022 if(!defined($include_duplicates) or !$include_duplicates){
1023 # test if this slice *could* have a duplicate (exception) region
1024 $self->_build_exception_cache() if(!exists $self->{'asm_exc_cache'});
1025 if(exists $self->{asm_exc_cache}->{$seq_region_id}) {
1026
1027 # Dereference symlinked assembly regions. Take out
1028 # any regions which are symlinked because these are duplicates
1029 my @projection = @{$self->fetch_normalized_slice_projection($slice)};
1030 foreach my $segment ( @projection) {
1031 if($segment->[2]->seq_region_name() eq $slice->seq_region_name() &&
1032 $segment->[2]->coord_system->equals($slice->coord_system)) {
1033 push @out, $segment->[2];
1034 }
1035 }
1036 } else {
1037 # no duplicate regions
1038 push @out, $slice;
1039 }
1040 } else {
1041 # we want duplicates anyway so do not do any checks
1042 push @out, $slice;
1043 }
1044 }
1045 }
1046
1047 return \@out;
1048 }
1049
1050 =head2 is_toplevel
1051 Arg : int seq_region_id
1052 Example : my $top = $slice_adptor->is_toplevel($seq_region_id)
1053 Description: Returns 1 if slice is a toplevel slice else 0
1054 Returntype : int
1055 Caller : Slice method is_toplevel
1056 Status : At Risk
1057
1058 =cut
1059
1060 sub is_toplevel {
1061 my $self = shift;
1062 my $id = shift;
1063
1064 my $sth = $self->prepare(
1065 "SELECT at.code from seq_region_attrib sra, attrib_type at "
1066 . "WHERE sra.seq_region_id = ? "
1067 . "AND at.attrib_type_id = sra.attrib_type_id "
1068 . "AND at.code = 'toplevel'" );
1069
1070 $sth->bind_param( 1, $id, SQL_INTEGER );
1071 $sth->execute();
1072
1073 my $code;
1074 $sth->bind_columns( \$code );
1075
1076 while ( $sth->fetch ) {
1077 $sth->finish;
1078 return 1;
1079 }
1080
1081 $sth->finish;
1082 return 0;
1083 }
1084
1085 =head2 is_reference
1086 Arg : int seq_region_id
1087 Example : my $reference = $slice_adptor->is_reference($seq_region_id)
1088 Description: Returns 1 if slice is a reference slice else 0
1089 Returntype : int
1090 Caller : Slice method is_reference
1091 Status : At Risk
1092
1093 =cut
1094
1095 sub is_reference {
1096 my $self = shift;
1097 my $id = shift;
1098
1099 my $sth = $self->prepare(
1100 "SELECT at.code from seq_region_attrib sra, attrib_type at "
1101 . "WHERE sra.seq_region_id = ? "
1102 . "AND at.attrib_type_id = sra.attrib_type_id "
1103 . "AND at.code = 'non_ref'" );
1104
1105 $sth->bind_param( 1, $id, SQL_INTEGER );
1106 $sth->execute();
1107
1108 my $code;
1109 $sth->bind_columns( \$code );
1110
1111 while ( $sth->fetch ) {
1112 $sth->finish;
1113 return 0;
1114 }
1115
1116 $sth->finish;
1117 return 1;
1118 }
1119
1120 =head2 is_circular
1121
1122 Arg[1] : int seq_region_id
1123 Example : my $circular = $slice_adptor->is_circular($seq_region_id);
1124 Description : Indicates if the sequence region was circular or not
1125 Returntype : Boolean
1126
1127 =cut
1128
1129 sub is_circular {
1130 my ($self, $id) = @_;
1131
1132 if (! defined $self->{is_circular}) {
1133 $self->_build_circular_slice_cache();
1134 }
1135
1136 return 0 if $self->{is_circular} == 0;
1137 return (exists $self->{circular_sr_id_cache}->{$id}) ? 1 : 0;
1138 }
1139
1140 =head2 fetch_by_band
1141
1142 Title : fetch_by_band
1143 Usage :
1144 Function: Does not work please use fetch_by_chr_band
1145 Example :
1146 Returns : Bio::EnsEMBL::Slice
1147 Args : the band name
1148 Status : AT RISK
1149
1150 =cut
1151
1152 sub fetch_by_band {
1153 my ($self,$band) = @_;
1154
1155 my $sth = $self->dbc->prepare
1156 ("select s.name,max(k.seq_region_id)-min(k.seq_region_id, min(k.seq_region_start), max(k.seq_region_id) " .
1157 "from karyotype as k " .
1158 "where k.band like ? and k.seq_region_id = s.seq_region_id");
1159
1160 $sth->bind_param(1,"$band%",SQL_VARCHAR);
1161 $sth->execute();
1162 my ( $seq_region_name, $discrepancy, $seq_region_start, $seq_region_end) = $sth->fetchrow_array;
1163
1164 if($seq_region_name && $discrepancy>0) {
1165 throw("Band maps to multiple seq_regions");
1166 } else {
1167 return $self->fetch_by_region('toplevel',$seq_region_name,$seq_region_start,$seq_region_end);
1168 }
1169 throw("Band not recognised in database");
1170 }
1171
1172 =head2 fetch_by_chr_band
1173
1174 Title : fetch_by_chr_band
1175 Usage :
1176 Function: create a Slice representing a series of bands
1177 Example :
1178 Returns : Bio::EnsEMBL::Slice
1179 Args : the band name
1180 Status : Stable
1181
1182 =cut
1183
1184 sub fetch_by_chr_band {
1185 my ( $self, $chr, $band ) = @_;
1186
1187 my $chr_slice = $self->fetch_by_region( 'toplevel', $chr );
1188 my $seq_region_id = $self->get_seq_region_id($chr_slice);
1189
1190 my $sth =
1191 $self->prepare( 'SELECT MIN(k.seq_region_start), '
1192 . 'MAX(k.seq_region_end) '
1193 . 'FROM karyotype k '
1194 . 'WHERE k.seq_region_id = ? '
1195 . 'AND k.band LIKE ?' );
1196
1197 $sth->bind_param( 1, $seq_region_id, SQL_INTEGER );
1198 $sth->bind_param( 2, "$band%", SQL_VARCHAR );
1199 $sth->execute();
1200
1201 my ( $slice_start, $slice_end ) = $sth->fetchrow_array;
1202
1203 if ( defined $slice_start ) {
1204 return
1205 $self->fetch_by_region( 'toplevel', $chr,
1206 $slice_start, $slice_end );
1207 }
1208
1209 throw("Band not recognised in database");
1210 } ## end sub fetch_by_chr_band
1211
1212
1213
1214 =head2 fetch_by_exon_stable_id
1215
1216 Arg [1] : string $exonid
1217 The stable id of the exon around which the slice is
1218 desired
1219 Arg [2] : (optional) int $size
1220 The length of the flanking regions the slice should encompass
1221 on either side of the exon (0 by default)
1222 Example : $slc = $sa->fetch_by_exon_stable_id('ENSE00000302930',10);
1223 Description: Creates a slice around the region of the specified exon.
1224 If a context size is given, the slice is extended by that
1225 number of basepairs on either side of the exon.
1226
1227 The slice will be created in the exon's native coordinate system
1228 and in the forward orientation.
1229 Returntype : Bio::EnsEMBL::Slice
1230 Exceptions : Thrown if the exon is not in the database.
1231 Caller : general
1232 Status : Stable
1233
1234 =cut
1235
1236 sub fetch_by_exon_stable_id{
1237 my ($self,$exonid,$size) = @_;
1238
1239 throw('Exon argument is required.') if(!$exonid);
1240
1241 my $ea = $self->db->get_ExonAdaptor;
1242 my $exon = $ea->fetch_by_stable_id($exonid);
1243
1244 throw("Exon [$exonid] does not exist in DB.") if(!$exon);
1245
1246 return $self->fetch_by_Feature($exon, $size);
1247 }
1248
1249 =head2 fetch_by_transcript_stable_id
1250
1251 Arg [1] : string $transcriptid
1252 The stable id of the transcript around which the slice is
1253 desired
1254 Arg [2] : (optional) int $size
1255 The length of the flanking regions the slice should encompass
1256 on either side of the transcript (0 by default)
1257 Example : $slc = $sa->fetch_by_transcript_stable_id('ENST00000302930',10);
1258 Description: Creates a slice around the region of the specified transcript.
1259 If a context size is given, the slice is extended by that
1260 number of basepairs on either side of the
1261 transcript.
1262
1263 The slice will be created in the transcript's native coordinate
1264 system and in the forward orientation.
1265 Returntype : Bio::EnsEMBL::Slice
1266 Exceptions : Thrown if the transcript is not in the database.
1267 Caller : general
1268 Status : Stable
1269
1270 =cut
1271
1272 sub fetch_by_transcript_stable_id{
1273 my ($self,$transcriptid,$size) = @_;
1274
1275 throw('Transcript argument is required.') if(!$transcriptid);
1276
1277 my $ta = $self->db->get_TranscriptAdaptor;
1278 my $transcript = $ta->fetch_by_stable_id($transcriptid);
1279
1280 throw("Transcript [$transcriptid] does not exist in DB.") if(!$transcript);
1281
1282 return $self->fetch_by_Feature($transcript, $size);
1283 }
1284
1285
1286
1287 =head2 fetch_by_transcript_id
1288
1289 Arg [1] : int $transcriptid
1290 The unique database identifier of the transcript around which
1291 the slice is desired
1292 Arg [2] : (optional) int $size
1293 The length of the flanking regions the slice should encompass
1294 on either side of the transcript (0 by default)
1295 Example : $slc = $sa->fetch_by_transcript_id(24, 1000);
1296 Description: Creates a slice around the region of the specified transcript.
1297 If a context size is given, the slice is extended by that
1298 number of basepairs on either side of the
1299 transcript.
1300
1301 The slice will be created in the transcript's native coordinate
1302 system and in the forward orientation.
1303 Returntype : Bio::EnsEMBL::Slice
1304 Exceptions : throw on incorrect args
1305 throw if transcript is not in database
1306 Caller : general
1307 Status : Stable
1308
1309 =cut
1310
1311 sub fetch_by_transcript_id {
1312 my ($self,$transcriptid,$size) = @_;
1313
1314 throw('Transcript id argument is required.') if(!$transcriptid);
1315
1316 my $transcript_adaptor = $self->db()->get_TranscriptAdaptor();
1317 my $transcript = $transcript_adaptor->fetch_by_dbID($transcriptid);
1318
1319 throw("Transcript [$transcriptid] does not exist in DB.") if(!$transcript);
1320
1321 return $self->fetch_by_Feature($transcript, $size);
1322 }
1323
1324
1325
1326 =head2 fetch_by_gene_stable_id
1327
1328 Arg [1] : string $geneid
1329 The stable id of the gene around which the slice is
1330 desired
1331 Arg [2] : (optional) int $size
1332 The length of the flanking regions the slice should encompass
1333 on either side of the gene (0 by default)
1334 Example : $slc = $sa->fetch_by_gene_stable_id('ENSG00000012123',10);
1335 Description: Creates a slice around the region of the specified gene.
1336 If a context size is given, the slice is extended by that
1337 number of basepairs on either side of the gene.
1338
1339 The slice will be created in the gene's native coordinate system
1340 and in the forward orientation.
1341 Returntype : Bio::EnsEMBL::Slice
1342 Exceptions : throw on incorrect args
1343 throw if transcript does not exist
1344 Caller : general
1345 Status : Stable
1346
1347 =cut
1348
1349 sub fetch_by_gene_stable_id {
1350 my ($self,$geneid,$size) = @_;
1351
1352 throw('Gene argument is required.') if(!$geneid);
1353
1354 my $gene_adaptor = $self->db->get_GeneAdaptor();
1355 my $gene = $gene_adaptor->fetch_by_stable_id($geneid);
1356
1357 throw("Gene [$geneid] does not exist in DB.") if(!$gene);
1358
1359 return $self->fetch_by_Feature($gene, $size);
1360 }
1361
1362
1363
1364 =head2 fetch_by_Feature
1365
1366 Arg [1] : Bio::EnsEMBL::Feature $feat
1367 The feature to fetch the slice around
1368 Arg [2] : int size (optional)
1369 The desired number of flanking basepairs around the feature.
1370 The size may also be provided as a percentage of the feature
1371 size such as 200% or 80.5%.
1372 Example : $slice = $slice_adaptor->fetch_by_Feature($feat, 100);
1373 Description: Retrieves a slice around a specific feature. All this really
1374 does is return a resized version of the slice that the feature
1375 is already on. Note that slices returned from this method
1376 are always on the forward strand of the seq_region regardless of
1377 the strandedness of the feature passed in.
1378 Returntype : Bio::EnsEMBL::Slice
1379 Exceptions : throw if the feature does not have an attached slice
1380 throw if feature argument is not provided
1381 Caller : fetch_by_gene_stable_id, fetch_by_transcript_stable_id,
1382 fetch_by_gene_id, fetch_by_transcript_id
1383 Status : Stable
1384
1385 =cut
1386
1387 sub fetch_by_Feature{
1388 my ($self, $feature, $size) = @_;
1389
1390 $size ||= 0;
1391
1392 if(!ref($feature) || !$feature->isa('Bio::EnsEMBL::Feature')) {
1393 throw('Feature argument expected.');
1394 }
1395
1396 my $slice = $feature->slice();
1397 if(!$slice || !($slice->isa('Bio::EnsEMBL::Slice') or $slice->isa('Bio::EnsEMBL::LRGSlice') )) {
1398 throw('Feature must be attached to a valid slice.');
1399 }
1400
1401
1402 my $fstart = $feature->start();
1403 my $fend = $feature->end();
1404 if(!defined($fstart) || !defined($fend)) {
1405 throw('Feature must have defined start and end.');
1406 }
1407
1408 #convert the feature slice coordinates to seq_region coordinates
1409 my $slice_start = $slice->start();
1410 my $slice_end = $slice->end();
1411 my $slice_strand = $slice->strand();
1412 if($slice_start != 1 || $slice_strand != 1) {
1413 if($slice_strand == 1) {
1414 $fstart = $fstart + $slice_start - 1;
1415 $fend = $fend + $slice_start - 1;
1416 } else {
1417 my $tmp_start = $fstart;
1418 $fstart = $slice_end - $fend + 1;
1419 $fend = $slice_end - $tmp_start + 1;
1420 }
1421 }
1422
1423 ## Size may be stored as a %age of the length of the feature
1424 ## Size = 100% gives no context
1425 ## Size = 200% gives context - 50% the size of the feature either side of
1426 ## feature
1427
1428 $size = int( ($1-100)/200 * ($fend-$fstart+1) ) if( $size =~/([\d+\.]+)%/ );
1429
1430 #return a new slice covering the region of the feature
1431 my $S = Bio::EnsEMBL::Slice->new_fast({
1432 'seq_region_name' => $slice->seq_region_name,
1433 'seq_region_length' => $slice->seq_region_length,
1434 'coord_system' => $slice->coord_system,
1435 'start' => $fstart - $size,
1436 'end' => $fend + $size,
1437 'strand' => 1,
1438 'adaptor' => $self});
1439 $S->{'_raw_feature_strand'} = $feature->strand * $slice_strand if $feature->can('strand');
1440 return $S;
1441 }
1442
1443
1444
1445 =head2 fetch_by_misc_feature_attribute
1446
1447 Arg [1] : string $attribute_type
1448 The code of the attribute type
1449 Arg [2] : (optional) string $attribute_value
1450 The value of the attribute to fetch by
1451 Arg [3] : (optional) int $size
1452 The amount of flanking region around the misc feature desired.
1453 Example : $slice = $sa->fetch_by_misc_feature_attribute('superctg',
1454 'NT_030871');
1455 $slice = $sa->fetch_by_misc_feature_attribute('synonym',
1456 'AL00012311',
1457 $flanking);
1458 Description: Fetches a slice around a MiscFeature with a particular
1459 attribute type and value. If no value is specified then
1460 the feature with the particular attribute is used.
1461 If no size is specified then 0 is used.
1462 Returntype : Bio::EnsEMBL::Slice
1463 Exceptions : Throw if no feature with the specified attribute type and value
1464 exists in the database
1465 Warning if multiple features with the specified attribute type
1466 and value exist in the database.
1467 Caller : webcode
1468 Status : Stable
1469
1470 =cut
1471
1472 sub fetch_by_misc_feature_attribute {
1473 my ($self, $attrib_type_code, $attrib_value, $size) = @_;
1474
1475 my $mfa = $self->db()->get_MiscFeatureAdaptor();
1476
1477 my $feats = $mfa->fetch_all_by_attribute_type_value($attrib_type_code,
1478 $attrib_value);
1479
1480 if(@$feats == 0) {
1481 throw("MiscFeature with $attrib_type_code=$attrib_value does " .
1482 "not exist in DB.");
1483 }
1484
1485 if(@$feats > 1) {
1486 warning("MiscFeature with $attrib_type_code=$attrib_value is " .
1487 "ambiguous - using first one found.");
1488 }
1489
1490 my ($feat) = @$feats;
1491
1492 return $self->fetch_by_Feature($feat, $size);
1493 }
1494
1495 =head2 fetch_normalized_slice_projection
1496
1497 Arg [1] : Bio::EnsEMBL::Slice $slice
1498 Example : ( optional )
1499 Description: gives back a project style result. The returned slices
1500 represent the areas to which there are symlinks for the
1501 given slice. start, end show which area on given slice is
1502 symlinked
1503 Returntype : [[start,end,$slice][]]
1504 Exceptions : none
1505 Caller : BaseFeatureAdaptor
1506 Status : Stable
1507
1508 =cut
1509
1510
1511 sub fetch_normalized_slice_projection {
1512 my $self = shift;
1513 my $slice = shift;
1514
1515 my $slice_seq_region_id = $self->get_seq_region_id( $slice );
1516
1517 $self->_build_exception_cache() if(!exists($self->{'asm_exc_cache'}));
1518
1519 my $result = $self->{'asm_exc_cache'}->{$slice_seq_region_id};
1520
1521 $result ||= [];
1522
1523 my (@haps, @pars);
1524
1525 foreach my $row (@$result) {
1526 my ( $seq_region_id, $seq_region_start, $seq_region_end,
1527 $exc_type, $exc_seq_region_id, $exc_seq_region_start,
1528 $exc_seq_region_end ) = @$row;
1529
1530 # need overlapping PAR and all HAPs if any
1531 if( $exc_type eq "PAR" ) {
1532 if( $seq_region_start <= $slice->end() &&
1533 $seq_region_end >= $slice->start() ) {
1534 push( @pars, [ $seq_region_start, $seq_region_end, $exc_seq_region_id,
1535 $exc_seq_region_start, $exc_seq_region_end ] );
1536 }
1537 } else {
1538 push( @haps, [ $seq_region_start, $seq_region_end, $exc_seq_region_id,
1539 $exc_seq_region_start, $exc_seq_region_end ] );
1540 }
1541 }
1542
1543 if(!@pars && !@haps) {
1544 #just return this slice, there were no haps or pars
1545 return [bless ( [1,$slice->length, $slice], "Bio::EnsEMBL::ProjectionSegment")];
1546 }
1547
1548 my @syms;
1549 if( @haps >= 1 ) {
1550 my @sort_haps = sort { $a->[1] <=> $b->[1] } @haps;
1551
1552 my $count =0;
1553 my $chr_start = 1;
1554 my $hap_start = 1;
1555 my $last = 0;
1556
1557 my $seq_reg_slice = $self->fetch_by_seq_region_id($slice_seq_region_id);
1558 my $exc_slice = $self->fetch_by_seq_region_id( $sort_haps[0][2] );
1559 my $len1 = $seq_reg_slice->length();
1560 my $len2 = $exc_slice->length();
1561 my $max_len = ($len1 > $len2) ? $len1 : $len2;
1562
1563 while($count <= scalar(@sort_haps) and !$last){
1564 my $chr_end;
1565 my $hap_end;
1566 if(defined($sort_haps[$count]) and defined($sort_haps[$count][0]) ){
1567 $hap_end = $sort_haps[$count][0]-1;
1568 $chr_end = $sort_haps[$count][3]-1
1569 }
1570 else{
1571 $last = 1;
1572 $hap_end = $len1;
1573 $chr_end = $len2;
1574 my $diff = ($hap_end-$hap_start)-($chr_end-$chr_start);
1575 if($diff > 0){
1576 push( @syms, [ $hap_start, $hap_end, $sort_haps[0][2], $chr_start, $chr_end+$diff] );
1577 }
1578 elsif($diff < 0){
1579 push( @syms, [ $hap_start, $hap_end - $diff, $sort_haps[0][2], $chr_start, $chr_end] );
1580 }
1581 else{
1582 push( @syms, [ $hap_start, $hap_end, $sort_haps[0][2], $chr_start, $chr_end] );
1583 }
1584 next;
1585 }
1586 if($hap_end and $hap_start < $len1){ # if hap at start or end of chromosome
1587 push( @syms, [ $hap_start, $hap_end, $sort_haps[0][2], $chr_start, $chr_end] );
1588 }
1589 $chr_start = $chr_end + ($sort_haps[$count][4]-$sort_haps[$count][3]) + 2;
1590 $hap_start = $hap_end + ($sort_haps[$count][1]-$sort_haps[$count][0]) + 2;
1591 $count++;
1592 }
1593
1594
1595 }
1596
1597
1598 # for now haps and pars should not be both there, but in theory we
1599 # could handle it here by cleverly merging the pars into the existing syms,
1600 # for now just:
1601 push( @syms, @pars );
1602
1603 my $mapper = Bio::EnsEMBL::Mapper->new( "sym", "org" );
1604 my $count = 0;
1605 for my $sym ( @syms ) {
1606 $mapper->add_map_coordinates( $slice_seq_region_id, $sym->[0], $sym->[1],
1607 1, $sym->[2], $sym->[3], $sym->[4] );
1608 }
1609
1610 my @linked = $mapper->map_coordinates( $slice_seq_region_id,
1611 $slice->start(), $slice->end(),
1612 $slice->strand(), "sym" );
1613
1614 # gaps are regions where there is no mapping to another region
1615 my $rel_start = 1;
1616
1617 #if there was only one coord and it is a gap, we know it is just the
1618 #same slice with no overlapping symlinks
1619 if(@linked == 1 && $linked[0]->isa('Bio::EnsEMBL::Mapper::Gap')) {
1620 return [bless( [1,$slice->length, $slice], "Bio::EnsEMBL::ProjectionSegment" )];
1621 }
1622
1623 my @out;
1624 for my $coord ( @linked ) {
1625 if( $coord->isa( "Bio::EnsEMBL::Mapper::Gap" )) {
1626 my $exc_slice = Bio::EnsEMBL::Slice->new_fast({
1627 'start' => $coord->start(),
1628 'end' => $coord->end(),
1629 'strand' => $slice->strand(),
1630 'coord_system' => $slice->coord_system(),
1631 'adaptor' => $self,
1632 'seq_region_name' => $slice->seq_region_name(),
1633 'seq_region_length' => $slice->seq_region_length()});
1634 push( @out, bless ( [ $rel_start, $coord->length()+$rel_start-1,
1635 $exc_slice ], "Bio::EnsEMBL::ProjectionSegment") );
1636 } else {
1637 my $exc_slice = $self->fetch_by_seq_region_id( $coord->id() );
1638 my $exc2_slice = Bio::EnsEMBL::Slice->new_fast({
1639
1640 'start' => $coord->start(),
1641 'end' => $coord->end(),
1642 'strand' => $coord->strand(),
1643 'seq_region_name' => $exc_slice->seq_region_name(),
1644 'seq_region_length' => $exc_slice->seq_region_length(),
1645 'coord_system' => $exc_slice->coord_system(),
1646 'adaptor' => $self
1647 });
1648
1649 push( @out, bless( [ $rel_start, $coord->length() + $rel_start - 1,
1650 $exc2_slice ], "Bio::EnsEMBL::ProjectionSegment") );
1651 }
1652 $rel_start += $coord->length();
1653 }
1654
1655 return \@out;
1656 }
1657
1658
1659
1660
1661 =head2 store
1662
1663 Arg [1] : Bio::EnsEMBL::Slice $slice
1664 Arg [2] : (optional) $seqref reference to a string
1665 The sequence associated with the slice to be stored.
1666 Example : $slice = Bio::EnsEMBL::Slice->new(...);
1667 $seq_region_id = $slice_adaptor->store($slice, \$sequence);
1668 Description: This stores a slice as a sequence region in the database
1669 and returns the seq region id. The passed in slice must
1670 start at 1, and must have a valid seq_region name and coordinate
1671 system. The attached coordinate system must already be stored in
1672 the database. The sequence region is assumed to start at 1 and
1673 to have a length equalling the length of the slice. The end of
1674 the slice must equal the seq_region_length.
1675 If the slice coordinate system is the sequence level coordinate
1676 system then the seqref argument must also be passed. If the
1677 slice coordinate system is NOT a sequence level coordinate
1678 system then the sequence argument cannot be passed.
1679 Returntype : int
1680 Exceptions : throw if slice has no coord system.
1681 throw if slice coord system is not already stored.
1682 throw if slice coord system is seqlevel and no sequence is
1683 provided.
1684 throw if slice coord system is not seqlevel and sequence is
1685 provided.
1686 throw if slice does not start at 1
1687 throw if sequence is provided and the sequence length does not
1688 match the slice length.
1689 throw if the SQL insert fails (e.g. on duplicate seq region)
1690 throw if slice argument is not passed
1691 throw if the slice end is not equal to seq_region_length
1692 Caller : database loading scripts
1693 Status : Stable
1694
1695 =cut
1696
1697
1698
1699 sub store {
1700 my $self = shift;
1701 my $slice = shift;
1702 my $seqref = shift;
1703
1704 #
1705 # Get all of the sanity checks out of the way before storing anything
1706 #
1707
1708 if(!ref($slice) || !($slice->isa('Bio::EnsEMBL::Slice') or $slice->isa('Bio::EnsEMBL::LRGSlice'))) {
1709 throw('Slice argument is required');
1710 }
1711
1712 my $cs = $slice->coord_system();
1713 throw("Slice must have attached CoordSystem.") if(!$cs);
1714
1715 my $db = $self->db();
1716 if(!$cs->is_stored($db)) {
1717 throw("Slice CoordSystem must already be stored in DB.")
1718 }
1719
1720 if($slice->start != 1 || $slice->strand != 1) {
1721 throw("Slice must have start==1 and strand==1.");
1722 }
1723
1724 if($slice->end() != $slice->seq_region_length()) {
1725 throw("Slice must have end==seq_region_length");
1726 }
1727
1728 my $sr_len = $slice->length();
1729 my $sr_name = $slice->seq_region_name();
1730
1731 if(!$sr_name) {
1732 throw("Slice must have valid seq region name.");
1733 }
1734
1735 if($cs->is_sequence_level()) {
1736 if(!$seqref) {
1737 throw("Must provide sequence for sequence level coord system.");
1738 }
1739 if(ref($seqref) ne 'SCALAR') {
1740 throw("Sequence must be a scalar reference.");
1741 }
1742 my $seq_len = length($$seqref);
1743
1744 if($seq_len != $sr_len) {
1745 throw("Sequence length ($seq_len) must match slice length ($sr_len).");
1746 }
1747 } else {
1748 if($seqref) {
1749 throw("Cannot provide sequence for non-sequence level seq regions.");
1750 }
1751 }
1752
1753 #store the seq_region
1754
1755 my $sth = $db->dbc->prepare("INSERT INTO seq_region " .
1756 "SET name = ?, " .
1757 " length = ?, " .
1758 " coord_system_id = ?" );
1759
1760 $sth->bind_param(1,$sr_name,SQL_VARCHAR);
1761 $sth->bind_param(2,$sr_len,SQL_INTEGER);
1762 $sth->bind_param(3,$cs->dbID,SQL_INTEGER);
1763
1764 $sth->execute();
1765
1766 my $seq_region_id = $sth->{'mysql_insertid'};
1767
1768 if(!$seq_region_id) {
1769 throw("Database seq_region insertion failed.");
1770 }
1771
1772 if($cs->is_sequence_level()) {
1773 #store sequence if it was provided
1774 my $seq_adaptor = $db->get_SequenceAdaptor();
1775 $seq_adaptor->store($seq_region_id, $$seqref);
1776 }
1777
1778 #synonyms
1779 if(defined($slice->{'synonym'})){
1780 foreach my $syn (@{$slice->{'synonym'}} ){
1781 $syn->seq_region_id($seq_region_id); # set the seq_region_id
1782 $syn->adaptor->store($syn);
1783 }
1784 }
1785
1786
1787 $slice->adaptor($self);
1788
1789 return $seq_region_id;
1790 }
1791
1792
1793 =head2 store_assembly
1794
1795 Arg [1] : Bio::EnsEMBL::Slice $asm_slice
1796 Arg [2] : Bio::EnsEMBL::Slice $cmp_slice
1797 Example : $asm = $slice_adaptor->store_assembly( $slice1, $slice2 );
1798 Description: Creates an entry in the analysis table based on the
1799 coordinates of the two slices supplied. Returns a string
1800 representation of the assembly that gets created.
1801 Returntype : string
1802 Exceptions : throw if either slice has no coord system (cs).
1803 throw unless the cs rank of the asm_slice is lower than the
1804 cmp_slice.
1805 throw if there is no mapping path between coord systems
1806 throw if the lengths of each slice are not equal
1807 throw if there are existing mappings between either slice
1808 and the oposite cs
1809 Caller : database loading scripts
1810 Status : Experimental
1811
1812 =cut
1813
1814 sub store_assembly{
1815 my $self = shift;
1816 my $asm_slice = shift;
1817 my $cmp_slice = shift;
1818
1819 #
1820 # Get all of the sanity checks out of the way before storing anything
1821 #
1822
1823 if(!ref($asm_slice) || !($asm_slice->isa('Bio::EnsEMBL::Slice') or $asm_slice->isa('Bio::EnsEMBL::LRGSlice'))) {
1824 throw('Assembled Slice argument is required');
1825 }
1826 if(!ref($cmp_slice) || !($cmp_slice->isa('Bio::EnsEMBL::Slice') or $cmp_slice->isa('Bio::EnsEMBL::LRGSlice')) ) {
1827 throw('Assembled Slice argument is required');
1828 }
1829
1830 my $asm_cs = $asm_slice->coord_system();
1831 throw("Slice must have attached CoordSystem.") if(!$asm_cs);
1832 my $cmp_cs = $cmp_slice->coord_system();
1833 throw("Slice must have attached CoordSystem.") if(!$cmp_cs);
1834
1835 unless( $asm_cs->rank < $cmp_cs->rank ){
1836 throw("Assembled Slice CoordSystem->rank must be lower than ".
1837 "the component Slice Coord_system" );
1838 }
1839
1840 my @path =
1841 @{ $asm_cs->adaptor()->get_mapping_path( $asm_cs, $cmp_cs ) };
1842
1843 if ( !@path ) {
1844 throw("No mapping path defined between "
1845 . $asm_cs->name() . " and "
1846 . $cmp_cs->name() );
1847 }
1848
1849 if( $asm_slice->length != $cmp_slice->length ){
1850 throw("The lengths of the assembled and component slices are not equal" );
1851 }
1852
1853 # For now we disallow any existing mappings between the asm slice and cmp
1854 # CoordSystem and vice-versa.
1855 # Some cases of multiple mappings may be allowable by the API, but their
1856 # logic needs to be coded below.
1857
1858 my $asm_proj = $asm_slice->project( $cmp_cs->name, $cmp_cs->version );
1859 if( @$asm_proj ){
1860 throw("Regions of the assembled slice are already assembled ".
1861 "into the component CoordSystem" );
1862 }
1863 my $cmp_proj = $cmp_slice->project( $asm_cs->name, $asm_cs->version );
1864 if( @$cmp_proj ){
1865 throw("Regions of the component slice are already assembled ".
1866 "into the assembled CoordSystem" );
1867 }
1868
1869 #
1870 # Checks complete. Store the data
1871 #
1872 my $sth = $self->db->dbc->prepare
1873 ("INSERT INTO assembly " .
1874 "SET asm_seq_region_id = ?, " .
1875 " cmp_seq_region_id = ?, " .
1876 " asm_start = ?, " .
1877 " asm_end = ?, " .
1878 " cmp_start = ?, " .
1879 " cmp_end = ?, " .
1880 " ori = ?" );
1881
1882 my $asm_seq_region_id = $self->get_seq_region_id( $asm_slice );
1883 my $cmp_seq_region_id = $self->get_seq_region_id( $cmp_slice );
1884 my $ori = $asm_slice->strand * $cmp_slice->strand;
1885
1886 $sth->bind_param(1,$asm_seq_region_id,SQL_INTEGER);
1887 $sth->bind_param(2,$cmp_seq_region_id,SQL_INTEGER);
1888 $sth->bind_param(3,$asm_slice->start,SQL_INTEGER);
1889 $sth->bind_param(4,$asm_slice->end,SQL_INTEGER);
1890 $sth->bind_param(5,$cmp_slice->start,SQL_INTEGER);
1891 $sth->bind_param(6,$cmp_slice->end,SQL_INTEGER);
1892 $sth->bind_param(7,$ori,SQL_INTEGER);
1893
1894 $sth->execute();
1895
1896 #use Data::Dumper qw( Dumper );
1897 #warn Dumper( $self->db->{seq_region_cache} );
1898 #$self->db->{seq_region_cache} = undef;
1899 #$self->_cache_seq_regions();
1900
1901 my $ama = $self->db->get_AssemblyMapperAdaptor();
1902 $ama->delete_cache();
1903
1904
1905 return $asm_slice->name . "<>" . $cmp_slice->name;
1906
1907 }
1908
1909
1910 =head2 prepare
1911
1912 Arg [1] : String $sql
1913 Example : ( optional )
1914 Description: overrides the default adaptor prepare method.
1915 All slice sql will usually use the dna_db.
1916 Returntype : DBD::sth
1917 Exceptions : none
1918 Caller : internal, convenience method
1919 Status : Stable
1920
1921 =cut
1922
1923 sub prepare {
1924 my ( $self, $sql ) = @_;
1925 return $self->db()->dnadb()->dbc->prepare($sql);
1926 }
1927
1928 sub _build_exception_cache {
1929 my $self = shift;
1930
1931 # build up a cache of the entire assembly exception table
1932 # it should be small anyway
1933 my $sth =
1934 $self->prepare( 'SELECT ae.seq_region_id, ae.seq_region_start, '
1935 . 'ae.seq_region_end, ae.exc_type, ae.exc_seq_region_id, '
1936 . 'ae.exc_seq_region_start, ae.exc_seq_region_end '
1937 . 'FROM assembly_exception ae, '
1938 . 'seq_region sr, coord_system cs '
1939 . 'WHERE sr.seq_region_id = ae.seq_region_id '
1940 . 'AND sr.coord_system_id = cs.coord_system_id '
1941 . 'AND cs.species_id = ?' );
1942
1943 $sth->bind_param( 1, $self->species_id(), SQL_INTEGER );
1944 $sth->execute();
1945
1946 my %hash;
1947 $self->{'asm_exc_cache'} = \%hash;
1948
1949 my $row;
1950 while ( $row = $sth->fetchrow_arrayref() ) {
1951 my @result = @$row;
1952 $hash{ $result[0] } ||= [];
1953 push( @{ $hash{ $result[0] } }, \@result );
1954 }
1955 $sth->finish();
1956 } ## end sub _build_exception_cache
1957
1958 =head2 cache_toplevel_seq_mappings
1959
1960 Args : none
1961 Example : $slice_adaptor->cache_toplevel_seq_mappings();
1962 Description: caches all the assembly mappings needed for genes
1963 Returntype : None
1964 Exceptions : None
1965 Caller : general
1966 Status : At Risk
1967 : New experimental code
1968
1969 =cut
1970
1971 sub cache_toplevel_seq_mappings {
1972 my ($self) = @_;
1973
1974 # Get the sequence level to map too
1975
1976 my $sql = (<<SSQL);
1977 SELECT name
1978 FROM coord_system
1979 WHERE attrib like "%sequence_level%"
1980 AND species_id = ?
1981 SSQL
1982
1983 my $sth = $self->prepare($sql);
1984 $sth->bind_param( 1, $self->species_id(), SQL_INTEGER );
1985 $sth->execute();
1986
1987 my $sequence_level = $sth->fetchrow_array();
1988
1989 $sth->finish();
1990
1991 my $csa = $self->db->get_CoordSystemAdaptor();
1992 my $ama = $self->db->get_AssemblyMapperAdaptor();
1993
1994 my $cs1 = $csa->fetch_by_name($sequence_level);
1995
1996 #get level to map too.
1997
1998 $sql = (<<LSQL);
1999 SELECT DISTINCT(cs.name)
2000 FROM seq_region sr,
2001 seq_region_attrib sra,
2002 attrib_type at,
2003 coord_system cs
2004 WHERE sra.seq_region_id = sr.seq_region_id
2005 AND sra.attrib_type_id = at.attrib_type_id
2006 AND at.code = "toplevel"
2007 AND cs.coord_system_id = sr.coord_system_id
2008 AND cs.species_id = ?
2009 LSQL
2010
2011 $sth = $self->prepare($sql);
2012 $sth->bind_param( 1, $self->species_id(), SQL_INTEGER );
2013 $sth->execute();
2014
2015 while ( my $csn = $sth->fetchrow_array() ) {
2016 if ( $csn eq $sequence_level ) { next }
2017 my $cs2 = $csa->fetch_by_name($csn);
2018 my $am = $ama->fetch_by_CoordSystems( $cs1, $cs2 );
2019 $am->register_all();
2020 }
2021
2022 } ## end sub cache_toplevel_seq_mappings
2023
2024
2025 sub _build_circular_slice_cache {
2026 my $self = shift;
2027
2028 # build up a cache of circular sequence region ids
2029 my $sth =
2030 $self->prepare( "SELECT sra.seq_region_id FROM seq_region_attrib sra "
2031 . "INNER JOIN attrib_type at ON sra.attrib_type_id = at.attrib_type_id "
2032 . "INNER JOIN seq_region sr ON sra.seq_region_id = sr.seq_region_id "
2033 . "INNER JOIN coord_system cs ON sr.coord_system_id = cs.coord_system_id "
2034 . "WHERE code = 'circular_seq' and cs.species_id = ?");
2035
2036 $sth->bind_param( 1, $self->species_id(), SQL_INTEGER );
2037 $sth->execute();
2038
2039 my $id;
2040 my %hash;
2041 if ( ($id) = $sth->fetchrow_array() ) {
2042 $self->{'circular_sr_id_cache'} = \%hash;
2043 $self->{'is_circular'} = 1;
2044 $hash{ $id } = $id;
2045 while ( ($id) = $sth->fetchrow_array() ) {
2046 $hash{ $id } = $id;
2047 }
2048 } else {
2049 $self->{'is_circular'} = 0;
2050 }
2051 $sth->finish();
2052 } ## end _build_circular_slice_cache
2053
2054
2055 #####################################
2056 # sub DEPRECATED METHODs
2057 #####################################
2058
2059 =head2 fetch_by_mapfrag
2060
2061 Function: DEPRECATED use fetch_by_misc_feature_attribute('synonym',$mapfrag)
2062
2063 =cut
2064
2065 sub fetch_by_mapfrag{
2066 my ($self,$mymapfrag,$flag,$size) = @_;
2067 deprecate('Use fetch_by_misc_feature_attribute instead');
2068 $flag ||= 'fixed-width'; # alt.. 'context'
2069 $size ||= $flag eq 'fixed-width' ? 100000 : 0;
2070 return $self->fetch_by_misc_feature_attribute('synonym',$mymapfrag,$size);
2071 }
2072
2073
2074
2075 =head2 fetch_by_chr_start_end
2076
2077 Description: DEPRECATED use fetch_by_region instead
2078
2079 =cut
2080
2081 sub fetch_by_chr_start_end {
2082 my ($self,$chr,$start,$end) = @_;
2083 deprecate('Use fetch_by_region() instead');
2084
2085 #assume that by chromosome the user actually meant top-level coord
2086 #system since this is the old behaviour of this deprecated method
2087 my $csa = $self->db->get_CoordSystemAdaptor();
2088 my ($cs) = @{$csa->fetch_all()}; # get the highest coord system
2089
2090 return $self->fetch_by_region($cs->name,$chr,$start,$end,1,$cs->version);
2091 }
2092
2093
2094
2095 =head2 fetch_by_contig_name
2096
2097 Description: Deprecated. Use fetch_by_region(), Slice::project(),
2098 Slice::expand() instead
2099
2100 =cut
2101
2102 sub fetch_by_contig_name {
2103 my ($self, $name, $size) = @_;
2104
2105 deprecate('Use fetch_by_region(), Slice::project() and Slice::expand().');
2106
2107 #previously wanted chromosomal slice on a given contig. Assume this means
2108 #a top-level slice on a given seq_region in the seq_level coord system
2109 my $csa = $self->db()->get_CoordSystemAdaptor();
2110 my $seq_level = $csa->fetch_sequence_level();
2111
2112 my $seq_lvl_slice = $self->fetch_by_region($seq_level->name(), $name);
2113
2114 if(!$seq_lvl_slice) {
2115 return undef;
2116 }
2117
2118 my @projection = @{$seq_lvl_slice->project('toplevel')};
2119
2120 if(@projection != 1) {
2121 warning("$name is mapped to multiple toplevel locations.");
2122 }
2123
2124 return $projection[0]->[2]->expand($size, $size);
2125 }
2126
2127
2128 =head2 fetch_by_clone_accession
2129
2130 Description: DEPRECATED. Use fetch_by_region, Slice::project, Slice::expand
2131 instead.
2132
2133 =cut
2134
2135 sub fetch_by_clone_accession{
2136 my ($self,$name,$size) = @_;
2137
2138 deprecate('Use fetch_by_region(), Slice::project() and Slice::expand().');
2139
2140 my $csa = $self->db()->get_CoordSystemAdaptor();
2141 my $clone_cs = $csa->fetch_by_name('clone');
2142
2143 if(!$clone_cs) {
2144 warning('Clone coordinate system does not exist for this species');
2145 return undef;
2146 }
2147
2148 #this unfortunately needs a version on the end to work
2149 if(! ($name =~ /\./)) {
2150 my $sth =
2151 $self->prepare( "SELECT sr.name "
2152 . "FROM seq_region sr, coord_system cs "
2153 . "WHERE cs.name = 'clone' "
2154 . "AND cs.coord_system_id = sr.coord_system_id "
2155 . "AND sr.name LIKE '$name.%'"
2156 . "AND cs.species_id = ?" );
2157
2158 $sth->bind_param( 1, $self->species_id(), SQL_INTEGER );
2159 $sth->execute();
2160
2161 if(!$sth->rows()) {
2162 $sth->finish();
2163 throw("Clone $name not found in database");
2164 }
2165
2166 ($name) = $sth->fetchrow_array();
2167
2168 $sth->finish();
2169 }
2170
2171 my $clone = $self->fetch_by_region($clone_cs->name(), $name);
2172 return undef if(!$clone);
2173
2174 my @projection = @{$clone->project('toplevel')};
2175
2176 if(@projection != 1) {
2177 warning("$name is mapped to multiple locations.");
2178 }
2179
2180 return $projection[0]->[2]->expand($size, $size);
2181 }
2182
2183
2184 =head2 fetch_by_supercontig_name
2185
2186 Description: DEPRECATED. Use fetch_by_region(), Slice::project() and
2187 Slice::expand() instead
2188
2189 =cut
2190
2191 sub fetch_by_supercontig_name {
2192 my ($self,$name, $size) = @_;
2193
2194 deprecate('Use fetch_by_region(), Slice::project() and Slice::expand().');
2195
2196 my $csa = $self->db()->get_CoordSystemAdaptor();
2197 my $sc_level = $csa->fetch_by_name('supercontig');
2198
2199 if(!$sc_level) {
2200 warning('No supercontig coordinate system exists for this species.');
2201 return undef;
2202 }
2203
2204 my $sc_slice = $self->fetch_by_region($sc_level->name(),$name);
2205
2206 return undef if(!$sc_slice);
2207
2208 my @projection = @{$sc_slice->project('toplevel')};
2209
2210 if(@projection > 1) {
2211 warning("$name is mapped to multiple locations in toplevel");
2212 }
2213
2214 return $projection[0]->[2]->expand($size, $size);
2215 }
2216
2217
2218
2219
2220 =head2 list_overlapping_supercontigs
2221
2222 Description: DEPRECATED use Slice::project instead
2223
2224 =cut
2225
2226 sub list_overlapping_supercontigs {
2227 my ($self,$slice) = @_;
2228
2229 deprecate('Use Slice::project() instead.');
2230
2231 my $csa = $self->db()->get_CoordSystemAdaptor();
2232 my $sc_level = $csa->fetch_by_name('supercontig');
2233
2234 if(!$sc_level) {
2235 warning('No supercontig coordinate system exists for this species.');
2236 return undef;
2237 }
2238
2239 my @out;
2240 foreach my $seg (@{$slice->project($sc_level->name(), $sc_level->version)}){
2241 push @out, $seg->[2]->seq_region_name();
2242 }
2243
2244 return \@out;
2245 }
2246
2247
2248
2249 =head2 fetch_by_chr_name
2250
2251 Description: DEPRECATED. Use fetch by region instead
2252
2253 =cut
2254
2255 sub fetch_by_chr_name{
2256 my ($self,$chr_name) = @_;
2257 deprecate('Use fetch_by_region() instead.');
2258
2259 my $csa = $self->db->get_CoordSystemAdaptor();
2260
2261 my $top_cs = @{$csa->fetch_all()};
2262
2263 return $self->fetch_by_region($top_cs->name(),$chr_name,
2264 undef,undef,undef,$top_cs->version);
2265 }
2266
2267
2268
2269
2270
2271
2272 1;