Added Therm. Para. SantaLucia2004 and fixed SC Owczarzy

primer3-org · Oct 10, 2024 · c5b092b · c5b092b
1 parent 23fee1d
commit c5b092b
Show file tree

Hide file tree

Showing 7 changed files with 749 additions and 43 deletions.
diff --git a/src/oligotm.c b/src/oligotm.c
@@ -39,7 +39,7 @@ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 #include <math.h>
 #include <string.h>
 #include "oligotm.h"
-#include <stdio.h>
+/* #include <stdio.h> */ /* Only for testing */
 
 #ifndef MAX_PRIMER_LENGTH
 #define MAX_PRIMER_LENGTH 36
@@ -117,6 +117,32 @@ static const int SantaLucia_1998_dG[5][5] =  {
      { 580, 1300, 1280,  880,  580}}; /* NA, NC, NG, NT, NN; */
 
 
+/* Table 3, updated parameters:
+ * Tables of nearest-neighbor thermodynamics for DNA bases, from the
+ * paper [SantaLucia JR and Hicks (2004) "The thermodynamics of DNA
+ * Structuralmotifs", PAnnu. Rev. Biophys. Biomol. Struct. 2004. 33:
+ * 415–40 http://dx.doi.org/10.1146/annurev.biophys.32.110601.141800]
+ */
+
+/* dH *-100 cal/mol */
+static const int SantaLucia_2004_dH[5][5] =  {
+     {76, 84,  78, 72, 72},  /* AA, AC, AG, AT, AN; */
+     {85, 80, 106, 78, 78},  /* CA, CC, CG, CT, CN; */
+     {82, 98,  80, 84, 80},  /* GA, GC, GG, GT, GN; */
+     {72, 82,  85, 76, 72},  /* TA, TC, TG, TT, TN; */
+     {72, 80,  78, 72, 72}}; /* NA, NC, NG, NT, NN; */
+
+/* dS *-0.1 cal/k*mol */
+static const int SantaLucia_2004_dS[5][5] =  {
+     {213, 224, 210, 204, 224},  /* AA, AC, AG, AT, AN; */
+     {227, 199, 272, 210, 272},  /* CA, CC, CG, CT, CN; */
+     {222, 244, 199, 224, 244},  /* GA, GC, GG, GT, GN; */
+     {213, 222, 227, 213, 227},  /* TA, TC, TG, TT, TN; */
+     {168, 210, 220, 215, 220}}; /* NA, NC, NG, NT, NN; */
+
+/* dG *-0.001 cal/mol */
+/* unmodified to SantaLucia_1998_dG */
+
 /* Calculate the melting temperature of oligo s.  See
    oligotm.h for documentation of arguments.
 */
@@ -149,8 +175,8 @@ oligotm(const  char *s,
     return ret;
   }
   /** K_mM = K_mM + divalent_to_monovalent(divalent_conc, dntp_conc); **/
