ViennaRNA-bindings-0.233.1.1: C/ViennaRNA/interior_loops.h
#ifndef VIENNA_RNA_PACKAGE_INTERIOR_LOOPS_H
#define VIENNA_RNA_PACKAGE_INTERIOR_LOOPS_H
#include <ViennaRNA/utils.h>
#include "ViennaRNA/energy_par.h"
#include <ViennaRNA/data_structures.h>
#include <ViennaRNA/params.h>
#include <ViennaRNA/constraints.h>
#ifdef __GNUC__
# define INLINE inline
#else
# define INLINE
#endif
#ifdef ON_SAME_STRAND
#undef ON_SAME_STRAND
#endif
#define ON_SAME_STRAND(I,J,C) (((I)>=(C))||((J)<(C)))
/**
* @file interior_loops.h
* @ingroup loops
* @brief Energy evaluation of interior loops for MFE and partition function calculations
*/
/**
* @{
* @ingroup loops
*/
/**
* <H2>Compute the Energy of an interior-loop</H2>
* This function computes the free energy @f$\Delta G@f$ of an interior-loop with the
* following structure: <BR>
* <PRE>
* 3' 5'
* | |
* U - V
* a_n b_1
* . .
* . .
* . .
* a_1 b_m
* X - Y
* | |
* 5' 3'
* </PRE>
* This general structure depicts an interior-loop that is closed by the base pair (X,Y).
* The enclosed base pair is (V,U) which leaves the unpaired bases a_1-a_n and b_1-b_n
* that constitute the loop. In this example, the length of the interior-loop is @f$(n+m)@f$
* where n or m may be 0 resulting in a bulge-loop or base pair stack.
* The mismatching nucleotides for the closing pair (X,Y) are:<BR>
* 5'-mismatch: a_1<BR>
* 3'-mismatch: b_m<BR>
* and for the enclosed base pair (V,U):<BR>
* 5'-mismatch: b_1<BR>
* 3'-mismatch: a_n<BR>
* @note Base pairs are always denoted in 5'->3' direction. Thus the enclosed base pair
* must be 'turned arround' when evaluating the free energy of the interior-loop
* @see scale_parameters()
* @see vrna_param_t
* @note This function is threadsafe
*
* @param n1 The size of the 'left'-loop (number of unpaired nucleotides)
* @param n2 The size of the 'right'-loop (number of unpaired nucleotides)
* @param type The pair type of the base pair closing the interior loop
* @param type_2 The pair type of the enclosed base pair
* @param si1 The 5'-mismatching nucleotide of the closing pair
* @param sj1 The 3'-mismatching nucleotide of the closing pair
* @param sp1 The 3'-mismatching nucleotide of the enclosed pair
* @param sq1 The 5'-mismatching nucleotide of the enclosed pair
* @param P The datastructure containing scaled energy parameters
* @return The Free energy of the Interior-loop in dcal/mol
*/
PRIVATE INLINE int E_IntLoop(int n1,
int n2,
int type,
int type_2,
int si1,
int sj1,
int sp1,
int sq1,
vrna_param_t *P);
/**
* <H2>Compute Boltzmann weight @f$e^{-\Delta G/kT} @f$ of interior loop</H2>
* multiply by scale[u1+u2+2] for scaling
* @see get_scaled_pf_parameters()
* @see vrna_exp_param_t
* @see E_IntLoop()
* @note This function is threadsafe
*
* @param u1 The size of the 'left'-loop (number of unpaired nucleotides)
* @param u2 The size of the 'right'-loop (number of unpaired nucleotides)
* @param type The pair type of the base pair closing the interior loop
* @param type2 The pair type of the enclosed base pair
* @param si1 The 5'-mismatching nucleotide of the closing pair
* @param sj1 The 3'-mismatching nucleotide of the closing pair
* @param sp1 The 3'-mismatching nucleotide of the enclosed pair
* @param sq1 The 5'-mismatching nucleotide of the enclosed pair
* @param P The datastructure containing scaled Boltzmann weights of the energy parameters
* @return The Boltzmann weight of the Interior-loop
*/
PRIVATE INLINE FLT_OR_DBL exp_E_IntLoop(int u1,
int u2,
int type,
int type2,
short si1,
short sj1,
short sp1,
short sq1,
vrna_exp_param_t *P);
PRIVATE INLINE int E_IntLoop_Co(int type,
int type_2,
int i,
int j,
int p,
int q,
int cutpoint,
short si1,
short sj1,
short sp1,
short sq1,
int dangles,
vrna_param_t *P);
/*
#################################
# BEGIN OF FUNCTION DEFINITIONS #
#################################
*/
/*
* ugly but fast interior loop evaluation
*
* Avoid including this function in your own code. It only serves
* as a fast inline block internally re-used throughout the RNAlib. It
* evalutes the free energy of interior loops in single sequences or sequence
* hybrids. Soft constraints are also applied if available.
*
* NOTE: do not include into doxygen reference manual!
*/
PRIVATE INLINE int
ubf_eval_int_loop( int i,
int j,
int p,
int q,
int i1,
int j1,
int p1,
int q1,
short si,
short sj,
short sp,
short sq,
unsigned char type,
unsigned char type_2,
int *rtype,
int ij,
int cp,
vrna_param_t *P,
vrna_sc_t *sc){
int energy, u1, u2;
u1 = p1 - i;
u2 = j1 - q;
if((cp < 0) || (ON_SAME_STRAND(i, p, cp) && ON_SAME_STRAND(q, j, cp))){ /* regular interior loop */
energy = E_IntLoop(u1, u2, type, type_2, si, sj, sp, sq, P);
} else { /* interior loop like cofold structure */
short Si, Sj;
Si = ON_SAME_STRAND(i, i1, cp) ? si : -1;
Sj = ON_SAME_STRAND(j1, j, cp) ? sj : -1;
energy = E_IntLoop_Co(rtype[type], rtype[type_2],
i, j, p, q,
cp,
Si, Sj,
sp, sq,
P->model_details.dangles,
P);
}
/* add soft constraints */
if(sc){
if(sc->energy_up)
energy += sc->energy_up[i1][u1]
+ sc->energy_up[q1][u2];
if(sc->energy_bp)
energy += sc->energy_bp[ij];
if(sc->energy_stack)
if(u1 + u2 == 0){
int a = sc->energy_stack[i]
+ sc->energy_stack[p]
+ sc->energy_stack[q]
+ sc->energy_stack[j];
energy += a;
}
if(sc->f)
energy += sc->f(i, j, p, q, VRNA_DECOMP_PAIR_IL, sc->data);
}
return energy;
}
/*
* ugly but fast exterior interior loop evaluation
*
* Avoid including this function in your own code. It only serves
* as a fast inline block internally re-used throughout the RNAlib. It
* evalutes the free energy of interior loops in single sequences or sequence
* hybrids. Soft constraints are also applied if available.
*
* NOTE: do not include into doxygen reference manual!
*/
PRIVATE INLINE int
ubf_eval_ext_int_loop(int i,
int j,
int p,
int q,
int i1,
int j1,
int p1,
int q1,
short si,
short sj,
short sp,
short sq,
unsigned char type,
unsigned char type_2,
int length,
vrna_param_t *P,
vrna_sc_t *sc){
int energy, u1, u2, u3;
u1 = i1;
u2 = p1 - j;
u3 = length - q;
energy = E_IntLoop(u2, u1 + u3, type, type_2, si, sj, sp, sq, P);
/* add soft constraints */
if(sc){
if(sc->energy_up){
energy += sc->energy_up[j1][u2]
+ ((u3 > 0) ? sc->energy_up[q1][u3] : 0)
+ ((u1 > 0) ? sc->energy_up[1][u1] : 0);
}
if(sc->energy_stack)
if(u1 + u2 + u3 == 0)
energy += sc->energy_stack[i]
+ sc->energy_stack[p]
+ sc->energy_stack[q]
+ sc->energy_stack[j];
if(sc->f)
energy += sc->f(i, j, p, q, VRNA_DECOMP_PAIR_IL, sc->data);
}
return energy;
}
PRIVATE INLINE int
E_IntLoop(int n1,
int n2,
int type,
int type_2,
int si1,
int sj1,
int sp1,
int sq1,
vrna_param_t *P){
/* compute energy of degree 2 loop (stack bulge or interior) */
int nl, ns, u, energy;
energy = INF;
if (n1>n2) { nl=n1; ns=n2;}
else {nl=n2; ns=n1;}
if (nl == 0)
return P->stack[type][type_2]; /* stack */
if (ns==0) { /* bulge */
energy = (nl<=MAXLOOP)?P->bulge[nl]:
(P->bulge[30]+(int)(P->lxc*log(nl/30.)));
if (nl==1) energy += P->stack[type][type_2];
else {
if (type>2) energy += P->TerminalAU;
if (type_2>2) energy += P->TerminalAU;
}
return energy;
}
else { /* interior loop */
if (ns==1) {
if (nl==1) /* 1x1 loop */
return P->int11[type][type_2][si1][sj1];
if (nl==2) { /* 2x1 loop */
if (n1==1)
energy = P->int21[type][type_2][si1][sq1][sj1];
else
energy = P->int21[type_2][type][sq1][si1][sp1];
return energy;
}
else { /* 1xn loop */
energy = (nl+1<=MAXLOOP)?(P->internal_loop[nl+1]) : (P->internal_loop[30]+(int)(P->lxc*log((nl+1)/30.)));
energy += MIN2(MAX_NINIO, (nl-ns)*P->ninio[2]);
energy += P->mismatch1nI[type][si1][sj1] + P->mismatch1nI[type_2][sq1][sp1];
return energy;
}
}
else if (ns==2) {
if(nl==2) { /* 2x2 loop */
return P->int22[type][type_2][si1][sp1][sq1][sj1];}
else if (nl==3){ /* 2x3 loop */
energy = P->internal_loop[5]+P->ninio[2];
energy += P->mismatch23I[type][si1][sj1] + P->mismatch23I[type_2][sq1][sp1];
return energy;
}
}
{ /* generic interior loop (no else here!)*/
u = nl + ns;
energy = (u <= MAXLOOP) ? (P->internal_loop[u]) : (P->internal_loop[30]+(int)(P->lxc*log((u)/30.)));
energy += MIN2(MAX_NINIO, (nl-ns)*P->ninio[2]);
energy += P->mismatchI[type][si1][sj1] + P->mismatchI[type_2][sq1][sp1];
}
}
return energy;
}
PRIVATE INLINE FLT_OR_DBL
exp_E_IntLoop(int u1,
int u2,
int type,
int type2,
short si1,
short sj1,
short sp1,
short sq1,
vrna_exp_param_t *P){
int ul, us, no_close = 0;
double z = 0.;
int noGUclosure = P->model_details.noGUclosure;
if ((noGUclosure) && ((type2==3)||(type2==4)||(type==3)||(type==4)))
no_close = 1;
if (u1>u2) { ul=u1; us=u2;}
else {ul=u2; us=u1;}
if (ul==0) /* stack */
z = P->expstack[type][type2];
else if(!no_close){
if (us==0) { /* bulge */
z = P->expbulge[ul];
if (ul==1) z *= P->expstack[type][type2];
else {
if (type>2) z *= P->expTermAU;
if (type2>2) z *= P->expTermAU;
}
return (FLT_OR_DBL)z;
}
else if (us==1) {
if (ul==1){ /* 1x1 loop */
return (FLT_OR_DBL)(P->expint11[type][type2][si1][sj1]);
}
if (ul==2) { /* 2x1 loop */
if (u1==1)
return (FLT_OR_DBL)(P->expint21[type][type2][si1][sq1][sj1]);
else
return (FLT_OR_DBL)(P->expint21[type2][type][sq1][si1][sp1]);
}
else { /* 1xn loop */
z = P->expinternal[ul+us] * P->expmismatch1nI[type][si1][sj1] * P->expmismatch1nI[type2][sq1][sp1];
return (FLT_OR_DBL)(z * P->expninio[2][ul-us]);
}
}
else if (us==2) {
if(ul==2) /* 2x2 loop */
return (FLT_OR_DBL)(P->expint22[type][type2][si1][sp1][sq1][sj1]);
else if(ul==3){ /* 2x3 loop */
z = P->expinternal[5]*P->expmismatch23I[type][si1][sj1]*P->expmismatch23I[type2][sq1][sp1];
return (FLT_OR_DBL)(z * P->expninio[2][1]);
}
}
/* generic interior loop (no else here!)*/
z = P->expinternal[ul+us] * P->expmismatchI[type][si1][sj1] * P->expmismatchI[type2][sq1][sp1];
return (FLT_OR_DBL)(z * P->expninio[2][ul-us]);
}
return (FLT_OR_DBL)z;
}
PRIVATE INLINE int
E_IntLoop_Co( int type,
int type_2,
int i,
int j,
int p,
int q,
int cutpoint,
short si1,
short sj1,
short sp1,
short sq1,
int dangles,
vrna_param_t *P){
int energy, ci, cj, cp, cq, d3, d5, d5_2, d3_2, tmm, tmm_2;
energy = 0;
if(type > 2) energy += P->TerminalAU;
if(type_2 > 2) energy += P->TerminalAU;
if(!dangles) return energy;
ci = ON_SAME_STRAND(i, i + 1, cutpoint);
cj = ON_SAME_STRAND(j - 1, j, cutpoint);
cp = ON_SAME_STRAND(p - 1, p, cutpoint);
cq = ON_SAME_STRAND(q, q + 1, cutpoint);
d3 = ci ? P->dangle3[type][si1] : 0;
d5 = cj ? P->dangle5[type][sj1] : 0;
d5_2 = cp ? P->dangle5[type_2][sp1] : 0;
d3_2 = cq ? P->dangle3[type_2][sq1] : 0;
tmm = (cj && ci) ? P->mismatchExt[type][sj1][si1] : d5 + d3;
tmm_2 = (cp && cq) ? P->mismatchExt[type_2][sp1][sq1] : d5_2 + d3_2;
if(dangles == 2) return energy + tmm + tmm_2;
/* now we may have non-double dangles only */
if(i+2 < p){
if(q+2 < j){ energy += tmm + tmm_2;}
else if(q+2 == j){ energy += (cj && cq) ? MIN2(tmm + d5_2, tmm_2 + d3) : tmm + tmm_2;}
else energy += d3 + d5_2;
}
else if(i+2 == p){
if(q+2 < j){ energy += (ci && cp) ? MIN2(tmm + d3_2, tmm_2 + d5) : tmm + tmm_2;}
else if(q+2 == j){
energy += MIN2(tmm, MIN2(tmm_2, MIN2(d5 + d5_2, d3 + d3_2)));
}
else energy += MIN2(d3, d5_2);
}
else{
if(q+2 < j){ energy += d5 + d3_2;}
else if(q+2 == j){ energy += MIN2(d5, d3_2);}
}
return energy;
}
int
vrna_E_int_loop(vrna_fold_compound_t *vc,
int i,
int j);
int
vrna_eval_int_loop( vrna_fold_compound_t *vc,
int i,
int j,
int k,
int l);
FLT_OR_DBL
vrna_exp_E_int_loop(vrna_fold_compound_t *vc,
int i,
int j);
FLT_OR_DBL
vrna_exp_E_interior_loop( vrna_fold_compound_t *vc,
int i,
int j,
int k,
int l);
int
vrna_E_ext_int_loop(vrna_fold_compound_t *vc,
int i,
int j,
int *ip,
int *iq);
int
vrna_E_stack( vrna_fold_compound_t *vc,
int i,
int j);
/**
* @brief Backtrack a stacked pair closed by @f$ (i,j) @f$
*
*/
int
vrna_BT_stack(vrna_fold_compound_t *vc,
int *i,
int *j,
int *en,
vrna_bp_stack_t *bp_stack,
int *stack_count);
/**
* @brief Backtrack an interior loop closed by @f$ (i,j) @f$
*
*/
int
vrna_BT_int_loop( vrna_fold_compound_t *vc,
int *i,
int *j,
int en,
vrna_bp_stack_t *bp_stack,
int *stack_count);
/**
* @}
*/
#endif