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