--- /dev/null
+#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_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);
+
+/**
+ * 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_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;
+ 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;
+ 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;
+ q += qe * G[index[k]-l] * expintern[u];
+ }
+
+ return q;
+}
+
+#endif
git://source.jalview.org / jabaws.git / blobdiff