Mercurial > repos > mahtabm > ensembl

comparison variant_effect_predictor/Bio/EnsEMBL/DBSQL/BaseFeatureAdaptor.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::BaseFeatureAdaptor - An Abstract Base class for all
24 FeatureAdaptors
25
26 =head1 SYNOPSIS
27
28 Abstract class - should not be instantiated. Implementation of
29 abstract methods must be performed by subclasses.
30
31 =head1 DESCRIPTION
32
33 This is a base adaptor for feature adaptors. This base class is simply a way
34 of eliminating code duplication through the implementation of methods
35 common to all feature adaptors.
36
37 =head1 METHODS
38
39 =cut
40
41 package Bio::EnsEMBL::DBSQL::BaseFeatureAdaptor;
42 use vars qw(@ISA @EXPORT);
43 use strict;
44
45 use Bio::EnsEMBL::DBSQL::BaseAdaptor;
46 use Bio::EnsEMBL::Utils::Cache;
47 use Bio::EnsEMBL::Utils::Exception qw(warning throw deprecate stack_trace_dump);
48 use Bio::EnsEMBL::Utils::Argument qw(rearrange);
49 use Bio::EnsEMBL::Utils::Iterator;
50
51 @ISA = qw(Bio::EnsEMBL::DBSQL::BaseAdaptor);
52
53 @EXPORT = (@{$DBI::EXPORT_TAGS{'sql_types'}});
54
55 our $SLICE_FEATURE_CACHE_SIZE = 4;
56 our $MAX_SPLIT_QUERY_SEQ_REGIONS = 3;
57 our $SILENCE_CACHE_WARNINGS = 0;
58
59 =head2 new
60
61 Arg [1] : list of args @args
62 Superclass constructor arguments
63 Example : none
64 Description: Constructor which warns if caching has been switched off
65 Returntype : Bio::EnsEMBL::BaseFeatureAdaptor
66 Exceptions : none
67 Caller : implementing subclass constructors
68 Status : Stable
69
70 =cut
71
72 sub new {
73 my ($class, @args) = @_;
74 my $self = $class->SUPER::new(@args);
75 if ( defined $self->db->no_cache() && $self->db->no_cache() && ! $SILENCE_CACHE_WARNINGS) {
76 warning( "You are using the API without caching most recent features. "
77 . "Performance might be affected." );
78 }
79 return $self;
80 }
81
82 =head2 start_equals_end
83
84 Arg [1] : (optional) boolean $newval
85 Example : $bfa->start_equals_end(1);
86 Description: Getter/Setter for the start_equals_end flag. If set
87 to true sub _slice_fetch will use a simplified sql to retrieve 1bp slices.
88 Returntype : boolean
89 Exceptions : none
90 Caller : EnsemblGenomes variation DB build
91 Status : Stable
92
93 =cut
94
95 sub start_equals_end {
96 my ( $self, $value ) = @_;
97
98 if ( defined($value) ) {
99 $self->{'start_equals_end'} = $value;
100 }
101 return $self->{'start_equals_end'};
102 }
103
104
105 =head2 clear_cache
106
107 Args : None
108 Example : my $sa =
109 $registry->get_adaptor( 'Mus musculus', 'Core',
110 'Slice' );
111 my $ga =
112 $registry->get_adaptor( 'Mus musculus', 'Core',
113 'Gene' );
114
115 my $slice =
116 $sa->fetch_by_region( 'Chromosome', '1', 1e8,
117 1.05e8 );
118
119 my $genes = $ga->fetch_all_by_Slice($slice);
120
121 $ga->clear_cache();
122
123 Description : Empties the feature cache associated with this
124 feature adaptor.
125 Return type : None
126 Exceptions : None
127 Caller : General
128 Status : At risk (under development)
129
130 =cut
131
132 sub clear_cache {
133 my ($self) = @_;
134 %{$self->{_slice_feature_cache}} = ();
135 return;
136 }
137
138 =head2 _slice_feature_cache
139
140 Description : Returns the feature cache if we are allowed to cache and
141 will build it if we need to. We will never return a reference
142 to the hash to avoid unintentional auto-vivfying caching
143 Returntype : Bio::EnsEMBL::Utils::Cache
144 Exceptions : None
145 Caller : Internal
146
147 =cut
148
149 sub _slice_feature_cache {
150 my ($self) = @_;
151 return if $self->db()->no_cache();
152 if(! exists $self->{_slice_feature_cache}) {
153 tie my %cache, 'Bio::EnsEMBL::Utils::Cache', $SLICE_FEATURE_CACHE_SIZE;
154 $self->{_slice_feature_cache} = \%cache;
155 }
156 return $self->{_slice_feature_cache};
157 }
158
159 =head2 fetch_all_by_Slice
160
161 Arg [1] : Bio::EnsEMBL::Slice $slice
162 the slice from which to obtain features
163 Arg [2] : (optional) string $logic_name
164 the logic name of the type of features to obtain
165 Example : $fts = $a->fetch_all_by_Slice($slice, 'Swall');
166 Description: Returns a listref of features created from the database
167 which are on the Slice defined by $slice. If $logic_name is
168 defined only features with an analysis of type $logic_name
169 will be returned.
170 NOTE: only features that are entirely on the slice's seq_region
171 will be returned (i.e. if they hang off the start/end of a
172 seq_region they will be discarded). Features can extend over the
173 slice boundaries though (in cases where you have a slice that
174 doesn't span the whole seq_region).
175 Returntype : listref of Bio::EnsEMBL::SeqFeatures in Slice coordinates
176 Exceptions : none
177 Caller : Bio::EnsEMBL::Slice
178 Status : Stable
179
180 =cut
181
182 sub fetch_all_by_Slice {
183 my ($self, $slice, $logic_name) = @_;
184 #fetch by constraint with empty constraint
185 return $self->fetch_all_by_Slice_constraint($slice, '', $logic_name);
186 }
187
188
189
190 =head2 fetch_Iterator_by_Slice_method
191
192 Arg [1] : CODE ref of Slice fetch method
193 Arg [2] : ARRAY ref of parameters for Slice fetch method
194 Arg [3] : Optional int: Slice index in parameters array
195 Arg [4] : Optional int: Slice chunk size. Default=500000
196 Example : my $slice_iter = $feature_adaptor->fetch_Iterator_by_Slice_method
197 ($feature_adaptor->can('fetch_all_by_Slice_Arrays'),
198 \@fetch_method_params,
199 0,#Slice idx
200 );
201
202 while(my $feature = $slice_iter->next && defined $feature){
203 #Do something here
204 }
205
206 Description: Creates an Iterator which chunks the query Slice to facilitate
207 large Slice queries which would have previously run out of memory
208 Returntype : Bio::EnsEMBL::Utils::Iterator
209 Exceptions : Throws if mandatory params not valid
210 Caller : general
211 Status : at risk
212
213 =cut
214
215 #Does not support Collections. See Funcgen ResultFeatureAdaptor::fetch_collection_Iterator_by_Slice_method
216
217 sub fetch_Iterator_by_Slice_method{
218 my ($self, $slice_method_ref, $params_ref, $slice_idx, $chunk_size) = @_;
219
220 if(! ( defined $slice_method_ref &&
221 ref($slice_method_ref) eq 'CODE')
222 ){
223 throw('Must pass a valid Slice fetch method CODE ref');
224 }
225
226 if (! ($params_ref &&
227 ref($params_ref) eq 'ARRAY')) {
228 #Don't need to check size here so long as we have valid Slice
229 throw('You must pass a method params ARRAYREF');
230 }
231
232 $slice_idx = 0 if(! defined $slice_idx);
233 my $slice = $params_ref->[$slice_idx];
234 $chunk_size ||= 1000000;
235
236 my @feat_cache;
237 my $finished = 0;
238 my $start = 1; #local coord for sub slice
239 my $end = $slice->length;
240 my $num_overlaps = 0;
241
242 my $coderef =
243 sub {
244
245 while (scalar(@feat_cache) == 0 &&
246 ! $finished) {
247
248 my $new_end = ($start + $chunk_size - 1);
249
250 if ($new_end >= $end) {
251 # this is our last chunk
252 $new_end = $end;
253 $finished = 1;
254 }
255
256 #Chunk by sub slicing
257 my $sub_slice = $slice->sub_Slice($start, $new_end);
258 $params_ref->[$slice_idx] = $sub_slice;
259 @feat_cache = @{ $slice_method_ref->($self, @$params_ref)};
260
261 #Remove & count overlapping features
262 splice(@feat_cache, 0, $num_overlaps) if($num_overlaps);
263 my $i;
264
265 if (scalar(@feat_cache) > 0) {
266
267 my $feat_end = $feat_cache[$#feat_cache]->seq_region_end;
268 my $slice_end = $sub_slice->end;
269 $num_overlaps = 0;
270
271 for ($i = $#feat_cache; $i >=0; $i--) {
272
273 if ($feat_end > $slice_end) {
274 $feat_end = $feat_cache[$i]->end;
275 $num_overlaps ++;
276 } else {
277 last;
278 }
279
280 }
281 }
282
283 # update the start coordinate
284 $start = $new_end + 1;
285 }
286
287 #this maybe returning from an undef cache
288 #Need to sub this out even more?
289 return shift @feat_cache;
290 };
291
292 return Bio::EnsEMBL::Utils::Iterator->new($coderef);
293 }
294
295
296 =head2 fetch_Iterator_by_Slice
297
298 Arg [1] : Bio::EnsEMBL::Slice
299 Arg [2] : Optional string: logic name of analysis
300 Arg [3] : Optional int: Chunk size to iterate over. Default is 500000
301 Example : my $slice_iter = $feature_adaptor->fetch_Iterator_by_Slice($slice);
302
303 while(my $feature = $slice_iter->next && defined $feature){
304 #Do something here
305 }
306
307 Description: Creates an Iterator which chunks the query Slice to facilitate
308 large Slice queries which would have previously run out of memory
309 Returntype : Bio::EnsEMBL::Utils::Iterator
310 Exceptions : None
311 Caller : general
312 Status : at risk
313
314 =cut
315
316 sub fetch_Iterator_by_Slice{
317 my ($self, $slice, $logic_name, $chunk_size) = @_;
318
319 my $method_ref = $self->can('fetch_all_by_Slice');
320
321 return $self->fetch_Iterator_by_Slice_method($method_ref, [$slice, $logic_name], 0, $chunk_size);
322 }
323
324
325 =head2 fetch_all_by_Slice_and_score
326
327 Arg [1] : Bio::EnsEMBL::Slice $slice
328 the slice from which to obtain features
329 Arg [2] : (optional) float $score
330 lower bound of the the score of the features retrieved
331 Arg [3] : (optional) string $logic_name
332 the logic name of the type of features to obtain
333 Example : $fts = $a->fetch_all_by_Slice_and_score($slice,90,'Swall');
334 Description: Returns a list of features created from the database which are
335 are on the Slice defined by $slice and which have a score
336 greater than $score. If $logic_name is defined,
337 only features with an analysis of type $logic_name will be
338 returned.
339 Returntype : listref of Bio::EnsEMBL::SeqFeatures in Slice coordinates
340 Exceptions : none
341 Caller : Bio::EnsEMBL::Slice
342 Status : Stable
343
344 =cut
345
346 sub fetch_all_by_Slice_and_score {
347 my ( $self, $slice, $score, $logic_name ) = @_;
348
349 my $constraint;
350 if ( defined($score) ) {
351 # Get the synonym of the primary_table
352 my @tabs = $self->_tables();
353 my $syn = $tabs[0]->[1];
354
355 $constraint = sprintf( "%s.score > %s",
356 $syn,
357 $self->dbc()->db_handle()->quote( $score, SQL_FLOAT ) );
358 }
359
360 return
361 $self->fetch_all_by_Slice_constraint( $slice, $constraint,
362 $logic_name );
363 }
364
365
366 =head2 fetch_all_by_Slice_constraint
367
368 Arg [1] : Bio::EnsEMBL::Slice $slice
369 the slice from which to obtain features
370 Arg [2] : (optional) string $constraint
371 An SQL query constraint (i.e. part of the WHERE clause)
372 Arg [3] : (optional) string $logic_name
373 the logic name of the type of features to obtain
374 Example : $fs = $a->fetch_all_by_Slice_constraint($slc, 'perc_ident > 5');
375 Description: Returns a listref of features created from the database which
376 are on the Slice defined by $slice and fulfill the SQL
377 constraint defined by $constraint. If logic name is defined,
378 only features with an analysis of type $logic_name will be
379 returned.
380 Returntype : listref of Bio::EnsEMBL::SeqFeatures in Slice coordinates
381 Exceptions : thrown if $slice is not defined
382 Caller : Bio::EnsEMBL::Slice
383 Status : Stable
384
385 =cut
386
387 sub fetch_all_by_Slice_constraint {
388 my ( $self, $slice, $constraint, $logic_name ) = @_;
389
390
391 my @result = ();
392
393 if ( !ref($slice)
394 || !( $slice->isa('Bio::EnsEMBL::Slice')
395 or $slice->isa('Bio::EnsEMBL::LRGSlice') ) )
396 {
397 throw("Bio::EnsEMBL::Slice argument expected.");
398 }
399
400 $constraint ||= '';
401 $constraint =
402 $self->_logic_name_to_constraint( $constraint, $logic_name );
403
404 # If the logic name was invalid, undef was returned
405 if ( !defined($constraint) ) { return [] }
406
407 my $key;
408 my $cache;
409
410 # Will only use feature_cache if hasn't been no_cache attribute set
411 if (
412 !( defined( $self->db()->no_cache() ) && $self->db()->no_cache() ) )
413 {
414
415 #strain test and add to constraint if so to stop caching.
416 if ( $slice->isa('Bio::EnsEMBL::StrainSlice') ) {
417 my $string =
418 $self->dbc()->db_handle()->quote( $slice->strain_name() );
419
420 if ( $constraint ne "" ) {
421 $constraint .= " AND $string = $string ";
422 } else {
423 $constraint .= " $string = $string ";
424 }
425 }
426
427 # Check the cache and return the cached results if we have already
428 # done this query. The cache key is the made up from the slice
429 # name, the constraint, and the bound parameters (if there are any).
430 $key = uc( join( ':', $slice->name(), $constraint ) );
431
432 my $bind_params = $self->bind_param_generic_fetch();
433
434 if ( defined($bind_params) ) {
435 $key .= ':'
436 . join( ':', map { $_->[0] . '/' . $_->[1] } @{$bind_params} );
437 }
438
439 $cache = $self->_slice_feature_cache();
440 if ( exists( $cache->{$key} ) ) {
441 # Clear the bound parameters and return the cached data.
442 $self->{'_bind_param_generic_fetch'} = ();
443 return $cache->{$key};
444 }
445 } ## end if ( !( defined( $self...)))
446
447 my $sa = $slice->adaptor();
448
449 # Hap/PAR support: retrieve normalized 'non-symlinked' slices.
450 my @proj = @{ $sa->fetch_normalized_slice_projection($slice) };
451
452
453
454 if ( !@proj ) {
455 throw( 'Could not retrieve normalized Slices. '
456 . 'Database contains incorrect assembly_exception information.'
457 );
458 }
459
460 # Want to get features on the FULL original slice as well as any
461 # symlinked slices.
462
463 # Filter out partial slices from projection that are on same
464 # seq_region as original slice.
465
466 my $sr_id = $slice->get_seq_region_id();
467
468 @proj = grep { $_->to_Slice->get_seq_region_id() != $sr_id } @proj;
469
470 my $segment = bless( [ 1, $slice->length(), $slice ],
471 'Bio::EnsEMBL::ProjectionSegment' );
472 push( @proj, $segment );
473
474 # construct list of Hap/PAR boundaries for entire seq region
475 my @bounds;
476
477 my $ent_slice = $sa->fetch_by_seq_region_id($sr_id);
478 if ( $slice->strand() == -1 ) {
479 $ent_slice = $ent_slice->invert();
480 }
481
482 my @ent_proj =
483 @{ $sa->fetch_normalized_slice_projection($ent_slice) };
484 shift(@ent_proj); # skip first
485
486 @bounds = map { $_->from_start() - $slice->start() + 1 } @ent_proj;
487
488
489 # fetch features for the primary slice AND all symlinked slices
490 foreach my $seg (@proj) {
491
492
493 my $offset = $seg->from_start();
494 my $seg_slice = $seg->to_Slice();
495 my $features =
496 $self->_slice_fetch( $seg_slice, $constraint );
497
498 # If this was a symlinked slice offset the feature coordinates as
499 # needed.
500 if ( $seg_slice->name() ne $slice->name() ) {
501
502 FEATURE:
503 foreach my $f ( @{$features} ) {
504 if ( $offset != 1 ) {
505 $f->{'start'} += $offset - 1;
506 $f->{'end'} += $offset - 1;
507 }
508
509 # discard boundary crossing features from symlinked regions
510 foreach my $bound (@bounds) {
511 if ( $f->{'start'} < $bound && $f->{'end'} >= $bound ) {
512 next FEATURE;
513 }
514 }
515
516 $f->{'slice'} = $slice;
517 push( @result, $f );
518 }
519 } else {
520 push( @result, @{$features} );
521 }
522 } ## end foreach my $seg (@proj)
523
524 # Will only use feature_cache when set attribute no_cache in DBAdaptor
525 if ( defined($key) ) {
526 $cache->{$key} = \@result;
527 }
528
529 return \@result;
530 } ## end sub fetch_all_by_Slice_constraint
531
532
533 =head2 fetch_all_by_logic_name
534
535 Arg [3] : string $logic_name
536 the logic name of the type of features to obtain
537 Example : $fs = $a->fetch_all_by_logic_name('foobar');
538 Description: Returns a listref of features created from the database.
539 only features with an analysis of type $logic_name will
540 be returned. If the logic name is invalid (not in the
541 analysis table), a reference to an empty list will be
542 returned.
543 Returntype : listref of Bio::EnsEMBL::SeqFeatures
544 Exceptions : thrown if no $logic_name
545 Caller : General
546 Status : Stable
547
548 =cut
549
550 sub fetch_all_by_logic_name {
551 my ( $self, $logic_name ) = @_;
552
553 if ( !defined($logic_name) ) {
554 throw("Need a logic_name");
555 }
556
557 my $constraint = $self->_logic_name_to_constraint( '', $logic_name );
558
559 if ( !defined($constraint) ) {
560 warning("Invalid logic name: $logic_name");
561 return [];
562 }
563
564 return $self->generic_fetch($constraint);
565 }
566
567 # Method that creates an object. Called by the _objs_from_sth() method
568 # in the sub-classes (the various feature adaptors). Overridden by the
569 # feature collection classes.
570
571 sub _create_feature {
572 my ( $self, $feature_type, $args ) = @_;
573 return $feature_type->new( %{$args} );
574 }
575
576 # This is the same as the above, but calls the new_fast() constructor of
577 # the feature type.
578
579 sub _create_feature_fast {
580 my ( $self, $feature_type, $args ) = @_;
581 return $feature_type->new_fast($args);
582 }
583
584 #
585 # helper function used by fetch_all_by_Slice_constraint method
586 #
587 sub _slice_fetch {
588 my ( $self, $slice, $orig_constraint ) = @_;
589
590 my $slice_start = $slice->start();
591 my $slice_end = $slice->end();
592 my $slice_strand = $slice->strand();
593 my $slice_cs = $slice->coord_system();
594 my $slice_seq_region = $slice->seq_region_name();
595 my $slice_seq_region_id = $slice->get_seq_region_id();
596
597 #get the synonym and name of the primary_table
598 my @tabs = $self->_tables;
599 my ( $tab_name, $tab_syn ) = @{ $tabs[0] };
600
601 #find out what coordinate systems the features are in
602 my $mcc = $self->db->get_MetaCoordContainer();
603 my @feat_css = ();
604
605 my $mca = $self->db->get_MetaContainer();
606 my $value_list = $mca->list_value_by_key( $tab_name . "build.level" );
607 if ( @$value_list and $slice->is_toplevel() ) {
608 push @feat_css, $slice_cs;
609 } else {
610 @feat_css =
611 @{ $mcc->fetch_all_CoordSystems_by_feature_type($tab_name) };
612 }
613
614 my $asma = $self->db->get_AssemblyMapperAdaptor();
615 my @features;
616
617 # fetch the features from each coordinate system they are stored in
618 COORD_SYSTEM: foreach my $feat_cs (@feat_css) {
619 my $mapper;
620 my @coords;
621 my @ids;
622
623 if ( $feat_cs->equals($slice_cs) ) {
624 # no mapping is required if this is the same coord system
625
626 my $max_len = $self->_max_feature_length()
627 || $mcc->fetch_max_length_by_CoordSystem_feature_type( $feat_cs,
628 $tab_name );
629
630 my $constraint = $orig_constraint;
631
632 my $sr_id;
633 if ( $slice->adaptor() ) {
634 $sr_id = $slice->adaptor()->get_seq_region_id($slice);
635 } else {
636 $sr_id =
637 $self->db()->get_SliceAdaptor()->get_seq_region_id($slice);
638 }
639
640 # If there is mapping information, use the external_seq_region_id
641 # to get features.
642
643 my @sr_ids = ($sr_id);
644
645 while (1) {
646 my $ext_sr_id = $self->get_seq_region_id_external($sr_id);
647
648 if ( $ext_sr_id == $sr_id ) { last }
649
650 push( @sr_ids, $ext_sr_id );
651 $sr_id = $ext_sr_id;
652 }
653
654 $constraint .= " AND " if ($constraint);
655
656
657 $constraint .= "${tab_syn}.seq_region_id IN ("
658 . join( ',', @sr_ids ) . ") AND";
659
660 #faster query for 1bp slices where SNP data is not compressed
661 if ( $self->start_equals_end && $slice_start == $slice_end ) {
662 $constraint .=
663 " AND ${tab_syn}.seq_region_start = $slice_end" .
664 " AND ${tab_syn}.seq_region_end = $slice_start";
665
666 } else {
667
668 if ( !$slice->is_circular() ) {
669 # Deal with the default case of a non-circular chromosome.
670 $constraint .= " ${tab_syn}.seq_region_start <= $slice_end AND "
671 . "${tab_syn}.seq_region_end >= $slice_start";
672
673 if ( $max_len ) {
674 my $min_start = $slice_start - $max_len;
675 $constraint .= " AND ${tab_syn}.seq_region_start >= $min_start";
676 }
677
678 } else {
679 # Deal with the case of a circular chromosome.
680 if ( $slice_start > $slice_end ) {
681 $constraint .= " ( ${tab_syn}.seq_region_start >= $slice_start "
682 . "OR ${tab_syn}.seq_region_start <= $slice_end "
683 . "OR ${tab_syn}.seq_region_end >= $slice_start "
684 . "OR ${tab_syn}.seq_region_end <= $slice_end "
685 . "OR ${tab_syn}.seq_region_start > ${tab_syn}.seq_region_end)";
686
687 } else {
688 $constraint .= " ((${tab_syn}.seq_region_start <= $slice_end "
689 . "AND ${tab_syn}.seq_region_end >= $slice_start) "
690 . "OR (${tab_syn}.seq_region_start > ${tab_syn}.seq_region_end "
691 . "AND (${tab_syn}.seq_region_start <= $slice_end "
692 . "OR ${tab_syn}.seq_region_end >= $slice_start)))";
693 }
694 }
695
696 }
697
698 my $fs = $self->generic_fetch( $constraint, undef, $slice );
699
700 # features may still have to have coordinates made relative to slice
701 # start
702 $fs = $self->_remap( $fs, $mapper, $slice );
703
704 push @features, @$fs;
705 } else {
706 $mapper = $asma->fetch_by_CoordSystems( $slice_cs, $feat_cs );
707
708 next unless defined $mapper;
709
710 # Get list of coordinates and corresponding internal ids for
711 # regions the slice spans
712 @coords =
713 $mapper->map( $slice_seq_region, $slice_start, $slice_end,
714 $slice_strand, $slice_cs );
715
716 @coords = grep { !$_->isa('Bio::EnsEMBL::Mapper::Gap') } @coords;
717
718 next COORD_SYSTEM if ( !@coords );
719
720 @ids = map { $_->id() } @coords;
721 #coords are now id rather than name
722 # @ids = @{$asma->seq_regions_to_ids($feat_cs, \@ids)};
723
724 # When regions are large and only partially spanned
725 # by slice it is faster to to limit the query with
726 # start and end constraints. Take simple approach:
727 # use regional constraints if there are less than a
728 # specific number of regions covered.
729
730 if ( @coords > $MAX_SPLIT_QUERY_SEQ_REGIONS ) {
731 my $constraint = $orig_constraint;
732 my $id_str = join( ',', @ids );
733 $constraint .= " AND " if ($constraint);
734 $constraint .= "${tab_syn}.seq_region_id IN ($id_str)";
735 my $fs = $self->generic_fetch( $constraint, $mapper, $slice );
736
737 $fs = $self->_remap( $fs, $mapper, $slice );
738
739 push @features, @$fs;
740
741 } else {
742 # do multiple split queries using start / end constraints
743
744 my $max_len = (
745 $self->_max_feature_length()
746 || $mcc->fetch_max_length_by_CoordSystem_feature_type(
747 $feat_cs, $tab_name
748 ) );
749
750 my $len = @coords;
751 for ( my $i = 0; $i < $len; $i++ ) {
752 my $constraint = $orig_constraint;
753 $constraint .= " AND " if ($constraint);
754 $constraint .=
755 "${tab_syn}.seq_region_id = "
756 . $ids[$i] . " AND "
757 . "${tab_syn}.seq_region_start <= "
758 . $coords[$i]->end() . " AND "
759 . "${tab_syn}.seq_region_end >= "
760 . $coords[$i]->start();
761
762 if ($max_len) {
763 my $min_start = $coords[$i]->start() - $max_len;
764 $constraint .=
765 " AND ${tab_syn}.seq_region_start >= $min_start";
766 }
767 my $fs = $self->generic_fetch( $constraint, $mapper, $slice );
768
769 $fs = $self->_remap( $fs, $mapper, $slice );
770
771 push @features, @$fs;
772 }
773 } ## end else [ if ( @coords > $MAX_SPLIT_QUERY_SEQ_REGIONS)]
774 } ## end else [ if ( $feat_cs->equals(...))]
775 } ## end foreach my $feat_cs (@feat_css)
776
777 return \@features;
778 } ## end sub _slice_fetch
779
780
781 #for a given seq_region_id, gets the one used in an external database, if present, otherwise, returns the internal one
782 sub get_seq_region_id_external {
783 my ( $self, $sr_id ) = @_;
784 my $cs_a = $self->db()->get_CoordSystemAdaptor();
785 return ( exists( $cs_a->{'_internal_seq_region_mapping'}->{$sr_id} )
786 ? $cs_a->{'_internal_seq_region_mapping'}->{$sr_id}
787 : $sr_id );
788 }
789
790 #for a given seq_region_id and coord_system, gets the one used in the internal (core) database
791 sub get_seq_region_id_internal{
792 my ( $self, $sr_id ) = @_;
793 my $cs_a = $self->db()->get_CoordSystemAdaptor();
794 return ( exists $cs_a->{'_external_seq_region_mapping'}->{$sr_id}
795 ? $cs_a->{'_external_seq_region_mapping'}->{$sr_id}
796 : $sr_id);
797 }
798
799 #
800 # Helper function containing some common feature storing functionality
801 #
802 # Given a Feature this will return a copy (or the same feature if no changes
803 # to the feature are needed) of the feature which is relative to the start
804 # of the seq_region it is on. The seq_region_id of the seq_region it is on
805 # is also returned.
806 #
807 # This method will also ensure that the database knows which coordinate
808 # systems that this feature is stored in.
809 #
810
811 sub _pre_store {
812 my $self = shift;
813 my $feature = shift;
814
815 if(!ref($feature) || !$feature->isa('Bio::EnsEMBL::Feature')) {
816 throw('Expected Feature argument.');
817 }
818 my $slice = $feature->slice();
819
820 $self->_check_start_end_strand($feature->start(),$feature->end(),
821 $feature->strand(), $slice);
822
823
824 my $db = $self->db();
825
826 my $slice_adaptor = $db->get_SliceAdaptor();
827
828 if(!ref($slice) || !($slice->isa('Bio::EnsEMBL::Slice') or $slice->isa('Bio::EnsEMBL::LRGSlice')) ) {
829 throw('Feature must be attached to Slice to be stored.');
830 }
831
832 # make sure feature coords are relative to start of entire seq_region
833
834 if($slice->start != 1 || $slice->strand != 1) {
835 #move feature onto a slice of the entire seq_region
836 $slice = $slice_adaptor->fetch_by_region($slice->coord_system->name(),
837 $slice->seq_region_name(),
838 undef, #start
839 undef, #end
840 undef, #strand
841 $slice->coord_system->version());
842
843 $feature = $feature->transfer($slice);
844
845 if(!$feature) {
846 throw('Could not transfer Feature to slice of ' .
847 'entire seq_region prior to storing');
848 }
849 }
850
851 # Ensure this type of feature is known to be stored in this coord system.
852 my $cs = $slice->coord_system;
853
854 my ($tab) = $self->_tables();
855 my $tabname = $tab->[0];
856
857 my $mcc = $db->get_MetaCoordContainer();
858
859 $mcc->add_feature_type($cs, $tabname, $feature->length);
860
861 my $seq_region_id = $slice_adaptor->get_seq_region_id($slice);
862
863 if(!$seq_region_id) {
864 throw('Feature is associated with seq_region which is not in this DB.');
865 }
866
867 return ($feature, $seq_region_id);
868 }
869
870
871 # The same function as _pre_store
872 # This one is used to store user uploaded features in XXX_userdata db
873
874 sub _pre_store_userdata {
875 my $self = shift;
876 my $feature = shift;
877
878 if(!ref($feature) || !$feature->isa('Bio::EnsEMBL::Feature')) {
879 throw('Expected Feature argument.');
880 }
881
882 my $slice = $feature->slice();
883 my $slice_adaptor = $slice->adaptor;
884
885 $self->_check_start_end_strand($feature->start(),$feature->end(),
886 $feature->strand(), $slice);
887
888
889 if(!ref($slice) || !($slice->isa('Bio::EnsEMBL::Slice') or $slice->isa('Bio::EnsEMBL::LRGSlice')) ) {
890 throw('Feature must be attached to Slice to be stored.');
891 }
892
893 # make sure feature coords are relative to start of entire seq_region
894
895 if($slice->start != 1 || $slice->strand != 1) {
896 #move feature onto a slice of the entire seq_region
897 $slice = $slice_adaptor->fetch_by_region($slice->coord_system->name(),
898 $slice->seq_region_name(),
899 undef, #start
900 undef, #end
901 undef, #strand
902 $slice->coord_system->version());
903
904 $feature = $feature->transfer($slice);
905
906 if(!$feature) {
907 throw('Could not transfer Feature to slice of ' .
908 'entire seq_region prior to storing');
909 }
910 }
911
912 # Ensure this type of feature is known to be stored in this coord system.
913 my $cs = $slice->coord_system;
914
915 my ($tab) = $self->_tables();
916 my $tabname = $tab->[0];
917
918 my $db = $self->db;
919 my $mcc = $db->get_MetaCoordContainer();
920
921 $mcc->add_feature_type($cs, $tabname, $feature->length);
922
923 my $seq_region_id = $slice_adaptor->get_seq_region_id($slice);
924
925 if(!$seq_region_id) {
926 throw('Feature is associated with seq_region which is not in this DB.');
927 }
928
929 return ($feature, $seq_region_id);
930 }
931
932
933 #
934 # helper function used to validate start/end/strand and
935 # hstart/hend/hstrand etc.
936 #
937 sub _check_start_end_strand {
938 my $self = shift;
939 my $start = shift;
940 my $end = shift;
941 my $strand = shift;
942 my $slice = shift;
943
944 #
945 # Make sure that the start, end, strand are valid
946 #
947 if(int($start) != $start) {
948 throw("Invalid Feature start [$start]. Must be integer.");
949 }
950 if(int($end) != $end) {
951 throw("Invalid Feature end [$end]. Must be integer.");
952 }
953 if(int($strand) != $strand || $strand < -1 || $strand > 1) {
954 throw("Invalid Feature strand [$strand]. Must be -1, 0 or 1.");
955 }
956 if($end < $start && !$slice->is_circular()) {
957 throw("Invalid Feature start/end [$start/$end]. Start must be less " .
958 "than or equal to end.");
959 }
960
961 return 1;
962 }
963
964
965 #
966 # Given a list of features checks if they are in the correct coord system
967 # by looking at the first features slice. If they are not then they are
968 # converted and placed on the slice.
969 #
970 sub _remap {
971 my ( $self, $features, $mapper, $slice ) = @_;
972
973 #check if any remapping is actually needed
974 if(@$features && (!$features->[0]->isa('Bio::EnsEMBL::Feature') ||
975 $features->[0]->slice == $slice)) {
976 return $features;
977 }
978
979 #remapping has not been done, we have to do our own conversion from
980 #to slice coords
981
982 my @out;
983
984 my $slice_start = $slice->start();
985 my $slice_end = $slice->end();
986 my $slice_strand = $slice->strand();
987 my $slice_cs = $slice->coord_system();
988
989 my ($seq_region, $start, $end, $strand);
990
991 my $slice_seq_region_id = $slice->get_seq_region_id();
992 my $slice_seq_region = $slice->seq_region_name();
993
994 foreach my $f (@$features) {
995 #since feats were obtained in contig coords, attached seq is a contig
996 my $fslice = $f->slice();
997 if(!$fslice) {
998 throw("Feature does not have attached slice.\n");
999 }
1000 my $fseq_region = $fslice->seq_region_name();
1001 my $fseq_region_id = $fslice->get_seq_region_id();
1002 my $fcs = $fslice->coord_system();
1003
1004 if(!$slice_cs->equals($fcs)) {
1005 #slice of feature in different coord system, mapping required
1006
1007 ($seq_region, $start, $end, $strand) =
1008 $mapper->fastmap($fseq_region_id,$f->start(),$f->end(),$f->strand(),$fcs);
1009
1010 # undefined start means gap
1011 next if(!defined $start);
1012 } else {
1013 $start = $f->start();
1014 $end = $f->end();
1015 $strand = $f->strand();
1016 $seq_region = $f->slice->seq_region_name();
1017 }
1018
1019 # maps to region outside desired area
1020 next if ($start > $slice_end) || ($end < $slice_start) ||
1021 ($slice_seq_region ne $seq_region);
1022
1023 #shift the feature start, end and strand in one call
1024 if($slice_strand == -1) {
1025 $f->move( $slice_end - $end + 1, $slice_end - $start + 1, $strand * -1 );
1026 } else {
1027 $f->move( $start - $slice_start + 1, $end - $slice_start + 1, $strand );
1028 }
1029
1030 $f->slice($slice);
1031
1032 push @out,$f;
1033 }
1034
1035 return \@out;
1036 }
1037
1038
1039 #
1040 # Given a logic name and an existing constraint this will
1041 # add an analysis table constraint to the feature. Note that if no
1042 # analysis_id exists in the columns of the primary table then no
1043 # constraint is added at all
1044 #
1045 sub _logic_name_to_constraint {
1046 my $self = shift;
1047 my $constraint = shift;
1048 my $logic_name = shift;
1049
1050 return $constraint if(!$logic_name);
1051
1052 #make sure that an analysis_id exists in the primary table
1053 my ($prim_tab) = $self->_tables();
1054 my $prim_synonym = $prim_tab->[1];
1055
1056 my $found_analysis=0;
1057 foreach my $col ($self->_columns) {
1058 my ($syn,$col_name) = split(/\./,$col);
1059 next if($syn ne $prim_synonym);
1060 if($col_name eq 'analysis_id') {
1061 $found_analysis = 1;
1062 last;
1063 }
1064 }
1065
1066 if(!$found_analysis) {
1067 warning("This feature is not associated with an analysis.\n" .
1068 "Ignoring logic_name argument = [$logic_name].\n");
1069 return $constraint;
1070 }
1071
1072 my $aa = $self->db->get_AnalysisAdaptor();
1073 my $an = $aa->fetch_by_logic_name($logic_name);
1074
1075 if ( !defined($an) ) {
1076 return undef;
1077 }
1078
1079 my $an_id = $an->dbID();
1080
1081 $constraint .= ' AND' if($constraint);
1082 $constraint .= " ${prim_synonym}.analysis_id = $an_id";
1083 return $constraint;
1084 }
1085
1086
1087 =head2 store
1088
1089 Arg [1] : list of Bio::EnsEMBL::SeqFeature
1090 Example : $adaptor->store(@feats);
1091 Description: ABSTRACT Subclasses are responsible for implementing this
1092 method. It should take a list of features and store them in
1093 the database.
1094 Returntype : none
1095 Exceptions : thrown method is not implemented by subclass
1096 Caller : general
1097 Status : At Risk
1098 : throws if called.
1099
1100 =cut
1101
1102 sub store{
1103 my $self = @_;
1104
1105 throw("Abstract method store not defined by implementing subclass\n");
1106 }
1107
1108
1109 =head2 remove
1110
1111 Arg [1] : A feature $feature
1112 Example : $feature_adaptor->remove($feature);
1113 Description: This removes a feature from the database. The table the
1114 feature is removed from is defined by the abstract method
1115 _tablename, and the primary key of the table is assumed
1116 to be _tablename() . '_id'. The feature argument must
1117 be an object implementing the dbID method, and for the
1118 feature to be removed from the database a dbID value must
1119 be returned.
1120 Returntype : none
1121 Exceptions : thrown if $feature arg does not implement dbID(), or if
1122 $feature->dbID is not a true value
1123 Caller : general
1124 Status : Stable
1125
1126 =cut
1127
1128
1129 sub remove {
1130 my ($self, $feature) = @_;
1131
1132 if(!$feature || !ref($feature) || !$feature->isa('Bio::EnsEMBL::Feature')) {
1133 throw('Feature argument is required');
1134 }
1135
1136 if(!$feature->is_stored($self->db)) {
1137 throw("This feature is not stored in this database");
1138 }
1139
1140 my @tabs = $self->_tables;
1141 my ($table) = @{$tabs[0]};
1142
1143 my $sth = $self->prepare("DELETE FROM $table WHERE ${table}_id = ?");
1144 $sth->bind_param(1,$feature->dbID,SQL_INTEGER);
1145 $sth->execute();
1146
1147 #unset the feature dbID ad adaptor
1148 $feature->dbID(undef);
1149 $feature->adaptor(undef);
1150
1151 return;
1152 }
1153
1154
1155 =head2 remove_by_Slice
1156
1157 Arg [1] : Bio::Ensembl::Slice $slice
1158 Example : $feature_adaptor->remove_by_Slice($slice);
1159 Description: This removes features from the database which lie on a region
1160 represented by the passed in slice. Only features which are
1161 fully contained by the slice are deleted; features which overlap
1162 the edge of the slice are not removed.
1163 The table the features are removed from is defined by
1164 the abstract method_tablename.
1165 Returntype : none
1166 Exceptions : thrown if no slice is supplied
1167 Caller : general
1168 Status : Stable
1169
1170 =cut
1171
1172 sub remove_by_Slice {
1173 my ($self, $slice) = @_;
1174
1175 if(!$slice || !ref($slice) || !($slice->isa('Bio::EnsEMBL::Slice') or $slice->isa('Bio::EnsEMBL::LRGSlice')) ) {
1176 throw("Slice argument is required");
1177 }
1178
1179 my @tabs = $self->_tables;
1180 my ($table_name) = @{$tabs[0]};
1181
1182 my $seq_region_id = $self->db->get_SliceAdaptor->get_seq_region_id($slice);
1183 my $start = $slice->start();
1184 my $end = $slice->end();
1185
1186 #
1187 # Delete only features fully on the slice, not overlapping ones
1188 #
1189 my $sth = $self->prepare("DELETE FROM $table_name " .
1190 "WHERE seq_region_id = ? " .
1191 "AND seq_region_start >= ? " .
1192 "AND seq_region_end <= ?");
1193
1194 $sth->bind_param(1,$seq_region_id,SQL_INTEGER);
1195 $sth->bind_param(2,$start,SQL_INTEGER);
1196 $sth->bind_param(3,$end,SQL_INTEGER);
1197 $sth->execute();
1198 $sth->finish();
1199 }
1200
1201
1202 #
1203 # Internal function. Allows the max feature length which is normally
1204 # retrieved from the meta_coord table to be overridden. This allows
1205 # for some significant optimizations to be put in when it is known
1206 # that requested features will not be over a certain size.
1207 #
1208 sub _max_feature_length {
1209 my $self = shift;
1210 return $self->{'_max_feature_length'} = shift if(@_);
1211 return $self->{'_max_feature_length'};
1212 }
1213
1214
1215 #
1216 # Lists all seq_region_ids that a particular feature type is found on.
1217 # Useful e.g. for finding out which seq_regions have genes.
1218 # Returns a listref of seq_region_ids.
1219 #
1220 sub _list_seq_region_ids {
1221 my ($self, $table) = @_;
1222
1223 my @out;
1224
1225 my $sql = qq(
1226 SELECT DISTINCT
1227 sr.seq_region_id
1228 FROM seq_region sr,
1229 $table a,
1230 coord_system cs
1231 WHERE sr.seq_region_id = a.seq_region_id
1232 AND sr.coord_system_id = cs.coord_system_id
1233 AND cs.species_id = ?);
1234
1235 my $sth = $self->prepare($sql);
1236
1237 $sth->bind_param( 1, $self->species_id(), SQL_INTEGER );
1238
1239 $sth->execute();
1240
1241 while (my ($id) = $sth->fetchrow) {
1242 push(@out, $id);
1243 }
1244
1245 $sth->finish;
1246
1247 return \@out;
1248 }
1249
1250
1251 =head1 DEPRECATED METHODS
1252
1253 =cut
1254
1255
1256 =head2 fetch_all_by_RawContig_constraint
1257
1258 Description: DEPRECATED use fetch_all_by_RawContig_constraint instead
1259
1260 =cut
1261
1262 sub fetch_all_by_RawContig_constraint {
1263 my $self = shift;
1264 deprecate('Use fetch_all_by_Slice_constraint() instead.');
1265 return $self->fetch_all_by_slice_constraint(@_);
1266 }
1267
1268 =head2 fetch_all_by_RawContig
1269
1270 Description: DEPRECATED use fetch_all_by_Slice instead
1271
1272 =cut
1273
1274 sub fetch_all_by_RawContig {
1275 my $self = shift;
1276 deprecate('Use fetch_all_by_Slice() instead.');
1277 return $self->fetch_all_by_Slice(@_);
1278 }
1279
1280 =head2 fetch_all_by_RawContig_and_score
1281
1282 Description: DEPRECATED use fetch_all_by_Slice_and_score instead
1283
1284 =cut
1285
1286 sub fetch_all_by_RawContig_and_score{
1287 my $self = shift;
1288 deprecate('Use fetch_all_by_Slice_and_score() instead.');
1289 return $self->fetch_all_by_Slice_and_score(@_);
1290 }
1291
1292 =head2 remove_by_RawContig
1293
1294 Description: DEPRECATED use remove_by_Slice instead
1295
1296 =cut
1297
1298 sub remove_by_RawContig {
1299 my $self = shift;
1300 deprecate("Use remove_by_Slice instead");
1301 return $self->remove_by_Slice(@_);
1302 }
1303
1304
1305 sub remove_by_analysis_id {
1306 my ($self, $analysis_id) = @_;
1307
1308 $analysis_id or throw("Must call with analysis id");
1309
1310 my @tabs = $self->_tables;
1311 my ($tablename) = @{$tabs[0]};
1312
1313 my $sql = "DELETE FROM $tablename WHERE analysis_id = $analysis_id";
1314 # warn "SQL : $sql";
1315
1316 my $sth = $self->prepare($sql);
1317 $sth->execute();
1318 $sth->finish();
1319 }
1320
1321 sub remove_by_feature_id {
1322 my ($self, $features_list) = @_;
1323
1324 my @feats = @$features_list or throw("Must call store with features");
1325
1326 my @tabs = $self->_tables;
1327 my ($tablename) = @{$tabs[0]};
1328
1329 my $sql = sprintf "DELETE FROM $tablename WHERE ${tablename}_id IN (%s)", join ', ', @feats;
1330 # warn "SQL : $sql";
1331
1332 my $sth = $self->prepare($sql);
1333 $sth->execute();
1334 $sth->finish();
1335 }
1336
1337
1338 1;