rf-fold

#!/usr/bin/env perl

##
# RF Fold
# RNA Framework [http://www.rnaframework.com]
#
# Author:  Danny Incarnato (dincarnato[at]rnaframework.com)
# Summary: Produces secondary structures for analyzed transcripts using structural
#          probing data to guide folding
#
# This program is free software, and can be redistribute  and/or modified
# under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# any later version.
#
# Please see <http://www.gnu.org/licenses/> for more informations.
##

use strict;
use Config;
use Cwd;
use File::Basename;
use File::Spec;
use FindBin qw($Bin);
use Getopt::Long qw(:config no_ignore_case);
use POSIX qw(floor);

use lib $Bin . "/lib";

use Core::Mathematics;
use Core::Statistics;
use Core::Process::Queue;
use Core::Utils;
use Data::IO;
use Data::IO::Sequence;
use Data::Sequence;
use Data::Sequence::Structure;
use Data::Sequence::Utils;
use Graphics::Chart::Arcs;
use Graphics::Chart::Area;
use Graphics::Chart::Barplot;
use Graphics::Image;
use Graphics::RNA;
use Interface::ViennaRNA;
use RNA::Utils;
use RF::Data::IO::XML;
use RF::Utils;
use Term::Constants qw(:screen);
use Term::Progress;

use constant MINWINLENGTH => 50; # Minimum window length for windowed folding

$|++;

my ($output, $ct, $method, $installed, $maxdist, 
    $threads, $img, $t,
    $cutoff, $nlp, $ngu, $rnastructure, $R,
    $viennarna, $data, $slope, $intercept,
    $overwrite, $help, $tmpdir, $ignore, $hc, 
    $progressBar, $processManager, $windowed, $cwd, 
    $partition, $probplot, $partwin, $partstep, $wintrim, 
    $foldwin, $foldstep, $dp, $sh, $pk, $suboptDeltaEnergy,
    $pkwin, $pkstep, $p1, $p2, $dotplotOnly,
    $pkcutoff, $rnasubopt, $pkmax, $pksubopt,
    $shapeknots, $pkmethod, $constraintDir, 
    $pkpercentmfe, $nozuker, $zukersubopt, $ignoreSeq, 
    $onlyCommon, $gradient, $vrna, $useAlifold, 
    @skparams, @vrfoldparams, @rsfoldparams, @slope,
    @intercept, @samples, %xml, %lambda, %constraints, 
    %results);

%results = ( folded    => 0,
             parseErr  => 0,
             foldErr   => 0,
             pkErr     => 0,
             diffSeq   => 0,
             diffConst => 0 );
%lambda = ( "2"  => 0,
            "3"  => 24,
            "4"  => 24,
            "5"  => 34,
            "6"  => 24,
            "7"  => 35,
            "8"  => 58,
            "9"  => 82,
            "10" => 65,
            "11" => 527,
            "12" => 2447,
            "13" => 4199,
            "14" => 6564,
            "15" => 12540 );

do {

    local $SIG{__WARN__} = sub { };

    GetOptions( "h|help"                     => \$help,
                "o|output-dir=s"             => \$output,
                "ow|overwrite"               => \$overwrite,
                "ct|connectivity-table"      => \$ct,
                "m|folding-method=i"         => \$method,
                "p|processors=i"             => \$threads,
                "g|img"                      => \$img,
                "c|constraints=s"            => \$constraintDir,
                "t|temperature=s"            => \$t,
                "sl|slope=s"                 => \$slope,
                "in|intercept=s"             => \$intercept,
                "md|maximum-distance=i"      => \$maxdist,
                "vrf|vienna-rnafold=s"       => \$viennarna,
                "f|cutoff=s"                 => \$cutoff,
                "nlp|no-lonelypairs"         => \$nlp,
                "ua|use-alifold"             => \$useAlifold,
                "ngu|no-closing-gu"          => \$ngu,
                "rs|rnastructure=s"          => \$rnastructure,
                "d|data-path=s"              => \$data,
                "i|ignore-reactivity"        => \$ignore,
                "is|ignore-sequence"         => \$ignoreSeq,
                "oc|only-common=s"           => \$onlyCommon,
                "hc|hard-constraint"         => \$hc,
                "w|windowed"                 => \$windowed,
                "fw|fold-window=i"           => \$foldwin,
                "fo|fold-offset=i"           => \$foldstep,
                "wt|window-trim=i"           => \$wintrim,
                "pw|partition-window=i"      => \$partwin,
                "po|partition-offset=i"      => \$partstep,
                "dp|dotplot"                 => \$dp,
                "dpo|dotplotOnly"            => \$dotplotOnly,
                "kt|pseudoknot-tollerance=s" => \$pkpercentmfe,
                "sh|shannon-entropy"         => \$sh,
                "pk|pseudoknots"             => \$pk,
                "pt|partition=s"             => \$partition,
                "pp|probabilityplot=s"       => \$probplot,
                "vrs|vienna-rnasubopt=s"     => \$rnasubopt,
                "kw|pseudoknot-window=i"     => \$pkwin,
                "ko|pseudoknot-offset=i"     => \$pkstep,
                "kp1|pseudoknot-penalty1=s"  => \$p1,
                "kp2|pseudoknot-penalty2=s"  => \$p2,
                "kc|pseudoknot-cutoff=s"     => \$pkcutoff,
                "ks|pseudoknot-suboptimal=i" => \$pksubopt,
                "kh|pseudoknot-helices=i"    => \$pkmax,
                "km|pseudoknot-method=i"     => \$pkmethod,
                "sk|shapeknots=s"            => \$shapeknots,
                "nz|no-zuker"                => \$nozuker,
                "zs|zuker-suboptimal=i"      => \$zukersubopt,
                "ke|subopt-delta-energy=s"   => \$suboptDeltaEnergy,
                "R|R-path=s"                 => \$R,
                "gr|gradient"                => \$gradient ) or help(1);

    @samples = uniq(@ARGV);

};

help() if ($help);

# Default
$cwd = cwd();
$output ||= "rf_fold/";
$method ||= 1;
$threads ||= 1;
$t //= 37.0;
$cutoff //= 0.7;
$slope //= 1.8;
$intercept //= -0.6;
$maxdist //= 0;
$partwin ||= 600;
$partstep ||= 200;
$wintrim //= 100;
$foldwin ||= 600;
$foldstep ||= 200;
$pkwin ||= 600;
$pkstep ||= 200;
$p1 //= 0.35;
$p2 //= 0.65;
$pkcutoff //= 0.5;
$pkmax ||= 100;
$pksubopt ||= 1000;
$zukersubopt ||= 1000;
$pkmethod ||= 1;
$pkpercentmfe ||= 0.5;
$suboptDeltaEnergy ||= 1;
$viennarna = $viennarna || which("RNAfold");
$rnasubopt = $rnasubopt || which("RNAsubopt");
$rnastructure = $rnastructure || which("Fold-smp") || which("Fold");
$shapeknots = $shapeknots || which("ShapeKnots-smp") || which("ShapeKnots");
$partition = $partition || which("partition-smp") || which("partition");
$probplot = $probplot || which("ProbabilityPlot");
$data ||= $ENV{DATAPATH};
$R = checkRinstall($R) if ($img);
$installed = eval { require RNA; 1; };

$output =~ s/\/?$/\//;
$tmpdir = $output . "/tmp/";
$constraintDir =~ s/\/?$/\// if (defined $constraintDir);

for ($viennarna, $rnasubopt, $rnastructure, 
     $shapeknots, $partition, $probplot) { $_ = File::Spec->rel2abs($_) if (defined $_); }

##
# Input validation
##

die "\n  [!] Error: No output directory specified\n\n" unless(defined $output);
die "\n  [!] Error: No XML directory\/file specified\n\n" if (!@samples);
die "\n  [!] Error: Provided constraint directory doesn't exist\n\n" if (defined $constraintDir && !-d $constraintDir);
die "\n  [!] Error: Number of processors must be an integer greater than 0\n\n" if ($threads < 1);
die "\n  [!] Error: Invalid folding method\n\n" unless($method =~ m/^[12]$/);
die "\n  [!] Error: Temperature must be a positive value\n\n" unless(ispositive($t));
die "\n  [!] Error: Temperature must be comprised between 0 and 100 degree Celsius\n\n" if ($t > 100);
die "\n  [!] Error: Maximum distance value must be a positive integer\n\n" if (!ispositive($maxdist));
die "\n  [!] Error: No closing GU parameter requires folding method 1 (ViennaRNA)\n\n" if ($ngu && $method != 1);
die "\n  [!] Error: Reactivity cutoff must be a value >0 and <= 1\n\n" if (!inrange($cutoff, [0, 1]));

if ($onlyCommon) {

    die "\n  [!] Error: Parameter -oc value must be a positive INT >= 1\n\n" if (!isint($onlyCommon) || $onlyCommon < 1);
    die "\n  [!] Error: Parameter -oc value cannot exceed the number of samples (" . scalar(@samples) . ")\n\n" if ($onlyCommon > scalar(@samples));

}

if (!$ignore) {

    @slope = split /,/, $slope;
    @intercept = split /,/, $intercept;
    @slope = ($slope[0]) x scalar(@samples) if (@slope == 1 && @samples > 1);
    @intercept = ($intercept[0]) x scalar(@samples) if (@intercept == 1 && @samples > 1);

    die "\n  [!] Error: Number of provided slope values differs from number of samples\n\n" if (@slope != @samples);
    die "\n  [!] Error: Number of provided intercept values differs from number of samples\n\n" if (@intercept != @samples);

    for (@slope) { die "\n  [!] Error: Slope must be numeric\n\n" if (!isnumeric($_)); }
    for (@intercept) { die "\n  [!] Error: Intercept must be numeric\n\n" if (!isnumeric($_)); }

}

if ($useAlifold) {

    die "\n  [!] Error: RNAalifold-based folding requires folding method 1 (ViennaRNA)\n\n" if ($method != 1);
    die "\n  [!] Error: RNAalifold does not support sample-specific slopes\n\n" if (uniq(@slope) > 1);
    die "\n  [!] Error: RNAalifold does not support sample-specific intercepts\n\n" if (uniq(@intercept) > 1);

}

