#ifndef __VIENNA_RNA_PACKAGE_GQUAD_H__
#define __VIENNA_RNA_PACKAGE_GQUAD_H__
#include "data_structures.h"
#ifndef INLINE
#ifdef __GNUC__
# define INLINE inline
#else
# define INLINE
#endif
#endif
/**
* \file gquad.h
* \brief Various functions related to G-quadruplex computations
*/
int E_gquad(int L,
int l[3],
paramT *P);
FLT_OR_DBL exp_E_gquad( int L,
int l[3],
pf_paramT *pf);
int E_gquad_ali(int i,
int L,
int l[3],
const short **S,
int n_seq,
paramT *P);
void E_gquad_ali_en( int i,
int L,
int l[3],
const short **S,
int n_seq,
int en[2],
paramT *P);
/**
* \brief Get a triangular matrix prefilled with minimum free energy
* contributions of G-quadruplexes.
*
* At each position ij in the matrix, the minimum free energy of any
* G-quadruplex delimited by i and j is stored. If no G-quadruplex formation
* is possible, the matrix element is set to INF.
* Access the elements in the matrix via matrix[indx[j]+i]. To get
* the integer array indx see get_jindx().
*
* \see get_jindx(), encode_sequence()
*
* \param S The encoded sequence
* \param P A pointer to the data structure containing the precomputed energy contributions
* \return A pointer to the G-quadruplex contribution matrix
*/
int *get_gquad_matrix(short *S, paramT *P);
int *get_gquad_ali_matrix(short *S_cons,
short **S,
int n_seq,
paramT *P);
FLT_OR_DBL *get_gquad_pf_matrix( short *S,
FLT_OR_DBL *scale,
pf_paramT *pf);
int **get_gquad_L_matrix( short *S,
int start,
int maxdist,
int **g,
paramT *P);
void get_gquad_pattern_mfe(short *S,
int i,
int j,
paramT *P,
int *L,
int l[3]);
void
get_gquad_pattern_exhaustive( short *S,
int i,
int j,
paramT *P,
int *L,
int *l,
int threshold);
void get_gquad_pattern_pf( short *S,
int i,
int j,
pf_paramT *pf,
int *L,
int l[3]);
plist *get_plist_gquad_from_pr( short *S,
int gi,
int gj,
FLT_OR_DBL *G,
FLT_OR_DBL *probs,
FLT_OR_DBL *scale,
pf_paramT *pf);
plist *get_plist_gquad_from_pr_max(short *S,
int gi,
int gj,
FLT_OR_DBL *G,
FLT_OR_DBL *probs,
FLT_OR_DBL *scale,
int *L,
int l[3],
pf_paramT *pf);
plist *get_plist_gquad_from_db( const char *structure,
float pr);
int get_gquad_count(short *S,
int i,
int j);
int get_gquad_layer_count(short *S,
int i,
int j);
/**
* given a dot-bracket structure (possibly) containing gquads encoded
* by '+' signs, find first gquad, return end position or 0 if none found
* Upon return L and l[] contain the number of stacked layers, as well as
* the lengths of the linker regions.
* To parse a string with many gquads, call parse_gquad repeatedly e.g.
* end1 = parse_gquad(struc, &L, l); ... ;
* end2 = parse_gquad(struc+end1, &L, l); end2+=end1; ... ;
* end3 = parse_gquad(struc+end2, &L, l); end3+=end2; ... ;
*/
int parse_gquad(const char *struc, int *L, int l[3]);
/**
* backtrack an interior loop like enclosed g-quadruplex
* with closing pair (i,j)
*
* \param c The total contribution the loop should resemble
* \param i position i of enclosing pair
* \param j position j of enclosing pair
* \param type base pair type of enclosing pair (must be reverse type)
* \param S integer encoded sequence
* \param ggg triangular matrix containing g-quadruplex contributions
* \param index the index for accessing the triangular matrix
* \param p here the 5' position of the gquad is stored
* \param q here the 3' position of the gquad is stored
* \param P the datastructure containing the precalculated contibutions
*
* \return 1 on success, 0 if no gquad found
*/
INLINE PRIVATE int backtrack_GQuad_IntLoop(int c,
int i,
int j,
int type,
short *S,
int *ggg,
int *index,
int *p,
int *q,
paramT *P){
int energy, dangles, k, l, maxl, minl, c0, l1;
short si, sj;
dangles = P->model_details.dangles;
si = S[i + 1];
sj = S[j - 1];
energy = 0;
if(dangles == 2)
energy += P->mismatchI[type][si][sj];
if(type > 2)
energy += P->TerminalAU;
k = i + 1;
if(S[k] == 3){
if(k < j - VRNA_GQUAD_MIN_BOX_SIZE){
minl = j - i + k - MAXLOOP - 2;
c0 = k + VRNA_GQUAD_MIN_BOX_SIZE - 1;
minl = MAX2(c0, minl);
c0 = j - 3;
maxl = k + VRNA_GQUAD_MAX_BOX_SIZE + 1;
maxl = MIN2(c0, maxl);
for(l = minl; l < maxl; l++){
if(S[l] != 3) continue;
if(c == energy + ggg[index[l] + k] + P->internal_loop[j - l - 1]){
*p = k; *q = l;
return 1;
}
}
}
}
for(k = i + 2;
k < j - VRNA_GQUAD_MIN_BOX_SIZE;
k++){
l1 = k - i - 1;
if(l1>MAXLOOP) break;
if(S[k] != 3) continue;
minl = j - i + k - MAXLOOP - 2;
c0 = k + VRNA_GQUAD_MIN_BOX_SIZE - 1;
minl = MAX2(c0, minl);
c0 = j - 1;
maxl = k + VRNA_GQUAD_MAX_BOX_SIZE + 1;
maxl = MIN2(c0, maxl);
for(l = minl; l < maxl; l++){
if(S[l] != 3) continue;
if(c == energy + ggg[index[l] + k] + P->internal_loop[l1 + j - l - 1]){
*p = k; *q = l;
return 1;
}
}
}
l = j - 1;
if(S[l] == 3)
for(k = i + 4;
k < j - VRNA_GQUAD_MIN_BOX_SIZE;
k++){
l1 = k - i - 1;
if(l1>MAXLOOP) break;
if(S[k] != 3) continue;
if(c == energy + ggg[index[l] + k] + P->internal_loop[l1]){
*p = k; *q = l;
return 1;
}
}
return 0;
}
/**
* backtrack an interior loop like enclosed g-quadruplex
* with closing pair (i,j) with underlying Lfold matrix
*
* \param c The total contribution the loop should resemble
* \param i position i of enclosing pair
* \param j position j of enclosing pair
* \param type base pair type of enclosing pair (must be reverse type)
* \param S integer encoded sequence
* \param ggg triangular matrix containing g-quadruplex contributions
* \param p here the 5' position of the gquad is stored
* \param q here the 3' position of the gquad is stored
* \param P the datastructure containing the precalculated contibutions
*
* \return 1 on success, 0 if no gquad found
*/
INLINE PRIVATE int backtrack_GQuad_IntLoop_L(int c,
int i,
int j,
int type,
short *S,
int **ggg,
int maxdist,
int *p,
int *q,
paramT *P){
int energy, dangles, k, l, maxl, minl, c0, l1;
short si, sj;
dangles = P->model_details.dangles;
si = S[i + 1];
sj = S[j - 1];
energy = 0;
if(dangles == 2)
energy += P->mismatchI[type][si][sj];
if(type > 2)
energy += P->TerminalAU;
k = i + 1;
if(S[k] == 3){
if(k < j - VRNA_GQUAD_MIN_BOX_SIZE){
minl = j - i + k - MAXLOOP - 2;
c0 = k + VRNA_GQUAD_MIN_BOX_SIZE - 1;
minl = MAX2(c0, minl);
c0 = j - 3;
maxl = k + VRNA_GQUAD_MAX_BOX_SIZE + 1;
maxl = MIN2(c0, maxl);
for(l = minl; l < maxl; l++){
if(S[l] != 3) continue;
if(c == energy + ggg[k][l - k] + P->internal_loop[j - l - 1]){
*p = k; *q = l;
return 1;
}
}
}
}
for(k = i + 2;
k < j - VRNA_GQUAD_MIN_BOX_SIZE;
k++){
l1 = k - i - 1;
if(l1>MAXLOOP) break;
if(S[k] != 3) continue;
minl = j - i + k - MAXLOOP - 2;
c0 = k + VRNA_GQUAD_MIN_BOX_SIZE - 1;
minl = MAX2(c0, minl);
c0 = j - 1;
maxl = k + VRNA_GQUAD_MAX_BOX_SIZE + 1;
maxl = MIN2(c0, maxl);
for(l = minl; l < maxl; l++){
if(S[l] != 3) continue;
if(c == energy + ggg[k][l - k] + P->internal_loop[l1 + j - l - 1]){
*p = k; *q = l;
return 1;
}
}
}
l = j - 1;
if(S[l] == 3)
for(k = i + 4;
k < j - VRNA_GQUAD_MIN_BOX_SIZE;
k++){
l1 = k - i - 1;
if(l1>MAXLOOP) break;
if(S[k] != 3) continue;
if(c == energy + ggg[k][l - k] + P->internal_loop[l1]){
*p = k; *q = l;
return 1;
}
}
return 0;
}
INLINE PRIVATE
int
E_GQuad_IntLoop(int i,
int j,
int type,
short *S,
int *ggg,
int *index,
paramT *P){
int energy, ge, en1, en2, dangles, p, q, l1, minq, maxq;
int c0, c1, c2, c3, up, d53, d5, d3;
short si, sj;
dangles = P->model_details.dangles;
si = S[i + 1];
sj = S[j - 1];
energy = 0;
if(dangles == 2)
energy += P->mismatchI[type][si][sj];
if(type > 2)
energy += P->TerminalAU;
ge = INF;
p = i + 1;
if(S[p] == 3){
if(p < j - VRNA_GQUAD_MIN_BOX_SIZE){
minq = j - i + p - MAXLOOP - 2;
c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
minq = MAX2(c0, minq);
c0 = j - 3;
maxq = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
maxq = MIN2(c0, maxq);
for(q = minq; q < maxq; q++){
if(S[q] != 3) continue;
c0 = energy + ggg[index[q] + p] + P->internal_loop[j - q - 1];
ge = MIN2(ge, c0);
}
}
}
for(p = i + 2;
p < j - VRNA_GQUAD_MIN_BOX_SIZE;
p++){
l1 = p - i - 1;
if(l1>MAXLOOP) break;
if(S[p] != 3) continue;
minq = j - i + p - MAXLOOP - 2;
c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
minq = MAX2(c0, minq);
c0 = j - 1;
maxq = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
maxq = MIN2(c0, maxq);
for(q = minq; q < maxq; q++){
if(S[q] != 3) continue;
c0 = energy + ggg[index[q] + p] + P->internal_loop[l1 + j - q - 1];
ge = MIN2(ge, c0);
}
}
q = j - 1;
if(S[q] == 3)
for(p = i + 4;
p < j - VRNA_GQUAD_MIN_BOX_SIZE;
p++){
l1 = p - i - 1;
if(l1>MAXLOOP) break;
if(S[p] != 3) continue;
c0 = energy + ggg[index[q] + p] + P->internal_loop[l1];
ge = MIN2(ge, c0);
}
#if 0
/* here comes the additional stuff for the odd dangle models */
if(dangles % 1){
en1 = energy + P->dangle5[type][si];
en2 = energy + P->dangle5[type][sj];
en3 = energy + P->mismatchI[type][si][sj];
/* first case with 5' dangle (i.e. j-1) onto enclosing pair */
p = i + 1;
if(S[p] == 3){
if(p < j - VRNA_GQUAD_MIN_BOX_SIZE){
minq = j - i + p - MAXLOOP - 2;
c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
minq = MAX2(c0, minq);
c0 = j - 4;
maxq = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
maxq = MIN2(c0, maxq);
for(q = minq; q < maxq; q++){
if(S[q] != 3) continue;
c0 = en1 + ggg[index[q] + p] + P->internal_loop[j - q - 1];
ge = MIN2(ge, c0);
}
}
}
for(p = i + 2; p < j - VRNA_GQUAD_MIN_BOX_SIZE; p++){
l1 = p - i - 1;
if(l1>MAXLOOP) break;
if(S[p] != 3) continue;
minq = j - i + p - MAXLOOP - 2;
c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
minq = MAX2(c0, minq);
c0 = j - 2;
maxq = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
maxq = MIN2(c0, maxq);
for(q = minq; q < maxq; q++){
if(S[q] != 3) continue;
c0 = en1 + ggg[index[q] + p] + P->internal_loop[l1 + j - q - 1];
ge = MIN2(ge, c0);
}
}
q = j - 2;
if(S[q] == 3)
for(p = i + 4; p < j - VRNA_GQUAD_MIN_BOX_SIZE; p++){
l1 = p - i - 1;
if(l1>MAXLOOP) break;
if(S[p] != 3) continue;
c0 = en1 + ggg[index[q] + p] + P->internal_loop[l1 + 1];
ge = MIN2(ge, c0);
}
/* second case with 3' dangle (i.e. i+1) onto enclosing pair */
}
#endif
return ge;
}
INLINE PRIVATE
int *
E_GQuad_IntLoop_exhaustive( int i,
int j,
int **p_p,
int **q_p,
int type,
short *S,
int *ggg,
int threshold,
int *index,
paramT *P){
int energy, *ge, en1, en2, dangles, p, q, l1, minq, maxq;
int c0, c1, c2, c3, up, d53, d5, d3;
short si, sj;
int cnt = 0;
dangles = P->model_details.dangles;
si = S[i + 1];
sj = S[j - 1];
energy = 0;
if(dangles == 2)
energy += P->mismatchI[type][si][sj];
if(type > 2)
energy += P->TerminalAU;
/* guess how many gquads are possible in interval [i+1,j-1] */
*p_p = (int *)space(sizeof(int) * 256);
*q_p = (int *)space(sizeof(int) * 256);
ge = (int *)space(sizeof(int) * 256);
p = i + 1;
if(S[p] == 3){
if(p < j - VRNA_GQUAD_MIN_BOX_SIZE){
minq = j - i + p - MAXLOOP - 2;
c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
minq = MAX2(c0, minq);
c0 = j - 3;
maxq = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
maxq = MIN2(c0, maxq);
for(q = minq; q < maxq; q++){
if(S[q] != 3) continue;
c0 = energy + ggg[index[q] + p] + P->internal_loop[j - q - 1];
if(c0 <= threshold){
ge[cnt] = energy + P->internal_loop[j - q - 1];
(*p_p)[cnt] = p;
(*q_p)[cnt++] = q;
}
}
}
}
for(p = i + 2;
p < j - VRNA_GQUAD_MIN_BOX_SIZE;
p++){
l1 = p - i - 1;
if(l1>MAXLOOP) break;
if(S[p] != 3) continue;
minq = j - i + p - MAXLOOP - 2;
c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
minq = MAX2(c0, minq);
c0 = j - 1;
maxq = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
maxq = MIN2(c0, maxq);
for(q = minq; q < maxq; q++){
if(S[q] != 3) continue;
c0 = energy + ggg[index[q] + p] + P->internal_loop[l1 + j - q - 1];
if(c0 <= threshold){
ge[cnt] = energy + P->internal_loop[l1 + j - q - 1];
(*p_p)[cnt] = p;
(*q_p)[cnt++] = q;
}
}
}
q = j - 1;
if(S[q] == 3)
for(p = i + 4;
p < j - VRNA_GQUAD_MIN_BOX_SIZE;
p++){
l1 = p - i - 1;
if(l1>MAXLOOP) break;
if(S[p] != 3) continue;
c0 = energy + ggg[index[q] + p] + P->internal_loop[l1];
if(c0 <= threshold){
ge[cnt] = energy + P->internal_loop[l1];
(*p_p)[cnt] = p;
(*q_p)[cnt++] = q;
}
}
(*p_p)[cnt] = -1;
return ge;
}
INLINE PRIVATE
int
E_GQuad_IntLoop_L(int i,
int j,
int type,
short *S,
int **ggg,
int maxdist,
paramT *P){
int energy, ge, en1, en2, dangles, p, q, l1, minq, maxq;
int c0, c1, c2, c3, up, d53, d5, d3;
short si, sj;
dangles = P->model_details.dangles;
si = S[i + 1];
sj = S[j - 1];
energy = 0;
if(dangles == 2)
energy += P->mismatchI[type][si][sj];
if(type > 2)
energy += P->TerminalAU;
ge = INF;
p = i + 1;
if(S[p] == 3){
if(p < j - VRNA_GQUAD_MIN_BOX_SIZE){
minq = j - i + p - MAXLOOP - 2;
c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
minq = MAX2(c0, minq);
c0 = j - 3;
maxq = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
maxq = MIN2(c0, maxq);
for(q = minq; q < maxq; q++){
if(S[q] != 3) continue;
c0 = energy + ggg[p][q-p] + P->internal_loop[j - q - 1];
ge = MIN2(ge, c0);
}
}
}
for(p = i + 2;
p < j - VRNA_GQUAD_MIN_BOX_SIZE;
p++){
l1 = p - i - 1;
if(l1>MAXLOOP) break;
if(S[p] != 3) continue;
minq = j - i + p - MAXLOOP - 2;
c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
minq = MAX2(c0, minq);
c0 = j - 1;
maxq = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
maxq = MIN2(c0, maxq);
for(q = minq; q < maxq; q++){
if(S[q] != 3) continue;
c0 = energy + ggg[p][q - p] + P->internal_loop[l1 + j - q - 1];
ge = MIN2(ge, c0);
}
}
q = j - 1;
if(S[q] == 3)
for(p = i + 4;
p < j - VRNA_GQUAD_MIN_BOX_SIZE;
p++){
l1 = p - i - 1;
if(l1>MAXLOOP) break;
if(S[p] != 3) continue;
c0 = energy + ggg[p][q - p] + P->internal_loop[l1];
ge = MIN2(ge, c0);
}
return ge;
}
INLINE PRIVATE
FLT_OR_DBL
exp_E_GQuad_IntLoop(int i,
int j,
int type,
short *S,
FLT_OR_DBL *G,
int *index,
pf_paramT *pf){
int k, l, minl, maxl, u, r;
FLT_OR_DBL q, qe, *expintern;
short si, sj;
q = 0;
si = S[i + 1];
sj = S[j - 1];
qe = pf->expmismatchI[type][si][sj];
expintern = pf->expinternal;
if(type > 2)
qe *= pf->expTermAU;
k = i + 1;
if(S[k] == 3){
if(k < j - VRNA_GQUAD_MIN_BOX_SIZE){
minl = j - i + k - MAXLOOP - 2;
u = k + VRNA_GQUAD_MIN_BOX_SIZE - 1;
minl = MAX2(u, minl);
u = j - 3;
maxl = k + VRNA_GQUAD_MAX_BOX_SIZE + 1;
maxl = MIN2(u, maxl);
for(l = minl; l < maxl; l++){
if(S[l] != 3) continue;
if(G[index[k]-l] == 0.) continue;
q += qe * G[index[k]-l] * expintern[j - l - 1];
}
}
}
for(k = i + 2;
k <= j - VRNA_GQUAD_MIN_BOX_SIZE;
k++){
u = k - i - 1;
if(u > MAXLOOP) break;
if(S[k] != 3) continue;
minl = j - i + k - MAXLOOP - 2;
r = k + VRNA_GQUAD_MIN_BOX_SIZE - 1;
minl = MAX2(r, minl);
maxl = k + VRNA_GQUAD_MAX_BOX_SIZE + 1;
r = j - 1;
maxl = MIN2(r, maxl);
for(l = minl; l < maxl; l++){
if(S[l] != 3) continue;
if(G[index[k]-l] == 0.) continue;
q += qe * G[index[k]-l] * expintern[u + j - l - 1];
}
}
l = j - 1;
if(S[l] == 3)
for(k = i + 4; k < j - VRNA_GQUAD_MIN_BOX_SIZE; k++){
u = k - i - 1;
if(u>MAXLOOP) break;
if(S[k] != 3) continue;
if(G[index[k]-l] == 0.) continue;
q += qe * G[index[k]-l] * expintern[u];
}
return q;
}
#endif