IndividualSlice.pm 23.7 KB
Newer Older
1
=head1 LICENSE
Daniel Rios's avatar
Daniel Rios committed
2

3
  Copyright (c) 1999-2013 The European Bioinformatics Institute and
4
  Genome Research Limited.  All rights reserved.
Daniel Rios's avatar
Daniel Rios committed
5

6 7
  This software is distributed under a modified Apache license.
  For license details, please see
Daniel Rios's avatar
Daniel Rios committed
8

9
    http://www.ensembl.org/info/about/code_licence.html
Daniel Rios's avatar
Daniel Rios committed
10

11
=head1 CONTACT
Daniel Rios's avatar
Daniel Rios committed
12

13
  Please email comments or questions to the public Ensembl
14
  developers list at <dev@ensembl.org>.
Daniel Rios's avatar
Daniel Rios committed
15

16 17
  Questions may also be sent to the Ensembl help desk at
  <helpdesk@ensembl.org>.
Daniel Rios's avatar
Daniel Rios committed
18

19
=cut
Daniel Rios's avatar
Daniel Rios committed
20

21
=head1 NAME
Daniel Rios's avatar
Daniel Rios committed
22

23 24 25
Bio::EnsEMBL::IndividualSlice - SubClass of the Slice. Represents the
slice of the genome for a certain individual (applying the alleles for
this individual)
Daniel Rios's avatar
Daniel Rios committed
26

27
=head1 SYNOPSIS
Daniel Rios's avatar
Daniel Rios committed
28

29
  $sa = $db->get_SliceAdaptor;
Daniel Rios's avatar
Daniel Rios committed
30

31 32
  $slice =
    $sa->fetch_by_region( 'chromosome', 'X', 1_000_000, 2_000_000 );
Daniel Rios's avatar
Daniel Rios committed
33

34 35 36 37 38 39 40 41 42 43
  $individualSlice = $slice->get_by_Individual($individual_name);

  # Get the sequence from the Individual Slice: will contain IUPAC codes
  # for SNPs and Ensembl ambiguity codes for indels
  my $seq = $individualSlice->seq();
  print $seq;

  # Get a subSlice of the Strain
  my $subSlice_individual =
    $individualSlice->sub_Slice( 5_000, 8_000, 1 )
Daniel Rios's avatar
Daniel Rios committed
44

45 46 47 48 49 50 51 52 53 54 55 56 57 58 59
  # Compare two different individuals in the same Slice
  my $sliceIndividual2 = $slice->get_by_Individual($individual_name2);
  my $differences =
    $individualSlice->get_all_differences_IndividualSlice(
    $sliceIndividual2);

  foreach my $af ( @{$differences} ) {
    print
      "There is a difference between $individual_name "
      . "and $individual_name2 at ",
      $af->start, "-", $af->end,
      " with allele ", $af->allele_string(), "\n";
  }

=head1 DESCRIPTION
Daniel Rios's avatar
Daniel Rios committed
60

61 62 63
A IndividualSlice object represents a region of a genome for a certain
individual.  It can be used to retrieve sequence or features from a
individual.
Daniel Rios's avatar
Daniel Rios committed
64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114

=head1 METHODS

=cut

package Bio::EnsEMBL::IndividualSlice;
use vars qw(@ISA);
use strict;

use Bio::EnsEMBL::Utils::Argument qw(rearrange);
use Bio::EnsEMBL::Utils::Sequence qw(reverse_comp);
use Bio::EnsEMBL::Slice;
use Bio::EnsEMBL::Mapper;
use Bio::EnsEMBL::Utils::Exception qw(throw deprecate warning);

@ISA = qw(Bio::EnsEMBL::Slice);


=head2 new

    Arg [1..N]  : List of named arguments
                  Bio::EnsEMBL::CoordSystem COORD_SYSTEM
                  string SEQ_REGION_NAME,
                  int    START,
                  int    END,
                  string VERSION (optional, defaults to '')
                  int    STRAND, (optional, defaults to 1)
                  Bio::EnsEMBL::DBSQL::SliceAdaptor ADAPTOR (optional)
    Arg[N+1]      : string $individual_name
    Example     : $individualSlice = Bio::EnsEMBL::IndividualSlice->new(-coord_system => $cs,
									-start => 1,
									-end => 10000,
									-strand => 1,
									-seq_region_name => 'X',
									-seq_region_length => 12e6,
									-individual_name => $individual_name);
    Description : Creates a new Bio::EnsEMBL::IndividualSlice object that will contain a shallow copy of the
                  Slice object, plus additional information such as the individual this Slice refers to
                  and listref of Bio::EnsEMBL::Variation::AlleleFeatures of differences with the
                  reference sequence
    ReturnType  : Bio::EnsEMBL::IndividualSlice
    Exceptions  : none
    Caller      : general

