diff --git a/ChangeLog b/ChangeLog index 2458a32..1e330af 100644 --- a/ChangeLog +++ b/ChangeLog @@ -14,6 +14,8 @@ ## Interface and handling +- IntaRNAeval / --rri evaluates predefined RNA-RNA interactions (issue #184) + - compressed binary .agz accessibility caches for repeated screens (issue #245) - too short sequences are skipped with warning if they are shorter than the @@ -61,6 +63,24 @@ ################################################################################ ################################################################################ +261001 Alexander Mitrofanov + * IntaRNA/PredictorEvalOnly, src/IntaRNA/Makefile.am : + + parse colon-separated hybridDB structures with sequence/index/pair validation + + evaluate explicit base pairs with the selected energy and accessibility model + + report distinct structures, energies, restricted partition sums and trackers + * bin/CommandLineParsing, bin/IntaRNA : + + --rri selects evaluation; IntaRNAeval personality requires interaction input + + require one query/target and report ignored prediction and output constraints + * bypass prediction result filtering to preserve distinct structures with equal + energies, boundaries and base-pair counts; retain energy/accessibility settings + * tests/PredictorEvalOnly_test.cpp, tests/runIntaRNAeval.sh, tests/Makefile.am : + + check malformed structures, shifted indices, ignored filters, duplicate/tied + structures, tracker coordinates, partition reset and energy contributions + + reproduce ViennaRNA predictions through API and CLI evaluation round trips + * README.md and CLI help : + + document input format, retained settings and supplied-structure ensemble scope + * resolves https://github.com/BackofenLab/IntaRNA/issues/184 + 260930 Alexander Mitrofanov * IntaRNA/Matrix.h : + const/mutable upper-band row views that exclude structural zeros and padding diff --git a/README.md b/README.md index 4db61ce..746ed49 100644 --- a/README.md +++ b/README.md @@ -100,6 +100,7 @@ The following topics are covered by this documentation: - [IntaRNAsTar - optimized for sRNA-target prediction](#IntaRNAsTar) - [IntaRNAseed - identifys and reports seed interactions only](#IntaRNAseed) - [IntaRNAens - ensemble-based prediction and partition function computation](#IntaRNAens) + - [IntaRNAeval - evaluate predefined interactions](#IntaRNAeval) - [How to constrain predicted interactions](#constraintSetup) - [Interaction restrictions](#interConstr) - [Seed constraints](#seed) @@ -999,6 +1000,62 @@ IntaRNA --mode=S ... ``` +[![up](doc/figures/icon-up.28.png) back to overview](#overview) + + +### IntaRNAeval + +**IntaRNAeval** evaluates predefined RNA-RNA interactions with the selected +energy and accessibility model. It requires `--rri`; supplying `--rri` to +`IntaRNA` also enables evaluation automatically. + +Use the CSV `hybridDB` format, `startTdotbarT&startQdotbarQ`. Both strands are +written in their original 5'-3' direction: `|` marks a paired nucleotide and +`.` an unpaired nucleotide. Pairing is antiparallel, so the first target bar +pairs with the last query bar. Starts follow `--tIdxPos0` and `--qIdxPos0` +(default: 1), including the usual skipped zero when indexing starts negative. +Full-length `hybridDBfull` encodings with flanking dots are also accepted. +Separate multiple interactions with `:` and quote the argument in the shell. + +```sh +IntaRNAeval -t GGGG -q CCCC --rri='1||||&1||||:2||&2||' +# Equivalent calls: +IntaRNA --personality=IntaRNAeval -t GGGG -q CCCC --rri='1||||&1||||:2||&2||' +IntaRNA -t GGGG -q CCCC --rri='1||||&1||||:2||&2||' +``` + +Exactly one selected query and target are required. Each encoding must contain +at least one base pair, equal numbers of bars, in-range positions, and only +complementary AU, GC or GU pairs. Intramolecular pairs and crossing interactions +are not represented by this format. Every distinct supplied structure is +reported in input order, including single pairs and positive energies; repeated +structures count once. `--outCsvSort` can change the reporting order. + +Evaluation ignores prediction modes/models, seed and helix restrictions, +interaction length/loop limits, query/target search regions, prediction windows, +and output filters (including `--outNumber`, overlap, energy/accessibility +cutoffs, `--outNoLP` and `--outNoGUend`). Explicitly supplied ignored options +are listed in the INFO log. Sequence selection, output formats and columns, +energy parameters, temperature, `--energyAdd`, dangling ends, and accessibility +settings (including SHAPE and accessibility files) still apply. Accessibility +windows and input matrices must cover the supplied sites; structures with no +finite energy under the selected model cause an error. + +The reported energy includes initiation, all intervening loops/stackings, +accessibility penalties, dangling ends, terminal-pair penalties and the energy +shift. CSV energy contributions and ordinary text output use the same reporting +interfaces as predictions. No seed is assigned, so seed-related columns report +unavailable values. `Zall`, `Eall`, `P_E` and probability/minimum-energy trackers +refer **only to the distinct supplied structures**, not all possible +interactions of the sequences. Intramolecular ensemble quantities retain their +usual meaning. + +For library users, `IntaRNA/PredictorEvalOnly.h` provides `parseInteractions()` +and a `Predictor` subclass accepting a list of `Interaction` objects. The energy +model and its sequences must outlive the predictor. Library callers control +loop/accessibility limits through their energy model and storage/filtering +through their chosen output handler. + [![up](doc/figures/icon-up.28.png) back to overview](#overview) ### IntaRNAens diff --git a/src/IntaRNA/Makefile.am b/src/IntaRNA/Makefile.am index 22291e6..9828816 100644 --- a/src/IntaRNA/Makefile.am +++ b/src/IntaRNA/Makefile.am @@ -65,6 +65,7 @@ libIntaRNA_a_HEADERS = \ PredictionTrackerSpotProbAll.h \ PredictionTrackerProfileSpotProb.h \ Predictor.h \ + PredictorEvalOnly.h \ PredictorMfe.h \ PredictorMfeSeedOnly.h \ PredictorMfe2d.h \ @@ -126,6 +127,7 @@ libIntaRNA_a_SOURCES = \ PredictionTrackerSpotProb.cpp \ PredictionTrackerSpotProbAll.cpp \ PredictionTrackerProfileSpotProb.cpp \ + PredictorEvalOnly.cpp \ PredictorMfe.cpp \ PredictorMfeSeedOnly.cpp \ PredictorMfe2d.cpp \ diff --git a/src/IntaRNA/PredictorEvalOnly.cpp b/src/IntaRNA/PredictorEvalOnly.cpp new file mode 100644 index 0000000..8c5acde --- /dev/null +++ b/src/IntaRNA/PredictorEvalOnly.cpp @@ -0,0 +1,162 @@ +#include "IntaRNA/PredictorEvalOnly.h" + +#include +#include +#include +#include + +namespace IntaRNA { + +namespace { + +// Parse one strand without performing arithmetic on unvalidated external indices. +std::vector parseStrand( const std::string_view strand, const RnaSequence & sequence ) +{ + const size_t structureStart = strand.find_first_of(".|"); + if (structureStart == std::string_view::npos || structureStart == 0) { + throw std::invalid_argument("--rri: expected start1dotbar1&start2dotbar2, e.g. 1|||&1|||"); + } + long start = 0; + const auto parsed = std::from_chars(strand.data(), strand.data()+structureStart, start); + if (parsed.ec != std::errc() || parsed.ptr != strand.data()+structureStart + || start < sequence.getInOutIndex(0) + || start > sequence.getInOutIndex(sequence.size()-1) + || (start == 0 && sequence.getInOutIndex(0) < 0)) { + throw std::invalid_argument("--rri: invalid or out-of-range start index"); + } + const size_t offset = sequence.getIndex(start); + const auto structure = strand.substr(structureStart); + if (structure.size() > sequence.size()-offset || structure.find_first_not_of(".|") != std::string_view::npos) { + throw std::invalid_argument("--rri: dot-bar strand exceeds sequence length or contains invalid symbols"); + } + std::vector paired; + for (size_t p=0; p +PredictorEvalOnly::parseInteractions( const std::string & encoding, + const RnaSequence & target, const RnaSequence & query ) +{ + std::vector result; + size_t start = 0; + do { + const size_t end = encoding.find(':', start); + const auto entry = std::string_view(encoding).substr(start, + end == std::string::npos ? end : end-start); + const size_t separator = entry.find('&'); + if (separator == std::string_view::npos || entry.find('&', separator+1) != std::string_view::npos) { + throw std::invalid_argument("--rri: each interaction requires exactly one '&' separator"); + } + const auto paired1 = parseStrand(entry.substr(0, separator), target); + const auto paired2 = parseStrand(entry.substr(separator+1), query); + if (paired1.empty() || paired1.size() != paired2.size()) { + throw std::invalid_argument("--rri: strands must contain the same nonzero number of pairing bars"); + } + Interaction interaction(target, query); + for (size_t p=0; p & input ) + : Predictor(energy, output, predTracker) +{ + const auto & target = energy.getAccessibility1().getSequence(); + const auto & query = energy.getAccessibility2().getAccessibilityOrigin().getSequence(); + if (input.empty()) throw std::invalid_argument("PredictorEvalOnly: no interactions provided"); + std::set seen; + for (const auto & interaction : input) { + if (!interaction.s1 || !interaction.s2 || !interaction.isValid() + || interaction.s1->asString() != target.asString() + || interaction.s2->asString() != query.asString()) { + throw std::invalid_argument("PredictorEvalOnly: invalid interaction or incompatible sequences"); + } + for (const auto & bp : interaction.basePairs) { + if (bp.first >= target.size() || bp.second >= query.size() + || !RnaSequence::areComplementary(target, query, bp.first, bp.second)) { + throw std::invalid_argument("PredictorEvalOnly: out-of-range or non-complementary base pair"); + } + } + if (seen.insert(interaction.basePairs).second) { + interactions.emplace_back(target, query); + interactions.back().basePairs = interaction.basePairs; + } + } +} + +void +PredictorEvalOnly::predict( const IndexRange &, const IndexRange & ) +{ + initOptima(); + // Finish evaluation before reporting anything, so an unevaluable structure + // cannot leave a partially reported list or partition function. + for (auto & interaction : interactions) { + E_type hybridE = energy.getE_init(); + for (size_t p=1; pupdateOptimumCalled(i1, j1, i2, j2, totalE); +} + +void +PredictorEvalOnly::reportOptima() +{ + output.incrementZ(Zall); + for (const auto & interaction : interactions) output.add(interaction); +} + +} // namespace IntaRNA diff --git a/src/IntaRNA/PredictorEvalOnly.h b/src/IntaRNA/PredictorEvalOnly.h new file mode 100644 index 0000000..cf9c1df --- /dev/null +++ b/src/IntaRNA/PredictorEvalOnly.h @@ -0,0 +1,72 @@ +#ifndef INTARNA_PREDICTOREVALONLY_H_ +#define INTARNA_PREDICTOREVALONLY_H_ + +#include "IntaRNA/Predictor.h" + +#include +#include + +namespace IntaRNA { + +/** + * Evaluates predefined, nested intermolecular base pairs without searching. + * Prediction ranges and output filters are ignored. Identical structures are + * counted once; Zall and trackers describe only the supplied structures. + */ +class PredictorEvalOnly : public Predictor { +public: + + /** + * Copies and validates the structures, discarding input energies and seeds. + * @param energy the energy model; it and its sequences must outlive this object + * @param output the reporting destination, which must outlive this object + * @param predTracker optional tracker owned by this predictor + * @param interactions nonempty list with ascending target/descending query + * base-pair indices in the original, zero-based sequence coordinates + * @throws std::invalid_argument for empty, incompatible or invalid structures + */ + PredictorEvalOnly( const InteractionEnergy & energy, OutputHandler & output, + PredictionTracker * predTracker, const std::vector & interactions ); + + /** + * Evaluates all structures using initiation, loop, accessibility, dangling-end, + * terminal-pair and energy-shift contributions from the selected model. + * @param r1 ignored; every supplied interaction is evaluated + * @param r2 ignored; every supplied interaction is evaluated + * @throws std::runtime_error if the model cannot assign a finite energy + */ + void predict( const IndexRange & r1 = IndexRange(0,RnaSequence::lastPos), + const IndexRange & r2 = IndexRange(0,RnaSequence::lastPos) ) override; + + /** + * Parses a colon-separated list of hybridDB (start1dotbar1&start2dotbar2) + * encodings, including full-length encodings with flanking dots. Pairing bars + * are matched antiparallel. Starts use each sequence's input/output indexing. + * @param encoding one or more nonempty encodings, each containing a base pair + * @param target the first sequence; must outlive the returned interactions + * @param query the second sequence in its original 5'-3' orientation + * @return validated interactions referencing the provided sequences + * @throws std::invalid_argument for malformed, out-of-bounds, unbalanced or + * non-complementary structures + */ + static std::vector parseInteractions( const std::string & encoding, + const RnaSequence & target, const RnaSequence & query ); + +protected: + //! reset the partition function before each evaluation + void initOptima() override; + //! record an evaluated energy for the partition function and tracker + void updateOptima( const size_t i1, const size_t j1, + const size_t i2, const size_t j2, const E_type energy, + const bool isHybridE, const bool incrementZall ) override; + //! forward the partition function and all evaluated structures to output + void reportOptima() override; + +private: + //! unique structures referencing the energy model's original sequences + std::vector interactions; +}; + +} // namespace IntaRNA + +#endif diff --git a/src/bin/CommandLineParsing.cpp b/src/bin/CommandLineParsing.cpp index e8b0a64..8da0809 100644 --- a/src/bin/CommandLineParsing.cpp +++ b/src/bin/CommandLineParsing.cpp @@ -8,6 +8,7 @@ #include #include #include +#include #if INTARNA_MULITHREADING #include @@ -41,6 +42,7 @@ extern "C" { #include "IntaRNA/PredictorMfe2dHeuristic.h" #include "IntaRNA/PredictorMfe2dHelixBlockHeuristic.h" #include "IntaRNA/PredictorMfe2d.h" +#include "IntaRNA/PredictorEvalOnly.h" #include "IntaRNA/PredictorMfeSeedOnly.h" #include "IntaRNA/PredictorMfe2dHeuristicSeed.h" @@ -90,6 +92,8 @@ const std::string CommandLineParsing::outCsvLstSep = ":"; CommandLineParsing::CommandLineParsing( const Personality personality ) : personality( personality ), + rri(), + rriInteractions(), personalityParamValue(""), stdinUsed(false), opts_query("Query"), @@ -226,6 +230,7 @@ CommandLineParsing::CommandLineParsing( const Personality personality ) switch (personality) { case IntaRNA : case IntaRNA3 : + case IntaRNAeval : // no changes break; case IntaRNA1 : @@ -786,6 +791,11 @@ CommandLineParsing::CommandLineParsing( const Personality personality ) //// INTERACTION/ENERGY OPTIONS //////////////////////// opts_inter.add_options() + ("rri", value(&rri), + "evaluate predefined interactions in hybridDB format: startTdotbarT&startQdotbarQ" + " (e.g. '1|||&1|||'); separate interactions with ':'. Requires one query and target." + " Ignores prediction/seed/helix/region/window and output-filter constraints." + " Energy and accessibility settings remain in effect; ensemble output covers only the supplied structures.") ((mode.name+",m").c_str() , value(&(mode.val)) ->default_value(mode.def) @@ -1187,6 +1197,12 @@ parse(int argc, char** argv) // if parsing was successful, continue with additional checks if (parsingCode == ReturnCode::KEEP_GOING) { try { + if (vm.count("rri") || personality == IntaRNAeval) { + if (!vm.count("rri") || vm.at("rri").as().empty()) { + throw error("IntaRNAeval requires a nonempty --rri interaction input"); + } + prepareEvaluation(vm); + } // run all notifier checks notify(vm); } catch (required_option& e) { @@ -1221,6 +1237,12 @@ parse(int argc, char** argv) // parse the sequences parseSequences("query",qId,queryArg,query,qSet,qIdxPos0.val); parseSequences("target",tId,targetArg,target,tSet,tIdxPos0.val); + if (!rri.empty()) { + if (query.size() != 1 || target.size() != 1) { + throw error("--rri requires exactly one query and one target sequence"); + } + rriInteractions = PredictorEvalOnly::parseInteractions(rri, target.front(), query.front()); + } // check if same number if pairwise mode if (outPairwise && query.size() != target.size()) { @@ -1641,6 +1663,41 @@ parse(int argc, char** argv) //////////////////////////////////////////////////////////////////////////// +void +CommandLineParsing::prepareEvaluation( boost::program_options::variables_map & vm ) +{ + // Erase before notify(), so even out-of-sequence regions and seed encodings + // cannot constrain or invalidate evaluation. Energy/accessibility options stay. + const std::set ignored = { + "model", "mode", "noSeed", "intLenMax", "qIntLenMax", "tIntLenMax", + "intLoopMax", "qIntLoopMax", "tIntLoopMax", "qRegion", "tRegion", + "qRegionLenMax", "tRegionLenMax", "windowWidth", "windowOverlap", + "outNumber", "outOverlap", "outMaxE", "outDeltaE", "outMinPu", + "outNoLP", "outNoGUend", "outBestSeedOnly", "outPerRegion", "outPairwise" + }; + for (auto it=vm.begin(); it!=vm.end();) { + if (ignored.count(it->first) || it->first.starts_with("seed") || it->first.starts_with("helix")) { + if (!it->second.defaulted()) LOG(INFO) <<"--rri evaluation: ignoring --"<first; + it = vm.erase(it); + } else { + ++it; + } + } + // Override prediction defaults, including those of another personality. + model.val = 'S'; + mode.val = 'H'; + noSeedRequired = true; + intLenMax.val = intLenMax.def = 0; + qIntLenMax.val = tIntLenMax.val = 0; + windowWidth.val = 0; + qRegionLenMax.val = tRegionLenMax.val = 0; + outNoLP = outNoGUend = outPerRegion = outPairwise = false; + LOG(INFO) <<"Evaluating predefined interactions; prediction constraints are ignored. " + <<"Energy/accessibility settings are retained; ensemble statistics cover only the supplied structures."; +} + +//////////////////////////////////////////////////////////////////////////// + void CommandLineParsing:: validate_charArgument(const CommandLineParsing::CharParameter& param, const char & value) @@ -2076,13 +2133,15 @@ getEnergyHandler( const Accessibility& accTarget, const ReverseAccessibility& ac // check whether to compute ES values (for multi-site predictions) const bool initES = std::string("M").find(model.val) != std::string::npos; + const size_t loopMax1 = rri.empty() ? tIntLoopMax.val : accTarget.getSequence().size(); + const size_t loopMax2 = rri.empty() ? qIntLoopMax.val : accQuery.getSequence().size(); switch( energy.val ) { case 'B' : return new InteractionEnergyBasePair( accTarget, accQuery - , tIntLoopMax.val, qIntLoopMax.val + , loopMax1, loopMax2 , initES, Z_type(1.0), Ekcal_2_E(-1), 3 , Ekcal_2_E(energyAdd.val), !energyNoDangles, !outNoGUend ); - case 'V' : return new InteractionEnergyVrna( accTarget, accQuery, vrnaHandler, tIntLoopMax.val, qIntLoopMax.val, initES, Ekcal_2_E(energyAdd.val), !energyNoDangles, !outNoGUend ); + case 'V' : return new InteractionEnergyVrna( accTarget, accQuery, vrnaHandler, loopMax1, loopMax2, initES, Ekcal_2_E(energyAdd.val), !energyNoDangles, !outNoGUend ); default : INTARNA_NOT_IMPLEMENTED("CommandLineParsing::getEnergyHandler : energy = '"+toString(energy.val)+"' is not supported"); } @@ -2097,6 +2156,10 @@ CommandLineParsing:: getOutputConstraint( const InteractionEnergy & energy ) const { checkIfParsed(); + if (!rri.empty()) { + return OutputConstraint(rriInteractions.size(), OutputConstraint::OVERLAP_BOTH, + E_INF, E_INF, false, false, false, outNeedsZall, true); + } OutputConstraint::ReportOverlap overlap = OutputConstraint::ReportOverlap::OVERLAP_BOTH; switch(outOverlap.val) { case 'N' : overlap = OutputConstraint::ReportOverlap::OVERLAP_NONE; break; @@ -2425,6 +2488,10 @@ getPredictor( const InteractionEnergy & energy, OutputHandler & output ) const predTracker = NULL; } + if (!rri.empty()) { + return new PredictorEvalOnly(energy, output, predTracker, rriInteractions); + } + if (noSeedRequired) { // predictors without seed constraint switch( model.val ) { @@ -2829,6 +2896,9 @@ getPersonality( int argc, char ** argv ) } // parse personality + if (value == "IntaRNAeval") { + return Personality::IntaRNAeval; + } if (boost::regex_match(value,boost::regex("IntaRNAexact"), boost::match_perl)) { return Personality::IntaRNAexact; } diff --git a/src/bin/CommandLineParsing.h b/src/bin/CommandLineParsing.h index 8594c50..3731406 100644 --- a/src/bin/CommandLineParsing.h +++ b/src/bin/CommandLineParsing.h @@ -4,6 +4,7 @@ #include "IntaRNA/general.h" #include "IntaRNA/RnaSequence.h" +#include "IntaRNA/Interaction.h" #include #include @@ -14,6 +15,7 @@ #include #include #include +#include #include "IntaRNA/Accessibility.h" #include "IntaRNA/InteractionEnergy.h" @@ -57,6 +59,7 @@ class CommandLineParsing { IntaRNA2, // IntaRNA v2 like setup IntaRNA3, // default IntaRNA v3 setup IntaRNAens, // ensemble-based prediction + IntaRNAeval, // evaluate predefined interactions IntaRNAsTar, // sRNA-target prediction (optimized parameter) IntaRNAseed, // seed-only predictions IntaRNAhelix, // helix-block-based predictions @@ -80,6 +83,7 @@ class CommandLineParsing { case IntaRNA2 : return "IntaRNA2"; case IntaRNA3 : return "IntaRNA3"; case IntaRNAens : return "IntaRNAens"; + case IntaRNAeval : return "IntaRNAeval"; case IntaRNAsTar : return "IntaRNAsTar"; case IntaRNAseed : return "IntaRNAseed"; case IntaRNAhelix : return "IntaRNAhelix"; @@ -141,6 +145,12 @@ class CommandLineParsing { CommandLineParsing( const Personality personality ); virtual ~CommandLineParsing(); + /** + * Whether predefined interactions are evaluated instead of predicted. + * @return true if a nonempty --rri input was parsed + */ + bool isEvaluation() const; + /** * Parses the commandline arguments as passed to the 'main' method * @param argc the number of arguments @@ -515,6 +525,11 @@ class CommandLineParsing { //! what is the requested personality for which we parse the parameters Personality personality; + //! colon-separated predefined interactions; nonempty enables evaluation + std::string rri; + //! validated structures referencing the parsed target and query sequences + std::vector rriInteractions; + //! might hold the personality string after parsing if given via parameter std::string personalityParamValue; @@ -1131,6 +1146,13 @@ class CommandLineParsing { */ void initOutputHandler(); + /** + * Removes prediction-only arguments before their notifiers run and reports + * explicitly supplied arguments that evaluation ignores. + * @param vm parsed options from the command line and parameter file + */ + void prepareEvaluation( boost::program_options::variables_map & vm ); + /** * Writes the accessibility to file or stream if requested by the user * @param acc the accessibility data assigned @@ -1160,6 +1182,13 @@ class CommandLineParsing { +inline +bool +CommandLineParsing::isEvaluation() const +{ + return !rri.empty(); +} + //////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////// diff --git a/src/bin/IntaRNA.cpp b/src/bin/IntaRNA.cpp index ab88dfe..b7bdbd4 100644 --- a/src/bin/IntaRNA.cpp +++ b/src/bin/IntaRNA.cpp @@ -285,7 +285,11 @@ int main(int argc, char **argv){ <<" ..."; } // get interaction prediction handler - std::unique_ptr predictor(parameters.getPredictor( *energy, bestInteractions )); + // Evaluation reports every distinct supplied structure directly. + // The prediction collector merges equal-energy structures with + // identical boundaries and pair counts, even if inner pairs differ. + std::unique_ptr predictor(parameters.getPredictor( *energy, + parameters.isEvaluation() ? *output : bestInteractions )); INTARNA_CHECK_NOT_NULL(predictor.get(),"predictor initialization failed"); // run prediction for this window combination @@ -448,4 +452,3 @@ int main(int argc, char **argv){ el::Loggers::flushAll(); return 0; } - diff --git a/tests/Makefile.am b/tests/Makefile.am index 4a4b95d..571deef 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -15,7 +15,7 @@ TEST_EXTENSIONS = $(EXEEXT) .sh SH_LOG_COMPILER = $(SHELL) # the script needed for tests -dist_check_SCRIPTS = runIntaRNA.sh runAccessibilityBinary.sh +dist_check_SCRIPTS = runIntaRNA.sh runAccessibilityBinary.sh runIntaRNAeval.sh # the program to build check_PROGRAMS = runApiTests @@ -50,6 +50,7 @@ runApiTests_SOURCES = \ PredictorMfeHeuristicCellState_test.cpp \ PredictorMfeEnsRegression_test.cpp \ PredictorTinyOracle_test.cpp \ + PredictorEvalOnly_test.cpp \ PredictorSeedOracle_test.cpp \ Matrix_test.cpp \ NussinovHandler_test.cpp \ diff --git a/tests/PredictorEvalOnly_test.cpp b/tests/PredictorEvalOnly_test.cpp new file mode 100644 index 0000000..96c92ae --- /dev/null +++ b/tests/PredictorEvalOnly_test.cpp @@ -0,0 +1,182 @@ +#include "catch.hpp" + +#undef NDEBUG + +#include "IntaRNA/AccessibilityDisabled.h" +#include "IntaRNA/InteractionEnergyBasePair.h" +#include "IntaRNA/InteractionEnergyVrna.h" +#include "IntaRNA/OutputHandlerInteractionList.h" +#include "IntaRNA/PredictorEvalOnly.h" +#include "IntaRNA/PredictorMfe2d.h" + +#include +#include +#include +#include + +using namespace IntaRNA; + +namespace { + +class EvalAccessibility : public Accessibility { +public: + EvalAccessibility(const RnaSequence & sequence, const E_type perBase); + E_type getED(const size_t from, const size_t to) const override; +private: + E_type perBase; +}; + +inline EvalAccessibility::EvalAccessibility(const RnaSequence & sequence, const E_type perBase) + : Accessibility(sequence, 0, NULL), perBase(perBase) +{} + +inline E_type EvalAccessibility::getED(const size_t from, const size_t to) const +{ + checkIndices(from, to); + return (to-from+1)*perBase; +} + +using EvalUpdate = std::tuple; + +class EvalTracker : public PredictionTracker { +public: + EvalTracker(std::vector & updates); + void updateOptimumCalled(const size_t i1, const size_t j1, + const size_t i2, const size_t j2, const E_type energy) override; +private: + std::vector & updates; +}; + +inline EvalTracker::EvalTracker(std::vector & updates) : updates(updates) {} + +inline void EvalTracker::updateOptimumCalled(const size_t i1, const size_t j1, + const size_t i2, const size_t j2, const E_type energy) +{ + updates.emplace_back(i1, j1, i2, j2, energy); +} + +} // namespace + +TEST_CASE("Evaluation parses hybridDB in original sequence coordinates", "[PredictorEvalOnly]") +{ + #include "testEasyLoggingSetup.icc" + RnaSequence target("t", "GGAGG", -2), query("q", "CCACC", 8); + auto input = PredictorEvalOnly::parseInteractions("-2||.||&8||.||:-1|&9|", target, query); + REQUIRE(input.size() == 2); + const Interaction::PairingVec expected = {{0,4}, {1,3}, {3,1}, {4,0}}; + REQUIRE(input[0].basePairs == expected); + REQUIRE(Interaction::dotBar(input[0]) == "-2||.||&8||.||"); + REQUIRE(input[1].basePairs.front() == Interaction::BasePair(1,1)); + REQUIRE_THROWS_AS(PredictorEvalOnly::parseInteractions("0|&8|", target, query), std::invalid_argument); + + RnaSequence g("g", "GGGG"), c("c", "CCCC"); + auto full = PredictorEvalOnly::parseInteractions("1.||.&1.||.", g, c); + REQUIRE(Interaction::dotBar(full[0], true) == "1.||.&1.||."); + for (const auto & bad : {"", "1|||", "1|&1|&1|", "1|&1||", "1...&1...", + "1x|&1||", "1(&1)", "0|&1|", "-1|&1|", "4||&1||", "1|&5|", + "9223372036854775808|&1|", "-9223372036854775809|&1|", "1|&1|:", ":1|&1|", "1|&1|::1|&1|"}) { + CAPTURE(bad); + REQUIRE_THROWS_AS(PredictorEvalOnly::parseInteractions(bad, g, c), std::invalid_argument); + } + REQUIRE_THROWS_AS(PredictorEvalOnly::parseInteractions("1|&1|", g, g), std::invalid_argument); + RnaSequence n("n", "NCCC"); + REQUIRE_THROWS_AS(PredictorEvalOnly::parseInteractions("1|&1|", g, n), std::invalid_argument); + RnaSequence zero("zero", "GGGG", 0); + REQUIRE(PredictorEvalOnly::parseInteractions("0|&1|", zero, c)[0].basePairs.front().first == 0); +} + +TEST_CASE("Evaluation reports energies, structures and restricted ensemble without filtering", "[PredictorEvalOnly]") +{ + #include "testEasyLoggingSetup.icc" + RnaSequence target("t", "GGGGGG"), query("q", "UUUUUUU"); + EvalAccessibility acc1(target, 20), acc2(query, 30); + ReverseAccessibility reversed(acc2); + InteractionEnergyBasePair energy(acc1, reversed, 20, 20, false, 1, -100, 3, 500); + // Deliberately exclude these positive-energy, lonely GU pairs with filters. + OutputConstraint constraints(0, OutputConstraint::OVERLAP_NONE, -1000, 0, true, true, true, true, false, 0); + OutputHandlerInteractionList output(constraints, 10); + std::vector updates; + auto input = PredictorEvalOnly::parseInteractions("2||.||&2|..|||:1|&1|:2||.||&2|..|||", target, query); + input[0].energy = -12345; + input[0].setSeedRange(input[0].basePairs.front(), input[0].basePairs.back(), -999); + PredictorEvalOnly predictor(energy, output, new EvalTracker(updates), input); + // Ranges do not trim or shift the supplied coordinates. + predictor.predict(IndexRange(0,0), IndexRange(0,0)); + REQUIRE(output.reported() == 2); + auto it = output.begin(); + REQUIRE((*it)->energy == 380); // -4 bp + 5*0.2 ED1 + 6*0.3 ED2 + 5 shift + REQUIRE((*it)->basePairs == input[0].basePairs); + REQUIRE((*it)->seed == NULL); + REQUIRE((*++it)->energy == 450); + REQUIRE(updates.size() == 2); + REQUIRE(updates[0] == EvalUpdate(1,5,0,5,380)); + const double partition = std::exp(-3.8) + std::exp(-4.5); + REQUIRE(static_cast(predictor.getZall()) == Approx(partition)); + REQUIRE(static_cast(output.getZ()) == Approx(partition)); + predictor.predict(); + REQUIRE(static_cast(predictor.getZall()) == Approx(partition)); + REQUIRE(static_cast(output.getZ()) == Approx(2*partition)); + REQUIRE(updates.size() == 4); +} + +TEST_CASE("Evaluation validates API inputs and reports unavailable energies before output", "[PredictorEvalOnly]") +{ + #include "testEasyLoggingSetup.icc" + RnaSequence target("t", "GGGG"), query("q", "CCCC"); + AccessibilityDisabled acc1(target, 0, NULL), acc2(query, 0, NULL); + ReverseAccessibility reversed(acc2); + InteractionEnergyBasePair energy(acc1, reversed, 0, 0); + OutputHandlerInteractionList output(OutputConstraint(), 10); + REQUIRE_THROWS_AS(PredictorEvalOnly(energy, output, NULL, {}), std::invalid_argument); + Interaction invalid(target, query); + REQUIRE_THROWS_AS(PredictorEvalOnly(energy, output, NULL, {invalid}), std::invalid_argument); + invalid.basePairs = {{0,0}, {1,1}}; + REQUIRE_THROWS_AS(PredictorEvalOnly(energy, output, NULL, {invalid}), std::invalid_argument); + invalid.basePairs = {{4,0}}; + REQUIRE_THROWS_AS(PredictorEvalOnly(energy, output, NULL, {invalid}), std::invalid_argument); + Interaction mismatch(query, target); + mismatch.basePairs = {{0,0}}; + REQUIRE_THROWS_AS(PredictorEvalOnly(energy, output, NULL, {mismatch}), std::invalid_argument); + auto input = PredictorEvalOnly::parseInteractions("1|&1|:1|..|&1|..|", target, query); + PredictorEvalOnly predictor(energy, output, NULL, input); + REQUIRE_THROWS_AS(predictor.predict(), std::runtime_error); + REQUIRE(output.empty()); + REQUIRE(output.getZ() == 0); + + AccessibilityDisabled narrow(target, 1, NULL); + InteractionEnergyBasePair limited(narrow, reversed, 10, 10); + PredictorEvalOnly inaccessible(limited, output, NULL, input); + REQUIRE_THROWS_AS(inaccessible.predict(), std::runtime_error); + REQUIRE(output.empty()); +} + +TEST_CASE("Evaluation reproduces ViennaRNA prediction energies and contributions", "[PredictorEvalOnly]") +{ + #include "testEasyLoggingSetup.icc" + RnaSequence target("t", "AGCGACGCA"), query("q", "UGCGUCGCU"); + EvalAccessibility acc1(target, 7), acc2(query, 13); + ReverseAccessibility reversed(acc2); + VrnaHandler vrna; + for (const bool dangles : {false, true}) { + InteractionEnergyVrna energy(acc1, reversed, vrna, 20, 20, false, 37, dangles); + OutputConstraint constraints(20, OutputConstraint::OVERLAP_BOTH, E_INF, E_INF); + OutputHandlerInteractionList predicted(constraints, 20), evaluated(constraints, 20); + PredictorMfe2d search(energy, predicted, NULL); + search.predict(); + REQUIRE_FALSE(predicted.empty()); + std::vector input; + for (const auto * interaction : predicted) input.push_back(*interaction); + PredictorEvalOnly evaluation(energy, evaluated, NULL, input); + evaluation.predict(); + REQUIRE(evaluated.reported() == predicted.reported()); + auto result = evaluated.begin(); + for (const auto * reference : predicted) { + REQUIRE((*result)->basePairs == reference->basePairs); + REQUIRE((*result)->energy == reference->energy); + const auto parts = energy.getE_contributions(**result); + REQUIRE(parts.init + parts.loops + parts.ED1 + parts.ED2 + parts.dangleLeft + + parts.dangleRight + parts.endLeft + parts.endRight + parts.energyAdd == reference->energy); + ++result; + } + } +} diff --git a/tests/runIntaRNAeval.sh b/tests/runIntaRNAeval.sh new file mode 100644 index 0000000..71a1eae --- /dev/null +++ b/tests/runIntaRNAeval.sh @@ -0,0 +1,111 @@ +#!/usr/bin/env bash +# Exercise evaluation dispatch, validation, energy round trips and personalities. +set -euo pipefail +bin="$INTARNABINPATH/src/bin/IntaRNA" +tmp=$(mktemp -d) +trap 'rm -rf "$tmp"' EXIT +common=(--target=GGGGGG --query=UUUUUUU --energy=B --acc=N --threads=1 --default-log-file="$tmp/info.log") +csv=(--outMode=C --outCsvCols=hybridDB,E) +rri='2||.||&2|..|||:1|&1|' +"$bin" "${common[@]}" "${csv[@]}" --rri="$rri" > "$tmp/reference" +printf 'hybridDB;E\n2||.||&2|..|||;-4\n1|&1|;-1\n' > "$tmp/expected" +cmp "$tmp/expected" "$tmp/reference" + +# Distinct inner pairs must survive even with the same energy and boundaries. +"$bin" --target=GGGG --query=CCCC --energy=B --acc=N "${csv[@]}" \ + --rri='1||.|&1||.|:1|.||&1|.||' --default-log-file="$tmp/info.log" > "$tmp/tied" +printf 'hybridDB;E\n1||.|&1||.|;-3\n1|.||&1|.||;-3\n' > "$tmp/expected" +cmp "$tmp/expected" "$tmp/tied" + +# Prediction filters, invalid regions/seeds, length limits and windows are ignored. +"$bin" "${common[@]}" "${csv[@]}" --rri="$rri" --model=B --mode=S --seedBP=20 --seedTQ=invalid \ + --seedQRange=999-1000 --seedNoGU --helixMinBP=4 --helixMaxBP=2 \ + --qRegion=999-1000 --tRegion=999-1000 --qRegionLenMax=1 --tRegionLenMax=1 \ + --intLenMax=1 --qIntLenMax=2 --tIntLenMax=3 --intLoopMax=0 --qIntLoopMax=1 \ + --windowWidth=5 --windowOverlap=100 --outNumber=0 --outOverlap=N --outDeltaE=0 \ + --outMaxE=-999 --outMinPu=1 --outNoLP --outNoGUend --outBestSeedOnly --outPerRegion \ + > "$tmp/ignored" 2> "$tmp/ignored.log" +cmp "$tmp/reference" "$tmp/ignored" +grep -q 'ignoring --outNumber' "$tmp/info.log" +grep -q 'ignoring --seedTQ' "$tmp/info.log" + +# Both personality dispatch paths, plus overriding another personality with --rri. +ln -s "$bin" "$tmp/IntaRNAeval" +"$tmp/IntaRNAeval" "${common[@]}" "${csv[@]}" --rri="$rri" > "$tmp/personality" +cmp "$tmp/reference" "$tmp/personality" +for personality in IntaRNAeval IntaRNAseed IntaRNAhelix IntaRNAsTar; do + "$bin" "${common[@]}" "${csv[@]}" --personality="$personality" --rri="$rri" > "$tmp/personality" + cmp "$tmp/reference" "$tmp/personality" +done +printf 'rri=%s\nseedTQ=invalid\noutNumber=0\n' "$rri" > "$tmp/parameters" +"$bin" "${common[@]}" "${csv[@]}" --parameterFile="$tmp/parameters" > "$tmp/configured" +cmp "$tmp/reference" "$tmp/configured" + +# Energy shifts and positive energies remain visible, even with default outMaxE. +"$bin" "${common[@]}" "${csv[@]}" --rri="$rri" --energyAdd=10 > "$tmp/positive" +printf 'hybridDB;E\n2||.||&2|..|||;6\n1|&1|;9\n' > "$tmp/expected" +cmp "$tmp/expected" "$tmp/positive" +"$bin" "${common[@]}" "${csv[@]}" --tIdxPos0=-2 --qIdxPos0=8 \ + --rri='-1||.||&9|..|||:-2|&8|' > "$tmp/shifted" +printf 'hybridDB;E\n-1||.||&9|..|||;-4\n-2|&8|;-1\n' > "$tmp/expected" +cmp "$tmp/expected" "$tmp/shifted" + +# Repeated structures do not inflate Zall or the number of results. +for structures in "$rri" "$rri:$rri"; do + "$bin" "${common[@]}" --outMode=C --rri="$structures" --outCsvCols=hybridDB,E,Zall,Eall,P_E \ + > "$tmp/ensemble-current" + if test -f "$tmp/ensemble"; then cmp "$tmp/ensemble" "$tmp/ensemble-current"; fi + cp "$tmp/ensemble-current" "$tmp/ensemble" +done +"$bin" "${common[@]}" --rri="$rri" --outMode=E > "$tmp/ensemble-mode" +test -s "$tmp/ensemble-mode" + +# Loops beyond the normal prediction limits can be evaluated. +printf -v dots '%40s' '' +dots=${dots// /.} +printf -v gs '%42s' '' +gs=${gs// /G} +cs=${gs//G/C} +"$bin" --target="$gs" --query="$cs" --acc=N --energy=B --intLoopMax=0 \ + --rri="1|${dots}|&1|${dots}|" --outMode=C --outCsvCols=E --default-log-file="$tmp/info.log" > "$tmp/long-loop" +printf 'E\n-2\n' > "$tmp/expected" +cmp "$tmp/expected" "$tmp/long-loop" + +# Re-evaluate ViennaRNA predictions with the same energy/accessibility settings. +thermo=(--target=AGCGACGCA --query=UGCGUCGCU --accW=0 --accL=0 --temperature=25 + --energyAdd=1.2 --outMode=C --outCsvCols=hybridDB,E,ED1,ED2,E_init,E_loops,E_dangleL,E_dangleR,E_endL,E_endR,E_add + --threads=1 --default-log-file=/dev/null) +for dangles in '' --energyNoDangles; do + extra=() + if test -n "$dangles"; then extra+=("$dangles"); fi + "$bin" "${thermo[@]}" "${extra[@]}" --noSeed --model=S --mode=M -n 5 > "$tmp/predicted" + test "$(wc -l < "$tmp/predicted")" -gt 1 + structures=$(awk -F';' 'NR>1 {printf "%s%s", sep, $1; sep=":"}' "$tmp/predicted") + "$bin" "${thermo[@]}" "${extra[@]}" --rri="$structures" > "$tmp/evaluated" + cmp "$tmp/predicted" "$tmp/evaluated" + # All ordinary report formats must accept the evaluated structures. + for mode in N D; do + "$bin" --target=AGCGACGCA --query=UGCGUCGCU --rri="$structures" --outMode="$mode" > "$tmp/text" + test -s "$tmp/text" + done +done + +expect_error() { + local status=0 + "$bin" "$@" > "$tmp/bad.out" 2> "$tmp/bad.err" || status=$? + test "$status" -eq 1 || test "$status" -eq 255 +} +expect_error "${common[@]}" "${csv[@]}" --personality=IntaRNAeval +for bad in '' '1|&1||' '1...&1...' '0|&1|' '6||&1||' '1|&1|:' '1|&1|::1|&1|' '9223372036854775808|&1|' '1((&1))'; do + expect_error "${common[@]}" "${csv[@]}" --rri "$bad" +done +expect_error --target=AAAA --query=CCCC --rri='1|&1|' +expect_error --target=GGGG --query=NCCC --rri='1|&1|' +expect_error --target=GGGG --query=CCCC --tIdxPos0=-2 --rri='0|&1|' +printf '>one\nGGGGGG\n>two\nGGGGGG\n' > "$tmp/multiple.fa" +expect_error "${common[@]:1}" --target="$tmp/multiple.fa" --rri="$rri" +expect_error --target=CCCCCC --query="$tmp/multiple.fa" --rri="$rri" +# Retain the selected accessibility model: unavailable intervals fail clearly. +expect_error --target=GGGGGGGG --query=CCCCCCCC --accW=3 --accL=3 --rri='1||||||||&1||||||||' +grep -q 'accessibility' "$tmp/bad.err" "$tmp/bad.out" +echo 'IntaRNAeval CLI checks passed'