]> git.donarmstrong.com Git - biopieces.git/blob - code_perl/Maasha/Plot.pm
disabled quality score conversion in FASTQ and Solexa
[biopieces.git] / 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 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<