-  if (tm_method != breslauer_auto
-      && tm_method != santalucia_auto) {
+  if (tm_method != breslauer_auto && tm_method != santalucia_auto
+      && tm_method != santalucia_2004) {
     ret.Tm = OLIGOTM_ERROR;
     return ret;
   }
@@ -225,6 +251,35 @@ oligotm(const  char *s,
         dh += -1;
         break;
     }
+  } else if(tm_method == santalucia_2004) {
+    ds += 57;
+    dh += -2;
+
+    if(sym == 1) {
+      ds += 14;
+    }
+
+    /** Terminal penalty **/
+    switch (s[0]) {
+      case 'A':
+      case 'T':
+        ds += -69;
+        dh += -22;
+        break;
+    }
+    /* Sum the pairs up */
+    for (i = 0; i < len; i++) {
+      ds += SantaLucia_2004_dS[oligo_int[i]][oligo_int[i+1]];
+      dh += SantaLucia_2004_dH[oligo_int[i]][oligo_int[i+1]];
+    }
+    /** Terminal penalty **/
+    switch (s[len]) {
+      case 'A':
+      case 'T':
+        ds += -69;
+        dh += -22;
+        break;
+    }
   } else {
     ds = 108;
     /* Sum the pairs up */
@@ -245,6 +300,9 @@ oligotm(const  char *s,
   ret.Tm=0;  /* Melting temperature */
   len=len+1;
 
+  /* For testing */
+  /* printf("dH: %.2f kcal/mol\ndS: %.2f cal/k*mol\n", delta_H / 1000.0, delta_S); */
+
    /**********************************************/
   if (salt_corrections == schildkraut) {
     K_mM = K_mM + divalent_to_monovalent(divalent_conc, dntp_conc);
@@ -282,9 +340,9 @@ oligotm(const  char *s,
       }
     }
   } else if (salt_corrections == owczarzy) {
-    double gcPercent;
+    double gcFract;
     double free_divalent; /* conc of divalent cations minus dNTP conc */
-    gcPercent = (double)GC_count/((double)len);
+    gcFract = (double)GC_count/((double)len);
     /**** BEGIN: UPDATED SALT BY OWCZARZY *****/
     /* different salt corrections for monovalent (Owczarzy et al.,2004)
        and divalent cations (Owczarzy et al.,2008)
@@ -310,17 +368,17 @@ oligotm(const  char *s,
     }
     if (div_monov_ratio < crossover_point) {
       /* use only monovalent salt correction, Eq 22 (Owczarzy et al., 2004) */
-      correction = (((4.29 * gcPercent) - 3.95) * pow(10,-5) * log(K_mM / 1000.0))
+      correction = (((4.29 * gcFract) - 3.95) * pow(10,-5) * log(K_mM / 1000.0))
                    + (9.40 * pow(10,-6) * (pow(log(K_mM / 1000.0),2)));
     } else {
       /* magnesium effects are dominant, Eq 16 (Owczarzy et al., 2008) is used */
-      b =- 9.11 * pow(10,-6);
-      c = 6.26 * pow(10,-5);
-      e =- 4.82 * pow(10,-4);
-      f = 5.25 * pow(10,-4);
-      a = 3.92 * pow(10,-5);
-      d = 1.42 * pow(10,-5);
-      g = 8.31 * pow(10,-5);
+      a =  3.92 * pow(10,-5);
+      b = -9.11 * pow(10,-6);
+      c =  6.26 * pow(10,-5);
+      d =  1.42 * pow(10,-5);
+      e = -4.82 * pow(10,-4);
+      f =  5.25 * pow(10,-4);
+      g =  8.31 * pow(10,-5);
       if(div_monov_ratio < 6.0) {
         /* in particular ratio of conc of monov and div cations
          *             some parameters of Eq 16 must be corrected (a,d,g) */
@@ -330,9 +388,8 @@ oligotm(const  char *s,
       }
 
       correction = a + (b * log(free_divalent))
-                     + gcPercent * (c + (d * log(free_divalent)))
-                     + (1/(2 * (len - 1))) * (e + (f * log(free_divalent))
-                     + g * (pow((log(free_divalent)),2)));
+                     + gcFract * (c + (d * log(free_divalent)))
+                     + (1/(2 * (((double) len)- 1))) * (e + (f * log(free_divalent)) + g * (pow((log(free_divalent)),2)));
     }
     /**** END: UPDATED SALT BY OWCZARZY *****/
     if (sym == 1) {
@@ -348,7 +405,7 @@ oligotm(const  char *s,
                       (delta_S + 1.9872 * log(DNA_nM/4000000000.0)))) + correction) - t_kelvin;
     }
     ret.Tm -= dmso_conc * dmso_fact;
-    ret.Tm += (0.453 * ((double) GC_count) / len - 2.88) * formamide_conc;
+    ret.Tm += (0.453 * ((double) GC_count) / ((double) len) - 2.88) * formamide_conc;
   } /* END else if (salt_corrections == owczarzy) { */
   /***************************************/
   return ret;
@@ -390,7 +447,7 @@ oligodg(const char *s,      /* The sequence. */
     }
   }
 
