1 #ifdef HAVE_CONFIG_H
2 #include "config.h"
3 #endif
4
5 #include <stdio.h>
6 #include <stdlib.h>
7
8 #include <svm.h>
9
10 #ifdef __GNUC__
11 # define INLINE inline
12 #else
13 # define INLINE
14 #endif
15
16 #include "ViennaRNA/utils/basic.h"
17 #include "ViennaRNA/utils/svm.h"
18
19 #include "ViennaRNA/zscore_dat.inc"
20
21
22 PRIVATE INLINE double
23 get_zscore(vrna_fold_compound_t *fc,
24 int i,
25 int j,
26 int e);
27
28
29 PUBLIC int
vrna_zsc_filter_init(vrna_fold_compound_t * fc,double min_z,unsigned int options)30 vrna_zsc_filter_init(vrna_fold_compound_t *fc,
31 double min_z,
32 unsigned int options)
33 {
34 if (fc) {
35 vrna_zsc_filter_free(fc);
36
37 fc->zscore_data = (vrna_zsc_dat_t)vrna_alloc(sizeof(struct vrna_zsc_dat_s));
38 fc->zscore_data->filter_on = (options & VRNA_ZSCORE_FILTER_ON) ? 1 : 0;
39 fc->zscore_data->pre_filter = (options & VRNA_ZSCORE_PRE_FILTER) ? 1 : 0;
40 fc->zscore_data->report_subsumed = (options & VRNA_ZSCORE_REPORT_SUBSUMED) ? 1 : 0;
41 fc->zscore_data->min_z = min_z;
42 fc->zscore_data->avg_model = svm_load_model_string(avg_model_string);
43 fc->zscore_data->sd_model = svm_load_model_string(sd_model_string);
44
45 if (fc->zscore_data->pre_filter)
46 fc->zscore_data->current_z = (double *)vrna_alloc(sizeof(double) * (fc->window_size + 2));
47 else
48 fc->zscore_data->current_z = NULL;
49
50 fc->zscore_data->current_i = 0;
51
52 return 1;
53 }
54
55 return 0;
56 }
57
58
59 PUBLIC int
vrna_zsc_filter_update(vrna_fold_compound_t * fc,double min_z,unsigned int options)60 vrna_zsc_filter_update(vrna_fold_compound_t *fc,
61 double min_z,
62 unsigned int options)
63 {
64 if (fc) {
65 if (!fc->zscore_data) {
66 vrna_zsc_filter_init(fc, min_z, options);
67 } else {
68 fc->zscore_data->min_z = min_z;
69
70 if (!(options & VRNA_ZSCORE_OPTIONS_NONE)) {
71 fc->zscore_data->filter_on = (options & VRNA_ZSCORE_FILTER_ON) ? 1 : 0;
72 fc->zscore_data->pre_filter = (options & VRNA_ZSCORE_PRE_FILTER) ? 1 : 0;
73 fc->zscore_data->report_subsumed = (options & VRNA_ZSCORE_REPORT_SUBSUMED) ? 1 : 0;
74 }
75
76 if (fc->zscore_data->pre_filter) {
77 if (fc->zscore_data->current_z) {
78 fc->zscore_data->current_z += fc->zscore_data->current_i;
79 free(fc->zscore_data->current_z);
80 }
81
82 fc->zscore_data->current_z = (double *)vrna_alloc(sizeof(double) * (fc->window_size + 2));
83 fc->zscore_data->current_i = 0;
84 } else if (fc->zscore_data->current_z) {
85 fc->zscore_data->current_z += fc->zscore_data->current_i;
86 free(fc->zscore_data->current_z);
87 fc->zscore_data->current_z = NULL;
88 fc->zscore_data->current_i = 0;
89 }
90 }
91
92 return 1;
93 }
94
95 return 0;
96 }
97
98
99 PUBLIC void
vrna_zsc_filter_free(vrna_fold_compound_t * fc)100 vrna_zsc_filter_free(vrna_fold_compound_t *fc)
101 {
102 if ((fc) && (fc->zscore_data)) {
103 vrna_zsc_dat_t zsc_data = fc->zscore_data;
104
105 zsc_data->current_z += zsc_data->current_i;
106 free(zsc_data->current_z);
107 svm_free_model_content(zsc_data->avg_model);
108 svm_free_model_content(zsc_data->sd_model);
109 free(zsc_data);
110
111 fc->zscore_data = NULL;
112 }
113 }
114
115
116 PUBLIC int
vrna_zsc_filter_on(vrna_fold_compound_t * fc)117 vrna_zsc_filter_on(vrna_fold_compound_t *fc)
118 {
119 if ((fc) &&
120 (fc->zscore_data) &&
121 (fc->zscore_data->filter_on))
122 return 1;
123
124 return 0;
125 }
126
127
128 PUBLIC double
vrna_zsc_filter_threshold(vrna_fold_compound_t * fc)129 vrna_zsc_filter_threshold(vrna_fold_compound_t *fc)
130 {
131 if ((fc) &&
132 (fc->zscore_data) &&
133 (fc->zscore_data->filter_on))
134 return fc->zscore_data->min_z;
135
136 return (double)INF;
137 }
138
139
140 PUBLIC double
vrna_zsc_compute(vrna_fold_compound_t * fc,unsigned int i,unsigned int j,int e)141 vrna_zsc_compute(vrna_fold_compound_t *fc,
142 unsigned int i,
143 unsigned int j,
144 int e)
145 {
146 if ((fc) &&
147 (fc->zscore_data) &&
148 (fc->zscore_data->filter_on))
149 return get_zscore(fc, i, j, e);
150
151 return (double)INF;
152 }
153
154
155 PRIVATE INLINE double
get_zscore(vrna_fold_compound_t * fc,int i,int j,int e)156 get_zscore(vrna_fold_compound_t *fc,
157 int i,
158 int j,
159 int e)
160 {
161 short *S;
162 int info_avg, start, end, dangle_model, length;
163 double average_free_energy;
164 double sd_free_energy;
165 double z;
166 vrna_zsc_dat_t d;
167
168 length = fc->length;
169 S = fc->sequence_encoding2;
170 dangle_model = fc->params->model_details.dangles;
171 z = (double)INF;
172 d = fc->zscore_data;
173
174 start = (dangle_model) ? MAX2(1, i - 1) : i;
175 end = (dangle_model) ? MIN2(length, j + 1) : j;
176
177 int *AUGC = get_seq_composition(S, start, end, length);
178
179 /*\svm*/
180 average_free_energy = avg_regression(AUGC[0],
181 AUGC[1],
182 AUGC[2],
183 AUGC[3],
184 AUGC[4],
185 d->avg_model,
186 &info_avg);
187
188 if (info_avg == 0) {
189 double min_sd = minimal_sd(AUGC[0], AUGC[1], AUGC[2], AUGC[3], AUGC[4]);
190 double difference = ((double)e / 100.) - average_free_energy;
191
192 if (difference - (d->min_z * min_sd) <= 0.0001) {
193 sd_free_energy = sd_regression(AUGC[0], AUGC[1], AUGC[2], AUGC[3], AUGC[4], d->sd_model);
194 z = difference / sd_free_energy;
195 }
196 }
197
198 free(AUGC);
199
200 return z;
201 }
202