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