PairAlign.pm

Go to the documentation of this file.

=head1 LICENSE
 
Copyright [1999-2015] Wellcome Trust Sanger Institute and the EMBL-European Bioinformatics Institute
Copyright [2016-2024] EMBL-European Bioinformatics Institute
 
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
 
     http://www.apache.org/licenses/LICENSE-2.0
 
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
 
=cut
 
 
=head1 CONTACT
 
  Please email comments or questions to the public Ensembl
  developers list at <http://lists.ensembl.org/mailman/listinfo/dev>.
 
  Questions may also be sent to the Ensembl help desk at
  <http://www.ensembl.org/Help/Contact>.
 
=cut
 
=head1 NAME
 
PairAlign - Dna pairwise alignment module
 
=head1 SYNOPSIS
 
#To convert between coordinates:
 
  my $cdna_coord = $pair->genomic2cDNA($gen_coord);
  my $gen_coord  = $pair->cDNA2genomic($cdna_coord);
 
=head1 DESCRIPTION
 
Contains list of sub alignments making up a dna-dna alignment
 
Creation:
 
  my $pair = new Bio::EnsEMBL::FeaturePair(
    -start  => $qstart,
    -end    => $qend,
    -strand => $qstrand,
    -hstart => $hstart,
    -hend   => $hend,
    -hend   => $hstrand,
  );
 
  my $pairaln = new Bio::EnsEMBL::Analysis::PairAlign;
  $pairaln->addFeaturePair($pair);
 
Any number of pair alignments can be added to the PairAlign object
 
=cut
 
package Bio::EnsEMBL::Analysis::PairAlign;
 
use vars qw(@ISA);
use strict;
 
 
sub new {
    my($class,@args) = @_;
    my $self = {};
    bless $self, $class;
 
    $self->{'_homol'} = [];
    
    return $self; # success - we hope!
}
 
sub addFeaturePair {
    my ($self,$pair) = @_;
 
    $self->throw("Not a Bio::EnsEMBL::FeaturePair object") unless ($pair->isa("Bio::EnsEMBL::FeaturePair"));
 
    push(@{$self->{'_pairs'}},$pair);
    
}
 
 
=head2 eachFeaturePair
 
 Title   : eachFeaturePair 
 Example : my @pairs = $pair->eachFeaturePair
 Returns : Array of Bio::SeqFeature::FeaturePair
 Args    : none
=cut
 
sub eachFeaturePair {
    my ($self) = @_;
 
    if (defined($self->{'_pairs'})) {
    return @{$self->{'_pairs'}};
    }
}
 
sub get_hstrand {
    my ($self) = @_;
 
    my @features = $self->eachFeaturePair;
 
    return $features[0]->hstrand;
}
 
=head2 genomic2cDNA
 
 Title   : genomic2cDNA
 Usage   : my $cdna_coord = $pair->genomic2cDNA($gen_coord)
 Function: Converts a genomic coordinate to a cdna coordinate
 Returns : int
 Args    : int
=cut
 
sub genomic2cDNA {
  my ($self,$coord) = @_;
  my @pairs = $self->eachFeaturePair;
  
  @pairs = sort {$a->start <=> $b->start} @pairs;
  
  my $newcoord;
  
  HOMOL: while (my $sf1 = shift(@pairs)) {
    next HOMOL unless ($coord >= $sf1->start && $coord <= $sf1->end);
    
    if ($sf1->strand == 1 && $sf1->hstrand == 1) {
      $newcoord = $sf1->hstart + ($coord - $sf1->start);
      last HOMOL;
    } elsif ($sf1->strand == 1 && $sf1->hstrand == -1) {
      $newcoord = $sf1->hend   - ($coord - $sf1->start);
      last HOMOL;
    } elsif ($sf1->strand == -1 && $sf1->hstrand == 1) {
      $newcoord = $sf1->hstart + ($sf1->end - $coord);
      last HOMOL;
    } elsif ($sf1->strand == -1 && $sf1->hstrand == -1) {
      $newcoord = $sf1->hend   - ($sf1->end - $coord);
      last HOMOL;
    } else {
      $self->throw("ERROR: Wrong strand value in FeaturePair (" . $sf1->strand . "/" . $sf1->hstrand . "\n");
    }
  }    
  
  if (defined($newcoord)) {
    
    return $newcoord;
  } else {
    $self->throw("Couldn't convert $coord");
  }
}
 
