binaries/src/ViennaRNA/Utils/cmount.pl

   1 #!/usr/bin/perl -w
   2 # -*-Perl-*-
   3 # Last changed Time-stamp: <2005-11-07 13:31:34 ivo>
   4 # Produce coloured Hogeweg mountain representation in PostScript.
   5 # Input is a colour _dp.ps file from alidot or RNAalifold
   6 # definition: mm[i],mp[i]=number of base pairs enclosing base i
   7
   8 use strict;
   9 use vars qw(@mm @mp @hue @sat @pair @pr);
  10 sub usage {
  11     die "Usage: $0 alidot.ps\n";
  12 }
  13
  14 usage() if eof();
  15
  16 print "%!PS-Adobe-2.0 EPSF-1.2
  17 %%Title: Mount Ali
  18 %%Creator: in hiding
  19 %%BoundingBox: 66 209 518 686
  20 %%Pages: 1
  21 %%EndComments\n";
  22
  23
  24 print <<EOF; # PS macros
  25 /trapez {   % use as i j height prob hue sat
  26   dup 0.3 mul 1 exch sub sethsbcolor
  27   newpath
  28   3 index 0.5 sub 2 index moveto % i-0.5 h moveto
  29   dup 1 exch rlineto
  30   4 2 roll exch sub 2 sub 0 rlineto
  31   neg 1 exch rlineto
  32   pop closepath
  33   gsave fill grestore
  34   0 setgray stroke
  35 } def
  36
  37 /centershow {
  38   gsave 1 xs div 1 ys div scale 60 rotate
  39 % dup stringwidth pop 2 div neg 0 rmoveto
  40   show
  41   grestore
  42 } def
  43 EOF
  44
  45 my ($from, $to)=(1,0);
  46 my ($seq, $length);
  47 while (<>) {
  48     chomp;
  49     if (/^% Subsequence from (\d+) to (\d+)/) { # get start and end
  50         $from=$1; $to=$2;
  51     }
  52     if (/\/sequence \{ \((\S*)[\\\)]/) {
  53         print "$_\n";
  54         $seq = $1;              # empty for new version
  55         while (!/\) \} def/) {  # read until end of definition
  56             $_ = <>;
  57             print $_;
  58             /(\S*)[\\\)]/;      # ends either with `)' or `\'
  59             $seq .= $1;
  60         }
  61         print "/len { sequence length } def\n";
  62         $length = length($seq);
  63         next;
  64     }
  65
  66     next if /^%/;
  67     next unless /[ul]box$/;
  68     # Damn! alidot and alifold use different postscript macros.
  69     # Try to recognize both
  70     my ($h, $s, $blah, $i, $j, $p, $tok);
  71     if (/ hsb /) { # alifold version
  72         ($h, $s, $blah, $i, $j, $p, $tok) = split;
  73
  74     } else {             # alidot version
  75         ($i, $j, $p, $h, $s, $tok) = split;
  76     }
  77
  78     if ($tok eq "lbox") { # only read lbox entries
  79         $mp[$i+1]+=$p*$p;
  80         $mp[$j]  -=$p*$p;
  81         $pair[$i] = $j;
  82         $hue[$i]  = $h;
  83         $sat[$i]  = $s;
  84         $pr[$i]   = $p*$p;
  85     }
  86 }
  87 my $max=0;
  88 $mp[0]=$mm[0]=0;
  89 for (my $i=1; $i<=$length; $i++) { #find maximum for scaling
  90     $mp[$i]+=$mp[$i-1];
  91     $max = $mp[$i] if ($mp[$i]>$max);
  92     $mm[$i]+=$mm[$i-1];
  93     $max = $mp[$i] if ($mp[$i]>$max);
  94 }
  95
  96 # postscript scaleing etc
  97 print "72 216 translate\n";
  98 print "/xs {72 6 mul len div} def /ys {72 6 mul $max div} def xs ys scale\n";
  99 print "0.03 setlinewidth
 100 /Times-Roman findfont 1.8 scalefont setfont
 101 0 1 len 1 sub {
 102     dup
 103     0.8 add -0.5 moveto
 104     sequence exch 1 getinterval
 105     gsave 1 $max len div scale show grestore
 106 } for\n\n";
 107
 108 print "/Times-Roman findfont 10 scalefont setfont
 109 0.01 setlinewidth
 110 len log 0.7 sub cvi 10 exch exp  % grid spacing
 111 gsave 0.5 0 translate
 112 /temp 12 string def
 113 0 exch len {
 114    dup dup
 115    0 moveto
 116    $max 1.03 mul lineto
 117    cvi $from 1 sub add temp cvs centershow
 118 } for
 119 stroke
 120 grestore\n";
 121
 122
 123 for (my $i=1; $i<=$length; $i++) { # print pairs as coloured trapezes
 124     next unless ($pair[$i]);
 125     print "$i $pair[$i] ";
 126     printf "%6.4f %6.4f $hue[$i] $sat[$i] trapez\n", $mp[$i], $pr[$i] ;
 127 }
 128
 129 print "showpage\n";
 130
 131 =head1 NAME
 132
 133 cmount.pl - produce a colored mountain plot of a consensus structure
 134
 135 =head1 SYNOPSIS
 136
 137   cmount.pl alidot.ps
 138
 139 =head1 DESCRIPTION
 140
 141 cmount.pl reads a color dot plot as produced by either C<RNAalifold
 142 -p> or C<alidot -p> and writes a postscript colored mountain plot to
 143 stdout.
 144
 145 In the mountain plot a colored trapez with baseline from position i to
 146 j represents the pair (i,j) in the consensus structure. The color of
 147 the trapez encodes the sequence variation at that pair: Red marks
 148 pairs with no sequence variation; ochre, green, turquoise, blue, and
 149 violet mark pairs with 2,3,4,5,6 different tpyes of pairs,
 150 respectively.
 151
 152 =cut