#ifdef HAVE_CONFIG_H #include "config.h" #endif #include #include #include #include #include #include "ViennaRNA/utils.h" #include "ViennaRNA/fold_vars.h" #include "ViennaRNA/energy_par.h" #include "ViennaRNA/constraints.h" #include "ViennaRNA/exterior_loops.h" #include "ViennaRNA/gquad.h" #include "ViennaRNA/structured_domains.h" #include "ViennaRNA/unstructured_domains.h" #include "ViennaRNA/multibranch_loops.h" struct default_data { int *idx; char *mx; int cp; int *hc_up; void *hc_dat; vrna_callback_hc_evaluate *hc_f; }; /* ################################# # PRIVATE FUNCTION DECLARATIONS # ################################# */ PRIVATE int E_mb_loop_fast(vrna_fold_compound_t *vc, int i, int j, int *dmli1, int *dmli2); PRIVATE int E_mb_loop_fast_comparative(vrna_fold_compound_t *vc, int i, int j, int *dmli1, int *dmli2); PRIVATE FLT_OR_DBL exp_E_mb_loop_fast(vrna_fold_compound_t *vc, int i, int j, FLT_OR_DBL *qqm1); PRIVATE FLT_OR_DBL exp_E_mb_loop_fast_comparative(vrna_fold_compound_t *vc, int i, int j, FLT_OR_DBL *qqm1); PRIVATE int E_ml_stems_fast(vrna_fold_compound_t *vc, int i, int j, int *fmi, int *dmli); PRIVATE int E_ml_stems_fast_comparative(vrna_fold_compound_t *vc, int i, int j, int *fmi, int *dmli); PRIVATE int extend_fm_3p(int i, int j, int *fm, vrna_fold_compound_t *vc); PRIVATE char hc_default(int i, int j, int k, int l, char d, void *data); PRIVATE char hc_default_user(int i, int j, int k, int l, char d, void *data); PRIVATE FLT_OR_DBL exp_E_ml_fast(vrna_fold_compound_t *vc, int i, int j, vrna_mx_pf_aux_ml_t *aux_mx); PRIVATE FLT_OR_DBL exp_E_ml_fast_comparative(vrna_fold_compound_t *vc, int i, int j, vrna_mx_pf_aux_ml_t *aux_mx); /* ################################# # BEGIN OF FUNCTION DEFINITIONS # ################################# */ PUBLIC int vrna_E_mb_loop_fast(vrna_fold_compound_t *vc, int i, int j, int *dmli1, int *dmli2) { int e = INF; if (vc) { switch (vc->type) { case VRNA_FC_TYPE_SINGLE: e = E_mb_loop_fast(vc, i, j, dmli1, dmli2); break; case VRNA_FC_TYPE_COMPARATIVE: e = E_mb_loop_fast_comparative(vc, i, j, dmli1, dmli2); break; } } return e; } PUBLIC int vrna_E_ml_stems_fast(vrna_fold_compound_t *vc, int i, int j, int *fmi, int *dmli) { int e = INF; if (vc) { switch (vc->type) { case VRNA_FC_TYPE_SINGLE: e = E_ml_stems_fast(vc, i, j, fmi, dmli); break; case VRNA_FC_TYPE_COMPARATIVE: e = E_ml_stems_fast_comparative(vc, i, j, fmi, dmli); break; } } return e; } PUBLIC FLT_OR_DBL vrna_exp_E_mb_loop_fast(vrna_fold_compound_t *vc, int i, int j, FLT_OR_DBL *qqm1) { FLT_OR_DBL q = 0.; if (vc) { switch (vc->type) { case VRNA_FC_TYPE_SINGLE: q = exp_E_mb_loop_fast(vc, i, j, qqm1); break; case VRNA_FC_TYPE_COMPARATIVE: q = exp_E_mb_loop_fast_comparative(vc, i, j, qqm1); break; } } return q; } PRIVATE int E_mb_loop_fast_comparative(vrna_fold_compound_t *vc, int i, int j, int *dmli1, int *dmli2) { short **S, **S5, **S3; int *indx, e, decomp, s, n_seq, dangle_model, *type; vrna_param_t *P; vrna_md_t *md; vrna_hc_t *hc; vrna_sc_t **scs; vrna_callback_hc_evaluate *evaluate; struct default_data hc_dat_local; n_seq = vc->n_seq; indx = vc->jindx; P = vc->params; md = &(P->model_details); hc = vc->hc; scs = vc->scs; dangle_model = md->dangles; e = INF; hc_dat_local.idx = vc->jindx; hc_dat_local.mx = hc->matrix; hc_dat_local.hc_up = hc->up_ml; hc_dat_local.cp = vc->cutpoint; if (hc->f) { evaluate = &hc_default_user; hc_dat_local.hc_f = hc->f; hc_dat_local.hc_dat = hc->data; } else { evaluate = &hc_default; } /* multi-loop decomposition ------------------------*/ if (evaluate(i, j, i + 1, j - 1, VRNA_DECOMP_PAIR_ML, &hc_dat_local)) { decomp = dmli1[j - 1]; type = (int *)vrna_alloc(n_seq * sizeof(int)); S = vc->S; S5 = vc->S5; /* S5[s][i] holds next base 5' of i in sequence s */ S3 = vc->S3; /* Sl[s][i] holds next base 3' of i in sequence s */ for (s = 0; s < n_seq; s++) { type[s] = md->pair[S[s][j]][S[s][i]]; if (type[s] == 0) type[s] = 7; } if (dangle_model) for (s = 0; s < n_seq; s++) decomp += E_MLstem(type[s], S5[s][j], S3[s][i], P); else for (s = 0; s < n_seq; s++) decomp += E_MLstem(type[s], -1, -1, P); if (scs) { for (s = 0; s < n_seq; s++) { if (scs[s]) if (scs[s]->energy_bp) decomp += scs[s]->energy_bp[indx[j] + i]; } } free(type); e = decomp + n_seq * P->MLclosing; } return e; } PRIVATE int E_mb_loop_fast(vrna_fold_compound_t *vc, int i, int j, int *dmli1, int *dmli2) { unsigned char type, tt; char *ptype; short S_i1, S_j1, *S; unsigned int *sn; int decomp, en, e, cp, *indx, *fc, ij, dangle_model, *rtype; vrna_hc_t *hc; vrna_sc_t *sc; vrna_param_t *P; vrna_callback_hc_evaluate *evaluate; struct default_data hc_dat_local; cp = vc->cutpoint; ptype = vc->ptype; S = vc->sequence_encoding; indx = vc->jindx; sn = vc->strand_number; hc = vc->hc; sc = vc->sc; fc = vc->matrices->fc; P = vc->params; ij = indx[j] + i; dangle_model = P->model_details.dangles; rtype = &(P->model_details.rtype[0]); type = (unsigned char)ptype[ij]; /* init values */ e = INF; decomp = INF; hc_dat_local.idx = vc->jindx; hc_dat_local.mx = hc->matrix; hc_dat_local.hc_up = hc->up_ml; hc_dat_local.cp = vc->cutpoint; if (hc->f) { evaluate = &hc_default_user; hc_dat_local.hc_f = hc->f; hc_dat_local.hc_dat = hc->data; } else { evaluate = &hc_default; } if (cp < 0) { S_i1 = S[i + 1]; S_j1 = S[j - 1]; } else { S_i1 = (sn[i] == sn[i + 1]) ? S[i + 1] : -1; S_j1 = (sn[j - 1] == sn[j]) ? S[j - 1] : -1; } if ((S_i1 >= 0) && (S_j1 >= 0)) { /* regular multi branch loop */ /* new closing pair (i,j) with mb part [i+1,j-1] */ if (evaluate(i, j, i + 1, j - 1, VRNA_DECOMP_PAIR_ML, &hc_dat_local)) { decomp = dmli1[j - 1]; tt = rtype[type]; if (tt == 0) tt = 7; if (decomp != INF) { switch (dangle_model) { /* no dangles */ case 0: decomp += E_MLstem(tt, -1, -1, P); if (sc) { if (sc->energy_bp) decomp += sc->energy_bp[ij]; if (sc->f) decomp += sc->f(i, j, i + 1, j - 1, VRNA_DECOMP_PAIR_ML, sc->data); } break; /* double dangles */ case 2: decomp += E_MLstem(tt, S_j1, S_i1, P); if (sc) { if (sc->energy_bp) decomp += sc->energy_bp[ij]; if (sc->f) decomp += sc->f(i, j, i + 1, j - 1, VRNA_DECOMP_PAIR_ML, sc->data); } break; /* normal dangles, aka dangles = 1 || 3 */ default: decomp += E_MLstem(tt, -1, -1, P); if (sc) { if (sc->energy_bp) decomp += sc->energy_bp[ij]; if (sc->f) decomp += sc->f(i, j, i + 1, j - 1, VRNA_DECOMP_PAIR_ML, sc->data); } break; } } } if (dangle_model % 2) { /* dangles == 1 || dangles == 3 */ /* new closing pair (i,j) with mb part [i+2,j-1] */ if (evaluate(i, j, i + 2, j - 1, VRNA_DECOMP_PAIR_ML, &hc_dat_local)) { if (dmli2[j - 1] != INF) { tt = rtype[type]; if (tt == 0) tt = 7; en = dmli2[j - 1] + E_MLstem(tt, -1, S_i1, P) + P->MLbase; if (sc) { if (sc->energy_up) en += sc->energy_up[i + 1][1]; if (sc->energy_bp) en += sc->energy_bp[ij]; if (sc->f) en += sc->f(i, j, i + 2, j - 1, VRNA_DECOMP_PAIR_ML, sc->data); } decomp = MIN2(decomp, en); } } /* new closing pair (i,j) with mb part [i+2.j-2] */ if (evaluate(i, j, i + 2, j - 2, VRNA_DECOMP_PAIR_ML, &hc_dat_local)) { if (dmli2[j - 2] != INF) { tt = rtype[type]; if (tt == 0) tt = 7; en = dmli2[j - 2] + E_MLstem(tt, S_j1, S_i1, P) + 2 * P->MLbase; if (sc) { if (sc->energy_up) en += sc->energy_up[i + 1][1] + sc->energy_up[j - 1][1]; if (sc->energy_bp) en += sc->energy_bp[ij]; if (sc->f) en += sc->f(i, j, i + 2, j - 2, VRNA_DECOMP_PAIR_ML, sc->data); } decomp = MIN2(decomp, en); } } /* new closing pair (i,j) with mb part [i+1, j-2] */ if (evaluate(i, j, i + 1, j - 2, VRNA_DECOMP_PAIR_ML, &hc_dat_local)) { if (dmli1[j - 2] != INF) { tt = rtype[type]; if (tt == 0) tt = 7; en = dmli1[j - 2] + E_MLstem(tt, S_j1, -1, P) + P->MLbase; if (sc) { if (sc->energy_up) en += sc->energy_up[j - 1][1]; if (sc->energy_bp) en += sc->energy_bp[ij]; if (sc->f) en += sc->f(i, j, i + 1, j - 2, VRNA_DECOMP_PAIR_ML, sc->data); } decomp = MIN2(decomp, en); } } } /* end if dangles % 2 */ if (decomp != INF) e = decomp + P->MLclosing; } /* end regular multibranch loop */ if (sn[i] != sn[j]) { /* multibrach like cofold structure with cut somewhere between i and j */ if (evaluate(i, j, i + 1, j - 1, VRNA_DECOMP_PAIR_ML, &hc_dat_local)) { if ((fc[i + 1] != INF) && (fc[j - 1] != INF)) { decomp = fc[i + 1] + fc[j - 1]; tt = rtype[type]; if (tt == 0) tt = 7; switch (dangle_model) { case 0: decomp += E_ExtLoop(tt, -1, -1, P); break; case 2: decomp += E_ExtLoop(tt, S_j1, S_i1, P); break; default: decomp += E_ExtLoop(tt, -1, -1, P); break; } } } if (dangle_model % 2) { /* dangles == 1 || dangles == 3 */ if (evaluate(i, j, i + 2, j - 1, VRNA_DECOMP_PAIR_ML, &hc_dat_local)) { if ((fc[i + 2] != INF) && (fc[j - 1] != INF)) { tt = rtype[type]; if (tt == 0) tt = 7; en = fc[i + 2] + fc[j - 1] + E_ExtLoop(tt, -1, S_i1, P); decomp = MIN2(decomp, en); } } if (evaluate(i, j, i + 1, j - 2, VRNA_DECOMP_PAIR_ML, &hc_dat_local)) { if ((fc[i + 1] != INF) && (fc[j - 2] != INF)) { tt = rtype[type]; if (tt == 0) tt = 7; en = fc[i + 1] + fc[j - 2] + E_ExtLoop(tt, S_j1, -1, P); decomp = MIN2(decomp, en); } } if (evaluate(i, j, i + 2, j - 2, VRNA_DECOMP_PAIR_ML, &hc_dat_local)) { if ((fc[i + 2] != INF) && (fc[j - 2] != INF)) { tt = rtype[type]; if (tt == 0) tt = 7; en = fc[i + 2] + fc[j - 2] + E_ExtLoop(tt, S_j1, S_i1, P); decomp = MIN2(decomp, en); } } } e = MIN2(e, decomp); } return e; } PUBLIC int E_mb_loop_stack(int i, int j, vrna_fold_compound_t *vc) { unsigned char type, type_2; char *ptype; int e, decomp, en, i1k, k1j1, ij, k, *indx, *c, *fML, turn, *rtype; vrna_param_t *P; vrna_md_t *md; vrna_hc_t *hc; vrna_sc_t *sc; vrna_callback_hc_evaluate *evaluate; struct default_data hc_dat_local; indx = vc->jindx; hc = vc->hc; c = vc->matrices->c; fML = vc->matrices->fML; P = vc->params; md = &(P->model_details); turn = md->min_loop_size; ptype = vc->ptype; rtype = &(md->rtype[0]); sc = vc->sc; e = INF; ij = indx[j] + i; type = ptype[ij]; hc_dat_local.idx = vc->jindx; hc_dat_local.mx = hc->matrix; hc_dat_local.hc_up = hc->up_ml; hc_dat_local.cp = vc->cutpoint; if (hc->f) { evaluate = &hc_default_user; hc_dat_local.hc_f = hc->f; hc_dat_local.hc_dat = hc->data; } else { evaluate = &hc_default; } if (evaluate(i, j, i + 1, j - 1, VRNA_DECOMP_PAIR_ML, &hc_dat_local)) { if (type == 0) type = 7; decomp = INF; k1j1 = indx[j - 1] + i + 2 + turn + 1; for (k = i + 2 + turn; k < j - 2 - turn; k++, k1j1++) { i1k = indx[k] + i + 1; if (evaluate(i, j, i + 1, k, VRNA_DECOMP_ML_COAXIAL, &hc_dat_local)) { type_2 = rtype[(unsigned char)ptype[i1k]]; if (type_2 == 0) type_2 = 7; en = c[i1k] + P->stack[type][type_2] + fML[k1j1]; if (sc) if (sc->f) en += sc->f(i, j, i + 1, k, VRNA_DECOMP_ML_COAXIAL, sc->data); decomp = MIN2(decomp, en); } if (evaluate(i, j, k + 1, j - 1, VRNA_DECOMP_ML_COAXIAL, &hc_dat_local)) { type_2 = rtype[(unsigned char)ptype[k1j1]]; if (type_2 == 0) type_2 = 7; en = c[k1j1] + P->stack[type][type_2] + fML[i1k]; if (sc) if (sc->f) en += sc->f(i, j, k + 1, j - 1, VRNA_DECOMP_ML_COAXIAL, sc->data); decomp = MIN2(decomp, en); } } /* no TermAU penalty if coax stack */ decomp += 2 * P->MLintern[1] + P->MLclosing; if (sc) { if (sc->energy_bp) decomp += sc->energy_bp[ij]; if (sc->f) decomp += sc->f(i, j, i + 1, j - 1, VRNA_DECOMP_PAIR_ML, sc->data); } e = decomp; } return e; } PUBLIC int E_ml_rightmost_stem(int i, int j, vrna_fold_compound_t *vc) { if ((vc) && (vc->matrices) && (vc->matrices->fM1)) return extend_fm_3p(i, j, vc->matrices->fM1, vc); return INF; } /* * compose a multibranch loop part fm[i:j] * by either c[i,j]/ggg[i,j] or fm[i:j-1] * * This function can be used for fM and fM1 */ PRIVATE int extend_fm_3p(int i, int j, int *fm, vrna_fold_compound_t *vc) { short *S; unsigned int *sn; int en, length, *indx, *c, *ggg, ij, type, dangle_model, with_gquad, e, u, k, cnt, with_ud; vrna_param_t *P; vrna_hc_t *hc; vrna_sc_t *sc; vrna_ud_t *domains_up; vrna_callback_hc_evaluate *evaluate; struct default_data hc_dat_local; P = vc->params; length = vc->length; S = vc->sequence_encoding; indx = vc->jindx; sn = vc->strand_number; hc = vc->hc; sc = vc->sc; c = vc->matrices->c; ggg = vc->matrices->ggg; ij = indx[j] + i; type = vc->ptype[ij]; dangle_model = P->model_details.dangles; with_gquad = P->model_details.gquad; domains_up = vc->domains_up; with_ud = (domains_up && domains_up->energy_cb) ? 1 : 0; e = INF; hc_dat_local.idx = vc->jindx; hc_dat_local.mx = hc->matrix; hc_dat_local.hc_up = hc->up_ml; hc_dat_local.cp = vc->cutpoint; if (hc->f) { evaluate = &hc_default_user; hc_dat_local.hc_f = hc->f; hc_dat_local.hc_dat = hc->data; } else { evaluate = &hc_default; } if (sn[i - 1] == sn[i]) { if (sn[j] == sn[j + 1]) { if (evaluate(i, j, i, j, VRNA_DECOMP_ML_STEM, &hc_dat_local)) { if (type == 0) type = 7; e = c[ij]; if (e != INF) { switch (dangle_model) { case 2: e += E_MLstem(type, (i == 1) ? S[length] : S[i - 1], S[j + 1], P); break; default: e += E_MLstem(type, -1, -1, P); break; } if (sc) if (sc->f) e += sc->f(i, j, i, j, VRNA_DECOMP_ML_STEM, sc->data); } } if (with_gquad) { if (sn[i] == sn[j]) { en = ggg[ij] + E_MLstem(0, -1, -1, P); e = MIN2(e, en); } } } if (sn[j - 1] == sn[j]) { if (evaluate(i, j, i, j - 1, VRNA_DECOMP_ML_ML, &hc_dat_local)) { if (fm[indx[j - 1] + i] != INF) { en = fm[indx[j - 1] + i] + P->MLbase; if (sc) { if (sc->energy_up) en += sc->energy_up[j][1]; if (sc->f) en += sc->f(i, j, i, j - 1, VRNA_DECOMP_ML_ML, sc->data); } e = MIN2(e, en); } } } if (with_ud) { for (cnt = 0; cnt < domains_up->uniq_motif_count; cnt++) { u = domains_up->uniq_motif_size[cnt]; k = j - u + 1; if ((k > i) && (sn[j - u] == sn[j])) { if (evaluate(i, j, i, k - 1, VRNA_DECOMP_ML_ML, &hc_dat_local)) { if (fm[indx[k - 1] + i] != INF) { en = domains_up->energy_cb(vc, k, j, VRNA_UNSTRUCTURED_DOMAIN_MB_LOOP | VRNA_UNSTRUCTURED_DOMAIN_MOTIF, domains_up->data); if (en != INF) { en += fm[indx[k - 1] + i] + u * P->MLbase; if (sc) { if (sc->energy_up) en += sc->energy_up[k][u]; if (sc->f) en += sc->f(i, j, i, k - 1, VRNA_DECOMP_ML_ML, sc->data); } e = MIN2(e, en); } } } } } } } return e; } PRIVATE int E_ml_stems_fast(vrna_fold_compound_t *vc, int i, int j, int *fmi, int *dmli) { char *ptype; short *S; unsigned int *sn; int k, en, decomp, mm5, mm3, type_2, k1j, stop, length, *indx, *c, *fm, ij, dangle_model, turn, type, *rtype, circular, cp, e, u, cnt, with_ud; vrna_hc_t *hc; vrna_sc_t *sc; vrna_param_t *P; vrna_ud_t *domains_up; vrna_callback_hc_evaluate *evaluate; struct default_data hc_dat_local; length = (int)vc->length; ptype = vc->ptype; S = vc->sequence_encoding; indx = vc->jindx; sn = vc->strand_number; hc = vc->hc; sc = vc->sc; c = vc->matrices->c; fm = vc->matrices->fML; P = vc->params; ij = indx[j] + i; dangle_model = P->model_details.dangles; turn = P->model_details.min_loop_size; type = ptype[ij]; rtype = &(P->model_details.rtype[0]); circular = P->model_details.circ; cp = vc->cutpoint; domains_up = vc->domains_up; with_ud = (domains_up && domains_up->energy_cb) ? 1 : 0; e = INF; hc_dat_local.idx = vc->jindx; hc_dat_local.mx = hc->matrix; hc_dat_local.hc_up = hc->up_ml; hc_dat_local.cp = vc->cutpoint; if (hc->f) { evaluate = &hc_default_user; hc_dat_local.hc_f = hc->f; hc_dat_local.hc_dat = hc->data; } else { evaluate = &hc_default; } /* * extension with one unpaired nucleotide at the right (3' site) * or full branch of (i,j) */ e = extend_fm_3p(i, j, fm, vc); /* * extension with one unpaired nucleotide at 5' site * and all other variants which are needed for odd * dangle models */ if (sn[i - 1] == sn[i]) { if (sn[i] == sn[i + 1]) { if (evaluate(i, j, i + 1, j, VRNA_DECOMP_ML_ML, &hc_dat_local)) { if (fm[ij + 1] != INF) { en = fm[ij + 1] + P->MLbase; if (sc) { if (sc->energy_up) en += sc->energy_up[i][1]; if (sc->f) en += sc->f(i, j, i + 1, j, VRNA_DECOMP_ML_ML, sc->data); } e = MIN2(e, en); } } } /* extension with bound ligand on 5'site */ if (with_ud) { for (cnt = 0; cnt < domains_up->uniq_motif_count; cnt++) { u = domains_up->uniq_motif_size[cnt]; k = i + u - 1; if ((k < j) && (sn[i] == sn[k + 1])) { if (evaluate(i, j, k + 1, j, VRNA_DECOMP_ML_ML, &hc_dat_local)) { if (fm[ij + u] != INF) { en = domains_up->energy_cb(vc, i, k, VRNA_UNSTRUCTURED_DOMAIN_MB_LOOP | VRNA_UNSTRUCTURED_DOMAIN_MOTIF, domains_up->data); if (en != INF) { en += fm[ij + u] + u * P->MLbase; if (sc) { if (sc->energy_up) en += sc->energy_up[i][u]; if (sc->f) en += sc->f(i, j, k + 1, j, VRNA_DECOMP_ML_ML, sc->data); } e = MIN2(e, en); } } } } } } if (dangle_model % 2) { /* dangle_model = 1 || 3 */ mm5 = ((i > 1) || circular) ? S[i] : -1; mm3 = ((j < length) || circular) ? S[j] : -1; if (sn[i] == sn[i + 1]) { if (evaluate(i, j, i + 1, j, VRNA_DECOMP_ML_STEM, &hc_dat_local)) { if (c[ij + 1] != INF) { type = ptype[ij + 1]; if (type == 0) type = 7; en = c[ij + 1] + E_MLstem(type, mm5, -1, P) + P->MLbase; if (sc) { if (sc->energy_up) en += sc->energy_up[i][1]; if (sc->f) en += sc->f(i, j, i + 1, j, VRNA_DECOMP_ML_STEM, sc->data); } e = MIN2(e, en); } } } if (sn[j - 1] == sn[j]) { if (evaluate(i, j, i, j - 1, VRNA_DECOMP_ML_STEM, &hc_dat_local)) { if (c[indx[j - 1] + i] != INF) { type = ptype[indx[j - 1] + i]; if (type == 0) type = 7; en = c[indx[j - 1] + i] + E_MLstem(type, -1, mm3, P) + P->MLbase; if (sc) { if (sc->energy_up) en += sc->energy_up[j][1]; if (sc->f) en += sc->f(i, j, i, j - 1, VRNA_DECOMP_ML_STEM, sc->data); } e = MIN2(e, en); } } } if ((sn[j - 1] == sn[j]) && (sn[i] == sn[i + 1])) { if (evaluate(i, j, i + 1, j - 1, VRNA_DECOMP_ML_STEM, &hc_dat_local)) { if (c[indx[j - 1] + i + 1] != INF) { type = ptype[indx[j - 1] + i + 1]; if (type == 0) type = 7; en = c[indx[j - 1] + i + 1] + E_MLstem(type, mm5, mm3, P) + 2 * P->MLbase; if (sc) { if (sc->energy_up) en += sc->energy_up[j][1] + sc->energy_up[i][1]; if (sc->f) en += sc->f(i, j, i + 1, j - 1, VRNA_DECOMP_ML_STEM, sc->data); } e = MIN2(e, en); } } } } /* end special cases for dangles == 1 || dangles == 3 */ } /* modular decomposition -------------------------------*/ k1j = indx[j] + i + turn + 2; stop = (cp > 0) ? (cp - 1) : (j - 2 - turn); /* duplicated code is faster than conditions in loop */ if (hc->f) { if (sc && sc->f) { for (decomp = INF, k = i + 1 + turn; k <= stop; k++, k1j++) { if ((fmi[k] != INF) && (fm[k1j] != INF) && hc->f(i, j, k, k + 1, VRNA_DECOMP_ML_ML_ML, &hc_dat_local)) { en = fmi[k] + fm[k1j]; en += sc->f(i, j, k, k + 1, VRNA_DECOMP_ML_ML_ML, sc->data); decomp = MIN2(decomp, en); } } k++; k1j++; for (; k <= j - 2 - turn; k++, k1j++) { if ((fmi[k] != INF) && (fm[k1j] != INF) && hc->f(i, j, k, k + 1, VRNA_DECOMP_ML_ML_ML, &hc_dat_local)) { en = fmi[k] + fm[k1j]; en += sc->f(i, j, k, k + 1, VRNA_DECOMP_ML_ML_ML, sc->data); decomp = MIN2(decomp, en); } } } else { for (decomp = INF, k = i + 1 + turn; k <= stop; k++, k1j++) { if ((fmi[k] != INF) && (fm[k1j] != INF) && hc->f(i, j, k, k + 1, VRNA_DECOMP_ML_ML_ML, &hc_dat_local)) { en = fmi[k] + fm[k1j]; decomp = MIN2(decomp, en); } } k++; k1j++; for (; k <= j - 2 - turn; k++, k1j++) { if ((fmi[k] != INF) && (fm[k1j] != INF) && hc->f(i, j, k, k + 1, VRNA_DECOMP_ML_ML_ML, &hc_dat_local)) { en = fmi[k] + fm[k1j]; decomp = MIN2(decomp, en); } } } } else { if (sc && sc->f) { for (decomp = INF, k = i + 1 + turn; k <= stop; k++, k1j++) { if ((fmi[k] != INF) && (fm[k1j] != INF)) { en = fmi[k] + fm[k1j]; en += sc->f(i, j, k, k + 1, VRNA_DECOMP_ML_ML_ML, sc->data); decomp = MIN2(decomp, en); } } k++; k1j++; for (; k <= j - 2 - turn; k++, k1j++) { if ((fmi[k] != INF) && (fm[k1j] != INF)) { en = fmi[k] + fm[k1j]; en += sc->f(i, j, k, k + 1, VRNA_DECOMP_ML_ML_ML, sc->data); decomp = MIN2(decomp, en); } } } else { for (decomp = INF, k = i + 1 + turn; k <= stop; k++, k1j++) { if ((fmi[k] != INF) && (fm[k1j] != INF)) { en = fmi[k] + fm[k1j]; decomp = MIN2(decomp, en); } } k++; k1j++; for (; k <= j - 2 - turn; k++, k1j++) { if ((fmi[k] != INF) && (fm[k1j] != INF)) { en = fmi[k] + fm[k1j]; decomp = MIN2(decomp, en); } } } } dmli[j] = decomp; /* store for use in fast ML decompositon */ e = MIN2(e, decomp); /* coaxial stacking */ if (dangle_model == 3) { /* additional ML decomposition as two coaxially stacked helices */ int ik; k1j = indx[j] + i + turn + 2; for (decomp = INF, k = i + 1 + turn; k <= stop; k++, k1j++) { ik = indx[k] + i; if (evaluate(i, k, k + 1, j, VRNA_DECOMP_ML_COAXIAL_ENC, &hc_dat_local)) { type = rtype[(unsigned char)ptype[ik]]; type_2 = rtype[(unsigned char)ptype[k1j]]; if (type == 0) type = 7; if (type_2 == 0) type_2 = 7; en = c[ik] + c[k1j] + P->stack[type][type_2]; if (sc) if (sc->f) en += sc->f(i, k, k + 1, j, VRNA_DECOMP_ML_COAXIAL_ENC, sc->data); decomp = MIN2(decomp, en); } } k++; k1j++; for (; k <= j - 2 - turn; k++, k1j++) { ik = indx[k] + i; if (evaluate(i, k, k + 1, j, VRNA_DECOMP_ML_COAXIAL_ENC, &hc_dat_local)) { type = rtype[(unsigned char)ptype[ik]]; type_2 = rtype[(unsigned char)ptype[k1j]]; if (type == 0) type = 7; if (type_2 == 0) type_2 = 7; en = c[ik] + c[k1j] + P->stack[type][type_2]; if (sc) if (sc->f) en += sc->f(i, k, k + 1, j, VRNA_DECOMP_ML_COAXIAL, sc->data); decomp = MIN2(decomp, en); } } decomp += 2 * P->MLintern[1]; /* no TermAU penalty if coax stack */ #if 0 /* * This is needed for Y shaped ML loops with coax stacking of * interior pairts, but backtracking will fail if activated */ DMLi[j] = MIN2(DMLi[j], decomp); DMLi[j] = MIN2(DMLi[j], DMLi[j - 1] + P->MLbase); DMLi[j] = MIN2(DMLi[j], DMLi1[j] + P->MLbase); new_fML = MIN2(new_fML, DMLi[j]); #endif e = MIN2(e, decomp); } fmi[j] = e; return e; } PRIVATE int E_ml_stems_fast_comparative(vrna_fold_compound_t *vc, int i, int j, int *fmi, int *dmli) { char *hard_constraints; short **S, **S5, **S3; unsigned short **a2s; int e, energy, *c, *fML, *ggg, ij, *indx, s, n_seq, k, dangle_model, decomp, turn, *type; vrna_param_t *P; vrna_md_t *md; vrna_mx_mfe_t *matrices; vrna_hc_t *hc; vrna_sc_t **scs; n_seq = vc->n_seq; matrices = vc->matrices; P = vc->params; md = &(P->model_details); c = matrices->c; fML = matrices->fML; ggg = matrices->ggg; indx = vc->jindx; hc = vc->hc; scs = vc->scs; hard_constraints = hc->matrix; dangle_model = md->dangles; turn = md->min_loop_size; a2s = vc->a2s; ij = indx[j] + i; e = INF; if (hc->up_ml[i]) { energy = fML[ij + 1] + n_seq * P->MLbase; if (scs) { for (s = 0; s < n_seq; s++) { if (scs[s]) if (scs[s]->energy_up) energy += scs[s]->energy_up[a2s[s][i]][1]; } } e = MIN2(e, energy); } if (hc->up_ml[j]) { energy = fML[indx[j - 1] + i] + n_seq * P->MLbase; if (scs) { for (s = 0; s < n_seq; s++) { if (scs[s]) if (scs[s]->energy_up) energy += scs[s]->energy_up[a2s[s][j]][1]; } } e = MIN2(e, energy); } if (hard_constraints[ij] & VRNA_CONSTRAINT_CONTEXT_MB_LOOP_ENC) { energy = c[ij]; type = (int *)vrna_alloc(n_seq * sizeof(int)); S = vc->S; S5 = vc->S5; /* S5[s][i] holds next base 5' of i in sequence s */ S3 = vc->S3; /* Sl[s][i] holds next base 3' of i in sequence s */ for (s = 0; s < n_seq; s++) { type[s] = md->pair[S[s][i]][S[s][j]]; if (type[s] == 0) type[s] = 7; } if (dangle_model) for (s = 0; s < n_seq; s++) energy += E_MLstem(type[s], S5[s][i], S3[s][j], P); else for (s = 0; s < n_seq; s++) energy += E_MLstem(type[s], -1, -1, P); e = MIN2(e, energy); if (md->gquad) { decomp = ggg[indx[j] + i] + n_seq * E_MLstem(0, -1, -1, P); e = MIN2(e, decomp); } free(type); } /* modular decomposition -------------------------------*/ for (decomp = INF, k = i + 1 + turn; k <= j - 2 - turn; k++) decomp = MIN2(decomp, fmi[k] + fML[indx[j] + k + 1]); dmli[j] = decomp; /* store for later use in ML decompositon */ e = MIN2(e, decomp); fmi[j] = e; /* store for later use in ML decompositon */ return e; } PRIVATE FLT_OR_DBL exp_E_mb_loop_fast(vrna_fold_compound_t *vc, int i, int j, FLT_OR_DBL *qqm1) { unsigned char type, tt; char *ptype; short *S1; unsigned int *sn; int ij, k, kl, *my_iindx, *jindx, *rtype; FLT_OR_DBL qbt1, temp, qqqmmm, *qm, *scale, expMLclosing; vrna_hc_t *hc; vrna_sc_t *sc; vrna_exp_param_t *pf_params; vrna_md_t *md; vrna_callback_hc_evaluate *evaluate; struct default_data hc_dat_local; my_iindx = vc->iindx; jindx = vc->jindx; sc = vc->sc; ptype = vc->ptype; S1 = vc->sequence_encoding; qm = vc->exp_matrices->qm; scale = vc->exp_matrices->scale; pf_params = vc->exp_params; md = &(pf_params->model_details); ij = jindx[j] + i; sn = vc->strand_number; hc = vc->hc; expMLclosing = pf_params->expMLclosing; qbt1 = 0.; hc_dat_local.idx = vc->jindx; hc_dat_local.mx = hc->matrix; hc_dat_local.hc_up = hc->up_ml; hc_dat_local.cp = vc->cutpoint; if (hc->f) { evaluate = &hc_default_user; hc_dat_local.hc_f = hc->f; hc_dat_local.hc_dat = hc->data; } else { evaluate = &hc_default; } /* multiple stem loop contribution */ if (evaluate(i, j, i + 1, j - 1, VRNA_DECOMP_PAIR_ML, &hc_dat_local) && (sn[i] == sn[i + 1]) && (sn[j - 1] == sn[j])) { type = (unsigned char)ptype[ij]; rtype = &(md->rtype[0]); tt = rtype[type]; if (tt == 0) tt = 7; qqqmmm = expMLclosing * exp_E_MLstem(tt, S1[j - 1], S1[i + 1], pf_params) * scale[2]; temp = 0.0; kl = my_iindx[i + 1] - (i + 1); if (sc) { if (sc->exp_energy_bp) qqqmmm *= sc->exp_energy_bp[my_iindx[i] - j]; if (sc->exp_f) { qqqmmm *= sc->exp_f(i, j, i, j, VRNA_DECOMP_PAIR_ML, sc->data); for (k = i + 2; k <= j - 1; k++, kl--) { if (sn[k - 1] == sn[k]) { temp += qm[kl] * qqm1[k] * sc->exp_f(i + 1, j - 1, k - 1, k, VRNA_DECOMP_ML_ML_ML, sc->data); } } } else { for (k = i + 2; k <= j - 1; k++, kl--) { if (sn[k - 1] == sn[k]) temp += qm[kl] * qqm1[k]; } } } else { for (k = i + 2; k <= j - 1; k++, kl--) { if (sn[k - 1] == sn[k]) temp += qm[kl] * qqm1[k]; } } qbt1 += temp * qqqmmm; } return qbt1; } PRIVATE FLT_OR_DBL exp_E_mb_loop_fast_comparative(vrna_fold_compound_t *vc, int i, int j, FLT_OR_DBL *qqm1) { short **S, **S5, **S3; int k, kl, *my_iindx, *types, n_seq, s; FLT_OR_DBL qbt1, temp, qqqmmm, *qm, *scale, expMLclosing; vrna_hc_t *hc; vrna_sc_t **scs; vrna_exp_param_t *pf_params; vrna_md_t *md; vrna_callback_hc_evaluate *evaluate; struct default_data hc_dat_local; my_iindx = vc->iindx; qm = vc->exp_matrices->qm; scale = vc->exp_matrices->scale; pf_params = vc->exp_params; md = &(pf_params->model_details); hc = vc->hc; expMLclosing = pf_params->expMLclosing; qbt1 = 0.; types = NULL; hc_dat_local.idx = vc->jindx; hc_dat_local.mx = hc->matrix; hc_dat_local.hc_up = hc->up_ml; hc_dat_local.cp = vc->cutpoint; if (hc->f) { evaluate = &hc_default_user; hc_dat_local.hc_f = hc->f; hc_dat_local.hc_dat = hc->data; } else { evaluate = &hc_default; } /* multiple stem loop contribution */ if (evaluate(i, j, i + 1, j - 1, VRNA_DECOMP_PAIR_ML, &hc_dat_local)) { S = vc->S; S5 = vc->S5; /* S5[s][i] holds next base 5' of i in sequence s */ S3 = vc->S3; /* Sl[s][i] holds next base 3' of i in sequence s */ scs = vc->scs; n_seq = vc->n_seq; types = (int *)vrna_alloc(sizeof(int) * n_seq); qqqmmm = 1.; for (s = 0; s < n_seq; s++) { types[s] = md->pair[S[s][j]][S[s][i]]; if (types[s] == 0) types[s] = 7; } for (s = 0; s < n_seq; s++) qqqmmm *= exp_E_MLstem(types[s], S5[s][j], S3[s][i], pf_params) * expMLclosing; if (scs) { for (s = 0; s < n_seq; s++) { if (scs[s]) if (scs[s]->exp_energy_bp) qqqmmm *= scs[s]->exp_energy_bp[my_iindx[i] - j]; } } /* multi-loop loop contribution */ temp = 0.; kl = my_iindx[i + 1] - (i + 1); for (k = i + 2; k <= j - 1; k++, kl--) temp += qm[kl] * qqm1[k]; temp *= scale[2]; qbt1 = temp * qqqmmm; } /* cleanup */ free(types); return qbt1; } /* ################################# # Backtracking functions below # ################################# */ PUBLIC int vrna_BT_mb_loop_fake(vrna_fold_compound_t *vc, int *u, int *i, int *j, vrna_bp_stack_t *bp_stack, int *stack_count) { unsigned char type; char *ptype; short mm5, mm3, *S1; unsigned int *sn; int length, ii, jj, k, en, cp, fij, fi, *my_c, *my_fc, *my_ggg, *idx, with_gquad, dangle_model, turn; vrna_param_t *P; vrna_md_t *md; vrna_hc_t *hc; vrna_sc_t *sc; cp = vc->cutpoint; length = vc->length; P = vc->params; md = &(P->model_details); sn = vc->strand_number; hc = vc->hc; sc = vc->sc; S1 = vc->sequence_encoding; ptype = vc->ptype; idx = vc->jindx; my_c = vc->matrices->c; my_fc = vc->matrices->fc; my_ggg = vc->matrices->ggg; turn = md->min_loop_size; with_gquad = md->gquad; dangle_model = md->dangles; ii = *i; jj = *j; if (ii < cp) { /* 'lower' part (fc[iup_ext[ii]) ? my_fc[ii + 1] : INF; if (sc) if (sc->energy_up) fi += sc->energy_up[ii][1]; if (++ii == jj) break; } while (fij == fi); ii--; if (jj < ii + turn + 2) { /* no more pairs */ *u = *i = *j = -1; return 1; } mm5 = (ii > 1 && (sn[ii - 1] == sn[ii])) ? S1[ii - 1] : -1; /* i or i+1 is paired. Find pairing partner */ switch (dangle_model) { case 0: for (k = ii + turn + 1; k <= jj; k++) { if (hc->matrix[idx[k] + ii] & VRNA_CONSTRAINT_CONTEXT_EXT_LOOP) { type = (unsigned char)ptype[idx[k] + ii]; if (type == 0) type = 7; if (fij == my_fc[k + 1] + my_c[idx[k] + ii] + E_ExtLoop(type, -1, -1, P)) { bp_stack[++(*stack_count)].i = ii; bp_stack[(*stack_count)].j = k; *u = k + 1; *i = ii; *j = k; return 1; } } if (with_gquad) { if (fij == my_fc[k + 1] + my_ggg[idx[k] + ii]) { *u = k + 1; *i = *j = -1; vrna_BT_gquad_mfe(vc, ii, k, bp_stack, stack_count); return 1; } } } break; case 2: for (k = ii + turn + 1; k <= jj; k++) { if (hc->matrix[idx[k] + ii] & VRNA_CONSTRAINT_CONTEXT_EXT_LOOP) { mm3 = (sn[k] == sn[k + 1]) ? S1[k + 1] : -1; type = (unsigned char)ptype[idx[k] + ii]; if (type == 0) type = 7; if (fij == my_fc[k + 1] + my_c[idx[k] + ii] + E_ExtLoop(type, mm5, mm3, P)) { bp_stack[++(*stack_count)].i = ii; bp_stack[(*stack_count)].j = k; *u = k + 1; *i = ii; *j = k; return 1; } } if (with_gquad) { if (fij == my_fc[k + 1] + my_ggg[idx[k] + ii]) { *u = k + 1; *i = *j = -1; vrna_BT_gquad_mfe(vc, ii, k, bp_stack, stack_count); return 1; } } } break; default: for (k = ii + turn + 1; k <= jj; k++) { if (hc->matrix[idx[k] + ii] & VRNA_CONSTRAINT_CONTEXT_EXT_LOOP) { type = (unsigned char)ptype[idx[k] + ii]; if (type == 0) type = 7; if (fij == my_fc[k + 1] + my_c[idx[k] + ii] + E_ExtLoop(type, -1, -1, P)) { bp_stack[++(*stack_count)].i = ii; bp_stack[(*stack_count)].j = k; *u = k + 1; *i = ii; *j = k; return 1; } if (hc->up_ext[k + 1]) { mm3 = (sn[k] == sn[k + 1]) ? S1[k + 1] : -1; en = my_c[idx[k] + ii]; if (sc) if (sc->energy_up) en += sc->energy_up[k + 1][1]; if (fij == my_fc[k + 2] + en + E_ExtLoop(type, -1, mm3, P)) { bp_stack[++(*stack_count)].i = ii; bp_stack[(*stack_count)].j = k; *u = k + 2; *i = ii; *j = k; return 1; } } } if (with_gquad) { if (fij == my_fc[k + 1] + my_ggg[idx[k] + ii]) { *u = k + 1; *i = *j = -1; vrna_BT_gquad_mfe(vc, ii, k, bp_stack, stack_count); return 1; } } if (hc->matrix[idx[k] + ii + 1] & VRNA_CONSTRAINT_CONTEXT_EXT_LOOP) { if (hc->up_ext[ii]) { mm5 = (sn[ii] == sn[ii + 1]) ? S1[ii] : -1; mm3 = (sn[k] == sn[k + 1]) ? S1[k + 1] : -1; type = ptype[idx[k] + ii + 1]; if (type == 0) type = 7; en = my_c[idx[k] + ii + 1]; if (sc) if (sc->energy_up) en += sc->energy_up[ii][1]; if (fij == en + my_fc[k + 1] + E_ExtLoop(type, mm5, -1, P)) { bp_stack[++(*stack_count)].i = ii + 1; bp_stack[(*stack_count)].j = k; *u = k + 1; *i = ii + 1; *j = k; return 1; } if (k < jj) { if (hc->up_ext[k + 1]) { if (sc) if (sc->energy_up) en += sc->energy_up[k + 1][1]; if (fij == en + my_fc[k + 2] + E_ExtLoop(type, mm5, mm3, P)) { bp_stack[++(*stack_count)].i = ii + 1; bp_stack[(*stack_count)].j = k; *u = k + 2; *i = ii + 1; *j = k; return 1; } } } } } } break; } } else { /* 'upper' part (fc[i=cut,j>cut]) */ /* nibble off unpaired 3' bases */ do { fij = my_fc[jj]; fi = (hc->up_ext[jj]) ? my_fc[jj - 1] : INF; if (sc) if (sc->energy_up) fi += sc->energy_up[jj][1]; if (--jj == ii) break; } while (fij == fi); jj++; if (jj < ii + turn + 2) { /* no more pairs */ *u = *i = *j = -1; return 1; } /* j or j-1 is paired. Find pairing partner */ mm3 = ((jj < length) && (sn[jj] == sn[jj + 1])) ? S1[jj + 1] : -1; switch (dangle_model) { case 0: for (k = jj - turn - 1; k >= ii; k--) { if (with_gquad) { if (fij == my_fc[k - 1] + my_ggg[idx[jj] + k]) { *u = k - 1; *i = *j = -1; vrna_BT_gquad_mfe(vc, k, jj, bp_stack, stack_count); return 1; } } if (hc->matrix[idx[jj] + k] & VRNA_CONSTRAINT_CONTEXT_EXT_LOOP) { type = (unsigned char)ptype[idx[jj] + k]; if (type == 0) type = 7; en = my_c[idx[jj] + k]; if (sn[k] != sn[jj]) en += P->DuplexInit; if (fij == my_fc[k - 1] + en + E_ExtLoop(type, -1, -1, P)) { bp_stack[++(*stack_count)].i = k; bp_stack[(*stack_count)].j = jj; *u = k - 1; *i = k; *j = jj; return 1; } } } break; case 2: for (k = jj - turn - 1; k >= ii; k--) { if (with_gquad) { if (fij == my_fc[k - 1] + my_ggg[idx[jj] + k]) { *u = k - 1; *i = *j = -1; vrna_BT_gquad_mfe(vc, k, jj, bp_stack, stack_count); return 1; } } if (hc->matrix[idx[jj] + k] & VRNA_CONSTRAINT_CONTEXT_EXT_LOOP) { mm5 = ((k > 1) && (sn[k - 1] == sn[k])) ? S1[k - 1] : -1; type = (unsigned char)ptype[idx[jj] + k]; if (type == 0) type = 7; en = my_c[idx[jj] + k]; if (sn[k] != sn[jj]) en += P->DuplexInit; if (fij == my_fc[k - 1] + en + E_ExtLoop(type, mm5, mm3, P)) { bp_stack[++(*stack_count)].i = k; bp_stack[(*stack_count)].j = jj; *u = k - 1; *i = k; *j = jj; return 1; } } } break; default: for (k = jj - turn - 1; k >= ii; k--) { if (with_gquad) { if (fij == my_fc[k - 1] + my_ggg[idx[jj] + k]) { *u = k - 1; *i = *j = -1; vrna_BT_gquad_mfe(vc, k, jj, bp_stack, stack_count); return 1; } } if (hc->matrix[idx[jj] + k] & VRNA_CONSTRAINT_CONTEXT_EXT_LOOP) { type = (unsigned char)ptype[idx[jj] + k]; if (type == 0) type = 7; en = my_c[idx[jj] + k]; if (sn[k] != sn[jj]) en += P->DuplexInit; if (fij == my_fc[k - 1] + en + E_ExtLoop(type, -1, -1, P)) { bp_stack[++(*stack_count)].i = k; bp_stack[(*stack_count)].j = jj; *u = k - 1; *i = k; *j = jj; return 1; } if (hc->up_ext[k - 1]) { if ((k > 1) && (sn[k - 1] == sn[k])) { mm5 = S1[k - 1]; if (sc) if (sc->energy_up) en += sc->energy_up[k - 1][1]; if (fij == my_fc[k - 2] + en + E_ExtLoop(type, mm5, -1, P)) { bp_stack[++(*stack_count)].i = k; bp_stack[(*stack_count)].j = jj; *u = k - 2; *i = k; *j = jj; return 1; } } } } if (hc->matrix[idx[jj - 1] + k] & VRNA_CONSTRAINT_CONTEXT_EXT_LOOP) { type = (unsigned char)ptype[idx[jj - 1] + k]; if (type == 0) type = 7; if (hc->up_ext[jj]) { if (sn[jj - 1] == sn[jj]) { mm3 = S1[jj]; en = my_c[idx[jj - 1] + k]; if (sn[k] != sn[jj - 1]) en += P->DuplexInit; /* ??? */ if (sc) if (sc->energy_up) en += sc->energy_up[jj][1]; if (fij == en + my_fc[k - 1] + E_ExtLoop(type, -1, mm3, P)) { bp_stack[++(*stack_count)].i = k; bp_stack[(*stack_count)].j = jj - 1; *u = k - 1; *i = k; *j = jj - 1; return 1; } if (k > ii) { if (hc->up_ext[k - 1]) { mm5 = (sn[k - 1] == sn[k]) ? S1[k - 1] : -1; if (sc) if (sc->energy_up) en += sc->energy_up[k - 1][1]; if (fij == my_fc[k - 2] + en + E_ExtLoop(type, mm5, mm3, P)) { bp_stack[++(*stack_count)].i = k; bp_stack[(*stack_count)].j = jj - 1; *u = k - 2; *i = k; *j = jj - 1; return 1; } } } } } } } break; } } return 0; } PUBLIC int vrna_BT_mb_loop_split(vrna_fold_compound_t *vc, int *i, int *j, int *k, int *l, int *component1, int *component2, vrna_bp_stack_t *bp_stack, int *stack_count) { unsigned char type, type_2; char *ptype; short *S1; int ij, ii, jj, fij, fi, u, en, *my_c, *my_fML, *my_ggg, turn, *idx, with_gquad, dangle_model, *rtype, kk, cnt, with_ud; vrna_param_t *P; vrna_md_t *md; vrna_hc_t *hc; vrna_sc_t *sc; vrna_ud_t *domains_up; vrna_callback_hc_evaluate *evaluate; struct default_data hc_dat_local; P = vc->params; md = &(P->model_details); hc = vc->hc; sc = vc->sc; idx = vc->jindx; ptype = vc->ptype; rtype = &(md->rtype[0]); S1 = vc->sequence_encoding; domains_up = vc->domains_up; my_c = vc->matrices->c; my_fML = vc->matrices->fML; my_ggg = vc->matrices->ggg; turn = md->min_loop_size; with_gquad = md->gquad; with_ud = (domains_up && domains_up->energy_cb) ? 1 : 0; dangle_model = md->dangles; hc_dat_local.idx = vc->jindx; hc_dat_local.mx = hc->matrix; hc_dat_local.hc_up = hc->up_ml; hc_dat_local.cp = vc->cutpoint; if (hc->f) { evaluate = &hc_default_user; hc_dat_local.hc_f = hc->f; hc_dat_local.hc_dat = hc->data; } else { evaluate = &hc_default; } ii = *i; jj = *j; if (with_ud) { /* nibble off unpaired stretches at 3' site */ do { fij = my_fML[idx[jj] + ii]; fi = INF; /* process regular unpaired nucleotides (unbound by ligand) first */ if (evaluate(ii, jj, ii, jj - 1, VRNA_DECOMP_ML_ML, &hc_dat_local)) { fi = my_fML[idx[jj - 1] + ii] + P->MLbase; if (sc) { if (sc->energy_up) fi += sc->energy_up[jj][1]; if (sc->f) fi += sc->f(ii, jj, ii, jj - 1, VRNA_DECOMP_ML_ML, sc->data); } if (jj == ii) return 0; /* no more pairs */ if (fij == fi) { jj--; continue; } } /* next try to nibble off ligand */ for (cnt = 0; cnt < domains_up->uniq_motif_count; cnt++) { u = domains_up->uniq_motif_size[cnt]; kk = jj - u + 1; if ((kk >= ii) && evaluate(ii, jj, ii, jj - u, VRNA_DECOMP_ML_ML, &hc_dat_local)) { en = domains_up->energy_cb(vc, kk, jj, VRNA_UNSTRUCTURED_DOMAIN_MB_LOOP | VRNA_UNSTRUCTURED_DOMAIN_MOTIF, domains_up->data); if (sc) { if (sc->energy_up) en += sc->energy_up[kk][u]; if (sc->f) en += sc->f(ii, jj, ii, jj - u, VRNA_DECOMP_ML_ML, sc->data); } fi = my_fML[idx[kk - 1] + ii] + u * P->MLbase; fi += en; if (fij == fi) { /* skip remaining motifs after first hit */ jj = kk - 1; break; } } } if (jj < ii) return 0; /* no more pairs */ } while (fij == fi); /* nibble off unpaired stretches at 5' site */ do { fij = my_fML[idx[jj] + ii]; fi = INF; /* again, process regular unpaired nucleotides (unbound by ligand) first */ if (evaluate(ii, jj, ii + 1, jj, VRNA_DECOMP_ML_ML, &hc_dat_local)) { fi = my_fML[idx[jj] + ii + 1] + P->MLbase; if (sc) { if (sc->energy_up) fi += sc->energy_up[ii][1]; if (sc->f) fi += sc->f(ii, jj, ii + 1, jj, VRNA_DECOMP_ML_ML, sc->data); } if (ii + 1 == jj) return 0; /* no more pairs */ if (fij == fi) { ii++; continue; } } /* next try to nibble off ligand again */ for (cnt = 0; cnt < domains_up->uniq_motif_count; cnt++) { u = domains_up->uniq_motif_size[cnt]; kk = ii + u - 1; if ((kk <= jj) && evaluate(ii, jj, ii + u, jj, VRNA_DECOMP_ML_ML, &hc_dat_local)) { en = domains_up->energy_cb(vc, ii, kk, VRNA_UNSTRUCTURED_DOMAIN_MB_LOOP | VRNA_UNSTRUCTURED_DOMAIN_MOTIF, domains_up->data); if (sc) { if (sc->energy_up) en += sc->energy_up[ii][u]; if (sc->f) en += sc->f(ii, jj, ii + u, jj, VRNA_DECOMP_ML_ML, sc->data); } fi = my_fML[idx[jj] + kk + 1] + u * P->MLbase; fi += en; if (fij == fi) { /* skip remaining motifs after first hit */ ii = kk + 1; break; } } } if (ii > jj) return 0; /* no more pairs */ } while (fij == fi); } else { /* nibble off unpaired 3' bases */ do { fij = my_fML[idx[jj] + ii]; fi = INF; if (evaluate(ii, jj, ii, jj - 1, VRNA_DECOMP_ML_ML, &hc_dat_local)) { fi = my_fML[idx[jj - 1] + ii] + P->MLbase; if (sc) { if (sc->energy_up) fi += sc->energy_up[jj][1]; if (sc->f) fi += sc->f(ii, jj, ii, jj - 1, VRNA_DECOMP_ML_ML, sc->data); } } if (--jj == 0) break; } while (fij == fi); jj++; /* nibble off unpaired 5' bases */ do { fij = my_fML[idx[jj] + ii]; fi = INF; if (evaluate(ii, jj, ii + 1, jj, VRNA_DECOMP_ML_ML, &hc_dat_local)) { fi = my_fML[idx[jj] + ii + 1] + P->MLbase; if (sc) { if (sc->energy_up) fi += sc->energy_up[ii][1]; if (sc->f) fi += sc->f(ii, jj, ii + 1, jj, VRNA_DECOMP_ML_ML, sc->data); } } if (++ii == jj) break; } while (fij == fi); ii--; if (jj < ii + turn + 1) /* no more pairs */ return 0; } ij = idx[jj] + ii; *component1 = *component2 = 1; /* split into two multi loop parts by default */ /* 1. test for single component */ if (with_gquad) { if (fij == my_ggg[ij] + E_MLstem(0, -1, -1, P)) { *i = *j = -1; *k = *l = -1; vrna_BT_gquad_mfe(vc, ii, jj, bp_stack, stack_count); return 1; } } type = (unsigned char)ptype[ij]; en = my_c[ij]; if (sc) if (sc->f) en += sc->f(ii, jj, ii, jj, VRNA_DECOMP_ML_STEM, sc->data); switch (dangle_model) { case 0: if (evaluate(ii, jj, ii, jj, VRNA_DECOMP_ML_STEM, &hc_dat_local)) { if (type == 0) type = 7; if (fij == en + E_MLstem(type, -1, -1, P)) { *i = *j = -1; *k = ii; *l = jj; *component2 = 2; /* 2nd part is structure enclosed by base pair */ return 1; } } break; case 2: if (evaluate(ii, jj, ii, jj, VRNA_DECOMP_ML_STEM, &hc_dat_local)) { if (type == 0) type = 7; if (fij == en + E_MLstem(type, S1[ii - 1], S1[jj + 1], P)) { *i = *j = -1; *k = ii; *l = jj; *component2 = 2; return 1; } } break; default: if (evaluate(ii, jj, ii, jj, VRNA_DECOMP_ML_STEM, &hc_dat_local)) { if (type == 0) type = 7; if (fij == en + E_MLstem(type, -1, -1, P)) { *i = *j = -1; *k = ii; *l = jj; *component2 = 2; return 1; } } if (evaluate(ii, jj, ii + 1, jj, VRNA_DECOMP_ML_STEM, &hc_dat_local)) { int tmp_en = fij; if (sc) { if (sc->energy_up) tmp_en -= sc->energy_up[ii][1]; if (sc->f) tmp_en -= sc->f(ii, jj, ii + 1, jj, VRNA_DECOMP_ML_STEM, sc->data); } type = (unsigned char)ptype[ij + 1]; if (type == 0) type = 7; if (tmp_en == my_c[ij + 1] + E_MLstem(type, S1[ii], -1, P) + P->MLbase) { *i = *j = -1; *k = ii + 1; *l = jj; *component2 = 2; return 1; } } if (evaluate(ii, jj, ii, jj - 1, VRNA_DECOMP_ML_STEM, &hc_dat_local)) { int tmp_en = fij; if (sc) { if (sc->energy_up) tmp_en -= sc->energy_up[jj][1]; if (sc->f) tmp_en -= sc->f(ii, jj, ii, jj - 1, VRNA_DECOMP_ML_STEM, sc->data); } type = (unsigned char)ptype[idx[jj - 1] + ii]; if (type == 0) type = 7; if (tmp_en == my_c[idx[jj - 1] + ii] + E_MLstem(type, -1, S1[jj], P) + P->MLbase) { *i = *j = -1; *k = ii; *l = jj - 1; *component2 = 2; return 1; } } if (evaluate(ii, jj, ii + 1, jj - 1, VRNA_DECOMP_ML_STEM, &hc_dat_local)) { int tmp_en = fij; if (sc) { if (sc->energy_up) tmp_en -= sc->energy_up[ii][1] + sc->energy_up[jj][1]; if (sc->f) tmp_en -= sc->f(ii, jj, ii + 1, jj - 1, VRNA_DECOMP_ML_STEM, sc->data); } type = (unsigned char)ptype[idx[jj - 1] + ii + 1]; if (type == 0) type = 7; if (tmp_en == my_c[idx[jj - 1] + ii + 1] + E_MLstem(type, S1[ii], S1[jj], P) + 2 * P->MLbase) { *i = *j = -1; *k = ii + 1; *l = jj - 1; *component2 = 2; return 1; } } break; } /* 2. Test for possible split point */ for (u = ii + 1 + turn; u <= jj - 2 - turn; u++) { en = my_fML[idx[u] + ii] + my_fML[idx[jj] + u + 1]; if (sc) if (sc->f) en += sc->f(ii, jj, u, u + 1, VRNA_DECOMP_ML_ML_ML, sc->data); if (fij == en) { *i = ii; *j = u; *k = u + 1; *l = jj; return 1; } } /* 3. last chance! Maybe coax stack */ if (dangle_model == 3) { int ik, k1j, tmp_en; for (k1j = idx[jj] + ii + turn + 2, u = ii + 1 + turn; u <= jj - 2 - turn; u++, k1j++) { ik = idx[u] + ii; if (evaluate(ii, u, u + 1, jj, VRNA_DECOMP_ML_COAXIAL_ENC, &hc_dat_local)) { type = rtype[(unsigned char)ptype[ik]]; type_2 = rtype[(unsigned char)ptype[k1j]]; if (type == 0) type = 7; if (type_2 == 0) type_2 = 7; tmp_en = my_c[ik] + my_c[k1j] + P->stack[type][type_2] + 2 * P->MLintern[1]; if (sc) if (sc->f) tmp_en += sc->f(ii, u, u + 1, jj, VRNA_DECOMP_ML_COAXIAL, sc->data); if (fij == tmp_en) { *i = ii; *j = u; *k = u + 1; *l = jj; *component1 = *component2 = 2; return 1; } } } } return 0; } PUBLIC int vrna_BT_mb_loop(vrna_fold_compound_t *vc, int *i, int *j, int *k, int en, int *component1, int *component2) { unsigned char type, type_2, tt; char *ptype; short s5, s3, *S1; unsigned int *sn; int ij, p, q, r, e, tmp_en, cp, *idx, turn, dangle_model, *my_c, *my_fML, *my_fc, *rtype; vrna_param_t *P; vrna_md_t *md; vrna_hc_t *hc; vrna_sc_t *sc; vrna_callback_hc_evaluate *evaluate; struct default_data hc_dat_local; cp = vc->cutpoint; idx = vc->jindx; ij = idx[*j] + *i; S1 = vc->sequence_encoding; P = vc->params; md = &(P->model_details); sn = vc->strand_number; hc = vc->hc; sc = vc->sc; my_c = vc->matrices->c; my_fML = vc->matrices->fML; my_fc = vc->matrices->fc; turn = md->min_loop_size; ptype = vc->ptype; rtype = &(md->rtype[0]); type = (unsigned char)ptype[ij]; tt = type; type = rtype[type]; dangle_model = md->dangles; hc_dat_local.idx = vc->jindx; hc_dat_local.mx = hc->matrix; hc_dat_local.hc_up = hc->up_ml; hc_dat_local.cp = vc->cutpoint; if (hc->f) { evaluate = &hc_default_user; hc_dat_local.hc_f = hc->f; hc_dat_local.hc_dat = hc->data; } else { evaluate = &hc_default; } p = *i + 1; q = *j - 1; r = q - turn - 1; if (evaluate(*i, *j, p, q, VRNA_DECOMP_PAIR_ML, &hc_dat_local)) { if (type == 0) type = 7; if (tt == 0) tt = 7; /* is it a fake multi-loop? */ /* NOTE: do we really want to evaluate it hard-constraint-wise as a multibranch loop? */ if (sn[*i] != sn[*j]) { int ii, jj; ii = jj = 0; e = my_fc[p] + my_fc[q]; if (sc) if (sc->energy_bp) e += sc->energy_bp[ij]; s5 = (sn[q] == sn[*j]) ? S1[q] : -1; s3 = (sn[*i] == sn[p]) ? S1[p] : -1; switch (dangle_model) { case 0: if (en == e + E_ExtLoop(type, -1, -1, P)) ii = p, jj = q; break; case 2: if (en == e + E_ExtLoop(type, s5, s3, P)) ii = p, jj = q; break; default: if (en == e + E_ExtLoop(type, -1, -1, P)) { ii = p, jj = q; break; } if (hc->up_ext[p]) { e = my_fc[p + 1] + my_fc[q]; if (sc) { if (sc->energy_up) e += sc->energy_up[p][1]; if (sc->energy_bp) e += sc->energy_bp[ij]; } if (en == e + E_ExtLoop(type, -1, s3, P)) { ii = p + 1; jj = q; break; } } if (hc->up_ext[q]) { e = my_fc[p] + my_fc[q - 1]; if (sc) { if (sc->energy_up) e += sc->energy_up[q][1]; if (sc->energy_bp) e += sc->energy_bp[ij]; } if (en == e + E_ExtLoop(type, s5, -1, P)) { ii = p; jj = q - 1; break; } } if ((hc->up_ext[q]) && (hc->up_ext[p])) { e = my_fc[p + 1] + my_fc[q - 1]; if (sc) { if (sc->energy_up) e += sc->energy_up[p][1] + sc->energy_up[q][1]; if (sc->energy_bp) e += sc->energy_bp[ij]; } if (en == e + E_ExtLoop(type, s5, s3, P)) { ii = p + 1; jj = q - 1; break; } } break; } if (ii) { /* found a decomposition */ *component1 = 3; *i = ii; *k = cp - 1; *j = jj; *component2 = 4; return 1; } } /* true multi loop? */ *component1 = *component2 = 1; /* both components are MB loop parts by default */ s5 = (sn[q] == sn[*j]) ? S1[q] : -1; s3 = (sn[*i] == sn[p]) ? S1[p] : -1; switch (dangle_model) { case 0: e = en - E_MLstem(type, -1, -1, P) - P->MLclosing; if (sc) { if (sc->energy_bp) e -= sc->energy_bp[ij]; if (sc->f) e -= sc->f(*i, *j, p, q, VRNA_DECOMP_PAIR_ML, sc->data); } for (r = *i + 2 + turn; r < *j - 2 - turn; ++r) { if (evaluate(p, q, r, r + 1, VRNA_DECOMP_ML_ML_ML, &hc_dat_local)) { tmp_en = my_fML[idx[r] + p] + my_fML[idx[q] + r + 1]; if (sc) if (sc->f) tmp_en += sc->f(p, q, r, r + 1, VRNA_DECOMP_ML_ML_ML, sc->data); if (e == tmp_en) break; } } break; case 2: e = en - E_MLstem(type, s5, s3, P) - P->MLclosing; if (sc) { if (sc->energy_bp) e -= sc->energy_bp[ij]; if (sc->f) e -= sc->f(*i, *j, p, q, VRNA_DECOMP_PAIR_ML, sc->data); } for (r = p + turn + 1; r < q - turn - 1; ++r) { if (evaluate(p, q, r, r + 1, VRNA_DECOMP_ML_ML_ML, &hc_dat_local)) { tmp_en = my_fML[idx[r] + p] + my_fML[idx[q] + r + 1]; if (sc) if (sc->f) tmp_en += sc->f(p, q, r, r + 1, VRNA_DECOMP_ML_ML_ML, sc->data); if (e == tmp_en) break; } } break; default: e = en - P->MLclosing; if (sc) if (sc->energy_bp) e -= sc->energy_bp[ij]; for (r = p + turn + 1; r < q - turn - 1; ++r) { if (evaluate(p, q, r, r + 1, VRNA_DECOMP_ML_ML_ML, &hc_dat_local)) { tmp_en = my_fML[idx[r] + p] + my_fML[idx[q] + r + 1] + E_MLstem(type, -1, -1, P); if (sc) { if (sc->f) { tmp_en += sc->f(*i, *j, p, q, VRNA_DECOMP_PAIR_ML, sc->data); tmp_en += sc->f(p, q, r, r + 1, VRNA_DECOMP_ML_ML_ML, sc->data); } } if (e == tmp_en) break; } if (evaluate(*i, *j, p + 1, q, VRNA_DECOMP_PAIR_ML, &hc_dat_local)) { if (evaluate(p + 1, q, r, r + 1, VRNA_DECOMP_ML_ML_ML, &hc_dat_local)) { tmp_en = e; if (sc) { if (sc->energy_up) tmp_en -= sc->energy_up[p][1]; if (sc->f) { tmp_en -= sc->f(*i, *j, p + 1, q, VRNA_DECOMP_PAIR_ML, sc->data); tmp_en -= sc->f(p + 1, q, r, r + 1, VRNA_DECOMP_ML_ML_ML, sc->data); } } if (tmp_en == my_fML[idx[r] + p + 1] + my_fML[idx[q] + r + 1] + E_MLstem(type, -1, s3, P) + P->MLbase) { p += 1; break; } } } if (evaluate(*i, *j, p, q - 1, VRNA_DECOMP_PAIR_ML, &hc_dat_local)) { if (evaluate(p, q - 1, r, r + 1, VRNA_DECOMP_ML_ML_ML, &hc_dat_local)) { tmp_en = e; if (sc) { if (sc->energy_up) tmp_en -= sc->energy_up[q][1]; if (sc->f) { tmp_en -= sc->f(*i, *j, p, q - 1, VRNA_DECOMP_PAIR_ML, sc->data); tmp_en -= sc->f(p, q - 1, r, r + 1, VRNA_DECOMP_ML_ML_ML, sc->data); } } if (tmp_en == my_fML[idx[r] + p] + my_fML[idx[q - 1] + r + 1] + E_MLstem(type, s5, -1, P) + P->MLbase) { q -= 1; break; } } } if (evaluate(*i, *j, p + 1, q - 1, VRNA_DECOMP_PAIR_ML, &hc_dat_local)) { if (evaluate(p + 1, q - 1, r, r + 1, VRNA_DECOMP_ML_ML_ML, &hc_dat_local)) { tmp_en = e; if (sc) { if (sc->energy_up) tmp_en -= sc->energy_up[p][1] + sc->energy_up[q][1]; if (sc->f) { tmp_en -= sc->f(*i, *j, p + 1, q - 1, VRNA_DECOMP_PAIR_ML, sc->data); tmp_en -= sc->f(p + 1, q - 1, r, r + 1, VRNA_DECOMP_ML_ML_ML, sc->data); } } if (tmp_en == my_fML[idx[r] + p + 1] + my_fML[idx[q - 1] + r + 1] + E_MLstem(type, s5, s3, P) + 2 * P->MLbase) { p += 1; q -= 1; break; } } } /* coaxial stacking of (i.j) with (i+1.r) or (r.j-1) */ /* use MLintern[1] since coax stacked pairs don't get TerminalAU */ if (dangle_model == 3) { tmp_en = e; if (evaluate(*i, *j, p, r, VRNA_DECOMP_ML_COAXIAL, &hc_dat_local)) { type_2 = rtype[(unsigned char)ptype[idx[r] + p]]; if (type_2 == 0) type_2 = 7; tmp_en = my_c[idx[r] + p] + P->stack[tt][type_2] + my_fML[idx[q] + r + 1]; if (sc) { if (sc->f) { tmp_en += sc->f(*i, *j, p, q, VRNA_DECOMP_PAIR_ML, sc->data); tmp_en += sc->f(*i, *j, p, r, VRNA_DECOMP_ML_COAXIAL, sc->data); } } if (e == tmp_en + 2 * P->MLintern[1]) { *component1 = 2; break; } } if (evaluate(*i, *j, r + 1, q, VRNA_DECOMP_ML_COAXIAL, &hc_dat_local)) { type_2 = rtype[(unsigned char)ptype[idx[q] + r + 1]]; if (type_2 == 0) type_2 = 7; tmp_en = my_c[idx[q] + r + 1] + P->stack[tt][type_2] + my_fML[idx[r] + p]; if (sc) { if (sc->f) { tmp_en += sc->f(*i, *j, p, q, VRNA_DECOMP_PAIR_ML, sc->data); tmp_en += sc->f(*i, *j, r + 1, q, VRNA_DECOMP_ML_COAXIAL, sc->data); } } if (e == tmp_en + 2 * P->MLintern[1]) { *component2 = 2; break; } } } } break; } } if (r <= *j - turn - 3) { *i = p; *k = r; *j = q; return 1; } else { #if 0 /* Y shaped ML loops fon't work yet */ if (dangle_model == 3) { d5 = P->dangle5[tt][S1[j - 1]]; d3 = P->dangle3[tt][S1[i + 1]]; /* (i,j) must close a Y shaped ML loop with coax stacking */ if (cij == fML[indx[j - 2] + i + 2] + mm + d3 + d5 + P->MLbase + P->MLbase) { i1 = i + 2; j1 = j - 2; } else if (cij == fML[indx[j - 2] + i + 1] + mm + d5 + P->MLbase) { j1 = j - 2; } else if (cij == fML[indx[j - 1] + i + 2] + mm + d3 + P->MLbase) { i1 = i + 2; } else /* last chance */ if (cij != fML[indx[j - 1] + i + 1] + mm + P->MLbase) { fprintf(stderr, "backtracking failed in repeat"); } /* if we arrive here we can express cij via fML[i1,j1]+dangles */ bt_stack[++s].i = i1; bt_stack[s].j = j1; } #endif } return 0; } PUBLIC vrna_mx_pf_aux_ml_t * vrna_exp_E_ml_fast_init(vrna_fold_compound_t *vc) { vrna_mx_pf_aux_ml_t *aux_mx = NULL; if (vc) { int i, j, d, n, u, turn, ij, *iidx; FLT_OR_DBL *qm; n = (int)vc->length; iidx = vc->iindx; turn = vc->exp_params->model_details.min_loop_size; qm = vc->exp_matrices->qm; /* allocate memory for helper arrays */ aux_mx = (vrna_mx_pf_aux_ml_t *)vrna_alloc(sizeof(vrna_mx_pf_aux_ml_t)); aux_mx->qqm = (FLT_OR_DBL *)vrna_alloc(sizeof(FLT_OR_DBL) * (n + 2)); aux_mx->qqm1 = (FLT_OR_DBL *)vrna_alloc(sizeof(FLT_OR_DBL) * (n + 2)); aux_mx->qqmu_size = 0; aux_mx->qqmu = NULL; if (vc->type == VRNA_FC_TYPE_SINGLE) { vrna_ud_t *domains_up = vc->domains_up; int with_ud = (domains_up && domains_up->exp_energy_cb); int ud_max_size = 0; /* pre-processing ligand binding production rule(s) and auxiliary memory */ if (with_ud) { for (u = 0; u < domains_up->uniq_motif_count; u++) if (ud_max_size < domains_up->uniq_motif_size[u]) ud_max_size = domains_up->uniq_motif_size[u]; aux_mx->qqmu_size = ud_max_size; aux_mx->qqmu = (FLT_OR_DBL **)vrna_alloc(sizeof(FLT_OR_DBL *) * (ud_max_size + 1)); for (u = 0; u <= ud_max_size; u++) aux_mx->qqmu[u] = (FLT_OR_DBL *)vrna_alloc(sizeof(FLT_OR_DBL) * (n + 2)); } } for (d = 0; d <= turn; d++) for (i = 1; i <= n - d; i++) { j = i + d; ij = iidx[i] - j; if (j > n) continue; qm[ij] = 0.; } } return aux_mx; } PUBLIC void vrna_exp_E_ml_fast_rotate(vrna_fold_compound_t *vc, vrna_mx_pf_aux_ml_t *aux_mx) { if (vc && aux_mx) { int u; FLT_OR_DBL *tmp; tmp = aux_mx->qqm1; aux_mx->qqm1 = aux_mx->qqm; aux_mx->qqm = tmp; /* rotate auxiliary arrays for unstructured domains */ if (aux_mx->qqmu) { tmp = aux_mx->qqmu[aux_mx->qqmu_size]; for (u = aux_mx->qqmu_size; u > 0; u--) aux_mx->qqmu[u] = aux_mx->qqmu[u - 1]; aux_mx->qqmu[0] = tmp; } } } PUBLIC void vrna_exp_E_ml_fast_free(vrna_fold_compound_t *vc, vrna_mx_pf_aux_ml_t *aux_mx) { if (vc && aux_mx) { int u; free(aux_mx->qqm); free(aux_mx->qqm1); if (aux_mx->qqmu) { for (u = 0; u <= aux_mx->qqmu_size; u++) free(aux_mx->qqmu[u]); free(aux_mx->qqmu); } free(aux_mx); } } PUBLIC FLT_OR_DBL vrna_exp_E_ml_fast(vrna_fold_compound_t *vc, int i, int j, vrna_mx_pf_aux_ml_t *aux_mx) { if (vc) { switch (vc->type) { case VRNA_FC_TYPE_SINGLE: return exp_E_ml_fast(vc, i, j, aux_mx); break; case VRNA_FC_TYPE_COMPARATIVE: return exp_E_ml_fast_comparative(vc, i, j, aux_mx); break; default: vrna_message_warning("vrna_exp_E_ml_fast@multibranch_loops.c: Unknown fold_compound type"); return 0.; break; } } else { return 0.; } } PRIVATE FLT_OR_DBL exp_E_ml_fast(vrna_fold_compound_t *vc, int i, int j, vrna_mx_pf_aux_ml_t *aux_mx) { short *S1; unsigned char type; int n, *iidx, k, ij, kl, maxk, ii, with_ud, u, circular, with_gquad, *hc_up_ml; FLT_OR_DBL qbt1, temp, *qm, *qb, *qqm, *qqm1, **qqmu, q_temp, q_temp2, *G, *expMLbase, expMLstem; vrna_md_t *md; vrna_exp_param_t *pf_params; vrna_ud_t *domains_up; vrna_hc_t *hc; vrna_sc_t *sc; vrna_callback_hc_evaluate *evaluate; struct default_data hc_dat_local; n = (int)vc->length; iidx = vc->iindx; ij = iidx[i] - j; qqm = aux_mx->qqm; qqm1 = aux_mx->qqm1; qqmu = aux_mx->qqmu; qm = vc->exp_matrices->qm; qb = vc->exp_matrices->qb; G = vc->exp_matrices->G; expMLbase = vc->exp_matrices->expMLbase; pf_params = vc->exp_params; md = &(pf_params->model_details); hc = vc->hc; sc = vc->sc; domains_up = vc->domains_up; circular = md->circ; with_gquad = md->gquad; with_ud = (domains_up && domains_up->exp_energy_cb); hc_up_ml = hc->up_ml; hc_dat_local.idx = vc->jindx; hc_dat_local.mx = hc->matrix; hc_dat_local.hc_up = hc->up_ml; hc_dat_local.cp = vc->cutpoint; if (hc->f) { evaluate = &hc_default_user; hc_dat_local.hc_f = hc->f; hc_dat_local.hc_dat = hc->data; } else { evaluate = &hc_default; } qbt1 = 0; q_temp = 0.; qqm[i] = 0.; if (with_ud) qqmu[0][i] = 0.; if (with_gquad) expMLstem = exp_E_MLstem(0, -1, -1, pf_params); if (evaluate(i, j, i, j - 1, VRNA_DECOMP_ML_ML, &hc_dat_local)) { q_temp = qqm1[i] * expMLbase[1]; if (sc) { if (sc->exp_energy_up) q_temp *= sc->exp_energy_up[j][1]; if (sc->exp_f) q_temp *= sc->exp_f(i, j, i, j - 1, VRNA_DECOMP_ML_ML, sc->data); } if (with_ud) { int cnt; for (cnt = 0; cnt < domains_up->uniq_motif_count; cnt++) { u = domains_up->uniq_motif_size[cnt]; if (j - u >= i) { if (evaluate(i, j, i, j - u, VRNA_DECOMP_ML_ML, &hc_dat_local)) { q_temp2 = qqmu[u][i] * domains_up->exp_energy_cb(vc, j - u + 1, j, VRNA_UNSTRUCTURED_DOMAIN_MB_LOOP | VRNA_UNSTRUCTURED_DOMAIN_MOTIF, domains_up->data) * expMLbase[u]; if (sc) { if (sc->exp_energy_up) q_temp2 *= sc->exp_energy_up[j - u + 1][u]; if (sc->exp_f) q_temp2 *= sc->exp_f(i, j, i, j - 1, VRNA_DECOMP_ML_ML, sc->data); } q_temp += q_temp2; } } } qqmu[0][i] += q_temp; } qqm[i] += q_temp; } if (evaluate(i, j, i, j, VRNA_DECOMP_ML_STEM, &hc_dat_local)) { S1 = vc->sequence_encoding; type = md->pair[S1[i]][S1[j]]; if (type == 0) type = 7; qbt1 = qb[ij] * exp_E_MLstem(type, ((i > 1) || circular) ? S1[i - 1] : -1, ((j < n) || circular) ? S1[j + 1] : -1, pf_params); if (sc) if (sc->exp_f) qbt1 *= sc->exp_f(i, j, i, j, VRNA_DECOMP_ML_STEM, sc->data); qqm[i] += qbt1; if (with_ud) qqmu[0][i] += qbt1; } if (with_gquad) { /*include gquads into qqm*/ qqm[i] += G[ij] * expMLstem; if (with_ud) qqmu[0][i] += G[ij] * expMLstem; } /* * construction of qm matrix containing multiple loop * partition function contributions from segment i,j */ temp = 0.0; kl = iidx[i] - j + 1; /* ii-k=[i,k-1] */ if (sc && sc->exp_f) { for (k = j; k > i; k--, kl++) { q_temp = qm[kl] * qqm[k]; q_temp *= sc->exp_f(i, j, k - 1, k, VRNA_DECOMP_ML_ML_ML, sc->data); temp += q_temp; } } else { for (k = j; k > i; k--, kl++) temp += qm[kl] * qqm[k]; } maxk = MIN2(i + hc_up_ml[i], j); ii = maxk - i; /* length of unpaired stretch */ if (with_ud) { if (sc) { for (k = maxk; k > i; k--, ii--) { q_temp = expMLbase[ii] * qqm[k]; if (sc->exp_energy_up) q_temp *= sc->exp_energy_up[i][ii]; if (sc->exp_f) q_temp *= sc->exp_f(i, j, k, j, VRNA_DECOMP_ML_ML, sc->data); temp += q_temp; temp += q_temp * domains_up->exp_energy_cb(vc, i, k - 1, VRNA_UNSTRUCTURED_DOMAIN_MB_LOOP, domains_up->data); } } else { for (k = maxk; k > i; k--, ii--) { q_temp = expMLbase[ii] * qqm[k]; temp += q_temp; temp += q_temp * domains_up->exp_energy_cb(vc, i, k - 1, VRNA_UNSTRUCTURED_DOMAIN_MB_LOOP, domains_up->data); } } } else { if (sc) { for (k = maxk; k > i; k--, ii--) { q_temp = expMLbase[ii] * qqm[k]; if (sc->exp_energy_up) q_temp *= sc->exp_energy_up[i][ii]; if (sc->exp_f) q_temp *= sc->exp_f(i, j, k, j, VRNA_DECOMP_ML_ML, sc->data); temp += q_temp; } } else { for (k = maxk; k > i; k--, ii--) temp += expMLbase[ii] * qqm[k]; } } return temp + qqm[i]; } PRIVATE FLT_OR_DBL exp_E_ml_fast_comparative(vrna_fold_compound_t *vc, int i, int j, vrna_mx_pf_aux_ml_t *aux_mx) { unsigned char type; unsigned short **a2s; short **S, **S5, **S3; int n, s, n_seq, *iidx, k, ij, kl, maxk, ii, circular, *hc_up_ml; FLT_OR_DBL qbt1, temp, *qm, *qb, *qqm, *qqm1, q_temp, *expMLbase; vrna_md_t *md; vrna_exp_param_t *pf_params; vrna_hc_t *hc; vrna_sc_t **scs; vrna_callback_hc_evaluate *evaluate; struct default_data hc_dat_local; n = (int)vc->length; n_seq = vc->n_seq; iidx = vc->iindx; ij = iidx[i] - j; S = vc->S; S5 = vc->S5; /* S5[s][i] holds next base 5' of i in sequence s */ S3 = vc->S3; /* Sl[s][i] holds next base 3' of i in sequence s */ a2s = vc->a2s; qqm = aux_mx->qqm; qqm1 = aux_mx->qqm1; qm = vc->exp_matrices->qm; qb = vc->exp_matrices->qb; expMLbase = vc->exp_matrices->expMLbase; pf_params = vc->exp_params; md = &(pf_params->model_details); hc = vc->hc; scs = vc->scs; circular = md->circ; hc_up_ml = hc->up_ml; hc_dat_local.idx = vc->jindx; hc_dat_local.mx = hc->matrix; hc_dat_local.hc_up = hc->up_ml; hc_dat_local.cp = vc->cutpoint; if (hc->f) { evaluate = &hc_default_user; hc_dat_local.hc_f = hc->f; hc_dat_local.hc_dat = hc->data; } else { evaluate = &hc_default; } qbt1 = 0; q_temp = 0.; qqm[i] = 0.; if (evaluate(i, j, i, j - 1, VRNA_DECOMP_ML_ML, &hc_dat_local)) { q_temp = qqm1[i] * expMLbase[1]; if (scs) { for (s = 0; s < n_seq; s++) { if (scs[s]) if (scs[s]->exp_energy_up) q_temp *= scs[s]->exp_energy_up[a2s[s][j]][1]; } } qqm[i] += q_temp; } if (evaluate(i, j, i, j, VRNA_DECOMP_ML_STEM, &hc_dat_local)) { q_temp = qb[ij]; for (s = 0; s < n_seq; s++) { type = md->pair[S[s][i]][S[s][j]]; if (type == 0) type = 7; q_temp *= exp_E_MLstem(type, ((i > 1) || circular) ? S5[s][i] : -1, ((j < n) || circular) ? S3[s][j] : -1, pf_params); } qqm[i] += q_temp; } /* * construction of qm matrix containing multiple loop * partition function contributions from segment i,j */ temp = 0.0; kl = iidx[i] - j + 1; /* ii-k=[i,k-1] */ for (k = j; k > i; k--, kl++) temp += qm[kl] * qqm[k]; maxk = MIN2(i + hc_up_ml[i], j); ii = maxk - i; /* length of unpaired stretch */ if (scs) { for (k = maxk; k > i; k--, ii--) { q_temp = expMLbase[ii] * qqm[k]; for (s = 0; s < n_seq; s++) { if (scs[s]) if (scs[s]->exp_energy_up) q_temp *= scs[s]->exp_energy_up[a2s[s][i]][a2s[s][k] - a2s[s][i]]; } temp += q_temp; } } else { for (k = maxk; k > i; k--, ii--) temp += expMLbase[ii] * qqm[k]; } return temp + qqm[i]; } PRIVATE char hc_default(int i, int j, int k, int l, char d, void *data) { int ij, kl, di, dj, u; char eval; struct default_data *dat = (struct default_data *)data; eval = (char)0; di = k - i; dj = j - l; switch (d) { case VRNA_DECOMP_ML_ML_ML: u = l - k - 1; eval = (char)1; if ((u != 0) && (dat->hc_up[k + 1] < u)) eval = (char)0; break; case VRNA_DECOMP_ML_ML: eval = (char)1; if ((di != 0) && (dat->hc_up[i] < di)) eval = (char)0; if ((dj != 0) && (dat->hc_up[l + 1] < dj)) eval = (char)0; break; case VRNA_DECOMP_ML_STEM: kl = dat->idx[l] + k; if (dat->mx[kl] & VRNA_CONSTRAINT_CONTEXT_MB_LOOP_ENC) { eval = (char)1; if ((di != 0) && (dat->hc_up[i] < di)) eval = (char)0; if ((dj != 0) && (dat->hc_up[l + 1] < dj)) eval = (char)0; } break; case VRNA_DECOMP_PAIR_ML: ij = dat->idx[j] + i; if (dat->mx[ij] & VRNA_CONSTRAINT_CONTEXT_MB_LOOP) { eval = (char)1; di--; dj--; if ((di != 0) && (dat->hc_up[i + 1] < di)) eval = (char)0; if ((dj != 0) && (dat->hc_up[l + 1] < dj)) eval = (char)0; } break; case VRNA_DECOMP_ML_COAXIAL: kl = dat->idx[l] + k; if (dat->mx[kl] & VRNA_CONSTRAINT_CONTEXT_MB_LOOP_ENC) eval = (char)1; break; case VRNA_DECOMP_ML_COAXIAL_ENC: ij = dat->idx[j] + i; kl = dat->idx[l] + k; if ((dat->mx[ij] & VRNA_CONSTRAINT_CONTEXT_MB_LOOP_ENC) && (dat->mx[kl] & VRNA_CONSTRAINT_CONTEXT_MB_LOOP_ENC)) eval = (char)1; break; default: nrerror("wtf"); } return eval; } PRIVATE char hc_default_user(int i, int j, int k, int l, char d, void *data) { char eval; struct default_data *dat = (struct default_data *)data; eval = hc_default(i, j, k, l, d, data); eval = (dat->hc_f(i, j, k, l, d, dat->hc_dat)) ? eval : (char)0; return eval; }