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         $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";
 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";
 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";
 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     @nums = sort { $a <=> $b } @nums;
 546
 547     $cent = ( $nums[ 1 ] + $nums[ 2 ] ) / 2;
 548
 549     return $cent;
 550 }
 551
 552
 553 sub draw_chr
 554 {
 555     # Martin A. Hansen, December 2003.
 556
 557     # draws a whole cromosome with or without centromeric region
 558
 559     my ( $svg,         # image object
 560          $x,           # x position
 561          $y,           # y position
 562          $chr_len,     # lenght of chromosome
 563          $chr_width,   # width of chromosome
 564          $chr_cent,    # position of centromeric region
 565          $chr,         # chromosome
 566          $karyo_list,  # hashref with karyo data
 567          $feat_list,   # hashref with features
 568        ) = @_;
 569
 570     # returns image object
 571
 572     my ( $chr_obj, $clip_obj, $gr_obj );
 573
 574     $chr_obj = $svg->group(
 575         id => $chr,
 576     );
 577
 578     if ( exists $feat_list->{ $chr } ) {
 579         draw_chr_feat( $chr_obj, $x, $y, $chr_width, $feat_list->{ $chr } );
 580     }
 581
 582     $clip_obj = $chr_obj->clipPath(
 583         id => $chr . "_clipPath",
 584     );
 585
 586     $clip_obj->rectangle(
 587         x      => sprintf( "%.3f", $x ),
 588         y      => sprintf( "%.3f", $y ),
 589         width  => sprintf( "%.3f", $chr_width ),
 590         height => sprintf( "%.3f", $chr_cent ),
 591         rx     => 10,
 592         ry     => 10,
 593     );
 594
 595     $clip_obj->rectangle(
 596         x      => sprintf( "%.3f", $x ),
 597         y      => sprintf( "%.3f", $y + $chr_cent ),
 598         width  => sprintf( "%.3f", $chr_width ),
 599         height => sprintf( "%.3f", $chr_len - $chr_cent ),
 600         rx     => 10,
 601         ry     => 10,
 602     );
 603
 604     $gr_obj = $chr_obj->group(
 605         "clip-path"  => "url(#$chr" . "_clipPath)",
 606     );
 607
 608     if ( exists $karyo_list->{ $chr } ) {
 609         draw_karyo_data( $gr_obj, $x, $y, $chr_width, $karyo_list->{ $chr } );
 610     }
 611
 612     $gr_obj->rectangle(
 613         x      => sprintf( "%.3f", $x ),
 614         y      => sprintf( "%.3f", $y ),
 615         width  => sprintf( "%.3f", $chr_width ),
 616         height => sprintf( "%.3f", $chr_cent ),
 617         fill   => 'none',
 618         rx     => 10,
 619         ry     => 10,
 620     );
 621
 622     $gr_obj->rectangle(
 623         x      => sprintf( "%.3f", $x ),
 624         y      => sprintf( "%.3f", $y + $chr_cent ),
 625         width  => sprintf( "%.3f", $chr_width ),
 626         height => sprintf( "%.3f", $chr_len - $chr_cent ),
 627         fill   => 'none',
 628         rx     => 10,
 629         ry     => 10,
 630     );
 631
 632     draw_chr_num( $chr_obj, $x, $y, $chr_len, $chr_width, $chr );
 633 }
 634
 635
 636 sub draw_chr_num
 637 {
 638     # Martin A. Hansen, December 2003.
 639
 640     # draws a cromosome number
 641
 642     my ( $svg,        # image object
 643          $x,         # x position
 644          $y,         # y position
 645          $chr_len,   # lenght of chromosome
 646          $chr_width, # width of chromosome
 647          $chr,       # chromosome number
 648        ) = @_;
 649
 650     # returns image object
 651
 652     my ( $chr_num, $chars, @a, $word_width );
 653
 654     $chr_num = $chr;
 655     $chr_num =~ s/chr//;
 656
 657     $chars = @a = split "", $chr_num;
 658
 659     $word_width = ( $chars * 8 ) / 2;
 660
 661     $svg->text(
 662         x  => sprintf("%.3f", $x + ( $chr_width / 2 ) - $word_width ),
 663         y  => sprintf("%.3f", $y + $chr_len + 15 ),
 664     )->cdata( $chr_num );
 665 }
 666
 667
 668 sub draw_karyo_data
 669 {
 670     # Martin A. Hansen, February 2004.
 671
 672     # Plots chromosome features
 673
 674     my ( $svg,
 675          $x,
 676          $y,
 677          $chr_width,
 678          $list,
 679        ) = @_;
 680
 681     # returns an image object
 682
 683     my ( $feat_beg, $feat_end, $feat_height, $i, $color, $label );
 684
 685     for ( $i = 0; $i < @{ $list }; $i++ )
 686     {
 687         ( $label, $feat_beg, $feat_end, $color ) = @{ $list->[ $i ] };
 688
 689         $feat_height = $feat_end - $feat_beg;
 690
 691         $svg->rectangle(
 692             x      => sprintf("%.3f", $x ),
 693             y      => sprintf("%.3f", $y + $feat_beg ),
 694             width  => sprintf("%.3f", $chr_width ),
 695             height => sprintf("%.3f", $feat_height ),
 696             'stroke-width' => 0,
 697             fill   => $color,
 698         );
 699     }
 700 }
 701
 702
 703 sub draw_chr_feat
 704 {
 705     # Martin A. Hansen, February 2004.
 706
 707     # Plots chromosome features
 708
 709     my ( $svg,
 710          $x,
 711          $y,
 712          $chr_width,
 713          $list,
 714        ) = @_;
 715
 716     # returns an image object
 717
 718     my ( $feat_beg, $feat_end, $feat_height, $i, $color, $height, $width, $x1, $y1, %lookup );
 719
 720     for ( $i = 0; $i < @{ $list }; $i++ )
 721     {
 722         ( $feat_beg, $feat_end, $color ) = @{ $list->[ $i ] };
 723
 724         $feat_height = $feat_end - $feat_beg;
 725
 726         $x1     = sprintf("%.0f", $x + ( $chr_width / 2 ) ),
 727         $y1     = sprintf("%.0f", $y + $feat_beg ),
 728         $width  = sprintf("%.0f", ( $chr_width / 2 ) + 5 ),
 729         $height = sprintf("%.0f", $feat_height );
 730
 731         if ( $height < 1 )
 732         {
 733             $height = 1;
 734
 735             if ( exists $lookup{ $x1 . $y1 } ) {
 736                 next;
 737             } else {
 738                 $lookup{ $x1 . $y1 } = 1;
 739             }
 740         }
 741
 742         $svg->rectangle(
 743             x      => $x1,
 744             y      => $y1,
 745             width  => $width,
 746             height => $height,
 747             stroke => $color,
 748             fill   => $color,
 749         );
 750     }
 751 }
 752
 753
 754 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> SEQUENCE LOGO <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 755
 756
 757 sub seq_logo
 758 {
 759     # Martin A. Hansen, August 2007.
 760
 761     # Calculates and renders a sequence logo in SVG format.
 762
 763     my ( $entries,   # aligned sequence entries - list of tuples
 764        ) = @_;
 765
 766     # Returns string.
 767
 768     my ( $type, $bit_max, $logo_data, $svg );
 769
 770     $type = Maasha::Seq::seq_guess_type( $entries->[ 0 ]->[ 1 ] );
 771
 772     if ( $type =~ /^p/i ) {
 773         $bit_max = 4;
 774     } else {
 775         $bit_max = 2;
 776     }
 777
 778     $logo_data = Maasha::Seq::seqlogo_calc( $bit_max, $entries );
 779
 780     $svg = Maasha::Plot::svg_init();
 781
 782     svg_draw_logo( $svg, $logo_data, $bit_max, $type );
 783     svg_draw_logo_scale( $svg, $bit_max );
 784
 785     return $svg->xmlify;
 786 }
 787
 788
 789 sub svg_init
 790 {
 791     # Martin A. Hansen, October 2005.
 792
 793     # inititalizes SVG object, which is returned.
 794
 795     my $svg;
 796
 797     $svg = SVG->new(
 798         style  => {
 799             'font-weight'  => 'normal',
 800             'font-family'  => 'Courier New',
 801             'font-size'    => 10,
 802         },
 803     );
 804
 805     return $svg;
 806 }
 807
 808
 809 sub svg_draw_logo
 810 {
 811     # Martin A. Hansen, January 2007.
 812
 813     # Renders a sequence logo in SVG using a
 814     # given data structure with logo details.
 815
 816     my ( $svg,         # SVG object,
 817          $logo_data,   # data structure
 818          $bit_max,     # maximum bit height
 819          $type,        # sequence type
 820          $nocolor,     # render black and white - OPTIONAL
 821        ) = @_;
 822
 823     my ( $pos, $elem, $char, $char_height_bit, $char_height_px, $block, $x, $y, $scale_factor, $color );
 824
 825     $x = 0;
 826
 827     foreach $pos ( @{ $logo_data } )
 828     {
 829         $y = 30;
 830
 831         foreach $elem ( @{ $pos } )
 832         {
 833             ( $char, $char_height_bit ) = @{ $elem };
 834
 835             $char_height_px = $char_height_bit * ( 30 / $bit_max );
 836
 837             $block = $svg->group(
 838                 transform => "translate($x,$y)",
 839             );
 840
 841             $scale_factor = $char_height_px / 7;
 842
 843             if ( $nocolor ) {
 844                 $color = "black";
 845             } elsif ( $type eq "dna" or $type eq "rna" ) {
 846                 $color = Maasha::Seq::color_nuc( $char );
 847             } else {
 848                 $color = Maasha::Seq::color_pep( $char );
 849             }
 850
 851             $block->text(
 852                 transform => "scale(1,$scale_factor)",
 853                 x  => 0,
 854                 y  => 0,
 855                 style => {
 856                     'font-weight' => 'bold',
 857                     fill          => Maasha::Seq::color_palette( $color ),
 858                 }
 859             )->cdata( $char );
 860
 861             $y -= $char_height_px;
 862         }
 863
 864         $x += 7;
 865     }
 866 }
 867
 868
 869 sub svg_draw_logo_scale
 870 {
 871     # Martin A. Hansen, January 2007.
 872
 873     # draws the bit scale for the sequence logo
 874
 875     my ( $svg,         # SVG object,
 876          $bit_max,     # maximum bit height
 877        ) = @_;
 878
 879     my ( $scale, $i );
 880
 881     $scale = $svg->group(
 882         transform => "translate(-10)",
 883         style => {
 884             stroke      => 'black',
 885             'font-size' => '8px',
 886         }
 887     );
 888
 889     $svg->text(
 890 #        transform => "translate(0,$logo_y)",
 891         transform => "rotate(-90)",
 892         x  => -26,
 893         y  => -30,
 894         style => {
 895             stroke => 'none',
 896         }
 897     )->cdata( "bits" );
 898
 899     $scale->line(
 900         x1 => 0,
 901         x2 => 0,
 902         y1 => 0,
 903         y2 => 30,
 904     );
 905
 906     for ( $i = 0; $i <= $bit_max; $i++ )
 907     {
 908         $scale->line(
 909             x1 => -5,
 910             x2 => 0,
 911             y1 => ( 30 / $bit_max ) * $i,
 912             y2 => ( 30 / $bit_max ) * $i,
 913         );
 914
 915         $scale->text(
 916             x  => -13,
 917             y  => ( 30 / $bit_max ) * $i + 2,
 918             style => {
 919                 stroke => 'none',
 920             }
 921         )->cdata( $bit_max - $i );
 922     }
 923 }
 924
 925
 926 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<