← Back to branch summary

~ubuntu-branches/ubuntu/intrepid/primer3/intrepid

« back to all changes in this revision

Viewing changes to src/oligotm.c

  • Committer: Bazaar Package Importer
  • Author(s): Charles Plessy
  • Date: 2006-09-28 20:18:54 UTC
  • Revision ID: james.westby@ubuntu.com-20060928201854-45pwapz5e3a6d684
Tags: upstream-1.0b
Import upstream version 1.0b

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/oligotm.c
 
1
/*
 
2
Copyright (c) 1996,1997,1998,1999,2000,2001,2004,2006
 
3
Whitehead Institute for Biomedical Research, Steve Rozen
 
4
(http://jura.wi.mit.edu/rozen), and Helen Skaletsky
 
5
All rights reserved.
 
6
 
 
7
Redistribution and use in source and binary forms, with or without
 
8
modification, are permitted provided that the following conditions are
 
9
met:
 
10
 
 
11
   * Redistributions of source code must retain the above copyright
 
12
notice, this list of conditions and the following disclaimer.
 
13
   * Redistributions in binary form must reproduce the above
 
14
copyright notice, this list of conditions and the following disclaimer
 
15
in the documentation and/or other materials provided with the
 
16
distribution.
 
17
   * Neither the names of the copyright holders nor contributors may
 
18
be used to endorse or promote products derived from this software
 
19
without specific prior written permission.
 
20
 
 
21
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 
22
"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 
23
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
 
24
A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
 
25
OWNERS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
 
26
SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
 
27
LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
 
28
DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
 
29
THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 
30
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 
31
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
32
*/
 
33
 
 
34
#include <limits.h>
 
35
#include <math.h>
 
36
#include <string.h>
 
37
#include "oligotm.h"
 
38
#include "primer3_release.h"
 
39
 
 
40
/* 
 
41
 * Tables of nearest-neighbor thermodynamics for DNA bases.
 
42
 * See Breslauer, Frank, Blocker, and Markey, 
 
43
 * "Predicting DNA duplex stability from the base sequence."
 
44
 * Proc. Natl. Acad. Sci. USA, vol 83, page 3746 (1986).
 
45
 * Article free at
 
46
 * http://www.pubmedcentral.nih.gov/picrender.fcgi?artid=323600&blobtype=pdf
 
47
 * See table 2.
 
48
 */
 
49
#define S_A_A 240
 
50
#define S_A_C 173
 
51
#define S_A_G 208
 
52
#define S_A_T 239
 
53
#define S_A_N 215
 
54
  
 
55
#define S_C_A 129
 
56
#define S_C_C 266
 
57
#define S_C_G 278
 
58
#define S_C_T 208
 
59
#define S_C_N 220  
 
60
  
 
61
#define S_G_A 135
 
62
#define S_G_C 267
 
63
#define S_G_G 266
 
64
#define S_G_T 173
 
65
#define S_G_N 210
 
66
  
 
67
#define S_T_A 169
 
68
#define S_T_C 135
 
69
#define S_T_G 129
 
70
#define S_T_T 240
 
71
#define S_T_N 168
 
72
  
 
73
#define S_N_A 168
 
74
#define S_N_C 210
 
75
#define S_N_G 220
 
76
#define S_N_T 215
 
77
#define S_N_N 203
 
78
 
 
79
 
 
80
#define H_A_A  91
 
81
#define H_A_C  65
 
82
#define H_A_G  78
 
83
#define H_A_T  86
 
84
#define H_A_N  80
 
85
 
 
86
#define H_C_A  58
 
87
#define H_C_C 110
 
88
#define H_C_G 119
 
89
#define H_C_T  78
 
90
#define H_C_N  91
 
91
 
 
92
#define H_G_A  56
 
93
#define H_G_C 111
 
94
#define H_G_G 110
 
95
#define H_G_T  65
 
96
#define H_G_N  85
 
97
 
 
98
#define H_T_A  60
 
99
#define H_T_C  56
 
100
#define H_T_G  58
 
101
#define H_T_T  91
 
102
#define H_T_N  66
 