=cut

sub new{
    my $caller = shift;
    my $class = ref($caller) || $caller;

    #create the IndividualSlice object as the Slice, plus the individual attribute
Daniel Rios's avatar
Daniel Rios committed
115
    my ($individual_name, $sample_id) = rearrange(['INDIVIDUAL', 'SAMPLE_ID'],@_);
Daniel Rios's avatar
Daniel Rios committed
116 117 118 119

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

    $self->{'individual_name'} = $individual_name;
Daniel Rios's avatar
Daniel Rios committed
120
    $self->{'sample_id'} = $sample_id;
Daniel Rios's avatar
Daniel Rios committed
121

Daniel Rios's avatar
Daniel Rios committed
122
    return $self;
Daniel Rios's avatar
Daniel Rios committed
123 124 125 126 127 128 129 130 131 132 133 134 135 136

}

=head2 individual_name

    Arg [1]     : (optional) string $individual_name
    Example     : my $individual_name = $individualSlice->individual_name();
    Description : Getter/Setter for the name of the individual in the slice
    ReturnType  : string
    Exceptions  : none
    Caller      : general

=cut

Daniel Rios's avatar
Daniel Rios committed
137
sub individual_name{
Daniel Rios's avatar
Daniel Rios committed
138 139 140 141 142 143
   my $self = shift;
   if (@_){
       $self->{'individual_name'} = shift @_;
   }
   return $self->{'individual_name'};
}
Daniel Rios's avatar
Daniel Rios committed
144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171

=head2 seq

  Arg [1]    : none
  Example    : print "SEQUENCE = ", $strainSlice->seq();
  Description: Returns the sequence of the region represented by this
               StrainSlice formatted as a string.
  Returntype : string
  Exceptions : none
  Caller     : general

=cut

sub seq {
  my $self = shift;

  # special case for in-between (insert) coordinates
  return '' if($self->start() == $self->end() + 1);

  return $self->{'seq'} if($self->{'seq'});

  if($self->adaptor()) {
    my $seqAdaptor = $self->adaptor()->db()->get_SequenceAdaptor();
    my $reference_sequence = $seqAdaptor->fetch_by_Slice_start_end_strand($self,1,undef,1); #get the reference sequence for that slice
    #apply all differences to the reference sequence

    # sort edits in reverse order to remove complication of
    # adjusting downstream edits
172 173
    my @allele_features_ordered;
    @allele_features_ordered = sort {$b->start() <=> $a->start() || $b->end() <=> $a->end()} @{$self->{'alleleFeatures'}} if (defined $self->{'alleleFeatures'});
Daniel Rios's avatar
Daniel Rios committed
174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346

    foreach my $af (@allele_features_ordered){
	$af->apply_edit($reference_sequence); #change, in the reference sequence, the af
    }
#    return substr(${$reference_sequence},0,1) if ($self->length == 1); 
    return ${$reference_sequence}; #returns the reference sequence, applying the alleleFeatures
  }

  # no attached sequence, and no db, so just return Ns
  return 'N' x $self->length();
}

=head2 get_all_differences_Slice

    Args        : none
    Example     : my $differences = $individualSlice->get_all_differences_Slice()
    Description : Gets all differences between the IndividualSlice object and the Slice is defined
    ReturnType  : listref of Bio::EnsEMBL::Variation::AlleleFeature
    Exceptions  : none
    Caller      : general

=cut

sub get_all_differences_Slice{
    my $self = shift;
    my $differences; #reference to the array with the differences between Slice and StrainSlice
    my $ref_allele;
    foreach my $difference (@{$self->{'alleleFeatures'}}){
	if ($difference->length_diff == 0){
	    #the difference is a SNP, check if it is the same as the reference allele
	    $ref_allele = $self->SUPER::subseq($difference->start,$difference->end,$difference->strand);
	    $ref_allele = '-' if ($ref_allele eq '');
	    if ($ref_allele ne $difference->allele_string){
		#when the alleleFeature is different from the reference allele, add to the differences list
		push @{$differences},$difference;
	    }
	}
	else{
	    push @{$differences},$difference;
	}
    }

    return $differences;

}

