view release on metacpan or  search on metacpan

lib/Bio/Graphics/Pictogram.pm  view on Meta::CPAN

    my @seq;
    while (my $seq = $sio->next_seq){
      push @seq, $seq;
    }
    @pwm = @{$self->_make_pwm(\@seq)};
  }


  my $svg = $self->svg_obj;
  my $seq_length = scalar(@pwm + 1) * $width + $x + $x;
  my $seq_grp;

  #scale the svg if length greater than svg width
  if($seq_length > $svg->{-document}->{'width'}){
    my $ratio = $svg->{-document}->{'width'}/($seq_length);
    $seq_grp = $svg->group(transform=>"scale($ratio,1)");
  }
  else {
    $seq_grp= $svg->group();
  }

  #do the drawing, each set is a base position
  foreach my $set(@pwm){
    my ($A,$C,$G,$T,$bits) = @$set;
    my @array;
    push @array,  ['a',($A)];
    push @array, ['g',($G)];
    push @array, ['c',($C)];
    push @array, ['t',($T)];
    @array = sort {$b->[1]<=>$a->[1]}@array;
    my $count = 1;
    my $pos_group = $seq_grp->group(id=>"bp $bp");
    my $prev_size;
    my $y_trans;

    #draw each letter at each position
    foreach my $letter(@array){
	  my $scale;
	  if($self->normalize){
		$scale = $letter->[1];
	  } else {
		$scale = $letter->[1] * ($bits / MAXBITS);
	  }

      if($count == 1){
		if($self->normalize){
		  $y_trans = 0;
		} else {
		  $y_trans = (1 - ($bits / MAXBITS)) * $size;
		}
      }
      else {
        $y_trans += $prev_size;
      }
      $pos_group->text('id'=> uc($letter->[0]).$bp,height=>$height,
                      'width'=>$width,x=>$x,y=>$size,
                      'transform'=>"translate(0,$y_trans),scale(1,$scale)",
                      'style'=>{"font-size"=>$fontsize,
                      'text-anchor'=>'middle',
                      'font-family'=>'Verdana',
                      'fill'=>$color->{uc $letter->[0]}})->cdata(uc $letter->[0]) if $scale > 0;

     $prev_size = $scale * $size;
     $count++;
    }
    #plot the bit if required
    if($self->plot_bits){
         $seq_grp->text('x'=>$x,
                        'y'=>$bit_y,
                        'style'=>{"font-size"=>'10',
                                'text-anchor'=>'middle',
                                'font-family'=>'Verdana',
                                'fill'=>'black'})->cdata($bits);
    }
    $bp++;
    $x+=$width;
  }

  #plot the tags
  $seq_grp->text(x=>int($width/2),y=>$bit_y,style=>{"font-size"=>'10','text-anchor'=>'middle','font-family'=>'Verdana','fill'=>'black'})->cdata("Bits:") if $self->plot_bits;

 $seq_grp->text(x=>int($width/2),y=>$pos_y,style=>{"font-size"=>'10','text-anchor'=>'middle','font-family'=>'Verdana','fill'=>'black'})->cdata("Position:");

  #plot the base positions
  $x = 45+$size/2-int($width/2);
  foreach my $nbr(1..($bp-1)){
    $seq_grp->text(x=>$x+int($width/2),y=>$pos_y,style=>{"font-size"=>'10','text-anchor'=>'left','font-family'=>'Verdana','fill'=>'black'})->cdata($nbr);
    $x+=$width;
  }


#  $seq_grp->transform("scale(2,2)");

  return $self->svg_obj($svg);
}

sub _make_pwm_from_site_matrix{
  my ($self,$matrix) = @_;
  my $bgd = $self->background;
  my @pwm;
  my $consensus = $matrix->consensus;
  foreach my $i(1..length($consensus)){
    my %base = $matrix->next_pos;
    my $bits;
    $bits+=($base{pA} * log2($base{pA}/$bgd->{'A'}));
    $bits+=($base{pC} * log2($base{pC}/$bgd->{'C'}));
    $bits+=($base{pG} * log2($base{pG}/$bgd->{'G'}));
    $bits+=($base{pT} * log2($base{pT}/$bgd->{'T'}));
    push @pwm, [$base{pA},$base{pC},$base{pG},$base{pT},abs(sprintf("%.3f",$bits))];
  }
  return @pwm;
}

sub _make_pwm {
  my ($self,$input) = @_;
  my $count = 1;
  my %hash;
  my $bgd = $self->background;
  #sum up the frequencies at each base pair
  foreach my $seq(@$input){
    my $string = $seq->seq;
    $string =  uc $string;
    my @motif = split('',$string);
    my $pos = 1;
    foreach my $t(@motif){
      $hash{$pos}{$t}++;
      $pos++;
    }
    $count++;
  }

  #calculate relative freq
  my @pwm;

  #decrement last count
  $count--;
  foreach my $pos(sort{$a<=>$b} keys %hash){
    my @array;
    push @array,($hash{$pos}{'A'}||0)/$count;
    push @array,($hash{$pos}{'C'}||0)/$count;
    push @array,($hash{$pos}{'G'}||0)/$count;
    push @array,($hash{$pos}{'T'}||0)/$count;

    #calculate bits
    # relative entropy (RelEnt) or Kullback-Liebler distance
    # relent = sum fk * log2(fk/gk) where fk is frequency of nucleotide k and
    # gk the background frequency of nucleotide k