scripts/seqstats/bp_chaos_plot.PLS

   1 #!perl
   2
   3 use strict;
   4 use warnings;
   5
   6 use Bio::SeqIO;
   7 use Getopt::Long;
   8 use GD;
   9
  10 use vars qw( $USAGE %VALIDFORMATS);
  11
  12 %VALIDFORMATS = ( 'png'  => 1,
  13                   'jpeg' => 1,
  14                   'gd2'  => 1,
  15                   'gd'   => 1,
  16                   'gif'  => 1,
  17                   'wbmp' => 1 );
  18
  19 $USAGE = "usage:\tchaos_plot -i/--input=INPUTFILE -f/--format=SEQFORMAT \n".
  20     "\t-o/--output=OUTPUTFILE -g/--graphics=GRAPHIC TYPE\n".
  21     "\t-w/--width=600 -h/--height=400\n";
  22
  23 $USAGE .= "\tValid graphics formats: (" . join(",", ( keys %VALIDFORMATS )) .")\n";
  24 $USAGE .= "\tImage size defaults to 600x400, SEQFORMAT to fasta\n";
  25 $USAGE .= "\tINPUTFILE can also be read from STDIN\n";
  26
  27 my ($format,$graph,$width,$height,$seqfile,$output) = ('fasta', 'png', 600, 400);
  28 GetOptions( "i|input:s"           => \$seqfile,
  29             "f|format:s"          => \$format,
  30             "o|output:s"          => \$output,
  31             "g|graph|graphics:s"  => \$graph,
  32             "width:i"             => \$width,
  33             "height:i"            => \$height
  34             );
  35
  36 if( ! $output || ! $VALIDFORMATS{$graph} ) {
  37     die $USAGE ;
  38 }
  39 my $seqin;
  40 $seqfile = shift unless $seqfile;
  41 if( defined $seqfile ) {
  42     print "Could not open file [$seqfile]\n$USAGE" and exit unless -e $seqfile;
  43     $seqin = new Bio::SeqIO(-format => $format,
  44                             -file   => $seqfile);
  45 } else {
  46     $seqin = new Bio::SeqIO(-format => $format,
  47                             -fh     => \*STDIN);
  48 }
  49
  50 my $img = new GD::Image($width,$height);
  51 my $white = $img->colorAllocate(255,255,255);
  52 my $black = $img->colorAllocate(0,0,0);
  53
  54 my $seq = $seqin->next_seq;
  55 die("Sequence type must be DNA not " . $seq->alphabet())
  56     unless $seq->alphabet ne 'dna' or $seq->alphabet ne 'rna';
  57 my %nmerdata;
  58 my $len = $seq->length();
  59 my $max = 0;
  60
  61 my ($x,$y) = ( 0.5, 0.5);
  62 $img->string(gdGiantFont, 1,1, 'A', $black);
  63 $img->string(gdGiantFont, 0,$height - 15, 'C', $black);
  64 $img->string(gdGiantFont, $width - 15,1, 'T', $black);
  65 $img->string(gdGiantFont, $width - 15,$height -20, 'G', $black);
  66
  67 for( my $i = 1; $i <= $len; $i++ ) {
  68
  69     my $base = lc $seq->subseq($i,$i);
  70     if( $base eq 'a' ) {
  71         $x *= 0.5;
  72         $y  *= 0.5;
  73     } elsif ( $base eq 'g' ) {
  74         $x = ( $x + 1.0 ) * 0.5;
  75         $y  = ( $y + 1.0  ) * 0.5;
  76     } elsif ( $base eq 'c' ) {
  77         $x *= 0.5;
  78         $y  = ( $y + 1.0  ) * 0.5;
  79     } elsif ( $base eq 't' or $base eq 'u' ) {
  80         $x = ( $x + 1.0 ) * 0.5;
  81         $y  *= 0.5;
  82     }
  83
  84     $img->setPixel($x * $width,$y * $height, $black);
  85 }
  86 open(OUT, ">$output");
  87 binmode OUT;
  88 $graph =~ s/jpg/jpeg/;
  89
  90 print OUT $img->$graph();
  91
  92
  93 __END__
  94
  95 =head1 NAME
  96
  97 bp_chaos_plot - a chaos plot from DNA and RNA sequences
  98
  99 =head1 SYNOPSIS
 100
 101   bp_chaos_plot.pl -i/--input=INPUTFILE -f/--format=SEQFORMAT
 102         -o/--output=OUTPUTFILE -g/--graphics=GRAPHIC FORMAT
 103         -w/--width=WIGHT -h/--height=HEIGHT
 104
 105 =head1 DESCRIPTION
 106
 107 This scripts generates image files using GD image library to visualize
 108 nucleotide sequences using chaos plot.
 109
 110 =head1 OPTIONS
 111
 112 Valid graphics formats are currently gd, gd2, png, wbmp, jpeg and gif.
 113
 114 The default size of the image file is 600x400.
 115
 116 The sequence input can be provided using any of the three methods:
 117
 118 =over 3
 119
 120 =item unnamed argument
 121
 122   bp_chaos_plot filename
 123
 124 =item named argument
 125
 126   bp_chaos_plot -i filename
 127
 128 =item standard input
 129
 130   bp_chaos_plot < filename
 131
 132 =back
 133
 134 =head1 FEEDBACK
 135
 136 =head2 Mailing Lists
 137
 138 User feedback is an integral part of the evolution of this and other
 139 Bioperl modules. Send your comments and suggestions preferably to
 140 the Bioperl mailing list.  Your participation is much appreciated.
 141
 142   bioperl-l@bioperl.org                  - General discussion
 143   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
 144
 145 =head2 Reporting Bugs
 146
 147 Report bugs to the Bioperl bug tracking system to help us keep track
 148 of the bugs and their resolution. Bug reports can be submitted via the
 149 web:
 150
 151   https://redmine.open-bio.org/projects/bioperl/
 152
 153 =head1 AUTHOR - Jason Stajich
 154
 155 Email jason@bioperl.org
 156
 157 =head1 HISTORY
 158
 159 This code is based on EMBOSS C code for chaos.c by Ian Longden.
 160 Included are documentation from EMBOSS code:
 161
 162 Chaos produces a chaos plot.  The original application is part of the
 163 ACEDB genome database package, written by ** Richard Durbin (MRC LMB,
 164 UK) rd@mrc-lmba.cam.ac.uk, and Jean Thierry-Mieg (CRBM du CNRS,
 165 France) mieg@crbm1.cnusc.fr
 166
 167 =cut
 168