=head2 get_all_differences_IndividualSlice

    Arg[1]      : Bio::EnsEMBL::IndividualSlice $is
    Example     : my $differences = $individualSlice->get_all_differences_IndividualSlice($individualslice)
    Description : Gets differences between 2 IndividualSlice objects
    ReturnType  : listref of Bio::EnsEMBL::Variation::AlleleFeature
    Exceptions  : thrown on bad argument
    Caller      : general

=cut

sub get_all_differences_IndividualSlice{
    my $self = shift;
    my $individualSlice = shift;

    if (!ref($individualSlice) || !$individualSlice->isa('Bio::EnsEMBL::IndividualSlice')){
	throw('Bio::EnsEMBL::IndividualSlice arg expected');
    }
    if ( @{$self->{'alleleFeatures'}} == 0 && @{$individualSlice->{'alleleFeatures'}} == 0){
	return undef; #there are no differences in any of the Individuals
	
    }
    my $differences; #differences between individuals
    if (@{$individualSlice->{'alleleFeatures'}} == 0){
	#need to create a copy of alleleFeature for the first Individual
	foreach my $difference (@{$self->{'alleleFeatures'}}){
	    my %vf = %$difference;
	    push @{$differences},bless \%vf,ref($difference);
	}
    }
    elsif (@{$self->{'alleleFeatures'}} == 0){
	#need to create a copy of AlleleFeature, but changing the allele by the allele in the reference sequence
	foreach my $difference (@{$individualSlice->{'alleleFeatures'}}){
	    push @{$differences}, $individualSlice->_convert_difference($difference);
	}
    }
    else{
	#both individuals have differences
	#create a hash with the differences in the first slice
	my %allele_features_self = map {$_->start.'-'.$_->end => $_} @{$self->{'alleleFeatures'}};
	foreach my $difference (@{$individualSlice->{'alleleFeatures'}}){
	    #there is no difference in the other individual slice, convert the allele
	  if (!defined $allele_features_self{$difference->start.'-'.$difference->end}){
	      push @{$differences},$individualSlice->_convert_difference($difference);
	    }
	    else{
		#if it is defined and have the same allele, delete from the hash since it is not a difference
		#between the individuals
		if ($allele_features_self{$difference->start.'-'.$difference->end}->allele_string eq $difference->allele_string){
		  delete $allele_features_self{$difference->start.'-'.$difference->end};
		}
	    }
	}	
	#and finally, make a shallow copy of the differences in the first individual
	foreach my $difference (values %allele_features_self){
	  my %vf = %$difference;
	  push @{$differences},bless \%vf,ref($difference);
	}
	
    }
    #need to map differences to the first individual, self, since the coordinates are in the Slice coordinate system
    my $mapper = $self->mapper(); #now that we have the differences, map them in the IndividualSlice
    my @results;
    foreach my $difference (@{$differences}){
      @results = $mapper->map_coordinates('Slice',$difference->start,$difference->end,$difference->strand,'Slice');
      #we can have 3 possibilities:
      #the difference is an insertion and when mapping returns the boundaries of the insertion in the IndividualSlice
      if (@results == 2){
	#the first position in the result is the beginning of the insertion
	if($results[0]->start < $results[1]->start){
	  $difference->start($results[0]->end+1);
	  $difference->end($results[1]->start-1);
	}
	else{
	  #it is the second position the beginning of the insertion
	  $difference->start($results[1]->end+1);
	  $difference->end($results[0]->start-1);
	}
	$difference->strand($results[0]->strand);
      }
      else{
	#it can be either a SNP or a deletion, and we have the coordinates in the result, etither a Bio::EnsEMBL::Mapper::Coordinate
	# or a Bio::EnsEMBL::Mapper::IndelCoordinate
	$difference->start($results[0]->start);
	$difference->end($results[0]->end);
	$difference->strand($results[0]->strand);
      }
    }
    
    return $differences;
}

#for a given AlleleFeature, converts the allele into the reference allele and returns
#the converted AlleleFeature

sub _convert_difference{
    my $self = shift;
    my $difference = shift;
    my %new_af = %$difference; #make a copy of the alleleFeature
    #and change the allele with the one from the reference Slice
    $new_af{'allele_string'} = $self->SUPER::subseq($difference->start,$difference->end,$difference->strand);
    return bless \%new_af,ref($difference);
}

