package Bio::Graphics::Glyph::vista_plot;

use strict;
use base qw(Bio::Graphics::Glyph::wiggle_xyplot 
            Bio::Graphics::Glyph::wiggle_density
            Bio::Graphics::Glyph::wiggle_whiskers
            Bio::Graphics::Glyph::heat_map 
            Bio::Graphics::Glyph::smoothing);


our $VERSION = '1.0';

sub my_options {
    {
        start_color =>  [
            'color',
            'white',
            'Beginning of the color gradient, expressed as a named color or',
            'RGB hex string.'],
        end_color   => [
            'color',
            'red',
            'End of the color gradient.'],
        min_peak => [
            'integer',
            1,
            "Minimum value of the peak feature's \"score\" attribute."],
        max_peak => [
            'integer',
            255,
            "Maximum value of the peak feature's \"score\" attribute."],
        min_score => [
            'integer',
            undef,
            "Minimum value of the signal graph feature's \"score\" attribute."],
        max_score => [
            'integer',
            undef,
            "Maximum value of the signal graph feature's \"score\" attribute."],
        peakwidth => [
            'integer',
            3,
            "Line width determine the thickness of the line representing a peak."],
        glyph_subtype => [
	    ['peaks+signal','peaks','signal','density'],
            'vista',
            "What to show, peaks or signal, both (vista plot) or density graph."],
        graph_type => [
	     ['whiskers','histogram','boxes','line','points','linepoints'],
            'boxes',
            "Type of signal graph to show."],
	alpha  => [
	    'integer',
	    100,
	    "Alpha transparency of peak regions",
	],
    };
}

