00001 #ifndef __VIENNA_RNA_PACKAGE_GQUAD_H__
00002 #define __VIENNA_RNA_PACKAGE_GQUAD_H__
00003
00004 #include "data_structures.h"
00005
00006 #ifndef INLINE
00007 #ifdef __GNUC__
00008 # define INLINE inline
00009 #else
00010 # define INLINE
00011 #endif
00012 #endif
00013
00020 int E_gquad(int L,
00021 int l[3],
00022 paramT *P);
00023
00024 FLT_OR_DBL exp_E_gquad( int L,
00025 int l[3],
00026 pf_paramT *pf);
00027
00028 int E_gquad_ali(int i,
00029 int L,
00030 int l[3],
00031 const short **S,
00032 int n_seq,
00033 paramT *P);
00034
00035
00036 void E_gquad_ali_en( int i,
00037 int L,
00038 int l[3],
00039 const short **S,
00040 int n_seq,
00041 int en[2],
00042 paramT *P);
00043
00060 int *get_gquad_matrix(short *S, paramT *P);
00061
00062 int *get_gquad_ali_matrix(short *S_cons,
00063 short **S,
00064 int n_seq,
00065 paramT *P);
00066
00067 FLT_OR_DBL *get_gquad_pf_matrix( short *S,
00068 FLT_OR_DBL *scale,
00069 pf_paramT *pf);
00070
00071 int **get_gquad_L_matrix( short *S,
00072 int start,
00073 int maxdist,
00074 int **g,
00075 paramT *P);
00076
00077 void get_gquad_pattern_mfe(short *S,
00078 int i,
00079 int j,
00080 paramT *P,
00081 int *L,
00082 int l[3]);
00083
00084 void get_gquad_pattern_pf( short *S,
00085 int i,
00086 int j,
00087 pf_paramT *pf,
00088 int *L,
00089 int l[3]);
00090
00091 plist *get_plist_gquad_from_pr( short *S,
00092 int gi,
00093 int gj,
00094 FLT_OR_DBL *G,
00095 FLT_OR_DBL *probs,
00096 FLT_OR_DBL *scale,
00097 pf_paramT *pf);
00098 plist *get_plist_gquad_from_pr_max(short *S,
00099 int gi,
00100 int gj,
00101 FLT_OR_DBL *G,
00102 FLT_OR_DBL *probs,
00103 FLT_OR_DBL *scale,
00104 int *L,
00105 int l[3],
00106 pf_paramT *pf);
00107
00108 plist *get_plist_gquad_from_db( const char *structure,
00109 float pr);
00110
00121 int parse_gquad(const char *struc, int *L, int l[3]);
00122
00123
00124
00142 INLINE PRIVATE int backtrack_GQuad_IntLoop(int c,
00143 int i,
00144 int j,
00145 int type,
00146 short *S,
00147 int *ggg,
00148 int *index,
00149 int *p,
00150 int *q,
00151 paramT *P){
00152
00153 int energy, dangles, k, l, maxl, minl, c0, l1;
00154 short si, sj;
00155
00156 dangles = P->model_details.dangles;
00157 si = S[i + 1];
00158 sj = S[j - 1];
00159 energy = 0;
00160
00161 if(dangles == 2)
00162 energy += P->mismatchI[type][si][sj];
00163
00164 if(type > 2)
00165 energy += P->TerminalAU;
00166
00167 k = i + 1;
00168 if(S[k] == 3){
00169 if(k < j - VRNA_GQUAD_MIN_BOX_SIZE){
00170 minl = j - i + k - MAXLOOP - 2;
00171 c0 = k + VRNA_GQUAD_MIN_BOX_SIZE - 1;
00172 minl = MAX2(c0, minl);
00173 c0 = j - 3;
00174 maxl = k + VRNA_GQUAD_MAX_BOX_SIZE + 1;
00175 maxl = MIN2(c0, maxl);
00176 for(l = minl; l < maxl; l++){
00177 if(S[l] != 3) continue;
00178 if(c == energy + ggg[index[l] + k] + P->internal_loop[j - l - 1]){
00179 *p = k; *q = l;
00180 return 1;
00181 }
00182 }
00183 }
00184 }
00185
00186 for(k = i + 2;
00187 k < j - VRNA_GQUAD_MIN_BOX_SIZE;
00188 k++){
00189 l1 = k - i - 1;
00190 if(l1>MAXLOOP) break;
00191 if(S[k] != 3) continue;
00192 minl = j - i + k - MAXLOOP - 2;
00193 c0 = k + VRNA_GQUAD_MIN_BOX_SIZE - 1;
00194 minl = MAX2(c0, minl);
00195 c0 = j - 1;
00196 maxl = k + VRNA_GQUAD_MAX_BOX_SIZE + 1;
00197 maxl = MIN2(c0, maxl);
00198 for(l = minl; l < maxl; l++){
00199 if(S[l] != 3) continue;
00200 if(c == energy + ggg[index[l] + k] + P->internal_loop[l1 + j - l - 1]){
00201 *p = k; *q = l;
00202 return 1;
00203 }
00204 }
00205 }
00206
00207 l = j - 1;
00208 if(S[l] == 3)
00209 for(k = i + 4;
00210 k < j - VRNA_GQUAD_MIN_BOX_SIZE;
00211 k++){
00212 l1 = k - i - 1;
00213 if(l1>MAXLOOP) break;
00214 if(S[k] != 3) continue;
00215 if(c == energy + ggg[index[l] + k] + P->internal_loop[l1]){
00216 *p = k; *q = l;
00217 return 1;
00218 }
00219 }
00220
00221 return 0;
00222 }
00223
00240 INLINE PRIVATE int backtrack_GQuad_IntLoop_L(int c,
00241 int i,
00242 int j,
00243 int type,
00244 short *S,
00245 int **ggg,
00246 int maxdist,
00247 int *p,
00248 int *q,
00249 paramT *P){
00250
00251 int energy, dangles, k, l, maxl, minl, c0, l1;
00252 short si, sj;
00253
00254 dangles = P->model_details.dangles;
00255 si = S[i + 1];
00256 sj = S[j - 1];
00257 energy = 0;
00258
00259 if(dangles == 2)
00260 energy += P->mismatchI[type][si][sj];
00261
00262 if(type > 2)
00263 energy += P->TerminalAU;
00264
00265 k = i + 1;
00266 if(S[k] == 3){
00267 if(k < j - VRNA_GQUAD_MIN_BOX_SIZE){
00268 minl = j - i + k - MAXLOOP - 2;
00269 c0 = k + VRNA_GQUAD_MIN_BOX_SIZE - 1;
00270 minl = MAX2(c0, minl);
00271 c0 = j - 3;
00272 maxl = k + VRNA_GQUAD_MAX_BOX_SIZE + 1;
00273 maxl = MIN2(c0, maxl);
00274 for(l = minl; l < maxl; l++){
00275 if(S[l] != 3) continue;
00276 if(c == energy + ggg[k][l - k] + P->internal_loop[j - l - 1]){
00277 *p = k; *q = l;
00278 return 1;
00279 }
00280 }
00281 }
00282 }
00283
00284 for(k = i + 2;
00285 k < j - VRNA_GQUAD_MIN_BOX_SIZE;
00286 k++){
00287 l1 = k - i - 1;
00288 if(l1>MAXLOOP) break;
00289 if(S[k] != 3) continue;
00290 minl = j - i + k - MAXLOOP - 2;
00291 c0 = k + VRNA_GQUAD_MIN_BOX_SIZE - 1;
00292 minl = MAX2(c0, minl);
00293 c0 = j - 1;
00294 maxl = k + VRNA_GQUAD_MAX_BOX_SIZE + 1;
00295 maxl = MIN2(c0, maxl);
00296 for(l = minl; l < maxl; l++){
00297 if(S[l] != 3) continue;
00298 if(c == energy + ggg[k][l - k] + P->internal_loop[l1 + j - l - 1]){
00299 *p = k; *q = l;
00300 return 1;
00301 }
00302 }
00303 }
00304
00305 l = j - 1;
00306 if(S[l] == 3)
00307 for(k = i + 4;
00308 k < j - VRNA_GQUAD_MIN_BOX_SIZE;
00309 k++){
00310 l1 = k - i - 1;
00311 if(l1>MAXLOOP) break;
00312 if(S[k] != 3) continue;
00313 if(c == energy + ggg[k][l - k] + P->internal_loop[l1]){
00314 *p = k; *q = l;
00315 return 1;
00316 }
00317 }
00318
00319 return 0;
00320 }
00321
00322 INLINE PRIVATE
00323 int
00324 E_GQuad_IntLoop(int i,
00325 int j,
00326 int type,
00327 short *S,
00328 int *ggg,
00329 int *index,
00330 paramT *P){
00331
00332 int energy, ge, en1, en2, dangles, p, q, l1, minq, maxq;
00333 int c0, c1, c2, c3, up, d53, d5, d3;
00334 short si, sj;
00335
00336 dangles = P->model_details.dangles;
00337 si = S[i + 1];
00338 sj = S[j - 1];
00339 energy = 0;
00340
00341 if(dangles == 2)
00342 energy += P->mismatchI[type][si][sj];
00343
00344 if(type > 2)
00345 energy += P->TerminalAU;
00346
00347 ge = INF;
00348
00349 p = i + 1;
00350 if(S[p] == 3){
00351 if(p < j - VRNA_GQUAD_MIN_BOX_SIZE){
00352 minq = j - i + p - MAXLOOP - 2;
00353 c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
00354 minq = MAX2(c0, minq);
00355 c0 = j - 3;
00356 maxq = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
00357 maxq = MIN2(c0, maxq);
00358 for(q = minq; q < maxq; q++){
00359 if(S[q] != 3) continue;
00360 c0 = energy + ggg[index[q] + p] + P->internal_loop[j - q - 1];
00361 ge = MIN2(ge, c0);
00362 }
00363 }
00364 }
00365
00366 for(p = i + 2;
00367 p < j - VRNA_GQUAD_MIN_BOX_SIZE;
00368 p++){
00369 l1 = p - i - 1;
00370 if(l1>MAXLOOP) break;
00371 if(S[p] != 3) continue;
00372 minq = j - i + p - MAXLOOP - 2;
00373 c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
00374 minq = MAX2(c0, minq);
00375 c0 = j - 1;
00376 maxq = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
00377 maxq = MIN2(c0, maxq);
00378 for(q = minq; q < maxq; q++){
00379 if(S[q] != 3) continue;
00380 c0 = energy + ggg[index[q] + p] + P->internal_loop[l1 + j - q - 1];
00381 ge = MIN2(ge, c0);
00382 }
00383 }
00384
00385 q = j - 1;
00386 if(S[q] == 3)
00387 for(p = i + 4;
00388 p < j - VRNA_GQUAD_MIN_BOX_SIZE;
00389 p++){
00390 l1 = p - i - 1;
00391 if(l1>MAXLOOP) break;
00392 if(S[p] != 3) continue;
00393 c0 = energy + ggg[index[q] + p] + P->internal_loop[l1];
00394 ge = MIN2(ge, c0);
00395 }
00396
00397 #if 0
00398
00399 if(dangles % 1){
00400 en1 = energy + P->dangle5[type][si];
00401 en2 = energy + P->dangle5[type][sj];
00402 en3 = energy + P->mismatchI[type][si][sj];
00403
00404
00405 p = i + 1;
00406 if(S[p] == 3){
00407 if(p < j - VRNA_GQUAD_MIN_BOX_SIZE){
00408 minq = j - i + p - MAXLOOP - 2;
00409 c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
00410 minq = MAX2(c0, minq);
00411 c0 = j - 4;
00412 maxq = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
00413 maxq = MIN2(c0, maxq);
00414 for(q = minq; q < maxq; q++){
00415 if(S[q] != 3) continue;
00416 c0 = en1 + ggg[index[q] + p] + P->internal_loop[j - q - 1];
00417 ge = MIN2(ge, c0);
00418 }
00419 }
00420 }
00421
00422 for(p = i + 2; p < j - VRNA_GQUAD_MIN_BOX_SIZE; p++){
00423 l1 = p - i - 1;
00424 if(l1>MAXLOOP) break;
00425 if(S[p] != 3) continue;
00426 minq = j - i + p - MAXLOOP - 2;
00427 c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
00428 minq = MAX2(c0, minq);
00429 c0 = j - 2;
00430 maxq = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
00431 maxq = MIN2(c0, maxq);
00432 for(q = minq; q < maxq; q++){
00433 if(S[q] != 3) continue;
00434 c0 = en1 + ggg[index[q] + p] + P->internal_loop[l1 + j - q - 1];
00435 ge = MIN2(ge, c0);
00436 }
00437 }
00438
00439 q = j - 2;
00440 if(S[q] == 3)
00441 for(p = i + 4; p < j - VRNA_GQUAD_MIN_BOX_SIZE; p++){
00442 l1 = p - i - 1;
00443 if(l1>MAXLOOP) break;
00444 if(S[p] != 3) continue;
00445 c0 = en1 + ggg[index[q] + p] + P->internal_loop[l1 + 1];
00446 ge = MIN2(ge, c0);
00447 }
00448
00449
00450
00451 }
00452 #endif
00453 return ge;
00454 }
00455
00456 INLINE PRIVATE
00457 int
00458 E_GQuad_IntLoop_L(int i,
00459 int j,
00460 int type,
00461 short *S,
00462 int **ggg,
00463 int maxdist,
00464 paramT *P){
00465
00466 int energy, ge, en1, en2, dangles, p, q, l1, minq, maxq;
00467 int c0, c1, c2, c3, up, d53, d5, d3;
00468 short si, sj;
00469
00470 dangles = P->model_details.dangles;
00471 si = S[i + 1];
00472 sj = S[j - 1];
00473 energy = 0;
00474
00475 if(dangles == 2)
00476 energy += P->mismatchI[type][si][sj];
00477
00478 if(type > 2)
00479 energy += P->TerminalAU;
00480
00481 ge = INF;
00482
00483 p = i + 1;
00484 if(S[p] == 3){
00485 if(p < j - VRNA_GQUAD_MIN_BOX_SIZE){
00486 minq = j - i + p - MAXLOOP - 2;
00487 c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
00488 minq = MAX2(c0, minq);
00489 c0 = j - 3;
00490 maxq = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
00491 maxq = MIN2(c0, maxq);
00492 for(q = minq; q < maxq; q++){
00493 if(S[q] != 3) continue;
00494 c0 = energy + ggg[p][q-p] + P->internal_loop[j - q - 1];
00495 ge = MIN2(ge, c0);
00496 }
00497 }
00498 }
00499
00500 for(p = i + 2;
00501 p < j - VRNA_GQUAD_MIN_BOX_SIZE;
00502 p++){
00503 l1 = p - i - 1;
00504 if(l1>MAXLOOP) break;
00505 if(S[p] != 3) continue;
00506 minq = j - i + p - MAXLOOP - 2;
00507 c0 = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
00508 minq = MAX2(c0, minq);
00509 c0 = j - 1;
00510 maxq = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
00511 maxq = MIN2(c0, maxq);
00512 for(q = minq; q < maxq; q++){
00513 if(S[q] != 3) continue;
00514 c0 = energy + ggg[p][q - p] + P->internal_loop[l1 + j - q - 1];
00515 ge = MIN2(ge, c0);
00516 }
00517 }
00518
00519 q = j - 1;
00520 if(S[q] == 3)
00521 for(p = i + 4;
00522 p < j - VRNA_GQUAD_MIN_BOX_SIZE;
00523 p++){
00524 l1 = p - i - 1;
00525 if(l1>MAXLOOP) break;
00526 if(S[p] != 3) continue;
00527 c0 = energy + ggg[p][q - p] + P->internal_loop[l1];
00528 ge = MIN2(ge, c0);
00529 }
00530
00531 return ge;
00532 }
00533
00534 INLINE PRIVATE
00535 FLT_OR_DBL
00536 exp_E_GQuad_IntLoop(int i,
00537 int j,
00538 int type,
00539 short *S,
00540 FLT_OR_DBL *G,
00541 int *index,
00542 pf_paramT *pf){
00543
00544 int k, l, minl, maxl, u, r;
00545 FLT_OR_DBL q, qe, *expintern;
00546 short si, sj;
00547
00548 q = 0;
00549 si = S[i + 1];
00550 sj = S[j - 1];
00551 qe = pf->expmismatchI[type][si][sj];
00552 expintern = pf->expinternal;
00553
00554 if(type > 2)
00555 qe *= pf->expTermAU;
00556
00557 k = i + 1;
00558 if(S[k] == 3){
00559 if(k < j - VRNA_GQUAD_MIN_BOX_SIZE){
00560 minl = j - i + k - MAXLOOP - 2;
00561 u = k + VRNA_GQUAD_MIN_BOX_SIZE - 1;
00562 minl = MAX2(u, minl);
00563 u = j - 3;
00564 maxl = k + VRNA_GQUAD_MAX_BOX_SIZE + 1;
00565 maxl = MIN2(u, maxl);
00566 for(l = minl; l < maxl; l++){
00567 if(S[l] != 3) continue;
00568 q += qe * G[index[k]-l] * expintern[j - l - 1];
00569 }
00570 }
00571 }
00572
00573
00574 for(k = i + 2;
00575 k <= j - VRNA_GQUAD_MIN_BOX_SIZE;
00576 k++){
00577 u = k - i - 1;
00578 if(u > MAXLOOP) break;
00579 if(S[k] != 3) continue;
00580 minl = j - i + k - MAXLOOP - 2;
00581 r = k + VRNA_GQUAD_MIN_BOX_SIZE - 1;
00582 minl = MAX2(r, minl);
00583 maxl = k + VRNA_GQUAD_MAX_BOX_SIZE + 1;
00584 r = j - 1;
00585 maxl = MIN2(r, maxl);
00586 for(l = minl; l < maxl; l++){
00587 if(S[l] != 3) continue;
00588 q += qe * G[index[k]-l] * expintern[u + j - l - 1];
00589 }
00590 }
00591
00592 l = j - 1;
00593 if(S[l] == 3)
00594 for(k = i + 4; k < j - VRNA_GQUAD_MIN_BOX_SIZE; k++){
00595 u = k - i - 1;
00596 if(u>MAXLOOP) break;
00597 if(S[k] != 3) continue;
00598 q += qe * G[index[k]-l] * expintern[u];
00599 }
00600
00601 return q;
00602 }
00603
00604 #endif
