Bio::Graphics::Glyph::wiggle_xyplot - An xyplot plot compatible with dense "wig"data

  See <Bio::Graphics::Panel> and <Bio::Graphics::Glyph>.

 reference = chr1
 ChipCHIP Feature1 1..10000 wigfile=./test.wig
 ChipCHIP Feature2 10001..20000 wigfile=./test.wig
 ChipCHIP Feature3 25001..35000 wigfile=./test.wig

   Name        Value        Description
   ----        -----        -----------

   basedir     path         Path to be used to resolve "wigfile" and "densefile"
                                tags giving relative paths. Default is to use the
                                current working directory. Absolute wigfile &
                                densefile paths will not be changed.

   autoscale   "local" or "global"
                             If one or more of min_score and max_score options 
                             are absent, then these values will be calculated 
                             automatically. The "autoscale" option controls how
                             the calculation is done. The "local" value will
                             scale values according to the minimum and maximum
                             values present in the window being graphed. "global"   
                             will use chromosome-wide statistics for the entire
                             wiggle or dense file to find min and max values.




   smoothing   method name  Smoothing method: one of "mean", "max", "min" or "none"

   smoothing_window 
               integer      Number of values across which data should be smoothed.

   variance_band boolean    If true, draw a grey band across entire plot showing mean
                               and +/- 1 standard deviation (for wig files only).

   bicolor_pivot
               name         Where to pivot the two colors when drawing bicolor plots.
                               Options are "mean" and "zero". A numeric value can
                               also be provided.

   pos_color   color        When drawing bicolor plots, the fill color to use for values
                              that are above the pivot point.

   neg_color   color        When drawing bicolor plots, the fill color to use for values
                              that are below the pivot point.

   Name        Value        Description
   ----        -----        -----------

   wigfile     path name    Path to the Bio::Graphics::Wiggle file or object
                            for quantitative values.

   wigdata     string       Data exported from a Bio::Graphics::Wiggle in WIF
                            format using its export_to_wif64() method.

   densefile   path name    Path to a Bio::Graphics::DenseFeature object
                               (deprecated)

   denseoffset integer      Integer offset to where the data begins in the
                               Bio::Graphics::DenseFeature file (deprecated)

   densesize   integer      Integer size of the data in the Bio::Graphics::DenseFeature
                               file (deprecated)

package Bio::Graphics::Glyph::wiggle_xyplot;

use strict;
use base qw(Bio::Graphics::Glyph::wiggle_minmax 
            Bio::Graphics::Glyph::xyplot 
            Bio::Graphics::Glyph::smoothing);
use IO::File;
use File::Spec;


sub my_description {
    return <<END;
This glyph draws quantitative data as an xyplot. It is designed to be
used in conjunction with features in "wiggle" format as generated by
Bio::Graphics::Wiggle, or base pair coverage data generated by the
Bio::DB::Sam module.

For this glyph to work, the feature must define one of the following tags:

  wigfile -- a path to a Bio::Graphics::Wiggle file

  wigdata -- Wiggle data in the Bio::Graphics::Wiggle "wif" format, as created
             by \$wig->export_to_wif().

  coverage-- a simple comma-delimited string containing the quantitative values,
             assumed to be one value per pixel.

END
}

sub my_options {
    {
	basedir => [
	    'string',
	    undef,
	    'If a relative path is used for "wigfile", then this option provides',
	    'the base directory on which to resolve the path.'
	    ],
	variance_band => [
	    'boolean',
	    0,
	    'If true, draw a semi-transparent band across the image that indicates',
	    'the mean and standard deviation of the data set. Only of use when a wig',
	    'file is provided.'
        ],
	z_score_bounds => [
	    'integer',
            4,
	    'When using z_score autoscaling, this option controls how many standard deviations',
	    'above and below the mean to show.'
	],
	autoscale => [
	    ['local','chromosome','global','z_score','clipped_global'],
            'clipped_global',
	    'If set to "global" , then the minimum and maximum values of the XY plot',
	    'will be taken from the wiggle file as a whole. If set to "chromosome", then',
            'scaling will be to minimum and maximum on the current chromosome.',
	    '"clipped_global" is similar to "global", but clips the top and bottom values',
	    'to the multiples of standard deviations indicated by "z_score_bounds"',
	    'If set to "z_score", then the whole plot will be rescaled to z-scores in which',
	    'the "0" value corresponds to the mean across the genome, and the units correspond',
	    'to standard deviations above and below the mean. The number of SDs to show are',
	    'controlled by the "z_score_bound" option.',
	    'Otherwise, the plot will be',
	    'scaled to the minimum and maximum values of the region currently on display.',
	    'min_score and max_score override autoscaling if one or both are defined'
        ],
    };
}

