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