Bio::Graphics::Glyph::phylo_align - The "phylogenetic alignment" glyph

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

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

  -fgcolor      Foreground color	       black

  -outlinecolor	Synonym for -fgcolor

  -bgcolor      Background color               turquoise

  -fillcolor    Synonym for -bgcolor

  -linewidth    Line width                     1

  -height       Height of glyph		       10

  -font         Glyph font		       gdSmallFont

  -connector    Connector type                 0 (false)

  -connector_color
                Connector color                black

  -label        Whether to draw a label	       0 (false)

  -description  Whether to draw a description  0 (false)

  -hilite       Highlight color                undef (no color)

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

  -draw_clado_left
              Draws the Cladogram on left 0

  -species_spacing
              Spacing of species in DNA   1
              mode in units of font height

  -species_spacing_score
              Spacing of spcies in        5
              conservation view in units
              of font height

  -hide_label Whether to label spcies     0

  -tree_file  Path of file containing     undef
              cladogram tree information
  -tree_format Format of tree file	   newick

  -axis_color Color of the vertical axes  fgcolor
              in the GC content graph

  -errcolor   Color of all misalignment   fgcolor
              indicators

  -mid_axis_color
              Color of the middle axis of
              the conservation score graph axis_color

  -clado_bg  Color of the clado bg       bgcolor
              indicators

  -ref_color  Color of base pair bg for   bgcolor
              the reference sequence

  -targ_color Color of base pair bg for   bgcolor
              all base pairs that match
              reference

package Bio::Graphics::Glyph::phylo_align;

use strict;
use base qw(Bio::Graphics::Glyph::generic Bio::Graphics::Glyph::xyplot);
use Bio::TreeIO;
use Bio::Graphics::Wiggle;
use POSIX qw(log10);

use Carp 'croak','cluck';
use Data::Dumper;

my %complement = (g=>'c',a=>'t',t=>'a',c=>'g',n=>'n',
		  G=>'C',A=>'T',T=>'A',C=>'G',N=>'N');


# turn off description
sub description { 0 }

# turn off label
# sub label { 1 }

sub height {
  my $self = shift;
  my $font = $self->font;
  
  #adjust the space to take if conservation scores are drawn instead
  if (! $self->dna_fits) {
    my $species_spacing_score = $self->option('species_spacing_score') || 5;
    $self->factory->set_option('species_spacing', $species_spacing_score);
  }
  
  my $species_spacing = $self->option('species_spacing') || 1;
  
  #Height = NumSpecies x Spacing/species x FontHeight
  my $height = ($self->known_species + $self->unknown_species + 1)
            * $species_spacing
            * $font->height;
  
  #use this if you want to show only those species that have alignments in the viewing window
  #$height = ($self->extract_features + 1) * 2 * $font->height;
  
  $self->factory->set_option('height', $height);
  
  
  return $height;
  
}

#####
# TODO: extract the wigfiles covering the range as well
#####
# get all features within the viewing window
sub extract_features {
  my $self = shift;
  #my $segment = $self->feature->{'factory'}->segment($self->feature->refseq,
  #						      $self->feature->start => $self->feature->stop);
  #my @match = $segment->features('submatch:pa'); 
  my @match = $self->feature->features('submatch:pa');
  warn $self->feature->features;
#print "Match has ",$#match,"<p>\n";
  
  # exract wifiles here too:
  my @wmatch = $self->feature->features('wfile:pa');
  #push @match, $self->feature->features('wfile:pa');
#  print "xxxxxx<pre>",Dumper(@wmatch),"</pre>cccccc";
#print "WMatch has ",$#wmatch,"<p>\n";

#print "Halfwaypoint<p>";  
#for my $feature (@match,@wmatch) {
#my %attributes = $feature->attributes;
#my $species = $attributes{'species'};
#print "<pre>Feature $species:\n",Dumper(%attributes),"</pre>\n";
#}

  my %alignments;
#  for my $feature (@match) {
  for my $feature (@match,@wmatch) {
    my %attributes = $feature->attributes;
    my $species = $attributes{'species'};
    
    push @{$alignments{$species}}, $feature;
  }
  
  %alignments;
}

#known species (all those that are in the Phylo tree)
sub known_species {
  my $self = shift;
  
  my $tree = shift;
  
  if ($tree) {
    my @leaves = $tree->get_leaf_nodes;
    my @allspecies = map {$_->id} @leaves;
    return @allspecies
    
  } else {
    #this may be too simple of an assumption, especially for non newick files
    my $tree_file = $self->option('tree_file');
    
    open (FH, $tree_file);
    my $newick = <FH>;
    close FH;
    
    my @allspecies = $newick =~ /([a-zA-Z]\w*)/g;
    return @allspecies;
  }
}

