[BioPython] Tm: First working version
Sebastian Bassi
sbassi at asalup.org
Tue Feb 24 09:53:46 EST 2004
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
More information about the BioPython
mailing list