binaries/src/ViennaRNA/lib/gquad.c

   1 /*
   2   gquad.c
   3
   4   Ronny Lorenz 2012
   5
   6   Vienna RNA package
   7 */
   8
   9 #include <stdio.h>
  10 #include <stdlib.h>
  11 #include <math.h>
  12 #include <string.h>
  13
  14 #include "../config.h"
  15 #include "fold_vars.h"
  16 #include "data_structures.h"
  17 #include "energy_const.h"
  18 #include "utils.h"
  19 #include "aln_util.h"
  20 #include "gquad.h"
  21
  22 #ifndef INLINE
  23 #ifdef __GNUC__
  24 # define INLINE inline
  25 #else
  26 # define INLINE
  27 #endif
  28 #endif
  29
  30 /**
  31  *  Use this macro to loop over each G-quadruplex
  32  *  delimited by a and b within the subsequence [c,d]
  33  */
  34 #define FOR_EACH_GQUAD(a, b, c, d)  \
  35           for((a) = (d) - VRNA_GQUAD_MIN_BOX_SIZE + 1; (a) >= (c); (a)--)\
  36             for((b) = (a) + VRNA_GQUAD_MIN_BOX_SIZE - 1;\
  37                 (b) <= MIN2((d), (a) + VRNA_GQUAD_MAX_BOX_SIZE - 1);\
  38                 (b)++)
  39
  40 /**
  41  *  This macro does almost the same as FOR_EACH_GQUAD() but keeps
  42  *  the 5' delimiter fixed. 'b' is the 3' delimiter of the gquad,
  43  *  for gquads within subsequence [a,c] that have 5' delimiter 'a'
  44  */
  45 #define FOR_EACH_GQUAD_AT(a, b, c)  \
  46           for((b) = (a) + VRNA_GQUAD_MIN_BOX_SIZE - 1;\
  47               (b) <= MIN2((c), (a) + VRNA_GQUAD_MAX_BOX_SIZE - 1);\
  48               (b)++)
  49
  50
  51 /*
  52 #################################
  53 # PRIVATE FUNCTION DECLARATIONS #
  54 #################################
  55 */
  56
  57 PRIVATE INLINE
  58 int *
  59 get_g_islands(short *S);
  60
  61 PRIVATE INLINE
  62 int *
  63 get_g_islands_sub(short *S, int i, int j);
  64
  65 /**
  66  *  IMPORTANT:
  67  *  If you don't know how to use this function, DONT'T USE IT!
  68  *
  69  *  The function pointer this function takes as argument is
  70  *  used for individual calculations with each g-quadruplex
  71  *  delimited by [i,j].
  72  *  The function it points to always receives as first 3 arguments
  73  *  position i, the stack size L and an array l[3] containing the
  74  *  individual linker sizes.
  75  *  The remaining 4 (void *) pointers of the callback function receive
  76  *  the parameters 'data', 'P', 'aux1' and 'aux2' and thus may be
  77  *  used to pass whatever data you like to.
  78  *  As the names of those parameters suggest the convention is that
  79  *  'data' should be used as a pointer where data is stored into,
  80  *  e.g the MFE or PF and the 'P' parameter should actually be a
  81  *  'paramT *' or 'pf_paramT *' type.
  82  *  However, what you actually pass obviously depends on the
  83  *  function the pointer is pointing to.
  84  *
  85  *  Although all of this may look like an overkill, it is found
  86  *  to be almost as fast as implementing g-quadruplex enumeration
  87  *  in each individual scenario, i.e. code duplication.
  88  *  Using this function, however, ensures that all g-quadruplex
  89  *  enumerations are absolutely identical.
  90  */
  91 PRIVATE
  92 void
  93 process_gquad_enumeration(int *gg,
  94                           int i,
  95                           int j,
  96                           void (*f)(int, int, int *,
  97                                     void *, void *, void *, void *),
  98                           void *data,
  99                           void *P,
 100                           void *aux1,
 101                           void *aux2);
 102
 103 /**
 104  *  MFE callback for process_gquad_enumeration()
 105  */
 106 PRIVATE
 107 void
 108 gquad_mfe(int i,
 109           int L,
 110           int *l,
 111           void *data,
 112           void *P,
 113           void *NA,
 114           void *NA2);
 115
 116 PRIVATE
 117 void
 118 gquad_mfe_pos(int i,
 119               int L,
 120               int *l,
 121               void *data,
 122               void *P,
 123               void *Lmfe,
 124               void *lmfe);
 125
 126 /**
 127  * Partition function callback for process_gquad_enumeration()
 128  */
 129 PRIVATE
 130 void
 131 gquad_pf( int i,
 132           int L,
 133           int *l,
 134           void *data,
 135           void *P,
 136           void *NA,
 137           void *NA2);
 138
 139 /**
 140  * Partition function callback for process_gquad_enumeration()
 141  * in contrast to gquad_pf() it stores the stack size L and
 142  * the linker lengths l[3] of the g-quadruplex that dominates
 143  * the interval [i,j]
 144  * (FLT_OR_DBL *)data must be 0. on entry
 145  */
 146 PRIVATE
 147 void
 148 gquad_pf_pos( int i,
 149               int L,
 150               int *l,
 151               void *data,
 152               void *pf,
 153               void *Lmax,
 154               void *lmax);
 155
 156 /**
 157  * MFE (alifold) callback for process_gquad_enumeration()
 158  */
 159 PRIVATE
 160 void
 161 gquad_mfe_ali(int i,
 162               int L,
 163               int *l,
 164               void *data,
 165               void *P,
 166               void *S,
 167               void *n_seq);
 168
 169 /**
 170  * MFE (alifold) callback for process_gquad_enumeration()
 171  * with seperation of free energy and penalty contribution
 172  */
 173 PRIVATE
 174 void
 175 gquad_mfe_ali_en( int i,
 176                   int L,
 177                   int *l,
 178                   void *data,
 179                   void *P,
 180                   void *S,
 181                   void *n_seq);
 182
 183 PRIVATE
 184 void
 185 gquad_interact( int i,
 186                 int L,
 187                 int *l,
 188                 void *data,
 189                 void *pf,
 190                 void *index,
 191                 void *NA2);
 192
 193 /* other useful static functions */
 194
 195 PRIVATE
 196 int
 197 gquad_ali_penalty(int i,
 198                   int L,
 199                   int l[3],
 200                   const short **S,
 201                   paramT *P);
 202
 203 /*
 204 #########################################
 205 # BEGIN OF PUBLIC FUNCTION DEFINITIONS  #
 206 #      (all available in RNAlib)        #
 207 #########################################
 208 */
 209
 210 /********************************
 211   Here are the G-quadruplex energy
 212   contribution functions
 213 *********************************/
 214
 215 PUBLIC int E_gquad( int L,
 216                     int l[3],
 217                     paramT *P){
 218
 219   int i, c = INF;
 220
 221   for(i=0;i<3;i++){
 222     if(l[i] > VRNA_GQUAD_MAX_LINKER_LENGTH) return c;
 223     if(l[i] < VRNA_GQUAD_MIN_LINKER_LENGTH) return c;
 224   }
 225   if(L > VRNA_GQUAD_MAX_STACK_SIZE) return c;
 226   if(L < VRNA_GQUAD_MIN_STACK_SIZE) return c;
 227
 228   gquad_mfe(0, L, l,
 229             (void *)(&c),
 230             (void *)P,
 231             NULL,
 232             NULL);
 233   return c;
 234 }
 235
 236 PUBLIC FLT_OR_DBL exp_E_gquad(int L,
 237                               int l[3],
 238                               pf_paramT *pf){
 239
 240   int i;
 241   FLT_OR_DBL q = 0.;
 242
 243   for(i=0;i<3;i++){
 244     if(l[i] > VRNA_GQUAD_MAX_LINKER_LENGTH) return q;
 245     if(l[i] < VRNA_GQUAD_MIN_LINKER_LENGTH) return q;
 246   }
 247   if(L > VRNA_GQUAD_MAX_STACK_SIZE) return q;
 248   if(L < VRNA_GQUAD_MIN_STACK_SIZE) return q;
 249
 250   gquad_pf( 0, L, l,
 251             (void *)(&q),
 252             (void *)pf,
 253             NULL,
 254             NULL);
 255   return q;
 256 }
 257
 258 PUBLIC int E_gquad_ali( int i,
 259                         int L,
 260                         int l[3],
 261                         const short **S,
 262                         int n_seq,
 263                         paramT *P){
 264
 265   int en[2];
 266   E_gquad_ali_en(i, L, l, S, n_seq, en, P);
 267   return en[0] + en[1];
 268 }
 269
 270
 271 PUBLIC void E_gquad_ali_en( int i,
 272                             int L,
 273                             int l[3],
 274                             const short **S,
 275                             int n_seq,
 276                             int en[2],
 277                             paramT *P){
 278
 279   int j;
 280   en[0] = en[1] = INF;
 281
 282   for(j=0;j<3;j++){
 283     if(l[j] > VRNA_GQUAD_MAX_LINKER_LENGTH) return;
 284     if(l[j] < VRNA_GQUAD_MIN_LINKER_LENGTH) return;
 285   }
 286   if(L > VRNA_GQUAD_MAX_STACK_SIZE) return;
 287   if(L < VRNA_GQUAD_MIN_STACK_SIZE) return;
 288
 289   gquad_mfe_ali_en( i, L, l,
 290                     (void *)(&(en[0])),
 291                     (void *)P,
 292                     (void *)S,
 293                     (void *)(&n_seq));
 294 }
 295
 296 /********************************
 297   Now, the triangular matrix
 298   generators for the G-quadruplex
 299   contributions are following
 300 *********************************/
 301
 302 PUBLIC int *get_gquad_matrix(short *S, paramT *P){
 303
 304   int n, size, i, j, *gg, *my_index, *data;
 305
 306   n         = S[0];
 307   my_index  = get_indx(n);
 308   gg        = get_g_islands(S);
 309   size      = (n * (n+1))/2 + 2;
 310   data      = (int *)space(sizeof(int) * size);
 311
 312   /* prefill the upper triangular matrix with INF */
 313   for(i = 0; i < size; i++) data[i] = INF;
 314
 315   FOR_EACH_GQUAD(i, j, 1, n){
 316     process_gquad_enumeration(gg, i, j,
 317                               &gquad_mfe,
 318                               (void *)(&(data[my_index[j]+i])),
 319                               (void *)P,
 320                               NULL,
 321                               NULL);
 322   }
 323
 324   free(my_index);
 325   free(gg);
 326   return data;
 327 }
 328
 329 PUBLIC FLT_OR_DBL *get_gquad_pf_matrix( short *S,
 330                                         FLT_OR_DBL *scale,
 331                                         pf_paramT *pf){
 332
 333   int n, size, *gg, i, j, *my_index;
 334   FLT_OR_DBL *data;
 335
 336
 337   n         = S[0];
 338   size      = (n * (n+1))/2 + 2;
 339   data      = (FLT_OR_DBL *)space(sizeof(FLT_OR_DBL) * size);
 340   gg        = get_g_islands(S);
 341   my_index  = get_iindx(n);
 342
 343   FOR_EACH_GQUAD(i, j, 1, n){
 344     process_gquad_enumeration(gg, i, j,
 345                               &gquad_pf,
 346                               (void *)(&(data[my_index[i]-j])),
 347                               (void *)pf,
 348                               NULL,
 349                               NULL);
 350     data[my_index[i]-j] *= scale[j-i+1];
 351   }
 352
 353   free(my_index);
 354   free(gg);
 355   return data;
 356 }
 357
 358 PUBLIC int *get_gquad_ali_matrix( short *S_cons,
 359                                   short **S,
 360                                   int n_seq,
 361                                   paramT *P){
 362
 363   int n, size, *data, *gg;
 364   int i, j, *my_index;
 365
 366
 367   n         = S[0][0];
 368   size      = (n * (n+1))/2 + 2;
 369   data      = (int *)space(sizeof(int) * size);
 370   gg        = get_g_islands(S_cons);
 371   my_index  = get_indx(n);
 372
 373   /* prefill the upper triangular matrix with INF */
 374   for(i=0;i<size;i++) data[i] = INF;
 375
 376   FOR_EACH_GQUAD(i, j, 1, n){
 377     process_gquad_enumeration(gg, i, j,
 378                               &gquad_mfe_ali,
 379                               (void *)(&(data[my_index[j]+i])),
 380                               (void *)P,
 381                               (void *)S,
 382                               (void *)(&n_seq));
 383   }
 384
 385   free(my_index);
 386   free(gg);
 387   return data;
 388 }
 389
 390 PUBLIC int **get_gquad_L_matrix(short *S,
 391                                 int start,
 392                                 int maxdist,
 393                                 int **g,
 394                                 paramT *P){
 395
 396   int **data;
 397   int n, i, j, k, l, *gg;
 398
 399   n   = S[0];
 400   gg  = get_g_islands_sub(S, start, MIN2(n, start + maxdist + 4));
 401
 402   if(g){ /* we just update the gquadruplex contribution for the current
 403             start and rotate the rest */
 404     data = g;
 405     /* we re-use the memory allocated previously */
 406     data[start] = data[start + maxdist + 5];
 407     data[start + maxdist + 5] = NULL;
 408
 409     /* prefill with INF */
 410     for(i = 0; i < maxdist + 5; i++)
 411       data[start][i] = INF;
 412
 413     /*  now we compute contributions for all gquads with 5' delimiter at
 414         position 'start'
 415     */
 416     FOR_EACH_GQUAD_AT(start, j, start + maxdist + 4){
 417       process_gquad_enumeration(gg, start, j,
 418                                 &gquad_mfe,
 419                                 (void *)(&(data[start][j-start])),
 420                                 (void *)P,
 421                                 NULL,
 422                                 NULL);
 423     }
 424
 425   } else { /* create a new matrix from scratch since this is the first
 426               call to this function */
 427
 428     /* allocate memory and prefill with INF */
 429     data = (int **) space(sizeof(int *) * (n+1));
 430     for(k = n; (k>n-maxdist-5) && (k>=0); k--){
 431       data[k] = (int *) space(sizeof(int)*(maxdist+5));
 432       for(i = 0; i < maxdist+5; i++) data[k][i] = INF;
 433     }
 434
 435     /* compute all contributions for the gquads in this interval */
 436     FOR_EACH_GQUAD(i, j, n - maxdist - 4, n){
 437       process_gquad_enumeration(gg, i, j,
 438                                 &gquad_mfe,
 439                                 (void *)(&(data[i][j-i])),
 440                                 (void *)P,
 441                                 NULL,
 442                                 NULL);
 443     }
 444   }
 445
 446   gg += start - 1;
 447   free(gg);
 448   return data;
 449 }
 450
 451 PUBLIC plist *get_plist_gquad_from_db(const char *structure, float pr){
 452   int x, size, actual_size, L, n, ge, ee, gb, l[3];
 453   plist *pl;
 454
 455   actual_size = 0;
 456   ge          = 0;
 457   n           = 2;
 458   size        = strlen(structure);
 459   pl          = (plist *)space(n*size*sizeof(plist));
 460
 461   while((ee = parse_gquad(structure + ge, &L, l)) > 0){
 462     ge += ee;
 463     gb = ge - L*4 - l[0] - l[1] - l[2] + 1;
 464     /* add pseudo-base pair encloding gquad */
 465     for(x = 0; x < L; x++){
 466       if (actual_size >= n * size - 5){
 467         n *= 2;
 468         pl = (plist *)xrealloc(pl, n * size * sizeof(plist));
 469       }
 470       pl[actual_size].i = gb + x;
 471       pl[actual_size].j = ge + x - L + 1;
 472       pl[actual_size].p = pr;
 473       pl[actual_size++].type = 0;
 474
 475       pl[actual_size].i = gb + x;
 476       pl[actual_size].j = gb + x + l[0] + L;
 477       pl[actual_size].p = pr;
 478       pl[actual_size++].type = 0;
 479
 480       pl[actual_size].i = gb + x + l[0] + L;
 481       pl[actual_size].j = ge + x - 2*L - l[2] + 1;
 482       pl[actual_size].p = pr;
 483       pl[actual_size++].type = 0;
 484
 485       pl[actual_size].i = ge + x - 2*L - l[2] + 1;
 486       pl[actual_size].j = ge + x - L + 1;
 487       pl[actual_size].p = pr;
 488       pl[actual_size++].type = 0;
 489     }
 490   }
 491
 492   pl[actual_size].i = pl[actual_size].j = 0;
 493   pl[actual_size++].p = 0;
 494   pl = (plist *)xrealloc(pl, actual_size * sizeof(plist));
 495   return pl;
 496 }
 497
 498 PUBLIC void get_gquad_pattern_mfe(short *S,
 499                                   int i,
 500                                   int j,
 501                                   paramT *P,
 502                                   int *L,
 503                                   int l[3]){
 504
 505   int *gg = get_g_islands_sub(S, i, j);
 506   int c = INF;
 507
 508   process_gquad_enumeration(gg, i, j,
 509                             &gquad_mfe_pos,
 510                             (void *)(&c),
 511                             (void *)P,
 512                             (void *)L,
 513                             (void *)l);
 514
 515   gg += i - 1;
 516   free(gg);
 517 }
 518
 519 PUBLIC void get_gquad_pattern_pf( short *S,
 520                                   int i,
 521                                   int j,
 522                                   pf_paramT *pf,
 523                                   int *L,
 524                                   int l[3]){
 525
 526   int *gg = get_g_islands_sub(S, i, j);
 527   FLT_OR_DBL q = 0.;
 528
 529   process_gquad_enumeration(gg, i, j,
 530                             &gquad_pf_pos,
 531                             (void *)(&q),
 532                             (void *)pf,
 533                             (void *)L,
 534                             (void *)l);
 535
 536   gg += i - 1;
 537   free(gg);
 538 }
 539
 540 PUBLIC plist *get_plist_gquad_from_pr(short *S,
 541                                       int gi,
 542                                       int gj,
 543                                       FLT_OR_DBL *G,
 544                                       FLT_OR_DBL *probs,
 545                                       FLT_OR_DBL *scale,
 546                                       pf_paramT *pf){
 547
 548   int L, l[3];
 549   return  get_plist_gquad_from_pr_max(S, gi, gj, G, probs, scale, &L, l, pf);
 550 }
 551
 552
 553 PUBLIC plist *get_plist_gquad_from_pr_max(short *S,
 554                                       int gi,
 555                                       int gj,
 556                                       FLT_OR_DBL *G,
 557                                       FLT_OR_DBL *probs,
 558                                       FLT_OR_DBL *scale,
 559                                       int *Lmax,
 560                                       int lmax[3],
 561                                       pf_paramT *pf){
 562
 563   int n, size, *gg, counter, i, j, *my_index;
 564   FLT_OR_DBL pp, *tempprobs;
 565   plist *pl;
 566
 567   n         = S[0];
 568   size      = (n * (n + 1))/2 + 2;
 569   tempprobs = (FLT_OR_DBL *)space(sizeof(FLT_OR_DBL) * size);
 570   pl        = (plist *)space((S[0]*S[0])*sizeof(plist));
 571   gg        = get_g_islands_sub(S, gi, gj);
 572   counter   = 0;
 573   my_index  = get_iindx(n);
 574
 575   process_gquad_enumeration(gg, gi, gj,
 576                             &gquad_interact,
 577                             (void *)tempprobs,
 578                             (void *)pf,
 579                             (void *)my_index,
 580                             NULL);
 581
 582   pp = 0.;
 583   process_gquad_enumeration(gg, gi, gj,
 584                             &gquad_pf_pos,
 585                             (void *)(&pp),
 586                             (void *)pf,
 587                             (void *)Lmax,
 588                             (void *)lmax);
 589
 590   pp = probs[my_index[gi]-gj] * scale[gj-gi+1] / G[my_index[gi]-gj];
 591   for (i=gi;i<gj; i++) {
 592     for (j=i; j<=gj; j++) {
 593       if (tempprobs[my_index[i]-j]>0.) {
 594         pl[counter].i=i;
 595         pl[counter].j=j;
 596         pl[counter++].p = pp * tempprobs[my_index[i]-j];
 597       }
 598     }
 599   }
 600   pl[counter].i = pl[counter].j = 0;
 601   pl[counter++].p = 0.;
 602   /* shrink memory to actual size needed */
 603   pl = (plist *) xrealloc(pl, counter * sizeof(plist));
 604
 605   gg += gi - 1; free(gg);
 606   free(my_index);
 607   free (tempprobs);
 608   return pl;
 609 }
 610
 611 PUBLIC int parse_gquad(const char *struc, int *L, int l[3]) {
 612   int i, il, start, end, len;
 613
 614   for (i=0; struc[i] && struc[i]!='+'; i++);
 615   if (struc[i] == '+') { /* start of gquad */
 616     for (il=0; il<=3; il++) {
 617       start=i; /* pos of first '+' */
 618       while (struc[++i] == '+'){
 619         if((il) && (i-start == *L))
 620           break;
 621       }
 622       end=i; len=end-start;
 623       if (il==0) *L=len;
 624       else if (len!=*L)
 625         nrerror("unequal stack lengths in gquad");
 626       if (il==3) break;
 627       while (struc[++i] == '.'); /* linker */
 628       l[il] = i-end;
 629       if (struc[i] != '+')
 630         nrerror("illegal character in gquad linker region");
 631     }
 632   }
 633   else return 0;
 634   /* printf("gquad at %d %d %d %d %d\n", end, *L, l[0], l[1], l[2]); */
 635   return end;
 636 }
 637
 638
 639
 640 /*
 641 #########################################
 642 # BEGIN OF PRIVATE FUNCTION DEFINITIONS #
 643 #          (internal use only)          #
 644 #########################################
 645 */
 646
 647 PRIVATE int gquad_ali_penalty(int i,
 648                               int L,
 649                               int l[3],
 650                               const short **S,
 651                               paramT *P){
 652
 653   int s, cnt;
 654   int penalty     = 0;
 655   int gg_mismatch = 0;
 656
 657   /* check for compatibility in the alignment */
 658   for(s = 0; S[s]; s++){
 659     unsigned int  ld  = 0; /* !=0 if layer destruction was detected */
 660     int           pen = 0;
 661
 662     /* check bottom layer */
 663     if(S[s][i] != 3)                            ld |= 1U;
 664     if(S[s][i + L + l[0]] != 3)                 ld |= 2U;
 665     if(S[s][i + 2*L + l[0] + l[1]] != 3)        ld |= 4U;
 666     if(S[s][i + 3*L + l[0] + l[1] + l[2]] != 3) ld |= 8U;
 667      /* add 1x penalty for missing bottom layer */
 668     if(ld) pen += VRNA_GQUAD_MISMATCH_PENALTY;
 669
 670     /* check top layer */
 671     ld = 0;
 672     if(S[s][i + L - 1] != 3)                        ld |= 1U;
 673     if(S[s][i + 2*L + l[0] - 1] != 3)               ld |= 2U;
 674     if(S[s][i + 3*L + l[0] + l[1] - 1] != 3)        ld |= 4U;
 675     if(S[s][i + 4*L + l[0] + l[1] + l[2] - 1] != 3) ld |= 8U;
 676      /* add 1x penalty for missing top layer */
 677     if(ld) pen += VRNA_GQUAD_MISMATCH_PENALTY;
 678
 679     /* check inner layers */
 680     for(cnt=1;cnt<L-1;cnt++){
 681       if(S[s][i + cnt] != 3)                            ld |= 1U;
 682       if(S[s][i + L + l[0] + cnt] != 3)                 ld |= 2U;
 683       if(S[s][i + 2*L + l[0] + l[1] + cnt] != 3)        ld |= 4U;
 684       if(S[s][i + 3*L + l[0] + l[1] + l[2] + cnt] != 3) ld |= 8U;
 685       /* add 2x penalty for missing inner layer */
 686       if(ld) pen += 2*VRNA_GQUAD_MISMATCH_PENALTY;
 687     }
 688
 689     /* if all layers are missing, we have a complete gg mismatch */
 690     if(pen >= (2*VRNA_GQUAD_MISMATCH_PENALTY * (L-1)))
 691       gg_mismatch++;
 692
 693     /* add the penalty to the score */
 694     penalty += pen;
 695   }
 696   /* if gg_mismatch exceeds maximum allowed, this g-quadruplex is forbidden */
 697   if(gg_mismatch > VRNA_GQUAD_MISMATCH_NUM_ALI) return INF;
 698   else return penalty;
 699 }
 700
 701
 702 PRIVATE void gquad_mfe( int i,
 703                         int L,
 704                         int *l,
 705                         void *data,
 706                         void *P,
 707                         void *NA,
 708                         void *NA2){
 709
 710   int cc = ((paramT *)P)->gquad[L][l[0] + l[1] + l[2]];
 711   if(cc < *((int *)data))
 712     *((int *)data) = cc;
 713 }
 714
 715 PRIVATE void gquad_mfe_pos( int i,
 716                             int L,
 717                             int *l,
 718                             void *data,
 719                             void *P,
 720                             void *Lmfe,
 721                             void *lmfe){
 722
 723   int cc = ((paramT *)P)->gquad[L][l[0] + l[1] + l[2]];
 724   if(cc < *((int *)data)){
 725     *((int *)data)        = cc;
 726     *((int *)Lmfe)        = L;
 727     *((int *)lmfe)        = l[0];
 728     *(((int *)lmfe) + 1)  = l[1];
 729     *(((int *)lmfe) + 2)  = l[2];
 730   }
 731 }
 732
 733 PRIVATE void gquad_pf(int i,
 734                       int L,
 735                       int *l,
 736                       void *data,
 737                       void *pf,
 738                       void *NA,
 739                       void *NA2){
 740
 741   *((FLT_OR_DBL *)data) += ((pf_paramT *)pf)->expgquad[L][l[0] + l[1] + l[2]];
 742 }
 743
 744 PRIVATE void gquad_pf_pos(int i,
 745                           int L,
 746                           int *l,
 747                           void *data,
 748                           void *pf,
 749                           void *Lmax,
 750                           void *lmax){
 751
 752   FLT_OR_DBL gq = ((pf_paramT *)pf)->expgquad[L][l[0] + l[1] + l[2]];
 753   if(gq > *((FLT_OR_DBL *)data)){
 754     *((FLT_OR_DBL *)data) = gq;
 755     *((int *)Lmax)        = L;
 756     *((int *)lmax)        = l[0];
 757     *(((int *)lmax) + 1)  = l[1];
 758     *(((int *)lmax) + 2)  = l[2];
 759   }
 760 }
 761
 762 PRIVATE void gquad_mfe_ali( int i,
 763                             int L,
 764                             int *l,
 765                             void *data,
 766                             void *P,
 767                             void *S,
 768                             void *n_seq){
 769
 770   int en[2], cc;
 771   en[0] = en[1] = INF;
 772   gquad_mfe_ali_en(i, L, l, (void *)(&(en[0])), P, S, n_seq);
 773   if(en[1] != INF){
 774     cc  = en[0] + en[1];
 775     if(cc < *((int *)data)) *((int *)data) = cc;
 776   }
 777 }
 778
 779 PRIVATE void gquad_mfe_ali_en(int i,
 780                               int L,
 781                               int *l,
 782                               void *data,
 783                               void *P,
 784                               void *S,
 785                               void *n_seq){
 786
 787   int en[2], cc, dd;
 788   en[0] = ((paramT *)P)->gquad[L][l[0] + l[1] + l[2]] * (*(int *)n_seq);
 789   en[1] = gquad_ali_penalty(i, L, l, (const short **)S, (paramT *)P);
 790   if(en[1] != INF){
 791     cc = en[0] + en[1];
 792     dd = ((int *)data)[0] + ((int *)data)[1];
 793     if(cc < dd){
 794       ((int *)data)[0] = en[0];
 795       ((int *)data)[1] = en[1];
 796     }
 797   }
 798 }
 799
 800 PRIVATE void gquad_interact(int i,
 801                       int L,
 802                       int *l,
 803                       void *data,
 804                       void *pf,
 805                       void *index,
 806                       void *NA2){
 807
 808   int x, *idx;
 809   FLT_OR_DBL gq, *pp;
 810
 811   idx = (int *)index;
 812   pp  = (FLT_OR_DBL *)data;
 813   gq  = exp_E_gquad(L, l, (pf_paramT *)pf);
 814
 815   for(x = 0; x < L; x++){
 816     pp[idx[i + x] - (i + x + 3*L + l[0] + l[1] + l[2])] += gq;
 817     pp[idx[i + x] - (i + x + L + l[0])] += gq;
 818     pp[idx[i + x + L + l[0]] - (i + x + 2*L + l[0] + l[1])] += gq;
 819     pp[idx[i + x + 2*L + l[0] + l[1]] - (i + x + 3*L + l[0] + l[1] + l[2])] += gq;
 820   }
 821
 822 }
 823
 824 PRIVATE INLINE int *get_g_islands(short *S){
 825   return get_g_islands_sub(S, 1, S[0]);
 826 }
 827
 828 PRIVATE INLINE int *get_g_islands_sub(short *S, int i, int j){
 829   int x, *gg;
 830
 831   gg = (int *)space(sizeof(int)*(j-i+2));
 832   gg -= i - 1;
 833
 834   if(S[j]==3) gg[j] = 1;
 835   for(x = j - 1; x >= i; x--)
 836     if(S[x] == 3)
 837       gg[x] = gg[x+1]+1;
 838
 839   return gg;
 840 }
 841
 842 /**
 843  *  We could've also created a macro that loops over all G-quadruplexes
 844  *  delimited by i and j. However, for the fun of it we use this function
 845  *  that receives a pointer to a callback function which in turn does the
 846  *  actual computation for each quadruplex found.
 847  */
 848 PRIVATE
 849 void
 850 process_gquad_enumeration(int *gg,
 851                           int i,
 852                           int j,
 853                           void (*f)(int, int, int *,
 854                                     void *, void *, void *, void *),
 855                           void *data,
 856                           void *P,
 857                           void *aux1,
 858                           void *aux2){
 859
 860   int L, l[3], n, max_linker, maxl0, maxl1;
 861
 862   n = j - i + 1;
 863
 864   if((n >= VRNA_GQUAD_MIN_BOX_SIZE) && (n <= VRNA_GQUAD_MAX_BOX_SIZE))
 865     for(L = MIN2(gg[i], VRNA_GQUAD_MAX_STACK_SIZE);
 866         L >= VRNA_GQUAD_MIN_STACK_SIZE;
 867         L--)
 868       if(gg[j-L+1] >= L){
 869         max_linker = n-4*L;
 870         if(     (max_linker >= 3*VRNA_GQUAD_MIN_LINKER_LENGTH)
 871             &&  (max_linker <= 3*VRNA_GQUAD_MAX_LINKER_LENGTH)){
 872           maxl0 = MIN2( VRNA_GQUAD_MAX_LINKER_LENGTH,
 873                         max_linker - 2*VRNA_GQUAD_MIN_LINKER_LENGTH
 874                       );
 875           for(l[0] = VRNA_GQUAD_MIN_LINKER_LENGTH;
 876               l[0] <= maxl0;
 877               l[0]++)
 878             if(gg[i+L+l[0]] >= L){
 879               maxl1 = MIN2( VRNA_GQUAD_MAX_LINKER_LENGTH,
 880                             max_linker - l[0] - VRNA_GQUAD_MIN_LINKER_LENGTH
 881                           );
 882               for(l[1] = VRNA_GQUAD_MIN_LINKER_LENGTH;
 883                   l[1] <= maxl1;
 884                   l[1]++)
 885                 if(gg[i + 2*L + l[0] + l[1]] >= L){
 886                   l[2] = max_linker - l[0] - l[1];
 887                   f(i, L, &(l[0]), data, P, aux1, aux2);
 888                 }
 889             }
 890         }
 891       }
 892 }
 893