binaries/src/ViennaRNA/H/gquad.h

   1 #ifndef __VIENNA_RNA_PACKAGE_GQUAD_H__
   2 #define __VIENNA_RNA_PACKAGE_GQUAD_H__
   3
   4 #include "data_structures.h"
   5
   6 #ifndef INLINE
   7 #ifdef __GNUC__
   8 # define INLINE inline
   9 #else
  10 # define INLINE
  11 #endif
  12 #endif
  13
  14 /**
  15  *  \file gquad.h
  16  *  \brief Various functions related to G-quadruplex computations
  17  */
  18
  19
  20 int         E_gquad(int L,
  21                     int l[3],
  22                     paramT *P);
  23
  24 FLT_OR_DBL exp_E_gquad( int L,
  25                         int l[3],
  26                         pf_paramT *pf);
  27
  28 int         E_gquad_ali(int i,
  29                         int L,
  30                         int l[3],
  31                         const short **S,
  32                         int n_seq,
  33                         paramT *P);
  34
  35
  36 void        E_gquad_ali_en( int i,
  37                             int L,
  38                             int l[3],
  39                             const short **S,
  40                             int n_seq,
  41                             int en[2],
  42                             paramT *P);
  43
  44 /**
  45  *  \brief Get a triangular matrix prefilled with minimum free energy
  46  *  contributions of G-quadruplexes.
  47  *
  48  *  At each position ij in the matrix, the minimum free energy of any
  49  *  G-quadruplex delimited by i and j is stored. If no G-quadruplex formation
  50  *  is possible, the matrix element is set to INF.
  51  *  Access the elements in the matrix via matrix[indx[j]+i]. To get
  52  *  the integer array indx see get_jindx().
  53  *
  54  *  \see get_jindx(), encode_sequence()
  55  *
  56  *  \param S  The encoded sequence
  57  *  \param P  A pointer to the data structure containing the precomputed energy contributions
  58  *  \return   A pointer to the G-quadruplex contribution matrix
  59 */
  60 int         *get_gquad_matrix(short *S, paramT *P);
  61
  62 int         *get_gquad_ali_matrix(short *S_cons,
  63                                   short **S,
  64                                   int n_seq,
  65                                   paramT *P);
  66
  67 FLT_OR_DBL  *get_gquad_pf_matrix( short *S,
  68                                   FLT_OR_DBL *scale,
  69                                   pf_paramT *pf);
  70
  71 int         **get_gquad_L_matrix( short *S,
  72                                   int start,
  73                                   int maxdist,
  74                                   int **g,
  75                                   paramT *P);
  76
  77 void        get_gquad_pattern_mfe(short *S,
  78                                   int i,
  79                                   int j,
  80                                   paramT *P,
  81                                   int *L,
  82                                   int l[3]);
  83
  84 void        get_gquad_pattern_pf( short *S,
  85                                   int i,
  86                                   int j,
  87                                   pf_paramT *pf,
  88                                   int *L,
  89                                   int l[3]);
  90
  91 plist       *get_plist_gquad_from_pr( short *S,
  92                                       int gi,
  93                                       int gj,
  94                                       FLT_OR_DBL *G,
  95                                       FLT_OR_DBL *probs,
  96                                       FLT_OR_DBL *scale,
  97                                       pf_paramT *pf);
  98 plist       *get_plist_gquad_from_pr_max(short *S,
  99                                       int gi,
 100                                       int gj,
 101                                       FLT_OR_DBL *G,
 102                                       FLT_OR_DBL *probs,
 103                                       FLT_OR_DBL *scale,
 104                                       int *L,
 105                                       int l[3],
 106                                       pf_paramT *pf);
 107
 108 plist       *get_plist_gquad_from_db( const char *structure,
 109                                       float pr);
 110
 111 /**
 112  *  given a dot-bracket structure (possibly) containing gquads encoded
 113  *  by '+' signs, find first gquad, return end position or 0 if none found
 114  *  Upon return L and l[] contain the number of stacked layers, as well as
 115  *  the lengths of the linker regions.
 116  *  To parse a string with many gquads, call parse_gquad repeatedly e.g.
 117  *  end1 = parse_gquad(struc, &L, l); ... ;
 118  *  end2 = parse_gquad(struc+end1, &L, l); end2+=end1; ... ;
 119  *  end3 = parse_gquad(struc+end2, &L, l); end3+=end2; ... ;
 120  */
 121 int         parse_gquad(const char *struc, int *L, int l[3]);
 122
 123
 124
 125 /**
 126  *  backtrack an interior loop like enclosed g-quadruplex
 127  *  with closing pair (i,j)
 128  *
 129  *  \param c      The total contribution the loop should resemble
 130  *  \param i      position i of enclosing pair
 131  *  \param j      position j of enclosing pair
 132  *  \param type   base pair type of enclosing pair (must be reverse type)
 133  *  \param S      integer encoded sequence
 134  *  \param ggg    triangular matrix containing g-quadruplex contributions
 135  *  \param index  the index for accessing the triangular matrix
 136  *  \param p      here the 5' position of the gquad is stored
 137  *  \param q      here the 3' position of the gquad is stored
 138  *  \param P      the datastructure containing the precalculated contibutions
 139  *
 140  *  \return       1 on success, 0 if no gquad found
 141  */
 142 INLINE  PRIVATE int backtrack_GQuad_IntLoop(int c,
 143                                             int i,
 144                                             int j,
 145                                             int type,
 146                                             short *S,
 147                                             int *ggg,
 148                                             int *index,
 149                                             int *p,
 150                                             int *q,
 151                                             paramT *P){
 152
 153   int energy, dangles, k, l, maxl, minl, c0, l1;
 154   short si, sj;
 155
 156   dangles = P->model_details.dangles;
 157   si      = S[i + 1];
 158   sj      = S[j - 1];
 159   energy  = 0;
 160
 161   if(dangles == 2)
 162     energy += P->mismatchI[type][si][sj];
 163
 164   if(type > 2)
 165     energy += P->TerminalAU;
 166
 167   k = i + 1;
 168   if(S[k] == 3){
 169     if(k < j - VRNA_GQUAD_MIN_BOX_SIZE){
 170       minl  = j - i + k - MAXLOOP - 2;
 171       c0    = k + VRNA_GQUAD_MIN_BOX_SIZE - 1;
 172       minl  = MAX2(c0, minl);
 173       c0    = j - 3;
 174       maxl  = k + VRNA_GQUAD_MAX_BOX_SIZE + 1;
 175       maxl  = MIN2(c0, maxl);
 176       for(l = minl; l < maxl; l++){
 177         if(S[l] != 3) continue;
 178         if(c == energy + ggg[index[l] + k] + P->internal_loop[j - l - 1]){
 179           *p = k; *q = l;
 180           return 1;
 181         }
 182       }
 183     }
 184   }
 185
 186   for(k = i + 2;
 187       k < j - VRNA_GQUAD_MIN_BOX_SIZE;
 188       k++){
 189     l1    = k - i - 1;
 190     if(l1>MAXLOOP) break;
 191     if(S[k] != 3) continue;
 192     minl  = j - i + k - MAXLOOP - 2;
 193     c0    = k + VRNA_GQUAD_MIN_BOX_SIZE - 1;
 194     minl  = MAX2(c0, minl);
 195     c0    = j - 1;
 196     maxl  = k + VRNA_GQUAD_MAX_BOX_SIZE + 1;
 197     maxl  = MIN2(c0, maxl);
 198     for(l = minl; l < maxl; l++){
 199       if(S[l] != 3) continue;
 200       if(c == energy + ggg[index[l] + k] + P->internal_loop[l1 + j - l - 1]){
 201         *p = k; *q = l;
 202         return 1;
 203       }
 204     }
 205   }
 206
 207   l = j - 1;
 208   if(S[l] == 3)
 209     for(k = i + 4;
 210         k < j - VRNA_GQUAD_MIN_BOX_SIZE;
 211         k++){
 212       l1    = k - i - 1;
 213       if(l1>MAXLOOP) break;
 214       if(S[k] != 3) continue;
 215       if(c == energy + ggg[index[l] + k] + P->internal_loop[l1]){
 216         *p = k; *q = l;
 217         return 1;
 218       }
 219     }
 220
 221   return 0;
 222 }
 223
 224 /**
 225  *  backtrack an interior loop like enclosed g-quadruplex
 226  *  with closing pair (i,j) with underlying Lfold matrix
 227  *
 228  *  \param c      The total contribution the loop should resemble
 229  *  \param i      position i of enclosing pair
 230  *  \param j      position j of enclosing pair
 231  *  \param type   base pair type of enclosing pair (must be reverse type)
 232  *  \param S      integer encoded sequence
 233  *  \param ggg    triangular matrix containing g-quadruplex contributions
 234  *  \param p      here the 5' position of the gquad is stored
 235  *  \param q      here the 3' position of the gquad is stored
 236  *  \param P      the datastructure containing the precalculated contibutions
 237  *
 238  *  \return       1 on success, 0 if no gquad found
 239  */
 240 INLINE  PRIVATE int backtrack_GQuad_IntLoop_L(int c,
 241                                               int i,
 242                                               int j,
 243                                               int type,
 244                                               short *S,
 245                                               int **ggg,
 246                                               int maxdist,
 247                                               int *p,
 248                                               int *q,
 249                                               paramT *P){
 250
 251   int energy, dangles, k, l, maxl, minl, c0, l1;
 252   short si, sj;
 253
 254   dangles = P->model_details.dangles;
 255   si      = S[i + 1];
 256   sj      = S[j - 1];
 257   energy  = 0;
 258
 259   if(dangles == 2)
 260     energy += P->mismatchI[type][si][sj];
 261
 262   if(type > 2)
 263     energy += P->TerminalAU;
 264
 265   k = i + 1;
 266   if(S[k] == 3){
 267     if(k < j - VRNA_GQUAD_MIN_BOX_SIZE){
 268       minl  = j - i + k - MAXLOOP - 2;
 269       c0    = k + VRNA_GQUAD_MIN_BOX_SIZE - 1;
 270       minl  = MAX2(c0, minl);
 271       c0    = j - 3;
 272       maxl  = k + VRNA_GQUAD_MAX_BOX_SIZE + 1;
 273       maxl  = MIN2(c0, maxl);
 274       for(l = minl; l < maxl; l++){
 275         if(S[l] != 3) continue;
 276         if(c == energy + ggg[k][l - k] + P->internal_loop[j - l - 1]){
 277           *p = k; *q = l;
 278           return 1;
 279         }
 280       }
 281     }
 282   }
 283
 284   for(k = i + 2;
 285       k < j - VRNA_GQUAD_MIN_BOX_SIZE;
 286       k++){
 287     l1    = k - i - 1;
 288     if(l1>MAXLOOP) break;
 289     if(S[k] != 3) continue;
 290     minl  = j - i + k - MAXLOOP - 2;
 291     c0    = k + VRNA_GQUAD_MIN_BOX_SIZE - 1;
 292     minl  = MAX2(c0, minl);
 293     c0    = j - 1;
 294     maxl  = k + VRNA_GQUAD_MAX_BOX_SIZE + 1;
 295     maxl  = MIN2(c0, maxl);
 296     for(l = minl; l < maxl; l++){
 297       if(S[l] != 3) continue;
 298       if(c == energy + ggg[k][l - k] + P->internal_loop[l1 + j - l - 1]){
 299         *p = k; *q = l;
 300         return 1;
 301       }
 302     }
 303   }
 304
 305   l = j - 1;
 306   if(S[l] == 3)
 307     for(k = i + 4;
 308         k < j - VRNA_GQUAD_MIN_BOX_SIZE;
 309         k++){
 310       l1    = k - i - 1;
 311       if(l1>MAXLOOP) break;
 312       if(S[k] != 3) continue;
 313       if(c == energy + ggg[k][l - k] + P->internal_loop[l1]){
 314         *p = k; *q = l;
 315         return 1;
 316       }
 317     }
 318
 319   return 0;
 320 }
 321
 322 INLINE PRIVATE
 323 int
 324 E_GQuad_IntLoop(int i,
 325                 int j,
 326                 int type,
 327                 short *S,
 328                 int *ggg,
 329                 int *index,
 330                 paramT *P){
 331
 332   int energy, ge, en1, en2, dangles, p, q, l1, minq, maxq;
 333   int c0, c1, c2, c3, up, d53, d5, d3;
 334   short si, sj;
 335
 336   dangles = P->model_details.dangles;
 337   si      = S[i + 1];
 338   sj      = S[j - 1];
 339   energy  = 0;
 340
 341   if(dangles == 2)
 342     energy += P->mismatchI[type][si][sj];
 343
 344   if(type > 2)
 345     energy += P->TerminalAU;
 346
 347   ge = INF;
 348
 349   p = i + 1;
 350   if(S[p] == 3){
 351     if(p < j - VRNA_GQUAD_MIN_BOX_SIZE){
 352       minq  = j - i + p - MAXLOOP - 2;
 353       c0    = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
 354       minq  = MAX2(c0, minq);
 355       c0    = j - 3;
 356       maxq  = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
 357       maxq  = MIN2(c0, maxq);
 358       for(q = minq; q < maxq; q++){
 359         if(S[q] != 3) continue;
 360         c0  = energy + ggg[index[q] + p] + P->internal_loop[j - q - 1];
 361         ge  = MIN2(ge, c0);
 362       }
 363     }
 364   }
 365
 366   for(p = i + 2;
 367       p < j - VRNA_GQUAD_MIN_BOX_SIZE;
 368       p++){
 369     l1    = p - i - 1;
 370     if(l1>MAXLOOP) break;
 371     if(S[p] != 3) continue;
 372     minq  = j - i + p - MAXLOOP - 2;
 373     c0    = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
 374     minq  = MAX2(c0, minq);
 375     c0    = j - 1;
 376     maxq  = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
 377     maxq  = MIN2(c0, maxq);
 378     for(q = minq; q < maxq; q++){
 379       if(S[q] != 3) continue;
 380       c0  = energy + ggg[index[q] + p] + P->internal_loop[l1 + j - q - 1];
 381       ge   = MIN2(ge, c0);
 382     }
 383   }
 384
 385   q = j - 1;
 386   if(S[q] == 3)
 387     for(p = i + 4;
 388         p < j - VRNA_GQUAD_MIN_BOX_SIZE;
 389         p++){
 390       l1    = p - i - 1;
 391       if(l1>MAXLOOP) break;
 392       if(S[p] != 3) continue;
 393       c0  = energy + ggg[index[q] + p] + P->internal_loop[l1];
 394       ge  = MIN2(ge, c0);
 395     }
 396
 397 #if 0
 398   /* here comes the additional stuff for the odd dangle models */
 399   if(dangles % 1){
 400     en1 = energy + P->dangle5[type][si];
 401     en2 = energy + P->dangle5[type][sj];
 402     en3 = energy + P->mismatchI[type][si][sj];
 403
 404     /* first case with 5' dangle (i.e. j-1) onto enclosing pair */
 405     p = i + 1;
 406     if(S[p] == 3){
 407       if(p < j - VRNA_GQUAD_MIN_BOX_SIZE){
 408         minq  = j - i + p - MAXLOOP - 2;
 409         c0    = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
 410         minq  = MAX2(c0, minq);
 411         c0    = j - 4;
 412         maxq  = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
 413         maxq  = MIN2(c0, maxq);
 414         for(q = minq; q < maxq; q++){
 415           if(S[q] != 3) continue;
 416           c0  = en1 + ggg[index[q] + p] + P->internal_loop[j - q - 1];
 417           ge  = MIN2(ge, c0);
 418         }
 419       }
 420     }
 421
 422     for(p = i + 2; p < j - VRNA_GQUAD_MIN_BOX_SIZE; p++){
 423       l1    = p - i - 1;
 424       if(l1>MAXLOOP) break;
 425       if(S[p] != 3) continue;
 426       minq  = j - i + p - MAXLOOP - 2;
 427       c0    = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
 428       minq  = MAX2(c0, minq);
 429       c0    = j - 2;
 430       maxq  = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
 431       maxq  = MIN2(c0, maxq);
 432       for(q = minq; q < maxq; q++){
 433         if(S[q] != 3) continue;
 434         c0  = en1 + ggg[index[q] + p] + P->internal_loop[l1 + j - q - 1];
 435         ge   = MIN2(ge, c0);
 436       }
 437     }
 438
 439     q = j - 2;
 440     if(S[q] == 3)
 441       for(p = i + 4; p < j - VRNA_GQUAD_MIN_BOX_SIZE; p++){
 442         l1    = p - i - 1;
 443         if(l1>MAXLOOP) break;
 444         if(S[p] != 3) continue;
 445         c0  = en1 + ggg[index[q] + p] + P->internal_loop[l1 + 1];
 446         ge  = MIN2(ge, c0);
 447       }
 448
 449     /* second case with 3' dangle (i.e. i+1) onto enclosing pair */
 450
 451   }
 452 #endif
 453   return ge;
 454 }
 455
 456 INLINE PRIVATE
 457 int
 458 E_GQuad_IntLoop_L(int i,
 459                   int j,
 460                   int type,
 461                   short *S,
 462                   int **ggg,
 463                   int maxdist,
 464                   paramT *P){
 465
 466   int energy, ge, en1, en2, dangles, p, q, l1, minq, maxq;
 467   int c0, c1, c2, c3, up, d53, d5, d3;
 468   short si, sj;
 469
 470   dangles = P->model_details.dangles;
 471   si      = S[i + 1];
 472   sj      = S[j - 1];
 473   energy  = 0;
 474
 475   if(dangles == 2)
 476     energy += P->mismatchI[type][si][sj];
 477
 478   if(type > 2)
 479     energy += P->TerminalAU;
 480
 481   ge = INF;
 482
 483   p = i + 1;
 484   if(S[p] == 3){
 485     if(p < j - VRNA_GQUAD_MIN_BOX_SIZE){
 486       minq  = j - i + p - MAXLOOP - 2;
 487       c0    = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
 488       minq  = MAX2(c0, minq);
 489       c0    = j - 3;
 490       maxq  = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
 491       maxq  = MIN2(c0, maxq);
 492       for(q = minq; q < maxq; q++){
 493         if(S[q] != 3) continue;
 494         c0  = energy + ggg[p][q-p] + P->internal_loop[j - q - 1];
 495         ge  = MIN2(ge, c0);
 496       }
 497     }
 498   }
 499
 500   for(p = i + 2;
 501       p < j - VRNA_GQUAD_MIN_BOX_SIZE;
 502       p++){
 503     l1    = p - i - 1;
 504     if(l1>MAXLOOP) break;
 505     if(S[p] != 3) continue;
 506     minq  = j - i + p - MAXLOOP - 2;
 507     c0    = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
 508     minq  = MAX2(c0, minq);
 509     c0    = j - 1;
 510     maxq  = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
 511     maxq  = MIN2(c0, maxq);
 512     for(q = minq; q < maxq; q++){
 513       if(S[q] != 3) continue;
 514       c0  = energy + ggg[p][q - p] + P->internal_loop[l1 + j - q - 1];
 515       ge   = MIN2(ge, c0);
 516     }
 517   }
 518
 519   q = j - 1;
 520   if(S[q] == 3)
 521     for(p = i + 4;
 522         p < j - VRNA_GQUAD_MIN_BOX_SIZE;
 523         p++){
 524       l1    = p - i - 1;
 525       if(l1>MAXLOOP) break;
 526       if(S[p] != 3) continue;
 527       c0  = energy + ggg[p][q - p] + P->internal_loop[l1];
 528       ge  = MIN2(ge, c0);
 529     }
 530
 531   return ge;
 532 }
 533
 534 INLINE PRIVATE
 535 FLT_OR_DBL
 536 exp_E_GQuad_IntLoop(int i,
 537                     int j,
 538                     int type,
 539                     short *S,
 540                     FLT_OR_DBL *G,
 541                     int *index,
 542                     pf_paramT *pf){
 543
 544   int k, l, minl, maxl, u, r;
 545   FLT_OR_DBL q, qe, *expintern;
 546   short si, sj;
 547
 548   q         = 0;
 549   si        = S[i + 1];
 550   sj        = S[j - 1];
 551   qe        = pf->expmismatchI[type][si][sj];
 552   expintern = pf->expinternal;
 553
 554   if(type > 2)
 555     qe *= pf->expTermAU;
 556
 557   k = i + 1;
 558   if(S[k] == 3){
 559     if(k < j - VRNA_GQUAD_MIN_BOX_SIZE){
 560       minl  = j - i + k - MAXLOOP - 2;
 561       u     = k + VRNA_GQUAD_MIN_BOX_SIZE - 1;
 562       minl  = MAX2(u, minl);
 563       u     = j - 3;
 564       maxl  = k + VRNA_GQUAD_MAX_BOX_SIZE + 1;
 565       maxl  = MIN2(u, maxl);
 566       for(l = minl; l < maxl; l++){
 567         if(S[l] != 3) continue;
 568         q += qe * G[index[k]-l] * expintern[j - l - 1];
 569       }
 570     }
 571   }
 572
 573
 574   for(k = i + 2;
 575       k <= j - VRNA_GQUAD_MIN_BOX_SIZE;
 576       k++){
 577     u = k - i - 1;
 578     if(u > MAXLOOP) break;
 579     if(S[k] != 3) continue;
 580     minl  = j - i + k - MAXLOOP - 2;
 581     r     = k + VRNA_GQUAD_MIN_BOX_SIZE - 1;
 582     minl  = MAX2(r, minl);
 583     maxl  = k + VRNA_GQUAD_MAX_BOX_SIZE + 1;
 584     r     = j - 1;
 585     maxl  = MIN2(r, maxl);
 586     for(l = minl; l < maxl; l++){
 587       if(S[l] != 3) continue;
 588       q += qe * G[index[k]-l] * expintern[u + j - l - 1];
 589     }
 590   }
 591
 592   l = j - 1;
 593   if(S[l] == 3)
 594     for(k = i + 4; k < j - VRNA_GQUAD_MIN_BOX_SIZE; k++){
 595       u    = k - i - 1;
 596       if(u>MAXLOOP) break;
 597       if(S[k] != 3) continue;
 598       q += qe * G[index[k]-l] * expintern[u];
 599     }
 600
 601   return q;
 602 }
 603
 604 #endif