Vol. 21 no.

6 2005, pages 711–722

BIOINFORMATICS ORIGINAL PAPER doi:10.1093/bioinformatics/bti066

Structural bioinformatics

Comparison of different melting temperature calculation

methods for short DNA sequences
Alejandro Panjkovich and Francisco Melo∗
Departamento de Genética Molecular y Microbiología, Facultad de Ciencias Biológicas, Pontificia
Universidad Católica de Chile, Alameda 340, Santiago, Chile
Received on June 23, 2004; revised on September 10, 2004; accepted on October 5, 2004

Downloaded from https://academic.oup.com/bioinformatics/article/21/6/711/199347 by guest on 14 March 2021

Advance Access publication October 22, 2004

ABSTRACT to the authors or from our website http://protein.bio.puc.cl/melting-

Motivation: The overall performance of several molecular biology temperatures.html
techniques involving DNA/DNA hybridization depends on the accur- Contact: fmelo@bio.puc.cl
ate prediction of the experimental value of a critical parameter: the Supplementary information: The large set of oligonucleotides, the
melting temperature Tm . Till date, many computer software programs detailed results of the comparative and accuracy benchmarks, and
based on different methods and/or parameterizations are available for hundreds of comparative graphs generated during this work are avail-
the theoretical estimation of the experimental Tm value of any given able at our website http://protein.bio.puc.cl/melting-temperatures.html.
short oligonucleotide sequence. However, in most cases, large and
significant differences in the estimations of Tm were obtained while INTRODUCTION
using different methods. Thus, it is difficult to decide which Tm value is
The experimental performance and the outcome of several molecular
the accurate one. In addition, it seems that most people who use these
biology techniques depend on the accurate prediction of the DNA
methods are unaware about the limitations, which are well described
melting temperature (Tm ). This is particularly critical in the case
in the literature but not stated properly or restricted the inputs of
of those techniques that involve several oligonucleotides simultan-
most of the web servers and standalone software programs that
eously such as DNA microarrays of short sequences with fixed length
implement them.
or ‘Affymetrix chips’ (Pease et al., 1994), quantitative polymerase
Results: A quantitative comparison on the similarities and differences
chain reaction (PCR) (Buck et al., 1991) and multiplex PCR, where
among some of the published DNA/DNA Tm calculation methods is
large errors in the Tm estimation can lead to the amplification of non-
reported. The comparison was carried out for a large set of short oli-
specific products or to an inappropriate hybridization performance
gonucleotide sequences ranging from 16 to 30 nt long, which span
in general (Steger, 1994).
the whole range of CG-content. The results showed that significant
A variety of methods are available to calculate the Tm of short oli-
differences were observed in all the methods, which in some cases
gonucleotide sequences. The first approach described by Marmur and
depend on the oligonucleotide length and CG-content in a non-trivial
Doty (1962) uses a rather simple formula where the Tm depends only
manner. Based on these results, the regions of consensus and dis-
on the relative content of cytosine and guanine. This formula was later
agreement for the methods in the oligonucleotide feature space were
improved by adding a correction factor which is also responsible for
reported. Owing to the lack of sufficient experimental data, a fair and
salt concentration, thus adjusting the Tm value for different exper-
complete assessment of accuracy for the different methods is not
imental conditions (Wetmur, 1991). An indepth analysis of DNA
yet possible. Inspite of this limitation, a consensus Tm with minimal
oligonucleotides and their corresponding experimental Tm s has led
error probability was calculated by averaging the values obtained from
to the conclusion that not only the relative amounts of cytosine and
two or more methods that exhibit similar behavior to each particular
guanine concentrations determine the thermal denaturation of DNA,
combination of oligonucleotide length and CG-content class. Using a
but also the sequential arrangement of different nucleotides in DNA
total of 348 DNA sequences in the size range between 16mer and
sequences were found to play a major role in the experimental value
30mer, for which the experimental Tm data are available, we demon-
of Tm s. Hence the ‘Nearest Neighbor’ (NN) model was adopted for
strated that the consensus Tm is a robust and accurate measure.
the calculation of sequence related Tm s (Borer et al., 1974). The NN
It is expected that the results of this work would be constituted as
model postulates that the free energy for duplex formation depends
a useful set of guidelines to be followed for the successful experi-
mostly on two factors: first, the initiation-free energy given by an
mental implementation of various molecular biology techniques, such
unfavorable entropy as a consequence of a loss of translational free-
as quantitative PCR, multiplex PCR and the design of optimal DNA
dom after the first DNA/DNA pair is formed; and second, the sum
microarrays.
of the complementary pairwise terms between the oligonucleotide
Availability: A binary software distribution to calculate the consensus
sequences (propagation terms), which are based on dinucleotide
Tm described in this work for thousands of oligonucleotides simultan-
entities. In addition to these two terms, an entropic penalty is also
eously for the LINUX operating system is freely available upon request
considered, which takes into account the maintenance of symmetry
in self-complementary sequences. The calculation of Tm values
∗ To whom correspondence should be addressed. by the NN method requires a set of experimental thermodynamic

© The Author 2004. Published by Oxford University Press. All rights reserved. For Permissions, please email: journals.permissions@oupjournals.org 711
A.Panjkovich and F.Melo

parameters as an input. Till date, several tables with DNA/DNA ther- none of these web servers limits the input to be entered by the user.
modynamic parameters have been published (Gotoh and Tagashira, Thus, we believe that it is important to be aware of the current limit-
1981; Vologodskii et al., 1984; Breslauer et al., 1986; Delcourt and ations that these methods have and also to be aware of the magnitude
Blake, 1991; Doktycz et al., 1992; SantaLucia et al., 1996; Sugimoto of errors that could arise when using them for practical applications.
et al., 1996; Allawi and SantaLucia, 1997). Although a detailed Unfortunately, most people use these softwares with no apparent
comparison among these tables is difficult because several inde- conscience about the possible costs and risks, given the magnitude
pendent variables are involved in their derivation, attempts to clarify of errors that could arise from these calculations. It must be advised
the similarities and differences among the various NN parameters that depending on the method that is used, for a given specific oli-
have been carried out (Doktycz et al., 1992; SantaLucia, 1998). One gonucleotide sequence, the absolute differences in the calculated Tm
major conclusion of the SantaLucia’s (1998) work has been that a values could be large. Then, in practical terms, the question of which
remarkable consensus exists among most of the different thermo- Tm value should be considered arises. The major aim of this work
dynamic parameters, leading to the proposition and derivation of a is to compare the Tm values calculated by the different methods in

Downloaded from https://academic.oup.com/bioinformatics/article/21/6/711/199347 by guest on 14 March 2021

