binaries/src/ViennaRNA/Utils/mountain.pl

   1 #!/usr/bin/perl -w
   2 # -*-Perl-*-
   3 # Last changed Time-stamp: <2008-08-26 16:04:00 ivo>
   4 # produce Pauline Hogeweg's mountain representation *_dp.ps files
   5 # writes 3 sets of x y data separated by a "&"
   6 # the first two sets are mountain representations from base pair probabilities
   7 # and mfe structure, respectively.
   8 # definition: mm[i],mp[i] = (mean) number of base pairs enclosing base i
   9 # third set a measure of well-definedness: the entropy of the pair probs of
  10 # base i, sp[i] = -Sum p_i * ln(p_i). Well-defined regions have low entropy.
  11 #
  12 # use e.g. as  mountain.pl dot.ps | xmgrace -pipe
  13
  14 use strict;
  15 our (@mm, @mp, @pp, @sp, @p0, $i, $max, $length, $do_png);  # perl5 only
  16
  17 my $sep = "&";   # xmgr uses & to separate data sets
  18
  19 if (@ARGV && ($ARGV[0] eq '-png')) {
  20     eval "use Chart::Lines";
  21     die($@,
  22         "\nCould not load the Chart::Lines module required with -png option\n")
  23         if $@;
  24     $do_png=1;
  25     shift;
  26 }
  27
  28
  29 while (<>) {
  30     chomp;
  31     if (/\/sequence \{ \((\S*)[\\\)]/) {
  32         my $seq = $1;           # empty for new version
  33         while (!/\) \} def/) {  # read until end of definition
  34             $_ = <>;
  35             /(\S*)[\\\)]/;      # ends either with `)' or `\'
  36             $seq .= $1;
  37         }
  38         $length = length($seq);
  39         next;
  40     }
  41
  42     next unless /(\d+) (\d+) (\d+\.\d+) (.box)$/;
  43     my ($i, $j, $p, $id) = ($1,$2,$3,$4);
  44     if ($id eq "ubox") {
  45         $p *= $p;           # square it to probability
  46         $mp[$i+1] += $p;
  47         $mp[$j]   -= $p;
  48         my $ss = $p>0 ? $p*log($p) : 0;
  49         $sp[$i] += $ss;
  50         $sp[$j] += $ss;
  51         $pp[$i] += $p;
  52         $pp[$j] += $p;
  53     }
  54     if ($id eq "lbox") {
  55         $mm[$i+1]++;
  56         $mm[$j]--;
  57     }
  58 }
  59 $mp[0] = $mm[0] = $max = 0;
  60 for ($i=1; $i<=$length; $i++) {
  61     no warnings;
  62     $mp[$i]+=$mp[$i-1];
  63     $max = $mp[$i] if ($mp[$i]>$max);
  64     $mm[$i]+=$mm[$i-1];
  65     $max = $mp[$i] if ($mp[$i]>$max);
  66     $sp[$i] += (1-$pp[$i])*log(1-$pp[$i]);
  67 }
  68
  69 if ($do_png) {
  70     my $width =  800;
  71     my $height = 600;
  72
  73     # FIXME: legend_lables when doing mfe only
  74     my $skip = 10**(int (log($length)/log(10.) - 0.5));
  75     my $obj = Chart::Lines->new( $width, $height );
  76     $obj->set ('title' => $ARGV,
  77                'x_label' => 'Position',
  78                'y_label' => 'Height',
  79                'min_val' => 0,
  80                'precision' => 0,
  81                'legend_labels' => ['mfe', 'pf'],
  82                'skip_x_ticks' => $skip);
  83
  84     $obj->add_dataset ((0..$length));
  85
  86     $obj->add_dataset (@mp);
  87     $obj->add_dataset (@mm);
  88     $obj->png("mountain.png");
  89
  90 } else {
  91     # print the results for plotting
  92     for ($i=1; $i<=$length; $i++) {
  93         printf("%4d  %7.5g\n", $i, $mp[$i]);
  94     }
  95     print "$sep\n";
  96
  97     for ($i=1; $i<=$length; $i++) {
  98         printf("%4d  %4d\n", $i, $mm[$i]);
  99     }
 100     print "&\n";
 101     my $log2 = log(2);
 102     for ($i=1; $i<=$length; $i++) {
 103         printf("%4d  %7.5g\n", $i, -$sp[$i]/$log2);
 104     }
 105 }
 106
 107 =head1 NAME
 108
 109     mountain - produce coordinates for a mountain plot from a dot plot
 110
 111 =head1 SYNOPSIS
 112
 113     mountain.pl myseq_dp.ps > mountain.dat
 114
 115 =head1 DESCRIPTION
 116
 117     reads pair proabilities and MFE structure from a probability dot
 118     plot as produced by C<RNAfold -p>, and produces x-y data suitable
 119     for producing a mountain plot using standard xy-plotting programs.
 120
 121     Output consists of 3 data sets separated by a line containing only
 122     the C<&> character. The first two sets are mountain representations
 123     computed from base pair probabilities and mfe structure, respectively.
 124     For the mfe case the moutain plot graphs the number base pairs
 125     enclosing a position k, in case of pair probabilities we use the average
 126     number of base pairs computed as m_k = \Sum_i<k<j p_ij.
 127     The third set contains the positional entropy, which provides a measure
 128     of local structural welldefinedness, s_i = -\Sum_j p_ij * ln(p_ij).
 129
 130     The output is suitable for graphing with xmgrace, e.g.:
 131     C< RNAfold -p < foo.seq; mountain.pl foo_dp.ps | xmgrace -pipe>
 132
 133 =head1 AUTHOR
 134
 135 Ivo L. Hofacker <ivo@tbi.univie.ac.at>