sub unknown_species {
  my $self = shift;

  my %alignments;        #all species in viewing window
  my $refspecies;        #all species from cladogram info
  my @current_species;   #all species in viewing window
  my @known_species;     #species in GFF but and clado
  my @unknown_species;   #species in GFF but not in clado
  # current - known = unknown
  
  
  if (@_) {
    %alignments = %{$_[0]};
    $refspecies = $_[1];
    @current_species = @{$_[2]};
    @known_species = @{$_[3]};
    @unknown_species;
  } else {
    %alignments = $self->extract_features;
    $refspecies = $self->option('reference_species');
    @current_species =  keys %alignments;   #all species in viewing window
    @known_species = $self->known_species;  #all species from cladogram info
    @unknown_species;                       #species in GFF but not in clado
  } #would have combined the two cases into one line using || but Perl will treat the arrays as num of elem
  
  #do set subtraction to see which species in viewing range but not in tree
  my %seen;  # build lookup table
  @seen{@known_species} = ();
  foreach my $item (@current_species, $refspecies) {
    push(@unknown_species, $item) unless exists $seen{$item};
  }
  
  return @unknown_species;
  
}

sub set_tree {
  my $self = shift;
  
  #warn"My species are ".Dumper(@species);
  
  my $tree_file   = $self->option('tree_file');
  my $tree_format = $self->option('tree_format') || 'newick';
  
  my $treeio = new Bio::TreeIO(-file   => $tree_file,
                               -format => $tree_format);
  
  
  
  my $tree = $treeio->next_tree;
  my $root = $tree->get_root_node;
  
  # would be ideal to remove all species that don't have features (alignments) within
  # viewing window but there is a bug in Bio::Tree library where you can't remove the
  # first leaf node.
  
  #  $tree->remove_Node('dog'); # etc...
  
  #set the leaf x coodinate (make all evenly spaced)
  my @leaves = $tree->get_leaf_nodes;
  for (my $i=0; $i<@leaves; $i++) {
    my $leaf = $leaves[$i];
  
    #note that leaves can use "description" functions while intermediate nodes cannot
    #thus objects must be handled directly
    $leaf->description({'x'=>$i});
  
  }
  
  
  #set root height to 0
  $root->{'description'}{'y'} = 0;
  
  #set the x and y coordinates of all intermediate nodes
  get_n_set_next_treenode($root, 0);
  
  flip_xy($tree);

  
  $tree;
  
}

sub get_max_height {
  my $tree = shift;
  my $max_height;
  
  #get the max height
  for my $child ($tree->get_leaf_nodes) {
    my $x = $child->{'_description'}{'x'};
    $max_height = $x if $max_height < $x;
  }
  
  $max_height;
}


sub draw_clado {
  my $self = shift;
  my $tree = shift;
  my $gd = shift;
  my ($x1, $y1, $x2, $y2, $color, $xscale, $yscale, $xoffset, $yoffset, $start_x, $draw_clado_left) = @_;
  
  my @bounds = $gd->getBounds;
  
  my $root = $tree->get_root_node;
  my @nodes = $root->get_all_Descendents;
  
  #draw bg for cladogram
  my $clado_bg = $self->color('clado_bg') || $self->bgcolor;
  my @coords = (0, $y1, $start_x+$xoffset+$self->font->width-1, $y2+1);
  my @coords2 = ($x1, $y1, $start_x+$xoffset/2, $y2);
  if ($draw_clado_left) {
    $gd->filledRectangle(@coords, $clado_bg);
    $gd->filledRectangle(@coords2, $self->color('clado_bg'));
    $gd->filledRectangle($x2, $x1, $bounds[0], $bounds[1], $self->color('bg_color')) if $self->dna_fits;
  } else {
    $gd->filledRectangle($bounds[0]-$coords[2], $coords[1], $bounds[0]-$coords[0], $coords[3],
  			 $self->color('bg_color'));
    $gd->filledRectangle($bounds[0]-$coords2[2], $coords2[1], $bounds[0]-$coords2[0], $coords2[3],
			 $clado_bg);  
    $gd->filledRectangle(0, $y1, $x1, $y2+1, $self->color('bg_color')) if $self->dna_fits;
  }

  
  
  #draw the lines of the tree
  for my $node ($root,@nodes) {
    next if $node->is_Leaf;
    my $x = $node->{'_description'}{'x'} * $xscale;
    
    #draw vertical line covering all children
    my $topx = $node->{'_description'}{'childmin'} * $yscale;
    my $botx = $node->{'_description'}{'childmax'} * $yscale;
    
    @coords = ($x+$xoffset, $topx+$yoffset, $x+$xoffset, $botx+$yoffset);
    if ($draw_clado_left) {
      $gd->line(@coords, $self->fgcolor);
    } else {
      $gd->line($bounds[0]-$coords[2], $coords[1], $bounds[0]-$coords[0], $coords[3], $self->fgcolor);
    }
    
    #draw a line connecting the bar to each child
    my @children = $node->each_Descendent;
    for my $child (@children) {
      my $cx = $child->{'_description'}{'x'} * $xscale;
      my $cy = $child->{'_description'}{'y'} * $yscale;
      $cx = $start_x if $child->is_Leaf;
      
      #print"($cx, $cy)";
      
      @coords = ($x+$xoffset, $cy+$yoffset, $cx+$xoffset, $cy+$yoffset);
      if ($draw_clado_left) {
        $gd->line(@coords, $self->fgcolor);
      } else {
          $gd->line($bounds[0]-$coords[2], $coords[1], $bounds[0]-$coords[0], $coords[3], $self->fgcolor);
      }
      
    }
    
  }
  
  
  
  
  #my $tree = $treeio->next_tree;
  #my @nodes = grep { $_->bootstrap > 70 } $tree->get_nodes;
  #warn"my nodes are:\n".Dumper(@nodes);
  print"</pre>";
  
  $start_x + $xscale;
  
  }

