[Bioperl-l] how to create Graphics object from home-made SQL table with parsed blast data?

Sean Davis sdavis2 at mail.nih.gov
Fri Oct 29 06:18:59 EDT 2004


Marcus,

The following was NOT!!! meant for public consumption, but it does show  
a script that I use for generating graphics from BLAT output (it is in  
a CGI content, so sorry about that part of things).  You could likely  
munge it to something useful to you.  In particular, look at how I set  
up the Panel and how I generate seqfeatures on the fly from DB query  
results.  Perl Gurus can use this as fodder for the critique mill.

Sean


sub return_graphical_results {
     my $self = shift;

     # Get Query Object
     my $q        = $self->query();
     my $template = $self->load_tmpl('mode4.html');
     my $oligoid  = $q->param('name');
     # Get database handles
     my $probedbh = $self->param('probedbh');
     my $ensdbh   = $self->param('ensdbh');

     $template->param(OLIGOID => $oligoid);

     my  
($fh,$filename)=tempfile(DIR=>'/Library/WebServer/Documents/ 
oligodir',SUFFIX=>'.png');

     my $sth=$probedbh->prepare("select length from oligo where  
oligo='$oligoid'");
     eval {
	$sth->execute();
     };
     if ($@) {
	carp "Database error: $DBI::errstr\n";
	$probedbh->rollback;
     }
     my ($Qsize)=$sth->fetchrow_array;

     my $panel=Bio::Graphics::Panel->new(-key_style=>'between',
					-length => $Qsize,
					-width  => 800,
					-pad_left => 10,
					-pad_right => 10,
				       );
     my $full_length =  
Bio::SeqFeature::Generic->new(-start=>1,-end=>$Qsize);
     my $track=$panel->add_track($full_length,
				-glyph   => 'arrow',
				-tick    => 2,
				-fgcolor => 'black',
				-double  => 1,
			       );
     my %track;
     my @db=qw/ ensGene knownGene H-inv RefSeq genomic dbEST /;
     foreach my $dbname (@db) {
	$track{$dbname}=$panel->add_track(-key    => $dbname,
					  -glyph   => 'graded_segments',
					  -label   => 1,
					  -min_score => 100,
					  -max_score => 140,
					  -bgcolor   => 'blue',
					  -font2color=> 'red',
					  -description => sub {
					      my $feat=shift;
					      return("score=" . $feat->score);
					  },
					 );
     }

     my $sql=q{select db,matches,mismatch,strand,Qstart,Qend,
	      target_id,Tstart,Tend,block_sizes,block_count,Qstarts
	      from hit,oligo,db,analysis
	      where oligo.oligo=?
	      and oligo.oligo_id=hit.oligo_id
	      AND hit.analysis_id=analysis.analysis_id
	      AND db.db_id=analysis.db_id
	     };

     $sth=$probedbh->prepare($sql);
     eval {
	$sth->execute($oligoid);
     };
     if ($@) {
	carp "Database error: $DBI::errstr\n";
	$probedbh->rollback;
     }


# LOOP OVER THE HITS AND GENERATE GRAPHICS

     while (my  
($db,$match,$mismatch,$strand,$Qstart,$Qend,$target_id,$Tstart,
	       $Tend,$block_sizes,$block_count,$Qstarts) =  
$sth->fetchrow_array) {
	my $score = 2*$match-$mismatch;
	my $feature = Bio::SeqFeature::Generic->new(-display_name =>  
$target_id,
						    -start        => $Qstart+1,
						    -end          => $Qend,
						    -tstart       => $Tstart,
						    -tend         => $Tend,
						    -score        => $score,
						   );
	my @Qstarts=split(',',$Qstarts);
	my @block_sizes=split(',',$block_sizes);
	for (my $i=0;$i<=$#Qstarts;$i++) {
	    if ($strand eq '+') {
		my $subfeature = Bio::SeqFeature::Generic->new(
							       -score        => $score,
							       -start        => $Qstarts[$i]+1,
							       -end          => $Qstarts[$i]+$block_sizes[$i],
							      );
		$feature->add_sub_SeqFeature($subfeature);
	    } else {		# negative strand!
		my $subfeature = Bio::SeqFeature::Generic->new(
							       -score        => $score,
							       -end          => $Qsize-$Qstarts[$i],
							       -start        => $Qsize-$block_sizes[$i]-$Qstarts[$i]+1
							      );
		$feature->add_sub_SeqFeature($subfeature);
	    }
	}
	$track{$db}->add_feature($feature);

     }
     print $fh $panel->png;
     $panel->finished;
     $fh->close;
     $filename=~s/\/Library\/WebServer\/Documents//;
     $template->param(FILENAME => $filename);
     return $template->output;
}



More information about the Bioperl-l mailing list