/* partiton function for RNA secondary structures Ivo L Hofacker + Ronny Lorenz Vienna RNA package */ #ifdef HAVE_CONFIG_H #include "config.h" #endif #include #include #include #include #include /* #defines FLT_MAX ... */ #include #include "ViennaRNA/utils.h" #include "ViennaRNA/energy_par.h" #include "ViennaRNA/fold_vars.h" #include "ViennaRNA/loop_energies.h" #include "ViennaRNA/gquad.h" #include "ViennaRNA/constraints.h" #include "ViennaRNA/mfe.h" #include "ViennaRNA/part_func.h" #ifdef _OPENMP #include #endif /* ################################# # GLOBAL VARIABLES # ################################# */ PUBLIC int st_back = 0; /* ################################# # PRIVATE VARIABLES # ################################# */ #ifdef VRNA_BACKWARD_COMPAT /* some backward compatibility stuff */ PRIVATE vrna_fold_compound_t *backward_compat_compound = NULL; PRIVATE int backward_compat = 0; #ifdef _OPENMP #pragma omp threadprivate(backward_compat_compound, backward_compat) #endif #endif /* ################################# # PRIVATE FUNCTION DECLARATIONS # ################################# */ PRIVATE void pf_circ(vrna_fold_compound_t *vc); PRIVATE void pf_linear(vrna_fold_compound_t *vc); PRIVATE void alipf_linear(vrna_fold_compound_t *vc); PRIVATE void wrap_alipf_circ(vrna_fold_compound_t *vc, char *structure); #ifdef VRNA_BACKWARD_COMPAT PRIVATE float wrap_pf_fold( const char *sequence, char *structure, vrna_exp_param_t *parameters, int calculate_bppm, int is_constrained, int is_circular); #endif PRIVATE double wrap_mean_bp_distance(FLT_OR_DBL *p, int length, int *index, int turn); /* ################################# # BEGIN OF FUNCTION DEFINITIONS # ################################# */ PUBLIC float vrna_pf_fold( const char *seq, char *structure, vrna_plist_t **pl){ float free_energy; double mfe; vrna_fold_compound_t *vc; vrna_md_t md; vrna_md_set_default(&md); /* no need to backtrack MFE structure */ md.backtrack = 0; if(!pl){ /* no need for pair probability computations if we do not store them somewhere */ md.compute_bpp = 0; } vc = vrna_fold_compound(seq, &md, 0); mfe = (double)vrna_pf(vc, NULL); vrna_exp_params_rescale(vc, &mfe); free_energy = vrna_pf(vc, structure); /* fill plist */ if(pl){ *pl = vrna_plist_from_probs(vc, /*cut_off:*/ 1e-6); } vrna_fold_compound_free(vc); return free_energy; } PUBLIC float vrna_pf_circfold( const char *seq, char *structure, vrna_plist_t **pl){ float free_energy; double mfe; vrna_fold_compound_t *vc; vrna_md_t md; vrna_md_set_default(&md); md.circ = 1; /* no need to backtrack MFE structure */ md.backtrack = 0; if(!pl){ /* no need for pair probability computations if we do not store them somewhere */ md.compute_bpp = 0; } vc = vrna_fold_compound(seq, &md, 0); mfe = (double)vrna_mfe(vc, NULL); vrna_exp_params_rescale(vc, &mfe); free_energy = vrna_pf(vc, structure); /* fill plist */ if(pl){ *pl = vrna_plist_from_probs(vc, /*cut_off:*/ 1e-6); } vrna_fold_compound_free(vc); return free_energy; } PUBLIC float vrna_pf(vrna_fold_compound_t *vc, char *structure){ int n; FLT_OR_DBL Q; double free_energy; vrna_md_t *md; vrna_exp_param_t *params; vrna_mx_pf_t *matrices; free_energy = (float)(INF/100.); if(vc){ /* make sure, everything is set up properly to start partition function computations */ if(!vrna_fold_compound_prepare(vc, VRNA_OPTION_PF)){ vrna_message_warning("vrna_pf@part_func.c: Failed to prepare vrna_fold_compound"); return free_energy; } n = vc->length; params = vc->exp_params; matrices = vc->exp_matrices; md = &(params->model_details); #ifdef _OPENMP /* Explicitly turn off dynamic threads */ omp_set_dynamic(0); #endif #ifdef SUN4 nonstandard_arithmetic(); #else #ifdef HP9 fpsetfastmode(1); #endif #endif /* call user-defined recursion status callback function */ if(vc->stat_cb) vc->stat_cb(VRNA_STATUS_PF_PRE, vc->auxdata); switch(vc->type){ case VRNA_FC_TYPE_SINGLE: /* do the linear pf fold and fill all matrices */ pf_linear(vc); if(md->circ) pf_circ(vc); /* do post processing step for circular RNAs */ break; case VRNA_FC_TYPE_COMPARATIVE: /* do the linear pf fold and fill all matrices */ alipf_linear(vc); /* calculate post processing step for circular */ /* RNAs */ if(md->circ) wrap_alipf_circ(vc, structure); break; default: vrna_message_warning("vrna_pf@part_func.c: Unrecognized fold compound type"); return free_energy; break; } /* call user-defined recursion status callback function */ if(vc->stat_cb) vc->stat_cb(VRNA_STATUS_PF_POST, vc->auxdata); /* calculate base pairing probability matrix (bppm) */ if(md->compute_bpp){ vrna_pairing_probs(vc, structure); #ifdef VRNA_BACKWARD_COMPAT /* * Backward compatibility: * This block may be removed if deprecated functions * relying on the global variable "pr" vanish from within the package! */ pr = matrices->probs; /* { if(pr) free(pr); pr = (FLT_OR_DBL *) vrna_alloc(sizeof(FLT_OR_DBL) * ((n+1)*(n+2)/2)); memcpy(pr, probs, sizeof(FLT_OR_DBL) * ((n+1)*(n+2)/2)); } */ #endif } if (md->backtrack_type=='C') Q = matrices->qb[vc->iindx[1]-n]; else if (md->backtrack_type=='M') Q = matrices->qm[vc->iindx[1]-n]; else Q = (md->circ) ? matrices->qo : matrices->q[vc->iindx[1]-n]; /* ensemble free energy in Kcal/mol */ if (Q<=FLT_MIN) vrna_message_warning("pf_scale too large"); switch(vc->type){ case VRNA_FC_TYPE_COMPARATIVE: free_energy = (-log(Q)-n*log(params->pf_scale))*params->kT/(1000.0 * vc->n_seq); break; case VRNA_FC_TYPE_SINGLE: /* fall through */ default: free_energy = (-log(Q)-n*log(params->pf_scale))*params->kT/1000.0; break; } #ifdef SUN4 standard_arithmetic(); #else #ifdef HP9 fpsetfastmode(0); #endif #endif } return free_energy; } PRIVATE void pf_linear(vrna_fold_compound_t *vc){ char *hard_constraints; int n, i,j, k, ij, d, *my_iindx, *jindx, with_gquad, turn, with_ud, hc_decompose; FLT_OR_DBL temp, Qmax, qbt1, *q, *qb, *qm, *qm1, *q1k, *qln; double max_real; vrna_ud_t *domains_up; vrna_md_t *md; vrna_hc_t *hc; vrna_mx_pf_t *matrices; vrna_mx_pf_aux_el_t *aux_mx_el; vrna_mx_pf_aux_ml_t *aux_mx_ml; vrna_exp_param_t *pf_params; n = vc->length; my_iindx = vc->iindx; jindx = vc->jindx; matrices = vc->exp_matrices; pf_params = vc->exp_params; hc = vc->hc; domains_up = vc->domains_up; q = matrices->q; qb = matrices->qb; qm = matrices->qm; qm1 = matrices->qm1; q1k = matrices->q1k; qln = matrices->qln; md = &(pf_params->model_details); with_gquad = md->gquad; turn = md->min_loop_size; hard_constraints = hc->matrix; with_ud = (domains_up && domains_up->exp_energy_cb); Qmax = 0; max_real = (sizeof(FLT_OR_DBL) == sizeof(float)) ? FLT_MAX : DBL_MAX; if (with_ud && domains_up->exp_prod_cb) domains_up->exp_prod_cb(vc, domains_up->data); if(with_gquad){ free(vc->exp_matrices->G); vc->exp_matrices->G = get_gquad_pf_matrix(vc->sequence_encoding2, vc->exp_matrices->scale, vc->exp_params); } /* init auxiliary arrays for fast exterior/multibranch loops */ aux_mx_el = vrna_exp_E_ext_fast_init(vc); aux_mx_ml = vrna_exp_E_ml_fast_init(vc); /*array initialization ; qb,qm,q qb,qm,q (i,j) are stored as ((n+1-i)*(n-i) div 2 + n+1-j */ for (d=0; d<=turn; d++) for (i=1; i<=n-d; i++) { j=i+d; ij = my_iindx[i]-j; qb[ij] = 0.0; } for (j = turn + 2; j <= n; j++) { for (i = j - turn - 1; i >= 1; i--) { /* construction of partition function of segment i,j */ /* firstly that given i binds j : qb(i,j) */ ij = my_iindx[i] - j; hc_decompose = hard_constraints[jindx[j] + i]; qbt1 = 0; if(hc_decompose){ /* process hairpin loop(s) */ qbt1 += vrna_exp_E_hp_loop(vc, i, j); /* process interior loop(s) */ qbt1 += vrna_exp_E_int_loop(vc, i, j); /* process multibranch loop(s) */ qbt1 += vrna_exp_E_mb_loop_fast(vc, i, j, aux_mx_ml->qqm1); } qb[ij] = qbt1; /* Multibranch loop */ qm[ij] = vrna_exp_E_ml_fast(vc, i, j, aux_mx_ml); if (qm1) qm1[jindx[j] + i] = aux_mx_ml->qqm[i]; /* for stochastic backtracking and circfold */ /* Exterior loop */ q[ij] = temp = vrna_exp_E_ext_fast(vc, i, j, aux_mx_el); if (temp>Qmax) { Qmax = temp; if (Qmax>max_real/10.) vrna_message_warning("Q close to overflow: %d %d %g", i,j,temp); } if (temp>=max_real) { vrna_message_error("overflow in pf_fold while calculating q[%d,%d]\n" "use larger pf_scale", i,j); } } /* rotate auxiliary arrays */ vrna_exp_E_ext_fast_rotate(vc, aux_mx_el); vrna_exp_E_ml_fast_rotate(vc, aux_mx_ml); } /* prefill linear qln, q1k arrays */ if(q1k && qln){ for (k=1; k<=n; k++) { q1k[k] = q[my_iindx[1] - k]; qln[k] = q[my_iindx[k] - n]; } q1k[0] = 1.0; qln[n+1] = 1.0; } /* free memory occupied by auxiliary arrays for fast exterior/multibranch loops */ vrna_exp_E_ml_fast_free(vc, aux_mx_ml); vrna_exp_E_ext_fast_free(vc, aux_mx_el); } /* calculate partition function for circular case */ /* NOTE: this is the postprocessing step ONLY */ /* You have to call pf_linear first to calculate */ /* complete circular case!!! */ PRIVATE void pf_circ(vrna_fold_compound_t *vc){ char *ptype; short *S1; int u, p, q, k, l, turn, n, *my_iindx, *jindx, *rtype; FLT_OR_DBL *scale, *qb, *qm, *qm1, *qm2, qo, qho, qio, qmo, qbt1, qot, expMLclosing; vrna_exp_param_t *pf_params; vrna_mx_pf_t *matrices; n = vc->length; matrices = vc->exp_matrices; my_iindx = vc->iindx; jindx = vc->jindx; ptype = vc->ptype; S1 = vc->sequence_encoding; pf_params = vc->exp_params; qb = matrices->qb; qm = matrices->qm; qm1 = matrices->qm1; qm2 = matrices->qm2; scale = matrices->scale; expMLclosing = pf_params->expMLclosing; turn = pf_params->model_details.min_loop_size; rtype = &(pf_params->model_details.rtype[0]); qo = qho = qio = qmo = 0.; /* construct qm2 matrix from qm1 entries */ for(k=1; kMAXLOOP) break; lstart = ln1+p-1+n-MAXLOOP; if(lstart MAXLOOP) continue; type2 = ptype[jindx[l] + k]; if(!type2) continue; qio += qb[my_iindx[p]-q] * qb[my_iindx[k]-l] * exp_E_IntLoop(ln2, ln1, rtype[type2], type, S1[l+1], S1[k-1], S1[p-1], S1[q+1], pf_params) * scale[ln1+ln2]; } } /* end of kl double loop */ } } /* end of pq double loop */ /* 3. Multiloops */ for(k=turn+2; kqo = qo; matrices->qho = qho; matrices->qio = qio; matrices->qmo = qmo; } PUBLIC int vrna_pf_float_precision(void){ return (sizeof(FLT_OR_DBL) == sizeof(float)); } PRIVATE void alipf_linear( vrna_fold_compound_t *vc){ char *hard_constraints; int i,j, ij, jij, d, turn, n, *my_iindx, *jindx, *pscore; FLT_OR_DBL temp, Qmax, qbt1, *q, *qb, *qm, *qm1; double kTn, max_real; vrna_exp_param_t *pf_params; vrna_mx_pf_t *matrices; vrna_mx_pf_aux_el_t *aux_mx_el; vrna_mx_pf_aux_ml_t *aux_mx_ml; vrna_md_t *md; vrna_hc_t *hc; n = vc->length; pf_params = vc->exp_params; matrices = vc->exp_matrices; hc = vc->hc; my_iindx = vc->iindx; jindx = vc->jindx; pscore = vc->pscore; /* precomputed array of pair types */ md = &(pf_params->model_details); q = matrices->q; qb = matrices->qb; qm = matrices->qm; qm1 = matrices->qm1; hard_constraints = hc->matrix; turn = md->min_loop_size; kTn = pf_params->kT/10.; /* kT in cal/mol */ Qmax = 0.; max_real = (sizeof(FLT_OR_DBL) == sizeof(float)) ? FLT_MAX : DBL_MAX; /* init auxiliary arrays for fast exterior/multibranch loops */ aux_mx_el = vrna_exp_E_ext_fast_init(vc); aux_mx_ml = vrna_exp_E_ml_fast_init(vc); /* array initialization ; qb,qm,q qb,qm,q (i,j) are stored as ((n+1-i)*(n-i) div 2 + n+1-j */ for (d = 0; d <= turn; d++) for (i = 1; i <= n - d; i++) { j = i + d; ij = my_iindx[i]-j; qb[ij] = 0.0; } for (j = turn + 2; j <= n; j++) { for (i = j - turn - 1; i >= 1; i--) { int psc; /* construction of partition function for segment i,j */ /* calculate pf given that i and j pair: qb(i,j) */ ij = my_iindx[i] - j; jij = jindx[j] + i; psc = pscore[jij]; qbt1 = 0.; if (hard_constraints[jij]) { /* process hairpin loop(s) */ qbt1 += vrna_exp_E_hp_loop(vc, i, j); /* process interior loop(s) */ qbt1 += vrna_exp_E_int_loop(vc, i, j); /* process multibranch loop(s) */ qbt1 += vrna_exp_E_mb_loop_fast(vc, i, j, aux_mx_ml->qqm1); qbt1 *= exp(psc/kTn); } qb[ij] = qbt1; /* Multibranch loop */ qm[ij] = vrna_exp_E_ml_fast(vc, i, j, aux_mx_ml); if (qm1) qm1[jindx[j] + i] = aux_mx_ml->qqm[i]; /* for stochastic backtracking and circfold */ /* Exterior loop */ q[ij] = temp = vrna_exp_E_ext_fast(vc, i, j, aux_mx_el); if (temp > Qmax) { Qmax = temp; if (Qmax > max_real/10.) vrna_message_warning("Q close to overflow: %d %d %g", i,j,temp); } if (temp >= max_real) { vrna_message_error("overflow in pf_fold while calculating q[%d,%d]\n" "use larger pf_scale", i,j); } } /* rotate auxiliary arrays */ vrna_exp_E_ext_fast_rotate(vc, aux_mx_el); vrna_exp_E_ml_fast_rotate(vc, aux_mx_ml); } /* free memory occupied by auxiliary arrays for fast exterior/multibranch loops */ vrna_exp_E_ml_fast_free(vc, aux_mx_ml); vrna_exp_E_ext_fast_free(vc, aux_mx_el); } /* calculate partition function for circular case */ /* NOTE: this is the postprocessing step ONLY */ /* You have to call alipf_linear first to calculate */ /* circular case!!! */ PRIVATE void wrap_alipf_circ(vrna_fold_compound_t *vc, char *structure){ char **Ss, *hard_constraints; unsigned short **a2s; short **S, **S5, **S3; int u, p, q, pq, k, l, s, *type, n_seq, n, *my_iindx, *jindx, *rtype; FLT_OR_DBL qbt1, qot, qo, qho, qio, qmo, *qb, *qm, *qm1, *qm2, *scale, expMLclosing; vrna_exp_param_t *pf_params; vrna_mx_pf_t *matrices; vrna_md_t *md; vrna_hc_t *hc; vrna_sc_t **sc; n_seq = vc->n_seq; n = vc->length; 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 */ Ss = vc->Ss; a2s = vc->a2s; pf_params = vc->exp_params; matrices = vc->exp_matrices; my_iindx = vc->iindx; jindx = vc->jindx; hc = vc->hc; sc = vc->scs; qb = matrices->qb; qm = matrices->qm; qm1 = matrices->qm1; qm2 = matrices->qm2; scale = matrices->scale; expMLclosing = pf_params->expMLclosing; md = &(pf_params->model_details); hard_constraints = hc->matrix; rtype = &(md->rtype[0]); qo = qho = qio = qmo = 0.; type = (int *)vrna_alloc(sizeof(int) * n_seq); /* calculate the qm2 matrix */ for(k=1; kpair[S[s][p]][S[s][q]]; if (type[s]==0) type[s]=7; } /* 1. exterior hairpin contribution */ /* Note, that we do not scale Hairpin Energy by u+2 but by u cause the scale */ /* for the closing pair was already done in the forward recursion */ if(hard_constraints[pq] & VRNA_CONSTRAINT_CONTEXT_HP_LOOP){ if(hc->up_hp[q+1] > u){ for (qbt1=1,s=0; sexp_energy_up){ qbt1 *= ((p > 1) ? sc[s]->exp_energy_up[1][a2s[s][p]-1] : 1.) * ((q < n) ? sc[s]->exp_energy_up[a2s[s][q]+1][a2s[s][n] - a2s[s][q]] : 1.); } } } qho += qb[my_iindx[p]-q] * qbt1 * scale[u]; } } /* 2. exterior interior loop contribution*/ if(hard_constraints[pq] & VRNA_CONSTRAINT_CONTEXT_INT_LOOP){ for(k=q+1; k < n; k++){ int ln1, lstart; ln1 = k - q - 1; if(ln1 + p - 1 > MAXLOOP) break; if(hc->up_int[q+1] < ln1) break; lstart = ln1+p-1+n-MAXLOOP; if(lstart MAXLOOP) continue; if(hc->up_int[l+1] < ln2) continue; FLT_OR_DBL qloop=1.; if (qb[my_iindx[k]-l]==0.){ qloop=0.; continue;} for (s=0; spair[S[s][l]][S[s][k]]; if (type_2 == 0) type_2 = 7; qloop *= exp_E_IntLoop(ln1a, ln2a, rt, type_2, S3[s][q], S5[s][p], S5[s][k], S3[s][l], pf_params); } if(sc) for(s = 0; s < n_seq; s++){ int ln1a = a2s[s][k] - 1 - a2s[s][q]; int ln2a = a2s[s][n] - a2s[s][l] + a2s[s][p] - 1; if(sc[s]){ if((ln1a+ln2a == 0) && (sc[s]->exp_energy_stack)){ if(S[s][p] && S[s][q] && S[s][k] && S[s][l]){ /* don't allow gaps in stack */ qloop *= sc[s]->exp_energy_stack[a2s[s][p]] * sc[s]->exp_energy_stack[a2s[s][q]] * sc[s]->exp_energy_stack[a2s[s][k]] * sc[s]->exp_energy_stack[a2s[s][l]]; } } if(sc[s]->exp_energy_up) qloop *= sc[s]->exp_energy_up[a2s[s][q] + 1][ln1a] * ((l < n) ? sc[s]->exp_energy_up[a2s[s][l]+1][a2s[s][n] - a2s[s][l]] : 1.) * ((p > 1) ? sc[s]->exp_energy_up[1][a2s[s][p]-1] : 1.); } } qio += qb[my_iindx[p]-q] * qb[my_iindx[k]-l] * qloop * scale[ln1+ln2]; } } /* end of kl double loop */ } } } /* end of pq double loop */ /* 3. exterior multiloop contribution */ for(k=TURN+2; kup_ext[1] >= n) qo += 1.0 * scale[n]; matrices->qo = qo; matrices->qho = qho; matrices->qio = qio; matrices->qmo = qmo; free(type); } /*###########################################*/ /*# deprecated functions below #*/ /*###########################################*/ #ifdef VRNA_BACKWARD_COMPAT PRIVATE double wrap_mean_bp_distance(FLT_OR_DBL *p, int length, int *index, int turn){ int i,j; double d = 0.; /* compute the mean base pair distance in the thermodynamic ensemble */ /* = \sum_{a,b} p_a p_b d(S_a,S_b) this can be computed from the pair probs p_ij as = \sum_{ij} p_{ij}(1-p_{ij}) */ for (i=1; i<=length; i++) for (j=i+turn+1; j<=length; j++) d += p[index[i]-j] * (1-p[index[i]-j]); return 2*d; } PRIVATE float wrap_pf_fold( const char *sequence, char *structure, vrna_exp_param_t *parameters, int calculate_bppm, int is_constrained, int is_circular){ vrna_fold_compound_t *vc; vrna_md_t md; vc = NULL; /* we need vrna_exp_param_t datastructure to correctly init default hard constraints */ if(parameters) md = parameters->model_details; else{ set_model_details(&md); /* get global default parameters */ } md.circ = is_circular; md.compute_bpp = calculate_bppm; vc = vrna_fold_compound(sequence, &md, VRNA_OPTION_DEFAULT); /* prepare exp_params and set global pf_scale */ vc->exp_params = vrna_exp_params(&md); vc->exp_params->pf_scale = pf_scale; if(is_constrained && structure){ unsigned int constraint_options = 0; constraint_options |= VRNA_CONSTRAINT_DB | VRNA_CONSTRAINT_DB_PIPE | VRNA_CONSTRAINT_DB_DOT | VRNA_CONSTRAINT_DB_X | VRNA_CONSTRAINT_DB_ANG_BRACK | VRNA_CONSTRAINT_DB_RND_BRACK; vrna_constraints_add(vc, (const char *)structure, constraint_options); } if(backward_compat_compound && backward_compat) vrna_fold_compound_free(backward_compat_compound); backward_compat_compound = vc; backward_compat = 1; iindx = backward_compat_compound->iindx; return vrna_pf(vc, structure); } PUBLIC vrna_plist_t * stackProb(double cutoff){ if(!(backward_compat_compound && backward_compat)){ vrna_message_error("stackProb: run pf_fold() first!"); } else if( !backward_compat_compound->exp_matrices->probs){ vrna_message_error("stackProb: probs==NULL!"); } return vrna_stack_prob(backward_compat_compound, cutoff); } PUBLIC char * centroid( int length, double *dist) { if (pr==NULL) vrna_message_error("pr==NULL. You need to call pf_fold() before centroid()"); return vrna_centroid_from_probs(length, dist, pr); } PUBLIC double mean_bp_dist(int length) { /* compute the mean base pair distance in the thermodynamic ensemble */ /* = \sum_{a,b} p_a p_b d(S_a,S_b) this can be computed from the pair probs p_ij as = \sum_{ij} p_{ij}(1-p_{ij}) */ int i, j, *my_iindx; double d = 0; if (pr==NULL) vrna_message_error("pr==NULL. You need to call pf_fold() before mean_bp_dist()"); my_iindx = vrna_idx_row_wise(length); for (i=1; i<=length; i++) for (j=i+TURN+1; j<=length; j++) d += pr[my_iindx[i]-j] * (1-pr[my_iindx[i]-j]); free(my_iindx); return 2*d; } /* get the free energy of a subsequence from the q[] array */ PUBLIC double get_subseq_F( int i, int j){ if(backward_compat_compound) if(backward_compat_compound->exp_matrices) if(backward_compat_compound->exp_matrices->q){ int *my_iindx = backward_compat_compound->iindx; vrna_exp_param_t *pf_params = backward_compat_compound->exp_params; FLT_OR_DBL *q = backward_compat_compound->exp_matrices->q; return ((-log(q[my_iindx[i]-j])-(j-i+1)*log(pf_params->pf_scale))*pf_params->kT/1000.0); } vrna_message_error("call pf_fold() to fill q[] array before calling get_subseq_F()"); return 0.; /* we will never get to this point */ } /*----------------------------------------------------------------------*/ PUBLIC double expHairpinEnergy( int u, int type, short si1, short sj1, const char *string) { /* compute Boltzmann weight of a hairpin loop, multiply by scale[u+2] */ vrna_exp_param_t *pf_params = backward_compat_compound->exp_params; double q, kT; kT = pf_params->kT; /* kT in cal/mol */ if(u <= 30) q = pf_params->exphairpin[u]; else q = pf_params->exphairpin[30] * exp( -(pf_params->lxc*log( u/30.))*10./kT); if ((tetra_loop)&&(u==4)) { char tl[7]={0}, *ts; strncpy(tl, string, 6); if ((ts=strstr(pf_params->Tetraloops, tl))) return (pf_params->exptetra[(ts-pf_params->Tetraloops)/7]); } if ((tetra_loop)&&(u==6)) { char tl[9]={0}, *ts; strncpy(tl, string, 6); if ((ts=strstr(pf_params->Hexaloops, tl))) return (pf_params->exphex[(ts-pf_params->Hexaloops)/9]); } if (u==3) { char tl[6]={0}, *ts; strncpy(tl, string, 5); if ((ts=strstr(pf_params->Triloops, tl))) return (pf_params->exptri[(ts-pf_params->Triloops)/6]); if (type>2) q *= pf_params->expTermAU; } else /* no mismatches for tri-loops */ q *= pf_params->expmismatchH[type][si1][sj1]; return q; } PUBLIC double expLoopEnergy(int u1, int u2, int type, int type2, short si1, short sj1, short sp1, short sq1) { /* compute Boltzmann weight of interior loop, multiply by scale[u1+u2+2] for scaling */ double z=0; int no_close = 0; vrna_exp_param_t *pf_params = backward_compat_compound->exp_params; if ((no_closingGU) && ((type2==3)||(type2==4)||(type==2)||(type==4))) no_close = 1; if ((u1==0) && (u2==0)) /* stack */ z = pf_params->expstack[type][type2]; else if (no_close==0) { if ((u1==0)||(u2==0)) { /* bulge */ int u; u = (u1==0)?u2:u1; z = pf_params->expbulge[u]; if (u2+u1==1) z *= pf_params->expstack[type][type2]; else { if (type>2) z *= pf_params->expTermAU; if (type2>2) z *= pf_params->expTermAU; } } else { /* interior loop */ if (u1+u2==2) /* size 2 is special */ z = pf_params->expint11[type][type2][si1][sj1]; else if ((u1==1) && (u2==2)) z = pf_params->expint21[type][type2][si1][sq1][sj1]; else if ((u1==2) && (u2==1)) z = pf_params->expint21[type2][type][sq1][si1][sp1]; else if ((u1==2) && (u2==2)) z = pf_params->expint22[type][type2][si1][sp1][sq1][sj1]; else if (((u1==2)&&(u2==3))||((u1==3)&&(u2==2))){ /*2-3 is special*/ z = pf_params->expinternal[5]* pf_params->expmismatch23I[type][si1][sj1]* pf_params->expmismatch23I[type2][sq1][sp1]; z *= pf_params->expninio[2][1]; } else if ((u1==1)||(u2==1)) { /*1-n is special*/ z = pf_params->expinternal[u1+u2]* pf_params->expmismatch1nI[type][si1][sj1]* pf_params->expmismatch1nI[type2][sq1][sp1]; z *= pf_params->expninio[2][abs(u1-u2)]; } else { z = pf_params->expinternal[u1+u2]* pf_params->expmismatchI[type][si1][sj1]* pf_params->expmismatchI[type2][sq1][sp1]; z *= pf_params->expninio[2][abs(u1-u2)]; } } } return z; } PUBLIC void init_pf_circ_fold(int length){ /* DO NOTHING */ } PUBLIC void init_pf_fold(int length){ /* DO NOTHING */ } /** *** Allocate memory for all matrices and other stuff **/ PUBLIC void free_pf_arrays(void){ if(backward_compat_compound && backward_compat){ vrna_fold_compound_free(backward_compat_compound); backward_compat_compound = NULL; backward_compat = 0; iindx = NULL; } } PUBLIC FLT_OR_DBL * export_bppm(void){ if(backward_compat_compound) if(backward_compat_compound->exp_matrices) if(backward_compat_compound->exp_matrices->probs) return backward_compat_compound->exp_matrices->probs; return NULL; } /*-------------------------------------------------------------------------*/ /* make arrays used for pf_fold available to other routines */ PUBLIC int get_pf_arrays(short **S_p, short **S1_p, char **ptype_p, FLT_OR_DBL **qb_p, FLT_OR_DBL **qm_p, FLT_OR_DBL **q1k_p, FLT_OR_DBL **qln_p){ if(backward_compat_compound){ if(backward_compat_compound->exp_matrices) if(backward_compat_compound->exp_matrices->qb){ *S_p = backward_compat_compound->sequence_encoding2; *S1_p = backward_compat_compound->sequence_encoding; *ptype_p = backward_compat_compound->ptype_pf_compat; *qb_p = backward_compat_compound->exp_matrices->qb; *qm_p = backward_compat_compound->exp_matrices->qm; *q1k_p = backward_compat_compound->exp_matrices->q1k; *qln_p = backward_compat_compound->exp_matrices->qln; return 1; } } return 0; } /*-----------------------------------------------------------------*/ PUBLIC float pf_fold(const char *sequence, char *structure){ return wrap_pf_fold(sequence, structure, NULL, do_backtrack, fold_constrained, 0); } PUBLIC float pf_circ_fold( const char *sequence, char *structure){ return wrap_pf_fold(sequence, structure, NULL, do_backtrack, fold_constrained, 1); } PUBLIC float pf_fold_par(const char *sequence, char *structure, vrna_exp_param_t *parameters, int calculate_bppm, int is_constrained, int is_circular){ return wrap_pf_fold(sequence, structure, parameters, calculate_bppm, is_constrained, is_circular); } PUBLIC char * pbacktrack(char *seq){ int n = (int)strlen(seq); return vrna_pbacktrack5(backward_compat_compound, n); } PUBLIC char * pbacktrack5(char *seq, int length){ /* the seq parameter must no differ to the one stored globally anyway, so we just ignore it */ return vrna_pbacktrack5(backward_compat_compound, length); } PUBLIC char * pbacktrack_circ(char *seq){ char *structure; vrna_md_t *md; structure = NULL; if(backward_compat_compound){ md = &(backward_compat_compound->exp_params->model_details); if(md->circ && backward_compat_compound->exp_matrices->qm2){ structure = vrna_pbacktrack(backward_compat_compound); } } return structure; } PUBLIC void update_pf_params(int length){ if(backward_compat_compound && backward_compat){ vrna_md_t md; set_model_details(&md); vrna_exp_params_reset(backward_compat_compound, &md); /* compatibility with RNAup, may be removed sometime */ pf_scale = backward_compat_compound->exp_params->pf_scale; } } PUBLIC void update_pf_params_par( int length, vrna_exp_param_t *parameters){ if(backward_compat_compound && backward_compat){ vrna_md_t md; if(parameters){ vrna_exp_params_subst(backward_compat_compound, parameters); } else { set_model_details(&md); vrna_exp_params_reset(backward_compat_compound, &md); } /* compatibility with RNAup, may be removed sometime */ pf_scale = backward_compat_compound->exp_params->pf_scale; } } PUBLIC char * get_centroid_struct_gquad_pr( int length, double *dist){ return vrna_centroid(backward_compat_compound, dist); } PUBLIC void assign_plist_gquad_from_pr( vrna_plist_t **pl, int length, /* ignored */ double cut_off){ if(!backward_compat_compound){ *pl = NULL; } else if( !backward_compat_compound->exp_matrices->probs){ *pl = NULL; } else { *pl = vrna_plist_from_probs(backward_compat_compound, cut_off); } } PUBLIC double mean_bp_distance(int length){ if(backward_compat_compound) if(backward_compat_compound->exp_matrices) if(backward_compat_compound->exp_matrices->probs) return vrna_mean_bp_distance(backward_compat_compound); vrna_message_error("mean_bp_distance: you need to call vrna_pf_fold first"); return 0.; /* we will never get to this point */ } PUBLIC double mean_bp_distance_pr(int length, FLT_OR_DBL *p){ double d=0; int *index = vrna_idx_row_wise((unsigned int) length); if (p==NULL) vrna_message_error("p==NULL. You need to supply a valid probability matrix for mean_bp_distance_pr()"); d = wrap_mean_bp_distance(p, length, index, TURN); free(index); return d; } #endif