1 /*
2 Copyright (c) 1996,1997,1998,1999,2000,2001,2004,2006,2007
3 Whitehead Institute for Biomedical Research, Steve Rozen
4 (http://purl.com/STEVEROZEN/), Andreas Untergasser and Helen Skaletsky
5 All rights reserved.
6
7 This file is part of the oligotm library.
8
9 The oligotm library is free software; you can redistribute it and/or modify
10 it under the terms of the GNU General Public License as published by
11 the Free Software Foundation; either version 2 of the License, or
12 (at your option) any later version.
13
14 The oligotm library is distributed in the hope that it will be useful,
15 but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 GNU General Public License for more details.
18
19 You should have received a copy of the GNU General Public License
20 along with the oligtm library (file gpl-2.0.txt in the source
21 distribution); if not, write to the Free Software
22 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
23
24 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
25 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
26 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
27 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
28 OWNERS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
29 SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
30 LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
31 DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
32 THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
33 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
34 OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
35
36 */
37
38 #include <limits.h>
39 #include <math.h>
40 #include <string.h>
41 #include "oligotm.h"
42
43 #define T_KELVIN 273.15;
44 #define A_CHAR 'A'
45 #define G_CHAR 'G'
46 #define T_CHAR 'T'
47 #define C_CHAR 'C'
48 #define N_CHAR 'N'
49
50 #define CATID5(A,B,C,D,E) A##B##C##D##E
51 #define CATID2(A,B) A##B
52
53 #define DO_PAIR(LAST,THIS) \
54 if (CATID2(THIS,_CHAR) == c) { \
55 dh += CATID5(H,_,LAST,_,THIS); \
56 ds += CATID5(S,_,LAST,_,THIS); \
57 goto CATID2(THIS,_STATE); \
58 }
59
60 #define STATE(LAST) \
61 CATID2(LAST,_STATE): \
62 c = *s; s++; \
63 DO_PAIR(LAST,A) \
64 else DO_PAIR(LAST,T) \
65 else DO_PAIR(LAST,G) \
66 else DO_PAIR(LAST,C) \
67 else DO_PAIR(LAST,N) \
68 else if ('\0' == c) \
69 goto DONE; \
70 else goto ERROR \
71
72 #define DO_PAIR2(LAST,THIS) \
73 if (CATID2(THIS,_CHAR) == c) { \
74 dh += CATID5(DH,_,LAST,_,THIS); \
75 ds += CATID5(DS,_,LAST,_,THIS); \
76 goto CATID2(THIS,_STATE2); \
77 }
78 #define STATE2(LAST) \
79 CATID2(LAST,_STATE2): \
80 c = *s; s++; \
81 DO_PAIR2(LAST,A) \
82 else DO_PAIR2(LAST,T) \
83 else DO_PAIR2(LAST,G) \
84 else DO_PAIR2(LAST,C) \
85 else DO_PAIR2(LAST,N) \
86 else if ('\0' == c) \
87 goto DONE; \
88 else goto ERROR \
89
90 /*
91 * Two tables of nearest-neighbor parameters for di-nucleotide
92 * base pairs.
93 *
94 * These are included in this file because they are not needed by
95 * clients (callers) of oligtm().
96 */
97
98 /* Table 1 (old parameters):
99 * See table 2 in the paper [Breslauer KJ, Frank R, Bl�cker H and
100 * Marky LA (1986) "Predicting DNA duplex stability from the base
101 * sequence" Proc Natl Acad Sci 83:4746-50
102 * http://dx.doi.org/10.1073/pnas.83.11.3746]
103 */
104
105 #define S_A_A 240
106 #define S_A_C 173
107 #define S_A_G 208
108 #define S_A_T 239
109 #define S_A_N 215
110
111 #define S_C_A 129
112 #define S_C_C 266
113 #define S_C_G 278
114 #define S_C_T 208
115 #define S_C_N 220
116
117 #define S_G_A 135
118 #define S_G_C 267
119 #define S_G_G 266
120 #define S_G_T 173
121 #define S_G_N 210
122
123 #define S_T_A 169
124 #define S_T_C 135
125 #define S_T_G 129
126 #define S_T_T 240
127 #define S_T_N 168
128
129 #define S_N_A 168
130 #define S_N_C 210
131 #define S_N_G 220
132 #define S_N_T 215
133 #define S_N_N 203
134
135
136 #define H_A_A 91
137 #define H_A_C 65
138 #define H_A_G 78
139 #define H_A_T 86
140 #define H_A_N 80
141
142 #define H_C_A 58
143 #define H_C_C 110
144 #define H_C_G 119
145 #define H_C_T 78
146 #define H_C_N 91
147
148 #define H_G_A 56
149 #define H_G_C 111
150 #define H_G_G 110
151 #define H_G_T 65
152 #define H_G_N 85
153
154 #define H_T_A 60
155 #define H_T_C 56
156 #define H_T_G 58
157 #define H_T_T 91
158 #define H_T_N 66
159
160 #define H_N_A 66
161 #define H_N_C 85
162 #define H_N_G 91
163 #define H_N_T 80
164 #define H_N_N 80
165
166 /* Delta G's of disruption * 1000. */
167 #define G_A_A 1900
168 #define G_A_C 1300
169 #define G_A_G 1600
170 #define G_A_T 1500
171 #define G_A_N 1575
172
173 #define G_C_A 1900
174 #define G_C_C 3100
175 #define G_C_G 3600
176 #define G_C_T 1600
177 #define G_C_N 2550
178
179 #define G_G_A 1600
180 #define G_G_C 3100
181 #define G_G_G 3100
182 #define G_G_T 1300
183 #define G_G_N 2275
184
185 #define G_T_A 900
186 #define G_T_C 1600
187 #define G_T_G 1900
188 #define G_T_T 1900
189 #define G_T_N 1575
190
191 #define G_N_A 1575
192 #define G_N_C 2275
193 #define G_N_G 2550
194 #define G_N_T 1575
195 #define G_N_N 1994
196
197 /* Table 2, new parameters:
198 * Tables of nearest-neighbor thermodynamics for DNA bases, from the
199 * paper [SantaLucia JR (1998) "A unified view of polymer, dumbbell
200 * and oligonucleotide DNA nearest-neighbor thermodynamics", Proc Natl
201 * Acad Sci 95:1460-65 http://dx.doi.org/10.1073/pnas.95.4.1460]
202 */
203
204 #define DS_A_A 222
205 #define DS_A_C 224
206 #define DS_A_G 210
207 #define DS_A_T 204
208 #define DS_A_N 224
209
210 #define DS_C_A 227
211 #define DS_C_C 199
212 #define DS_C_G 272
213 #define DS_C_T 210
214 #define DS_C_N 272
215
216 #define DS_G_A 222
217 #define DS_G_C 244
218 #define DS_G_G 199
219 #define DS_G_T 224
220 #define DS_G_N 244
221
222 #define DS_T_A 213
223 #define DS_T_C 222
224 #define DS_T_G 227
225 #define DS_T_T 222
226 #define DS_T_N 227
227
228 #define DS_N_A 168
229 #define DS_N_C 210
230 #define DS_N_G 220
231 #define DS_N_T 215
232 #define DS_N_N 220
233
234
235 #define DH_A_A 79
236 #define DH_A_C 84
237 #define DH_A_G 78
238 #define DH_A_T 72
239 #define DH_A_N 72
240
241 #define DH_C_A 85
242 #define DH_C_C 80
243 #define DH_C_G 106
244 #define DH_C_T 78
245 #define DH_C_N 78
246
247 #define DH_G_A 82
248 #define DH_G_C 98
249 #define DH_G_G 80
250 #define DH_G_T 84
251 #define DH_G_N 80
252
253 #define DH_T_A 72
254 #define DH_T_C 82
255 #define DH_T_G 85
256 #define DH_T_T 79
257 #define DH_T_N 72
258
259 #define DH_N_A 72
260 #define DH_N_C 80
261 #define DH_N_G 78
262 #define DH_N_T 72
263 #define DH_N_N 72
264
265 /* Delta G's of disruption * 1000. */
266 #define DG_A_A 1000
267 #define DG_A_C 1440
268 #define DG_A_G 1280
269 #define DG_A_T 880
270 #define DG_A_N 880
271
272 #define DG_C_A 1450
273 #define DG_C_C 1840
274 #define DG_C_G 2170
275 #define DG_C_T 1280
276 #define DG_C_N 1450
277
278 #define DG_G_A 1300
279 #define DG_G_C 2240
280 #define DG_G_G 1840
281 #define DG_G_T 1440
282 #define DG_G_N 1300
283
284 #define DG_T_A 580
285 #define DG_T_C 1300
286 #define DG_T_G 1450
287 #define DG_T_T 1000
288 #define DG_T_N 580
289
290 #define DG_N_A 580
291 #define DG_N_C 1300
292 #define DG_N_G 1280
293 #define DG_N_T 880
294 #define DG_N_N 580
295
296 /* End of tables nearest-neighbor parameter. */
297
298
299 /* Calculate the melting temperature of oligo s. See
300 oligotm.h for documentation of arguments.
301 */
302 double
oligotm(const char * s,double DNA_nM,double K_mM,double divalent_conc,double dntp_conc,tm_method_type tm_method,salt_correction_type salt_corrections)303 oligotm(const char *s,
304 double DNA_nM,
305 double K_mM,
306 double divalent_conc,
307 double dntp_conc,
308 tm_method_type tm_method,
309 salt_correction_type salt_corrections)
310 {
311 register int dh = 0, ds = 0;
312 register char c;
313 double delta_H, delta_S;
314 double Tm; /* Melting temperature */
315 double correction;
316 int len, sym;
317 const char* d = s;
318 if(divalent_to_monovalent(divalent_conc, dntp_conc) == OLIGOTM_ERROR)
319 return OLIGOTM_ERROR;
320
321 /** K_mM = K_mM + divalent_to_monovalent(divalent_conc, dntp_conc); **/
322 if (tm_method != breslauer_auto
323 && tm_method != santalucia_auto)
324 return OLIGOTM_ERROR;
325 if (salt_corrections != schildkraut
326 && salt_corrections != santalucia
327 && salt_corrections != owczarzy)
328 return OLIGOTM_ERROR;
329
330 len = (strlen(s)-1);
331
332 sym = symmetry(s); /*Add symmetry correction if seq is symmetrical*/
333 if( tm_method == breslauer_auto ) {
334 ds=108;
335 }
336 else {
337 if(sym == 1) {
338 ds+=14;
339 }
340
341 /** Terminal AT penalty **/
342
343 if(strncmp("A", s, 1)==0
344 || strncmp("T", s, 1)==0) {
345 ds += -41;
346 dh += -23;
347 } else if (strncmp("C", s, 1)==0
348 || strncmp("G", s, 1)==0) {
349 ds += 28;
350 dh += -1;
351 }
352 s+=len;
353 if(strncmp("T", s, 1)==0
354 || strncmp("A", s, 1)==0) {
355 ds += -41;
356 dh += -23;
357 } else if (strncmp("C", s, 1)==0
358 || strncmp("G", s, 1)==0) {
359 ds += 28;
360 dh += -1;
361 }
362 s-=len;
363 }
364 /* Use a finite-state machine (DFA) to calucluate dh and ds for s. */
365 c = *s; s++;
366 if (tm_method == breslauer_auto) {
367 if (c == 'A') goto A_STATE;
368 else if (c == 'G') goto G_STATE;
369 else if (c == 'T') goto T_STATE;
370 else if (c == 'C') goto C_STATE;
371 else if (c == 'N') goto N_STATE;
372 else goto ERROR;
373 STATE(A);
374 STATE(T);
375 STATE(G);
376 STATE(C);
377 STATE(N);
378 } else {
379 if (c == 'A') goto A_STATE2;
380 else if (c == 'G') goto G_STATE2;
381 else if (c == 'T') goto T_STATE2;
382 else if (c == 'C') goto C_STATE2;
383 else if (c == 'N') goto N_STATE2;
384 else goto ERROR;
385 STATE2(A);
386 STATE2(T);
387 STATE2(G);
388 STATE2(C);
389 STATE2(N);
390 }
391
392
393 DONE: /* dh and ds are now computed for the given sequence. */
394 delta_H = dh * -100.0; /*
395 * Nearest-neighbor thermodynamic values for dh
396 * are given in 100 cal/mol of interaction.
397 */
398 delta_S = ds * -0.1; /*
399 * Nearest-neighbor thermodynamic values for ds
400 * are in in .1 cal/K per mol of interaction.
401 */
402 Tm=0; /* Melting temperature */
403 len=len+1;
404
405 /**********************************************/
406 if (salt_corrections == schildkraut) {
407 K_mM = K_mM + divalent_to_monovalent(divalent_conc, dntp_conc);
408 correction = 16.6 * log10(K_mM/1000.0) - T_KELVIN;
409 Tm = delta_H / (delta_S + 1.987 * log(DNA_nM/4000000000.0)) + correction;
410 } else if (salt_corrections== santalucia) {
411 K_mM = K_mM + divalent_to_monovalent(divalent_conc, dntp_conc);
412 delta_S = delta_S + 0.368 * (len - 1) * log(K_mM / 1000.0 );
413 if(sym == 1) { /* primer is symmetrical */
414 /* Equation A */
415 Tm = delta_H / (delta_S + 1.987 * log(DNA_nM/1000000000.0)) - T_KELVIN;
416 } else {
417 /* Equation B */
418 Tm = delta_H / (delta_S + 1.987 * log(DNA_nM/4000000000.0)) - T_KELVIN;
419 }
420 } else if (salt_corrections == owczarzy) {
421 double gcPercent=0;
422 double free_divalent; /* conc of divalent cations minus dNTP conc */
423 int i;
424 for(i = 0; i <= len && d != NULL && d != '\0';) {
425 if(*d == 'C' || *d == 'G') {
426 gcPercent++;
427 }
428 d++;
429 i++;
430 }
431 gcPercent = (double)gcPercent/((double)len);
432 /**** BEGIN: UPDATED SALT BY OWCZARZY *****/
433 /* different salt corrections for monovalent (Owczarzy et al.,2004)
434 and divalent cations (Owczarzy et al.,2008)
435 */
436 /* competition bw magnesium and monovalent cations, see Owczarzy et al., 2008 Figure 9 and Equation 16 */
437
438 static const double crossover_point = 0.22; /* depending on the value of div_monov_ratio respect
439 to value of crossover_point Eq 16 (divalent corr, Owczarzy et al., 2008)
440 or Eq 22 (monovalent corr, Owczarzy et al., 2004) should be used */
441 double div_monov_ratio;
442 if(dntp_conc >= divalent_conc) {
443 free_divalent = 0.00000000001; /* to not to get log(0) */
444 } else {
445 free_divalent = (divalent_conc - dntp_conc)/1000.0;
446 }
447 static double a = 0,b = 0,c = 0,d = 0,e = 0,f = 0,g = 0;
448 if(K_mM==0) {
449 div_monov_ratio = 6.0;
450 } else {
451 div_monov_ratio = (sqrt(free_divalent))/(K_mM/1000); /* if conc of monov cations is provided
452 a ratio is calculated to further calculate
453 the _correct_ correction */
454 }
455 if (div_monov_ratio < crossover_point) {
456 /* use only monovalent salt correction, Eq 22 (Owczarzy et al., 2004) */
457 correction
458 = (((4.29 * gcPercent) - 3.95) * pow(10,-5) * log(K_mM / 1000.0))
459 + (9.40 * pow(10,-6) * (pow(log(K_mM / 1000.0),2)));
460 } else {
461 /* magnesium effects are dominant, Eq 16 (Owczarzy et al., 2008) is used */
462 b =- 9.11 * pow(10,-6);
463 c = 6.26 * pow(10,-5);
464 e =- 4.82 * pow(10,-4);
465 f = 5.25 * pow(10,-4);
466 a = 3.92 * pow(10,-5);
467 d = 1.42 * pow(10,-5);
468 g = 8.31 * pow(10,-5);
469 if(div_monov_ratio < 6.0) {
470 /* in particular ratio of conc of monov and div cations
471 * some parameters of Eq 16 must be corrected (a,d,g) */
472 a = 3.92 * pow(10,-5) * (0.843 - (0.352 * sqrt(K_mM/1000.0) * log(K_mM/1000.0)));
473 d = 1.42 * pow(10,-5) * (1.279 - 4.03 * pow(10,-3) * log(K_mM/1000.0) - 8.03 * pow(10,-3) * pow(log(K_mM/1000.0),2));
474 g = 8.31 * pow(10,-5) * (0.486 - 0.258 * log(K_mM/1000.0) + 5.25 * pow(10,-3) * pow(log(K_mM/1000.0),3));
475 }
476
477 correction = a + (b * log(free_divalent))
478 + gcPercent * (c + (d * log(free_divalent)))
479 + (1/(2 * (len - 1))) * (e + (f * log(free_divalent))
480 + g * (pow((log(free_divalent)),2)));
481 }
482 /**** END: UPDATED SALT BY OWCZARZY *****/
483 if (sym == 1) {
484 /* primer is symmetrical */
485 /* Equation A */
486 Tm = 1/((1/(delta_H
487 /
488 (delta_S + 1.9872 * log(DNA_nM/1000000000.0)))) + correction) - T_KELVIN;
489 } else {
490 /* Equation B */
491 Tm = 1/((1/(delta_H
492 /
493 (delta_S + 1.9872 * log(DNA_nM/4000000000.0)))) + correction) - T_KELVIN;
494 }
495 } /* END else if (salt_corrections == owczarzy) { */
496
497
498 /***************************************/
499 return Tm;
500 ERROR: /*
501 * length of s was less than 2 or there was an illegal character in
502 * s.
503 */
504 return OLIGOTM_ERROR;
505 }
506 #undef DO_PAIR
507 #undef DO_PAIR2
508
509 #define DO_PAIR(LAST,THIS) \
510 if (CATID2(THIS,_CHAR) == c) { \
511 dg += CATID5(G,_,LAST,_,THIS); \
512 goto CATID2(THIS,_STATE); \
513 }
514
515 #define DO_PAIR2(LAST,THIS) \
516 if (CATID2(THIS,_CHAR) == c) { \
517 dg += CATID5(DG,_,LAST,_,THIS); \
518 goto CATID2(THIS,_STATE2); \
519 }
520
521 double
oligodg(const char * s,int tm_method)522 oligodg(const char *s, /* The sequence. */
523 int tm_method)
524 {
525 register int dg = 0;
526 register char c;
527
528 if (tm_method != breslauer_auto
529 && tm_method != santalucia_auto)
530 return OLIGOTM_ERROR;
531
532 /* Use a finite-state machine (DFA) to calucluate dg s. */
533 c = *s; s++;
534 if(tm_method != breslauer_auto) {
535 dg=-1960; /* Initial dG */
536 if(c == 'A' || c == 'T') {
537 dg += -50; /* terminal AT penalty */
538 }
539 if (c == 'A') goto A_STATE2;
540 else if (c == 'G') goto G_STATE2;
541 else if (c == 'T') goto T_STATE2;
542 else if (c == 'C') goto C_STATE2;
543 else if (c == 'N') goto N_STATE2;
544 else goto ERROR;
545 STATE2(A);
546 STATE2(T);
547 STATE2(G);
548 STATE2(C);
549 STATE2(N);
550 } else {
551 if (c == 'A') goto A_STATE;
552 else if (c == 'G') goto G_STATE;
553 else if (c == 'T') goto T_STATE;
554 else if (c == 'C') goto C_STATE;
555 else if (c == 'N') goto N_STATE;
556 else goto ERROR;
557 STATE(A);
558 STATE(T);
559 STATE(G);
560 STATE(C);
561 STATE(N);
562
563 }
564 DONE: /* dg is now computed for the given sequence. */
565 if(tm_method != breslauer_auto) {
566 int sym;
567 --s; --s; c = *s;
568 if(c == 'A' || c == 'T') {
569 dg += -50; /* terminal AT penalty */
570 }
571 sym = symmetry(s);
572 if(sym==1) {
573 dg +=-430; /* symmetry correction for dG */
574 }
575 }
576 return dg / 1000.0;
577
578 ERROR: /*
579 * length of s was less than 2 or there was an illegal character in
580 * s.
581 */
582 return OLIGOTM_ERROR;
583 }
584
end_oligodg(const char * s,int len,int tm_method)585 double end_oligodg(const char *s,
586 int len, /* The number of characters to return. */
587 int tm_method)
588 {
589 int x = strlen(s);
590
591 if (tm_method != breslauer_auto
592 && tm_method != santalucia_auto)
593 return OLIGOTM_ERROR;
594
595 return
596 x < len
597 ? oligodg(s,tm_method) :
598 oligodg(s + (x - len),tm_method);
599 }
600
601 /* See oligotm.h for documentation of arguments. */
seqtm(const char * seq,double dna_conc,double salt_conc,double divalent_conc,double dntp_conc,int nn_max_len,tm_method_type tm_method,salt_correction_type salt_corrections)602 double seqtm(const char *seq,
603 double dna_conc,
604 double salt_conc,
605 double divalent_conc,
606 double dntp_conc,
607 int nn_max_len,
608 tm_method_type tm_method,
609 salt_correction_type salt_corrections)
610 {
611 int len = strlen(seq);
612 if (tm_method != breslauer_auto
613 && tm_method != santalucia_auto)
614 return OLIGOTM_ERROR;
615 if (salt_corrections != schildkraut
616 && salt_corrections != santalucia
617 && salt_corrections != owczarzy)
618 return OLIGOTM_ERROR;
619
620 if (len > nn_max_len) {
621 return long_seq_tm(seq, 0, len, salt_conc, divalent_conc, dntp_conc);
622 } else {
623 return oligotm(seq, dna_conc, salt_conc,
624 divalent_conc, dntp_conc, tm_method, salt_corrections);
625 }
626 }
627
628 /* See oligotm.h for documentation on this function and the formula it
629 uses. */
630 double
long_seq_tm(const char * s,int start,int len,double salt_conc,double divalent_conc,double dntp_conc)631 long_seq_tm(const char *s,
632 int start,
633 int len,
634 double salt_conc,
635 double divalent_conc,
636 double dntp_conc)
637 {
638 int GC_count = 0;
639 const char *p, *end;
640
641 if (divalent_to_monovalent(divalent_conc, dntp_conc) == OLIGOTM_ERROR)
642 return OLIGOTM_ERROR;
643
644 salt_conc = salt_conc + divalent_to_monovalent(divalent_conc, dntp_conc);
645
646 if ((unsigned) (start + len) > strlen(s) || start < 0 || len <= 0)
647 return OLIGOTM_ERROR;
648 end = &s[start + len];
649 /* Length <= 0 is nonsensical. */
650 for (p = &s[start]; p < end; p++) {
651 if ('G' == *p || 'C' == *p)
652 GC_count++;
653 }
654
655 return
656 81.5
657 + (16.6 * log10(salt_conc / 1000.0))
658 + (41.0 * (((double) GC_count) / len))
659 - (600.0 / len);
660 }
661
662 /* Return 1 if string is symmetrical, 0 otherwise. */
symmetry(const char * seq)663 int symmetry(const char* seq) {
664 register char s;
665 register char e;
666 const char *seq_end=seq;
667 int i = 0;
668 int seq_len=strlen(seq);
669 int mp = seq_len/2;
670 if(seq_len%2==1) {
671 return 0;
672 }
673 seq_end+=seq_len;
674 seq_end--;
675 while(i<mp) {
676 i++;
677 s=*seq;
678 e=*seq_end;
679 if ((s=='A' && e!='T')
680 || (s=='T' && e!='A')
681 || (e=='A' && s!='T')
682 || (e=='T' && s!='A')) {
683 return 0;
684 }
685 if ((s=='C' && e!='G')
686 || (s=='G' && e!='C')
687 || (e=='C' && s!='G')
688 || (e=='G' && s!='C')) {
689 return 0;
690 }
691 seq++;
692 seq_end--;
693 }
694 return 1;
695 }
696
697 /* Convert divalent salt concentration to monovalent */
divalent_to_monovalent(double divalent,double dntp)698 double divalent_to_monovalent(double divalent,
699 double dntp)
700 {
701 if(divalent==0) dntp=0;
702 if(divalent<0 || dntp<0) return OLIGOTM_ERROR;
703 if(divalent<dntp)
704 /* According to theory, melting temperature does not depend on
705 divalent cations */
706 divalent=dntp;
707 return 120*(sqrt(divalent-dntp));
708 }
709