[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