=head2 mapper

  Args       : none
  Description: Getter for the mapper between the between the IndividualSlice and the Slice it refers to. 
               It is done automatically when necessary to create subSlice or to get the differences between individuals
  Returntype : Bio::EnsEMBL::Mapper
  Exceptions : none
  Caller     : Internal function

=cut

sub mapper{
    my $self = shift;
   
    if (@_) {
	#allow to create again the mapper
	delete $self->{'mapper'};
    }
    if(!defined $self->{'mapper'}){
	#create the mapper between the Slice and StrainSlice
	my $mapper = Bio::EnsEMBL::Mapper->new('Slice','IndividualSlice');
	#align with Slice
	#get all the VariationFeatures in the Individual Slice, from start to end in the Slice
347 348
  my @allele_features_ordered;
	@allele_features_ordered = sort {$a->start() <=> $b->start() || $b->end() <=> $a->end()} @{$self->{'alleleFeatures'}} if (defined $self->{'alleleFeatures'});
Daniel Rios's avatar
Daniel Rios committed
349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469
	
	my $start_slice = 1;
	my $end_slice;
	my $start_individual = 1;
	my $end_individual;
	my $length_allele;
	my $total_length_diff = 0;
	#we will walk from left to right in the slice object, updating the start and end individual every time
	#there is a new alleleFeature in the Individual
	foreach my $allele_feature (@allele_features_ordered){
	    #we have a insertion/deletion: marks the beginning of new slice move coordinates
	    if ($allele_feature->length_diff != 0){
	        $total_length_diff += $allele_feature->length_diff;
		$length_allele = $allele_feature->length + $allele_feature->length_diff(); #length of the allele in the Individual
		$end_slice = $allele_feature->start() - 1; #set the end of the slice before the alleleFeature
		if ($end_slice >= $start_slice){
		    #normal cases (not with gaps)
		    $end_individual = $end_slice - $start_slice + $start_individual; #set the end of the individual from the beginning plus the offset
		    #add the sequence that maps
		    $mapper->add_map_coordinates('Slice',$start_slice,$end_slice,1,'IndividualSlice',$start_individual,$end_individual);
		    #and add the indel
		    $mapper->add_indel_coordinates('Slice',$end_slice+1,$end_slice + $allele_feature->length,1,'IndividualSlice',$end_individual+1,$end_individual + $length_allele);
		    $start_individual = $end_individual + $length_allele + 1; #set the beginning of the individual after the allele
		}
		else{
		    #add the indel
		    $mapper->add_indel_coordinates('Slice',$end_slice+1,$end_slice + $allele_feature->length,1,'IndividualSlice',$end_individual+1,$end_individual + $length_allele);
		    $start_individual += $length_allele;
		}
		$start_slice = $end_slice + $allele_feature->length+ 1; #set the beginning of the slice after the variation feature
	    }
	}
	if ($start_slice <= $self->length){
	    #if we haven't reached the end of the IndividualSlice, add the final map coordinates between the individual and the slice
	    $mapper->add_map_coordinates('Slice',$start_slice,$self->length,1,'IndividualSlice',$start_individual,$start_individual + $self->length - $start_slice);
	}

	$mapper->add_map_coordinates('Slice', -$self->start+1, 0,1, 'IndividualSlice', -$self->start +1,0) if ($self->start > 0); #before individualSlice
	$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
	$self->{'mapper'} = $mapper;
    }
    return $self->{'mapper'};
}

=head2 sub_Slice

  Arg   1    : int $start
  Arg   2    : int $end
  Arge [3]   : int $strand
  Example    : none
  Description: Makes another IndividualSlice that covers only part of this IndividualSlice
               with the appropriate differences to the reference Slice
               If a slice is requested which lies outside of the boundaries
               of this function will return undef.  This means that
               behaviour will be consistant whether or not the slice is
               attached to the database (i.e. if there is attached sequence
               to the slice).  Alternatively the expand() method or the
               SliceAdaptor::fetch_by_region method can be used instead.
  Returntype : Bio::EnsEMBL::IndividualSlice or undef if arguments are wrong
  Exceptions : thrown when trying to get the subSlice in the middle of a
               insertion
  Caller     : general

=cut

