#!/usr/local/bin/perl -w $| =1; use strict; use sigtrap; use GD; use BPlite; use Getopt::Std; use vars qw($opt_x $opt_y); getopts('x:y:'); my $usage = " plotBlast - a simple BLAST report plotter using GD. usage: $0 [options] <blast report> <query> options: -x <int> width of image (scales) -y <int> height of image (does not scale) "; die $usage unless @ARGV == 2; my ($blast_file, $fasta_file) = @ARGV; my ($qdef, $qseq) = getFasta($fasta_file); # image setup my $BEGIN = 0; my $END = length($qseq) -1; my $LENGTH = length($qseq); my $WIDTH = $opt_x ? $opt_x : 1000; my $HEIGHT = $opt_y ? $opt_y : 500; my $STEP = 10; my $TEXTAREA = 100; my $MARGIN = $WIDTH - 10; my $bpp = ($MARGIN-$TEXTAREA)/($END-$BEGIN); my $B = $BEGIN * $bpp - $TEXTAREA; my $img = new GD::Image($WIDTH, $HEIGHT); # color allocation my $black = $img->colorAllocate(0,0,0); my $white = $img->colorAllocate(255,255,255); my $gray = $img->colorAllocate(128,128,128); my $red = $img->colorAllocate(255,0,0); my $yellow = $img->colorAllocate(255,255,0); my $green = $img->colorAllocate(0,255,0); my $cyan = $img->colorAllocate(0,255,255); my $blue = $img->colorAllocate(0,0,255); # header $img->string(gdSmallFont, 10,10, "plotBlast", $white); # tick marks and numbers my $MARK = $LENGTH / 10000; if ($MARK < 1) {$MARK = 1} else {$MARK = int($MARK)} $MARK *= 1000; $img->line($TEXTAREA, 10, $MARGIN, 10, $white); for(my $i=$BEGIN;$i<=$END;$i++) { if ($i % $MARK == 0 and $i != $BEGIN) { $img->line($i*$bpp - $B, 10, $i*$bpp - $B, 15, $white); $img->string(gdSmallFont, $i*$bpp-2-$B, 17, $i, $white); } } # plot blast report my $Y = 40; # spacing from top open(BLAST, $blast_file) or die; my $blast = new BPlite(\*BLAST); while (my $sbjct = $blast->nextSbjct) { my ($name) = $sbjct->name =~ /^>(\S+)/; $img->string(gdSmallFont, 10, $Y, $name, $white); $Y += $STEP; while (my $hsp = $sbjct->nextHSP) { # color ramp my $color; if ($hsp->bits > 200) {$color = $white} elsif ($hsp->bits > 150) {$color = $blue} elsif ($hsp->bits > 100) {$color = $cyan} elsif ($hsp->bits > 80) {$color = $green} elsif ($hsp->bits > 60) {$color = $yellow} elsif ($hsp->bits > 40) {$color = $red} else {$color = $gray} my $x1 = $hsp->qb * $bpp; my $x2 = $hsp->qe * $bpp; my $y1 = $Y; my $y2 = $Y + 5; my $method = 'filledRectangle'; if ($x1 > $x2) { ($x1, $x2) = ($x2, $x1); # swqp coordiates $method = 'rectangle'; # open rectangle on negative strand } $img->$method($x1 + $TEXTAREA, $y1, $x2 + $TEXTAREA, $y2, $color); } } print $img->png; sub getFasta { my ($file) = @_; open(FILE, $file) or die; my ($def, @seq) = <FILE>; chomp @seq; my $seq = join("", @seq); close FILE; return ($def, $seq); }