#tree is made with root on top but this will switch x and y coords so root is on left
sub flip_xy {
  my $tree = shift;
  my $root = $tree->get_root_node;
  
  my @nodes = $root->get_all_Descendents;
  
  for my $node ($root, @nodes) {
    if ($node->is_Leaf) {
      $node->{'_description'} = {
        'x'       => $node->{'_description'}{'y'},
        'y'       => $node->{'_description'}{'x'}
      }
    } else {
      $node->{'_description'} = {
        'x'         => $node->{'_description'}{'y'},
        'y'         => $node->{'_description'}{'x'},
        'child_pos' => $node->{'_description'}{'child_pos'},
        'childmin'  => $node->{'_description'}{'childmin'},
        'childmax'  => $node->{'_description'}{'childmax'}
      }
    }
  }
  
  
}


#recursive function that sets the x and y coordinates of tree nodes
sub get_n_set_next_treenode {
  my $node = shift;
  my $height = $node->{'_description'}{'y'};
  
  my @children = $node->each_Descendent;
  
  
  my $x = 0;
  my $min_child_x = -1;
  my $max_child_x = -1;
  
  #iterate through children to find the x's and set the y's (height's)
  for my $child (@children) {
    #set the y coordinate as parent's height + 1
    $child->{'_description'}{'y'} = $height + 1;
    get_n_set_next_treenode($child);
    
    #retrieve the child's x coordinate
    my $child_x = $child->{'_description'}{'x'} || 0;
    $x += $child_x;
    
    $min_child_x = $child_x if $min_child_x==-1 || $child_x < $min_child_x;
    $max_child_x = $child_x if $max_child_x==-1 || $max_child_x < $child_x;
    #$x += $child->discription->{'x'};  #cannot do this for intermediate nodes
    
    #store the x values of all children
    push @{$node->{'_description'}{'child_pos'}}, $child_x;
    
    
  }
  
  #set the current x coordinate as the average of all children's x's
  if (@children) {
    $x = $x / @children;
    $node->{'_description'}{'x'} = $x;
    $node->{'_description'}{'childmin'} = $min_child_x;
    $node->{'_description'}{'childmax'} = $max_child_x;

  }
  
  $node->{'_description'}{'y'} = $height;
   
}

sub get_legend_and_scale {
  my $yscale = shift;
  my $height = shift;
  
  if ($yscale < 2*$height - 1) {
    $height = 0;
  }
  
  #########
  # chage scale later so that the base can  be anything and not just 1!!
  
  # scale legend goes in order from min, axis, max, if either min or max = 1, then will only have min & max
  my @order = sort {$a <=> $b} (1, @_);
  my $graph_scale = - ($yscale - $height) / (log10($order[2]) - log10($order[0]));
  my $graph_legend = {1 => $graph_scale * (log10(1) - log10($order[2])),
  		   $order[0] => $graph_scale * (log10($order[0]) - log10($order[2])),
    		   $order[2] => 0};
    
  #print "order is @order and the yscale is $yscale and height is $height<br>";
  
  return ($graph_legend, $graph_scale);
}

#main method that draws everything
sub draw {
  my $self = shift;
  my $height = $self->font->height;
  my $scale = $self->scale;
  
  my $gd = shift;
  my ($left,$top,$partno,$total_parts) = @_;
  my ($x1,$y1,$x2,$y2) = $self->bounds($left, $top);
  
  
  my @bounds = $gd->getBounds;
  my $draw_clado_left = $self->option('draw_clado_left');
  
  
  #spacing of either DNA alignments or score histograms in units of font height
  my $species_spacing = $self->option('species_spacing') || 1;
  
  my $xscale = $self->font->width;
  my $yscale = $height * $species_spacing;
  
  
  my $xoffset = $x1;
  my $yoffset = $y1 + 0.5*$self->font->height;
  
  #method that reads the tree file to create the tree objects
  my $tree = $self->set_tree;
  my $max_height = get_max_height($tree);
  my $start_x =($max_height-1) * $xscale +$xoffset;
  
  
  my $connector = $self->connector;
  
  #all species having alignments in viewing window (key=name, val=feat obj)
  my %alignments = $self->extract_features;
  #print "Species are:",keys %alignments,"<br>\n";
  
  my ($min_score, $max_score) = $self->get_score_bounds(%alignments);
  #$min_score = 0 unless $min_score;
  my ($graph_legend, $graph_scale) = get_legend_and_scale($yscale, $height, $min_score, $max_score);
#print "min/max scores: $min_score, $max_score<br>\n",
#"graph legend and scale: $graph_legend, $graph_scale";
# TODO: Gap entries give an undef for the min values for some reason
  
  
  my $refspecies = $self->option('reference_species');
  
  my @current_species = keys %alignments;    #all species in viewing window
  my @known_species = $self->known_species($tree);  #all species from cladogram info
  my @unknown_species = $self->unknown_species(\%alignments, 
    						 $refspecies,
    						\@current_species,
    						\@known_species);
                                              #species in GFF but not in clado

  ########is this even used?
  #my @allfeats;
  #for my $species (keys %alignments) {
  #  push @allfeats, @{$alignments{$species}};
  #}
  
  
  #this y value is the base for the next species' alignment/histogram and is incremented at each step 
  my $y = $y1;
  
  
  for my $species (@known_species,@unknown_species) {
    my $y_track_top = $y + $height;
    my $y_track_bottom = $y + $yscale;
    
    
    if ($yscale < 2*$height-1) {
      #small scale
      $y_track_top = $y;# + $height;
      my $y_track_bottom = $y + $height;
    }
    
    
    #process the reference sequence differently
    if ($species eq $refspecies) {
      #draw DNA alignments if zoomed close enough
      my ($fx1,$fy1) = ($x1, $y_track_top);
      my ($fx2,$fy2) = ($x2,$y_track_bottom);
      

      if ($self->dna_fits) {
	my $dna = eval { $self->feature->seq };
        $dna    = $dna->seq if ref($dna) and $dna->can('seq'); # to catch Bio::PrimarySeqI objects
        my $bg_color = $self->color('ref_color') || $self->bgcolor;
        
        $fy2 = $fy1 + $self->font->height || $y2;
  
        
	$self->_draw_dna($gd,$dna,$fx1,$fy1,$fx2,$fy2, $self->fgcolor, $bg_color);
      } else {
      }
      
      my $x_label_start = $start_x + $xoffset + $self->font->width;
      $self->species_label($gd, $draw_clado_left, $x_label_start, $y, $species) unless ($self->option('hide_label'));
      
      $y += $yscale;
      next;
    }
    
    
    #skip if the there is no alignments for this species in this window
    unless ($alignments{$species}) {
      my $x_label_start = $start_x + $xoffset + $self->font->width;
      $self->species_label($gd, $draw_clado_left, $x_label_start, $y, $species) unless ($self->option('hide_label'));
      
      $y += $yscale;
      next;
    }
    
    
    
    my @features = @{$alignments{$species}};
    
    
    
    
    #draw the axis for the plots
    $self->draw_pairwisegraph_axis($gd,
    				    $graph_legend,
    				    $x1,
    				    $x2,
    				    $y_track_top,
    				    $y_track_bottom,
    				    $draw_clado_left,
    				    @bounds) unless $self->dna_fits;
      
    
    #iterate through the wigfiles and put them on the graph
    ###
    # todo
    ###
    
    
    #iterate through features, and put them on the graph
    for my $feat (@features) {
      my ($start, $stop, %attributes) = ($feat->start, $feat->stop, $feat->attributes);
      
      my ($fx1,$fy1) = ($x1 + ($start-$self->start)*$scale, $y_track_top);
      my ($fx2,$fy2) = ($x1 + ($stop-$self->start)*$scale,$y_track_bottom);
      
      my $gapstr = $attributes{'Gap'} || "M".($stop-$start+1);
      my @gapstr = split " ", $gapstr;
      my @gaps;
      for my $gap (@gapstr) {
        my ($type, $num) = $gap =~ /^(.)(\d+)/; 
      	push @gaps, [$type, $num+0];
      }
      
      
      
      #draw DNA alignments if zoomed close enough
      if ($self->dna_fits) {
	
	
	my $test = 0;
	
	my $hit = $feat->hit;
	if (!defined $hit) {
	  warn "No hit for feature $feat, skipping drawing DNA";
	  next;
	}

	my $hit_seq = $hit->seq || print "No seq for hit<br>\n";
	my $targ_dna = $hit_seq->seq || print "No seq for hit_seq<br>\n";
	
	
	my $seq = $feat->seq || print "No sequence object for feature<br>\n";
	my $ref_dna = $seq->seq || print "No ref dna";
	
	#doesn't work as planned
	# my $ref_dna = $feat->seq->seq || print "No ref dna";
	# my $targ_dna = $feat->hit->seq->seq || print "No targ dna";
	
	
	next if !defined $ref_dna || !defined $targ_dna;
	
	$self->draw_dna($gd,$ref_dna, $targ_dna,$fx1,$fy1,$fx2,$fy2,\@gaps);
      } else {
      	my $wigfile = $attributes{'wigfile'};
      	if ($wigfile) {
      	  if (-e $wigfile) {
      	    $self->pairwise_draw_wig_graph($gd, $feat, $x1, $scale, \@gaps, $graph_legend->{1}, $graph_scale, $fx1, $fy1, $fx2, $fy2,$wigfile);
      	  } else {
      	    warn "Wigfile $wigfile does not exist, skipping ...";
      	  }
      	  
      	} else {
      	  $self->pairwise_draw_graph($gd, $feat, $x1, $scale, \@gaps, $graph_legend->{1}, $graph_scale, $fx1, $fy1, $fx2, $fy2);
      	}
      	
      }
    }
    
    
    #label the species in the cladogram
    my $x_label_start = $start_x + $xoffset + $self->font->width;
    $self->species_label($gd, $draw_clado_left, $x_label_start, $y, $species) unless ($self->option('hide_label'));
    
    $y += $yscale;
  }
  
    $self->draw_clado($tree, $gd, $x1, $y1, $x2, $y2, $self->fgcolor,
  		     $xscale, $yscale, $xoffset, $yoffset, $start_x, $draw_clado_left);

}



sub species_label {
  my $self = shift;
  my $gd = shift;
  my $draw_clado_left = shift;
  my $x_start = shift;
  my $y_start = shift;
  my $species = shift;
  
  $x_start += 2;
  my $text_width = $self->font->width * length($species);
  my $bgcolor = $self->color('bg_color');
  
  #make label
  if ($draw_clado_left) {
    
    $gd->filledRectangle($x_start-2, $y_start, $x_start + $text_width, $y_start+$self->font->height, $bgcolor);
    $gd->rectangle($x_start-2, $y_start, $x_start + $text_width, $y_start+$self->font->height, $self->fgcolor);
    $gd->string($self->font, $x_start, $y_start, $species, $self->fgcolor);
    
  } else {
    my ($x_max, $y_max) = $gd->getBounds;
    my $write_pos = $x_max - $x_start - $text_width;
    
    $gd->filledRectangle($write_pos, $y_start, $write_pos + $text_width+2, $y_start+$self->font->height, $bgcolor);
    $gd->rectangle($write_pos, $y_start, $write_pos + $text_width+2, $y_start+$self->font->height, $self->fgcolor);
    $gd->string($self->font, $write_pos+2, $y_start, $species, $self->fgcolor);
    
  }
}


# draws the legends on the conservation scale
sub draw_pairwisegraph_axis {
  my $self = shift;
  my ($gd, $graph_legend, $x1, $x2, $y_track_top, $y_track_bottom, $draw_clado_left, @bounds) = @_;
  
  
  my $axis_color = $self->color('axis_color') || $self->fgcolor;
  my $mid_axis_color = $self->color('mid_axis_color') || $axis_color;
  
  for my $label (keys %$graph_legend) {
    my $y_label = $graph_legend->{$label} + $y_track_top;

    
    my $col = $axis_color;
    $col = $mid_axis_color if ($y_label != $y_track_top && $y_label != $y_track_bottom);
    $gd->line($x1,$y_label,$x2,$y_label,$col);
    
    my @coords = (0, $y_label, $x1, $y_label);
    
    
    if ($draw_clado_left) {
      #draw the legend on the right
      $coords[0] = $bounds[0] - $coords[0];
      $coords[2] = $bounds[0] - $coords[2];
      
      my $x_text_offset = length($label) * $self->font->width;
      
      $gd->string($self->font, $coords[0]-$x_text_offset, $coords[1], $label, $self->fgcolor);
      $gd->line(@coords, $self->fgcolor);
      
      $gd->line($x2,$y_track_top,$x2,$y_track_bottom,$self->fgcolor);
    } else {
      #draw the legned on the left
      $gd->string($self->font, @coords[0..1], $label, $self->fgcolor);
      $gd->line(@coords, $self->fgcolor);
      
      $gd->line($x1,$y_track_top,$x1,$y_track_bottom,$self->fgcolor);
    }
  
  }
  
}


#find min and max from features within the bounds
sub get_score_bounds {
  my $self = shift;
  my %alignments = @_;
  
  my $min = -1;
  my $max = -1;
  
  for my $species (keys %alignments) {
    for my $feature (@{$alignments{$species}}) {
      my $score = $feature->score;
      
      #check to see if wigfile
      if ($score == undef) {
      	my %attributes = $feature->attributes;
      	if (-e $attributes{'wigfile'}) {
      	  ($min, $max) = $self->get_score_bounds_wigfile($feature,$min,$max,$attributes{'wigfile'});
      	}
      	next;
      }
      
      $min = $score if $min == -1 || $score < $min;
      $max = $score if $max == -1 || $max < $score;
    }
  }
  
  
  my @parts = $self->parts;
  
  return ($min, $max)
}

#find min and max of sampled wigfile
sub get_score_bounds_wigfile {
  my $self = shift;
  my $feature = shift;
  my ($min,$max,$wigfile) = @_;
  
  my $start = $feature->start < $self->start ? $self->start : $feature->start;
  my $stop  = $self->stop < $feature->stop   ? $self->stop  : $feature->stop;
  
  #print $self->stop, "-", $feature->stop, "checking @_ $start - $stop\n";
  
  #extract wig file contents  
  my $wig = Bio::Graphics::Wiggle->new($wigfile, 0, {step => 1}) or die;
  
  #todo: make step configurable
  my $step = 100;
  
  for (my $i=$start; $i<$stop; $i += $step) {
    my $v = $wig->value($i);
    
    next unless defined $v;
    
    $min = $v if !defined $min or $v < $min;
    $max = $v if !defined $max or $max < $v;
    
  }
  
  #print "min and max are $min , $max<br>\n";
  
  return ($min,$max);
  
}



sub pairwise_draw_graph {
  my $self = shift;
  my $gd = shift;
  my $feat = shift;		# current feature object
  my $x_edge = shift;		# x start position of the track
  my $scale = shift;		# pixels / bp
  my $gaps = shift;		# gap data for insertions, deletions and matches
  my $zero_y = shift;		# y coordinate of 0 position
  my $graph_scale = shift;	# scale for the graph. y_coord = graph_scale x log(score)
  
  my ($x1,$y1,$x2,$y2) = @_;
  my $fgcolor = $self->fgcolor;
  my $errcolor  = $self->color('errcolor') || $fgcolor;
  
  my $score = $feat->score;
  my %attributes = $feat->attributes;
  
  
  my $log_y = log10($score);
  my $y_bottom = log10($score) * $graph_scale + $zero_y + $y1;
  my $y_top = $zero_y+$y1;
  
  my @y = sort {$a <=> $b} ($y_bottom, $y_top);
  
  
  
  #missing gap data
  unless ($gaps) {
    $x1 = $x_edge if $x1 < $x_edge;
    return if $x2 < $x_edge;
    #$gd->filledRectangle($x1,$y[0],$x2,$y[1],$fgcolor);
    return;
  }
  
  my $bp = 0;
  
  #draw a bar representing the score for the span of base pairs
  for my $tuple (@$gaps) {
    my ($type, $num) = @$tuple;
    
    #warn"$type and $num";
    if ($type eq "M") {
      my $x_left  = $x1 + ($bp*$scale);
      my $x_right = $x_left + $num*$scale;
      
      $bp += $num;
      
      $x_left = $x_edge if $x_left < $x_edge;
      next if $x_right < $x_edge;
      $gd->filledRectangle($x_left,$y[0],$x_right,$y[1],$fgcolor);
    } elsif ($type eq "D") {
      
      $bp += $num;
      
    } elsif ($type eq "I") {
      my $x_left  = $x1 + ($bp*$scale);
      $gd->line($x_left-2, $y1-4, $x_left, $y1, $errcolor);
      $gd->line($x_left, $y1, $x_left+2, $y1-4, $errcolor);
    }
    

  }
  
  
}



sub pairwise_draw_wig_graph {
  my $self = shift;
  my $gd = shift;
  my $feat = shift;		# current feature object
  my $x_edge = shift;		# x start position of the track
  my $scale = shift;		# pixels / bp
  my $gaps = shift;		# gap data for insertions, deletions and matches
  my $zero_y = shift;		# y coordinate of 0 position
  my $graph_scale = shift;	# scale for the graph. y_coord = graph_scale x log(score)
  
  my ($x1,$y1,$x2,$y2,$wigfile) = @_;
  my $fgcolor = $self->fgcolor;
  my $errcolor  = $self->color('errcolor') || $fgcolor;
  
  my $score = $feat->score;
  my %attributes = $feat->attributes;
  
  
#  my $log_y = log10($score);
#  my $y_bottom = log10($score) * $graph_scale + $zero_y + $y1;
#  my $y_top = $zero_y+$y1;
#  
#  my @y = sort {$a <=> $b} ($y_bottom, $y_top);
  
  
  
  #print "checking wigfile $wigfile<br>\n";
  
  
  #todo: make step variable
  my $start = $feat->start < $self->start ? $self->start : $feat->start;
  my $stop  = $self->stop < $feat->stop   ? $self->stop  : $feat->stop;
#  my $wig = Bio::Graphics::Wiggle->new($wigfile, 0) or die;
  my $wig = Bio::Graphics::Wiggle->new($wigfile, 0, {step => 1}) or die;


##### not sure why this was here
#  my $vals = $wig->values($start,$stop);
#  my ($vmin,$vmax);
#  for my $v (@$vals) {
#    next unless defined $v;
#    $vmin = $v if !defined $vmin or $v < $vmin;
#    $vmax = $v if !defined $vmax or $vmax < $v;
#  }
  
  
#  print "min and max are $vmin , $vmax\n";
#  $min = $min < $vmin ? $min : $vmin;
#  $max = $vmax < $max ? $max : $vmax;
#  print "min and max are $min , $max\n";
#  return ($min,$max);
  
  
#  print Dumper($vals);
  
  
#  my $pos = $start;
  
  
#  $gd->rectangle($x1,$y1,$x2,$y2,$fgcolor);
#  $gd->line($x1,$y1,$x2,$y2,$fgcolor);
#  $gd->line($x2,$y1,$x1,$y2,$fgcolor);
#  print "Start and stop: $start and $stop<br>\n";
#  print "<br>\n$x1, $x2, $x_edge, $start, $stop<br>\n";
  
  
  
  if ($scale < 1) {
    #### Algorithm 1, when zoomed at scale 1bp / pixel or more
    #### sample 10 values across the pixel (prevents redrawing same pixel over
    #### and over).  The sample values are averaged.
    
    my $bp = $start;
    for (my $pix=$x1; $pix < $x2; $pix++) {
      
      $bp = ($pix - $x1)/$scale + $start;
      
      #todo: make samplesize an option
      my $samplesize = 10;
      my $score = 0;
      my $trial;
      for ($trial=0; $trial < $samplesize; $trial++) {
        
        my $samp_bp = $bp + ($trial/$samplesize)/$scale;
        last if $samp_bp > $stop;
  
        $score += $wig->value($samp_bp);
        
        #print "trial $trial: $samp_bp\tPixel $pix\tpos $x1<br>\n" if ($bp > 1700 && $bp < 1750);
        
        
      }
      
      #print "Pixel: $pix : Total $score and trial $trial<br>\n";
      next if $trial == 0 || $score == 0;
      $score = $score / $trial;
      
      $self->draw_log10_rectangle($score, $graph_scale, $zero_y, $y1,
      				$zero_y, $pix, $pix, $gd, $fgcolor);
      
      
    }
    
  } else {
    #### Algorithm 2, when zoomed in close at less than 1bp / pixel (1 bp spans
    #### 1 pixel or more), draw each value directly
    
    my $step = 1;
    
    for (my $i=$start; $i<$stop; $i += $step) {
      my $val = $wig->value($i);
      next if !defined $val;
      
      my $bp = $i - $start;
      
      my $x_left  = $x1 + ($bp*$scale);
      my $x_right = $x_left + 1*$scale;
      
      $score = $val;
      
      
      $self->draw_log10_rectangle($score, $graph_scale, $zero_y, $y1, 
      				$zero_y, $x_left, $x_right, $gd, $fgcolor);
      
      
    }
    
    
  }
  
  return;

  
}

sub draw_log10_rectangle {
  my $self = shift;
  my $score = shift;
  my $graph_scale = shift;
  my $zero_y = shift;
  my $y1 = shift;
  my $zero_y = shift;
  my $x_left = shift;
  my $x_right = shift;
  my $gd = shift;
  my $fgcolor = shift;
  
  my $log_y = log10($score);
  my $y_bottom = log10($score) * $graph_scale + $zero_y + $y1;
  my $y_top = $zero_y+$y1;
  
  my @y = sort {$a <=> $b} ($y_bottom, $y_top);
  
  #print "$x_left,$y[0],$x_right,$y[1]<br>\n";
  
  $gd->filledRectangle($x_left,$y[0],$x_right,$y[1],$fgcolor);
}



sub draw_dna {
  my $self = shift;
  my ($gd,$ref_dna, $dna,$x1,$y1,$x2,$y2, $gaps) = @_;
  my $pixels_per_base = $self->scale;
  
  my $fgcolor = $self->fgcolor;
  my $bg_color = $self->color('targ_color') || $self->bgcolor;
  my $errcolor  = $self->color('errcolor') || $fgcolor;
  
  $y2 = $y1 + $self->font->height || $y2;
  
  
  
  #missing gap data, draw as is
  unless ($gaps) {
    warn"no gap data for DNA sequence $dna";
    $self->_draw_dna($gd, $dna, $x1, $y1, $x2, $y2);
    return;
  }
  
  #parse the DNA segments by the gaps
  for my $tuple (@$gaps) {
    my ($type, $num) = @$tuple;
    
    if ($type eq "M") {
      my $dnaseg = substr($dna, 0, $num);
      my $ref_dnaseg = substr($ref_dna, 0, $num);
      
      $self->_draw_dna($gd,$dnaseg, $x1, $y1, $x2, $y2, $fgcolor, $bg_color,$ref_dnaseg);
      $dna = substr($dna, $num);
      $ref_dna = substr($ref_dna, $num);
      $x1 += $num * $pixels_per_base;
    } elsif ($type eq "D") {
      my $dnaseg = '-' x $num;
      $self->_draw_dna($gd, $dnaseg, $x1, $y1, $x2, $y2, $fgcolor); 
      $gd->rectangle($x1, $y1-1, $x1+$num * $pixels_per_base, $y2+1, $errcolor);
      $ref_dna = substr($ref_dna, $num);
      $x1 += $num * $pixels_per_base;
    } elsif ($type eq "I") {
      $dna = substr($dna, $num);
      $gd->line($x1-2, $y1-2, $x1+2, $y1-2, $errcolor);
      $gd->line($x1, $y1-1, $x1, $y2+1, $errcolor);
      $gd->line($x1-2, $y2+1, $x1+2, $y2+1, $errcolor);
    }
    

  }
  
  
}

sub _draw_dna {
  my $self = shift;

  #the last argument is optional.  If the reference seq is given, it will check it with target
  my ($gd,$dna,$x1,$y1,$x2,$y2, $color, $bg_color, $ref_dna) = @_;
  
  my $pixels_per_base = $self->scale;
  my $feature = $self->feature;
  
  unless ($ref_dna) {
    $gd->filledRectangle($x1+1, $y1, $x2, $y2, $bg_color);
  }
  
  
  my $feature = $self->feature;
  
  
  my $strand = $feature->strand || 1;
  $strand *= -1 if $self->{flip};

  my @bases = split '',$strand >= 0 ? $dna : $self->reversec($dna);
  my @refbases = split '',$strand >= 0 ? $ref_dna : $self->reversec($ref_dna);
  
  
  
  $color = $self->fgcolor unless $color;
  $bg_color = 0 unless $bg_color;
  my $font  = $self->font;
  my $lineheight = $font->height;
#  $y1 -= $lineheight/2 - 3;          ##################NOT SURE WHY THIS WAS HERE BEFORE
  my $strands = $self->option('strand') || 'auto';

  my ($forward,$reverse);
  if ($strands eq 'auto') {
    $forward = $feature->strand >= 0;
    $reverse = $feature->strand <= 0;
  } elsif ($strands eq 'both') {
    $forward = $reverse = 1;
  } elsif ($strands eq 'reverse') {
    $reverse = 1;
  } else {
    $forward = 1;
  }
  # minus strand features align right, not left
  $x1 += $pixels_per_base - $font->width - 1 if $strand < 0;
  for (my $i=0;$i<@bases;$i++) {
    my $x = $x1 + $i * $pixels_per_base;
    
    my $x_next = $x + $pixels_per_base;

    
    #draw background if DNA base aligns with reference (if ref given)
    $gd->filledRectangle($x+1, $y1, $x_next, $y2, $bg_color) 
    			if ( ($forward && $bases[$i] eq $refbases[$i]) ||
    			     ($reverse && $complement{$bases[$i]} eq $refbases[$i]) );
    
    $gd->char($font,$x+2,$y1,$bases[$i],$color)                                   if $forward;
    $gd->char($font,$x+2,$y1+($forward ? $lineheight:0),
	      $complement{$bases[$i]}||$bases[$i],$color)                         if $reverse;
  }
}


1;

__END__