# Added pad_top subroutine (pad_top of Glyph.pm, which is called when executing $self->pad_top
# returns 0, so we need to override it here)
sub pad_top {
  my $self = shift;
  my $pad = $self->Bio::Graphics::Glyph::generic::pad_top(@_);
  if ($pad < ($self->font('gdTinyFont')->height)) {
    $pad = $self->font('gdTinyFont')->height;  # extra room for the scale
  }
  $pad;
}

sub pad_left {
    my $self = shift;
    my $pad  = $self->SUPER::pad_left(@_);
    return $pad unless $self->option('variance_band');
    $pad    += length('+1sd')/2 * $self->font('gdTinyFont')->width+3;
    return $pad;
}

# we override the draw method so that it dynamically creates the parts needed
# from the wig file rather than trying to fetch them from the database
sub draw {
  my $self = shift;
  my ($gd,$dx,$dy) = @_;

  my $feature     = $self->feature;
  my ($wigfile)   = eval{$feature->get_tag_values('wigfile')};
  return $self->draw_wigfile($feature,$self->rel2abs($wigfile),@_) if $wigfile;

  my ($wigdata) = eval{$feature->get_tag_values('wigdata')};
  return $self->draw_wigdata($feature,$wigdata,@_) if $wigdata;

  my ($densefile) = eval{$feature->get_tag_values('densefile')};
  return $self->draw_densefile($feature,$self->rel2abs($densefile),@_) if $densefile;

  my ($coverage)  = eval{$feature->get_tag_values('coverage')};
  return $self->draw_coverage($feature,$coverage,@_) if $coverage;

  # support for BigWig/BigBed
  if ($feature->can('statistical_summary')) {
      my $stats = $feature->statistical_summary($self->width);
      $stats   ||= [];
      my @vals  = map {$_->{validCount} ? $_->{sumData}/$_->{validCount}:0} @$stats;
      return $self->draw_coverage($feature,\@vals,@_);
  }

  return $self->SUPER::draw(@_);
}

sub draw_wigfile {
  my $self = shift;
  my $feature = shift;
  my $wigfile = shift;

  eval "require Bio::Graphics::Wiggle" unless Bio::Graphics::Wiggle->can('new');
  my $wig = ref $wigfile && $wigfile->isa('Bio::Graphics::Wiggle') 
      ? $wigfile
      : eval { Bio::Graphics::Wiggle->new($wigfile) };
  unless ($wig) {
      warn $@;
      return $self->SUPER::draw(@_);
  }
  $self->_draw_wigfile($feature,$wig,@_);
}

sub draw_wigdata {
    my $self    = shift;
    my $feature = shift;
    my $data    = shift;

    if (ref $data eq 'ARRAY') {
	my ($start,$end) = $self->effective_bounds($feature);
	my $parts = $self->subsample($data,$start,$end);
	$self->draw_plot($parts,@_);
    }

    else {
	my $wig = eval { Bio::Graphics::Wiggle->new() };
	unless ($wig) {
	    warn $@;
	    return $self->SUPER::draw(@_);
	}

	$wig->import_from_wif64($data);
	$self->_draw_wigfile($feature,$wig,@_);
    }
}

