Add IntaRNAeval for predefined RNA–RNA interaction energies

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

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

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 \

src/IntaRNA/PredictorEvalOnly.cpp

@@ -0,0 +1,162 @@
+#include "IntaRNA/PredictorEvalOnly.h"
+
+#include <charconv>
+#include <set>
+#include <stdexcept>
+#include <string_view>
+
+namespace IntaRNA {
+
+namespace {
+
+// Parse one strand without performing arithmetic on unvalidated external indices.
+std::vector<size_t> 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<size_t> paired;
+	for (size_t p=0; p<structure.size(); ++p) {
+		if (structure[p] == '|') paired.push_back(offset+p);
+	}
+	return paired;
+}
+
+} // namespace
+
+std::vector<Interaction>
+PredictorEvalOnly::parseInteractions( const std::string & encoding,
+		const RnaSequence & target, const RnaSequence & query )
+{
+	std::vector<Interaction> 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<paired1.size(); ++p) {
+			const size_t q = paired2[paired2.size()-1-p];
+			if (!RnaSequence::areComplementary(target, query, paired1[p], q)) {
+				throw std::invalid_argument("--rri: non-complementary base pair at target "
+						+toString(target.getInOutIndex(paired1[p]))+", query "+toString(query.getInOutIndex(q)));
+			}
+			interaction.basePairs.emplace_back(paired1[p], q);
+		}
+		result.push_back(interaction);
+		if (end == std::string::npos) break;
+		start = end+1;
+	} while (true);
+	return result;
+}
+
+PredictorEvalOnly::PredictorEvalOnly( const InteractionEnergy & energy,
+		OutputHandler & output, PredictionTracker * predTracker,
+		const std::vector<Interaction> & 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<Interaction::PairingVec> 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; p<interaction.basePairs.size(); ++p) {
+			const auto & left = interaction.basePairs[p-1];
+			const auto & right = interaction.basePairs[p];
+			const E_type loopE = energy.getE_interLeft(energy.getIndex1(left), energy.getIndex1(right),
+					energy.getIndex2(left), energy.getIndex2(right));
+			if (E_isINF(loopE)) {
+				throw std::runtime_error("PredictorEvalOnly: interaction loop cannot be evaluated by the selected energy model");
+			}
+			hybridE += loopE;
+		}
+		const auto & left = interaction.basePairs.front();
+		const auto & right = interaction.basePairs.back();
+		interaction.energy = energy.getE(energy.getIndex1(left), energy.getIndex1(right),
+				energy.getIndex2(left), energy.getIndex2(right), hybridE);
+		if (E_isINF(interaction.energy)) {
+			throw std::runtime_error("PredictorEvalOnly: interaction cannot be evaluated with the selected accessibility model; check accessibility windows, constraints and input coverage");
+		}
+	}
+	for (const auto & interaction : interactions) {
+		const auto & left = interaction.basePairs.front();
+		const auto & right = interaction.basePairs.back();
+		updateOptima(energy.getIndex1(left), energy.getIndex1(right),
+				energy.getIndex2(left), energy.getIndex2(right), interaction.energy, false, true);
+	}
+	reportOptima();
+}
+
+void
+PredictorEvalOnly::initOptima()
+{
+	Zall = 0;
+}
+
+void
+PredictorEvalOnly::updateOptima( const size_t i1, const size_t j1,
+		const size_t i2, const size_t j2, const E_type interactionE,
+		const bool isHybridE, const bool incrementZ )
+{
+	const E_type totalE = isHybridE ? energy.getE(i1, j1, i2, j2, interactionE) : interactionE;
+	if (incrementZ) incrementZall(energy.getBoltzmannWeight(totalE));
+	if (predTracker != NULL) predTracker->updateOptimumCalled(i1, j1, i2, j2, totalE);
+}
+
+void
+PredictorEvalOnly::reportOptima()
+{
+	output.incrementZ(Zall);
+	for (const auto & interaction : interactions) output.add(interaction);
+}
+
+} // namespace IntaRNA

src/IntaRNA/PredictorEvalOnly.h

@@ -0,0 +1,72 @@
+#ifndef INTARNA_PREDICTOREVALONLY_H_
+#define INTARNA_PREDICTOREVALONLY_H_
+
+#include "IntaRNA/Predictor.h"
+
+#include <string>
+#include <vector>
+
+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<Interaction> & 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<Interaction> 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<Interaction> interactions;
+};
+
+} // namespace IntaRNA
+
+#endif