[BioPython] Tm: First working version

[Sebastian Bassi]

Hello,

I've just finished a working version of the Tm function.
I contacted Nicolas Le Novere, author of "MELTING, computing the melting 
temperature of nucleic acid duplex" (Bioinformatics.17:1226-1227,2001) 
and help me understand some issues. He made a program (MELTING) that 
calculates Tm at various conditions, taking into account SantaLucia 
work. It performs better than EMBOSS and any other open source software. 
I dare to say the only better software is Santalucia's one, but it is 
commercial closed source.
What I am posting here is a working version of the function, IT IS NOT 
READY YET for submission into CVS, because I have to adjust some values 
to make it easy to use (like convenient units) and remove some stuff I 
used for developemnt that are not needed for daily use.

import string
import math
def Tm_staluc(s,dnac=50,saltc=50,mgc=15,rna=0):
     """Returns DNA/RNA tm using nearest neighbor thermodynamics. dnac is
     DNA/RNA concentration [nM] and saltc is salt concentration [mM]
     rna=0 is for DNA (default), for RNA, rna should be <>0.
     Overcount function thanks to Greg Singer <singerg at tcd.ie>, rest by
     Sebastian Bassi <sbassi at genesdigitales.com>"""
     dh=0 #DeltaH. Enthalpy
     ds=0 #deltaS Entropy

     def tercorr(stri):
         deltah=0
         deltas=0
         if stri[0]=="G" or stri[0]=="C":
             deltah=deltah-0.1
             deltas=deltas+2.8
         elif stri[0]=="A" or stri[0]=="T":
             deltah=deltah-2.3
             deltas=deltas-4.1
         if stri[-1]=="G" or stri[-1]=="C":
             deltah=deltah-0.1
             deltas=deltas+2.8
         elif stri[-1]=="A" or stri[-1]=="T":
             deltah=deltah-2.3
             deltas=deltas-4.1
         dhL=dh+deltah
	dsL=ds+deltas
         return dsL,dhL

     def saltcorr(s,saltc):
         saltc=(saltc/1E3)
         deltas=0
         deltas=deltas-0.368*(len(s)-1)*math.log(saltc)
         dsL=deltas
         return dsL #ver si esto va en el return o no!

     def overcount(st,p):
         """Returns how many p are on st, works even for overlapping"""
         ocu=0
         x=0
         while 1:
             try:
                 i=st.index(p,x)
             except ValueError:
                 break
             ocu=ocu+1
             x=i+1
         return ocu

     sup=string.upper(s)
     R=1.987 # universal gas constant in Cal/degrees C*Mol
     vsTC,vh=tercorr(sup)
     vs=vsTC

     k=(dnac/4.0) #* 1E-9 #; //converts from nanomolar to Molar. Note 
this ignores the contribution of the target since this is << than primer 
concentration.
     if rna==0:
 
vh=vh+(overcount(sup,"AA"))*7.9+(overcount(sup,"TT"))*7.9+(overcount(sup,"AT"))*7.2+(overcount(sup,"TA"))*7.2+(overcount(sup,"CA"))*8.5+(overcount(sup,"TG"))*8.5+(overcount(sup,"GT"))*8.4+(overcount(sup,"AC"))*8.4
 
vh=vh+(overcount(sup,"CT"))*7.8+(overcount(sup,"AG"))*7.8+(overcount(sup,"GA"))*8.2+(overcount(sup,"TC"))*8.2
 
vh=vh+(overcount(sup,"CG"))*10.6+(overcount(sup,"GC"))*10.6+(overcount(sup,"GG"))*8+(overcount(sup,"CC"))*8
 
vs=vs+(overcount(sup,"AA"))*22.2+(overcount(sup,"TT"))*22.2+(overcount(sup,"AT"))*20.4+(overcount(sup,"TA"))*21.3
 
vs=vs+(overcount(sup,"CA"))*22.7+(overcount(sup,"TG"))*22.7+(overcount(sup,"GT"))*22.4+(overcount(sup,"AC"))*22.4
 
vs=vs+(overcount(sup,"CT"))*21.0+(overcount(sup,"AG"))*21.0+(overcount(sup,"GA"))*22.2+(overcount(sup,"TC"))*22.2
 
vs=vs+(overcount(sup,"CG"))*27.2+(overcount(sup,"GC"))*27.2+(overcount(sup,"GG"))*19.9+(overcount(sup,"CC"))*19.9
         ds=vs
         dh=vh

     else:
 
vh=(overcount(sup,"AA"))*6.6+(overcount(sup,"TT"))*6.6+(overcount(sup,"AT"))*5.7+(overcount(sup,"TA"))*8.1+(overcount(sup,"CA"))*10.5+(overcount(sup,"TG"))*10.5+(overcount(sup,"GT"))*10.2+(overcount(sup,"AC"))*10.2
 
vh=vh+(overcount(sup,"CT"))*7.6+(overcount(sup,"AG"))*7.6+(overcount(sup,"GA"))*13.3+(overcount(sup,"TC"))*13.3
 
vh=vh+(overcount(sup,"CG"))*8.0+(overcount(sup,"GC"))*14.2+(overcount(sup,"GG"))*12.2+(overcount(sup,"CC"))*12.2
 
vs=(overcount(sup,"AA"))*18.4+(overcount(sup,"TT"))*18.4+(overcount(sup,"AT"))*15.5+(overcount(sup,"TA"))*16.9
 
vs=vs+(overcount(sup,"CA"))*27.8+(overcount(sup,"TG"))*27.8+(overcount(sup,"GT"))*26.2+(overcount(sup,"AC"))*26.2
 
vs=vs+(overcount(sup,"CT"))*19.2+(overcount(sup,"AG"))*19.2+(overcount(sup,"GA"))*35.5+(overcount(sup,"TC"))*35.5
 
vs=vs+(overcount(sup,"CG"))*19.4+(overcount(sup,"GC"))*34.9+(overcount(sup,"GG"))*29.7+(overcount(sup,"CC"))*29.7
         ds=vs
         dh=vh

     ds=ds-0.368*(len(s)-1)*math.log(saltc)
     tm=((1000* (-dh))/(-ds+(R * (math.log(k)))))-273.15
     return tm


-- 
Best regards,

//=\ Sebastian Bassi - Diplomado en Ciencia y Tecnologia, UNQ   //=\
\=// IT Manager Advanta Seeds - Balcarce Research Center -      \=//
//=\ Pro secretario ASALUP - www.asalup.org - PGP key available //=\
\=// E-mail: sbassi at genesdigitales.com - ICQ UIN: 3356556 -     \=//

                 http://Bioinformatica.info