hw6/genbank_CD_draw.pl

#!/usr/local/bin/perl -w

# Extract sequence features from a GenBank report with BioPerl.
# Once the features are parsed, they can be printed in any format.
# WI Bioinformatics course - Feb 2002 - Answers to homework6.2

use Bio::SeqIO; 
use GD;         
use integer;    # Do integer division - helpful since pixels are the quanta of graphics

########################  User defined variables  #########################

$imagefile = "bmp7.png";                          # output file

# $div shows the number of nt or aa per pixel - requires a larger integer for longer sequences
$div = 3;
$image_width = 800; 
$image_height = 600;
$seq_width = 20;
$margin_left = 20;
$margin_top = 50;

#########################  Set up the picture and colors  ##########################

# Create a new image
$img = new GD::Image($image_width, $image_height);

# Describe colors you want to use: each integer is 0-255 for Red, Green, and Blue
$white = $img->colorAllocate(255,255,255);
$black = $img->colorAllocate(0,0,0);
$green = $img->colorAllocate(0, 255, 0);
$blue = $img->colorAllocate(0, 0, 255);
$red = $img->colorAllocate(255, 0, 0); 
  
########################## Main script #########################
 
# Formats: Fasta, EMBL. GenBank, Swissprot (swiss), PIR and GCG
$seqin = Bio::SeqIO->new( '-format' => 'Genbank' , -file => 'bmp7.txt'); # create an object for input file

while(my $seqobj = $seqin->next_seq()) {
    foreach $feat ($seqobj->all_SeqFeatures()) {
   if ($feat->primary_tag eq "gene") {       # gene feature
       $geneStart = $feat->start;            # start position of the gene
       $geneEnd = $feat->end;                # end position of the gene
       foreach $tag($feat->all_tags()) {
      if ($tag eq "gene") {
          $geneName = join(' ',$feat->each_tag_value($tag)); # gene name:BMP7

# Draw the name of the gene in black:          
          $img->string(gdGiantFont, $margin_left, 5, $geneName, $black);

# Draw a bar in red for the sequence of the gene:
          my $x1 = $margin_left + ($geneStart / $div);
          my $x2 = $margin_left + ($geneEnd / $div);
          my $y1 = $margin_top;
          my $y2 = $y1 + $seq_width;
          $img->filledRectangle($x1, $y1, $x2, $y2, $red);

# Draw nt(red) along length of sequence, with vertical lines(blue):
          for ($i = 1; $i < ($geneEnd - $geneStart); $i+= 200) {
         my $x1 = $margin_left + ($geneStart / $div) + (($i - 1) / $div);
         $img->string(gdGiantFont, $x1, 20, $i, $red); #draw nt in red along length of the sequence
         $img->line($x1, 50, $x1, 30, $blue); #
          }
          
      }
       }
   }
    if ($feat->primary_tag eq "CDS") {       # CDS feature
       $CDstart = $feat->start;             # start positionof the CD
       $CDend = $feat->end;                 # end position of the CD

#Draw a bar in green for the sequence of the CD:
       my $x1 = $margin_left + ($CDstart / $div);
       my $x2 = $margin_left + ($CDend / $div);
       my $y1 = $margin_top + $seq_width + 30;
       $y2 = $y1 + 20;
       $img->filledRectangle($x1, $y1, $x2, $y2, $green);
   }

   if ($feat->primary_tag eq "variation") { # variation feature
       $variationStart = $feat->start;      # start position
       $variationEnd = $feat->end;          # end position

#Draw a line in black for the location of the SNPs:
       my $x1 = $margin_left + ($variationStart / $div);
       my $x2 = $margin_left + ($variationEnd / $div);
       my $y1 = $margin_top;
       my $y2 = $y1 + $seq_width;
       $img->line($x1, $y1, $x2, $y2, $black);
   }
   
    }
}

open(OUT, ">$imagefile") || die "Cannot write to $imagefile: $!\n";
    
# Print this image to a file as a PNG graphic
   print OUT $img->png;
   close OUT;