single and unified set of parameters. In this study, however, we will a large set of oligonucleotides that are representative of real-world
demonstrate that when different sets of NN parameters are used to applications in experimental biology. Thus, on the one hand, all the
predict the Tm for a large and diverse number of DNA sequences of oligonucleotides used in this study fall in the size range between
practical application value large differences are observed in many 16mer and 30mer, which is the widely used size range for the design
cases, which suggest that the current consensus parameters that have of PCR primers and microarrays of short-length oligonucleotides.
been proposed must be biased toward a reduced subset of possible Although it can be argued that in this size range the two-state trans-
DNA sequences. itions are not guaranteed and therefore the NN model should not be
When all the experimental conditions are fixed, the Tm value calcu- used, it must be stated that it is not the major aim of this work to
lated for any particular DNA sequence using the NN model depends assess or to develop a more accurate method, but to highlight the
on the thermodynamic table that is used. Till date, a single large-scale magnitude of the existing differences and variations among methods
benchmark comparing the predictive accuracy of 11 different ther- when they are applied to a typical DNA sequence that is used in the
modynamic tables has been published previously (Owczarzy et al., laboratory. On the other hand, the set of oligonucleotides used in
1998). In that work, a high-predictive accuracy with a small error was this study covers the complete space of CG-content, in contrast to
found only for three independent parameter sets: the set of Doktycz the oligonucleotides that have been used to derive the existing ther-
et al. (1992); the set of Allawi and SantaLucia (1997) and the set of modynamic parameters. Therefore, the present study constitutes an
Sugimoto et al. (1996). Although this large-scale accuracy assess- extensive and representative set of benchmark sequences to assess
ment constitutes a significant contribution, it must be mentioned that how similar the different predictions are.
the experimental data that are currently available to perform such The present study is divided into two major sections: the first sec-
benchmarks have several strong biases. First, there is a bias toward tion containing a comparative assessment of Tm s calculated using
oligonucleotide length, because some short-length oligonucleotides different methods and a second section containing an accuracy
(∼10mers long) are clearly overrepresented in the datasets. Unfor- assessment of different methods for experimental Tm prediction. It
tunately, most of the practical applications in molecular biology use must be noted that, in both cases, only DNA sequences with practical
oligonucleotides greater than this length. Second, most of the oligo- application values have been used. In the first section, we compare
nucleotides in the set fall in the CG-content range between 40 and the three methods: the basic, salt adjusted and thermodynamic. The
60%, with a few representations at the most extreme cases. Third, and first objective was to evaluate if simple models such as the basic and
most importantly, a large fraction of the oligonucleotides for which salt adjusted methods could give similar Tm predictions compared
experimentally determined Tm s are available have been already used to the ones obtained by the more complex thermodynamic method.
to obtain or derive the same thermodynamic parameters with the help The second objective was to compare in detail the similarities and
of which performance and accuracy are going to be assessed. Finally, differences of the Tm predictions obtained by the thermodynamic
not all of these oligonucleotides available are known to melt in a two- method while using three different parameter tables. Out of the dif-
state manner, which is a primary condition in order to apply the NN ferent thermodynamic tables that have been published, three tables
model to calculate the Tm . These existing biases can be explained were derived from the analysis of optical melting curves for a vari-
on the basis that the larger an oligonucleotide is, greater the chances ety of short synthetic DNA duplexes and are the most commonly
that it will not melt in a two-state manner, and hence it will be of used by the scientific community as expressed in the implementa-
marginal interest for the derivation or validation of thermodynamic tions of several standalone software programs and web browsers on
parameters to be used in the NN model. the Internet (for a list of web servers along with their parameter-
Therefore, there are still several limitations that need to be con- izations, see Supplementary material). These three NN tables were
sidered when using these methods to estimate the experimental Tm used in this comparative study and include the original set of values
of any given oligonucleotide for a practical application. Surprisingly, reported by Breslauer et al. (1986), and the putatively improved sets
the majority of the current web servers and software programs that are of values reported by SantaLucia et al. (1996) and Sugimoto et al.
available to calculate the Tm of oligonucleotides do not inform about (1996). It could be argued that Breslauer set is not as accurate as
these limitations and do not limit the input provided by the user, thus those of SantaLucia and Sugimoto sets based on the results of the
any valid DNA/RNA sequence is allowed. Out of the 17 web servers previous accuracy benchmarks mentioned above. However, till date,
that we have found as freely available on the Internet (for a list of not a single comparison has been performed among these thermo-
these web servers see the Supplementary material), only 4 of them dynamic sets using a large, representative and unbiased set of DNA
provide some sort of warning about the limitations of the Tm calcu- sequences of practical application in molecular biology. Moreover,
lation methods. However, at the time of performing the calculations, we will demonstrate here that all these thermodynamic sets, when

712
 Comparison of melting temperature calculation methods

compared in a pairwise fashion, exhibit some differences and share buffered solution of 50 mM Na+ and 50 nM of oligonucleotide concentration,
a similar behavior at different regions of the oligonucleotide fea- with a pH close to 7.0, but the Tm of DNA is unaffected within a significant
ture space (represented by length and CG-content). In the second range of pH around 7.0 due to the lack of titratable groups close to this pH
section of this study, the third and last objective of this work is in the Watson–Crick paired DNA. The salt adjusted Tm calculations were
addressed, which consists of the accuracy assessment of these meth- performed using the following equation (Howley et al., 1979):


ods in predicting the experimental Tm of several DNA sequences. yG + zC − 16.4
Tm = 100.5 + 41.0 ×
This accuracy benchmark set contained all DNA sequences in the wA + xT + yG + zC
length range between 16mer and 30mer for which the experimental

820.0
Tm , salt concentration and oligonucleotide concentration values were − + 16.6 log([Na+ ])
wA + xT + yG + zC
available. Also, a consensus Tm value proposed in this study was also
where x, y, w and z are the number of the bases of T, G, A and C, respectively.
assessed, thereby giving the lowest average error. The consensus Tm
In the above equation, the second term adjusts for the GC-content and the
value calculation is based on the consensus Tm map among different
third term adjusts for the length of the sequence. The equation also assumes

Downloaded from https://academic.oup.com/bioinformatics/article/21/6/711/199347 by guest on 14 March 2021

