/usr/share/pyshared/Bio/SeqUtils/MeltingTemp.py is in python-biopython 1.58-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 | # Copyright 2004-2008 by Sebastian Bassi.
# All rights reserved.
# This code is part of the Biopython distribution and governed by its
# license. Please see the LICENSE file that should have been included
# as part of this package.
"""Calculate the thermodynamic melting temperatures of nucleotide sequences."""
import math
def Tm_staluc(s,dnac=50,saltc=50,rna=0):
"""Returns DNA/DNA tm using nearest neighbor thermodynamics.
dnac is DNA concentration [nM]
saltc is salt concentration [mM].
rna=0 is for DNA/DNA (default), for RNA, rna should be 1.
Sebastian Bassi <sbassi@genesdigitales.com>"""
#Credits:
#Main author: Sebastian Bassi <sbassi@genesdigitales.com>
#Overcount function: Greg Singer <singerg@tcd.ie>
#Based on the work of Nicolas Le Novere <lenov@ebi.ac.uk> Bioinformatics.
#17:1226-1227(2001)
#This function returns better results than EMBOSS DAN because it uses
#updated thermodynamics values and takes into account inicialization
#parameters from the work of SantaLucia (1998).
#Things to do:
#+Detect complementary sequences. Change K according to result.
#+Add support for heteroduplex (see Sugimoto et al. 1995).
#+Correction for Mg2+. Now supports only monovalent ions.
#+Put thermodinamics table in a external file for users to change at will
#+Add support for danglings ends (see Le Novele. 2001) and mismatches.
dh = 0 #DeltaH. Enthalpy
ds = 0 #deltaS Entropy
def tercorr(stri):
deltah = 0
deltas = 0
if rna==0:
#DNA/DNA
#Allawi and SantaLucia (1997). Biochemistry 36 : 10581-10594
if stri.startswith('G') or stri.startswith('C'):
deltah -= 0.1
deltas += 2.8
elif stri.startswith('A') or stri.startswith('T'):
deltah -= 2.3
deltas -= 4.1
if stri.endswith('G') or stri.endswith('C'):
deltah -= 0.1
deltas += 2.8
elif stri.endswith('A') or stri.endswith('T'):
deltah -= 2.3
deltas -= 4.1
dhL = dh + deltah
dsL = ds + deltas
return dsL,dhL
elif rna==1:
#RNA
if stri.startswith('G') or stri.startswith('C'):
deltah -= 3.61
deltas -= 1.5
elif stri.startswith('A') or stri.startswith('T') or \
stri.startswith('U'):
deltah -= 3.72
deltas += 10.5
if stri.endswith('G') or stri.endswith('C'):
deltah -= 3.61
deltas -= 1.5
elif stri.endswith('A') or stri.endswith('T') or \
stri.endswith('U'):
deltah -= 3.72
deltas += 10.5
dhL = dh + deltah
dsL = ds + deltas
# print "delta h=",dhL
return dsL,dhL
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 += 1
x = i + 1
return ocu
R = 1.987 # universal gas constant in Cal/degrees C*Mol
sup = s.upper()
vsTC,vh = tercorr(sup)
vs = vsTC
k = (dnac/4.0)*1e-9
#With complementary check on, the 4.0 should be changed to a variable.
if rna==0:
#DNA/DNA
#Allawi and SantaLucia (1997). Biochemistry 36 : 10581-10594
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"))*\
9.8 + (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"))*\
24.4 + (overcount(sup,"GG"))*19.9 + (overcount(sup,"CC"))*19.9
ds = vs
dh = vh
else:
#RNA/RNA hybridisation of Xia et al (1998)
#Biochemistry 37: 14719-14735
vh = vh+(overcount(sup,"AA"))*6.82+(overcount(sup,"TT"))*6.6+\
(overcount(sup,"AT"))*9.38 + (overcount(sup,"TA"))*7.69+\
(overcount(sup,"CA"))*10.44 + (overcount(sup,"TG"))*10.5+\
(overcount(sup,"GT"))*11.4 + (overcount(sup,"AC"))*10.2
vh = vh + (overcount(sup,"CT"))*10.48 + (overcount(sup,"AG"))\
*7.6+(overcount(sup,"GA"))*12.44+(overcount(sup,"TC"))*13.3
vh = vh + (overcount(sup,"CG"))*10.64 + (overcount(sup,"GC"))\
*14.88+(overcount(sup,"GG"))*13.39+(overcount(sup,"CC"))*12.2
vs = vs + (overcount(sup,"AA"))*19.0 + (overcount(sup,"TT"))*\
18.4+(overcount(sup,"AT"))*26.7+(overcount(sup,"TA"))*20.5
vs = vs + (overcount(sup,"CA"))*26.9 + (overcount(sup,"TG"))*\
27.8 + (overcount(sup,"GT"))*29.5 + (overcount(sup,"AC"))*26.2
vs = vs + (overcount(sup,"CT"))*27.1 + (overcount(sup,"AG"))*\
19.2 + (overcount(sup,"GA"))*32.5 + (overcount(sup,"TC"))*35.5
vs = vs + (overcount(sup,"CG"))*26.7 + (overcount(sup,"GC"))\
*36.9 + (overcount(sup,"GG"))*32.7 + (overcount(sup,"CC"))*29.7
ds = vs
dh = vh
ds = ds-0.368*(len(s)-1)*math.log(saltc/1e3)
tm = ((1000* (-dh))/(-ds+(R * (math.log(k)))))-273.15
# print "ds="+str(ds)
# print "dh="+str(dh)
return tm
if __name__ == "__main__":
print "Quick self test"
assert Tm_staluc('CAGTCAGTACGTACGTGTACTGCCGTA') == 59.865612727457972
assert Tm_staluc('CAGTCAGTACGTACGTGTACTGCCGTA',rna=1) == 68.141611264576682
print "Done"
|