sub draw_coverage {
    my $self    = shift;
    my $feature = shift;
    my $array   = shift;

    $array      = [split ',',$array] unless ref $array;
    return unless @$array;

    my ($start,$end)    = $self->effective_bounds($feature);
    my $bases_per_bin   = ($end-$start)/@$array;
    my $pixels_per_base = $self->scale;
    my @parts;
    for (my $pixel=0;$pixel<$self->width;$pixel++) {
	my $offset = $pixel/$pixels_per_base;
	my $s      = $start + $offset;
	my $e      = $s+1;  # fill in gaps
	my $v      = $array->[$offset/$bases_per_bin];
	push @parts,[$s,$s,$v];
    }
    $self->draw_plot(\@parts,@_);
}

sub _draw_wigfile {
    my $self    = shift;
    my $feature = shift;
    my $wig     = shift;

    $wig->smoothing($self->get_smoothing);
    $wig->window($self->smooth_window);

    my ($start,$end) = $self->effective_bounds($feature);

    $self->wig($wig);
    my $parts = $self->create_parts_for_dense_feature($wig,$start,$end);
    $self->draw_plot($parts,@_);
}

sub effective_bounds {
    my $self    = shift;
    my $feature = shift;
    my $panel_start = $self->panel->start;
    my $panel_end   = $self->panel->end;
    my $start       = $feature->start>$panel_start 
                         ? $feature->start 
                         : $panel_start;
    my $end         = $feature->end<$panel_end   
                         ? $feature->end   
                         : $panel_end;
    return ($start,$end);
}

