code_perl/Maasha/Plot.pm

   1 package Maasha::Plot;
   2
   3 # Copyright (C) 2007 Martin A. Hansen.
   4
   5 # This program is free software; you can redistribute it and/or
   6 # modify it under the terms of the GNU General Public License
   7 # as published by the Free Software Foundation; either version 2
   8 # of the License, or (at your option) any later version.
   9
  10 # This program is distributed in the hope that it will be useful,
  11 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13 # GNU General Public License for more details.
  14
  15 # You should have received a copy of the GNU General Public License
  16 # along with this program; if not, write to the Free Software
  17 # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18
  19 # http://www.gnu.org/copyleft/gpl.html
  20
  21
  22 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> DESCRIPTION <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  23
  24
  25 # Routines to plot stuff with Gnuplot and SVG.
  26
  27
  28 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  29
  30
  31 use strict;
  32 use Data::Dumper;
  33 use SVG;
  34 use IPC::Open2;
  35 use Maasha::Common;
  36 use Maasha::Filesys;
  37 use Maasha::Seq;
  38 use vars qw ( @ISA @EXPORT );
  39
  40 use constant {
  41     WIDTH  => 800,
  42     HEIGHT => 600,
  43     MARGIN => 40,
  44 };
  45
  46 @ISA = qw( Exporter );
  47
  48
  49 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  50
  51
  52 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> LINEPLOTS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  53
  54
  55 sub lineplot_simple
  56 {
  57     # Martin A. Hansen, January 2008.
  58
  59     # Plots a simple lineplot using Gnuplot.
  60
  61     my ( $data,       # data table - each column will be plottet as one line.
  62          $options,    # options hash
  63          $tmp_dir,    # temporary directory
  64        ) = @_;
  65
  66     # Returns list.
  67
  68     my ( $tmp_file, $pid, $fh_in, $fh_out, $cmd, $i, $line, @lines, $xtic_space, @plot_cmd );
  69
  70     $tmp_dir ||= $ENV{ 'BP_TMP' };
  71
  72     $tmp_file = "$tmp_dir/lineplot_simple.tab";
  73
  74     $fh_out = Maasha::Filesys::file_write_open( $tmp_file );
  75
  76     map { print $fh_out join( "\t", @{ $_ } ), "\n" } @{ $data };
  77
  78     close $fh_out;
  79
  80     $options->{ "terminal" } ||= "dumb";
  81
  82     $cmd  = "gnuplot";
  83
  84     $pid = open2( $fh_out, $fh_in, $cmd );
  85
  86     # $fh_in = \*STDERR;
  87
  88     print $fh_in "set terminal $options->{ 'terminal' }\n";
  89     print $fh_in "set title \"$options->{ 'title' }\"\n"   if $options->{ "title" };
  90     print $fh_in "set xlabel \"$options->{ 'xlabel' }\"\n" if $options->{ "xlabel" };
  91     print $fh_in "set ylabel \"$options->{ 'ylabel' }\"\n" if $options->{ "ylabel" };
  92     print $fh_in "set grid\n"                              if not $options->{ "terminal" } eq "dumb";
  93     print $fh_in "set autoscale\n";
  94     print $fh_in "unset key\n";
  95     print $fh_in "set xtics border in scale 0 nomirror rotate by 90 offset character 0, 0, 0\n";
  96
  97     for ( $i = 1; $i < scalar @{ $data->[ 0 ] } + 1; $i++ ) {
  98         push @plot_cmd, qq("$tmp_file" using $i with lines ls 1);
  99     }
 100
 101     print $fh_in "plot " . join( ", ", @plot_cmd ) . "\n";
 102
 103     close $fh_in;
 104
 105     while ( $line = <$fh_out> )
 106     {
 107         chomp $line;
 108
 109         push @lines, $line;
 110     }
 111
 112     close $fh_out;
 113
 114     waitpid $pid, 0;
 115
 116     unlink $tmp_file;
 117
 118     return wantarray ? @lines : \@lines;
 119 }
 120
 121
 122 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> HISTOGRAMS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 123
 124
 125 sub histogram_simple
 126 {
 127     # Martin A. Hansen, August 2007.
 128
 129     # Plots a simple histogram using Gnuplot.
 130
 131     my ( $data,       # list of [ xlabel, data value ] tuples
 132          $options,    # options hash
 133        ) = @_;
 134
 135     # Returns list.
 136
 137     my ( $pid, $fh_in, $fh_out, $cmd, $i, $line, @lines );
 138
 139     $options->{ "terminal" } ||= "dumb";
 140
 141     $cmd  = "gnuplot";
 142
 143     $pid = open2( $fh_out, $fh_in, $cmd );
 144
 145 #    $fh_in = \*STDERR;
 146
 147     # print $fh_in "set terminal $options->{ 'terminal' } 10 \n"; # adsjust fontsize to 10 - find some other way to do this, because it don't work with SVG.
 148     print $fh_in "set terminal $options->{ 'terminal' }\n";
 149     print $fh_in "set title \"$options->{ 'title' }\"\n"   if $options->{ "title" };
 150     print $fh_in "set xlabel \"$options->{ 'xlabel' }\"\n" if $options->{ "xlabel" };
 151     print $fh_in "set ylabel \"$options->{ 'ylabel' }\"\n" if $options->{ "ylabel" };
 152     print $fh_in "set autoscale\n";
 153     print $fh_in "unset key\n";
 154     print $fh_in "set style fill solid\n";
 155     print $fh_in "set style histogram title offset character 0, 0, 0\n";
 156     print $fh_in "set style data histograms\n";
 157     print $fh_in "set xtics border in scale 0 nomirror rotate by 90 offset character 0, 0, 0\n";
 158
 159     print $fh_in "plot '-' using 2:xticlabels(1)\n";
 160
 161     for ( $i = 0; $i < @{ $data }; $i++ )
 162     {
 163         print $fh_in join( "\t", "\"$data->[ $i ]->[ 0 ]\"", $data->[ $i ]->[ 1 ] ), "\n";
 164     }
 165
 166     close $fh_in;
 167
 168     while ( $line = <$fh_out> )
 169     {
 170         chomp $line;
 171
 172         push @lines, $line;
 173     }
 174
 175     close $fh_out;
 176
 177     waitpid $pid, 0;
 178
 179     return wantarray ? @lines : \@lines;
 180 }
 181
 182
 183 sub histogram_lendist
 184 {
 185     # Martin A. Hansen, August 2007.
 186
 187     # Plots a histogram using Gnuplot.
 188
 189     my ( $data,       # list of [ xlabel, data value ] tuples
 190          $options,    # options hash
 191        ) = @_;
 192
 193     # Returns list.
 194
 195     my ( $pid, $fh_in, $fh_out, $cmd, $i, $line, @lines, $xtic_space );
 196
 197     $options->{ "terminal" } ||= "dumb";
 198
 199     if ( $data->[ -1 ]->[ 0 ] <= 10 ) {
 200         $xtic_space = 1;
 201     } elsif ( $data->[ -1 ]->[ 0 ] <= 100 ) {
 202         $xtic_space = 5;
 203     } elsif ( $data->[ -1 ]->[ 0 ] <= 250 ) {
 204         $xtic_space = 10;
 205     } elsif ( $data->[ -1 ]->[ 0 ] <= 500 ) {
 206         $xtic_space = 20;
 207     } elsif ( $data->[ -1 ]->[ 0 ] <= 1000 ) {
 208         $xtic_space = 50;
 209     } elsif ( $data->[ -1 ]->[ 0 ] <= 2500 ) {
 210         $xtic_space = 100;
 211     } elsif ( $data->[ -1 ]->[ 0 ] <= 5000 ) {
 212         $xtic_space = 250;
 213     } elsif ( $data->[ -1 ]->[ 0 ] <= 10000 ) {
 214         $xtic_space = 500;
 215     } elsif ( $data->[ -1 ]->[ 0 ] <= 50000 ) {
 216         $xtic_space = 1000;
 217     } elsif ( $data->[ -1 ]->[ 0 ] <= 100000 ) {
 218         $xtic_space = 5000;
 219     }
 220
 221     $cmd  = "gnuplot";
 222
 223     $pid = open2( $fh_out, $fh_in, $cmd );
 224
 225     print $fh_in "set terminal $options->{ 'terminal' }\n";
 226     print $fh_in "set title \"$options->{ 'title' }\"\n"   if $options->{ "title" };
 227     print $fh_in "set xlabel \"$options->{ 'xlabel' }\"\n" if $options->{ "xlabel" };
 228     print $fh_in "set ylabel \"$options->{ 'ylabel' }\"\n" if $options->{ "ylabel" };
 229     print $fh_in "set autoscale\n";
 230     print $fh_in "unset key\n";
 231     print $fh_in "set style fill solid\n";
 232     print $fh_in "set style histogram clustered gap 1 title offset character 0, 0, 0\n";
 233     print $fh_in "set style data histograms\n";
 234     print $fh_in "set xtics 0,$xtic_space border out nomirror\n";
 235
 236     print $fh_in "plot '-' using 1\n";
 237
 238     for ( $i = 0; $i < @{ $data }; $i++ )
 239     {
 240         $data->[ $i ]->[ 0 ] = "." if $data->[ $i ]->[ 0 ] % 10 != 0;
 241
 242         print $fh_in join( "\t", $data->[ $i ]->[ 1 ] ), "\n";
 243     }
 244
 245     close $fh_in;
 246
 247     while ( $line = <$fh_out> )
 248     {
 249         chomp $line;
 250
 251         push @lines, $line;
 252     }
 253
 254     close $fh_out;
 255
 256     waitpid $pid, 0;
 257
 258     return wantarray ? @lines : \@lines;
 259 }
 260
 261
 262 sub histogram_chrdist
 263 {
 264     # Martin A. Hansen, August 2007.
 265
 266     # Plots a histogram using Gnuplot.
 267
 268     my ( $data,       # list of [ xlabel, data value ] tuples
 269          $options,    # options hash
 270        ) = @_;
 271
 272     # Returns list.
 273
 274     my ( $pid, $fh_in, $fh_out, $cmd, $i, $line, @lines );
 275
 276     $options->{ "terminal" } ||= "dumb";
 277
 278     $cmd  = "gnuplot";
 279
 280     $pid = open2( $fh_out, $fh_in, $cmd );
 281
 282     print $fh_in "set terminal $options->{ 'terminal' }\n";
 283     print $fh_in "set title \"$options->{ 'title' }\"\n"   if $options->{ "title" };
 284     print $fh_in "set xlabel \"$options->{ 'xlabel' }\"\n" if $options->{ "xlabel" };
 285     print $fh_in "set ylabel \"$options->{ 'ylabel' }\"\n" if $options->{ "ylabel" };
 286     print $fh_in "set autoscale\n";
 287     print $fh_in "unset key\n";
 288     print $fh_in "set style fill solid\n";
 289     print $fh_in "set style histogram title offset character 0, 0, 0\n";
 290     print $fh_in "set style data histograms\n";
 291     print $fh_in "set xtics border in scale 0 nomirror rotate by 90 offset character 0, 0, 0\n";
 292
 293     print $fh_in "plot '-' using 2:xticlabels(1)\n";
 294
 295     for ( $i = 0; $i < @{ $data }; $i++ ) {
 296         print $fh_in join( "\t", "\"$data->[ $i ]->[ 0 ]\"", $data->[ $i ]->[ 1 ] ), "\n";
 297     }
 298
 299     close $fh_in;
 300
 301     while ( $line = <$fh_out> )
 302     {
 303         chomp $line;
 304
 305         push @lines, $line;
 306     }
 307
 308     close $fh_out;
 309
 310     waitpid $pid, 0;
 311
 312     return wantarray ? @lines : \@lines;
 313 }
 314
 315
 316 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> KARYOGRAM <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 317
 318
 319 sub karyogram
 320 {
 321     # Martin A. Hansen, August 2007.
 322
 323     # Plot hits on a karyogram for a given genome.
 324
 325     my ( $data,      # list of [ chr, beg, end ] triples
 326          $options,   # hashref with options
 327        ) = @_;
 328
 329     # Returns string
 330
 331     my ( $karyo_file, $svg, $features, $karyo );
 332
 333     if ( $options->{ "genome" } eq "hg18" )
 334     {
 335         $karyo_file = $ENV{ 'BP_DIR' } . "/bp_data/human_cytobands.txt";
 336     }
 337     elsif( $options->{ "genome" } eq "mm9" )
 338     {
 339         $karyo_file = $ENV{ 'BP_DIR' } . "/bp_data/mouse_cytobands.txt";
 340     }
 341
 342     $karyo = parse_karyo_data( $karyo_file );
 343
 344     $svg = init_svg();
 345
 346     chromosome_layout( $svg, $karyo, $data );
 347
 348     return $svg->xmlify;
 349 }
 350
 351
 352 sub parse_karyo_data
 353 {
 354     # X       q26.1   129700001       130200000       gneg
 355
 356     # color: /etc/X11/rgb.txt
 357
 358     my ( $file,
 359        ) = @_;
 360
 361     my ( $fh, $chr, $line, $name, $beg, $end, $color, %features, %color_hash );
 362
 363     %color_hash = (
 364         acen    => "DarkGrey",
 365         gneg    => "white",
 366         gpos100 => "black",
 367         gpos75  => "DarkGrey",
 368         gpos66  => "DarkGrey",
 369         gpos50  => "grey",
 370         gpos33  => "LightGrey",
 371         gpos25  => "LightGrey",
 372         gvar    => "LightGrey",
 373         stalk    => "DarkGrey",
 374 #        gpos75  => "rgb(169,169,169)",
 375 #        gpos66  => "gray66",
 376 #        gpos66  => "#8e8e8e",
 377 #        gpos50  => "gray50",
 378 #        gpos33  => "#e3e3e3",
 379 #        gpos33  => "gray33",
 380 #        gpos25  => "gray25",
 381 #        stalk   => "rgb(169,169,169)",
 382 #        stalk    => "gray66",
 383     );
 384
 385     $fh = Maasha::Filesys::file_read_open( $file );
 386
 387     while ( $line = <$fh> )
 388     {
 389         chomp $line;
 390
 391         next if $line =~ /^#/;
 392
 393         # ( $chr, $name, $beg, $end, $color ) = split "\t", $line;
 394         ( $chr, $beg, $end, $name, $color ) = split "\t", $line;
 395
 396 #        if ( $color =~ /^gpos(\d+)/ ) {
 397 #            $color = color_intensity( $1 );
 398 #        } elsif ( exists $color_hash{ $color } ) {
 399             $color = $color_hash{ $color };
 400 #        } else {
 401 #            die qq(ERROR: Unknown color->$color\n);
 402 #        }
 403
 404         if ( exists $features{ $chr } )
 405         {
 406             push @{ $features{ $chr } }, [ $name, $beg, $end, $color ];
 407         }
 408         else
 409         {
 410             $features{ $chr } = [ [ $name, $beg, $end, $color ] ];
 411         }
 412     }
 413
 414     close $fh;
 415
 416     return wantarray ? %features : \%features;
 417 }
 418
 419
 420 sub color_intensity
 421 {
 422     # Martin A. Hansen, September 2007.
 423
 424     # Converts a gray scale intensity in percent to rgb.
 425
 426     my ( $percent,   # color intensity
 427        ) = @_;
 428
 429     # Returns string
 430
 431     my ( $num, $hex );
 432
 433     $num = int( $percent * 256 / 100 );
 434
 435     $num--;
 436
 437     $hex = sprintf "%x", $num;
 438
 439     return "#$hex$hex$hex";
 440
 441 #    return "rgb($num,$num,$num)";
 442 }
 443
 444
 445 sub init_svg
 446 {
 447     # Martin A. Hansen, September 2005.
 448
 449     # initializes svg image.
 450
 451     # returns an image object
 452
 453     my $svg = SVG->new(
 454         width  => WIDTH,
 455         height => HEIGHT,
 456         style  => {
 457             'stroke-width' => 1,
 458             stroke => "black",
 459             font  => 'Helvetica',
 460         },
 461     );
 462
 463     return $svg;
 464 }
 465
 466
 467 sub chromosome_layout
 468 {
 469     # Martin A. Hansen, January 2004 - August 2007.
 470
 471     # Plots all chromosomes in a single
 472
 473     my ( $svg,          # image object
 474          $karyo_list,   # hashref with karyo data
 475          $feat_list,    # hashref with features
 476        ) = @_;
 477
 478     # returns an image object
 479
 480     my ( $layout_obj, $i, $x, $y, $max, $factor, $chr_len, $chr_width, $chr_cent, $chr, $feat, $karyo, @list, $A, $B );
 481
 482     $layout_obj = $svg->group(
 483         id => "chr_layout",
 484     );
 485
 486     $max       = $karyo_list->{ "chr1" }->[ -1 ]->[ 2 ];
 487     $factor    = ( HEIGHT / 2 ) / $max;
 488     $chr_width = ( HEIGHT / 4 ) / 13;
 489
 490     foreach $karyo ( keys %{ $karyo_list } ) {
 491         map { $_->[ 1 ] *= $factor; $_->[ 2 ] *= $factor } @{ $karyo_list->{ $karyo } };
 492     }
 493
 494     foreach $feat ( keys %{ $feat_list } ) {
 495         map { $_->[ 0 ] *= $factor; $_->[ 1 ] *= $factor } @{ $feat_list->{ $feat } };
 496     }
 497
 498 #    @list = sort { $A = $a; $B = $b; $A =~ s/chr//; $B =~ s/chr//; $A <=> $B } keys %{ $karyo_list };
 499
 500 #    splice @list, 0, 2;
 501 #    push @list, "chrX", "chrY";
 502
 503     $i = 0;
 504
 505     while ( $i < @list )
 506     {
 507         $chr      = $list[ $i ];
 508         $chr_len  = $karyo_list->{ $chr }->[ -1 ]->[ 2 ];
 509         $chr_cent = find_cent( $karyo_list->{ $list[ $i ] } );
 510
 511         $y = HEIGHT / 2 - $chr_len;
 512         $x = ( WIDTH / ( @list + 2 ) ) * ( $i + 1 );
 513
 514         draw_chr( $layout_obj, $x, $y, $chr_len, $chr_width, $chr_cent, $chr, $karyo_list, $feat_list );
 515
 516         $i++;
 517     }
 518 }
 519
 520
 521 sub find_cent
 522 {
 523     # Martin A. Hansen, December 2003.
 524
 525     # Finds the centromeric region in the karyo data.
 526
 527     my ( $list ) = @_;
 528
 529     my ( $acen, @nums, $cent );
 530
 531     @{ $acen } = grep { grep { /^DarkGrey$/ } @{ $_ } } @{ $list };
 532
 533     push @nums, $acen->[ 0 ]->[ 1 ];
 534     push @nums, $acen->[ 0 ]->[ 2 ];
 535     push @nums, $acen->[ 1 ]->[ 1 ];
 536     push @nums, $acen->[ 1 ]->[ 2 ];
 537
 538     @nums = sort { $a <=> $b } @nums;
 539
 540     $cent = ( $nums[ 1 ] + $nums[ 2 ] ) / 2;
 541
 542     return $cent;
 543 }
 544
 545
 546 sub draw_chr
 547 {
 548     # Martin A. Hansen, December 2003.
 549
 550     # draws a whole cromosome with or without centromeric region
 551
 552     my ( $svg,         # image object
 553          $x,           # x position
 554          $y,           # y position
 555          $chr_len,     # lenght of chromosome
 556          $chr_width,   # width of chromosome
 557          $chr_cent,    # position of centromeric region
 558          $chr,         # chromosome
 559          $karyo_list,  # hashref with karyo data
 560          $feat_list,   # hashref with features
 561        ) = @_;
 562
 563     # returns image object
 564
 565     my ( $chr_obj, $clip_obj, $gr_obj );
 566
 567     $chr_obj = $svg->group(
 568         id => $chr,
 569     );
 570
 571     if ( exists $feat_list->{ $chr } ) {
 572         draw_chr_feat( $chr_obj, $x, $y, $chr_width, $feat_list->{ $chr } );
 573     }
 574
 575     $clip_obj = $chr_obj->clipPath(
 576         id => $chr . "_clipPath",
 577     );
 578
 579     $clip_obj->rectangle(
 580         x      => sprintf( "%.3f", $x ),
 581         y      => sprintf( "%.3f", $y ),
 582         width  => sprintf( "%.3f", $chr_width ),
 583         height => sprintf( "%.3f", $chr_cent ),
 584         rx     => 10,
 585         ry     => 10,
 586     );
 587
 588     $clip_obj->rectangle(
 589         x      => sprintf( "%.3f", $x ),
 590         y      => sprintf( "%.3f", $y + $chr_cent ),
 591         width  => sprintf( "%.3f", $chr_width ),
 592         height => sprintf( "%.3f", $chr_len - $chr_cent ),
 593         rx     => 10,
 594         ry     => 10,
 595     );
 596
 597     $gr_obj = $chr_obj->group(
 598         "clip-path"  => "url(#$chr" . "_clipPath)",
 599     );
 600
 601     if ( exists $karyo_list->{ $chr } ) {
 602         draw_karyo_data( $gr_obj, $x, $y, $chr_width, $karyo_list->{ $chr } );
 603     }
 604
 605     $gr_obj->rectangle(
 606         x      => sprintf( "%.3f", $x ),
 607         y      => sprintf( "%.3f", $y ),
 608         width  => sprintf( "%.3f", $chr_width ),
 609         height => sprintf( "%.3f", $chr_cent ),
 610         fill   => 'none',
 611         rx     => 10,
 612         ry     => 10,
 613     );
 614
 615     $gr_obj->rectangle(
 616         x      => sprintf( "%.3f", $x ),
 617         y      => sprintf( "%.3f", $y + $chr_cent ),
 618         width  => sprintf( "%.3f", $chr_width ),
 619         height => sprintf( "%.3f", $chr_len - $chr_cent ),
 620         fill   => 'none',
 621         rx     => 10,
 622         ry     => 10,
 623     );
 624
 625     draw_chr_num( $chr_obj, $x, $y, $chr_len, $chr_width, $chr );
 626 }
 627
 628
 629 sub draw_chr_num
 630 {
 631     # Martin A. Hansen, December 2003.
 632
 633     # draws a cromosome number
 634
 635     my ( $svg,        # image object
 636          $x,         # x position
 637          $y,         # y position
 638          $chr_len,   # lenght of chromosome
 639          $chr_width, # width of chromosome
 640          $chr,       # chromosome number
 641        ) = @_;
 642
 643     # returns image object
 644
 645     my ( $chr_num, $chars, @a, $word_width );
 646
 647     $chr_num = $chr;
 648     $chr_num =~ s/chr//;
 649
 650     $chars = @a = split "", $chr_num;
 651
 652     $word_width = ( $chars * 8 ) / 2;
 653
 654     $svg->text(
 655         x  => sprintf("%.3f", $x + ( $chr_width / 2 ) - $word_width ),
 656         y  => sprintf("%.3f", $y + $chr_len + 15 ),
 657     )->cdata( $chr_num );
 658 }
 659
 660
 661 sub draw_karyo_data
 662 {
 663     # Martin A. Hansen, February 2004.
 664
 665     # Plots chromosome features
 666
 667     my ( $svg,
 668          $x,
 669          $y,
 670          $chr_width,
 671          $list,
 672        ) = @_;
 673
 674     # returns an image object
 675
 676     my ( $feat_beg, $feat_end, $feat_height, $i, $color, $label );
 677
 678     for ( $i = 0; $i < @{ $list }; $i++ )
 679     {
 680         ( $label, $feat_beg, $feat_end, $color ) = @{ $list->[ $i ] };
 681
 682         $feat_height = $feat_end - $feat_beg;
 683
 684         $svg->rectangle(
 685             x      => sprintf("%.3f", $x ),
 686             y      => sprintf("%.3f", $y + $feat_beg ),
 687             width  => sprintf("%.3f", $chr_width ),
 688             height => sprintf("%.3f", $feat_height ),
 689             'stroke-width' => 0,
 690             fill   => $color,
 691         );
 692     }
 693 }
 694
 695
 696 sub draw_chr_feat
 697 {
 698     # Martin A. Hansen, February 2004.
 699
 700     # Plots chromosome features
 701
 702     my ( $svg,
 703          $x,
 704          $y,
 705          $chr_width,
 706          $list,
 707        ) = @_;
 708
 709     # returns an image object
 710
 711     my ( $feat_beg, $feat_end, $feat_height, $i, $color, $height, $width, $x1, $y1, %lookup );
 712
 713     for ( $i = 0; $i < @{ $list }; $i++ )
 714     {
 715         ( $feat_beg, $feat_end, $color ) = @{ $list->[ $i ] };
 716
 717         $feat_height = $feat_end - $feat_beg;
 718
 719         $x1     = sprintf("%.0f", $x + ( $chr_width / 2 ) ),
 720         $y1     = sprintf("%.0f", $y + $feat_beg ),
 721         $width  = sprintf("%.0f", ( $chr_width / 2 ) + 5 ),
 722         $height = sprintf("%.0f", $feat_height );
 723
 724         if ( $height < 1 )
 725         {
 726             $height = 1;
 727
 728             if ( exists $lookup{ $x1 . $y1 } ) {
 729                 next;
 730             } else {
 731                 $lookup{ $x1 . $y1 } = 1;
 732             }
 733         }
 734
 735         $svg->rectangle(
 736             x      => $x1,
 737             y      => $y1,
 738             width  => $width,
 739             height => $height,
 740             stroke => $color,
 741             fill   => $color,
 742         );
 743     }
 744 }
 745
 746
 747 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> SEQUENCE LOGO <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 748
 749
 750 sub seq_logo
 751 {
 752     # Martin A. Hansen, August 2007.
 753
 754     # Calculates and renders a sequence logo in SVG format.
 755
 756     my ( $entries,   # aligned sequence entries - list of tuples
 757        ) = @_;
 758
 759     # Returns string.
 760
 761     my ( $type, $bit_max, $logo_data, $svg );
 762
 763     $type = Maasha::Seq::seq_guess_type( $entries->[ 0 ]->[ 1 ] );
 764
 765     if ( $type =~ /^p/i ) {
 766         $bit_max = 4;
 767     } else {
 768         $bit_max = 2;
 769     }
 770
 771     $logo_data = Maasha::Seq::seqlogo_calc( $bit_max, $entries );
 772
 773     $svg = Maasha::Plot::svg_init();
 774
 775     svg_draw_logo( $svg, $logo_data, $bit_max, $type );
 776     svg_draw_logo_scale( $svg, $bit_max );
 777
 778     return $svg->xmlify;
 779 }
 780
 781
 782 sub svg_init
 783 {
 784     # Martin A. Hansen, October 2005.
 785
 786     # inititalizes SVG object, which is returned.
 787
 788     my $svg;
 789
 790     $svg = SVG->new(
 791         style  => {
 792             'font-weight'  => 'normal',
 793             'font-family'  => 'Courier New',
 794             'font-size'    => 10,
 795         },
 796     );
 797
 798     return $svg;
 799 }
 800
 801
 802 sub svg_draw_logo
 803 {
 804     # Martin A. Hansen, January 2007.
 805
 806     # Renders a sequence logo in SVG using a
 807     # given data structure with logo details.
 808
 809     my ( $svg,         # SVG object,
 810          $logo_data,   # data structure
 811          $bit_max,     # maximum bit height
 812          $type,        # sequence type
 813          $nocolor,     # render black and white - OPTIONAL
 814        ) = @_;
 815
 816     my ( $pos, $elem, $char, $char_height_bit, $char_height_px, $block, $x, $y, $scale_factor, $color );
 817
 818     $x = 0;
 819
 820     foreach $pos ( @{ $logo_data } )
 821     {
 822         $y = 30;
 823
 824         foreach $elem ( @{ $pos } )
 825         {
 826             ( $char, $char_height_bit ) = @{ $elem };
 827
 828             $char_height_px = $char_height_bit * ( 30 / $bit_max );
 829
 830             $block = $svg->group(
 831                 transform => "translate($x,$y)",
 832             );
 833
 834             $scale_factor = $char_height_px / 7;
 835
 836             if ( $nocolor ) {
 837                 $color = "black";
 838             } elsif ( $type eq "dna" or $type eq "rna" ) {
 839                 $color = Maasha::Seq::color_nuc( $char );
 840             } else {
 841                 $color = Maasha::Seq::color_pep( $char );
 842             }
 843
 844             $block->text(
 845                 transform => "scale(1,$scale_factor)",
 846                 x  => 0,
 847                 y  => 0,
 848                 style => {
 849                     'font-weight' => 'bold',
 850                     fill          => Maasha::Seq::color_palette( $color ),
 851                 }
 852             )->cdata( $char );
 853
 854             $y -= $char_height_px;
 855         }
 856
 857         $x += 7;
 858     }
 859 }
 860
 861
 862 sub svg_draw_logo_scale
 863 {
 864     # Martin A. Hansen, January 2007.
 865
 866     # draws the bit scale for the sequence logo
 867
 868     my ( $svg,         # SVG object,
 869          $bit_max,     # maximum bit height
 870        ) = @_;
 871
 872     my ( $scale, $i );
 873
 874     $scale = $svg->group(
 875         transform => "translate(-10)",
 876         style => {
 877             stroke      => 'black',
 878             'font-size' => '8px',
 879         }
 880     );
 881
 882     $svg->text(
 883 #        transform => "translate(0,$logo_y)",
 884         transform => "rotate(-90)",
 885         x  => -26,
 886         y  => -30,
 887         style => {
 888             stroke => 'none',
 889         }
 890     )->cdata( "bits" );
 891
 892     $scale->line(
 893         x1 => 0,
 894         x2 => 0,
 895         y1 => 0,
 896         y2 => 30,
 897     );
 898
 899     for ( $i = 0; $i <= $bit_max; $i++ )
 900     {
 901         $scale->line(
 902             x1 => -5,
 903             x2 => 0,
 904             y1 => ( 30 / $bit_max ) * $i,
 905             y2 => ( 30 / $bit_max ) * $i,
 906         );
 907
 908         $scale->text(
 909             x  => -13,
 910             y  => ( 30 / $bit_max ) * $i + 2,
 911             style => {
 912                 stroke => 'none',
 913             }
 914         )->cdata( $bit_max - $i );
 915     }
 916 }
 917
 918
 919 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<