Bio::Graphics::Glyph::pairplot - The "pairwise plot" glyph

 use Bio::Graphics;

 # create the panel, etc.  See Bio::Graphics::Panel
 # for the synopsis

 # Create one big feature using the PairFeature
 # glyph (see end of synopsis for an implementation)
 my $block = PairFeature->new(-start=>  2001,
 			      -end  => 10000);

 # It will contain a series of subfeatures.
 my $start = 2001;
 while ($start < 10000) {
   my $end = $start+120;
   $block->add_SeqFeature($bsg->new(-start=>$start,
				    -end  =>$end
				   ),'EXPAND');
   $start += 200;
 }

 $panel->add_track($block,
 		   -glyph => 'pairplot',
		   -angle => 45,
		   -bgcolor => 'red',
		   -point => 1,
		  );

 print $panel->png;

 package PairFeature;
 use base 'Bio::SeqFeature::Generic';

 sub pair_score {
   my $self = shift;
   my ($sf1,$sf2) = @_;
   # simple distance function
   my $dist  = $sf2->end    - $sf1->start;
   my $total = $self->end   - $self->start;
   return sprintf('%2.2f',1-$dist/$total);
 }

  Option      Description                  Default
  ------      -----------                  -------

  -point      If true, the plot will be         0
              drawn relative to the
              midpoint between each adjacent
              subfeature.  This is appropriate
              for point-like subfeatures, such
              as SNPs.

  -angle      Angle to draw the plot.  Values   45
              between 1 degree and 89 degrees
              are valid.  Higher angles give
              a more vertical plot.

  -bgcolor    The color of the plot.            cyan

package Bio::Graphics::Glyph::pairplot;

# $Id: pairplot.pm,v 1.1 2009-06-04 22:00:49 lstein Exp $

# Triangle plot for showing pairwise quantitative relationships.
# Similar to pairwise_plot, which was originally distributed in GBrowse,
# except that the signal intensity is handled using a callback that takes
# the two features to be compared.

sub my_description {
    return <<END;
This is a triangle plot that displays quantitative relationships between point-like features, 
such as LD between SNPs. The quantitative information is generated by the "score" callback,
which receives two feature objects and returns a floating point number between "min_score"
and "max_score". The bgcolor intensity will be scaled from min to max like graded_segments.

NOTE (2 April 2008): This glyph is not yet fully functional.
END
}

sub my_options {
    {
	score => [
	    'CODEREF',
	    undef,
	    'Pass a coderef to a subroutine with the following signature: sub ($$). The',
	    'two arguments are the first and second feature to compare. Return a floating point',
	    'number to indicate the score at the intersection of these two features.'
	],
	min_score => [
	    'float',
	    0.0,
	    'Minimum possible pairwise score.'
	],
	max_score => [
	    'float',
	    1.0,
	    'Maximum possible pairwise score.'
	],
        angle => [
	    'float',
	    45,
	    'Angle between the side of the triangle and the base, in degrees.'
        ],
    }
}



use strict;
use Math::Trig;

use base 'Bio::Graphics::Glyph::generic';

sub maxdepth {
  my $self = shift;
  my $md   = $self->Bio::Graphics::Glyph::maxdepth;
  return $md if defined $md;
  return 1;
}

# return angle in radians
sub angle {
  my $self  = shift;
  my $angle = $self->{angle} ||= $self->option('angle') || 45;
  $self->{angle} = shift if @_;
  deg2rad($angle);
}

sub slope {
  my $self = shift;
  return $self->{slope} if exists $self->{slope};
  return $self->{slope} = tan($self->angle);
}

sub x2y {
  my $self = shift;
  shift() * $self->slope;
}

sub intercept {
  my $self = shift;
  my ($x1,$x2) = @_;
  my $mid = ($x1+$x2)/2;
  my $y   = $self->x2y($mid-$x1);
  return (int($mid+0.5),int($y+0.5));
}

# height calculated from width
sub layout_height {
  my $self = shift;
  return $self->{height} if exists $self->{height};
  return $self->{height} = $self->x2y($self->width)/2;
}

sub min_score {
    my $self = shift;
    my $min  = $self->option('min_score');
    return   0 unless defined $min;
}
sub max_score {
    my $self = shift;
    my $max  = $self->option('max_score');
    return 1.0 unless defined $max;
}

sub calculate_color {
  my $self = shift;
  my ($s,$rgb) = @_;
  return $self->{colors}{$s} if exists $self->{colors}{$s};
  my $max_score = $self->max_score;
  my $min_score = $self->min_score;
  my $scale     = 255/($max_score-$min_score);
  my $value     = ($s-$min_score) * $scale; # will range from 0 to 255
  return $self->{colors}{$s} = 
    $self->panel->translate_color(map { 255 - $value} @$rgb);
}

sub draw {
  my $self = shift;
  my $gd   = shift;
  my ($left,$top,$partno,$total_parts) = @_;
  my $fgcolor = $self->fgcolor;

  my ($red,$green,$blue) = $self->panel->rgb($self->bgcolor);

  my @points = $self->get_points();
  $gd->line($self->left+$left, $top+1,
	    $self->right+$left,$top+1,
	    $fgcolor);

  my $points = $self->option('point');

  my @parts = sort {$a->left<=>$b->left} $self->parts;
  $_->draw_component($gd,$left,$top-10) foreach @parts;

  # assumption: parts are not overlapping
  if ($points) {
    @points = map { int (($parts[$_]->right+$parts[$_+1]->left)/2)} (0..$#parts-1);
    unshift @points,int($parts[0]->left);
    push @points,int($parts[-1]->right);
  }

  for (my $ia=0;$ia<@parts-1;$ia++) {
    for (my $ib=$ia+1;$ib<@parts;$ib++) {
      my ($l1,$r1,$l2,$r2);
      if (@points) {
	($l1,$r1) = ($points[$ia]+1,$points[$ia+1]-1);
	($l2,$r2) = ($points[$ib]+1,$points[$ib+1]-1);
      } else {
	($l1,$r1) = ($parts[$ia]->left,$parts[$ia]->right);
	($l2,$r2) = ($parts[$ib]->left,$parts[$ib]->right);
      }

      my $intensity = eval{$self->feature->pair_score($parts[$ia],$parts[$ib])};
      warn $@ if $@;
      $intensity    = 1.0 unless defined $intensity;
      my $c         = $self->calculate_color($intensity,[$red,$green,$blue]);

      # left corner
      my ($lcx,$lcy) = $self->intercept($l1,$l2);
      my ($tcx,$tcy) = $self->intercept($r1,$l2);
      my ($rcx,$rcy) = $self->intercept($r1,$r2);
      my ($bcx,$bcy) = $self->intercept($l1,$r2);
      my $poly = GD::Polygon->new();
      $poly->addPt($lcx+$left,$lcy+$top);
      $poly->addPt($tcx+$left,$tcy+$top);
      $poly->addPt($rcx+$left,$rcy+$top);
      $poly->addPt($bcx+$left,$bcy+$top);
      $gd->filledPolygon($poly,$c);
    }
  }
}

sub get_points {
  my $self = shift;
  my @points;
  my @parts = $self->parts;
  return unless @parts;

  for my $g (@parts) {
    push @points,$g->left;
    push @points,$g->right;
  }
  @points;
}

# never allow our internal parts to bump;
sub bump { 0 }

1;

__END__