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
10 $Getopt::Std::STANDARD_HELP_VERSION=1;
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";
20 HELP_MESSAGE() unless $#ARGV >0;
22 my %mfe = (); # hash of mfe pairs
23 my @ss_ps = ('',''); # head and tail of the ss.ps file
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;
30 $Smax = $_ if $_>$Smax;
32 $Smax = ($Smax>0.2) ? sprintf("%3.1f", $Smax) : 0.2 ;
35 print $ss_ps[0]; # print head
46 invert {range exch sub} if
54 /colorbar { % xloc yloc colorbar -> []
57 xmin xmax add size sub 2 div
58 ymin ymax add size sub 2 div translate
69 invert {range exch sub} if
71 1 0 rlineto 0 1 rlineto -1 0 rlineto closepath fill
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
84 printf " %7.5f\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
93 # read the secondary structure plot
97 $macro_seen=1 if /drawreliability /;
98 $length ++ if /^\/coor/ .. /^\] def/;
99 if (/^\/pairs/ .. /^\] def/) {
100 $mfe{$1,$2}=1 if /(\d+)\s+(\d+)/;
102 $tail=1 if /^drawoutline/;
110 # compute positional entropy from pair probs in the dot plot file
111 # or, with $opt_p, find pair probs corresponding to mfe pairs
115 next unless /(\d+)\s+(\d+)\s+([0-9.Ee-]+)\s+ubox/;
116 my ($i, $j, $p) = ($1, $2, $3);
120 $sp[$i] = $sp[$j] = $p if exists $mfe{$i,$j};
123 $sp[$i] = $sp[$j] = 1-$p if exists $mfe{$i,$j};
125 my $ss = ($p>0)?$p*log($p):0;
135 no warnings; # $p[$i] may be undef
136 if ($opt_p || $opt_a) {
137 $sp[$i] = 1-$pp[$i] if !defined $sp[$i];
139 $sp[$i] += ($pp[$i]<1) ? (1-$pp[$i])*log(1-$pp[$i]) : 0;
143 shift @sp; # get rid of [0] entry
149 relplot - annotate a secdonary structure plot with reliability information
153 relplot [-p] [-a] file_ss.ps file_dp.ps > file_rss.ps
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).
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.
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.
173 If the -a switch is given, the script colors base pairs by their accessibility, i.e. the probability of being unpaired.
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).
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.
186 Ivo L. Hofacker <ivo@tbi.univie.ac.at>