binaries/src/ViennaRNA/interfaces/Perl/RNAfold.pl

   1 #!/usr/bin/perl
   2 # -*-Perl-*-
   3 # Last changed Time-stamp: <2006-04-04 14:51:06 ivo>
   4
   5 # Perl implementation of RNAfold
   6 # This is an example script demonstrating how to use the RNA Perl module
   7
   8 use RNA;
   9 use Getopt::Long;
  10 use strict;
  11 use warnings;
  12
  13  Getopt::Long::config("no_ignore_case");
  14
  15 use vars qw/$opt_debug $opt_v $ParamFile $pf $ns_bases/;
  16 my $sfact=1.07;
  17 &usage() unless GetOptions("p|p1" => \$pf,
  18                            "p0" => sub {$pf=1; $RNA::do_backtrack=0},
  19                            "C"   => \$RNA::fold_constrained,
  20                            "T=f" => \$RNA::temperature,
  21                            "4" => sub {$RNA::tetra_loop = 0},
  22                            "d|d0" => sub {$RNA::dangles=0},
  23                            "d2" => sub {$RNA::dangles=2},,
  24                            "noGU" => \$RNA::noGU,
  25                            "noCloseGU" => \$RNA::no_closingGU,
  26                            "noLP" => \$RNA::noLonelyPairs,
  27                            "e=i" => \$RNA::energy_set,
  28                            "P=s" => \$ParamFile,
  29                            "nsp=s" => \$ns_bases,
  30                            "S=f" => \$sfact);
  31
  32  RNA::read_parameter_file($ParamFile) if ($ParamFile);
  33
  34 if ($ns_bases) {
  35    $RNA::nonstandards = "";
  36    foreach my $p ( split(/,/, $ns_bases) ) {
  37         if ($p =~ s/^-//) {
  38             $RNA::nonstandards .= reverse($p)
  39             }
  40         $RNA::nonstandards .= $p;
  41     }
  42     print "$RNA::nonstandards\n";
  43 }
  44
  45 my $istty = (-t STDIN) && (-t STDOUT);
  46 if (($RNA::fold_constrained)&&($istty)) {
  47 print <<END
  48 Input constraints using the following notation:
  49  | : paired with another base
  50  . : no constraint at all
  51  x : base must not pair
  52  < : base i is paired with a base j<i
  53  > : base i is paired with a base j>i
  54 matching brackets ( ): base i pairs base j
  55 END
  56 }
  57
  58 if ($istty) {
  59    print "\nInput string (upper or lower case); @ to quit\n";
  60    for (1..8) { print "....,....$_";}
  61    print "\n";
  62 }
  63 my $fname;
  64 while (<>) {    # main loop: continue until end of file
  65    my ($string, $cstruc, $structure, $min_en);
  66    # skip comment lines and get filenames
  67    if (/^>\s*(\S*)/) {
  68       $fname = $1;
  69       next;
  70    }
  71    last if (/^@/);
  72
  73    if (/(\S+)/) {
  74       $string = $1;
  75    } else {
  76       next;
  77    }
  78
  79    $string = uc($string);
  80    my $length = length($string);
  81    printf("length = %d\n", $length) if ($istty);
  82
  83    if ($RNA::fold_constrained) {
  84       $_ = <>;
  85       $cstruc = $1 if (/(\S+)/);
  86       die("constraint string has wrong length")
  87           if (length($cstruc)!=$length);
  88       ($structure, $min_en) = RNA::fold($string, $structure);
  89    } else {
  90      ($structure, $min_en) = RNA::fold($string);
  91    }
  92    print "$string\n$structure";
  93    if ($istty) {
  94       printf("\n minimum free energy = %6.2f kcal/mol\n", $min_en);
  95    } else {
  96       printf(" (%6.2f)\n", $min_en);
  97    }
  98    my $ffname = ($fname) ? ($fname . '_ss.ps') : 'rna.ps';
  99    RNA::PS_rna_plot($string, $structure, $ffname);
 100
 101    if ($pf) {
 102
 103       # recompute with dangles as in pf_fold()
 104       $RNA::dangles=2 if ($RNA::dangles);
 105       $min_en = RNA::energy_of_struct($string, $structure);
 106
 107       my $kT = ($RNA::temperature+273.15)*1.98717/1000.; # in Kcal
 108       $RNA::pf_scale = exp(-($sfact*$min_en)/$kT/$length);
 109       print STDERR "scaling factor $RNA::pf_scale\n" if ($length>2000);
 110
 111       $structure = $cstruc if ($RNA::fold_constrained);
 112       my $energy = RNA::pf_fold($string, $structure);
 113
 114       if ($RNA::do_backtrack) {
 115          print $structure;
 116          printf(" [%6.2f]\n", $energy) if (!$istty);
 117          print "\n";
 118       }
 119       if (($istty)||(!$RNA::do_backtrack)) {
 120          printf(" free energy of ensemble = %6.2f kcal/mol\n", $energy);
 121          printf(" frequency of mfe structure in ensemble %g\n",
 122                 exp(($energy-$min_en)/$kT));
 123       }
 124
 125       if ($RNA::do_backtrack) {
 126          $ffname = ($fname)?($fname . "_dp.ps"):"dot.ps";
 127          &RNA::PS_dot_plot($string, $ffname);
 128       }
 129    }
 130    undef $fname;
 131 }
 132
 133  RNA::free_pf_arrays() if ($pf);
 134  RNA::free_arrays();
 135
 136 sub usage()
 137 {
 138    die("usage: " .
 139        "RNAfold [-p[0]] [-C] [-T temp] [-4] [-d] [-noGU] [-noCloseGU]\n" .
 140        "               [-e e_set] [-P paramfile] [-nsp pairs] [-S scale]");
 141 }
 142
 143
 144 # End of file