if ($windowed) {

    die "\n  [!] Error: Folding window's size must be an integer >= 50\n\n" unless($foldwin >= MINWINLENGTH);
    die "\n  [!] Error: Folding window's sliding offset cannot be greater than window's size\n\n" if ($foldstep > $foldwin);
    die "\n  [!] Error: Parition function window's size must be an integer >= 50\n\n" unless($partwin >= MINWINLENGTH);
    die "\n  [!] Error: Partition function window's sliding offset cannot be greater than window's size\n\n" if ($partstep > $partwin);
    die "\n  [!] Error: Partition function window's trimming value cannot exceed half of" .
        "\n             the difference between window's size and sliding offset size\n\n" if ($wintrim > int(($partwin - $partstep) / 2));

}

if ($pk) {

    die "\n  [!] Error: Invalid pseudoknots prediction method\n\n" unless($pkmethod =~ m/^[12]$/);
    die "\n  [!] Error: Pseudoknot window's size must be an integer >= 50\n\n" unless($pkwin >= MINWINLENGTH);
    die "\n  [!] Error: Pseudoknot window's sliding offset cannot be greater than window's size\n\n" if ($pkstep > $pkwin);
    die "\n  [!] Error: Pseudoknot penalty 1 must be numeric\n\n" unless(isnumeric($p1));
    die "\n  [!] Error: Pseudoknot penalty 2 must be numeric\n\n" unless(isnumeric($p2));
    die "\n  [!] Error: Pseudoknot reactivity cutoff must be > 0 and <= 1\n\n" if (!inrange($pkcutoff, [0, 1]));
    die "\n  [!] Error: Pseudoknot suboptimal structures energy tollerance must be > 0 and <= 1\n\n" if (!isnumeric($pkpercentmfe) || !inrange($pkpercentmfe, [0, 1]));
    die "\n  [!] Error: Pseudoknot suboptimal structure energy delta must be >= 0\n\n" if (!ispositive($suboptDeltaEnergy));

}

die "\n  [!] Error: ViennaRNA package's Perl bindings are required." .
    "\n             Please ensure that ViennaRNA package v2.7.0 (or greater) is installed and try again\n\n" unless($installed);

print "\n[+] Checking method's requirements...";

if ($method == 1) { # ViennaRNA

    if (!defined $viennarna) { die "\n\n  [!] Error: RNAfold is not in PATH\n\n"; }
    elsif (!-e $viennarna) { die "\n\n  [!] Error: RNAfold doesn't exist\n\n"; }
    elsif (!-x $viennarna) { die "\n\n  [!] Error: RNAfold is not executable\n\n"; }

    my $ret = `$viennarna --version`;

    if ($ret =~ /RNAfold (\d+)\.(\d+)/) {

        my ($v1, $v2) = ($1, $2);

        die "\n\n  [!] Error: RF Fold requires ViennaRNA package v2.5.0 or greater (Detected: v" . $v1 . "." . $v2 . ")\n\n" if ($v1 < 2 || $v2 < 5);

    }
    else { warn "\n\n  [!] Warning: Unable to detect ViennaRNA package version\n"; }

}
else { # RNAstructure

    if (!defined $data) { die "\n\n  [!] Error: Environment variable DATAPATH is not set\n\n"; }
    elsif (!-d $data) { die "\n\n  [!] Error: Provided DATAPATH directory doesn't exist\n\n"; }

    foreach my $prog (["Fold", $rnastructure], ["partition", $partition], ["ProbabilityPlot", $probplot]) {

        if (!defined $prog->[1]) { die "\n\n  [!] Error: RNAstructure $prog->[0] is not in PATH\n\n"; }
        elsif (!-e $prog->[1]) { die "\n\n  [!] Error: RNAstructure $prog->[0] does not exist\n\n"; }
        elsif (!-x $prog->[1]) { die "\n\n  [!] Error: RNAstructure $prog->[0] is not executable\n\n"; }
        else {

            my $ret = `$prog->[1] --version`;
            
            if ($ret =~ /Mathews Lab/i) {

                if ($ret =~ /Version (\d+)\.(\d+)/) {         
                    
                    my ($v1, $v2) = ($1, $2);

                    die "\n\n  [!] Error: RF Fold requires RNAstructure $prog->[0] v6.5 or greater (Detected: v" . $v1 . "." . $v2 . ")\n\n" if ($v1 < 6 || ($v1 == 6 && $v2 < 5));

                }
                else { warn "\n\n  [!] Warning: Unable to detect RNAstructure $prog->[0] version\n"; }

            }
            else { warn "\n\n  [!] Warning: $prog->[0] (path: $prog->[1]) does not look like an RNAstructure tool\n"; }

        }

    }

    $ENV{DATAPATH} = $data;

}

if ($pk) {

    if ($pkmethod == 1) {

        if (!defined $rnasubopt) { die "\n\n  [!] Error: RNAsubopt is not in PATH\n\n"; }
        elsif (!-e $rnasubopt) { die "\n\n  [!] Error: RNAsubopt doesn't exist\n\n"; }
        elsif (!-x $rnasubopt) { die "\n\n  [!] Error: RNAsubopt is not executable\n\n"; }

    }
    else {

        if (!defined $shapeknots) { die "\n\n  [!] Error: ShapeKnots is not in PATH\n\n"; }
        elsif (!-e $shapeknots) { die "\n\n  [!] Error: ShapeKnots doesn't exist\n\n"; }
        elsif (!-x $shapeknots) { die "\n\n  [!] Error: ShapeKnots is not executable\n\n"; }

    }

}

# Setting folding parameters
if ($method == 1 || $pkmethod == 1) { # ViennaRNA

    my %vrfoldparams = ( temperature   => $t,
                         noLonelyPairs => $nlp ? 1 : 0,
                         noClosingGU   => $ngu ? 1 : 0,
                         maxBPspan     => $maxdist );

    for my $i (0 .. $#samples) {

        my %tmpParams = %vrfoldparams;

        if (!$ignore) {

            $tmpParams{slope} = $slope[$i];
            $tmpParams{intercept} = $intercept[$i];

        }

        push(@vrfoldparams, \%tmpParams);

    }

}

if ($method == 2 || $pkmethod == 2) { # RNAstructure

    my $rsfoldparams = " -t " . ($t + 273.15);
    $rsfoldparams .= " -md " . $maxdist if ($maxdist);

    for my $i (0 .. $#samples) {

        my ($tmpParams);
        $tmpParams = " -si " . $intercept[$i] . " -sm " . $slope[$i] if (!$ignore);
        push(@rsfoldparams, $rsfoldparams . $tmpParams);
        push(@skparams, $tmpParams);
        
    }

}

$RNA::temperature = $t;

# Output directory tree
print "\n[+] Making output directory tree...";

if (-e $output) {

    if ($overwrite) {

        my $error = rmtree($output);

        die "\n\n  [!] Error: " . $error . "\n\n" if ($error);

    }
    else { die "\n\n  [!] Error: Output directory already exists." .
               "\n      Please use -ow (or --overwrite) to overwrite output directory\n\n"; }

}

if (my $error = mktree($tmpdir, "755")) { die "\n\n  [!] Error: Unable to create output directory ($error)\n\n"; }

if (!$dotplotOnly) {

    mktree($output . "structures", "755");
    mktree($output . "shannon", "755") if ($sh);

}

if ($img) {

    mktree($output . "plots/summaries", "755");
    mktree($output . "plots/structures", "755");                                     

}

mktree($output . "dotplot", "755") if ($dp);

# Sets all paths to absolute, before changing working directory to temp dir
$output = File::Spec->rel2abs($output) . "/";
$tmpdir = File::Spec->rel2abs($tmpdir) . "/";

print "\n[+] Importing XML file(s)...";

foreach my $sample (@samples) {

    die "\n\n  [!] Error: Sample XML file/directory \"$sample\" does not exist\n\n" if (!-e $sample);

    if (-d $sample) {

        opendir(my $dh, $sample) or die "\n\n  [!] Error: Unable to read XML files from directory \"$sample\"\n\n";
        while(my $file = readdir($dh)) {

            next if (substr($file, -4) ne ".xml");

            push(@{$xml{(fileparse($file, ".xml"))[0]}}, File::Spec->rel2abs($sample . "/" . $file));

        }
        closedir($dh);

    }
    else { # Single xml file

        die "\n\n  [!] Error: Provided file \"$sample\" lacks XML extension\n\n" if (substr($sample, -4) ne ".xml");

        push(@{$xml{(fileparse($sample, ".xml"))[0]}}, File::Spec->rel2abs($sample));

    }

}

%xml = map { $_ => $xml{$_} } grep { @{$xml{$_}} >= $onlyCommon } keys %xml if ($onlyCommon);

print " " . scalar(keys %xml) . " imported (of which " . scalar(grep { @{$xml{$_}} >= ($onlyCommon ? $onlyCommon : @samples) } keys %xml) . " common to " . ($onlyCommon ? "at least $onlyCommon" : "all") . " replicates)";

die "\n\n  [!] Error: Specified folder doesn't contain any XML file\n\n" unless (keys %xml);

warn "\n\n  [i] Ignoring XML reactivity data files. Predicting MFE unconstrained structures...\n" if ($ignore);

if ($constraintDir) {

    my $err = 0;

    print "\n[+] Importing constraint file(s)...";

    opendir(my $dh, $constraintDir) or die "\n\n  [!] Error: Unable to read constraint files from directory\n\n";
    while(my $file = readdir($dh)) {

        next if (substr($file, 0, 1) eq "." || substr($file, -3) !~ /^\.(?:db|ct)$/);

        my $eval = do { 

            local $@;

            eval {

                my ($structIO, $entry, $structure);
                $structIO = Data::IO::Sequence->new( file        => $constraintDir . $file,
                                                     lonelypairs => 1 );
                $entry = $structIO->read();
                $structure = $entry->structure();
                $constraints{$entry->id()} = $structure;

            };

            $@; 

        };

        $err++ if ($eval);

    }
    closedir($dh);

    print " " . scalar(keys %constraints) . " imported";
    print ", $err discarded" if ($err);

}

chdir($tmpdir);  # Change the working directory before starting threads

$processManager = Core::Process::Queue->new( processors => $threads,
                                             stderr     => $output . "error.out",
                                             tmpDir     => $tmpdir,
                                             verbosity  => 1 );
$processManager->enqueue( command   => \&fold,
                          arguments => $xml{$_},
                          id        => $_ ) for (keys %xml);