sub my_description {
    return <<END;
This glyph draws peak calls (features with discreet boundaries,
i.e. putative transcription sites, over signal graph (wiggle_xyplot)
requires a special load gff file that uses attributes 'wigfile' and 'peak_type'
BigWig support is available also, see POD documentation for more details

Example:

2L   chip_seq  vista    5407   23011573    .     .     .     Name=Chip-Seq Experiment 1;wigfile=SomeWigFile.wigdb;peak_type=transcript_region:exp1

END
}


BEGIN {
  no strict 'refs';

  my @subs = qw/ h_start   s_start   v_start h_range s_range  v_range
                 min_peak_score max_peak_score low_rgb low_hsv high_rgb peak_score_range/;

  for my $sub ( @subs ) {
    *{$sub} = sub {
      my ($self, $v) = @_;
      my $k = "_$sub";

      if (defined $v) {
        $self->{$k} = $v;
      }

      return $self->{$k};
    }
  }
}

sub peakwidth {
  shift->option('peakwidth') || 3;
}

sub alpha_c {
    my $self = shift;
    return $self->option('alpha') || 100;
}

# Need to override wiggle_xyplot padding function to enable adequate height control in density mode
sub pad_top {
  my $self = shift;
  return 0 if $self->glyph_subtype eq 'density';
  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 bigwig_summary {
    my $self = shift;
    my $d    = $self->{bigwig_summary};
    $self->{bigwig_summary} = shift if @_;
    $d;
}

# Need to override this too b/c we need unconventional mean and stdev calculation
sub global_mean_and_variance {
    my $self = shift;
    if (my $wig = $self->wig) {
	my @result = eval{($wig->mean,$wig->stdev)};
        return @result if @result;
    } 
    if (my $sum = $self->bigwig_summary){
        use Bio::DB::BigWig qw(binMean binStdev);
        my $stats = $sum->statistical_summary(1);
        return eval{(binMean($stats->[0]),binStdev($stats->[0]))};
    }
    return;
}

sub glyph_subtype {
    my $self = shift;
    my $only_show = $self->option('only_show') || $self->option('glyph_subtype') || 'vista';
    $only_show    = 'vista' if $only_show eq 'both' || $only_show eq 'peaks+signal';
    return $only_show;
}

sub graph_type {
    my $self = shift;
    return $self->option('graph_type') || $self->options->{graph_type}[1];
}

# 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 $only_show = $self->glyph_subtype;
    my $feature   = $self->feature;
 
    # Draw dual graph if we have both types of attributes, BigWig and wiggle format supported
    my %features = (wig  => (eval{$feature->get_tag_values('wigfile')})[0],
		    peak => (eval{$feature->get_tag_values('peak_type')})[0],
		    fasta=> (eval{$feature->get_tag_values('fasta')})[0]);
    $self->panel->startGroup($gd);

    $self->draw_signal($only_show,\%features,@_) if $only_show =~ /signal|density|vista/;
    $self->draw_peaks(\%features,@_)             if $features{peak} && $only_show =~ /peaks|vista|both/;
    $self->Bio::Graphics::Glyph::xyplot::draw_label(@_)       if $self->option('label');
    $self->draw_description(@_) if $self->option('description');
    $self->panel->endGroup($gd);
}

# this should be refactored from wiggle_xyplot and wiggle_density
sub draw_signal {
    my $self        = shift;
    my $signal_type = shift;
    my $paths       = shift;

    my $feature  = $self->feature;

    # Signal Graph drawing:
    if ($paths->{wig} && $paths->{wig}=~/\.wi\w{1,3}$/) {
	eval "require Bio::Graphics::Wiggle" unless Bio::Graphics::Wiggle->can('new');
	my $wig = eval { Bio::Graphics::Wiggle->new($paths->{wig}) };
	$self->wig($wig);
	$self->_draw_wigfile($feature,$wig,@_);
    } elsif ($paths->{wig} && $paths->{wig}=~/\.bw$/i) { 
	eval "use Bio::DB::BigWig 'binMean'" unless Bio::DB::BigWig->can('new');
	my @args = (-bigwig => "$paths->{wig}");
	if ($paths->{fasta}) {
	    eval "use Bio::DB::Sam" unless Bio::DB::Sam::Fai->can('open');
	    my $fasta_accessor = Bio::DB::Sam::Fai->can('open') ? Bio::DB::Sam::Fai->open("$paths->{fasta}")
		                                                : Bio::DB::Fasta->new("$paths->{fasta}");
	    push @args,(-fasta  => $fasta_accessor);
	}
	my $bigwig = Bio::DB::BigWig->new(@args);
	$self->wig($bigwig);
	my ($summary) = $bigwig->features(-seq_id => $feature->segment->ref,
					  -start  => $self->panel->start,
					  -end    => $self->panel->end,
					  -type   => 'summary');

	local $self->{feature} = $summary;
	if ($signal_type ne 'density' and  $self->graph_type eq 'whiskers') {
	    $self->Bio::Graphics::Glyph::wiggle_whiskers::draw(@_);
	} else {
	    my $stats = $summary->statistical_summary($self->width);
		my $interval_method = $self->option('interval_method') || 'mean';
		my @vals;
		if ($interval_method eq 'mean') {
			@vals  = map {$_->{validCount} ? Bio::DB::BigWig::binMean($_) : undef} @$stats;
		}
		elsif ($interval_method eq 'sum') {
			@vals  = map {$_->{validCount} ? $_->{sumData} : undef} @$stats;
		}
		elsif ($interval_method eq 'min') {
			@vals  = map {$_->{validCount} ? $_->{minVal} : undef} @$stats;
		}
		elsif ($interval_method eq 'max') {
			@vals  = map {$_->{validCount} ? $_->{maxVal} : undef} @$stats;
		}
		else {
			warn "unrecognized interval method $interval_method!";
		}
	    $self->bigwig_summary($summary);
	    if ($signal_type eq 'density') {
		$self->Bio::Graphics::Glyph::wiggle_density::_draw_coverage($summary,\@vals,@_);
            } else {
		$self->Bio::Graphics::Glyph::wiggle_data::_draw_coverage($summary,\@vals,@_);
            }
	}
    }
}

sub draw_peaks {
    my $self = shift;
    my $paths = shift;
    my($gd,$dx,$dy) = @_;
    my($left,$top,$right,$bottom) = $self->calculate_boundaries($dx,$dy);

    # Peak drawing:
    my $alpha_c = $self->alpha_c;
    my $feature = $self->feature;
    
    my $p_type = $paths->{peak};
    my @peaks = $self->peaks();
    my $x_scale     = $self->scale;
    my $panel_start = $self->panel->start;
    my $f_start     = $feature->start > $panel_start
	? $feature->start
	: $panel_start;  
    my $lw  = $self->peakwidth;
    my($max_s,$min_s) = ($self->option('max_peak'),$self->option('min_peak'));
    $max_s          = 255 if !defined $max_s;
    $min_s          = 1   if !defined $min_s;
    my $grad_ok = 0;
    if (defined $max_s && defined $min_s) {
	$grad_ok = $self->calculate_gradient($min_s,$max_s);
    }

    my $flip     = $self->{flip};

    $self->{peak_cache} = [];  # remember coordinates of the peaks

    foreach my $peak (@peaks) {
	my $x1     = $left    + ($peak->{start} - $f_start) * $x_scale;
	my $x2     = $left    + ($peak->{stop}  - $f_start) * $x_scale;
	if ($x2 >= $left and $x1 <= $right) {
	    my $y1     = $top;
	    my $y2     = $bottom; 
	    $x1        = $left   if $x1 < $left;
	    $x2        = $right  if $x2 > $right;
	    $alpha_c = $alpha_c <=127 ? $alpha_c : 0; # Reset to zero if illegal value is passed
	    my $score = $peak->{score};
	    if ($score eq "."){$score = 255;} # Set score to 255 if peak is unscored 
	    my $color;
	    if ($grad_ok && defined $score && $score!=255) {
		my @rgb = $self->Bio::Graphics::Glyph::heat_map::calculate_color($score,
										 $self->min_peak_score,
										 $self->max_peak_score,
										 $self->peak_score_range);
		$color = $self->color_index(@rgb);
	    }else{
		$color = $self->fgcolor;
	    }

	    my $bgcolor = $self->bgcolor;
		
	    if ($alpha_c > 0){
		$gd->alphaBlending(1);
		$bgcolor = $self->add_alpha($gd,$bgcolor,$alpha_c);
	    }

	    if ($flip) {
		$x1 = $right - ($x1-$left);
		$x2 = $right - ($x2-$left);
		($x1,$x2) = ($x2,$x1);
	    }
	    
	    my @rect = (int($x1+0.5),int($y1+0.5),int($x2+0.5),int($y2+0.5));
	    $self->filled_box($gd,@rect,$bgcolor,$bgcolor,0.5) if abs($y2-$y1) > 0;
	    $gd->setThickness($lw);
	    $gd->line(int($x1+0.5),int($y1+0.5),int($x2+0.5),int($y1+0.5),$color);
	    $gd->setThickness(1);
	    push @{$self->{peak_cache}},[$peak,@rect];
	}
    }
}

# Adding alpha channel to a color:
sub add_alpha {
 my($self,$im,$color,$alpha) = @_;
 my($r,$g,$b) = $im->rgb($color);
 return $im->colorAllocateAlpha($r,$g,$b,$alpha);
}

# Slightly modified function from heat_map.pm
sub calculate_gradient {
  my($self, $min, $max) = @_;
  my $start_color = lc $self->option('start_color') || 'white';
  my $stop_color  = lc $self->option('end_color')   || 'red';
  my $hsv_start   = $self->color2hsv($start_color);
  my $hsv_stop    = $self->color2hsv($stop_color);

  my ($h_start,$s_start,$v_start) = @$hsv_start;
  my ($h_stop,$s_stop,$v_stop )   = @$hsv_stop;

  my $s_range = abs($s_stop - $s_start);
  my $v_range = abs($v_stop - $v_start);

  my $h_range;
  # special case: if start hue = end hue, we want to go round
  # the whole wheel once. Otherwise round the wheel clockwise
  # or counterclockwise depending on start and end coordinate
  if ($h_start != $h_stop) {
   my $direction = abs($h_stop - $h_start)/($h_stop - $h_start);
   my ($sstart,$sstop) = sort {$a <=> $b} ($h_start,$h_stop);
   $direction *= -1 if $sstop - $sstart > 256/2; #reverse the direction if we cross 0
   $h_range = ($sstop - $sstart) <= 256/2 ? ($sstop - $sstart)*$direction : (256 - $sstop + $sstart)*$direction;
  }
  else {
   $h_range = 256;
  }
 # darkness or monochrome gradient?
  if ( !_isa_color($start_color) || !_isa_color($stop_color) ) {
    # hue (H) is fixed
    $h_range = 0;

    #    gradient         S       V    
    # white -> color    0->255   255
    # color -> white    255->0   255
    # white -> black    0        255->0
    # black -> white    0        0->255
    # black -> color    0->255   0->255
    # color -> black    255->0   255->0
    if ( $start_color eq 'white' && _isa_color($stop_color) ) {
      $s_range = 255;
      $s_start = 0;
      $v_range = 0;
      $v_start = 255;
      $h_start = $h_stop;
    }
    elsif ( _isa_color($start_color) && $stop_color eq 'white' ) {
      $s_range = -255;
      $s_start = 255;
      $v_range = 0;
      $v_start = 255;
    }
    elsif ( $start_color eq 'white' ) { # end black
      $s_range = 0;
      $s_start = 0;
      $v_range = -255;
      $v_start = 255;
    }
    elsif ( $stop_color eq 'white' ) { # start black
      $s_range = 0;
      $s_start = 0;
      $v_range = 255;
      $v_start = 0;
    }
    elsif ( _isa_color($start_color) ) { # end black
      $s_range = 255;
      $s_start = 0;
      $v_range = 255;
      $v_start = 0;
    }
    elsif ( _isa_color($stop_color) ) { # start black
      $s_range = -255;
      $s_start = 255;
      $v_range = -255;
      $v_start = 255;
    }

  }

  # store gradient info
  $self->h_range($h_range);
  $self->h_start($h_start);
  $self->s_start($s_start);
  $self->v_start($v_start);
  $self->s_range($s_range);
  $self->v_range($v_range);

  # store score info
  $self->peak_score_range($max - $min);
  $self->min_peak_score($min);
  $self->max_peak_score($max);

  # store color extremes
  my @low_rgb  = $self->HSVtoRGB(@$hsv_start);
  my @high_rgb = $self->HSVtoRGB(@$hsv_stop);
  $self->low_hsv($hsv_start);
  $self->high_rgb(\@high_rgb);
  $self->low_rgb(\@low_rgb);
  return 1;
}


sub _isa_color {
  my $color = shift;
  return $color =~ /white|black|FFFFFF|000000/i ? 0 : 1;
}

sub level { -1 }

# Need to override this so we have a nice image map for overlayed peaks
sub boxes {
    my $self = shift;
    my($left,$top,$parent) = @_;

    return if $self->glyph_subtype eq 'density'; # No boxes for density plot
    my @boxes = $self->SUPER::boxes(@_);
    
    if (my $rects = $self->{peak_cache}) {
	push @boxes,[@$_,$parent] foreach @$rects;
    }
    return wantarray ? @boxes : \@boxes;
}


# Modified and fused functions from wiggle_density.pm and wiggle_xyplot.pm
sub _draw_wigfile {
    my $self = shift;
    my $feature = shift;
    my $wig     = shift;

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

    my ($gd,$left,$top) = @_;
    my ($start,$end) = $self->effective_bounds($feature); 

    if ($self->glyph_subtype eq 'density') {
     my ($x1,$y1,$x2,$y2) = $self->bounds($left,$top);
     $self->draw_segment($gd,
                         $start,$end,
                         $wig,$start,$end,
                         1,1,
                         $x1,$y1,$x2,$y2);
     $self->Bio::Graphics::Glyph::xyplot::draw_label(@_)       if $self->option('label');
     $self->draw_description(@_) if $self->option('description');
    } else {
     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 peaks {
    my $self = shift;
    return @{$self->{_peaks}} if $self->{_peaks};

    my $feature = $self->feature;
    my $db = $feature->object_store;
    my ($p_type) = eval{$feature->get_tag_values('peak_type')};

    unless ($db && $p_type) {
	$self->{_peaks}	 = [];
	return;
    }

    my @peaks = $db->features(-seq_id => $feature->segment->ref,
			      -start  => $self->panel->start,
			      -end    => $self->panel->end,
			      -type   => $p_type); 

    $self->{_peaks} = \@peaks;
    return @{$self->{_peaks}};
}

1;

=head1 NAME

Bio::Graphics::Glyph::vista_plot - The "vista_plot" glyph

=head1 SYNOPSIS

See 
L<Bio::Graphics::Glyph>, L<Bio::Graphics::Glyph::wiggle_xyplot> and L<Bio::Graphics::Glyph::heat_map>.

=head1 DESCRIPTION

This glyph draws peak calls (features with discreet boundaries,
i.e. putative transcription sites, over signal graph (wiggle_xyplot)
requires a special load gff file that uses attributes 'wigfile' and 'peak_type'

B<Example:>

2L  chip_seq  vista  5407   23011573  .  .  .  Name=ChipSeq Exp 1;wigfile=SomeWigFile.wigdb;peak_type=binding_site:exp1

The glyph will draw the wiggle file first, than overlay the peaks (if there are any)
over signal graph. Elsewhere in the GFF3 file, there should be one or more features 
of type "binding_site:exp1", e.g.:

2L  exp1  binding_site  91934  92005  .  .  .

Options like 'balloon hover' and 'link' are available to customize
interaction with peaks in detail view.

B<BigWig support:>

Supported bigwig format also requires another attribute to be supplied
in load gff file (fasta) which specifies sequence index file for the
organism in use. The data file should have the 'bw' extension - it is
used to detect the BigWig format by vista_plot

3L  chip_seq  vista   1    24543530  .  .  .   Name=ChipSeq Exp 2;wigfile=SomeBigWigFile.bw;peak_type=binding_site:exp2;fasta=YourOrganism.fasta

Note that all attributes should be present in load gff, as the code currently does not handle situation when
only some of the attributes are in gff. To omit peak or signal drawing use "" (i.e. peak_type="")
In both cases, the stanza code will look the same (only essential parameters shown):

 [VISTA_PLOT]
 feature         = vista:chip_seq
 glyph           = vista_plot
 label           = 1 
 smoothing       = mean
 smoothing_window = 10
 bump density    = 250
 autoscale       = local
 variance_band   = 1
 max_peak        = 255
 min_peak        = 1
 peakwidth       = 3
 start_color     = lightgray
 end_color       = black
 pos_color       = blue
 neg_color       = orange
 bgcolor         = orange
 alpha           = 80
 fgcolor         = black
 database        = database_with_load_gff_data
 box_subparts    = 1
 bicolor_pivot   = min
 key             = VISTA plot 

=head1 OPTIONS

Options are the same as for wiggle_xyplot and heat_map

B<Additional parameters:>

B<alpha>
set transparency for peak area.

B<glyph_subtype>
Display only 'peaks', 'signal', 'density' or 'peaks+signal'. 
Aliases for 'peaks+signal' include "both" and "vista".

B<Recommended global settings:>

for proper peak drawing transparency should be enabled
by setting 
B<truecolors=1> 
in I<GBrowse.conf> file

=head1 BUGS

Please report them.

=head1 SEE ALSO

L<Bio::Graphics::Panel>
L<Bio::Graphics::Glyph>
L<Bio::Graphics::Glyph::wiggle_xyplot>
L<Bio::Graphics::Glyph::heat_map>
L<GD>        

=head1 AUTHOR

Peter Ruzanov pruzanov@oicr.on.ca

Copyright (c) 2010 Ontario Institute for Cancer Research

 This package and its accompanying libraries is free software; you can
 redistribute it and/or modify it under the terms of the GPL (either
 version 1, or at your option, any later version) or the Artistic
 License 2.0.  Refer to LICENSE for the full license text. In addition,
 please see DISCLAIMER.txt for disclaimers of warranty.

=cut