-  if( tm_method == santalucia_auto ) {
+  if ((tm_method == santalucia_auto) || (tm_method == santalucia_2004)) {
     dg  = -1960; /* Initial dG */
 
     sym = symmetry(s);
@@ -461,7 +518,8 @@ tm_ret seqtm(const  char *seq,
   ret.bound = OLIGOTM_ERROR;
   ret.Tm = OLIGOTM_ERROR;
   if (tm_method != breslauer_auto
-      && tm_method != santalucia_auto)
+      && tm_method != santalucia_auto
+      && tm_method != santalucia_2004)
     return ret;
   if (salt_corrections != schildkraut
       && salt_corrections != santalucia
@@ -511,10 +569,10 @@ long_seq_tm(const char *s,
   }
 
   ret.Tm = 81.5 - dmso_conc * dmso_fact
-                + (0.453 * ((double) GC_count) / len - 2.88) * formamide_conc
+                + (0.453 * ((double) GC_count) / ((double) len) - 2.88) * formamide_conc
                 + (16.6 * log10(salt_conc / 1000.0))
-                + (41.0 * (((double) GC_count) / len))
-                - (600.0 / len);
+                + (41.0 * (((double) GC_count) / ((double) len)))
+                - (600.0 / ((double) len));
 
   return ret;
 }

diff --git a/src/oligotm.h b/src/oligotm.h
@@ -111,6 +111,7 @@ tm_ret long_seq_tm(const char *seq,
 typedef enum tm_method_type {
         breslauer_auto      = 0,
         santalucia_auto     = 1,
+        santalucia_2004     = 2,
 } tm_method_type;
 
 typedef enum salt_correction_type {

diff --git a/src/oligotm_main.c b/src/oligotm_main.c
@@ -79,6 +79,7 @@ main(int argc, char **argv)
     "                    (used by primer3 up to and including release 1.1.0).\n"
     "                 1  Use nearest neighbor parameters from SantaLucia 1998\n"
     "                    *This is the default and recommended value*\n"
+    "                 2  Use nearest neighbor parameters from SantaLucia 2004\n"
     "\n"
     "-sc [0..2]    - Specifies salt correction formula for the melting \n"
     "                 temperature calculation\n"
@@ -89,8 +90,8 @@ main(int argc, char **argv)
     "                  2  Owczarzy et al., 2004\n\n"
     "\n\n"
     "Prints oligo's melting temperature on stdout.\n";
-   
-   const char *copyright = 
+
+   const char *copyright =
 "Copyright (c) 1996,1997,1998,1999,2000,2001,2004,2006\n"
 "Whitehead Institute for Biomedical Research, Steve Rozen\n"
 "(http://purl.com/STEVEROZEN/), Andreas Untergasser and Helen Skaletsky\n"
@@ -122,7 +123,7 @@ main(int argc, char **argv)
    int tm_santalucia=1, salt_corrections=1;
    int i, j, len;
    if (argc < 2 || argc > 14) {
-     fprintf(stderr, msg, argv[0]);       
+     fprintf(stderr, msg, argv[0]);
      fprintf(stderr, "%s", copyright);
      return -1;
    }
@@ -219,7 +220,7 @@ main(int argc, char **argv)
          exit(-1);
        }
        tm_santalucia = (int)strtol(argv[i+1], &endptr, 10);
-       if ('\0' != *endptr || tm_santalucia<0 || tm_santalucia>1) {          
+       if ('\0' != *endptr || tm_santalucia<0 || tm_santalucia>2) {
          fprintf(stderr, msg, argv[0]);
          exit(-1);
        }
@@ -243,7 +244,7 @@ main(int argc, char **argv)
      } else
        break;                /* all args processed. go on to sequences. */
    }
-   
+
   if(!argv[i]) { /* if no oligonucleotide sequence is specified */
     fprintf(stderr, msg, argv[0]);
     exit(-1);
@@ -252,15 +253,15 @@ main(int argc, char **argv)
   seq = argv[i];
   len=strlen(seq);
   for(j=0;j<len;j++) seq[j]=toupper(seq[j]);
-   
+
   tm_calc = oligotm(seq, d, mv, dv, n, dmso, dmso_fact, formamide,
                     (tm_method_type) tm_santalucia, (salt_correction_type) salt_corrections, -10.0);
   tm = tm_calc.Tm;
 
   if (OLIGOTM_ERROR == tm) {
     fprintf(stderr,
             "%s ERROR: length of sequence %s is less than 2 or\n"
-            "             the sequence contains an illegal character or\n" 
+            "             the sequence contains an illegal character or\n"
             "             you have specified incorrect value for concentration of divalent cations or\n"
             "             you have specified incorrect value for concentration of dNTPs\n",
             argv[0], argv[i]);

diff --git a/src/release_notes.txt b/src/release_notes.txt
@@ -19,6 +19,12 @@ Probes can be quenched if a G is on the 5' end. Setting
 PRIMER_INTERNAL_NO_5_PRIME_G=1 will force Primer3 to pick only internal 
 oligos with no 5' G.
 
+A bug was fixed in Owczarzy salt correction which changes the Tm 
+results.
+
+The nearest-neighbor thermodynamics from SantaLucia JR and Hicks (2004)
+were added.
+
 Minor fixes:
 1. Compiler complains were fixed.
 
@@ -28,6 +34,8 @@ Minor fixes:
 
 4. oligotm was rewritten to avoid makros.
 
+5. A bug was fixed in Owczarzy salt correction.
+
 
 release 2.6.1 2022-01-26  ===============================================
 This release fixes a small memory bug.

diff --git a/test/oligotm.txt b/test/oligotm.txt
@@ -1,7 +1,7 @@
 SEQ	TM_MODEL	SALT_CORRECTION	SALT_CONC	DIVALENT_CONC	DNTP_CONC	TM
-AAGGCGAGTCAGGCTCAGTG	1	2	5	1.5	0	62.123401
-AAGGCGAGTCAGGCTCAGTG	1	2	0	1.5	0	62.123401
-AAGGCGAGTCAGGCTCAGTG	1	2	10	1.5	0	62.895198
+AAGGCGAGTCAGGCTCAGTG	1	2	5	1.5	0	63.257974
+AAGGCGAGTCAGGCTCAGTG	1	2	0	1.5	0	63.257974
+AAGGCGAGTCAGGCTCAGTG	1	2	10	1.5	0	62.346915
 AAAGAAAGACCCAGTGATATGCTGTCT	0	0	50	0	0	63.684650
 ACAACCCAGGCCCTGTGGAC	0	0	50	0	0.8	66.783376
 ACTGGGATGACAGCAAGGCCC	0	0	50	0	1	67.947254
@@ -20,9 +20,9 @@ CCCGCGTGTGGTTTCTCTATTTG	0	1	50	0	0	71.826149
 CCTCCCGGGTTCACACCATT	0	1	50	1.5	0	79.028727
 CCTTCTGACAAGGGATTAATAACCAGA	0	1	50	2.5	0.8	75.999617
 CGCAAGCCACCACTGTCCAT	0	2	50	0	2	72.669877
-CTCGGTTATCCCTGGAATCTGGA	0	2	50	1	0	77.619980
-CTGCCTGCACCCGCTTGTTT	0	2	50	1.5	0.9	80.472263
-CTTCATGATCCACTGGCTATTACGG	0	2	50	2.5	1.0	76.728968
+CTCGGTTATCCCTGGAATCTGGA	0	2	50	1	0	76.717289
+CTGCCTGCACCCGCTTGTTT	0	2	50	1.5	0.9	78.069281
+CTTCATGATCCACTGGCTATTACGG	0	2	50	2.5	1.0	76.648737
 GAAACCTCACCCGGAGCCAGTA	1	0	50	0	2	51.910898
 GACAACACCTGGACACCGGTTA	1	0	50	0	10	50.261568
 GAGGCAGAGAGAGAGGAAGAGACAGA	1	0	50	0	0	52.569759
@@ -41,11 +41,37 @@ GTGAAACTCACGGTCACGCAAA	1	1	50	2.5	0	64.004323
 TCACGCTCACAAACGAGAGGAAG	1	2	50	0	1	56.755143
 TCCCTTGCCTCGTTCAGTGC	1	2	50	0	0	57.675190
 TCGGCCTCCCAAAGTGCTG	1	2	50	0	1.5	57.961561
-TCTCGCCATTTTGGTCTTGCAG	1	2	50	1	0	62.162400
-TGATCAGCACAGTGGCTCACG	1	2	50	1.5	0	63.699372
-TGCCGTTCACGTCGTAGCTG	1	2	50	1.5	0.8	63.844611
-TGGGCCTCCCAAAGTGCTG	1	2	50	2.5	0.1	63.716743
-TTGCATGCCCAGAGTTGACG	1	2	50	2.5	1	61.900322
-AAGGCGAGTCAGGCTCAGTG	1	2	5.5	1.5	0	62.123401
-AAGGCGAGTCAGGCTCAGTG	1	2	0.5	1.5	0	62.123401
-AAGGCGAGTCAGGCTCAGTG	1	2	10.5	1.5	0	62.906493
+TCTCGCCATTTTGGTCTTGCAG	1	2	50	1	0	61.298237
+TGATCAGCACAGTGGCTCACG	1	2	50	1.5	0	63.610135
+TGCCGTTCACGTCGTAGCTG	1	2	50	1.5	0.8	62.043487
+TGGGCCTCCCAAAGTGCTG	1	2	50	2.5	0.1	64.566813
+TTGCATGCCCAGAGTTGACG	1	2	50	2.5	1	61.807390
+AAGGCGAGTCAGGCTCAGTG	1	2	5.5	1.5	0	63.257974
+AAGGCGAGTCAGGCTCAGTG	1	2	0.5	1.5	0	63.257974
+AAGGCGAGTCAGGCTCAGTG	2	2	10.5	1.5	0	62.270853
+GAAACCTCACCCGGAGCCAGTA	2	0	50	0	2	51.889195
+GACAACACCTGGACACCGGTTA	2	0	50	0	10	50.234394
+GAGGCAGAGAGAGAGGAAGAGACAGA	2	0	50	0	0	52.531380
+GATCCTCGGGCAGGCAGTTG	2	0	50	1	0	59.632470
+GCAAGACAAAGGCATGCTTACAAGA	2	0	50	1.5	1	57.546566
+GCACAGTGGCCCTTGTGACTG	2	0	50	2.5	1.0	61.551715
+GCAGCCTGGACATAGCCCTGTA	2	1	50	0	0	58.044193
+GCAGGCCAGAAGGCAGGATT	2	1	50	0	5	56.833896
+GCCATGGACGTGGTCTACGC	2	1	50	0	0	57.733702
+GCGACCCAAGCGAGTTACCA	2	1	50	1	0	62.897052
+GCTTCATGTCAGATTTGTTATACCACTGGG	2	1	50	1	1	58.033837
+GGCAGGCAGTTGGTCGGTGT	2	1	50	1.5	1.5	59.481855
+GGGCCATGGGATTAGTCTGC	2	1	50	2.5	0	62.395268
+GGTATGGTGTTTTAATACCTGTGGCAGA	1	1	50	2.5	1	64.204218
+GTGAAACTCACGGTCACGCAAA	2	1	50	2.5	0	63.965960
+TCACGCTCACAAACGAGAGGAAG	2	2	50	0	1	56.751572
+TCCCTTGCCTCGTTCAGTGC	2	2	50	0	0	57.648545
+TCGGCCTCCCAAAGTGCTG	2	2	50	0	1.5	57.932687
+TCTCGCCATTTTGGTCTTGCAG	2	2	50	1	0	61.306914
+TGATCAGCACAGTGGCTCACG	2	2	50	1.5	0	63.542322
+TGCCGTTCACGTCGTAGCTG	2	2	50	1.5	0.8	61.995870
+TGGGCCTCCCAAAGTGCTG	2	2	50	2.5	0.1	64.536728
+TTGCATGCCCAGAGTTGACG	2	2	50	2.5	1	61.774378
+AAGGCGAGTCAGGCTCAGTG	2	2	5.5	1.5	0	63.207179
+AAGGCGAGTCAGGCTCAGTG	2	2	0.5	1.5	0	63.207179
+AAGGCGAGTCAGGCTCAGTG	2	2	10.5	1.5	0	62.270853