1 /*
2 * partiton function and base pair probabilities
3 * for RNA secvondary structures
4 * of a set of aligned sequences
5 *
6 * Ivo L Hofacker
7 * Vienna RNA package
8 */
9
10 /**
11 *** \file alipfold.c
12 **/
13
14 #ifdef HAVE_CONFIG_H
15 #include "config.h"
16 #endif
17
18 #ifndef VRNA_DISABLE_BACKWARD_COMPATIBILITY
19
20 #include <stdio.h>
21 #include <stdlib.h>
22 #include <string.h>
23 #include <math.h>
24 #include <float.h> /* #defines FLT_MIN */
25 #include <limits.h>
26
27 #include "ViennaRNA/utils/basic.h"
28 #include "ViennaRNA/params/default.h"
29 #include "ViennaRNA/fold_vars.h"
30 #include "ViennaRNA/plotting/probabilities.h"
31 #include "ViennaRNA/ribo.h"
32 #include "ViennaRNA/params/basic.h"
33 #include "ViennaRNA/loops/all.h"
34 #include "ViennaRNA/eval.h"
35 #include "ViennaRNA/mfe.h"
36 #include "ViennaRNA/part_func.h"
37 #include "ViennaRNA/utils/structures.h"
38 #include "ViennaRNA/alifold.h"
39
40 #ifdef _OPENMP
41 #include <omp.h>
42 #endif
43
44 /*
45 #################################
46 # PUBLIC GLOBAL VARIABLES #
47 #################################
48 */
49
50 /*
51 #################################
52 # PRIVATE GLOBAL VARIABLES #
53 #################################
54 */
55
56 /* some backward compatibility stuff */
57 PRIVATE vrna_fold_compound_t *backward_compat_compound = NULL;
58 PRIVATE int backward_compat = 0;
59 PRIVATE unsigned short **backward_compat_a2s = NULL;
60
61 #ifdef _OPENMP
62
63 #pragma omp threadprivate(backward_compat_compound, backward_compat, backward_compat_a2s)
64
65 #endif
66
67 /*
68 #################################
69 # PRIVATE FUNCTION DECLARATIONS #
70 #################################
71 */
72
73 PRIVATE float
74 wrap_alipf_fold(const char **sequences,
75 char *structure,
76 plist **pl,
77 vrna_exp_param_t *parameters,
78 int calculate_bppm,
79 int is_constrained,
80 int is_circular);
81
82
83 /*
84 #################################
85 # BEGIN OF FUNCTION DEFINITIONS #
86 #################################
87 */
88 /*-----------------------------------------------------------------*/
89 PRIVATE float
wrap_alipf_fold(const char ** sequences,char * structure,plist ** pl,vrna_exp_param_t * parameters,int calculate_bppm,int is_constrained,int is_circular)90 wrap_alipf_fold(const char **sequences,
91 char *structure,
92 plist **pl,
93 vrna_exp_param_t *parameters,
94 int calculate_bppm,
95 int is_constrained,
96 int is_circular)
97 {
98 int i, n_seq;
99 float free_energy;
100 vrna_fold_compound_t *vc;
101 vrna_md_t md;
102
103 if (sequences == NULL)
104 return 0.;
105
106 for (n_seq = 0; sequences[n_seq]; n_seq++) ; /* count the sequences */
107
108 vc = NULL;
109
110 /*
111 * if present, extract model details from provided parameters variable,
112 * to properly initialize the fold compound. Otherwise use default
113 * settings taken from deprecated global variables
114 */
115 if (parameters)
116 vrna_md_copy(&md, &(parameters->model_details));
117 else
118 set_model_details(&md);
119
120 /* set circular and backtracing options */
121 md.circ = is_circular;
122 md.compute_bpp = calculate_bppm;
123
124 vc = vrna_fold_compound_comparative(sequences, &md, VRNA_OPTION_DEFAULT);
125
126 /*
127 * if present, attach a copy of the parameters structure instead of the
128 * default parameters but take care of re-setting it to (initialized)
129 * model details
130 */
131 free(vc->exp_params);
132 if (parameters) {
133 vrna_md_copy(&(parameters->model_details), &(vc->params->model_details));
134 vc->exp_params = vrna_exp_params_copy(parameters);
135 } else {
136 vc->exp_params = vrna_exp_params_comparative(n_seq, &(vc->params->model_details));
137 }
138
139 /* propagate global pf_scale into vc->exp_params */
140 vc->exp_params->pf_scale = pf_scale;
141
142 if (is_constrained && structure) {
143 unsigned int constraint_options = 0;
144 constraint_options |= VRNA_CONSTRAINT_DB
145 | VRNA_CONSTRAINT_DB_PIPE
146 | VRNA_CONSTRAINT_DB_DOT
147 | VRNA_CONSTRAINT_DB_X
148 | VRNA_CONSTRAINT_DB_ANG_BRACK
149 | VRNA_CONSTRAINT_DB_RND_BRACK;
150
151 vrna_constraints_add(vc, (const char *)structure, constraint_options);
152 }
153
154 if (backward_compat && backward_compat_compound) {
155 for (n_seq = 0; n_seq < backward_compat_compound->n_seq; n_seq++)
156 free(backward_compat_a2s[n_seq]);
157 free(backward_compat_a2s);
158 vrna_fold_compound_free(backward_compat_compound);
159 }
160
161 backward_compat_compound = vc;
162 iindx = backward_compat_compound->iindx;
163
164 /* create alignment-column to sequence position mapping compatibility array */
165 backward_compat_a2s = (unsigned short **)vrna_alloc(sizeof(unsigned short *) * (vc->n_seq + 1));
166 for (n_seq = 0; n_seq < vc->n_seq; n_seq++) {
167 backward_compat_a2s[n_seq] =
168 (unsigned short *)vrna_alloc(sizeof(unsigned short) * (vc->length + 2));
169 for (i = 1; i <= vc->length; i++)
170 backward_compat_a2s[n_seq][i] = (unsigned short)vc->a2s[n_seq][i];
171 }
172
173 backward_compat = 1;
174
175 free_energy = vrna_pf(vc, structure);
176
177 /* fill plist */
178 if (pl && calculate_bppm)
179 *pl = vrna_plist_from_probs(vc, /*cut_off:*/ 1e-6);
180
181 return free_energy;
182 }
183
184
185 /*###########################################*/
186 /*# deprecated functions below #*/
187 /*###########################################*/
188
189 PUBLIC float
alipf_fold(const char ** sequences,char * structure,plist ** pl)190 alipf_fold(const char **sequences,
191 char *structure,
192 plist **pl)
193 {
194 return wrap_alipf_fold(sequences, structure, pl, NULL, do_backtrack, fold_constrained, 0);
195 }
196
197
198 PUBLIC float
alipf_circ_fold(const char ** sequences,char * structure,plist ** pl)199 alipf_circ_fold(const char **sequences,
200 char *structure,
201 plist **pl)
202 {
203 return wrap_alipf_fold(sequences, structure, pl, NULL, do_backtrack, fold_constrained, 1);
204 }
205
206
207 PUBLIC float
alipf_fold_par(const char ** sequences,char * structure,plist ** pl,vrna_exp_param_t * parameters,int calculate_bppm,int is_constrained,int is_circular)208 alipf_fold_par(const char **sequences,
209 char *structure,
210 plist **pl,
211 vrna_exp_param_t *parameters,
212 int calculate_bppm,
213 int is_constrained,
214 int is_circular)
215 {
216 return wrap_alipf_fold(sequences,
217 structure,
218 pl,
219 parameters,
220 calculate_bppm,
221 is_constrained,
222 is_circular);
223 }
224
225
226 PUBLIC FLT_OR_DBL *
alipf_export_bppm(void)227 alipf_export_bppm(void)
228 {
229 if (backward_compat_compound)
230 if (backward_compat_compound->exp_matrices)
231 if (backward_compat_compound->exp_matrices->probs)
232 return backward_compat_compound->exp_matrices->probs;
233
234 return NULL;
235 }
236
237
238 PUBLIC FLT_OR_DBL *
export_ali_bppm(void)239 export_ali_bppm(void)
240 {
241 if (backward_compat_compound)
242 if (backward_compat_compound->exp_matrices)
243 if (backward_compat_compound->exp_matrices->probs)
244 return backward_compat_compound->exp_matrices->probs;
245
246 return NULL;
247 }
248
249
250 /*brauch ma nurnoch pscores!*/
251 PUBLIC char *
alipbacktrack(double * prob)252 alipbacktrack(double *prob)
253 {
254 if (backward_compat_compound) {
255 if (backward_compat_compound->exp_matrices) {
256 vrna_exp_param_t *params = backward_compat_compound->exp_params;
257 int n = backward_compat_compound->length;
258 int n_seq = backward_compat_compound->n_seq;
259 int *idx = backward_compat_compound->iindx;
260 double Q = (double)backward_compat_compound->exp_matrices->q[idx[1] - n];
261 char *s = vrna_pbacktrack(backward_compat_compound);
262 double e = (double)vrna_eval_structure(backward_compat_compound, s);
263 e -= (double)vrna_eval_covar_structure(backward_compat_compound, s);
264 double fe = (-log(Q) - n * log(params->pf_scale)) * params->kT / (1000.0 * n_seq);
265 *prob = exp((fe - e) / params->kT);
266 return s;
267 }
268 }
269
270 return NULL;
271 }
272
273
274 /*-------------------------------------------------------------------------*/
275 /* make arrays used for alipf_fold available to other routines */
276 PUBLIC int
get_alipf_arrays(short *** S_p,short *** S5_p,short *** S3_p,unsigned short *** a2s_p,char *** Ss_p,FLT_OR_DBL ** qb_p,FLT_OR_DBL ** qm_p,FLT_OR_DBL ** q1k_p,FLT_OR_DBL ** qln_p,short ** pscore_p)277 get_alipf_arrays(short ***S_p,
278 short ***S5_p,
279 short ***S3_p,
280 unsigned short ***a2s_p,
281 char ***Ss_p,
282 FLT_OR_DBL **qb_p,
283 FLT_OR_DBL **qm_p,
284 FLT_OR_DBL **q1k_p,
285 FLT_OR_DBL **qln_p,
286 short **pscore_p)
287 {
288 if (backward_compat_compound) {
289 if (backward_compat_compound->exp_matrices) {
290 if (backward_compat_compound->exp_matrices->qb) {
291 *S_p = backward_compat_compound->S;
292 *S5_p = backward_compat_compound->S5;
293 *S3_p = backward_compat_compound->S3;
294 *Ss_p = backward_compat_compound->Ss;
295 *qb_p = backward_compat_compound->exp_matrices->qb;
296 *qm_p = backward_compat_compound->exp_matrices->qm;
297 *q1k_p = backward_compat_compound->exp_matrices->q1k;
298 *qln_p = backward_compat_compound->exp_matrices->qln;
299 *pscore_p = backward_compat_compound->pscore_pf_compat;
300 *a2s_p = backward_compat_a2s;
301 return 1;
302 }
303 }
304 }
305
306 return 0;
307 }
308
309
310 PUBLIC void
free_alipf_arrays(void)311 free_alipf_arrays(void)
312 {
313 if (backward_compat_compound && backward_compat) {
314 vrna_fold_compound_free(backward_compat_compound);
315 backward_compat_compound = NULL;
316 backward_compat = 0;
317 iindx = NULL;
318 }
319 }
320
321 #endif
322