1package Bio::Graphics::Glyph::hybrid_plot;
2
3use strict;
4use base qw(Bio::Graphics::Glyph::wiggle_xyplot);
5use constant DEBUG=>0;
6use constant NEGCOL=>"orange";
7use constant POSCOL=>"blue";
8our $VERSION = '1.0';
9
10sub my_options {
11    {
12        min_score => [
13            'integer',
14            undef,
15            "Minimum value of the signal graph feature's \"score\" attribute."],
16        max_score => [
17            'integer',
18            undef,
19            "Maximum value of the signal graph feature's \"score\" attribute."],
20        flip_sign => [
21            'boolean',
22            0,
23            "Optionally flip the signal for wigfileB (if the scores are positive but we wish to paint the signal along negative y-axis)"]
24    };
25}
26
27sub my_description {
28    return <<END;
29This glyph draws signal graph (wiggle_xyplot) using wiggle or BigWig files
30requires a special load gff file that uses attributes 'wigfileA' and 'wigfileB'
31For BigWig support it is also required to have 'fasta' attribute set to point
32to a fasta file for your organism of interest
33
34Example:
35
362L  test   hybrid    5407   23011573    .     .     .     Name=Experiment 1;wigfileA=SomeWigFile1.wigdb;wigfileB=SomeWigFile2.wigdb
37
38END
39}
40
41#Checking the method for individual features (RNA-Seq reads)
42sub _check_uni {
43 return shift->option('u_method') || 'match';
44}
45
46
47# Override height and pad functions (needed to correctly space features with different sources):
48sub height {
49  my $self = shift;
50  my $h    = $self->SUPER::height;
51  return $self->feature->method eq $self->_check_uni ? 3 : $h;
52}
53
54sub pad_top {
55  my $self = shift;
56  return $self->feature->method eq $self->_check_uni ? 0 : 4;
57}
58
59sub pad_bottom {
60  my $self = shift;
61  return $self->feature->method eq $self->_check_uni ? 0 : 4;
62}
63
64# we override the draw method so that it dynamically creates the parts needed
65# from the wig file rather than trying to fetch them from the database
66sub draw {
67
68 my $self = shift;
69 my ($gd,$dx,$dy) = @_;
70 my ($left,$top,$right,$bottom) = $self->calculate_boundaries($dx,$dy);
71 my $height   = $bottom - $top;
72 my $feature  = $self->feature;
73 my $set_flip = $self->option('flip_sign') || 0;
74
75 #Draw individual features for reads (unlike wiggle features reads will have scores)
76 my $t_id = $feature->method;
77 if($t_id && $t_id eq $self->_check_uni){return Bio::Graphics::Glyph::generic::draw_component($self,@_);}
78
79 #Draw multiple graph if we don't have a score
80 my @wiggles = $self->get_wiggles($feature);
81
82 my ($fasta)   = $feature->get_tag_values('fasta');
83 my($scale,$y_origin,$min_score,$max_score);
84
85 $self->panel->startGroup($gd);
86
87 #Depending on what we have (wiggle or BigWig) pick the way to paint the signal graph
88 for(my $w = 0; $w < @wiggles; $w++){
89     if ($w > 0) {
90	 $self->configure('bgcolor', NEGCOL);
91	 $self->configure('no_grid', 1);
92     } else {
93	 $self->configure('bgcolor', POSCOL);
94     }
95     if ($wiggles[$w] =~ /\.wi\w{1,3}$/) {
96	 $self->draw_wigfile($feature,$wiggles[$w],@_);
97     } elsif ($wiggles[$w] =~ /\.bw$/) {
98	 my $flip = ($w > 0 && $set_flip) ? -1 : 1;
99	 eval "require Bio::DB::BigWig;1" or die $@;
100	 eval "require Bio::DB::Sam; 1"   or die $@;
101	 my @args = (-bigwig => "$wiggles[$w]");
102	 push @args,(-fasta  => Bio::DB::Sam::Fai->open($fasta)) if $fasta;
103	 my $wig = Bio::DB::BigWig->new(@args);
104	 my ($summary) = $wig->features(-seq_id => $feature->segment->ref,
105					-start  => $self->panel->start,
106					-end    => $self->panel->end,
107					-type   => 'summary');
108	 my $stats = $summary->statistical_summary($self->width);
109	 my $interval_method = $self->option('interval_method') || 'mean';
110	 my @vals;
111	 if ($interval_method eq 'mean') {
112		@vals  = map {$_->{validCount} ? $_->{sumData}/$_->{validCount} * $flip : undef} @$stats;
113	 }
114	 elsif ($interval_method eq 'sum') {
115		@vals  = map {$_->{validCount} ? $_->{sumData} * $flip : undef} @$stats;
116	 }
117	 elsif ($interval_method eq 'min') {
118		@vals  = map {$_->{validCount} ? $_->{minVal} * $flip : undef} @$stats;
119	 }
120	 elsif ($interval_method eq 'max') {
121		@vals  = map {$_->{validCount} ? $_->{maxVal} * $flip : undef} @$stats;
122	 }
123	 else {
124		warn "unrecognized interval method $interval_method!";
125	 }
126	 $self->_draw_coverage($summary,\@vals,@_);
127     }
128 }
129}
130
131sub get_wiggles {
132    my $self = shift;
133    my $feature = shift;
134    my @wiggles;
135    foreach ('A'..'Z') {
136	my $filename = 'wigfile'.$_;
137	my ($wiggle) = $feature->get_tag_values('wigfile'.$_);
138	push (@wiggles, $wiggle) if $wiggle;
139    }
140    return @wiggles;
141}
142
143sub minmax {
144    my $self   = shift;
145    my $parts  = shift;
146
147    my $autoscale  = $self->option('autoscale') || 'local';
148    my $set_flip = $self->option('flip_sign') || 0;
149
150    my $min_score  = $self->min_score  unless $autoscale eq 'z_score';
151    my $max_score  = $self->max_score  unless $autoscale eq 'z_score';
152
153    my $do_min     = !defined $min_score;
154    my $do_max     = !defined $max_score;
155
156    my @wiggles = $self->get_wiggles($self->feature);
157    my ($min,$max,$mean,$stdev);
158    my @args = (-seq_id => (eval{$self->feature->segment->ref}||''),
159		-start  => $self->panel->start,
160		-end    => $self->panel->end,
161		-type   => 'summary');
162
163    for my $w (@wiggles) {
164	my ($a,$b,$c,$d);
165	if ($w =~ /\.bw$/) {
166	    eval "require Bio::DB::BigWig;1" or die $@;
167	    my $wig = Bio::DB::BigWig->new(-bigwig=>$w) or next;
168	    ($a,$b,$c,$d) = $self->bigwig_stats($autoscale,$wig->features(@args));
169	} elsif ($w =~ /\.wi\w{1,3}$/) {
170	    eval "require Bio::Graphics::Wiggle;1" or die $@;
171	    my $wig = Bio::Graphics::Wiggle->new($w);
172	    ($a,$b,$c,$d) = $self->wig_stats($autoscale,$wig);
173	}
174	$min    = $a if !defined $min || $min > $a;
175	$max    = $b if !defined $max || $max < $b;
176	$mean  += $c;
177	$stdev += $d**2;
178    }
179    $stdev = sqrt($stdev);
180    $min = $max * -1 if ($set_flip);
181
182    $min_score = $min if $do_min;
183    $max_score = $max if $do_max;
184    return $self->sanity_check($min_score,$max_score,$mean,$stdev);
185}
186
1871;
188__END__
189
190=head1 NAME
191
192Bio::Graphics::Glyph::hybrid_plot - An xyplot plot drawing dual graph using data from two or more wiggle files per track
193
194=head1 SYNOPSIS
195
196See <Bio::Graphics::Panel> <Bio::Graphics::Glyph> and <Bio::Graphics::Glyph::wiggle_xyplot>.
197
198=head1 DESCRIPTION
199
200Note that for full functionality this glyph requires Bio::Graphics::Glyph::generic (generic glyph is used for drawing individual
201matches for small RNA alignments at a high zoom level, specified by semantic zooming in GBrowse conf file)
202Unlike the regular xyplot, this glyph draws two overlapping graphs
203using value data in Bio::Graphics::Wiggle file format:
204
205track type=wiggle_0 name="Experiment" description="snRNA seq data" visibility=pack viewLimits=-2:2 color=255,0,0 altColor=0,0,255 windowingFunction=mean smoothingWindow=16
206
207 2L 400 500 0.5
208 2L 501 600 0.5
209 2L 601 700 0.4
210 2L 701 800 0.1
211 2L 800 900 0.1
212
213##gff-version 3
214
2152L      Sample_rnaseq  rnaseq_wiggle 41   3009 . . . ID=Samlpe_2L;Name=Sample;Note=YourNoteHere;wigfileA=/datadir/track_001.2L.wig;wigfileB=/datadir/track_002.2L.wig
216
217
218The "wigfileA" and "wigfileB" attributes give a relative or absolute pathname to
219Bio::Graphics::Wiggle format files for two concurrent sets of data. Basically,
220these wigfiles contain the data on signal intensity (counts) for sequences
221aligned with genomic regions. In wigfileA these data are additive, so for each
222sequence region the signal is calculated as a sum of signals from overlapping
223matches (signal). In wigfileB the signal represents the maximum value among all
224sequences (signal quality) aligned with the current region so the user can see
225the difference between accumulated signal from overlapping multiple matches
226(which may likely be just a noise from products of degradation) and high-quality
227signal from unique sequences.
228
229For a third wiggle file use the attribute "wigfileC" and so forth.
230
231It is essential that wigfile entries in gff file do not have score, because
232score used to differentiate between data for dual graph and data for matches
233(individual features visible at higher magnification). After an update to
234wiggle_xyplot code colors for dual plot are now hard-coded (blue for signal and
235orange for signal quality). Alpha channel is also handled by wiggle_xyplot code now.
236
237=head2 OPTIONS
238
239In addition to some of the wiggle_xyplot glyph options, the following options are
240recognized:
241
242 Name        Value        Description
243 ----        -----        -----------
244
245 wigfileA    path name    Path to a Bio::Graphics::Wiggle file for accumulated vales in 10-base bins
246
247 wigfileB    path name    Path to a Bio::Graphics::Wiggle file for max values in 10-base bins
248
249 fasta       path name    Path to fasta file to enable BigWig drawing
250
251 u_method    method name  Use method of [method name] to identify individual features (like alignment matches)
252                          to show at high zoom level. By default it is set to 'match'
253
254=head1 BUGS
255
256 Please report them.
257
258=head1 SEE ALSO
259
260L<Bio::Graphics::Panel>,
261L<Bio::Graphics::Glyph>,
262L<Bio::Graphics::Glyph::arrow>,
263L<Bio::Graphics::Glyph::cds>,
264L<Bio::Graphics::Glyph::crossbox>,
265L<Bio::Graphics::Glyph::diamond>,
266L<Bio::Graphics::Glyph::dna>,
267L<Bio::Graphics::Glyph::dot>,
268L<Bio::Graphics::Glyph::ellipse>,
269L<Bio::Graphics::Glyph::extending_arrow>,
270L<Bio::Graphics::Glyph::generic>,
271L<Bio::Graphics::Glyph::graded_segments>,
272L<Bio::Graphics::Glyph::heterogeneous_segments>,
273L<Bio::Graphics::Glyph::line>,
274L<Bio::Graphics::Glyph::pinsertion>,
275L<Bio::Graphics::Glyph::primers>,
276L<Bio::Graphics::Glyph::rndrect>,
277L<Bio::Graphics::Glyph::segments>,
278L<Bio::Graphics::Glyph::ruler_arrow>,
279L<Bio::Graphics::Glyph::toomany>,
280L<Bio::Graphics::Glyph::transcript>,
281L<Bio::Graphics::Glyph::transcript2>,
282L<Bio::Graphics::Glyph::translation>,
283L<Bio::Graphics::Glyph::allele_tower>,
284L<Bio::DB::GFF>,
285L<Bio::SeqI>,
286L<Bio::SeqFeatureI>,
287L<Bio::Das>,
288L<GD>
289
290=head1 AUTHOR
291
292Peter Ruzanov E<lt>pruzanov@oicr.on.caE<gt>.
293
294Copyright (c) 2008 Ontario Institute for Cancer Research
295
296This package and its accompanying libraries is free software; you can
297redistribute it and/or modify it under the terms of the GPL (either
298version 1, or at your option, any later version) or the Artistic
299License 2.0.  Refer to LICENSE for the full license text. In addition,
300please see DISCLAIMER.txt for disclaimers of warranty.
301
302=cut
303