binaries/src/ViennaRNA/Utils/relplot.pl

   1 #!/usr/bin/perl -w
   2 # -*-CPerl-*-
   3 # Last changed Time-stamp: <2012-11-19 17:22:03 ivo>
   4 # $Id: relplot.pl,v 1.10 2008/10/09 07:11:21 ivo Exp $
   5 # colorize a secondary structure plot with reliability annotation
   6 # from positional entropy
   7 use strict;
   8 use Getopt::Std;
   9 $main::VERSION = 1.3;
  10 $Getopt::Std::STANDARD_HELP_VERSION=1;
  11
  12 our ($opt_p, $opt_a);
  13 getopts('pa');
  14
  15 sub HELP_MESSAGE {
  16   print STDERR "\nusage: $0 [-p] FOO_ss.ps FOO_dp.ps > FOO_rss.ps\n";
  17   print STDERR "For more details run\n\tperldoc -F $0\n";
  18 }
  19
  20 HELP_MESSAGE() unless $#ARGV >0;
  21 my $macro_seen= 0;
  22 my %mfe = ();        # hash of mfe pairs
  23 my @ss_ps = ('',''); # head and tail of the ss.ps file
  24
  25 my $n = swallow_ss_ps();    # read ss plot
  26 my @sp = posent();   # read dot plot and compute entropies
  27 my $Smax = ($opt_p || $opt_a) ? 1 : 0;
  28 if (!$opt_p) {
  29   foreach (@sp) {
  30     $Smax = $_ if $_>$Smax;
  31   }
  32   $Smax = ($Smax>0.2) ? sprintf("%3.1f", $Smax) : 0.2 ;
  33
  34 }
  35 print $ss_ps[0];     # print head
  36 if (!$macro_seen) {
  37   print <<_E_O_F_
  38 /range 0.8 def
  39 /drawreliability {
  40   /Smax $Smax def
  41   0
  42   coor {
  43     aload pop
  44     S 3 index get
  45     Smax div range mul
  46     invert {range exch sub} if
  47     1 1 sethsbcolor
  48     newpath
  49     fsize 2 div 0 360 arc
  50     fill
  51     1 add
  52   } forall
  53 } bind def
  54 /colorbar { % xloc yloc colorbar -> []
  55   /STR 8 string def
  56   gsave
  57     xmin xmax add size sub 2 div
  58     ymin ymax add size sub 2 div translate
  59     size dup scale
  60     translate
  61     0.015 dup scale
  62     /tics 64 def
  63     gsave
  64       10 tics div 1 scale
  65       0 1 tics
  66       {
  67           dup 0 moveto 0.5 add
  68           tics div range mul
  69           invert {range exch sub} if
  70           1 1 sethsbcolor
  71           1 0 rlineto 0 1 rlineto -1 0 rlineto closepath fill
  72       } for
  73     grestore
  74     0 setgray
  75     -0.1 1.01 moveto (0) gsave 0.1 dup scale show grestore
  76     10 1.01 moveto Smax STR cvs
  77     gsave 0.1 dup scale dup stringwidth pop -2 div 0 rmoveto show grestore
  78   grestore
  79 } bind def
  80 _E_O_F_
  81 }
  82 print "/S [\n";
  83 foreach (@sp) {
  84   printf "  %7.5f\n", $_;
  85 }
  86 print "] def\n\n";
  87 print "/invert ", $opt_p||$opt_a ? 'true' : 'false', " def\n";
  88 print "drawreliability\n";
  89 print "0.1 0.1 colorbar\n";
  90 print $ss_ps[1];  # print tail
  91
  92 sub swallow_ss_ps {
  93   # read the secondary structure plot
  94   my $length=0;
  95   my $tail=0;
  96   while (<>) {
  97     $macro_seen=1 if /drawreliability /;
  98     $length ++ if /^\/coor/ .. /^\] def/;
  99     if (/^\/pairs/ .. /^\] def/) {
 100       $mfe{$1,$2}=1 if /(\d+)\s+(\d+)/;
 101     }
 102     $tail=1 if /^drawoutline/;
 103     $ss_ps[$tail] .= $_;
 104     last if eof;
 105   }
 106   return $length-2;
 107 }
 108
 109 sub posent {
 110   # compute positional entropy from pair probs in the dot plot file
 111   # or, with $opt_p, find pair probs corresponding to mfe pairs
 112   my @pp;
 113   my @sp;
 114   while (<>) {
 115     next unless /(\d+)\s+(\d+)\s+([0-9.Ee-]+)\s+ubox/;
 116     my ($i, $j, $p) = ($1, $2, $3);
 117
 118     $p *= $p;
 119     if ($opt_p) {
 120       $sp[$i] = $sp[$j] = $p if exists  $mfe{$i,$j};
 121     } else {
 122       if ($opt_a) {
 123         $sp[$i] = $sp[$j] = 1-$p if exists  $mfe{$i,$j};
 124       } else {
 125         my $ss = ($p>0)?$p*log($p):0;
 126         $sp[$i] += $ss;
 127         $sp[$j] += $ss;
 128       }
 129     }
 130     $pp[$i] += $p;
 131     $pp[$j] += $p;
 132   }
 133   my $log2 = log(2);
 134   for my $i (1..$n) {
 135     no warnings;  # $p[$i] may be undef
 136     if ($opt_p || $opt_a) {
 137       $sp[$i] = 1-$pp[$i] if !defined $sp[$i];
 138     } else {
 139       $sp[$i] += ($pp[$i]<1) ? (1-$pp[$i])*log(1-$pp[$i]) : 0;
 140       $sp[$i] /= -$log2;
 141     }
 142   }
 143   shift @sp; # get rid of [0] entry
 144   return @sp;
 145 }
 146
 147 =head1 NAME
 148
 149 relplot - annotate a secdonary structure plot with reliability information
 150
 151 =head1 SYNOPSIS
 152
 153    relplot [-p] [-a] file_ss.ps file_dp.ps > file_rss.ps
 154
 155 =head1 DESCRIPTION
 156
 157 relplot reads an RNA secondary structure plot and a dot plot
 158 containing pair probabilities, as produces by C<RNAfold -p>, and
 159 writes a new secondary structure with reliability annotation to
 160 stdout.  The anotation is used to colorize the plot and can use either
 161 "positional entropy" (default), or pair probabilities (with -p).
 162
 163 Positional entropies are computed from the pair probabilities as
 164 C<S(i) = - Sum_i p(ij) log(p(ij))>. Low entropy regions have little
 165 structural flexibility and the reliability of the predicted structure
 166 is high. High entropy implies many structural alternatives. While
 167 these alternatives may be functionally important, they make structure
 168 prediction more difficult and thus less reliable.
 169
 170 If the -p switch is given, the script colors base pairs by their pair
 171 probability, unpaired bases use the probability of being unpaired.
 172
 173 If the -a switch is given, the script colors base pairs by their accessibility, i.e. the probability of being unpaired.
 174
 175 Entropy (repsectively probability) is encoded as color hue, ranging
 176 from red for low entropy, well-defined regions, (high probability
 177 pairs) via yellow and green to blue and violet for regions with very
 178 high entropy (low probability).
 179
 180 You may have to manually move the color legend to a convenient
 181 position. Just edit the postscript file and change the two numbers in the
 182 line reading C<0.1 0.1 colobar>. Or delete the line to remove the legend.
 183
 184 =head1 AUTHOR
 185
 186 Ivo L. Hofacker <ivo@tbi.univie.ac.at>
 187
 188 =cut
 189
 190 #  End of file