--- /dev/null
+#!/usr/bin/perl -w
+# -*-CPerl-*-
+# Last changed Time-stamp: <2012-11-19 17:22:03 ivo>
+# $Id: relplot.pl,v 1.10 2008/10/09 07:11:21 ivo Exp $
+# colorize a secondary structure plot with reliability annotation
+# from positional entropy
+use strict;
+use Getopt::Std;
+$main::VERSION = 1.3;
+$Getopt::Std::STANDARD_HELP_VERSION=1;
+
+our ($opt_p, $opt_a);
+getopts('pa');
+
+sub HELP_MESSAGE {
+ print STDERR "\nusage: $0 [-p] FOO_ss.ps FOO_dp.ps > FOO_rss.ps\n";
+ print STDERR "For more details run\n\tperldoc -F $0\n";
+}
+
+HELP_MESSAGE() unless $#ARGV >0;
+my $macro_seen= 0;
+my %mfe = (); # hash of mfe pairs
+my @ss_ps = ('',''); # head and tail of the ss.ps file
+
+my $n = swallow_ss_ps(); # read ss plot
+my @sp = posent(); # read dot plot and compute entropies
+my $Smax = ($opt_p || $opt_a) ? 1 : 0;
+if (!$opt_p) {
+ foreach (@sp) {
+ $Smax = $_ if $_>$Smax;
+ }
+ $Smax = ($Smax>0.2) ? sprintf("%3.1f", $Smax) : 0.2 ;
+
+}
+print $ss_ps[0]; # print head
+if (!$macro_seen) {
+ print <<_E_O_F_
+/range 0.8 def
+/drawreliability {
+ /Smax $Smax def
+ 0
+ coor {
+ aload pop
+ S 3 index get
+ Smax div range mul
+ invert {range exch sub} if
+ 1 1 sethsbcolor
+ newpath
+ fsize 2 div 0 360 arc
+ fill
+ 1 add
+ } forall
+} bind def
+/colorbar { % xloc yloc colorbar -> []
+ /STR 8 string def
+ gsave
+ xmin xmax add size sub 2 div
+ ymin ymax add size sub 2 div translate
+ size dup scale
+ translate
+ 0.015 dup scale
+ /tics 64 def
+ gsave
+ 10 tics div 1 scale
+ 0 1 tics
+ {
+ dup 0 moveto 0.5 add
+ tics div range mul
+ invert {range exch sub} if
+ 1 1 sethsbcolor
+ 1 0 rlineto 0 1 rlineto -1 0 rlineto closepath fill
+ } for
+ grestore
+ 0 setgray
+ -0.1 1.01 moveto (0) gsave 0.1 dup scale show grestore
+ 10 1.01 moveto Smax STR cvs
+ gsave 0.1 dup scale dup stringwidth pop -2 div 0 rmoveto show grestore
+ grestore
+} bind def
+_E_O_F_
+}
+print "/S [\n";
+foreach (@sp) {
+ printf " %7.5f\n", $_;
+}
+print "] def\n\n";
+print "/invert ", $opt_p||$opt_a ? 'true' : 'false', " def\n";
+print "drawreliability\n";
+print "0.1 0.1 colorbar\n";
+print $ss_ps[1]; # print tail
+
+sub swallow_ss_ps {
+ # read the secondary structure plot
+ my $length=0;
+ my $tail=0;
+ while (<>) {
+ $macro_seen=1 if /drawreliability /;
+ $length ++ if /^\/coor/ .. /^\] def/;
+ if (/^\/pairs/ .. /^\] def/) {
+ $mfe{$1,$2}=1 if /(\d+)\s+(\d+)/;
+ }
+ $tail=1 if /^drawoutline/;
+ $ss_ps[$tail] .= $_;
+ last if eof;
+ }
+ return $length-2;
+}
+
+sub posent {
+ # compute positional entropy from pair probs in the dot plot file
+ # or, with $opt_p, find pair probs corresponding to mfe pairs
+ my @pp;
+ my @sp;
+ while (<>) {
+ next unless /(\d+)\s+(\d+)\s+([0-9.Ee-]+)\s+ubox/;
+ my ($i, $j, $p) = ($1, $2, $3);
+
+ $p *= $p;
+ if ($opt_p) {
+ $sp[$i] = $sp[$j] = $p if exists $mfe{$i,$j};
+ } else {
+ if ($opt_a) {
+ $sp[$i] = $sp[$j] = 1-$p if exists $mfe{$i,$j};
+ } else {
+ my $ss = ($p>0)?$p*log($p):0;
+ $sp[$i] += $ss;
+ $sp[$j] += $ss;
+ }
+ }
+ $pp[$i] += $p;
+ $pp[$j] += $p;
+ }
+ my $log2 = log(2);
+ for my $i (1..$n) {
+ no warnings; # $p[$i] may be undef
+ if ($opt_p || $opt_a) {
+ $sp[$i] = 1-$pp[$i] if !defined $sp[$i];
+ } else {
+ $sp[$i] += ($pp[$i]<1) ? (1-$pp[$i])*log(1-$pp[$i]) : 0;
+ $sp[$i] /= -$log2;
+ }
+ }
+ shift @sp; # get rid of [0] entry
+ return @sp;
+}
+
+=head1 NAME
+
+relplot - annotate a secdonary structure plot with reliability information
+
+=head1 SYNOPSIS
+
+ relplot [-p] [-a] file_ss.ps file_dp.ps > file_rss.ps
+
+=head1 DESCRIPTION
+
+relplot reads an RNA secondary structure plot and a dot plot
+containing pair probabilities, as produces by C<RNAfold -p>, and
+writes a new secondary structure with reliability annotation to
+stdout. The anotation is used to colorize the plot and can use either
+"positional entropy" (default), or pair probabilities (with -p).
+
+Positional entropies are computed from the pair probabilities as
+C<S(i) = - Sum_i p(ij) log(p(ij))>. Low entropy regions have little
+structural flexibility and the reliability of the predicted structure
+is high. High entropy implies many structural alternatives. While
+these alternatives may be functionally important, they make structure
+prediction more difficult and thus less reliable.
+
+If the -p switch is given, the script colors base pairs by their pair
+probability, unpaired bases use the probability of being unpaired.
+
+If the -a switch is given, the script colors base pairs by their accessibility, i.e. the probability of being unpaired.
+
+Entropy (repsectively probability) is encoded as color hue, ranging
+from red for low entropy, well-defined regions, (high probability
+pairs) via yellow and green to blue and violet for regions with very
+high entropy (low probability).
+
+You may have to manually move the color legend to a convenient
+position. Just edit the postscript file and change the two numbers in the
+line reading C<0.1 0.1 colobar>. Or delete the line to remove the legend.
+
+=head1 AUTHOR
+
+Ivo L. Hofacker <ivo@tbi.univie.ac.at>
+
+=cut
+
+# End of file
git://source.jalview.org / jabaws.git / blobdiff