[Bioperl-l] How can I get add_sub_SeqFeature to work as I want?

[Jason Stajich]

The problem is that subfeatures are never added to the same top-level 
feature.

Try this (note that the 'score' and 'description' [if you were to 
assign description] will be whatever the 1st instance of the subject 
is.  See also Lincoln's example here
http://bioperl.org/HOWTOs/Graphics-HOWTO/parsing-blast.html if you 
haven't seen this already.

We're using a hash (%sbj) to figure out how to group things by a key 
($sbj_count).

my %sbj;
while(<>) {
   my  ($sbj_count,$hsp_count,$score,$qbegin,$qend,$strand) = split;
     my $feature = Bio::SeqFeature::Generic->new(-start   => $qbegin,
						-end     => $qend,
						-strand  => $strand,
						-score   => $score);

   if( defined $sbj{$sbj_count} ) {
    # add this feature as a subfeature instead since we've already seen 
one for this subject
     $sbj{$sbj_count}->add_sub_SeqFeature($feature,"EXPAND");
  } else {
    # first time we've seen this subject so we'll just store this feature
     $sbj{$sbj_count} = $feature;
  }
}

# now lets add all those features to the track
for my $f ( values %sbj ) {
  $track->add_feature($f);
}

On Nov 24, 2004, at 1:53 PM, Allen Day wrote:

> someone correct me if i'm wrong, but i think you need to add all the
> features for a track at once.  maybe you should load this file into a
> database (even SQLite) so you don't have to keep track of the
> previous_row/next_row-ness of the text file.
>
> -allen
>
> On Wed, 24 Nov 2004, Marcus Claesson wrote:
>
>> 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
>>
>> _______________________________________________
>> Bioperl-l mailing list
>> Bioperl-l at portal.open-bio.org
>> http://portal.open-bio.org/mailman/listinfo/bioperl-l
>>
> _______________________________________________
> Bioperl-l mailing list
> Bioperl-l at portal.open-bio.org
> http://portal.open-bio.org/mailman/listinfo/bioperl-l
>
>
--
Jason Stajich
jason.stajich at duke.edu
http://www.duke.edu/~jes12/