methods observed from the comparative assessment. that the annealing occurs under standard conditions of pH close to 7.0 and 50
To achieve the comparative assessment, we have generated a large nM of oligonucleotide concentration.
set of oligonucleotide sequences in a computer, calculated the Tm The melting temperatures are calculated using the nearest-neighbor model
value for each of them by using different methods, and assessed and thermodynamic data as previously described by SantaLucia et al. (1996).
the observed similarities, differences and correlations among pre- The equation used is as follows:
dictions. These comparisons were performed independently for each 
(
Hd ) +
Hi
combination of oligonucleotide length and CG-content, for a given Tm =  + 16.6 log[Na+ ],
(
Sd ) +
Si +
Sself + R × ln CbT
range of oligonucleotide length that is of practical value. The regions
of consensus and disagreement are highlighted. The results show that where sums of enthalpy (
Hd ) and entropy (
Sd ) are calculated over all
significant differences could be obtained while using different meth- internal nearest-neighbor doublets,
Sself is the entropic penalty for self-
ods and that a consensus Tm value with a minimal error probability complementary sequences, and
Hi and
Si are the sums of initiation
enthalpies and entropies, respectively (Table 1). R is the gas constant (fixed
should be defined. Owing to the lack of experimental data available
at 1.987 cal/K · mol), CT is the total strand concentration in molar units and
till date, which is not extensive or representative enough, we address Tm is the melting temperature given in Kelvin units. Constant b adopts the
this point by averaging the Tm values of these methods that consist- value of 4 for non-self-complementary sequences or equal to 1 for duplexes
ently exhibited a similar behaviour for each particular combination of self-complementary strands or for duplexes when one of the strands is in
of oligonucleotide length and CG-content. To support our model, significant excess. The thermodynamic calculations assume that the anneal-
an accuracy benchmark using all DNA sequences that fall within ing occurs in a buffered solution at pH near 7.0 and that a two-state transition
the size range covered in this study for which experimental Tm data occurs.
are available was performed, demonstrating that the consensus Tm For the NN Tm calculations, three different thermodynamic tables were
proposed here constitutes a robust and accurate measure. Finally, a used in this work (summarized in Table 1): the first one is the original table
list of guidelines to calculate the Tm of short DNA sequences with a published by Breslauer et al. (1986), the second table was published by
SantaLucia et al. (1996) and the last table included in this study was published
minimal error probability is provided at the end of this paper.
by Sugimoto et al. (1996). Therefore, five different Tm values were calculated
and compared for each DNA oligonucleotide: one basic term (denominated
METHODS ‘bas’), one salt adjusted term (denominated ‘sal’) and three NN thermody-
Set of artificial oligonucleotide sequences used in the namic sets [denominated by date of publication: ‘Th1’ for Breslauer et al.
(1986), ‘Th2’ for SantaLucia et al. (1996) and ‘Th3’ for Sugimoto et al.
comparative benchmark (1996)]. To perform a comparison of all the methods, the Tm calculations
A total of 300 000 DNA oligonucleotide sequences were randomly generated were carried out by fixing the oligonucleotide and sodium concentrations
in a computer. The length of these random DNA sequences was restricted at 50 nM and 50 mM, respectively. In the case of the accuracy benchmark
to be in the range of 16–30 nt. The random sequences were also generated reported at the end of this work, the Tm calculations were carried out using
and selected in such a way that they span a homogeneous distribution in the the specific experimental oligonucleotide and salt concentrations reported for
defined size range and additionally in the possible CG-content range. For each particular case. The detailed experimental and calculated data used in
this purpose, 10 uniform CG-content classes were defined (fixed intervals this accuracy benchmark is available as Supplementary material.
of magnitude equal to 10, ranging between 0 and 100% of CG-content).
Therefore, each particular combination of sequence length and CG-content Salt correction in melting temperature calculations
class was populated with the same total number of 2000 oligonucleotide As described in the previous section, a unique salt correction term was used in
sequences (i.e. 15 length classes×10 CG-content classes×2000 sequences = this comparative benchmark for the thermodynamic and salt adjusted meth-
300 000 total sequences). ods. This salt correction factor was reported by Schildkraut and Lifson (1965)
Melting temperature calculations for large DNA polymers and is still the most frequently used correction term
in the current implementations available for Tm calculations. Little is known
In this study, three different methods were used to calculate and compare the about the dependence of Tm s for short DNA sequences on salt concentration,
Tm s of short DNA oligonucleotides: basic, salt adjusted and NN thermody- particularly for bivalent ions. A recent and complete study has addressed
namic calculations. The basic Tm calculations were performed according to the effects of sodium ions on experimental Tm s for short DNA oligomers
the following equation (Marmur and Doty, 1962): (Owczarzy et al., 2004). That work contains an exhaustive and complete

 review about salt correction formulas that have been published in the literat-
yG + zC − 16.4
Tm = 64.9 + 41.0 × , ure, and dispatches a new and more accurate salt correction factor. However,
wA + xT + yG + zC
it must be stated that the present comparative study is not influenced by the
where x, y, w and z are the number of the bases of T, G, A and C, respectively. above salt correction factor, because all the calculated relationships of Tm s
This equation assumes that the annealing occurs under standard conditions in a among these methods are relative and not absolute. Thus, the differences

713
A.Panjkovich and F.Melo

Table 1. Thermodynamic parameters for DNA helix initiation and propagation in 1 M NaCl

Author Breslauer et al. (1986) SantaLucia et al. (1996) Sugimoto et al. (1996)
Abbreviation TH1 TH2 TH3
Propagation sequence
H
S
G
H
S
G
H
S
G
(kcal/mol) (cal/◦ K mol) (kcal/mol) (kcal/mol) (cal/◦ K mol) (kcal/mol) (kcal/mol) (cal/◦ K mol) (kcal/mol)

AA/TT −9.1 −24.0 −1.9 −8.4 −23.6 −1.02 −8.0 −21.9 −1.2
AT/TA −8.6 −23.9 −1.5 −6.5 −18.8 −0.73 −5.6 −15.2 −0.9
TA/AT −6.0 −16.9 −0.9 −6.3 −18.5 −0.60 −6.6 −18.4 −0.9
CA/GT −5.8 −12.9 −1.9 −7.4 −19.3 −1.38 −8.2 −21.0 −1.7
GT/CA −6.5 −17.3 −1.3 −8.6 −23.0 −1.43 −9.4 −25.5 −1.5

Downloaded from https://academic.oup.com/bioinformatics/article/21/6/711/199347 by guest on 14 March 2021

CT/GA −7.8 −20.8 −1.6 −6.1 −16.1 −1.16 −6.6 −16.4 −1.5
GA/CT −5.6 −13.5 −1.6 −7.7 −20.3 −1.46 −8.8 −23.5 −1.5
CG/GC −11.9 −27.8 −3.6 −10.1 −25.5 −2.09 −11.8 −29.0 −2.8
GC/CG −11.1 −26.7 −3.1 −11.1 −28.4 −2.28 −10.5 −26.4 −2.3
GG/CC −11.0 −26.6 −3.1 −6.7 −15.6 −1.77 −10.9 −28.4 −2.1
Any G·C pair? 0.0 −16.77 +5.0 0.0 −5.9 +1.82 +0.6 −9.0 +3.4
Only A·T pairs? 0.0 −20.13 +6.0 0.0 −9.0 +2.8 +0.6 −9.0 +3.4
Symmetry correction 0.0 −1.34 +0.4 0.0 −1.4 +0.4 0.0 −1.4 +0.4
5 -terminal-T·A-3 bp 0.0 0.0 0.0 +0.4 0.0 +0.4 0.0 0.0 0.0

or correlations among Tm s calculated by any two methods are not sensitive concentrations, accounting for a total of 348 data points. All these sequences
to this factor, as it is a constant value on both sides of the equation that is have a length ranging between 16mer and 30mer. A total of 37 unique
subtracted and therefore eliminated. sequences were extracted from Owczarzy et al. (1998); a total of 11 unique
In the case of accuracy benchmark reported at the end of this work, the sequences were extracted from the NTDB database (Chiu et al., 2003); and
salt correction factor of course becomes relevant, because the calculated Tm finally, a total of 60 unique sequences were extracted from Owczarzy et al.
is then compared to the experimentally observed temperature. Thus, for the (2004). The experimental Tm for each one of these unique 60 sequences
accuracy benchmark we have used not only the correction factor described was measured at five different salt concentrations, thus constituting a total
above, but also the new and improved correction factor recently published of 300 different experimental points (Owczarzy et al., 2004). Therefore,
by Owczarzy et al. (2004). In addition, we also used the correction factor a total of 348 data points were used in the accuracy benchmark reported
published by SantaLucia et al. (1996), which is slightly smaller than the one at the end of this work. The detailed dataset information is available as
described above. The best-performing salt corrections were used to assess Supplementary material at our website http://protein.bio.puc.cl/melting-
each method. temperatures.html.

Comparative measures Electronic availability of software, raw and calculated

