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