1 /*
2 * minimum free energy folding
3 * for a set of aligned sequences
4 *
5 * c Ivo Hofacker
6 *
7 * Vienna RNA package
8 */
9
10 /**
11 *** \file alifold.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 <math.h>
23 #include <ctype.h>
24 #include <string.h>
25 #include <limits.h>
26
27 #include "ViennaRNA/fold_vars.h"
28 #include "ViennaRNA/datastructures/basic.h"
29 #include "ViennaRNA/mfe.h"
30 #include "ViennaRNA/fold.h"
31 #include "ViennaRNA/eval.h"
32 #include "ViennaRNA/utils/basic.h"
33 #include "ViennaRNA/params/default.h"
34 #include "ViennaRNA/params/basic.h"
35 #include "ViennaRNA/ribo.h"
36 #include "ViennaRNA/gquad.h"
37 #include "ViennaRNA/alifold.h"
38 #include "ViennaRNA/utils/alignments.h"
39 #include "ViennaRNA/loops/all.h"
40
41 #ifdef _OPENMP
42 #include <omp.h>
43 #endif
44
45 /*
46 #################################
47 # GLOBAL VARIABLES #
48 #################################
49 */
50
51 /*
52 #################################
53 # PRIVATE VARIABLES #
54 #################################
55 */
56
57
58 #define MAXSECTORS 500 /* dimension for a backtrack array */
59
60 /* some backward compatibility stuff */
61 PRIVATE vrna_fold_compound_t *backward_compat_compound = NULL;
62 PRIVATE int backward_compat = 0;
63
64 #ifdef _OPENMP
65
66 #pragma omp threadprivate(backward_compat_compound, backward_compat)
67
68 #endif
69
70 /*
71 #################################
72 # PRIVATE FUNCTION DECLARATIONS #
73 #################################
74 */
75
76 PRIVATE float
77 wrap_alifold(const char **strings,
78 char *structure,
79 vrna_param_t *parameters,
80 int is_constrained,
81 int is_circular);
82
83
84 /*
85 #################################
86 # BEGIN OF FUNCTION DEFINITIONS #
87 #################################
88 */
89
90 /*###########################################*/
91 /*# deprecated functions below #*/
92 /*###########################################*/
93
94 PRIVATE float
wrap_alifold(const char ** strings,char * structure,vrna_param_t * parameters,int is_constrained,int is_circular)95 wrap_alifold(const char **strings,
96 char *structure,
97 vrna_param_t *parameters,
98 int is_constrained,
99 int is_circular)
100 {
101 vrna_fold_compound_t *vc;
102 vrna_param_t *P;
103 float mfe;
104
105 #ifdef _OPENMP
106 /* Explicitly turn off dynamic threads */
107 omp_set_dynamic(0);
108 #endif
109
110 /* we need the parameter structure for hard constraints */
111 if (parameters) {
112 P = vrna_params_copy(parameters);
113 } else {
114 vrna_md_t md;
115 set_model_details(&md);
116 md.temperature = temperature;
117 P = vrna_params(&md);
118 }
119
120 P->model_details.circ = is_circular;
121
122 vc = vrna_fold_compound_comparative(strings, &(P->model_details), VRNA_OPTION_DEFAULT);
123
124 if (parameters) {
125 /* replace params if necessary */
126 free(vc->params);
127 vc->params = P;
128 } else {
129 free(P);
130 }
131
132 /* handle hard constraints in pseudo dot-bracket format if passed via simple interface */
133 if (is_constrained && structure)
134 vrna_constraints_add(vc, (const char *)structure, VRNA_CONSTRAINT_DB_DEFAULT);
135
136 if (backward_compat_compound && backward_compat)
137 vrna_fold_compound_free(backward_compat_compound);
138
139 backward_compat_compound = vc;
140 backward_compat = 1;
141
142 /* call mfe() function without backtracking */
143 mfe = vrna_mfe(vc, NULL);
144
145 /* backtrack structure */
146 if (structure && vc->params->model_details.backtrack) {
147 char *ss;
148 int length;
149 sect bt_stack[MAXSECTORS];
150 vrna_bp_stack_t *bp;
151
152 length = vc->length;
153 bp = (vrna_bp_stack_t *)vrna_alloc(sizeof(vrna_bp_stack_t) * (4 * (1 + length / 2))); /* add a guess of how many G's may be involved in a G quadruplex */
154
155 vrna_backtrack_from_intervals(vc, bp, bt_stack, 0);
156
157 ss = vrna_db_from_bp_stack(bp, length);
158 strncpy(structure, ss, length + 1);
159 free(ss);
160
161 if (base_pair)
162 free(base_pair);
163
164 base_pair = bp;
165 }
166
167 return mfe;
168 }
169
170
171 PUBLIC void
free_alifold_arrays(void)172 free_alifold_arrays(void)
173 {
174 if (backward_compat_compound && backward_compat) {
175 vrna_fold_compound_free(backward_compat_compound);
176 backward_compat_compound = NULL;
177 backward_compat = 0;
178 }
179 }
180
181
182 PUBLIC float
alifold(const char ** strings,char * structure)183 alifold(const char **strings,
184 char *structure)
185 {
186 return wrap_alifold(strings, structure, NULL, fold_constrained, 0);
187 }
188
189
190 PUBLIC float
circalifold(const char ** strings,char * structure)191 circalifold(const char **strings,
192 char *structure)
193 {
194 return wrap_alifold(strings, structure, NULL, fold_constrained, 1);
195 }
196
197
198 PUBLIC void
update_alifold_params(void)199 update_alifold_params(void)
200 {
201 vrna_fold_compound_t *v;
202
203 if (backward_compat_compound && backward_compat) {
204 v = backward_compat_compound;
205
206 if (v->params)
207 free(v->params);
208
209 vrna_md_t md;
210 set_model_details(&md);
211 v->params = vrna_params(&md);
212 }
213 }
214
215
216 PUBLIC float
energy_of_ali_gquad_structure(const char ** sequences,const char * structure,int n_seq,float * energy)217 energy_of_ali_gquad_structure(const char **sequences,
218 const char *structure,
219 int n_seq,
220 float *energy)
221 {
222 if (sequences[0] != NULL) {
223 vrna_fold_compound_t *vc;
224 vrna_md_t md;
225
226 set_model_details(&md);
227 md.gquad = 1;
228
229 vc = vrna_fold_compound_comparative(sequences, &md, VRNA_OPTION_EVAL_ONLY);
230
231 energy[0] = vrna_eval_structure(vc, structure);
232 energy[1] = vrna_eval_covar_structure(vc, structure);
233
234 vrna_fold_compound_free(vc);
235 } else {
236 vrna_message_warning("energy_of_ali_gquad_structure: "
237 "no sequences in alignment!");
238 return (float)(INF / 100.);
239 }
240
241 return energy[0];
242 }
243
244
245 PUBLIC float
energy_of_alistruct(const char ** sequences,const char * structure,int n_seq,float * energy)246 energy_of_alistruct(const char **sequences,
247 const char *structure,
248 int n_seq,
249 float *energy)
250 {
251 if (sequences[0] != NULL) {
252 vrna_fold_compound_t *vc;
253
254 vrna_md_t md;
255 set_model_details(&md);
256
257 vc = vrna_fold_compound_comparative(sequences, &md, VRNA_OPTION_EVAL_ONLY);
258
259 energy[0] = vrna_eval_structure(vc, structure);
260 energy[1] = vrna_eval_covar_structure(vc, structure);
261
262 vrna_fold_compound_free(vc);
263 } else {
264 vrna_message_warning("energy_of_alistruct(): "
265 "no sequences in alignment!");
266 return (float)(INF / 100.);
267 }
268
269 return energy[0];
270 }
271
272
273 #endif
274