103
 
 
104
#define H_N_A  66
 
105
#define H_N_C  85
 
106
#define H_N_G  91
 
107
#define H_N_T  80
 
108
#define H_N_N  80
 
109
 
 
110
/* Delta G's of disruption * 1000. */
 
111
#define G_A_A  1900
 
112
#define G_A_C  1300
 
113
#define G_A_G  1600
 
114
#define G_A_T  1500
 
115
#define G_A_N  1575
 
116
 
 
117
#define G_C_A  1900 
 
118
#define G_C_C  3100
 
119
#define G_C_G  3600
 
120
#define G_C_T  1600
 
121
#define G_C_N  2550
 
122
 
 
123
#define G_G_A  1600
 
124
#define G_G_C  3100
 
125
#define G_G_G  3100
 
126
#define G_G_T  1300
 
127
#define G_G_N  2275
 
128
 
 
129
#define G_T_A   900
 
130
#define G_T_C  1600
 
131
#define G_T_G  1900
 
132
#define G_T_T  1900
 
133
#define G_T_N  1575
 
134
 
 
135
#define G_N_A  1575
 
136
#define G_N_C  2275
 
137
#define G_N_G  2550
 
138
#define G_N_T  1575
 
139
#define G_N_N  1994
 
140
 
 
141
#define A_CHAR 'A'
 
142
#define G_CHAR 'G'
 
143
#define T_CHAR 'T'
 
144
#define C_CHAR 'C'
 
145
#define N_CHAR 'N'
 
146
 
 
147
#define CATID5(A,B,C,D,E) A##B##C##D##E
 
148
#define CATID2(A,B) A##B
 
149
#define DO_PAIR(LAST,THIS)          \
 
150
  if (CATID2(THIS,_CHAR) == c) {    \
 
151
     dh += CATID5(H,_,LAST,_,THIS); \
 
152
     ds += CATID5(S,_,LAST,_,THIS); \
 
153
     goto CATID2(THIS,_STATE);      \
 
154
  }
 
155
 
 
156
#define STATE(LAST)     \
 
157
   CATID2(LAST,_STATE): \
 
158
   c = *s; s++;         \
 
159
   DO_PAIR(LAST,A)      \
 
160
   else DO_PAIR(LAST,T) \
 
161
   else DO_PAIR(LAST,G) \
 
162
   else DO_PAIR(LAST,C) \
 
163
   else DO_PAIR(LAST,N) \
 
164
   else if ('\0' == c)  \
 
165
             goto DONE; \
 
166
   else goto ERROR \
 
167
 
 
168
double 
 
169
oligotm(s, DNA_nM, K_mM)
 
170
     const  char *s;
 
171
     double DNA_nM;
 
172
     double K_mM;
 
173
{
 
174
    register int dh = 0, ds = 108;
 
175
    register char c;
 
176
    double delta_H, delta_S;
 
177
 
 
178
    /* Use a finite-state machine (DFA) to calucluate dh and ds for s. */
 
179
    c = *s; s++;
 
180
    if (c == 'A') goto A_STATE;
 
181
    else if (c == 'G') goto G_STATE;
 
182
    else if (c == 'T') goto T_STATE;
 
183
    else if (c == 'C') goto C_STATE;
 
184
    else if (c == 'N') goto N_STATE;
 
185
    else goto ERROR;
 
186
    STATE(A);
 
187
    STATE(T);
 
188
    STATE(G);
 
189
    STATE(C);
 
190
    STATE(N);
 
191
 
 
192
 DONE:  /* dh and ds are now computed for the given sequence. */
 
193
    delta_H = dh * -100.0;  /* 
 
194
                             * Nearest-neighbor thermodynamic values for dh
 
195
                             * are given in 100 cal/mol of interaction.
 
196
                             */
 
197
    delta_S = ds * -0.1;     /*
 
198
                              * Nearest-neighbor thermodynamic values for ds
 
199
                              * are in in .1 cal/K per mol of interaction.
 
200
                              */
 
201
 
 
202
    /* 
 
203
     * See Rychlik, Spencer, Rhoads,
 
204
     * "Optimization of the annealing temperature for
 
205
     * DNA amplification in vitro."
 
206
     * Nucleic Acids Research, vol 18, no 21, page 6409 (1990).
 
207
     * Article free at 
 
208
     * http://www.pubmedcentral.nih.gov/articlerender.fcgi?tool=pubmed&pubmedid=2243783
 
209
     * See eqn (ii).
 
210
     */
 
211
    return delta_H / (delta_S + 1.987 * log(DNA_nM/4000000000.0))
 
212
        - 273.15 + 16.6 * log10(K_mM/1000.0);
 
213
 
 
214
 ERROR:  /* 
 
215
          * length of s was less than 2 or there was an illegal character in
 
216
          * s.
 
217
          */
 
218
    return OLIGOTM_ERROR;
 
219
}
 