$progressBar = Term::Progress->new( max     => $processManager->queueSize(),
                                    colored => 1 );

print "\n[+] Folding RNA structures. Please wait...\n\n";
                                
$progressBar->init();
$processManager->parentOnExit(sub { 
    
    my ($file, $pid, $uid, $ret);
    ($file, $pid, $uid) = @_[0..2];
    $ret = ($processManager->dequeue($uid)->exitcode())[0];

    $progressBar->update(1); 
    
    if (ref($ret) eq "HASH") { 
        
        my $err = (keys %$ret)[0];
        $results{$err} += $ret->{$err};
        
    }
    else { $results{"folded"}++; }
    
});
$processManager->start();
$processManager->waitall();

chdir($cwd);  # Reset to original working directory

my $otherErr = scalar(keys %xml);
$otherErr -= $_ for (values %results);

print "\n\n[+] Folding statistics:\n" .
      "\n  [*] Folded transcripts:    " . $results{"folded"} .
      "\n  [*] Discarded transcripts: " . (scalar(keys %xml) - $results{"folded"}) . " total" .
      "\n                             " . $results{"parseErr"} . " XML parsing failed" .
      "\n                             " . $results{"foldErr"} . " folding failed";
print "\n                             " . $results{"diffConst"} . " wrong constraint length" if (keys %constraints);
print "\n                             " . $results{"diffSeq"} . " different sequence between replicates" if (@samples > 1);
print "\n                             " . $otherErr . " other" if ($otherErr);
print "\n                             " . $results{"pkErr"} . " pseudoknotted folding failed" if ($pk);

rmtree($tmpdir);

if (!-s $output . "error.out") { unlink($output . "error.out"); }
else { print "\n\n  [!] Warning: Execution completed with error(s)/warning(s). Please check the \"${output}error.out\" file"; }

print "\n\n[+] All done.\n\n";

