code_perl/Maasha/Plot.pm

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