220
#undef DO_PAIR
 
221
 
 
222
#define DO_PAIR(LAST,THIS)          \
 
223
  if (CATID2(THIS,_CHAR) == c) {    \
 
224
     dg += CATID5(G,_,LAST,_,THIS); \
 
225
     goto CATID2(THIS,_STATE);      \
 
226
  }
 
227
 
 
228
double 
 
229
oligodg(s)
 
230
    const char *s;       /* The sequence. */
 
231
{
 
232
    register int dg = 0;
 
233
    register char c;
 
234
 
 
235
    /* Use a finite-state machine (DFA) to calucluate dg s. */
 
236
    c = *s; s++;
 
237
    if (c == 'A') goto A_STATE;
 
238
    else if (c == 'G') goto G_STATE;
 
239
    else if (c == 'T') goto T_STATE;
 
240
    else if (c == 'C') goto C_STATE;
 
241
    else if (c == 'N') goto N_STATE;
 
242
    else goto ERROR;
 
243
    STATE(A);
 
244
    STATE(T);
 
245
    STATE(G);
 
246
    STATE(C);
 
247
    STATE(N);
 
248
 
 
249
 DONE:  /* dg is now computed for the given sequence. */
 
250
    return dg / 1000.0;
 
251
 
 
252
 ERROR:  /* 
 
253
          * length of s was less than 2 or there was an illegal character in
 
254
          * s.
 
255
          */
 
256
    return OLIGOTM_ERROR;
 
257
}
 
258
 
 
259
double end_oligodg(s, len)
 
260
  const char *s;
 
261
  int len; /* The number of characters to return. */
 
262
{
 
263
  int x = strlen(s);
 
264
  return x < len ? oligodg(s) : oligodg(s + (x - len));
 
265
}
 
266
 
 
267
double seqtm(seq, dna_conc, salt_conc, nn_max_len)
 
268
  const  char *seq;
 
269
  double dna_conc;
 
270
  double salt_conc;
 
271
  int    nn_max_len;
 
272
{
 
273
  int len = strlen(seq);
 
274
  return (len > nn_max_len)
 
275
    ? long_seq_tm(seq, 0, len, salt_conc) : oligotm(seq, dna_conc, salt_conc);
 
276
}
 
277
 
 
278
/* See oligotm.h for documentation on this function and the formula it
 
279
   uses. */
 
280
double
 
281
long_seq_tm(s, start, len, salt_conc)
 
282
  const char *s;
 
283
  int start, len;
 
284
  double salt_conc;
 
285
{
 
286
  int GC_count = 0;
 
287
  const char *p, *end;
 
288
 
 
289
  if(start + len > strlen(s) || start < 0 || len <= 0) return OLIGOTM_ERROR;
 
290
  end = &s[start + len];
 
291
  /* Length <= 0 is nonsensical. */
 
292
  for (p = &s[start]; p < end; p++) {
 
293
    if ('G' == *p || 'g' == *p || 'C' == *p || 'c' == *p)
 
294
      GC_count++;
 
295
  }
 
296
 
 
297
  return
 
298
    81.5
 
299
    + (16.6 * log10(salt_conc / 1000.0))
 
300
    + (41.0 * (((double) GC_count) / len))
 
301
    - (600.0 / len);
 
302
 
 
303
}