sub sub_Slice {
  my ( $self, $start, $end, $strand ) = @_;
  my $mapper = $self->mapper();
  #map from the Individual to the Slice to get the sub_Slice, and then, apply the differences in the subSlice
  my @results = $mapper->map_coordinates('IndividualSlice',$start,$end,$strand,'IndividualSlice');
  my $new_start;
  my $new_end;
  my $new_strand;
  my $new_seq;
  #Get need start and end for the subSlice of the IndividualSlice
  my @results_ordered = sort {$a->start <=> $b->start} grep {ref($_) eq 'Bio::EnsEMBL::Mapper::Coordinate'} @results;
  $new_start = $results_ordered[0]->start();
  $new_strand = $results_ordered[0]->strand() if (ref($results_ordered[0]) eq 'Bio::EnsEMBL::Mapper::Coordinate');
#  $new_strand = $results_ordered[-1]->strand() if (ref($results_ordered[-1]) eq 'Bio::EnsEMBL::Mapper::Coordinate');
  $new_end = $results_ordered[-1]->end();  #get last element of the array, the end of the slice

  my $subSlice = $self->SUPER::sub_Slice($new_start,$new_end,$new_strand);
  $subSlice->{'individual_name'} = $self->{'individual_name'};

  my $new_alleles; #reference to an array that will contain the variationFeatures in the new subSlice
  #update the VariationFeatures in the sub_Slice of the Individual
  my %af;
  my $new_allele_feature;
  foreach my $alleleFeature (@{$self->{'alleleFeatures'}}){      
      $new_allele_feature = $alleleFeature->transfer($subSlice);
      #only transfer the coordinates to the SubSlice that are within the boundaries
      if ($new_allele_feature->start >= 1 && $new_allele_feature->end <= $subSlice->length){
	  push @{$new_alleles}, $new_allele_feature;
      }
  }
  $subSlice->{'alleleFeatures'} = $new_alleles;
  return $subSlice;

}

=head2 subseq

  Arg  [1]   : int $startBasePair
               relative to start of slice, which is 1.
  Arg  [2]   : int $endBasePair
               relative to start of slice.
  Arg  [3]   : (optional) int $strand
               The strand of the individual slice to obtain sequence from. Default
               value is 1.
  Description: returns string of dna sequence
  Returntype : txt
  Exceptions : end should be at least as big as start
               strand must be set
  Caller     : general

=cut

sub subseq {
  my ( $self, $start, $end, $strand ) = @_;

  if ( $end+1 < $start ) {
470
    throw("End coord + 1 is less than start coord");
Daniel Rios's avatar
Daniel Rios committed
471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489
  }

  # handle 'between' case for insertions
  return '' if( $start == $end + 1);

  $strand = 1 unless(defined $strand);

  if ( $strand != -1 && $strand != 1 ) {
    throw("Invalid strand [$strand] in call to Slice::subseq.");
  }

  my $subseq;
  my $seq;
  if($self->adaptor){
    my $seqAdaptor = $self->adaptor()->db()->get_SequenceAdaptor();
    $subseq = ${$seqAdaptor->fetch_by_Slice_start_end_strand($self,$start,$end,$strand)}; #get the reference sequence for that slice
    #apply all differences to the reference sequence
    # sort edits in reverse order to remove complication of
    # adjusting downstream edits
490 491
    my @allele_features_ordered;
    @allele_features_ordered = sort {$b->start() <=> $a->start() || $b->end() <=> $a->end()} @{$self->{'alleleFeatures'}} if (defined $self->{'alleleFeatures'});
Daniel Rios's avatar
Daniel Rios committed
492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658
    my $af_start;
    my $af_end;
    foreach my $af (@allele_features_ordered){
	if (($af->start - $start +1 > 0) && ($end - $af->end > 0)){
	    #save the current start and end of the alleleFeature before changing for apply_edit
	    $af_start = $af->start;
	    $af_end = $af->end;
	    #apply the difference if the feature is in the new slice
	    $af->start($af->start - $start +1);
	    $af->end($af->end - $start +1);
	    $af->apply_edit(\$subseq); #change, in the reference sequence, the af
	    #restore the initial values of alleleFeature start and end
	    $af->start($af_start);
	    $af->end($af_end);
	    
	}
    }
  } 
  else {
      ## check for gap at the beginning and pad it with Ns
      if ($start < 1) {
	  $subseq = "N" x (1 - $start);
	  $start = 1;
      }
      $subseq .= substr ($self->seq(), $start-1, $end - $start + 1);
      ## check for gap at the end and pad it with Ns
    if ($end > $self->length()) {
	$subseq .= "N" x ($end - $self->length());
    }
      reverse_comp(\$subseq) if($strand == -1);
  }
  return $subseq;
  
}