Several measures of similarity between the Tm values reported for any
data
two methods were used in this study. All comparisons were done within The software used in this work was fully written by us in ANSI C lan-
each grid point (i.e. for each fixed combination of oligonucleotide length guage and the executable LINUX version is freely available upon request
and percentage of CG-content class), thus involving a total of 2000 oligo- or from our website. The complete set of oligonucleotides used in this work,
nucleotide sequences. The calculated measures include the maximal observed along with the experimental and calculated data, and the complete set of
absolute difference (MaxAD) of the 2000 pairwise comparisons between the generated comparative graphs in full color are available at our website
any two methods, the minimal observed absolute difference (MinAD), the http://protein.bio.puc.cl/melting-temperatures.html.
average absolute difference (AveAD), the standard deviation of absolute dif-
ferences (DevAD), the maximal observed difference (MaxD), the minimal
observed difference (MinD), the average observed difference (AveD), the RESULTS
correlation coefficient (CC) and the percentage of cases, where the abso- In this study, thousands of short DNA sequences were generated
lute difference between Tm values was equal to or less than 10 (Per10C), 5 in a computer, and the Tm for each of them was calculated using
(Per5C) and 3 (Per3C)◦ C. The correlation coefficient was only calculated
five different approaches. Several comparative measures were used
between thermodynamic methods because the other two methods have a
to assess the differences and similarities of the calculated Tm s. The
variance equal to zero for oligonucleotides of fixed length and CG-content
(with the salt concentration being fixed, they depend only on these two approaches included the basic Tm calculation (denominated ‘bas’),
variables, which were of course identical at each grid point where the the salt adjusted Tm calculation (denominated ‘sal’) and the NN ther-
analysis was carried out). Because of space constraints, only some of the modynamic method based on three of the most commonly used
comparative measures are reported in this paper. The complete comparat- parameter sets. The thermodynamic parameter sets included the
ive data is available graphically as Supplementary material at our web site Breslauer table (Breslauer et al., 1986) denominated here as Th1,
http://protein.bio.puc.cl/melting-temperatures.html. the SantaLucia table (SantaLucia et al., 1996) called Th2 and the
Sugimoto table (Sugimoto et al., 1996) denominated as Th3. These
Accuracy benchmark dataset parameters are listed in Table 1 and the mathematical expressions
The experimental data used here were extracted from different sources and used to calculate the Tm s are described in detail in the Methods
consisted of a total of 108 unique oligonucleotide sequences in various salt section.

714
 Comparison of melting temperature calculation methods

The length of the DNA sequences was limited between 16 and based on the observed differences of Tm predictions offered by these
30 nt, which is the most commonly used length range for PCR primer two methods, where Breslauer Tm values were consistently higher
design and in situ synthesized oligonucleotide microarrays. For each than SantaLucia estimations, in the whole range of sequence length
length, 10 different CG-content classes ranging between 0 and 100 and CG-content (data not shown). The comparison of Breslauer
were defined, thus covering the complete CG-content range. Finally, (Breslauer et al., 1986) or Th1 and Sugimoto (Sugimoto et al., 1996)
a total of 2000 DNA sequences were randomly generated for each or Th3 calculations exhibit high similarity in the CG-content range
particular combination of length and CG-content class. For each of 10–40%, irrespective of the sequence length (Fig. 2C). These two
sequence, the Tm s were calculated using the methods described above approaches presented a good correlation (Fig. 2D), but it is puzz-
and several comparisons were carried out (see Methods section). ling that the highest correlation values are not expressed at the same
The first comparison involved the maximal observed absolute points where the highest similarities are observed (Fig. 2C and D).
difference of Tm s among all the methods at each combination of Finally, the comparison of SantaLucia (SantaLucia et al., 1996) or
sequence length and percentage of CG-content. The results are shown Th2 and Sugimoto (Sugimoto et al., 1996) or Th3 showed simil-

Downloaded from https://academic.oup.com/bioinformatics/article/21/6/711/199347 by guest on 14 March 2021

in Figure 1A. It is clear that large differences are observed, ran- arities only for medium to rich CG sequences, irrespective of the
ging between 16 and 24◦ C. The sequences in the middle range of length (Fig. 2E). However, the same previously described effect was
CG-content present the smallest differences, irrespective of their observed, the correlation is high only for the low-similarity regions
length. When the same comparison is carried out, but only for the (Fig. 2F).
thermodynamic methods (Fig. 1B), the absolute differences remain To illustrate the magnitude of the observed disparity between
in a similar range, but with the maximal values slightly lower at similarities and correlation while comparing two thermodynamic
20◦ C. However, the difference dependencies are totally different. In methods, two particular grid points of the comparison between
this case, sequences with a larger fraction of C and G nucleotides SantaLucia (SantaLucia et al., 1996) and Sugimoto (Sugimoto et al.,
exhibit larger differences in the Tm estimations, almost irrespective 1996) calculations were chosen. It should be noted that these two
of sequence length. When the average absolute differences were plot- grid points are not the most extreme points in the graph. On the one
ted, the previous trends were perfectly conserved, but the variation hand, it is clearly demonstrated that the correlation is high for a grid
range decreased as expected (Fig. 1C and D). The average differ- point where the similarity is low (Fig. 3A). On the other hand, it is
ences range between 5 and 9◦ C for all the methods and between demonstrated that the correlation is low for a grid point where the
4 and 12◦ C for the thermodynamic methods. Average differences similarity is high (Fig. 3B).
were smaller for sequences in the CG-content range of 40–60%, but In practical terms, irrespective of the correlation and similarities
still above 4–5◦ C. The average differences slightly increased with among different methods, what is desired is to have the possible most
sequence length, irrespective of the CG-content. accurate Tm estimation. This simple study has clearly shown that the
Another comparison was a similarity measure that consisted in methods do not correlate nor exhibit a clear common pattern of sim-
the percentage of cases or oligonucleotide sequences where all the ilarity. Thus, it is quite difficult to address the question of which
methods shared Tm estimations with a maximal difference of a fixed method or approach one should follow in order to be successful at
figure in Celsius degrees. When all the methods are compared simul- obtaining an accurate Tm prediction for a wide and diverse range
taneously, only a small fraction of the plotted area exhibit significant of oligonucleotide sequences. The best solution available, given the
similarities within 10◦ C. This area was restricted to short sequences biased selection of oligonucleotide sequences that has been used to
of 16mer with a CG-content between 40 and 80% (Fig. 1E). When obtain and parameterize the current existing thermodynamic tables,
the thermodynamic methods are compared using this measure, sim- is to minimize the possible error that can arise when estimating the
ilar results are observed, but the similarity area increases if the same Tm of a given DNA oligonucleotide. Owing to the lack of suffi-
figure value is considered (Fig. 1F). These results demonstrate that cient representative experimental data, this can only be achieved by
large and significant differences exist among different methods in following the average Tm of such methods that exhibited a similar
the Tm estimations. High similarities are only observed for short behavior at a particular combination of sequence length and CG-
oligonucleotide sequences with a medium and restricted CG-content. content. The observed consensus among thermodynamic sets that
To assess the possibility that a particular method could be con- has been obtained in this study is illustrated in Figure 4A. It can
cealing some additional existing similarities among the methods, a be seen that consensus among all the NN sets is only observed for
pairwise comparison of the thermodynamic methods was carried out. short sequences in the middle CG-content range (Fig. 4A, white
The results of these comparisons are shown in Figure 2. The com- color). This is not surprising because almost all the sequences that
parison of the Breslauer set (Breslauer et al., 1986) or Th1 and the have been used to obtain the thermodynamic parameters used in
SantaLucia set (SantaLucia et al., 1996) or Th2 (Fig. 2A) demon- this study lie exactly in this area (Breslauer et al., 1986; SantaLucia
strates that the overall similarities observed earlier among all the et al., 1996; Sugimoto et al., 1996). However, two additional
thermodynamic methods (Fig. 1F) were a simple consequence of regions where only two NN sets overlap are also highlighted. These
the existing similarities among these two methods, which present the showed the similarity regions of Breslauer (Breslauer et al., 1986)
largest differences in the all-against-all pairwise comparisons among and Sugimoto (Sugimoto et al., 1996) on the one hand (Fig. 4A,
the three thermodynamic methods (Fig. 2). However, it is quite light gray color), and that of SantaLucia (SantaLucia et al., 1996)
surprising that despite the low similarity observed, the Tm estim- and Sugimoto (Sugimoto et al., 1996) on the other hand (Fig. 4A,
ations of these two methods were found to be highly correlated in a dark gray color). In addition to these zones, some dangerous regions
large fraction of the oligonucleotide feature space (Fig. 2B). In fact, where none of these methods exhibits a similar behaviour are also
these two thermodynamic methods are the most dissimilar ones, but illustrated (Fig. 4A, black color). These regions mostly occur at
indeed the most highly correlated for a large fraction of the oligonuc- extreme CG-content values, almost irrespective of the sequence
leotide sequences tested in this study. This can be explained partially length.