sub fold { 

    my @xml = @_;

    my ($sequence, $id, $escapedId, $structure, $cmd, $plot,
        $ret, $entry, $length, $nwins, $file, $io, $energy,
        $constraint, @foldwins, @windows, @constraint,
        @pkpairs, @keptpairs, @dotplot, @shannon,
        @reactivity, @sequences, @meanReact, %pairprobs,
        %npairs, %seen);

    if ($img) {
 
        $plot = Graphics::Image->new( width  => 6 * 3,
                                      height => ($windowed || $dp || $sh ? 7 : 5) * 3,
                                      R      => $R,
                                      tmpdir => $tmpdir );

    }

    # The interface is instantiated here (and not in the parent) so that each
    # forked process gets its own set of temporary files
    if ($method == 1 || $pkmethod == 1) {

        eval { $vrna = Interface::ViennaRNA->new( RNAfold   => $viennarna,
                                                  RNAsubopt => $rnasubopt,
                                                  tmpdir    => $tmpdir,
                                                  verbosity => $ENV{"RF_VERBOSITY"} // -1 ); };

        if ($@) {

            _logError($@, $xml[0]);

            return({ foldErr => 1 });

        }

    }

    foreach my $xml (@xml) {

        my ($xmlref);
    
        do {
            
            local $@;

            eval { $xmlref = RF::Data::IO::XML->new(file => $xml); };

            if ($@) { 

                _logError($@, $xml);

                return({ parseErr => 1 }); 
                
            } 
            
        };

        if (!defined $id) { 
            
            $id = $xmlref->id();
            $escapedId = $id;
            $escapedId =~ s/[^\w_\-]/_/g; # To handle temporary files for ids with non-canonical characters
            $sequence = $xmlref->sequence(); 
            $length = $xmlref->length();
            
        }
        else { 

            if (($sequence ne $xmlref->sequence() && !$ignoreSeq) ||
                ($length != $xmlref->length() && $ignoreSeq)) {

                _logError("Sequence mismatch between input XML files", $id);

                return({ diffSeq => 1 });

            }

        }

        push(@reactivity, [ $xmlref->reactivity() ]);
        push(@sequences, $xmlref->sequence());

    }

    @meanReact = map { my $i = $_; [ map { $_->[$i] } @reactivity ] } 0 .. $#{$reactivity[0]};
    @meanReact = map { isnumeric(@$_) ? mean(@$_) : "NaN" } @meanReact;
    $constraint = exists $constraints{$id} ? $constraints{$id} : "." x $length;

    if (length($constraint) != $length) {

        _logError("Constraint length != sequence length", $id);

        return({ diffConst => 1 });

    }

    ######### Partition function block ##########

    # We enter here only if we need to do windowed analysis, or if
    # the user requests the bp dotplot or the shannon entropy
    if ($windowed || $sh || $dp) {

        $nwins = $windowed ? max(1, floor(($length - $partwin) / $partstep) + 1) : 1;

        # Define windows to be analyzed
        # Windows are defined as: [start, end, 5'-trim, 3'-trim]
        if ($nwins > 1) {

            for (50, 100) {

                my $end = $partwin - $_ - 1;

                last if ($end < MINWINLENGTH - 1);

                push(@foldwins, [0, $end, 0, $wintrim]);

            }

            for (0 .. $nwins - 1) {

                my ($start, $end);
                $start = $partstep * $_;
                $end = $_ == $nwins - 1 ? $length - 1 : $start + $partwin - 1;
                push(@foldwins, [$start, $end, $wintrim, $wintrim]);

            }

            for (50, 100) {

                my $start = $length - $partwin + $_;

                last if ($start > $length - MINWINLENGTH);

                push(@foldwins, [$start, $length - 1, $wintrim, 0]);

            }

        }
        else { push(@foldwins, [0, $length - 1, 0, 0]); }

        # This fixes the cases in which the first or last windows are < MINWINLENGTH
        $foldwins[0]->[2] = 0;
        $foldwins[-1]->[3] = 0;

        foreach my $sampleSet (samplesets(scalar(@reactivity))) {

            # Calculate window partition function
            foreach my $win (@foldwins) {

                my ($start, $end, $trim5, $trim3, $winconst,
                    $subseq, $tmpfile, $alignment, $alnReact, @winreact);
                ($start, $end, $trim5, $trim3) = @{$win};

                $subseq = substr($sequence, $start, $end - $start + 1);
                $winconst = fixdotbracket(substr($constraint, $start, $end - $start + 1));
                $tmpfile = "$escapedId\_$start-$end\_$$\_" . ($useAlifold ? "ali" : $sampleSet->[0]);

                if ($useAlifold) {

                    ($alignment, $alnReact) = winalignment($start, $end, \@sequences, \@reactivity);
                    @winreact = @meanReact[$start .. $end]; # Hard constraints are derived from the mean reactivity profile

                }
                else { @winreact = @{$reactivity[$sampleSet->[0]]}[$start .. $end]; }

                my ($win, $shannon) = partition( file          => $tmpfile,
                                                 sample        => $sampleSet->[0],
                                                 sequence      => $subseq,
                                                 reactivity    => \@winreact,
                                                 alignment     => $alignment,
                                                 alnReactivity => $alnReact,
                                                 pos           => $start,
                                                 trim5         => $trim5,
                                                 trim3         => $trim3,
                                                 constraint    => $winconst );

                if (!$win) {

                    unlink(glob($tmpfile . "*"));

                    return({ foldErr => 1 });

                }

                push(@{$pairprobs{$_}->{pairs}}, $win->{$_}) for (keys %{$win});
                push(@{$shannon[$_]}, $shannon->{$_}) for (keys %{$shannon});
                push(@windows, [$start + $trim5, $end - $trim3]);

            }

        }

        @shannon = map { mean(@{$_}) } @shannon; # Shannon entropy vector

        # Just do this for windowed analysis, otherwise just fold
        if ($windowed) {

            # Build constraint with pairs > 0.99 probability
            foreach my $pair (sort keys %pairprobs) {

                my ($i, $j, $prob);
                ($i, $j) = split(/-/, $pair);

                next if (exists $seen{$i} || exists $seen{$j});

                for (@windows) { $pairprobs{$pair}->{n}++ if (inrange($i, $_) && inrange($j, $_)); }

                $prob = sum(@{$pairprobs{$pair}->{pairs}}) / $pairprobs{$pair}->{n};

                # If average p(i,j) is >= 0.99
                if ($prob >= 0.99) {

                    push(@constraint, [$i, $j, $prob]);
                    $seen{$i} = 1;
                    $seen{$j} = 1;

                }

            }

        }
        else { @constraint = map { [@$_, 1] } listpairs($constraint); }

        undef($nwins);
        undef(@foldwins);
        undef(@windows);

        # We can take advantage of the rmpseudoknots() function to discard incompatible pairs
        # The rmpseudoknots() function takes the sequence, an array ref of base-pairs, and
        # optionally a scoring function to be applied (in this case the average of bp probs)
        $constraint = (rmpseudoknots($sequence, \@constraint, \&Core::Mathematics::mean))[0];

    }

    ######### End of partition function block ##########

    # Build dotplot of pairing probabilities (only if necessary)
    if ($windowed || $dp || $sh) {

        foreach my $pair (keys %pairprobs) {

            my ($i, $j) = split(/-/, $pair);
            push(@dotplot, [$i, $j, mean(@{$pairprobs{$pair}->{pairs}})]);

        }

        @dotplot = sort {$a->[0] <=> $b->[0] || $a->[1] <=> $b->[1]} @dotplot;

    }

    if ($dp) {

        open(my $wh, ">", $output . "dotplot/" . $id . ".dp") or _logError("Unable to write output dotplot file ($!)", $id) and next;
        select((select($wh), $|=1)[0]);
        print $wh length($sequence) . "\ni\tj\t-log10(Probability)\n";
        print $wh join("\t", $_->[0] + 1, $_->[1] + 1, -logarithm($_->[2], 10)) . "\n" for (@dotplot);
        close($wh);

    }

    return if ($dotplotOnly);

    ######### MFE folding block ##########

    $nwins = $windowed ? max(1, floor(($length - $foldwin) / $foldstep) + 1) : 1;

    # Define windows to be analyzed
    # Windows are defined as: [start, end]
    if ($windowed && $nwins > 1) {

        for (100, 50, -50, -100) {

            my $end = $foldwin - $_ - 1;
            $end = $length - 1 if ($end > $length - 1);

            next if ($end < MINWINLENGTH - 1);

            push(@foldwins, [0, $end]);

            last if ($end == $length - 1);

        }

        for (0 .. $nwins - 1) {

            my ($start, $end);
            $start = $foldstep * $_;
            $end = $_ == $nwins - 1 ? $length - 1 : $start + $foldwin - 1;
            push(@foldwins, [$start, $end]);

        }

        for (100, 50, -50, -100) {

            my $start = $length - $foldwin + $_;
            $start = 0 if (!ispositive($start));

            next if ($start > $length - MINWINLENGTH);

            push(@foldwins, [$start, $length - 1]);

            last if (!$start);

        }

    }
    else { push(@foldwins, [0, $length - 1]); }

    foreach my $sampleSet (samplesets(scalar(@reactivity))) {

        # Calculate MFE folding in window
        foreach my $win (@foldwins) {

            my ($start, $end, $subseq, $winconst, $tmpfile,
                $alignment, $alnReact, @winreact);
            ($start, $end) = @{$win};

            $subseq = substr($sequence, $start, $end - $start + 1);
            $winconst = rmlonelypairs(fixdotbracket(substr($constraint, $start, $end - $start + 1)));
            $tmpfile = "$escapedId\_$start-$end\_$$\_" . ($useAlifold ? "ali" : $sampleSet->[0]);

            if ($useAlifold) {

                ($alignment, $alnReact) = winalignment($start, $end, \@sequences, \@reactivity);
                @winreact = @meanReact[$start .. $end]; # Hard constraints are derived from the mean reactivity profile

            }
            else { @winreact = @{$reactivity[$sampleSet->[0]]}[$start .. $end]; }

            my $win = winfold( file          => $tmpfile,
                               sample        => $sampleSet->[0],
                               sequence      => $subseq,
                               reactivity    => \@winreact,
                               alignment     => $alignment,
                               alnReactivity => $alnReact,
                               pos           => $start,
                               constraint    => $winconst );

            if (!$win) {

                unlink(glob($tmpfile . "*"));

                return({ foldErr => 1 });

            }

            $npairs{$_}->{pairs}++ for (keys %{$win});
            push(@windows, [$start, $end]);

        }

    }

    # Build structure with pairs in > 0.50 of windows
    foreach my $pair (sort {(split(/-/, $a))[0] <=> (split(/-/, $b))[0]} keys %npairs) {

        my ($i, $j, $n1, $n2);
        ($i, $j) = split(/-/, $pair);

        for (@windows) { 
            
            $n1++ if (inrange($i, $_));
            $n2++ if (inrange($j, $_)); 
            
        }

        # If pair(i,j) is present in > 0.50 windows
        if ($npairs{$pair}->{pairs} / $n1 > 0.5 &&
            $npairs{$pair}->{pairs} / $n2 > 0.5) { push(@keptpairs, [$i, $j]); }

    }

    $entry = Data::Sequence::Structure->new( id          => $id,
                                             sequence    => $sequence,
                                             basepairs   => \@keptpairs,
                                             lonelypairs => $nlp ? 0 : 1 );

    undef(@foldwins);

    ######### End of MFE folding block ##########

    ######### Pseudoknots detection block ##########

    if ($pk) { # Pseudoknot detection

        my (@keptpkpairs, %pkpairs, %basepairs, %pkbases);
        $nwins = $windowed ? max(1, floor(($length - $pkwin) / $pkstep) + 1) : 1;

        # Define windows to be analyzed
        # Windows are defined as: [start, end]
        if ($windowed && $nwins > 1) {

            for (50, 100) {

                my $end = $pkwin - $_ - 1;

                last if ($end < MINWINLENGTH - 1);

                push(@foldwins, [0, $end]);

            }

            for (0 .. $nwins - 1) {

                my ($start, $end);
                $start = $pkstep * $_;
                $end = $_ == $nwins - 1 ? $length - 1 : $start + $pkwin - 1;
                push(@foldwins, [$start, $end]);

            }

            for (50, 100) {

                my $start = $length - $pkwin + $_;

                last if ($start > $length - MINWINLENGTH);

                push(@foldwins, [$start, $length - 1]);

            }

        }
        else { push(@foldwins, [0, $length - 1]); }

        for my $i (0 .. $#reactivity) {

            # Predict pseudoknot in window
            foreach my $win (@foldwins) {

                my ($start, $end, $subseq, $tmpfile,
                    $winconst, @winreact);
                ($start, $end) = @{$win};

                $subseq = substr($sequence, $start, $end - $start + 1);
                $winconst = fixdotbracket(substr($entry->structure(), $start, $end - $start + 1));
                @winreact = @{$reactivity[$i]}[$start .. $end];
                $tmpfile = "$escapedId\_$start-$end\_$$\_$i";

                my ($winpk, $winbp) = pkdetect( file       => $tmpfile,
                                                sample     => $i,
                                                sequence   => $subseq,
                                                reactivity => \@winreact,
                                                pos        => $start,
                                                constraint => $winconst );

                if (!$winpk) {

                    unlink(glob($tmpfile . "*"));

                    return({ pkErr => 1 });

                }

                $pkpairs{$_}++ for (keys %{$winpk});
                $basepairs{$_}++ for (keys %{$winbp});
                push(@windows, [$start, $end]);

            }

        }

        # Discard pseudoknotted pairs in < 0.50 of windows
        foreach my $pair (keys %pkpairs) {

            my ($i, $j, $n1, $n2);
            ($i, $j) = split(/-/, $pair);

            for (@windows) { 
                
                $n1++ if (inrange($i, $_));
                $n2++ if (inrange($j, $_)); 
                
            }

            # We essentially check whether the pseudoknot base-pair exists in other windows,
            # even if it has not been marked as a pseudoknot
            $pkpairs{$pair} += $basepairs{$pair} if (exists $basepairs{$pair});

            # If pseudoknotted pair is present in > 0.50 windows
            push(@pkpairs, [$i, $j]) if ($pkpairs{$pair} / $n1 > 0.5 && $pkpairs{$pair} / $n2 > 0.5);

        }

        # In the pkdetect() function, we have already discarded the pk helices that overlap by more than 50%
        # with helices in the MFE structure. Now, we remove the pairs from the MFE structure that still overlap (if any)
        %pkbases = map { $_ => 1 } map { @$_ } @pkpairs;
        @keptpairs = grep { !exists $pkbases{$_->[0]} && !exists $pkbases{$_->[1]} } @keptpairs;

        # Shouldn't be necessary anymore
        # In case we have incompatible base-pairs, we let rmpseudoknots() resolve that
        #@keptpkpairs = rmpseudoknots($sequence, [ @keptpairs, @pkpairs ]);
        #@pkpairs = @{$keptpkpairs[1]} ? @{$keptpkpairs[1]} : listpairs($keptpkpairs[0]);

        undef($nwins);
        undef(@foldwins);
        undef(@windows);

    }

    ######### End of pseudoknots detection block ##########

    $entry = Data::Sequence::Structure->new( id          => $id,
                                             sequence    => $sequence,
                                             basepairs   => [ @keptpairs, @pkpairs ],
                                             pseudoknots => 1,
                                             lonelypairs => $nlp ? 0 : 1 );
    $energy = RNA::energy_of_struct($sequence, $entry->structure());
    $energy += pseudodG($entry->structure(), \@meanReact) if (!$ignore);
    $entry->energy(sprintf("%.2f", min(0, $energy)));
    $file = $output . "structures/" . $id . "." . ($ct ? "ct" : "db");

    # Fix the ID changed by ViennaRNA
    $entry->id($id) if ($method == 1);

    eval { $io = Data::IO::Sequence->new( file   => $file,
                                          mode   => "w",
                                          format => $ct ? "CT" : "Vienna" ); };

    if ($@) {

        _logError($@, $id);

        return({ foldErr => 1 });

        next;

    }

    $io->write($entry);
    $io->close();

    if ($sh) {

        my $i = 0;

        open(my $wh, ">", $output . "shannon/" . $id . ".wig") or _logError("Unable to write Shannon entropy WIG file ($!)", $id) and next;
        select((select($wh), $|=1)[0]);

        print $wh "track type=wiggle_0 name=\"" . $id . " Shannon entropy\"\n" .
                    "variableStep chrom=" . $id . "\n";

        for (@shannon) {

            $i++;
            print $wh $i . " " . $_ . "\n";

        }

        close($wh);

    }

    if ($img) {

        my ($hasReact, $winSize, $rnaPlot, @plots, @heights, %legend);
        $winSize = 20;
        $plot->file("${output}plots/summaries/$id.pdf");

        $hasReact = scalar(grep { isnumeric($_) } @meanReact);

        if ($hasReact) {

            my ($reactPlot, %dataLabels);
            %dataLabels = ("Reactivity" => [ map { $_ <= 0.4 || isnan($_) ? "0-0.4" : ($_ < 0.7 ? "0.4-0.7" : "0.7+") } @meanReact ]);
            
            if (@reactivity > 1) {
                
                my @meanStdev = map { my $i = $_; [ map { $_->[$i] } @reactivity ] } 0 .. $#{$reactivity[0]};
                @meanStdev = map { isnumeric(@$_) ? stdev(@$_) : "NaN" } @meanStdev;
                $dataLabels{"sd"} = \@meanStdev;

            }

            $reactPlot = Graphics::Chart::Barplot->new( data            => \@meanReact,
                                                        dataLabels      => \%dataLabels,
                                                        fill            => "Reactivity",
                                                        legendTitle     => "Reactivity",
                                                        legendColors    => { "0-0.4"   => "#000000",
                                                                             "0.4-0.7" => "#fdb934",
                                                                             "0.7+"    => "#aa1a1f" },
                                                        legendKeyWidth  => 14,
                                                        legendKeyHeight => 14,
                                                        legendTextSize  => 12,
                                                        legendTitleSize => 14,
                                                        axisTitleSize   => 14,
                                                        labelTextSize   => 12,
                                                        legendSort      => [ "0.7+", "0.4-0.7", "0-0.4" ],
                                                        yTitle          => "Reactivity",
                                                        yLimit          => [0, 2],
                                                        xLimit          => [0, $length + 1],
                                                        stdev           => @reactivity > 1 ? "sd" : undef );

            push(@plots, $reactPlot);
            push(@heights, 1);

            if (2 * $winSize + 1 < $length) {

                my ($medianReact, $medianReactPlot, @winReact);
                @winReact = ("NaN") x $length;
                $medianReact = median(grep { isnumeric($_) } @meanReact);

                for(my $i = $winSize; $i < @meanReact - $winSize; $i += 1) {

                    my (@win, @i);
                    @win = @meanReact[$i - $winSize .. $i + $winSize];
                    @i = grep { isnumeric($win[$_]) } 0 .. $#win;

                    if (@i / @win >= 0.25) {

                        @win = @win[@i];
                        $winReact[$i] = median(@win) - $medianReact;

                    }

                }

                $medianReactPlot = Graphics::Chart::Area->new( data            => \@winReact,
                                                               dataLabels      => { "fill" => [ ("black") x @winReact ] },
                                                               fill            => "fill",
                                                               legendColors    => { "black" => "#000000" },
                                                               lineThickness   => 0,
                                                               legend          => 0,
                                                               legendTextSize  => 12,
                                                               legendTitleSize => 14,
                                                               axisTitleSize   => 14,
                                                               labelTextSize   => 12,
                                                               xLimit          => [0, $length + 1],
                                                               yTitle          => "Reactivity\n(Window median - Transcript median)" );

                push(@plots, $medianReactPlot);
                push(@heights, 1);

            }

        }

        if (@keptpairs) {

            my ($maxBp, $structPlot);
            $maxBp = max(map { abs(diff(@$_)) } (@keptpairs, @pkpairs));
            $structPlot = Graphics::Chart::Arcs->new( data            => [ map { [ map { $_ + 1 } @$_ ] } (@keptpairs, @pkpairs) ],
                                                      dataLabels      => { "fill" => [ ("black") x (@keptpairs + @pkpairs) ] },
                                                      flip            => "down",
                                                      fill            => "fill",
                                                      legendColors    => { "black" => "#000000" },
                                                      xLimit          => [0, $length + 1],
                                                      lineThickness   => 0.2,
                                                      legend          => 0,
                                                      background      => 0,
                                                      grid            => 0,
                                                      xTicks          => 0,
                                                      yTicks          => 0,
                                                      xLabels         => 0,
                                                      yLabels         => 0,
                                                      legendTextSize  => 12,
                                                      legendTitleSize => 14,
                                                      axisTitleSize   => 14,
                                                      labelTextSize   => 12 );

            push(@plots, $structPlot);
            push(@heights, $maxBp / $length * ($windowed || $dp || $sh ? 1.75 : 1.5));

        }

        if ($windowed || $dp || $sh) {

            if (2 * $winSize + 1 < $length) {

                my ($medianShan, $shanPlot, @winShannon);
                @winShannon = ("NaN") x $length;
                $medianShan = median(grep { isnumeric($_) } @shannon);

                for(my $i = $winSize; $i < @shannon - $winSize; $i += 1) {

                    my (@win, @i);
                    @win = @shannon[$i - $winSize .. $i + $winSize];
                    @i = grep { isnumeric($win[$_]) } 0 .. $#win;

                    if (@i / @win >= 0.25) {

                        @win = @win[@i];
                        $winShannon[$i] = median(@win) - $medianShan;

                    }

                }

                $shanPlot = Graphics::Chart::Area->new( data            => \@winShannon,
                                                        dataLabels      => { "fill" => [ ("brown") x @winShannon ] },
                                                        fill            => "fill",
                                                        legendColors    => { "brown" => "#916A28" },
                                                        lineThickness   => 0.2,
                                                        legend          => 0,
                                                        xLimit          => [0, $length + 1],
                                                        yTitle          => "Shannon entropy\n(Window median - Transcript median)",
                                                        legendTextSize  => 12,
                                                        legendTitleSize => 14,
                                                        axisTitleSize   => 14,
                                                        labelTextSize   => 12 );

                push(@plots, $shanPlot);
                push(@heights, 1);

            }

            @dotplot = grep { $_->[2] >= 0.05 } sort { $a->[2] <=> $b->[2] } @dotplot;

            if (@dotplot) {

                my ($maxBp, $probsPlot, @data);
                @data = map { [ @$_[0, 1] ] } @dotplot;
                $maxBp = max(map { abs(diff(@$_)) } @data);
                $probsPlot = Graphics::Chart::Arcs->new( data            => \@data,
                                                         dataLabels      => { "prob" => [ map { $_->[2] < 0.1 ? "5-10%" : ($_->[2] < 0.4 ? "10-40%" : ($_->[2] < 0.7 ? "40-70%" : "70-100%")) } @dotplot ] },
                                                         flip            => "down",
                                                         fill            => "prob",
                                                         legendKeyWidth  => 14,
                                                         legendKeyHeight => 4,
                                                         legendColors    => { "5-10%"   => "#979797",
                                                                             "10-40%"  => "#ffcc07",
                                                                             "40-70%"  => "#478ecc",
                                                                             "70-100%" => "#58b14a" },
                                                         legendSort      => [ "5-10%", "10-40%", "40-70%", "70-100%" ],
                                                         xLimit          => [0, $length + 1],
                                                         lineThickness   => 0.2,
                                                         legendTitle     => "Probability",
                                                         legend          => 1,
                                                         background      => 0,
                                                         grid            => 0,
                                                         xTicks          => 0,
                                                         yTicks          => 0,
                                                         xLabels         => 0,
                                                         yLabels         => 0,
                                                         legendTextSize  => 12,
                                                         legendTitleSize => 14,
                                                         axisTitleSize   => 14,
                                                         labelTextSize   => 12 );

                push(@plots, $probsPlot);
                push(@heights, $maxBp / $length * 1.75);

            }

        }

        $plot->plot([ map { [$_] } @plots ], { heights => \@heights});

        if (!$gradient) {

            @meanReact = map { isnan($_) ? $_ : ($_ <= 0.4 ? "0 - 0.4" : ($_ <= 0.8 ? "0.4 - 0.8" : "0.8+")) } @meanReact;
            %legend = ( "0 - 0.4"   => "#000000", 
                        "0.4 - 0.8" => "#FFCD2F", 
                        "0.8+"      => "#9A2322" );

        }
        else {

            %legend = ( "0"   => "#FFFFFF", 
                        "0.4" => "#FFCD2F", 
                        "0.8" => "#9A2322" );

        }

        # Generate secondary structure plot
        $rnaPlot = Graphics::RNA->new( file        => "${output}plots/structures/$id.svg",
                                       sequence    => $entry->sequence(),
                                       structure   => $entry->structure(), 
                                       values      => \@meanReact,
                                       plotLegend  => 1, 
                                       legend      => \%legend,
                                       legendScale => 1.5,
                                       legendTitle => "Reactivity",
                                       colorScale  => $gradient ? "gradient" : "discrete" );
        $rnaPlot->plot();

    }

}

##
#  Pseudoknots detection
#
#  pkdetect( file       => <temporary file prefix>,
#            sequence   => <window sequence>,
#            reactivity => <ARRAY ref of win reactivities>,
#            pos        => <true position> );
##

sub pkdetect {

    my %parameters = @_;

    my ($ret, $fastafile, $cmd, $energy, 
        $pktollerance, $mfetollerance, $io, $count,
        @subopt, @results, @hardconstraint, @pkpairs,
        @basepairs, @pkreact, %pkpairs, %mfepaired,
        %candidatepk, %basepairs);
    $count = 0;
    $fastafile = $parameters{file} . ".fa";
    @hardconstraint = @{hardbases(@{$parameters{reactivity}})} if (!$ignore && $hc);

    foreach my $helix (@{(listhelices($parameters{constraint}, 1))[0]}) { %mfepaired = (%mfepaired, map { $_ => 1 } (@{$helix->{h5bases}}, @{$helix->{h3bases}})); }

    if ($pkmethod == 1) {

        my (@sampled);

        $energy = RNA::energy_of_struct($parameters{sequence}, $parameters{constraint});
        $energy += pseudodG($parameters{constraint}, $parameters{reactivity}) if (!$ignore);
        #$mfetollerance = $energy - ($energy * ($entry->length() > 100 ? 0.25 : 0.2));
        #$mfetollerance = $energy - ($energy * 0.25); # We use increased tollerance due to some differences between RNAfold and ShapeKnots energy calculation
        $mfetollerance = $energy - ($energy * (length($parameters{constraint}) > 100 ? max($pkpercentmfe, 0.25) : $pkpercentmfe));
        $pktollerance = $energy - ($energy * 0.1);

        # Note: suboptimal structures are sampled without SHAPE-derived pseudo-energies,
        # as the pseudo free energy contribution is added afterwards by pseudodG()
        eval { @sampled = $vrna->subopt($parameters{sequence}, { %{$vrfoldparams[$parameters{sample}]},
                                                                 constraint  => $ignore ? undef : vrconstraint(\@hardconstraint, $parameters{constraint}),
                                                                 reactivity  => undef,
                                                                 sampleSize  => $pksubopt,
                                                                 deltaEnergy => $suboptDeltaEnergy }); };

        if ($@) {

            _logError($@, $parameters{file});

            return;

        }

        foreach my $struct (@sampled) {

            my $suboptenergy = RNA::energy_of_struct($parameters{sequence}, $struct->{structure});
            $suboptenergy += pseudodG($struct->{structure}, $parameters{reactivity}) if (!$ignore);

            # Here we compare the suboptimal energy and the mfetollerance down to the 1st decimal digit only
            # to avoid being too stringent (in some cases the optimal helix is within less than 0.1kcal/mol from the tollerance)
            push(@subopt, [$struct->{structure}, $suboptenergy]) if (sprintf("%.1f", $suboptenergy) <= sprintf("%.1f", $mfetollerance));

        }

        if (!$nozuker) {

            my $structs = 0;

            undef(@sampled);

            # I added this to sample suboptimal structures according to Zuker (increases the folding space searched)
            eval { @sampled = $vrna->subopt($parameters{sequence}, { temperature => $t,
                                                                     sampleSize  => 0,
                                                                     zuker       => 1 }); };

            if ($@) {

                _logError($@, $parameters{file});

                return;

            }

            foreach my $struct (@sampled) {

                last if ($structs == $zukersubopt);

                my $suboptenergy = RNA::energy_of_struct($parameters{sequence}, $struct->{structure});
                $suboptenergy += pseudodG($struct->{structure}, $parameters{reactivity}) if (!$ignore);

                # Here we compare the suboptimal energy and the mfetollerance down to the 1st decimal digit only
                # to avoid being too stringent (in some cases the optimal helix is within less than 0.1kcal/mol from the tollerance)
                push(@subopt, [$struct->{structure}, $suboptenergy]) if (sprintf("%.1f", $suboptenergy) <= sprintf("%.1f", $mfetollerance));

                $structs++;

            }

        }

        @subopt = uniq(@subopt);

        foreach my $subopt (@subopt) {

            foreach my $helix (@{(listhelices($subopt->[0], $nlp ? 0 : 1))[0]}) {

                my ($pairedinmfe, $h5, $h3, $tmpsubopt,
                    $tmpenergy, @tmpsubopt, @pkbases);
                @tmpsubopt = split(//, $subopt->[0]);
                @pkbases = (@{$helix->{h5bases}}, @{$helix->{h3bases}});
                $pairedinmfe += $mfepaired{$_} for (@pkbases);

                next if ($pairedinmfe / @pkbases > 0.5);

                $tmpsubopt[$_] = "." for (@pkbases);
                $tmpsubopt = join("", @tmpsubopt);
                $tmpenergy = RNA::energy_of_struct($parameters{sequence}, $tmpsubopt);
                $tmpenergy += pseudodG($tmpsubopt, $parameters{reactivity}) if (!$ignore);

                $candidatepk{join("-", @pkbases)} = { bases  => [$helix->{h5bases}, $helix->{h3bases}],
                                                      energy => $subopt->[1] - $tmpenergy };

                # Patched (instead of calculating the energy gain by ourselves, now we just let ViennaRNA do the dirty job):
                #$h5 = substr($parameters{sequence}, $helix->{h5bases}->[-1], $helix->{h5bases}->[0] - $helix->{h5bases}->[-1] + 1);
                #$h3 = substr($parameters{sequence}, $helix->{h3bases}->[0], $helix->{h3bases}->[-1] - $helix->{h3bases}->[0] + 1);

                #$candidatepk{join("-", @pkbases)} = { bases  => [$helix->{h5bases}, $helix->{h3bases}],
                #								      energy => RNA::duplexfold($h5, $h3)->{energy} };

            }

        }

        undef(@subopt);

        foreach my $pk (sort {$candidatepk{$a}->{energy} <=> $candidatepk{$b}->{energy}} keys %candidatepk) {

            last if ($count == $pkmax);

            my (@pkbases, @sampled, %seen);
            @pkbases = grep { !exists $mfepaired{$_} } (@{$candidatepk{$pk}->{bases}->[0]}, @{$candidatepk{$pk}->{bases}->[1]});

            eval { @sampled = $vrna->subopt($parameters{sequence}, { %{$vrfoldparams[$parameters{sample}]},
                                                                     constraint => $ignore ? undef : vrconstraint([@hardconstraint, @pkbases], $parameters{constraint}),
                                                                     reactivity => $ignore ? undef : $parameters{reactivity},
                                                                     sampleSize => 100 }); };

            if ($@) {

                _logError($@, $parameters{file});

                return;

            }

            foreach my $struct (@sampled) {

                my ($line, $suboptenergy);
                $line = $struct->{structure};

                next if (exists $seen{$line});

                $seen{$line} = 1;
                $suboptenergy = RNA::energy_of_struct($parameters{sequence}, $line);
                $suboptenergy += pseudodG($line, $parameters{reactivity}) if (!$ignore);

                if ($suboptenergy < $mfetollerance) {

                    my @pairs = listpairs($line);

                    for (0 .. $#{$candidatepk{$pk}->{bases}->[0]}) {

                        my ($i, $j) = ($candidatepk{$pk}->{bases}->[0]->[$_], $candidatepk{$pk}->{bases}->[1]->[$_]);
                        push(@pairs, [$i, $j]);

                    }

                    push(@subopt, { pairs   => \@pairs,
                                    energy  => $suboptenergy + $candidatepk{$pk}->{energy} });

                }

            }

            $count++;

        }

        foreach my $struct (sort {$a->{energy} <=> $b->{energy}} @subopt) {

            my ($dotbracket, $pkpairs, $pkregion, $noilregion,
                $ss, $ne, $il, $pkenergy, $ildb, @il);
            ($dotbracket, $pkpairs) = rmpseudoknots($parameters{sequence}, [ @{$struct->{pairs}} ]);
            ($ss, $ne, $il) = (0, 0, 0);

            next if (@{$pkpairs} <= 1);

            $pkregion = substr($dotbracket, $pkpairs->[0]->[0], $pkpairs->[0]->[1] - $pkpairs->[0]->[0] + 1);
            ($ss) = $pkregion =~ tr/././;
            $noilregion = fixdotbracket($pkregion);
            $ne = scalar(@{(listhelices($noilregion))[0]});
            $ildb = join("", map { substr($pkregion, $_, 1) eq substr($noilregion, $_, 1) ? "." : substr($pkregion, $_, 1) } 0 .. length($pkregion) - 1);

            while ($ildb =~ m/((?:\(+.?\(*)|(?:\)+.?\)*))/g) {

                my ($helix, $length);
                $helix = $1;
                ($length) = $helix =~ tr/()/()/;
                push(@il, $length);

            }

            $il += $_ * ($lambda{($_ < 2 ? 2 : ($_ > 15 ? 15 : $_))} ** 2) for (@il);
            $il ||= 1;

            $pkenergy = $struct->{energy} + ($p1 * ($ss == 0 && $ne == 0 ? 1 : logarithm((6.5 ** 2) * $ss + (15 ** 2) * $ne, Core::Mathematics::e)) + $p2 * logarithm($il, Core::Mathematics::e));

            push(@results, { pkpairs => $pkpairs,
                             energy  => $pkenergy }) if ($pkenergy < $pktollerance);

        }

        if (@results) {

            my $result = (sort {$a->{energy} <=> $b->{energy}} @results)[0];
            @pkpairs = @{$result->{pkpairs}};

        }


    }
    else {

        my $skparams = $skparams[$parameters{sample}];

        # Writes temporary constraint file and sequence
        return unless(seq2fasta($parameters{file}, $parameters{file}, $parameters{sequence}));
        if (!$ignore) { return unless(xml2shape($parameters{file}, $parameters{reactivity})); }
        return unless(applyhardconstraint($parameters{file}, \@hardconstraint, $parameters{constraint}));

        $cmd = "'" . $fastafile . "' '" . $parameters{file} . ".ct' -p1 " . $p1 . " -p2 " . $p2 . " -ip " . ($pkpercentmfe * 100) . " -ph " . $pkmax;
        $cmd .= " -sh '" . $parameters{file} . ".shape'" if (!$ignore);
        #$cmd .= " -C '" . $parameters{file} . ".constraint'" if ($hc); # ShapeKnots takes -C parameter, but it is ignored
        $ret = `$shapeknots $cmd $skparams 2>&1`;

        if ($ret !~ m/#+\s+DONE\s+#+/m) {

            _logError("Exception in RNAstructure ShapeKnots\n\n$ret", $parameters{file});

            return;

        }

        # If the sequence is all single-stranded, ShapeKnots may return an empty CT file
        if (-s $parameters{file} . ".ct") {

            eval { $io = Data::IO::Sequence->new( file        => $parameters{file} . ".ct",
                                                  pseudoknots => 1,
                                                  lonelypairs => $nlp ? 0 : 1 ); };

            if ($@) { 

                _logError($@, $parameters{file});

                return; 
                
            }

            if (my $entry = $io->read()) {

                @pkpairs = $entry->pkpairs();
                @basepairs = $entry->basepairs();

            }
            else { 

                _logError("ShapeKnots returned an empty structure", $parameters{file});

                return; 
                
            }

        }

    }

    if (my @mapPkPairs = map {@$_} @pkpairs) {

        my $pairedinmfe = 0; 
        $pairedinmfe += $mfepaired{$_} for (@mapPkPairs);

        return({}, {}) if ($pairedinmfe / scalar(@mapPkPairs) > 0.5);

    }

    @pkreact = grep { !isnan($_) } @{$parameters{reactivity}}[sort {$a <=> $b} map {@$_} @pkpairs];

    if ($ignore || !@pkreact || mean(@pkreact) <= $pkcutoff) {

        $pkpairs{($_->[0] + $parameters{pos}) . "-" . ($_->[1] + $parameters{pos})} = 1 for (@pkpairs);
        $basepairs{($_->[0] + $parameters{pos}) . "-" . ($_->[1] + $parameters{pos})} = 1 for (@basepairs);

    }

    unlink(glob($parameters{file} . "*"));

    return(\%pkpairs, \%basepairs);

}

##
#  Windowed folding function
#
#  winfold( file       => <temporary file prefix>,
#           sequence   => <window sequence>,
#           reactivity => <ARRAY ref of window reactivities>,
#           pos        => <true position>,
#           constraint => <base-pair constraint from windowed partition>,
#           pseudoknot => <ARRAY of indexes of pseudoknotted bases> );
##

sub winfold {

    my %parameters = @_;

    my ($ret, $cmd, $io, $entry,
        @hardconstraint, %pairs);
    @hardconstraint = @{hardbases(@{$parameters{reactivity}})} if (!$ignore && $hc);

    if ($method == 1) { # ViennaRNA

        my ($foldParams);
        $foldParams = { %{$vrfoldparams[$parameters{sample}]},
                        constraint        => vrconstraint(\@hardconstraint, $parameters{constraint}),
                        enforceConstraint => 1 };

        if ($useAlifold) { $foldParams->{reactivity} = $ignore ? undef : $parameters{alnReactivity}; }
        else { $foldParams->{reactivity} = $ignore ? undef : $parameters{reactivity}; }

        eval { $entry = $useAlifold ? $vrna->alifold($parameters{alignment}, $foldParams)
                                    : $vrna->fold($parameters{sequence}, $foldParams); };

        # RNAfold/RNAalifold fail to backtrack when hard constraints cannot be
        # enforced, so in that case we retry without enforcing them
        if ($@ || !defined $entry) {

            $foldParams->{enforceConstraint} = 0;

            eval { $entry = $useAlifold ? $vrna->alifold($parameters{alignment}, $foldParams)
                                        : $vrna->fold($parameters{sequence}, $foldParams); };

            if ($@) {

                _logError($@, $parameters{file});

                return;

            }

        }

    }
    else { # RNAstructure

        my $rsfoldparams = $rsfoldparams[$parameters{sample}];

        # Writes temporary constraint file and sequence
        return unless(seq2fasta($parameters{file}, $parameters{file}, $parameters{sequence}));
        if (!$ignore) { return unless(xml2shape($parameters{file}, $parameters{reactivity})); }
        return unless(applyhardconstraint($parameters{file}, \@hardconstraint, $parameters{constraint}));

        $cmd = $rnastructure . " '" . $parameters{file} . ".fa' '" . $parameters{file} . ".ct' -C '" . $parameters{file} . ".constraint'";
        $cmd .= " -sh '" . $parameters{file} . ".shape'" if (!$ignore);

        $ret = `$cmd $rsfoldparams 2>&1`;

        if ($ret !~ m/Single strand folding complete\./) {

            _logError("Exception in RNAstructure Fold\n\n$ret", $parameters{file});

            return;

        }

        eval { $io = Data::IO::Sequence->new( file        => $parameters{file} . ".ct",
                                              lonelypairs => $nlp ? 0 : 1 ); };

        if ($@) { 

            _logError($@, $parameters{file});

            return; 
            
        }

        $entry = $io->read();

    }

    if (defined $entry) {

        foreach my $pair ($entry->basepairs()) { # i,j numbering is 0-based

            $_ += $parameters{pos} for (@{$pair});
            $pairs{$pair->[0] . "-" . $pair->[1]} = 1;

        }

    }
    else { 

        _logError("Fold returned an empty structure", $parameters{file});

        return; 
        
    }

    unlink(glob($parameters{file} . "*"));

    return(\%pairs);

}

##
#  Windowed partition function calculation
#
#  partition( file       => <temporary file prefix>,
#             sequence   => <window sequence>,
#             reactivity => <ARRAY ref of window reactivities>,
#             pos        => <true position>,
#             trim5      => <n bases to trim from win 5'>,
#             trim3      => <n bases to trim from win 3'> );
##

sub partition {

    my %parameters = @_;

    my ($ret, $cmd, @hardconstraint, %pairs,
        %shannon);
    @hardconstraint = @{hardbases(@{$parameters{reactivity}})} if (!$ignore && $hc);

    # Fill window's Shannon entropy with 0
    %shannon = map { $_ => 0 } ($parameters{pos} + $parameters{trim5} .. $parameters{pos} + length($parameters{sequence}) - $parameters{trim3} - 1);

    if ($method == 1) { # ViennaRNA

        my ($entry, $foldParams, $bpprobs);
        $foldParams = { %{$vrfoldparams[$parameters{sample}]},
                        constraint        => $ignore ? undef : vrconstraint(\@hardconstraint, $parameters{constraint}),
                        enforceConstraint => $ignore ? 0 : 1,
                        partitionFunction => 1 };

        if ($useAlifold) { $foldParams->{reactivity} = $ignore ? undef : $parameters{alnReactivity}; }
        else { $foldParams->{reactivity} = $ignore ? undef : $parameters{reactivity}; }

        eval { $entry = $useAlifold ? $vrna->alifold($parameters{alignment}, $foldParams)
                                    : $vrna->fold($parameters{sequence}, $foldParams); };

        if ($@ || !defined $entry) {

            _logError($@ || ($useAlifold ? "RNAalifold" : "RNAfold") . " returned an empty structure", $parameters{file});

            return;

        }

        $bpprobs = $entry->bpprobability();

        foreach my $i0 (sort {$a <=> $b} keys %{$bpprobs}) {

            foreach my $j0 (sort {$a <=> $b} keys %{$bpprobs->{$i0}}) {

                my ($i, $j, $p, $entropy);

                next if ($j0 <= $i0); # Base-pairing probabilities are stored symmetrically

                ($i, $j) = ($i0 + 1, $j0 + 1); # Base-pairing probabilities are 0-based
                $p = $bpprobs->{$i0}->{$j0};

                next if ($i <= $parameters{trim5} || $j >= length($parameters{sequence}) - $parameters{trim3});

                $i += $parameters{pos} - 1;         # Base numbering is 1-based
                $j += $parameters{pos} - 1;

                $pairs{"$i-$j"} = $p;
                $entropy = $p ? $p * log($p) / log(10) : 0;
                $shannon{$i} += $entropy;
                $shannon{$j} += $entropy;

            }

        }

    }
    else { # RNAstructure

        my $rsfoldparams = $rsfoldparams[$parameters{sample}];

        # Writes temporary constraint file and sequence
        return unless(seq2fasta($parameters{file}, $parameters{file}, $parameters{sequence}));
        if (!$ignore) { return unless(xml2shape($parameters{file}, $parameters{reactivity})); }
        return unless(applyhardconstraint($parameters{file}, \@hardconstraint, $parameters{constraint}));

        $cmd = $partition . " '" . $parameters{file} . ".fa' '" . $parameters{file} . ".pfs'";
        $cmd .= " -sh '" . $parameters{file} . ".shape' -C '" . $parameters{file} . ".constraint'" if (!$ignore);

        $ret = `$cmd $rsfoldparams 2>&1`;

        if ($ret !~ /Single strand partition function complete\./) {

            _logError("Exception in RNAstructure partition\n\n$ret", $parameters{file});

            return;

        }

        $cmd = $probplot . " -t '" . $parameters{file} . ".pfs' '" . $parameters{file} . ".pfs.txt'";
        $ret = `$cmd 2>&1`;

        if ($ret !~ /Probability dot plot complete\./) {

            _logError("Exception in RNAstructure ProbabilityPlot\n\n$ret", $parameters{file});

            return;

        }

        open(my $fh, "<", $parameters{file} . ".pfs.txt") or _logError("Unable to read from temporary PFS file ($!)", $parameters{file}) and return;
        while(<$fh>) {

            if ($_ =~ m/^(\d+)\t(\d+)\t([\d\.]+(?:e-\d+)?)$/) {

                my ($i, $j, $p, $logp, $entropy);
                ($i, $j, $logp) = ($1, $2, $3);

                next if ($i <= $parameters{trim5} || $j >= length($parameters{sequence}) - $parameters{trim3});

                $i += $parameters{pos} - 1;         # Base numbering is 1-based
                $j += $parameters{pos} - 1;
                $p = 10 ** (-$logp);


                $pairs{"$i-$j"} = $p;
                $entropy = -$logp * $p;
                $shannon{$i} += $entropy;
                $shannon{$j} += $entropy;

            }

        }
        close($fh);

    }

    $_ = -$_ for (values %shannon);

    unlink(glob($parameters{file} . "*"));

    return(\%pairs, \%shannon);

}

sub xml2shape {

    my ($file, $reactivity) = @_;

    my $i = 0;

    open(my $wh, ">", $file . ".shape") or _logError("Unable to write temporary SHAPE file ($!)") and return();
    select((select($wh), $|=1)[0]);
    for (@{$reactivity}) {

        $i++;

        #next if (isnan($_));

        print $wh $i . " " . (isnan($_) ? -999 : $_) . "\n";

    }
    close($wh);

    return(1);

}

##
#  Returns the sets of samples to be folded together. The canonical RNAfold
#  approach folds (and calculates the partition function for) one sample at a
#  time, then aggregates the results; the RNAalifold approach instead takes all
#  the reactivity profiles at once, so it yields a single set
##

sub samplesets {

    my $samples = shift;

    return($useAlifold ? [0 .. $samples - 1] : map { [$_] } 0 .. $samples - 1);

}

##
#  Builds the input for Interface::ViennaRNA's alifold(): an "alignment" made of
#  each sample's window sequence, and the matching reactivity profiles. Sample IDs
#  are zero-padded, as sequences and SHAPE files are paired by sorting the IDs
##

sub winalignment {

    my ($start, $end, $sequences, $reactivity) = @_;

    my (%alignment, %alnReactivity);

    for my $i (0 .. $#{$sequences}) {

        my $sampleId = sprintf("%05d", $i + 1);
        $alignment{$sampleId} = substr($sequences->[$i], $start, $end - $start + 1);
        $alnReactivity{$sampleId} = [ @{$reactivity->[$i]}[$start .. $end] ];

    }

    return(\%alignment, \%alnReactivity);

}

##
#  Builds a ViennaRNA hard constraint string, in which bases to be forced
#  single-stranded are marked with an "x"
##

sub vrconstraint {

    my ($hardbases, $constraint) = @_;

    my (@pairs, %partner);
    @pairs = listpairs($constraint);
    %partner = map { $_->[0] => $_->[1],
                     $_->[1] => $_->[0] } @pairs;

    # A hard-constrained base inside a user-forced helix would leave its partner
    # orphaned, and the resulting unbalanced string throws in _nestedConstraint
    for (@{$hardbases}) {

        substr($constraint, $_, 1) = "x";
        substr($constraint, $partner{$_}, 1) = "x" if (exists $partner{$_});

    }
 
     return($constraint);
 
 }

##
#  Writes an RNAstructure hard constraint file
##

sub applyhardconstraint {

    my ($file, $hardbases, $constraint) = @_;

    $hardbases = [ sort {$a <=> $b} @{$hardbases} ];

    open(my $wh, ">", $file . ".constraint") or _logError("Unable to write temporary constraint file ($!)") and return();
    select((select($wh), $|=1)[0]);

    print $wh "DS:\n-1\nSS:\n";

    print $wh ++$_ . "\n" for (@{$hardbases});

    print $wh "-1\nMod:\n-1\nPairs:\n";

    if ($constraint) { print $wh ++$_->[0] . " " . ++$_->[1] . "\n" for (listpairs($constraint)); } # RNAstructure numbering is 1-based

    print $wh "-1 -1\nFMS:\n-1\nForbids:\n-1 -1";

    close($wh);

    return(1);

}

sub hardbases { # Returns 0-based coords of bases to be hard-constrained

    my (@positions);

    if ($hc &&
        !$ignore) {

        foreach my $i (0 .. $#_) { push(@positions, $i) if (isnumeric($_[$i]) &&
                                                            $_[$i] >= $cutoff); }

    }

    return(\@positions);

}

sub seq2fasta {

    my ($file, $id, $sequence) = @_;

    eval {

        my ($io, $fasta);
        $io = Data::IO::Sequence->new( file      => $file . ".fa",
                                       mode      => "w",
                                       overwrite => 1,
                                       flush     => 1,
                                       format    => "fasta" );
        $fasta = Data::Sequence->new( id       => $id,
                                      sequence => $sequence );

        $io->write($fasta);

    };

    if ($@) { 

        _logError($@, $id);

        return(); 
        
    }

    return(1);

}

##
# Pseudo dG calculation
#
# Calculates pseudo free energy contribution given structure and reactivity
# Essentially, the pseudodG calculated for each base, is counted twice if the base
# resides within a stack of bases, once if the base is at the terminal position of
# an helix, and is ignored if it's a lonely pair.
# When a bulge is encountered, the 2 bases immediately surrounding the bulge are
# counted once, as well as for the 2 bases to which they are paired.
##

sub pseudodG {

    my ($structure, $reactivity) = @_;

    my $dG = 0;

    foreach my $helix (@{(listhelices($structure))[0]}) {

        next if (@{$helix->{h5bases}} == 1);

        my (@h5bases, @h3bases, @diff, @external, @internal, %allbases);
        @h5bases = @{$helix->{h5bases}};
        @h3bases = @{$helix->{h3bases}};
        @external = (shift(@h5bases), shift(@h3bases), pop(@h5bases), pop(@h3bases));
        %allbases = map { $_ => 1 } (@h5bases, @h3bases);

        for ([\@h5bases, \@h3bases], [\@h3bases, \@h5bases]) {

            my (@h1, @h2, @diff);
            @h1 = @{$_->[0]};
            @h2 = @{$_->[1]};

            $diff[$_] = abs($h1[$_ + 1] - $h1[$_]) for(0 .. $#h1 - 1);

            for (my $i = 0; $i < $#h1; $i++) {

                if ($diff[$i] > 1) {

                    splice(@diff, $i, 1);
                    @external = (@external, splice(@h1, $i, 2), splice(@h2, $i, 2));
                    $i--;

                }

            }

        }

        @external = uniq(@external);
        delete($allbases{$_}) for (@external);
        @internal = map { $_ } keys %allbases;

        $dG += isnegative($reactivity->[$_]) || isnan($reactivity->[$_]) ? 0 : ($slope * log($reactivity->[$_] + 1) + $intercept) for (@external);
        $dG += isnegative($reactivity->[$_]) || isnan($reactivity->[$_]) ? 0 : 2 * ($slope * log($reactivity->[$_] + 1) + $intercept) for (@internal);

    }

    return($dG);

}

sub _logError {

    my ($error, $transcript) = @_;
    my $caller = (caller(1))[3];
    $caller =~ s/^main:://;
    $error =~ s/^/    /gm;

    print STDERR "[!] Exception in $caller()" . ($transcript ? ", transcript/file $transcript:" : ":") . "\n$error\n"; 

}

sub help {

    print "\n  [!] Error: Invalid option. Please check the help\n" if ($_[0]);

    die <<HELP;

 RF Fold (v$Core::Utils::VERSION)
 RNA Framework [http://www.rnaframework.com]

 Author:  Danny Incarnato (dincarnato[at]rnaframework.com)
 Summary: Produces RNA secondary structures using structural probing data to guide folding

 Usage:   rf-fold [Options] XML_folder/                                      # Whole transcriptome
          rf-fold [Options] transcript.xml                                   # Single transcript
          rf-fold [Options] XML_folder_1/ XML_folder_2/ .. XML_folder_n/     # Whole transcriptome, with replicates

 Options                                     Description
 -o   or --output-dir            <string>    Output directory (Default: rf_fold/)
 -ow  or --overwrite                         Overwrites output directory (if the specified path already exists)
 -ct  or --connectivity-table                Writes predicted structures in CT format (Default: Dot-bracket notation)
 -m   or --folding-method        <int>       Specifies the folding method (1-2, Default: 1): [1] ViennaRNA
                                                                                             [2] RNAstructure

 -p   or --processors            <int>       Number of processors to use for the analysis (Default: 1)
 -oc  or --only-common           <int>       In case of multiple experiments, only transcripts covered across at least this number of experiments
                                             will be folded
 -g   or --img                               Enables the generation of graphical reports (requires R)
 -gr  or --gradient                          Reactivity coloring is performed using a continuos gradient color scale rather than discrete
                                             reactivity binds (requires -g)
 -R   or --R-path                <string>    Path to R executable (Default: assumes R is in PATH)
 -t   or --temperature           <float>     Temperature in Celsius degrees (Default: 37.0)
 -sl  or --slope                 <float>     Sets slope used with structural probing data restraints (Default: 1.8 [kcal/mol])
 -in  or --intercept             <float>     Sets intercept used with structural probing data restraints (Default: -0.6 [kcal/mol])
 -md  or --maximum-distance      <int>       Sets the maximum pairing distance in nucleotides between transcript's residues (Default: 0 [No limit])
 -nlp or --no-lonelypairs                    Disallows lonely base-pairs (1 bp helices) inside predicted structure
 -i   or --ignore-reactivity                 Ignores XML reactivity data when performing folding (MFE unconstrained prediction)
 -is  or --ignore-sequence                   In case of multiple experiments, nucleotide differences (e.g. SNVs) between XML files are ignored
 -hc  or --hard-constraint                   Besides performing soft-constraint folding, allows specifying a reactivity cutoff (specified by -f)
                                             for hard-constraining a base to be single-stranded
 -c   or --constraints           <string>    Path to a directory containing constraint files (in dot-bracket notation), that will be used to
                                             enforce specific base-pairs in the structure models
                                             Note: pseudoknots are currently not supported
 -f   or --cutoff                <float>     Reactivity cutoff for constraining a position as unpaired (>0, Default: 0.7)
 -w   or --windowed                          Enables windowed folding
 -pt  or --partition             <string>    Path to partition executable (Default: assumes partition is in PATH)
                                             Note: by default, partition-smp will be used (if available).
 -pp  or --probabilityplot       <string>    Path to ProbabilityPlot executable (Default: assumes ProbabilityPlot is in PATH)
 -fw  or --fold-window           <int>       Window size (in nt) for performing MFE folding (>=50, Default: 600)
 -fo  or --fold-offset           <int>       Offset (in nt) for MFE folding window sliding (Default: 200)
 -pw  or --partition-window      <int>       Window size (in nt) for computing partition function (>=50, Default: 600)
 -po  or --partition-offset      <int>       Offset (in nt) for partition function window sliding (Default: 200)
 -wt  or --window-trim           <int>       Number of bases to trim from both ends of the partition windows to avoid end biases (Default: 100)
 -dp  or --dotplot                           Enables generation of dot-plots of base-pairing probabilities
 -dpo or --dotplotOnly                       Only performs the partition function step and generates the dot-plots,
                                             skipping MFE folding (implies -dp)
 -sh  or --shannon-entropy                   Enables generation of a WIGGLE track file with per-base Shannon entropies

 Pseudoknot detection
 -pk  or --pseudoknots                       Enables detection of pseudoknots (computationally intensive)
 -km  or --pseudoknot-method     <int>       Algorithm for pseudoknots prediction (1-2, Default: 1): [1] RNA Framework
                                                                                                     [2] ShapeKnots
                                             Note: the chosen folding method (specified by -m) affects the algorithm used by RNA Framework
                                                   (pseudoknot detection method #1) to define the initial MFE structure
 -kp1 or --pseudoknot-penalty1   <float>     Pseudoknot penalty P1 (Default: 0.35)
 -kp2 or --pseudoknot-penalty2   <float>     Pseudoknot penalty P2 (Default: 0.65)
 -kt  or --pseudoknot-tollerance <float>     Maximum tollerated deviation of suboptimal structures energy from MFE (>0-1, Default: 0.5 [50%])
 -kh  or --pseudoknot-helices    <int>       Number of candidate pseudoknotted helices to evaluate (>0, Default: 100)
 -kw  or --pseudoknot-window     <int>       Window size (in nt) for performing pseudoknots detection (Default: 600)
 -ko  or --pseudoknot-offset     <int>       Offset (in nt) for pseudoknots detection window sliding (Default: 200)
 -kc  or --pseudoknot-cutoff     <float>     Reactivity cutoff for retaining a pseudoknotted helix (0-1, Default: 0.5)
 
 |
 +- RNA Framework pseudoknots detection algorithm options
    -vrs or --vienna-rnasubopt      <string>    Path to ViennaRNA RNAsubopt executable (Default: assumes RNAsubopt is in PATH)
    -ke  or --subopt-delta-energy   <float>     Computes suboptimal structures with energy in a certain range of the optimum (>=0, Default: 1 [kcal/mol])
    -ks  or --pseudoknot-suboptimal <int>       Number of suboptimal structures to evaluate for pseudoknots prediction (>0, Default: 1000)
    -nz  or --no-zuker                          Disables the inclusion of Zuker suboptimal structures (reduces the sampled folding space)
    -zs  or --zuker-suboptimal                  Number of Zuker suboptimal structures to include (>0, Default: 1000)

 |
 +- ShapeKnots pseudoknots detection algorithm options
    -sk or --shapeknots             <string>    Path to ShapeKnots executable (Default: assumes ShapeKnots is in PATH)
                                                Note: by default, ShapeKnots-smp will be used (if available).

 Folding method #1 options (ViennaRNA)
 -vrf or --vienna-rnafold     <string>    Path to ViennaRNA RNAfold executable (Default: assumes RNAfold is in PATH)
 -ngu or --no-closing-gu                  Disallows G:U wobbles at the end of helices
 -ua  or --use-alifold                    Uses RNAalifold to restrain the folding with the reactivity profiles of all the
                                          samples at once, yielding a consensus structure, instead of folding each sample
                                          individually with RNAfold and aggregating the resulting base-pairs
                                          Note: this affects both the partition function and the MFE folding steps, but
                                                not pseudoknots detection. A single slope/intercept pair is allowed.

 Folding method #2 options (RNAstructure)
 -rs or --rnastructure        <string>    Path to RNAstructure Fold executable (Default: assumes Fold is in PATH)
                                          Note: by default, Fold-smp will be used (if available).
 -d  or --data-path           <string>    Path to RNAstructure data tables (Default: assumes DATAPATH environment variable is already set)

HELP

}