=head2 get_all_Transcripts

  Args       : None
  Example    : @transcripts = @{$individualslice->get_all_Transcripts)};
  Description: Gets all transcripts which overlap this Individual Slice.  If you want to
               specify a particular analysis or type, then you are better off
               using get_all_Genes or get_all_Genes_by_type and iterating
               through the transcripts of each gene.
  Returntype : reference to a list of Bio::EnsEMBL::Transcripts
  Exceptions : none
  Caller     : general

=cut

sub get_all_Transcripts {
  my $self = shift;

  my $transcripts = $self->SUPER::get_all_Transcripts(1);
  $self->map_to_Individual($transcripts);

  return $transcripts;
}


=head2 get_all_Exons

  Arg [1]    : (optional) string $dbtype
               The dbtype of exons to obtain.  This assumes that the db has
               been added to the DBAdaptor under this name (using the
               DBConnection::add_db_adaptor method).
  Example    : @exons = @{$individualSlice->get_all_Exons};
  Description: Gets all exons which overlap this IndividualSlice.  Note that these exons
               will not be associated with any transcripts, so this may not
               be terribly useful.
  Returntype : reference to a list of Bio::EnsEMBL::Exons
  Exceptions : none
  Caller     : general

=cut

sub get_all_Exons {
  my $self = shift;
  my $dbtype = shift;

  my $exons = $self->SUPER::get_all_Exons($dbtype);
  $self->map_to_Individual($exons); #map the exons to the Individual

  return $exons;
}

=head2 get_all_Genes

  Arg [1]    : (optional) string $logic_name
               The name of the analysis used to generate the genes to retrieve
  Arg [2]    : (optional) string $dbtype
               The dbtype of genes to obtain.  This assumes that the db has
               been added to the DBAdaptor under this name (using the
               DBConnection::add_db_adaptor method).
  Example    : @genes = @{$individualSlice->get_all_Genes};
  Description: Retrieves all genes that overlap this slice.
  Returntype : listref of Bio::EnsEMBL::Genes
  Exceptions : none
  Caller     : none

=cut

sub get_all_Genes{
  my ($self, $logic_name, $dbtype) = @_;

  my $genes = $self->SUPER::get_all_Genes($logic_name, $dbtype, 1);

  $self->map_to_Individual($genes);

  foreach my $gene (@{$genes}){
      $self->map_to_Individual($gene->get_all_Exons); #map the Exons to the Individual
      $self->map_to_Individual($gene->get_all_Transcripts); #map the Transcripts to the Individual
  }

  return $genes;
}

=head2 map_to_Individual

    Arg[1]      : ref $features
    Example     : $individualSlice->map_to_Individual($exons);
    Description : Gets the features from the Slice and maps it in the IndividualSlice, using the mapper
                  between Slice and IndividualSlice
    ReturnType  : None
    Exceptions  : None
    Caller      : general
    
=cut

sub map_to_Individual{
    my $self = shift;
    my $features = shift;

    my $mapper = $self->mapper();
    my (@results, @results_ordered, $new_start, $new_end, $new_strand);
    #foreach of the transcripts, map them to the IndividualSlice and replace the Slice with the IndividualSlice
    foreach my $feature (@{$features}){
	$feature->slice($self); #replace the IndividualSlice as the Slice for this feature (the Slice plus the AlleleFeatures)
	#map from the Slice to the Individual Slice
	my @results = $mapper->map_coordinates('Slice',$feature->start,$feature->end,$feature->strand,'Slice');
	#from the results, order them but filter out those that are not coordinates
	@results_ordered = sort {$a->start <=> $b->start} grep {ref($_) eq 'Bio::EnsEMBL::Mapper::Coordinate'} @results;
	$new_start = $results_ordered[0]->start();
	$new_strand = $results_ordered[0]->strand();
	$new_end = $results_ordered[-1]->end();  #get last element of the array, the end of the slice
	$feature->start($new_start);  #update new coordinates
	$feature->end($new_end);
	$feature->strand($new_strand);
  }
}

sub alleleFeatures{
    my $self = shift;
    return $self->{'alleleFeatures'};
}

sub add_AlleleFeature{
    my $self = shift;

    if (@_){
	if(!ref($_[0]) || !$_[0]->isa('Bio::EnsEMBL::Variation::AlleleFeature')) {
	    throw("Bio::EnsEMBL::Variation::AlleleFeature argument expected");
	}
	#add the alleleFeature to the individualSlice
	push @{$self->{'alleleFeatures'}},shift;
    }
}
1;