Raw content of Bio::EnsEMBL::Analysis::PairAlign
<![CDATA[=head1 LICENSE
Copyright (c) 1999-2009 The European Bioinformatics Institute and
Genome Research Limited. All rights reserved.
This software is distributed under a modified Apache license.
For license details, please see
/info/about/code_licence.html
=head1 CONTACT
Please email comments or questions to the public Ensembl
developers list at <ensembl-dev@ebi.ac.uk>.
Questions may also be sent to the Ensembl help desk at
<helpdesk@ensembl.org>.
=cut
# POD documentation - main docs before the code
=head1 NAME
PairAlign - Dna pairwise alignment module
=head1 SYNOPSIS
Give standard usage here
=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
Manipulation:
To convert between coordinates:
my $cdna_coord = $pair->genomic2cDNA($gen_coord);
my $gen_coord = $pair->cDNA2genomic($cdna_coord);
=head1 METHODS
=cut
package Bio::EnsEMBL::Analysis::PairAlign;
use vars qw(@ISA);
use strict;
@ISA = qw(Bio::EnsEMBL::Root);
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 : eachFeaturePait
Usage : my @pairs = $pair->eachFeaturePair
Function:
Example :
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
Example :
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);
Example :
Returns : @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;
]]>