sub draw_plot {
    my $self            = shift;
    my $parts           = shift;
    my ($gd,$dx,$dy)    = @_;
    my $pivot = $self->bicolor_pivot;

    $self->panel->startGroup($gd);

    my ($left,$top,$right,$bottom) = $self->calculate_boundaries($dx,$dy);

    # There is a minmax inherited from xyplot as well as wiggle_minmax, and I don't want to
    # rely on Perl's multiple inheritance DFS to find the right one.
    my ($min_score,$max_score,$mean,$stdev)     = $self->Bio::Graphics::Glyph::wiggle_minmax::minmax($parts);
    my $rescale  = $self->option('autoscale') eq 'z_score';
    my $side    = $self->_determine_side();

    my ($scaled_min,$scaled_max);
    if ($rescale) {
	my $bound  = $self->z_score_bound;
	$scaled_min = -$bound;
	$scaled_max = +$bound;
    }
    elsif ($side) {
	$scaled_min = int($min_score - 0.5);
	$scaled_max = int($max_score + 0.5);
    } else {
	($scaled_min,$scaled_max) = ($min_score,$max_score);
    }

    my $height = $bottom - $top;
    my $y_scale  = $scaled_max > $scaled_min ? $height/($scaled_max-$scaled_min)
	                                   : 1;
    my $x = $left;
    my $y = $top;
    
    my $x_scale     = $self->scale;
    my $panel_start = $self->panel->start;
    my $feature     = $self->feature;
    my $f_start     = $feature->start > $panel_start 
	                  ? $feature->start 
			  : $panel_start;

    $y += $self->pad_top;

    # position of "0" on the scale
    my $y_origin = $scaled_min <= 0 && $pivot ne 'min' ? $bottom - (0 - $scaled_min) * $y_scale : $bottom;
    $y_origin    = int($y_origin+0.5);

    $self->panel->startGroup($gd);
    $self->_draw_grid($gd,$x_scale,$scaled_min,$scaled_max,$dx,$dy,$y_origin) unless ($self->option('no_grid') == 1);
    $self->panel->endGroup($gd);

    return unless $scaled_max > $scaled_min;

    my $lw       = $self->linewidth;
    my $positive = $self->pos_color;
    my $negative = $self->neg_color;
    my $midpoint = $self->midpoint;
    my $flip     = $self->{flip};
    $midpoint    = ($midpoint - $mean)/$stdev if $rescale;

    my @points = map {
	my ($start,$end,$score) = @$_;
	$score     = ($score-$mean)/$stdev if $rescale;
	my $x1     = $left    + ($start - $f_start) * $x_scale;
	my $x2     = $left    + ($end   - $f_start) * $x_scale;
	if ($x2 >= $left and $x1 <= $right) {
	    my $y1     = $bottom  - ($score - $scaled_min) * $y_scale;
	    my $y2     = $y_origin;
	    $y1        = $top    if $y1 < $top;
	    $y1        = $bottom if $y1 > $bottom;
	    $x1        = $left   if $x1 < $left;
	    $x2        = $right  if $x2 > $right;

	    $x1        = $right - ($x1-$left) if $flip;
	    $x2        = $right - ($x2-$left) if $flip;
 
	    my $color = $score > $midpoint ? $positive : $negative;
	    [int($x1+0.5),int($y1+0.5),int($x2+0.5),int($y2+0.5),$color,$lw];
	} else {
	    ();
	}
    } @$parts;

    my $svg_workaround = $gd->isa('GD::SVG::Image');

    $self->panel->startGroup($gd);
    my $type           = $self->graph_type;
    if ($type eq 'boxes') {
	for (@points) {
	    my ($x1,$y1,$x2,$y2,$color,$lw) = @$_;
	    next unless abs($y2-$y1) > 0;
	    if ($svg_workaround && $x1==$x2) {
		$gd->line($x1,$y1,$x2,$y2,$color);
	    } else {
		$gd->filledRectangle($x1,$y1,$x2,$y2,$color);
	    }
	}
    }

    if ($type eq 'line' or $type eq 'linepoints') {
	my $current = shift @points;
	my $lw      = $self->option('linewidth');
	$gd->setThickness($lw) if $lw > 1;
	for (@points) {
	    my ($x1,$y1,$x2,$y2,$color,$lw) = @$_;
	    $gd->line(@{$current}[0,1],@{$_}[0,1],$color);
	    $current = $_;
	}
	$gd->setThickness(1);
    }

    if ($type eq 'points' or $type eq 'linepoints') {
	my $symbol_name = $self->option('point_symbol') || 'point';
	my $filled      = $symbol_name =~ s/^filled_//;
	my $symbol_ref  = $self->symbols->{$symbol_name};
	my $pr          = $self->point_radius;
	for (@points) {
	    my ($x1,$y1,$x2,$y2,$color,$lw) = @$_;
	    $symbol_ref->($gd,$x1,$y1,$pr,$color,$filled);
	}
    }

    if ($type eq 'histogram') {
	my $current = shift @points;
	for (@points) {
	    my ($x1, $y1, $x2, $y2, $color, $lw)  = @$_;
	    my ($y_start,$y_end) = $y1 < $y_origin ? ($y1,$y_origin) : ($y_origin,$y1);
	    if ($y1-$y2) {
		my $delta = abs($x2-$current->[0]);
		$gd->filledRectangle($current->[0],$y_start,$x2,$y_end,$color) if $delta > 1;
		$gd->line($current->[0],$y_start,$current->[0],$y_end,$color)  if $delta == 1;
		$current = $_;
	    }
	}	
    }

    if ($self->option('variance_band') && 
	(my ($mean,$variance) = $self->global_mean_and_variance())) {
	if ($rescale) {
	    $mean     = 0;
	    $variance = 1;
	}
	my $y1             = $bottom - ($mean+$variance   - $scaled_min) * $y_scale;
	my $y2             = $bottom - ($mean-$variance   - $scaled_min) * $y_scale;
	my $yy1            = $bottom - ($mean+$variance*2 - $scaled_min) * $y_scale;
	my $yy2            = $bottom - ($mean-$variance*2 - $scaled_min) * $y_scale;
	my ($clip_top,$clip_bottom);
	if ($y1 < $top) {
	    $y1                = $top;
	    $clip_top++;
	}
	if ($y2 > $bottom) {
	    $y2                = $bottom;
	    $clip_bottom++;
	}
	my $y              = $bottom - ($mean - $scaled_min) * $y_scale;
	my $mean_color     = $self->panel->translate_color('yellow:0.80');
	my $onesd_color = $self->panel->translate_color('grey:0.30');
	my $twosd_color = $self->panel->translate_color('grey:0.20');
	$gd->filledRectangle($left,$y1,$right,$y2,$onesd_color);
	$gd->filledRectangle($left,$yy1,$right,$yy2,$twosd_color);
	$gd->line($left,$y,$right,$y,$mean_color);

	my $side = $self->_determine_side();
	my $fcolor=$self->panel->translate_color('grey:0.50');
	my $font  = $self->font('gdTinyFont');
	my $x1    = $left - length('+2sd') * $font->width - ($side=~/left|three/ ? 15 : 0);
	my $x2    = $left - length('mn')   * $font->width - ($side=~/left|three/ ? 15 : 0);
	$gd->string($font,$x1,$yy1-$font->height/2,'+2sd',$fcolor) unless $clip_top;
	$gd->string($font,$x1,$yy2-$font->height/2,'-2sd',$fcolor) unless $clip_bottom;
	$gd->string($font,$x2,$y -$font->height/2,'mn',  $fcolor);
    }
    $self->panel->endGroup($gd);

    $self->panel->startGroup($gd);
    $self->_draw_scale($gd,$x_scale,$scaled_min,$scaled_max,$dx,$dy,$y_origin);
    $self->panel->endGroup($gd);

    $self->Bio::Graphics::Glyph::xyplot::draw_label(@_)       if $self->option('label');
    $self->draw_description(@_) if $self->option('description');

  $self->panel->endGroup($gd);
}

