primer3-org / primer3

Primer3 is a command line tool to select primers for polymerase chain reaction (PCR).
GNU General Public License v2.0
218 stars 62 forks source link

The salt ion correction formula in file thal.c is inconsistent with that in file oligotm.c #36

Closed Wy2160640 closed 2 years ago

Wy2160640 commented 4 years ago

@brantfaircloth @untergasser @zamaudio @ioanacut @bioinfo-ut Thank you very much for your outstanding contribution to this project. About salt concentration correction formula, I checked file thal.c and file oligotm.c.
There is only one function in file thal.c that implements method santalucia, and three implemented method in file oligotm.c, schildkraut, santalucia, owczarzy. Can you add the other two methods to file thal.c.

file thal.c

static double
saltCorrectS (double mv, double dv, double dntp)
{
  if(dv <= 0) dntp = dv;
  return 0.368*((log((mv + 120*(sqrt(fmax(0.0, dv - dntp)))) / 1000)));
}

file oligotm.c

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)
{
  register int dh = 0, ds = 0;
  register char c;
  double delta_H, delta_S;
  double Tm; /* Melting temperature */
  double correction;
  int len, sym;
  const char* seq = s;
  if(divalent_to_monovalent(divalent_conc, dntp_conc) == OLIGOTM_ERROR) {
    return OLIGOTM_ERROR;
  }

  /** K_mM = K_mM + divalent_to_monovalent(divalent_conc, dntp_conc); **/
  if (tm_method != breslauer_auto
      && tm_method != santalucia_auto) {
    return OLIGOTM_ERROR;
  }
  if (salt_corrections != schildkraut
      && salt_corrections != santalucia
      && salt_corrections != owczarzy) {
    return OLIGOTM_ERROR;
  }

  len = ((int)strlen(s)-1);

  sym = symmetry(s); /*Add symmetry correction if seq is symmetrical*/
  if( tm_method == breslauer_auto ) {
    ds=108;
  }
  else {
    if(sym == 1) {
      ds+=14;
    }

    /** Terminal AT penalty **/

    if (    (strncmp("A", s, 1) == 0)
        ||  (strncmp("T", s, 1) == 0) )  {
      ds += -41;
      dh += -23;
    } else if (   (strncmp("C", s, 1) == 0)
              ||  (strncmp("G", s, 1) == 0) ) {
      ds += 28;
      dh += -1;
    }
    s += len;
    if (  (strncmp("T", s, 1) == 0)
       || (strncmp("A", s, 1) == 0) ) {
      ds += -41;
      dh += -23;
    } else if ( (strncmp("C", s, 1) == 0)
              ||(strncmp("G", s, 1) == 0) ) {
      ds += 28;
      dh += -1;
    }
    s -= len;
  }
  /* Use a finite-state machine (DFA) to calucluate dh and ds for s. */
  c = *s; s++;
  if (tm_method == breslauer_auto) {
    if (c == 'A') goto A_STATE;
    else if (c == 'G') goto G_STATE;
    else if (c == 'T') goto T_STATE;
    else if (c == 'C') goto C_STATE;
    else if (c == 'N') goto N_STATE;
    else goto ERROR;
    STATE(A);
    STATE(T);
    STATE(G);
    STATE(C);
    STATE(N);
  } else {
    if (c == 'A') goto A_STATE2;
    else if (c == 'G') goto G_STATE2;
    else if (c == 'T') goto T_STATE2;
    else if (c == 'C') goto C_STATE2;
    else if (c == 'N') goto N_STATE2;
    else goto ERROR;
    STATE2(A);
    STATE2(T);
    STATE2(G);
    STATE2(C);
    STATE2(N);
  }

 DONE:  /* dh and ds are now computed for the given sequence. */
  delta_H = dh * -100.0;  /*
         * Nearest-neighbor thermodynamic values for dh
         * are given in 100 cal/mol of interaction.
         */
  delta_S = ds * -0.1;     /*
          * Nearest-neighbor thermodynamic values for ds
          * are in in .1 cal/K per mol of interaction.
          */
  Tm=0;  /* Melting temperature */
  len=len+1;

   /**********************************************/
  if (salt_corrections == schildkraut) {
     K_mM = K_mM + divalent_to_monovalent(divalent_conc, dntp_conc);
     correction = 16.6 * log10(K_mM/1000.0) - T_KELVIN;
     Tm = delta_H / (delta_S + 1.987 * log(DNA_nM/4000000000.0)) + correction;
  } else if (salt_corrections== santalucia) {
    K_mM = K_mM + divalent_to_monovalent(divalent_conc, dntp_conc);
    delta_S = delta_S + 0.368 * (len - 1) * log(K_mM / 1000.0 );
    if(sym == 1) { /* primer is symmetrical */
      /* Equation A */
      Tm = delta_H / (delta_S + 1.987 * log(DNA_nM/1000000000.0)) - T_KELVIN;
    } else {
      /* Equation B */
      Tm = delta_H / (delta_S + 1.987 * log(DNA_nM/4000000000.0)) - T_KELVIN;
    }
  } else if (salt_corrections == owczarzy) {
    double gcPercent=0;
    double free_divalent; /* conc of divalent cations minus dNTP conc */
    int i;

    /**** BEGIN: UPDATED SALT BY OWCZARZY *****/
    /* different salt corrections for monovalent (Owczarzy et al.,2004)
      and divalent cations (Owczarzy et al.,2008) */
    /* competition bw magnesium and monovalent cations, see Owczarzy et al., 2008 Figure 9 and Equation 16 */

    static const double crossover_point = 0.22; /* depending on the value of div_monov_ratio respect
               to value of crossover_point Eq 16 (divalent corr, Owczarzy et al., 2008)
               or Eq 22 (monovalent corr, Owczarzy et al., 2004) should be used */
    double div_monov_ratio;
    static double a = 0, b = 0, c = 0, d = 0, e = 0, f = 0, g = 0;

    for(i = 0; i <= len && seq != NULL && seq != '\0';) {
      if(*seq == 'C' || *seq == 'G') {
        gcPercent++;
      }
      seq++;
      i++;
    }
    gcPercent = (double)gcPercent/((double)len);

    if(dntp_conc >= divalent_conc) {
      free_divalent = 0.00000000001; /* to not to get log(0) */
    } else {
      free_divalent = (divalent_conc - dntp_conc)/1000.0;
    }

    if(K_mM==0) {
      div_monov_ratio = 6.0;
    } else {
      div_monov_ratio = (sqrt(free_divalent))/(K_mM/1000); /* if conc of monov cations is provided
                    a ratio is calculated to further calculate
                    the _correct_ correction */
    }
    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))
          + (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);
      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) */
        a = 3.92 * pow(10,-5) * (0.843 - (0.352 * sqrt(K_mM/1000.0) * log(K_mM/1000.0)));
        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));
        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));
      }

      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)));
    }
    /**** END: UPDATED SALT BY OWCZARZY *****/
    if (sym == 1) {
      /* primer is symmetrical */
      /* Equation A */
      Tm = 1/((1/(delta_H
                    /
              (delta_S + 1.9872 * log(DNA_nM/1000000000.0)))) + correction) - T_KELVIN;
    } else {
      /* Equation B */
      Tm = 1/((1/(delta_H
          /
        (delta_S + 1.9872 * log(DNA_nM/4000000000.0)))) + correction) - T_KELVIN;
    }
  } /* END else if (salt_corrections == owczarzy) { */

   /***************************************/
   return Tm;
ERROR:  /*
    * length of s was less than 2 or there was an illegal character in
    * s.
    */
  return OLIGOTM_ERROR;
}
untergasser commented 2 years ago

Yes, that is true, but it is intentionally. A melting temperature is calculated for primers and subsequent corrected for salt concentration. This is the case the correction formulas were made for. thal.c calculates the thermodynamic parameters of dimers and hairpins. This calculations are adapted to primer design and not as complete as for example in the program UnaFold http://www.unafold.org/. If you need accurate melting temperatures of secondary structures please use UnaFold. Primer3 only uses this parameters to discriminate between primers and choose the better.

Best, Andreas