Package Bio :: Package SeqUtils :: Module MeltingTemp
[hide private]
[frames] | no frames]

Source Code for Module Bio.SeqUtils.MeltingTemp

  1  # Copyright 2004-2008 by Sebastian Bassi. 
  2  # All rights reserved. 
  3  # This code is part of the Biopython distribution and governed by its 
  4  # license.  Please see the LICENSE file that should have been included 
  5  # as part of this package. 
  6   
  7  """Calculate the thermodynamic melting temperatures of nucleotide sequences.""" 
  8   
  9  from __future__ import print_function 
 10   
 11  import math 
 12   
 13   
14 -def Tm_staluc(s, dnac=50, saltc=50, rna=0):
15 """Returns DNA/DNA tm using nearest neighbor thermodynamics. 16 17 dnac is DNA concentration [nM] 18 saltc is salt concentration [mM]. 19 rna=0 is for DNA/DNA (default), use 1 for RNA/RNA hybridisation. 20 21 For DNA/DNA, see Allawi & SantaLucia (1997), Biochemistry 36: 10581-10594 22 For RNA/RNA, see Xia et al (1998), Biochemistry 37: 14719-14735 23 24 Example: 25 26 >>> print("%0.2f" % Tm_staluc('CAGTCAGTACGTACGTGTACTGCCGTA')) 27 59.87 28 >>> print("%0.2f" % Tm_staluc('CAGTCAGTACGTACGTGTACTGCCGTA', rna=True)) 29 68.14 30 31 You can also use a Seq object instead of a string, 32 33 >>> from Bio.Seq import Seq 34 >>> from Bio.Alphabet import generic_nucleotide 35 >>> s = Seq('CAGTCAGTACGTACGTGTACTGCCGTA', generic_nucleotide) 36 >>> print("%0.2f" % Tm_staluc(s)) 37 59.87 38 >>> print("%0.2f" % Tm_staluc(s, rna=True)) 39 68.14 40 41 """ 42 43 #Credits: 44 #Main author: Sebastian Bassi <sbassi@genesdigitales.com> 45 #Overcount function: Greg Singer <singerg@tcd.ie> 46 #Based on the work of Nicolas Le Novere <lenov@ebi.ac.uk> Bioinformatics. 47 #17:1226-1227(2001) 48 49 #This function returns better results than EMBOSS DAN because it uses 50 #updated thermodynamics values and takes into account inicialization 51 #parameters from the work of SantaLucia (1998). 52 53 #Things to do: 54 #+Detect complementary sequences. Change K according to result. 55 #+Add support for heteroduplex (see Sugimoto et al. 1995). 56 #+Correction for Mg2+. Now supports only monovalent ions. 57 #+Put thermodinamics table in a external file for users to change at will 58 #+Add support for danglings ends (see Le Novele. 2001) and mismatches. 59 60 dh = 0 # DeltaH. Enthalpy 61 ds = 0 # deltaS Entropy 62 63 def tercorr(stri): 64 deltah = 0 65 deltas = 0 66 if rna == 0: 67 #DNA/DNA 68 #Allawi and SantaLucia (1997). Biochemistry 36 : 10581-10594 69 if stri.startswith('G') or stri.startswith('C'): 70 deltah -= 0.1 71 deltas += 2.8 72 elif stri.startswith('A') or stri.startswith('T'): 73 deltah -= 2.3 74 deltas -= 4.1 75 if stri.endswith('G') or stri.endswith('C'): 76 deltah -= 0.1 77 deltas += 2.8 78 elif stri.endswith('A') or stri.endswith('T'): 79 deltah -= 2.3 80 deltas -= 4.1 81 dhL = dh + deltah 82 dsL = ds + deltas 83 return dsL, dhL 84 elif rna == 1: 85 #RNA 86 if stri.startswith('G') or stri.startswith('C'): 87 deltah -= 3.61 88 deltas -= 1.5 89 elif stri.startswith('A') or stri.startswith('T') or \ 90 stri.startswith('U'): 91 deltah -= 3.72 92 deltas += 10.5 93 if stri.endswith('G') or stri.endswith('C'): 94 deltah -= 3.61 95 deltas -= 1.5 96 elif stri.endswith('A') or stri.endswith('T') or \ 97 stri.endswith('U'): 98 deltah -= 3.72 99 deltas += 10.5 100 dhL = dh + deltah 101 dsL = ds + deltas 102 # print("delta h=%f" % dhL) 103 return dsL, dhL 104 else: 105 raise ValueError("rna = %r not supported" % rna)
106 107 def overcount(st, p): 108 """Returns how many p are on st, works even for overlapping""" 109 ocu = 0 110 x = 0 111 while True: 112 try: 113 i = st.index(p, x) 114 except ValueError: 115 break 116 ocu += 1 117 x = i + 1 118 return ocu 119 120 R = 1.987 # universal gas constant in Cal/degrees C*Mol 121 sup = str(s).upper() # turn any Seq object into a string (need index method) 122 vsTC, vh = tercorr(sup) 123 vs = vsTC 124 125 k = (dnac/4.0)*1e-9 126 #With complementary check on, the 4.0 should be changed to a variable. 127 128 if rna == 0: 129 #DNA/DNA 130 #Allawi and SantaLucia (1997). Biochemistry 36 : 10581-10594 131 vh = vh + (overcount(sup, "AA"))*7.9 + (overcount(sup, "TT"))*\ 132 7.9 + (overcount(sup, "AT"))*7.2 + (overcount(sup, "TA"))*7.2 \ 133 + (overcount(sup, "CA"))*8.5 + (overcount(sup, "TG"))*8.5 + \ 134 (overcount(sup, "GT"))*8.4 + (overcount(sup, "AC"))*8.4 135 vh = vh + (overcount(sup, "CT"))*7.8+(overcount(sup, "AG"))*\ 136 7.8 + (overcount(sup, "GA"))*8.2 + (overcount(sup, "TC"))*8.2 137 vh = vh + (overcount(sup, "CG"))*10.6+(overcount(sup, "GC"))*\ 138 9.8 + (overcount(sup, "GG"))*8 + (overcount(sup, "CC"))*8 139 vs = vs + (overcount(sup, "AA"))*22.2+(overcount(sup, "TT"))*\ 140 22.2 + (overcount(sup, "AT"))*20.4 + (overcount(sup, "TA"))*21.3 141 vs = vs + (overcount(sup, "CA"))*22.7+(overcount(sup, "TG"))*\ 142 22.7 + (overcount(sup, "GT"))*22.4 + (overcount(sup, "AC"))*22.4 143 vs = vs + (overcount(sup, "CT"))*21.0+(overcount(sup, "AG"))*\ 144 21.0 + (overcount(sup, "GA"))*22.2 + (overcount(sup, "TC"))*22.2 145 vs = vs + (overcount(sup, "CG"))*27.2+(overcount(sup, "GC"))*\ 146 24.4 + (overcount(sup, "GG"))*19.9 + (overcount(sup, "CC"))*19.9 147 ds = vs 148 dh = vh 149 elif rna == 1: 150 #RNA/RNA hybridisation of Xia et al (1998) 151 #Biochemistry 37: 14719-14735 152 vh = vh+(overcount(sup, "AA"))*6.82+(overcount(sup, "TT"))*6.6+\ 153 (overcount(sup, "AT"))*9.38 + (overcount(sup, "TA"))*7.69+\ 154 (overcount(sup, "CA"))*10.44 + (overcount(sup, "TG"))*10.5+\ 155 (overcount(sup, "GT"))*11.4 + (overcount(sup, "AC"))*10.2 156 vh = vh + (overcount(sup, "CT"))*10.48 + (overcount(sup, "AG"))\ 157 *7.6+(overcount(sup, "GA"))*12.44+(overcount(sup, "TC"))*13.3 158 vh = vh + (overcount(sup, "CG"))*10.64 + (overcount(sup, "GC"))\ 159 *14.88+(overcount(sup, "GG"))*13.39+(overcount(sup, "CC"))*12.2 160 vs = vs + (overcount(sup, "AA"))*19.0 + (overcount(sup, "TT"))*\ 161 18.4+(overcount(sup, "AT"))*26.7+(overcount(sup, "TA"))*20.5 162 vs = vs + (overcount(sup, "CA"))*26.9 + (overcount(sup, "TG"))*\ 163 27.8 + (overcount(sup, "GT"))*29.5 + (overcount(sup, "AC"))*26.2 164 vs = vs + (overcount(sup, "CT"))*27.1 + (overcount(sup, "AG"))*\ 165 19.2 + (overcount(sup, "GA"))*32.5 + (overcount(sup, "TC"))*35.5 166 vs = vs + (overcount(sup, "CG"))*26.7 + (overcount(sup, "GC"))\ 167 *36.9 + (overcount(sup, "GG"))*32.7 + (overcount(sup, "CC"))*29.7 168 ds = vs 169 dh = vh 170 else: 171 raise ValueError("rna = %r not supported" %rna) 172 173 ds = ds-0.368*(len(s)-1)*math.log(saltc/1e3) 174 tm = ((1000* (-dh))/(-ds+(R * (math.log(k)))))-273.15 175 # print("ds=%f" % ds) 176 # print("dh=%f" % dh) 177 return tm 178 179
180 -def _test():
181 """Run the module's doctests (PRIVATE).""" 182 import doctest 183 print("Running doctests...") 184 doctest.testmod() 185 print("Done")
186 187 if __name__ == "__main__": 188 _test() 189