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