715
A.Panjkovich and F.Melo

Downloaded from https://academic.oup.com/bioinformatics/article/21/6/711/199347 by guest on 14 March 2021

Fig. 1. Simultaneous comparison of Tm s calculated using all the methods. For each class of oligonucleotide length and percentage of CG-content, 2000
DNA sequences were randomly generated in a computer. The Tm of each sequence was calculated using different methods: basic, salt adjusted and the NN
thermodynamic method based on three published tables (Table 1). Subsequently, the similarities and differences of the calculated values were evaluated. (A, C
and E) The simultaneous comparison of all the methods is shown. (B, D and F) The simultaneous comparison of the three thermodynamic methods is shown.
The maximum observed absolute differences are shown in A and B; the average absolute differences in C and D; and finally, the percentage of cases where the
absolute difference is ≤10 or 5◦ C is shown in E and F. In the case of the simultaneous comparison of all the thermodynamic methods, a threshold of 5◦ C was
used to define similarity, because in our judgement this figure represents a reasonable error estimation of Tm values, as it has been suggested previously. In the
case of the simultaneous comparison of all methods, a threshold value of 10◦ C was used because no similarity was observed below that value (i.e. a flat graph
in the XY plane was generated).

716
 Comparison of melting temperature calculation methods

Downloaded from https://academic.oup.com/bioinformatics/article/21/6/711/199347 by guest on 14 March 2021

Fig. 2. Correlation and percentage of similarities among thermodynamic parameter sets. The pairwise similarities and correlation of the Tm values calculated
by the thermodynamic sets using the parameters described in Table 1 were assessed. The procedure adopted was the same as described in the legend of Figure 1.
(A, C and E) The percentage of oligonucleotide sequences where the absolute Tm difference is ≤5◦ C is shown as a function of sequence length and CG-content.
(B, D and F) The correlation coefficient of the calculated Tm s is shown as a function of sequence length and percentage of CG-content. The corresponding
pairwise comparisons are as follows: (A and B) Th1 versus Th2; (C and D) Th1 versus Th3 and (E and F) Th2 versus Th3. Th1 stands for Breslauer et al.
(1986), Th2 for SantaLucia et al. (1996) and Th3 for Sugimoto et al. (1996).

717
A.Panjkovich and F.Melo

Downloaded from https://academic.oup.com/bioinformatics/article/21/6/711/199347 by guest on 14 March 2021

Fig. 3. Correlation of Tm estimations. The correlation of the Tm estimations by two thermodynamic parameter sets are shown at two specific grid points of
Figure 2F. The NN sets compared were that of SantaLucia (SantaLucia et al., 1996) and Sugimoto (Sugimoto et al., 1996) as shown in Figure 2E and F. (A)
Scatter plots of Tm values estimated by SantaLucia and Sugimoto for 2000 sequences of length equal to 18mer and CG-content ranging between 20 and 30%.
The observed correlation coefficient value is 0.821994. (B) Scatter plots of Tm values estimated by SantaLucia and Sugimoto for 2000 sequences of length
equal to 23mer and CG-content ranging between 70 and 80%. The observed correlation coefficient value is 0.006366.

In summary, four different regions or zones in oligonucleotide giving a total of 348 experimental data points to support the accur-
feature space were obtained: (1) Zone 3, where Th1, Th2 and Th3 acy benchmark. The frequency of occurrence for all these sequences
exhibit similar Tm values (Fig. 4A, white color); (2) Zone 1, where mapped onto oligonucleotide feature space is shown in Figure 4C.
Th1 and Th3 exhibit similar Tm values (Fig. 4A, light gray color); The accuracy benchmark set was performed by using the salt cor-
(3) Zone 2, where Th2 and Th3 exhibit similar Tm values (Fig. 4A, rection factor that yielded the best results for each method, except
dark gray color); and (4) Zone 0, where none of the NN parameter for the consensus method, which used the Tm from each method
sets share similar Tm values (Fig. 4A, black color). Based on these calculated with the salt correction factor reported in the Methods
comparative results, a consensus Tm is then defined as the average section. Although this is unfair for the consensus method, the reason
of the Tm values calculated using those NN sets that exhibited a for doing so was to obtain values consequent with the comparat-
similar behavior at a given length and percentage of CG-content. ive benchmark that generated the consensus map, which used this
Thus, depending on the grid point where a given oligonucleotide unique and fixed salt correction factor for all the methods (salt adjus-
maps, the Tm values that are considered and averaged. Hence, the ted, Th1, Th2 and Th3). Three different salt correction factors were
consensus Tm value of oligonucleotides mapping in Zone 3 would be calculated and applied to the Tm predictions of each individual Tm
the average of the Tm values calculated by the Th1, Th2 and Th3 sets. calculation method. First, the salt correction described in the Meth-
In the case of a DNA sequence falling in Zone 1, the consensus Tm ods section (Schildkraut and Lifson, 1965), where the salt adjusted
value would be the average of the Tm values obtained by Th1 and Th3 method and the NN model with the Th1 and Th3 parameter sets
sets. Finally, the consensus Tm value of a sequence falling in Zone 2 showed the best performance. Second, the salt correction factor sug-
would be the average of the Tm values obtained by Th2 and Th3 sets. gested by SantaLucia et al. (1996), where the NN model with the
When a sequence maps into Zone 0 (black regions of Fig. 4A), it is Th2 parameter set showed the best performance. None of the meth-
not clear which parameter sets should be used. Although an average ods performed well when the new salt correction factor recently
of all Tm values could be used for these regions where no similarities suggested by Owczarzy et al. (2004) was used. In fact, when this cor-
are observed, we suggest that oligonucleotides falling in the black rection factor was considered, the accuracy of different predictions
regions should be avoided, because it is not clear which NN set could was significantly reduced in all the cases.
be considered. The results of the benchmark are summarized in Table 2. The
A final objective of this work was to assess the accuracy of dif- best-performing method in terms of giving the closest value to the
ferent methods at predicting the experimental Tm values for DNA experimental Tm most of the times is Th2 (SantaLucia et al., 1996).
sequences with a length between 16mer and 30mer. For that pur- In 40% of the cases, this method gives the most accurate experimental
pose, we recollected all the unique DNA sequences in that size range Tm prediction. However, when a more useful accuracy measure for
for which the experimental Tm s were available. A total of 108 unique the experimental biologist is used to assess the performance of these
oligonucleotide sequences that fulfill these requirements have been methods, the overall picture changes. When the percentage of cases
considered, for which frequencies of occurrence were mapped on oli- where the predictions containing a maximal error of 3 or 5◦ C is
gonucleotide feature space and are displayed in Figure 4B. It must be used, the best performance is achieved by the consensus method
mentioned that 60 of these 108 sequences had a total of five exper- proposed in this study, closely followed by the NN model using the
imental Tm values determined at different salt concentrations, thus thermodynamic parameters of Th3 (Sugimoto et al., 1996) and Th2

718
 Comparison of melting temperature calculation methods

