1 /*
2 * Utitlities for plot functions
3 *
4 * c Ronny Lorenz
5 * ViennaRNA package
6 */
7
8 #ifdef HAVE_CONFIG_H
9 #include "config.h"
10 #endif
11
12 #include <stdio.h>
13 #include <stdlib.h>
14 #include <string.h>
15 #include <ctype.h>
16 #include "ViennaRNA/utils/basic.h"
17 #include "ViennaRNA/utils/structures.h"
18 #include "ViennaRNA/alphabet.h"
19 #include "ViennaRNA/plotting/utils.h"
20
21 /*
22 #################################
23 # PRIVATE MACROS #
24 #################################
25 */
26
27 /*
28 #################################
29 # GLOBAL VARIABLES #
30 #################################
31 */
32
33 /*
34 #################################
35 # PRIVATE VARIABLES #
36 #################################
37 */
38
39 /*
40 #################################
41 # PRIVATE FUNCTION DECLARATIONS #
42 #################################
43 */
44
45 /*
46 #################################
47 # BEGIN OF FUNCTION DEFINITIONS #
48 #################################
49 */
50 PUBLIC char **
vrna_annotate_covar_db(const char ** alignment,const char * structure,vrna_md_t * md_p)51 vrna_annotate_covar_db(const char **alignment,
52 const char *structure,
53 vrna_md_t *md_p)
54 {
55 return vrna_annotate_covar_db_extended(alignment,
56 structure,
57 md_p,
58 VRNA_BRACKETS_RND);
59 }
60
61
62 PUBLIC char **
vrna_annotate_covar_db_extended(const char ** alignment,const char * structure,vrna_md_t * md_p,unsigned int options)63 vrna_annotate_covar_db_extended(const char **alignment,
64 const char *structure,
65 vrna_md_t *md_p,
66 unsigned int options)
67 {
68 /* produce annotation for colored drawings from vrna_file_PS_rnaplot_a() */
69 char *ps, *colorps, **A;
70 int i, n, s, pairings, maxl, a, b;
71 short *ptable;
72 char *colorMatrix[6][3] = {
73 { "0.0 1", "0.0 0.6", "0.0 0.2" }, /* red */
74 { "0.16 1", "0.16 0.6", "0.16 0.2" }, /* ochre */
75 { "0.32 1", "0.32 0.6", "0.32 0.2" }, /* turquoise */
76 { "0.48 1", "0.48 0.6", "0.48 0.2" }, /* green */
77 { "0.65 1", "0.65 0.6", "0.65 0.2" }, /* blue */
78 { "0.81 1", "0.81 0.6", "0.81 0.2" } /* violet */
79 };
80 vrna_md_t md;
81
82 if ((!alignment) || (!structure))
83 return NULL;
84
85 if (md_p)
86 vrna_md_copy(&md, md_p);
87 else
88 vrna_md_set_default(&md);
89
90 n = strlen(alignment[0]);
91 maxl = 1024;
92
93 A = (char **)vrna_alloc(sizeof(char *) * 2);
94 ps = (char *)vrna_alloc(maxl);
95 colorps = (char *)vrna_alloc(maxl);
96
97 ptable = vrna_ptable_from_string(structure, options);
98
99 for (i = 1; i <= n; i++) {
100 char pps[64], ci = '\0', cj = '\0';
101 int j, type, pfreq[8] = {
102 0, 0, 0, 0, 0, 0, 0, 0
103 }, vi = 0, vj = 0;
104
105 if ((j = ptable[i]) < i)
106 continue;
107
108 for (s = 0; alignment[s] != NULL; s++) {
109 a = vrna_nucleotide_encode(alignment[s][i - 1], &md);
110 b = vrna_nucleotide_encode(alignment[s][j - 1], &md);
111 type = md.pair[a][b];
112 pfreq[type]++;
113 if (type) {
114 if (alignment[s][i - 1] != ci) {
115 ci = alignment[s][i - 1];
116 vi++;
117 }
118
119 if (alignment[s][j - 1] != cj) {
120 cj = alignment[s][j - 1];
121 vj++;
122 }
123 }
124 }
125 for (pairings = 0, s = 1; s <= 7; s++)
126 if (pfreq[s])
127 pairings++;
128
129 if ((maxl - strlen(ps) < 192) || ((maxl - strlen(colorps)) < 64)) {
130 maxl *= 2;
131 ps = (char *)vrna_realloc(ps, sizeof(char) * maxl);
132 colorps = (char *)vrna_realloc(colorps, sizeof(char) * maxl);
133 if ((ps == NULL) || (colorps == NULL))
134 vrna_message_error("out of memory in realloc");
135 }
136
137 if (pfreq[0] <= 2 && pairings > 0) {
138 snprintf(pps, 64, "%d %d %s colorpair\n",
139 i, j, colorMatrix[pairings - 1][pfreq[0]]);
140 strcat(colorps, pps);
141 }
142
143 if (pfreq[0] > 0) {
144 snprintf(pps, 64, "%d %d %d gmark\n", i, j, pfreq[0]);
145 strcat(ps, pps);
146 }
147
148 if (vi > 1) {
149 snprintf(pps, 64, "%d cmark\n", i);
150 strcat(ps, pps);
151 }
152
153 if (vj > 1) {
154 snprintf(pps, 64, "%d cmark\n", j);
155 strcat(ps, pps);
156 }
157 }
158 free(ptable);
159 A[0] = colorps;
160 A[1] = ps;
161 return A;
162 }
163
164
165 /* produce info for PS_color_dot_plot */
166 PUBLIC vrna_cpair_t *
vrna_annotate_covar_pairs(const char ** alignment,vrna_ep_t * pl,vrna_ep_t * mfel,double threshold,vrna_md_t * md_p)167 vrna_annotate_covar_pairs(const char **alignment,
168 vrna_ep_t *pl,
169 vrna_ep_t *mfel,
170 double threshold,
171 vrna_md_t *md_p)
172 {
173 unsigned int n_seq;
174 int i, n, s, a, b, z, j, c, pfreq[7];
175 vrna_cpair_t *cp;
176 vrna_ep_t *ptr;
177 vrna_md_t md;
178
179 if ((!alignment) || (!pl))
180 return NULL;
181
182 if (md_p)
183 vrna_md_copy(&md, md_p);
184 else
185 vrna_md_set_default(&md);
186
187 /* count number of sequences */
188 for (n_seq = 0; alignment[n_seq] != NULL; n_seq++);
189
190 /* count number of entries in pl */
191 for (n = 0; pl[n].i > 0; n++);
192
193 c = 0;
194 cp = (vrna_cpair_t *)vrna_alloc(sizeof(vrna_cpair_t) * (n + 1));
195
196 for (i = 0; i < n; i++) {
197 int ncomp = 0;
198 if (pl[i].p > threshold) {
199 cp[c].i = pl[i].i;
200 cp[c].j = pl[i].j;
201 cp[c].p = pl[i].p;
202 cp[c].type = pl[i].type;
203 for (z = 0; z < 7; z++)
204 pfreq[z] = 0;
205 for (s = 0; s < n_seq; s++) {
206 a = vrna_nucleotide_encode(alignment[s][cp[c].i - 1], &md);
207 b = vrna_nucleotide_encode(alignment[s][cp[c].j - 1], &md);
208 if ((alignment[s][cp[c].j - 1] == '~') || (alignment[s][cp[c].i - 1] == '~'))
209 continue;
210
211 if ((md.gquad) && (a == 3) && (b == 3))
212 continue;
213
214 pfreq[md.pair[a][b]]++;
215 }
216 for (z = 1; z < 7; z++)
217 if (pfreq[z] > 0)
218 ncomp++;
219
220 cp[c].hue = MAX2(0.0, (ncomp - 1.0) / 6.2); /* hue<6/6.9 (hue=1 == hue=0) */
221 cp[c].sat = 1 - MIN2(1.0, (float)(pfreq[0] * 2. /*pi[i].bp[0]*/ / (n_seq)));
222 c++;
223 }
224 }
225
226 /*
227 * check whether all MFE base pairs are present in above list, and
228 * insert them if they are missing
229 */
230 if (mfel) {
231 for (ptr = mfel; ptr->i > 0; ptr++) {
232 int nofound = 1;
233 for (j = 0; j < c; j++) {
234 if ((cp[j].i == ptr->i) && (cp[j].j == ptr->j)) {
235 cp[j].mfe = 1;
236 nofound = 0;
237 break;
238 }
239 }
240 if (nofound) {
241 vrna_message_warning("mfe base pair with very low prob in pf: %d %d",
242 ptr->i,
243 ptr->j);
244
245 cp = (vrna_cpair_t *)vrna_realloc(cp, sizeof(vrna_cpair_t) * (c + 2));
246 cp[c].i = ptr->i;
247 cp[c].j = ptr->j;
248 cp[c].p = 0.;
249 cp[c].type = VRNA_PLIST_TYPE_BASEPAIR;
250 cp[c].hue = 0;
251 cp[c].sat = 0;
252 cp[c].mfe = 1;
253 c++;
254 cp[c].i = cp[c].j = 0;
255 }
256 }
257 }
258
259 return cp;
260 }
261