Mercurial > repos > mahtabm > ensembl
comparison variant_effect_predictor/Bio/EnsEMBL/IndividualSlice.pm @ 0:1f6dce3d34e0
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::IndividualSlice - SubClass of the Slice. Represents the | |
24 slice of the genome for a certain individual (applying the alleles for | |
25 this individual) | |
26 | |
27 =head1 SYNOPSIS | |
28 | |
29 $sa = $db->get_SliceAdaptor; | |
30 | |
31 $slice = | |
32 $sa->fetch_by_region( 'chromosome', 'X', 1_000_000, 2_000_000 ); | |
33 | |
34 $individualSlice = $slice->get_by_Individual($individual_name); | |
35 | |
36 # Get the sequence from the Individual Slice: will contain IUPAC codes | |
37 # for SNPs and Ensembl ambiguity codes for indels | |
38 my $seq = $individualSlice->seq(); | |
39 print $seq; | |
40 | |
41 # Get a subSlice of the Strain | |
42 my $subSlice_individual = | |
43 $individualSlice->sub_Slice( 5_000, 8_000, 1 ) | |
44 | |
45 # Compare two different individuals in the same Slice | |
46 my $sliceIndividual2 = $slice->get_by_Individual($individual_name2); | |
47 my $differences = | |
48 $individualSlice->get_all_differences_IndividualSlice( | |
49 $sliceIndividual2); | |
50 | |
51 foreach my $af ( @{$differences} ) { | |
52 print | |
53 "There is a difference between $individual_name " | |
54 . "and $individual_name2 at ", | |
55 $af->start, "-", $af->end, | |
56 " with allele ", $af->allele_string(), "\n"; | |
57 } | |
58 | |
59 =head1 DESCRIPTION | |
60 | |
61 A IndividualSlice object represents a region of a genome for a certain | |
62 individual. It can be used to retrieve sequence or features from a | |
63 individual. | |
64 | |
65 =head1 METHODS | |
66 | |
67 =cut | |
68 | |
69 package Bio::EnsEMBL::IndividualSlice; | |
70 use vars qw(@ISA); | |
71 use strict; | |
72 | |
73 use Bio::EnsEMBL::Utils::Argument qw(rearrange); | |
74 use Bio::EnsEMBL::Utils::Sequence qw(reverse_comp); | |
75 use Bio::EnsEMBL::Slice; | |
76 use Bio::EnsEMBL::Mapper; | |
77 use Bio::EnsEMBL::Utils::Exception qw(throw deprecate warning); | |
78 | |
79 @ISA = qw(Bio::EnsEMBL::Slice); | |
80 | |
81 | |
82 =head2 new | |
83 | |
84 Arg [1..N] : List of named arguments | |
85 Bio::EnsEMBL::CoordSystem COORD_SYSTEM | |
86 string SEQ_REGION_NAME, | |
87 int START, | |
88 int END, | |
89 string VERSION (optional, defaults to '') | |
90 int STRAND, (optional, defaults to 1) | |
91 Bio::EnsEMBL::DBSQL::SliceAdaptor ADAPTOR (optional) | |
92 Arg[N+1] : string $individual_name | |
93 Example : $individualSlice = Bio::EnsEMBL::IndividualSlice->new(-coord_system => $cs, | |
94 -start => 1, | |
95 -end => 10000, | |
96 -strand => 1, | |
97 -seq_region_name => 'X', | |
98 -seq_region_length => 12e6, | |
99 -individual_name => $individual_name); | |
100 Description : Creates a new Bio::EnsEMBL::IndividualSlice object that will contain a shallow copy of the | |
101 Slice object, plus additional information such as the individual this Slice refers to | |
102 and listref of Bio::EnsEMBL::Variation::AlleleFeatures of differences with the | |
103 reference sequence | |
104 ReturnType : Bio::EnsEMBL::IndividualSlice | |
105 Exceptions : none | |
106 Caller : general | |
107 | |
108 =cut | |
109 | |
110 sub new{ | |
111 my $caller = shift; | |
112 my $class = ref($caller) || $caller; | |
113 | |
114 #create the IndividualSlice object as the Slice, plus the individual attribute | |
115 my ($individual_name, $sample_id) = rearrange(['INDIVIDUAL', 'SAMPLE_ID'],@_); | |
116 | |
117 my $self = $class->SUPER::new(@_); | |
118 | |
119 $self->{'individual_name'} = $individual_name; | |
120 $self->{'sample_id'} = $sample_id; | |
121 | |
122 return $self; | |
123 | |
124 } | |
125 | |
126 =head2 individual_name | |
127 | |
128 Arg [1] : (optional) string $individual_name | |
129 Example : my $individual_name = $individualSlice->individual_name(); | |
130 Description : Getter/Setter for the name of the individual in the slice | |
131 ReturnType : string | |
132 Exceptions : none | |
133 Caller : general | |
134 | |
135 =cut | |
136 | |
137 sub individual_name{ | |
138 my $self = shift; | |
139 if (@_){ | |
140 $self->{'individual_name'} = shift @_; | |
141 } | |
142 return $self->{'individual_name'}; | |
143 } | |
144 | |
145 =head2 seq | |
146 | |
147 Arg [1] : none | |
148 Example : print "SEQUENCE = ", $strainSlice->seq(); | |
149 Description: Returns the sequence of the region represented by this | |
150 StrainSlice formatted as a string. | |
151 Returntype : string | |
152 Exceptions : none | |
153 Caller : general | |
154 | |
155 =cut | |
156 | |
157 sub seq { | |
158 my $self = shift; | |
159 | |
160 # special case for in-between (insert) coordinates | |
161 return '' if($self->start() == $self->end() + 1); | |
162 | |
163 return $self->{'seq'} if($self->{'seq'}); | |
164 | |
165 if($self->adaptor()) { | |
166 my $seqAdaptor = $self->adaptor()->db()->get_SequenceAdaptor(); | |
167 my $reference_sequence = $seqAdaptor->fetch_by_Slice_start_end_strand($self,1,undef,1); #get the reference sequence for that slice | |
168 #apply all differences to the reference sequence | |
169 | |
170 # sort edits in reverse order to remove complication of | |
171 # adjusting downstream edits | |
172 my @allele_features_ordered = sort {$b->start() <=> $a->start() || $b->end() <=> $a->end()} @{$self->{'alleleFeatures'}} if (defined $self->{'alleleFeatures'}); | |
173 | |
174 foreach my $af (@allele_features_ordered){ | |
175 $af->apply_edit($reference_sequence); #change, in the reference sequence, the af | |
176 } | |
177 # return substr(${$reference_sequence},0,1) if ($self->length == 1); | |
178 return ${$reference_sequence}; #returns the reference sequence, applying the alleleFeatures | |
179 } | |
180 | |
181 # no attached sequence, and no db, so just return Ns | |
182 return 'N' x $self->length(); | |
183 } | |
184 | |
185 =head2 get_all_differences_Slice | |
186 | |
187 Args : none | |
188 Example : my $differences = $individualSlice->get_all_differences_Slice() | |
189 Description : Gets all differences between the IndividualSlice object and the Slice is defined | |
190 ReturnType : listref of Bio::EnsEMBL::Variation::AlleleFeature | |
191 Exceptions : none | |
192 Caller : general | |
193 | |
194 =cut | |
195 | |
196 sub get_all_differences_Slice{ | |
197 my $self = shift; | |
198 my $differences; #reference to the array with the differences between Slice and StrainSlice | |
199 my $ref_allele; | |
200 foreach my $difference (@{$self->{'alleleFeatures'}}){ | |
201 if ($difference->length_diff == 0){ | |
202 #the difference is a SNP, check if it is the same as the reference allele | |
203 $ref_allele = $self->SUPER::subseq($difference->start,$difference->end,$difference->strand); | |
204 $ref_allele = '-' if ($ref_allele eq ''); | |
205 if ($ref_allele ne $difference->allele_string){ | |
206 #when the alleleFeature is different from the reference allele, add to the differences list | |
207 push @{$differences},$difference; | |
208 } | |
209 } | |
210 else{ | |
211 push @{$differences},$difference; | |
212 } | |
213 } | |
214 | |
215 return $differences; | |
216 | |
217 } | |
218 | |
219 =head2 get_all_differences_IndividualSlice | |
220 | |
221 Arg[1] : Bio::EnsEMBL::IndividualSlice $is | |
222 Example : my $differences = $individualSlice->get_all_differences_IndividualSlice($individualslice) | |
223 Description : Gets differences between 2 IndividualSlice objects | |
224 ReturnType : listref of Bio::EnsEMBL::Variation::AlleleFeature | |
225 Exceptions : thrown on bad argument | |
226 Caller : general | |
227 | |
228 =cut | |
229 | |
230 sub get_all_differences_IndividualSlice{ | |
231 my $self = shift; | |
232 my $individualSlice = shift; | |
233 | |
234 if (!ref($individualSlice) || !$individualSlice->isa('Bio::EnsEMBL::IndividualSlice')){ | |
235 throw('Bio::EnsEMBL::IndividualSlice arg expected'); | |
236 } | |
237 if ( @{$self->{'alleleFeatures'}} == 0 && @{$individualSlice->{'alleleFeatures'}} == 0){ | |
238 return undef; #there are no differences in any of the Individuals | |
239 | |
240 } | |
241 my $differences; #differences between individuals | |
242 if (@{$individualSlice->{'alleleFeatures'}} == 0){ | |
243 #need to create a copy of alleleFeature for the first Individual | |
244 foreach my $difference (@{$self->{'alleleFeatures'}}){ | |
245 my %vf = %$difference; | |
246 push @{$differences},bless \%vf,ref($difference); | |
247 } | |
248 } | |
249 elsif (@{$self->{'alleleFeatures'}} == 0){ | |
250 #need to create a copy of AlleleFeature, but changing the allele by the allele in the reference sequence | |
251 foreach my $difference (@{$individualSlice->{'alleleFeatures'}}){ | |
252 push @{$differences}, $individualSlice->_convert_difference($difference); | |
253 } | |
254 } | |
255 else{ | |
256 #both individuals have differences | |
257 #create a hash with the differences in the first slice | |
258 my %allele_features_self = map {$_->start.'-'.$_->end => $_} @{$self->{'alleleFeatures'}}; | |
259 foreach my $difference (@{$individualSlice->{'alleleFeatures'}}){ | |
260 #there is no difference in the other individual slice, convert the allele | |
261 if (!defined $allele_features_self{$difference->start.'-'.$difference->end}){ | |
262 push @{$differences},$individualSlice->_convert_difference($difference); | |
263 } | |
264 else{ | |
265 #if it is defined and have the same allele, delete from the hash since it is not a difference | |
266 #between the individuals | |
267 if ($allele_features_self{$difference->start.'-'.$difference->end}->allele_string eq $difference->allele_string){ | |
268 delete $allele_features_self{$difference->start.'-'.$difference->end}; | |
269 } | |
270 } | |
271 } | |
272 #and finally, make a shallow copy of the differences in the first individual | |
273 foreach my $difference (values %allele_features_self){ | |
274 my %vf = %$difference; | |
275 push @{$differences},bless \%vf,ref($difference); | |
276 } | |
277 | |
278 } | |
279 #need to map differences to the first individual, self, since the coordinates are in the Slice coordinate system | |
280 my $mapper = $self->mapper(); #now that we have the differences, map them in the IndividualSlice | |
281 my @results; | |
282 foreach my $difference (@{$differences}){ | |
283 @results = $mapper->map_coordinates('Slice',$difference->start,$difference->end,$difference->strand,'Slice'); | |
284 #we can have 3 possibilities: | |
285 #the difference is an insertion and when mapping returns the boundaries of the insertion in the IndividualSlice | |
286 if (@results == 2){ | |
287 #the first position in the result is the beginning of the insertion | |
288 if($results[0]->start < $results[1]->start){ | |
289 $difference->start($results[0]->end+1); | |
290 $difference->end($results[1]->start-1); | |
291 } | |
292 else{ | |
293 #it is the second position the beginning of the insertion | |
294 $difference->start($results[1]->end+1); | |
295 $difference->end($results[0]->start-1); | |
296 } | |
297 $difference->strand($results[0]->strand); | |
298 } | |
299 else{ | |
300 #it can be either a SNP or a deletion, and we have the coordinates in the result, etither a Bio::EnsEMBL::Mapper::Coordinate | |
301 # or a Bio::EnsEMBL::Mapper::IndelCoordinate | |
302 $difference->start($results[0]->start); | |
303 $difference->end($results[0]->end); | |
304 $difference->strand($results[0]->strand); | |
305 } | |
306 } | |
307 | |
308 return $differences; | |
309 } | |
310 | |
311 #for a given AlleleFeature, converts the allele into the reference allele and returns | |
312 #the converted AlleleFeature | |
313 | |
314 sub _convert_difference{ | |
315 my $self = shift; | |
316 my $difference = shift; | |
317 my %new_af = %$difference; #make a copy of the alleleFeature | |
318 #and change the allele with the one from the reference Slice | |
319 $new_af{'allele_string'} = $self->SUPER::subseq($difference->start,$difference->end,$difference->strand); | |
320 return bless \%new_af,ref($difference); | |
321 } | |
322 | |
323 =head2 mapper | |
324 | |
325 Args : none | |
326 Description: Getter for the mapper between the between the IndividualSlice and the Slice it refers to. | |
327 It is done automatically when necessary to create subSlice or to get the differences between individuals | |
328 Returntype : Bio::EnsEMBL::Mapper | |
329 Exceptions : none | |
330 Caller : Internal function | |
331 | |
332 =cut | |
333 | |
334 sub mapper{ | |
335 my $self = shift; | |
336 | |
337 if (@_) { | |
338 #allow to create again the mapper | |
339 delete $self->{'mapper'}; | |
340 } | |
341 if(!defined $self->{'mapper'}){ | |
342 #create the mapper between the Slice and StrainSlice | |
343 my $mapper = Bio::EnsEMBL::Mapper->new('Slice','IndividualSlice'); | |
344 #align with Slice | |
345 #get all the VariationFeatures in the Individual Slice, from start to end in the Slice | |
346 my @allele_features_ordered = sort {$a->start() <=> $b->start() || $b->end() <=> $a->end()} @{$self->{'alleleFeatures'}} if (defined $self->{'alleleFeatures'}); | |
347 | |
348 my $start_slice = 1; | |
349 my $end_slice; | |
350 my $start_individual = 1; | |
351 my $end_individual; | |
352 my $length_allele; | |
353 my $total_length_diff = 0; | |
354 #we will walk from left to right in the slice object, updating the start and end individual every time | |
355 #there is a new alleleFeature in the Individual | |
356 foreach my $allele_feature (@allele_features_ordered){ | |
357 #we have a insertion/deletion: marks the beginning of new slice move coordinates | |
358 if ($allele_feature->length_diff != 0){ | |
359 $total_length_diff += $allele_feature->length_diff; | |
360 $length_allele = $allele_feature->length + $allele_feature->length_diff(); #length of the allele in the Individual | |
361 $end_slice = $allele_feature->start() - 1; #set the end of the slice before the alleleFeature | |
362 if ($end_slice >= $start_slice){ | |
363 #normal cases (not with gaps) | |
364 $end_individual = $end_slice - $start_slice + $start_individual; #set the end of the individual from the beginning plus the offset | |
365 #add the sequence that maps | |
366 $mapper->add_map_coordinates('Slice',$start_slice,$end_slice,1,'IndividualSlice',$start_individual,$end_individual); | |
367 #and add the indel | |
368 $mapper->add_indel_coordinates('Slice',$end_slice+1,$end_slice + $allele_feature->length,1,'IndividualSlice',$end_individual+1,$end_individual + $length_allele); | |
369 $start_individual = $end_individual + $length_allele + 1; #set the beginning of the individual after the allele | |
370 } | |
371 else{ | |
372 #add the indel | |
373 $mapper->add_indel_coordinates('Slice',$end_slice+1,$end_slice + $allele_feature->length,1,'IndividualSlice',$end_individual+1,$end_individual + $length_allele); | |
374 $start_individual += $length_allele; | |
375 } | |
376 $start_slice = $end_slice + $allele_feature->length+ 1; #set the beginning of the slice after the variation feature | |
377 } | |
378 } | |
379 if ($start_slice <= $self->length){ | |
380 #if we haven't reached the end of the IndividualSlice, add the final map coordinates between the individual and the slice | |
381 $mapper->add_map_coordinates('Slice',$start_slice,$self->length,1,'IndividualSlice',$start_individual,$start_individual + $self->length - $start_slice); | |
382 } | |
383 | |
384 $mapper->add_map_coordinates('Slice', -$self->start+1, 0,1, 'IndividualSlice', -$self->start +1,0) if ($self->start > 0); #before individualSlice | |
385 $mapper->add_map_coordinates('Slice', $self->length + 1,$self->seq_region_length - ($self->length +1),1, 'IndividualSlice', $self->length + 1 + $total_length_diff,$self->seq_region_length + $total_length_diff - ($self->length +1) ) if ($self->length <= $self->seq_region_length); #after strainSlice | |
386 $self->{'mapper'} = $mapper; | |
387 } | |
388 return $self->{'mapper'}; | |
389 } | |
390 | |
391 =head2 sub_Slice | |
392 | |
393 Arg 1 : int $start | |
394 Arg 2 : int $end | |
395 Arge [3] : int $strand | |
396 Example : none | |
397 Description: Makes another IndividualSlice that covers only part of this IndividualSlice | |
398 with the appropriate differences to the reference Slice | |
399 If a slice is requested which lies outside of the boundaries | |
400 of this function will return undef. This means that | |
401 behaviour will be consistant whether or not the slice is | |
402 attached to the database (i.e. if there is attached sequence | |
403 to the slice). Alternatively the expand() method or the | |
404 SliceAdaptor::fetch_by_region method can be used instead. | |
405 Returntype : Bio::EnsEMBL::IndividualSlice or undef if arguments are wrong | |
406 Exceptions : thrown when trying to get the subSlice in the middle of a | |
407 insertion | |
408 Caller : general | |
409 | |
410 =cut | |
411 | |
412 sub sub_Slice { | |
413 my ( $self, $start, $end, $strand ) = @_; | |
414 my $mapper = $self->mapper(); | |
415 #map from the Individual to the Slice to get the sub_Slice, and then, apply the differences in the subSlice | |
416 my @results = $mapper->map_coordinates('IndividualSlice',$start,$end,$strand,'IndividualSlice'); | |
417 my $new_start; | |
418 my $new_end; | |
419 my $new_strand; | |
420 my $new_seq; | |
421 #Get need start and end for the subSlice of the IndividualSlice | |
422 my @results_ordered = sort {$a->start <=> $b->start} grep {ref($_) eq 'Bio::EnsEMBL::Mapper::Coordinate'} @results; | |
423 $new_start = $results_ordered[0]->start(); | |
424 $new_strand = $results_ordered[0]->strand() if (ref($results_ordered[0]) eq 'Bio::EnsEMBL::Mapper::Coordinate'); | |
425 # $new_strand = $results_ordered[-1]->strand() if (ref($results_ordered[-1]) eq 'Bio::EnsEMBL::Mapper::Coordinate'); | |
426 $new_end = $results_ordered[-1]->end(); #get last element of the array, the end of the slice | |
427 | |
428 my $subSlice = $self->SUPER::sub_Slice($new_start,$new_end,$new_strand); | |
429 $subSlice->{'individual_name'} = $self->{'individual_name'}; | |
430 | |
431 my $new_alleles; #reference to an array that will contain the variationFeatures in the new subSlice | |
432 #update the VariationFeatures in the sub_Slice of the Individual | |
433 my %af; | |
434 my $new_allele_feature; | |
435 foreach my $alleleFeature (@{$self->{'alleleFeatures'}}){ | |
436 $new_allele_feature = $alleleFeature->transfer($subSlice); | |
437 #only transfer the coordinates to the SubSlice that are within the boundaries | |
438 if ($new_allele_feature->start >= 1 && $new_allele_feature->end <= $subSlice->length){ | |
439 push @{$new_alleles}, $new_allele_feature; | |
440 } | |
441 } | |
442 $subSlice->{'alleleFeatures'} = $new_alleles; | |
443 return $subSlice; | |
444 | |
445 } | |
446 | |
447 =head2 subseq | |
448 | |
449 Arg [1] : int $startBasePair | |
450 relative to start of slice, which is 1. | |
451 Arg [2] : int $endBasePair | |
452 relative to start of slice. | |
453 Arg [3] : (optional) int $strand | |
454 The strand of the individual slice to obtain sequence from. Default | |
455 value is 1. | |
456 Description: returns string of dna sequence | |
457 Returntype : txt | |
458 Exceptions : end should be at least as big as start | |
459 strand must be set | |
460 Caller : general | |
461 | |
462 =cut | |
463 | |
464 sub subseq { | |
465 my ( $self, $start, $end, $strand ) = @_; | |
466 | |
467 if ( $end+1 < $start ) { | |
468 throw("End coord + 1 is less than start coord"); | |
469 } | |
470 | |
471 # handle 'between' case for insertions | |
472 return '' if( $start == $end + 1); | |
473 | |
474 $strand = 1 unless(defined $strand); | |
475 | |
476 if ( $strand != -1 && $strand != 1 ) { | |
477 throw("Invalid strand [$strand] in call to Slice::subseq."); | |
478 } | |
479 | |
480 my $subseq; | |
481 my $seq; | |
482 if($self->adaptor){ | |
483 my $seqAdaptor = $self->adaptor()->db()->get_SequenceAdaptor(); | |
484 $subseq = ${$seqAdaptor->fetch_by_Slice_start_end_strand($self,$start,$end,$strand)}; #get the reference sequence for that slice | |
485 #apply all differences to the reference sequence | |
486 # sort edits in reverse order to remove complication of | |
487 # adjusting downstream edits | |
488 my @allele_features_ordered = sort {$b->start() <=> $a->start() || $b->end() <=> $a->end()} @{$self->{'alleleFeatures'}} if (defined $self->{'alleleFeatures'}); | |
489 my $af_start; | |
490 my $af_end; | |
491 foreach my $af (@allele_features_ordered){ | |
492 if (($af->start - $start +1 > 0) && ($end - $af->end > 0)){ | |
493 #save the current start and end of the alleleFeature before changing for apply_edit | |
494 $af_start = $af->start; | |
495 $af_end = $af->end; | |
496 #apply the difference if the feature is in the new slice | |
497 $af->start($af->start - $start +1); | |
498 $af->end($af->end - $start +1); | |
499 $af->apply_edit(\$subseq); #change, in the reference sequence, the af | |
500 #restore the initial values of alleleFeature start and end | |
501 $af->start($af_start); | |
502 $af->end($af_end); | |
503 | |
504 } | |
505 } | |
506 } | |
507 else { | |
508 ## check for gap at the beginning and pad it with Ns | |
509 if ($start < 1) { | |
510 $subseq = "N" x (1 - $start); | |
511 $start = 1; | |
512 } | |
513 $subseq .= substr ($self->seq(), $start-1, $end - $start + 1); | |
514 ## check for gap at the end and pad it with Ns | |
515 if ($end > $self->length()) { | |
516 $subseq .= "N" x ($end - $self->length()); | |
517 } | |
518 reverse_comp(\$subseq) if($strand == -1); | |
519 } | |
520 return $subseq; | |
521 | |
522 } | |
523 | |
524 =head2 get_all_Transcripts | |
525 | |
526 Args : None | |
527 Example : @transcripts = @{$individualslice->get_all_Transcripts)}; | |
528 Description: Gets all transcripts which overlap this Individual Slice. If you want to | |
529 specify a particular analysis or type, then you are better off | |
530 using get_all_Genes or get_all_Genes_by_type and iterating | |
531 through the transcripts of each gene. | |
532 Returntype : reference to a list of Bio::EnsEMBL::Transcripts | |
533 Exceptions : none | |
534 Caller : general | |
535 | |
536 =cut | |
537 | |
538 sub get_all_Transcripts { | |
539 my $self = shift; | |
540 | |
541 my $transcripts = $self->SUPER::get_all_Transcripts(1); | |
542 $self->map_to_Individual($transcripts); | |
543 | |
544 return $transcripts; | |
545 } | |
546 | |
547 | |
548 =head2 get_all_Exons | |
549 | |
550 Arg [1] : (optional) string $dbtype | |
551 The dbtype of exons to obtain. This assumes that the db has | |
552 been added to the DBAdaptor under this name (using the | |
553 DBConnection::add_db_adaptor method). | |
554 Example : @exons = @{$individualSlice->get_all_Exons}; | |
555 Description: Gets all exons which overlap this IndividualSlice. Note that these exons | |
556 will not be associated with any transcripts, so this may not | |
557 be terribly useful. | |
558 Returntype : reference to a list of Bio::EnsEMBL::Exons | |
559 Exceptions : none | |
560 Caller : general | |
561 | |
562 =cut | |
563 | |
564 sub get_all_Exons { | |
565 my $self = shift; | |
566 my $dbtype = shift; | |
567 | |
568 my $exons = $self->SUPER::get_all_Exons($dbtype); | |
569 $self->map_to_Individual($exons); #map the exons to the Individual | |
570 | |
571 return $exons; | |
572 } | |
573 | |
574 =head2 get_all_Genes | |
575 | |
576 Arg [1] : (optional) string $logic_name | |
577 The name of the analysis used to generate the genes to retrieve | |
578 Arg [2] : (optional) string $dbtype | |
579 The dbtype of genes to obtain. This assumes that the db has | |
580 been added to the DBAdaptor under this name (using the | |
581 DBConnection::add_db_adaptor method). | |
582 Example : @genes = @{$individualSlice->get_all_Genes}; | |
583 Description: Retrieves all genes that overlap this slice. | |
584 Returntype : listref of Bio::EnsEMBL::Genes | |
585 Exceptions : none | |
586 Caller : none | |
587 | |
588 =cut | |
589 | |
590 sub get_all_Genes{ | |
591 my ($self, $logic_name, $dbtype) = @_; | |
592 | |
593 my $genes = $self->SUPER::get_all_Genes($logic_name, $dbtype, 1); | |
594 | |
595 $self->map_to_Individual($genes); | |
596 | |
597 foreach my $gene (@{$genes}){ | |
598 $self->map_to_Individual($gene->get_all_Exons); #map the Exons to the Individual | |
599 $self->map_to_Individual($gene->get_all_Transcripts); #map the Transcripts to the Individual | |
600 } | |
601 | |
602 return $genes; | |
603 } | |
604 | |
605 =head2 map_to_Individual | |
606 | |
607 Arg[1] : ref $features | |
608 Example : $individualSlice->map_to_Individual($exons); | |
609 Description : Gets the features from the Slice and maps it in the IndividualSlice, using the mapper | |
610 between Slice and IndividualSlice | |
611 ReturnType : None | |
612 Exceptions : None | |
613 Caller : general | |
614 | |
615 =cut | |
616 | |
617 sub map_to_Individual{ | |
618 my $self = shift; | |
619 my $features = shift; | |
620 | |
621 my $mapper = $self->mapper(); | |
622 my (@results, @results_ordered, $new_start, $new_end, $new_strand); | |
623 #foreach of the transcripts, map them to the IndividualSlice and replace the Slice with the IndividualSlice | |
624 foreach my $feature (@{$features}){ | |
625 $feature->slice($self); #replace the IndividualSlice as the Slice for this feature (the Slice plus the AlleleFeatures) | |
626 #map from the Slice to the Individual Slice | |
627 my @results = $mapper->map_coordinates('Slice',$feature->start,$feature->end,$feature->strand,'Slice'); | |
628 #from the results, order them but filter out those that are not coordinates | |
629 @results_ordered = sort {$a->start <=> $b->start} grep {ref($_) eq 'Bio::EnsEMBL::Mapper::Coordinate'} @results; | |
630 $new_start = $results_ordered[0]->start(); | |
631 $new_strand = $results_ordered[0]->strand(); | |
632 $new_end = $results_ordered[-1]->end(); #get last element of the array, the end of the slice | |
633 $feature->start($new_start); #update new coordinates | |
634 $feature->end($new_end); | |
635 $feature->strand($new_strand); | |
636 } | |
637 } | |
638 | |
639 sub alleleFeatures{ | |
640 my $self = shift; | |
641 return $self->{'alleleFeatures'}; | |
642 } | |
643 | |
644 sub add_AlleleFeature{ | |
645 my $self = shift; | |
646 | |
647 if (@_){ | |
648 if(!ref($_[0]) || !$_[0]->isa('Bio::EnsEMBL::Variation::AlleleFeature')) { | |
649 throw("Bio::EnsEMBL::Variation::AlleleFeature argument expected"); | |
650 } | |
651 #add the alleleFeature to the individualSlice | |
652 push @{$self->{'alleleFeatures'}},shift; | |
653 } | |
654 } | |
655 1; |