(SantaLucia et al., 1996). The same trend is observed when the aver-
age error is used as the accuracy measure. It must be noted that when
the benchmark is carried out considering only those sequences that
map into a consensus area of the oligonucleotide feature space (i.e.
excluding those sequences that fall into a black region of Fig. 4A),
the observed trend is not only confirmed, but also reinforced (Table 2,
bottom values within each cell).
Statistical significance tests for different methods’ average
accuracies were performed (Table 3). No significant average accur-
acy differences were observed among Th2, Th3 and consensus
methods when all the sequences from the benchmark set are used.
These three methods are however more accurate than the basic, salt-

Downloaded from https://academic.oup.com/bioinformatics/article/21/6/711/199347 by guest on 14 March 2021

adjusted and Th1, with statistical significance. When only those 281
sequences that fall in the consensus regions of the consensus map
are used (Fig. 4C), the situation is different. In this case, Th3 and
consensus methods are more accurate than the others, with statistical
significance. These results demonstrated that the Tm s estimated by
the Th3 NN set and the consensus method are the most accurate.
Thus, the consensus Tm is a robust measure that exhibits a minimum
average error and the largest fraction of cases predicted within an
acceptable experimental error.

DISCUSSION
In this study, we have compared the similarities and correlations of
the Tm values calculated using different methods. For this purpose,
we have used a large and representative benchmark set of short oli-
gonucleotide sequences. We did not address completely the problem
of judging which method gives the closest value to the experimental
Tm . We have only performed an accuracy benchmark using all the
relevant sequences for which experimental Tm data are available to
support the idea of using a consensus Tm with a minimal error prob-
ability. But we still believe that a benchmark based on the currently
available sequences is not sufficient to validate or discard the per-
formance of any existing method for Tm estimation. Although this last
issue has been addressed by other authors previously and an effort to
reconcile the existing differences has been made (SantaLucia, 1998;
Rouzina and Bloomfield, 1999; Owczarzy et al., 1998), we believe
that this kind of assessment is biased because most of the experi-
mental data available correspond to very short DNA sequences and

Fig. 4. Consensus of Tm values among thermodynamic parameter

sets. The consensus among two or three parameter sets is defined
when at least 80% of the sequences exhibit an absolute difference
between the calculated Tm values <5◦ C. All possible pairwise com-
parisons were carried out, as well as simultaneous comparison of the
three thermodynamic sets. Th1 stands for Breslauer et al. (1986);
Th2 stands for SantaLucia et al. (1996) and Th3 stands for
Sugimoto et al. (1996). Th1 and Th2 did not show similar behavior
in the whole range of sequence length and percentage of CG-content.
(A) The observed consensus among the methods is as follows: Sim-
ultaneously, Th1 and Th3, Th2 and Th3, exhibit similar values (white
color); only Th1 and Th3 exhibit similar values (light gray color);
only Th2 and Th3 exhibit similar values (dark gray color) and finally,
no consensus is observed among any of the methods (black color).
(B) At each grid point, the total number of unique oligonucleotide
sequences for which existing experimental Tm data available are also
displayed (Legend to Table 2). (C) At each grid point, the total
number of different cases (i.e. a unique combination of oligonuc-
leotide sequence, salt and oligonucleotide concentration) for which
the existing experimental Tm data available are also displayed.

719
A.Panjkovich and F.Melo

Table 2. Accuracy benchmark of methods of the experimental Tm s (5◦ C) for several oligonucleotides that were
not following a two-state transition. Although it is still not clear if
the NN model could be a valid approximation for larger sequences,
Accuracy measure BAS SAL BRE SAN SUG CON
where long range interactions and salt dependence could have a com-
plex effect, most of the computational implementations currently
BEST (%) 0.6 5.2 3.5 40.8 26.2 23.9 available on the Internet or standalone software use this method for
(0.0) (6.4) (2.1) (38.1) (27.1) (26.3)
sequences that fall in the length range covered in this study. Thus, we
ERROR WITHIN 5◦ C (%) 11.2 31.0 26.2 83.3 83.6 83.9
(9.6) (35.2) (24.6) (84.0) (84.3) (86.1)
believe that it is important to be aware of the potential errors and/or
ERROR WITHIN 3◦ C (%) 3.7 14.9 14.4 60.9 60.1 61.5 existing variations in the Tm predictions that will be generated by
(2.9) (17.4) (12.5) (57.3) (62.6) (64.1) using those softwares without care.
AVERAGE ERROR (◦ C) 12.3 7.1 8.5 2.9 2.9 2.8 Our results showed that complex relationships exist among meth-
(12.6) (6.7) (8.6) (3.0) (2.7) (2.6) ods. For instance, as it was shown in Figure 2, the similarity

Downloaded from https://academic.oup.com/bioinformatics/article/21/6/711/199347 by guest on 14 March 2021

and correlation of Tm values among the methods are usually not
A total of 348 DNA sequences 16–30mers long with experimental Tm , salt and oligo- aligned or in phase as may be expected. It was observed quite
nucleotide concentrations available were used in this benchmark: 37 sequences were often that when two methods give similar predictions, they do not
obtained from the work of Owczarzy et al. (1998); 11 sequences were obtained from
the NTDB database (Chiu et al., 2003); and the remaining 300 sequences were obtained
correlate. In addition, the pairwise comparison of average Tm differ-
from Owczarzy et al. (2004). The complete table containing all the experimental val- ences between SantaLucia (SantaLucia et al., 1996) and Sugimoto
ues and the theoretical predictions made by using the various methods is available as (Sugimoto et al., 1996) revealed that SantaLucia predicts higher Tm
Supplementary material. The Tm s were predicted with the basic method (BAS), the values for sequences in the middle to high-CG-content range, irre-
salt adjusted method (SAL), and the NN model with the thermodynamic parameters
spective of the sequence length (Supplementary material). In the
of Breslauer et al. (1986) (BRE), SantaLucia et al. (1996) (SAN) and Sugimoto et al.
(1996) (SUG). The Tm was also predicted using the consensus method (CON) proposed low-CG-content range, the opposite result was observed. On the
in this study, which is based on the results obtained and shown in Figure 4. The consensus other hand, Breslauer (Breslauer et al., 1986) Tm predictions were
Tm corresponds to the average Tm of those methods that exhibit similar results at a given higher than SantaLucia and Sugimoto Tm estimations, in the whole
grid point of the oligonucleotide feature space. In those cases where no similarities are range of sequence length and CG-content, but with a different mag-
observed among methods (black regions of Fig. 4A), the average of all melting temper-
ature values was used (top values within each cell). The results of the benchmark using
nitude and behavior (Supplementary material). When the standard
the 281 sequences that are mapped in Zones 1, 2 and 3 (excluding the 67 sequences from deviations of the Tm predictions were calculated for each grid point, it
the black regions in Fig. 4C) are shown in this table within parentheses. Four different was surprising to note that Breslauer presented a large figure (>4◦ C)
accuracy measures are reported here. First, the percentage of cases where the method in the middle range of the CG-content and for short length sequences
gives the closest prediction to the experimental Tm (BEST); second, the percentage of
(Supplementary material). Perhaps this is the reason why Breslauer
cases where the method gives a prediction within 5 and 3◦ C from the experimental Tm
(ERROR WITHIN); and finally, the average of the absolute differences between the Tm predictions do not agree very well with the predictions based
prediction method and the experimental Tm for all the cases considered (AVERAGE on other parameter sets (SantaLucia et al., 1996). SantaLucia and
ERROR). Sugimoto Tm predictions showed low-standard deviations at each
grid point (the maximal values were ∼2◦ C). However, when the
minimum absolute difference Tm value was plotted at each grid point
have also been used to optimize and parameterize the existing meth- for all pairwise comparisons among thermodynamic methods, large
ods. In addition, it is also true that the currently available data are differences were observed in all cases (Supplementary material). As
not representative of the oligonucleotide sequence space either. Most an example of these observations, Sugimoto and Breslauer showed
of the experimental data include oligonucleotide sequences <16mer large differences in the high-CG-content range (>12◦ C). Breslauer
and with a CG-content in the range of 40–60%. The aim of this and SantaLucia showed a similar trend, but in this case the differ-
study has not been to disqualify any of the existing methods, but ences were >4◦ C. SantaLucia and Sugimoto showed the opposite
to demonstrate that significant differences in the Tm predictions of trend, with the largest differences occurring in the low-CG-content
short DNA sequences are observed among them when a large num- range. In this case, the minimum differences were >3◦ C.
ber of sequences of practical application value are tested. In practical As a complement of the comparative assessment, we finally per-
terms, a large ‘error’ in the estimation of the Tm could easily arise, formed an accuracy benchmark to evaluate the ability of different
irrespective of which method is used. Thus, we believe that this com- methods to predict the experimental Tm of oligonucleotide sequences
parative analysis will provide some guidelines to be followed in order with a practical application value, under varying conditions of salt
to avoid or minimize large and frequent errors in the estimation of and oligonucleotide concentrations. It must be mentioned that this
the experimental Tm of short oligonucleotide sequences. Also, we benchmark was not fair in one major respect: all the methods were
encourage people working in the experimental determination of Tm s compared under such conditions which did not apply to some of them,
to cover more extensively the practical oligonucleotide feature space such as the basic method under varying conditions of salt and oligo-
when deriving new and improved thermodynamic tables. nucleotide concentration, and the salt adjusted method under varying
It has been suggested that NN thermodynamics apply to duplexes conditions of oligonucleotide concentration. This is one of the reas-
ranging from 4 to 20 bp, because beyond 20 bp the transitions are ons why these methods performed very poorly in this benchmark. In
unlikely to be two-state (SantaLucia et al., 1996). Thus, it could be this respect, it is important to note that in our comparative assess-
argued that most of the oligonucleotides used in this study would not ment, the basic method and the NN model using the Th3 parameters
be following a two-state transition. However, in the same work, the (Sugimoto et al., 1996) exhibit very similar Tm values in the complete
authors suggested that the NN model can also provide reasonable oligonucleotide feature space (data not shown, see Supplementary
approximations for oligonucleotide sequences that do not have two- material). This means that under certain experimental conditions of
state transitions. Using the NN model, they obtained good predictions salt and oligonucleotide concentrations, a very simple method that

