1 #ifndef VIENNA_RNA_PACKAGE_PAIR_MAT_H
2 #define VIENNA_RNA_PACKAGE_PAIR_MAT_H
3 
4 #include <ctype.h>
5 #include <ViennaRNA/utils/basic.h>
6 #include <ViennaRNA/fold_vars.h>
7 
8 #define NBASES 8
9 /*@notnull@*/
10 
11 #ifndef INLINE
12 # ifdef __GNUC__
13 #  define INLINE inline
14 # else
15 #  define INLINE
16 # endif
17 #endif
18 
19 static const char Law_and_Order[]         = "_ACGUTXKI";
20 static int        BP_pair[NBASES][NBASES] =
21   /* _  A  C  G  U  X  K  I */
22 { { 0, 0, 0, 0, 0, 0, 0, 0 },
23   { 0, 0, 0, 0, 5, 0, 0, 5 },
24   { 0, 0, 0, 1, 0, 0, 0, 0 },
25   { 0, 0, 2, 0, 3, 0, 0, 0 },
26   { 0, 6, 0, 4, 0, 0, 0, 6 },
27   { 0, 0, 0, 0, 0, 0, 2, 0 },
28   { 0, 0, 0, 0, 0, 1, 0, 0 },
29   { 0, 6, 0, 0, 5, 0, 0, 0 } };
30 
31 #define MAXALPHA 20       /* maximal length of alphabet */
32 
33 static short  alias[MAXALPHA + 1];
34 static int    pair[MAXALPHA + 1][MAXALPHA + 1];
35 /* rtype[pair[i][j]]:=pair[j][i] */
36 static int    rtype[8] = {
37   0, 2, 1, 4, 3, 6, 5, 7
38 };
39 
40 #ifdef _OPENMP
41 #pragma omp threadprivate(Law_and_Order, BP_pair, alias, pair, rtype)
42 #endif
43 
44 /* for backward compatibility */
45 #define ENCODE(c) encode_char(c)
46 
47 static INLINE int
encode_char(char c)48 encode_char(char c)
49 {
50   /* return numerical representation of base used e.g. in pair[][] */
51   int code;
52 
53   c = toupper(c);
54 
55   if (energy_set > 0) {
56     code = (int)(c - 'A') + 1;
57   } else {
58     const char *pos;
59     pos = strchr(Law_and_Order, c);
60     if (pos == NULL)
61       code = 0;
62     else
63       code = (int)(pos - Law_and_Order);
64 
65     if (code > 5)
66       code = 0;
67 
68     if (code > 4)
69       code--;           /* make T and U equivalent */
70   }
71 
72   return code;
73 }
74 
75 
76 /*@+boolint +charint@*/
77 /*@null@*/
78 extern char *nonstandards;
79 
80 static INLINE void
make_pair_matrix(void)81 make_pair_matrix(void)
82 {
83   int i, j;
84 
85   if (energy_set == 0) {
86     for (i = 0; i < 5; i++)
87       alias[i] = (short)i;
88     alias[5]  = 3;  /* X <-> G */
89     alias[6]  = 2;  /* K <-> C */
90     alias[7]  = 0;  /* I <-> default base '@' */
91     for (i = 0; i < NBASES; i++)
92       for (j = 0; j < NBASES; j++)
93         pair[i][j] = BP_pair[i][j];
94     if (noGU)
95       pair[3][4] = pair[4][3] = 0;
96 
97     if (nonstandards != NULL) {
98       /* allow nonstandard bp's */
99       for (i = 0; i < (int)strlen(nonstandards); i += 2)
100         pair[encode_char(nonstandards[i])]
101         [encode_char(nonstandards[i + 1])] = 7;
102     }
103 
104     for (i = 0; i < NBASES; i++)
105       for (j = 0; j < NBASES; j++)
106         rtype[pair[i][j]] = pair[j][i];
107   } else {
108     for (i = 0; i <= MAXALPHA; i++)
109       for (j = 0; j <= MAXALPHA; j++)
110         pair[i][j] = 0;
111     if (energy_set == 1) {
112       for (i = 1; i < MAXALPHA; ) {
113         alias[i++]  = 3;      /* A <-> G */
114         alias[i++]  = 2;      /* B <-> C */
115       }
116       for (i = 1; i < MAXALPHA; i++) {
117         pair[i][i + 1] = 2;       /* AB <-> GC */
118         i++;
119         pair[i][i - 1] = 1;       /* BA <-> CG */
120       }
121     } else if (energy_set == 2) {
122       for (i = 1; i < MAXALPHA; ) {
123         alias[i++]  = 1;      /* A <-> A*/
124         alias[i++]  = 4;      /* B <-> U */
125       }
126       for (i = 1; i < MAXALPHA; i++) {
127         pair[i][i + 1] = 5;       /* AB <-> AU */
128         i++;
129         pair[i][i - 1] = 6;       /* BA <-> UA */
130       }
131     } else if (energy_set == 3) {
132       for (i = 1; i < MAXALPHA - 2; ) {
133         alias[i++]  = 3;    /* A <-> G */
134         alias[i++]  = 2;    /* B <-> C */
135         alias[i++]  = 1;    /* C <-> A */
136         alias[i++]  = 4;    /* D <-> U */
137       }
138       for (i = 1; i < MAXALPHA - 2; i++) {
139         pair[i][i + 1] = 2;     /* AB <-> GC */
140         i++;
141         pair[i][i - 1] = 1;     /* BA <-> CG */
142         i++;
143         pair[i][i + 1] = 5;     /* CD <-> AU */
144         i++;
145         pair[i][i - 1] = 6;     /* DC <-> UA */
146       }
147     } else {
148       vrna_message_error("What energy_set are YOU using??");
149     }
150 
151     for (i = 0; i <= MAXALPHA; i++)
152       for (j = 0; j <= MAXALPHA; j++)
153         rtype[pair[i][j]] = pair[j][i];
154   }
155 }
156 
157 
158 static INLINE short *
encode_sequence(const char * sequence,short how)159 encode_sequence(const char  *sequence,
160                 short       how)
161 {
162   unsigned int  i, l = (unsigned int)strlen(sequence);
163   short         *S = (short *)vrna_alloc(sizeof(short) * (l + 2));
164 
165   switch (how) {
166     /* standard encoding as always used for S */
167     case 0:
168       for (i = 1; i <= l; i++)    /* make numerical encoding of sequence */
169         S[i] = (short)encode_char(sequence[i - 1]);
170       S[l + 1]  = S[1];
171       S[0]      = (short)l;
172       break;
173     /* encoding for mismatches of nostandard bases (normally used for S1) */
174     case 1:
175       for (i = 1; i <= l; i++)
176         S[i] = alias[(short)encode_char(sequence[i - 1])];
177       S[l + 1]  = S[1];
178       S[0]      = S[l];
179       break;
180   }
181 
182   return S;
183 }
184 
185 
186 #endif /* VIENNA_RNA_PACKAGE_PAIR_MAT_H */
187