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