720
 Comparison of melting temperature calculation methods

Table 3. Statistical significance of differences in average Tm prediction errors among the methods

BAS SAL BRE SAN SUG CON

BAS — 9.1E – 47 2.8E – 23 7.3E – 83 1.3E – 112 1.6E – 107
SAL 2.9E – 46 — 2.0E – 20 3.0E – 53 5.2E – 67 1.1E – 56
BRE 6.7E – 26 1.0E – 10 — 3.2E – 73 1.2E – 75 7.6E – 71
SAN 2.7E – 99 7.9E – 68 3.1E – 85 — 0.050 0.009
SUG 5.1E – 129 2.1E – 92 3.0E – 81 0.479 — 0.1702
CON 5.0E – 122 8.3E – 75 4.0E – 83 0.477 0.432 —

Paired one-tail distribution Student’s t-tests were carried out to assess the statistical significance of average error or accuracy differences between any two methods, as described by

Downloaded from https://academic.oup.com/bioinformatics/article/21/6/711/199347 by guest on 14 March 2021

Press et al. (1997). The null hypothesis stated that there is no difference in the average error of melting temperatures predicted by any two methods (the observed average differences
could arise only by chance). The upper right triangle of the table (dark gray shaded cells) contains the obtained P -values for the pairwise statistical tests of all the methods based on
the 281 sequences used in the accuracy benchmark mapped in Zones 1, 2 and 3 of the consensus map shown in Figure 4. This set does not contain the sequences falling in the black
regions (zones of no consensus among any two methods), which were excluded from the analysis. The lower left triangle of the table (light gray shaded cells) contains the obtained
P -values for the same analysis mentioned above, but considering all 348 sequences used in the accuracy benchmark. The abbreviation of methods is as follows: BAS for the basic,
SAL for the salt adjusted, BRE for Th1 NN set (Breslauer et al., 1986), SAN for Th2 NN set (SantaLucia et al., 1996), SUG for Th3 NN set (Sugimoto et al., 1996) and CON for the
consensus method described in this study. For more details about the sequences used in the accuracy benchmark, see the legend of Table 2 or the section Methods.

does not take into account these parameters can give results similar two methods are highly similar, as shown by the comparative bench-
to the more complex methods that indeed considered these variables. mark results. Out of the 348 experimental data points used in this
However, under varying conditions of salt and oligonucleotide con- benchmark, a total of 282 sequences fall in Zone 1 (dark gray region,
centration, the NN model with proper thermodynamic parameters Fig. 4C). On the other hand, it was also shown that the poor per-
clearly outperforms the simple methods. The NN model using the formance in this benchmark by the Th1 set (Breslauer et al., 1986)
Th1 thermodynamic parameters (Breslauer et al., 1986) showed a did not degrade the overall performance of the consensus method in
very low performance in our accuracy benchmark, when compared Zone 2, although the Tm value obtained with the Th1 set is averaged
to the other parameter sets. This was in agreement with what was in this zone with values obtained from the Th3 set. This demon-
observed in previous works for shorter oligonucleotide sequences strates that the consensus Tm is a robust measure and validates the
(SantaLucia, 1998; Owczarzy et al., 1998). The accuracy of the NN usefulness of this comparative study. Although the accuracy differ-
model using the Th1 set was too low, that even the salt adjusted ences of Tm predictions obtained in this benchmark for Th3 and
method had a better performance, although highly unsatisfactory. consensus method were not statistically significant, we suggest that
The NN model using the Th2 (SantaLucia et al., 1996) and Th3 in a large-scale application, the consensus Tm will turn out to be
(Sugimoto et al., 1996) thermodynamic sets exhibited a good per- significantly more accurate than the Tm estimated by the Th3 NN
formance in this accuracy benchmark. This result suggests that both set alone. Unfortunately, not enough experimental data are available
sets could be successfully used to predict the experimental Tm of yet to perform such a large-scale accuracy benchmark. Thus, we
oligonucleotide sequences in the range of 16–30mers. However, the finally suggest that additional experimental data, covering a larger
best result was achieved by the consensus method described in this fraction of the oligonucleotide sequence space, are required to derive
study, which is based on the three thermodynamic sets and in the con- more accurate and robust thermodynamic parameters. Also, a large,
sensus map generated in the comparative assessment carried out in heterogeneous, representative and unbiased set of sequences should
this study. The consensus Tm value is a more robust measure, which be used to carry out a complete assessment of accuracy for the exist-
is less sensitive to large errors that could arise while using a single ing methods. The consensus method proposed in this study does not
parameter table. Irrespective of the method that is used, oligonuc- guarantee the best accuracy for any possible sequence. However, it
leotide sequences falling in regions where no consensus was observed minimizes the chances of error when using the existing methods for
are more prone to large errors in the experimental Tm estimation, as a diverse and large number of sequences, which is the case in the
it was demonstrated in our accuracy benchmark. Thus, we recom- currently used practical molecular biology applications.
mend the use of the consensus Tm value for sequences in the range
of 16–30mers, avoiding those sequences that fall in those regions
where no consensus was observed (black regions of Fig. 4A). CONCLUSIONS
The consensus Tm suggested in this work will minimize the error in Significant differences are observed for the Tm values of short
the long run. Owing to the Lack of enough experimental data covering DNA oligonucleotides calculated by different Tm prediction meth-
the complete oligonucleotide feature space under varying conditions ods. Additional experimental data covering a larger fraction of the
of salt and oligonucleotide concentration, this is the safest way to oligonucleotide feature space are required in order to evaluate the
proceed. The consensus Tm measure does not certainly guarantee accuracy of the current methods or to obtain a more precise estima-
the minimal error in all individual cases, but none of the methods tion of the experimental Tm for any short oligonucleotide sequence.
can do that either. It must be mentioned that the accuracy bench- Meanwhile, the use of a consensus Tm calculation with a minimal
mark proposed here favors those Tm values obtained with the Th2 error probability is suggested, which should be derived from the com-
(SantaLucia et al., 1996) and Th3 (Sugimoto et al., 1996) parameter parison of existing methods in a large benchmark set of sequences,
sets, because most of the sequences fall in the regions where these as it was the case in this study. The guidelines to follow in order

