Module bussilab.rna2d
Module containing an RNA secondary-structure model with continuous pairing penalties.
See Molecule.
Functions
def reset_default_parameters()-
Expand source code
def reset_default_parameters(): """Restore the original thermodynamic defaults for new molecules.""" with _DEFAULT_PARAMETERS_LOCK: _default_parameters.clear() _default_parameters.update(_INITIAL_DEFAULT_PARAMETERS)Restore the original thermodynamic defaults for new molecules.
def sample_to_numpy(samples, *, deduplicate=True)-
Expand source code
def sample_to_numpy(samples, *, deduplicate=True): """ Convert sampled dot-bracket structures to NumPy arrays. The input may be either a sequence of structure strings, as returned by ``Molecule.sample(weights=False)``, or a sequence of ``(structure, log_weight)`` pairs, as returned by ``Molecule.sample(weights=True)``. Mixed inputs are rejected. Structures are represented by zero-based pair tables: ``states[k, i]`` is the index paired with nucleotide ``i`` in structure ``k``, or ``-1`` when it is unpaired. This is a NumPy-oriented adaptation of ViennaRNA's pair table format. A row ``table`` can be converted to ViennaRNA's convention with ``np.concatenate(([len(table)], table + 1))``: ViennaRNA prepends the sequence length, uses one-based partner indices, and uses zero for an unpaired nucleotide. Parameters ---------- samples : sequence of str or sequence of (str, float) Unweighted structures or structures with unnormalized log weights. deduplicate : bool, default=True If True, merge identical structures. Unweighted occurrences are combined through their counts; supplied log weights are combined using logarithmic addition. If False, retain every occurrence. Returns ------- states : ndarray of int16, shape (n_structures, sequence_length) Zero-based pair tables, with ``-1`` denoting an unpaired nucleotide. logweights : ndarray of float, shape (n_structures,) Normalized log weights, satisfying ``sum(exp(logweights)) == 1`` up to floating-point precision. """ if not isinstance(deduplicate, (bool, np.bool_)): raise ValueError("deduplicate must be a boolean") samples = list(samples) if not samples: raise ValueError("samples cannot be empty") unweighted = all(isinstance(item, str) for item in samples) weighted = all( not isinstance(item, str) and hasattr(item, "__len__") and len(item) == 2 and isinstance(item[0], str) for item in samples ) if not (unweighted or weighted): raise ValueError( "samples must contain either strings or (structure, log_weight) " "pairs, without mixing the two forms" ) if unweighted: structures = samples input_logweights = np.zeros(len(samples), dtype=float) else: structures = [item[0] for item in samples] try: input_logweights = np.asarray( [item[1] for item in samples], dtype=float, ) except (TypeError, ValueError) as error: raise ValueError("log weights must be real numbers") from error if input_logweights.ndim != 1: raise ValueError("log weights must be scalar") if not np.all(np.isfinite(input_logweights)): raise ValueError("log weights must be finite") if deduplicate: unique_structures = [] unique_logweights = [] indices = {} for structure, logweight in zip(structures, input_logweights): if structure in indices: index = indices[structure] unique_logweights[index] = np.logaddexp( unique_logweights[index], logweight, ) else: indices[structure] = len(unique_structures) unique_structures.append(structure) unique_logweights.append(float(logweight)) structures = unique_structures input_logweights = np.asarray(unique_logweights, dtype=float) states = _structures_to_pairtables(structures) return states, _normalize_logweights(input_logweights)Convert sampled dot-bracket structures to NumPy arrays.
The input may be either a sequence of structure strings, as returned by
Molecule.sample(weights=False), or a sequence of(structure, log_weight)<code> pairs, as returned by </code>Molecule.sample(weights=True). Mixed inputs are rejected.Structures are represented by zero-based pair tables:
states[k, i]is the index paired with nucleotideiin structurek, or-1when it is unpaired. This is a NumPy-oriented adaptation of ViennaRNA's pair table format. A rowtablecan be converted to ViennaRNA's convention withnp.concatenate(([len(table)], table + 1)): ViennaRNA prepends the sequence length, uses one-based partner indices, and uses zero for an unpaired nucleotide.Parameters
samples:sequenceofstrorsequenceof(str, float)- Unweighted structures or structures with unnormalized log weights.
deduplicate:bool, default=True- If True, merge identical structures. Unweighted occurrences are combined through their counts; supplied log weights are combined using logarithmic addition. If False, retain every occurrence.
Returns
states:ndarrayofint16, shape (n_structures, sequence_length)- Zero-based pair tables, with
-1denoting an unpaired nucleotide. logweights:ndarrayoffloat, shape (n_structures,)- Normalized log weights, satisfying
sum(exp(logweights)) == 1up to floating-point precision.
def set_default_parameters(*, temperature=None, no_lonely_pair=None, pf_smooth=None, NaCl=None, parameters=None)-
Expand source code
def set_default_parameters( *, temperature=None, no_lonely_pair=None, pf_smooth=None, NaCl=None, parameters=None, ): """ Update the default thermodynamic parameters for new molecules. Parameters set to ``None`` are left unchanged. Existing molecules are not affected. Use :func:`reset_default_parameters` to restore all original defaults, including ViennaRNA's default salt concentration. """ updates = {} if temperature is not None: if not temperature >= 0.0: raise ValueError( f"Temperature {temperature} should be positive" ) updates["temperature"] = temperature if no_lonely_pair is not None: if not isinstance(no_lonely_pair, (bool, np.bool_)): raise ValueError("no_lonely_pair must be a boolean") updates["no_lonely_pair"] = bool(no_lonely_pair) if pf_smooth is not None: if not isinstance(pf_smooth, (bool, np.bool_)): raise ValueError("pf_smooth must be a boolean") updates["pf_smooth"] = bool(pf_smooth) if NaCl is not None: if not NaCl >= 0.0: raise ValueError( f"Salt concentration {NaCl} should be positive" ) updates["NaCl"] = NaCl if parameters is not None: parameters = str(parameters).lower() if parameters not in _THERMODYNAMIC_PARAMETERS: raise ValueError( f"Thermodynamic parameters {parameters} not known" ) updates["parameters"] = parameters with _DEFAULT_PARAMETERS_LOCK: _default_parameters.update(updates)Update the default thermodynamic parameters for new molecules.
Parameters set to
Noneare left unchanged. Existing molecules are not affected. Use :func:reset_default_parameters()to restore all original defaults, including ViennaRNA's default salt concentration. def suboptimal_to_numpy(suboptimal, temperature)-
Expand source code
def suboptimal_to_numpy(suboptimal, temperature): """ Convert suboptimal dot-bracket structures and energies to NumPy arrays. Structures use the same zero-based pair-table representation documented by :func:`sample_to_numpy`. Energies are converted to Boltzmann log weights at the explicitly supplied temperature. Parameters ---------- suboptimal : sequence of (str, float) Structure and energy pairs, in kcal/mol, as returned by ``Molecule.suboptimal_structures()``. temperature : float Temperature in kelvin. Returns ------- states : ndarray of int16, shape (n_structures, sequence_length) Zero-based pair tables, with ``-1`` denoting an unpaired nucleotide. logweights : ndarray of float, shape (n_structures,) Normalized Boltzmann log weights, satisfying ``sum(exp(logweights)) == 1`` up to floating-point precision. """ suboptimal = list(suboptimal) if not suboptimal: raise ValueError("suboptimal cannot be empty") if not all( not isinstance(item, str) and hasattr(item, "__len__") and len(item) == 2 and isinstance(item[0], str) for item in suboptimal ): raise ValueError( "suboptimal must contain (structure, energy) pairs" ) temperature = float(temperature) if not np.isfinite(temperature) or temperature <= 0.0: raise ValueError("temperature must be finite and positive") structures = [item[0] for item in suboptimal] try: energies = np.asarray([item[1] for item in suboptimal], dtype=float) except (TypeError, ValueError) as error: raise ValueError("energies must be real numbers") from error if energies.ndim != 1: raise ValueError("energies must be scalar") if not np.all(np.isfinite(energies)): raise ValueError("energies must be finite") states = _structures_to_pairtables(structures) logweights = -energies / (_KB * temperature) if not np.all(np.isfinite(logweights)): raise ValueError("energies are too large to convert to log weights") return states, _normalize_logweights(logweights)Convert suboptimal dot-bracket structures and energies to NumPy arrays.
Structures use the same zero-based pair-table representation documented by :func:
sample_to_numpy(). Energies are converted to Boltzmann log weights at the explicitly supplied temperature.Parameters
suboptimal:sequenceof(str, float)- Structure and energy pairs, in kcal/mol, as returned by
Molecule.suboptimal_structures(). temperature:float- Temperature in kelvin.
Returns
states:ndarrayofint16, shape (n_structures, sequence_length)- Zero-based pair tables, with
-1denoting an unpaired nucleotide. logweights:ndarrayoffloat, shape (n_structures,)- Normalized Boltzmann log weights, satisfying
sum(exp(logweights)) == 1up to floating-point precision.
Classes
class Molecule (seq: str,
*,
lambdas1d=None,
temperature=None,
force_paired=None,
force_unpaired=None,
state_positions=None,
state_biases=None,
reduce_state_space=True,
no_lonely_pair=None,
pf_smooth=None,
NaCl=None,
parameters=None)-
Expand source code
class Molecule: """ RNA secondary-structure model with continuous pairing penalties. The class wraps one or more ViennaRNA dynamic-programming ensembles and supports continuous per-nucleotide pairing penalties. Multiple ensembles may be used to assign arbitrary energy biases to the paired/unpaired states of selected nucleotides. A penalty λᵢ is added whenever nucleotide *i* is paired. Internally, these penalties are automatically represented as an equivalent hybrid combination of unpaired and pair soft constraints. This representation avoids numerical overflows in partition-function calculations while preserving the requested thermodynamic model. Partition-function calculations use the exact continuous penalties. Minimum- free-energy and suboptimal structure prediction use ViennaRNA's rounded soft constraints to generate candidate structures, which are then rescored using the exact continuous penalties. Parameters ---------- seq : str RNA sequence. lambdas1d : array-like, optional Per-nucleotide pairing penalties (kcal/mol). Positive values penalize pairing, whereas negative values favor pairing. If omitted, all penalties are zero. force_paired : array-like of int, optional Zero-based indices of nucleotides that are required to be paired, without specifying their pairing partners. force_unpaired : array-like of int, optional Zero-based indices of nucleotides that are required to be unpaired. state_positions : array-like of int or sequence of array-like, optional Zero-based indices defining binary paired/unpaired states. A flat sequence defines one state set. A sequence of sequences defines multiple disjoint state sets whose energy biases are additive. If `state_biases` is omitted, every state is assigned zero bias. state_biases : array-like or sequence of array-like, optional Energy biases (kcal/mol) for the states defined by `state_positions`. For one state set, its shape must be `(2,) * len(state_positions)`. For multiple sets, provide one tensor per set, with shape `(2,) * len(positions)`. Index zero denotes an unpaired nucleotide and index one a paired nucleotide. It cannot be provided without `state_positions`. reduce_state_space : bool, default=True If True, represent the final selected state position as an equivalent 1D pairing penalty, reducing the number of dynamic-programming ensembles from 2**N to 2**(N-1). If False, use one hard-conditioned ensemble for every state. no_lonely_pair : bool or None, default=None If True, exclude structures containing isolated base pairs using ViennaRNA's `noLP` model option. If None, use the current module default. pf_smooth : bool or None, default=None Whether ViennaRNA should smooth energies in partition-function Boltzmann factors. If None, use the current module default, which is initially False. If True, partition functions and sampling use ViennaRNA's smoothed model, while MFE, evaluation, and suboptimal structures continue to use its unsmoothed energy model. temperature : float or None, default=None Temperature in kelvin. If None, use the current module default, which is initially 310.15 K. NaCl : float or None, default=None Sodium concentration (M). If None, use the current module default, which initially selects ViennaRNA's default value. parameters : {"turner1999", "turner2004", "andronescu2007", "langdon2018"} or None Thermodynamic parameter set. If None, use the current module default. Notes ----- By default, a state set containing N positions uses 2**(N-1) dynamic-programming ensembles. For multiple sets of sizes N_k, the number is 2**(sum(N_k)-K). The final position of each set is represented within each ensemble by an equivalent constant energy shift and 1D pairing penalty. Set `reduce_state_space=False` to explicitly condition every selected position. Default ViennaRNA builds round soft constraints to the nearest 0.01 kcal/mol in partition-function calculations. When continuous soft constraints are not natively supported, this class automatically applies a lightweight Python callback to recover the exact continuous model. """ def __init__( self, seq: str, *, lambdas1d=None, temperature=None, force_paired = None, force_unpaired = None, state_positions=None, state_biases=None, reduce_state_space=True, no_lonely_pair=None, pf_smooth=None, NaCl=None, parameters=None, ): ( temperature, no_lonely_pair, pf_smooth, NaCl, parameters, ) = _resolve_default_parameters( temperature, no_lonely_pair, pf_smooth, NaCl, parameters, ) self._has_state_biases = state_positions is not None if not isinstance(reduce_state_space, (bool, np.bool_)): raise ValueError("reduce_state_space must be a boolean") self._reduce_state_space = bool(reduce_state_space) if not isinstance(no_lonely_pair, (bool, np.bool_)): raise ValueError("no_lonely_pair must be a boolean") self._no_lonely_pair = bool(no_lonely_pair) if not isinstance(pf_smooth, (bool, np.bool_)): raise ValueError("pf_smooth must be a boolean") self._pf_smooth = bool(pf_smooth) if state_positions is None and state_biases is not None: raise ValueError( "state_biases cannot be provided without state_positions" ) base_force_paired = ( [] if force_paired is None else list(force_paired) ) base_force_unpaired = ( [] if force_unpaired is None else list(force_unpaired) ) if state_positions is None: self._state_sets_were_nested = False self._state_position_sets = () self._state_bias_sets = () else: try: raw_positions = list(state_positions) except TypeError as error: raise ValueError( "state_positions must be a sequence" ) from error # A flat sequence retains the original single-state-set syntax. # Nesting is detected from state_positions rather than state_biases, # since a two-element bias vector is inherently ambiguous. if all(np.isscalar(value) for value in raw_positions): self._state_sets_were_nested = False raw_position_sets = [raw_positions] else: self._state_sets_were_nested = True raw_position_sets = raw_positions position_sets = [] for positions in raw_position_sets: positions = np.asarray(positions) if positions.ndim != 1: raise ValueError( "each state_positions set must be one-dimensional" ) if ( positions.size and not np.issubdtype(positions.dtype, np.integer) ): raise ValueError( "state_positions must contain integers" ) position_set = tuple(int(i) for i in positions) if len(set(position_set)) != len(position_set): raise ValueError( "state_positions sets must not contain duplicates" ) if any( i < 0 or i >= len(str(seq)) for i in position_set ): raise ValueError( "state_positions must contain valid nucleotide indices" ) position_sets.append(position_set) all_state_positions = [ position for position_set in position_sets for position in position_set ] if len(set(all_state_positions)) != len(all_state_positions): raise ValueError("state_positions sets must be disjoint") fixed_positions = ( set(base_force_paired) | set(base_force_unpaired) ) if fixed_positions.intersection(all_state_positions): raise ValueError( "state_positions must not overlap force_paired or " "force_unpaired" ) expected_shapes = [ (2,) * len(position_set) for position_set in position_sets ] if state_biases is None: bias_sets = [ np.zeros(shape, dtype=float) for shape in expected_shapes ] elif len(position_sets) == 1: expected_shape = expected_shapes[0] try: candidate = np.asarray(state_biases, dtype=float) except (TypeError, ValueError): candidate = np.asarray([], dtype=float) if candidate.shape == expected_shape: bias_sets = [candidate.copy()] else: try: if len(state_biases) != 1: raise ValueError candidate = np.asarray( state_biases[0], dtype=float, ) except (TypeError, ValueError, IndexError) as error: raise ValueError( f"state_biases must have shape {expected_shape}" ) from error if candidate.shape != expected_shape: raise ValueError( f"state_biases must have shape {expected_shape}" ) bias_sets = [candidate.copy()] else: try: if len(state_biases) != len(position_sets): raise ValueError( "state_biases must contain one tensor per " "state_positions set" ) except TypeError as error: raise ValueError( "state_biases must contain one tensor per " "state_positions set" ) from error bias_sets = [] for biases, expected_shape in zip( state_biases, expected_shapes, ): biases = np.asarray(biases, dtype=float) if biases.shape != expected_shape: raise ValueError( "each state_biases tensor must have shape " f"{expected_shape}" ) bias_sets.append(biases.copy()) if any( not np.all(np.isfinite(biases)) for biases in bias_sets ): raise ValueError( "state_biases must contain only finite values" ) self._state_position_sets = tuple(position_sets) self._state_bias_sets = tuple(bias_sets) # Preserve the input representation for compatibility. New code uses # the normalized plural attributes above. if not self._state_position_sets: self._state_positions = () self._state_biases = np.zeros((), dtype=float) elif self._state_sets_were_nested: self._state_positions = self._state_position_sets self._state_biases = self._state_bias_sets else: self._state_positions = self._state_position_sets[0] self._state_biases = self._state_bias_sets[0] if lambdas1d is None: base_lambdas1d = np.zeros(len(str(seq))) else: base_lambdas1d = np.asarray( lambdas1d, dtype=float, ).copy() if base_lambdas1d.ndim != 1: raise ValueError("lambdas1d must be one-dimensional") if len(base_lambdas1d) != len(str(seq)): raise ValueError( "lambdas1d must contain one value per nucleotide" ) if not np.all(np.isfinite(base_lambdas1d)): raise ValueError( "lambdas1d must contain only finite values" ) # The final variable of every state set does not require explicit # branching. For each assignment of the preceding variables, its two # biases are equivalent to a constant plus a 1D pairing penalty. explicit_position_sets = [] implicit_positions = [] state_slices = [] offset = 0 for position_set in self._state_position_sets: if position_set and self._reduce_state_space: explicit_positions = position_set[:-1] implicit_position = position_set[-1] else: explicit_positions = position_set implicit_position = None explicit_position_sets.append(explicit_positions) implicit_positions.append(implicit_position) state_slices.append(slice( offset, offset + len(explicit_positions), )) offset += len(explicit_positions) self._explicit_state_position_sets = tuple(explicit_position_sets) self._implicit_state_positions = tuple(implicit_positions) self._state_slices = tuple(state_slices) explicit_state_positions = tuple( position for position_set in explicit_position_sets for position in position_set ) self._explicit_state_positions = explicit_state_positions component_shape = (2,) * len(explicit_state_positions) # Preserve the legacy singular attribute for one state set. self._implicit_state_position = ( implicit_positions[0] if len(implicit_positions) == 1 else None ) self._states = list(np.ndindex(component_shape)) self._component_biases = np.empty(component_shape, dtype=float) self._dp_molecules = [] for state in self._states: state_paired = [ position for position, value in zip( explicit_state_positions, state, ) if value ] state_unpaired = [ position for position, value in zip( explicit_state_positions, state, ) if not value ] component_lambdas1d = base_lambdas1d.copy() component_bias = 0.0 for biases, state_slice, implicit_position in zip( self._state_bias_sets, self._state_slices, self._implicit_state_positions, ): set_state = state[state_slice] if implicit_position is None: component_bias += float(biases[set_state]) else: unpaired_bias = float(biases[set_state + (0,)]) paired_bias = float(biases[set_state + (1,)]) component_bias += unpaired_bias component_lambdas1d[implicit_position] += ( paired_bias - unpaired_bias ) self._component_biases[state] = component_bias self._dp_molecules.append( _DPMolecule( seq, lambdas1d=component_lambdas1d, temperature=temperature, NaCl=NaCl, force_paired=base_force_paired + state_paired, force_unpaired=base_force_unpaired + state_unpaired, no_lonely_pair=self._no_lonely_pair, pf_smooth=self._pf_smooth, parameters=parameters, ) ) first_molecule = self._dp_molecules[0] self._seq = first_molecule._seq self._lambdas1d = base_lambdas1d self._temperature = first_molecule._temperature self._salt = first_molecule._salt self._parameters = first_molecule._parameters self._no_lonely_pair = first_molecule._no_lonely_pair self._pf_smooth = first_molecule._pf_smooth self._force_paired = tuple(base_force_paired) self._force_unpaired = tuple(base_force_unpaired) def _condition_unpaired(self, position): """ Return an equivalent molecule conditioned on one nucleotide being unpaired. """ if position in self._force_paired: return None if position in self._force_unpaired: return self force_unpaired = self._force_unpaired + (position,) state_position_sets = list(self._state_position_sets) state_bias_sets = list(self._state_bias_sets) for set_index, position_set in enumerate(state_position_sets): if position not in position_set: continue axis = position_set.index(position) state_position_sets[set_index] = ( position_set[:axis] + position_set[axis + 1:] ) state_bias_sets[set_index] = np.take( state_bias_sets[set_index], 0, axis=axis, ) break if not state_position_sets: state_positions = None state_biases = None elif not self._state_sets_were_nested: state_positions = state_position_sets[0] state_biases = state_bias_sets[0] else: state_positions = state_position_sets state_biases = state_bias_sets return Molecule( self._seq, lambdas1d=self._lambdas1d, temperature=self._temperature, force_paired=self._force_paired, force_unpaired=force_unpaired, state_positions=state_positions, state_biases=state_biases, reduce_state_space=self._reduce_state_space, no_lonely_pair=self._no_lonely_pair, pf_smooth=self._pf_smooth, NaCl=self._salt, parameters=self._parameters, ) def _component_probabilities(self): """ Return the total free energy and normalized component probabilities. """ free_energies = np.array([ molecule.total_free_energy() for molecule in self._dp_molecules ]) biased_free_energies = ( free_energies + self._component_biases.ravel() ) finite = np.isfinite(biased_free_energies) if not np.any(finite): raise RuntimeError( "No state has a finite partition function" ) reference = np.min(biased_free_energies[finite]) kT = _KB * self._dp_molecules[0]._temperature relative_weights = np.zeros(len(self._dp_molecules)) relative_weights[finite] = np.exp( -(biased_free_energies[finite] - reference) / kT ) normalization = math.fsum(relative_weights) probabilities = relative_weights / normalization total_free_energy = reference - kT * math.log(normalization) return float(total_free_energy), probabilities def _component_mfes(self): """ Return component MFE structures and their biased energies. """ structures = [] energies = [] for molecule, bias in zip( self._dp_molecules, self._component_biases.ravel(), ): structure, energy = molecule.mfe() structures.append(structure) if molecule._satisfies_hard_constraints(structure): energies.append(energy + bias) else: energies.append(math.inf) energies = np.asarray(energies) if not np.any(np.isfinite(energies)): raise RuntimeError( "No state has a feasible MFE structure" ) return structures, energies def mfe(self): """ Return the minimum-free-energy structure. The returned energy always corresponds to the exact continuous pairing penalties, even when ViennaRNA internally rounds soft constraints. Returns ------- structure : str Dot-bracket representation of the MFE structure. energy : float Exact free energy (kcal/mol). """ # A single DP component needs no mixture-level search. if len(self._dp_molecules) == 1: structure, energy = self._dp_molecules[0].mfe() return ( structure, float(energy + self._component_biases.item()), ) structures, energies = self._component_mfes() index = int(np.argmin(energies)) return structures[index], float(energies[index]) def evaluate(self, structure): """ Return the exact free energy of a secondary structure. Parameters ---------- structure : str Dot-bracket representation of a secondary structure. Returns ------- float Exact free energy (kcal/mol), including continuous pairing penalties and state biases. """ if not isinstance(structure, str): raise ValueError("structure must be a string") if len(structure) != len(self._seq): raise ValueError( "structure must contain one symbol per nucleotide" ) if structure.strip(".()"): raise ValueError( "structure must use '.', '(', and ')' dot-bracket symbols" ) # Components follow np.ndindex order, so the paired/unpaired state is # also the component index interpreted as a binary integer. Avoid # constructing a state tuple on this latency-sensitive path. component_index = 0 for position in self._explicit_state_positions: component_index = ( 2 * component_index + (structure[position] != ".") ) energy = self._dp_molecules[component_index].evaluate(structure) return float( energy + self._component_biases.flat[component_index] ) def base_pairing_probability(self): """ Return the base-pairing probability matrix. Returns ------- ndarray Symmetric NxN matrix whose element (i,j) is the equilibrium probability that nucleotides i and j form a base pair. """ # Compute BPPs before mixture weights: a BPP calculation also produces # the component free energy, while doing this in the opposite order # would require rebuilding each PF without probability backtracking. for molecule in self._dp_molecules: molecule._ensure_pf(compute_bpp=True) if len(self._dp_molecules) == 1: return self._dp_molecules[0]._base_pairing_probability.copy() _, probabilities = self._component_probabilities() matrix = np.zeros_like( self._dp_molecules[0]._base_pairing_probability ) for probability, molecule in zip( probabilities, self._dp_molecules, ): if probability == 0.0: continue matrix += ( probability * molecule._base_pairing_probability ) return matrix def total_free_energy(self): """ Return the ensemble free energy. Returns ------- float Ensemble free energy (kcal/mol) corresponding to the exact continuous pairing penalties. """ total_free_energy, _ = self._component_probabilities() return total_free_energy def d_free_energy_d_lambdas1d(self): """ Return the derivatives of the free energy with respect to `lambdas1d`. The derivative for nucleotide ``i`` is its equilibrium pairing probability, ``dF / d lambda_i = <s_i>``. Returns ------- ndarray One-dimensional array containing one derivative per nucleotide. """ return np.sum(self.base_pairing_probability(), axis=1) def d_free_energy_d_state_biases(self): """ Return the derivatives of the free energy with respect to state biases. The derivative with respect to a state's energy bias is the equilibrium probability of that state. Returns ------- ndarray or list of ndarray With the flat, single-set `state_positions` syntax, an array with the same shape as `state_biases`. With the list-of-sets syntax, one array per set, including when the list contains only one set. Raises ------ ValueError If this molecule was not initialized with state biases. """ if not self._has_state_biases: raise ValueError( "This molecule was not initialized with state biases" ) needs_pairing_probabilities = any( position is not None for position in self._implicit_state_positions ) if needs_pairing_probabilities: component_pairing_probabilities = [ np.sum(molecule.base_pairing_probability(), axis=1) for molecule in self._dp_molecules ] else: component_pairing_probabilities = [None] * len( self._dp_molecules ) _, component_probabilities = self._component_probabilities() probability_sets = [ np.zeros_like(biases) for biases in self._state_bias_sets ] for ( state, component_probability, pairing_probabilities, ) in zip( self._states, component_probabilities, component_pairing_probabilities, ): if component_probability == 0.0: continue for probabilities, state_slice, implicit_position in zip( probability_sets, self._state_slices, self._implicit_state_positions, ): set_state = state[state_slice] if implicit_position is None: probabilities[set_state] += component_probability continue paired_probability = float(np.clip( pairing_probabilities[implicit_position], 0.0, 1.0, )) probabilities[set_state + (0,)] += ( component_probability * (1.0 - paired_probability) ) probabilities[set_state + (1,)] += ( component_probability * paired_probability ) if not self._state_sets_were_nested: return probability_sets[0] return probability_sets def suboptimal_structures(self, delta): """ Enumerate suboptimal secondary structures. Candidate structures are generated using ViennaRNA's rounded soft constraints, rescored with the exact continuous penalties, and returned sorted by exact energy. For state mixtures, the MFE of every feasible state is calculated first. Each component is then enumerated with the local energy window required by the common, globally biased cutoff. Parameters ---------- delta : float Maximum energy difference (kcal/mol) above the exact MFE. Returns ------- list of (str, float) List of (structure, energy) pairs sorted by increasing exact energy. """ # A single DP component needs no mixture-level enumeration. if len(self._dp_molecules) == 1: suboptimal = self._dp_molecules[0].suboptimal_structures( delta ) bias = float(self._component_biases.item()) if bias == 0.0: return suboptimal return [ (structure, float(energy + bias)) for structure, energy in suboptimal ] delta = float(delta) if not np.isfinite(delta) or delta < 0.0: raise ValueError("delta must be finite and non-negative") _, component_mfe_energies = self._component_mfes() global_mfe_energy = float(np.min(component_mfe_energies)) cutoff = global_mfe_energy + delta merged = [] for molecule, bias, component_mfe_energy in zip( self._dp_molecules, self._component_biases.ravel(), component_mfe_energies, ): if not np.isfinite(component_mfe_energy): continue local_delta = cutoff - component_mfe_energy if local_delta < 0.0: continue for structure, energy in molecule.suboptimal_structures( local_delta ): biased_energy = float(energy + bias) if biased_energy <= cutoff: merged.append((structure, biased_energy)) merged.sort(key=lambda item: item[1]) return merged def sample(self, number, weights=False): """ Generate Boltzmann-distributed secondary structures. By default, rejection sampling corrects ViennaRNA's rounded soft constraints and returns unweighted structures from the exact continuous model. Weighted samples from the rounded model can be requested instead. Parameters ---------- number : int Number of structures to sample. weights : bool, default=False If False, return unweighted structures sampled from the exact model. If True, return rounded-model samples and their log-weight corrections. The weighted path is faster when rounding corrections are present. Returns ------- list of str or list of (str, float) With ``weights=False``, a list of dot-bracket structures. With ``weights=True``, each element contains a structure and its unnormalized log-weight correction log(w) = -(E_exact - E_rounded) / (k_B T). For state mixtures, the log-weight also includes the component-specific exact/rounded normalization correction. This correction is caused only by lambda rounding, not by the state bias. When all lambdas are multiples of 0.01 kcal/mol, every returned log-weight is zero. """ number = int(number) if number <= 0: raise ValueError("number must be a positive integer") if not isinstance(weights, (bool, np.bool_)): raise ValueError("weights must be a boolean") # A single DP component needs no mixture-level sampling. if len(self._dp_molecules) == 1: return self._dp_molecules[0].sample( number, weights=weights, ) _, probabilities = self._component_probabilities() component_indices = np.random.choice( len(self._dp_molecules), size=number, p=probabilities, ) result = [None] * number for component_index, molecule in enumerate(self._dp_molecules): output_indices = np.flatnonzero( component_indices == component_index ) if not len(output_indices): continue samples = molecule.sample( len(output_indices), weights=weights, ) if weights: rounding_correction = ( molecule.sample_rounding_correction() ) for output_index, (structure, log_weight) in zip( output_indices, samples, ): result[output_index] = ( structure, float(log_weight + rounding_correction), ) else: for output_index, structure in zip( output_indices, samples, ): result[output_index] = structure return result def suboptimal_coverage(self, delta): """ Return the fraction of the partition function represented by the suboptimal ensemble. The suboptimal ensemble contains all structures whose exact energy is within ``delta`` kcal/mol of the minimum-free-energy structure. Parameters ---------- delta : float Maximum energy difference (kcal/mol) above the exact MFE. Returns ------- float Fraction of the total partition function represented by the enumerated structures. The result lies between zero and one, apart from possible small numerical errors. """ suboptimal = self.suboptimal_structures(delta) reference_energy = suboptimal[0][1] inverse_kT = 1.0 / (_KB * self._temperature) relative_partition_function = math.fsum( math.exp(-(energy - reference_energy) * inverse_kT) for _, energy in suboptimal ) log_coverage = ( math.log(relative_partition_function) + (self.total_free_energy() - reference_energy) * inverse_kT ) coverage = float(math.exp(log_coverage)) if coverage > 1.0: if coverage <= 1.0 + 1e-5: coverage = 1.0 else: raise RuntimeError( f"Suboptimal coverage is unexpectedly larger than one: " f"{coverage}" ) return coverage def _pairing_correlation_matrix_pf(self): """ Return the exact joint pairing-probability matrix. """ p = np.sum(self.base_pairing_probability(), axis=1) p = np.clip(p, 0.0, 1.0) p[np.asarray(self._force_paired, dtype=int)] = 1.0 p[np.asarray(self._force_unpaired, dtype=int)] = 0.0 n = len(self._seq) conditional = np.empty((n, n)) for i in range(n): if p[i] == 0.0: conditional[i, :] = p continue if p[i] == 1.0: conditional[i, :] = 0.0 continue conditioned = self._condition_unpaired(i) conditional[i, :] = np.sum( conditioned.base_pairing_probability(), axis=1, ) matrix = np.empty((n, n)) for i in range(n): matrix[i, :] = p - (1.0 - p[i]) * conditional[i, :] # Since s_i**2 = s_i, the diagonal must equal P(s_i = 1). np.fill_diagonal(matrix, p) # Remove small numerical asymmetries from independent constrained PF runs. return 0.5 * (matrix + matrix.T) def pairing_correlation_matrix(self): """ Return the joint pairing-probability matrix. Element ``(i, j)`` is the probability that nucleotides ``i`` and ``j`` are simultaneously paired, irrespective of their pairing partners. The matrix is computed exactly using constrained partition-function calculations. Returns ------- ndarray Symmetric NxN matrix whose diagonal contains the pairing probabilities. """ return self._pairing_correlation_matrix_pf()RNA secondary-structure model with continuous pairing penalties.
The class wraps one or more ViennaRNA dynamic-programming ensembles and supports continuous per-nucleotide pairing penalties. Multiple ensembles may be used to assign arbitrary energy biases to the paired/unpaired states of selected nucleotides.
A penalty λᵢ is added whenever nucleotide i is paired. Internally, these penalties are automatically represented as an equivalent hybrid combination of unpaired and pair soft constraints. This representation avoids numerical overflows in partition-function calculations while preserving the requested thermodynamic model.
Partition-function calculations use the exact continuous penalties. Minimum- free-energy and suboptimal structure prediction use ViennaRNA's rounded soft constraints to generate candidate structures, which are then rescored using the exact continuous penalties.
Parameters
seq:str- RNA sequence.
lambdas1d:array-like, optional- Per-nucleotide pairing penalties (kcal/mol). Positive values penalize pairing, whereas negative values favor pairing. If omitted, all penalties are zero.
force_paired:array-likeofint, optional- Zero-based indices of nucleotides that are required to be paired, without specifying their pairing partners.
force_unpaired:array-likeofint, optional- Zero-based indices of nucleotides that are required to be unpaired.
state_positions:array-likeofintorsequenceofarray-like, optional- Zero-based indices defining binary paired/unpaired states. A flat
sequence defines one state set. A sequence of sequences defines
multiple disjoint state sets whose energy biases are additive. If
state_biasesis omitted, every state is assigned zero bias. state_biases:array-likeorsequenceofarray-like, optional- Energy biases (kcal/mol) for the states defined by
state_positions. For one state set, its shape must be(2,) * len(state_positions). For multiple sets, provide one tensor per set, with shape(2,) * len(positions). Index zero denotes an unpaired nucleotide and index one a paired nucleotide. It cannot be provided withoutstate_positions. reduce_state_space:bool, default=True- If True, represent the final selected state position as an equivalent 1D pairing penalty, reducing the number of dynamic-programming ensembles from 2N to 2(N-1). If False, use one hard-conditioned ensemble for every state.
no_lonely_pair:boolorNone, default=None- If True, exclude structures containing isolated base pairs using
ViennaRNA's
noLPmodel option. If None, use the current module default. pf_smooth:boolorNone, default=None- Whether ViennaRNA should smooth energies in partition-function Boltzmann factors. If None, use the current module default, which is initially False. If True, partition functions and sampling use ViennaRNA's smoothed model, while MFE, evaluation, and suboptimal structures continue to use its unsmoothed energy model.
temperature:floatorNone, default=None- Temperature in kelvin. If None, use the current module default, which is initially 310.15 K.
NaCl:floatorNone, default=None- Sodium concentration (M). If None, use the current module default, which initially selects ViennaRNA's default value.
parameters:{"turner1999", "turner2004", "andronescu2007", "langdon2018"}orNone- Thermodynamic parameter set. If None, use the current module default.
Notes
By default, a state set containing N positions uses 2(N-1) dynamic-programming ensembles. For multiple sets of sizes N_k, the number is 2(sum(N_k)-K). The final position of each set is represented within each ensemble by an equivalent constant energy shift and 1D pairing penalty. Set
reduce_state_space=Falseto explicitly condition every selected position.Default ViennaRNA builds round soft constraints to the nearest 0.01 kcal/mol in partition-function calculations. When continuous soft constraints are not natively supported, this class automatically applies a lightweight Python callback to recover the exact continuous model.
Methods
def base_pairing_probability(self)-
Expand source code
def base_pairing_probability(self): """ Return the base-pairing probability matrix. Returns ------- ndarray Symmetric NxN matrix whose element (i,j) is the equilibrium probability that nucleotides i and j form a base pair. """ # Compute BPPs before mixture weights: a BPP calculation also produces # the component free energy, while doing this in the opposite order # would require rebuilding each PF without probability backtracking. for molecule in self._dp_molecules: molecule._ensure_pf(compute_bpp=True) if len(self._dp_molecules) == 1: return self._dp_molecules[0]._base_pairing_probability.copy() _, probabilities = self._component_probabilities() matrix = np.zeros_like( self._dp_molecules[0]._base_pairing_probability ) for probability, molecule in zip( probabilities, self._dp_molecules, ): if probability == 0.0: continue matrix += ( probability * molecule._base_pairing_probability ) return matrixReturn the base-pairing probability matrix.
Returns
ndarray- Symmetric NxN matrix whose element (i,j) is the equilibrium probability that nucleotides i and j form a base pair.
def d_free_energy_d_lambdas1d(self)-
Expand source code
def d_free_energy_d_lambdas1d(self): """ Return the derivatives of the free energy with respect to `lambdas1d`. The derivative for nucleotide ``i`` is its equilibrium pairing probability, ``dF / d lambda_i = <s_i>``. Returns ------- ndarray One-dimensional array containing one derivative per nucleotide. """ return np.sum(self.base_pairing_probability(), axis=1)Return the derivatives of the free energy with respect to
lambdas1d.The derivative for nucleotide
iis its equilibrium pairing probability,dF / d lambda_i = <s_i>.Returns
ndarray- One-dimensional array containing one derivative per nucleotide.
def d_free_energy_d_state_biases(self)-
Expand source code
def d_free_energy_d_state_biases(self): """ Return the derivatives of the free energy with respect to state biases. The derivative with respect to a state's energy bias is the equilibrium probability of that state. Returns ------- ndarray or list of ndarray With the flat, single-set `state_positions` syntax, an array with the same shape as `state_biases`. With the list-of-sets syntax, one array per set, including when the list contains only one set. Raises ------ ValueError If this molecule was not initialized with state biases. """ if not self._has_state_biases: raise ValueError( "This molecule was not initialized with state biases" ) needs_pairing_probabilities = any( position is not None for position in self._implicit_state_positions ) if needs_pairing_probabilities: component_pairing_probabilities = [ np.sum(molecule.base_pairing_probability(), axis=1) for molecule in self._dp_molecules ] else: component_pairing_probabilities = [None] * len( self._dp_molecules ) _, component_probabilities = self._component_probabilities() probability_sets = [ np.zeros_like(biases) for biases in self._state_bias_sets ] for ( state, component_probability, pairing_probabilities, ) in zip( self._states, component_probabilities, component_pairing_probabilities, ): if component_probability == 0.0: continue for probabilities, state_slice, implicit_position in zip( probability_sets, self._state_slices, self._implicit_state_positions, ): set_state = state[state_slice] if implicit_position is None: probabilities[set_state] += component_probability continue paired_probability = float(np.clip( pairing_probabilities[implicit_position], 0.0, 1.0, )) probabilities[set_state + (0,)] += ( component_probability * (1.0 - paired_probability) ) probabilities[set_state + (1,)] += ( component_probability * paired_probability ) if not self._state_sets_were_nested: return probability_sets[0] return probability_setsReturn the derivatives of the free energy with respect to state biases.
The derivative with respect to a state's energy bias is the equilibrium probability of that state.
Returns
ndarrayorlistofndarray- With the flat, single-set
state_positionssyntax, an array with the same shape asstate_biases. With the list-of-sets syntax, one array per set, including when the list contains only one set.
Raises
ValueError- If this molecule was not initialized with state biases.
def evaluate(self, structure)-
Expand source code
def evaluate(self, structure): """ Return the exact free energy of a secondary structure. Parameters ---------- structure : str Dot-bracket representation of a secondary structure. Returns ------- float Exact free energy (kcal/mol), including continuous pairing penalties and state biases. """ if not isinstance(structure, str): raise ValueError("structure must be a string") if len(structure) != len(self._seq): raise ValueError( "structure must contain one symbol per nucleotide" ) if structure.strip(".()"): raise ValueError( "structure must use '.', '(', and ')' dot-bracket symbols" ) # Components follow np.ndindex order, so the paired/unpaired state is # also the component index interpreted as a binary integer. Avoid # constructing a state tuple on this latency-sensitive path. component_index = 0 for position in self._explicit_state_positions: component_index = ( 2 * component_index + (structure[position] != ".") ) energy = self._dp_molecules[component_index].evaluate(structure) return float( energy + self._component_biases.flat[component_index] )Return the exact free energy of a secondary structure.
Parameters
structure:str- Dot-bracket representation of a secondary structure.
Returns
float- Exact free energy (kcal/mol), including continuous pairing penalties and state biases.
def mfe(self)-
Expand source code
def mfe(self): """ Return the minimum-free-energy structure. The returned energy always corresponds to the exact continuous pairing penalties, even when ViennaRNA internally rounds soft constraints. Returns ------- structure : str Dot-bracket representation of the MFE structure. energy : float Exact free energy (kcal/mol). """ # A single DP component needs no mixture-level search. if len(self._dp_molecules) == 1: structure, energy = self._dp_molecules[0].mfe() return ( structure, float(energy + self._component_biases.item()), ) structures, energies = self._component_mfes() index = int(np.argmin(energies)) return structures[index], float(energies[index])Return the minimum-free-energy structure.
The returned energy always corresponds to the exact continuous pairing penalties, even when ViennaRNA internally rounds soft constraints.
Returns
structure:str- Dot-bracket representation of the MFE structure.
energy:float- Exact free energy (kcal/mol).
def pairing_correlation_matrix(self)-
Expand source code
def pairing_correlation_matrix(self): """ Return the joint pairing-probability matrix. Element ``(i, j)`` is the probability that nucleotides ``i`` and ``j`` are simultaneously paired, irrespective of their pairing partners. The matrix is computed exactly using constrained partition-function calculations. Returns ------- ndarray Symmetric NxN matrix whose diagonal contains the pairing probabilities. """ return self._pairing_correlation_matrix_pf()Return the joint pairing-probability matrix.
Element
(i, j)is the probability that nucleotidesiandjare simultaneously paired, irrespective of their pairing partners.The matrix is computed exactly using constrained partition-function calculations.
Returns
ndarray- Symmetric NxN matrix whose diagonal contains the pairing probabilities.
def sample(self, number, weights=False)-
Expand source code
def sample(self, number, weights=False): """ Generate Boltzmann-distributed secondary structures. By default, rejection sampling corrects ViennaRNA's rounded soft constraints and returns unweighted structures from the exact continuous model. Weighted samples from the rounded model can be requested instead. Parameters ---------- number : int Number of structures to sample. weights : bool, default=False If False, return unweighted structures sampled from the exact model. If True, return rounded-model samples and their log-weight corrections. The weighted path is faster when rounding corrections are present. Returns ------- list of str or list of (str, float) With ``weights=False``, a list of dot-bracket structures. With ``weights=True``, each element contains a structure and its unnormalized log-weight correction log(w) = -(E_exact - E_rounded) / (k_B T). For state mixtures, the log-weight also includes the component-specific exact/rounded normalization correction. This correction is caused only by lambda rounding, not by the state bias. When all lambdas are multiples of 0.01 kcal/mol, every returned log-weight is zero. """ number = int(number) if number <= 0: raise ValueError("number must be a positive integer") if not isinstance(weights, (bool, np.bool_)): raise ValueError("weights must be a boolean") # A single DP component needs no mixture-level sampling. if len(self._dp_molecules) == 1: return self._dp_molecules[0].sample( number, weights=weights, ) _, probabilities = self._component_probabilities() component_indices = np.random.choice( len(self._dp_molecules), size=number, p=probabilities, ) result = [None] * number for component_index, molecule in enumerate(self._dp_molecules): output_indices = np.flatnonzero( component_indices == component_index ) if not len(output_indices): continue samples = molecule.sample( len(output_indices), weights=weights, ) if weights: rounding_correction = ( molecule.sample_rounding_correction() ) for output_index, (structure, log_weight) in zip( output_indices, samples, ): result[output_index] = ( structure, float(log_weight + rounding_correction), ) else: for output_index, structure in zip( output_indices, samples, ): result[output_index] = structure return resultGenerate Boltzmann-distributed secondary structures.
By default, rejection sampling corrects ViennaRNA's rounded soft constraints and returns unweighted structures from the exact continuous model. Weighted samples from the rounded model can be requested instead.
Parameters
number:int- Number of structures to sample.
weights:bool, default=False- If False, return unweighted structures sampled from the exact model. If True, return rounded-model samples and their log-weight corrections. The weighted path is faster when rounding corrections are present.
Returns
listofstrorlistof(str, float)-
With
weights=False, a list of dot-bracket structures. Withweights=True, each element contains a structure and its unnormalized log-weight correctionlog(w) = -(E_exact - E_rounded) / (k_B T).For state mixtures, the log-weight also includes the component-specific exact/rounded normalization correction. This correction is caused only by lambda rounding, not by the state bias.
When all lambdas are multiples of 0.01 kcal/mol, every returned log-weight is zero.
def suboptimal_coverage(self, delta)-
Expand source code
def suboptimal_coverage(self, delta): """ Return the fraction of the partition function represented by the suboptimal ensemble. The suboptimal ensemble contains all structures whose exact energy is within ``delta`` kcal/mol of the minimum-free-energy structure. Parameters ---------- delta : float Maximum energy difference (kcal/mol) above the exact MFE. Returns ------- float Fraction of the total partition function represented by the enumerated structures. The result lies between zero and one, apart from possible small numerical errors. """ suboptimal = self.suboptimal_structures(delta) reference_energy = suboptimal[0][1] inverse_kT = 1.0 / (_KB * self._temperature) relative_partition_function = math.fsum( math.exp(-(energy - reference_energy) * inverse_kT) for _, energy in suboptimal ) log_coverage = ( math.log(relative_partition_function) + (self.total_free_energy() - reference_energy) * inverse_kT ) coverage = float(math.exp(log_coverage)) if coverage > 1.0: if coverage <= 1.0 + 1e-5: coverage = 1.0 else: raise RuntimeError( f"Suboptimal coverage is unexpectedly larger than one: " f"{coverage}" ) return coverageReturn the fraction of the partition function represented by the suboptimal ensemble.
The suboptimal ensemble contains all structures whose exact energy is within
deltakcal/mol of the minimum-free-energy structure.Parameters
delta:float- Maximum energy difference (kcal/mol) above the exact MFE.
Returns
float- Fraction of the total partition function represented by the enumerated structures. The result lies between zero and one, apart from possible small numerical errors.
def suboptimal_structures(self, delta)-
Expand source code
def suboptimal_structures(self, delta): """ Enumerate suboptimal secondary structures. Candidate structures are generated using ViennaRNA's rounded soft constraints, rescored with the exact continuous penalties, and returned sorted by exact energy. For state mixtures, the MFE of every feasible state is calculated first. Each component is then enumerated with the local energy window required by the common, globally biased cutoff. Parameters ---------- delta : float Maximum energy difference (kcal/mol) above the exact MFE. Returns ------- list of (str, float) List of (structure, energy) pairs sorted by increasing exact energy. """ # A single DP component needs no mixture-level enumeration. if len(self._dp_molecules) == 1: suboptimal = self._dp_molecules[0].suboptimal_structures( delta ) bias = float(self._component_biases.item()) if bias == 0.0: return suboptimal return [ (structure, float(energy + bias)) for structure, energy in suboptimal ] delta = float(delta) if not np.isfinite(delta) or delta < 0.0: raise ValueError("delta must be finite and non-negative") _, component_mfe_energies = self._component_mfes() global_mfe_energy = float(np.min(component_mfe_energies)) cutoff = global_mfe_energy + delta merged = [] for molecule, bias, component_mfe_energy in zip( self._dp_molecules, self._component_biases.ravel(), component_mfe_energies, ): if not np.isfinite(component_mfe_energy): continue local_delta = cutoff - component_mfe_energy if local_delta < 0.0: continue for structure, energy in molecule.suboptimal_structures( local_delta ): biased_energy = float(energy + bias) if biased_energy <= cutoff: merged.append((structure, biased_energy)) merged.sort(key=lambda item: item[1]) return mergedEnumerate suboptimal secondary structures.
Candidate structures are generated using ViennaRNA's rounded soft constraints, rescored with the exact continuous penalties, and returned sorted by exact energy.
For state mixtures, the MFE of every feasible state is calculated first. Each component is then enumerated with the local energy window required by the common, globally biased cutoff.
Parameters
delta:float- Maximum energy difference (kcal/mol) above the exact MFE.
Returns
listof(str, float)- List of (structure, energy) pairs sorted by increasing exact energy.
def total_free_energy(self)-
Expand source code
def total_free_energy(self): """ Return the ensemble free energy. Returns ------- float Ensemble free energy (kcal/mol) corresponding to the exact continuous pairing penalties. """ total_free_energy, _ = self._component_probabilities() return total_free_energyReturn the ensemble free energy.
Returns
float- Ensemble free energy (kcal/mol) corresponding to the exact continuous pairing penalties.