3 # Copyright (C) 2007 Martin A. Hansen.
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.
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.
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.
19 # http://www.gnu.org/copyleft/gpl.html
22 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> DESCRIPTION <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
25 # Routines to plot stuff with Gnuplot and SVG.
28 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
39 use vars qw ( @ISA @EXPORT );
47 @ISA = qw( Exporter );
50 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
53 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> LINEPLOTS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
58 # Martin A. Hansen, January 2008.
60 # Plots a simple lineplot using Gnuplot.
62 my ( $data, # data table - each column will be plottet as one line.
63 $options, # options hash
64 $tmp_dir, # temporary directory
69 my ( $tmp_file, $pid, $fh_in, $fh_out, $cmd, $i, $line, @lines, $xtic_space, @plot_cmd, $title );
71 $tmp_dir ||= $ENV{ 'BP_TMP' };
73 $tmp_file = "$tmp_dir/lineplot_simple.tab";
75 $fh_out = Maasha::Filesys::file_write_open( $tmp_file );
77 map { print $fh_out join( "\t", @{ $_ } ), "\n" } @{ $data };
81 $options->{ "terminal" } ||= "dumb";
83 $cmd = "gnuplot -persist";
85 $pid = open2( $fh_out, $fh_in, $cmd );
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";
98 for ( $i = 1; $i < scalar @{ $data->[ 0 ] } + 1; $i++ )
100 $title = $options->{ 'keys' }->[ $i - 1 ] || $options->{ 'list' } || "";
101 push @plot_cmd, qq("$tmp_file" using $i with lines title "$title");
104 print $fh_in "plot " . join( ", ", @plot_cmd ) . "\n";
108 while ( $line = <$fh_out> )
121 return wantarray ? @lines : \@lines;
125 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> HISTOGRAMS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
130 # Martin A. Hansen, August 2007.
132 # Plots a simple histogram using Gnuplot.
134 my ( $data, # list of [ xlabel, data value ] tuples
135 $options, # options hash
140 my ( $pid, $fh_in, $fh_out, $cmd, $i, $line, @lines );
142 $options->{ "terminal" } ||= "dumb";
144 $cmd = "gnuplot -persist";
146 $pid = open2( $fh_out, $fh_in, $cmd );
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";
164 for ( $i = 0; $i < @{ $data }; $i++ )
166 print $fh_in join( "\t", "\"$data->[ $i ]->[ 0 ]\"", $data->[ $i ]->[ 1 ] ), "\n";
171 while ( $line = <$fh_out> )
182 return wantarray ? @lines : \@lines;
186 sub histogram_lendist
188 # Martin A. Hansen, August 2007.
190 # Plots a histogram using Gnuplot.
192 my ( $data, # list of [ xlabel, data value ] tuples
193 $options, # options hash
198 my ( $pid, $fh_in, $fh_out, $cmd, $i, $line, @lines, $xtic_space );
200 $options->{ "terminal" } ||= "dumb";
202 if ( $data->[ -1 ]->[ 0 ] <= 10 ) { # FIXME: some day figure the formula for this!
204 } elsif ( $data->[ -1 ]->[ 0 ] <= 100 ) {
206 } elsif ( $data->[ -1 ]->[ 0 ] <= 250 ) {
208 } elsif ( $data->[ -1 ]->[ 0 ] <= 500 ) {
210 } elsif ( $data->[ -1 ]->[ 0 ] <= 1000 ) {
212 } elsif ( $data->[ -1 ]->[ 0 ] <= 2500 ) {
214 } elsif ( $data->[ -1 ]->[ 0 ] <= 5000 ) {
216 } elsif ( $data->[ -1 ]->[ 0 ] <= 10000 ) {
218 } elsif ( $data->[ -1 ]->[ 0 ] <= 50000 ) {
220 } elsif ( $data->[ -1 ]->[ 0 ] <= 100000 ) {
222 } elsif ( $data->[ -1 ]->[ 0 ] <= 1000000 ) {
226 $cmd = "gnuplot -persist";
228 $pid = open2( $fh_out, $fh_in, $cmd );
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";
242 print $fh_in "plot '-' using 1\n";
244 for ( $i = 0; $i < @{ $data }; $i++ )
246 $data->[ $i ]->[ 0 ] = "." if $data->[ $i ]->[ 0 ] % 10 != 0;
248 print $fh_in join( "\t", $data->[ $i ]->[ 1 ] ), "\n";
253 while ( $line = <$fh_out> )
264 return wantarray ? @lines : \@lines;
268 sub histogram_chrdist
270 # Martin A. Hansen, August 2007.
272 # Plots a histogram using Gnuplot.
274 my ( $data, # list of [ xlabel, data value ] tuples
275 $options, # options hash
280 my ( $pid, $fh_in, $fh_out, $cmd, $i, $line, @lines );
282 $options->{ "terminal" } ||= "dumb";
284 $cmd = "gnuplot -persist";
286 $pid = open2( $fh_out, $fh_in, $cmd );
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";
300 print $fh_in "plot '-' using 2:xticlabels(1)\n";
302 for ( $i = 0; $i < @{ $data }; $i++ ) {
303 print $fh_in join( "\t", "\"$data->[ $i ]->[ 0 ]\"", $data->[ $i ]->[ 1 ] ), "\n";
308 while ( $line = <$fh_out> )
319 return wantarray ? @lines : \@lines;
323 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> KARYOGRAM <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
328 # Martin A. Hansen, August 2007.
330 # Plot hits on a karyogram for a given genome.
332 my ( $data, # list of [ chr, beg, end ] triples
333 $options, # hashref with options
338 my ( $karyo_file, $svg, $features, $karyo );
340 if ( $options->{ "genome" } eq "hg18" )
342 $karyo_file = $ENV{ 'BP_DIR' } . "/bp_data/human_cytobands.txt";
344 elsif( $options->{ "genome" } eq "mm9" )
346 $karyo_file = $ENV{ 'BP_DIR' } . "/bp_data/mouse_cytobands.txt";
349 $karyo = parse_karyo_data( $karyo_file );
353 chromosome_layout( $svg, $karyo, $data );
361 # X q26.1 129700001 130200000 gneg
363 # color: /etc/X11/rgb.txt
368 my ( $fh, $chr, $line, $name, $beg, $end, $color, %features, %color_hash );
374 gpos75 => "DarkGrey",
375 gpos66 => "DarkGrey",
377 gpos33 => "LightGrey",
378 gpos25 => "LightGrey",
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)",
392 $fh = Maasha::Filesys::file_read_open( $file );
394 while ( $line = <$fh> )
398 next if $line =~ /^#/;
400 # ( $chr, $name, $beg, $end, $color ) = split "\t", $line;
401 ( $chr, $beg, $end, $name, $color ) = split "\t", $line;
403 # if ( $color =~ /^gpos(\d+)/ ) {
404 # $color = color_intensity( $1 );
405 # } elsif ( exists $color_hash{ $color } ) {
406 $color = $color_hash{ $color };
408 # die qq(ERROR: Unknown color->$color\n);
411 if ( exists $features{ $chr } )
413 push @{ $features{ $chr } }, [ $name, $beg, $end, $color ];
417 $features{ $chr } = [ [ $name, $beg, $end, $color ] ];
423 return wantarray ? %features : \%features;
429 # Martin A. Hansen, September 2007.
431 # Converts a gray scale intensity in percent to rgb.
433 my ( $percent, # color intensity
440 $num = int( $percent * 256 / 100 );
444 $hex = sprintf "%x", $num;
446 return "#$hex$hex$hex";
448 # return "rgb($num,$num,$num)";
454 # Martin A. Hansen, September 2005.
456 # initializes svg image.
458 # returns an image object
474 sub chromosome_layout
476 # Martin A. Hansen, January 2004 - August 2007.
478 # Plots all chromosomes in a single
480 my ( $svg, # image object
481 $karyo_list, # hashref with karyo data
482 $feat_list, # hashref with features
485 # returns an image object
487 my ( $layout_obj, $i, $x, $y, $max, $factor, $chr_len, $chr_width, $chr_cent, $chr, $feat, $karyo, @list, $A, $B );
489 $layout_obj = $svg->group(
493 $max = $karyo_list->{ "chr1" }->[ -1 ]->[ 2 ];
494 $factor = ( HEIGHT / 2 ) / $max;
495 $chr_width = ( HEIGHT / 4 ) / 13;
497 foreach $karyo ( keys %{ $karyo_list } ) {
498 map { $_->[ 1 ] *= $factor; $_->[ 2 ] *= $factor } @{ $karyo_list->{ $karyo } };
501 foreach $feat ( keys %{ $feat_list } ) {
502 map { $_->[ 0 ] *= $factor; $_->[ 1 ] *= $factor } @{ $feat_list->{ $feat } };
505 @list = sort { $A = $a; $B = $b; $A =~ s/chr//; $B =~ s/chr//; $A <=> $B } keys %{ $karyo_list };
508 push @list, "chrX", "chrY";
515 $chr_len = $karyo_list->{ $chr }->[ -1 ]->[ 2 ];
516 $chr_cent = find_cent( $karyo_list->{ $list[ $i ] } );
518 $y = HEIGHT / 2 - $chr_len;
519 $x = ( WIDTH / ( @list + 2 ) ) * ( $i + 1 );
521 draw_chr( $layout_obj, $x, $y, $chr_len, $chr_width, $chr_cent, $chr, $karyo_list, $feat_list );
530 # Martin A. Hansen, December 2003.
532 # Finds the centromeric region in the karyo data.
536 my ( $acen, @nums, $cent );
538 @{ $acen } = grep { grep { /^DarkGrey$/ } @{ $_ } } @{ $list };
540 push @nums, $acen->[ 0 ]->[ 1 ];
541 push @nums, $acen->[ 0 ]->[ 2 ];
542 push @nums, $acen->[ 1 ]->[ 1 ];
543 push @nums, $acen->[ 1 ]->[ 2 ];
546 @nums = grep { defined $_ } @nums; # FIXME
547 @nums = sort { $a <=> $b } @nums;
549 $cent = ( $nums[ 1 ] + $nums[ 2 ] ) / 2;
557 # Martin A. Hansen, December 2003.
559 # draws a whole cromosome with or without centromeric region
561 my ( $svg, # image object
564 $chr_len, # lenght of chromosome
565 $chr_width, # width of chromosome
566 $chr_cent, # position of centromeric region
568 $karyo_list, # hashref with karyo data
569 $feat_list, # hashref with features
572 # returns image object
574 my ( $chr_obj, $clip_obj, $gr_obj );
576 $chr_obj = $svg->group(
580 if ( exists $feat_list->{ $chr } ) {
581 draw_chr_feat( $chr_obj, $x, $y, $chr_width, $feat_list->{ $chr } );
584 $clip_obj = $chr_obj->clipPath(
585 id => $chr . "_clipPath",
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 ),
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 ),
606 $gr_obj = $chr_obj->group(
607 "clip-path" => "url(#$chr" . "_clipPath)",
610 if ( exists $karyo_list->{ $chr } ) {
611 draw_karyo_data( $gr_obj, $x, $y, $chr_width, $karyo_list->{ $chr } );
615 x => sprintf( "%.3f", $x ),
616 y => sprintf( "%.3f", $y ),
617 width => sprintf( "%.3f", $chr_width ),
618 height => sprintf( "%.3f", $chr_cent ),
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 ),
634 draw_chr_num( $chr_obj, $x, $y, $chr_len, $chr_width, $chr );
640 # Martin A. Hansen, December 2003.
642 # draws a cromosome number
644 my ( $svg, # image object
647 $chr_len, # lenght of chromosome
648 $chr_width, # width of chromosome
649 $chr, # chromosome number
652 # returns image object
654 my ( $chr_num, $chars, @a, $word_width );
659 $chars = @a = split "", $chr_num;
661 $word_width = ( $chars * 8 ) / 2;
664 x => sprintf("%.3f", $x + ( $chr_width / 2 ) - $word_width ),
665 y => sprintf("%.3f", $y + $chr_len + 15 ),
666 )->cdata( $chr_num );
672 # Martin A. Hansen, February 2004.
674 # Plots chromosome features
683 # returns an image object
685 my ( $feat_beg, $feat_end, $feat_height, $i, $color, $label );
687 for ( $i = 0; $i < @{ $list }; $i++ )
689 ( $label, $feat_beg, $feat_end, $color ) = @{ $list->[ $i ] };
691 $feat_height = $feat_end - $feat_beg;
694 x => sprintf("%.3f", $x ),
695 y => sprintf("%.3f", $y + $feat_beg ),
696 width => sprintf("%.3f", $chr_width ),
697 height => sprintf("%.3f", $feat_height ),
707 # Martin A. Hansen, February 2004.
709 # Plots chromosome features
718 # returns an image object
720 my ( $feat_beg, $feat_end, $feat_height, $i, $color, $height, $width, $x1, $y1, %lookup );
722 for ( $i = 0; $i < @{ $list }; $i++ )
724 ( $feat_beg, $feat_end, $color ) = @{ $list->[ $i ] };
726 $feat_height = $feat_end - $feat_beg;
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 );
737 if ( exists $lookup{ $x1 . $y1 } ) {
740 $lookup{ $x1 . $y1 } = 1;
756 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> SEQUENCE LOGO <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
761 # Martin A. Hansen, August 2007.
763 # Calculates and renders a sequence logo in SVG format.
765 my ( $entries, # aligned sequence entries - list of tuples
770 my ( $type, $bit_max, $logo_data, $svg );
772 $type = Maasha::Seq::seq_guess_type( $entries->[ 0 ]->[ 1 ] );
774 if ( $type eq "PROTEIN" ) {
780 $logo_data = Maasha::Seq::seqlogo_calc( $bit_max, $entries );
782 $svg = Maasha::Plot::svg_init();
784 svg_draw_logo( $svg, $logo_data, $bit_max, $type );
785 svg_draw_logo_scale( $svg, $bit_max );
793 # Martin A. Hansen, October 2005.
795 # inititalizes SVG object, which is returned.
801 'font-weight' => 'normal',
802 'font-family' => 'Courier New',
813 # Martin A. Hansen, January 2007.
815 # Renders a sequence logo in SVG using a
816 # given data structure with logo details.
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
825 my ( $pos, $elem, $char, $char_height_bit, $char_height_px, $block, $x, $y, $scale_factor, $color );
829 foreach $pos ( @{ $logo_data } )
833 foreach $elem ( @{ $pos } )
835 ( $char, $char_height_bit ) = @{ $elem };
837 $char_height_px = $char_height_bit * ( 30 / $bit_max );
839 $block = $svg->group(
840 transform => "translate($x,$y)",
843 $scale_factor = $char_height_px / 7;
847 } elsif ( $type eq "dna" or $type eq "rna" ) {
848 $color = Maasha::Seq::color_nuc( $char );
850 $color = Maasha::Seq::color_pep( $char );
854 transform => "scale(1,$scale_factor)",
858 'font-weight' => 'bold',
859 fill => Maasha::Seq::color_palette( $color ),
863 $y -= $char_height_px;
871 sub svg_draw_logo_scale
873 # Martin A. Hansen, January 2007.
875 # draws the bit scale for the sequence logo
877 my ( $svg, # SVG object,
878 $bit_max, # maximum bit height
883 $scale = $svg->group(
884 transform => "translate(-10)",
887 'font-size' => '8px',
892 # transform => "translate(0,$logo_y)",
893 transform => "rotate(-90)",
908 for ( $i = 0; $i <= $bit_max; $i++ )
913 y1 => ( 30 / $bit_max ) * $i,
914 y2 => ( 30 / $bit_max ) * $i,
919 y => ( 30 / $bit_max ) * $i + 2,
923 )->cdata( $bit_max - $i );
928 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<