721
A.Panjkovich and F.Melo

to increase the success for practical molecular biology applications, Doktycz,M.J., Goldstein,R.F., Paner,T.M., Gallo,F.J. and Benight,A.S. (1992) Studies
from top to bottom priority, are as follows: (1) apply safely the cur- of DNA dumbbells. I. Melting curves of 17 DNA dumbbells with different duplex
stem sequences linked by T4 endloops: evaluation of the nearest-neighbor stacking
rent methods after considering the restrictions or limitations they
interactions in DNA. Biopolymers, 32, 849–864.
have (i.e. avoid sequences that form stable alternative secondary Gotoh,O. and Tagashira,Y. (1981) Stabilities of nearest-neighbor doublets in double-
structures, because such sequences do not follow a two-state trans- helical DNA determined by fitting calculated melting profiles to observed profiles.
ition); (2) if possible, use oligonucleotide sequences that fall in the Biopolymers, 20, 1033–1042.
middle range of CG-content and are shorter than 20–22mers (i.e. Howley,P.M., Israel,M.F., Law,M.-F. and Martin,M.A. (1979) A rapid method for detect-
ing and mapping homology between heterologous DNAs. Evaluation of polyomavirus
where most of the current Tm prediction methods agree); (3) avoid genomes. J. Biol. Chem., 254, 4876–4883.
the use of sequences that fall in those regions of oligonucleotide fea- Marmur,J. and Doty,P. (1962) Determination of the base composition of deoxyribo-
ture space where none of the current methods agrees (black regions nucleic acid from its thermal denaturation temperature. J. Mol. Biol., 5,
of Fig. 4A); (4) for large-scale applications with sequences where 109–118.
Owczarzy,R., Vallone,P.M., Gallo,F.J., Paner,T.M., Lane,M.J. and Benight,A.S. (1998)
a two-state transition is not known to occur, use a consensus Tm

Downloaded from https://academic.oup.com/bioinformatics/article/21/6/711/199347 by guest on 14 March 2021

Predicting sequence-dependent melting stability of short duplex DNA oligomers.
calculation method like the one suggested in this study, and then Biopolymers, 44, 217–239.
maximizes the rate of success and (5) refer to the upcoming literat- Owczarzy,R., You,Y., Moreira,B.G., Manthey,J.A., Huang,L., Behlke,M.A. and
ure, wherein new and improved methods for Tm prediction will be Walder,J.A. (2004) Effects of sodium ions on DNA duplex oligomers: improved
developed. predictions of melting temperatures. Biochemistry, 43, 3537–3554.
Pease,A.C., Solas,D., Sullivan,E.J., Cronin,M.T., Holmes,C.P. and Fodor,S.P.A. (1994)
Light-generated oligonucleotide arrays for rapid DNA sequence analysis. Proc. Natl
ACKNOWLEDGEMENTS Acad. Sci. USA, 91, 5022–5026.
Press,W.H., Teukolsky,S.A., Vetterling,W.T. and Flannery,B.P. (1997) Numerical
We gratefully acknowledge the helpful comments and suggestions
Recipes in C. The Art of Scientific Computing, 2nd edn. Cambridge University
made by the two anonymous reviewers of this manuscript. This Press, NY.
work was funded by grants from Fundación Andes (#13600/4), Rouzina,I. and Bloomfield,V.A. (1999) Heat capacity effects on the melting of DNA. 2.
FONDECYT (#1010959) and DIPUC (#2004/01PF). Analysis of nearest-neighbor base pair effects. Biophys. J., 77, 3252–3255.
SantaLucia,J.,Jr. (1998) A unified view of polymer, dumbbell, and oligonucleotide DNA
nearest-neighbor thermodynamics. Proc. Natl Acad. Sci. USA, 95, 1460–1465.
REFERENCES SantaLucia,J.,Jr, Allawi,H.T. and Seneviratne,P.A. (1996) Improved nearest-neighbor
Allawi,H.T. and SantaLucia,J.,Jr. (1997) Thermodynamics and NMR of internal G-T parameters for predicting DNA duplex stability. Biochemistry, 35, 3555–3562.
mismatches in DNA. Biochemistry, 36, 10581–10594. Schildkraut,C. and Lifson,S. (1965) Dependence of the melting temperature of DNA on
Borer,P.N., Dengler,B., Tinoco,I.,Jr. and Uhlenbeck,O.C. (1974) Stability of ribonucleic salt concentration. Biopolymers, 3, 195–208.
acid double-stranded helices. J. Mol. Biol., 86, 843–853. Steger,G. (1994) Thermal denaturation of double-stranded nucleic acids: prediction of
Breslauer,K.J., Frank,R., Blöcker,H. and Marky,L.A. (1986) Predicting DNA duplex temperatures critical for gradient gel electrophoresis and polymerase chain reaction.
stability from the base sequence. Proc. Natl Acad. Sci. USA, 83, 3746–3750. Nucleic Acids Res., 22, 2760–2768.
Buck,K.J., Harris,R.A. and Sikela,J.M. (1991) A general method for quantitative PCR Sugimoto,N., Nakano,S., Yoneyama,M. and Honda,K. (1996) Improved thermodynamic
analysis of mRNA levels for members of gene families: application to GABAA parameters and helix initiation factor to predict stability of DNA duplexes. Nucleic
receptor subunits. Biotechniques, 11, 636–641. Acids Res., 24, 4501–4505.
Chiu,W.L., Sze,C.N., Ma,N.T., Chiu,L.F., Leung,C.W. and Au-Yeung,S.C. (2003) Vologodskii,A.V., Amirikyan,B.R., Lyubchenko,Y.L. and Frank-Kamenetskii, M.D.
NTDB: Thermodynamic Database for Nucleic Acids, Version 2.0. Nucleic Acids (1984) Allowance for heterogeneous stacking in the DNA helix-coil transition theory.
Res., 31, 483–485. J. Biomol. Struct. Dyn., 2, 131–148.
Delcourt,S.G. and Blake,R.D. (1991) Stacking energies in DNA. J. Biol. Chem., 266, Wetmur,J.G. (1991) DNA probes: applications of the principles of nucleic acid
15160–15169. hybridization. Crit. Rev. Biochem. Mol. Biol., 26, 227–259.

722