[Bioperl-l] How can I get add_sub_SeqFeature to work as I want?
Marcus Claesson
m.claesson at student.ucc.ie
Wed Nov 24 12:43:51 EST 2004
Hi!
I've had a problem some time with getting the add_sub_SeqFeature to work
for my own parsed blast data. Last time I posted it to the list nobody
was be able to help me, but I have now simplified the problem a bit.
I have blast parsed data like this in a text file (the columns are
sbj_count,hsp_count,score,query_begin,query_end,strand),(I left out
sbj_name for simplicity):
1 1 445 1148 375 -1
1 2 341 1717 1151 -1
2 1 364 1148 378 -1
2 2 344 1690 1151 -1
3 1 283 1145 381 -1
3 2 233 1714 1151 -1
4 1 182 1154 375 -1
4 2 124 1702 1160 -1
The only way I know the hit has more than one hsp is if the next line
has the same sbj_count. Based on this I'd like to construct a Graphics
Panel as an graphical overview of the blast hits. Below is the script
that shows every hsp as separate hits (I want hsps belonging to the same
hit on one line):
#!/usr/bin/perl
use Bio::Graphics;
use Bio::SeqFeature::Generic;
my $panel = Bio::Graphics::Panel->new(-length => 2000,
-width => 500,
-pad_left => 10,
-pad_right => 10);
my $track = $panel->add_track(-glyph => 'graded_segments',
-label => 1,
-strand_arrow => 1,
-connector => 'dashed',
-bgcolor => 'blue',
-font2color => 'red',
-sort_order => 'score',
-description => sub {
my $feature = shift;
my ($score) = $feature->score;
"Score=$score"});
$old_sbj_count = 0;
while (<>) {
($sbj_count,$hsp_count,$score,$qbegin,$qend,$strand) = split;
my $feature = Bio::SeqFeature::Generic->new(-start => $qbegin,
-end => $qend,
-strand => $strand,
-score => $score);
if ($old_sbj_count eq $sbj_count) {
my $subfeature = Bio::SeqFeature::Generic->new(
-start => $qbegin,
-end => $qend,
-strand=> $strand,
-score => $score);
$feature->add_sub_SeqFeature($subfeature,"EXPAND");
}
$old_sbj_count = $sbj_count;
$track->add_feature($feature);
}
print $panel->png;
I would be extremely grateful for any help on this since I've been
struggling for a long time...
Best regards,
Marcus
More information about the Bioperl-l
mailing list