sub draw_label {
    my $self = shift;
    my ($gd,$left,$top,$partno,$total_parts) = @_;
    return $self->Bio::Graphics::Glyph::xyplot::draw_label(@_) unless $self->option('variance_band');
    return $self->Bio::Graphics::Glyph::xyplot::draw_label($gd,$left,$top,$partno,$total_parts);
}


sub draw_densefile {
    my $self = shift;
    my $feature = shift;
    my $densefile = shift;
    
    my ($denseoffset) = eval{$feature->get_tag_values('denseoffset')};
    my ($densesize)   = eval{$feature->get_tag_values('densesize')};
    $denseoffset ||= 0;
    $densesize   ||= 1;
    
    my $smoothing      = $self->get_smoothing;
    my $smooth_window  = $self->smooth_window;
    my $start          = $self->smooth_start;
    my $end            = $self->smooth_end;

    my $fh         = IO::File->new($densefile) or die "can't open $densefile: $!";
    eval "require Bio::Graphics::DenseFeature" unless Bio::Graphics::DenseFeature->can('new');
    my $dense = Bio::Graphics::DenseFeature->new(-fh=>$fh,
						 -fh_offset => $denseoffset,
						 -start     => $feature->start,
						 -smooth    => $smoothing,
						 -recsize   => $densesize,
						 -window    => $smooth_window,
	) or die "Can't initialize DenseFeature: $!";
    my $parts = $self->create_parts_for_dense_feature($dense,$start,$end);
    $self->draw_plot($parts);
}

# BUG: the next two subroutines should be merged
sub create_parts_for_dense_feature {
    my $self = shift;
    my ($dense,$start,$end) = @_;

    my $span = $self->scale> 1 ? $end - $start : $self->width;
    my $data = $dense->values($start,$end,$span);
    my $points_per_span = ($end-$start+1)/$span;
    my @parts;

    for (my $i=0; $i<$span;$i++) {
	my $offset = $i * $points_per_span;
	my $value  = shift @$data;
	next unless defined $value;
	push @parts,[$start + int($i * $points_per_span),
		     $start + int($i * $points_per_span),
		     $value];
    }
    return \@parts;
}

sub subsample {
  my $self = shift;
  my ($data,$start,$end) = @_;
  my $span = $self->scale > 1 ? $end - $start 
                              : $self->width;
  my $points_per_span = ($end-$start+1)/$span;
  my @parts;
  for (my $i=0; $i<$span;$i++) {
    my $offset = $i * $points_per_span;
    my $value  = $data->[$offset + $points_per_span/2];
    push @parts,[$start + int($i*$points_per_span),
		 $start + int($i*$points_per_span),
		 $value];
  }
  return \@parts;
}

sub rel2abs {
    my $self = shift;
    my $wig  = shift;
    return $wig if ref $wig;
    my $path = $self->option('basedir');
    return File::Spec->rel2abs($wig,$path);
}


1;

__END__