=head2 cDNA2genomic
 
 Title   : cDNA2genomic
 Usage   : my $gen_coord = $pair->genomic2cDNA($cdna_coord)
 Function: Converts a cdna coordinate to a genomic coordinate
 Example : 
 Returns : int
 Args    : int
 
 
=cut
 
sub cDNA2genomic {
  my ($self,$coord) = @_;
  
  my @pairs = $self->eachFeaturePair;
  
  my $newcoord;
  
  HOMOL: while (my $sf1 = shift(@pairs)) {
    next HOMOL unless ($coord >= $sf1->hstart && $coord <= $sf1->hend);
    
    if ($sf1->strand == 1 && $sf1->hstrand == 1) {
      $newcoord = $sf1->start + ($coord - $sf1->hstart);
      last HOMOL;
    } elsif ($sf1->strand == 1 && $sf1->hstrand == -1) {
      $newcoord = $sf1->start  +($sf1->hend - $coord);
      last HOMOL;
    } elsif ($sf1->strand == -1 && $sf1->hstrand == 1) {
      $newcoord = $sf1->end   - ($coord - $sf1->hstart);
      last HOMOL;
    } elsif ($sf1->strand == -1 && $sf1->hstrand == -1) {
      $newcoord = $sf1->end   - ($sf1->hend - $coord);
      last HOMOL; 
    } else {
      $self->throw("ERROR: Wrong strand value in homol (" . $sf1->strand . "/" . $sf1->hstrand . "\n");
    }
  }
  
  if (defined ($newcoord)) {
    return $newcoord;
  } else {
    $self->throw("Couldn't convert $coord\n");
  }
}
 
sub find_Pair {
  my ($self,$coord) = @_;
  
  foreach my $p ($self->eachFeaturePair) {
    if ($coord >= $p->hstart && $coord <= $p->hend) {
      return $p;
    }
  }
}
 
=head2 convert_cDNA_feature
 
 Title   : convert_cDNA_feature
 Usage   : my @newfeatures = $self->convert_cDNA_feature($f);
 Function: Converts a feature on the cDNA into an array of 
           features on the genomic (for features that span across introns);
 Returns : Array of Bio::EnsEMBL::FeaturePair
 Args    : Bio::EnsEMBL::FeaturePair
 
=cut
 
