Comparing different methods of calculating DNA melting temperature​
​by Prachi Kharade

Introduction

During replication, DNA strands separate and this loss of helical structure is referred to as denaturation or DNA melting. Heating the DNA disrupts the hydrogen bonds so
T
m
is the temperature at which half the DNA molecules are denatured. On a graph of absorbance of UV light vs temperature, the
T
m
would be the point at which the sigmoid graph changes concavity. As DNA denatures, the absorbance at 260nm increases. On the graph before the
T
m
point, the DNA is dsDNA (double stranded in helical structure) which has a lower absorbance and after the
T
m
point, the DNA is ssDNA (single stranded and no longer in helical structure) which has greater absorbance. The one big factor that affects the
T
m
is the CG content because G and C nucleotides have one more hydrogen bond than A and T nucleotides, so higher CG content can contribute to a higher
T
m
because more energy is needed to break the bonds. Other factors include ionic strength (total ion concentration in the DNA solution) and the length of the sequence. A higher ionic strength means there are more or stronger ionic interaction and therefore results in higher stability so more heat is required.
There are four common methods of calculating Melting Temperature (
T
m
) of DNA: Wallace equation, Marmur-Doty equation, Salt-adjusted method, and Nearest Neighbor method. Each method accounts for different factors that affect the melting temperature of the DNA strand. During the Wolfram Summer camp with the help of mentors and TAs, I used the Wolfram Language to implement these four methods and a condensed function to apply each of these methods on a given sequence. Then, the methods were applied to experimental data and compared to the experimental results to analyze the accuracy and effectiveness of each method.

Implementation of the Different Algorithms

I. Wallace Equation

The Wallace equation is the most basic method of calculating the
T
m
and works best for shorter oligonucleotide sequences. It only takes into account the content of each nucleotide base.
Method to count the number of times each base occurs in the strand:
In[]:=
CountBases[strand_]:=LetterCounts[strand]
Depending on the base, multiply it by a factor and add them together to get the
T
m
:
In[]:=
Wallace[strandS_]:=​​Module[{count,Tm},​​count=CountBases[strandS];​​Tm=(2*(Lookup[count,"A",0]+Lookup[count,"T",0])+4*(Lookup[count,"C",0]+Lookup[count,"G",0]));​​Return[Tm]​​]
An example of the Wallace Equation:
In[]:=
Wallace["ATGCCGAT"]
Out[]=
24

II. Marmur-Doty Equation

The Marmur-Doty equation is best used for longer sequences. Instead of taking the counts for each nucleotide base, it only looks at the percentage of CG bases present and the total length of the sequence.
Use the CountBases function and then calculate
T
m
using the formula:
In[]:=
MarmurDoty[strandD_]:=​​Module[{count,Tm},​​count=CountBases[strandD];​​Tm=69.3+(0.41*((Lookup[count,"C",0]+Lookup[count,"G",0])/StringLength[strandD])*100)-(650/StringLength[strandD]);​​Return[Tm]​​]
An example of the Marmur-Doty Equation:
In[]:=
MarmurDoty["ATCGATGCGATAAACCGATTAG"]
Out[]=
56.5273

III. Salt-Adjusted Method

The salt-adjusted method accounts for the ionic strength of the sequence so it takes into account the sodium ion content. The standard amount is 50 mM so that is what is used in this project. For different sequence length subgroups there are slightly different formulas used. The sequences are categorized into three groups, sequences under 14 base pairs (bps), sequences between 14 and 50 bps and sequences longer than 50 bps.
Check the length of the given strand and based on that, use a different formula:
In[]:=
SaltAdusted[StrandS_]:=​​Module[{count,Tm1,TmFinal},​​count=CountBases[StrandS];​​If[StringLength[StrandS]<14,​​Tm1=(2*(Lookup[count,"A",0]+Lookup[count,"T",0])+4*(Lookup[count,"C",0]+Lookup[count,"G",0]));​​TmFinal=Tm1-(16.6*Log[10,0.050])];​​If[(StringLength[StrandS]>13)&&(StringLength[StrandS]<50),​​TmFinal=100.5+(41*(Lookup[count,"G",0]+Lookup[count,"C",0])/(StringLength[StrandS]))-(820/(StringLength[StrandS]))+16.6*Log[10,0.05]];​​If[StringLength[StrandS]>50,TmFinal=81.5+(41*(Lookup[count,"G",0]+Lookup[count,"C",0])/(StringLength[StrandS]))-(500/(StringLength[StrandS]))+Log[10,0.05]-(0.62*0.05)];​​Return[TmFinal]​​]
An example of the Salt-Adjusted Equation:
In[]:=
SaltAdusted["ACGATCAGAGCTATTTACG"]
Out[]=
53.0082

