Revision 548c817f948d808f6b2b9857773dc05f01fa16b7 authored by Ivo Hofacker on 03 November 2009, 15:29:09 UTC, committed by Ivo Hofacker on 03 November 2009, 15:29:09 UTC
1 parent 4528ead
RNAfold.pl
#!/usr/bin/perl
# -*-Perl-*-
# Last changed Time-stamp: <2006-04-04 14:51:06 ivo>
# Perl implementation of RNAfold
# This is an example script demonstrating how to use the RNA Perl module
use RNA;
use Getopt::Long;
use strict;
use warnings;
Getopt::Long::config("no_ignore_case");
use vars qw/$opt_debug $opt_v $ParamFile $pf $ns_bases/;
my $sfact=1.07;
&usage() unless GetOptions("p|p1" => \$pf,
"p0" => sub {$pf=1; $RNA::do_backtrack=0},
"C" => \$RNA::fold_constrained,
"T=f" => \$RNA::temperature,
"4" => sub {$RNA::tetra_loop = 0},
"d|d0" => sub {$RNA::dangles=0},
"d2" => sub {$RNA::dangles=2},,
"noGU" => \$RNA::noGU,
"noCloseGU" => \$RNA::no_closingGU,
"noLP" => \$RNA::noLonelyPairs,
"e=i" => \$RNA::energy_set,
"P=s" => \$ParamFile,
"nsp=s" => \$ns_bases,
"S=f" => \$sfact);
RNA::read_parameter_file($ParamFile) if ($ParamFile);
if ($ns_bases) {
$RNA::nonstandards = "";
foreach my $p ( split(/,/, $ns_bases) ) {
if ($p =~ s/^-//) {
$RNA::nonstandards .= reverse($p)
}
$RNA::nonstandards .= $p;
}
print "$RNA::nonstandards\n";
}
my $istty = (-t STDIN) && (-t STDOUT);
if (($RNA::fold_constrained)&&($istty)) {
print <<END
Input constraints using the following notation:
| : paired with another base
. : no constraint at all
x : base must not pair
< : base i is paired with a base j<i
> : base i is paired with a base j>i
matching brackets ( ): base i pairs base j
END
}
if ($istty) {
print "\nInput string (upper or lower case); @ to quit\n";
for (1..8) { print "....,....$_";}
print "\n";
}
my $fname;
while (<>) { # main loop: continue until end of file
my ($string, $cstruc, $structure, $min_en);
# skip comment lines and get filenames
if (/^>\s*(\S*)/) {
$fname = $1;
next;
}
last if (/^@/);
if (/(\S+)/) {
$string = $1;
} else {
next;
}
$string = uc($string);
my $length = length($string);
printf("length = %d\n", $length) if ($istty);
if ($RNA::fold_constrained) {
$_ = <>;
$cstruc = $1 if (/(\S+)/);
die("constraint string has wrong length")
if (length($cstruc)!=$length);
($structure, $min_en) = RNA::fold($string, $structure);
} else {
($structure, $min_en) = RNA::fold($string);
}
print "$string\n$structure";
if ($istty) {
printf("\n minimum free energy = %6.2f kcal/mol\n", $min_en);
} else {
printf(" (%6.2f)\n", $min_en);
}
my $ffname = ($fname) ? ($fname . '_ss.ps') : 'rna.ps';
RNA::PS_rna_plot($string, $structure, $ffname);
if ($pf) {
# recompute with dangles as in pf_fold()
$RNA::dangles=2 if ($RNA::dangles);
$min_en = RNA::energy_of_struct($string, $structure);
my $kT = ($RNA::temperature+273.15)*1.98717/1000.; # in Kcal
$RNA::pf_scale = exp(-($sfact*$min_en)/$kT/$length);
print STDERR "scaling factor $RNA::pf_scale\n" if ($length>2000);
$structure = $cstruc if ($RNA::fold_constrained);
my $energy = RNA::pf_fold($string, $structure);
if ($RNA::do_backtrack) {
print $structure;
printf(" [%6.2f]\n", $energy) if (!$istty);
print "\n";
}
if (($istty)||(!$RNA::do_backtrack)) {
printf(" free energy of ensemble = %6.2f kcal/mol\n", $energy);
printf(" frequency of mfe structure in ensemble %g\n",
exp(($energy-$min_en)/$kT));
}
if ($RNA::do_backtrack) {
$ffname = ($fname)?($fname . "_dp.ps"):"dot.ps";
&RNA::PS_dot_plot($string, $ffname);
}
}
undef $fname;
}
RNA::free_pf_arrays() if ($pf);
RNA::free_arrays();
sub usage()
{
die("usage: " .
"RNAfold [-p[0]] [-C] [-T temp] [-4] [-d] [-noGU] [-noCloseGU]\n" .
" [-e e_set] [-P paramfile] [-nsp pairs] [-S scale]");
}
# End of file
![swh spinner](/static/img/swh-spinner.gif)
Computing file changes ...