--- /dev/null
+/*
+ gquad.c
+
+ Ronny Lorenz 2012
+
+ Vienna RNA package
+*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <string.h>
+
+#include "../config.h"
+#include "fold_vars.h"
+#include "data_structures.h"
+#include "energy_const.h"
+#include "utils.h"
+#include "aln_util.h"
+#include "gquad.h"
+
+#ifndef INLINE
+#ifdef __GNUC__
+# define INLINE inline
+#else
+# define INLINE
+#endif
+#endif
+
+/**
+ * Use this macro to loop over each G-quadruplex
+ * delimited by a and b within the subsequence [c,d]
+ */
+#define FOR_EACH_GQUAD(a, b, c, d) \
+ for((a) = (d) - VRNA_GQUAD_MIN_BOX_SIZE + 1; (a) >= (c); (a)--)\
+ for((b) = (a) + VRNA_GQUAD_MIN_BOX_SIZE - 1;\
+ (b) <= MIN2((d), (a) + VRNA_GQUAD_MAX_BOX_SIZE - 1);\
+ (b)++)
+
+/**
+ * This macro does almost the same as FOR_EACH_GQUAD() but keeps
+ * the 5' delimiter fixed. 'b' is the 3' delimiter of the gquad,
+ * for gquads within subsequence [a,c] that have 5' delimiter 'a'
+ */
+#define FOR_EACH_GQUAD_AT(a, b, c) \
+ for((b) = (a) + VRNA_GQUAD_MIN_BOX_SIZE - 1;\
+ (b) <= MIN2((c), (a) + VRNA_GQUAD_MAX_BOX_SIZE - 1);\
+ (b)++)
+
+
+/*
+#################################
+# PRIVATE FUNCTION DECLARATIONS #
+#################################
+*/
+
+PRIVATE INLINE
+int *
+get_g_islands(short *S);
+
+PRIVATE INLINE
+int *
+get_g_islands_sub(short *S, int i, int j);
+
+/**
+ * IMPORTANT:
+ * If you don't know how to use this function, DONT'T USE IT!
+ *
+ * The function pointer this function takes as argument is
+ * used for individual calculations with each g-quadruplex
+ * delimited by [i,j].
+ * The function it points to always receives as first 3 arguments
+ * position i, the stack size L and an array l[3] containing the
+ * individual linker sizes.
+ * The remaining 4 (void *) pointers of the callback function receive
+ * the parameters 'data', 'P', 'aux1' and 'aux2' and thus may be
+ * used to pass whatever data you like to.
+ * As the names of those parameters suggest the convention is that
+ * 'data' should be used as a pointer where data is stored into,
+ * e.g the MFE or PF and the 'P' parameter should actually be a
+ * 'paramT *' or 'pf_paramT *' type.
+ * However, what you actually pass obviously depends on the
+ * function the pointer is pointing to.
+ *
+ * Although all of this may look like an overkill, it is found
+ * to be almost as fast as implementing g-quadruplex enumeration
+ * in each individual scenario, i.e. code duplication.
+ * Using this function, however, ensures that all g-quadruplex
+ * enumerations are absolutely identical.
+ */
+PRIVATE
+void
+process_gquad_enumeration(int *gg,
+ int i,
+ int j,
+ void (*f)(int, int, int *,
+ void *, void *, void *, void *),
+ void *data,
+ void *P,
+ void *aux1,
+ void *aux2);
+
+/**
+ * MFE callback for process_gquad_enumeration()
+ */
+PRIVATE
+void
+gquad_mfe(int i,
+ int L,
+ int *l,
+ void *data,
+ void *P,
+ void *NA,
+ void *NA2);
+
+PRIVATE
+void
+gquad_mfe_pos(int i,
+ int L,
+ int *l,
+ void *data,
+ void *P,
+ void *Lmfe,
+ void *lmfe);
+
+/**
+ * Partition function callback for process_gquad_enumeration()
+ */
+PRIVATE
+void
+gquad_pf( int i,
+ int L,
+ int *l,
+ void *data,
+ void *P,
+ void *NA,
+ void *NA2);
+
+/**
+ * Partition function callback for process_gquad_enumeration()
+ * in contrast to gquad_pf() it stores the stack size L and
+ * the linker lengths l[3] of the g-quadruplex that dominates
+ * the interval [i,j]
+ * (FLT_OR_DBL *)data must be 0. on entry
+ */
+PRIVATE
+void
+gquad_pf_pos( int i,
+ int L,
+ int *l,
+ void *data,
+ void *pf,
+ void *Lmax,
+ void *lmax);
+
+/**
+ * MFE (alifold) callback for process_gquad_enumeration()
+ */
+PRIVATE
+void
+gquad_mfe_ali(int i,
+ int L,
+ int *l,
+ void *data,
+ void *P,
+ void *S,
+ void *n_seq);
+
+/**
+ * MFE (alifold) callback for process_gquad_enumeration()
+ * with seperation of free energy and penalty contribution
+ */
+PRIVATE
+void
+gquad_mfe_ali_en( int i,
+ int L,
+ int *l,
+ void *data,
+ void *P,
+ void *S,
+ void *n_seq);
+
+PRIVATE
+void
+gquad_interact( int i,
+ int L,
+ int *l,
+ void *data,
+ void *pf,
+ void *index,
+ void *NA2);
+
+/* other useful static functions */
+
+PRIVATE
+int
+gquad_ali_penalty(int i,
+ int L,
+ int l[3],
+ const short **S,
+ paramT *P);
+
+/*
+#########################################
+# BEGIN OF PUBLIC FUNCTION DEFINITIONS #
+# (all available in RNAlib) #
+#########################################
+*/
+
+/********************************
+ Here are the G-quadruplex energy
+ contribution functions
+*********************************/
+
+PUBLIC int E_gquad( int L,
+ int l[3],
+ paramT *P){
+
+ int i, c = INF;
+
+ for(i=0;i<3;i++){
+ if(l[i] > VRNA_GQUAD_MAX_LINKER_LENGTH) return c;
+ if(l[i] < VRNA_GQUAD_MIN_LINKER_LENGTH) return c;
+ }
+ if(L > VRNA_GQUAD_MAX_STACK_SIZE) return c;
+ if(L < VRNA_GQUAD_MIN_STACK_SIZE) return c;
+
+ gquad_mfe(0, L, l,
+ (void *)(&c),
+ (void *)P,
+ NULL,
+ NULL);
+ return c;
+}
+
+PUBLIC FLT_OR_DBL exp_E_gquad(int L,
+ int l[3],
+ pf_paramT *pf){
+
+ int i;
+ FLT_OR_DBL q = 0.;
+
+ for(i=0;i<3;i++){
+ if(l[i] > VRNA_GQUAD_MAX_LINKER_LENGTH) return q;
+ if(l[i] < VRNA_GQUAD_MIN_LINKER_LENGTH) return q;
+ }
+ if(L > VRNA_GQUAD_MAX_STACK_SIZE) return q;
+ if(L < VRNA_GQUAD_MIN_STACK_SIZE) return q;
+
+ gquad_pf( 0, L, l,
+ (void *)(&q),
+ (void *)pf,
+ NULL,
+ NULL);
+ return q;
+}
+
+PUBLIC int E_gquad_ali( int i,
+ int L,
+ int l[3],
+ const short **S,
+ int n_seq,
+ paramT *P){
+
+ int en[2];
+ E_gquad_ali_en(i, L, l, S, n_seq, en, P);
+ return en[0] + en[1];
+}
+
+
+PUBLIC void E_gquad_ali_en( int i,
+ int L,
+ int l[3],
+ const short **S,
+ int n_seq,
+ int en[2],
+ paramT *P){
+
+ int j;
+ en[0] = en[1] = INF;
+
+ for(j=0;j<3;j++){
+ if(l[j] > VRNA_GQUAD_MAX_LINKER_LENGTH) return;
+ if(l[j] < VRNA_GQUAD_MIN_LINKER_LENGTH) return;
+ }
+ if(L > VRNA_GQUAD_MAX_STACK_SIZE) return;
+ if(L < VRNA_GQUAD_MIN_STACK_SIZE) return;
+
+ gquad_mfe_ali_en( i, L, l,
+ (void *)(&(en[0])),
+ (void *)P,
+ (void *)S,
+ (void *)(&n_seq));
+}
+
+/********************************
+ Now, the triangular matrix
+ generators for the G-quadruplex
+ contributions are following
+*********************************/
+
+PUBLIC int *get_gquad_matrix(short *S, paramT *P){
+
+ int n, size, i, j, *gg, *my_index, *data;
+
+ n = S[0];
+ my_index = get_indx(n);
+ gg = get_g_islands(S);
+ size = (n * (n+1))/2 + 2;
+ data = (int *)space(sizeof(int) * size);
+
+ /* prefill the upper triangular matrix with INF */
+ for(i = 0; i < size; i++) data[i] = INF;
+
+ FOR_EACH_GQUAD(i, j, 1, n){
+ process_gquad_enumeration(gg, i, j,
+ &gquad_mfe,
+ (void *)(&(data[my_index[j]+i])),
+ (void *)P,
+ NULL,
+ NULL);
+ }
+
+ free(my_index);
+ free(gg);
+ return data;
+}
+
+PUBLIC FLT_OR_DBL *get_gquad_pf_matrix( short *S,
+ FLT_OR_DBL *scale,
+ pf_paramT *pf){
+
+ int n, size, *gg, i, j, *my_index;
+ FLT_OR_DBL *data;
+
+
+ n = S[0];
+ size = (n * (n+1))/2 + 2;
+ data = (FLT_OR_DBL *)space(sizeof(FLT_OR_DBL) * size);
+ gg = get_g_islands(S);
+ my_index = get_iindx(n);
+
+ FOR_EACH_GQUAD(i, j, 1, n){
+ process_gquad_enumeration(gg, i, j,
+ &gquad_pf,
+ (void *)(&(data[my_index[i]-j])),
+ (void *)pf,
+ NULL,
+ NULL);
+ data[my_index[i]-j] *= scale[j-i+1];
+ }
+
+ free(my_index);
+ free(gg);
+ return data;
+}
+
+PUBLIC int *get_gquad_ali_matrix( short *S_cons,
+ short **S,
+ int n_seq,
+ paramT *P){
+
+ int n, size, *data, *gg;
+ int i, j, *my_index;
+
+
+ n = S[0][0];
+ size = (n * (n+1))/2 + 2;
+ data = (int *)space(sizeof(int) * size);
+ gg = get_g_islands(S_cons);
+ my_index = get_indx(n);
+
+ /* prefill the upper triangular matrix with INF */
+ for(i=0;i<size;i++) data[i] = INF;
+
+ FOR_EACH_GQUAD(i, j, 1, n){
+ process_gquad_enumeration(gg, i, j,
+ &gquad_mfe_ali,
+ (void *)(&(data[my_index[j]+i])),
+ (void *)P,
+ (void *)S,
+ (void *)(&n_seq));
+ }
+
+ free(my_index);
+ free(gg);
+ return data;
+}
+
+PUBLIC int **get_gquad_L_matrix(short *S,
+ int start,
+ int maxdist,
+ int **g,
+ paramT *P){
+
+ int **data;
+ int n, i, j, k, l, *gg;
+
+ n = S[0];
+ gg = get_g_islands_sub(S, start, MIN2(n, start + maxdist + 4));
+
+ if(g){ /* we just update the gquadruplex contribution for the current
+ start and rotate the rest */
+ data = g;
+ /* we re-use the memory allocated previously */
+ data[start] = data[start + maxdist + 5];
+ data[start + maxdist + 5] = NULL;
+
+ /* prefill with INF */
+ for(i = 0; i < maxdist + 5; i++)
+ data[start][i] = INF;
+
+ /* now we compute contributions for all gquads with 5' delimiter at
+ position 'start'
+ */
+ FOR_EACH_GQUAD_AT(start, j, start + maxdist + 4){
+ process_gquad_enumeration(gg, start, j,
+ &gquad_mfe,
+ (void *)(&(data[start][j-start])),
+ (void *)P,
+ NULL,
+ NULL);
+ }
+
+ } else { /* create a new matrix from scratch since this is the first
+ call to this function */
+
+ /* allocate memory and prefill with INF */
+ data = (int **) space(sizeof(int *) * (n+1));
+ for(k = n; (k>n-maxdist-5) && (k>=0); k--){
+ data[k] = (int *) space(sizeof(int)*(maxdist+5));
+ for(i = 0; i < maxdist+5; i++) data[k][i] = INF;
+ }
+
+ /* compute all contributions for the gquads in this interval */
+ FOR_EACH_GQUAD(i, j, n - maxdist - 4, n){
+ process_gquad_enumeration(gg, i, j,
+ &gquad_mfe,
+ (void *)(&(data[i][j-i])),
+ (void *)P,
+ NULL,
+ NULL);
+ }
+ }
+
+ gg += start - 1;
+ free(gg);
+ return data;
+}
+
+PUBLIC plist *get_plist_gquad_from_db(const char *structure, float pr){
+ int x, size, actual_size, L, n, ge, ee, gb, l[3];
+ plist *pl;
+
+ actual_size = 0;
+ ge = 0;
+ n = 2;
+ size = strlen(structure);
+ pl = (plist *)space(n*size*sizeof(plist));
+
+ while((ee = parse_gquad(structure + ge, &L, l)) > 0){
+ ge += ee;
+ gb = ge - L*4 - l[0] - l[1] - l[2] + 1;
+ /* add pseudo-base pair encloding gquad */
+ for(x = 0; x < L; x++){
+ if (actual_size >= n * size - 5){
+ n *= 2;
+ pl = (plist *)xrealloc(pl, n * size * sizeof(plist));
+ }
+ pl[actual_size].i = gb + x;
+ pl[actual_size].j = ge + x - L + 1;
+ pl[actual_size].p = pr;
+ pl[actual_size++].type = 0;
+
+ pl[actual_size].i = gb + x;
+ pl[actual_size].j = gb + x + l[0] + L;
+ pl[actual_size].p = pr;
+ pl[actual_size++].type = 0;
+
+ pl[actual_size].i = gb + x + l[0] + L;
+ pl[actual_size].j = ge + x - 2*L - l[2] + 1;
+ pl[actual_size].p = pr;
+ pl[actual_size++].type = 0;
+
+ pl[actual_size].i = ge + x - 2*L - l[2] + 1;
+ pl[actual_size].j = ge + x - L + 1;
+ pl[actual_size].p = pr;
+ pl[actual_size++].type = 0;
+ }
+ }
+
+ pl[actual_size].i = pl[actual_size].j = 0;
+ pl[actual_size++].p = 0;
+ pl = (plist *)xrealloc(pl, actual_size * sizeof(plist));
+ return pl;
+}
+
+PUBLIC void get_gquad_pattern_mfe(short *S,
+ int i,
+ int j,
+ paramT *P,
+ int *L,
+ int l[3]){
+
+ int *gg = get_g_islands_sub(S, i, j);
+ int c = INF;
+
+ process_gquad_enumeration(gg, i, j,
+ &gquad_mfe_pos,
+ (void *)(&c),
+ (void *)P,
+ (void *)L,
+ (void *)l);
+
+ gg += i - 1;
+ free(gg);
+}
+
+PUBLIC void get_gquad_pattern_pf( short *S,
+ int i,
+ int j,
+ pf_paramT *pf,
+ int *L,
+ int l[3]){
+
+ int *gg = get_g_islands_sub(S, i, j);
+ FLT_OR_DBL q = 0.;
+
+ process_gquad_enumeration(gg, i, j,
+ &gquad_pf_pos,
+ (void *)(&q),
+ (void *)pf,
+ (void *)L,
+ (void *)l);
+
+ gg += i - 1;
+ free(gg);
+}
+
+PUBLIC 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){
+
+ int L, l[3];
+ return get_plist_gquad_from_pr_max(S, gi, gj, G, probs, scale, &L, l, pf);
+}
+
+
+PUBLIC 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 *Lmax,
+ int lmax[3],
+ pf_paramT *pf){
+
+ int n, size, *gg, counter, i, j, *my_index;
+ FLT_OR_DBL pp, *tempprobs;
+ plist *pl;
+
+ n = S[0];
+ size = (n * (n + 1))/2 + 2;
+ tempprobs = (FLT_OR_DBL *)space(sizeof(FLT_OR_DBL) * size);
+ pl = (plist *)space((S[0]*S[0])*sizeof(plist));
+ gg = get_g_islands_sub(S, gi, gj);
+ counter = 0;
+ my_index = get_iindx(n);
+
+ process_gquad_enumeration(gg, gi, gj,
+ &gquad_interact,
+ (void *)tempprobs,
+ (void *)pf,
+ (void *)my_index,
+ NULL);
+
+ pp = 0.;
+ process_gquad_enumeration(gg, gi, gj,
+ &gquad_pf_pos,
+ (void *)(&pp),
+ (void *)pf,
+ (void *)Lmax,
+ (void *)lmax);
+
+ pp = probs[my_index[gi]-gj] * scale[gj-gi+1] / G[my_index[gi]-gj];
+ for (i=gi;i<gj; i++) {
+ for (j=i; j<=gj; j++) {
+ if (tempprobs[my_index[i]-j]>0.) {
+ pl[counter].i=i;
+ pl[counter].j=j;
+ pl[counter++].p = pp * tempprobs[my_index[i]-j];
+ }
+ }
+ }
+ pl[counter].i = pl[counter].j = 0;
+ pl[counter++].p = 0.;
+ /* shrink memory to actual size needed */
+ pl = (plist *) xrealloc(pl, counter * sizeof(plist));
+
+ gg += gi - 1; free(gg);
+ free(my_index);
+ free (tempprobs);
+ return pl;
+}
+
+PUBLIC int parse_gquad(const char *struc, int *L, int l[3]) {
+ int i, il, start, end, len;
+
+ for (i=0; struc[i] && struc[i]!='+'; i++);
+ if (struc[i] == '+') { /* start of gquad */
+ for (il=0; il<=3; il++) {
+ start=i; /* pos of first '+' */
+ while (struc[++i] == '+'){
+ if((il) && (i-start == *L))
+ break;
+ }
+ end=i; len=end-start;
+ if (il==0) *L=len;
+ else if (len!=*L)
+ nrerror("unequal stack lengths in gquad");
+ if (il==3) break;
+ while (struc[++i] == '.'); /* linker */
+ l[il] = i-end;
+ if (struc[i] != '+')
+ nrerror("illegal character in gquad linker region");
+ }
+ }
+ else return 0;
+ /* printf("gquad at %d %d %d %d %d\n", end, *L, l[0], l[1], l[2]); */
+ return end;
+}
+
+
+
+/*
+#########################################
+# BEGIN OF PRIVATE FUNCTION DEFINITIONS #
+# (internal use only) #
+#########################################
+*/
+
+PRIVATE int gquad_ali_penalty(int i,
+ int L,
+ int l[3],
+ const short **S,
+ paramT *P){
+
+ int s, cnt;
+ int penalty = 0;
+ int gg_mismatch = 0;
+
+ /* check for compatibility in the alignment */
+ for(s = 0; S[s]; s++){
+ unsigned int ld = 0; /* !=0 if layer destruction was detected */
+ int pen = 0;
+
+ /* check bottom layer */
+ if(S[s][i] != 3) ld |= 1U;
+ if(S[s][i + L + l[0]] != 3) ld |= 2U;
+ if(S[s][i + 2*L + l[0] + l[1]] != 3) ld |= 4U;
+ if(S[s][i + 3*L + l[0] + l[1] + l[2]] != 3) ld |= 8U;
+ /* add 1x penalty for missing bottom layer */
+ if(ld) pen += VRNA_GQUAD_MISMATCH_PENALTY;
+
+ /* check top layer */
+ ld = 0;
+ if(S[s][i + L - 1] != 3) ld |= 1U;
+ if(S[s][i + 2*L + l[0] - 1] != 3) ld |= 2U;
+ if(S[s][i + 3*L + l[0] + l[1] - 1] != 3) ld |= 4U;
+ if(S[s][i + 4*L + l[0] + l[1] + l[2] - 1] != 3) ld |= 8U;
+ /* add 1x penalty for missing top layer */
+ if(ld) pen += VRNA_GQUAD_MISMATCH_PENALTY;
+
+ /* check inner layers */
+ for(cnt=1;cnt<L-1;cnt++){
+ if(S[s][i + cnt] != 3) ld |= 1U;
+ if(S[s][i + L + l[0] + cnt] != 3) ld |= 2U;
+ if(S[s][i + 2*L + l[0] + l[1] + cnt] != 3) ld |= 4U;
+ if(S[s][i + 3*L + l[0] + l[1] + l[2] + cnt] != 3) ld |= 8U;
+ /* add 2x penalty for missing inner layer */
+ if(ld) pen += 2*VRNA_GQUAD_MISMATCH_PENALTY;
+ }
+
+ /* if all layers are missing, we have a complete gg mismatch */
+ if(pen >= (2*VRNA_GQUAD_MISMATCH_PENALTY * (L-1)))
+ gg_mismatch++;
+
+ /* add the penalty to the score */
+ penalty += pen;
+ }
+ /* if gg_mismatch exceeds maximum allowed, this g-quadruplex is forbidden */
+ if(gg_mismatch > VRNA_GQUAD_MISMATCH_NUM_ALI) return INF;
+ else return penalty;
+}
+
+
+PRIVATE void gquad_mfe( int i,
+ int L,
+ int *l,
+ void *data,
+ void *P,
+ void *NA,
+ void *NA2){
+
+ int cc = ((paramT *)P)->gquad[L][l[0] + l[1] + l[2]];
+ if(cc < *((int *)data))
+ *((int *)data) = cc;
+}
+
+PRIVATE void gquad_mfe_pos( int i,
+ int L,
+ int *l,
+ void *data,
+ void *P,
+ void *Lmfe,
+ void *lmfe){
+
+ int cc = ((paramT *)P)->gquad[L][l[0] + l[1] + l[2]];
+ if(cc < *((int *)data)){
+ *((int *)data) = cc;
+ *((int *)Lmfe) = L;
+ *((int *)lmfe) = l[0];
+ *(((int *)lmfe) + 1) = l[1];
+ *(((int *)lmfe) + 2) = l[2];
+ }
+}
+
+PRIVATE void gquad_pf(int i,
+ int L,
+ int *l,
+ void *data,
+ void *pf,
+ void *NA,
+ void *NA2){
+
+ *((FLT_OR_DBL *)data) += ((pf_paramT *)pf)->expgquad[L][l[0] + l[1] + l[2]];
+}
+
+PRIVATE void gquad_pf_pos(int i,
+ int L,
+ int *l,
+ void *data,
+ void *pf,
+ void *Lmax,
+ void *lmax){
+
+ FLT_OR_DBL gq = ((pf_paramT *)pf)->expgquad[L][l[0] + l[1] + l[2]];
+ if(gq > *((FLT_OR_DBL *)data)){
+ *((FLT_OR_DBL *)data) = gq;
+ *((int *)Lmax) = L;
+ *((int *)lmax) = l[0];
+ *(((int *)lmax) + 1) = l[1];
+ *(((int *)lmax) + 2) = l[2];
+ }
+}
+
+PRIVATE void gquad_mfe_ali( int i,
+ int L,
+ int *l,
+ void *data,
+ void *P,
+ void *S,
+ void *n_seq){
+
+ int en[2], cc;
+ en[0] = en[1] = INF;
+ gquad_mfe_ali_en(i, L, l, (void *)(&(en[0])), P, S, n_seq);
+ if(en[1] != INF){
+ cc = en[0] + en[1];
+ if(cc < *((int *)data)) *((int *)data) = cc;
+ }
+}
+
+PRIVATE void gquad_mfe_ali_en(int i,
+ int L,
+ int *l,
+ void *data,
+ void *P,
+ void *S,
+ void *n_seq){
+
+ int en[2], cc, dd;
+ en[0] = ((paramT *)P)->gquad[L][l[0] + l[1] + l[2]] * (*(int *)n_seq);
+ en[1] = gquad_ali_penalty(i, L, l, (const short **)S, (paramT *)P);
+ if(en[1] != INF){
+ cc = en[0] + en[1];
+ dd = ((int *)data)[0] + ((int *)data)[1];
+ if(cc < dd){
+ ((int *)data)[0] = en[0];
+ ((int *)data)[1] = en[1];
+ }
+ }
+}
+
+PRIVATE void gquad_interact(int i,
+ int L,
+ int *l,
+ void *data,
+ void *pf,
+ void *index,
+ void *NA2){
+
+ int x, *idx;
+ FLT_OR_DBL gq, *pp;
+
+ idx = (int *)index;
+ pp = (FLT_OR_DBL *)data;
+ gq = exp_E_gquad(L, l, (pf_paramT *)pf);
+
+ for(x = 0; x < L; x++){
+ pp[idx[i + x] - (i + x + 3*L + l[0] + l[1] + l[2])] += gq;
+ pp[idx[i + x] - (i + x + L + l[0])] += gq;
+ pp[idx[i + x + L + l[0]] - (i + x + 2*L + l[0] + l[1])] += gq;
+ pp[idx[i + x + 2*L + l[0] + l[1]] - (i + x + 3*L + l[0] + l[1] + l[2])] += gq;
+ }
+
+}
+
+PRIVATE INLINE int *get_g_islands(short *S){
+ return get_g_islands_sub(S, 1, S[0]);
+}
+
+PRIVATE INLINE int *get_g_islands_sub(short *S, int i, int j){
+ int x, *gg;
+
+ gg = (int *)space(sizeof(int)*(j-i+2));
+ gg -= i - 1;
+
+ if(S[j]==3) gg[j] = 1;
+ for(x = j - 1; x >= i; x--)
+ if(S[x] == 3)
+ gg[x] = gg[x+1]+1;
+
+ return gg;
+}
+
+/**
+ * We could've also created a macro that loops over all G-quadruplexes
+ * delimited by i and j. However, for the fun of it we use this function
+ * that receives a pointer to a callback function which in turn does the
+ * actual computation for each quadruplex found.
+ */
+PRIVATE
+void
+process_gquad_enumeration(int *gg,
+ int i,
+ int j,
+ void (*f)(int, int, int *,
+ void *, void *, void *, void *),
+ void *data,
+ void *P,
+ void *aux1,
+ void *aux2){
+
+ int L, l[3], n, max_linker, maxl0, maxl1;
+
+ n = j - i + 1;
+
+ if((n >= VRNA_GQUAD_MIN_BOX_SIZE) && (n <= VRNA_GQUAD_MAX_BOX_SIZE))
+ for(L = MIN2(gg[i], VRNA_GQUAD_MAX_STACK_SIZE);
+ L >= VRNA_GQUAD_MIN_STACK_SIZE;
+ L--)
+ if(gg[j-L+1] >= L){
+ max_linker = n-4*L;
+ if( (max_linker >= 3*VRNA_GQUAD_MIN_LINKER_LENGTH)
+ && (max_linker <= 3*VRNA_GQUAD_MAX_LINKER_LENGTH)){
+ maxl0 = MIN2( VRNA_GQUAD_MAX_LINKER_LENGTH,
+ max_linker - 2*VRNA_GQUAD_MIN_LINKER_LENGTH
+ );
+ for(l[0] = VRNA_GQUAD_MIN_LINKER_LENGTH;
+ l[0] <= maxl0;
+ l[0]++)
+ if(gg[i+L+l[0]] >= L){
+ maxl1 = MIN2( VRNA_GQUAD_MAX_LINKER_LENGTH,
+ max_linker - l[0] - VRNA_GQUAD_MIN_LINKER_LENGTH
+ );
+ for(l[1] = VRNA_GQUAD_MIN_LINKER_LENGTH;
+ l[1] <= maxl1;
+ l[1]++)
+ if(gg[i + 2*L + l[0] + l[1]] >= L){
+ l[2] = max_linker - l[0] - l[1];
+ f(i, L, &(l[0]), data, P, aux1, aux2);
+ }
+ }
+ }
+ }
+}
+
git://source.jalview.org / jabaws.git / blobdiff