Skip to content
Snippets Groups Projects
Commit 048b9f36 authored by juguang's avatar juguang
Browse files

converter from GenericHSP to BaseAlignFeature's submodules, for blast report

GenericHPS::algorithm tells what program of blastall is used, and in response of this, the code try to find the specific alignment feature in EnsEMBL.
parent 1ec434ad
No related branches found
No related tags found
No related merge requests found
=head1
Sequence alignment hits were previously stored within the core database as
ungapped alignments. This imposed 2 major constraints on alignments:
a) alignments for a single hit record would require multiple rows in the
database, and
b) it was not possible to accurately retrieve the exact original alignment.
Therefore, in the new branch sequence alignments are now stored as ungapped
alignments in the cigar line format (where CIGAR stands for Concise
Idiosyncratic Gapped Alignment Report).
In the cigar line format alignments are sotred as follows:
M: Match
D: Deletino
I: Insertion
An example of an alignment for a hypthetical protein match is shown below:
Query: 42 PGPAGLP----GSVGLQGPRGLRGPLP-GPLGPPL...
PG P G GP R PLGP
Sbjct: 1672 PGTP*TPLVPLGPWVPLGPSSPR--LPSGPLGPTD...
protein_align_feature table as the following cigar line:
7M4D12M2I2MD7M
=cut
package Bio::EnsEMBL::Utils::Converter::bio_ens_hit;
use strict;
use vars qw(@ISA);
use Bio::EnsEMBL::Utils::Converter::bio_ens;
use Bio::EnsEMBL::DnaDnaAlignFeature;
use Bio::EnsEMBL::DnaPepAlignFeature;
use Bio::EnsEMBL::PepDnaAlignFeature;
@ISA = qw(Bio::EnsEMBL::Utils::Converter::bio_ens);
sub _initialize {
my ($self, @args) = @_;
$self->SUPER::_initialize(@args);
# After super initialized, analysis and contig are ready.
my $bio_ens_seqFeature_converter = new Bio::EnsEMBL::Utils::Converter(
-in => 'Bio::SeqFeature::Generic',
-out => 'Bio::EnsEMBL::SeqFeature',
-analysis => $self->analysis,
-contig => $self->contig
);
$self->_bio_ens_seqFeature_converter($bio_ens_seqFeature_converter);
}
sub _convert_single {
my ($self, $input) = @_;
my $in = $self->in;
my $out = $self->out;
if($in =~ /Bio::Search::Hit::GenericHit/){
return $self->_convert_single_hit($input);
}elsif($in =~ /Bio::Search::HSP::GenericHSP/){
return $self->_convert_single_hsp($input);
}else{
$self->throw("[$in]->[$out], not implemented");
}
}
sub _convert_single_hit {
}
sub _convert_single_hsp {
my ($self, $hsp) = @_;
unless(ref($hsp) && $hsp->isa('Bio::Search::HSP::GenericHSP')){
$self->throw("a GenericHSP object needed");
}
my $bio_ens_seqFeature_converter = $self->_bio_ens_seqFeature_converter;
my $ens_feature1 = $bio_ens_seqFeature_converter->_convert_single(
$hsp->feature1);
my $ens_feature2 = $bio_ens_seqFeature_converter->_convert_single(
$hsp->feature2);
$ens_feature1->p_value($hsp->evalue);
$ens_feature1->score($hsp->score);
$ens_feature1->percent_id($hsp->percent_identity);
$ens_feature2->p_value($hsp->evalue);
$ens_feature2->score($hsp->score);
$ens_feature2->percent_id($hsp->percent_identity);
my $cigar_string = $hsp->cigar_string;
my @args = (
-feature1 => $ens_feature1,
-feature2 => $ens_feature2,
-cigar_string => $cigar_string
);
my $contig = $self->contig;
# choose the AlignFeature based on the blast program
my $program = $hsp->algorithm;
$self->throw("HSP does not have algorithm value") unless(defined($program));
my $align_feature;
if($program =~ /blastn/i){
$align_feature = new Bio::EnsEMBL::DnaDnaAlignFeature(@args);
$align_feature->attach_seq($contig);
}elsif($program =~ /blastx/i){
$align_feature = new Bio::EnsEMBL::DnaPepAlignFeature(@args);
$align_feature->attach_seq($contig);
}else{
$self->throw("$program is not supported yet");
}
return $align_feature;
}
# an internal getter/setter for a converter used for seq feature conversion.
sub _bio_ens_seqFeature_converter {
my ($self, $arg) = @_;
if(defined $arg){
$self->{_bio_ens_seqFeature_converter} = $arg;
}
return $self->{_bio_ens_seqFeature_converter};
}
1;
This diff is collapsed.
......@@ -4,11 +4,9 @@ use strict;
BEGIN {
$| = 1;
use Test;
plan tests => 22;
plan tests => 36;
}
use Bio::EnsEMBL::Utils::Converter;
END {
}
......@@ -49,6 +47,8 @@ my $featurePair1 = new Bio::SeqFeature::FeaturePair(
-feature2 => $seqFeature2
);
use Bio::EnsEMBL::Utils::Converter;
&test_SeqFeature;
&test_FeaturePair;
&test_hit;
......@@ -125,9 +125,48 @@ sub test_FeaturePair{
}
# 14 OKs
sub test_hit {
my $converter = new Bio::EnsEMBL::Utils::Converter(
-in => 'Bio::Search::HSP::GenericHSP',
-out => 'Bio::EnsEMBL::BaseAlignFeature',
-analysis => $ens_analysis,
-contig => $ens_contig
);
use Bio::SearchIO;
my $searchio = new Bio::SearchIO(
-format => 'blast',
-file => 't/converter.blast'
);
my @hsps = ();
while(my $result = $searchio->next_result){
while(my $hit = $result->next_hit){
while(my $hsp = $hit->next_hsp){
push @hsps, $hsp;
}
}
}
my @align_features = @{$converter->convert(\@hsps)};
my $align_feature = shift @align_features;
ok $align_feature->isa('Bio::EnsEMBL::BaseAlignFeature');
ok $align_feature->start, 1;
ok $align_feature->end, 504;
ok $align_feature->strand, 0;
ok $align_feature->hstart, 21;
ok $align_feature->hend, 1532;
ok $align_feature->cigar_string, '504M';
$align_feature = shift @align_features;
ok $align_feature->isa('Bio::EnsEMBL::BaseAlignFeature');
ok $align_feature->start, 64;
ok $align_feature->end, 324;
ok $align_feature->strand, 0;
ok $align_feature->hstart, 182;
ok $align_feature->hend, 844;
ok $align_feature->cigar_string, '29M18I22M11I20MD33M4I22M3I25M5I21MI33MD14M';
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment