view release on metacpan or  search on metacpan

lib/Bio/Graphics/Glyph/phylo_align.pm  view on Meta::CPAN

    } 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 $font = $self->mono_font;
  my $height = $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 = $font->width;
  my $yscale = $height * $species_spacing;
  
  
  my $xoffset = $x1;
  my $yoffset = $y1 + 0.5*$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 + $font->height || $y2;
  
        
	$self->_draw_dna($gd,$dna,$fx1,$fy1,$fx2,$fy2, $self->fgcolor, $bg_color);
      } else {
      }
      
      my $x_label_start = $start_x + $xoffset + $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 + $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 + $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;
  
  my $font = $self->mono_font;

  $x_start += 2;
  my $text_width = $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+$font->height, $bgcolor);
    $gd->rectangle($x_start-2, $y_start, $x_start + $text_width, $y_start+$font->height, $self->fgcolor);
    $gd->string($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+$font->height, $bgcolor);
    $gd->rectangle($write_pos, $y_start, $write_pos + $text_width+2, $y_start+$font->height, $self->fgcolor);
    $gd->string($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 $font = $self->mono_font;
  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) * $font->width;
      
      $gd->string($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($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) = @_;