IV. Nearest Neighbor Method

The nearest neighbor method utilizes DNA thermodynamic data for each pair of consecutive nucleotides in a sequence. Instead of looking at the content of each individual base, it looks at the relationship between neighboring bases. Then, based on what bases are next to each other, the entropy and enthalpy data are different. This is because the entropy and enthalpy numbers give more data about the energy changes which are different based on which bases are next to each other.
Import data table of thermodynamic data:
In[]:=
data=Import["https://www.sigmaaldrich.com/technical-documents/articles/biology/oligos-melting-temp.html","Data"];
*Different researchers have found slightly varying thermodynamic data. This data is from an online source but the origin of the data is from Breslauer et al. Other data sets include SantaLucia *
In[]:=
neighbors=data[[2,1,1,2,{3,4,5,6,7,8,9,10,11,12}]]
Out[]=
{{AA/TT*,-9.1**,-0.024,-6.6,-0.0184},{AT/TA,-8.6,-0.0239,-5.7,-0.0155},{TA/AT,-6.,-0.0169,-8.1,-0.0226},{CA/GT,-5.8,-0.0129,-10.5,-0.0278},{GT/CA,-6.5,-0.0173,-10.2,-0.0262},{CT/GA,-7.8,-0.0208,-7.6,-0.0192},{GA/CT,-5.6,-0.0135,-13.3,-0.0355},{CG/GC,-11.9,-0.0278,-8.,-0.0194},{GC/CG,-11.1,-0.0267,-14.2,-0.0349},{GG/CC,-11.,-0.0266,-12.2,-0.0297}}
After getting the thermodynamic data, the data was cleaned and formatted to only include entropy and enthalpy data for DNA:
In[]:=
neighbors[[1,1]]="AA/TT";
In[]:=
neighbors[[1,2]]=-9.1;
The original data table included the thermodynamic values for RNA and DNA and since this project is only working on DNA, drop the last two columns of data:
In[]:=
neighbor={Drop[neighbors[[1]],-2],Drop[neighbors[[2]],-2],Drop[neighbors[[3]],-2],Drop[neighbors[[4]],-2],​​Drop[neighbors[[5]],-2],Drop[neighbors[[6]],-2],Drop[neighbors[[7]],-2],Drop[neighbors[[8]],-2],​​Drop[neighbors[[9]],-2],Drop[neighbors[[10]],-2]}
Out[]=
{{AA/TT,-9.1,-0.024},{AT/TA,-8.6,-0.0239},{TA/AT,-6.,-0.0169},{CA/GT,-5.8,-0.0129},{GT/CA,-6.5,-0.0173},{CT/GA,-7.8,-0.0208},{GA/CT,-5.6,-0.0135},{CG/GC,-11.9,-0.0278},{GC/CG,-11.1,-0.0267},{GG/CC,-11.,-0.0266}}
The following two functions are helper methods that, given a neighbor pair, get the dH and dS values from the list.
Search the data set for the given base pair and return the enthalpy value from the data table:
In[]:=
getH[neighbor_,pair_]:=​​Module[{a,b,H,a2},​​b=0;​​a=Table[StringSplit[neighbor[[n,1]],"/"][[1]],{n,1,10}];​​a2=Table[StringSplit[neighbor[[n,1]],"/"][[2]],{n,1,10}];​​(If[StringMatchQ[a[[#]],pair]True,b=#])&/@Range[10];​​(*For[i=1,i≤10,i++,If[StringMatchQ[a[[i]],pair]True,b=i]];*)​​If[b0,(If[StringMatchQ[a2[[#]],pair]True,b=#])&/@Range[10]];​​(*If[b0,For[j=1,j≤10,j++,If[StringMatchQ[a2[[j]],pair]True,b=j]]];*)​​If[b0,(If[StringMatchQ[StringReverse[a2[[#]]],pair]True,b=#])&/@Range[10]];​​If[b0,(If[StringMatchQ[StringReverse[a[[#]]],pair]True,b=#])&/@Range[10]];​​neighbor[[b,2]]​​]
Search the data set for the given base pair and return the entropy value from the data table:
In[]:=
getS[neighbor_,pair_]:=​​Module[{a,b,S,a2},​​b=0;​​a=Table[StringSplit[neighbor[[n,1]],"/"][[1]],{n,1,10}];​​a2=Table[StringSplit[neighbor[[n,1]],"/"][[2]],{n,1,10}];​​(If[StringMatchQ[a[[#]],pair]True,b=#])&/@Range[10];​​(*For[i=1,i≤10,i++,If[StringMatchQ[a[[i]],pair]True,b=i]];*)​​If[b0,(If[StringMatchQ[a2[[#]],pair],b=#])&/@Range[10]];​​(*If[b0,For[j=1,j≤10,j++,If[StringMatchQ[a2[[j]],pair],b=j]]];*)​​If[b0,(If[StringMatchQ[StringReverse[a2[[#]]],pair]True,b=#])&/@Range[10]];​​If[b0,(If[StringMatchQ[StringReverse[a[[#]]],pair]True,b=#])&/@Range[10]];​​neighbor[[b,3]]​​]
Create a list of all the neighboring pairs in the given sequence:
In[]:=
neighborPairs[strandP_]:=Table[StringTake[strandP,{n,n+1}],{n,1,StringLength[strandP]-1}]
In[]:=
neighborPairs["TACTAAC"]
Out[]=
{TA,AC,CT,TA,AA,AC}
Go through the list of neighbor pairs for the sequence and calculate the dS and dH. After looking through many papers and websites, I have come across a number of slightly varying formulas for
T
m
. This is one of the formulas that I found to work best with the information I have. The following are parameters in the formula:
◼
  • dS is the sum of all the entropy change values for each neighbor pair and dH is the sum of all the enthalpy change values for each neighbor pair.
  • ◼
  • 0.00199 is R, the gas constant.
  • ◼
  • 0.0000005 is the oligonucleotide concentration.
  • ◼
  • The last term is for salt adjusting.
  • Get the sum of the entropy and enthalpy values for each base pair and apply the formula on the parameters:
    In[]:=
    NNTm[strandN_,neighbor_]:=​​Module[{Nps,Tm,dH,dS,L},​​Nps=neighborPairs[strandN];​​L=StringLength[strandN];​​dH=0;​​dS=0;​​For[i=1,i≤Length[Nps],i++,dH+=getH[neighbor,Nps[[i]]];dS+=getS[neighbor,Nps[[i]]]];​​Tm=((dH)/((dS)+(0.00199*Log[0.0000005/4])-0.0108))-273.15+(16.6*Log[10,0.050]);​​Tm​​]
    Example of the Nearest Neighbor method:

    Condensed Function

    Call the individual functions on the given sequence and based on the method parameter, return the result:

    Putting Together the Data Sets

    A table that puts together the results of the each of the methods and the experimental results:
    A table of all the distances between the methods and the experimental data. The first column is the differences between the experimental data and the Wallace results for each sequence and so on for each of the methods:

    Plotting the Graphs

    Create a grid of the histograms for each method which shows that each succeeding method gives increasingly more accurate results:

    Results & Conclusions

    Future Extensions

    References

    Salt-Adjusted Formulas https://labcalculator.net/wiki/oligo-tm​
    Nearest Neighbor Formula https://www.sigmaaldrich.com/technical-documents/articles/biology/oligos-melting-temp.html

    Keywords