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