sub convert_cDNA_feature {
  my ($self,$feature) = @_;
    
  my $foundstart = 0;
  my $foundend   = 0;
  
  my @pairs = $self->eachFeaturePair;
  my @newfeatures;
  
  HOMOL: while (my $sf1 = shift(@pairs)) {
    my $skip = 0;
    #print STDERR "Looking at cDNA exon " . $sf1->hstart . "\t" . $sf1->hend . "\t" . $sf1->strand ."\n";
    
    $skip = 1 unless ($feature->start >= $sf1->hstart 
                      && $feature->start <= $sf1->hend);
    
    if($skip){
      #print STDERR "Skipping ".$sf1->hstart . "\t" . $sf1->hend . "\t" . $sf1->strand ."\n";
      next HOMOL;
    }
    if ($feature->end >= $sf1->hstart && $feature->end <= $sf1->hend) {
      $foundend = 1;
    }
    
    my $startcoord = $self->cDNA2genomic($feature->start);
    my $endcoord;
    
    if ($sf1->hstrand == 1) {
      $endcoord   = $sf1->end;
    } else {
      $endcoord   = $sf1->start;
    }
    
    if ($foundend) {
      $endcoord = $self->cDNA2genomic($feature->end);
    }
      
    #print STDERR "Making new genomic feature $startcoord\t$endcoord\n";
    
    my $tmpf = new Bio::EnsEMBL::Feature(-seqname => $feature->seqname,
                                            -start   => $startcoord,
                                            -end     => $endcoord,
                                            -strand  => $feature->strand);
    push(@newfeatures,$tmpf);
    last;
  }
    
  # Now the rest of the pairs until we find the endcoord
  
  while ((my $sf1 = shift(@pairs)) && ($foundend == 0)) {
    
    if ($feature->end >= $sf1->hstart && $feature->end <= $sf1->hend) {
      $foundend = 1;
    }
    
    my $startcoord;
    my $endcoord;
    
    if ($sf1->hstrand == 1) {
      $startcoord = $sf1->start;
      $endcoord   = $sf1->end;
    } else {
      $startcoord = $sf1->end;
      $endcoord   = $sf1->start;
    }
    
    if ($foundend) {
      $endcoord = $self->cDNA2genomic($feature->end);
    }
    
    #   #print STDERR "Making new genomic feature $startcoord\t$endcoord\n";
    
    my $tmpf = new Bio::EnsEMBL::Feature(-seqname => $feature->seqname,
                                         -start   => $startcoord,
                                         -end     => $endcoord,
                                         -strand  => $feature->strand);
    push(@newfeatures,$tmpf);
  }
  #print STDERR "Have ".@newfeatures." features from ".$feature."\n";
  return @newfeatures;
}
 
 
sub convert_FeaturePair {
  my ($self,$pair) = @_;
 
  my $hstrand = $self->get_hstrand;
  
  my $feat = $self->create_Feature($pair->start, $pair->end, $pair->strand,
                                   $pair->slice);
  my @newfeatures = $self->convert_cDNA_feature($feat);
  my @newpairs;
  
  my $hitpairaln  = new Bio::EnsEMBL::Analysis::PairAlign;
  $hitpairaln->addFeaturePair($pair);
  
  foreach my $new (@newfeatures) {
    
    # Now we want to convert these cDNA coords into hit coords
    
    my $hstart1 = $self->genomic2cDNA($new->start);
    my $hend1   = $self->genomic2cDNA($new->end);
    
    my $hstart2 = $hitpairaln->genomic2cDNA($hstart1);
    my $hend2   = $hitpairaln->genomic2cDNA($hend1);
    
    # We can now put the final feature together
    
    my $finalstrand = $hstrand * $pair->strand * $pair->hstrand;
    
    if ($hstart2 > $hend2) {
      my $tmp = $hstart2;
      $hstart2 = $hend2;
      $hend2   = $tmp;
    }
    
    my $finalpair = $self->create_FeaturePair($new->start, $new->end, 
                                              $new->strand,
                                              $hstart2, $hend2, 
                                              $finalstrand, $pair->score);
    
    push(@newpairs,$finalpair);
    
  }
  
  return @newpairs;
}
 
sub create_FeaturePair {
    my ($self, $start, $end, $strand, $hstart, $hend, 
        $hstrand, $score) = @_;
   
    my $fp = Bio::EnsEMBL::FeaturePair->new(
                                            -start    => $start,
                                            -end      => $end,
                                            -strand   => $strand,
                                            -hstart   => $hstart,
                                            -hend     => $hend,
                                            -hstrand  => $hstrand,
                                            -score    => $score,
                                           );
 
   
    return $fp;
}
 
sub create_Feature{
  my ($self, $start, $end, $strand,  $slice) = @_;
 
  my $feat = new Bio::EnsEMBL::Feature(-start   => $start,
                                       -end     => $end,
                                       -strand  => $strand,
                                       -slice   => $slice,
                                      );
  return $feat;
}
 
1;